Home

Imperial College London Optimal Control Software User Guide

1. Plot the solutions figure 1 hold on plot time states title States vs time xlabel Time figure 2 hold on plot time inputs title Inputs vs Time xlabel Time figure 3 hold on plot states 1 states 2 T States vs time 1 5 T T 0 57 0 2 4 6 8 10 12 14 16 18 20 Time Figure 25 State trajectories for the continuous time nonlinear Model Predictive Control example title ptimal state trajectory xlabel x 1 ylabel x 2 Notice that mainMPC m calls the following two additional files testPlant m It returns the ODE right hand side of the model describing the plant applyControl m It applies the optimal control stored in the variable solution to the system described in the file testPlant m The state and control variables solution to this problem using the Optimal Control Toolbox are shown in Figs 25 26 and 27 4 10 Example 10 Singular Arc problem Find the control input u R over t in 0 ty solving the following optimisation problem 7 5 p 57 71 nputs vs Time 1 5 T T T 0 2 4 6 8 10 12 14 16 18 20 Figure 26 Input trajectory for the Continuous time nonlinear Model Predictive Control example Optimal state trajectory T e 0 5r C N 4 1 1 1 1 1 1 1 1 Eri 0 2 0 0 2 0 4 0 6 0 8 1 12 1 4 Figure 27 Phase plane for the continuous time nonlinear Model Predictive
2. dE dx0 dE du0 dE dxf dE duf The gradient of the stage cost and the boundary cost are supplied and so dL flag 1 and dE flag 1 The final time ty and p are not decision variables and then the deriva tives dL dp dL dt dE dtf and dE dp are set to be empty matrices The derivative of the boundary cost does not depend on any variable and then all variables containing the derivatives are set to be empty matrices jacConst m function df dg db jacConst f g X U P t b xO xf uO uf p tf tO data Lt ones size t coefedata data T1 1 data data mdrt df flag 1 df dp 1 1 Of a u p t Op df dt 1 f z u p t t Of a u p t Ox1 42 4 11 10 16 X 3 72 X 1 72 42 X 1 73 X 3 xX 2 X 1 72 1 Of x u p t Ox df dx 2 1 Lt O Lt X 3 X 1 Of x u p t Ox3 df dx 3 0 Lt 2 X 3 X 1 X 2 X 1 1 df du 1 0 Lt coef cos U 1 coef sin U 1 dg flag 0 db flag 1 52 db dtf db dp db dx0 db du0 db duf Yo Derivatives with respect to z1 ty va tr va ty db dxf 1 sqrt xf 1 xf 1 2 0 1 0 1 01 The path constraints are not present and then their derivatives have not be sup plied dg flag 0 Instead the derivatives of the boundary constraints are supplied in the structured variable db as shown in the illustrated example The boundary constraints b tf p o uo
3. 0z30u 07 fz Ox20u 8 fa Ox2du Hf 2 5 0 t 1 500 X 4 O t O t 0 fi 0130u 0 fo Org0u 0 fs 8x30u 0 fa Ox30u Hf 3 5 0 t O t 1 500 X 4 O t 10 fi Oxs0u 0 f2 0v40u 0 fs za u 9 fa Oxs0u Hf 4 5 X 1 500 X 4 72 X 2 500 X 4 72 1 X 4 72 1 X 3 500 O t 0 fi 0w 0 fo Ou 0 fs Ou 0 fa Ou Hf 5 5 0 t O t O t O t Lt ones size t HL cell nfz nfz HL 1 1 0 t 96 9 L 0x HL 1 2 0 t O L 02021 HL 2 2 0 t 96 O L 8x3 HL 1 3 0 t 96 O L 021023 32 HL 2 3 0 t 96 0 L Ox20z3 HL 3 3 0 t O L 923 HL 1 4 0 t O L 0210c4 HL 2 4 0 t O L 020x4 HL 3 4 0 t 96 O L 813024 HL 4 4 0 t 96 9 L 822 HL 1 5 0 t 96 0 L Ox10u HL 2 5 0 t 96 0 L Ox20u HL 3 5 0 t 96 O L 8x38u HL 4 5 0 t O L 8x40u HL 5 5 2 0 00001 Lt 0 L 8u nE n Ez zeros nE nE HE num2cell Ez HE 2 nE 1 Hg 1 Hb Each Hessian is specified in a two dimensional cell array of dimension given by the size of y z ul or y xy Each entry i j of the cell array corresponds to the derivative with respect to y i and y j The order to follow is given by the order of the vector variable
4. https computation llnl gov casc sundials main html It is highly recommended that the free NLP solver IPOPT be installed since this will dramatically improve functionality and performance of the code If you would like to use a multiple shooting method to solve your problem the sensitivity solver CVODES it comes with the SUNDIALS suite also has to be installed Before running the code you will need to include the IPOPT libraries in the path by using setenv or possibly export commands and include the file Tpopi Contrib MatlabInterface ipopt mexglz in your MATLAB path To compile the CVODES mex file simply add the sundialsTB directory to your MATLAB path and run startup STB m Since the optimal control code consists only of m files no installation is necessary but don t forget to add ICLOCS src to your MATLAB path The current version of ICLOCS has been tested under Linux Red Hat Enterprise Linux Version 5 2 and Ubuntu 9 0 4 with MATLAB 7 6 0 The compilation of IPoPT has been performed using the compilers gcc 4 1 and gfortran 3 Solving Optimal Control Problems This section details the procedure for defining and solving optimal control problems using the MATLAB code provided In general the following steps have to be performed 1 Copy main m settings m myProblem m from ICLOCS usr to your working directory The other files in this directory callback m gradCost m jacConst m and hessianLagrangian m should only be cop
5. 0 35 0 15 0 3 0 1 0 25 0 05 o 2t 1 of Qe x x 0 15 0 05 4 0 1 0 05 N 1 oast 0 1 ii 0 2 Res 1 1 1 1 0 20 40 60 80 100 0 20 40 60 80 100 Time Time Figure 13 State trajectories for the discrete time optimization problem dx 1 x1 0 1 x2 dx 2 1 8 x2 0 0787 u function c g x u p t data c function bc b x0 xf u0 uf p tf data P data P bc xf P xf e The solution method and solver settings are set in settings Dis m See the file included in the directory ICLOCS examples Discrete TimeMPC The solution of the optimisation problem is computed running the file main m inside the directory ICLOCS examples Discrete TimeMPC The state control and adjoint variables solution to this problem using the Optimal Control Toolbox are shown in Figs 12 El and 3 4 6 Example 6 Minimum fuel orbit raising problem Find the control u t over t in to ty solving the following optimisation problem 3 47 MPC 0 5 se Input T 1 5 5 10 20 30 40 50 60 70 80 90 100 Time Figure 14 Input trajectory for the discrete time optimization problem Cost 0 8 0 4 4 0 2 N 1 1 10 20 30 40 50 60 70 80 90 100 Steps Figure 15 Cost for the discrete time optimization problem 48 ty min z t dt ul 0 subject to d 12 2 o gt T sin u r 22 1 mdt 12313 T cos u da 21 1
6. 3 0 4 1 j x NI 0 0 5 1 1 5 2 Figure 35 Input trajectory for the two link robot arm problem Adjoint variables 8 d Time Figure 36 Adjoint variables for the two link robot arm problem 85 2 A E Bryson Dynamic Optimization Addison Wesley Longman Menlo Park CA USA 1999 3 A E Bryson and Y C Ho Applied optimal control John Wiley amp Sons 1975 4 J E Cuthrell and L T Biegler Simultaneous optimization and solution methods for batch reactor control profiles Comput Chem Eng 13 49 62 1989 5 Wendell H Fleming and Raymond W Rishel Deterministic and stochastic optimal control Springer Verlag Berlin New York 1975 6 G Hicks and W Ray Approximation methods for optimal control synthesis The Canadian Journal of Chemical Engineering 49 522 528 1971 7 L S Jennings and M E Fisher MISER3 optimal control toolbox User Manual Department of Mathematics The University of Western Australia 8 S Kameswaran and L T Biegler Simultaneous dynamic optimization strategies recent advances and challenges Computers and Chemical Engineering 30 1560 1575 2006 9 Y Chen Liang M Meng and R Fullmer Solving tough optimal control problems by network enabled optimization server neos Technical report School of Engineering Utah State University USA Chinene University of Hong Kong China 2003 10 Rein Luus Iterative Dynamic Prog
7. HL 1 1 2 HL 2 2 2 HL 3 3 2 0 005 HE Gz zeros nfz nfz Hg num2cell Gz Hb e The solution of the optimisation problem is computed running the file main m inside the di rectory ICLOCS examples PathConstraint The state control and adjoint variables solution to this problem using the Optimal Control Toolbox are shown in Figs 20 and 21 4 8 Example 8 Continuously stirred tank reactor The following optimisation problem is based on the model of a continuously stirred tank reactor proposed in 6 and modified in 8 60 Optimal state traiectories 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 Time Figure 19 State trajectories for the path constrained optimal control problem Optimal input sequences T T T T il 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 Time Figure 20 Input trajectory for the path constrained optimal control problem 61 Adjoint variables 0 15 F i l X 1 2 we IN Qu i i s i i i l i 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 Time Figure 21 Adjoint variables for the path constrained optimal control problem min f ea walt s fi ult ult u tr subject to 1 En t kaye 23 Tf _ En t2 1 kriel ag au x2 Te q 0 10 98 0 39 ai tf 0 26 1 0 65 u t 0 76 0x u t lt 2 where a 0
8. The gradient of the stage cost and the boundary cost are supplied and so dL flag 1 and dE flag 1 The final time 2 is a decision variable and then the derivative of the stage cost with respect to the time t is set to 0 t It is zero at any time instant The same rule applies for the derivatives of the stage cost with respect to the state and the input The derivative of the stage cost with respect to the state x is stored in a matrix with n columns and a number of rows given by the different evaluation times Each row corresponds to the evaluation of the derivative at a given time instant The derivative of the boundary cost depends only on the final time and then the variables containing the derivatives with respect to the other decision variables are set to be empty matrices jacConst m function df dg db jacConst f g X U P t b xO xf uO uf p tf tO data 96 vector of ones with the same size of the number of point evaluation in time Lt ones size t df flag 1 df dp 1 96 Derivatives of f x u p wrt p df dt 1 O Lt O Lt 96 Derivatives of f a u p 2 wrt t df dx 1 0 Lt O Lt Derivatives of f z u p wrt 71 df dx 2 1 Lt 0 Lt Derivative of f x u p wrt z df du 1 0 Lt 1 Lt 96 Derivative of f x u p t wrt u dg flag 0 db flag 0 The path constraints and boundary constraints are not present and then their deriva tives have not be supplied dg flag 0 and db flag 0 Instead the deriv
9. problem constraints gu inf 10 96 Bounds for boundary constraints problem constraints bl problem constraints bu 96 store the necessary problem parameters used in the functions problem data problem functions L E Of g b function stageCost L x xr u ur p t data x1 x 1 x2 x 2 ul u 1 stageCost x1 72 x2 72 0 005 u1 72 57 function boundaryCost E x0 xf u0 uf p tf data boundaryCost 0 function dx f x u p t data x1 x 1 x2 x 2 ul u 1 dx 1 x2 dx 2 x2 u1 function c g x u p t data x2 x 2 c 8 t 0 5 2 0 5 x2 x 1 function bc b x0 xf u0 uf p tf data be e The solution method and solver settings are set in settings m See the file included in the directory ICLOCS examples PathConstraint e The files gradCost m jacConst m and hessianLagrangian m for this example are sup plied See inside the directory ICLOCS examples PathConstraint gradCost m function dL dE gradCost L X Xr U Ur P t E xO xf uO uf p tf data dL Gradient of the stage cost L wrt t p 7 u dL flag 1 dL dp dL dt dL dx 2 X 1 2 X 2 dL du 2 0 005 U 1 Yo dE Gradient of E with respect to tf p zo uo Uf vf dE flag 1 dE dtf dE dp 58 dE dx0 4 0 1 dE dxf dE duf jacConst m function df dg db jacConst f g X U P t b xO xf uO uf p tf tO data Lt ones size t
10. df flag 1 df dp 1 df dt 1 df dx 1 0 Lt O Lt df dx 2 1 Lt 1 Lt df du 1 0 Lt 1 Lt dg flag 1 dg dp 1 dg dt 1 16 t 0 5 O t dg dx 1 0 Lt 1 Lt dg dx 2 1 Lt O Lt dg du 1 0 Lt O Lt db flag 1 db dtf db dp L1 db dx0 L1 db du0 L1 db dxf db duf The path constraints are present and their derivatives have been supplied dg flag 1 The derivatives of g z u p gi v up t gn x u p t are supplied in the structured variable dg as shown in the illustrated example The derivative of g x u p t with respect to x are stored in dg dx which is a cell array where each entry corresponds to the derivative with respect to a state variable For instance the deriva tive of g x u p t with respect to x is stored in dg dx i dg dx i is a matrix with ng columns and a number of rows given by the different evaluation times Notice that the derivative of the path constraints with respect to the time has been specified but it was not necessary since the final time is not a decision variable hessianlagrangian m 59 function HL HE Hf Hg Hb hessianLagrangian X U P t E x0 xf u0 uf p tf data The Hessian of the different parts of the Lagrangian HL HE Hf Hg Hb must be supplied in cell arrays as follows nt np n m deal data sizes 1 4 nfz nt nptnt m Fz zeros nfz nfz Hf num2cell Fz Hz zeros nfz nfz HL num2cell Hz
11. T T T Figure 29 Input trajectory for the singular arc problem 76 Adjoint variables i t 0 5 4 X oL fee 4 0 5 mn 0 5 1 1 5 2 5 3 3 5 4 2 Time Figure 30 Adjoint variables for the singular arc problem 1 min x x3 x3 0 00001u7 dt tuis 0 subject to T xo io u 1 1 x1 0 1 x2 0 1000 lt u t lt 1000 39 s 28 ie vet a Problem setup The Optimal Control Problem is defined in the file Speyer m in the following way Initial time to lt ty problem time t0 0 Final time Let tf min tf max if ty is fixed problem time tf_min 1 problem time tf_max 1 77 guess tf 1 Parameters bounds pl lt p lt pu problem parameters pl problem parameters pu guess parameters Initial condition z t9 and its bounds problem states xO problem states x01 100 100 problem states xOu 100 100 96 State bounds xl x t xu problem states x1l 100 100 problem states xu 100 100 Yo Terminal state bounds xfl z ty lt xfu problem states xfl 100 100 problem states xfu 100 100 Yo Guess the state trajectories with x to x ty guess states 1 1 1 guess states 2 1 1 96 Number of control actions N 96 Set problem inputs N 0 if N is equal to the number of integration steps problem inputs N 0 Input bounds ul lt u t uu problem inputs ul 1000 p
12. lowerbound xn tf lowerbound problem states xfu xi tf upperbound xn tf upperbound If the final states are unbounded their bounds have to be set to inf or inf Guess state trajectories By default the initial guess for the state trajectories is auto matically generated by linearly interpolating between the expected initial and final value for each state which are provided as shown below guess states 1 x1 tO x1 tf guess states n xn tO xn tf If the variable guess states is the empty matrix the initial trajectories will be generated by linearly interpolating between random but feasible initial and final values for each state Note that the initial guess generated here can be overwritten and a user supplied guess can be assigned in main m by defining the variable infoNLP zO see Remark D Number of control actions The number of piecewise constant control actions can be defined here For direct collocation methods N 0 sets the number of control actions equal to the number of integration steps Note that the number of integration steps M 1 defined in settings m have to be divisible by the number of control actions For mul tiple shooting the number of control actions is equal to the number of integration steps N M 1 problem inputs N number of control actions 10 11 12 13 Control inputs Upper and lower bounds for the control inputs are defined as follows problem
13. m dx 1 f1 x1 xn ul um p t dx n fn x1 xn ul um p t If the i ODE right hand side does not depend on variables it is necessary to multiply the assigned value by a vector of ones with the same length of t in order to have a vector with the right dimension when called for the optimization Example dx i O ones size t 1 Path constraint function The path constraint function c g x u p t data is defined in a manner similar to the system equations It returns a row vector at each point in time Each entry correspond to the evaluation of a constraint for a given x u p for a point in time t Again this function should be vectorised as follows x1 ul x 1 xn x n u 1 um u m c 1 7gi x1 u1 p t c ng gag x1 xn ui um p t where ng is the number of constraints If the problem does not have any path constraints let 1 Boundary constraint function The boundary constraint function bc b x0 xf u0 uf p tf data returns a column vector corresponding to the evaluation of each boundary constraint If the problem does not have any path constraints let bc 3 2 Choosing Solver Settings Once the optimal control problem has been defined in myProblem m the solver methods and settings have to edited in the file settings m 1 Transcription Method As a first step the transcription method has to be chosen by setting the variable options transeription to eit
14. pi lowerbound problem parameters pu p1_upperbound 1 guess parameters pi guess p2 guess we Define all lower and upper bounds on the parameters as entries in a row vector If param eters are unbounded their bounds can be set to inf or inf As before an initial guess for the unknown parameters should be provided If there are no unknown parameters that can be optimised over set the bound and initial guess vectors to 4 Initial conditions The bounds for the initial condition xg of the system have to be defined in row vectors as shown in the following lines problem states x01 x1 t0 _lowerbound xn t0 _lowerbound problem states xOu x1 t0 _upperbound xn t0 _upperbound problem states x0 x1 60 xn t0 If the initial conditions are fixed let problem states xOl problem states xOu Note that there will be n bounds for a system with n states When the initial condition belongs to a box a value for zo can be assigned in problem states x0 it can be a guess or a desired value Otherwise problem states x0 can be the empty matrix State variables Box constraints for the state variables at to lt t lt ty are defined here problem states xl xi lowerbound xn lowerbound problem states xu xi upperbound xn upperbound If the states are unbounded their bounds have to be set to inf or inf Final state Bounds on the final state at t t are specified here problem states xfl xi tf
15. x2 x 2 ul u 1 dx 1 x2 19 dx 2 ul function c g x u p t data e 1 function bc b x0 xf u0 uf p tf data be The solution method and solver settings are set in settings m See the file included in the directory ICLOCS examples BangBang Notice that the multiple shooting option can not be used since the final time ty is a decision variable The files gradCost m jacConst m and hessianLagrangian m for this example are sup plied See inside the directory ICLOCS examples BangBang gradCost m function dL dE gradCost L X Xr U Ur P t E x0 xf u0 uf p tf data get the dimension of the state n and of the input m n m deal data sizes 3 4 vector of ones with the same size of the number of point evaluation in time Lt ones size T dL flag 1 dL dp 96 Derivatives of L x u p t wrt p dL dt 0 t 96 Derivatives of L z u p t wrt t Derivatives of L x u p t wrt x zi n dL dxekron zeros 1 n Lt Derivatives of L z u p wrt u u1 Um dL du kron zeros 1 m Lt dE flag 1 dE dtf 1 96 Derivatives of E xo z 7 uo uf p tf wrt ty dE dp Derivatives of E xo f uo ug p ty wrt p dE dx0 96 Derivatives of E xo f uo ug p tg wrt o dE du0 Derivatives of E zo f uo uf p t7 wrt uo dE dxf Derivatives of E xo x uo ug p ty wrt vy 20 dE duf 96 Derivatives of E xo f uo ug p tg wrt Uf
16. 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 Time Figure 12 Adjoint variables for the Optimal Mizing of a Catalyst 4 5 Example 5 Discrete time optimisation problem At each time instant measure the current state r and compute the control input sequence uz u 0 u 1 u N 1 R solving the following optimisation problem N 1 min z k Qz k u k Ru k NY Pz N k 0 subject to za k 1 a Kk 0 122 k g k 1 1 8z k 0 0787u k where x 0 0 3 0 1 N 20 R 1071 and Q is the identity matrix The parameters P and gm define a terminal invariant set where P is the positive definite solution of the algebraic Riccati equation for the selected Q and R Problem setup e The Optimal Control Problem is defined in the file Discrete_Sys m in the following way 44 Initial time fo lt For discrete time systems is the initial index problem time t0 0 Final time For discrete time systems tf_min tf_max problem time tf min problem time tf max guess tf 70 Parameter bounds pl p lt pu problem parameters pl problem parameters pu guess parameters 96 Initial conditions zo and its bounds problem states x0 0 3 0 1 problem states xOl problem states x0 problem states xOu problem states x01 Yo State bounds xl a t xu problem states xl inf infl problem states xu inf inf Yo Terminal state bounds xfl a t xfu problem states xfl
17. Control example 72 min tf subyect to di U 3 cos z 3 sin zi z 0 I2 4 0 73 tf z if 2 lt u 10 t 10 1 lt ts 100 t lt 2 10 a t 10 10 JE x 10 vt 0 tj x3 t Problem setup e The Optimal Control Problem is defined in the file SingularArc m in the following way Initial time to ty problem time t0 0 Yo Final time Let tf min tf max if ty is fixed problem time tf min 1 problem time tf 100 guess tf 4 70 Parameters bounds pl p lt pu problem parameters pl problem parameters pu guess parameters Initial condition z to and its bounds problem states xO pi 2 4 0 problem states xOl pi 2 4 0 problem states xOu pi 2 4 0 State bounds xl lt x t xu problem states xl 10 10 10 problem states xu 10 10 10 73 Terminal state bounds xfl x t lt xfu problem states xfl 10 O 0 problem states xfu 10 O 0 Yo Guess the state trajectories with x to x t guess states 1 pi 2 pi 2 guess states 2 4 0 guess states 3 0 01 96 Number of control actions N 96 Set problem inputs N 0 if N is equal to the number of integration steps problem inputs N 0 Input bounds ul lt u t uu problem inputs ul 2 problem inputs uu 2 Yo Guess the input sequences with u to u ty guess inputs 1 0 0 Choose the set points if required problem setp
18. Direct transcription settings e Number of integration nodes in the interval to tf The quantity steps N N num ber of control actions steps nodes 1 must be a positive integer 10 options nodes 1001 e Distribution of integration steps Set to tau 0 for equispaced steps Otherwise tau is a vector of length M 1 with 0 lt tau i lt 1 and sum tau 1 For discrete time system set tau 0 options tau 0 8 Multiple shooting settings e The ODE solver has to be set to cvodes options 0DEsolver cvodes e CVODES settings Refer to CVODES documentation for more details Method can be set to either Adams or BDF and solver to Newton or Functional options cvodes CVodeSet ptions RelTol 1 e 4 AbsTol 1 e 6 LinearSolver Dense MaxNumSteps 10000 LMM Method NonlinearSolver Solver options cvodesf CVodeSensSet ptions ErrControl true method Staggered 3 2 1 Derivative Definitions If some of the analytical gradients are supplied copy gradCost m jacConst m and hessianLagrangian m to the working directory and edit them as follows 1 Gradient of the cost In the file gradCost m it is possible to define the gradient of the stage cost function and the boundary cost by defining the partial derivatives dL dp dL dx dL du for the stage cost function and dE dtf dE dp dE dx0 dE duO dE dxf dE duf for the boundary cost function Whenever the gradient of the stage cost is
19. considers a continuous state constraint 7 1 min a t t 0 005u2 t dt ule 0 subject to ti xo Uu z 0 0 1 Problem setup e The Optimal Control Problem is defined in the file PathConstraint min the following way Initial time to ty problem time t0 0 Final time Let tf min tf max if tf is fixed problem time tf_min 1 problem time tf_max 1 guess tf 1 70 Parameters bounds pl p pu problem parameters pl problem parameters pu guess parameters Initial condition z t and its bounds problem states x0 0 1 problem states xOl O 1 problem states xOu O 1 State bounds xl z t xu problem states xl 10 10 problem states xu 10 10 56 Yo Terminal state bounds xfl z ty lt xfu problem states xfl 10 10 problem states xfu 10 10 Yo Guess the state trajectories with r to x t guess states 1 0 0 guess states 2 1 1 Number of control actions N Set problem inputs N 0 if N is equal to the number of integration steps problem inputs N 0 Input bounds ul lt u t lt uu problem inputs ul 20 problem inputs uu 20 Yo Guess the input sequences with u to u ty guess inputs 1 0 0 Choose the set points if required problem setpoints states problem setpoints inputs Bounds for path constraint function gl g x u p t lt gu problem constraints gl 0 10
20. its bounds problem states x0 0 1 0 problem states x01 0 1 0 problem states xOu O 1 0 36 96 State bounds xl z t lt xu problem states xl O 20 20 problem states xu 1 9 20 20 Yo Terminal state bounds xfl z ty xfu problem states xfl O 1 201 problem states xfu 0 1 20 Yo Guess the state trajectories with r to x t guess states 1 0 0 guess states 2 1 1 guess states 3 0 0 Number of control actions N If N is equal to the number of integration steps prob lem inputs N can be set to 0 problem inputs N 0 Input bounds ul lt u t lt uu problem inputs ul 10 problem inputs uu 10 Yo Guess the input sequences with u to u ty guess inputs 1 5 5 Choose the set points if required problem setpoints states problem setpoints inputs Bounds for path constraint function gl lt g x u p t lt gu problem constraints gl problem constraints gu Yo Bounds for boundary constraints bl lt b x to s 7 u to p to tr lt bu problem constraints bl problem constraints bu 96 store the necessary problem parameters used in the functions problem data problem functions L E Qf g b function stageCost L x xr u ur p t data 37 stageCost Oxt function boundaryCost E x0 xf u0 uf p tf data boundaryCost xf 3 function dx f x u p t data xl x 1 x2 x 2 ul u 1 d
21. state reference xr input u steady input reference ur parameters p and the time instant t The variables x xr u ur and p for a time instant are passed to the function as row vectors In general the arguments x xr u ur and p will be matrices whose rows correspond to the states inputs and parameters at different time instants For instance the it state variable xi can be obtained as follows xi x i Importantly the function should return a column vector if called with arguments that have more than one row Each entry of the output corresponds to the evaluation of the stage cost for a point in time If there is no stage cost let stageCost 0 t so that the output will have the right dimension 2 Boundary cost function The boundary cost function boundaryCost E x0 xf u0 uf p tf data returns a scalar cost as a function of its arguments If there is no boundary cost let boundaryCost 0 3 System equations The ODE right hand side of the dynamical system is defined in the function dx f x u p t data The input arguments follow the same rule as in the stage cost function The function returns a row containing the evaluation of each state equation for a given x u p for a point in time t The function should return a matrix if z u and p are matrices whose rows correspond to the states inputs and parameters at different points in time The adopted structure is shown here xi x 1 xnex n ul u 1 um u
22. u Hf 3 4 0 t 0 t 96 07 fi Ov Ou 0 f2 r u Hf 4 4 O t O t 96 O fi 0u 0 fo Ou URI III HL cell nfz nfz Hessian of L r u p t HL 1 1 0 t 0 L 0t nt 1 HL 1 2 0 t 9 L t z HL 2 2 0 t 197L z7 l HL 1 3 0 t 96 O L Ot Ox2 HL 2 3 0 t 96 9 L zi r HL 3 3 0 t 96 0 L 0x3 HL 1 4 0 t 96 O L HL 2 4 0 t 96 97L zi HL 3 4 0 t 96 97L z u HL 4 4 0 t O L 0u nE nt Ezezeros nE nE HE num2cell Ez Hessian of E ty Hg 1 96 Hessian of g x u p t Hb Hessian of b zo x p uo uf p tf If the Hessian of some component of the Lagrangian is not available set the corre sponding output term equal to the empty matrix For instance if the Hessian of the dynamical system is not available set Hf Each Hessian has to be specified in a two dimensional cell array of dimension given by the size of y defined as t p x ul or tf p o uo Uf vy or a subset of their variables Each entry 7 7 of the cell array corresponds to the derivative with respect to y i and y j The order to follow is given by the order of the vector variables inside 22 t p x ul or ty p zo uo uf zr If t or and p is not a decision variable the derivative with respect to time or and p must not be considered i e the Hessian is computed with respect to y p x u
23. uf ry depends only on the final state and then the variables containing the other derivatives are set to be empty matrices Each row of db dxf corresponds to the evaluated derivative of a terminal constraint with respect to zi 7 22087 va ty hessianlagrangian m function HL HE Hf Hg Hb hessianLagrangian X U P t E x0 xf u0 uf p tf data The Hessian of the different parts of the Lagrangian HL HE Hf Hg Hb must be supplied in cell arrays as follows nt np n m ng nb deal data sizes 1 6 nfz nt nptnt m Hf Lz zeros nfz nfz HL num2cell Lz HE He Bz zeros n n Hb num2ce11 Bz Hb 1 1 3 sqrt xf 1 xf 1 72 4 0l The Hessian of the dynamical system and the path constraints have not been specified and then the respective variables are set to be empty matrices Hf and Hg The Hessian of the stage cost is specified in a two dimensional cell array of dimen sion given by the size of y z u while the Hessian of the boundary constraints is given with respect to y vy Notice that dE flag 1 in gradCost m and does not depend on any variable so the Hessian of E is not evaluated and it has 53 Optimal state traiectories 152 Se mu zl E P x di m d ee t dem eps 4 E 0 5 i i i i J 2 d 3 E ie az p d ka 0 i s 0 0 5 1 1 5 2 2 5 3 Figure 16 State trajectories for minimum fuel orbit rai
24. variables options perturbation H and options perturbation J The perturbation size for sec ond derivatives can be set in options perturbation H The perturbation size for first derivatives can be set in options perturbation J It is possible to select default values for the perturbations by setting options perturbation H and options perturbation J to the empty matrix options perturbation H options perturbation J The default values for the Jacobian approximation is eps 2 1 3 while for the Hessian is 8 eps 1 3 5 NLP solver To choose the NLP solver the variable options NLPsolver can be set to either ipopt or fmincon For instance options NLPsolver ipopt If IPOPT has been chosen as the solver the following basic settings can be defined Desired convergence tolerance relative The default value is 1e 8 options ipopt tol 1e 9 Hessian computation can either be numeric to use the method selected in options derivatives or quasi newton for a limited memory quasi newton approximation options ipopt hessian approximation exact e Print level Check out the IPOPT documentation for details options ipopt print level b5 e Maximum number of iterations The default value is options ipopt max iter 3000 e Select the barrier parameter update strategy The default value for this string option is monotone Possible values monotone use the monotone Fiacco McCormick strategy adaptive use th
25. 000195 x 600 q 20 En 5 k 300 Te 0 38158 and T f 0 3947 The variable z is the product concentration z is the temperature and u is the coolant flow rate In the described optimisation problem the reactor undergoes a slow set point change from the stable steady state x 0 0 98 0 39 to the unstable steady state x t 10 26 0 65 62 Problem setup e The Optimal Control Problem is defined in the file RayHicksCSTR m in the following way lInitial time to tf problem time t0 0 Final time Let tf min tf max if ty is fixed problem time tf_min 10 problem time tf max inf guess tf 120 70 Parameters bounds pl p pu problem parameters pl problem parameters pu guess parameters Initial condition z t9 and its bounds problem states x0 0 9831 0 3918 problem states xOl problem states x0 problem states xOu problem states x01 State bounds xl lt x t xu problem states xl 0 01 problem states xu 1 1 Yo Terminal state bounds xfl lt z ty lt xfu problem states xfl 0 2632 0 6519 problem states xfu 0 2632 0 6519 Yo Guess the state trajectories with x to x t guess states 1 0 9831 0 2632 guess states 2 0 3918 0 6519 96 Number of control actions N Set problem inputs N 0 if N is equal to the number of integration steps problem inputs N 40 Input bounds ul lt u t uu problem inputs ul 0 problem inputs uu 2 Yo Guess the in
26. Imperial College London Optimal Control Software User Guide ICLOCS Paola Falugi Eric Kerrigan Eugene van Wyk Department of Electrical and Electronic Engineering Imperial College London London England UK iclocs imperial ac uk 6 May 2010 1 Introduction This document presents a brief user s guide to the optimal control software supplied The code allows users to define optimal control problems with general path and boundary constraints free or fixed final times and the ability to include constant design parameters as unknowns The following optimal control problems fall within the scope of the code mino pro J z u D ty subject to t f a t u t p t to zo Vt to trl gn g z t u t p t gu Vt to tr r X xo f U0 Us p tf v Vt to tr zy lt x t zu Vt to tr UL 2 u t lt uu Vte to tr PLS PS pu where uo u to x a ty and up u t Here the cost function is defined as JEO mt f Litt utt pdt E so 5 ap pty to where E is the cost associated with the boundary conditions and L the stage cost function The arguments over which the cost function can be minimised are the time varying control input signals u the initial state zo the final time t and a set of parameters p that are constant for the duration of the phase The function g describes general path constraints and imposes the boundary conditions at the beginning and
27. Optimal state trajectories Time Figure 31 State trajectories for Speyer s problem Optimal input sequences 0 4 0 5 0 6 Time 0 7 0 8 0 9 Figure 32 Input trajectory for Speyer s problem 80 Adioint variables 0 02 y SE 1 0 015 j N i 0 01 4 0 005 4 0 y 0 005r J 0 01 7 i K 0 0151 y LUN Mec 0 02 d F il il il il il il il il 0 0 1 02 038 04 05 06 07 08 039 1 Time Figure 33 Adjoint variables for Speyer s problem 4 12 Example 12 Two link Robot Arm Find 2 and u over t 0 l solving the following optimal control problem of a two link robot arm Section 12 4 2 Example 2 in 10 tr 0 01 ro uX dt min u tf subject to sin z3 cos 73 27 2x3 u1 3 cos x3 u2 T 3 sin z3 sin zs cos zs z2 521 cos ra un ua 31 z 2 33 a Sines T3 T2 Ly Ta z a 0 10 0 x t 10 0 0 5 0 522 1 lt u t lt 1 ty gt 0 1 81 Problem setup e The Optimal Control Problem is defined in the file TwoLinkRobotArm m in the following way Initial time to lt ty problem time t0 0 Final time Let tf min tf max if ty is fixed problem time tf_min 0 1 problem time tf max inf guess tf 3 1 70 Parameters bounds pl p pu problem parameters pl proble
28. a u p t Ox1 Ofi Ox1 t 9f 021 df dx 1 hi dhix1 X 1 U 1 500 X 4 h2 h1 0 47 dhix1 X 1 0 47 h2 1 2 0 029 X 3 0 00014X 3 O Lt 70 Of 3 u p t Ox3 E fi 022 wa 9f 023 df dx 2 0 Lt 0 01 U 1 500 X 4 O Lt O Lt 70 Of x u p t Ox3 Ofi Oxs aes Ofn Ox3 df dx 3 dhi1x3 X 1 dh2 X 1 dhix3 X 1 0 47 dh2 X 1 1 2 X 1 0 029 0 0001 0 0001 X 3 72 U 1 X 4 500 OxLtl 70 Of z u p t Ox4 Ofi Ox4 Of Ov4 df dx 4 UC 1 xX 1 500 X 4 72 UC 1 XC 2 500 X 4 72 U 1 X 4 72 1 X 3 500 O Lt df du 1 X 1 500 X 4 X 2 500 X 4 1 X 3 500 X 4 Lt 500 dg flag 0 db flag 0 The path constraints and boundary constraints are not present and then their deriva tives have not be supplied dg flag 0 and db flag 0 Instead the derivatives of f z u p t fi z u p t fu z u p t are supplied in the structured variable df as shown in the illustrated example The derivative of f r u p t with respect to x are stored in df dx which is a cell array where each entry corresponds to the derivative with respect to a state variable For instance the derivative of f x u p t with respect to x is stored in df dx i df dx i is a matrix with n columns and a number of rows given by the different evalu
29. aining the optimal solution for the inputs e solution T brings information about the evaluation time instants for solution X and solution U Remark 3 It is important to observe that the auxiliary data passed to may not change through the course of the IPOPT optimisation In order to store information changing over time the global variable so1 has been used Display output output m This function uses the display options in settings m and the solution generated by the call to solveNLP m in main m and plots the state input and adjoint variable trajectories If direct collocation methods are used output m can also estimate the relative local discretization error to estimate the accuracy with which the dynamics have been approximated 1 4 Examples The following sections contain all the optimal control examples included with the ICLOCS tool box 4 1 Example 1 A linear optimisation problem with bang bang con trol Find the final time ty and control input u R over t in 0 tr solving the following optimisation problem 9 17 x1 0 0 x2 0 0 zi ty 300 z tr uf lt 1 o 120 lt 10 e tz IA Problem setup e The Optimal Control Problem is defined in the file BangBang m in the following way problem time t0 0 Initial time to problem time tf_min 0 1 problem time tf_max 100 guess tf 1 Parameters bounds pl lt p lt pu problem parameters pl problem parame
30. and input reference trajectories The function transcribeOCP m generates constant references in to t expanding the set point values assigned in problem setpoints The user can overwrite the default references assigning other reference trajectories data references xr and data references ur are matrices with dimension M x n and M x m respectively If the transcription method is Hermite Simpson it is necessary to assign also the references at the interval midpoints M 2xoptions nodes 1 data sparsity contains sparsity information about the functions f L g and data tau contains information to evaluate the time instants in the interval to t see Remark l data map data map A data map B data map w and data map W contain the data used for the transcription with direct collocation methods while data map Vx data map xV data map Vu and data map uV are matrices for the mapping of the variables For the details see Remark 2 data FD contains information for the computation of the derivatives 14 e data acStruct brings the structure of the Jacobian for the constraints e data costStruct gives the structure of the Jacobian for the cost e data hessianStruct brings the structure of the Hessian for the Lagrangian e data multipliers In this field the initial value for the Lagrange multipliers can specified It is especially useful for warm starting the IPOPT solver They can be specified defining the fo
31. and ns 1 It is used to adjust sizes and values of some of the variables data x0 contains the initial condition for the state at the time fo If the vari able problem states x0 defined in myProblem n is empty transcribeOCP m set prob lem states x0 problem states x l data x0t contains the measured initial condition z ko The function transcribeOCP m sets data xOt data xO The variable data x0t needs to be updated to change the initial condition for receding horizon optimisation problems If the initial condition is fixed its bounds have to be updated in the following way infoNLP zl nt npti ntt nptn data x0t infoNLP zu nt npti nttnptn data x0t data cx0 transcribeOCP m sets data cx0 0 if problem states x0 defined in myProblem m is empty otherwise data cx0 1 When data cx0 1 the additional term Ao data x0t data x0 is introduced in the Lagrangian where Xo is the relative Lagrange multiplier data options contains information about the transcription method the derivative evaluation and the solver settings data functions contains the function definition of the stage cost boundary cost system dynamic path constraints and boundary constraints L E f g b deal data functions data data contains the data to be used inside the functions and defined in the variable problem data in the function myProblem m data references The variables data references ur and data references xr store respectively the state
32. ation times Each row of df dx i corre sponds to the evaluation of the derivative O f x u p t 0z Of x u p t Ox i at a given time instant They are expressed in term of the it state and input x and u evaluated at the time instants t to ty tf and taken as column vectors X i and U i 30 hessianLagrangian m function HL HE Hf Hg Hb hessianLagrangian X U P t E x0 xf u0 uf p tf data The Hessian of the different parts of the Lagrangian HL HE Hf Hg Hb must be supplied in cell arrays as follows nt np n m ng nb deal data sizes 1 6 nfz nt nptn m dhix1 0 11 0 006 X 3 0 006 X 1 X 3 72 96 h ri ra ri dhix3 0 11 0 006 X 1 0 006 X 1 4X 3 72 96 Ohi a1 x3 Ox3 Oho x3 Ox3 dh2 0 0055 0 0001 10 X 3 72 0 00014X 3 1 10 X 3 72 ha z3 0 3 d2h2 0 0055 2 30 X 3 0 0001 100 X 3 73 0 0001 0 0001 X 3 10 X 3 72 73 hi x1 03 0721 d2h1x1 2 0 11 0 00672 X 3 0 006 X 1 X 3 73 hi 11 23 021 0x3 d2hixix3 0 11 0 006 0 006 X 1 X 3 0 006 X 1 X 32 73 96 O hi z1 23 0 23 d2h x3 2 0 11 0 006 X 1 0 006 X 1 4X 3 3 Hf cell nfz nfz Hessian of f x u p t 07 fi 0x1 0 fo Oxj 0 fs Oxi O fa Oxy Hf 1 1 2 dhixi dOhixi O t 2 dhixi d2hix1 0 47 0 t 0 fi Ox10x2 0 f2 0102 0 f
33. atives of f z u p t fi z u p t fu x u p t are supplied in the structured variable df as shown in the illustrated example The derivative of f r u p t with respect to x are stored in df dx which is a cell array where each entry corresponds to the derivative with respect to a state variable For instance the derivative of f x u p t with respect to x is stored in df dx i df dxifi is a matrix with n columns and a number of rows given by the different evaluation times Each row of df dx i corre sponds to the evaluation of the derivative O f x u p t 0z Of x u p t Ox i at a given time instant The same rule holds for the other derivatives hessianLagrangian m function HL HE Hf Hg Hb hessianLagrangian X U P t E xO xf uO uf p tf data 21 The Hessian of the different parts of the Lagrangian must be supplied in cell arrays as follows nt np n m ng nb deal data sizes 1 4 nfz nt nptnt m Hf cell nfz nfz Hessian of f x u p t Hf 1 1 0 t 0 t 96 0 fi 0 0 f5 0t nt 1 Hf 1 2 0 t 0 t 96 O fi Ot zi 0 fo Ot Ox1 Hf 2 2 0 t 0 t 96 0 fi 9z2 0 f2 81 Hf 1 3 0 t 0 t 96 07 fi Ot zr 0 f2 Ot 8x3 Hf 2 3 0 t 0 t 96 O f zi zr 0 f3 Ox r Hf 3 3 0 t 0 t 96 19 fi Ox3 0 f2 82 Hf 1 4 0 t 0 t 96 07 fi t Ou 0 fo Ot Hf 2 4 0 t 0 t 96 O fi zi Ou 0 f2 0x4
34. cessary problem parameters used in the functions problem data problem functions L E Of g b function stageCost L x xr u ur p t data stageCost t 0 function boundaryCost E x0 xf u0 uf p tf data boundaryCost 1 xf 1 xf 2 function dx f x u p t data x1 x 1 x2 x 2 u u 1 dx 1 u 10 x2 x1 dx 2 u x1 10 x2 1 u x2 function c g x u p t data c function bc b x0 xf u0 uf p tf data be The solution method and solver settings are set in settings m See the file included in the directory ICLOCS examples CatalystMixing The solution of the optimisation problem is computed running the file main m inside the directory ICLOCS exzamples CatalystMixing The state control and adjoint variables solution to this problem using the Optimal Control Toolbox are shown in Figs 10 I1 and 42 0 9 0 8 0 7 0 6 0 5 0 4 0 3 0 2 0 1 Optimal state traiectories Figure 10 State trajectories for the Optimal Mixing of a Catalyst 0 3 0 4 Optimal input sequences 0 5 Time 0 6 0 97 0 8 0 7r 0 6 0 5r 0 4 0 3 0 2 pie Figure 11 Input trajectory for the Optimal Mixing of a Catalyst 0 1 0 2 0 3 0 4 0 5 Time 43 0 6 Adioint variables T t T 0 75 ee 7 TO 0 9 4 0 95 a 4 zi i 1 inque i 0 0 1 0
35. ction gl lt g x u p t lt gu problem constraints gl problem constraints gu Yo Bounds for boundary constraints bl lt b x to tf u to u t p to tr lt bu problem constraints bl O 0 problem constraints bu 0 0 96 Store the problem parameters used in the functions problem data Ti1 0 1405 problem data md 0 0749 problem functions L E Qf 6g 6b function stageCost L x xr u ur p t data 50 stageCost x 2 function boundaryCost E x0 xf u0 uf p tf data boundaryCost 0 function dx f x u p t data r x 1 q x 2 vex 3 phi u 1 Ti data T1 md data md dx q v v r r 2 Ti sin phi 1 md t q v r Ti cos phi 1 md t function c g x u p t data c function bc b x0 xf u0 uf p tf data bc xf 3 gqrt 1 xf 1 xf 2 1 e The solution method and solver settings are set in settings m See the file included in the directory ICLOCS examples OrbitRaising e The files gradCost m jacConst m and hessianLagrangian m for this example are sup plied See inside the directory ICLOCS examples OrbitRaising gradCost m function dL dE gradCost L X Xr U Ur P t E xO xf uO uf p tf data Lt ones size t dL Gradient of the stage cost L wrt t p 7 u dL flag 1 dL dp L1 dL dt dL dx O Lt Lt O Lt 51 dL du 0 Lt Yo dE Gradient of E with respect to tf p To uo Uf LF dE flag 1 dE dtf dE dp
36. e are no path constraints problem constraints gl gi lowerbound g2_lowerbound problem constraints gu gi upperbound g2_upperbound Bounds for boundary constraints Set the upper and lower bounds for the boundary constraint function as entries in a row vector if required Set the variables to if there are no boundary constraints problem constraints bl b1_lowerbound b2_lowerbound problem constraints bu b1_upperbound b2 upperbound 14 Function definition The stage cost function L the boundary cost function E the path constraints g system equations f and the boundary constraint function b have to be defined in func tions with the name L E g f and b respectively Their name has to be stored as follows and illustrated in the provided examples in the toolbox problem functions L GE f g b A detailed discussion of each function is carried out in the Section 3 1 2 15 User defined data problem It is possible to store constant parameters of the problem that are not optimisation variables for instance transition matrices reference trajectories cost weights etc in a structured variable as follows see also illustrative examples problem data a 2 problem data b 1 Set problem data if there are no data 3 1 2 Function Definitions 1 Stage cost function The stage cost function stageCost L x xr u ur p t data computes the stage cost for a given state x steady
37. e adaptive update strategy options ipopt mu strategy adaptive e Indicate which information for the Hessian of the Lagrangian function is used by the algorithm The default value is exact Possible values exact Use second derivatives provided by ICLOCS limited memory Perform a limited memory quasi Newton approximation implemented inside IPOPT options ipopt hessian_approximation exact There are many other options which are described in the IPOPT documentation These options can all be changed from the default settings in the file ICLOCS src solveNLP m if required If MATLAB s own NLP solver fmincon is used the options have to be set in options fmincon optimset Consult the MATLAB documentation for detailed information on the various options for fmincon 6 Output settings t is possible to use the function output m to display some information about the solved optimisation problem This function uses the display options defined in settings m The available options are described hereinafter Set to zero to disable options e Display computation time options print time 1 e Display relative local discretization error recommended for direct transcription options print relative local error 1 e Display optimal cost options print cost 1 e Plot states options plot states 1 e Plot inputs options plot inputs 1 e Plot Lagrange multipliers relative to the system equations options plot multipliers 1 7
38. e basis of the information stored in infoNLP and data and return the solution The function has two output arguments The first output solution is a structured variable containing the optimal final time parameters states controls and Lagrange multipliers Instead the second one status returns the exit condition of the solver Its values depend on the selected solver For its descrip tion see the help for the outputs of ipopt and fmincon In particular the output variable solution whenever IPOPT is employed returns e solution multipliers contains the Lagrange multipliers at the solution z see the IPOPT and MATLAB user guide for its structure solution computation_time Returns the CPU time in seconds that has been used by MATLAB to run the selected software e solution iterates gives the total number of iterations if fmincon is used The number of iterations for Ipopt can be obtained writing a function callback m see the example inside the folder src where its first input argment is stored in the global structured variable sol in sol iterates 16 solution z the computed solution z e solution tf the value of the final time solution p contains the computed values for the parameters p if any e solution X a matrix with dimension options nodes x n containing the optimal solution for the states solution x0 the initial value of the determined solution e solution U amatrix with dimension options nodes x mcont
39. end of the phase As a first step the user defined optimal control problem is transcribed to a static optimisation problem by either direct multiple shooting or direct collocation methods The direct multiple shooting formulation requires the solution of initial value problems that can be determined us ing the open source sensitivity solver package CVODES The entire Sundials suite includes the CVODES solver and the SUNDIALSTB MATLAB TOOLBOX The multiple shooting method does not yet support problems with the final time ty and the set of parameters p as variables The direct collocation formulations discretize the system dynamics using implicit Runge Kutta for mulae and can also be used to incorporate discrete time problems Once the optimal control problem has been transcribed it can be solved with a selection of nonlinear constrained optimi sation algorithms given by the open source code IPOPT or MATLAB s own NLP solver fmincon The derivatives of the ODE right hand side cost and constraint functions are also required for the optimisation and are either estimated numerically or supplied analytically 2 Installation The code is entirely MATLAB based and can be used without installing any additional software by making use of MATLAB s own built in functions The code can be used in conjunction with the following packages e POPT with mex interface highly recommended https projects coin or org Ipopt e SUNDIALS v 2 4 0 optional
40. he assignment of an initial guess for the variable z from guess defined on x and u can be easily performed by mean of data map Vx data map xV data map Vu and data map uV data map Vx maps z x data map xV maps x 2 data map Vu maps z u 15 data map uV maps u 2 Consider tf guess p guess x guess and u guess to be the guess for the final time if any the constant parameters if any the states and inputs respectively x guess and u guess are matrices with dimension M x n and N x m respectively where each row corresponds to the vector variables at some time instant It is possible to update the guess for z running the following lines u guess u guess u_guess u_guess x guess x guess x guess x guess infoNLP zO data map xV x guesstdata map uV u guess infoNLP zO 1 tf guess infoNLP zO 1 np p guess If the initial condition is fixed data xOt and its bounds have to be updated as previously explained e Once a solution z is available the variables x and u can be obtained in the following way x reshape data map Vx sol z n M u reshape data map Vu sol z m N U kron u ones M 1 N 1 x u and U are matrices of dimension M xn N xm and M 1 xm respectively Indeed the output of solveNLP m returns explicitly the variable x and u together with the optimisation variable z Function solveNLP m The function solveNLP m calls the selected solver to solve the optimisation problem on th
41. he gradient of the stage cost and the boundary cost are supplied and so dL flag 1 and dE flag 1 The final time ty and p are not decision variables and then the deriva tives dL dp dL dt dE dtf and dE dp are set to be empty matrices The derivative of the stage cost with respect to the state u is stored in a matrix with m columns and a number of rows given by the different evaluation times It is expressed in term of the input evaluated at different time instant The it input u evaluated at the time instants t to fk tf is a column vectors taken as U i The derivative of the stage cost with respect to the state z must be specified even if it is identically zero The derivative of the boundary cost depends only on the final state ry and then the variables containing the derivatives with respect to the other decision variables are set to be empty matrices jacConst m function df dg db jacConst f g X U P t b xO xf uO uf p tf tO data Lt ones size t df flag 1 hi h2 0 11 X 3 0 006 X 1 XC 3 5 96 hi zi 23 0 0055 X 3 0 00014X 3 1 10 X 3 ho x3 h a1 U3 O24 dhixi 0 11 0 006 X 3 0 006 X 1 X 32 72 Oh 21 3 Ox3 dh1x3 0 11 0 006 X 1 0 006 X 1 X 3 72 29 Yo h r3 x dh2 0 0055x 0 0001 10xX 3 72 0 0001 X 3 x 1410xX 3 72 df dp 1 L1 1 df dt 1 Of zx u p t dt Of
42. her discrete multiple_ shooting euler trapezoidal or hermite For instance options transcription trapezoidal Here discrete has to be chosen for discrete time systems The multiple shooting method can be used for continuous time systems whenever the final time and the constant parame ters are not decision variables of the optimisation problem In this case if IPOPT is used the quasi newton option options ipopt hessian_approximation limited memory for the Hessian computation has to be selected 2 Derivative generation This option selects how the derivatives are calculated The string derivative method can be set to the following values analytic or numeric options derivatives derivative_method For the analytic option the files gradCost m jacConst m and possibly hessianLagrangian m have to be defined see Section 3 2 1 For the numeric option the derivatives are com puted using finite differences and do not require definitions of any additional function If IPOPT is used the quasi newton option for the computation of the Hessian can be used 3 Whenever the numeric differentiation is enabled it is necessary to specify which kind of finite difference approximation to use between the following ones Central difference central forward difference forward For instance options hessianFD central 4 The perturbation size for numerical differentiation can be chosen by setting the
43. ied if required see Section 3 2 1 2 Define the optimal control problem by editing myProblem m see Section 5 1 Note that the file can be renamed as long as this change is reflected in the main m 3 Edit settings m to choose the solution method and solver settings see Section 3 2 4 Run main m Steps 2 and 3 are discussed in greater detail below 3 1 Defining the Optimal Control Problem The following section describes how the optimal control problem is to be defined in the problem definition file originally called myProblem m 3 1 1 General Problem Definition 1 Initial time The initial time fo has to be defined here problem time tO0 o For discrete time systems to is the initial index 2 Final time The final time can be a variable of the optimisation problem and the bounds for the final time have to be assigned problem time tf min final time min problem time tf max final time max guess tf final time guess If the final time is fixed set the minimum final time final time min equal to the maximum final time final time max Note that final time min fo and final time maz gt final time min For discrete time systems set the final time min final time max and final time guess as empty matrices In all other cases a final time guess has to be supplied 3 Parameters The bounds of any unknown constant parameters that are included in the optimisation should be defined here problem parameters pl
44. inal time Let tf min tf max if ty is fixed problem time tf_min 10 problem time tf_max 10 guess tf 10 70 Parameters bounds pl p pu problem parameters pl problem parameters pu guess parameters Initial condition z 0 and its bounds problem states xO 1 1 problem states xOl problem states x0 problem states xOu problem states x01 State bounds xl lt x t xu problem states xl 2 2 67 problem states xu 2 21 Yo Terminal state bounds xfl lt z f7 xfu problem states xfl O 0 problem states xfu 0 0 Yo Guess the state trajectories with r 0 x t guess states 1 1 0 guess states 2 1 0 96 Number of control actions N problem inputs N 10 Input bounds ul u t uu problem inputs ul 10 problem inputs uu 10 Yo Guess the input sequences with u 0 u ty guess inputs 1 2 0 96 Choose the set points if required problem setpoints states problem setpoints inputs Bounds for path constraint function gl lt g x u p t lt gu problem constraints gl problem constraints gu Yo Bounds for boundary constraints bl lt b x to t7 u to p to tr bu problem constraints bl problem constraints bu 96 Problem data problem data problem functions L E Of 6g 6b function stageCost L x xr u ur p t data xi x 1 x2 x 2 ui u 1 68 stageCost x1 x1 x2 x2tul ul function bounda
45. inf inf problem states xfu inf inf Yo Guess the state trajectories with r to x t guess states 1 problem states x0 1 01 guess states 2 problem states x0 2 0 Number of control actions N If N is equal to the number of integration steps prob lem inputs N can be set to 0 problem inputs N 0 Input bounds ul u t y lt uu problem inputs ul 2 problem inputs uu 2 Yo Guess the input sequences with u to wu ty guess inputs 1 2 0 Bounds for path constraint function gl lt g x u p t lt gu problem constraints gl 45 problem constraints gu Load information about the terminal set load termset Yo Bounds for boundary constraints bl lt b x to x t u to p to tr lt bu problem constraints bl 0 problem constraints bu gm 96 Choose the set points problem setpoints states 0 0 problem setpoints inputs 0 Rd 107 4 Qd 1 0 0 1 store the problem parameters used in the functions problem data P P problem data R Rd problem data Q Qd problem functions L E Of g b function stageCost L x xr u ur p t data RW data R QW data Q stagex x xr QW x xr stageu u ur RW u ur stageCost sum stagex 2 sum stageu 2 function boundaryCost E x0 xf u0 uf p tf data P data P boundaryCost xf P xf function dx f x u p t data xi x 1 x2 x 2 u u 1 46 MPC MPC
46. inputs ul ui_lowerbound lowerbound problem inputs uu ul_upperbound um_upperbound Note that there are m entries in the row vector if the problem is specified with m control inputs Guess input sequence Provide an initial guess for the optimal control sequence The initial guess is generated in a similar way to that of the initial state trajectory guess inputs 1 u1 t0 ui1 tf guess inputs m um t0 um tf If the variable guess inputs is the empty matrix the initial trajectories will be generated by linearly interpolating between random but feasible initial and final values for each state Note that the initial guess generated here can be overwritten and a user supplied guess can be supplied in main m by defining the variable infoNLP z0 see Remark D Choose set points Constant state and input setpoints can be defined here if required These will be formatted and passed to the stage cost function as xr and ur respectively to be used as reference trajectories along the optimisation horizon problem setpoints states x1 setpoint xn setpointl problem setpoints inputs lul setpoint um setpoint Alternatively time varying setpoints can also be passed to the stage cost through the struc tured variable data see Section B 3 Bounds for path constraint function Set the upper and lower bounds for the path constraint function as entries in a row vector if required Set the variables to if ther
47. l or y t x ul or y z u The same syntax has to be applied for y ty p zo uo us xf The Hessian with respect to the variable x and u must be specified If some variables do not contribute to the Hessian the corresponding entries can be set to zero otherwise the entries have to be vectors with the same length of t con taining the value of the Hessian at different time instants Notice that in the example the linear equations of dynamical system do not con tribute to the Hessian as a consequence all entries can be set to zero in the following way fz zeros nfz nfz Hf num2cell fz where nfz gives the dimension of the vector y t z ul The Hessian of E x xf uou p tr has to be considered on the basis of the specified analytic gradient when its evaluation is enabled dE flag 1 Since only dE dtf has been specified the other Jacobin are set to empty matrices in the file gradCost m the Hessian only with respect to t must be considered If dE flag 0 the Hessian with respect to the variables x0 u0 uf and xf must be always specified and it would be the following nE nt npt 2 n 2 m Ez zeros nE nE HE num2cell Ez Hessian of E zo zr tf Solution of the problem and results The solution of the optimisation problem is computed running the file main m The following lines are executed clear all format compact problem guess BangBang Fetch the problem definition options
48. l time t is fixed tf min tf max problem time tf_min 126 problem time tf_max 126 guess tf 126 Parameters bounds pl lt p lt pu problem parameters pl problem parameters pu guess parameters Initial condition zp and its bounds problem states x0 1 5 0 0 0 0 7 problem states x01 11 5 0 0 0 0 7 Bounds for zo 26 problem states xOu 1 5 0 0 0 0 71 State bounds zl x t x problem states xl 0 01 problem states xu 40 50 25 10 Yo Terminal state bounds zfl lt x t lt afu problem states xfl 0 0 01 problem states xfu 40 50 25 10 Yo Guess the state trajectories x to z t7 guess states 1 1 5 30 xi to xi ty guess states 2 0 8 5 ra to za tp guess states 3 0 0 za to va ty guess states 4 7 10 96 wa to c 27 Number of control actions problem inputs N 500 Input bounds ul lt u t lt uu problem inputs ul 0 problem inputs uu 50 Yo Guess the input sequences u to u t7 guess inputs 1 2 10 Choose the set points if required problem setpoints states problem setpoints inputs Bounds for path constraint function gl lt g x u p t lt gu problem constraints gl problem constraints gu Yo Bounds for boundary constraints bl lt b x to x tr u to p to tr bu problem constraints bl problem constraints bu 96 store the nece
49. liers lambda 2 ones 1 nc solution status solveNLP infoNLP data 96 Solve the problem output solution options data Output solutions The state control and adjoint variables solution to this problem using the Optimal Control Toolbox are shown in Figs 4 5 and 6 34 Optimal input sequences T 7 S T T Pox 40 35r 30r 25r 20r j 15 1 1 10r y Bp P d L 0 20 40 60 80 100 Time 120 Figure 5 Input trajectory for Fed batch fermentor problem Adjoint variables T ya 5 4 10 15 20 xu 25 J 1 1 1 1 0 20 40 60 80 100 Time Figure 6 Adjoint variables for Fed batch fermentor problem 35 4 3 Example 3 Bryson Denham optimal control problem Find the final time t and control input u IR over t in 0 t solving the following optimisation problem 2 min 223 1 i s ty subject to t 2 Lo 0 x1 t 1 9 20 lt z t lt 20 vt 0 ty 20 lt x3 tp lt 20 Problem setup e The Optimal Control Problem is defined in the file BrysonDenham m in the following way YInitial time to lt ty problem time t0 0 Final time Let tf min tf max if ty is fixed problem time tf_min 0 1 problem time tf_max 100 guess tf 1 Parameters bounds pl lt p lt pu problem parameters pl problem parameters pu guess parameters Initial condition zy and
50. llowing fields Multipliers corresponding to the lower bounds on the variables data multipliers zl Multipliers corresponding to the upper bounds on the variables data multipliers zu Multipliers corresponding to the constraints data multipliers lambda The fields data multipliers zl and data multipliers zu have to be a column vec tor with the same length of z The field data multipliers lambda have to be a col umn vector with the length given by the number of constraints size data jacStruct 1 Remark 1 The general formula to compute the vector t of the time instants ty to t for k 1 M is t data tf data t0 T data k0 where T is a vector taking values in the interval 0 1 such that t k data tf data t0 T k data k0 T depends on the distribution of integration steps stored in data tau ns and is given by T 0 cumsum data tau data Nm ns where ns 2ifthe transcription method is hermite and ns 1 otherwise so that sum data tau ns 1 for any transcription_ method The function transcribe0CP m set data k0 data tO initial time for the continuous time case while set data k0 initial time data t0 0 and data tf 1 for discrete time problems Then if the initial time has to be changed to solve a time varying op timisation problem in a receding horizon fashion the variable data k0 has to be set properly before to call the function solveNLP m Remark 2 The data mining of the state x and input u and t
51. m parameters pu guess parameters Initial condition z to and its bounds problem states x0 0 0 0 01 problem states x01 0 0 0 01 problem states xOu O 0 State bounds xl a t xu problem states xl inf inf inf inf problem states xu inf inf inf inf Terminal state bounds xfl x t lt xfu problem states xfl 0 0 5 0 522 problem states xfu 0 0 0 5 0 522 Yo Guess the state trajectories with r to x ty guess states 1 0 1 0 guess states 2 0 1 0 guess states 3 0 0 5 guess states 4 0 0 522 Number of control actions N Set problem inputs N 0 if N is equal to the number of integration steps problem inputs N 0 Input bounds ul lt u t lt uu problem inputs ul 1 1 82 problem inputs uu 1 1 Yo Guess the input sequences with u to u ty guess inputs Choose the set points if required problem setpoints states problem setpoints inputs Bounds for path constraint function gl lt g x u p t lt gu problem constraints gl problem constraints gu Yo Bounds for boundary constraints bl lt b x to x t u to p to tr bu problem constraints bl problem constraints bu 96 Store the necessary problem parameters used in the functions problem data problem functions 0L 0E f g b function stageCost L x xr u ur p t data stageCost 0 01 u 1 u 1 u 2 u 2 fu
52. md t x 0 1 0 J r t GA 1 0 ty 3 32 0 u t 27 0 5 x1 t 2 0 lt x t lt 1 Vt 0 tr with T 0 1405 and md 0 0749 Problem setup e The Optimal Control Problem is defined in the file OrbitRaising m in the following way Initial Time to ty problem time t0 0 Final time Let tf min tf max if tf is fixed problem time tf min 3 32 problem time tf max 3 32 guess tf 3 32 70 Parameters bounds pl p lt pu problem parameters pl problem parameters pu guess parameters Initial condition zp and its bounds problem states xO 1 O 1 problem states xOl 1 O 1 problem states xOu 1 O 1 49 State bounds xl a t x problem states xl 0 5 0 problem states xu 2 1 2 xu 0 51 13 Terminal state bounds xfl lt x t lt xfu problem states xfl O 0 problem states xfu 2 1 2 Yo Guess the state trajectories with r to x t guess states 1 1 1 4 guess states 2 0 0 8 guess states 3 1 0 05 Number of control actions N If N is equal to the number of integration steps prob lem inputs N can be set to 0 problem inputs N 0 Input bounds ul u t uu problem inputs ul 0 problem inputs uu 2 pi Yo Guess the input sequences with u to u ty guess inputs 1 0 2 pi Choose the set points if required problem setpoints states problem setpoints inputs Bounds for path constraint fun
53. nction boundaryCost E x0 xf u0 uf p tf data boundaryCost tf function dx f x u p t data x1 x 1 x2 x 2 x3 x 3 ul u 1 u2 u 2 dx 1 sin x3 x 9 4xcos x3 x1 72 2 x2 72 4 3x u1 u2 3 2 cos x3 u2 31 36 9 4 sin x3 72 dx 2 sin x3 9 4 cos x3 x2 72 7 2 x1 72 7 3xu2 3 2xcos x3 x u1 u2 31 36 9 4 sin x3 72 dx 3 x2 x1 83 Optimal state trajectories J S 0 8 Jg 0 6 d J y 0 4 1 i FAN oat ff 0 aed k i x ho 0 2E cw 7 0 4 b AV d A 0 6 X F 0 5 Figure 34 State trajectories for the two link robot arm problem dx 4 x1 function c g x u p t data c function bc b x0 xf u0 uf p tf data be e The solution method and solver settings are set in settings m the directory ICLOCS examples TwoLinkRobotArm See the file included in e The solution of the optimisation problem is computed running the file main m inside the directory ICLOCS eramples TuoLinkRobotArm The state control and adjoint variables solution to this problem using the Optimal Control Toolbox are shown in Figs B4 35 and B l References 1 John T Betts Practical Methods for Optimal Control Using Nonlinear Programming Ad vances in Design and Control SIAM Philadelphia PA 2001 84 Optimal input sequences T T T 0 6 1 i
54. ng cell array of appropriate dimensions otherwise setting the related variable to Refer to the directory ICLOCS examples for clear examples of how these analytic derivatives can be defined The derivatives that are not supplied the corresponding flag is set to zero or the Hessian is set to 11 are evaluated numerically Note that the limited memory quasi newton approximation of the Hessian can still be used if IPOPT is used 3 3 Solving the optimisation problem Once all the data for the definition and solution of the optimisation problem have been specified in myProblem m settings m and eventually in gradCost m jacConst m and hessianLagrangian m as described in the previous section the optimisation is performed running the following lines problem guess myProblem options settings infoNLP data transcribe0CP problem guess options solution status solveNLP infoNLP data The data inserted in myProblem m and settings m and contained in the variables problem guess and options have to be properly transcribed for the nonlinear solver The following subsections describes briefly the functions transcribe0CP m solveNLP m and their output ar guments Transcription function transcribe0CP m The function transcribe0CP m processes the information from problem guess and options for the nonlinear solver Mainly it defines bounds for the optimisation variable dynamic equa tions path and boundary constraints it format
55. nodes 1 For the trapezoidal method the required u M is imposed equal to u M 1 A detailed description of the output variable follows 1 Output variable infoNLP e infoNLP zl e infoNLP zu e infoNLP cl e infoNLP cu infoNLP zO Lower bound of the optimisation variable z Upper bound of the optimisation variable z Lower bound of the all set of constraints for the optimisation problem Upper bound of the all set of constraints for the optimisation problem Initial guess for the optimisation problem 2 Output variable data data t0 contains information to evaluate the time instants in the interval to tj see Remark l data k0 contains information to evaluate the time instants in the interval to tj see Remark l data Nm contains information to evaluate the time instants in the interval to tj see Remark 1 Nm 1 for continuous time systems Nm N for discrete time systems e data sizes contains information about dimensions involved in the problem nt np n m ng nb M N ns deal data sizes nt 1 if tf is a decision variable otherwise nt 0 np contains the number of free constant parameters n gives the number of states m gives the number of inputs ng gives the number of path constraints nb gives the number of boundary constraints N is the number of control actions M is the number of mesh nodes 13 ns 2 for Hermite Simpson transcription method
56. oints states problem setpoints inputs Bounds for path constraint function gl lt g x u p t lt gu problem constraints gl problem constraints gu Yo Bounds for boundary constraints bl lt b x to t7 u to p to tr lt bu problem constraints bl problem constraints bu 70 Problem data problem data problem functions L E Qf 6g 6b function stageCost L x xr u ur p t data stageCost Oxt function boundaryCost E x0 xf u0 uf p tf data 74 boundaryCost tf function dx f x u p t data x1 x 1 ul u 1 dx 1 ul dx 2 cos x1 dx 3 sin x1 function c g x u p t data c function bc b x0 xf u0 uf p tf data bc e The solution method and solver settings are set in settings m See the file included in the directory ICLOCS examples Singular Arc e The solution of the optimisation problem is computed running the file main m inside the directory ICLOCS examples SingularArc The state control and adjoint variables solution to this problem using the Optimal Control Toolbox are shown in Figs 28 29land 30 4 11 Example 11 Speyer s problem The Speyer s problem consists in the following periodic optimal control problem of a sailboat 111 75 3 5 2 5 0 5 0 8 0 6 0 4 0 2 Optimal state trajectories Figure 28 State trajectories for the singular arc problem Optimal input sequences
57. put sequences with u to u ty guess inputs 1 0 0 455 600 63 Bounds for path constraint function gl lt g x u p t lt gu problem constraints gl problem constraints gu Yo Bounds for boundary constraints bl lt b x to x t u to p to tr bu problem constraints bl problem constraints bu 96 Choose the set points problem setpoints states 0 2632 0 6519 problem setpoints inputs 455 600 96 store the necessary problem parameters used in the functions problem data problem functions L E Of g b function stageCost L x xr u ur p t data c x 1 T x 2 uzu 1 cr xr 1 Tr xr 2 stageCost 0 5 c cr c cr T Tr T Tr 0 5 u ur u ur function boundaryCost E x0 xf u0 uf p tf data boundaryCost 0 function dx f x u p t data Biegler s Model coefficient a 0 000195 600 q 20 En 5 k 300 Tc 0 38158 Tf 0 3947 c x 1 T x 2 uzu 1 dx 1 1 c q k c exp En T dx 2 T T qtk c exp En T ax u T Tc function c g x u p t data function bc b x0 xf u0 uf p tf data 64 0 7 0 6 0 5r QAL 0 3 Figure 22 c Time State trajectories for the continuous stirred tank reactor function bc b x0 xf u0 uf p tf data be e The solution method and solver settings are set in settings m See the file included in the directory ICLOCS examples RayHicksCSTR e The solu
58. ramming Chapman amp Hall CRC Monographs and Surveys in Pure and Applied Mathematics 2000 11 J L Speyer Periodic optimal flight Journal of Guidance Control and Dynamics 61 745 755 1996 12 J R Banga V S Vassiliadis and E Balsa Canto Second order sensitivities of general dy namic systems with application to optimal control problems Chemical Engineering Science 54 3851 3860 1999 86
59. roblem inputs uu 1000 Yo Guess the input sequences with u to u ty guess inputs 1 1 1 96 Choose the set points if required problem setpoints states problem setpoints inputs Bounds for path constraint function gl lt g x u p t lt gu problem constraints gl problem constraints gu 78 Yo Bounds for boundary constraints bl x b x to s 7 u to p to tr lt bu problem constraints bl O 0 problem constraints bu 0 0 96 Store the necessary problem parameters used in the functions problem data problem functions 10L 6E 6f 6g 6b function stageCosteL x xr u ur p t data xi x 1 x2 x 2 u u 1 a 1 stageCost 0 5 x1 72 a x2 72 x2 74 0 00001 u 72 function boundaryCost E x0 xf u0 uf p tf data boundaryCost 0 function dx f x u p t data xi x 1 x2 x 2 u u 1 dx 1 x2 dx 2 u function c g x u p t data c function bc b x0 xf u0 uf p tf data bc x0 xf e The solution method and solver settings are set in settings m See the file included in the directory ICLOCS examples Speyersproblem e The solution of the optimisation problem is computed running the file main m inside the di rectory ICLOCS examples Speyersproblem The state control and adjoint variables solution to this problem using the Optimal Control Toolbox are shown in Figs 32 and 33 79 0 6 0 4 0 2 150 100 50 100 150
60. rol Problem is defined in the file CatalystMixing min the following way Initial time to lt ty problem time t0 0 Final time Let tf min tf max if tf is fixed problem time tf min 1 40 problem time tf_max 1 guess tf 1 70 Parameters bounds pl p lt pu problem parameters pl problem parameters pu guess parameters Initial conditions zo and its bounds problem states x0 1 01 problem states xOl 1 0 problem states xOu 1 0 96 State bounds xl a t xu problem states xl 0 8 0 problem states xu 1 0 1 Yo Terminal state bounds xfl z ty xfu problem states xfl 0 8 0 problem states xfu 0 95 0 1 Yo Guess the state trajectories with x to x t guess states 1 1 0 915 guess states 2 0 0 05 Number of control actions N If N is equal to the number of integration steps prob lem inputs N can be set to 0 problem inputs N 0 Input bounds ul u t uu problem inputs ul 0 problem inputs uu 1 Yo Guess the input sequences with u to u ty guess inputs 1 1 0 96 Choose the set points if required problem setpoints states problem setpoints inputs Bounds for path constraint function gl lt g x u p t lt gu problem constraints gl problem constraints gu 41 Yo Bounds for boundary constraints bl x b z to s 7 u to U ty P to ty bu problem constraints bl problem constraints bu 96 store the ne
61. ryCost E x0 xf u0 uf p tf data boundaryCost 0 function dx f x u p t data xi x 1 x2 x 2 u u 1 dx 1 x2 dx 2 sin x1 u function c g x u p t data c function bc b x0 xf u0 uf p tf data be e The solution method and solver settings are set in settings m See the file included in the directory ICLOCS examples SimpleMPC e The solution of the optimisation problem is computed running the file mainMPC m inside the directory ICLOCS ezxamples SimpleMPC The following lines are executed clear all format compact problem guess testProblem Fetch the problem definition options settings Get options and solver settings plant testPlant Get function handle of plant model infoNLP data transcribeO0CP problem guess options Format for NLP solver nt np n m ng nb M N ns deal data sizes time states inputs P 20 Number of MPC iterations 69 Begin MPC loop for i 1 P disp Compute Control Action disp i solution solveNLP infoNLP data 96 Solve the NLP tc solution tf N Control horizon Apply control disp Apply Control x0 tv xv uv applyControl tc solution plant data i Store results time time tv states states xv inputs inputs uv Update initial condition infoNLP zl nttnp i nt4np tn x0 infoNLP zu nttnp i nt4np tn x0 data x0t x0 Update initial guess infoNLP zO solution z end End MPC loop
62. s 0z1022 0 fa 821023 Hf 1 2 0 t O t O t O t 0 fi 0x3 0 fo Ox3 O fs Ox3 0 fa Ox5 Hf 2 2 O t O t O t O t Yo 19 f Ox1 023 fe 0 04023 f5 0 0 023 8 f4 021023 Hf 1 3 dhix3 d2hixix3 X 1 dh2 dhix3 0 47 d2hixix3 X 1 0 47 dh2 1 2 0 029x0 0001 0 0001 X 3 72 0 xt1 107 0 zi za 0 fo Ox10x4 0 f3 02104 9 fa 021024 H 1 4 U 1 500 X 4 72 0 t 0 t 0 8t1 31 0 fi Or20x3 0 fa Ox20x3 0 fa Ox20x3 0 fa 8120x3 Hf 2 3 O t O t O t O t 107 fi Ox3 0 fo Ox3 O fa Ox3 0 fa Ox3 Hf 3 3 2 d2hi1x3 X 1 d2h2 X 1 d2hix3 X 1 0 47 d2h2 X 1 1 2 2 X 1 0 029 0 0001 0 0001 X 3 73 O t 07 fi Ov20x4 0 fo Ox20x4 0 fs 8120x4 8 fa Ox2024 HE 2 4 0 t UC 1 500 X 4 72 0 t 0 t 07 f Ovz0x4 0 fo Oxg0x4 0 f3 Ox30x4 0 fa Oxz0x4 HE 3 4 0 t O t UC 1 500 X 4 72 0 t 07 1 04 0 f2 0v 0 fs Oxg 0 fa zi HE 4 4 2 U 1 X 1 500 K 4 73 2XU 1 K 2 500 XC 4 73 2x0 1 XC 4 73 1 X 3 500 0 t 0 fi 0z10u 9 f2 0z10u 07 fz Ox19u 8 fa Ox10U Hf 1 5 1 500 X 4 O t O t O t 0 fi 0x30u 9 f2
63. s matrices for the direct transcription method if required it generates initial guesses for the optimisation variable and the structure of the Jacobian for the constraints and it constructs optimal finite difference perturbation sets The function has two output arguments that are structured variables The first output infoNLP con tains upper and lower bounds on the optimisation variable additional constraints and the initial guess Instead the second output variable data contains the data used in the functions evaluated during optimisation The optimisation variable z depends on the transcription method on the number of integration nodes options nodes M and on the number of control actions N If the transcription method is multiple shooting N M 1 and the optimisation variable is 12 The direct transcription methods present two different situations 1 If N M 1 the optimisation variable is in the most general case z tf p x 0 u 0 z 1 u 1 z M 1 M 1 z M 2 If N lt M 1 the optimisation variable is in the most general case es mant etos 4 4 1 u 0 2 AF 2 M 1 u N 2 M Notice that M 1 have to be divisible by the number of control actions N If tf and or p are not variables of the problem they are not introduced in z For discrete time systems N M 1 and ty cannot be a variable of the optimisation problem For the Hermite Simpson method it is imposed M 2xoptions
64. settings Get options and solver settings Format for the solver infoNLP data transcribe0CP problem guess options solution status solveNLP infoNLP data Solve the problem output solution options data Output solutions The state control and adjoint variables solution to this problem using the Optimal Control Toolbox are shown in Figs I 2 and 3 23 Optimal state traiectories 300 T T 250 150r 50 Figure 1 State trajectories for Bang Bang problem Optimal input sequences 1 T T 0 57 Time Figure 2 Input trajectory for Bang Bang problem 24 Adioint variables 1 T T T 0 8 0 6 iit 7 0 4 0 2 Me 4 Time Figure 3 Adjoint variables for Bang Bang problem 4 2 Example 2 Fed batch fermentor Find the control input u R over t in 0 ty solving the following optimisation problem 25 ts nn z t7 za t f 0 00001 u t de ut 0 subject to aa EN Ti med u hoax T 0 0129 il 2 500r4 5 21 21 0 029x3 u T3 h h 1 i an to mo 500 u t4 500 hadir m 0 006z x3 T3 ha cwm a 0 1 5 0 7 0 lt u t 50 lt 2 x vt 10 ty 0 za t 10 tf 126 Problem setup The Optimal Control Problem is defined in the file BatchFermentor m in the following way problem time t0 0 Initial time to Fina
65. sing problem not be specified If dE flag 0 whenever HE is not set to be empty the entries for the Hessian with respect to x to u to z tf are specified in the following way nt 0 np 0 nE nt npt 2 n 2 m Ez zeros nE nE HE num2cell Ez If db flag 1 it is necessary to check on which variables the boundary constraints b tr p x to u to u t depends on In this example it depends only on z and then it is necessary to specify only the Hessian with respect to z 5 If db flag 0 whenever Hb is not set to be empty it needs to specify at least the entries for the Hessian with respect to x to u to u t x t5 as follow here t and p are not de cision variables and they must not be considered Bzezeros nE nE Hb num2cell Bz Hb 1 6 3 sqrt xf 1 xf 1 72 4 01 e The solution of the optimisation problem is computed running the file main m inside the directory ICLOCS examples OrbitRaising The state control and adjoint variables solution to this problem using the Optimal Control Toolbox are shown in Figs 16 17 and L8 54 Optimal input sequences i 3 5 2 5 Figure 17 Input trajectory for minimum fuel orbit raising problem Adjoint variables 2 T T T Figure 18 Adjoint variables for minimum fuel orbit raising problem 55 4 7 Example 7 Path constrained optimal control problem The following optimisation problem
66. ssary problem parameters used in the functions problem data problem functions L E Qf 6g 6b 27 function stageCost L x xr u ur p t data stageCost 0 00001 u 1 u 1 function boundaryCost E x0 xf u0 uf p tf data boundaryCost xf 2 xf 4 function dx f x u p t data xi x 1 x2 x 2 x3 x 3 x4 x 4 ui u 1 h 20 11 x3 0 006 x1 4x3 h2 0 0055 x3 0 0001 x3 1 10 x3 dx 1 h1 x1 ul x1 500 x4 dx 2 h2 xx1 0 01xx2 u1 x x2 500 x4 dx 3 h1 x1 0 47 h2 x1 1 2 x1 0 029 x3 0 0001 x3 1 4 1 3 500 dx 4 u1 500 function c g x u p t data c function bc b x0 xf u0 uf p tf data bc e The solution method and solver settings are set in settings m See the file included in the directory ICLOCS examples BatchFermentor e The files gradCost m jacConst m and hessianLagrangian m for this example are sup plied See inside the directory ICLOCS examples BatchFermentor gradCost m function dL dE gradCost L X Xr U Ur P t E x0 xf u0 uf p tf data n deal data sizes 3 Lt ones size t 28 dL Gradient of the stage cost L wrt t p 7 u dL flag 1 dL dp L1 dL dt dL dx kron zeros 1 n Lt dL du 2 0 00001 U 1 Yo dE Gradient of E with respect to ty p zo uo Uf vf dE flag 1 dE dtf dE dp dE dx0 dE du0 L1 dE dxf O xf 4 0 xf 2 dE duf T
67. supplied set dL flag 1 otherwise set dL flag 0 Similarly if the gradient of the boundary cost is supplied set dE flag 1 otherwise set dE flag 0 The partial derivatives of the stage cost function must be expressed in vector form For instance the derivative of L with respect to the it state variable evaluated along all the horizon corresponds to the ith column of dL dx The same rule holds for dL du dL dp and dL dt For details refer to the sample files supplied in some of the examples 2 Jacobian of the constraint function In the file jacConst m it is possible to define the Jacobian of the ODE right hand side of the system equations path constraint and boundary constraint functions Whenever the gradient of the dynamics is supplied set df flag 1 otherwise set df flag 0 Similarly if the gradient of the path constraints is supplied set dg flag 1 otherwise set dg flag 0 The same rule holds for the boundary cost Set db flag 1 if the analytic expression of the gradient is available otherwise db flag 0 For details refer to the sample files supplied in some of the examples 11 3 Hessian of the Lagrangian The Hessian of the Lagrangian can be defined in the file hessianLagrangian m The toolbox allows to specify the Hessian of the following functions stage cost boundary cost path constraints boundary constraints are ODE right hand side of the system equations These can be supplied by defining the correspondi
68. ters pu guess parameters 96 Initial condition zo and its bounds problem states x0 0 01 x to problem states x01 0 01 Bounds for z fo problem states xOu 0 0 96 State bounds zl x t xu problem states xl 10 100 problem states xu 300 100 Yo Terminal state bounds v fl x t xfu problem states xfl 300 0 problem states xfu 300 0 Yo Guess the state trajectories with r to x t guess states 1 0 300 zi to x1 ty 18 guess states 2 0 01 96 z to va ty Number of control actions N If N is equal to the number of integration steps prob lem inputs N can be set to 0 problem inputs N 0 problem inputs ul 2 problem inputs uu 1 Yo Guess the input sequences with u to u ty guess inputs 1 2 1 Choose the set points if required problem setpoints states problem setpoints inputs Bounds for path constraint function gl lt g x u p t lt gu problem constraints gl problem constraints gu Yo Bounds for boundary constraints bl lt b x to tf u to p to tr lt bu problem constraints bl problem constraints bu 96 Store the necessary problem data used in the functions problem data problem functions L E Of g b function stageCost L x xr u ur p t data stageCost Oxt function boundaryCost E x0 xf u0 uf p tf data boundaryCost tf function dx f x u p t data x1 x 1
69. tion of the optimisation problem is computed running the file main m inside the di rectory ICLOCS ezamples RayHicksCSTR The state control and adjoint variables solution to this problem using the Optimal Control Toolbox are shown in Figs 23 and 241 4 9 Example 9 Continuous time Model Predictive Control example At any time instant t for k 0 1 2 control problem measure the state x t solve the following optimal 65 Optimal input sequences iE T T i T E 0 7 0 6 i i re 0 5 Tunel 0 4 el 0 3 i TE 1 1 1 0 20 40 60 80 100 120 140 Time Figure 23 Input trajectory for the continuous stirred tank reactor Adjoint variables 800 27 es 600 300 X 200 X 7 X 100r 2 4 0555550 r a 1 0 20 40 60 80 100 120 140 Time Figure 24 Adjoint variables for the continuous stirred tank reactor 66 10 min zi 0 s 0 u t dt u Jo subject to t 2 t sin z FU 10 ati 10 lt u t lt 10 LK and apply the first control action u 0 for the time window t t 1 The measures are collected from the plant described by the following ODE d 1 223 0 1sin t 9 0 2sin z1 u with z 0 1 1 Problem setup e The Optimal Control Problem is defined in the file testProblem m in the following way Initial time to lt ty problem time t0 0 F
70. x 1 x2 dx 2 ul dx 3 u1 u1 2 function c g x u p t data E Lis function bc b x0 xf u0 uf p tf data be e The solution method and solver settings are set in settings m See the file included in the directory ICLOCS examples BrysonDenham e The solution of the optimisation problem is computed running the file main m inside the directory ICLOCS ezamples BrysonDenham The state control and adjoint vari ables solution to this problem using the Optimal Control Toolbox are shown in Figs and D 4 4 Example 4 Optimal Mixing of a Catalyst This mixing problem considers a plug flow reactor packed with two catalysts involving the reactions 1 2 5 S3 The optimal mixing policy of the two catalysts has to be determined in order to maximise the production of species S3 12 The optimisation problem is formulated as follow 38 Optimal state traiectories 0 1 0 2 Time Figure 7 State trajectories for Bryson Denham problem Optimal input sequences i N N N Z 0 2 0 3 Time 0 4 Figure 8 Input trajectory for Bryson Denham problem 39 Adioint variables T eb 2 A 0 0 1 0 2 0 3 0 4 0 5 0 6 Time Figure 9 Adjoint variables for Bryson Denham problem 77 zi t x2 tf subject to T 1022 t2 u zi 1022 1 u zs 0 Il 0 iy 1 Problem setup e The Optimal Cont
71. y z ul or y xs If some variables do not contribute to the Hessian the corresponding entries can be set to zero otherwise the entries have to be vectors with the same length of t containing the value of the Hessian at different time instants For instance in the example the Hessian of the stage cost can be specified in a com pact way as follows Lz zeros nfz nfz Hf num2cell Lz HL 5 5 2 0 00001 Lt where nfz gives the dimension of the vector y z ul If dE f1ag 0 the Hessian with respect to the vector y r0 u0 uf x f must be spec ified as follows nE nt np 2 n42 m Ez zeros nE nE HE num2cell Ez Hessian of E zo vf uo ug p tg HE 8 nE 1 96 O E zo rf uo ug p tp Ova tf Ova tf 33 Optimal state traiectories 30 25 20r Figure 4 State trajectories for Fed batch fermentor problem Solution of the problem and results The solution of the optimisation problem is computed running the file main m The folloving lines are executed clear all format compact problem guess BangBang Fetch the problem definition options settings Get options and solver settings Format for the solver infoNLP data transcribe0CP problem guess options Initialize the dual point nc size data jacStruct 1 nt np n m M N deal data sizes 1 4 7 8 nz nt nptn M mN data multipliers zl 2 ones 1 nz data multipliers zu 2 ones 1 nz data multip

Imperial College London Optimal Control Software User Guide

Contents

Download Pdf Manuals

Related Search

Related Contents