In order to understand the basic modeling features in PROPT it is best to start with a simple example. Open the file called brysonDenham.m located in /tomlab/propt/examples.
There are also brysonDenhamShort.m, brysonDenhamDetailed.m, brysonDenhamTwoPhase.m that solve the same problem, utilizing different features of PROPT.
To solve the problem, simply enter the following in the Matlab command prompt:
>> brysonDenham
The following example can be found in vanDerPol.m in /tomlab/propt/examples.
The objective is to solve the following problem:
minimize:

subject to:

To solve the problem with PROPT the following compact code can be used:
toms t p = tomPhase('p', t, 0, 5, 60); setPhase(p); tomStates x1 x2 x3 tomControls u % Initial guess x0 = {icollocate({x1 == 0; x2 == 1; x3 == 0}) collocate(u == 0.01)}; % Box constraints cbox = {10 <= icollocate(x1) <= 10 10 <= icollocate(x2) <= 10 10 <= icollocate(x3) <= 10 0.3 <= collocate(u) <= 1}; % Boundary constraints cbnd = initial({x1 == 0; x2 == 1; x3 == 0}); % ODEs and path constraints ceq = collocate({dot(x1) == (1x2.^2).*x1x2+u dot(x2) == x1; dot(x3) == x1.^2+x2.^2+u.^2}); % Objective objective = final(x3); % Solve the problem options = struct; options.name = 'Van Der Pol'; solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
It is possible to compile permanently, in order to keep the autogenerated code:
>> Prob = sym2prob('lpcon',objective, {cbox, cbnd, ceq}, x0, options); >> Prob.FUNCS ans = f: 'lp_f' g: 'lp_g' H: 'lp_H' c: 'c_AFBHVT' dc: 'dc_AFBHVT' d2c: 'd2c_AFBHVT' r: [] J: [] d2r: [] fc: [] gdc: [] fg: [] cdc: [] rJ: [] >> edit c_AFBHVT
The code that was generated can be found in the $TEMP directory. The objective in this case is linear and can found in the Prob structure (Prob.QP.c).
Here is what the auto generated constraint code may look like:
function out = c_AFBHVT(tempX,Prob) % c_AFBHVT  Autogenerated file. tempD=Prob.tomSym.c; x3_p = reshape(tempX(183:243),61,1); u_p = reshape(tempX(1:60),60,1); x2_p = reshape(tempX(122:182),61,1); x1_p = reshape(tempX(61:121),61,1); tempC7 = tempD{2}*x2_p; tempC8 = tempC7.^2; tempC10 = tempD{2}*x1_p; out = [(tempD{3}tempC8).*tempC10tempC7+u_p0.2*(tempD{1}*x1_p);... tempC10.^2+tempC8+u_p.^20.2*(tempD{1}*x3_p)];
And the objective function to optimize (in this case a simple linear objective already available in TOMLAB (hence not auto generated, but defined by the Prob field used)):
function f = lp_f(x, Prob) f = Prob.QP.c'*x(:);
The PROPT system uses equations and expressions (collected in cells) to model optimal control problems.
Equations must be written either using (== <= >=) equality/inequality markers.
It is possible to include more than one equation on the same line.
For example:
toms a b c cnbd = {a == b; b == c};
does the same job as:
toms a b c cnbd = {a == b b == c};
The same is true for inequalities.
When working with optimal control problems, one typically work with expressions that are valid for all collocation points. The functions collocate and icollocate can be used to extend an arbitrary expressions to the necessary collocation points.
Consider for example the starting points:
tomStates x1 tomControls u1 x0 = {icollocate(x1==1); collocate(u1==1)};
Note: icollocate, as it is working with a state variable, also includes the starting point.
Scale the problem:
Proper scaling may speed up convergence, and conversely, improperly scaled problems may converge slowly or not at all. Both unknowns and equations can be scaled.
Don’t use inf in equations:
It is strongly discouraged to use inf/inf in equations. This is because the equation may be evaluated by subtracting its lefthand side from the righthand side and comparing the result to zero, which could have unexpected consequences.
Equations like x<=Inf are better left out completely (Although in many cases tomSym will know to ignore them anyway).
Avoid using nonsmooth functions:
Because the optimization routines rely on derivatives, nonsmooth
functions should be avoided. A common example of a nonsmooth function is
the Heaviside function H, defined as H(x) = 1 for x
> 0, and H(x) = 0 for x ≤ 0,
which in Matlab code can be written as (x>0). A smoother
expression, which has the same behavior for x >> a is:
H_{s} = 1/2 ( 1 + tanh( x/a ) )
When the problems are nonsmooth in the independent variable (time,
t), the problem should normally be separated into phases.
Use descriptive names for equations, states, controls, and variables:
The names used for states, controls, and variables make the code easier to read. Consider using names such as ”speed”, ”position”, ”concentration”, etc. instead of ”x1”, ”x2”, ”x3”.
The names used for equations do not matter in most cases. However,
they will be useful when accessing Lagrange multipliers.
Resolve on successively finer grids:
If a very fine grid is needed in order to obtain a good solution, it is usually a good idea to solve the problem on a less dense grid first, and then re solve, by using the obtained solution as a starting point. The following code will do the trick:
for n=[10 40] p = tomPhase('p', t, 0, 6, n); setPhase(p); tomStates x1 x2 % Initial guess if n == 10 x0 = {p1 == 0; p2 == 0}; else x0 = {p1 == p1opt; p2 == p2opt icollocate({x1 == x1opt; x2 == x2opt})}; end ... ... ... solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options); % Optimal x, p for starting point x1opt = subs(x1, solution); x2opt = subs(x2, solution); p1opt = subs(p1, solution); p2opt = subs(p2, solution); end
See for example denbighSystem.m and drugScheduling.m.
The examples included with the software in many cases just scratch the surface of the capabilities of PROPT. The software is designed for maximum user flexibility and with many "hidden" features. The following notes are worth keeping in mind when modeling optimal control problems.
Lefthand side treatment:
The left hand side of equations do not have to be a single timederivative. Things like collocate(m*dot(v) == F) work just fine, even if m is a state variable (A constraint that prevents m from becoming to small may improve convergence).
Second derivatives:
Second derivatives are allowed, as in collocate(dot(dot(x)) == a), although convergence is almost always better when including extra states to avoid this.
Equations:
Equations do not have to include any timederivative. For example collocate(0.5*m*v^{2} + m*g*h)==initial(0.5*m*v^{2} + m*g*h) will work just fine in the PROPT framework.
Fixed end times:
Problems with fixed endtimes enables the use of any form of expression involving the collocation points in time since the collocation points are predetermined.
The following illustrates how t can (should) be used.
toms t % Will cause an error if myfunc is undefined for tomSym myfunc(t); % This works since collocate(t) is a double vector and not a % tomSym object. myfunc(collocate(t))
It is also worth noting that feval bypasses tomSym, so it is possible to call any function, at the expense of slower derivatives. In the case of a fixed end time, there are no derivatives:
toms t % Works since tomSym is bypassed. collocate(feval('myfunc',t))
Segmented constraints:
With PROPT is it possible to define constraints that are valid for only a part of phase. In addition, it is possible to define constraints for other points than the collocation points which could assist in avoiding constraint violations in between the collocation points.
mcollocate automates this process for the entire phase and imposes twice the number of constraints, i.e. the solution will be constrained in between each collocation point as well.
The following example shows how to set segmented constraints (valid for a part of the phase) for timefree and timefixed problems.
% For variable end time (the number of constraints need to be fixed) % Relax all constraints below the cutoff, while 0 above con = {collocate((8*(t0.5).^20.5x2) >= 1e6*(t<0.4))}; % When the end time is fixed there are many possible solutions. % One can use custom points or a set of collocation points if needed. tsteps = collocate(t); n = length(tsteps); y = atPoints(p,tsteps(round(n/3):round(n/2)),(8*(t0.5).^20.5x2)); con = {y >= 0}; % Alternatively one could write for custom points. con = atPoints(p,linspace(0.4, 0.5, 10),(8*(t0.5).^20.5x2) >= 0);
Constant unknowns:
It is possible to mix constant unknowns (created with "toms") into the equations, and those will become part of the solution: icollocate(0.5*m*v^{2} + m*g*h == Etot).
Higherindex DAEs:
Higherindex DAEs usually converge nicely with PROPT. One can use dot() on entire expressions to get symbolic time derivatives to use when creating lower index DAEs.
Lagrangian equations:
Deducing Lagrangian equations manually is usually not necessary, but if needed, one can use tomSym for derivatives. TomSym allows for computation of symbolic derivatives with respect to constant scalars and matrices using the "derivative" function, but this does not work for tomStates or tomControls. The current implementation is therefore to create expressions using normal tomSym variables (created with the "toms" command), use the derivative function and then use "subs" function to replace the constants with states and controls before calling collocate.
Independent variables can be defined as follows (yes, it is that simple!):
toms t % Time toms tf % Final time % Phase name, time, start, end, number of nodes p = tomPhase('p', t, 0, tf, 25); cbox = {20 <= tf <= 30}; % Bounds on end time x0 = {tf == 22); % Starting point for end time (must be set before used % in other expressions)
It is also possible use scalar variables within the PROPT framework:
toms p1 cbox = {1 <= p1 <= 2}; x0 = {p1 == 1.3);
where the variables define lower, upper bound and a suitable guess.
A constant variable ki0 can be defined with the following statement (i.e. no difference from normal Matlab code):
ki0 = [1e3; 1e7; 10; 1e3];
The difference between state and control variables is that states are continuous between phases, while control variables can be discontinuous.
Examples:
Unconstrained state variable x_{1} with a starting guess from 0 to 1:
tomStates x1 x0 = icollocate(x1 == t/tf);
State variable x_{1} with bounds of 0.5 and 10. The starting guess should be set as well to avoid any potential singularities:
tomStates x1 cbox = {0.5 <= icollocate(x1) <= 10};
Control variable T with a "table" defining the starting point:
tomControls T x0 = collocate(T==273*(t<100)+415*(t>=100));
Boundary constraints are defined as expressions since the problem size is reduced. For example if state variable x_{1} has to start in 1 and end in 1 and x_{2} has to travel from 0 to 2 define:
cbnd1 = initial(x1 == 1); cbnd2 = initial(x2 == 0); cbnd3 = final(x1 == 1); cbnd4 = final(x2 == 2);
A variety of path, event and integral constraints are shown below:
x3min = icollocate(x3 >= 0.5); % Path constraint integral1 = {integrate(x3)  1 == 0}; % Integral constraint final3 = final(x3) >= 0.5; % Final event constraint init1 = initial(x1) <= 2.0; % Initial event constraint
NOTE: When defining integral constraints that span over several phase, use either of the following notations:
integrate(p1,x3p1) + integrate(p2,x3p1) + integrate(p3,x3p1) <= 2.1e3
or with expressions:
p1qx3 == integrate(p1,x3p1); p2qx3 == integrate(p2,x3p1); p3qx3 == integrate(p3,x3p1); qx2sum = {p1qx3 + p2qx3 + p3qx3 <= 2.1e3};