超敏感问题

背景

超敏感最优控制问题同时具有很长的时域和端点附近的快速变化。大部分轨迹接近缓慢变化的内部解，但狭窄的初始边界层和终端边界层决定了两个端点条件能否同时满足。若采用均匀细网格，大多数节点会被浪费在变化很小的区域。

问题定义

标准标量基准问题为

\[ \begin{aligned} \min_{x,u}\quad & \frac{1}{2}\int_0^{10000}\left(x^2+u^2\right)\,\mathrm{d}t \\ \text{约束}\quad & \dot{x}=-x^3+u, \\ & x(0)=1.5, \qquad x(10000)=1. \end{aligned} \]

二次型目标在整个时域上同时惩罚状态幅值和控制消耗。

变量与单位

该基准问题完全采用缩放后的无量纲变量，模型时域长度为 10,000。若未先定义参考尺度而直接把时间写成秒，归一化动力学在量纲上并不成立。

符号	含义	单位
\(t\)	模型时间	-
\(x\)	标量状态	-
\(u\)	标量控制	-

建模选择

实现采用单个固定时间的 Lobatto 阶段。控制初始猜测近似满足缓慢变化状态猜测对应的动力学。每次成功求解后，网格误差检查会判断是否需要继续更新；网格调整将分辨率集中在两个边界层附近，而不会在平缓的内部区段均匀堆积节点。

运行示例

在仓库根目录运行完整精度示例：

python -m examples.hyper_sensitive

这是一个有意设置得较困难的基准，运行时间可能超过一分钟。如需保存图片且不打开窗口：

python -m examples.hyper_sensitive --save hyper-sensitive.png --no-show

如需进行耗时更短的冒烟测试，可添加 --quick。快速模式使用较粗的误差容差，并最多更新三次网格；它用于演示流程，不是下文报告的验证结果。

关键实现

import numpy as np

from pockit.lobatto import System, linear_guess
from pockit.optimizer import ipopt

system = System(0)
phase = system.new_phase(["state"], ["control"])
(state,) = phase.x
(control,) = phase.u

phase.set_dynamics([-(state**3) + control])
phase.set_integral([(state**2 + control**2) / 2.0])
phase.set_boundary_condition([1.5], [1.0], 0.0, 10_000.0)
phase.set_discretization(10, 10)
system.set_phase([phase])
system.set_objective(phase.I[0])

guess = linear_guess(phase, 0.0)
state_at_control_nodes = np.interp(guess.t_u, guess.t_x, guess.x[0])
state_slope = (1.0 - 1.5) / 10_000.0
guess.u[0] = state_at_control_nodes**3 + state_slope
solution, info = ipopt.solve(system, guess)
for _ in range(20):
    if system.check(solution, absolute_tolerance_continuous=1e-8,
                    relative_tolerance_continuous=1e-8):
        break
    solution = system.refine(
        solution,
        absolute_tolerance_continuous=1e-8,
        relative_tolerance_continuous=1e-8,
        num_point_min=10,
        num_point_max=20,
        mesh_length_min=1e-8,
    )
    solution, info = ipopt.solve(system, solution)

由动力学构造的控制猜测使 \(\dot{x}\) 近似等于线性状态猜测的斜率；随后通过显式网格限制获得可靠收敛。

已验证结果

默认示例经过 18 次网格更新后，目标函数值达到 \(1.330806904944\)。全时域图显示内部区段几乎保持不变；下方局部图则展示了初始和终端的快速过渡，否则这些变化很难看清。

源码

完整可运行示例见 examples/hyper_sensitive.py。