Sparse NLP Structure and Time Normalization

Direct collocation produces nonlinear programs that are often large but sparse. This structure is one of the main reasons collocation methods are computationally practical.

This section develops:

local dependence of collocation constraints;
constraint-Jacobian sparsity;
banded structure;
nonzero counting;
sparsity scaling with mesh refinement;
implications for NLP solvers;
time normalization to a fixed interval;
transformation of dynamics and cost;
free initial and final times; and
extensions to large and multi-state systems.

Review of the Euler NLP¶

Consider the decision vector

\boldsymbol{z} = \begin{bmatrix} x_0 & x_1 & \cdots & x_K & u_0 & u_1 & \cdots & u_{K-1} \end{bmatrix}^{\mathsf{T}}.

(1)

The defect equations are

c_k = x_{k+1} - x_k - h(-x_k^3+u_k) = 0,

(2)

for

k=0,\ldots,K-1.

(3)

The endpoint constraints are

\begin{aligned} b_1&=x_0-a=0,\\ b_2&=x_K-b=0. \end{aligned}

(4)

Local Dependence¶

Each defect $c_k$ depends only on:

x_k,\qquad x_{k+1},\qquad u_k.

(5)

It does not depend on distant variables such as:

x_{k+4},\quad u_{k+6},\quad \text{etc.}

(6)

This local dependence creates sparsity.

Constraint Jacobian¶

Define

\boldsymbol{g}(\boldsymbol{z}) = \begin{bmatrix} c_0\\ c_1\\ \vdots\\ c_{K-1}\\ b_1\\ b_2 \end{bmatrix}.

(7)

The constraint Jacobian is

\boxed{ \boldsymbol{J}_c = \frac{\partial\boldsymbol{g}}{\partial\boldsymbol{z}}. }

(8)

Its dimensions are

(K+2)\times(2K+1).

(9)

Derivatives of One Defect¶

For

c_k = x_{k+1}-x_k+hx_k^3-hu_k,

(10)

the nonzero derivatives are

\begin{aligned} \frac{\partial c_k}{\partial x_k} &= -1+3hx_k^2,\\ \frac{\partial c_k}{\partial x_{k+1}} &= 1,\\ \frac{\partial c_k}{\partial u_k} &= -h. \end{aligned}

(11)

All other partial derivatives of $c_k$ are zero.

Sparsity Pattern¶

The structure is schematically:

\begin{array}{c|cccccc|ccccc} & x_0&x_1&x_2&\cdots&x_{K-1}&x_K & u_0&u_1&u_2&\cdots&u_{K-1} \\ \hline c_0&\times&\times&&&&&\times&&&&\\ c_1&&\times&\times&&&&&\times&&&\\ c_2&&&\times&\times&&&&&\times&&\\ \vdots&&&&\ddots&\ddots&&&&\ddots&\\ c_{K-1}&&&&&\times&\times&&&&\times\\ \hline b_1&\times&&&&&&&&&&\\ b_2&&&&&&\times&&&&& \end{array}

(12)

Only a narrow band of entries is nonzero.

Counting Nonzero Entries¶

Each defect contributes three nonzeros.

There are $K$ defects, so:

3K

(13)

nonzeros arise from the dynamics.

The two endpoint constraints contribute two more.

Therefore,

\boxed{ \operatorname{nnz}(\boldsymbol{J}_c) = 3K+2. }

(14)

Jacobian Density¶

The total number of entries is

(K+2)(2K+1).

(15)

The density is

\rho_K = \frac{3K+2}{(K+2)(2K+1)}.

(16)

K\rightarrow\infty,

(17)

\rho_K \sim \frac{3}{2K} \rightarrow0.

(18)

Thus, the problem becomes increasingly sparse as the mesh is refined.

Why Sparse Structure Matters¶

A dense treatment stores and processes all entries.

A sparse treatment stores and processes only nonzero entries.

This affects:

memory usage;
Jacobian assembly;
linear-system solution;
factorization cost;
scalability;
parallelization.

For large trajectory-optimization problems, sparse linear algebra is essential.

Sparse NLP Solvers¶

Modern sparse NLP solvers exploit derivative structure.

Examples include:

IPOPT;
SNOPT;
KNITRO;
WORHP;
filterSQP-type solvers;
sparse interior-point and SQP methods.

MATLAB’s fmincon may be useful for small instructional problems, but large collocation problems benefit from solvers designed explicitly for sparse NLPs.

Why Collocation Naturally Produces Sparsity¶

A collocation equation on interval $k$ depends only on variables in or near that interval.

Typical local dependencies include:

(\boldsymbol{x}_k,\boldsymbol{x}_{k+1},\boldsymbol{u}_k,\boldsymbol{u}_{k+1}).

(19)

Therefore, distant nodes do not interact directly in the Jacobian.

This remains true for many methods:

Euler forward;
Euler backward;
trapezoidal collocation;
Hermite–Simpson;
Runge–Kutta collocation;
finite-element transcription;
multiple shooting.

Multiple States and Controls¶

For

\boldsymbol{x}\in\mathbb{R}^{n_x}, \qquad \boldsymbol{u}\in\mathbb{R}^{n_u},

(20)

each defect block has dimension $n_x$ .

A defect typically depends on:

$n_x$ components of $\boldsymbol{x}_k$ ;
$n_x$ components of $\boldsymbol{x}_{k+1}$ ;
$n_u$ components of $\boldsymbol{u}_k$ .

The Jacobian becomes block sparse and block banded.

Time Normalization¶

Instead of discretizing directly on

t\in[t_0,t_f],

(21)

it is convenient to map the domain to

\tau\in[-1,1].

(22)

Define

\boxed{ \tau = 2\frac{t-t_0}{t_f-t_0}-1. }

(23)

The inverse transformation is

\boxed{ t(\tau) = \frac{t_f-t_0}{2}\tau + \frac{t_f+t_0}{2}. }

(24)

Derivation of the Mapping¶

Starting from

t_0\leq t\leq t_f,

(25)

subtract $t_0$ :

0\leq t-t_0\leq t_f-t_0.

(26)

Divide by $t_f-t_0$ :

0 \leq \frac{t-t_0}{t_f-t_0} \leq 1.

(27)

Multiply by 2 and subtract 1:

-1 \leq 2\frac{t-t_0}{t_f-t_0}-1 \leq 1.

(28)

Transformation of the Dynamics¶

The original dynamics are

\frac{\,\mathrm{d}\boldsymbol{x}}{\,\mathrm{d} t} = \boldsymbol{f}(t,\boldsymbol{x},\boldsymbol{u}).

(29)

By the chain rule,

\frac{\,\mathrm{d}\boldsymbol{x}}{\,\mathrm{d} \tau} = \frac{\,\mathrm{d}\boldsymbol{x}}{\,\mathrm{d} t} \frac{\,\mathrm{d} t}{\,\mathrm{d}\tau}.

(30)

From the inverse map,

\frac{\,\mathrm{d} t}{\,\mathrm{d}\tau} = \frac{t_f-t_0}{2}.

(31)

Therefore,

\boxed{ \frac{\,\mathrm{d}\boldsymbol{x}}{\,\mathrm{d}\tau} = \frac{t_f-t_0}{2} \boldsymbol{f} \left( t(\tau),\boldsymbol{x}(\tau),\boldsymbol{u}(\tau) \right). }

(32)

Transformation of the Cost¶

For

J = \Phi + \int_{t_0}^{t_f} L(t,\boldsymbol{x},\boldsymbol{u})\,\mathrm{d} t,

(33)

use

\,\mathrm{d} t = \frac{t_f-t_0}{2}\,\mathrm{d}\tau.

(34)

Then

\boxed{ J = \Phi + \frac{t_f-t_0}{2} \int_{-1}^{1} L \left( t(\tau),\boldsymbol{x}(\tau),\boldsymbol{u}(\tau) \right) \,\mathrm{d}\tau. }

(35)

Free Initial and Final Times¶

If $t_0$ or $t_f$ is free, they become NLP decision variables.

The normalized computational domain remains fixed:

\tau\in[-1,1].

(36)

However, the duration

t_f-t_0

(37)

appears explicitly in:

the transformed dynamics;
the transformed cost;
time-dependent path constraints;
endpoint conditions.

Why Normalize Time¶

Time normalization provides:

a fixed computational domain;
simpler mesh generation;
easier handling of free final time;
direct compatibility with orthogonal polynomials;
a natural foundation for pseudospectral methods;
improved software modularity.

Normalized Euler Defects¶

Let

-1=\tau_0<\tau_1<\cdots<\tau_K=1.

(38)

Define

\Delta\tau_k = \tau_{k+1}-\tau_k.

(39)

Euler forward in normalized time gives

\boldsymbol{x}_{k+1} = \boldsymbol{x}_k + \Delta\tau_k \frac{t_f-t_0}{2} \boldsymbol{f}(t_k,\boldsymbol{x}_k,\boldsymbol{u}_k).

(40)

Thus,

\boxed{ \boldsymbol{x}_{k+1} - \boldsymbol{x}_k - \Delta\tau_k \frac{t_f-t_0}{2} \boldsymbol{f}(t_k,\boldsymbol{x}_k,\boldsymbol{u}_k) = \boldsymbol{0}. }

(41)

Normalized Quadrature¶

A left-rectangle approximation becomes

J_a = \Phi + \frac{t_f-t_0}{2} \sum_{k=0}^{K-1} \Delta\tau_k L(t_k,\boldsymbol{x}_k,\boldsymbol{u}_k).

(42)

For a uniform normalized mesh,

\Delta\tau = \frac{2}{K}.

(43)

Sparsity with Free Final Time¶

When $t_f$ is a decision variable, every defect depends on $t_f$ through the factor

\frac{t_f-t_0}{2}.

(44)

This introduces one additional dense column in the Jacobian.

However, the rest of the Jacobian remains sparse.

Thus, free-time problems retain strong sparse structure.

Finite Versus Infinite Horizons¶

The affine map above applies to finite intervals.

An infinite horizon,

t\in[t_0,\infty),

(45)

requires a nonlinear transformation, for example:

t = t_0 + \frac{1+\tau}{1-\tau}, \qquad \tau\in[-1,1).

(46)

Such mappings require special treatment because:

the derivative scaling becomes singular near the mapped endpoint;
quadrature weights may become extreme;
asymptotic behavior must be represented accurately;
convergence analysis is more delicate.

Solver-Oriented Implementation Workflow¶

Normalize time.
Select the mesh.
Define the decision vector.
Define the discretized objective.
Define defect constraints.
Define endpoint and path constraints.
Provide sparse Jacobian information.
Provide variable and constraint bounds.
Solve the NLP.
Verify convergence under mesh refinement.

Derivative Information¶

NLP solvers may use:

finite differences;
analytic derivatives;
automatic differentiation;
complex-step derivatives.

For large collocation problems, exact or automatic sparse derivatives are usually preferable.

Common Errors¶

Treating the Jacobian as dense.
Counting $2K-1$ variables instead of $2K+1$ in the scalar Euler example.
Forgetting endpoint-constraint nonzeros.
Assuming sparsity disappears for nonlinear dynamics.
Omitting the factor $(t_f-t_0)/2$ after time normalization.
Transforming the dynamics but not the cost.
Forgetting that $t(\tau)$ depends on both endpoint times.
Using a finite-interval affine map for an infinite-horizon problem.
Supplying inconsistent derivative sparsity to the solver.
Failing to verify mesh convergence.

Summary¶

Collocation constraints depend locally on neighboring variables.
This produces sparse, banded Jacobians.
The scalar Euler example has $3K+2$ Jacobian nonzeros.
Jacobian density approaches zero as the mesh is refined.
Sparse NLP solvers exploit this structure.
Time normalization maps arbitrary finite intervals to $[-1,1]$ .
The dynamics scale by $(t_f-t_0)/2$ .
The cost scales by the same factor.
Free endpoint times become NLP decision variables.
Normalization is foundational for advanced collocation and pseudospectral methods.