次の線形最適化問題を Python で解いてみよう．
このような線形最適化問題はPuLP や SciPy いうモジュールを使えば，Python で簡単かつ高速に解くことができます．
PuLP は Anaconda をインストールしただけではインストールされていないので，まずは PuLP をインストールします．PuLP のインストールは「Anaconda Prompt」または「Anaconda Powershell Prompt」で pip install pulp を実行すれば良いでしょう．なお，学内の情報処理実習室では pip install pulp --user のように --user オプションを指定する必要があるかもしれませんし，Jupyter Lab の起動後に Anaconda Prompt で pip install pulp を実行した場合には Jupyter Lab を再起動する必要があるかもしれません．また，インストールされているモジュールの一覧は pip list で確認できます．Jupyter Lab や Jupyter Notebook のマジックコマンドを利用しても構いません．
まずは PuLP モジュールをインポートします．もしもエラーが表示されたら，準備を参考に，PuLP をインストールしよう．
モジュールのインポート
import pulp
次に問題オブジェクトを生成しよう．問題オブジェクトは pulp.LpProblem() 関数で生成できます．このとき sense 引数には，最大化問題を解きたい場合には sense=pulp.LpMaximize を，最小化問題を解きたい場合には sense=pulp.LpMinimize を指定します．なお，sense 引数を省略した場合は，最小化問題として生成されることに注意してください．
問題オブジェクトの生成
problem = pulp.LpProblem(name='例題1', sense=pulp.LpMaximize)
続いて，変数オブジェクトを生成します．x1 と x2 からなる2変数の問題を解きたいので，pulp.LpVariable 関数を2回呼び出します．このとき，x1 と x2 の最小値制約を lowBound = 0 のように指定することができます．
変数オブジェクトの生成
x1 = pulp.LpVariable('x1', lowBound=0)
x2 = pulp.LpVariable('x2', lowBound=0)
いよいよ目的関数を設定します．目的関数は次のとおり既に生成した問題オブジェクト problem に追加代入すれば設定できます．
目的関数の設定
problem += 3 * x1 + 4 * x2
さらに制約式を追加しよう．次のとおり目的関数と同様に問題オブジェクト problem に追加するだけで設定できてしまいます．
制約式の設定
problem += x1 + x2 <= 5
problem += x1 + 2 * x2 <= 7
problem += x1 <= 4
problem += x2 <= 3
なお，制約式には次のような演算子が利用できます．「等しい」は == のように = を2つ並べて書くことに注意してください．
| 演算子 | 意味 |
|---|---|
== | 等しい |
<= | 以下 |
>= | 以上 |
ここで設定した問題を表示し，その内容を確認しよう．
問題の確認
print(problem)
例題1: MAXIMIZE 3*x1 + 4*x2 + 0 SUBJECT TO _C1: x1 + x2 <= 5 _C2: x1 + 2 x2 <= 7 _C3: x1 <= 4 _C4: x2 <= 3 VARIABLES x1 Continuous x2 Continuous
問題を解く前に，現在の問題のステータスを確認しておきます．もちろん，状態はまだ解かれていない Not Solved です．
問題のステータスを確認
print("Status : ", pulp.LpStatus[problem.status])
Status : Not Solved
変数の一覧とその値も確認しておきます．変数 x1 と x2 が定義されているが，値はまだ設定されていないことがわかります．
変数の確認
for v in problem.variables():
print(v.name, "=", v.varValue)
x1 = None x2 = None
solve() 関数で問題を解きます．このとき，計算ログが多く表示されます．1行目の代わりに2行目のコードを利用すると計算ログを表示しないようにできます．
問題を解く
problem.solve()
# problem.solve(pulp.PULP_CBC_CMD(msg=0)) # ログを出力しない
1
問題が解けたようなので，再度ステータスを確認します．最適解 (Optimal solution) が得られたことがわかりました．
問題のステータス
print("Status : ", pulp.LpStatus[problem.status])
Status : Optimal
変数の値を確認します．これらの値が目的関数を最大化する変数の値です．
変数の値
for v in problem.variables():
print(v.name, "=", v.varValue)
x1 = 3.0 x2 = 2.0
最大化された目的関数の値を確認します．
目的関数の値
print("z = ", pulp.value(problem.objective))
z = 17.0
すべてのコードをまとめると次のようになります．
プログラム全体
# モジュールをインポート
import pulp
# 問題オブジェクトを作る
problem = pulp.LpProblem(name='例題1', sense=pulp.LpMaximize)
# 変数オブジェクトを作る
x1 = pulp.LpVariable('x1', lowBound=0)
x2 = pulp.LpVariable('x2', lowBound=0)
# 目的関数を設定する
problem += 3 * x1 + 4 * x2
# 制約条件を作る
problem += x1 + x2 <= 5
problem += x1 + 2 * x2 <= 7
problem += x1 <= 4
problem += x2 <= 3
# 問題を表示する
print(problem)
# 問題を解く
problem.solve()
# problem.solve(pulp.PULP_CBC_CMD(msg=0)) # ログを出力しない
# 問題のステータス
print("Status : ", pulp.LpStatus[problem.status])
# 結果を表示しよう
for v in problem.variables():
print(v.name, "=", v.varValue)
# 最適解を計算する
print("z = ", pulp.value(problem.objective))
例題1: MAXIMIZE 3*x1 + 4*x2 + 0 SUBJECT TO _C1: x1 + x2 <= 5 _C2: x1 + 2 x2 <= 7 _C3: x1 <= 4 _C4: x2 <= 3 VARIABLES x1 Continuous x2 Continuous Status : Optimal x1 = 3.0 x2 = 2.0 z = 17.0
なお，pulp.LpVariable で作成した変数は通常の Python 変数とは異なります．したがって，次のように print(x1) としても値を表示できません．
変数の表示
print(x1)
x1
pulp.LpVariable で作成した変数の型は pulp.pulp.LpVariable 型です．
変数の型
type(x1)
pulp.pulp.LpVariable
これら pulp.pulp.LpVariable 型変数の値を表示するには，次のような方法があります．
変数の値の取得方法
print(x1.varValue)
print(x1.value())
print(pulp.value(x1))
3.0 3.0 3.0
よって，3つ目の書き方を使う場合には，次のような方法で最適解を自分で計算することも可能です．
最適解の計算
# 最適解を計算する
z = 3 * pulp.value(x1) + 4 * pulp.value(x2)
print("z =", z)
z = 17.0
ここで説明しているように，非推奨ではありますが，ライブラリを from pulp import * という書き方でインポートする方法もあります．この場合は次のように pulp ライブラリで定義される関数や定数の全てで pulp. を付ける必要がなくなります．コードの入力は少なくて済みますが，他のモジュールで同じ名前の関数などが定義されていた場合には思わぬ結果になる可能性があります．
プログラムの全体（非推奨の書き方）
# モジュールをインポート
from pulp import *
# 問題オブジェクトを作る
problem = LpProblem(name='例題1', sense=LpMaximize)
# 変数オブジェクトを作る
x1 = LpVariable('x1', lowBound=0)
x2 = LpVariable('x2', lowBound=0)
# 目的関数を設定する
problem += 3 * x1 + 4 * x2
# 制約条件を作る
problem += x1 + x2 <= 5
problem += x1 + 2 * x2 <= 7
problem += x1 <= 4
problem += x2 <= 3
# 問題を表示する
print(problem)
# 問題を解く
problem.solve()
# problem.solve(PULP_CBC_CMD(msg=0)) # ログを出力しない
# 問題のステータス
print("Status : ", LpStatus[problem.status])
# 結果を表示しよう
for v in problem.variables():
print(v.name, "=", v.varValue)
# 最適解を計算する
print("z = ", value(problem.objective))
SciPyを使っても線形最適化問題を解くことができます．解きたい問題は次のような最大化問題でした．
現時点（2024年4月時点）では SciPy は最小化問題しか取り扱うことができないので，目的関数を \(-1\) 倍することで最大化問題を次のような最小化問題に書き換えます．
この問題を行列とベクトルを使った形式
に書き直します．
問題の定式化ができたので，モジュールをインポートします．
モジュールのインポート
from scipy.optimize import linprog
ベクトルや行列を順番に定義します．まず，目的関数のベクトル c を定義します．
c = [-3, -4]
制約式の左辺にある行列 A を定義します．
A = [[1, 1],[1,2]]
制約式の右辺にあるベクトル b を定義します．
b = [5, 7]
変数 x1 と x2 の上限・下限を定義します．なお，上限，または下限が必要ない場合は None を指定することに注意してください．
x1_bounds = (0, 4)
x2_bounds = (0, 3)
線形最適化問題を解き，その結果を result に代入します．
result = linprog(c, A_ub=A, b_ub=b, bounds=[x1_bounds, x2_bounds])
結果の詳細は次の通り取得できます．
result
message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
success: True
status: 0
fun: -17.0
x: [ 3.000e+00 2.000e+00]
nit: 2
lower: residual: [ 3.000e+00 2.000e+00]
marginals: [ 0.000e+00 0.000e+00]
upper: residual: [ 1.000e+00 1.000e+00]
marginals: [ 0.000e+00 0.000e+00]
eqlin: residual: []
marginals: []
ineqlin: residual: [ 0.000e+00 0.000e+00]
marginals: [-2.000e+00 -1.000e+00]
mip_node_count: 0
mip_dual_bound: 0.0
mip_gap: 0.0
上の結果から最適解の値だけを表示します．
result.fun
-17.0
最大化問題の目的関数を \(-1\) 倍して最小化問題に書き換えていたことを思い出して，得られた値を \(-1\) 倍すると最大化問題の最適解になります．
result.fun * (-1)
17.0
変数の値を確認します．
result.x
array([3., 2.])