Costate Estimation Using the Legendre–Gauss–Radau Method
This section develops costate estimation for an optimal-control problem transcribed with the Legendre–Gauss–Radau (LGR) pseudospectral method. The main objective is to derive a direct relationship between:
the Karush–Kuhn–Tucker (KKT) multipliers of the discretized nonlinear program; and
approximations of the continuous-time costate.
The LGR method is especially convenient for costate estimation because one endpoint is itself a collocation point. For the left-Radau grid considered here, the initial point is included among the LGR points, whereas the final point is a noncollocated interpolation point.
The final mapping is
λ 1 : N = W − 1 Λ 1 : N , \boxed{
\boldsymbol{\lambda}_{1:N}
=
\boldsymbol{W}^{-1}\boldsymbol{\Lambda}_{1:N},
} λ 1 : N = W − 1 Λ 1 : N , and
λ N + 1 = d N + 1 T Λ 1 : N , \boxed{
\boldsymbol{\lambda}_{N+1}
=
\boldsymbol{d}_{N+1}^{\mathsf{T}}\boldsymbol{\Lambda}_{1:N},
} λ N + 1 = d N + 1 T Λ 1 : N , where Λ 1 : N \boldsymbol{\Lambda}_{1:N} Λ 1 : N W \boldsymbol{W} W d N + 1 \boldsymbol{d}_{N+1} d N + 1
Continuous Mayer Optimal-Control Problem ¶ Consider the simplified Mayer problem
min x ( ⋅ ) , u ( ⋅ ) Φ ( x ( 1 ) ) , \min_{\boldsymbol{x}(\cdot),\boldsymbol{u}(\cdot)}
\Phi\bigl(\boldsymbol{x}(1)\bigr), x ( ⋅ ) , u ( ⋅ ) min Φ ( x ( 1 ) ) , subject to
x ˙ ( τ ) = f ( x ( τ ) , u ( τ ) ) , τ ∈ [ − 1 , 1 ] , \dot{\boldsymbol{x}}(\tau)
=
\boldsymbol{f}\bigl(\boldsymbol{x}(\tau),\boldsymbol{u}(\tau)\bigr),
\qquad
\tau\in[-1,1], x ˙ ( τ ) = f ( x ( τ ) , u ( τ ) ) , τ ∈ [ − 1 , 1 ] , and
x ( − 1 ) = x 0 . \boldsymbol{x}(-1)=\boldsymbol{x}_0. x ( − 1 ) = x 0 . Because the objective contains no running cost, the Hamiltonian is
H ( x , u , λ ) = λ T f ( x , u ) . H(\boldsymbol{x},\boldsymbol{u},\boldsymbol{\lambda})
=
\boldsymbol{\lambda}^{\mathsf{T}}\boldsymbol{f}(\boldsymbol{x},\boldsymbol{u}). H ( x , u , λ ) = λ T f ( x , u ) . Continuous First-Order Necessary Conditions ¶ The continuous first-order conditions are:
Terminal Transversality ¶ λ ( 1 ) = ∇ x Φ ( x ( 1 ) ) . \boxed{
\boldsymbol{\lambda}(1)
=
\nabla_{\boldsymbol{x}}
\Phi\bigl(\boldsymbol{x}(1)\bigr).
} λ ( 1 ) = ∇ x Φ ( x ( 1 ) ) . Initial Boundary Multiplier ¶ If μ \boldsymbol{\mu} μ
λ ( − 1 ) = μ . \boxed{
\boldsymbol{\lambda}(-1)=\boldsymbol{\mu}.
} λ ( − 1 ) = μ . Costate Equation ¶ λ ˙ ( τ ) = − ∇ x [ λ ( τ ) T f ( x ( τ ) , u ( τ ) ) ] . \boxed{
\dot{\boldsymbol{\lambda}}(\tau)
=
-
\nabla_{\boldsymbol{x}}
\left[
\boldsymbol{\lambda}(\tau)^{\mathsf{T}}
\boldsymbol{f}\bigl(\boldsymbol{x}(\tau),\boldsymbol{u}(\tau)\bigr)
\right].
} λ ˙ ( τ ) = − ∇ x [ λ ( τ ) T f ( x ( τ ) , u ( τ ) ) ] . Stationarity with Respect to the Control ¶ ∇ u [ λ ( τ ) T f ( x ( τ ) , u ( τ ) ) ] = 0 . \boxed{
\nabla_{\boldsymbol{u}}
\left[
\boldsymbol{\lambda}(\tau)^{\mathsf{T}}
\boldsymbol{f}\bigl(\boldsymbol{x}(\tau),\boldsymbol{u}(\tau)\bigr)
\right]
=
\boldsymbol{0}.
} ∇ u [ λ ( τ ) T f ( x ( τ ) , u ( τ ) ) ] = 0 . The discrete LGR KKT conditions will be transformed into approximations of these equations.
LGR Grid ¶ Let
τ 1 , … , τ N \tau_1,\ldots,\tau_N τ 1 , … , τ N denote the N N N
and
− 1 = τ 1 < τ 2 < ⋯ < τ N < 1. -1=\tau_1<\tau_2<\cdots<\tau_N<1. − 1 = τ 1 < τ 2 < ⋯ < τ N < 1. The noncollocated terminal interpolation point is
τ N + 1 = 1. \tau_{N+1}=1. τ N + 1 = 1. Thus:
the initial point is a collocation point;
the final point is included in the state interpolation;
the final point is not a collocation point.
State Approximation ¶ Approximate the state by a polynomial of degree at most N N N
X ( τ ) = ∑ i = 1 N + 1 X i L i ( τ ) , \boxed{
\boldsymbol{X}(\tau)
=
\sum_{i=1}^{N+1}
\boldsymbol{X}_iL_i(\tau),
} X ( τ ) = i = 1 ∑ N + 1 X i L i ( τ ) , where
L i ( τ ) = ∏ j = 1 j ≠ i N + 1 τ − τ j τ i − τ j . L_i(\tau)
=
\prod_{\substack{j=1\\j\neq i}}^{N+1}
\frac{\tau-\tau_j}{\tau_i-\tau_j}. L i ( τ ) = j = 1 j = i ∏ N + 1 τ i − τ j τ − τ j . The state values are
X i ≈ x ( τ i ) , i = 1 , … , N + 1. \boldsymbol{X}_i
\approx
\boldsymbol{x}(\tau_i),
\qquad
i=1,\ldots,N+1. X i ≈ x ( τ i ) , i = 1 , … , N + 1. The derivative of the approximation is
X ˙ ( τ ) = ∑ i = 1 N + 1 X i L i ′ ( τ ) . \dot{\boldsymbol{X}}(\tau)
=
\sum_{i=1}^{N+1}
\boldsymbol{X}_iL_i'(\tau). X ˙ ( τ ) = i = 1 ∑ N + 1 X i L i ′ ( τ ) . LGR Differentiation Matrix ¶ Define the LGR differentiation matrix by
D j i = L i ′ ( τ j ) , j = 1 , … , N , i = 1 , … , N + 1. D_{ji}
=
L_i'(\tau_j),
\qquad
j=1,\ldots,N,
\qquad
i=1,\ldots,N+1. D ji = L i ′ ( τ j ) , j = 1 , … , N , i = 1 , … , N + 1. Therefore,
D ∈ R N × ( N + 1 ) . \boxed{
\boldsymbol{D}\in\mathbb{R}^{N\times(N+1)}.
} D ∈ R N × ( N + 1 ) . Collocating at the N N N
∑ i = 1 N + 1 D j i X i = f ( X j , U j ) , j = 1 , … , N . \sum_{i=1}^{N+1}
D_{ji}\boldsymbol{X}_i
=
\boldsymbol{f}(\boldsymbol{X}_j,\boldsymbol{U}_j),
\qquad
j=1,\ldots,N. i = 1 ∑ N + 1 D ji X i = f ( X j , U j ) , j = 1 , … , N . In matrix form,
D X 1 : N + 1 = F 1 : N , \boxed{
\boldsymbol{D}\boldsymbol{X}_{1:N+1}
=
\boldsymbol{F}_{1:N},
} D X 1 : N + 1 = F 1 : N , where
F j = f ( X j , U j ) . \boldsymbol{F}_j
=
\boldsymbol{f}(\boldsymbol{X}_j,\boldsymbol{U}_j). F j = f ( X j , U j ) . LGR Nonlinear Program ¶ The LGR transcription is
min X , U Φ ( X N + 1 ) subject to F 1 : N − D X 1 : N + 1 = 0 , x 0 − X 1 = 0 . \boxed{
\begin{aligned}
\min_{\boldsymbol{X},\boldsymbol{U}}
\quad
&
\Phi(\boldsymbol{X}_{N+1})
\\
\text{subject to}
\quad
&
\boldsymbol{F}_{1:N}
-
\boldsymbol{D}\boldsymbol{X}_{1:N+1}
=
\boldsymbol{0},
\\
&
\boldsymbol{x}_0-\boldsymbol{X}_1
=
\boldsymbol{0}.
\end{aligned}
} X , U min subject to Φ ( X N + 1 ) F 1 : N − D X 1 : N + 1 = 0 , x 0 − X 1 = 0 . Unlike the LG method, no separate quadrature equation is needed to recover the terminal state, because X N + 1 \boldsymbol{X}_{N+1} X N + 1
LGR NLP Lagrangian ¶ Let
Λ 1 : N = [ Λ 1 ⋮ Λ N ] \boldsymbol{\Lambda}_{1:N}
=
\begin{bmatrix}
\boldsymbol{\Lambda}_1\\
\vdots\\
\boldsymbol{\Lambda}_N
\end{bmatrix} Λ 1 : N = ⎣ ⎡ Λ 1 ⋮ Λ N ⎦ ⎤ be the multipliers associated with the N N N
Let
be the multiplier associated with the initial condition.
The NLP Lagrangian is
L = Φ ( X N + 1 ) + ⟨ Λ 1 : N , F 1 : N − D X 1 : N + 1 ⟩ + ⟨ μ , x 0 − X 1 ⟩ . \begin{aligned}
\mathcal{L}
={}&
\Phi(\boldsymbol{X}_{N+1})
+
\left\langle
\boldsymbol{\Lambda}_{1:N},
\boldsymbol{F}_{1:N}
-
\boldsymbol{D}\boldsymbol{X}_{1:N+1}
\right\rangle
\nonumber\\
&+
\left\langle
\boldsymbol{\mu},
\boldsymbol{x}_0-\boldsymbol{X}_1
\right\rangle.
\end{aligned} L = Φ ( X N + 1 ) + ⟨ Λ 1 : N , F 1 : N − D X 1 : N + 1 ⟩ + ⟨ μ , x 0 − X 1 ⟩ . KKT Conditions for the LGR NLP ¶ The NLP variables are
X 1 , … , X N + 1 , U 1 , … , U N . \boldsymbol{X}_1,\ldots,\boldsymbol{X}_{N+1},
\qquad
\boldsymbol{U}_1,\ldots,\boldsymbol{U}_N. X 1 , … , X N + 1 , U 1 , … , U N . State Stationarity for Interior LGR Points ¶ For
j = 2 , … , N , j=2,\ldots,N, j = 2 , … , N , stationarity with respect to X j \boldsymbol{X}_j X j
∑ i = 1 N D i j Λ i = ( ∇ X f ( X j , U j ) ) T Λ j . \boxed{
\sum_{i=1}^{N}
D_{ij}\boldsymbol{\Lambda}_i
=
\left(
\nabla_{\boldsymbol{X}}
\boldsymbol{f}(\boldsymbol{X}_j,\boldsymbol{U}_j)
\right)^{\mathsf{T}}
\boldsymbol{\Lambda}_j.
} i = 1 ∑ N D ij Λ i = ( ∇ X f ( X j , U j ) ) T Λ j . In scalar shorthand,
∑ i = 1 N D i j Λ i = Λ j ∇ X f ( X j , U j ) . \sum_{i=1}^{N}
D_{ij}\Lambda_i
=
\Lambda_j
\nabla_{\boldsymbol{X}}f(\boldsymbol{X}_j,\boldsymbol{U}_j). i = 1 ∑ N D ij Λ i = Λ j ∇ X f ( X j , U j ) . State Stationarity at the Initial LGR Point ¶ Because X 1 \boldsymbol{X}_1 X 1
∑ i = 1 N D i 1 Λ i = ( ∇ X f ( X 1 , U 1 ) ) T Λ 1 − μ . \boxed{
\sum_{i=1}^{N}
D_{i1}\boldsymbol{\Lambda}_i
=
\left(
\nabla_{\boldsymbol{X}}
\boldsymbol{f}(\boldsymbol{X}_1,\boldsymbol{U}_1)
\right)^{\mathsf{T}}
\boldsymbol{\Lambda}_1
-
\boldsymbol{\mu}.
} i = 1 ∑ N D i 1 Λ i = ( ∇ X f ( X 1 , U 1 ) ) T Λ 1 − μ . Equivalently,
∑ i = 1 N D i 1 Λ i + μ = ( ∇ X f ( X 1 , U 1 ) ) T Λ 1 . \sum_{i=1}^{N}
D_{i1}\boldsymbol{\Lambda}_i
+
\boldsymbol{\mu}
=
\left(
\nabla_{\boldsymbol{X}}
\boldsymbol{f}(\boldsymbol{X}_1,\boldsymbol{U}_1)
\right)^{\mathsf{T}}
\boldsymbol{\Lambda}_1. i = 1 ∑ N D i 1 Λ i + μ = ( ∇ X f ( X 1 , U 1 ) ) T Λ 1 . Control Stationarity ¶ For
j = 1 , … , N , j=1,\ldots,N, j = 1 , … , N , stationarity with respect to U j \boldsymbol{U}_j U j
( ∇ U f ( X j , U j ) ) T Λ j = 0 . \boxed{
\left(
\nabla_{\boldsymbol{U}}
\boldsymbol{f}(\boldsymbol{X}_j,\boldsymbol{U}_j)
\right)^{\mathsf{T}}
\boldsymbol{\Lambda}_j
=
\boldsymbol{0}.
} ( ∇ U f ( X j , U j ) ) T Λ j = 0 . Terminal-State Stationarity ¶ Differentiation with respect to X N + 1 \boldsymbol{X}_{N+1} X N + 1
∇ X Φ ( X N + 1 ) = d N + 1 T Λ 1 : N , \boxed{
\nabla_{\boldsymbol{X}}
\Phi(\boldsymbol{X}_{N+1})
=
\boldsymbol{d}_{N+1}^{\mathsf{T}}
\boldsymbol{\Lambda}_{1:N},
} ∇ X Φ ( X N + 1 ) = d N + 1 T Λ 1 : N , where
d N + 1 \boldsymbol{d}_{N+1} d N + 1 is the final column of D \boldsymbol{D} D
This equation provides the terminal costate estimate.
Definition of the LGR Adjoint Differentiation Matrix ¶ Let
W = diag ( w 1 , … , w N ) , \boldsymbol{W}
=
\operatorname{diag}(w_1,\ldots,w_N), W = diag ( w 1 , … , w N ) , where w i w_i w i
Define
D † ∈ R N × N . \boldsymbol{D}^\dagger\in\mathbb{R}^{N\times N}. D † ∈ R N × N . A commonly used representation is
D 11 † = − D 11 − 1 w 1 , \boxed{
D_{11}^\dagger
=
-D_{11}-\frac{1}{w_1},
} D 11 † = − D 11 − w 1 1 , and, for the remaining entries,
D i j † = − w j w i D j i . \boxed{
D_{ij}^\dagger
=
-
\frac{w_j}{w_i}D_{ji}.
} D ij † = − w i w j D ji . This matrix plays the role of a differentiation matrix for the discrete costate approximation at the LGR points.
Define the costate estimate at the LGR points by
λ 1 : N = W − 1 Λ 1 : N . \boxed{
\boldsymbol{\lambda}_{1:N}
=
\boldsymbol{W}^{-1}\boldsymbol{\Lambda}_{1:N}.
} λ 1 : N = W − 1 Λ 1 : N . Componentwise,
λ i = Λ i w i , i = 1 , … , N . \boxed{
\boldsymbol{\lambda}_i
=
\frac{\boldsymbol{\Lambda}_i}{w_i},
\qquad
i=1,\ldots,N.
} λ i = w i Λ i , i = 1 , … , N . This mapping is simpler than the LG mapping because there is no additional terminal quadrature multiplier.
Substitute
Λ i = w i λ i \boldsymbol{\Lambda}_i=w_i\boldsymbol{\lambda}_i Λ i = w i λ i into the KKT state equations.
For the interior points, the transformed equation becomes
∑ j = 1 N D i j † λ j = − ∇ X [ λ i T f ( X i , U i ) ] , i = 2 , … , N . \boxed{
\sum_{j=1}^{N}
D_{ij}^\dagger\boldsymbol{\lambda}_j
=
-
\nabla_{\boldsymbol{X}}
\left[
\boldsymbol{\lambda}_i^{\mathsf{T}}
\boldsymbol{f}(\boldsymbol{X}_i,\boldsymbol{U}_i)
\right],
\qquad
i=2,\ldots,N.
} j = 1 ∑ N D ij † λ j = − ∇ X [ λ i T f ( X i , U i ) ] , i = 2 , … , N . This is the discrete LGR approximation of
λ ˙ = − ∇ x ( λ T f ) . \dot{\boldsymbol{\lambda}}
=
-
\nabla_{\boldsymbol{x}}
\left(
\boldsymbol{\lambda}^{\mathsf{T}}\boldsymbol{f}
\right). λ ˙ = − ∇ x ( λ T f ) . At the first LGR point,
∑ j = 1 N D 1 j † λ j = − ∇ X [ λ 1 T f ( X 1 , U 1 ) ] + 1 w 1 ( μ − λ 1 ) . \boxed{
\sum_{j=1}^{N}
D_{1j}^\dagger\boldsymbol{\lambda}_j
=
-
\nabla_{\boldsymbol{X}}
\left[
\boldsymbol{\lambda}_1^{\mathsf{T}}
\boldsymbol{f}(\boldsymbol{X}_1,\boldsymbol{U}_1)
\right]
+
\frac{1}{w_1}
\left(
\boldsymbol{\mu}-\boldsymbol{\lambda}_1
\right).
} j = 1 ∑ N D 1 j † λ j = − ∇ X [ λ 1 T f ( X 1 , U 1 ) ] + w 1 1 ( μ − λ 1 ) . The additional term reflects the initial-condition multiplier.
Asymptotic Endpoint Consistency ¶ The continuous boundary condition is
λ ( − 1 ) = μ . \boldsymbol{\lambda}(-1)=\boldsymbol{\mu}. λ ( − 1 ) = μ . Because
the first LGR costate estimate should approximate the initial costate:
λ 1 ≈ μ . \boldsymbol{\lambda}_1\approx\boldsymbol{\mu}. λ 1 ≈ μ . As the polynomial degree increases,
μ − λ 1 → 0 . \boxed{
\boldsymbol{\mu}-\boldsymbol{\lambda}_1\rightarrow\boldsymbol{0}.
} μ − λ 1 → 0 . Consequently, the correction term
1 w 1 ( μ − λ 1 ) \frac{1}{w_1}
\left(
\boldsymbol{\mu}-\boldsymbol{\lambda}_1
\right) w 1 1 ( μ − λ 1 ) vanishes asymptotically.
Terminal Costate Recovery ¶ The terminal-state KKT condition is
∇ X Φ ( X N + 1 ) = d N + 1 T Λ 1 : N . \nabla_{\boldsymbol{X}}
\Phi(\boldsymbol{X}_{N+1})
=
\boldsymbol{d}_{N+1}^{\mathsf{T}}
\boldsymbol{\Lambda}_{1:N}. ∇ X Φ ( X N + 1 ) = d N + 1 T Λ 1 : N . From the continuous transversality condition,
λ ( 1 ) = ∇ x Φ ( x ( 1 ) ) . \boldsymbol{\lambda}(1)
=
\nabla_{\boldsymbol{x}}
\Phi(\boldsymbol{x}(1)). λ ( 1 ) = ∇ x Φ ( x ( 1 )) . Therefore, define
λ N + 1 = d N + 1 T Λ 1 : N . \boxed{
\boldsymbol{\lambda}_{N+1}
=
\boldsymbol{d}_{N+1}^{\mathsf{T}}
\boldsymbol{\Lambda}_{1:N}.
} λ N + 1 = d N + 1 T Λ 1 : N . This gives the costate estimate at the noncollocated final point.
Alternative Interpretation Through Quadrature ¶ Since
λ ˙ = − ∇ x ( λ T f ) , \dot{\boldsymbol{\lambda}}
=
-
\nabla_{\boldsymbol{x}}
\left(
\boldsymbol{\lambda}^{\mathsf{T}}\boldsymbol{f}
\right), λ ˙ = − ∇ x ( λ T f ) , integration gives
λ ( 1 ) = λ ( − 1 ) + ∫ − 1 1 λ ˙ ( τ ) d τ . \boldsymbol{\lambda}(1)
=
\boldsymbol{\lambda}(-1)
+
\int_{-1}^{1}
\dot{\boldsymbol{\lambda}}(\tau)\,\mathrm{d}\tau. λ ( 1 ) = λ ( − 1 ) + ∫ − 1 1 λ ˙ ( τ ) d τ . Using LGR quadrature,
λ N + 1 ≈ λ 1 + ∑ i = 1 N w i [ − ∇ X ( λ i T F i ) ] . \boldsymbol{\lambda}_{N+1}
\approx
\boldsymbol{\lambda}_1
+
\sum_{i=1}^{N}
w_i
\left[
-
\nabla_{\boldsymbol{X}}
\left(
\boldsymbol{\lambda}_i^{\mathsf{T}}\boldsymbol{F}_i
\right)
\right]. λ N + 1 ≈ λ 1 + i = 1 ∑ N w i [ − ∇ X ( λ i T F i ) ] . When
λ 1 ≈ μ , \boldsymbol{\lambda}_1\approx\boldsymbol{\mu}, λ 1 ≈ μ , this quadrature relation is consistent with the terminal formula obtained directly from the KKT conditions.
Discrete Control Stationarity ¶ Using
Λ i = w i λ i , \boldsymbol{\Lambda}_i=w_i\boldsymbol{\lambda}_i, Λ i = w i λ i , the control stationarity conditions become
( ∇ U f ( X i , U i ) ) T w i λ i = 0 . \left(
\nabla_{\boldsymbol{U}}
\boldsymbol{f}(\boldsymbol{X}_i,\boldsymbol{U}_i)
\right)^{\mathsf{T}}
w_i\boldsymbol{\lambda}_i
=
\boldsymbol{0}. ( ∇ U f ( X i , U i ) ) T w i λ i = 0 . Since w i > 0 w_i>0 w i > 0
∇ U [ λ i T f ( X i , U i ) ] = 0 , i = 1 , … , N . \boxed{
\nabla_{\boldsymbol{U}}
\left[
\boldsymbol{\lambda}_i^{\mathsf{T}}
\boldsymbol{f}(\boldsymbol{X}_i,\boldsymbol{U}_i)
\right]
=
\boldsymbol{0},
\qquad
i=1,\ldots,N.
} ∇ U [ λ i T f ( X i , U i ) ] = 0 , i = 1 , … , N . This is the discrete equivalent of
∇ u H = 0 . \nabla_{\boldsymbol{u}}H=\boldsymbol{0}. ∇ u H = 0 . Complete LGR Costate Mapping Theorem ¶ Let Λ 1 : N \boldsymbol{\Lambda}_{1:N} Λ 1 : N W \boldsymbol{W} W d N + 1 \boldsymbol{d}_{N+1} d N + 1
λ 1 : N = W − 1 Λ 1 : N , \boldsymbol{\lambda}_{1:N}
=
\boldsymbol{W}^{-1}\boldsymbol{\Lambda}_{1:N}, λ 1 : N = W − 1 Λ 1 : N , and
λ N + 1 = d N + 1 T Λ 1 : N . \boldsymbol{\lambda}_{N+1}
=
\boldsymbol{d}_{N+1}^{\mathsf{T}}\boldsymbol{\Lambda}_{1:N}. λ N + 1 = d N + 1 T Λ 1 : N . The first point satisfies
λ 1 ≈ μ , \boldsymbol{\lambda}_1\approx\boldsymbol{\mu}, λ 1 ≈ μ , with the discrepancy tending to zero as the approximation is refined.
Computational Algorithm ¶ After solving the LGR NLP:
Extract the multipliers associated with the collocation constraints:
Λ 1 : N . \boldsymbol{\Lambda}_{1:N}. Λ 1 : N . Form the diagonal weight matrix:
W = diag ( w 1 , … , w N ) . \boldsymbol{W}
=
\operatorname{diag}(w_1,\ldots,w_N). W = diag ( w 1 , … , w N ) . Compute the costates at the LGR points:
λ 1 : N = W − 1 Λ 1 : N . \boldsymbol{\lambda}_{1:N}
=
\boldsymbol{W}^{-1}\boldsymbol{\Lambda}_{1:N}. λ 1 : N = W − 1 Λ 1 : N . Extract the final differentiation-matrix column:
d N + 1 = D ( : , N + 1 ) . \boldsymbol{d}_{N+1}
=
\boldsymbol{D}(:,N+1). d N + 1 = D ( : , N + 1 ) . Compute the terminal costate:
λ N + 1 = d N + 1 T Λ 1 : N . \boldsymbol{\lambda}_{N+1}
=
\boldsymbol{d}_{N+1}^{\mathsf{T}}\boldsymbol{\Lambda}_{1:N}. λ N + 1 = d N + 1 T Λ 1 : N . Compare
λ 1 and μ \boldsymbol{\lambda}_1
\quad\text{and}\quad
\boldsymbol{\mu} λ 1 and μ as an endpoint-consistency check.
Check the discrete adjoint residual:
D † λ 1 : N + ∇ X ⟨ λ 1 : N , F 1 : N ⟩ . \boldsymbol{D}^\dagger\boldsymbol{\lambda}_{1:N}
+
\nabla_{\boldsymbol{X}}
\left\langle
\boldsymbol{\lambda}_{1:N},
\boldsymbol{F}_{1:N}
\right\rangle. D † λ 1 : N + ∇ X ⟨ λ 1 : N , F 1 : N ⟩ . Verify convergence under polynomial-order or mesh refinement.
Comparison with the LG Costate Mapping ¶ Feature LG method LGR method Collocation points Interior Gauss points Includes one endpoint Initial point Noncollocated interpolation point Collocation point Final point Noncollocated and recovered by quadrature Noncollocated state variable Terminal quadrature constraint Required Not required Interior costate mapping W − 1 Λ + 1 Λ N + 1 \boldsymbol{W}^{-1}\boldsymbol{\Lambda}+\boldsymbol{1}\boldsymbol{\Lambda}_{N+1} W − 1 Λ + 1 Λ N + 1 W − 1 Λ \boldsymbol{W}^{-1}\boldsymbol{\Lambda} W − 1 Λ Terminal costate Terminal quadrature multiplier Final column of D \boldsymbol{D} D Algebraic complexity Higher Lower
Comparison of LG and LGR costate estimation.
Sign Conventions ¶ The formulas above assume the collocation constraints are written as
F − D X = 0 . \boldsymbol{F}-\boldsymbol{D}\boldsymbol{X}=\boldsymbol{0}. F − D X = 0 . If they are instead written as
D X − F = 0 , \boldsymbol{D}\boldsymbol{X}-\boldsymbol{F}=\boldsymbol{0}, D X − F = 0 , the reported KKT multipliers may have the opposite sign. Therefore, the implementation must be consistent with the exact NLP constraint convention.
Scaling and Multiplier Accuracy ¶ Accurate primal states do not automatically guarantee accurate dual variables.
Potential sources of poor multiplier accuracy include:
badly scaled state variables;
badly scaled dynamic defects;
loose NLP tolerances;
insufficient polynomial degree;
poor mesh placement;
inconsistent derivative information;
active path constraints near switching points.
A reliable implementation should check both primal and dual convergence.
Interpretation of the Costate ¶ The costate can be interpreted as a sensitivity of the optimal value to perturbations of the state. In the LGR transcription, the mapped KKT multipliers provide a discrete approximation of this sensitivity along the trajectory.
Costate estimates are useful for:
verifying first-order necessary conditions;
interpreting state importance;
identifying switching structure;
comparing direct and indirect methods;
constructing neighboring optimal-control laws;
sensitivity analysis;
mesh refinement diagnostics.
Summary ¶ The LGR state polynomial uses the N N N
The LGR NLP contains collocation constraints and an initial-condition constraint.
No terminal quadrature constraint is required.
The KKT multipliers of the collocation equations are scaled by the inverse Radau weights.
The interior costate mapping is
λ 1 : N = W − 1 Λ 1 : N . \boldsymbol{\lambda}_{1:N}
=
\boldsymbol{W}^{-1}\boldsymbol{\Lambda}_{1:N}. λ 1 : N = W − 1 Λ 1 : N . The terminal costate is recovered from
λ N + 1 = d N + 1 T Λ 1 : N . \boldsymbol{\lambda}_{N+1}
=
\boldsymbol{d}_{N+1}^{\mathsf{T}}\boldsymbol{\Lambda}_{1:N}. λ N + 1 = d N + 1 T Λ 1 : N . The initial costate estimate satisfies
λ 1 ≈ μ . \boldsymbol{\lambda}_1\approx\boldsymbol{\mu}. λ 1 ≈ μ . The transformed KKT equations reproduce the continuous adjoint and control-stationarity conditions in discrete form.
Connection. Completing the LGL case makes it possible to compare all three Legendre schemes and then restore the full Bolza, endpoint, and path-constraint terms.