GEKKO动态优化问题调整时间点、操纵变量或配置节点后不收敛的原因

使用Gekko求解动态优化问题时，调整时间点数量、操纵变量（manipulated variable）的初始值/上下限，或更改配置节点（collocation nodes）数后，部分参数组合下求解完全不收敛，出现如下错误：

raise Exception(response)

Exception:  @error: Solution Not Found

具体代码如下：

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt  

m = GEKKO()     
nt = 101        # no. of time steps
m.time = np.linspace(0,147,nt)

k17 = 0.0508
k26 = 1
k27 = 0.0577
k28 = 0.0104
k32 = 2
k33 = 2
k37 = 0.0016
k38 = 0.0107
k47 = 0.006
k48 = 0.072
k53 = 2
k57 = 0.0201
k58 = 0.082
k61 = 2
k62 = 2
k63 = 2
k72 = 2
k77 = 0.0133
k78 = 0.011
k87 = 0.081
k88 = 0.0148
kd = 0.06
w = 0.043

s1 = m.Var(value=0,lb=0) 
s2 = m.Var(value=3,lb=0)
s3 = m.Var(value=0.01,lb=0)
s4 = m.Var(value=0.46,lb=0)
s5 = m.Var(value=0.27,lb=0)
p1 = m.Var(value=0.51,lb=0)
p2 = m.Var(value=0.33,lb=0)
p3 = m.Var(value=0.3,lb=0)
p4 = m.Var(value=0.34e-4,lb=0)
p5 = m.Var(value=2.01,lb=0)
p6 = m.Var(value=0.05,lb=0)
X = m.Var(value=0.09,lb=0)
Xd = m.Var(value=0.02,lb=0)
V = m.Var(value=5,lb=0,ub=10)         
u = m.MV(value=0.05,lb=0,ub=0.1) # manipulated variable
u.STATUS = 1
u.DCOST = 0

f1 = (8.443e-4)*X*s1/(8.989e5 + s1)
f2 = (2.481e6)*X*s1*s3/((6.495e4 + s1)*(7.076e2 + s3))
f3 = (3.968e5)*X*s1*s3/((3.723e4 + s1)*(2.782e3 + s3))
f4 = (1.09e2)*X*s3/(0.019+s3)
f5 = 7.283*X*s4/(1.92e3 + s4)
f6 = (3.337e5)*X*s2*s5/((2.719e4 + s2)*(4.488e4 + s5))
f7 = (3.977e3)*X*s2/(9.324e3 + s2)
f8 = (6.697e-6)*X*s2/(0.537+s2)
f9 = (3.261e4)*X*s2/(6.683e5 + s2)

P = np.zeros(nt)
P[-1] = 1.0
final = m.Param(value=P)

# Equations
m.Equation(V.dt()==u)     
m.Equation(s1.dt()==-f1-f2-f3-k17*f7+(2-s1)*u/V)
m.Equation(s2.dt()==-f6-k27*f7-k28*f8-f9-s2*u/V)
m.Equation(s3.dt()==-k32*f2-k33*f3-f4+f6-k37*f7-k38*f8+f9-s3*u/V)
m.Equation(s4.dt()==-f5+f6-k47*f7-k48*f8-s4*u/V)
m.Equation(s5.dt()==k53*f3+f5-f6-k57*f7-k58*f8-s5*u/V)
m.Equation(p1.dt()==k61*f1+k62*f2+k63*f3-p1*u/V)
m.Equation(p2.dt()==k72*f2-k77*f7-k78*f8-p2*u/V)
m.Equation(p3.dt()==f4-k87*f7-k88*f8-p3*u/V)
m.Equation(p4.dt()==f8-p4*u/V)
m.Equation(p5.dt()==f7-p5*u/V)
m.Equation(p6.dt()==f5+f9-p6*u/V)
m.Equation(X.dt()==w*X-kd*X*Xd-X*u/V)
m.Equation(Xd.dt()==kd*X*Xd-Xd*u/V)

m.Obj(-final*V*p4) 
m.options.IMODE = 6
m.options.NODES = 4
#m.options.COLDSTART=2
#m.options.MAX_ITER=1000
m.solve(disp=True) 

p4_ = np.multiply(p4.value,1000)

plt.figure(1)
plt.subplot(2,1,1)
plt.plot(m.time,u.value,'r-')
plt.subplot(2,1,2)
plt.plot(m.time,p4_,'b--')

收敛问题的排查与解决方法 #

1. 冷启动与初始猜测优化 #

• 启用冷启动：取消注释m.options.COLDSTART=2，先求解稳态问题获得可行初始点，再启动动态优化，避免初始猜测不合理导致的求解失败。
• 预模拟获取初始值：先将IMODE设为4（动态模拟），固定操纵变量u为初始值运行一次，将得到的状态变量值作为动态优化的初始猜测，提升求解器起点合理性。

2. 数值稳定性优化 #

• 变量缩放：你的变量量级差异极大（如p4为1e-4量级，s1分母达1e5量级），易引发数值病态。可对小量级变量缩放（如p4乘以1000后定义变量），或在方程中统一量级，降低数值误差。
• 避免除零：V的下限设为0，当V趋近于0时u/V会出现数值爆炸。将V的下限改为lb=1e-3，规避除零风险。

3. 求解器参数调整 #

• 增加迭代次数：取消注释m.options.MAX_ITER=1000，甚至调至2000，给求解器足够的迭代空间。
• 调整配置节点数：NODES=4会增加问题复杂度，先尝试NODES=2或3，得到可行解后再逐步提高节点数提升精度。
• 切换求解器：默认APOPT求解器不收敛时，切换到IPOPT（m.options.SOLVER=3），不同求解器对非线性问题的处理能力不同，IPOPT在部分场景下收敛性更优。

4. 操纵变量约束与惩罚调整 #

• 添加控制平滑惩罚：当前u.DCOST=0，求解器可能生成剧烈控制动作导致数值不稳定。设置u.DCOST=1e-4，平滑控制曲线，帮助收敛。
• 暂时放宽约束：若u的0-0.1范围过紧，先放宽到0-0.2，确认收敛后再逐步缩小约束范围。

5. 问题规模优化 #

• 减少时间点数量nt（如先设为51），降低问题规模，得到可行解后再逐步增加时间点提高分辨率。

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