Costate time stamps: costate(sol)(t_i) approximates p(t_i + h/2), or p(t_{i+1}) with Gauss–Legendre

Aperta
  • Qualcosa non va?
#640 0 commenti 0 reazioni 1 assegnatario Vedi su GitHub

I maintainer di solito rispondono entro 1 giorno

@PierreMartinon ci sta già lavorando.

Dal 9/10/2026.

Valutazione

Questa issue non è ancora stata valutata.

Descrizione

Summary

The costate of a direct solve is built from the multipliers $\lambda_i$ of the step equations, one per step, and stamped at the left node $t_i$ (T_costate = T[1:end-1]). But $\lambda_i$ does not approximate $p(t_i)$:

  • with :euler, :euler_implicit, :midpoint and :trapeze, it approximates $p(t_i + h/2)$, the middle of the step;
  • with :gauss_legendre_2 and :gauss_legendre_3, it approximates $p(t_{i+1})$, the end of the step.

So costate(sol)(t) lags by half a step, or by a full step with Gauss–Legendre. This gives an $O(h)$ error, which spoils the order of the Gauss–Legendre schemes on the costate. Users notice it when they seed a shooting method with costate(sol)(t0): with gauss_legendre_2 and grid_size = 100, they get $p(h)$, not $p(t_0)$.

Measurements

Energy-minimal double integrator, $\dot q = v$, $\dot v = u$, $x(0) = (-1, 0)$, $x(1) = 0$, $\frac12\int u^2 \to \min$. Exact costate $p(t) = (12, 6 - 12t)$. The offset $\theta$ is defined by costate(sol)(t_i)[2] $= 6 - 12(t_i + \theta h)$:

Scheme N = 50, nodes 1, N/4, N/2, N N = 200
:euler 0.49, 0.495, 0.50, 0.51 0.498, 0.499, 0.50, 0.502
:euler_implicit 0.49, 0.495, 0.50, 0.51 0.498, 0.499, 0.50, 0.502
:midpoint 0.49, 0.495, 0.50, 0.51 0.498, 0.499, 0.50, 0.502
:trapeze 0.46, 0.48, 0.50, 0.54 0.49, 0.495, 0.50, 0.51
:gauss_legendre_2 1.0, 1.0, 1.0, 1.0 1.0, 1.0, 1.0, 1.0
:gauss_legendre_3 1.0, 1.0, 1.0, 1.0 1.0, 1.0, 1.0, 1.0

θ stays the same when the grid is 4 times finer: the error is a shift of the time stamps, not a discretisation error. For comparison, the multiplier of the initial condition gives $p(t_0)$ directly: dual(sol, ocp, :start) returns $(12.0, 6.0)$ with :midpoint and :gauss_legendre_2 (see #639 for the derivation).

using OptimalControl, NLPModelsIpopt
ocp = @def begin
    t ∈ [0, 1], time
    x = (q, v) ∈ R², state
    u ∈ R, control
    x(0) == [-1, 0], (start)
    x(1) == [0, 0]
    ẋ(t) == [v(t), u(t)]
    0.5∫(u(t)^2) → min
end
for scheme in (:midpoint, :gauss_legendre_2)
    sol = solve(ocp; scheme, grid_size=100, display=false)
    println(scheme, ": costate(sol)(0) = ", costate(sol)(0.0), ", dual(:start) = ", dual(sol, ocp, :start))
end
# midpoint: costate(sol)(0) = [12.0, 5.94], i.e. p(h/2); dual(:start) = [12.0, 6.0]
# gauss_legendre_2: costate(sol)(0) = [12.0, 5.88], i.e. p(h); dual(:start) = [12.0, 6.0]

The maths to do first

Where a multiplier lives in time follows from the discrete adjoint equation of each scheme. For the midpoint rule, with $x_{i+1} - x_i - h f\big(\tfrac{x_i + x_{i+1}}{2}, u_i\big) = 0$ and multiplier $\lambda_i$, stationarity in $x_i$ gives

$$
\frac{\lambda_i - \lambda_{i-1}}{h} = -\tfrac12\big(\partial_x f_{i}^\top \lambda_i + \partial_x f_{i-1}^\top \lambda_{i-1}\big) + \dots,
$$

a trapezoidal step of $\dot p = -\partial_x f^\top p$ centred at $t_i$, between $\lambda_{i-1}$ and $\lambda_i$. So $\lambda_i$ sits at $t_i + h/2$. For the Gauss–Legendre schemes (state at the nodes, stage equations, then $x_{i+1} = x_i + h \sum_j b_j k_j$), the multiplier of the last equation should be shown to approximate $p(t_{i+1})$, as measured.

I suggest writing this out for each scheme before changing anything, then stamping each costate value at the time it approximates: $t_i + h/2$ (or the stage times), and $t_{i+1}$ for Gauss–Legendre. Alternatively, interpolate it back onto the nodes. The step-to-step costate time grid is already separate from the state grid (T_costate), so only the stamps change.

Related: #405 and #639 (the scale and sign of the constraint multipliers).

Lingua principale
Julia
Stelle
12
Fork
5
Merge medio
6h 23m
PR unite (30g)
15

Preparare l'ambiente

Come iniziare

  1. Leggi tutta la issue e poi la guida ai contributi del progetto.
  2. Commenta sulla issue per dire che te ne occupi tu — evita che due persone facciano lo stesso lavoro.
  3. Fai un fork del repository e lavora su un branch.
  4. Apri una pull request che faccia riferimento al numero della issue.

Altre issue di control-toolbox/CTDirect.jl

Tutte le issue di control-toolbox/CTDirect.jl

Issue simili

Altre issue su Julia

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.