Julia中IPOPT求解最优控制问题的离散化N值异常排查

安全阿岳

2026-6-30

最优控制求解异常与动力学验证不匹配的排查方案

在Julia中基于IPOPT求解器实现最优控制问题：

采用梯形法则近似积分型目标函数，代码表达式：

@expression(sys, J1, 0.5 * Δt * sum((A * x[t] + B * u[t] + C * v[t] + A * i[t+1] + B * u[t+1] + C * v[t+1]) for t in 1:N))

遇到的异常现象：

核心变量校验：表达式中的i[t+1]疑似笔误，应为状态量x[t+1]？若误用非状态变量，会直接导致目标函数计算逻辑错误。
被积项形式校验：如果目标是常见的二次型性能指标（如LQR类问题），当前表达式缺少平方/内积运算，正确形式应为(A*x[t]+B*u[t]+C*v[t])'*(A*x[t]+B*u[t]+C*v[t])与(A*x[t+1]+B*u[t+1]+C*v[t+1])'*(A*x[t+1]+B*u[t+1]+C*v[t+1])之和，漏掉平方会完全扭曲优化目标。
梯形法则系数校验：标准梯形法系数为$\frac{\Delta t}{2}$，当前代码的0.5 * Δt是正确的，但需确认求和范围是否覆盖所有离散区间（$t=1$到$N$对应$N$个区间，终端状态$x_{N+1}$是否被正确包含）。

状态转移约束校验：严格遵循Crank-Nicolson标准形式：
$$x_{t+1} = x_t + \frac{\Delta t}{2}\left(f(x_t,u_t,v_t) + f(x_{t+1},u_{t+1},v_{t+1})\right)$$
检查代码中约束是否准确实现上述等式，例如：
```
@constraint(sys, x[t+1] == x[t] + 0.5*Δt*(f(x[t],u[t],v[t]) + f(x[t+1],u[t+1],v[t+1])))
```
边界条件校验：确认初始状态$x_1$是否正确固定，终端状态$x_{N+1}$是否符合问题要求（自由/固定/约束），边界条件错误会导致整个离散系统解偏离真实最优轨迹。
矩阵维度校验：检查$A,B,C$的维度是否与状态$x$、控制$u$、扰动$v$完全匹配，Julia的自动广播可能掩盖维度不匹配问题，导致计算结果不符合物理意义。

收敛容调优：将IPOPT的收敛容差tol从默认的$1e-6$调至$1e-8$，增大max_iter，避免求解器在大$N$场景下提前终止于局部最优解。
Hessian近似设置：当$N$较大时，精确Hessian计算可能出现数值病态，尝试开启有限记忆Hessian近似：
```
solver = Ipopt.Optimizer()
set_optimizer_attribute(solver, "hessian_approximation", "limited-memory")
```
初始值暖启动：用小$N$的最优解作为大$N$场景的初始值（warm start），避免求解器收敛到不同的局部最优解，验证解的突变是否由初始值导致。

动力学一致性校验：确保验证用的RK方法与原问题的连续动力学$\dot{x}=f(x,u,v)$完全一致，RK阶数、步长需与离散化的$\Delta t$匹配。
输入完整性校验：验证时严格使用求解器输出的所有$u[t],v[t]$（$t=1$到$N$），从初始状态$x_1$开始逐段积分，避免遗漏任意时刻的控制/扰动导致轨迹偏差。
约束残差校验：手动计算求解器输出$x$的Crank-Nicolson约束残差：
$$res_t = x_{t+1} - x_t - \frac{\Delta t}{2}\left(f(x_t,u_t,v_t) + f(x_{t+1},u_{t+1},v_{t+1})\right)$$
若残差远大于IPOPT的收敛容差，说明约束未被正确满足，是求解器设置或约束实现的核心问题。

Δt变化影响：计算不同$N$对应的$\Delta t=T/N$（$T$为总时间），检查$N=7000$到$7200$时$\Delta t$的变化是否触发数值方法的稳定性切换（尽管Crank-Nicolson是A-稳定，但刚性系统仍可能受步长影响）。
中间N值测试：测试$N=7100、7150$等中间值，观察解的突变是渐变还是阶跃，判断是数值阈值问题还是求解器收敛到不同局部最优。

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