最优控制问题求解求助：Python代码无解决方案求排查建议

问题背景 #

目标为在最终时刻最小化$M_u$和$F_u$，最大化$M_i$和$F_i$，对应目标函数：
$$J(u)=M_u^2+F_u2+ \int_0^T(M_i+F_i)u2dt$$
需控制的ODE系统如下：
$$
\begin{align}
\alpha &= \left( 1-\frac{A_u+A_i}{K_a}\right)\
\delta &= (\gamma + \mu_a )A_u \
\sigma &= M_uF_u+(1-\rho_1) M_iF_i +(1-\rho_2)M_uF_i+(1-\rho_3\lambda_m) M_iF_u \
\frac{dA_u}{dt}&= \phi\left( \frac{\sigma}{M_u+M_i}\right)\alpha-\delta \
\frac{dA_i}{dt}&= \phi\left( \frac{\rho_1M_iF_i+\rho_2M_uF_i+ \rho_3\lambda_m M_iF_u}{M_u+ M_i}\right) \alpha - \delta \
\frac{dF_u}{dt}&= b_f\gamma A_u -\beta_m \frac{M_i}{M_u+ M_i} F_u- \mu_f F_u \
\frac{dF_i}{dt}&= b_f\gamma A_i +\beta_m \frac{M_i}{M_u+ M_i} F_u- \mu_f F_i+uF_i \
\frac{dM_u}{dt}&= b_m\gamma A_u -\beta_f \frac{F_i}{F_u+F_i} M_u- \mu_m M_u \
\frac{dM_i}{dt}&= b_m\gamma A_i +\beta_f \frac{F_i}{F_u+F_i} M_u- \mu_m M_i +uM_i.
\end{align}
$$

用户提供的Python求解代码：

import numpy as np
from gekko import GEKKO
m = GEKKO(remote=True)

m.time = np.concatenate((np.array([0.0,0.01,0.02,0.05,0.1,0.2,0.5]),
                        np.linspace(1,500,500)))
nt = len(m.time)

Au = m.Var(10000)
Ai = m.Var(0)
Fu = m.Var(15000)
Fi = m.Var(200)
Mu = m.Var(15000)
Mi = m.Var(200)
M = m.Var(1,lb=1)
F = m.Var(1,lb=1)
u = m.MV(value=0.06,lb=0.02, ub=0.06)
u.STATUS=1
p = np.zeros(nt); p[-1] = 1.0; final = m.Param(value=p)


Ka = 20000
phi =5
bf = 0.5
bm = 0.5
ro1= 0.76
ro2 = 0.45
ro3 = 0.37
lm =0.51
gamma =0.2
betam = 0.035
betaf = 0.032
ua = 0.025

um = 0.12
uf = 0.1

# system
m.Equation(M==Mu+Mi)
m.Equation(F==Fu+Fi)
m.Equation(Au.dt() == phi*(Mu*Fu+(1-ro1)*Mi*Fi+(1-ro2)*Mu*Fi+(1-ro3*lm)*Mi*Fu)*(1-(Au+Ai)/Ka)/M-(gamma+ua)*Au)
m.Equation(Ai.dt() == phi*(ro1*Mi*Fi+ro2*Mu*Fi+ro3*lm*Mi*Fu)*(1-(Au+Ai)/Ka)/M-(gamma+ua)*Ai)
m.Equation(Fu.dt() == bf*gamma*Au-betam*Mi*Fu/M-uf*Fu)
m.Equation(Fi.dt() == bf*gamma*Ai+betam*Mi*Fu/M-uf*Fi+u*Fi)
m.Equation(Mu.dt() == bm*gamma*Au-betaf*Fi*Mu/F-um*Mu)
m.Equation(Mi.dt() == bm*gamma*Ai+betaf*Fi*Mu/F-um*Mi+u*Mi)
#m.Equation(final*(y-0.0)==0)
m.Minimize(final*((Fu-0.0)**2+(Mu-0.0)**2))

# cost function
m.Minimize(m.integral((Mi+Fi)*u**2)*final)
           
m.options.IMODE = 6
m.options.MAX_ITER = 2000
m.options.NODES = 2
m.options.MV_TYPE = 1
m.options.SOLVER = 2

m.options.COLDSTART = 1
m.solve()

m.options.COLDSTART = 0
m.options.TIME_SHIFT = 0
m.solve()

print('Objective: ' + str(m.options.OBJFCNVAL))

import matplotlib.pyplot as plt
plt.plot(m.time,Au.value+Fu.value+Mu.value,lw=2,label=r'$y$')
plt.plot(m.time,Ai.value+Fi.value+Mi.value,lw=2,label=r'$x$')
plt.grid()
plt.plot(m.time,u.value,'g',lw=2,label=r'$u$')
plt.grid()
plt.xlabel('Time')

关键修正与求解建议 #

1. 修复目标函数的核心错误 #

• 补充最大化目标的实现：原目标要求最大化$M_i$和$F_i$，需转化为求解器支持的最小化形式，在最终时刻代价项中加入加权负项，示例：

  # 权重可根据变量数量级调整，平衡最小化与最大化目标的优先级
  weight = 1e-3
  m.Minimize(final*((Fu)**2 + (Mu)**2 - weight*(Mi + Fi)))
  

• 修正积分项的时间范围：代码中多余的*final导致积分仅计算最后时刻的值，需移除以覆盖整个时间区间：

  m.Minimize(m.integral((Mi+Fi)*u**2))
  

  注：若积分项目的是辅助最大化$M_i+F_i$，需将其改为负号（最小化负积分等价于最大化原积分），需根据实际需求确认目标函数逻辑。

2. 简化变量定义，减少数值冗余 #

移除手动定义的M和F变量，改用Intermediate变量直接计算，避免额外等式约束带来的数值负担：

# 替换原M、F的Var定义
M = m.Intermediate(Mu + Mi)
F = m.Intermediate(Fu + Fi)

3. 优化求解器参数设置 #

• 提升动态优化精度，将NODES调整为3或5：

  m.options.NODES = 3
  

• 开启自动变量缩放，处理不同变量间的数量级差异：

  m.options.SCALING = 1
  

• 若收敛困难，可先放宽控制量u的上下界，或调整初始值，待收敛后再逐步收紧约束。

4. 增强调试与可视化 #

• 启用求解器调试日志，排查收敛失败的具体原因：

  m.solve(debug=True)
  

• 增加关键变量的单独可视化，直观观察状态变化趋势：

  plt.figure(figsize=(12,8))
  plt.subplot(2,2,1)
  plt.plot(m.time, Mu.value, label='Mu')
  plt.plot(m.time, Fu.value, label='Fu')
  plt.title('Minimized Variables')
  plt.legend(); plt.grid()
  plt.subplot(2,2,2)
  plt.plot(m.time, Mi.value, label='Mi')
  plt.plot(m.time, Fi.value, label='Fi')
  plt.title('Maximized Variables')
  plt.legend(); plt.grid()
  plt.subplot(2,2,3)
  plt.plot(m.time, Au.value, label='Au')
  plt.plot(m.time, Ai.value, label='Ai')
  plt.title('A States')
  plt.legend(); plt.grid()
  plt.subplot(2,2,4)
  plt.plot(m.time, u.value, label='u')
  plt.title('Control Input')
  plt.legend(); plt.grid()
  plt.tight_layout()
  plt.show()
  

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