Pyomo

Pyomo

PyomoPython線形計画法整数計画法などの数理最適化問題をモデル化、解決するための強力なツールです。

以下はPyomoの便利な使い方をいくつか紹介します。

1. 簡単な線形計画問題のモデル化

まず、Pyomoを使って簡単な線形計画問題をモデル化してみましょう。

問題

次の線形計画問題を解いてみます。
最大化:
$$
\text{maximize} \quad 2x + 3y
$$
制約:
$$
x + 2y \leq 20
$$
$$
3x + y \leq 15
$$
$$
x \geq 0, \quad y \geq 0
$$

解法

まず、Pyomoをインストールします。

1
2
!pip install pyomo
!apt-get install -y -qq glpk-utils

次に、Pyomoのコードを書きます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
from pyomo.environ import *

# モデルの作成
model = ConcreteModel()

# 変数の定義
model.x = Var(within=NonNegativeReals)
model.y = Var(within=NonNegativeReals)

# 目的関数の定義
model.objective = Objective(expr=2*model.x + 3*model.y, sense=maximize)

# 制約条件の定義
model.constraint1 = Constraint(expr=model.x + 2*model.y <= 20)
model.constraint2 = Constraint(expr=3*model.x + model.y <= 15)

# 求解器の呼び出し
solver = SolverFactory('glpk')
solver.solve(model)

# 結果の表示
print("Optimal value for x:", model.x())
print("Optimal value for y:", model.y())
print("Optimal objective value:", model.objective())

[実行結果]

Optimal value for x: 2.0
Optimal value for y: 9.0
Optimal objective value: 31.0

2. 制約付き最適化問題のモデル化

Pyomoは制約の数が多くても簡単にモデル化できます。

例えば、多くの変数や複雑な制約がある場合でも、Pythonのリストや辞書を使って効率的に定義できます。

問題

多変数の制約付き最適化問題を解く例です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
from pyomo.environ import *

# モデルの作成
model = ConcreteModel()

# データの定義
N = 5 # 変数の数
coefficients = [1, 2, 3, 4, 5]
constraints = [
[2, 1, 2, 1, 1, 10],
[1, 3, 2, 2, 1, 15]
]

# 変数の定義
model.x = Var(range(N), within=NonNegativeReals)

# 目的関数の定義
model.objective = Objective(expr=sum(coefficients[i] * model.x[i] for i in range(N)), sense=maximize)

# 制約条件の定義
model.constraints = ConstraintList()
for con in constraints:
model.constraints.add(sum(con[i] * model.x[i] for i in range(N)) <= con[-1])

# 求解器の呼び出し
solver = SolverFactory('glpk')
solver.solve(model)

# 結果の表示
for i in range(N):
print(f"Optimal value for x[{i}]:", model.x[i]())
print("Optimal objective value:", model.objective())

[実行結果]

Optimal value for x[0]: 0.0
Optimal value for x[1]: 0.0
Optimal value for x[2]: 0.0
Optimal value for x[3]: 0.0
Optimal value for x[4]: 10.0
Optimal objective value: 50.0

3. 整数計画問題のモデル化

Pyomo整数計画問題も扱えます。

例えば、整数変数を使ったモデルの作成です。

問題

整数計画問題の例です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
from pyomo.environ import *

# モデルの作成
model = ConcreteModel()

# 変数の定義
model.x = Var(within=NonNegativeIntegers)
model.y = Var(within=NonNegativeIntegers)

# 目的関数の定義
model.objective = Objective(expr=2*model.x + 3*model.y, sense=maximize)

# 制約条件の定義
model.constraint1 = Constraint(expr=model.x + 2*model.y <= 20)
model.constraint2 = Constraint(expr=3*model.x + model.y <= 15)

# 求解器の呼び出し
solver = SolverFactory('glpk')
solver.solve(model)

# 結果の表示
print("Optimal value for x:", model.x())
print("Optimal value for y:", model.y())
print("Optimal objective value:", model.objective())

[実行結果]

Optimal value for x: 2.0
Optimal value for y: 9.0
Optimal objective value: 31.0

4. 複雑な制約条件のモデル化

Pyomoを使えば、複雑な非線形制約条件もモデル化できます。

問題

非線形制約条件を含む最適化問題の例です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
!apt-get update
!apt-get install -y coinor-cbc
!pip install -q pyomo

from pyomo.environ import *

# モデルの作成
model = ConcreteModel()

# データの定義
N = 5 # 変数の数
coefficients = [1, 2, 3, 4, 5]
constraints = [
[2, 1, 2, 1, 1, 10],
[1, 3, 2, 2, 1, 15]
]

# 変数の定義
model.x = Var(range(N), within=NonNegativeReals)

# 目的関数の定義
model.objective = Objective(expr=sum(coefficients[i] * model.x[i] for i in range(N)), sense=maximize)

# 制約条件の定義
model.constraints = ConstraintList()
for con in constraints:
model.constraints.add(sum(con[i] * model.x[i] for i in range(N)) <= con[-1])

# CBCソルバーを使用する場合
solver = SolverFactory('cbc', executable='/usr/bin/cbc')
result = solver.solve(model)

# 結果の表示
print("CBC solver results:")
for i in range(N):
print(f"Optimal value for x[{i}]:", model.x[i]())
print("Optimal objective value:", model.objective())

[実行結果]

CBC solver results:
Optimal value for x[0]: 0.0
Optimal value for x[1]: 0.0
Optimal value for x[2]: 0.0
Optimal value for x[3]: 0.0
Optimal value for x[4]: 10.0
Optimal objective value: 50.0

Pyomoは非常に柔軟で強力なツールであり、さまざまな最適化問題に対応できます。

これらの例を参考に、実際の問題に応用してみてください。

scikit-learn

scikit-learn

scikit-learnは、Python機械学習モデル構築評価実行するためのライブラリです。

以下に、scikit-learnの便利な使い方をいくつか紹介します。

1. 基本的なデータの読み込みと前処理

1
2
3
4
5
6
7
8
9
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split

# データの読み込み
iris = load_iris()
X, y = iris.data, iris.target

# トレーニングセットとテストセットに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

2. モデルの構築とトレーニング

1
2
3
4
5
6
7
from sklearn.ensemble import RandomForestClassifier

# モデルの構築
model = RandomForestClassifier(n_estimators=100, random_state=42)

# モデルのトレーニング
model.fit(X_train, y_train)

3. モデルの評価

1
2
3
4
5
6
7
8
9
10
11
12
from sklearn.metrics import accuracy_score, classification_report

# 予測
y_pred = model.predict(X_test)

# 精度の評価
accuracy = accuracy_score(y_test, y_pred)
print(f'Accuracy: {accuracy}')

# 詳細な評価
report = classification_report(y_test, y_pred)
print(report)

[実行結果]

Accuracy: 1.0
              precision    recall  f1-score   support

           0       1.00      1.00      1.00        10
           1       1.00      1.00      1.00         9
           2       1.00      1.00      1.00        11

    accuracy                           1.00        30
   macro avg       1.00      1.00      1.00        30
weighted avg       1.00      1.00      1.00        30

4. クロスバリデーション

1
2
3
4
5
6
7
8
9
from sklearn.model_selection import cross_val_score

# クロスバリデーションの実行
scores = cross_val_score(model, X, y, cv=5)

# 平均スコアと標準偏差
print(f'Cross-validation scores: {scores}')
print(f'Mean score: {scores.mean()}')
print(f'Standard deviation: {scores.std()}')

[実行結果]

Cross-validation scores: [0.96666667 0.96666667 0.93333333 0.96666667 1.        ]
Mean score: 0.9666666666666668
Standard deviation: 0.02108185106778919

5. ハイパーパラメータチューニング

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
from sklearn.model_selection import GridSearchCV

# ハイパーパラメータの候補を設定
param_grid = {
'n_estimators': [50, 100, 200],
'max_depth': [None, 10, 20, 30]
}

# グリッドサーチの実行
grid_search = GridSearchCV(model, param_grid, cv=5)
grid_search.fit(X_train, y_train)

# 最良のパラメータとスコア
best_params = grid_search.best_params_
best_score = grid_search.best_score_

print(f'Best parameters: {best_params}')
print(f'Best cross-validation score: {best_score}')

[実行結果]

Best parameters: {'max_depth': None, 'n_estimators': 200}
Best cross-validation score: 0.95

6. パイプラインの作成

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

# パイプラインの構築
pipeline = Pipeline([
('scaler', StandardScaler()),
('classifier', RandomForestClassifier(n_estimators=100, random_state=42))
])

# パイプラインでのトレーニングと評価
pipeline.fit(X_train, y_train)
y_pred = pipeline.predict(X_test)

accuracy = accuracy_score(y_test, y_pred)
print(f'Accuracy: {accuracy}')

[実行結果]

Accuracy: 1.0

7. 主成分分析 (PCA)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from sklearn.decomposition import PCA

# PCAの適用
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)

# トレーニングセットとテストセットに分割
X_train, X_test, y_train, y_test = train_test_split(X_pca, y, test_size=0.2, random_state=42)

# モデルのトレーニングと評価
model.fit(X_train, y_train)
y_pred = model.predict(X_test)

accuracy = accuracy_score(y_test, y_pred)
print(f'Accuracy with PCA: {accuracy}')

[実行結果]

Accuracy with PCA: 1.0

scikit-learnは、シンプルなAPIで多くの機械学習タスクを簡単に実行できるように設計されています。

これらの基本的な使い方を理解することで、さらに高度な機械学習実装実験に役立てることができます。

SciPy

SciPy

SciPy科学計算エンジニアリングの分野で広く利用されているPythonライブラリで、便利な機能が多数揃っています。

ここでは、SciPyの便利な使い方をいくつか紹介します。

1. 数値積分

SciPyintegrateモジュールを使って数値積分を行うことができます。

以下は、関数の定積分を計算する例です。

1
2
3
4
5
6
7
8
9
import scipy.integrate as integrate

# 定義する関数
def f(x):
return x ** 2

# 数値積分
result, error = integrate.quad(f, 0, 1)
print("積分の結果:", result)

[実行結果]

積分の結果: 0.33333333333333337

2. 最適化

SciPyoptimizeモジュールを使って関数の最小化方程式の解を求めることができます。

関数の最小化

1
2
3
4
5
6
7
8
9
10
import numpy as np
import scipy.optimize as optimize

# 定義する関数
def f(x):
return x ** 2 + 10 * np.sin(x)

# 最小化
result = optimize.minimize(f, x0=0)
print("最適化の結果:", result.x)

[実行結果]

最適化の結果: [-1.30644012]

方程式の解を求める

1
2
3
4
5
6
7
# 定義する方程式
def equation(x):
return x ** 2 - 2

# 解を求める
solution = optimize.root(equation, x0=1)
print("方程式の解:", solution.x)

[実行結果]

方程式の解: [1.41421356]

3. 線形代数

SciPylinalgモジュールを使って行列演算を行うことができます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
import scipy.linalg as linalg
import numpy as np

# 行列の定義
A = np.array([[3, 2], [1, 4]])
B = np.array([1, 2])

# 線形方程式の解を求める
x = linalg.solve(A, B)
print("方程式の解:", x)

# 行列の逆行列を求める
A_inv = linalg.inv(A)
print("逆行列:", A_inv)

[実行結果]

方程式の解: [0.  0.5]
逆行列: [[ 0.4 -0.2]
 [-0.1  0.3]]

4. 信号処理

SciPysignalモジュールを使ってフィルタリング周波数解析を行うことができます。

信号のフィルタリング

1
2
3
4
5
6
7
8
9
10
11
12
from scipy.signal import butter, filtfilt
import numpy as np

# バターワースフィルタの設計
b, a = butter(N=4, Wn=0.2)

# 信号の定義(繰り返し)
signal = np.tile([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 3) # 繰り返して長さを増やす

# フィルタリング
filtered_signal = filtfilt(b, a, signal)
print("フィルタリング後の信号:", filtered_signal)

[実行結果]

フィルタリング後の信号: [ 0.97688963  2.12759398  3.32584065  4.55768608  5.73180471  6.68699207
  7.23854322  7.25256462  6.72119364  5.80208259  4.79233683  4.03518382
  3.79495244  4.15630622  4.99146245  6.00594272  6.84255021  7.20634692
  6.96867975  6.2130705   5.20271618  4.28015452  3.74274467  3.74971619
  4.29610177  5.25159104  6.43442239  7.6839036   8.90514189 10.07593652]

5. 統計

SciPystatsモジュールを使って統計解析を行うことができます。

基本的な統計量の計算

1
2
3
4
5
6
7
8
9
10
11
12
13
from scipy import stats

# データの定義
data = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

# 平均と標準偏差
mean = np.mean(data)
std_dev = np.std(data)
print("平均:", mean, "標準偏差:", std_dev)

# 正規分布のフィッティング
param = stats.norm.fit(data)
print("正規分布のパラメータ:", param)

[実行結果]

平均: 5.5 標準偏差: 2.8722813232690143
正規分布のパラメータ: (5.5, 2.8722813232690143)

6. 補間

SciPyinterpolateモジュールを使って補間を行うことができます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

# データの定義
x = np.linspace(0, 10, 10)
y = np.sin(x)

# 線形補間
f = interp1d(x, y, kind='linear')

# 新しいデータ点の生成
x_new = np.linspace(0, 10, 100)
y_new = f(x_new)

# プロット
plt.plot(x, y, 'o', label='Original data')
plt.plot(x_new, y_new, '-', label='Interpolated data')
plt.legend()
plt.show()

[実行結果]


これらはSciPyの便利な使い方のほんの一部です。

SciPyは多機能で多岐にわたる分野で使用できるため、目的に応じて必要なモジュールや関数を活用してみてください。

PuLP

PuLP

PuLPPython用の線形プログラミングライブラリで、最適化問題を簡単に定義し解くことができます。

以下にPuLPの基本的な使い方と便利な機能を紹介します。

基本的な使い方

1. インストール

まず、PuLPをインストールします。

1
pip install pulp

2. 問題の定義

次に、線形プログラミング問題を定義します。

例えば、以下は最大化問題の定義です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import pulp

# 問題の定義
model = pulp.LpProblem("Maximize_Profit", pulp.LpMaximize)

# 変数の定義
x = pulp.LpVariable("x", lowBound=0, cat='Continuous')
y = pulp.LpVariable("y", lowBound=0, cat='Continuous')

# 目的関数の定義
model += 40 * x + 30 * y, "Total Profit"

# 制約条件の定義
model += 2 * x + y <= 20, "Constraint 1"
model += 4 * x + 3 * y <= 45, "Constraint 2"

# 問題を解く
model.solve()

# 結果の表示
print(f"Status: {pulp.LpStatus[model.status]}")
print(f"x: {pulp.value(x)}")
print(f"y: {pulp.value(y)}")
print(f"Total Profit: {pulp.value(model.objective)}")

[実行結果]

Status: Optimal
x: 0.0
y: 15.0
Total Profit: 450.0

応用的な使い方

3. 整数計画問題

PuLPでは整数計画問題も扱うことができます。

以下は整数計画問題の例です。

1
2
3
4
5
# 変数の定義 (整数)
x = pulp.LpVariable("x", lowBound=0, cat='Integer')
y = pulp.LpVariable("y", lowBound=0, cat='Integer')

# 問題の定義、目的関数、制約条件は前述の通り

4. 複数の目的関数

PuLPでは複数の目的関数を扱うことができますが、ここでは重み付きの合成目的関数として実装します。

1
2
# 目的関数の定義 (例: 50% の重みを持つ2つの目的関数)
model += 0.5 * (40 * x + 30 * y) + 0.5 * (30 * x + 20 * y), "Combined Objective"

5. 制約条件の名前付きリスト

制約条件をリストで管理し、名前を付けて定義することができます。

1
2
3
4
5
6
7
constraints = [
(2 * x + y <= 20, "Raw Material Constraint"),
(4 * x + 3 * y <= 45, "Labor Constraint")
]

for constraint, name in constraints:
model += constraint, name

6. コマンドラインソルバーの利用

PuLPは多くのコマンドラインソルバー(例: CBC, Gurobi, CPLEX)と連携できます。

1
2
3
# CBCソルバーを使用
solver = pulp.PULP_CBC_CMD()
model.solve(solver)

実践的な例:供給と需要の最適化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
import pulp

# 問題の定義
model = pulp.LpProblem("Minimize_Transportation_Cost", pulp.LpMinimize)

# サプライチェーンの定義
warehouses = ["A", "B"]
stores = ["X", "Y", "Z"]
supply = {"A": 300, "B": 400}
demand = {"X": 250, "Y": 350, "Z": 100}

# 輸送コストの定義
costs = {"A": {"X": 2, "Y": 4, "Z": 5}, "B": {"X": 3, "Y": 1, "Z": 3}}

# 変数の定義
vars = pulp.LpVariable.dicts("Route", (warehouses, stores), lowBound=0, cat='Continuous')

# 目的関数の定義
model += pulp.lpSum([vars[w][s] * costs[w][s] for w in warehouses for s in stores]), "Total Cost"

# 制約条件の定義
for w in warehouses:
model += pulp.lpSum([vars[w][s] for s in stores]) <= supply[w], f"Supply_{w}"

for s in stores:
model += pulp.lpSum([vars[w][s] for w in warehouses]) >= demand[s], f"Demand_{s}"

# 問題を解く
model.solve()

# 結果の表示
print(f"Status: {pulp.LpStatus[model.status]}")
for w in warehouses:
for s in stores:
print(f"Route {w} to {s}: {pulp.value(vars[w][s])} units")
print(f"Total Cost: {pulp.value(model.objective)}")

[実行結果]

Status: Optimal
Route A to X: 250.0 units
Route A to Y: 0.0 units
Route A to Z: 50.0 units
Route B to X: 0.0 units
Route B to Y: 350.0 units
Route B to Z: 50.0 units
Total Cost: 1250.0

まとめ

PuLPを使うと、線形計画問題整数計画問題を簡単に定義し解くことができます。

複雑な制約条件や複数の目的関数も扱えるため、実際のビジネス研究最適化問題に広く応用可能です。

Polars

Polars

Polarsは、Pythonのデータフレームライブラリで、高速なデータ処理と分析が可能です。

Pandasと似た使い方をしながらも、特に大規模なデータセットに対して高いパフォーマンスを発揮します。

以下に、Polarsのインストール方法、基本的な使い方、そして現実的な利用例を説明します。

インストール

Polarsをインストールするには、以下のコマンドを使用します。

1
pip install polars

基本的な使い方

データフレームの作成

Polarsのデータフレームは、pl.DataFrameクラスを使用して作成します。

1
2
3
4
5
6
7
8
9
10
import polars as pl

data = {
"name": ["Alice", "Bob", "Charlie"],
"age": [25, 30, 35],
"city": ["New York", "Los Angeles", "Chicago"]
}

df = pl.DataFrame(data)
print(df)

CSVファイルの読み込み

Polarsは大規模なCSVファイルの読み込みにも優れています。

1
2
df = pl.read_csv("your_file.csv")
print(df)

データの操作

Polarsを使った基本的なデータ操作はPandasと似ています。

1
2
3
4
5
6
7
8
9
10
# 列の選択
df = df.select(["name", "age"])

# フィルタリング
df = df.filter(pl.col("age") > 30)

# 列の追加
df = df.with_column((pl.col("age") * 2).alias("age_doubled"))

print(df)

現実的な利用例

データのクレンジング

データのクレンジングは、データ分析の重要なステップです。

Polarsを使うと、大規模なデータセットに対して効率的にクレンジング処理が行えます。

1
2
3
4
5
6
7
# 欠損値の処理
df = df.fill_null("unknown")

# 重複の削除
df = df.unique(subset=["name"])

print(df)

集計処理

Polarsを使った集計処理も非常に効率的です。

1
2
3
4
5
6
7
# グループ化して集計
df = df.groupby("city").agg([
pl.sum("age").alias("total_age"),
pl.mean("age").alias("average_age")
])

print(df)

データの結合

複数のデータフレームを結合する操作もPolarsで簡単に行えます。

1
2
3
4
5
6
7
8
data1 = {"name": ["Alice", "Bob"], "age": [25, 30]}
data2 = {"name": ["Charlie", "David"], "city": ["Chicago", "San Francisco"]}

df1 = pl.DataFrame(data1)
df2 = pl.DataFrame(data2)

df = df1.join(df2, on="name", how="left")
print(df)

パフォーマンスの比較

Polarsは大規模データセットに対して特にパフォーマンスが優れています。

例えば、1億行のデータを生成し、PolarsPandasで処理時間を比較してみます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import time
import pandas as pd

data = {
"a": range(100000000),
"b": range(100000000, 0, -1)
}

# Pandasでの処理
start = time.time()
pdf = pd.DataFrame(data)
pdf["c"] = pdf["a"] + pdf["b"]
print("Pandas処理時間:", time.time() - start)

# Polarsでの処理
start = time.time()
pldf = pl.DataFrame(data)
pldf = pldf.with_column((pl.col("a") + pl.col("b")).alias("c"))
print("Polars処理時間:", time.time() - start)

まとめ

Polarsは、大規模データの処理や分析において非常に効率的であり、Pandasの使い慣れたインターフェースを保ちながらも、より高いパフォーマンスを提供します。

データのクレンジング集計結合など、日常的なデータ操作において非常に便利です。

大規模データセットを扱う際には、Polarsを試してみると良いでしょう。

python-constraint

python-constraint

python-constraintは、制約プログラミングを簡単に行うためのPythonライブラリです。

このライブラリを使うことで、制約条件を満たす組み合わせを効率的に探索することができます。

以下に、python-constraintの基本的な使い方と便利な機能について説明します。

インストール

まず、ライブラリをインストールします。

1
pip install python-constraint

基本的な使い方

以下は、基本的な制約プログラムの例です。

この例では、変数に対していくつかの制約を追加し、制約を満たす解を求めます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from constraint import *

# 問題を定義
problem = Problem()

# 変数の追加
problem.addVariable('x', range(10)) # xは0から9までの値をとる
problem.addVariable('y', range(10)) # yは0から9までの値をとる

# 制約条件の追加
def constraint_func(x, y):
return x + y == 5

problem.addConstraint(constraint_func, ['x', 'y'])

# 解を取得
solutions = problem.getSolutions()

# 解の表示
for solution in solutions:
print(solution)

便利な機能

1. 変数の範囲設定

変数の範囲をリストやジェネレータで指定できます。

1
problem.addVariable('z', [1, 2, 4, 8, 16])

2. 複数変数への同じ制約の追加

同じ制約を複数の変数に対して追加することができます。

1
2
# 全ての変数の値が異なるように制約を追加
problem.addConstraint(AllDifferentConstraint(), ['x', 'y', 'z'])

3. 事前定義された制約

python-constraintにはいくつかの事前定義された制約が用意されています。

  • AllDifferentConstraint(): すべての変数が異なる値を持つように制約する。
  • ExactSumConstraint(total): 変数の合計がtotalに等しくなるように制約する。
1
2
3
4
5
# すべての変数が異なる値を持つ
problem.addConstraint(AllDifferentConstraint(), ['x', 'y', 'z'])

# 変数の合計が10になるように制約する
problem.addConstraint(ExactSumConstraint(10), ['x', 'y', 'z'])

4. カスタム制約の追加

カスタム制約関数を定義して追加することができます。

1
2
3
4
def custom_constraint(a, b):
return a * b <= 20

problem.addConstraint(custom_constraint, ['x', 'y'])

5. 特定の値を除外

特定の値を変数から除外する制約も可能です。

1
problem.addConstraint(lambda x: x != 5, ['x'])  # xは5を取らない

実践例

以下に、より複雑な制約問題の例を示します。

この例では、簡単なナンプレ(数独)パズルの解を求めます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
from constraint import *

# 問題を定義
problem = Problem()

# 変数の追加
rows = range(1, 10)
cols = range(1, 10)

# 各セルに1から9の値を割り当てる
for row in rows:
for col in cols:
problem.addVariable((row, col), range(1, 10))

# 各行に対して全ての数が異なるように制約
for row in rows:
problem.addConstraint(AllDifferentConstraint(), [(row, col) for col in cols])

# 各列に対して全ての数が異なるように制約
for col in cols:
problem.addConstraint(AllDifferentConstraint(), [(row, col) for row in rows])

# 各3x3のボックスに対して全ての数が異なるように制約
for box_row in range(3):
for box_col in range(3):
cells = [(row, col) for row in range(box_row*3+1, box_row*3+4) for col in range(box_col*3+1, box_col*3+4)]
problem.addConstraint(AllDifferentConstraint(), cells)

# 解を取得
solutions = problem.getSolutions()

# 解の表示(最初の解のみ)
if solutions:
solution = solutions[0]
for row in rows:
print([solution[(row, col)] for col in cols])
else:
print("No solution found.")

この例では、ナンプレパズルの各セルに$1$から$9$の値を割り当て、行、列、および$3\times3$のボックス内で全ての数が異なるように制約しています。

python-constraintを使用することで、様々な制約付き問題を効率的に解くことができます。

制約を定義し、それに従って解を探索するこのライブラリの使い方を習得することで、複雑な問題にも対応できるようになります。

Google OR-Tools 4. ジョブショップスケジューリング問題(Job Shop Scheduling)

Google OR-Tools

Google OR-Toolsは、組合せ最適化問題を解決するための強力なツールキットです。

PythonC++, JavaC#などで利用でき、様々な最適化問題(線形最適化、整数最適化、車両ルーティング問題、スケジューリング問題など)に対応しています。

4. ジョブショップスケジューリング問題(Job Shop Scheduling)

ジョブショップスケジューリングは、複数のジョブとマシンがある状況で、ジョブを最適にスケジュールする問題です。

以下のソースコードは、ジョブショップスケジューリング問題をGoogle OR-ToolsCP-SATソルバーを使用して解く例です。

ジョブショップスケジューリング問題は、いくつかのジョブ(作業)があり、それぞれのジョブは特定の順序でマシンを必要とする複数のタスクで構成されています。

この問題は、すべてのジョブが完了するまでの最短時間(メイクスパン)を求めるものです。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
from ortools.sat.python import cp_model

def job_shop_scheduling():
# Create the model.
model = cp_model.CpModel()

jobs_data = [ # (machine_id, processing_time).
[(0, 3), (1, 2), (2, 2)],
[(0, 2), (2, 1), (1, 4)],
[(1, 4), (2, 3)]
]

machines_count = 3
all_tasks = {}
machine_to_intervals = [[] for _ in range(machines_count)]
horizon = sum(task[1] for job in jobs_data for task in job)

# Create job variables and constraints.
for job_id, job in enumerate(jobs_data):
for task_id, task in enumerate(job):
machine = task[0]
duration = task[1]
suffix = '_%i_%i' % (job_id, task_id)
start_var = model.NewIntVar(0, horizon, 'start' + suffix)
end_var = model.NewIntVar(0, horizon, 'end' + suffix)
interval_var = model.NewIntervalVar(start_var, duration, end_var, 'interval' + suffix)
all_tasks[job_id, task_id] = (start_var, end_var, interval_var)
machine_to_intervals[machine].append(interval_var)

# Create and add disjunctive constraints.
for machine_intervals in machine_to_intervals:
model.AddNoOverlap(machine_intervals)

# Precedences inside a job.
for job_id, job in enumerate(jobs_data):
for task_id in range(len(job) - 1):
model.Add(all_tasks[job_id, task_id + 1][0] >= all_tasks[job_id, task_id][1])

# Makespan objective.
obj_var = model.NewIntVar(0, horizon, 'makespan')
model.AddMaxEquality(obj_var, [all_tasks[job_id, len(job) - 1][1] for job_id, job in enumerate(jobs_data)])
model.Minimize(obj_var)

# Solve the model.
solver = cp_model.CpSolver()
status = solver.Solve(model)

# Print solution.
if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
print('Optimal Schedule Length: %i' % solver.ObjectiveValue())
print()
for job_id, job in enumerate(jobs_data):
print('Job %i:' % job_id)
for task_id, task in enumerate(job):
start_var, end_var, _ = all_tasks[job_id, task_id]
print('Task %i starts at %i and ends at %i' % (
task_id, solver.Value(start_var), solver.Value(end_var)))
else:
print('No solution found.')

job_shop_scheduling()

ソースコード解説

ソースコードの各部分を説明します。

1. インポートと関数定義

1
2
3
4
5
from ortools.sat.python import cp_model

def job_shop_scheduling():
# Create the model.
model = cp_model.CpModel()

まず、Google OR-ToolsCP-SATソルバーからcp_modelモジュールをインポートし、スケジューリングを行う関数job_shop_schedulingを定義します。

その後、制約プログラミングモデルを作成します。

2. ジョブデータの定義

1
2
3
4
5
6
7
8
9
10
jobs_data = [  # (machine_id, processing_time).
[(0, 3), (1, 2), (2, 2)],
[(0, 2), (2, 1), (1, 4)],
[(1, 4), (2, 3)]
]

machines_count = 3
all_tasks = {}
machine_to_intervals = [[] for _ in range(machines_count)]
horizon = sum(task[1] for job in jobs_data for task in job)

ジョブデータを定義します。

ここでは、各ジョブはタスクのリストで構成されており、各タスクは使用するマシンIDと処理時間を持っています。

また、マシンの数、すべてのタスクを保持する辞書、各マシンに対するタスクのリスト、およびスケジュールの全期間(ホライズン)を計算します。

3. 変数と制約の作成

1
2
3
4
5
6
7
8
9
10
11
# Create job variables and constraints.
for job_id, job in enumerate(jobs_data):
for task_id, task in enumerate(job):
machine = task[0]
duration = task[1]
suffix = '_%i_%i' % (job_id, task_id)
start_var = model.NewIntVar(0, horizon, 'start' + suffix)
end_var = model.NewIntVar(0, horizon, 'end' + suffix)
interval_var = model.NewIntervalVar(start_var, duration, end_var, 'interval' + suffix)
all_tasks[job_id, task_id] = (start_var, end_var, interval_var)
machine_to_intervals[machine].append(interval_var)

各タスクに対して、開始時間、終了時間、およびそれらを結びつけるインターバル変数を作成します。

これらの変数をall_tasks辞書に格納し、対応するマシンのインターバルリストに追加します。

4. 排他制約の追加

1
2
3
# Create and add disjunctive constraints.
for machine_intervals in machine_to_intervals:
model.AddNoOverlap(machine_intervals)

同じマシンで同時に複数のタスクが実行されないようにするため、AddNoOverlap制約を追加します。

この制約は、同じマシンに割り当てられたすべてのインターバル変数に対して適用されます。

5. ジョブ内の順序制約の追加

1
2
3
4
# Precedences inside a job.
for job_id, job in enumerate(jobs_data):
for task_id in range(len(job) - 1):
model.Add(all_tasks[job_id, task_id + 1][0] >= all_tasks[job_id, task_id][1])

各ジョブ内のタスクが定められた順序で実行されるようにするための制約を追加します。

具体的には、あるタスクの終了時間が次のタスクの開始時間よりも早くなるようにします。

6. メイクスパンの目標設定

1
2
3
4
# Makespan objective.
obj_var = model.NewIntVar(0, horizon, 'makespan')
model.AddMaxEquality(obj_var, [all_tasks[job_id, len(job) - 1][1] for job_id, job in enumerate(jobs_data)])
model.Minimize(obj_var)

メイクスパン(全ジョブが完了するまでの最短時間)を最小化することを目標として設定します。

AddMaxEquality制約を使用して、メイクスパン変数がすべてのジョブの最終タスクの終了時間の最大値に等しくなるようにします。

7. モデルの解決

1
2
3
# Solve the model.
solver = cp_model.CpSolver()
status = solver.Solve(model)

CP-SATソルバーを作成し、モデルを解決します。

8. 解の表示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
    # Print solution.
if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
print('Optimal Schedule Length: %i' % solver.ObjectiveValue())
print()
for job_id, job in enumerate(jobs_data):
print('Job %i:' % job_id)
for task_id, task in enumerate(job):
start_var, end_var, _ = all_tasks[job_id, task_id]
print('Task %i starts at %i and ends at %i' % (
task_id, solver.Value(start_var), solver.Value(end_var)))
else:
print('No solution found.')

job_shop_scheduling()

解が見つかった場合、そのメイクスパンと各タスクの開始時間と終了時間を表示します。

解が見つからなかった場合は、その旨を表示します。


このコード全体で、ジョブショップスケジューリング問題を解決するための制約プログラミングモデルを定義し、最適なスケジュールを求めています。

結果解説

[実行結果]

Optimal Schedule Length: 11

Job 0:
Task 0 starts at 2 and ends at 5
Task 1 starts at 5 and ends at 7
Task 2 starts at 7 and ends at 9
Job 1:
Task 0 starts at 0 and ends at 2
Task 1 starts at 2 and ends at 3
Task 2 starts at 7 and ends at 11
Job 2:
Task 0 starts at 0 and ends at 4
Task 1 starts at 4 and ends at 7

この結果は、ジョブショップスケジューリング問題に対して最適なスケジュールが見つかり、そのメイクスパン(全ジョブが完了するまでの最短時間)が$11$であることを示しています。

以下、解説します。

1. メイクスパンの解説

Optimal Schedule Length: 11

  • メイクスパンは$11$です。
    これは、すべてのジョブが完了するまでに最短で$11$単位時間かかることを意味します。

2. ジョブ0のスケジュール

Job 0:

  • Task 0 starts at 2 and ends at 5: このタスクはマシン$0$で処理され、処理時間は$3$単位時間です。
  • Task 1 starts at 5 and ends at 7: このタスクはマシン$1$で処理され、処理時間は$2$単位時間です。
  • Task 2 starts at 7 and ends at 9: このタスクはマシン$2$で処理され、処理時間は$2$単位時間です。

3. ジョブ1のスケジュール

Job 1:

  • Task 0 starts at 0 and ends at 2: このタスクはマシン$0$で処理され、処理時間は$2$単位時間です。
  • Task 1 starts at 2 and ends at 3: このタスクはマシン$2$で処理され、処理時間は$1$単位時間です。
  • Task 2 starts at 7 and ends at 11: このタスクはマシン$1$で処理され、処理時間は$4$単位時間です。

4. ジョブ2のスケジュール

Job 2:

  • Task 0 starts at 0 and ends at 4: このタスクはマシン$1$で処理され、処理時間は$4$単位時間です。
  • Task 1 starts at 4 and ends at 7: このタスクはマシン$2$で処理され、処理時間は$3$単位時間です。

5. スケジュール全体の概要

このスケジュールにより、以下のようなタスクの配置が見られます:

  • マシン0:

    • Job 1のTask 0: $0$から$2$
    • Job 0のTask 0: $2$から$5$
  • マシン1:

    • Job 2のTask 0: $0$から$4$
    • Job 0のTask 1: $5$から$7$
    • Job 1のTask 2: $7$から$11$
  • マシン2:

    • Job 1のTask 1: $2$から$3$
    • Job 2のTask 1: $4$から$7$
    • Job 0のTask 2: $7$から$9$

このように、各マシンでタスクが重ならないようにスケジュールされています。

また、ジョブ内のタスクが指定された順序で実行されています。

結果として、このスケジュールによりすべてのジョブが最短時間($11$単位時間)で完了することが確認できます。

Google OR-Tools 3. 車両ルーティング問題(Vehicle Routing Problem, VRP)

Google OR-Tools

Google OR-Toolsは、組合せ最適化問題を解決するための強力なツールキットです。

PythonC++, JavaC#などで利用でき、様々な最適化問題(線形最適化、整数最適化、車両ルーティング問題、スケジューリング問題など)に対応しています。

3. 車両ルーティング問題(Vehicle Routing Problem, VRP)

VRPは、指定された地点を巡回する車両の最適なルートを見つける問題です。

サンプルコード:VRP

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
from ortools.constraint_solver import pywrapcp
from ortools.constraint_solver import routing_enums_pb2

def create_data_model():
data = {}
data['distance_matrix'] = [
[0, 2, 9, 10],
[1, 0, 6, 4],
[15, 7, 0, 8],
[6, 3, 12, 0],
]
data['num_vehicles'] = 1
data['depot'] = 0
return data

def print_solution(manager, routing, solution):
print('Objective: {}'.format(solution.ObjectiveValue()))
index = routing.Start(0)
plan_output = 'Route for vehicle 0:\n'
route_distance = 0
while not routing.IsEnd(index):
plan_output += ' {} ->'.format(manager.IndexToNode(index))
previous_index = index
index = solution.Value(routing.NextVar(index))
route_distance += routing.GetArcCostForVehicle(previous_index, index, 0)
plan_output += ' {}\n'.format(manager.IndexToNode(index))
plan_output += 'Distance of the route: {}m\n'.format(route_distance)
print(plan_output)

def vehicle_routing_problem():
data = create_data_model()
manager = pywrapcp.RoutingIndexManager(len(data['distance_matrix']), data['num_vehicles'], data['depot'])
routing = pywrapcp.RoutingModel(manager)
def distance_callback(from_index, to_index):
from_node = manager.IndexToNode(from_index)
to_node = manager.IndexToNode(to_index)
return data['distance_matrix'][from_node][to_node]
transit_callback_index = routing.RegisterTransitCallback(distance_callback)
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)
solution = routing.SolveWithParameters(search_parameters)
if solution:
print_solution(manager, routing, solution)

vehicle_routing_problem()

ソースコード解説

ソースコードの各部分を説明します。

1. インポート

1
2
from ortools.constraint_solver import pywrapcp
from ortools.constraint_solver import routing_enums_pb2
  • pywrapcpOR-Tools制約ソルバーを提供します。
  • routing_enums_pb2 はルーティング用の定数を提供します。

2. データモデルの作成

1
2
3
4
5
6
7
8
9
10
11
def create_data_model():
data = {}
data['distance_matrix'] = [
[0, 2, 9, 10],
[1, 0, 6, 4],
[15, 7, 0, 8],
[6, 3, 12, 0],
]
data['num_vehicles'] = 1
data['depot'] = 0
return data
  • create_data_model 関数は、距離行列や車両の数、出発点(デポ)などのデータを定義します。
  • distance_matrix は各地点間の距離を示します。
  • num_vehicles は使用する車両の数を示します。
  • depot は出発点のインデックスを示します。

3. 解の表示

1
2
3
4
5
6
7
8
9
10
11
12
13
def print_solution(manager, routing, solution):
print('Objective: {}'.format(solution.ObjectiveValue()))
index = routing.Start(0)
plan_output = 'Route for vehicle 0:\n'
route_distance = 0
while not routing.IsEnd(index):
plan_output += ' {} ->'.format(manager.IndexToNode(index))
previous_index = index
index = solution.Value(routing.NextVar(index))
route_distance += routing.GetArcCostForVehicle(previous_index, index, 0)
plan_output += ' {}\n'.format(manager.IndexToNode(index))
plan_output += 'Distance of the route: {}m\n'.format(route_distance)
print(plan_output)
  • print_solution 関数は、解(ルートと距離)を出力します。
  • solution.ObjectiveValue()目的関数の値を取得します。
  • routing.Start(0) で車両$ 0 $のルートを開始します。
  • ループ内で、各地点のインデックスを取得し、次の地点へ進みます。
  • ルート全体の距離を計算し、最終的にルート距離を出力します。

4. 車両経路問題の定義と解の計算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def vehicle_routing_problem():
data = create_data_model()
manager = pywrapcp.RoutingIndexManager(len(data['distance_matrix']), data['num_vehicles'], data['depot'])
routing = pywrapcp.RoutingModel(manager)

def distance_callback(from_index, to_index):
from_node = manager.IndexToNode(from_index)
to_node = manager.IndexToNode(to_index)
return data['distance_matrix'][from_node][to_node]

transit_callback_index = routing.RegisterTransitCallback(distance_callback)
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)

solution = routing.SolveWithParameters(search_parameters)
if solution:
print_solution(manager, routing, solution)

vehicle_routing_problem()
  • vehicle_routing_problem 関数は、車両経路問題を定義し、解を計算します。
  • create_data_model を呼び出してデータを取得します。
  • RoutingIndexManager を使用してルーティングインデックスマネージャを作成します。
  • RoutingModel を使用してルーティングモデルを作成します。
  • distance_callback 関数を定義し、各地点間の距離を返します。
  • RegisterTransitCallback を使用してコールバックを登録し、距離行列を評価します。
  • SetArcCostEvaluatorOfAllVehicles を使用してすべての車両のアークコスト評価器を設定します。
  • DefaultRoutingSearchParameters を使用して検索パラメータを設定し、最初の解法戦略を指定します。
  • SolveWithParameters を使用して問題を解き、解が見つかった場合は print_solution を呼び出して結果を表示します。

5. 全体の流れ

  1. データモデルを作成し、距離行列車両情報を定義します。
  2. ルーティングインデックスマネージャとルーティングモデルを作成します。
  3. 距離評価用のコールバック関数を定義し、登録します。
  4. 検索パラメータを設定し、最適なルートを求めます。
  5. 解が見つかった場合、ルート距離を出力します。

以上がこのコードの詳細な説明です。

結果解説

[実行結果]

Objective: 21
Route for vehicle 0:
 0 -> 2 -> 3 -> 1 -> 0
Distance of the route: 21m

詳細な解説

  1. Objective: 21

    • 目的関数の値を示しています。
      この場合、車両が訪問するルートの総距離が$21$メートルであることを意味します。
      ここでは最小化された距離が$21$メートルという結果を得ています。
  2. Route for vehicle 0: 0 -> 2 -> 3 -> 1 -> 0

    • ルートの詳細を示しています。
      車両$0$が訪問する順番を示しており、出発点から始まり、各ノードを訪問し、最後に出発点に戻ります。
      • 0 は出発点(デポ)。
      • 2 は2番目の地点。
      • 3 は3番目の地点。
      • 1 は1番目の地点。
      • 0 は再び出発点に戻ります。
  3. Distance of the route: 21m

    • ルート全体の距離が$21$メートルであることを示しています。
      ルートの距離は、各ステップ間の距離の合計です。

Google OR-Tools 2. 整数最適化(Integer Optimization)

Google OR-Tools

Google OR-Toolsは、組合せ最適化問題を解決するための強力なツールキットです。

PythonC++, JavaC#などで利用でき、様々な最適化問題(線形最適化、整数最適化、車両ルーティング問題、スケジューリング問題など)に対応しています。

2. 整数最適化(Integer Optimization)

整数最適化では、変数が整数値を取るような最適化問題を解きます。

サンプルコード:整数最適化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
from ortools.linear_solver import pywraplp

def integer_optimization_example():
# Solverの作成
solver = pywraplp.Solver.CreateSolver('SCIP')
if not solver:
return

# 変数の作成
x = solver.IntVar(0, 10, 'x')
y = solver.IntVar(0, 10, 'y')

# 制約条件の追加
solver.Add(2 * x + 3 * y <= 12)

# 目的関数の設定
solver.Maximize(3 * x + 4 * y)

# 解を求める
status = solver.Solve()

# 結果を表示
if status == pywraplp.Solver.OPTIMAL:
print('Objective value =', solver.Objective().Value())
print('x =', x.solution_value())
print('y =', y.solution_value())
else:
print('The problem does not have an optimal solution.')

integer_optimization_example()

ソースコード解説

1. pywraplp モジュールのインポート

最初に、ortools.linear_solver モジュールから pywraplp をインポートします。
このモジュールは、Google OR-Toolsの線形および整数ソルバーを使用するためのインターフェースを提供します。

1
from ortools.linear_solver import pywraplp

2. ソルバーの作成

次に、整数最適化問題を解くためのソルバーを作成します。
CreateSolver('SCIP') メソッドは、SCIP(Solving Constraint Integer Programs)ソルバーを作成します。

1
2
3
solver = pywraplp.Solver.CreateSolver('SCIP')
if not solver:
return
  • solver が作成されない場合、None が返され、関数を終了します。

3. 変数の作成

ソルバーで使用する整数変数を作成します。
IntVar メソッドを使って変数を定義します。
ここでは、変数 xy がそれぞれ0から10の範囲に制限されています。

1
2
x = solver.IntVar(0, 10, 'x')
y = solver.IntVar(0, 10, 'y')
  • xy は整数変数で、それぞれ$0$から$10$の範囲の整数値を取ります。

4. 制約条件の追加

問題の制約条件を追加します。

ここでは、2 * x + 3 * y <= 12 という制約を追加しています。

1
solver.Add(2 * x + 3 * y <= 12)
  • この制約は、2倍の x と3倍の y の合計が$12$以下であることを要求します。

5. 目的関数の設定

最適化する目的関数を設定します。ここでは、3 * x + 4 * y を最大化するように設定しています。

1
solver.Maximize(3 * x + 4 * y)
  • 目的関数は 3 * x + 4 * y で、これを最大化します。

6. 解を求める

ソルバーを呼び出して、最適解を求めます。

Solve メソッドは、ソルバーを実行して最適化問題を解きます。

1
status = solver.Solve()
  • Solve() メソッドは、ソルバーを実行し、解が見つかったかどうかを示すステータスコードを返します。

7. 結果の表示

最適解が見つかったかどうかを status で確認し、結果を表示します。

pywraplp.Solver.OPTIMAL は、問題に最適解が見つかったことを示します。

1
2
3
4
5
6
if status == pywraplp.Solver.OPTIMAL:
print('Objective value =', solver.Objective().Value())
print('x =', x.solution_value())
print('y =', y.solution_value())
else:
print('The problem does not have an optimal solution.')
  • statuspywraplp.Solver.OPTIMAL であれば、最適解が見つかっています。
  • solver.Objective().Value()目的関数の値を取得します。
  • x.solution_value() および y.solution_value() で各変数の最適解を取得します。

8. 関数の実行

最後に、定義した関数を実行します。

1
integer_optimization_example()

このソースコード全体は、SCIPソルバーを使用して、整数最適化問題を解き、その結果を表示するプロセスを示しています。

結果解説

[実行結果]

Objective value = 18.0
x = 6.0
y = 0.0

この結果は、与えられた整数最適化問題の最適解を示しています。

結果の詳細

1
2
3
Objective value = 18.0
x = 6.0
y = 0.0

問題の設定の再確認

まず、ソースコードで設定した問題の内容を再確認します。

  • 変数の範囲:

    • x は$0$から$10$の範囲の整数
    • y は$0$から$10$の範囲の整数
  • 制約条件:

    • 2 * x + 3 * y <= 12
  • 目的関数:

    • 3 * x + 4 * y を最大化

結果の意味

  1. Objective value = 18.0

    • 目的関数 3 * x + 4 * y の最大値が$18$であることを示しています。
  2. x = 6.0

    • 変数 x の最適値が$6$であることを示しています。
  3. y = 0.0

    • 変数 y の最適値が$0$であることを示しています。

制約条件の確認

最適解 (x = 6, y = 0)制約条件を満たしていることを確認します。

  • 制約条件: 2 * x + 3 * y <= 12
    • 2 * 6 + 3 * 0 = 12 なので、この制約を満たしています。

目的関数の値の計算

最適解 (x = 6, y = 0) に対する目的関数の値を計算します。

  • 目的関数: 3 * x + 4 * y
    • 3 * 6 + 4 * 0 = 18

これにより、目的関数の値が$18$であることが確認できます。

結果の解釈

この結果から、以下のことがわかります:

  • 制約条件 2 * x + 3 * y <= 12 を満たしながら、目的関数 3 * x + 4 * y最大化するためには、x = 6y = 0 に設定するのが最適であり、そのときの目的関数の値は$18$である。
  • つまり、x を6に、y を$0$にすることで、与えられた制約条件の下で最大の利益(目的関数の値)を得ることができます。

このようにして、整数最適化問題の最適解を求め、その結果を評価することができます。

Google OR-Tools 1. 線形最適化(Linear Optimization)

Google OR-Tools

Google OR-Toolsは、組合せ最適化問題を解決するための強力なツールキットです。

PythonC++, JavaC#などで利用でき、様々な最適化問題(線形最適化、整数最適化、車両ルーティング問題、スケジューリング問題など)に対応しています。

1. 線形最適化(Linear Optimization)

線形最適化では、線形の目的関数最大化または最小化する問題を解きます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
!pip install ortools
from ortools.linear_solver import pywraplp

def linear_optimization_example():
# Solverの作成
solver = pywraplp.Solver.CreateSolver('GLOP')
if not solver:
return

# 変数の作成
x = solver.NumVar(0, 1, 'x')
y = solver.NumVar(0, 2, 'y')

# 制約条件の追加
solver.Add(x + y <= 2)

# 目的関数の設定
solver.Maximize(x + y)

# 解を求める
status = solver.Solve()

# 結果を表示
if status == pywraplp.Solver.OPTIMAL:
print('Objective value =', solver.Objective().Value())
print('x =', x.solution_value())
print('y =', y.solution_value())
else:
print('The problem does not have an optimal solution.')

linear_optimization_example()

コード解説

ソースコードを説明します。

1. インポート

1
from ortools.linear_solver import pywraplp

Google OR-Toolsの線形最適化ライブラリであるpywraplpをインポートします。

これにより、線形最適化問題を解くための機能を利用できます。

2. 関数定義

1
def linear_optimization_example():

線形最適化問題の例を示す関数linear_optimization_exampleを定義します。

3. ソルバーの作成

1
2
3
solver = pywraplp.Solver.CreateSolver('GLOP')
if not solver:
return

ここでは、GLOP(Google Linear Optimization Package)という線形最適化ソルバーを作成します。

CreateSolver('GLOP')は新しいソルバーオブジェクトを作成します。

ソルバーが作成できなかった場合、関数を終了します。

4. 変数の作成

1
2
x = solver.NumVar(0, 1, 'x')
y = solver.NumVar(0, 2, 'y')

線形最適化問題の変数を作成します。ここでは、変数xyを定義します。

  • xは$0$から$1$までの値を取ります。
  • yは$0$から$2$までの値を取ります。

5. 制約条件の追加

1
solver.Add(x + y <= 2)

変数に対する制約を追加します。ここでは、x + y <= 2という制約条件を設定しています。

これは、変数xyの和が2以下でなければならないことを意味します。

6. 目的関数の設定

1
solver.Maximize(x + y)

最適化の目的を設定します。ここでは、x + yを最大化するように設定しています。

7. 解を求める

1
status = solver.Solve()

設定した制約条件目的関数に基づいて、最適解を求めます。

このメソッドは、解が見つかったかどうかを示すステータスコードを返します。

8. 結果の表示

1
2
3
4
5
6
if status == pywraplp.Solver.OPTIMAL:
print('Objective value =', solver.Objective().Value())
print('x =', x.solution_value())
print('y =', y.solution_value())
else:
print('The problem does not have an optimal solution.')

最適解が見つかった場合、目的関数の値と変数xおよびyの値を表示します。

解が見つからなかった場合は、最適解が存在しないことを表示します。

9. 関数の呼び出し

1
linear_optimization_example()

定義した関数を呼び出し、線形最適化問題を実行します。

まとめ

このコードは、OR-Toolsを使用して簡単な線形最適化問題を解決する方法を示しています。

具体的には、2つの変数xyを定義し、制約条件x + y <= 2のもとでx + yを最大化する問題を解いています。

最適解が見つかった場合、その解を表示します。

結果解説

[実行結果]

Objective value = 2.0
x = 0.0
y = 2.0

この結果は、線形最適化問題の解として得られた最適解を示しています。
各部分の意味を以下に説明します。

結果の詳細

Objective value = 2.0

目的関数の値が$2.0$であることを示しています。
つまり、x + yを最大化する目的関数に対して得られた最大値が2.0です。

x = 0.0

変数xの最適解が$0.0$であることを示しています。

y = 2.0

変数yの最適解が$2.0$であることを示しています。

解の詳細

この結果に基づいて、以下の点を解説します。

  1. 目的関数の最大化:

    • 目的関数はx + yを最大化することを目指しています。
  2. 制約条件の確認:

    • 制約条件はx + y <= 2です。
      この条件により、xyの和が$2$以下でなければなりません。
  3. 最適解の妥当性:

    • 得られた解はx = 0.0y = 2.0です。
    • この解は制約条件x + y <= 2を満たしており、0.0 + 2.0 = 2.0となっています。
      制約条件の上限である$2$に一致しています。
  4. 最大化の達成:

    • 目的関数x + yを最大化した結果が$2.0$であることは、制約条件の下で可能な最大値であることを示しています。
    • 他の値の組み合わせ(例えば、x = 1y = 1など)では、x + yの値が2を超えることはありません。

結論

この結果は、与えられた線形最適化問題の最適解として、変数xが$0.0$、変数yが$2.0$であり、そのときの目的関数x + yの値が最大で$2.0$であることを示しています。

制約条件x + y <= 2を満たしながら、目的関数の値を最大化するための最適な値の組み合わせがx = 0.0y = 2.0であることが確認されました。