7  Solving optimal control problems

The preceding sections show how to setup and define an optimal control problem using PROPT. Many additional features have been implemented, such as advanced input error identification, to facilitate fast modeling and debugging.

7.1  Standard functions and operators

In the preceding sections several different Matlab functions and operators were used. This section lists the ones that are useful for the end user.

7.1.1   collocate — Expand a propt tomSym to all collocation points on a phase.

y = collocate(phase, x) for a m-by-n tomSym object x, on a phase with p collocation points, returns an p-by-m*n tomSym with values of x for each collocation point.

If x is a cell array of tomSym objects, then collocate is applied recursively to each element in the array.

7.1.2   tomSym/dot

Shortcut to overdot (alternatively dot product).

dot(p,x) gives the time derivative of x in the phase p.

dot(x) can be used to the same effect, if setPhase(p) has been called previously.

7.1.3   final — Evaluate a propt tomSym at the final point of a phase.

y = final(phase, x) for tomSym object x, returns an object of the same size as x, where the independent variable (usually t) has been replaced by its final value on the phase.

7.1.4   icollocate — Expand a propt tomSym to all interpolation points on a phase

y = icollocate(phase, x) is the same as y = collocate(phase, x), except that the interpolation points are used instead of the collocation points. This is typically useful when constructing an initial guess.

7.1.5   initial — Evaluate a propt tomSym at the initial point of a phase.

y = initial(phase, x) for tomSym object x, returns an object of the same size as x, where the independent variable (usually t) has been replaced by its initial value on the phase (often 0).

If x is a cell array of tomSym objects, then initial is applied recursively to each element in the array.

7.1.6   integrate — Evaluate the integral of an expression in a phase.

y = integrate(phase, x) for tomSym object x, returns an object which has the same size as x and is the integral of x in the given phase.

7.1.7   mcollocate — Expand to all collocation points, endpoints and midpoints on a phase.

y = mcollocate(phase, x) for a m-by-n tomSym object x, on a phase with p collocation points, returns an (2p+1)-by-m*n tomSym with values of x for each collocation point, the endpoints and all points that lie halfway inbetween these points.

The mcollocate function is useful in setting up inequalities that involve state variables. Because twice as many points are used, compared to collocate, the resulting problem is slightly slower to solve, but the obtained solution is often more correct, because overshoots in between collocation points are smaller.

Because it uses many more points than there are degrees of freedom, mcollocate should only be used for inequalities. Applying mcollocate to equalities will generally result in an optimization problem that has no solution. Care should also be taken to ensure that the mcollocated condition is not in conflict with any initial or final condition.

If x is a cell array of tomSym objects, then mcollocate is applied recursively to each element in the array.

If a finite element method is used, then mcollocate uses all points that are used in computing the numeric quadrature over elements.

7.1.8   setPhase — Set the active phase when modeling PROPT problem.

setPhase(p) sets the active phase to p.

It is not strictly necessary to use this command, but when doing so, it is possible to omit the phase argument to the commands tomState, tomControl, initial, final, integrate, etc.

7.1.9   tomControl — Generate a PROPT symbolic state.

``` x = tomControl creates a scalar PROPT control with an automatic name.
x = tomControl(phase,label) creates a scalar control with the provided name.
x = tomControl(phase,label,m,n) creates a m-by-n matrix of controls.
x = tomControl(phase,[],m,n) creates a matrix control with an automatic name.
x = tomControl(phase,label,m,n,'int') creates an integer matrix symbol.
x = tomControl(phase,label,m,n,'symmetric') creates a symmetric matrix symbol.
```

The tomControl symbols are different from tomState symbols in that the states are assumed to be continuous, but not the controls. This means that derivatives of tomControls should typically not be used in the differential equations, and no initial or final condition should be imposed on a tomControl.

If setPhase has been used previously, then the phase is stored in a global variable, and the phase argument can be omitted.

Constructs like "x = tomControl(’x’)" are very common. There is the shorthand notation "tomControls x".

7.1.10   tomControls — Create tomControl objects

Examples:

```   tomControls x y z
```

is equivalent to

```   x = tomControl('x');
y = tomControl('y');
z = tomControl('z');
```
```   tomControls 2x3 Q 3x3 -integer R -symmetric S
```

is equivalent to

```   Q = tomControl('Q', 2, 3);
R = tomControl('R', 3, 3, 'integer');
S = tomControl('S', 3, 3, 'integer', 'symmetric');
```

Note: While the "tomControls" shorthand is very convenient to use prototyping code, it is recommended to only use the longhand "tomControl" notation for production code. (See toms for the reason why.)

It is necessary to call setPhase before calling tomStates.

7.1.11   tomPhase — Create a phase struct

phase = tomPhase(label, t, tstart, tdelta, ipoints, cpoints) The ipoints (interpolation points) and cpoints (collocation points) input arguments must be vectors of unique sorted points on the interval 0 to 1.

phase = tomPhase(label, t, tstart, tdelta, n) automatically creates cpoints and ipoints using n Gauss points. (If n > 128 then Chebyshev points are used instead.)

phase = tomPhase(label, t, tstart, tdelta, n, [], ’cheb’) uses Chebyshev points instead of Gauss points. This yields better convergence for some problems, and worse for others, as compared to Gauss points.

7.1.12   tomState — Generate a PROPT symbolic state

``` x = tomState creates a scalar PROPT state with an automatic name.
x = tomState(phase,label) creates a scalar state with the provided name.
x = tomState(phase,label,m,n) creates a m-by-n matrix state.
x = tomState(phase,[],m,n) creates a matrix state with an automatic name.
x = tomState(phase,label,m,n,'int') creates an integer matrix symbol.
x = tomState(phase,label,m,n,'symmetric') creates a symmetric matrix symbol.
```

The tomState symbols are different from tomControl symbols in that the states are assumed to be continuous. This means that they have time derivatives, accessible via the dot() function, and that tomStates cannot be integers.

If setPhase has been used previously, then the phase is stored in a global variable, and the phase argument can be omitted.

Constructs like "x = tomState(’x’)" are very common. There is the shorthand notation "tomStates x".

7.1.13   tomStates — Create tomState objects as toms create tomSym objects

Examples:

```   tomStates x y z
```

is equivalent to

```   x = tomState('x');
y = tomState('y');
z = tomState('z');
```
```   tomStates 2x3 Q 3x3 R -symmetric S
is equivalent to
Q = tomState('Q', 2, 3);
R = tomState('R', 3, 3);
S = tomState('S', 3, 3, 'symmetric')
```
```   tomStates 3x1! v
is equivalent to
v1 = tomState('v1');
v2 = tomState('v2');
v3 = tomState('v3');
v = [v1;v2;v3];
```

See the help for "toms" for more info on how to use the different flags.

Note: While the "tomStates" shorthand is very convenient to use prototyping code, it is recommended to only use the longhand "tomState" notation for production code. (See toms for the reason why.)

It is necessary to call setPhase before calling tomStates.

Some user may wish to use more advanced function of PROPT as well. The following function can be useful in many situations.

7.2.1   atPoints — Expand a propt tomSym to a set of points on a phase.

y = atPoints(phase, t, x) for a m-by-n tomSym object x, and a vector t of length p, returns an p-by-m*n tomSym with values of x for each value of t (this is done via substitution).

If x is a cell array of tomSym objects, then atPoints is applied recursively to each element in the array.

If x is an equality or inequality, and the left-hand-side is divided by a symbolic phase.tdelta, then both sides are multiplied by phase.tdelta

```  Overloaded methods:
tomSym/atPoints
tomSym/atPoints
```

7.2.2   interp1p — Polynomial interpolation.

yi = interp1p(x,y,xi) fits a polynomial to the points (x,y) and evaluates that polynomial at xi to compute yi.

The behavior of interp1p is different from that of interp1 in that the same high-degree polynomial is fitted to all points. This is usually not a good idea, and for general interpolation interp1 is to be preferred. However, it works well in some cases, such as when the points x are the roots of a Legendre polynomial.

```  Overloaded methods:
tomSym/interp1p
tomSym/interp1p
```

7.2.3   proptGausspoints — Generate a list of gauss points.

[r, w] = gausspoints(nquad) = proptGausspoints(n)

Input: n - The number of points requested

Output r - The Gauss points (roots of a Legendre polynomial of order n.) w - The weights for Gaussian quadrature.

7.2.4   proptDiffMatrix — Differentiation matrix for interpPoly.

M = proptDiffMatrix(X,XI,N) creates a matrix M of the values of the N:th derivatives of the interpolating polynomials defined by the roots X, at the points XI. Each column of M corresponds to a specific root, and each row corresponds to a component of XI.

7.3  Screen output

The screen output depends on the functionality of the nonlinear programming solver. In most cases screen output can be generated with the following command:

`options.PriLevOpt = 1; % or higher for more output`

A print level of 3 is recommended when using the global NLP solver multiMin.

7.4  Solution structure

The solution structure from PROPT contains the relevant results.

It is possible to use the solution information to evaluate any custom expressions. For example:

` subs(x1.*x2,solution);`

7.5  Plotting

The best way to plot the data is to directly use the data fields as needed. There are also some automated plotting routines included, see for example help tomSym/ezplot.