Global Assembly of Multiple-Interval LGR Transcriptions
This section completes the derivation of the multiple-interval Legendre–Gauss–Radau (LGR) pseudospectral method. The preceding development established the local state approximation, local differentiation matrix, local quadrature rule, and collocation equations on each mesh interval. The present chapter assembles those local objects into a single global nonlinear programming problem.
The central result is that the multiple-interval method retains almost exactly the same mathematical form as the single-interval LGR method:
D X = t f − t 0 2 F , \boldsymbol{D}\boldsymbol{X}
=
\frac{t_f-t_0}{2}\boldsymbol{F}, D X = 2 t f − t 0 F , but the global differentiation matrix is sparse and is constructed from local interval matrices whose adjacent blocks share one state column.
The chapter develops:
global stacking of LGR nodes, weights, states, and controls;
shared state variables at mesh interfaces;
the unique control representation at interval interfaces;
assembly of the global differentiation matrix;
compact quadrature and cost expressions;
boundary and path constraints;
the final sparse nonlinear programming formulation;
implementation and verification procedures;
motivation for adaptive mesh refinement.
Transformed Time Domain and Mesh ¶ The physical time interval
t ∈ [ t 0 , t f ] t\in[t_0,t_f] t ∈ [ t 0 , t f ] is mapped to the fixed domain
s ∈ [ − 1 , 1 ] s\in[-1,1] s ∈ [ − 1 , 1 ] using
s = 2 ( t − t 0 ) t f − t 0 − 1. s
=
\frac{2(t-t_0)}{t_f-t_0}-1. s = t f − t 0 2 ( t − t 0 ) − 1. The inverse mapping is
t ( s ) = t f − t 0 2 s + t f + t 0 2 . t(s)
=
\frac{t_f-t_0}{2}s
+
\frac{t_f+t_0}{2}. t ( s ) = 2 t f − t 0 s + 2 t f + t 0 . The transformed dynamics are
d x d s = t f − t 0 2 f ( x ( s ) , u ( s ) , t ( s ) , p ) . \frac{\,\mathrm{d}\boldsymbol{x}}{\,\mathrm{d} s}
=
\frac{t_f-t_0}{2}
\boldsymbol{f}\left(
\boldsymbol{x}(s),\boldsymbol{u}(s),t(s),\boldsymbol{p}
\right). d s d x = 2 t f − t 0 f ( x ( s ) , u ( s ) , t ( s ) , p ) . The transformed domain is divided into K K K
− 1 = s 0 < s 1 < ⋯ < s K − 1 < s K = 1. -1=s_0<s_1<\cdots<s_{K-1}<s_K=1. − 1 = s 0 < s 1 < ⋯ < s K − 1 < s K = 1. The k k k
S k = [ s k − 1 , s k ] . \mathcal{S}_k=[s_{k-1},s_k]. S k = [ s k − 1 , s k ] . The mesh points are fixed during a single nonlinear programming solve. They may be changed between solves during mesh refinement.
Local LGR Approximation ¶ Let interval k k k N k N_k N k
s 1 ( k ) , s 2 ( k ) , … , s N k ( k ) . s_1^{(k)},s_2^{(k)},\ldots,s_{N_k}^{(k)}. s 1 ( k ) , s 2 ( k ) , … , s N k ( k ) . Forward LGR includes the left interval endpoint but excludes the right endpoint:
s 1 ( k ) = s k − 1 , s N k ( k ) < s k . s_1^{(k)}=s_{k-1},
\qquad
s_{N_k}^{(k)}<s_k. s 1 ( k ) = s k − 1 , s N k ( k ) < s k . The additional state point is
s N k + 1 ( k ) = s k . s_{N_k+1}^{(k)}=s_k. s N k + 1 ( k ) = s k . The state approximation is
x ( k ) ( s ) ≈ X ( k ) ( s ) = ∑ i = 1 N k + 1 X i ( k ) ℓ i ( k ) ( s ) , \boldsymbol{x}^{(k)}(s)
\approx
\boldsymbol{X}^{(k)}(s)
=
\sum_{i=1}^{N_k+1}
\boldsymbol{X}_i^{(k)}
\ell_i^{(k)}(s), x ( k ) ( s ) ≈ X ( k ) ( s ) = i = 1 ∑ N k + 1 X i ( k ) ℓ i ( k ) ( s ) , where
ℓ i ( k ) ( s ) = ∏ j = 1 j ≠ i N k + 1 s − s j ( k ) s i ( k ) − s j ( k ) . \ell_i^{(k)}(s)
=
\prod_{\substack{j=1\\j\neq i}}^{N_k+1}
\frac{s-s_j^{(k)}}{s_i^{(k)}-s_j^{(k)}}. ℓ i ( k ) ( s ) = j = 1 j = i ∏ N k + 1 s i ( k ) − s j ( k ) s − s j ( k ) . The control is represented at only the N k N_k N k
U i ( k ) ≈ u ( s i ( k ) ) , i = 1 , … , N k . \boldsymbol{U}_i^{(k)}
\approx
\boldsymbol{u}\left(s_i^{(k)}\right),
\qquad
i=1,\ldots,N_k. U i ( k ) ≈ u ( s i ( k ) ) , i = 1 , … , N k . Local Differentiation Matrix ¶ Differentiating the state approximation gives
d X ( k ) d s = ∑ i = 1 N k + 1 X i ( k ) d ℓ i ( k ) d s . \frac{\,\mathrm{d}\boldsymbol{X}^{(k)}}{\,\mathrm{d} s}
=
\sum_{i=1}^{N_k+1}
\boldsymbol{X}_i^{(k)}
\frac{\,\mathrm{d}\ell_i^{(k)}}{\,\mathrm{d} s}. d s d X ( k ) = i = 1 ∑ N k + 1 X i ( k ) d s d ℓ i ( k ) . At collocation node s j ( k ) s_j^{(k)} s j ( k )
d X ( k ) d s ∣ s = s j ( k ) = ∑ i = 1 N k + 1 D j i ( k ) X i ( k ) , \left.
\frac{\,\mathrm{d}\boldsymbol{X}^{(k)}}{\,\mathrm{d} s}
\right|_{s=s_j^{(k)}}
=
\sum_{i=1}^{N_k+1}
D_{ji}^{(k)}
\boldsymbol{X}_i^{(k)}, d s d X ( k ) ∣ ∣ s = s j ( k ) = i = 1 ∑ N k + 1 D ji ( k ) X i ( k ) , where
D j i ( k ) = d ℓ i ( k ) d s ∣ s = s j ( k ) . D_{ji}^{(k)}
=
\left.
\frac{\,\mathrm{d}\ell_i^{(k)}}{\,\mathrm{d} s}
\right|_{s=s_j^{(k)}}. D ji ( k ) = d s d ℓ i ( k ) ∣ ∣ s = s j ( k ) . Thus,
D ( k ) ∈ R N k × ( N k + 1 ) . \boldsymbol{D}^{(k)}
\in
\mathbb{R}^{N_k\times(N_k+1)}. D ( k ) ∈ R N k × ( N k + 1 ) . The local collocation equations are
∑ i = 1 N k + 1 D j i ( k ) X i ( k ) = t f − t 0 2 f ( X j ( k ) , U j ( k ) , t j ( k ) , p ) . \boxed{
\sum_{i=1}^{N_k+1}
D_{ji}^{(k)}
\boldsymbol{X}_i^{(k)}
=
\frac{t_f-t_0}{2}
\boldsymbol{f}\left(
\boldsymbol{X}_j^{(k)},
\boldsymbol{U}_j^{(k)},
t_j^{(k)},
\boldsymbol{p}
\right).
} i = 1 ∑ N k + 1 D ji ( k ) X i ( k ) = 2 t f − t 0 f ( X j ( k ) , U j ( k ) , t j ( k ) , p ) . These equations hold for
j = 1 , … , N k . j=1,\ldots,N_k. j = 1 , … , N k . State Continuity at Mesh Points ¶ The final state of interval k k k k + 1 k+1 k + 1
X N k + 1 ( k ) = X 1 ( k + 1 ) . \boxed{
\boldsymbol{X}_{N_k+1}^{(k)}
=
\boldsymbol{X}_1^{(k+1)}.
} X N k + 1 ( k ) = X 1 ( k + 1 ) . This relation can be implemented in two ways.
Explicit Continuity Constraints ¶ Retain separate variables and impose
X N k + 1 ( k ) − X 1 ( k + 1 ) = 0 . \boldsymbol{X}_{N_k+1}^{(k)}
-
\boldsymbol{X}_1^{(k+1)}
=
\boldsymbol{0}. X N k + 1 ( k ) − X 1 ( k + 1 ) = 0 . Shared Interface Variables ¶ Use one NLP variable for both quantities:
X N k + 1 ( k ) ≡ X 1 ( k + 1 ) . \boldsymbol{X}_{N_k+1}^{(k)}
\equiv
\boldsymbol{X}_1^{(k+1)}. X N k + 1 ( k ) ≡ X 1 ( k + 1 ) . The shared-variable formulation is particularly elegant because it creates the overlapping-column structure of the global differentiation matrix.
Control Representation at Mesh Points ¶ The forward LGR scheme does not contain a control at the right endpoint of interval k k k k + 1 k+1 k + 1
Therefore, the control at an interior mesh point has one representation:
U 1 ( k + 1 ) . \boldsymbol{U}_1^{(k+1)}. U 1 ( k + 1 ) . There is no second control variable
U N k + 1 ( k ) \boldsymbol{U}_{N_k+1}^{(k)} U N k + 1 ( k ) because the terminal state point is not a collocation point.
This property avoids an artificial double-valued control at the interface.
The only point on the full transformed domain without a control variable is the final endpoint:
Total Number of LGR Points ¶ Define the total number of collocation points:
N = ∑ k = 1 K N k . \boxed{
N
=
\sum_{k=1}^{K}N_k.
} N = k = 1 ∑ K N k . There are:
N N N
N + 1 N+1 N + 1
N N N
N N N
The additional state point is the terminal point of the final interval.
Global Collocation-Point Vector ¶ Stack the local LGR nodes in monotonically increasing order:
S ~ = [ s 1 ( 1 ) ⋮ s N 1 ( 1 ) s 1 ( 2 ) ⋮ s N 2 ( 2 ) ⋮ s 1 ( K ) ⋮ s N K ( K ) ] ∈ R N . \widetilde{\boldsymbol{S}}
=
\begin{bmatrix}
s_1^{(1)}\\
\vdots\\
s_{N_1}^{(1)}\\
s_1^{(2)}\\
\vdots\\
s_{N_2}^{(2)}\\
\vdots\\
s_1^{(K)}\\
\vdots\\
s_{N_K}^{(K)}
\end{bmatrix}
\in\mathbb{R}^N. S = ⎣ ⎡ s 1 ( 1 ) ⋮ s N 1 ( 1 ) s 1 ( 2 ) ⋮ s N 2 ( 2 ) ⋮ s 1 ( K ) ⋮ s N K ( K ) ⎦ ⎤ ∈ R N . The complete state-node vector includes the terminal point:
S = [ S ~ s N K + 1 ( K ) ] = [ S ~ 1 ] ∈ R N + 1 . \boldsymbol{S}
=
\begin{bmatrix}
\widetilde{\boldsymbol{S}}\\
s_{N_K+1}^{(K)}
\end{bmatrix}
=
\begin{bmatrix}
\widetilde{\boldsymbol{S}}\\
1
\end{bmatrix}
\in\mathbb{R}^{N+1}. S = [ S s N K + 1 ( K ) ] = [ S 1 ] ∈ R N + 1 . All entries in S ~ \widetilde{\boldsymbol{S}} S S \boldsymbol{S} S
Global Weight Vector ¶ For interval k k k
w 1 ( k ) , … , w N k ( k ) . w_1^{(k)},\ldots,w_{N_k}^{(k)}. w 1 ( k ) , … , w N k ( k ) . Stack them:
W = [ w 1 ( 1 ) ⋮ w N 1 ( 1 ) w 1 ( 2 ) ⋮ w N 2 ( 2 ) ⋮ w 1 ( K ) ⋮ w N K ( K ) ] ∈ R N . \boxed{
\boldsymbol{W}
=
\begin{bmatrix}
w_1^{(1)}\\
\vdots\\
w_{N_1}^{(1)}\\
w_1^{(2)}\\
\vdots\\
w_{N_2}^{(2)}\\
\vdots\\
w_1^{(K)}\\
\vdots\\
w_{N_K}^{(K)}
\end{bmatrix}
\in\mathbb{R}^N.
} W = ⎣ ⎡ w 1 ( 1 ) ⋮ w N 1 ( 1 ) w 1 ( 2 ) ⋮ w N 2 ( 2 ) ⋮ w 1 ( K ) ⋮ w N K ( K ) ⎦ ⎤ ∈ R N . There is no weight corresponding to the terminal state point.
Global State Matrix ¶ Assume there are n x n_x n x
X = [ X 1 ( 1 ) ⋮ X N 1 ( 1 ) X 1 ( 2 ) ⋮ X N 2 ( 2 ) ⋮ X 1 ( K ) ⋮ X N K ( K ) X N K + 1 ( K ) ] . \boldsymbol{X}
=
\begin{bmatrix}
\boldsymbol{X}_1^{(1)}\\
\vdots\\
\boldsymbol{X}_{N_1}^{(1)}\\
\boldsymbol{X}_1^{(2)}\\
\vdots\\
\boldsymbol{X}_{N_2}^{(2)}\\
\vdots\\
\boldsymbol{X}_1^{(K)}\\
\vdots\\
\boldsymbol{X}_{N_K}^{(K)}\\
\boldsymbol{X}_{N_K+1}^{(K)}
\end{bmatrix}. X = ⎣ ⎡ X 1 ( 1 ) ⋮ X N 1 ( 1 ) X 1 ( 2 ) ⋮ X N 2 ( 2 ) ⋮ X 1 ( K ) ⋮ X N K ( K ) X N K + 1 ( K ) ⎦ ⎤ . With one row per state node,
X ∈ R ( N + 1 ) × n x . \boldsymbol{X}\in\mathbb{R}^{(N+1)\times n_x}. X ∈ R ( N + 1 ) × n x . The interface states are not duplicated. For example,
X N 1 + 1 ( 1 ) = X 1 ( 2 ) \boldsymbol{X}_{N_1+1}^{(1)}
=
\boldsymbol{X}_1^{(2)} X N 1 + 1 ( 1 ) = X 1 ( 2 ) is represented by one shared row in the global state matrix.
Collocation-State Matrix ¶ Path constraints and dynamics are evaluated only at collocation nodes. Define
X ~ = X ( 1 : N , : ) ∈ R N × n x . \boxed{
\widetilde{\boldsymbol{X}}
=
\boldsymbol{X}(1:N,:)
\in\mathbb{R}^{N\times n_x}.
} X = X ( 1 : N , : ) ∈ R N × n x . Thus, X ~ \widetilde{\boldsymbol{X}} X
Global Control Matrix ¶ Assume there are n u n_u n u
U = [ U 1 ( 1 ) ⋮ U N 1 ( 1 ) U 1 ( 2 ) ⋮ U N 2 ( 2 ) ⋮ U 1 ( K ) ⋮ U N K ( K ) ] ∈ R N × n u . \boldsymbol{U}
=
\begin{bmatrix}
\boldsymbol{U}_1^{(1)}\\
\vdots\\
\boldsymbol{U}_{N_1}^{(1)}\\
\boldsymbol{U}_1^{(2)}\\
\vdots\\
\boldsymbol{U}_{N_2}^{(2)}\\
\vdots\\
\boldsymbol{U}_1^{(K)}\\
\vdots\\
\boldsymbol{U}_{N_K}^{(K)}
\end{bmatrix}
\in
\mathbb{R}^{N\times n_u}. U = ⎣ ⎡ U 1 ( 1 ) ⋮ U N 1 ( 1 ) U 1 ( 2 ) ⋮ U N 2 ( 2 ) ⋮ U 1 ( K ) ⋮ U N K ( K ) ⎦ ⎤ ∈ R N × n u . The number of control rows equals the number of global collocation points.
Global Dynamics Matrix ¶ Define
F = [ f ( X 1 ( 1 ) , U 1 ( 1 ) , t 1 ( 1 ) , p ) ⋮ f ( X N 1 ( 1 ) , U N 1 ( 1 ) , t N 1 ( 1 ) , p ) ⋮ f ( X N K ( K ) , U N K ( K ) , t N K ( K ) , p ) ] . \boldsymbol{F}
=
\begin{bmatrix}
\boldsymbol{f}\left(
\boldsymbol{X}_1^{(1)},\boldsymbol{U}_1^{(1)},t_1^{(1)},\boldsymbol{p}
\right)\\
\vdots\\
\boldsymbol{f}\left(
\boldsymbol{X}_{N_1}^{(1)},\boldsymbol{U}_{N_1}^{(1)},t_{N_1}^{(1)},\boldsymbol{p}
\right)\\
\vdots\\
\boldsymbol{f}\left(
\boldsymbol{X}_{N_K}^{(K)},\boldsymbol{U}_{N_K}^{(K)},t_{N_K}^{(K)},\boldsymbol{p}
\right)
\end{bmatrix}. F = ⎣ ⎡ f ( X 1 ( 1 ) , U 1 ( 1 ) , t 1 ( 1 ) , p ) ⋮ f ( X N 1 ( 1 ) , U N 1 ( 1 ) , t N 1 ( 1 ) , p ) ⋮ f ( X N K ( K ) , U N K ( K ) , t N K ( K ) , p ) ⎦ ⎤ . Therefore,
F ∈ R N × n x . \boldsymbol{F}
\in
\mathbb{R}^{N\times n_x}. F ∈ R N × n x . Assembly of the Global Differentiation Matrix ¶ Each local matrix is
D ( k ) ∈ R N k × ( N k + 1 ) . \boldsymbol{D}^{(k)}
\in
\mathbb{R}^{N_k\times(N_k+1)}. D ( k ) ∈ R N k × ( N k + 1 ) . The final state column of interval k k k k + 1 k+1 k + 1
The global matrix has dimension
D ∈ R N × ( N + 1 ) . \boxed{
\boldsymbol{D}
\in
\mathbb{R}^{N\times(N+1)}.
} D ∈ R N × ( N + 1 ) . It is not a conventional block-diagonal matrix because adjacent blocks share one column. It is more accurately described as a block-sparse matrix with one-column overlap.
Two-Interval Example ¶ For two intervals with N 1 N_1 N 1 N 2 N_2 N 2
D = [ D ( 1 ) 0 0 D ( 2 ) ] with shared interface column . \boldsymbol{D}
=
\begin{bmatrix}
\boldsymbol{D}^{(1)} & \boldsymbol{0}\\
\boldsymbol{0} & \boldsymbol{D}^{(2)}
\end{bmatrix}_{\text{with shared interface column}}. D = [ D ( 1 ) 0 0 D ( 2 ) ] with shared interface column . More explicitly, partition
D ( 1 ) = [ D A ( 1 ) d b ( 1 ) ] , \boldsymbol{D}^{(1)}
=
\begin{bmatrix}
\boldsymbol{D}_A^{(1)} & \boldsymbol{d}_b^{(1)}
\end{bmatrix}, D ( 1 ) = [ D A ( 1 ) d b ( 1 ) ] , where
D A ( 1 ) ∈ R N 1 × N 1 , d b ( 1 ) ∈ R N 1 . \boldsymbol{D}_A^{(1)}\in\mathbb{R}^{N_1\times N_1},
\qquad
\boldsymbol{d}_b^{(1)}\in\mathbb{R}^{N_1}. D A ( 1 ) ∈ R N 1 × N 1 , d b ( 1 ) ∈ R N 1 . Similarly,
D ( 2 ) = [ d a ( 2 ) D B ( 2 ) ] . \boldsymbol{D}^{(2)}
=
\begin{bmatrix}
\boldsymbol{d}_a^{(2)} & \boldsymbol{D}_B^{(2)}
\end{bmatrix}. D ( 2 ) = [ d a ( 2 ) D B ( 2 ) ] . Then
D = [ D A ( 1 ) d b ( 1 ) 0 0 d a ( 2 ) D B ( 2 ) ] . \boxed{
\boldsymbol{D}
=
\begin{bmatrix}
\boldsymbol{D}_A^{(1)}
&
\boldsymbol{d}_b^{(1)}
&
\boldsymbol{0}
\\
\boldsymbol{0}
&
\boldsymbol{d}_a^{(2)}
&
\boldsymbol{D}_B^{(2)}
\end{bmatrix}.
} D = [ D A ( 1 ) 0 d b ( 1 ) d a ( 2 ) 0 D B ( 2 ) ] . The central column corresponds to the shared interface state.
General K K K ¶ The global matrix has the schematic form
D = [ D ( 1 ) D ( 2 ) ⋱ D ( K ) ] , \boldsymbol{D}
=
\begin{bmatrix}
\boxed{\boldsymbol{D}^{(1)}} & & & \\
& \boxed{\boldsymbol{D}^{(2)}} & & \\
& & \ddots & \\
& & & \boxed{\boldsymbol{D}^{(K)}}
\end{bmatrix}, D = ⎣ ⎡ D ( 1 ) D ( 2 ) ⋱ D ( K ) ⎦ ⎤ , where neighboring boxes overlap by one state column.
This structure yields a sparse matrix even when the total number of collocation points is large.
Global Collocation Equations ¶ The local collocation equations from all mesh intervals combine into
D X = t f − t 0 2 F . \boxed{
\boldsymbol{D}\boldsymbol{X}
=
\frac{t_f-t_0}{2}\boldsymbol{F}.
} D X = 2 t f − t 0 F . Equivalently, define the global defect matrix:
Δ = D X − t f − t 0 2 F . \boxed{
\boldsymbol{\Delta}
=
\boldsymbol{D}\boldsymbol{X}
-
\frac{t_f-t_0}{2}\boldsymbol{F}.
} Δ = D X − 2 t f − t 0 F . The NLP enforces
Δ = 0 . \boldsymbol{\Delta}=\boldsymbol{0}. Δ = 0 . The dimensions are
D ∈ R N × ( N + 1 ) , X ∈ R ( N + 1 ) × n x , F ∈ R N × n x , Δ ∈ R N × n x . \begin{aligned}
\boldsymbol{D}&\in\mathbb{R}^{N\times(N+1)},\\
\boldsymbol{X}&\in\mathbb{R}^{(N+1)\times n_x},\\
\boldsymbol{F}&\in\mathbb{R}^{N\times n_x},\\
\boldsymbol{\Delta}&\in\mathbb{R}^{N\times n_x}.
\end{aligned} D X F Δ ∈ R N × ( N + 1 ) , ∈ R ( N + 1 ) × n x , ∈ R N × n x , ∈ R N × n x . The multiplication is therefore conformable.
Composite LGR Quadrature ¶ Let the transformed running cost at each collocation point be
L i ( k ) = L ( X i ( k ) , U i ( k ) , t i ( k ) , p ) . L_i^{(k)}
=
L\left(
\boldsymbol{X}_i^{(k)},
\boldsymbol{U}_i^{(k)},
t_i^{(k)},
\boldsymbol{p}
\right). L i ( k ) = L ( X i ( k ) , U i ( k ) , t i ( k ) , p ) . Stack the values:
L = [ L 1 ( 1 ) ⋮ L N 1 ( 1 ) ⋮ L 1 ( K ) ⋮ L N K ( K ) ] ∈ R N . \boldsymbol{L}
=
\begin{bmatrix}
L_1^{(1)}\\
\vdots\\
L_{N_1}^{(1)}\\
\vdots\\
L_1^{(K)}\\
\vdots\\
L_{N_K}^{(K)}
\end{bmatrix}
\in\mathbb{R}^N. L = ⎣ ⎡ L 1 ( 1 ) ⋮ L N 1 ( 1 ) ⋮ L 1 ( K ) ⋮ L N K ( K ) ⎦ ⎤ ∈ R N . The integral is approximated as
∫ − 1 1 L ( s ) d s ≈ ∑ k = 1 K ∑ i = 1 N k w i ( k ) L i ( k ) . \int_{-1}^{1}L(s)\,\mathrm{d} s
\approx
\sum_{k=1}^{K}
\sum_{i=1}^{N_k}
w_i^{(k)}L_i^{(k)}. ∫ − 1 1 L ( s ) d s ≈ k = 1 ∑ K i = 1 ∑ N k w i ( k ) L i ( k ) . Using the stacked vectors:
∫ − 1 1 L ( s ) d s ≈ W T L . \boxed{
\int_{-1}^{1}L(s)\,\mathrm{d} s
\approx
\boldsymbol{W}^\mathsf{T}\boldsymbol{L}.
} ∫ − 1 1 L ( s ) d s ≈ W T L . Discrete Cost Functional ¶ The Bolza objective is
J = Φ ( x ( t 0 ) , t 0 , x ( t f ) , t f , p ) + ∫ t 0 t f L ( x , u , t , p ) d t . \begin{aligned}
J
&=
\Phi\left(
\boldsymbol{x}(t_0),t_0,
\boldsymbol{x}(t_f),t_f,\boldsymbol{p}
\right)
+
\int_{t_0}^{t_f}
L\left(\boldsymbol{x},\boldsymbol{u},t,\boldsymbol{p}\right)\,\mathrm{d} t.
\end{aligned} J = Φ ( x ( t 0 ) , t 0 , x ( t f ) , t f , p ) + ∫ t 0 t f L ( x , u , t , p ) d t . The multiple-interval LGR approximation is
J h = Φ ( X 1 ( 1 ) , t 0 , X N K + 1 ( K ) , t f , p ) + t f − t 0 2 W T L . \boxed{
J_h
=
\Phi\left(
\boldsymbol{X}_1^{(1)},t_0,
\boldsymbol{X}_{N_K+1}^{(K)},t_f,\boldsymbol{p}
\right)
+
\frac{t_f-t_0}{2}
\boldsymbol{W}^\mathsf{T}\boldsymbol{L}.
} J h = Φ ( X 1 ( 1 ) , t 0 , X N K + 1 ( K ) , t f , p ) + 2 t f − t 0 W T L . This compact expression has exactly the same form as the single-interval LGR objective.
Boundary Constraints ¶ The discrete endpoint constraints are
ϕ ( X 1 ( 1 ) , t 0 , X N K + 1 ( K ) , t f , p ) = 0 . \boxed{
\boldsymbol{\phi}\left(
\boldsymbol{X}_1^{(1)},t_0,
\boldsymbol{X}_{N_K+1}^{(K)},t_f,
\boldsymbol{p}
\right)
=
\boldsymbol{0}.
} ϕ ( X 1 ( 1 ) , t 0 , X N K + 1 ( K ) , t f , p ) = 0 . Only the first state and final state enter the boundary constraints.
More general bounds may be written as
ϕ min ≤ ϕ ( X 1 ( 1 ) , t 0 , X N K + 1 ( K ) , t f , p ) ≤ ϕ max . \boldsymbol{\phi}_{\min}
\le
\boldsymbol{\phi}\left(
\boldsymbol{X}_1^{(1)},t_0,
\boldsymbol{X}_{N_K+1}^{(K)},t_f,\boldsymbol{p}
\right)
\le
\boldsymbol{\phi}_{\max}. ϕ m i n ≤ ϕ ( X 1 ( 1 ) , t 0 , X N K + 1 ( K ) , t f , p ) ≤ ϕ m a x . Path Constraints ¶ Path constraints are evaluated at all collocation points:
c min ≤ c ( X ~ , U , S ~ , t 0 , t f , p ) ≤ c max . \boxed{
\boldsymbol{c}_{\min}
\le
\boldsymbol{c}\left(
\widetilde{\boldsymbol{X}},
\boldsymbol{U},
\widetilde{\boldsymbol{S}},
t_0,t_f,\boldsymbol{p}
\right)
\le
\boldsymbol{c}_{\max}.
} c m i n ≤ c ( X , U , S , t 0 , t f , p ) ≤ c m a x . Here:
X ~ \widetilde{\boldsymbol{X}} X
U \boldsymbol{U} U
S ~ \widetilde{\boldsymbol{S}} S
physical time is recovered from S ~ \widetilde{\boldsymbol{S}} S
The final state point is excluded because no control or collocation equation is defined there.
Physical-Time Vector ¶ The global collocation-time vector is
t ~ = t f − t 0 2 S ~ + t f + t 0 2 1 N . \widetilde{\boldsymbol{t}}
=
\frac{t_f-t_0}{2}\widetilde{\boldsymbol{S}}
+
\frac{t_f+t_0}{2}\boldsymbol{1}_N. t = 2 t f − t 0 S + 2 t f + t 0 1 N . The full state-time vector is
t = t f − t 0 2 S + t f + t 0 2 1 N + 1 . \boldsymbol{t}
=
\frac{t_f-t_0}{2}\boldsymbol{S}
+
\frac{t_f+t_0}{2}\boldsymbol{1}_{N+1}. t = 2 t f − t 0 S + 2 t f + t 0 1 N + 1 . If t 0 t_0 t 0 t f t_f t f
Final Nonlinear Programming Problem ¶ The complete multiple-interval LGR transcription is
min X , U , t 0 , t f , p J h \boxed{
\min_{\boldsymbol{X},\boldsymbol{U},t_0,t_f,\boldsymbol{p}}
J_h
} X , U , t 0 , t f , p min J h subject to
D X − t f − t 0 2 F = 0 , ϕ ( X 1 , t 0 , X N + 1 , t f , p ) = 0 , c min ≤ c ( X ~ , U , t ~ , p ) ≤ c max , Z min ≤ Z ≤ Z max . \begin{aligned}
\boldsymbol{D}\boldsymbol{X}
-
\frac{t_f-t_0}{2}\boldsymbol{F}
&=
\boldsymbol{0},
\\
\boldsymbol{\phi}\left(
\boldsymbol{X}_1,t_0,
\boldsymbol{X}_{N+1},t_f,\boldsymbol{p}
\right)
&=
\boldsymbol{0},
\\
\boldsymbol{c}_{\min}
\le
\boldsymbol{c}\left(
\widetilde{\boldsymbol{X}},
\boldsymbol{U},
\widetilde{\boldsymbol{t}},
\boldsymbol{p}
\right)
&\le
\boldsymbol{c}_{\max},
\\
\boldsymbol{Z}_{\min}
\le
\boldsymbol{Z}
&\le
\boldsymbol{Z}_{\max}.
\end{aligned} D X − 2 t f − t 0 F ϕ ( X 1 , t 0 , X N + 1 , t f , p ) c m i n ≤ c ( X , U , t , p ) Z m i n ≤ Z = 0 , = 0 , ≤ c m a x , ≤ Z m a x . A convenient decision vector is
Z = [ vec ( X ) vec ( U ) t 0 t f p ] . \boldsymbol{Z}
=
\begin{bmatrix}
\operatorname{vec}(\boldsymbol{X})\\
\operatorname{vec}(\boldsymbol{U})\\
t_0\\
t_f\\
\boldsymbol{p}
\end{bmatrix}. Z = ⎣ ⎡ vec ( X ) vec ( U ) t 0 t f p ⎦ ⎤ . The single-interval LGR collocation equation is
D s i n g l e X = t f − t 0 2 F . \boldsymbol{D}_{\mathrm{single}}\boldsymbol{X}
=
\frac{t_f-t_0}{2}\boldsymbol{F}. D single X = 2 t f − t 0 F . The multiple-interval equation is
D m u l t i X = t f − t 0 2 F . \boldsymbol{D}_{\mathrm{multi}}\boldsymbol{X}
=
\frac{t_f-t_0}{2}\boldsymbol{F}. D multi X = 2 t f − t 0 F . The algebraic form is identical. The difference lies in the structure of D \boldsymbol{D} D
Likewise, the objective retains the form
J h = Φ h + t f − t 0 2 W T L . J_h
=
\Phi_h
+
\frac{t_f-t_0}{2}
\boldsymbol{W}^\mathsf{T}\boldsymbol{L}. J h = Φ h + 2 t f − t 0 W T L . This structural equivalence simplifies implementation.
Jacobian Sparsity ¶ The global defect constraints are
Δ = D X − t f − t 0 2 F ( X ~ , U , t ~ , p ) . \boldsymbol{\Delta}
=
\boldsymbol{D}\boldsymbol{X}
-
\frac{t_f-t_0}{2}
\boldsymbol{F}(\widetilde{\boldsymbol{X}},\boldsymbol{U},\widetilde{\boldsymbol{t}},\boldsymbol{p}). Δ = D X − 2 t f − t 0 F ( X , U , t , p ) . A defect row in interval k k k
Therefore, the Jacobian with respect to state and control variables is block sparse.
The principal Jacobian terms are
∂ Δ ∂ X = D − t f − t 0 2 ∂ F ∂ X , \frac{\partial\boldsymbol{\Delta}}{\partial\boldsymbol{X}}
=
\boldsymbol{D}
-
\frac{t_f-t_0}{2}
\frac{\partial\boldsymbol{F}}{\partial\boldsymbol{X}}, ∂ X ∂ Δ = D − 2 t f − t 0 ∂ X ∂ F , and
∂ Δ ∂ U = − t f − t 0 2 ∂ F ∂ U . \frac{\partial\boldsymbol{\Delta}}{\partial\boldsymbol{U}}
=
-
\frac{t_f-t_0}{2}
\frac{\partial\boldsymbol{F}}{\partial\boldsymbol{U}}. ∂ U ∂ Δ = − 2 t f − t 0 ∂ U ∂ F . Number of Variables and Constraints ¶ Assume:
n x n_x n x
n u n_u n u
n p n_p n p
N N N
The state variables contribute
( N + 1 ) n x . (N+1)n_x. ( N + 1 ) n x . The control variables contribute
Including time and parameters:
n z = ( N + 1 ) n x + N n u + 2 + n p . \boxed{
n_z
=
(N+1)n_x
+
Nn_u
+
2
+
n_p.
} n z = ( N + 1 ) n x + N n u + 2 + n p . The defect equations contribute
n d e f e c t = N n x . \boxed{
n_{\mathrm{defect}}
=
Nn_x.
} n defect = N n x . Additional constraints arise from endpoint and path conditions.
Implementation Strategy ¶ A robust implementation separates mesh construction from NLP evaluation.
Preprocessing ¶ Before solving:
choose mesh points s k s_k s k
choose degrees N k N_k N k
compute standard LGR nodes and weights;
map nodes into each interval;
scale interval quadrature weights;
scale interval differentiation matrices;
assemble S ~ \widetilde{\boldsymbol{S}} S S \boldsymbol{S} S W \boldsymbol{W} W D \boldsymbol{D} D
create state and control index maps.
NLP Evaluation ¶ At each NLP iteration:
unpack X \boldsymbol{X} X U \boldsymbol{U} U t 0 t_0 t 0 t f t_f t f p \boldsymbol{p} p
compute the physical-time vector;
evaluate F \boldsymbol{F} F L \boldsymbol{L} L
evaluate the cost;
evaluate defects;
evaluate endpoint constraints;
evaluate path constraints;
evaluate derivatives.
MATLAB Assembly of the Global Matrix ¶ function D = assembleGlobalLGR(Dlocal, degree)
K = numel(degree);
N = sum(degree);
D = sparse(N,N+1);
rowStart = 1;
colStart = 1;
for k = 1:K
Nk = degree(k);
rows = rowStart:(rowStart+Nk-1);
cols = colStart:(colStart+Nk);
D(rows,cols) = Dlocal{k};
rowStart = rowStart + Nk;
% Advance by Nk, not Nk+1, because the last
% state column is shared with the next interval.
colStart = colStart + Nk;
end
endMATLAB Stacking of Nodes and Weights ¶ function [Scol,Sstate,W] = stackLGRData(nodes,weights,mesh)
K = numel(nodes);
Scol = [];
W = [];
for k = 1:K
Scol = [Scol; nodes{k}(:)];
W = [W; weights{k}(:)];
end
Sstate = [Scol; mesh(end)];
endMATLAB Defect Function ¶ function defect = lgrDefect(X,U,t0,tf,Scol,D,p)
N = size(U,1);
t = 0.5*(tf-t0)*Scol ...
+ 0.5*(tf+t0)*ones(N,1);
F = dynamics(X(1:N,:),U,t,p);
defect = D*X - 0.5*(tf-t0)*F;
endMATLAB Objective Function ¶ function J = lgrObjective(X,U,t0,tf,Scol,W,p)
N = size(U,1);
t = 0.5*(tf-t0)*Scol ...
+ 0.5*(tf+t0)*ones(N,1);
L = runningCost(X(1:N,:),U,t,p);
M = endpointCost(X(1,:),t0,X(end,:),tf,p);
J = M + 0.5*(tf-t0)*(W.'*L);
endIndexing Example ¶ Suppose
K = 3 , N 1 = 3 , N 2 = 4 , N 3 = 2. K=3,
\qquad
N_1=3,
\qquad
N_2=4,
\qquad
N_3=2. K = 3 , N 1 = 3 , N 2 = 4 , N 3 = 2. Then
N = 3 + 4 + 2 = 9. N=3+4+2=9. N = 3 + 4 + 2 = 9. The global state rows are:
1 : 4 for interval 1 , 4 : 8 for interval 2 , 8 : 10 for interval 3 . \begin{aligned}
1\!:\!4
&\quad\text{for interval 1},\\
4\!:\!8
&\quad\text{for interval 2},\\
8\!:\!10
&\quad\text{for interval 3}.
\end{aligned} 1 : 4 4 : 8 8 : 10 for interval 1 , for interval 2 , for interval 3 . Rows 4 and 8 are shared interface states.
The global differentiation matrix has size
The controls have nine rows. The weight vector has nine entries.
Integral State Consistency ¶ The local state update can be checked using LGR quadrature:
X N k + 1 ( k ) − X 1 ( k ) ≈ t f − t 0 2 ∑ i = 1 N k w i ( k ) F i ( k ) . \boldsymbol{X}_{N_k+1}^{(k)}
-
\boldsymbol{X}_1^{(k)}
\approx
\frac{t_f-t_0}{2}
\sum_{i=1}^{N_k}
w_i^{(k)}
\boldsymbol{F}_i^{(k)}. X N k + 1 ( k ) − X 1 ( k ) ≈ 2 t f − t 0 i = 1 ∑ N k w i ( k ) F i ( k ) . For the full horizon:
X N + 1 − X 1 ≈ t f − t 0 2 W T F . \boldsymbol{X}_{N+1}
-
\boldsymbol{X}_1
\approx
\frac{t_f-t_0}{2}
\boldsymbol{W}^\mathsf{T}\boldsymbol{F}. X N + 1 − X 1 ≈ 2 t f − t 0 W T F . For vector-valued dynamics, the weighted sum is applied columnwise.
This identity is useful for debugging quadrature and indexing.
Why the Final Point Has No Control ¶ The LGR state polynomial requires N k + 1 N_k+1 N k + 1 N k N_k N k
Therefore:
state nodes = N + 1 , control nodes = N . \text{state nodes}=N+1,
\qquad
\text{control nodes}=N. state nodes = N + 1 , control nodes = N . A terminal control can be estimated after solving by extrapolation if needed, but it is not an independent NLP variable.
Path-Constraint Verification between Nodes ¶ Enforcing a path constraint only at collocation points does not guarantee that it is satisfied between nodes.
After solving, evaluate the state and control polynomials on a dense grid:
s c h e c k , m ∈ [ − 1 , 1 ] . s_{\mathrm{check},m}\in[-1,1]. s check , m ∈ [ − 1 , 1 ] . Then verify
c min ≤ c ( x i n t e r p , u i n t e r p , t i n t e r p ) ≤ c max . \boldsymbol{c}_{\min}
\le
\boldsymbol{c}\left(
\boldsymbol{x}_{\mathrm{interp}},
\boldsymbol{u}_{\mathrm{interp}},
t_{\mathrm{interp}}
\right)
\le
\boldsymbol{c}_{\max}. c m i n ≤ c ( x interp , u interp , t interp ) ≤ c m a x . Violations may require:
Motivation for Mesh Refinement ¶ The derivation above assumes that the mesh points and polynomial degrees are given. In practice, the important questions are:
Where should mesh points be placed?
How many intervals are required?
What polynomial degree should be used in each interval?
How should discontinuities be detected?
When should an interval be subdivided?
When should the polynomial degree be increased?
These questions motivate adaptive mesh-refinement algorithms.
h h h p p p h p hp h p ¶ h h h ¶ Decrease interval width by subdividing:
→ [ s k − 1 , s m ] ∪ [ s m , s k ] . \rightarrow
[s_{k-1},s_m]\cup[s_m,s_k]. → [ s k − 1 , s m ] ∪ [ s m , s k ] . This is suitable near discontinuities or sharp local features.
p p p ¶ Increase N k N_k N k
This is suitable when the solution is smooth.
h p hp h p ¶ Adapt both interval width and polynomial degree.
A practical strategy is:
use p p p
use h h h
use both when the appropriate classification is uncertain.
Error Estimation Concept ¶ Let X c o l l ( k ) ( s ) \boldsymbol{X}_{\mathrm{coll}}^{(k)}(s) X coll ( k ) ( s )
X i n t ( k ) ( s ) = X 1 ( k ) + t f − t 0 2 ∫ s k − 1 s f ( X c o l l ( k ) ( σ ) , U ( k ) ( σ ) , t ( σ ) ) d σ . \boldsymbol{X}_{\mathrm{int}}^{(k)}(s)
=
\boldsymbol{X}_1^{(k)}
+
\frac{t_f-t_0}{2}
\int_{s_{k-1}}^s
\boldsymbol{f}\left(
\boldsymbol{X}_{\mathrm{coll}}^{(k)}(\sigma),
\boldsymbol{U}^{(k)}(\sigma),
t(\sigma)
\right)
\,\mathrm{d}\sigma. X int ( k ) ( s ) = X 1 ( k ) + 2 t f − t 0 ∫ s k − 1 s f ( X coll ( k ) ( σ ) , U ( k ) ( σ ) , t ( σ ) ) d σ . Define
e k = max s ∈ S k ∥ X i n t ( k ) ( s ) − X c o l l ( k ) ( s ) ∥ ∞ 1 + ∥ X c o l l ( k ) ( s ) ∥ ∞ . e_k
=
\max_{s\in\mathcal{S}_k}
\frac{
\left\|
\boldsymbol{X}_{\mathrm{int}}^{(k)}(s)
-
\boldsymbol{X}_{\mathrm{coll}}^{(k)}(s)
\right\|_\infty
}{
1+\left\|\boldsymbol{X}_{\mathrm{coll}}^{(k)}(s)\right\|_\infty
}. e k = s ∈ S k max 1 + ∥ ∥ X coll ( k ) ( s ) ∥ ∥ ∞ ∥ ∥ X int ( k ) ( s ) − X coll ( k ) ( s ) ∥ ∥ ∞ . If
refine interval k k k
Warm Starting after Mesh Refinement ¶ After changing the mesh:
interpolate the previous state solution onto the new state nodes;
interpolate the previous control solution onto the new collocation nodes;
preserve the previous t 0 t_0 t 0 t f t_f t f
project variables onto bounds if necessary;
solve the refined NLP.
Warm starting greatly improves efficiency.
Numerical Conditioning ¶ Potential sources of poor conditioning include:
very high polynomial degree;
extremely small mesh intervals;
badly scaled states;
inconsistent derivative scaling;
large differences among state magnitudes;
nearly singular endpoint constraints.
A practical implementation should scale variables and avoid excessive polynomial degrees within a single interval.
Verification Checklist ¶ A multiple-interval LGR solution should be checked using:
defect residuals;
endpoint residuals;
path constraints at and between nodes;
mesh-refinement convergence;
objective convergence;
independent forward integration;
state continuity at interfaces;
control behavior near switches;
consistency of matrix dimensions;
sensitivity to the initial mesh.
Common Implementation Errors ¶ duplicating interface states;
failing to overlap one state column between adjacent D ( k ) \boldsymbol{D}^{(k)} D ( k )
assigning a control variable to the final state point;
attaching a quadrature weight to the final point;
using X \boldsymbol{X} X X ~ \widetilde{\boldsymbol{X}} X
forgetting the global factor ( t f − t 0 ) / 2 (t_f-t_0)/2 ( t f − t 0 ) /2
double-scaling local differentiation matrices;
using inconsistent node order;
mapping local nodes incorrectly;
constructing an N × N N\times N N × N N × ( N + 1 ) N\times(N+1) N × ( N + 1 )
imposing unnecessary control continuity;
failing to check constraints between nodes.
Summary ¶ The main results of this section are:
The total number of LGR collocation points is
N = ∑ k = 1 K N k . N=\sum_{k=1}^{K}N_k. N = k = 1 ∑ K N k . There are N + 1 N+1 N + 1 N N N
Adjacent intervals share one state variable at each mesh point.
The forward LGR control has one value at each interior mesh point.
The global weight vector is obtained by stacking local weights.
The global differentiation matrix has size
N × ( N + 1 ) . N\times(N+1). N × ( N + 1 ) . Adjacent local differentiation blocks overlap by one state column.
The global defect equation is
D X − t f − t 0 2 F = 0 . \boldsymbol{D}\boldsymbol{X}
-
\frac{t_f-t_0}{2}\boldsymbol{F}
=
\boldsymbol{0}. D X − 2 t f − t 0 F = 0 . The discrete objective is
J h = Φ h + t f − t 0 2 W T L . J_h
=
\Phi_h
+
\frac{t_f-t_0}{2}\boldsymbol{W}^\mathsf{T}\boldsymbol{L}. J h = Φ h + 2 t f − t 0 W T L . The multiple-interval and single-interval methods have the same compact algebraic form.
Their principal difference is the sparse structure of the differentiation matrix.
Mesh refinement determines where intervals and polynomial degrees should be placed.
Connection. Once the global method is available, interpolation and integration error estimates can drive an adaptive choice between mesh subdivision and polynomial enrichment.