Direct Collocation with Euler Discretization

The purpose of this section is to move from an abstract description of direct methods to a concrete nonlinear programming formulation. The focus is not on software syntax, but on the mathematical structure that must be created before an NLP solver can be used.

The main topics are:

Euler forward in collocation form;
simultaneous solution of discrete dynamics;
root finding for differential-equation residuals;
transcription of a nonlinear optimal-control problem;
discretization of time, state, control, and cost;
treatment of boundary conditions;
counting variables, constraints, and degrees of freedom; and
interpretation of the resulting nonlinear program.

Euler Forward: Time-Marching and Collocation Forms¶

For the ordinary differential equation

\dot{\boldsymbol{x}}(t)=\boldsymbol{f}(t,\boldsymbol{x}(t)),

(1)

Euler forward gives

\boldsymbol{x}_{k+1} = \boldsymbol{x}_k + h_k\boldsymbol{f}(t_k,\boldsymbol{x}_k).

(2)

This is the time-marching form.

The same equation can be written as a residual:

\boxed{ \boldsymbol{x}_{k+1} - \boldsymbol{x}_k - h_k\boldsymbol{f}(t_k,\boldsymbol{x}_k) = \boldsymbol{0}. }

(3)

This is the collocation or defect form.

Mesh and Discrete State Vector¶

Let

t_0<t_1<\cdots<t_K=t_f,

(4)

with

h_k=t_{k+1}-t_k.

(5)

Define

\boldsymbol{x}_k\approx\boldsymbol{x}(t_k).

(6)

Collect all state approximations into

\boldsymbol{y} = \begin{bmatrix} \boldsymbol{x}_0\\ \boldsymbol{x}_1\\ \vdots\\ \boldsymbol{x}_K \end{bmatrix}.

(7)

The Euler residuals form a nonlinear system:

\boldsymbol{G}(\boldsymbol{y})=\boldsymbol{0}.

(8)

Initial Conditions as Algebraic Constraints¶

Suppose

\boldsymbol{x}(t_0)=\boldsymbol{a},

(9)

where $\boldsymbol{a}$ is known.

One may either substitute

\boldsymbol{x}_0=\boldsymbol{a}

(10)

directly or add the residual equation

\boxed{ \boldsymbol{x}_0-\boldsymbol{a}=\boldsymbol{0}. }

(11)

The second form is common in collocation because it keeps all conditions in one residual system.

Newton’s Method for the Global Residual System¶

To solve

\boldsymbol{G}(\boldsymbol{y})=\boldsymbol{0},

(12)

consider

\boldsymbol{G}(\boldsymbol{y}+\Delta\boldsymbol{y})\approx \boldsymbol{G}(\boldsymbol{y}) + \frac{\partial\boldsymbol{G}}{\partial\boldsymbol{y}}\Delta\boldsymbol{y}.

(13)

Setting this approximation equal to zero gives

\frac{\partial\boldsymbol{G}}{\partial\boldsymbol{y}}\Delta\boldsymbol{y} = -\boldsymbol{G}(\boldsymbol{y}).

(14)

Hence,

\Delta\boldsymbol{y} = - \left( \frac{\partial\boldsymbol{G}}{\partial\boldsymbol{y}} \right)^{-1} \boldsymbol{G}(\boldsymbol{y}).

(15)

The Newton iteration is

\boxed{ \boldsymbol{y}^{(j+1)} = \boldsymbol{y}^{(j)} - \left[ \frac{\partial\boldsymbol{G}}{\partial\boldsymbol{y}} \left(\boldsymbol{y}^{(j)}\right) \right]^{-1} \boldsymbol{G}\left(\boldsymbol{y}^{(j)}\right). }

(16)

In implementation, one solves the linear system rather than explicitly forming the inverse.

A Nonlinear Optimal-Control Example¶

Consider

\boxed{ \min_{x(\cdot),u(\cdot)} J = \frac12 \int_0^{t_f} \left( x^2(t)+u^2(t) \right)\,\mathrm{d} t }

(17)

subject to

\dot{x}(t) = -x^3(t)+u(t),

(18)

with

x(0)=a, \qquad x(t_f)=b,

(19)

and fixed $t_f$ .

This problem is nonlinear because of the term $-x^3$ .

Uniform Time Mesh¶

Choose $K$ intervals:

h = \frac{t_f-t_0}{K}.

(20)

The mesh points are

t_k=t_0+kh, \qquad k=0,\ldots,K.

(21)

There are:

$K$ intervals;
$K+1$ mesh points.

State and Control Discretization¶

For Euler forward, the state is defined at all mesh points:

x_0,x_1,\ldots,x_K.

(22)

The control is required only at the beginning of each interval:

u_0,u_1,\ldots,u_{K-1}.

(23)

There is no $u_K$ in this transcription because the final node does not begin another interval.

Discrete Dynamics¶

Euler forward gives

x_{k+1} = x_k + h\left(-x_k^3+u_k\right).

(24)

The defect equation is

\boxed{ \Delta_k = x_{k+1} - x_k - h\left(-x_k^3+u_k\right) = 0, }

(25)

for

k=0,\ldots,K-1.

(26)

Equivalently,

\Delta_k = x_{k+1} - x_k + hx_k^3 - hu_k = 0.

(27)

Discretization of the Cost¶

The running cost is

L(x,u)=\frac12(x^2+u^2).

(28)

Using a left-rectangle rule consistent with Euler forward,

\int_0^{t_f}(x^2+u^2)\,\mathrm{d} t \approx h\sum_{k=0}^{K-1}(x_k^2+u_k^2).

(29)

Therefore,

\boxed{ J_a = \frac{h}{2} \sum_{k=0}^{K-1} \left( x_k^2+u_k^2 \right). }

(30)

The original functional has become an ordinary function of finitely many variables.

Boundary Conditions¶

The continuous boundary conditions map to

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

(31)

These are equality constraints in the NLP.

Decision Vector¶

Define the decision vector

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

(32)

The total number of variables is

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

(33)

Complete NLP Formulation¶

The Euler direct-collocation problem is

\boxed{ \begin{aligned} \min_{\boldsymbol{z}}\quad & \frac{h}{2} \sum_{k=0}^{K-1} \left( x_k^2+u_k^2 \right) \\[1mm] \text{subject to}\quad & x_{k+1}-x_k-h(-x_k^3+u_k)=0, \\ & k=0,\ldots,K-1, \\ & x_0-a=0, \\ & x_K-b=0. \end{aligned} }

(34)

This is a nonlinear programming problem.

Variables, Constraints, and Degrees of Freedom¶

The number of variables is

n_z=2K+1.

(35)

The number of equality constraints is

n_c=K+2.

(36)

Therefore,

\boxed{ n_{\mathrm{dof}} = n_z-n_c = K-1. }

(37)

The problem must have degrees of freedom because the optimizer must choose among many feasible trajectories.

Why the System Is Underdetermined Before Optimization¶

The defect equations and boundary conditions alone do not uniquely determine all state and control values.

There are many controls that may approximately transfer the state from $a$ to $b$ .

The objective selects the one that minimizes

\frac12\int_0^{t_f}(x^2+u^2)\,\mathrm{d} t.

(38)

Thus:

the constraints define feasibility;
the objective selects optimality.

Generalization to Multiple States and Controls¶

For

\dot{\boldsymbol{x}} = \boldsymbol{f}(t,\boldsymbol{x},\boldsymbol{u}),

(39)

with

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

(40)

Euler forward gives

\boldsymbol{x}_{k+1} - \boldsymbol{x}_k - h_k\boldsymbol{f}(t_k,\boldsymbol{x}_k,\boldsymbol{u}_k) = \boldsymbol{0}.

(41)

The number of variables becomes

(K+1)n_x+Kn_u.

(42)

The number of dynamic constraints becomes

Kn_x.

(43)

Endpoint and path constraints are added separately.

Forward Versus Backward Euler Indexing¶

Euler forward uses

\boldsymbol{f}(t_k,\boldsymbol{x}_k,\boldsymbol{u}_k).

(44)

Backward Euler uses

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

(45)

Hence:

forward Euler naturally uses control values $u_0,\ldots,u_{K-1}$ ;
backward Euler naturally uses control values $u_1,\ldots,u_K$ .

The quadrature indexing should be consistent with the dynamic discretization.

Accuracy of the Transcription¶

Euler forward is first-order accurate.

Therefore, the discrete optimal-control solution is expected to converge at first order under appropriate regularity and stability assumptions.

A valid numerical study should compare solutions for:

K,\quad 2K,\quad 4K,\ldots

(46)

and monitor:

objective convergence;
state convergence;
control convergence;
defect residuals;
endpoint residuals.

Implementation Workflow¶

Select $K$ .
Compute $h$ .
Create an initial guess for all states and controls.
Build the objective function.
Build the Euler defect constraints.
Add endpoint constraints.
Supply bounds if needed.
Call an NLP solver.
Reconstruct the trajectories.
Refine the mesh and repeat.

Common Errors¶

Using $K$ state points instead of $K+1$ .
Including $u_K$ in an Euler-forward transcription without a reason.
Summing the quadrature from 0 to $K$ instead of 0 to $K-1$ .
Forgetting the endpoint constraints.
Counting variables or constraints incorrectly.
Solving the defect equations without including the objective.
Calling the discrete solution exact.
Failing to refine the mesh.

Summary¶

Euler forward can be written as collocation defects.
The state and control values become NLP decision variables.
The running cost becomes a finite quadrature sum.
Endpoint conditions become algebraic constraints.
The resulting finite-dimensional problem is an NLP.
The example contains $2K+1$ variables and $K+2$ equality constraints.
The remaining $K-1$ degrees of freedom are optimized through the cost.
Mesh refinement is required to assess accuracy.

Connection. Once collocation has produced an NLP, its usefulness depends on exploiting local dependence, sparse derivatives, and a fixed normalized time interval.