#27 最適化計算の実行
宮崎 稔也
前回の投稿（#26）で、最適化計算のライブラリPyomoで、変数をベクトル化してみたことを報告しました。
https://note.com/toshiyamiyazaki/n/nd8dbf61b1b02
今回はその続き。
実際に最適化計算を行いました。
まずは、コード全体を掲載します。
import pyomo.environ as pyo
from pyomo.opt import SolverFactory
# ConcreteModelクラスを用いて問題のオブジェクトを作成
model = pyo.ConcreteModel()
# リストで指定した内容をキーワードとして保持
df_index = df_price_per_ton.index
NUM = len(df_index)
model.ITEMS = list(df_index)
uu = {df_index[i]:df_price_per_ton.iat[i, 0] for i in range(NUM)}
def initializeU(model, i):
'''変数の初期化'''
return uu[i]
model.u = pyo.Var(model.ITEMS, domain=pyo.NonNegativeReals, initialize=initializeU)
for x in range(len(df_index)):
# 投入量を計算
input = 0 # 投入分
output = 0 # 産出分
fuel = 0 # 燃料使用分
waste = 0 # 廃棄物発生量
for i in range(187):
for k in model.ITEMS:
if k == '0131': # 0131（農業サービス）からの投入は考慮しない
continue
if k == df_iotable.index[i]:
input += df_iotable.iat[i, x] * 1000000 / model.u[k]
# 燃料使用分を計算
if k == '2121': # 石炭製品の燃料使用分
fuel += df_iotable.iat[i, x] * 1000000 / model.u[k]
# 石油化学系基礎製品への石油製品の投入ではない場合
if k == '2111':
if df_iotable.columns[x] != '2031':
fuel += df_iotable.iat[i, x] * 1000000 / model.u[k] # 石油製品の燃料使用分
# 産出量を計算
for i in range(187):
if i == x:
continue
if model.ITEMS[x] == '0131': # 0131（農業サービス）の産出は考慮しない
output = 0
else:
output += df_iotable.iat[x, i] * 1000000 / model.u[k]
# 廃棄物発生量を計算
waste = input - output - fuel
# 計算結果をdf_estimateに格納
df_estimate.iat[x, 0] = waste
df_estimate.iat[x, 1] = input
df_estimate.iat[x, 2] = output
df_estimate.iat[x, 3] = fuel
# 目的関数を定義してruleに渡す
def ObjRule(model):
z = 0
slist = ['農業', '林業', '鉱業', '印刷', '皮革', '非鉄金属', '業務用機器', '情報通信機器','輸送用機器']
for sector, group in df_estimate.groupby('sector'):
if sector in slist:
continue
else:
z += (group['Wx'].sum() - df_statistics.loc[sector, 'Sy']) ** 2
return z
model.OBJ = pyo.Objective(rule = ObjRule,
sense = pyo.minimize)
# 制約条件を定義
#model.CONST = pyo.Constraint(expr = model.u['0131'] == 0 )
model.CONST = pyo.ConstraintList()
for i in model.ITEMS:
if i == '0131':
model.CONST.add(model.u[i] == 0)
else:
model.CONST.add(model.u[i] >= 100)
# ソルバを指定して最適化
opt = pyo.SolverFactory('ipopt')
opt.options['linear_solver'] = 'ma27'
# 最適化計算のログ出力の詳細度合い (print_level) を設定
opt.options["print_level"] = 4
res = opt.solve(model, tee=True) # 最適化計算を実行
print(f"評価関数：{model.OBJ()}")
for i in model.ITEMS:
print(f"u{i}: {model.u[i].value:.0f}") # 重量単価を1つずつ表示
print(res)問題のオブジェクを作成
まず、ConcreteModelクラスを用いて問題のオブジェクトを作成します。
# ConcreteModelクラスを用いて問題のオブジェクトを作成
model = pyo.ConcreteModel()変数の定義
Varクラスを用いて変数を定義します。
このとき、model.ITEMSというリストを引数に与えることで、ベクトル化された変数を定義できます。
model.ITEMSには、df_price_per_ton（重量単価初期値一覧表）のインデックス（code）がリスト形式で格納されています。
次に、uuに、キー（key）がcode、値（value）が重量単価初期値のペアを辞書形式で格納していきます。
続いて、関数InitializeUを定義し、変数の初期値を重量単価初期値で与えていきます。
# リストで指定した内容をキーワードとして保持
df_index = df_price_per_ton.index
NUM = len(df_index)
model.ITEMS = list(df_index)
uu = {df_index[i]:df_price_per_ton.iat[i, 0] for i in range(NUM)}
def initializeU(model, i):
'''変数の初期化'''
return uu[i]
model.u = pyo.Var(model.ITEMS, domain=pyo.NonNegativeReals, initialize=initializeU)廃棄物発生量の計算
各産業の廃棄物発生量を、先ほど定義した変数を用いて計算していきます。
for x in range(len(df_index)):
# 投入量を計算
input = 0 # 投入分
output = 0 # 産出分
fuel = 0 # 燃料使用分
waste = 0 # 廃棄物発生量
for i in range(187):
for k in model.ITEMS:
if k == '0131': # 0131（農業サービス）からの投入は考慮しない
continue
if k == df_iotable.index[i]:
input += df_iotable.iat[i, x] * 1000000 / model.u[k]
# 燃料使用分を計算
if k == '2121': # 石炭製品の燃料使用分
fuel += df_iotable.iat[i, x] * 1000000 / model.u[k]
# 石油化学系基礎製品への石油製品の投入ではない場合
if k == '2111':
if df_iotable.columns[x] != '2031':
fuel += df_iotable.iat[i, x] * 1000000 / model.u[k] # 石油製品の燃料使用分
# 産出量を計算
for i in range(187):
if i == x:
continue
if model.ITEMS[x] == '0131': # 0131（農業サービス）の産出は考慮しない
output = 0
else:
output += df_iotable.iat[x, i] * 1000000 / model.u[k]
# 廃棄物発生量を計算
waste = input - output - fuel
# 計算結果をdf_estimateに格納
df_estimate.iat[x, 0] = waste
df_estimate.iat[x, 1] = input
df_estimate.iat[x, 2] = output
df_estimate.iat[x, 3] = fuel目的関数の定義
Objectiveクラスを用いて目的関数を定義します。
ObjRule関数で目的関数を記述し、Objectiveクラスのruleに渡しています。
slistには、目的関数から除外するsector（部門分類）として、農業、林業、 鉱業、印刷、皮革、非鉄金属、業務用機器、情報通信機器、輸送用機器を指定します。
その理由として、
・農業部門では、動物の糞尿のみが統計値に計上されていること
・農業以外の部門では、廃棄物発生量が0となっていること
が挙げられます。
# 目的関数を定義してruleに渡す
def ObjRule(model):
z = 0
slist = ['農業', '林業', '鉱業', '印刷', '皮革', '非鉄金属', '業務用機器', '情報通信機器','輸送用機器']
for sector, group in df_estimate.groupby('sector'):
if sector in slist:
continue
else:
z += (group['Wx'].sum() - df_statistics.loc[sector, 'Sy']) ** 2
return z
model.OBJ = pyo.Objective(rule = ObjRule,
sense = pyo.minimize)制約条件の定義
Constraintクラスで制約条件を設定します。
# 制約条件を定義
model.CONST = pyo.ConstraintList()
for i in model.ITEMS:
if i == '0131':
model.CONST.add(model.u[i] == 0)
else:
model.CONST.add(model.u[i] >= 100)ソルバを指定
SolverFactoryでソルバを指定します。ここではIPOPTを用います。
オプションで、linearsolverをHSLのma27に指定します。
最初、linearsolverがデフォルトのMUMPS (MUltifrontal Massively Parallel sparse direct Solver)で最適化計算を実行したのですが、最適化計算がうまくいきませんでした。
そこで、オプションで、lineasolverをHSL MA27に指定しました。
# ソルバを指定して最適化
opt = pyo.SolverFactory('ipopt')
opt.options['linear_solver'] = 'ma27'HSLのインストールに四苦八苦しましたが、以下のサイトを参考にして、何とかできました^^;
Installation — do-mpc 4.0.0 documentation
1.Obtain the HSL shared library. Choose the personal licence.
2.Unpack the archive and copy its content to a destination of your choice. (e.g. /home/username/Documents/coinhsl/)
3.Rename libcoinhsl.so to libhsl.so. CasADi is searching for the shared libraries under a depreciated name.
4.Locate your .bashrc file on your home directory (e.g. /home/username/.bashrc)
5.Add the previously created directory to your LD_LIBRARY_PATH, by adding the following line to your .bashrc
export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:/home/ffiedler/Documents/coinhsl/lib"6.Install libgfortran with Anaconda:
conda install -c anaconda libgfortran
最適化計算のログ出力の詳細度合いを指定
最適化計算のログ出力の詳細の度合いを、print_levelで設定できます。
今回は4に設定。0～12の範囲で設定でき、数字が大きいほど詳細になります。
# 最適化計算のログ出力の詳細度合い (print_level) を設定
opt.options["print_level"] = 4最適化
SolverFactoryオブジェクトのsolveメソッドにmodelを指定して実行すると、最適化が行わます。resには最適化の結果が格納されます。
また、solveでtee=Trueとすると、ソルバの出力が表示されます。
res = opt.solve(model, tee=True) # 最適化計算を実行結果の表示
目的関数と重量単価【最適化後】を表示させます。
print(f"評価関数：{model.OBJ()}")
for i in model.ITEMS:
print(f"u{i}: {model.u[i].value:.0f}") # 重量単価を1つずつ表示
また、opt.solve()の戻り値をprintで表示させると、以下のような出力を得ます。
print(res)
Follow me!
