如何解读Gekko中影子价格数组形状并关联约束名称

我正在尝试为各类约束计算影子价格，并整理成c[i]:sp[i]形式的字典（c为约束名称，sp为影子价格数值）。目前已成功生成本地文件apm_lam.txt，但得到的数组形状为506，远多于模型中手动添加的约60个约束。请问有没有简便方法将该文件中的影子价格与Gekko模型里的实际命名约束对应起来？

解决方法 #

apm_lam.txt包含模型所有约束的影子价格，包括变量上下界、if3/max2等辅助函数自动展开的约束。要关联自定义约束，核心思路是给手动添加的约束命名并记录索引，再映射到影子价格数组。

修改后代码实现 #

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]*19
max_promo_price = [3.5]*19
base_srp = [3.48]*4 + [3.0799]*15
lnuc_min_promo_price = 1.99
lnuc_max_promo_price = 1.99
coeff_fedi = [0.022589]*19
coeff_feao = [0.02929995]*19
coeff_diso = [0.05292338]*19
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_ppi_coeff = [1.22293879]*19
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

# 存储命名约束的字典
named_constraints = {}

# 变量定义
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
    # 给价格约束命名并存入字典
    const_low = m.Equation(xi >= m.if3((x1[i]) - 0.5, min_promo_price[i], lnuc_min_promo_price))
    const_high = m.Equation(xi <= m.if3((x1[i]) - 0.5, max_promo_price[i], lnuc_max_promo_price))
    named_constraints[f"x2_low_week_{week[i]}"] = const_low
    named_constraints[f"x2_high_week_{week[i]}"] = const_high
    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

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)]

# 促销类型互斥约束
for i in range(n):
    const_promo_type = m.Equation(x3[i] + x4[i] + x5[i] + x6[i] == 1)
    named_constraints[f"promo_type_week_{week[i]}"] = const_promo_type

# 最大促销数量约束
const_max_promos = m.Equation(sum(x7) <= 10)
named_constraints["max_total_promos"] = const_max_promos

# 促销间隔约束
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:
            const_spacing = m.Equation(x7[i] + x7[f[0]] <= 1)
            named_constraints[f"promo_spacing_week_{week[i]}_and_{week[f[0]]}"] = const_spacing

m.Maximize(m.sum(simu_total_volume))

m.options.SOLVER = 3
m.options.DIAGLEVEL = 2
m.solve(disp=True, debug=True)

# 关联约束与影子价格
lam = np.loadtxt(m.path + '/apm_lam.txt')
# 获取每个约束在求解器中的索引
constraint_indices = {eq: idx for idx, eq in enumerate(m._equations)}

# 构建约束名称到影子价格的字典
shadow_price_dict = {}
for name, eq in named_constraints.items():
    idx = constraint_indices.get(eq)
    if idx is not None and idx < len(lam):
        shadow_price_dict[name] = lam[idx]

# 输出结果
print("约束名称: 影子价格")
for name, sp in shadow_price_dict.items():
    print(f"{name}: {sp}")

# 原数据分析部分
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']
print(df[df['x7']>0])

关键说明 #

1. 约束命名：所有手动添加的约束都被赋予清晰名称并存入named_constraints字典，避免匿名约束无法识别。
2. 索引映射：利用Gekko模型的_equations属性，获取每个约束在求解器中的索引，对应apm_lam.txt中的影子价格位置。
3. 过滤无关约束：最终生成的shadow_price_dict仅包含自定义命名约束的影子价格，自动排除变量上下界、辅助函数生成的冗余约束。

内容的提问来源于stack exchange，提问作者datadude558