Gekko小场景问题可行、大场景问题不可行的排查与解决

问题描述 #

使用Python的Gekko求解带指示变量的优化问题：指示变量I_s在θ_s>0时取1，θ_s=0时取0。当场景数N=10时，求解得到所有θ_s=0，符合预期；但设置N=100或200时求解器无法找到可行解。需要确认大N（如200）下θ_s是否仍全为0，并解决求解失败的问题。

原代码 #

# Import package
from gekko import GEKKO
import numpy as np

# Define parameters
P_CO = 600                      # $/tonCO
beta_CO2 = 1                    # no unit
P_CO2 = 80                      # $/tonCO2eq
E_ref = 3.1022616               # tonCO2eq/tonCO
E_dir = -1.600570692            # tonCO2eq/tonCO
E_indir_others = 0.3339226804   # tonCO2eq/tonCO
E_indir_elec_cons = 18.46607256 # GJ/tonCO
C1_CAPEX = 285695               # no unit
C2_CAPEX = 188.42               # no unit
C1_FOX = 82282                  # no unit
C2_FOX = 24.094                 # no unit
C1_ROX = 4471.5                 # no unit
C2_ROX = 96.034                 # no unit
C1_UOX = 7934.9                 # no unit
C2_UOX = 986.9                  # no unit
r = 0.08                        # discount rate
N = 10                         # number of scenarios
T = 30                          # total time period
GWP_init = 0.338723235          # 2020 Electricity GWP in EU 27 countries
theta_max = 1600000             # Max capacity

# Function to make GWP_EU matrix (TxN matrix)
def Electricity_GWP(GWP_init, n_years, num_episodes):

    GWP_mean = 0.36258224*np.exp(-0.16395611*np.arange(1, n_years+2)) + 0.03091272
    GWP_mean = GWP_mean.reshape(-1,1)
    GWP_Yearly = np.tile(GWP_mean, num_episodes) 

    noise = np.zeros((n_years+1, num_episodes))
    stdev2050 = GWP_mean[-1] * 0.25 
    stdev = np.arange(0, stdev2050 * (1 + 1/n_years), stdev2050/n_years)

    for i in range(n_years+1):
        noise[i,:] = np.random.normal(0, stdev[i], num_episodes) 

    GWP_forecast = GWP_Yearly + noise 

    return GWP_forecast

GWP_EU = Electricity_GWP(GWP_init, T, N) # (T+1)*N matrix
GWP_EU = GWP_EU[1:,:] # T*N matrix

print(np.shape(GWP_EU))

# Build Gekko model
m = GEKKO(remote=False)
theta = m.Array(m.Var, N, lb=0, ub=theta_max)
I = m.Array(m.Var, N, lb=0, ub=1, integer=True)
demand = np.ones((T,1))
demand[0] = 8031887.589
for k in range(1,11):
    demand[k] = demand[k-1] * 1.026 
for k in range(11,21):
    demand[k] = demand[k-1] * 1.016
for k in range(21,T):
    demand[k] = demand[k-1] * 1.011 
demand = 0.12 * demand
demand = np.tile(demand, N) # T*N matrix

print(np.shape(demand))

m3 = [[m.min3(demand[t,s],theta[s]) for t in range(T)] for s in range(N)]
obj = m.sum([sum([((1/(1+r))**(t+1))*((P_CO*m3[s][t]) \
            + (beta_CO2*P_CO2*m3[s][t]*(E_ref-E_dir-E_indir_others-E_indir_elec_cons*GWP_EU[t,s])) \
            - (C1_CAPEX*I[s]+C2_CAPEX*theta[s]+C1_FOX*I[s]+C2_FOX*theta[s])\
            - (C1_ROX*I[s]+C2_ROX*m3[s][t]+C1_UOX*I[s]+C2_UOX*m3[s][t])) for t in range(T)]) for s in range(N)]
for i in range(N):
    m.Equation(theta[i]<=1000000*I[i])
    m.Equation(-theta[i]<1000000*(1-I[i]))

# obj = m.sum([m.sum([((1/(1+r))**(t+1))*((P_CO*m.min3(demand[t,s], theta[s])) \
#             + (beta_CO2*P_CO2*m.min3(demand[t,s], theta[s])*(E_ref-E_dir-E_indir_others-E_indir_elec_cons*GWP_EU[t,s])) \
#             - (C1_CAPEX+C2_CAPEX*theta[s]+C1_FOX+C2_FOX*theta[s])-(C1_ROX+C2_ROX*m.min3(demand[t,s], theta[s])+C1_UOX+C2_UOX*m.min3(demand[t,s], theta[s]))) for t in range(T)]) for s in range(N)])
m.Maximize(obj/N)
m.solve(disp=True)   

# s = m.sum(m.sum(((1/(1+r))**(t+1))*((P_CO*m.min3(demand[t,s], theta[s])) \
#             + beta_CO2*P_CO2*m.min3(demand[t,s], theta[s])*(E_ref-E_dir-E_indir_others-E_indir_elec_cons*GWP_EU[t,s]) \
#             - (C1_CAPEX + C2_CAPEX*theta[s]) - (C1_FOX + C2_FOX*theta[s]) - (C1_ROX + C2_ROX*m.min3(demand[t,s], theta[s])) - (C1_UOX + C2_UOX*m.min3(demand[t,s], theta[s])))
#     for s in range(N)) for t in range(T))/N
         
print(theta)

解决步骤与验证方法 #

1. 先验证θ全0是否为大N下的最优解 #

手动固定所有θ_s=0、I_s=0，代入目标函数计算收益值。再针对单个场景，计算θ_s>0时的净收益：比较运营阶段累计折现收益与固定成本（C1_CAPEX+C1_FOX+C1_ROX+C1_UOX）。如果所有场景下运营收益都无法覆盖固定成本，那么θ_s=0确实是全局最优解，此时求解失败是数值层面的问题。

2. 修正指示变量约束的严谨性 #

原约束中-theta[i]<1000000*(1-I[i])等价于theta[i] > -1000000*(1-I[i])，但theta[i]已经设置了下限0，这个约束完全无效。替换为标准的大M约束，强制θ与I的逻辑关系：

M = theta_max  # 用theta的最大值作为大M参数
epsilon = 1e-3  # 避免θ极小值导致I不触发
for i in range(N):
    m.Equation(theta[i] <= M * I[i])       # θ>0时必须I=1
    m.Equation(theta[i] >= epsilon * I[i]) # I=1时θ至少取epsilon，避免数值歧义

3. 提升模型数值稳定性 #

变量量级差异过大（theta是1e6级别，固定成本是1e5级别）会导致求解器数值困难，对theta进行缩放：

scale_factor = 1e5
theta_scaled = m.Array(m.Var, N, lb=0, ub=theta_max/scale_factor)
# 后续所有涉及theta的地方替换为theta_scaled * scale_factor
# 约束修改为：
m.Equation(theta_scaled[i] * scale_factor <= M * I[i])
m.Equation(theta_scaled[i] * scale_factor >= epsilon * I[i])

4. 调整求解器参数 #

针对大规模MIP问题，增加迭代次数、调整容忍度：

m.options.SOLVER = 1  # 指定用APOPT求解器（MIP默认）
m.options.MAX_ITER = 10000  # 增加最大迭代次数
m.options.RTOL = 1e-6  # 相对容忍度
m.options.ATOL = 1e-6  # 绝对容忍度
m.options.MIP_GAP = 1e-4  # MIP最优性间隙，允许一定误差提升求解速度

5. 拆分独立场景并行求解 #

观察模型结构，每个场景的θ_s和I_s是独立的（目标函数是各场景的平均，约束无交叉），可以拆分成N个独立子问题，每个子问题求解单个场景的最优解，并行处理后汇总结果，大幅降低单问题规模：

# 示例：单个场景求解函数
def solve_single_scenario(s):
    m_single = GEKKO(remote=False)
    theta_s = m_single.Var(lb=0, ub=theta_max)
    I_s = m_single.Var(lb=0, ub=1, integer=True)
    # 构建单个场景的目标函数
    obj_single = 0
    for t in range(T):
        m3_t = m_single.min3(demand[t,s], theta_s)
        discount = (1/(1+r))**(t+1)
        revenue = P_CO * m3_t + beta_CO2 * P_CO2 * m3_t * (E_ref - E_dir - E_indir_others - E_indir_elec_cons * GWP_EU[t,s])
        cost_fixed = C1_CAPEX*I_s + C2_CAPEX*theta_s + C1_FOX*I_s + C2_FOX*theta_s
        cost_var = C1_ROX*I_s + C2_ROX*m3_t + C1_UOX*I_s + C2_UOX*m3_t
        obj_single += discount * (revenue - cost_fixed - cost_var)
    m_single.Maximize(obj_single)
    # 添加约束
    M = theta_max
    epsilon = 1e-3
    m_single.Equation(theta_s <= M * I_s)
    m_single.Equation(theta_s >= epsilon * I_s)
    m_single.solve(disp=False)
    return theta_s.value[0], I_s.value[0]

# 批量求解所有场景
theta_results = []
I_results = []
for s in range(N):
    theta_s, I_s = solve_single_scenario(s)
    theta_results.append(theta_s)
    I_results.append(I_s)

print("Theta results:", theta_results)

6. 验证随机场景的一致性 #

由于GWP_EU包含随机噪声，大N下可能存在少数场景θ_s>0的情况。可以先固定随机种子（np.random.seed(42)），确保每次生成的GWP_EU一致，便于排查差异。

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