Definitions

# PROPT

The PROPT MATLAB Optimal Control Software is a new generation platform for solving applied optimal control (with ODE or DAE formulation) and parameters estimation problems.

The platform was developed by MATLAB Programming Contest Winner, Per Rutquist in 2008.

## Description

PROPT is a combined modeling, compilation and solver engine for generation of highly complex optimal control problems. PROPT uses a pseudospectral collocation method for solving optimal control problems. This means that the solution takes the form of a polynomial, and this polynomial satisfies the DAE and the path constraints at the collocation points.

In general PROPT has five main functions:

• Computation of the constant matrices used for the differentiation and integration of the polynomials used to approximate the solution to the trajectory optimization problem.
• Text manipulation to turn user-supplied expressions into MATLAB code for the cost function $f$ and constraint function $c$ that are passed to a nonlinear programming solver in TOMLAB, ensuring that the code is compatible with MAD (TOMLAB package for automatic differentiation to achieve floating point precision for derivatives).
• Functionality for plotting and computing a variety of information for the solution to the problem.
• Partial automatic linearization for the following scenarios:
• Minimal (or maximal) time problem
• Problems with a linear objective (pure feasibility problems are also handled)
• Linear systems with a fixed end time
• Integrated support for non-smooth (hybrid) optimal control problems.

## Modeling

The PROPT system uses equations and expressions to model optimal control problems. It is possible to define independent variables, dependent functions, scalars and constant parameters:

t = proptIndependent('t', [0 0 0], [10 20 40]); problem.variables.z1 = proptScalar(0, 500 , 10); problem.parameters.ki0 = [1e3; 1e7; 10; 1e-3];

### States and Controls

States and controls only differ in the sense that controls need be continuous between phases.

x.x1 = proptState(-inf , inf , [0; 1]); x.u1 = proptControl(-2 , 1 , 0);

### Boundary, path, event and integral constraints

A variety of boundary, path, event and integral constraints are shown below:

expr.x1_i = '1'; % Starting point for x1 expr.x1_f = '1'; % End point for x1 expr.x2_f = '2'; % End point for x2 eq.x3min = proptEquation('x3 > 0.5'); % Path constraint for x3 eq.integral = proptEquation('quad_t *x2 - 1 = 0'); % Integral constraint for x2 eq.final3 = proptEquation('x3_f > 0.5'); % Final event constraint for x3 eq.init1 = proptEquation('x1_i < 2.0'); % Initial event constraint x1

## Single Phase Optimal Control Example

Van der Pol Oscillator

Minimize:

$begin\left\{matrix\right\} J_\left\{x,t\right\} & = & x_3\left(t_f\right) end\left\{matrix\right\}$

Subject to:

$begin\left\{cases\right\} frac\left\{dx_1\right\}\left\{dt\right\} = \left(1-x_2^2\right)*x_1-x_2+u frac\left\{dx_2\right\}\left\{dt\right\} = x_1 frac\left\{dx_3\right\}\left\{dt\right\} = x_1^2+x_2^2+u^2$

` x(t_0) = [0  1  0] `
` t_f = 5 `
` -0.3 le u le 1.0 `
end{cases}

To solve the problem with PROPT the following code can be used (with 60 collocation points):

tF = 5; t = proptIndependent('t', 0, tF);

clear x eq expr x.x1 = proptState(-10, 10, 0); x.x2 = proptState(-10, 10, 1); x.x3 = proptState(-10, 10, 0); x.u = proptControl(-0.3, 1, -0.01);

expr.x1_i = '0'; expr.x2_i = '1'; expr.x3_i = '0';

eq.eq1 = proptEquation('x1_t = (1-x2.^2).*x1-x2+u'); eq.eq2 = proptEquation('x2_t = x1'); eq.eq3 = proptEquation('x3_t = x1.^2+x2.^2+u.^2');

problem = proptPhase(t, x, [], eq, expr, 60); problem.cost = 'x3_f'; problem.name = 'Van Der Pol'; problem.language = 'Matlab, vectorized';

solution = proptSolve(problem);

## Multi Phase Optimal Control Example

One-dimensional rocket with free end time and undetermined phase shift

Minimize:

$begin\left\{matrix\right\} J_\left\{x,t\right\} & = & tCut end\left\{matrix\right\}$

Subject to:

$begin\left\{cases\right\} frac\left\{dx_1\right\}\left\{dt\right\} = x_2 frac\left\{dx_2\right\}\left\{dt\right\} = a-g \left(0 < t <= tCut\right) frac\left\{dx_2\right\}\left\{dt\right\} = -g \left(tCut < t < t_f\right)$

` x(t_0) = [0  0] `
` g = 1 `
` a = 2 `
` x_1(t_f) = 100 `
end{cases}

The problem is solved with PROPT by creating two phases and connecting them:

tF = 100; tCut = 10; t0 = 0; t1 = proptIndependent('t', 0.1, [tCut-9 tCut tCut+9]); t2 = proptIndependent('t', [tCut-9 tCut tCut+9], [tCut+9 tCut+9 tF]);

clear x1 x2 eq1 eq2 expr1 expr2 x1.x1 = proptState(0, inf, [0 50]); x1.x2 = proptState(0, inf, 0);

x2.x1 = proptState(0, inf, [50 100]); x2.x2 = proptState(0, inf, 0);

expr1.x1_i = '0'; expr1.x2_i = '0'; expr1.a = '2'; expr1.g = '1';

expr2.x1_f = '100'; expr2.a = '2'; expr2.g = '1';

eq1.eq1 = proptEquation('x1_t = x2'); eq2 = eq1; eq1.eq2 = proptEquation('x2_t = a-g'); eq2.eq2 = proptEquation('x2_t = -g');

problem = proptProblem; problem.name = 'One Dim Rocket'; problem.language = 'Matlab, vectorized';

problem.phases.p1 = proptPhase(t1, x1, [], eq1, expr1, 20); problem.phases.p2 = proptPhase(t2, x2, [], eq2, expr2, 20); problem.phases.p2.cost = 't_f_p1'; problem = proptAutoLink(problem, 'p1', 'p2');

solution = proptSolve(problem);

## Parameter Estimation Example

Parameter Estimation Problem

Minimize:

$begin\left\{matrix\right\} J_\left\{p\right\} & = & sum_\left\{i=1,2,3,5\right\}\left\{\left(x_1\left(t_i\right) - x_1^m\left(t_i\right)\right)^2\right\} end\left\{matrix\right\}$

Subject to:

$begin\left\{cases\right\} frac\left\{dx_1\right\}\left\{dt\right\} = x_2 frac\left\{dx_2\right\}\left\{dt\right\} = 1-2*x_2-x_1$

` x_0 = [p_1  p_2] `
` t_i = [1  2  3  5] `
` x_1^m(t_i) = [0.264  0.594  0.801  0.959] `
|p_{1:2}| <= 1.5 end{cases}

In the code below the problem is solved with a fine grid (10 collocation points). This solution is subsequently fine-tuned using 40 collocation points:

tF = 6; t = proptIndependent('t', 0, tF);

clear x eq expr vbl

x.x1 = proptState; x.x2 = proptState;

expr.x1meas = '[0.264;0.594;0.801;0.959]'; expr.tmeas = '[1;2;3;5]'; expr.x1err = 'sum((interp1(t,x1,tmeas) - x1meas).^2)';

eq.eq1 = proptEquation('x1_t = x2'); eq.eq2 = proptEquation('x2_t = 1-2*x2-x1');

expr.x1_i = 'p1'; expr.x2_i = 'p2';

vbl.p1 = proptScalar(-1.5,1.5,0); vbl.p2 = proptScalar(-1.5,1.5,0);

problem = proptPhase(t, x, vbl, eq, expr, 20); problem.cost = 'x1err'; problem.name = 'Parameter Estimation'; problem.language = 'Matlab, vectorized';

solution = proptSolveIter(problem,[10;40]);