#32  移輸出や移輸入、消費や固定資本、在庫を考慮に入れて最適化を実行

#30で、移輸出や移輸入、消費や固定資本、在庫を考慮に入れてやり直すことについて言及しました。

https://note.com/toshiyamiyazaki/n/n9303b5f4fa79

今回は、上記のことを考慮して、最適化計算をやり直しました。
なお、最適化計算のコードの詳細は、#27で触れています。

https://note.com/toshiyamiyazaki/n/n5857a61a2a36

手直したコードは以下のようになります。

import pyomo.environ as pyo
from pyomo.opt import SolverFactory
import numpy as np

# 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 # 廃棄物発生量
   from_other_industry = 0 # 他産業からの投入
   same_industry = 0 # 同産業
   to_other_industry = 0 # 他産業への産出
   iyunyu = 0 # 移輸入
   iyushutsu = 0 # 移輸出
   consume = 0 # 消費
   fixed_capital = 0 # 固定資本
   stock = 0 # 在庫
   Ux = model.u[model.ITEMS[x]]
   xcode = model.ITEMS[x]
   for i in range(187):
       for k in model.ITEMS:
           Ui = model.u[k]
           if k == '0131': # 0131（農業サービス）からの投入は考慮しない
               continue
           if k == df_iotable.index[i]:
               from_other_industry += df_iotable.iat[i, x] * 1000000 / Ui
               # 燃料使用分を計算
               if k == '2121': # 石炭製品の燃料使用分
                   fuel += df_iotable.iat[i, x] * 1000000 / Ui
               # 石油化学系基礎製品への石油製品の投入ではない場合
               if k == '2111':
                   if df_iotable.columns[x] != '2031':
                       fuel += df_iotable.iat[i, x] * 1000000 / Ui # 石油製品の燃料使用分

   # 産出量を計算
   if model.ITEMS[x] == '0131': # 0131（農業サービス）の産出は考慮しない
       same_industry = 0
       to_other_industry = 0
       consume = 0
       fixed_capital = 0
       stock = 0
       iyushutsu = 0
       iyunyu = 0
   else:
       for i in range(187):
           if i == x:
               same_industry = df_iotable.iat[x, i] * 1000000 / Ux
           else:
               to_other_industry += df_iotable.iat[x, i] * 1000000 / Ux

       ctuple = ('7111', '7211','7212', '7311', '7321')
       for c in ctuple:
           consume += df_iotable.loc[xcode, c] * 1000000 / Ux

       ftuple = ('7411', '7511')
       for f in ftuple:
           fixed_capital += df_iotable.loc[xcode, f] * 1000000 / Ux

       stock = df_iotable.loc[xcode, '7611'] * 1000000 / Ux

       iyushutsu = df_iotable.loc[xcode, '8100'] * 1000000 / Ux

       iyunyu = -(df_iotable.loc[xcode, '8700']) * 1000000 / Ux

   input = from_other_industry + iyunyu
   output = same_industry + to_other_industry + iyushutsu + consume + fixed_capital + stock
   waste = input - output - fuel

   # 計算結果をdf_under_estimationに格納
   df_under_estimation.iat[x, 0] = input
   df_under_estimation.iat[x, 1] = output
   df_under_estimation.iat[x, 2] = fuel
   df_under_estimation.iat[x, 3] = waste

# 目的関数を定義してruleに渡す
def ObjRule(model):
   z = 0
   slist = ['印刷', '皮革', '非鉄金属', '業務用機器', '情報通信機器','輸送用機器']
   for sector, group in df_under_estimation.groupby('sector'):
       if sector in slist or sector == np.NaN:
           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.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)

最適化後の重量単価の表示（一部）は、以下のようになります。

Follow me!