TOMLAB
Setting up a TOMLAB /SOCS optimal control problem involves creating a MATLAB structure with static information for the problem, and writing MATLAB m-files (callback functions) to define the functions in the problem.
The functions that can be defined through m-files are:
If the problem consists of functions depending linearly of any of the variables, this can be handled without writing m-files for all functions except the differential equation right-hand-side functions. Treating the linear term as a special case can improve the performance of the solution process. The coefficients of the linear terms are stored in the problem structure.
The mathematical problem definition given in section 2 does not illustrate how the problem is defined when using the software. A problem definition that directly connects to the problem structure and the callback functions follows:
Minimize
|
where N is the number of phases, yI(p) = [tI, x(p)(tI), q(p)]T , yF(p) = [tF, x(p)(tF), q(p)]T , δI(p) and δF(p) are real-valued vectors, while ΦI(p) , ΦF(p) and L(p) are arbitrary functions.
Subject to, for each phase p, with its set of state variables x, control variables u, parameter variables q and values of tI, tF the differential equations
| ẋ = f(x(t), u(t), t, q), t = [tI, tF], |
event (boundary) constraints
| ψlj ≤ ψIj(wI) + ψFj(wF) + γIjT wI + γFjT wF ≤ ψuj, j = 1,…,ne, |
and path constraints
| glk ≤ gk(x(t), u(t), t, q) + αkT v(t) + βkT ak(x(t), u(t), t, q) ≤ guk, |
| t = [tI, tF], k = 1,…,ng |
where v = [x(t), u(t), q]T , wI = [tI, x(tI), u(tI), q]T, wF = [tF, x(tF), u(tF), q]T , and αk, βk and γ are real-valued vectors. gk are arbitrary single-valued functions, while ak is a vector-valued function (auxiliary functions) defined for each constraint.
There are also link constraints, used to link quantities between phases. One can link any quantity (x, u, t, q) at any endpoint (initial or final) of any phase with any other quantity at any endpoint of any phase.
Simple bounds can and should be set on all variables.
If the bounds for initial and/or final time of a phase are not fixed, that is, the initial and/or final time are variables, then they will be available as parameters in the callback functions (m-files). They will be added last in the parameter vector, tI before tF if both are variables. So parameters will become either a 1 or 2 longer vector depending whether both or only one of the times are variables.
np (mod) is the number of parameters INCLUDING the free time values.
All static data of a problem is defined in a MATLAB problem structure, containing:
All fields in the problem structure needs to be set prior to calling the solver. This is easiest done automatically using the utility newSocsProblem.m:
problem = newSocsProblem();
The structure is as follows:
Table 1: TOMLAB /SOCS problem structure
problem structure costfun_endpoint string costfun_integrand string odefun string eventfun string pathfun string pathfun_auxiliary string phases structure array, N nodes double (integer values) lower structure time double matrix, 1 x 2 states double matrix, ns x 3 control double matrix, nc x 1 parameters double matrix, np x 1 events double matrix, necon x 1 path double matrix, npcon x 1 phaselength double scalar upper structure time double matrix, 1 x 2 states double matrix, ns x 3 control double matrix, nc x 1 parameters double matrix, np x 1 events double matrix, necon x 1 path double matrix, npcon x 1 phaselength double scalar A double matrix, npcon x ns+nc+np nBterms double matrix (integer values), npcon x 1 B double matrix, sum(nBterms) x 1 C_I double matrix, necon x 1+ns+nc+np C_F double matrix, necon x 1+ns+nc+np D_I double matrix, N x 1+ns+np D_F double matrix, N x 1+ns+np cost_pattern double matrix, 1 x 3 event_pattern double matrix, necon x 2 path_pattern double matrix, npcon x 1 scaling double matrix, 1+ns+nc+np x 2 guess structure array, N time double matrix, 1 x s states double matrix, ns x s controls double matrix, nc x s parameters double matrix, np x 1 links double matrix, nl x 6 options structure OC structure NLP structure PrintFile string ErrorFile string user any type
where
N Number of phases of the problem. Defined by the length of the phases array. nl Number of linkage constraints in the problem. Defined by the dimensions of links.
The rest of the variables used in the dimensions have phase specific values. There exists N number of sets of [ns, nc, np, necon, npcon, s ]:
ns Number of state variables in the phase. Defined by the dimensions of lower.states. nc Number of control variables in the phase. Defined by the dimensions of lower.controls. np Number of parameter variables in the phase. Defined by the dimensions of lower.parameters. necon Number of event constraints in the phase. Defined by the dimensions of lower.events. npcon Number of path constraints in the phase. Defined by the dimensions of lower.path. s Number of start guessing points. Defined by the length of guess(i).time for phase i.
Name of callback m-file that defines the endpoint functions, ΦI and ΦF, in the objective function. This m-file defines both the initial and final endpoint functions for all phases.
See more information about this callback in the callback section 3.4.
Name of callback m-file that defines the integrand functions, L, in the objective function. This m-file defines the integrand for all phases.
See more information about this callback in the callback section 3.4.
Name of callback m-file that defines the right hand side of the differential equation. This m-file defines the right hand side for all phases and all differental equations.
See more information about this callback in the callback section 3.4.
Name of callback m-file that defines the event constraints. This m-file defines both the initial and final endpoint functions of the event constraints for all phases.
See more information about this callback in the callback section 3.4.
Name of callback m-file that defines the functions gk of the path constraints. This m-file defines the functions gk for all the path constraints that are selected to have a g function, for all phases. Additional nonlinear terms can be added to a path constraints; the auxiliary functions, which are defined separately in pathfun_auxiliary.
See more information about this callback in the callback section 3.4.
Name of callback m-file that defines the auxiliary functions a of the path constraints. There can be several auxiliary function tied to one path constraint, and there can be different number of auxiliary functions tied to different constraints. This function should compute and return all auxiliary functions for a given phase, i.e. for all constraints at one time. This m-file defined the auxiliary functions for all constraints phases.
See more information about this callback in the callback section 3.4.
This is an array of N structures, where N is the number of phases in the problem. The data in a structure i in the array is related only to phase i. Information that is supplied by problem.phases is:
Must be an integer scalar given as a double. Gives the initial number of grid points. The grid will be refined through mesh refinement if it doesn't give sufficient precision. Note that the number of mesh refinements is limited by OC option MITODE.
If given, and a start guess is also given, s (the number of guess points) must be equal to nodes. An exception to this is when s is 2. Then the start guess will be linearly interpolated between the two guess points at nodes number of points (including the two given guess points).
Set nodes empty to let nodes take the value of s, if a start guess is given. If no start guess is given, nodes must not be empty.
See problem.guess (3.3.26).
The lower bounds for variables.
A 1 x 2 matrix giving the lower bounds of the initial and final time of the phase.
| lower.time(1) | lower bound for tI (initial time) |
| lower.time(2) | lower bound for tF (final time) |
An ns x 3 matrix giving the lower bounds of the state variables for initial phase time, during phase and for final phase time.
lower.states(i,1) lower bound for xi(tI) lower.states(i,2) lower bound for xi(t), tI ≤ t ≤ tF lower.states(i,3) lower bound for xi(tF)
NOTICE: The number of state variables does not need to be equal in all phases. The state variables in one phase are totally independent of the state variables in another phase. Two variables in different phases can be made dependent of each other through linkage constraints in the links matrix.
The number of state variables, ns, in the phase is determined from the dimensions of this matrix.
An nc x 1 matrix giving the lower bounds of the control variables during phase.
lower.controls(i) lower bound for ui(t), tI <= t <= tF
NOTICE: The number of control variables does not need to be equal in all phases. The control variables in one phase are totally independent of the control variables in another phase. Two variables in different phases can be made dependent of each other through linkage constraints in the links matrix.
The number of control variables, nc, in the phase is determined from the dimensions of this matrix.
An np x 1 matrix giving the lower bounds of the parameter variables.
lower.parameters(i) lower bound for qi
NOTICE: The number of parameter variables does not need to be equal in all phases. The parameter variables in one phase are totally independent of the control variables in another phase. Two variables in different phases can be made dependent of each other through linkage constraints in the links matrix.
The number of parameters, np, in the phase is determined from the dimensions of this matrix.
An necon x 1 matrix giving the lower bounds of the event constraints.
lower.events(i) lower bound for event constraint i, ψli
The number of event constraints, necon, in the phase is determined from the dimensions of this matrix.
An npcon x 1 matrix giving the lower bounds of the path constarints.
lower.path(i) lower bound for path constraint i, gli
The number of path constraints, npcon, in the phase is determined from the dimensions of this matrix.
A scalar giving the lower bound on the length of the path constraint, tF − tI, that is. phaselength may be negative.
lower.phaselenth(i) lower bound for phase length, tF − tI
The upper structure is like the lower structure, but change all occurrences of lower to upper in the text above.
An npcon x ns +nc +np matrix giving the coefficients of the linear terms in the path constraints. Each row represents one constraint, and the corresponding αk vector. The scalar product of αk and v = [x, u, q] is added to constraint k.
Let A be empty if there are no linear terms in any constraint.
It is recommended to let A be sparse.
EXAMPLE:
We have three states x(1−3), two controls u(1−2) and two parameters q(1−2). We want two linear path constraints, one nonlinear with linear terms, and one nonlinear without linear terms:
|
Then A would be the matrix:
A = [2 0 1 0 -10 1 0;
2 3 0 -1 0 3 0;
0 0 0 0 0 0 0;
0 0 0 -1 -2 0 0];
Off topic: A path_pattern needs to be set to handle the nonlinear functions g3 and g4.
path_pattern = [0;
0;
1;
1];
and g should be defined in callback pathfun.
nBterms and B give information about the auxiliary functions in the path constraints. One path constraint can contain an arbitrary number of auxiliary function, and no path constraints need to have the same number of auxiliary functions. Some may have no auxiliary functions, others do.
nBterms is a column vector of length npcon, i.e. the number of path constraints in the phase. nBterms(i) tells how many auxiliary functions there are for constraint i.
B is a vector of length sum(nBterms). It contains the coefficients
for each auxiliary function. Coefficient for auxiliary function j
in path constraint k is B(sum(nBterms(1:(k-1)))+j).
The auxiliary functions are defined in pathfun_auxiliary in the same order as the corresponding coefficient values in B.
EXAMPLE:
We have two path constraints. The first one has 2 auxiliary functions, while the second one has 4:
|
Then
nBterms = [2 4]';
B = [4 -2 1 -10 2 1]';
An necon x 1+ns +nc +np matrix giving the coefficients of the linear terms in the event constraints for the initial (I) and final (F) endpoints of the phase. Each row i in C_I/F corresponds to γ(I/F)i. The scalar product of γ(I/F)i and w(I/F) = [t(I/F), x(t(I/F)), u(t(I/F)), q] is added to constraint i.
Let C_I/F be empty if there are no linear terms in any event constraints.
It is recommended to let C_I/F be sparse.
EXAMPLE:
We have two states, x(1−2), one control, u1, and one parameter q1. We want to formulate the following four event constraints:
|
Then C_I and C_F will look like this:
C_I = [1 0 0 -1;
0 0 1 0;
0 0 1 1;
0 1 0 0];
C_F = [-1 0 0 0;
0 -1 1 0;
0 0 0 0;
0 1 0 0];
Notice that the parameter q1 does not depend on t, so one can give the coefficients of the parameters in any of C_I and C_F. In the example above, C_I was chosen to carry the coefficients of the parameters.
An N x 1+ns +np matrix giving the coefficients of the linear terms in the objective endpoint functions for the initial (I) and final (F) endpoints of the phase. Each row p in D_I/F corresponds to δ(I/F)(p). The scalar product of δ(I/F)(p) and y(I/F)(p) = [t(I/F), x(p)(t(I/F)), q(p)] is added to the objective function.
Let D_I/F be empty if there are no linear terms in the objective function.
It is recommended to let D_I/F be sparse.
See problem.phases.C_I/F above for an example. D_I and D_F works the same, but instead of the vector w(I/F), there is this vector y(I/F)(p), which does not hold the control variables.
A 1 x 3 matrix describing which endpoint and integrand functions are defined for the given phase.
If cost_pattern(1) is set to nonzero, an initial endpoint function is defined in the objective function for the given phase, otherwise not.
If cost_pattern(2) is set to nonzero, a final endpoint function is defined in the objective function for the given phase, otherwise not.
If cost_pattern(3) is set to nonzero, an integrand function is defined in the objective function for the given phase, otherwise not.
If a function is defined, TOMLAB /SOCS will call it and expect it to return a function value. cost_pattern(1:2) is related to the callback m-file given in costfun_endpoint, while cost_pattern(3) is related to costfun_integrand.
An undefined function is equal to zero.
Leaving cost_pattern empty automatically makes the pattern full. It might be handy, but might also waste performance and cause errors if the functions in reality are not defined in the callbacks.
An necon x 2 matrix describing which endpoint functions are defined in the event constraints.
If event_pattern(i, 1) is set to nonzero, an initial endpoint function is defined in event constraint i for the given phase, otherwise not.
If event_pattern(i, 2) is set to nonzero, a final endpoint function is defined in event constraint i for the given phase, otherwise not.
If a function is defined, TOMLAB /SOCS will call it and expect it to return a function value. event_pattern(1:2) is related to the callback m-file given in eventfun.
An undefined function is equal to zero.
Leaving event_pattern empty automatically makes the pattern full. It might be handy, but might also waste performance and cause errors if the functions in reality are not defined in the callbacks.
An npcon x 1 matrix describing which g functions are defined in the pathfun callback m-file.
If path_pattern(i) is set to nonzero, a g function is defined in path constraint i for the given phase, otherwise not.
If a function is defined, TOMLAB /SOCS will call it and expect it to return a function value. path_pattern is related to the callback m-file given in pathfun.
An undefined function is equal to zero.
Leaving path_pattern empty automatically makes the pattern full. It might be handy, but might also waste performance and cause errors if the functions in reality are not defined in the callbacks.
A 1+ns +nc +np x 2 matrix containing scale information. Each row corresponds to one variable. First row is the independent variable, time. The ns following rows are the state variables. The nc following are the control variables. And the np last rows are the parameter variables. Column two defines the variable scales for the NLP problem, overriding the automatic scaling procedure.
Leave empty for no scaling, which is equivalent to problem.phases.scaling is all ones.
This is an array of N structures, where N is the number of phases in the problem. The data in a structure i in the array is related only to the start guess of phase i. Information that is supplied by guess is:
Leave empty to generate a linearly interpolated start guess from the mean of lower and upper bounds given in phases for all phases.
A 1 x s matrix giving the time points at which a guess of the states and the controls is given. The number s needs to be equal to problem.phases.nodes for the corresponding phase or equal to 2. See problem.phases.nodes (3.3.8) for more info.
Leave empty to generate a linearly interpolated start guess from the lower and upper bounds given in phases for the given phase.
An ns x s matrix giving the guess of the state variables at the time points given in guess.time.
An nc x s matrix giving the guess of the control variables at the time points given in guess.time.
An np x 1 matrix giving the guess of the parameter variables.
An nl x 6 matrix where the linkage constraints of the problem are defined. nl is determined from the dimensions of this matrix.
One linkage constraint can connect any two variables (time, state, control or parameter) in either the beginning or the end of any phase.
For linkage constraint i, the following applies:
The first variable to link is given in links(i,2). The first variable belongs to phase links(i,1). The variable is given as a number j from 0 to ns +nc +np, where:
j = 0 means time 0 < j ≤ ns means state variable j ns < j ≤ ns+nc means control variable j-ns ns+nc < j ≤ ns+nc+np means parameter variable j-ns-nc
links(i,3) is equal to
-1 if we want to link the value of the first variable in the beginning of the phase to the second variable 1 if we want to link the value of the first variable in the end of the phase to the second variable.
The second variable to link with the first variable is given like the first variable in links(i,4:6) instead of links(i,1:3).
EXAMPLE:
We have two phases:
phase 1 two states, one control, variable ending time phase 2 one state and one parameter, variable starting time
We want to make three links:
(1) link state variable two in the end of phase 1 with the state variable in the beginning phase 2, (2) link the control variable in the end of phase 1 with the parameter variable in the beginning of phase 2, (3) link the ending time of phase 1 with the starting time of phase 2.
The links matrix would look like this:
links = [1 2 1 2 1 -1; % link (1)
1 3 1 2 2 -1; % link (2)
1 0 1 2 0 -1]; % link (3)
Even if it suggests so, this example does not say one can not connect the beginning of two variables, or different variable types with each other. Sometimes it can be practical to for example link two variable ending times together if a process is constrained to stop at the same time.
A structure with options information in it. It provides:
A structure containing name-value pairs for the options to set to the OC (optmal control) solver. One field makes one pair, with the field name as name and the field value as value. A list of available options and descriptions is available in section 4.
EXAMPLE:
options.OC.IPGRD = 20; % Sets the output level to 20.
A structure like options.OC, but with options for the NLP solver.
Name of print file where output from SOCS (OC + NLP) is printed. If not
given, or left empty, a default name: socsprint.txt will be
used.
A flag with two possible values: 0 or 1. If set to 0, no print file is created, regardless of other output flags. If set to 1, a print file is created and written according to the other output flags.
Name of print file where error output from SOCS (OC) is printed. Look
in this file for error messages, if the solution encountered errors. If
not given, or left empty, a default name: socserror.txt will be
used.
A flag with two possible values: 0 or 1. If set to 0, no error file is created, even if there are errors. If set to 1, an error file is created and written.
Discretization/integration method, used in conjunction with options.nstg (discretization stages). Valid combinations are:
method nstg Discretization 0 1 ANL: Analytic Transformation 1 1 EUL: Euler's Method 1 2 RK2: Runge-Kutta 2-stage 1 3 RK3: Runge-Kutta 3-stage 1 4 RK4: Runge-Kutta 4-stage 2 1 TRP: Trapezoidal 2 2 HSS: Hermite-Simpson (Separated) 3 1 HSC: Hermite-Simpson (Compressed) 4 1 LM1: Linear Multistep 1-step 4 2 LM2: Linear Multistep 2-step 5 |1| RKT: Runge-Kutta-Taylor Variable Step 5 |2| DOP: Dormand-Prince Variable Step 5 |3| CWG: Gear Stiff BDF, Variable Order and Step 5 |4| ADM: Adams Predictor-Corrector, Variable Order and Step
When method =5, if:
nstg > 0 Forward propagation (from the beginning of the phase to th end) nstg < 0 Backward propagation (from phase end to phase start)
Default: TRP: Trapezoidal (method = 2, nstg = 1)
Give method and nstg as scalars, or as vectors of length N (number of phases). In case of scalars, the same method is used for all phases. In case of vectors, (method(i), nstg(i)) is used for phase i.
This variable is passed unchanged to the callbacks as the last argument.
A callback function is implemented in an m-file and is the way to provide nonlinear functions to the problem. There are six callback functions that can be defined. Some callbacks have patterns. A pattern describes exactly which functions are defined in the callback function in question. In the table below, the callbacks are listed with the names of their respective patterns, and a description:
Callback generic name Pattern Description costfun_endpoint cost_pattern Defines ΦI and ΦF costfun_integrand cost_pattern Defines L odefun no pattern Defines f eventfun event_pattern Defines ψI and ψF pathfun path_pattern Defines g pathfun_auxiliary no pattern Defines a
The patterns describe the mathematical functions defined in the callbacks. This is useful when there are functions consisting of only linear terms, and the callback functions need not be called. If a problem consists of only linear terms, one can omit the callback by not giving an m-file name. In some cases though, there are mixed nonlinear and linear functions in a problem. Then the patterns are useful.
More information about patterns can be found in Sections 3.3.22, 3.3.23 and 3.3.24 and also below in the descriptions of the callbacks.
The name of the costfun_endpoint callback m-file is given in problem.costfun_endpoint (3.3.1).
Purpose
The purpose of the user supplied costfun_endpoint callback is to
compute the endpoint functions (ΦI and ΦF) of the objective
function for all phases. Exactly which functions to compute is determined by
cost_pattern.
Calling Syntax
[phi, iferr] = costfun_endpoint(phase_number, endpoint, states,
time, parameters, user)
Description of Input
| phase_number | Number of the phase, 1:N, for which to
evaluate the endpoint function. Length: 1 |
| endpoint | Flag telling whether to compute the initial
endpoint function or the final endpoint
function. If endpoint =-1, compute the
initial endpoint function, if endpoint =1,
compute the final endpoint function. 1 and
-1 are the only possible values of
endpoint. Length: 1 |
| states | The state variable values at time time. Length: ns |
| time | The time value. Is equal to tI or tF
depending on the value of endpoint. Length: 1 |
| parameters | The parameter variable values. Length: np |
| user | The user variable from problem.user. |
Description of Output
| phi | The computed function value. Length: 1 |
| iferr | A user defined error flag. Set it to 0 to indicate no error. SOCS will terminate if it is nonzero. |
cost_pattern(1:2) applies to this callback. A nonzero value in
cost_pattern(1) indicates the initial phi functions is defined in
costfun_endpoint and can be evaluated by calling it with endpoint
= -1. A zero value indicates it is not defined at all, and is
equal to zero. This function will then not be called with endpoint
= -1. It is analogues with cost_pattern(2), final phi functions
and endpoint = 1.
The name of the costfun_integrand callback m-file is given in problem.costfun_integrand (3.3.2).
Purpose
The purpose of the user supplied costfun_integrand callback is to
compute the integrand functions (L) of the objective
function for all phases. Use cost_pattern to tell what phases
have an integrand function defined.
Calling Syntax
[L, iferr] = costfun_integrand(phase_number, states, controls,
time, parameters, user)
Description of Input
| phase_number | Number of the phase, 1:N, for which to
evaluate the integrand function. Length: 1 |
| states | The state variable values at time time. Length: ns |
| controls | The control variable values at time time. Length: nc |
| time | The time value to evaluate at. Length: 1 |
| parameters | The parameter variable values. Length: np |
| user | The user variable from problem.user. |
Description of Output
| L | The computed function value. Length: 1 |
| iferr | See costfun_endpoint. |
cost_pattern(3) applies to this callback. A nonzero value
indicates there is an integrand defined in this callback, and the
callback will be called. A zero value indiciate there is no
integrand defined (equal to zero). No call to this function for
the phase in question will be done in that case.
The name of the odefun callback m-file is given in problem.odefun (3.3.3).
Purpose
The purpose of the user supplied odefun callback is to
compute the right-hand-side of the differential equations for
all phases. There is no pattern for odefun ; all differential
equations must be defined.
Calling Syntax
[f, iferr] = odefun(phase_number, states, controls, time,
parameters, user)
Description of Input
| phase_number | Number of the phase, 1:N, for which to
evaluate the right-hand-side functions. Length: 1 |
| states | The state variable values at time time. Length: ns |
| controls | The control variable values at time time. Length: nc |
| time | The time value to evaluate at. Length: 1 |
| parameters | The parameter variable values. Length: np |
| user | The user variable from problem.user. |
Description of Output
| f | The computed function values. Length: ns |
| iferr | See costfun_endpoint. |
The name of the eventfun callback m-file is given in problem.eventfun (3.3.4).
Purpose
The purpose of the user supplied eventfun callback is to
compute the endpoint functions of the event constraints
(ψI and ψF) for all phases.
Description
All selected nonlinear event constraint endpoint functions for a
given phase and endpoint must be evaluated in one call and returned.
The selected endpoint functions are determined by the pattern
event_pattern. The order of the function values to be
returned is determined by the order of the constraints: a function
value belonging to a constraint with lower number has a lower index
in the returned vector than a function value belonging to a
constraint with higher number.
Calling Syntax
[psi, iferr] = eventfun(phase_number, endpoint, states, controls,
time, parameters, user)
Description of Input
| phase_number | Number of the phase, 1:N, for which to
evaluate the event constraints. Length: 1 |
| endpoint | Flag telling whether to compute the initial
endpoint function or the final endpoint
function. If endpoint =-1, compute the
initial endpoint function, if endpoint =1,
compute the final endpoint function. 1 and
-1 are the only possible values of
endpoint. Length: 1 |
| states | The state variable values at time time. Length: ns |
| controls | The control variable values at time time. Length: nc |
| time | The time value. Is equal to tI or tF
depending on the value of endpoint. Length: 1 |
| parameters | The parameter variable values. Length: np |
| user | The user variable from problem.user. |
Description of Output
| psi | The computed function values. Each value
corresponds to one selected nonlinear
endpoint function. event_pattern is used to
determine what nonlinear endpoint functions
are selected. The length of psi is therefore
determined by event_pattern as well, as the
number of nonlinear endpoint function for a
given phase and endpoint, or equivalently,
as the number of nonzeros in respective
column of event_pattern. If no event_pattern is given, it is assumed that the pattern is full, and that all nonlinear endpoint functions are selected. Length: |
| iferr | See costfun_endpoint. |
Example
An example of how event_pattern affects the length of psi.
We have five event constraints:
|
Notice that ψI2, ψF2, ψF3 and ψI5 are undefined, and equal to zero, i.e. not selected. γ(I/F)k are the linear terms, and γ is set through the matrices problem.phases.C_I and problem.phases.C_F.
The event_pattern for this problem would look like:
event_pattern = [1 1; % psi_I1 and psi_F1 selected in
% constraint 1
0 0; % No nonlinear endpoint functions in
% constraint 2
1 0; % Only psi_I3 selected here.
1 1; % and so on...
0 1];
The output psi from eventfun should be:
if endpoint == -1
psi(1) = psi_I_1; % for constraint 1
psi(2) = psi_I_3; % for constraint 3
psi(3) = psi_I_4; % for constraint 4
else % if endpoint == 1
psi(1) = psi_F_1; % for constraint 1
psi(2) = psi_F_4; % for constraint 4
psi(3) = psi_F_5; % for constraint 5
end
In this example, it happened to be equal number of funcions (3) defined
for both endpoints. This is not neccesary in a general case.
The name of the pathfun callback m-file is given in problem.pathfun (3.3.5).
Purpose
The purpose of the user supplied pathfun callback is to
compute the nonlinear functions g of the path constraints
for all phases.
Description
All selected g functions of the path constraints for a given phase
must be evaluated in one call and returned. The selected g
functions of the path constraints are determined by the pattern
path_pattern. The order of the function values to be returned
is determined by the order of the constraints: a function value
belonging to a constraint with lower number has a lower index in the
returned vector than a function value belonging to a constraint with
higher number.
Calling Syntax
[g, iferr] = pathfun(phase_number, states, controls, time,
parameters, user)
Description of Input
| phase_number | Number of the phase, 1:N, for which to
evaluate the g functions. Length: 1 |
| states | The state variable values at time time. Length: ns |
| controls | The control variable values at time time. Length: nc |
| time | The time value to evaluate at. Length: 1 |
| parameters | The parameter variable values. Length: np |
| user | The user variable from problem.user. |
Description of Output
| g | The computed function values. Each value
corresponds to one selected g
function. path_pattern is used to determine
what g functions are selected. The
length of output g is therefore determined
by path_pattern as well, as the number of
selected g functions for a given phase,
or equivalently, as the number of nonzeros
in path_pattern. If no path_pattern is given, it is assumed that the pattern is full, and that all g functions are selected. Length: |
| iferr | See costfun_endpoint. |
Example
An example of how path_pattern affects the length of g.
We have five path constraints:
|
Notice that g1 and g4 are not selected, and equal to zero. αk are the linear terms set through the matrix problem.phases.A, and βi are the coefficients for some auxiliary functions ai. The β coefficients are defined in problem.phases.nBterms and problem.phases.B, while the functions a are defined in pathfun_auxiliary.
The path_pattern for this prolem would look like:
path_pattern = [0; % No nonlinear g function in constraint 1
1; % g_2 is selected
1; % g_3 as well
0; % No g function in constraint 4
1]; % But one in constraint 5.
The output from pathfun should be:
g(1) = g_2;
g(2) = g_3;
g(3) = g_5;
The name of the pathfun_auxiliary callback m-file is given in problem.pathfun_auxiliary (3.3.6).
Purpose
The purpose of the user supplied pathfun_auxiliary callback
is to compute the auxliliary functions a of the path constraints
for all phases.
Description
As one phase can consist of several path constraints, and each
path constraint can have several auxiliary functions the evaluated
values of them need to be ordered correctly when returned from
this callback. The major order is the constraints the auxiliary
functions belong to, while the inner order then is the order of
which they are tied to the coefficients of beta. The number of
returned function values is the total number of auxiliary
functions in the given phase. See below for an example.
There is no pattern for the auxiliary functions.
Calling Syntax
[a, iferr] = pathfun_auxiliary(phase_number, states, controls,
time, parameters, user)
Description of Input
| phase_number | Number of the phase, 1:N, for which to
evaluate auxiliary functions. Length: 1 |
| states | The state variable values at time time. Length: ns |
| controls | The control variable values at time time. Length: nc |
| time | The time value to evaluate at. Length: 1 |
| parameters | The parameter variable values. Length: np |
| user | The user variable from problem.user. |
Description of Output
| a | The computed function values. The ordering
and length is described above, and in an
example below. There is no pattern for the
auxiliary functions. Length: |
| iferr | See costfun_endpoint. |
Example
An example of how to order the auxiliary functions in a.
We have five path constraints:
|
We have auxiliary function terms in constraint one and four. Let β1 be a vector of length 5, and β4 a vector of length 3:
beta_1 = [1 2 3 4 5];
beta_4 = [9 8 7];
Consequently, a1 is a vector-valued function of length 5, and a4 a vector-valued function of length 3.
Constraint one comes before constraint four (1 < 4), hence the a1 function values should be returned before the a4 function values. The internal order within a1 is an order corresponding to the coefficients in β1, and the internal order within a4 corresponds to β4.
The user should return the following from pathfun_auxiliary:
a(1) = a_1(1);
a(2) = a_1(2);
a(3) = a_1(3);
a(4) = a_1(4);
a(5) = a_1(5);
a(6) = a_4(1);
a(7) = a_4(2);
a(8) = a_4(3);