第11章 数理最適化モジュール PuLPとOpenOpt

丼チェーン店長の悩み トンコケ丼，コケトン丼，ミックス丼の3種類 トンコケ丼を1杯作るには，200グラムの豚肉と100グラムの鶏肉コケトン丼を1杯作るには，100グラムの豚肉と200グラムの鶏肉ミックス丼は，豚肉，鶏肉，牛肉を100グラムずつミックス 豚，鶏，牛の肉は，最大1日あたり6キログラム，6キログラム，3キログラム 販売価格は，トンコケ丼1杯1500円，コケトン丼1杯1800円，ミックス丼1杯3000円 お店の利益を最大にするためには，あなたは丼を何杯ずつ作るように指示を出せばよいのだろうか？

丼チェーン店長の悩み 種類 トンコケ丼 コケトン丼 ミックス丼 使用可能量 （百グラム） 豚 2 1 60 鶏 牛 30 利益 （百円） 30 利益 （百円） 15 18

定式化（線形最適化問題） トンコケ丼の数：x1（杯） 変数 コケトン丼の数：x2 （杯） ミックス丼の数：x3 （杯） 目的関数 材料の使用可能量の上限 制約式 （百グラム） 丼の数は非負

Pythonによるモデル化 (1) モジュールの読み込み（mypulpは商用ソルバーGurobiと同じ文法） from mypulp import * モデルオブジェクトの生成 model = Model(“適当な名前") GuとRothbergとBixbyによる最速ソルバー http://www.gurobi.com/ 大学内でインストール可（無料）方法はノートPC持参でＴＡに聞いてね

Pythonによるモデル化 (2) 変数オブジェクトの追加 x1 = model.addVar(name='x1') x2 = model.addVar(name='x2') x3 = model.addVar(name='x3') モデルの更新 (Gurobiだと制約を追加する前には必須． PuLPだと省略可） model.update() 引数 name で名前をつけておく （省略可）

Pythonによるモデル化 (3) 目的関数の設定 model.setObjective(15*x1 + 18*x2 + 30*x3, GRB.MAXIMIZE) 制約の追加 model.addConstr(2*x1 + x2 + x3 <= 60) model.addConstr(x1 + 2*x2 + x3 <= 60) model.addConstr(x3 <= 30) 最適化 model.optimize() setObjectiveメソッドで 目的関数を設定 最大化を指定 GRB は Gurobi の略 addConstrメソッドで 制約を追加 optimizeメソッドで求解

Pythonによるモデル化 (4) 目的関数値の出力 print('Opt. Value=', model.ObjVal) >>> Opt. Value= 1230.0 最適解の出力 print(x1.VarName, x1.X) >>> x1 10.0 print(x2.VarName, x2.X) >>> x2 10.0 print(x3.VarName, x3.X) >>> x3 30.0 modelオブジェクトのObjVal属性 変数オブジェクト x1 のVarName属性で名前 X属性で最適解

Pythonによるモデル化 (５) 最適解の出力（別の方法） for v in model.getVars(): print(v.VarName,v.X) >>>出力 x1 10.0 x2 10.0 x3 30.0 modelオブジェクトのgetVars()メソッド vに変数が順番に入るforループ 変数オブジェクト v のVarName属性で名前 X属性で最適解

線形最適化の問題 テキスト 265ページからの 問題 17, 18, 19, 20 から１つを選んで解いてみよう！ （残りの問題は宿題）

0-1ナップサック問題 ナップサックの重量を超えないように，なるべく価値の高い クマさんをもって逃げるには，どれを持っていく？ （テキストの例題は，容量制約もある多制約ナップサック問題）

整数計画問題としての定式化 クマさんの集合 J クマさんjの価値 vj 変数 クマさんjの重量 aj クマさんiをもっていくとき xj =1 ナップサックの重量制限 b 変数 クマさんiをもっていくとき xj =1 そうでないとき xj =0

Pythonによるモデル化 (1) モジュールの読み込み（mypulpは商用ソルバーGurobiと同じ文法） from mypulp import * モデルオブジェクトの生成 model = Model(“適当な名前")

Pythonによるモデル化 (2) 変数オブジェクトの追加 x1 = model.addVar(vtype=‘B’, name='x1') x2 = model.addVar(vtype=‘B’, name='x2') x3 = model.addVar(vtype=‘B’, name='x3') x3 = model.addVar(vtype=‘B’, name='x4') モデルの更新 model.update() 引数 vtype で変数の種類を指定 C : Continuous 実数（規定値） I : Integer 整数 B : Binary 0-1（2値；バイナリ）

Pythonによるモデル化 (3) 目的関数の設定 model.setObjective(16*x1+19*x2+23*x3+28*x4, GRB.MAXIMIZE) 制約の追加 model.addConstr(2*x1+3*x2+4*x3+5*x4 <=7) 最適化 model.optimize() setObjectiveメソッドで 目的関数を設定 最大化を指定 GRB は Gurobi の略 addConstrメソッドで 制約を追加 optimizeメソッドで求解

Pythonによるモデル化 (4) 目的関数値の出力 print('Opt. Value=', model.ObjVal) >>> Optimal value= 44.0 最適解の出力 for x in model.getVars(): print( x.VarName,x.X) >>> x1 1.0 x2 0.0 x3 0.0 x4 1.0 modelオブジェクトのObjVal属性 modelオブジェクトのgetVars()メソッド vに変数が順番に入るforループ 変数オブジェクト v のVarName属性で名前 X属性で最適解

格好いいモデル化 (1) データの準備 J =[1,2,3,4] #クマさんの番号のリスト v = {1:16, 2:19, 3:23, 4:28} #クマさんの価値の辞書 a = {1:2, 2:3, 3:4, 4:5} #クマさんの重さの辞書 モデルと変数の準備 model = Model(‘mkp’) x = {} #空の辞書を準備 for j in J: #jに1,2,3,4を順に入れたループ 　 x[j] = model.addVar(vtype=‘B’, name=‘x({0})’.format(j)) #0-1変数 x[1],x[2],x[3],x[4] を生成して辞書に保管 model.update() “文字列{0}”.format(何々) で “文字列何々” になる！ 詳細はテキスト429ページ

格好いいモデル化 (2) 変数の和を計算するためのsumの高速化版 目的関数の設定 model.setObjective(quicksum(v[j]*x[j] for j in J), GRB.MAXIMIZE) 制約の追加 model.addConstr(quicksum(a[j]*x[j] for j in J) <=7) 内包表記の復習 リスト内包 [ 2*i for i in range(3) ] =>[0,2,4] のリスト sum関数 sum( [0,2,4] ) => 0+2+4 = 6　 ジェネレータ内包 sum( 2*i for i in range(3) ) => 6を返す

整数最適化の問題 テキスト 266ページの問題 22を解いてみよう！ 余力があれば他の問題も解いてみよう！加点します！！

0-1変数を用いた論理条件の定式化 (1) 𝑥1+𝑥4≤1 離接制約 （ＡかBの何れか） 0-1ナップサック問題 極小クマ（x1 もしくは x[1]）か大クマ（x4 もしくは x[4]）の何 れかしかもって逃げることができない！ 𝑥1+𝑥4≤1

0-1変数を用いた論理条件の定式化 (2) 𝑥1≤𝑥2 if A then B （もしＡならB） 0-1ナップサック問題 もし極小クマ（x1 もしくは x[1]）をもって逃げるなら， 小クマ（x2 もしくは x[2]）ももって逃げなければならない！ 𝑥1≤𝑥2

最大値の最小化 XとYの大きい方を最小化 min max { X,Y} 補助変数 Z を用いる！ 𝑍 最小化 𝑋≤𝑍 条件 最小化　 𝑍 𝑋≤𝑍 条件　 ZはX,Yの大きい方になる！ 𝑌≤𝑍

絶対値の最小化 𝑌=|𝑋| 𝑍 𝑍≥𝑋 𝑍≥−𝑋 𝑌=𝑋 𝑌=−𝑋 𝑌+𝑍 𝑌≥0 𝑍≥0 Xの絶対値の最小化 min |X| 方法1：補助変数 Z=|X| を導入 方法２：補助変数 Y,Z を導入 X=Y-Z 　　　 　 𝑍 最小化　 𝑍≥𝑋 𝑍≥−𝑋 条件　 𝑌=𝑋 𝑌=−𝑋 最小化　 𝑌+𝑍 Xが正のときは Y>0 Xが負のときは Z>0 条件　 𝑌≥0 𝑍≥0

半連続変数 10 𝑦 ≤𝑥 ≤20 𝑦 実数変数 x が 0もしくは[10,20]の範囲内 0-1変数 ｙを導入 yが0のときにはxは0 yが1のときにはxは 10以上20以下 10 𝑦 ≤𝑥 ≤20 𝑦

工夫が必要な定式化の問題 テキスト 272ページの問題 28, 29, 30を解いてみよう！ 余力があれば他の問題も解いてみよう！加点します！！

非線形最適化 (SciPyのoptimizeでも解ける） Rosenbrock関数の最小化 𝑓 𝑥 = 𝑛=1 𝑛−1 100( 𝑥 𝑖 2 − 𝑥 𝑖+1 2 +(1− 𝑥 𝑖 ) 2 ） ２次元(n=2)のときの図 x=(1,1)で最適値 0

Pythonによるモデル化 (1) 非線形最適化ソルバー OpenOpt (myopenoptモジュール）を用いる． モジュールの読み込み（myopenoptはmypulp, Gurobiと同じ文法） from myopenopt import * モデルと変数の準備 model=Model() n=5 x={} #空の辞書を準備 for i in range(n): x[i]=model.addVar(name='x[{0}]'.format(i), init =3.0 ) 非線形の場合には，引数 init で初期値を入れる！

Pythonによるモデル化 (2) 目的関数（最小化）の設定 model.setObjective( sum( 100*(x[i]**2-x[i+1])**2 + (1-x[i])**2 for i in range(n-1)), GRB.MINIMIZE) 最適化 model.optimize(plot=True) #引数plotをTrueにすると図が出る！ 解の出力 for v in model.getVars(): print(v.VarName,v.X)

最適化の結果 x[0] 1.00000791776 x[1] 1.00003710501 x[2] 1.00004699024 x[3] 1.00013621119 x[4] 1.00028715063 ------------------------- OpenOpt 0.5625 --------------- problem: unnamed type: NLP goal: minimum solver: ralg iter objFunVal log10(maxResidual) 0 1.442e+04 -100.00 10 4.678e+01 -100.00 20 1.181e-01 -100.00 30 4.780e-03 -100.00 40 4.019e-04 -100.00 50 8.077e-06 -100.00 58 3.415e-07 -100.00 istop: 4 (|| F[k] - F[k-1] || < ftol) Solver: Time Elapsed = 0.3600000000000003 CPU Time Elapsed = 0.6093610000000025 Plotting: Time Elapsed = 2.07 CPU Time Elapsed = 3.6906389999999973 objFuncValue: 3.4150106e-07 (feasible, MaxResidual = 0)

Weber問題（テキスト88ページ） ７つの家の住民が歩く距離の和が最小になる地点に井戸を掘りたい． x座標 x y座標 y 人数 w 1 24 54 2 60 63 3 84 4 23 100 5 48 6 15 64 7 52 74 ７つの家の住民が歩く距離の和が最小になる地点に井戸を掘りたい． ただし，中心 (50,50) から半径 40 の領域は砂漠なので，水は出ない．

定式化 (50−𝑋) 2 + 50−𝑌 2 ≥ 40 2 井戸の場所を変数 (X,Y) とする． 𝑖 𝑤 𝑖 ( 𝑥 𝑖 −𝑋) 2 + ( 𝑦 𝑖 −𝑌) 2 最小化　 条件　 (50−𝑋) 2 + 50−𝑌 2 ≥ 40 2

Pythonによるモデル化（Weber問題） from myopenopt import * n = 7 I = range(1,n+1)　 #点（家）の番号を1から7とする x = {1: 24, 2: 60, 3: 1, 4: 23, 5: 84, 6: 15, 7: 52} #家のx座標を表す辞書 y = {1: 54, 2: 63, 3: 84, 4: 100, 5: 48, 6: 64, 7: 74} #家のy座標を表す辞書 w = {1: 2, 2: 1, 3: 2, 4: 3, 5: 4, 6: 5, 7: 4} #家の人数を表す辞書 model=Model() X=model.addVar(name='X') Y=model.addVar(name='Y') model.addConstr( (X-50)**2+(Y-50)**2>=1600 ) model.setObjective(sum( w[i]* sqrt( (x[i]-X)**2 +(y[i]-Y)**2 ) for i in I), GRB.MINIMIZE ) model.optimize(plot=True) for v in model.getVars(): print(v.VarName,v.X)

最適化の結果（Weber問題） ------------------------- OpenOpt 0.5625 --------------- problem: unnamed type: NLP goal: minimum solver: ralg iter objFunVal log10(maxResidual) 0 1.758e+03 -100.00 10 6.846e+02 -100.00 20 6.779e+02 -100.00 30 6.777e+02 -100.00 40 6.777e+02 -100.00 50 6.777e+02 -100.00 52 6.777e+02 -100.00 istop: 3 (|| X[k] - X[k-1] || < xtol) Solver: Time Elapsed = 0.7300000000000002 CPU Time Elapsed = 0.9752789999999987 Plotting: Time Elapsed = 1.99 CPU Time Elapsed = 3.874721000000001 objFuncValue: 677.71943 (feasible, MaxResidual = 0) X 16.5375106626 Y 71.9148765872

非線形最適化の問題 テキスト 276ページからの問題 35-38から１つを選んで解いてみ よう！(問題36はちょっと難しい） 余力があれば他の問題も解いてみよう！加点します！！