Gekko本地模式下大规模场景优化的内存溢出问题求解

我正在求解上图所示的优化问题，并编写了如下代码。但当场景数N较大（如300）时，代码因内存溢出无法运行。尝试将m = GEKKO(remote=False)替换为m = GEKKO(remote=True)，但运行10-12小时仍未完成。现寻求在m = GEKKO(remote=False)模式下解决该问题的方法。

原始代码 #

# Import package
from gekko import GEKKO
import numpy as np

# Define parameters
P_CO = 600                      # $/tonCO
beta_CO2 = 1                    # no unit
P_CO2 = 60                      # $/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 = 1983.7                 # no unit
C2_UOX = 249.79                 # no unit
r = 0.08                        # discount rate
N = 300                         # 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)
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))

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

优化方案与修改后代码 #

核心优化思路 #

1. 减少重复变量创建：原代码每个时间步重复调用m.min3生成大量冗余中间变量，改为每个场景预定义产量变量q[t]复用结果。
2. 提取固定成本项：固定成本与时间无关，从时间循环中提取避免重复计算。
3. 调整求解器参数：本地模式切换到IPOPT求解器，限制内存使用并关闭冗余日志，提升内存效率。
4. 预计算常量：提前计算折现因子，减少循环内重复运算。

修改后代码 #

# Import package
from gekko import GEKKO
import numpy as np

# Define parameters
P_CO = 600                      # $/tonCO
beta_CO2 = 1                    # no unit
P_CO2 = 60                      # $/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 = 1983.7                 # no unit
C2_UOX = 249.79                 # no unit
r = 0.08                        # discount rate
N = 300                         # 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)[1:,:] # T*N matrix
print(np.shape(GWP_EU))

# Build Gekko model
m = GEKKO(remote=False)
m.options.SOLVER = 3  # 切换到IPOPT求解器，内存效率优于APOPT
m.options.MAX_MEMORY = 2048  # 限制内存使用（单位：MB），根据自身机器配置调整
m.options.DIAGLEVEL = 0  # 关闭诊断日志，减少内存占用与IO开销

theta = m.Array(m.Var, N, lb=0, ub=theta_max)

# 生成需求矩阵
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 * np.tile(demand, N) # T*N matrix
print(np.shape(demand))

# 预计算折现因子，避免循环内重复计算
discount_factors = [(1/(1+r))**(t+1) for t in range(T)]

# 构建目标函数
total_obj = 0
for s in range(N):
    # 为当前场景定义每个时间步的产量变量，复用min3结果
    q = [m.min3(demand[t,s], theta[s]) for t in range(T)]
    # 计算固定成本（与时间无关，仅计算一次）
    fixed_cost = C1_CAPEX + C2_CAPEX*theta[s] + C1_FOX + C2_FOX*theta[s]
    # 计算年度贡献总和
    yearly_total = 0
    for t in range(T):
        # 碳收益项
        carbon_revenue = beta_CO2 * P_CO2 * q[t] * (E_ref - E_dir - E_indir_others - E_indir_elec_cons * GWP_EU[t,s])
        # 可变成本项
        variable_cost = (C1_ROX + C1_UOX) + (C2_ROX + C2_UOX)*q[t]
        # 当前时间步的净贡献（折现后）
        period_contrib = (P_CO * q[t] + carbon_revenue - variable_cost) * discount_factors[t]
        yearly_total += period_contrib
    # 场景总贡献：年度折现贡献 - 固定成本（假设固定成本在初始时刻发生，折现因子为1）
    scenario_contrib = yearly_total - fixed_cost
    total_obj += scenario_contrib

# 最大化平均场景收益
m.Maximize(total_obj / N)
# 开启冷启动选项，加速求解
m.options.COLDSTART = 2
m.solve(disp=True)

# 输出结果
print("Optimal theta values:")
for s in range(N):
    print(f"Scenario {s+1}: theta = {theta[s].value[0]:.2f}")

额外优化建议 #

• 如果场景间存在统计相关性，可考虑场景聚合，减少场景数量后求解，再将结果推广到全场景。
• 若机器内存仍紧张，可尝试分批次求解，将N个场景分成若干组分别求解后取平均结果（需验证近似合理性）。

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