提取协态与约束乘子¶

pockit.optimizer.ipopt.solve 返回的第二个结果是 CyIpopt 的原始信息字典。其中包含 NLP 乘子，但这些数值还不是连续时间协态或路径约束乘子密度。转换公式取决于 Pockit 实际使用的缺陷方程、求积权重和约束排列。

本指南针对当前 Lobatto 和 Radau 离散严格推导该转换，并说明为什么不能把针对微分矩阵配点格式的公式直接套用到 Pockit。

Ipopt 的符号约定¶

设决策变量为 \(z\)，约束函数为 \(c(z)\)。Ipopt 使用以下驻值条件约定：

\[ \nabla F(z)+J_c(z)^\mathsf{T}\nu-z_L+z_U=0. \]

返回数组的含义如下：

字段	含义
`info["mult_g"]`	`system.constraints(z)` 的有符号乘子 \(\nu\)
`info["mult_x_L"]`	非负的变量下界乘子 \(z_L\)
`info["mult_x_U"]`	非负的变量上界乘子 \(z_U\)

对于按 \(c_L\le c(z)\le c_U\) 存储的约束，激活上界通常对应 \(\nu>0\)，激活下界通常对应 \(\nu<0\)；等式乘子没有预定符号。若把 \(c\) 改写为 \(-c\)，乘子也会反号。因此，保存对偶结果时必须同时记录原始约束表达式。

Pockit 的离散形式¶

定义

\[ \tau=\frac{t-t_0}{\Delta t},\qquad \Delta t=t_f-t_0, \]

并令 \(w_j=\texttt{phase.w_m[j]}\) 为归一化阶段区间 \([0,1]\) 上的求积权重。Pockit 将 Bolza 型目标离散为

\[ F_N=\Phi+\Delta t\sum_j w_j L_j. \]

对于第 \(i\) 个状态分量，Pockit 的积分缺陷方程为

\[ d_{i,r}=(T x_i)_r-\Delta t\sum_j I_{rj}f_{i,j}=0, \]

其中 \(I=\texttt{phase.I_m}\)。两种离散都从每个子区间的右端向左积分。正是 T*x - dt*I*f 这一方向决定了后续公式中的负号。

设 \(\alpha_{i,r}\) 是这些动态缺陷在 mult_g 中的乘子， \(\beta_{a,j}\) 是一般路径表达式 \(g_a(x_j,u_j,t_j)\) 的乘子。将离散 NLP 拉格朗日函数与

\[ H=L+\lambda^\mathsf{T}f+\mu^\mathsf{T}g \]

的求积近似逐项配平，可得

\[ \boxed{\lambda_{i,j}=-\frac{(I^\mathsf{T}\alpha_i)_j}{w_j}}, \qquad \boxed{\mu_{a,j}=\frac{\beta_{a,j}}{\Delta t\,w_j}}. \]

协态公式中没有 \(\Delta t\)：NLP 动态缺陷和连续 Hamiltonian 中的动态项都带有同一因子，它们相互抵消。一般路径约束在 system.constraints 中没有求积权重，因此其原始 NLP 乘子必须除以 \(\Delta t w_j\) 才是连续乘子密度。

应直接使用 phase.w_m。对于多区间 Lobatto 阶段，Pockit 已经在每个共享网格节点上累加了左右子区间的求积贡献。手工取某一侧参考单元的权重会使接口处的尺度错误。

约束排列¶

info["mult_g"] 与 system.constraints 的返回顺序完全一致：

所有有限维系统约束，长度为 system.n_c；
随后按照 system.p 的顺序逐阶段排列；
每个阶段先放动态缺陷，按状态分组，并由 phase.l_d / phase.r_d 索引；
再放一般路径约束，按表达式分组，每个表达式有 phase.L_m 个值。

第三部分的总长度为 phase.r_d[-1]，第四部分的长度为 phase.n_c * phase.L_m。这里 phase.n_c 只统计一般表达式。若 set_phase_constraint 中传入的是单独的状态、控制、时间或静态参数符号， Pockit 会将它转换为变量界；这种约束不会占用 mult_g 的行。

最前面的系统约束乘子是有限维对偶变量，应按原值返回，不能除以求积权重或阶段时长。

下面的辅助函数不硬编码 Lobatto 或 Radau 的矩阵尺寸，可以解码所有阶段：

import numpy as np


def extract_continuous_multipliers(system, solution, info):
    """提取系统乘子、缺陷乘子、协态和路径乘子密度。"""
    mult_g = np.asarray(info["mult_g"], dtype=float)
    phase_values = solution if isinstance(solution, (list, tuple)) else [solution]
    phase_values = phase_values[: system.n_p]
    if len(phase_values) != system.n_p:
        raise ValueError("solution does not contain every phase")

    result = {
        "system": mult_g[: system.n_c].copy(),
        "phases": [],
    }
    offset = system.n_c

    for phase, value in zip(system.p, phase_values):
        dynamic_size = int(phase.r_d[-1]) if phase.n_x else 0
        dynamic_raw = mult_g[offset : offset + dynamic_size]
        offset += dynamic_size

        alpha = [
            dynamic_raw[phase.l_d[i] : phase.r_d[i]].copy()
            for i in range(phase.n_x)
        ]
        defect_multiplier = (
            np.vstack(alpha) if alpha else np.empty((0, 0), dtype=float)
        )
        costate = (
            np.vstack(
                [
                    -np.asarray(phase.I_m.T @ alpha_i).ravel() / phase.w_m
                    for alpha_i in alpha
                ]
            )
            if alpha
            else np.empty((0, phase.L_m), dtype=float)
        )

        path_size = phase.n_c * phase.L_m
        path_nlp = mult_g[offset : offset + path_size].reshape(
            phase.n_c, phase.L_m
        )
        offset += path_size

        dt = value.t_f - value.t_0
        if dt <= 0.0:
            raise ValueError("phase duration must be positive")
        path_density = path_nlp / (dt * phase.w_m[None, :])
        time = value.t_0 + dt * phase.t_m

        result["phases"].append(
            {
                "time": time,
                "defect_multiplier": defect_multiplier,
                "costate": costate,
                "path_nlp": path_nlp.copy(),
                "path_density": path_density,
            }
        )

    if offset != len(mult_g):
        raise RuntimeError("multiplier layout does not match the configured system")
    return result

solution、info 和 system.p 必须来自同一次最终求解。网格调整会改变所有相关切片，因此上一轮 NLP 的乘子不能与新网格混用。

Lobatto 与 Radau 端点¶

Lobatto 的 phase.t_m 包含阶段的两个端点，因此上面的协态公式会直接返回端点协态。在内部网格点，phase.w_m 和 phase.I_m.T @ alpha 都会合并相邻两个子区间的贡献。

Radau 子区间的右端是状态节点，但不是控制、求积或一般路径约束节点，所以它不在 phase.t_m 中。若 \(\alpha_i^{(k)}\) 是第 \(k\) 个子区间的动态缺陷乘子，则其右端协态为

\[ \boxed{\lambda_i(t_{k+1}^{-})=\sum_r\alpha_{i,r}^{(k)}}. \]

Radau 的动态缺陷行分区与 phase.l_m / phase.r_m 一致：

alpha_i = decoded["phases"][phase_index]["defect_multiplier"][state_index]
right_costate = np.array(
    [alpha_i[left:right].sum() for left, right in zip(phase.l_m, phase.r_m)]
)
right_time = phase_value.t_0 + (
    phase_value.t_f - phase_value.t_0
) * phase.t_x[phase.r_x - 1]

在内部接口处，应将该右极限与下一子区间左侧 Radau 节点上的映射协态比较。明显跳变可能来自离散精度不足、真实内部事件或遗漏的阶段连接条件。

phase.set_boundary_condition 中的固定量由 Pockit 直接代入，并不会在 mult_g 中新增约束行。因此，不能从该数组中切出一个独立的“固定边界乘子”。需要使用端点协态与横截条件；若有限维关系的乘子本身就是研究对象，则应在模型允许时将该关系写成系统约束。

直接变量界¶

若路径表达式是单独的控制符号，Pockit 会把它提升为 NLP 变量界。此时控制上下界对应的连续非负乘子密度为

\[ \mu_{L,j}=\frac{z_{L,j}}{\Delta t w_j},\qquad \mu_{U,j}=\frac{z_{U,j}}{\Delta t w_j}. \]

原始有界表达式的有符号乘子为 \(\widetilde\mu=\mu_U-\mu_L\)：激活上界时为正，激活下界时为负。

def extract_control_bound_density(
    system, phase_value, info, phase_index, control_index
):
    phase = system.p[phase_index]
    phase_offset = sum(system.p[i].L for i in range(phase_index))
    local_left = phase.l_v[phase.n_x + control_index]
    local_right = phase.r_v[phase.n_x + control_index]
    selection = slice(phase_offset + local_left, phase_offset + local_right)

    dt = phase_value.t_f - phase_value.t_0
    scale = dt * phase.w_m
    lower = np.asarray(info["mult_x_L"])[selection] / scale
    upper = np.asarray(info["mult_x_U"])[selection] / scale
    return lower, upper, upper - lower

该函数有意只处理控制，因为两种离散的控制优化节点都与 phase.w_m 对齐。 Radau 状态还包含一个没有求积权重的右端节点；该节点上的变量界乘子是端点原子，不能当作路径密度除以 w_m。

可复现的协态验证¶

提出本指南问题的两个 notebook 研究了

\[ \min -y(2),\qquad \dot y=\frac{5}{2}(-y+yu-u^2),\qquad y(0)=1. \]

对于 \(H=\lambda f\)，其最优解和协态为

\[ y^*(t)=\frac{4}{1+3e^{5t/2}},\qquad u^*(t)=\frac{y^*(t)}{2}, \]

\[ \lambda^*(t)= -\frac{(1+3e^{5t/2})^2e^{-5t/2}}{e^{-5}+6+9e^5}. \]

以下脚本用 Lobatto 和 Radau 求解同一个多区间问题：

import numpy as np

from pockit.lobatto import System as LobattoSystem
from pockit.lobatto import linear_guess as lobatto_guess
from pockit.optimizer import ipopt
from pockit.radau import System as RadauSystem
from pockit.radau import linear_guess as radau_guess


def exact_y(t):
    return 4.0 / (1.0 + 3.0 * np.exp(2.5 * t))


def exact_u(t):
    return exact_y(t) / 2.0


def exact_costate(t):
    a = 1.0 + 3.0 * np.exp(2.5 * t)
    b = np.exp(-5.0) + 6.0 + 9.0 * np.exp(5.0)
    return -(a**2) * np.exp(-2.5 * t) / b


def run(System, make_guess):
    system = System(["y_final"])
    (y_final,) = system.s
    phase = system.new_phase(["y"], ["u"])
    (y,) = phase.x
    (u,) = phase.u

    phase.set_dynamics([2.5 * (-y + y * u - u**2)])
    phase.set_boundary_condition([1.0], [y_final], 0.0, 2.0)
    phase.set_discretization([0.0, 0.2, 0.55, 1.0], [8, 9, 10])
    system.set_phase([phase])
    system.set_objective(-y_final)

    guess = make_guess(phase, 0.0)
    guess.x[0] = exact_y(guess.t_x)
    guess.u[0] = exact_u(guess.t_u)
    solution, info = ipopt.solve(
        system,
        [guess, [exact_y(2.0)]],
        optimizer_options={
            "print_level": 0,
            "sb": "yes",
            "tol": 1.0e-11,
            "acceptable_tol": 1.0e-11,
            "bound_relax_factor": 0.0,
        },
    )
    if int(info["status"]) not in (0, 1):
        raise RuntimeError(info["status_msg"])

    decoded = extract_continuous_multipliers(system, solution, info)
    dual = decoded["phases"][0]
    error = np.max(np.abs(dual["costate"][0] - exact_costate(dual["time"])))
    assert error < 1.0e-6
    return error


print("Lobatto:", run(LobattoSystem, lobatto_guess))
print("Radau:  ", run(RadauSystem, radau_guess))

在上述网格上，Lobatto 和 Radau 的最大误差分别为 \(2.8\times10^{-8}\) 与 \(3.5\times10^{-9}\)。若反转框内协态公式的符号，误差约为 \(2\)，这直接验证了缺陷方向。Radau 区间右端公式与解析端点协态的误差为 \(1.4\times10^{-15}\)，最大接口左右差为 \(3.5\times10^{-9}\)。

可复现的路径乘子验证¶

为独立检查 \(\Delta t w_j\) 缩放，考虑

\[ \min\int_0^{3.7}\frac{1}{2}(u-2)^2\,\mathrm dt \]

且满足 \(g(u)=u+0.1u^3\le g(0.5)\)。最优解是 \(u^*=0.5\)，驻值条件给出的常数上界乘子密度为

\[ \mu^*=\frac{2-0.5}{1+0.3(0.5)^2}=1.395348837209\ldots. \]

from pockit.lobatto import System, constant_guess

system = System(0)
phase = system.new_phase(["x"], ["u"])
(x,) = phase.x
(u,) = phase.u
g = u + 0.1 * u**3
g_upper = 0.5 + 0.1 * 0.5**3

phase.set_dynamics([0.0 * x])
phase.set_integral([(u - 2.0) ** 2 / 2.0])
phase.set_phase_constraint([g], [-np.inf], [g_upper])
phase.set_boundary_condition([0.0], [0.0], 0.0, 3.7)
phase.set_discretization([0.0, 0.2, 0.55, 1.0], [7, 8, 9])
system.set_phase([phase])
system.set_objective(phase.I[0])

guess = constant_guess(phase, 0.0)
guess.u[0] = 0.5
solution, info = ipopt.solve(
    system,
    guess,
    optimizer_options={
        "print_level": 0,
        "sb": "yes",
        "tol": 1.0e-11,
        "acceptable_tol": 1.0e-11,
        "bound_relax_factor": 0.0,
    },
)
if int(info["status"]) not in (0, 1):
    raise RuntimeError(info["status_msg"])
decoded = extract_continuous_multipliers(system, solution, info)
mu = decoded["phases"][0]["path_density"][0]
expected = 1.5 / 1.075
np.testing.assert_allclose(mu, expected, rtol=0.0, atol=1.0e-8)

将其改为下界后，有符号 mult_g 密度会反号。对 Lobatto/Radau、上下界、一般约束和提升后的变量界进行组合测试时，恢复值与解析密度的误差均不超过 \(5.3\times10^{-11}\)。

可靠性检查¶

恢复出的乘子是数值结果，而不是无需检查的精确证书：

读取对偶变量前先确认 Ipopt 状态；
检查原始可行性和 KKT 驻值残差；
提高多项式次数或细化网格后比较乘子历史；
未激活不等式的乘子只会在求解器容差内接近零；
若存在冗余约束或激活约束梯度线性相关，不应单独解释某个非唯一乘子；
在切换、碰撞等真实不连续点保留左右极限；
若手工缩放了目标或约束，必须显式还原该缩放。

原始 NLP 乘子通常会随网格变化。应当收敛到与网格无关连续极限的是经过上述映射得到的协态和乘子密度。