Extracting Costates and Constraint Multipliers

pockit.optimizer.ipopt.solve returns the raw dictionary produced by CyIpopt as its second result. The dictionary includes NLP multipliers, but these values are not yet continuous-time costates or path-constraint multiplier densities. The conversion depends on Pockit's actual defect equations, quadrature weights, and constraint order.

This guide derives that conversion for the current Lobatto and Radau transcriptions. It also explains why formulas written for a differentiation-matrix transcription cannot be copied directly into Pockit.

Ipopt sign convention

For decision variables \(z\), constraint functions \(c(z)\), and variable bounds, Ipopt uses the stationarity convention

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

The returned arrays have the following meanings:

Entry	Meaning
info["mult_g"]	Signed multipliers \(\nu\) for system.constraints(z)
info["mult_x_L"]	Nonnegative lower-variable-bound multipliers \(z_L\)
info["mult_x_U"]	Nonnegative upper-variable-bound multipliers \(z_U\)

For a constraint stored as \(c_L\le c(z)\le c_U\), an active upper bound normally has \(\nu>0\), while an active lower bound normally has \(\nu<0\). An equality multiplier has no prescribed sign. Replacing \(c\) by \(-c\) reverses its multiplier, so the modeled expression must always be recorded with the dual.

Pockit's transcription

Let

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

and let \(w_j=\texttt{phase.w_m[j]}\) be a quadrature weight on the normalized phase interval \([0,1]\). Pockit approximates a Bolza objective as

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

For state component \(i\), Pockit stores the integral defect equations as

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

where \(I=\texttt{phase.I_m}\). Both schemes integrate backward from each subinterval's right endpoint. This T*x - dt*I*f orientation determines the minus sign below.

Let \(\alpha_{i,r}\) be the corresponding entries of mult_g, and let \(\beta_{a,j}\) multiply a general path expression \(g_a(x_j,u_j,t_j)\). Matching the discrete NLP Lagrangian to the quadrature approximation of

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

gives

\[ \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}}. \]

There is no \(\Delta t\) in the costate formula: the same factor multiplies the dynamics term in both the NLP defect and the continuous Hamiltonian. A general path constraint is not quadrature-weighted in system.constraints, so its raw NLP multiplier must be divided by \(\Delta t w_j\).

Use phase.w_m directly. For a multi-interval Lobatto phase, Pockit has already added the left- and right-interval quadrature contributions at every shared mesh node. Reconstructing a reference-element weight by hand gives the wrong scale at interfaces.

Constraint order

info["mult_g"] follows the exact order returned by system.constraints:

1. all finite-dimensional system constraints, with length system.n_c;
2. for each phase in system.p order:
3. dynamic defects, grouped by state and indexed by phase.l_d / phase.r_d;
4. general path constraints, grouped by expression, with phase.L_m values per expression.

The third item has total length phase.r_d[-1]. The fourth has length phase.n_c * phase.L_m. Here, phase.n_c counts only general expressions. Pockit recognizes a bare state, control, time, or static-parameter symbol in set_phase_constraint and converts it to a variable bound instead; such a bound does not occupy any mult_g rows.

The leading system-constraint multipliers are finite-dimensional duals. Return them as stored; do not divide them by a quadrature weight or phase duration.

The following helper decodes all phases without hard-coding Lobatto or Radau matrix dimensions:

import numpy as np


def extract_continuous_multipliers(system, solution, info):
    """Decode system duals, defect multipliers, costates, and path densities."""
    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

Always pass solution, info, and system.p from the same, final solve. Mesh refinement changes every affected slice and invalidates multipliers from the previous NLP.

Lobatto and Radau endpoints

For Lobatto, phase.t_m contains both phase endpoints. The boxed costate formula therefore returns endpoint costates directly. At an internal mesh point, phase.w_m and phase.I_m.T @ alpha both combine the adjacent subintervals.

For Radau, the right endpoint of a subinterval is a state node but not a control, quadrature, or general-path node. It is absent from phase.t_m. If \(\alpha_i^{(k)}\) denotes the defect multipliers of subinterval \(k\), its right endpoint costate is

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

In code, Radau's defect-row partition matches 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]

At an internal interface, compare this right-limit value with the mapped costate at the next interval's left Radau node. A significant jump can indicate an under-resolved solution, an actual interior event, or a missing linkage condition.

Fixed entries supplied to phase.set_boundary_condition are substituted by Pockit; they are not extra rows in mult_g. Consequently, there is no separate \"fixed boundary multiplier\" to slice from that array. Use the endpoint costate and transversality conditions, or model the required finite-dimensional relation as a system constraint when its multiplier is itself needed.

Direct variable bounds

When the path expression is a bare control symbol, Pockit promotes it to NLP variable bounds. For a control, the continuous nonnegative lower and upper multiplier densities are

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

The signed multiplier for the raw bounded expression is \(\widetilde\mu=\mu_U-\mu_L\). This is positive on an active upper bound and negative on an active lower bound.

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

This helper is intentionally limited to controls, whose optimization nodes match phase.w_m in both schemes. A Radau state has an additional right-endpoint node with no quadrature weight; a bound multiplier at that node is an endpoint atom, not a path-density value to divide by w_m.

Reproducible costate check

The notebooks that motivated this guide study

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

Its solution and costate for \(H=\lambda f\) are

\[ 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}. \]

The following script checks the same multi-interval problem with both schemes:

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))

With the displayed mesh, the maximum errors were \(2.8\times10^{-8}\) for Lobatto and \(3.5\times10^{-9}\) for Radau. Reversing the boxed costate sign gave an error of about \(2\), which directly checks the defect orientation. Radau's interval-right formula agreed with the analytical endpoint costate to \(1.4\times10^{-15}\); the largest left/right interface mismatch was \(3.5\times10^{-9}\).

Reproducible path-multiplier check

To check the \(\Delta t w_j\) scaling independently, minimize

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

subject to \(g(u)=u+0.1u^3\le g(0.5)\). The optimum is \(u^*=0.5\), and stationarity gives the constant upper-bound density

\[ \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)

Repeating this test with a lower bound reverses the signed mult_g density. Tests with Lobatto and Radau, upper and lower bounds, and both general constraints and promoted variable bounds all recovered the analytical density within \(5.3\times10^{-11}\).

Reliability checks

Treat a reconstructed multiplier as a numerical result, not as an exact certificate:

• verify Ipopt's status before reading duals;
• check primal feasibility and the KKT stationarity residual;
• compare multiplier histories after increasing the polynomial degree or refining the mesh;
• expect inactive inequality multipliers to be zero only to solver tolerance;
• avoid interpreting individual duals when constraints are redundant or the active-set gradients are linearly dependent;
• retain separate left and right limits near switches, impacts, and other genuine discontinuities;
• account explicitly for any manual objective or constraint rescaling.

The raw NLP multipliers generally change as the mesh changes. The reconstructed costates and multiplier densities are the quantities that should converge to a mesh-independent continuous-time limit.