9. Highlevel Interface¶
The FORCESPRO highlevel interface gives optimization users a familiar easytouse way to define an optimization problem. The interface also gives the advanced user full flexibility when importing external Ccoded functions to evaluate the quantities involved in the optimization problem.
This interface is provided with Variant L and partially with Variant M of FORCESPRO.
Important
Starting with FORCESPRO 1.8.0
, the solver generated from the highlevel interface supports nonlinear and convex
decision making problems with integer variables.
Note
The highlevel Python interface expects 0based indices in the model formulation, such as for the indices in lbidx, ubidx, hlidx, huidx, xinitidx and xfinalidx, as is usual in Python programs. Note that this is contrary to the lowlevel interface, which requires 1based indices for these fields.
9.1. Supported problems¶
9.1.1. Canonical problem for discretetime dynamics¶
The FORCESPRO NLP solver solves (potentially) nonconvex, finitetime nonlinear optimal control problems with horizon length \(N\) of the form:
for \(k = 1,\dots,N\), where \(z_k \in \mathbb{R}^{n_k}\) are the optimization variables, for example a collection of inputs, states or outputs in an MPC problem; \({\color{red} p_{\color{red} k}} \in \mathbb{R}^{l_k}\) are realtime data, which are not necessarily present in all problems; the functions \(f_k : \mathbb{R}^{n_k} \times \mathbb{R}^{l_k} \rightarrow \mathbb{R}\) are stage cost functions; the functions \(c_k : \mathbb{R}^{n_k} \times \mathbb{R}^{l_k} \rightarrow \mathbb{R}^{w_k}\) represents (potentially nonlinear) equality constraints, such as a state transition function; the matrices \(E_k\) are used to couple variables from the \((k+1)\)th stage to those of stage \(k\) through the function \(c_k\); and the functions \(h_k : \mathbb{R}^{n_k} \times \mathbb{R}^{l_k} \rightarrow \mathbb{R}^{m_k}\) are used to express potentially nonlinear, nonconvex inequality constraints. The index sets \(\mathcal{I}\) and \(\mathcal{N}\) are used to determine which variables are fixed to initial and final values, respectively. The initial and final values \(\color{red} z_{\mathrm{init}}\) and \(\color{red} z_{\mathrm{final}}\) can also be changed in realtime. At every stage \(k\), the matrix \(F_k\) is a selection matrix that picks some coordinates in vector \(z_k\).
All realtime data is coloured in red. Additionally, when integer variables are modelled, they need to be declared as realtime parameters. See Section Mixedinteger nonlinear solver.
To obtain a solver for this optimization problem using the FORCESPRO client, you need to define all functions involved \((f_k, c_k, h_k)\) along with the corresponding dimensions.
9.1.2. Continuoustime dynamics¶
Instead of having discretetime dynamics as can be seen in Section 9.1.1, we also support using continuoustime dynamics of the form:
and then discretizing this equation by one of the standard integration methods. See Section 9.2.4.2 for more details.
9.1.3. Other variants¶
Not all elements in Section 9.1.1 have to be necessarily present. Possible variants include problems:
where all functions are fixed at code generation time and do not need extra realtime data \(\color{red} p\);
with no lower (upper) bounds for variable \(z_{k,i}\), then \(\underline{z}_i \equiv \infty (\bar{z}_i \equiv + \infty)\);
without nonlinear inequalities \(h\);
with \(N = 1\) (single stage problem), then the interstage equality can be omitted;
with final equality, then the final stage cost function f_N can be omitted;
that optimize over the initial value \(\color{red} z_{\mathrm{init}}\) and do not include the initial equality;
that optimize over the final value \(\color{red} z_{\mathrm{final}}\) and do not include the final equality;
mixedinteger nonlinear programs, where some variables are declared as integers. See Section Mixedinteger nonlinear solver for more information about the MINLP solver.
9.1.4. Function evaluations¶
The FORCESPRO NLP solver requires external functions to evaluate:
the cost function terms \(f_k(z_k)\) and their gradients \(\nabla f_k(z_k)\),
the dynamics \(c_k(z_k)\) and their Jacobians \(\nabla c_k (z_k)\), and
the inequality constraints \(h_k(z_k)\) and their Jacobians \(\nabla h_k(z_k)\).
The FORCESPRO code generator supports the following ways of supplying these functions:
1. Automatic Ccode generation of these functions from MATLAB using a supported automatic differentiation (AD) tool, such as CasADi. This happens automatically in the background, as long as the AD tool is found on your system. By doing so, the user does not need to adhere to any toolspecific syntax but can use standard MATLAB commands to define the necessary functions instead (which are then automatically converted to match the speficics of the chosen AD tool). This is the recommended way of getting started with FORCESPRO NLP. See Section 9.2 to learn how to use this approach.
2. Cfunctions (source files). These can be handcoded, or generated by any automatic differentiation tool. See Section 9.5 for details on how to provide own function evaluations and derivatives to FORCESPRO.
9.2. Expressing the optimization problem in code¶
When solving nonlinear programs of the type in Section 9.1.1, FORCESPRO requires the functions \(f,c,h\) and their derivatives (Jacobians) to be evaluated in each iteration. There are two ways for accomplishing this: either implement all function evaluations in C by some other method (by hand or by another automatic differentiation tool), or use an AD tool integrated with FORCESPRO, such as the opensource package CasADi (see Automatic differentiation tool for a list of all supported tools). This is the easiest option to quickly get started with solving NLPs, and it generates efficient code. However, if you prefer other AD tools, see Section 9.5 to learn how to provide your own derivatives to FORCESPRO NLP solvers. This section will describe the CasADibased approach in detail, using either the MATLAB or the Python client of FORCESPRO. Please note that even though both the MATLAB and the Python client are intended to behave largely identical, there are some differences between the two clients. For details, refer to Differences between the MATLAB and the Python client.
9.2.1. Model Initialization¶
9.2.1.1. Model Initialization in Matlab¶
In the MATLAB highlevel interface, the formulation of the optimization problem is given through a simple structure array. In the following, we will describe the problem in such an array named model. It is advisable to zeroinitialize this variable at the beginning of your script such that no values set in previous iterations of your script interfere with this run:
model = {}
9.2.1.2. Model Initialization in Python¶
In the highlevel Python interface, optimization problems are described through objects of different types, depending on the problem. The following classes are available:
SymbolicModel  Allows you to describe your optimization problem using regular Python functions. These functions will be evaluated symbolically by CasADi to create optimized C code. Note that this model is meant to be used for nonlinear models. If you wish to express a convex model symbolically, consider using the ConvexSymbolicModel or forcing generation of a nonconvex solver by seting the option forcenonconvex to True.
ExternalFunctionModel  Enables more flexibility in describing nonlinear problems by allowing any external function to be used as objective function and constraints. This requires C code or already compiled code (object files or shared libraries) from any language. The approach using external function evaluations for your objective function and constraints is described in External function evaluations in C, including the required call signature of the callback function.
ConvexSymbolicModel  FORCESPRO can generate different solvers for convex problems.
Whichever model you choose, it can be initialized with no arguments, or with a single argument denoting the number of stages \(N\) in the problem:
import forcespro.nlp
model = forcespro.nlp.SymbolicModel(50)
Note that most symbolic problem descriptions will also require the Numpy and CasADi packages, so it is a good idea to import them at the beginning:
import numpy as np
import casadi
9.2.2. Dimensions¶
In order to define the dimensions of the stage variables \(z_i\), the number of equality and inequality constraints and the number of realtime parameters use the following fields (properties) in the client:
model.N = 50; % length of multistage problem
model.nvar = 6; % number of stage variables
model.neq = 4; % number of equality constraints
model.nh = 2; % number of nonlinear inequality constraints
model.npar = 0; % number of runtime parameters
model.N = 50 # not required if already specified in initializer
model.nvar = 6 # number of stage variables
model.neq = 4 # number of equality constraints
model.nh = 2 # number of nonlinear inequality constraints
model.npar = 0 # number of runtime parameters
If the dimensions vary for different stages use arrays of length \(N\). See Section 9.2.7.1 for an example.
9.2.3. Objective¶
The highlevel interface allows you to define the objective function using a handle to a MATLAB or Python function that evaluates the objective. This function is called with the variables of one stage as its first argument, i.e. a vector of model.nvar entries. FORCESPRO will process the given function symbolically and generate the necessary C code to be included in the solver.
model.objective = @eval_obj; % handle to objective function
model.objective = eval_obj # eval_obj is a Python function
For instance, the function could be:
function f = eval_obj ( z )
F = z(1);
s = z(2);
y = z(4);
f = 100*y + 0.1*F^2 + 0.01* s^2;
end
def eval_obj(z):
F = z[0]
s = z[1]
y = z[3]
return 100*y + 0.1*F**2 + 0.01*s**2
If the cost function varies for different stages use a cell array of function handles of length \(N\) in MATLAB, or a list of function handles in Python. See Section 9.2.7.1 for an example.
Note
Python and MATLAB use different indexing bases. The first element of any variable has index 1 in MATLAB, whereas it is accessed at offset 0 in Python!
The objective evaluation function can optionally accept an additional argument p which serves as a runtime parameter. In order to be able to change the terms in the cost function during runtime, one can define the objective function as:
function f = eval_obj ( z, p )
F = z(1);
s = z(2);
y = z(4);
f = 100*y + p(1)*F^2 + p(2)* s^2;
end
def eval_obj(z, p):
F = z[0]
s = z[1]
y = z[3]
return 100*y + p[0]*F**2 + p[1]*s**2
The length of this additional parameter vector in each stage is given by model.npar.
9.2.4. Equalities¶
9.2.4.1. Discretetime¶
For discretetime dynamics, one can define a handle to a function evaluating \(c\) as shown below. The selection matrix \(E\) that determines which variables are affected by the interstage equality must also be filled. For performance reasons, it is recommended to order variables such that the selection matrix has the following structure:
model.eq = @eval_dynamics; % handle to interstage function
model.E = [zeros(4,2), eye(4)]; % selection matrix
model.eq = eval_dynamics # handle to interstage function
model.E = np.concatenate([np.zeros((4, 2)), np.eye(4)], axis=1) # selection matrix
If the equality constraint function varies for different stages use a cell array (or list in Python) of function handles of length \(N1\), and similarly for \(E_k\). See Section 9.2.7.1 for an example.
9.2.4.2. Continuoustime¶
For continuoustime dynamics, FORCESPRO requires you to describe the dynamics of the system in the following form:
where \(x\) are the states of the system, \(u\) are the inputs and \(\color{red} p\) a vector of parameters, e.g. the mass or inertia. The selection matrix \(E\) determines which components of the stage variable \(z_i\) are to be considered state \(x\) or input \(u\) in this representation.
For example, let’s assume that the system to be controlled has the dynamics:
In order to discretize the system for use with FORCESPRO we have to:
Implement the continuoustime dynamics as a function:
function xdot = continous_dynamics(x, u, p)
xdot = p(1)*x(1)*x(2) + p(2)*u;
end
def continuous_dynamics(x, u, p):
return p[0]*x[0]*x[1] + p[1]*u[0]
Note that in general the parameter vector p can be omitted if there are no parameters. You can also implement short functions as anonymous function handles:
continous_dynamics_anonymous = @(x,u,p) p(1)*x(1)*x(2) + p(2)*u;
continuous_dynamics_anonymous = lambda x, u, p: p[0]*x[0]*x[1] + p[1]*u[0]
2. Tell FORCESPRO that you are using continuoustime dynamics by setting the
continuous_dynamics
field of the model
to a function handle to one of
the functions above:
model.continuous_dynamics = @continuous_dynamics;
model.continuous_dynamics = continuous_dynamics
or, if you are using anonymous functions:
model.continuous_dynamics = @continuous_dynamics_anonymous;
model.continuous_dynamics = continuous_dynamics_anonymous
3. Use the selection matrix \(E\) to link the stage variables \(z_i\) with the states \(x\) and inputs \(u\) of the continuous dynamics function:
model.E = [zeros(2, 1), eye(2)]
model.E = np.concatenate([np.zeros((2, 1)), np.eye(2)], axis=1)
Components of \(z_i\) are considered as state variables \(x\) according to the order prescribed by the selection matrix. If an entire \(k\)th column of the selection matrix is zero, the \(k\)th component of \(z_i\) is not governed by a dynamic equation and thus considered as input \(u\).
Choose one of the integrator functions from the Integrators section (the default is ERK4):
codeoptions.nlp.integrator.type = 'ERK2';
codeoptions.nlp.integrator.Ts = 0.1;
codeoptions.nlp.integrator.nodes = 5;
codeoptions.nlp.integrator.type = 'ERK2'
codeoptions.nlp.integrator.Ts = 0.1
codeoptions.nlp.integrator.nodes = 5
where the integrator type is set using the type field of the options struct
codeoptions.nlp.integrator
. The field Ts
determines the absolute time
between two integration intervals, while nodes
defines the number of
intermediate integration nodes within that integration interval. In the example
above, we use 5 steps to integrate for 0.1 seconds, i.e. each integration step
covers an interval of 0.02 seconds.
9.2.5. Initial and final conditions¶
The indices affected by the initial and final conditions can be set as follows:
model.xinitidx = 3:6; % indices affected by initial condition
model.xfinalidx = 5:6; % indices affected by final condition
model.xinitidx = range(2, 6) # indices affected by the initial condition
model.xfinalidx = range(4, 6) # indices affected by the final condition
Note
Python and MATLAB use different indexing bases. The first variable in a stage has index 1 in MATLAB, whereas it is accessed at offset 0 in Python! Further note that Python’s range does not include the upper limit, thus:
list(range(2, 6)) == [2, 3, 4, 5] # does not include upper limit
9.2.6. Inequalities¶
A function evaluating nonlinear inequalities can be provided in a similar way, for example:
function h = eval_const(z)
x = z(3);
y = z(4);
h = [x^2 + y^2;
(x+2)^2 + (y2.5)^2 ];
end
def eval_const(z):
x = z[2]
y = z[3]
return np.array([x**2 + y**2;
(x+2)**2 + (y2.5)**2])
Note
For Python installations with Numpy version 1.20 onwards it is advised to use CasADi arrays and CasADi functions (where available) for the implementation of the functions assigned to: model.objective, model.eq, model.ineq, model.continous_dynamics for the problem formulation in order to ensure maximum compatibility between CasADi and the FORCESPRO Python client.
The simple bounds and the nonlinear inequality bounds can have
+inf
and inf
elements, but must be the same length as the field
nvar
and nh
, respectively:
model.ineq = @eval_const; % handle to nonlinear constraints
model.hu = [9, +inf]; % upper bound for nonlinear constraints
model.hl = [1, 0.95^2]; % lower bound for nonlinear constraints
model.ub = [+5, +1, 0, 3, 2, +pi]; % simple upper bounds
model.lb = [5, 1, 3, inf, 0, 0]; % simple lower bounds
model.ineq = eval_const # handle to nonlinear constraints
model.hu = [9, +float('inf')] # upper bound for nonlinear constraints
model.hl = [1, 0.95**2] # lower bound for nonlinear constraints
model.ub = [+5, +1, 0, 3, 2, +np.pi] # simple upper bounds
model.lb = [5, 1, 3, float('inf'), 0, 0] # simple lower bounds
Note
While the FORCESPRO Python client does not require you to use numpy arrays, we encourage their use for vector and matrixvalued properties of the model, as it simplifies calculations for the user. Therefore, any of the above properties can also be set to Numpy arrays instead of lists. If lists are given, these are converted to Numpy arrays internally.
If the constraints vary for different stages, use cell arrays of length \(N\) for any of the quantities defined above. See Varying dimensions, parameters, constraints, or functions section for an example.
Bounds model.lb
, model.ub
, model.hl
and model.hu
can be made
parametric by leaving said fields empty and using the model.lbidx
,
model.ubidx
, model.hlidx
and model.huidx
fields respectively to
indicate on which variables/inequalities lower and upper bounds are present. The numerical
values will then be expected at runtime. The runtime parameters will be
created stagewise for the above bounds and will have the names
lb<n>
, ub<n>
, hl<n>
, hu<n>
where <n>
will be the 1based stage number they
belong it (padded with enough 0 based on the largest stage)
For example, to set parametric lower bounds on states 1 and 2,
and parametric upper bounds on states 2 and 3, use:
% Lower bounds are parametric (indices not mentioned here are inf)
model.lbidx = [1 2]';
% Upper bounds are parametric (indices not mentioned here are +inf)
model.ubidx = [2 3]';
% lb and ub have to be empty when using parametric bounds
model.lb = [];
model.ub = [];
# Lower bounds are parametric (indices not mentioned here are inf)
model.lbidx = [0, 1]
# Upper bounds are parametric (indices not mentioned here are +inf)
model.ubidx = [1, 2]
# There is no need to specify model.lb or model.ub to empty lists if
# model.lbidx or model.ubidx are set, and any nonempty value is disallowed.
and then specify the exact values at runtime through the fields
lb01
–lbN
and ub01
–ubN
:
% Specify lower bounds
problem.lb01 = [0 0]';
problem.lb02 = [0 0]';
% ...
% Specify upper bounds
problem.ub01 = [3 2]';
problem.ub02 = [3 2]';
% ...
# Specify lower bounds
problem["lb01"] = [0, 0]
problem["lb02"] = [0, 0]
# Specify upper bounds
problem["ub01"] = [3, 2]
problem["ub02"] = [3, 2]
Tip
One could use problem.(sprintf('lb%02u',i))
in an i
indexed loop to
set the parametric bounds more easily in the MATLAB client. Similarly, the
parametric bounds for the stages can be set using
problem["{:02d}".format(i+1)]
in Python. Alternatively, consider using
the option stack_parambounds, described below.
If the model.lbidx
and model.ubidx
fields vary for different stages use
cell arrays of length \(N\).
From Release 3.0.1
, the parametric bounds can be stacked on one same array covering all stages.
To enable this behaviour, users need to set the following codegeneration option:
codeoptions.nlp.stack_parambounds = 1;
codeoptions.nlp.stack_parambounds = True
This option is effective for both the PDIP_NLP
and SQP_NLP
solve methods and works with bounds on variables and inequalities. At runtime, users can then write
% Lower and upper bounds stacked over all stages
problem.lb = [0 0 0 0 ...];
problem.ub = [3 2 3 2 ...];
# Lower and upper bounds stacked over all stages
problem["lb"] = [0, 0, 0, 0, ...]
problem["ub"] = [3, 2, 3, 2, ...]
Alternatively, if you want to use the same bounds across all stages:
problem.lb = repmat([0, 0], 1, model.N);
problem.ub = repmat([3, 2], 1, model.N);
problem["lb"] = np.tile([0, 0], (model.N,))
problem["ub"] = np.tile([3, 2], (model.N,))
9.2.7. Variations¶
9.2.7.1. Varying dimensions, parameters, constraints, or functions¶
The example described above has the same dimensions, bounds and functions for the whole horizon. One can define varying dimensions using arrays and varying bounds and functions using MATLAB cell arrays or Python lists. For instance, to remove the first and second variables from the last stage one could write the following:
for i = 1:model.N1
model.nvar(i) = 6;
model.objective{i} = @(z) 100*z(4) + 0.1*z(1)^2 + 0.01*z(2)^2;
model.lb{i} = [5, 1, 3, 0, 0, 0];
model.ub{i} = [+5 , +1, 0, 3, 2, +pi];
if i < model.N1
model.E{i} = [zeros(4, 2), eye(4)];
else
model.E{i} = eye(4);
end
end
model.nvar(model.N) = 4;
model.objective{model.N} = @(z) 100*z(2);
model.lb{model.N} = [3, 0, 0, 0];
model.ub{model.N} = [ 0, 3, 2, +pi];
model = forcespro.nlp.SymbolicModel(50) # to set values stagewise, the model must be initialized this way
for i in range(0,model.N1):
model.nvar[i] = 6
model.objective[i] = lambda z: 100*z[3] + 0.1*z[0]**2 + 0.01*z[1]**2
model.lb[i] = [5, 1, 3, 0, 0, 0]
model.ub[i] = [+5 , +1, 0, 3, 2, +np.pi]
if i < model.N2:
model.E[i] = np.concatenate([np.zeros(4, 2), np.eye(4)], axis=1)
else:
model.E[i] = np.eye(4)
model.nvar[1] = 4
model.objective[1] = lambda z: 100*z[1]
model.lb[1] = [3, 0, 0, 0]
model.ub[1] = [ 0, 3, 2, +np.pi]
It is also typical for model predictive control problems (MPC) that only the last stage differs from the others (excluding the initial condition, which is handled separately). Instead of defining cell arrays as above for all stages, FORCESPRO offers the following shorthand notations that alter the last stage:
nvarN
: number of variables in last stagenparN
: number of parameters in last stageobjectiveN
: objective function for last stageEN
: selection matrix \(E\) for last stage updatenhN
: number of inequalities in last stageineqN
: inequalities for last stage
Add any of these fields to the model
struct/object to override the default
values, which is to make everything the same along the horizon. For example, to
add a terminal cost that is a factor 10 higher than the stage cost:
model.objectiveN = @(z) 10*model.objective(z);
model.objectiveN = lambda z: 10*model.objective(z)
9.2.7.2. Providing analytic derivatives¶
The algorithms inside FORCESPRO need the derivatives of the functions
describing the objective, equality and inequality constraints. The code
generation engine uses algorithmic differentiation (AD) to compute these
quantities. Instead, when analytic derivatives are available, the user can
provide them using the fields model.dobjective
, model.deq
, and
model.dineq
.
Note that the user must be particularly careful to make sure that the provided functions and derivatives are consistent, for example:
model.objective = @(z) z(3)^2;
model.dobjective = @(z) 2*z(3);
model.objective = lambda z: z[2]**2
model.dobjective = lambda z: 2*z[2]
The code generation system will not check the correctness of the provided derivatives.
9.2.8. Single precision callbacks¶
Evaluating objective function, dynamics and constraints as well as their respective derivatives may take a significant part of the overall solution time (both total and per iteration). In such situations solution time and memory consumption may be sped up by evaluating those functions in single, rather than double precision arithmetic. This can be done by specifying
codeoptions.callback_floattype = 'float';
# not yet supported
Note that this will allow to run the NLP solver in mixedprecision arithmetic,
where the callbacks are evaluated in single precision, but the overall algorithm in
double precision. In order for this to work well, all callbacks functions need to
be numerically wellposed and overall accuracy requirements of the solution must not
be too high. In particular, when using that feature, you may need to reduce some of the
accuracy settings (such as codeoptions.nlp.TolStat
) by one or two orders of magnitude,
see Accuracy requirements.
Note
Single precision callbacks are currently supported for legacy and chainrule integrators, but not yet for VDE integrators. Also, this features is currently only available via the Matlab client.
9.3. Generating a solver¶
In addition to the definition of the NLP, solver generation requires an (optional) set of options for customization (see the Solver Options section for more information). Using the default solver options we generate a solver using:
% Get the default solver options
codeoptions = getOptions('FORCESNLPsolver');
% Generate solver
FORCES_NLP(model, codeoptions);
# Get the default solver options
options = forcespro.CodeOptions('FORCESNLPsolver')
# Generate solver for previously initialized model
solver = model.generate_solver(options)
As the solver is generated, several files are downloaded into the current working directory of the calling script, including the compiled solver itself and MATLAB/Python interfaces for calling it.
Note
In the Python client, generate_solver() returns a solver object. This object can be used to call the solver. To get a solver object for a previously generated solver in some directory /path/to/solver, use:
import forcespro.nlp
solver = forcespro.nlp.Solver.from_directory('/path/to/solver')
9.3.1. Declaring outputs¶
By default, the solver returns the solution vector for all stages as multiple
outputs. Alternatively, the user can pass a third argument to the function
FORCES_NLP
with an array that specifies what the solver should output. For
instance, to define an output, named u0
, to be the first two elements of the
solution vector at stage 1, use the following commands:
output1 = newOutput('u0', 1, 1:2);
FORCES_NLP(model, codeoptions, output1);
output_1 = ("u0", 0, [0, 1], "")
model.generate_solver(options, [output_1])
Additionally, you can request that the solver returns the solution vectors for all different stages “stacked” together into a single vector, say called sol
, by using the following commands:
output1 = newOutput('sol');
FORCES_NLP(model, codeoptions, output1);
output1 = ("sol", [], [])
model.generate_solver(options, [output1])
Important
When using the MINLP solver and defining outputs, all integer variables need to be specified as custom outputs.
The dual variables at the solution returned by FORCESPRO provide useful information on the problem sensitivity. They can be exported from the nonlinear solver as well by giving the maps2const
field one of the following values:
‘nl_eq_dual’ for the dual variables associated to equality constraints
‘nl_lb_var_dual’ for the dual variables associated to lower bounds on variables
‘nl_ub_var_dual’ for the dual variables associated to upper bounds on variables
‘nl_ip_ineq_dual’ for the dual variables associated to nonlinear inequalities
‘nl_ineq_slack’ for the dual variables associated to slacks on nonlinear inequalities.
An example of exporting the marginals associated to nonlinear equalities is shown in the code snippet below.
outputs(4) = newOutput('dual_eq0', 1:model.N, 1:2, 'nl_eq_dual');
9.4. Calling the solver¶
After code generation has been successful, one can obtain information about the realtime data needed to call the generated solver by typing:
help FORCESNLPsolver
# Assuming `solver` is the return value of a `model.generate_solver()` call
solver.help()
In Python, a previously generated solver can be loaded as follows:
import forcespro.nlp
solver = forcespro.nlp.Solver.from_directory("/path/to/generated/solver/")
solver.help()
9.4.1. Initial guess¶
The FORCESPRO NLP solver solves NLPs to local optimality, hence the resulting optimal solution depends on the initialization of the solver. One can also choose another initialization point when a better guess is available. The following code sets the initial point to be in the middle of all bounds:
x0i = model.lb +(model.ub  model.lb)/2;
x0 = repmat(x0i', model.N, 1);
problem.x0 = x0;
xi = (model.lb + model.ub) / 2 # assuming lb and ub are numpy arrays
x0 = np.tile(xi, (model.N,))
problem = {"x0": x0}
9.4.2. Initial and final conditions¶
If there are initial and/or final conditions on the optimization variables, the solver will expect the corresponding runtime fields:
problem.xinit = model.xinit;
problem.xfinal = model.xfinal;
problem = {"xinit": np.array([1, 2, 3]),
"xfinal": np.array([4, 5, 6])}
Note that the Python client does not allow setting model.xinit or model.xfinal properties, as those are runtime parameters not needed at solver generation time.
9.4.3. Realtime parameters¶
Whenever there are any runtime parameters defined in the problem, i.e. the field
npar
is not zero, the solver will expect the following field containing the
parameters for all the \(N\) stages stacked in a single vector:
problem.all_parameters = repmat(1.0, model.N, 1);
problem["all_parameters"] = np.tile(1.0, (model.N,))
9.4.4. Tolerances as realtime parameters¶
From FORCESPRO 2.0 onwards, the NLP solver tolerances can be made realtime parameters, meaning that they do not need to be set when generating the solver but can be changed at runtime when calling the generated solver. The codesnippet below shows how to make the tolerances on the gradient of the Lagrangian, the equalities, the inequalities and the complementarity condition parametric. Essentially, when the tolerances are declared nonpositive at codegeneration, the corresponding runtime parameter is created in the solver.
codeoptions.nlp.TolStat = 1; % Tolerance on gradient of Lagrangian
codeoptions.nlp.TolEq = 1; % Tolerance on equality constraints
codeoptions.nlp.TolIneq = 1; % Tolerance on inequality constraints
codeoptions.nlp.TolComp = 1; % Tolerance on complementarity
codeoptions.nlp.TolStat = 1 # Tolerance on gradient of Lagrangian
codeoptions.nlp.TolEq = 1 # Tolerance on equality constraints
codeoptions.nlp.TolIneq = 1 # Tolerance on inequality constraints
codeoptions.nlp.TolComp = 1 # Tolerance on complementarity
Once the tolerance has been declared nonpositive and the solver has been generated, the corresponding parameter can be set at runtime as follows:
problem.ToleranceStationarity = 1e1;
problem.ToleranceEqualities = 1e1;
problem.ToleranceInequalities = 1e1;
problem.ToleranceComplementarity = 1e1;
problem["ToleranceStationarity"] = 1e1
problem["ToleranceEqualities"] = 1e1
problem["ToleranceInequalities"] = 1e1
problem["ToleranceComplementarity"] = 1e1
Tip
We do not recommend changing the tolerance on the complementarity condition since it is used internally to update the barrier parameter. Hence loosening it may hamper the solver convergence.
9.4.5. Exitflags and quality of the result¶
Once all parameters have been populated, the MEX interface of the solver can be used to invoke it:
[output, exitflag, info] = FORCESNLPsolver(problem);
output, exitflag, info = solver.solve(problem)
The possible exitflags are documented in Table 9.1. The exitflag should always be checked before continuing with program execution to avoid using spurious solutions later in the code. Check whether the solver has exited without an error before using the solution. For example, in MATLAB, we suggest to use an assert statement:
assert(exitflag == 1, 'Some issue with FORCESPRO solver');
assert exitflag == 1, "Some issue with FORCESPRO solver"
Exitflag 
Description 


Local optimal solution found (i.e. the point satisfies the KKT optimality conditions to the requested accuracy). 

Maximum number of iterations reached. You can examine the value of optimality
conditions returned inside the 

Wrong number of inequalities input to solver. 

Error occured during matrix factorization. 



The solver could not proceed. Most likely cause is that the problem is infeasible.Try formulating a problem with slack variables (soft constraints) to avoid this error. 

The internal QP solver could not proceed. This exitflag can only occur when
using the Sequential quadratic programming algorithm. The most likely cause
is that an infeasible QP or a numerical unstable QP was encountered.
Try increasing the hessian regularization parameter 



Invalid values in problem parameters. 

License error. This typically happens if you are trying to execute code that has been generated with a Simulation license of FORCESPRO on another machine. Regenerate the solver using your machine. 
Besides the exitflag, the solver also returns an information structure containing detailed KPIs on the solver performance. All its entries are listed and explained in
Table 9.2. Those values may also be useful in case the solver
exitflag has been 0
, which means the solver did not fail but reached the maximum number
of iterations. In that case, one should at least check whether the “solution” returned is
sufficiently feasible. This can be done by examining res_eq
and res_ineq
, respectively,
to check whether equality and inequality constraints are satisfied up to a sufficiently small
tolerance. The exact tolerances to check may be strongly application dependent.
Note
Applying a premature “solution” returned along with an exitflag of 0
to control your
system may have undesired effects to the behaviour of that system. Only do so if you fully
understand what you are doing.
Fieldname 
Description 


Number of solver iterations that led to this result 

Maximum norm of equality constraint residual 

Maximum norm of inequality constraint residual 

Maximum norm of stationarity condition 

Maximum norm of violations of the complementarity conditions 

Primal objective value (as provided by the user) 

Duality measure 

Time needed to run the solver (wall clock time) 

Time needed just for function evaluations of all user callbacks inside the solver (wall clock time) 

Solver ID of generated solver 
9.5. External function evaluations in C¶
This approach allows the user to integrate existing efficient C implementations to evaluate the required functions and their derivatives with respect to the stage variable. This gives the user full flexibility in defining the optimization problem. In this case, the functions do not necessarily have to be differentiable, although the convergence of the algorithm is not guaranteed if they are not. When following this route the user does not have to provide MATLAB code to evaluate the objective or constraint functions. However, the user is responsible for making sure that the provided derivatives and function evaluations are coherent. The FORCESPRO NLP code generator will not check this.
9.5.1. Interface¶
9.5.1.1. Expected function signature¶
There are two ways in which a user can supply custom functions written in C:
Either she can supply all functions (
model.objective
,model.eq
,model.ineq
etc.)Or she can supply one or a few of these together with its differential/Jacobian.
Depending on the case, the user will have to supply different information when generating a FORCESPRO solver.
When supplying all functions, the user will have supply a single C function having the following signature:
int myfunctions (
double *x, /* primal vars */
double *y, /* eq. constraint multiplers */
double *l, /* ineq . constraint multipliers */
double *p, /* runtime parameters */
double *f, /* objective function ( incremented in this function ) */
double *nabla_f , /* gradient of objective function */
double *c, /* dynamics */
double *nabla_c , /* Jacobian of the dynamics ( column major ) */
double *h, /* inequality constraints */
double *nabla_h , /* Jacobian of inequality constraints ( column major ) */
double *H, /* Hessian ( column major ) */
int stage, /* stage number (0 indexed ) */
int iteration, /* Solver iteration count */
int threadID /* Thread id */
)
If instead the user wants to supply just a single function she will have to supply a single C function having the following signature:
void function (
const double * const x, /* primal vars */
const double * const p, /* runtime parameters */
double * const zeroOrderFcn, /* Zero order function */
double * const firstOrderFcn , /* first order Fcn */
const int stage, /* stage number (0 indexed ) */
const int threadID /* thread number */
)
where zeroOrderFcn
and firstOrderFcn
denote the function which the user wants to supply, together with its differential/Jacobian respectively. E.g. if the user were to add the objective function and its differential externally, the function might look as follows:
void objective (
const double * const x, /* primal vars */
const double * const p, /* runtime parameters */
double * const obj, /* objective function */
double * const nabla_obj, /* jacobian of objective fcn */
const int stage, /* stage number (0 indexed ) */
const int threadID /* thread number */
)
Note
External Cfunctions should have the same name as the file it is contained in, minus the file extension. E.g. in the above example the source file containing the definition of the function objective
would have to have the name objective.c
.
If all functions are provided as external C functions through the FORCESPRO Python client, then one can provide a different name for the function and the file.
9.5.1.2. Custom data structures as parameters¶
If you have an advanced data structure that holds the userdefined runtime parameters, and you do not want to serialize it into an array of doubles to use the interface above, you can invoke the option:
codeoptions.customParams = 1;
options.customParams = 1
When doing this, it is important to note that runtime parameters can only be passed to externally provided functions. In particular, if some but not all function evaluations are provided externally, one will have to set model.npar = 0
.
When using custom parameters, if all functions are provided externally, the expected function signature is:
int myfunctions (
double *x, /* primal vars */
double *y, /* eq. constraint multiplers */
double *l, /* ineq . constraint multipliers */
void *p, /* runtime parameters (passed as void pointer) */
double *f, /* objective function ( incremented in this function ) */
double *nabla_f , /* gradient of objective function */
double *c, /* dynamics */
double *nabla_c , /* Jacobian of the dynamics ( column major ) */
double *h, /* inequality constraints */
double *nabla_h , /* Jacobian of inequality constraints ( column major ) */
double *H, /* Hessian ( column major ) */
int stage, /* stage number (0 indexed ) */
int iteration, /* Solver iteration count */
int threadID /* Thread id */
)
If instead only some of the functions are provided, the expected function signature of these is
void function (
const double * const x, /* primal vars */
void * const p, /* runtime parameters passed as void pointer */
double * const zeroOrderFcn, /* Zero order function */
double * const firstOrderFcn , /* first order Fcn */
const int stage, /* stage number (0 indexed ) */
const int threadID /* thread number */
)
At run time you can then pass arbitrary data structures to your own function by setting the pointer in the params struct:
myData p; /* define your own parameter structure */
/* ... */ /* fill it with data */
/* Set parameter pointer to your data structure */
mySolver_params params; /* Define solver parameters */
params.customParams = &p;
/* Call solver (assuming everything else is defined) */
mySolver_solv(¶ms, &output, &info, stdout, &external_func);
Note
Setting customParams
to 1
will disable building highlevel interfaces. Only C header and source files will be generated. As a consequence,
if CasADi is used to generate some of the function evaluations, the generated source code will have to be compiled by the user.
Note
When using custom parameters, generating callback code automatically is only supported for CasADi 3.5.x only.
9.5.2. Supplying function evaluation information¶
In MATLAB, if the user wants to supply all functions externally she can let the code generator know about the path to the C files implementing the necessary function as follows:
model.extfuncs = 'C/myfunctions.c';
Alternatively she could supply either one of the functions model.objective
, model.eq
or model.ineq
together with its differential by specifying the path to the corresponding C source file in the corresponding field of model.extfuncs
as follows:
model.extfuncs.objective = 'C/myobjective.c'; % adding model.objective externally %
model.extfuncs.dynamics = 'C/mydynamics.c'; % adding model.eq externally %
model.extfuncs.inequalities = 'C/myinequalities.c'; % adding model.ineq externally %
As noted above, FORCESPRO derives the function name used for the callback from the file name; the function must therefore have the same name as the file in which it is contained.
In Python, if the user wishes to add ALL functions externally she would use a ExternalFunctionModel as follows:
model = forcespro.nlp.ExternalFunctionModel(50)
model.add_auxiliary(["helper_functions.c", "compiled_helper_functions.obj"])
model.set_main_callback("myfunctions.c", function="myfunctions")
Herein, the add_auxiliary() method is used to add any helper C source files or object files that should be compiled and linked against, and the set_main_callback() function is used to define the path to a C source file or compiled object file, as well as the name of an exported function that conforms to the call signature given above. This function will be used to evaluate any nonlinear constraints and the objective function.
Alternatively, in order to add a single function externally the user would use the SymbolicModel and add the C source files containing code for the external functions through model.extfuncs as follows:
model = forcespro.nlp.SymbolicModel(50)
model.extfuncs.objective = "myobjective.c"
model.extfuncs.dynamics = "mydynamics.c"
model.extfuncs.inequalities = "myinequalities.c"
model.add_c_source("extra_source.c")
where the extra_source.c are additional C source files needed for evaluating some or all of the external functions.
9.5.3. Rules for function evaluation code¶
The contents of the function have to follow certain rules. We will use the following example to illustrate them:
/* cost */
if (f)
{ /* notice the increment of f */
(*f) += 100*x[3] + 0.1* x[0]*x[0] + 0.01*x [1]*x [1];
}
/* gradient  only nonzero elements have to be filled in */
if ( nabla_f )
{
nabla_f [0] = 0.2*x[0];
nabla_f [1] = 0.02*x[1];
nabla_f [3] = 100;
}
/* eq constr */
if (c)
{
vehicle_dyanmics (x, c);
}
/* jacobian equalities ( column major ) */
if ( nabla_c )
{
vehicle_dyanmics_jacobian (x, nabla_c );
}
/* ineq constr */
if (h)
{
h[0] = x [2]*x[2] + x[3]*x [3];
h[1] = (x[2]+2)*(x[2]+2) + (x[3] 2.5)*(x[3] 2.5);
}
/* jacobian inequalities ( column major )
 only non  zero elements to be filled in */
if ( nabla_h )
{
/* column 3 */
nabla_h [4] = 2*x [2];
nabla_h [5] = 2*x[2] + 4;
/* column 4 */
nabla_h [6] = 2*x [3];
nabla_h [7] = 2*x[3]  5;
}
Notice that every function evaluation is only carried out if the corresponding pointer is not null. This is used by the FORCESPRO NLP solver to call the same interface with different pointers depending on the functions that it requires.
9.5.4. Matrix format¶
Matrices are assumed to be stored in dense column major format. However, only the nonzero components need to be populated, as FORCESPRO NLP makes sure that the arrays are initialized to zero before calling this interface.
9.5.5. Multiple source files¶
The use of multiple C files is also supported. In the example above, the
functions dynamics
and dynamics_jacobian
are defined in another file
and included as external functions using:
extern void dynamics ( double *x, double *c);
extern void dynamics_jacobian ( double *x, double *J);
In MATLAB, to let the code generator know about the location of these other files use a string with spaces separating the different files. In Python, use the add_auxiliary() method:
codeoptions.nlp.other_srcs = 'C/dynamics.c';
model.add_auxiliary('C/dynamics.c')
9.5.6. Stagedependent functions¶
Whenever the cost function in one of the stages is different from the standard cost function \(f\), one can make use of the argument stage to evaluate different functions depending on the stage number. The same applies to all other quantities.
9.6. Calling the nonlinear functions from Matlab or Python¶
From FORCESPRO version 5.0.1 onwards it is possible to evaluate the nonlinear functions and their derivatives directly in Matlab or Python. This can be very useful for debugging a FORCESPRO solver. This in particular can be useful when calling the generated solver returns an exitflag \(6\), \(8\) or \(10\). For help on how to do this, see Realtime SQP Solver: Robotic Arm Manipulator (MATLAB & Python) and Issues when running the solver.
9.6.1. Calling the nonlinear functions from MATLAB¶
When generating a FORCESPRO solver one can request the FORCESPRO client to also generate a MEX entry point to any of the nonlinear functions (equality constraints, inequality constraints and objective function) so that one can evaluate these along with their derivatives directly in MATLAB. This is done via the following three options:
dynamics
: In order to generate a MEX interface for the equality constraints, setcodeoptions.MEXinterface.dynamics = 1
.inequalities
: In order to generate a MEX interface for the inequality constraints, setcodeoptions.MEXinterface.inequalities = 1
.objective
: In order to generate a MEX interface for the objective function, setcodeoptions.MEXinterface.objective = 1
.
The generated MEX entry point carries the name <solver name>_dynamics
, <solver name>_inequalities
or <solver name>_objective
depending on the case.
The first argument to each of these generated MEX entry points is the primal variable \(z_k\) and the second argument is the realtime parameter \(p_k\) (see Canonical problem for discretetime dynamics),
both passed as column vectors. Additionally a third argument can be passed, which indicates the stage (1based) at which the function should be evaluated. Each of the generated MEX entry points output
two variables: The first is the constraint/objective evaluated at the supplied point and the second is its jacobian with repect to \(z_k\). E.g. in order to generate a solver named FORCESsolver
and generate a MEX entry point for the inequality constraints, one would proceed as follows:
model = getModel();
codeoptions = getOptions('FORCESsolver');
codeoptions.MEXinterface.inequalities = 1;
FORCES_NLP(model, codeoptions);
Now one can evaluate the inequality constraints on stage \(8\) at the point \(z_8 = (4.0, 3.0, 2.4)^T\) with realtime parameters given by \(p_8 = (3.4, 2.1)^T\) by doing
z = [4.0; 3.0; 2.4];
p = [3.4; 2.1];
ineq, jacineq = FORCESsolver_inequalities(z,p,8);
which stores the values of the inequality constraints in ineq
and their jacobian in jacineq
.
Note
If there are no realtime parameters one simply passes an empty array: ineq, jacineq = FORCESsolver_inequalities(z,[],8)
9.6.2. Calling the nonlinear functions from Python¶
When generating a solver one obtains a solver object
solver = model.generate_solver(options)
or
import forcespro.nlp
solver = forcespro.nlp.Solver.from_directory('/path/to/solver')
and one can directly evaluate any of the nonlinear functions (equality constraints, inequalities or objective function) via the solver object’s methods. The methods for doing so are
dynamics
, ineq
and objective
. Each of these functions take as a first input the primal variable \(z_k\) and the realtime parameter \(p_k\) as a second input, if the problem
formulation has realtime parameters (see Canonical problem for discretetime dynamics). All of these inputs are expected to be passed as numpy arrays. Recall that the way in which states and inputs are
“organized” into \(z_k\) is determined by model.E
. Additionally a stage
input can be given in order to compute the dynamics at a given stage (it is not necessary to supply this argument
if the same dynamics are used throughout the entire horizon). This stage
input should give the 0based number of the corresponding stage. The first output of each of these function is the evaluated constraint/objective and the second output is the derivative with respect to
\(z_k\). E.g. in order to evaluate the equality constraints at stage number \(7\) on the primal variable \(z_8 = (4.0, 3.0, 2.4)^T\) in a formulation where there are no realtime
parameters, one would do as follows:
import numpy as np
import forcespro.nlp
solver = forcespro.nlp.Solver.from_directory('/path/to/solver')
z = np.array([4.0, 3.0, 2.4])
c, jacc = solver.dynamics(z, stage=7)
This stores the value of the equality constraints in the numpy array c
and its jacobian in the numpy array jacc
.
9.7. Mixedinteger nonlinear solver¶
From FORCESPRO 1.8.0
, mixedinteger nonlinear programs (MINLPs) are supported. This broad class of problems encompasses all nonlinear programs with some integer decision variables.
This interface is provided with Variant L of FORCESPRO.
9.7.1. Writing a mixedinteger model¶
In order to use this feature, the user has to declare lower and upper bounds on all variables as parametric, as shown in the code below.
model.lb = [];
model.ub = [];
model.lbidx = range(0, model.nvar)
model.ubidx = range(0, model.nvar)
The user is then expected to provide lower and upper bounds as runtime parameters. FORCESPRO switches to the MINLP solver as soon as some variables are
declared as integers in any stage. This information can be provided to FORCESPRO via the intidx
array at every stage. An example is shown below.
%% Add integer variables to existing nonlinear model
for s = 1:5
model.intidx{s} = [4, 5, 6];
end
# Add integer variables to existing nonlinear model
for s in range(0, 5):
model.intidx[s] = [3, 4, 5]
In the above code snippet, the user declares variables 4
, 5
and 6
(3
, 4
and 5
in Python’s zerobased indexing) as integers from stage 1
to 5
(stages 0
to 4
in Python’s zerobased indexing). The values that can be taken by an integer variable are derived from its lower and upper bounds.
For instance, if the variable lies between 1
and 1
, then it can take integer values 1
, 0
or 1
. If a variable has been declared as integer and does not have lower or upper bounds, FORCESPRO raises an exception during code generation. Stating
that a variable has lower and upper bounds should be done via the arrays lbidx
and ubidx
. For instance, in the code below, variables 1
to 6
(0
to 5
in Python) in stage 1
(0
) have lower and upper bounds, which are expected to be provided at runtime.
model.lbidx{1} = 1 : 6;
model.ubidx{1} = 1 : 6;
model.lbidx[0] = range(0, 6)
model.ubidx[0] = range(0, 6)
The FORCESPRO MINLP algorithm is based on the wellknown branchandbound algorithm but comes with several customization features which generally help for improving performance on some models by enabling the user to provide application specific knowledge into the search process. At every node of the search tree, the FORCESPRO nonlinear solver is called in order to compute a solution of a relaxed problem. The generated MINLP solver code can be customized via the options described in Table 9.3, which can be changed before running the code generation.
One of the salient features of the MINLP solver is that the branchandbound search can be run in parallel on several threads. Therefore the search is split in two phases. It starts with a sequential branchandbound and switches to a parallelizable process when the number of nodes in the queue is sufficiently high. The node selection strategy can be customized in both phases, as described in Table 9.3.
Code generation setting 
Values 
Default 


Any value \(\geq 0\) 


Any value \(\geq 0\) 








Any nonnegative value preferably smaller than 8 


\(0\) or \(1\) 
\(0\) 
The
minlp.int_gap_tol
setting corresponds to the final optimality tolerance below which the solver is claimed to have converged. It is the difference between the objective incumbent, which is the best integer feasible solution found so far and the lowest lower bound. As the node problems are generally not convex, it can be expected to become negative. FORCESPRO claims convergence to a local minimum only when the integrality gap is nonnegative and below the toleranceminlp.int_gap_tol
.The
minlp.max_num_nodes
setting is the maximum number of nodes which can be explored during the search.The
minlp.seq_search_strat
setting is the search strategy which is used to select candidate nodes during the sequential search phase.The
minlp.par_search_strat
setting is the search strategy which is used to select candidate nodes during the parallelizable search phase.The
minlp.max_num_threads
setting is the maximum number of threads allowed for a parallel search. The actual number of threads on which the branchandbound algorithm can be run can be set as a runtime parameter, as described below.The
minlp.output_relaxation
setting enables users to export the primal outputs of the root relaxation. With this option set to1
, the server automatically generates one additional output for every defined output. The name of the root relaxation output is the name of the output followed by_relax
.
Note
The MINLP solver is currently constrained to run on one thread on MacOS, meaning that minlp.max_num_threads is automatically set to 1 on MacOS.
Important
When generating a MINLP solver for MacOS the thread local feature (codeoptions.threadSafeStorage) is automatically set to 0 so if a dynamic library is used for a MINLP solver in a MacOS environment then one should not run at the same time more than one solvers linked to that library. A workaround for this would be to use the static library which is not bound by this restriction.
The FORCESPRO MINLP solver also features settings which can be set at runtime. These are the following:
minlp.numThreadsBnB
, the number of threads used to parallelize the search. Its default value is 1, if not provided by the user.minlp.solver_timeout
, the maximum amount of time allowed for completing the search. Its default value is 1.0 seconds, if not set by the user.minlp.parallelStrategy
, the method used for parallelizing the mixedinteger search (from FORCESPRO 1.9.0). Value 0 (default) corresponds to a single priority queue shared between threads. Value 1 corresponds to having each thread managing its own priority queue.
9.7.2. Mixedinteger solver customization via user callbacks¶
For advanced users, the mixedinteger branchandbound search can be customized after the rounding and the branching phases. In the rounding phase, an integer feasible solution is computed after each relaxed problem solve. The user is allowed to modify the rounded solution
according to some modelling requirements and constraints. This can be accomplished via the postRoundCallback_template.c
file provided in the FORCESPRO client. This callback is applied at every stage in a loop and updates the relaxed solution stagewise. It needs
to be provided before code generation, as shown in the following code snippet.
%% Add postrounding callback to existing model
postRndCall = fileread('postRoundCallback_template.c'); % The file name can be changed by the user
model.minlpPostRounding = postRndCall;
with open('postroundCallback_template.c') as f:
model.minlpPostRounding = f.read()
The branching process can be customized in order to discard some nodes during the search. To do so, the user is expected to overwrite the file postBranchCallback_template.c
and pass it to FORCESPRO before generating the MINLP solver code.
%% Add as postbranching callbacks as you want
postBranchCall_1 = fileread('postBranchCallback_template_1.c');
postBranchCall_2 = fileread('postBranchCallback_template_2.c');
postBranchCall_3 = fileread('postBranchCallback_template_3.c');
model.minlpPostBranching{1} = postBranchCall_1;
model.minlpPostBranching{2} = postBranchCall_2;
model.minlpPostBranching{3} = postBranchCall_3;
# Add as postbranching callbacks as you want
with open('postBranchCallback_template_1.c') as f:
model.minlpPostBranching[0] = f.read()
with open('postBranchCallback_template_2.c') as f:
model.minlpPostBranching[1] = f.read()
with open('postBranchCallback_template_3.c') as f:
model.minlpPostBranching[2] = f.read()
In each of those callbacks, the user is expected to update the lower and upper bounds of the sons computed during branching given the index of the stage in which the branched variables lies, the index of this variable inside the stage and the relaxed solution at the parent node.
9.7.3. Providing a guess for the incumbent¶
Internally, the mixedinteger branchandbound computes an integer feasible solution by rounding. Moreover, since version 1.9.0
, users are allowed to provide an initial guess for the incumbent. At codegeneration, the following options need to be set:
minlp.int_guess
, which tells whether an integer feasible guess is provided by the user (value 1). Its default value is 0.minlp.int_guess_stage_vars
, which specifies the indices of the integer variables that are userinitialized within one stage (MATLAB based indexing). Ifminlp.int_guess = 1
, a parameterint_guess
needs to be set at every stage. An example can be found there Mixedinteger nonlinear solver: F8 Crusader aircraft.
Another important related option is minlp.round_root
. If set to 1
, the solution of the root relaxation is rounded and set as incumbent if feasible. Its default value is 1
. The mixedinteger solver behaviour differs depending on the combinations of options. The different behaviours are listed below.
If
minlp.int_guess = 0
andminlp.round_root = 1
, then the solution of the root relaxation is taken as incumbent (if feasible). This is the default behaviour.If
minlp.int_guess = 1
andminlp.round_root = 0
, then the incumbent guess provided by the user is tested after the root solve. If feasible, it is taken as incumbent. Note that the user is allowed to provide guesses for a few integers per stage only. In this case, the other integer variables are rounded to the closest integer.If
minlp.int_guess = 1
andminlp.round_root = 1
, then the rounded solution of the root relaxation and the user guess are compared. The best integer feasible solution in terms of primal objective is then taken as incumbent.
This feature is illustrated in Example Mixedinteger nonlinear solver: F8 Crusader aircraft. The ability of providing an integer guess for the incumbent is a key feature to run the mixedinteger solver in a receding horizon setting.
9.8. Sequential quadratic programming algorithm¶
The FORCESPRO realtime sequential quadratic programming (SQP) algorithm allows one to solve problems of the type specified in the section Highlevel Interface. The algorithm iteratively solves a convex quadratic approximations of the (generally nonconvex) problem. Moreover, the solution is stored internally in the solver and used as an initial guess for the next time the solver is called. This and other features enables the solver to have fast solvetimes (compared to the interior point method), particularly suitable for MPC applications where the sampling time or the computational power of the hardware is small.
Important
The SQP algorithm currently only supports affine inequalities. This means that all the inequality functions \(h_k, k=1,...,N\) from (9.1) must be affine functions of the variable \(z_k\) (not necessarily of \(p_k\)).
Moreover, the SQP algorithm currently does not support problems comprising final equality constraints
(specified via model.xfinalidx
).
9.8.1. How to generate a SQP solver¶
To generate a FORCESPRO sequential quadratic programming realtime iteration solver one sets
codeoptions.solvemethod = 'SQP_NLP';
codeoptions.solvemethod = "SQP_NLP"
(see Generating a solver). In addition to the general code options specified in the previous section here are some of the important code options one can use to customize the generated SQP solver.
By default the FORCESPRO SQP solver solves a single convex quadratic approximation. This accomplishes a fast solvetime compared to a “full” sequential quadratic programming solver (which solves quadratic approximations to the nonlinear program until a KKT point is reached). The user might prefer to manually allow the SQP solver to solve multiple quadratic approximations: By setting
codeoptions.sqp_nlp.maxqps = k;
codeoptions.sqp_nlp.maxqps = k
for a positive integer k
one allows the solver to solve k
quadratic approximations at every call to the solver.
In general, the more quadratic approximations which are solved, the higher the control performance. The tradeoff is that the solvetime also increases.
9.8.2. Different SQP variants¶
In addition to the default FORCESPRO SQP solver (also refered to as SQP General here), FORCESPRO supports a new SQP Fast solver from version 6.0.0. An advantage of the fast SQP solver is that it is typically faster and can be made to run on more lowlevel hardware. An advantage of the general (default) SQP solver is that it is very robust. See Tuning the SQP Fast solver for a detailed discussion on how to generate a SQP Fast solver and see Highlevel interface: Path tracking example (MATLAB) for an example of how to use the SQP Fast solver.
Note
The SQP Fast solver is only supported when using a GaussNewton hessian approximation and disabling the line search (see The hessian approximation and line search settings). The SQP Fast solver is currently only supported via the MATLAB client of FORCESPRO.
9.8.3. Tuning the SQP Fast solver¶
One of the novel features of the fast SQP solver is that it can be tuned/tailored to achieve optimal performance on a specific application. FORCESPRO provides a tool for “autotuning” the solver
(see Autotuner). The main step in using this tool is to collect training problems which can be used to tune the solver. The standard way to do this is to first generate an SQP General
(codeoptions.sqp_nlp.qp_method = "general"
, which is also the default value) solver, perform a closedloop simulation with this solver while
caching/saving all problem instances (i.e. the inputs to the solver at runtime). In a second step a SQP Fast solver is generated and automatically tuned by calling
ForcesGenerateSqpFastSolver
. The complete workflow is as follows:
tuningoptions = ForcesAutotuneOptions(problems)
ForcesGenerateSqpFastSolver(model, codeoptions, tuningoptions);
where problems
is a cell array of problems collected from the simulation with the SQP General solver. The autotuning procedure allows for quite a bit of customization which can be specified
through the ForcesAutotuneOptions
object (see Autotuner Options), and one can also pass multiple problems for the tuning process which should result in better tuned parameters.
See Highlevel interface: Path tracking example (MATLAB) for an example of how to tune the SQP Fast solver.
9.8.4. The hessian approximation and line search settings¶
The SQP code generation currently supports two different types of hessian approximations. A good choice of hessian approximation can often improve the number of iterations required by the solver and thereby its solvetime. The default option for a SQP solver is the BFGS hessian approximation. When the objective function of the optimization problem is a least squares cost it is often benefitial to use the GaussNewton hessian approximation instead. To enable this option one proceeds as specified in the sections Hessian approximation and GaussNewton options. When the GaussNewton hessian approximation is chosen one can also disable the the internal linesearch by setting
codeoptions.sqp_nlp.use_line_search = 0;
options.sqp_nlp.use_line_search = False
A linesearch is required to ensure global convergence of an SQP method, but is not needed in a realtime context when a GaussNewton hessian approximation is used.
Note
One cannot disable the line search when using the BFGS hessian approximation.
9.8.5. Controlling the initial guess at runtime¶
Upon the first call to the generated FORCESPRO SQP solver one needs to specify a primal initial guess (problem.x0
, see also Initial guess). The default behaviour of the FORCESPRO SQP
solver is to use the solution from the previous call as initial guess in every subsequent call to it. However, one can also manually set an initial guess in subsequent calls to the solver. Wether a manual initial guess
(provided through problem.x0
) will be used or the internally stored solution from the previous call will be used can be controlled by the field problem.reinitialize
of the problem
struct which is passed as an argument to the solver when it is called.
The reinitialize
field can take two values: 0 or 1. For the default usage of the solver
problem.reinitialize = 0;
problem["reinitialize"] = False
should be used. This choice results in the solver using the solution from the previous call as initial guess. This feature is useful when running the realtime iteration scheme because it ensures that the initial guess is close to the optimal solution. If you want to specify an initial guess at runtime, you will need to set
problem.reinitialize = 1;
problem["reinitialize"] = True
So in summary: The first time the solver is called the initial guess the solver will use has to be provided by problem.x0
. In all subsequent calls the solver will only make use of problem.x0
as its initial guess if problem.reinitialize = 1
.
9.8.6. Additional code options specific to the SQPRTI solver¶
In addition to the above codeoptions, the following options are specific to the SQP algorithm. Each of these options can be supplied when generating a solver as a field of codeoptions.sqp_nlp
(e.g. codeoptions.sqp_nlp.TolStat
).

Possible values 
Default value 
Description 


positive 
\(10^{6}\) 
Set the stationarity tolerance required for terminating the algorithm (the tolerance required to claim convergence to a KKT point). 

positive 
\(10^{6}\) 
Set the feasibility tolerance required for terminating the algorithm (the tolerance required to claim convergence to a feasible point). 

positive 
\(5\cdot 10^{9}\) 
Set the level of regularization of the hessian approximation (often increasing this parameter can help if the SQP solver returns exitflag 8 for your problem) 

\(0\) or \(1\) 
\(0\) 
Set the initialization strategy for the internal QP solver. \(0\) = cold start and \(1\) = centered start. See also Solver Initialization (note however, that for the SQP solver 
In addition to these options one can also specify the maximum number of iterations the internal QP solver is allowed to run in order to solve the quadratic approximation. If one wishes the QP solver use no more than k
iterations to solve a problem one sets
codeoptions.maxit = k;
Note
The SQP algorithm currently does not support parallel execution, i.e. setting
codeoptions.parallel
will have no effect.
9.9. Differences between the MATLAB and the Python client¶
The Python NLP interface is largely similar to the MATLAB interface, but does come with some language and implementationspecific differences.
All indices in the problem formulation are expected to be 0based in Python, as is usual in this language. This does not include the indices of the generated solver, however, where outputs are named x01, x02, … as in MATLAB. Thus, the problem formulation before generation requires 0based indices, whereas the returned solver from the server uses 1based indices. This also does not apply to the lowlevel Python interface, where indices are 1based even in the model formulation.
In the Python client, different model objects must be used when using external functions or symbolic expressions, namely nlp.ExternalFunctionModel() and nlp.SymbolicModel(). Furthermore, if the highlevel interface is to be used for convex problems, this is only possible using the nlp.ConvexSymbolicModel(). This is different from the MATLAB client, where the FORCES_NLP function accepts problems of any kind and switches to the appropriate solver automatically.
When using the Python client with a nlp.SymbolicModel(), the C code generated for symbolic expressions is currently not entirely identical to the code generated by MATLAB. While the actual expression evaluation code generated by CasADi is the same, the structure of the files varies. Specifically, the MATLAB client creates individual C files for each problem stage with distinct symbolic expressions (leading to varying file names when changing the problem horizon) whereas all functions are gathered in one file in the Python client. Yet, the Python client does add one additional file for the FORCESPROCasADi glue code, which is not present when using the MATLAB client. Lastly, function names of the evaluation functions differ.
If you want to get the same code for MATLAB and Python, you must generate the CasADi C code from one of both clients and then supply this code as an external function in the other client.
9.10. Examples¶
Highlevel interface: Basic example: In this example, you learn the basics in how to use FORCESPRO to create an MPC regulation controllers.
Highlevel interface: Obstacle avoidance (MATLAB & Python): This example uses a simple nonlinear vehicle model to illustrate the use of FORCESPRO for realtime trajectory planning around nonconvex obstacles.
Highlevel interface: Indoor localization (MATLAB & Python): This examples describes a nonlinear optimization approach for the indoor localization problem.
Mixedinteger nonlinear solver: F8 Crusader aircraft: In this example, you learn the basics in how to use FORCESPRO MINLP solver to solve a mixedinteger optimal control problem.
Realtime SQP Solver: Robotic Arm Manipulator (MATLAB & Python): This example describes how to apply the FORCESPRO SQP solver to control a robotic arm.
Controlling a DC motor using a FORCESPRO SQP solver: This example describes how to apply the FORCESPRO SQP solver to control a DC motor.
Controlling a crane using a FORCESPRO NLP solver: This example describes how to apply the FORCESPRO interior point NLP solver to control a crane.
Highlevel interface: Optimal EV charging and speed profile example (MATLAB & PYTHON): This example describes how to apply the FORCESPRO interior point NLP solver to control the speed and charging profile of an electric vehicle.