OCP with Path Constraint

Pontryagin’s Maximum Principle

Pontryagin’s Maximum Principle (PMP) is one of the most important parts in OCP solving. Based on variational method, PMP allows us to solve OCPs whose objective functional stagnation point cannot fall in admissible set, which makes PMP more widely used in OCPs with constraints, especially control constraints.

The dynamic eqautions is given:

\[\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u},t)\]

Objective function

\[J(\mathbf{x},\mathbf{u},t)=\varphi(\mathbf{x}(t_f),t_f)+\int_{t_0}^{t_f}L(\mathbf{x},\mathbf{u},t)dt\]

We have Hamilton function:

\[H(\mathbf{x},\mathbf{u},t)=L(\mathbf{x},\mathbf{u},t)+\mathbf{\lambda}^\intercal \mathbf{f}(\mathbf{x},\mathbf{u},t)\]

Therefore, if there’s no control constraint, we can simply let our objective function take the stationary value, therefore we have control equation:

\[\frac{\partial H}{\partial \mathbf{u}}=\mathbf{0}\\\]

Costate equation:

\[\dot{\mathbf{\lambda}}=-\frac{\partial H}{\partial \mathbf{x}}\]

All the equations come from variational method. For terminal conditions, you can also work them out using variational method.

However, if there is control constraint (e.g. the amplitude of thrust):

\[C(\mathbf{u},t)\leq 0\]

We have

\[H=L+\mathbf{\lambda}^\intercal \mathbf{f}\]

In case that the stationary point cannot be taken in the admissible control set $U$, we modified the control equation into the one below:

\[\mathbf{u^*}=\arg \min_{u\in U}H(\mathbf{u})\]

This is PMP. However, in Bryson’s Applied Optimal Control, this is solved with variational method instead of PMP. For terminal conditions in different situations, there are a lot of summary tables on technology websties like CSDN or StackOverFlow.

Path Constraint and Solutions

Equality Constraint

\[\begin{align} &S(\mathbf{x},\mathbf{u},t)=0\\ &S(\mathbf{x},t)=0 \end{align}\]

For state-and-control-mixed constraint $S(\mathbf{x},\mathbf{u},t)=0$, we have Hamilton function:

\[H=L+\mathbf{\lambda}^\intercal \mathbf{f}+\mu S\]

For pure-state constraint $S(\mathbf{x},t)=0$, we have to take the time derivative of $S$, until $S$ includes $\mathbf{u}$ explicitly:

\[S^{(q)}_{\mathbf{u}}(\mathbf{x},\mathbf{u},t)\neq0, S^{(q)}(\mathbf{x},\mathbf{u},t)=0\]

Boundary conditions:

\[\begin{bmatrix} S^{(0)}(\mathbf{x}(t_0),t_0)\\ S^{(1)}(\mathbf{x}(t_0),t_0)\\ \cdots\\ S^{(q-1)}(\mathbf{x}(t_0),t_0) \end{bmatrix} =\mathbf{0}\]

If there’s no control constraint, we have

\[\frac{\partial H}{\partial \mathbf{u}}=\mathbf{0}\\ \dot{\mathbf{\lambda}}=-\frac{\partial H}{\partial \mathbf{x}}\]

Therefore we get $u$ including the multiplier $\mu$, substitute $\mathbf{u}$ into $S^{(q)}$ and then we can get the expression of $\mu$, and finally get the expression of $\mathbf{u}$.

Inequality Constraint

\[\begin{align} &S(\mathbf{x},\mathbf{u},t)\leq0\\ &S(\mathbf{x},t)\leq0 \end{align}\]

For state-and-control-mixed constraint $S(\mathbf{x},\mathbf{u},t)\leq0$, the Hamilton function is given by the last section, but we have:

\[\mu=\left\{ \begin{array} >\geq0,& S=0\\ =0,& S<0 \end{array} \right.\]

Get boundary control $\mathbf{u}_b$ first, if $\mathbf{u}$ is continuous, we have $\mu(t_k^+)=\mu(t_k^-)=0$.

For pure-state constraint $S(\mathbf{x},t)\leq0$, we have two ways:

• Use $S^{(q)}$: Hamilton function is given by the last section.

\[\mu=\left\{ \begin{array} >\geq0,& S=0\\ =0,& S<0 \end{array} \right.\]

Tangency condition (interior condition):

\[N= \begin{bmatrix} S^{(0)}(\mathbf{x}(t_k),t_k)\\ S^{(1)}(\mathbf{x}(t_k),t_k)\\ \cdots\\ S^{(q-1)}(\mathbf{x}(t_k),t_k) \end{bmatrix} =\mathbf{0}\] \[\left\{ \begin{array} \mathbf{\lambda}^\intercal (t_k^-)=\mathbf{\lambda}^\intercal (t_k^+)+\mathbf{\pi}^\intercal \frac{\partial N}{\partial \mathbf{x}(t_k)}\\ H(t_k^-)=H (t_k^+)+\mathbf{\pi}^\intercal \frac{\partial N}{\partial t_k}\\ \end{array} \right.\]

• Use $S$:

\[\mathcal{H}=L+\mathbf{\lambda}^\intercal \mathbf{f}+\eta S\]

Solve boundary control first, then solve the multiplier $\eta$. The tangency condition still exists, but the Hamilton function $H$ should be changed into $\mathcal{H}$.

Example: Van der Pol’s Oscillator

The example comes from Maurer’s Lecture - Part 3. The model is shown below:

Here we take $a=-0.3$, and we solved the OCP in two ways.

Source code: Van der Pol Oscillator

Result

• PMP: Objective $J=7.8728$.
• GPOPS-II: Objective $J=7.8736$.