我正在尝试计算各种约束的影子价格，并将它们按照 c[i]:sp[i] 的形式放入字典中，其中 c 是特定约束的名称，sp 是数字影子价格值。我能够在本地生成“apm_lam.txt”，但数组的形状是 506，这远远超过了应该添加到模型中的约束数量（我数了 60-某事）。
是否可以轻松将此输出文件中的影子价格与 Gekko 模型中的实际命名约束联系起来？
import numpy as np
import pandas as pd
from gekko import GEKKO
m = GEKKO(remote=False)
m.options.NODES = 3
m.options.IMODE = 3
m.options.MAX_ITER = 1000
lnuc_weeks = [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0]
min_promo_price = [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,3]
max_promo_price = [3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5,3.5, 3.5, 3.5, 3.5, 3.5, 3.5]
base_srp = [3.48, 3.48, 3.48, 3.48, 3.0799, 3.0799, 3.0799, 3.0799,3.0799, 3.0799, 3.0799, 3.0799, 3.0799, 3.0799, 3.0799, 3.0799, 3.0799, 3.0799, 3.0799]
lnuc_min_promo_price = 1.99
lnuc_max_promo_price = 1.99
coeff_fedi = [0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589,0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589]
coeff_feao = [0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995]
coeff_diso = [0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338]
sumproduct_base = [0.20560305, 0.24735297, 0.24957423, 0.23155435, 0.23424058,0.2368096 , 0.27567109, 0.27820648, 0.2826393 , 0.28660598, 0.28583971, 0.30238505, 0.31726649, 0.31428312, 0.31073792, 0.29036779, 0.32679041, 0.32156337, 0.24633734]
neg_ln = [[0.14842000515],[0.14842000512],[0.14842000515],[0.14842000512],[-0.10407483058],[0.43676249024],[0.43676249019],[0.43676249024],[0.43676249019],[0.43676249024],[0.43676249019], [0.026284840258],[0.026284840291],[0.026284840258],[0.026284840291], [0.026185109811],[0.026284840258],[0.026284840291],[0.026284840258]]
neg_ln_ppi_coeff = [1.22293879, 1.22293879, 1.22293879, 1.22293879, 1.22293879,1.22293879, 1.22293879, 1.22293879, 1.22293879, 1.22293879, 1.22293879, 1.22293879, 1.22293879, 1.22293879, 1.22293879,1.22293879, 1.22293879, 1.22293879, 1.22293879]
base_volume = [124.38, 193.2, 578.72, 183.88, 197.42, 559.01, 67.68, 110.01,60.38, 177.11, 102.65, 66.02, 209.83, 81.22, 250.44, 206.44, 87.99, 298.95, 71.07]
week = pd.Series([13, 14, 17, 18, 19, 26, 28, 33, 34, 35, 39, 42, 45, 46, 47, 48, 50, 51, 52])
n = 19
cons = []
shadow = []
x1 = m.Array(m.Var,(n), integer=True) #LNUC weeks
i = 0
for xi in x1:
xi.value = lnuc_weeks[i]
xi.lower = 0
xi.upper = lnuc_weeks[i]
i += 1
x2 = m.Array(m.Var,(n)) #Blended SRP
i = 0
for xi in x2:
xi.value = 5
m.Equation(xi >= m.if3((x1[i]) - 0.5, min_promo_price[i], lnuc_min_promo_price))
m.Equation(xi <= m.if3((x1[i]) - 0.5, max_promo_price[i], lnuc_max_promo_price))
i += 1
x3 = m.Array(m.Var,(n), integer=True) #F&D
x4 = m.Array(m.Var,(n), integer=True) #FO
x5 = m.Array(m.Var,(n), integer=True) #DO
x6 = m.Array(m.Var,(n), integer=True) #TPR
#Default to F&D
i = 0
for xi in x3:
xi.value = 1
xi.lower = 0
xi.upper = 1
i += 1
i = 0
for xi in x4:
xi.value = 0
xi.lower = 0
xi.upper = 1
i += 1
i = 0
for xi in x5:
xi.value = 0
xi.lower = 0
xi.upper = 1
i += 1
i = 0
for xi in x6:
xi.value = 0
xi.lower = 0
xi.upper = 1
i += 1
x7 = m.Array(m.Var,(n), integer=True) #Max promos
i = 0
for xi in x7:
xi.value = 1
xi.lower = 0
xi.upper = 1
i += 1
x = [x1,x2,x3,x4,x5,x6,x7]
neg_ln=[m.Intermediate(-m.log(x[1][i]/base_srp[i])) for i in range(n)]
total_vol_fedi =[m.Intermediate(coeff_fedi[0]+ sumproduct_base[i] + (neg_ln[i]*neg_ln_ppi_coeff[0])) for i in range(n)]
total_vol_feao =[m.Intermediate(coeff_feao[0]+ sumproduct_base[i] + (neg_ln[i]*neg_ln_ppi_coeff[0])) for i in range(n)]
total_vol_diso =[m.Intermediate(coeff_diso[0]+ sumproduct_base[i] + (neg_ln[i]*neg_ln_ppi_coeff[0])) for i in range(n)]
total_vol_tpro =[m.Intermediate(sumproduct_base[i] + (neg_ln[i]*neg_ln_ppi_coeff[0])) for i in range(n)]
simu_total_volume = [m.Intermediate((
(m.max2(0,base_volume[i]*(m.exp(total_vol_fedi[i])-1)) * x[2][i] +
m.max2(0,base_volume[i]*(m.exp(total_vol_feao[i])-1)) * x[3][i] +
m.max2(0,base_volume[i]*(m.exp(total_vol_diso[i])-1)) * x[4][i] +
m.max2(0,base_volume[i]*(m.exp(total_vol_tpro[i])-1)) * x[5][i]) + base_volume[i]) * x[6][i]) for i in range(n)]
[m.Equation(x3[i] + x4[i] + x5[i] + x6[i] == 1) for i in range(i)]
#Limit max promos
m.Equation(sum(x7)<=10)
#Enforce spacing and duration
s=1
for s2 in range(1, s+1):
for i in range(0, n-s2):
f = week[week == week[i] + s2].index
if len(f) > 0:
m.Equation(x7[i] + x7[f[0]]<=1)
m.Maximize(m.sum(simu_total_volume))
m.options.SOLVER=3
m.options.DIAGLEVEL = 2
m.solve(disp = True, debug=True)
df = pd.concat([pd.Series(week), pd.Series([i[0] for i in x7]), pd.Series([i[0] for i in simu_total_volume])], axis=1)
df.columns = ['week', 'x7', 'total_volume']
df[df['x7']>0]
m.open_folder()
lam = np.loadtxt(m.path+'/apm_lam.txt')
print(lam.shape)
有关约束的附加信息在运行目录中的
rto_4_eqn.txtremote=Falsem.solve()
对于不同的操作模式，文件的前3个字符会发生变化。
| 模式 | 方程式摘要文件名 | |
|---|---|---|
| 稳态模拟 | ss_4_eqn.txt | |
| 模型参数更新 | mpu_4_eqn.txt | |
| 实时优化 | rto_4_eqn.txt | |
| 动态模拟 | sim_4_eqn.txt | |
| 动态估计 | est_4_eqn.txt | |
| 动态控制 | ctl_4_eqn.txt | |
| 顺序模拟 | sqs_4_eqn.txt | |
| 顺序估计 | est_4_eqn.txt | |
| 顺序控制 | ctl_4_eqn.txt |
fid = open(m.path+'/rto_4_eqn.txt','r')
print(fid.read())
fid.close()
每个方程的更详细报告位于文件
************************************************
************* ACTIVE EQUATIONS *****************
************************************************
Number ID Node Horizon Unscaled Res Scaled Res Scaling Name
1 1 1 1 1.0000E-08 1.0000E-08 1.0000E+00 ss.max2_1.mpcc[1]: x[2] - x[1] = s[1] - s[2]
2 2 1 1 -1.0000E-08 -1.0000E-08 1.0000E+00 ss.max2_1.output_eqn: y = x[2] + s[2]
3 4 1 1 1.0000E-08 1.0000E-08 1.0000E+00 ss.max2_2.mpcc[1]: x[2] - x[1] = s[1] - s[2]
4 5 1 1 -1.0000E-08 -1.0000E-08 1.0000E+00 ss.max2_2.output_eqn: y = x[2] + s[2]
5 7 1 1 1.0000E-08 1.0000E-08 1.0000E+00 ss.max2_3.mpcc[1]: x[2] - x[1] = s[1] - s[2]
6 8 1 1 -1.0000E-08 -1.0000E-08 1.0000E+00 ss.max2_3.output_eqn: y = x[2] + s[2]
7 10 1 1 1.0000E-08 1.0000E-08 1.0000E+00 ss.max2_4.mpcc[1]: x[2] - x[1] = s[1] - s[2]
8 11 1 1 -1.0000E-08 -1.0000E-08 1.0000E+00 ss.max2_4.output_eqn: y = x[2] + s[2]
9 13 1 1 1.0000E-08 1.0000E-08 1.0000E+00 ss.max2_5.mpcc[1]: x[2] - x[1] = s[1] - s[2]
10 14 1 1 -1.0000E-08 -1.0000E-08 1.0000E+00 ss.max2_5.output_eqn: y = x[2] + s[2]
...
484 560 1 1 -1.1718E-12 -1.1718E-12 1.0000E+00 ss.Eqn(331): 0 = 1-((int_v204+int_v205))-slk_331
485 561 1 1 0.0000E+00 0.0000E+00 1.0000E+00 ss.Eqn(332): 0 = 1-((int_v205+int_v206))-slk_332
486 562 1 1 0.0000E+00 0.0000E+00 1.0000E+00 ss.Eqn(333): 0 = 1-((int_v207+int_v208))-slk_333
487 563 1 1 0.0000E+00 0.0000E+00 1.0000E+00 ss.Eqn(334): 0 = 1-((int_v208+int_v209))-slk_334
488 564 1 1 -1.8406E-06 -1.8406E-06 1.0000E+00 ss.Eqn(335): 0 = v438-(i209)
489 565 1 1 2.2181E-06 2.2181E-06 1.0000E+00 ss.Eqn(336): 0 = v439-(i210)
490 566 1 1 6.5732E-06 6.5732E-06 1.0000E+00 ss.Eqn(337): 0 = v440-(i211)
491 567 1 1 -5.0622E-09 -5.0622E-09 1.0000E+00 ss.Eqn(338): 0 = v441-(i212)
492 568 1 1 2.3812E-06 2.3812E-06 1.0000E+00 ss.Eqn(339): 0 = v442-(i213)
493 569 1 1 6.6581E-06 6.6581E-06 1.0000E+00 ss.Eqn(340): 0 = v443-(i214)
494 570 1 1 -1.5305E-06 -1.5305E-06 1.0000E+00 ss.Eqn(341): 0 = v444-(i215)
495 571 1 1 1.3645E-06 1.3645E-06 1.0000E+00 ss.Eqn(342): 0 = v445-(i216)
496 572 1 1 -1.3760E-06 -1.3760E-06 1.0000E+00 ss.Eqn(343): 0 = v446-(i217)
497 573 1 1 2.1610E-06 2.1610E-06 1.0000E+00 ss.Eqn(344): 0 = v447-(i218)
498 574 1 1 -1.1720E-06 -1.1720E-06 1.0000E+00 ss.Eqn(345): 0 = v448-(i219)
499 575 1 1 -9.3201E-07 -9.3201E-07 1.0000E+00 ss.Eqn(346): 0 = v449-(i220)
500 576 1 1 2.3803E-06 2.3803E-06 1.0000E+00 ss.Eqn(347): 0 = v450-(i221)
501 577 1 1 -1.1581E-06 -1.1581E-06 1.0000E+00 ss.Eqn(348): 0 = v451-(i222)
502 578 1 1 2.8314E-06 2.8314E-06 1.0000E+00 ss.Eqn(349): 0 = v452-(i223)
503 579 1 1 -1.1710E-09 -1.1710E-09 1.0000E+00 ss.Eqn(350): 0 = v453-(i224)
504 580 1 1 -1.2694E-06 -1.2694E-06 1.0000E+00 ss.Eqn(351): 0 = v454-(i225)
505 581 1 1 3.3784E-06 3.3784E-06 1.0000E+00 ss.Eqn(352): 0 = v455-(i226)
506 582 1 1 -9.4846E-07 -9.4846E-07 1.0000E+00 ss.Eqn(353): 0 = v456-(i227)
中，其中包括方程中每个变量的报告。不等式约束会转换为带有附加松弛变量的等式约束，因此
rto_4_eqn_var.txtgk0_model.apm