使用Python GEKKO求解优化问题时遇负自由度及目标值异常

我正在求解一个优化问题，编写了对应的Python GEKKO代码。代码可运行，但出现**负自由度（negative degrees of freedom）**问题，且得到的目标值为负数，与预期的正值不符，恳请提供解决建议。

代码 #

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

输出信息 #

(30, 10)
(30, 10)
 ----------------------------------------------------------------
 APMonitor, Version 1.0.0
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :  11
   Constants    :  0
   Variables    :  5121
   Intermediates:  0
   Connections  :  321
   Equations    :  3901
   Residuals    :  3901
 
 Number of state variables:    5121
 Number of total equations: -  3911
 Number of slack variables: -  2400
 ---------------------------------------
 Degrees of freedom       :    -1190
 
 * Warning: DOF <= 0
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:     18.61 NLPi:    5 Dpth:    0 Lvs:    0 Obj: -1.87E+09 Gap:  0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  18.619200000000003 sec
 Objective      :  -1.8677021320161405E+9
 Successful solution
 ---------------------------------------------------

问题原因及解决建议 #

1. 负自由度问题原因与修复

负自由度说明约束方程数远超决策变量数，模型被过度约束。核心问题是m.min3()函数会自动生成大量辅助变量和约束：

• 你仅定义了N=10个决策变量theta[s]，但每个m.min3(demand[t,s], theta[s])会生成至少3个约束，结合T=30、N=10的规模，总约束数直接爆炸，导致DOF为负。

修复方法：
显式定义产量变量并添加约束替代m.min3()，减少不必要的辅助变量：

# 显式定义产量变量x[t,s]
x = m.Array(m.Var, (T,N), lb=0)
# 添加产量约束：产量不超过需求，也不超过产能
for t in range(T):
    for s in range(N):
        m.Equation(x[t,s] <= demand[t,s])
        m.Equation(x[t,s] <= theta[s])

2. 目标值为负的原因与修复

目标值为负大概率是成本项的计算逻辑错误：

• 当前代码中每个时间步都重复扣除CAPEX和FOX，相当于把初始投资和年固定成本重复计算了30次，直接放大了成本规模；
• C1_ROX、C1_UOX如果是年度固定成本，不应该和产量绑定重复扣除。

修复方法：
拆分目标函数的盈利、固定成本、可变成本项，明确成本的计算周期：

obj = 0
for s in range(N):
    # CAPEX是初始投资，仅在第0年扣除一次
    capex = C1_CAPEX + C2_CAPEX*theta[s]
    # 年固定运营成本（FOX、C1_ROX、C1_UOX）
    annual_fixed_cost = (C1_FOX + C2_FOX*theta[s]) + (C1_ROX + C1_UOX)
    for t in range(T):
        discount = (1/(1+r))**(t+1)
        # 销售收入
        revenue = P_CO * x[t,s]
        # 碳收益（确认减排量计算逻辑正确）
        carbon_benefit = beta_CO2 * P_CO2 * x[t,s] * (E_ref - E_dir - E_indir_others - E_indir_elec_cons*GWP_EU[t,s])
        # 可变运营成本
        var_cost = (C2_ROX + C2_UOX) * x[t,s]
        
        # 单时间步净收益
        net = revenue + carbon_benefit - var_cost - annual_fixed_cost
        # 第0年额外扣除初始投资
        if t == 0:
            obj += discount * (net - capex)
        else:
            obj += discount * net
m.Maximize(obj/N)

3. 额外建议

• 先测试小规模模型：比如设置T=1、N=1，验证目标值是否为正，再逐步扩展到全规模；
• 打印中间变量：比如输出单个场景下的盈利、成本项，核对每部分的符号和数值是否符合预期。

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