8. Low-level interface

FORCESPRO supports designing solvers and controllers via MATLAB and Python scripts. When using the MATLAB client, a Simulink block is always created such that you can plug your advanced formulation directly into your simulation models, or download it to a real-time target platform.

The low-level interface gives advanced optimization users the full flexibility when designing custom optimization solvers and MPC controllers based on non-standard formulations.

This interface is provided with all variants of FORCESPRO, starting with Variant S.

8.1. Supported problem class

The FORCESPRO low-level interface supports the class of convex multistage quadratically constrained programs (QCQPs) of the form

\begin{align*} \text{minimize} \quad & \sum_{i=1}^N \frac{1}{2}z_i^{\top}H_iz_i + f_i^{\top}z_i && \quad \text{(separable objective)} \\ \text{subject to} \quad & D_1z_1 = c_1 && \quad \text{(initial equality)} \\ & C_{i-1}z_{i-1} + D_iz_i = c_i,~~i=2,\ldots,N && \quad \text{(inter-stage equality)} \\ & \underline{z}_i \leq z_i && \quad \text{(lower bound)} \\ & z_i \leq \bar{z}_i && \quad \text{(upper bound)} \\ & A_iz_i \leq b_i && \quad \text{(polytopic inequalities)} \\ & z_i^{\top}Q_{i,k}z_i + L_{i,k}^{\top}z_i \leq r_{i,k} && \quad \text{(quadratic inequalities)} \end{align*}

for \(i=1,...,N\) and \(k=1,...,M\). To obtain a solver for this optimization program using the FORCESPRO client, you need to define all data in the problem, that is the matrices \(H_i\), \(A_i\), \(Q_{i,j}\), \(D_i\), \(C_i\) and the vectors \(\underline{z}_i < \bar{z}_i\), \(b_i\), \(L_{i,k}\), \(r_{i,k}\), \(c_i\), in a MATLAB struct or Python dictionary, along with the corresponding dimensions. The following steps will take you through this process. Importantly, the matrices \(H_i\) and \(Q_{i,j}\) should all be positive definite.

Note

FORCESPRO supports all problem data to be parametric, i.e. to be unknown at code generation time. Read Section 12 to learn how to use parameters correctly.

In the following, we describe how to model a problem of the above form with FORCESPRO. First make sure that the FORCESPRO client is on the MATLAB/Python path. See Section 3 for more details on how to set up the MATLAB client and Section 3.3.

After the PYTHONPATH has been appropriately set up to include your FORCESPRO client directory (see Section 3.3.3), Python users have to import the FORCESPRO module and their user ID.

import forcespro
import get_userid

8.2. Multistage struct

First, an empty struct/class has to be initialized, which contains all fields needed and initialises matrices and vectors to empty matrices. The command

stages = MultistageProblem(N);
stages = forcespro.MultistagePoblem(N) # 0-indexed

creates such an empty structure/class of length \(N\). Once this structure/class has been created, the corresponding matrices, vectors and dimensions can be set for each element of stages.

8.3. Dimensions

In order to define the dimensions of the stage variables \(z_i\), the number of lower and upper bounds, the number of polytopic inequality constraints and the number of quadratic constraints use the following fields:

stages(i).dims.n = ...; % length of stage variable zi
stages(i).dims.r = ...; % number of equality constraints
stages(i).dims.l = ...; % number of lower bounds
stages(i).dims.u = ...; % number of upper bounds
stages(i).dims.p = ...; % number of polytopic constraints
stages(i).dims.q = ...; % number of quadratic constraints
# 0-indexed
stages.dims[i]['n'] = ... # length of stage variable zi
stages.dims[i]['r'] = ... # number of equality constraints
stages.dims[i]['l'] = ... # number of lower bounds
stages.dims[i]['u'] = ... # number of upper bounds
stages.dims[i]['p'] = ... # number of polytopic constraints
stages.dims[i]['q'] = ... # number of quadratic constraints

8.4. Cost function

The cost function is, for each stage, defined by the matrix \(H_i\) and the vector \(f_i\). These can be set by

stages(i).cost.H = ...; % Hessian
stages(i).cost.f = ...; % linear term
# 0-indexed
stages.cost[i]['H'] = ... # Hessian
stages.cost[i]['f'] = ... # linear term

Note: whenever one of these terms is zero, you have to set them to zero (otherwise the default of an empty matrix is assumed, which is different from a zero matrix).

H_i should be square of size dims[i]['n'] / stages(i).dims.n and f_i should also be of that length. It does not matter whether you use a column or row vector for f.

8.5. Equality constraints

The equality constraints for each stage, which are given by the matrices \(C_{i-1}\), \(D_i\) and the vector \(c_i\), have to be provided in the following form:

The matrices \(C_{i-1}, D_i\) should have dimension dims[i]['r'] by dims[i]['n'] and \(c_i\) should be of length dim[i]['r'].

Note the index shift in \(C_{i-1}\). In particular you should take care, that the vertical dimension of \(C_{i-1}\) matches dim[i]['r'] and not i-1.

stages(i).eq.C = ...;
stages(i).eq.c = ...;
stages(i).eq.D = ...;
# 0-indexed
stages.eq[i]['C'] = ...
stages.eq[i]['c'] = ...
stages.eq[i]['D'] = ...

In many parts of the literature a different notation is given for inter-stage equality, which places the next index on the right hand side of the equation like so:

\[A_{i}x_{i} + B_iu_i = x_{i+1}, \quad i=1,\ldots,N-1 \quad \text{(inter-stage equality)}\]

with states \(x_{i}\) and control inputs \(u_{i}\).

The correct way to translate this to FORCESPRO is as follows:

\[\left(\begin{array}{@{}c|c@{}} B_{i-1} & A_{i-1} \end{array}\right) z_{i-1} + \left(\begin{array}{@{}c|c@{}} 0 & -\mathbf{Id} \end{array}\right) z_i = 0, \quad i=2, \dots, N\]

where \(z_i = (u_i, x_i)\). Note that the number of columns of \(-\mathbf{Id}\) must match the size of \(x_i\).

There is the common case of a fixed initial value \(x_1 = x_{init}\). This can be implemented by adding an equality constraint for \(i=1\):

\[\left(\begin{array}{@{}c|c@{}} 0 & -\mathbf{Id} \end{array}\right) z_1 = -x_{init}\]

Alternatively, one stage can be eliminated by defining the state space as \(z_i = (u_i, x_{i+1})\) for \(i=1,\dots N-1\) (thus removing \(x_1\) and \(u_N\)). For that choice, the inter-stage equalities are given by:

\begin{align*} & \left(\begin{array}{@{}c|c@{}} B_{1} & -\mathbf{Id} \end{array}\right) z_1 = -A_1 x_{init} \\ & \left(\begin{array}{@{}c|c@{}} 0 & A_{i} \end{array}\right) z_{i-1} + \left(\begin{array}{@{}c|c@{}} B_{i} & -\mathbf{Id} \end{array}\right) z_i = 0 \quad i=2,\dots, N - 1 \end{align*}

8.6. Lower and upper bounds

Lower and upper bounds have to be set in sparse format, i.e. an index vector lbIdx/ubIdx that defines the elements of the stage variable \(z_i\) has to be provided, along with the corresponding upper/lower bound lb/ub:

These index vectors will be the same in both Matlab and Python, which means that the Python indices need to be adjusted to match Matlab’s 1-indexed style.

stages(i).ineq.b.lbidx = ...; % index vector for lower bounds
stages(i).ineq.b.lb = ...;    % lower bounds
stages(i).ineq.b.ubidx = ...; % index vector for upper bounds
stages(i).ineq.b.ub = ...;    % upper bounds
stages.ineq[ i ]['b']['lbidx'] = ... # index vector for lower bounds, 1-indexed
stages.ineq[ i ]['b']['lb'] = ...    # lower bounds
stages.ineq[ i ]['b']['ubidx'] = ... # index vector for upper bounds, 1-indexed
stages.ineq[ i ]['b']['ub'] = ...    # upper bounds

Both lb and lbIdx must have length stages(i).dims.l / stages.dims[ i ]['l'], and both ub and ubIdx must have length stages(i).dims.u / stages.dims[ i ]['u'].

8.7. Polytopic constraints

In order to define the inequality \(A_iz_i\leq b_i\), use

stages(i).ineq.p.A = ...; % Jacobian of linear inequality
stages(i).ineq.p.b = ...; % RHS of linear inequality
# 0-indexed
stages.ineq[i]['p']['A'] = ... # Jacobian of linear inequality
stages.ineq[i]['p']['b'] = ... # RHS of linear inequality

The matrix A must have stages(i).dims.p / stages.dims[ i ][‘p’] rows and stages(i).dims.n / stages.dims[ i ][‘n’] columns. The vector b must have stages(i).dims.p / stages.dims[ i ][‘p’] rows.

8.8. Quadratic constraints

Similar to lower and upper bounds, quadratic constraints are given in sparse form by means of an index vector, which determines on which variables the corresponding quadratic constraint acts.

stages(i).ineq.q.idx = { idx1, idx2,}; % index vectors
stages(i).ineq.q.Q = { Q1, Q2,};       % Hessians
stages(i).ineq.q.l = { L1, L2,};       % linear terms
stages(i).ineq.q.r = [ r1; r2;];      % RHSs
stages.ineq[i]['q']['idx'] = ... # index vectors, 1-indexed
stages.ineq[i]['q']['Q'] = ...   # Hessians
stages.ineq[i]['q']['l'] = ...   # linear terms
stages.ineq[i]['q']['r'] = ...   # RHSs

If the index vector idx1 has length \(m_1\), then the matrix Q must be square and of size \(m_1\times m_1\), the column vector l1 must be of size \(m_1\) and r_1 is a scalar. Of course this dimension rules apply to all further quadratic constraints that might be present. Note that L_1, L_2 etc. are column vectors in MATLAB!

Since multiple quadratic constraints can be present per stage, in MATLAB we make use of the cell notation for the Hessian, linear terms, and index vectors. In Python we make use of Python object arrays for the Hessians, linear terms, and index vectors.

8.8.1. Example

To express the two quadratic constraints

\begin{align*} & z_{3,3}^2 + 2z_{3,5}^2 - z_{3,5} \leq 3 \\ & 5z_{3,1}^2 \leq 1 \end{align*}

on the third stage variable, use the code

stages(3).ineq.q.idx = { [3 5], [1] } % index vectors
stages(3).ineq.q.Q = { [1 0; 0 2], [5] }; % Hessians
stages(3).ineq.q.l = { [0; -1], [0] }; % linear terms
stages(3).ineq.q.r = [ 3; 1 ]; % RHSs
stages.ineq[3-1]['q']['idx'] = np.zeros((2,), dtype=object) # index vectors
stages.ineq[3-1]['q']['idx'][0] = np.array([3,5])
stages.ineq[3-1]['q']['idx'][1] = np.array([1])
stages.ineq[3-1]['q']['Q'] = np.zeros((2,), dtype=object) # Hessians
stages.ineq[3-1]['q']['Q'][0] = np.array([1.0 0],[0 2.0])
stages.ineq[3-1]['q']['Q'][1] = np.array([5])
stages.ineq[3-1]['q']['l'] = np.zeros((2,), dtype=object) # linear terms
stages.ineq[3-1]['q']['l'][0] = np.array([0], [-1])
stages.ineq[3-1]['q']['l'][1] = np.array([0])
stages.ineq[3-1]['q']['r'] = np.array([3],[1]) # RHSs

8.9. Binary constraints

To declare binary variables, you can use the bidx field of the stages struct or object. For example, the following code declares variables 3 and 7 of stage 1 to be binary:

stages(1).bidx = [3 7]
stages.bidx[0] = np.array([3, 7])

That’s it! You can now generate a solver that will take into account the binary constraints on these variables. If binary variables are declared, FORCESPRO will add a branch-and-bound procedure to the standard convex solver it generates.

8.10. Declaring parameters

FORCESPRO is a parametric solver. As such, it is necessary to specify which parts of the model are to be parametric.

For a detailed introduction to setting parameters, see the section Section 12.

A common choice of parameter is an \(x_{init}\) value. To set this, use the following:

parameter = newParam('xinit', 1, 'eq.c')
stages.newParam("xinit", [1], 'eq.c') # 1-indexed

8.11. Declaring Solver Outputs

FORCESPRO gives you full control over the part of the solution that should be outputted by the solver. It is also possible to obtain the Lagrange multipliers of certain constraints. To define a standard output as a slice of the primal solution vector, use the function

output = newOutput(name, maps2stage, idxWithinStage)
stages.newOutput(name, maps2stage, idxWithinStage)

where name is the name you give to the output (you will need this to read it after calling the solver). The index vector (or integer) maps2stage defines to which stage this output maps to. The last argument, idxWithinStage allows the user to select which indices from the stage vector should be outputted by the solver.

To define an output as a slice of certain Lagrange multipliers, use the function

output = newOutput(name, maps2stage, idxWithinStage, maps2const)
stages.newOutput(name, maps2stage, idxWithinStage, maps2const)

where the remaining argument maps2const specifies the constraint associated with the Lagrange multipliers being requested.

Table 8.1 Possible string values for argument maps2const

maps2const

Constraint

r

Equalities

u

Upper bounds

l

Lower bounds

p

Polytopic bounds

8.11.1. Example

To define an output to be the first two elements of the primal solution vector, use the following command:

output1 = newOutput('u0', 1, 1:2)
stages.newOutput('u0', 1, range(1,3))

To define an output to be the first and third indices of the Lagrange multipliers for the equality constraints of the second stage, use the following command:

output2 = newOutput('dual_eq0', 2, [1 3], 'r')
stages.newOutput('dual_eq0', 2, [1,3], 'r')

8.12. Dumping and Loading Models

To dump a low-level formulation, simply run:

not implemented
save_all(stages)

Available options are listed in detail in Section 20.

Note that, unlike in the high-level implementation, specifying options or outputs does nothing and is ignored, as these should be included in the stages object.

Note

In the Python client, code generation and dumping are not interchangeable, as code generation mutates the CodeOptions object.

To get around this issue, use Python’s copy.deepcopy function to make a copy which will not be affected by code generation.

Note, that calling copy.copy makes a shallow copy and will contain references to the original object and will thus also mutate.

8.13. Generating the solver

After the optimization problem has been formulated into a structure stages, an optimized solver can be generated. To do so, the solver requires a name and a number of solver options, as described in Section 17.

codeoptions = getOptions('solver name');
generateCode(stages, params, codeoptions, outputs);
options = forcespro.CodeOptions('solver_name')
stages.codeoptions = options
stages.generateCode(get_userid.user_id)

8.14. Calling the generated low-level solver

After solver generation has completed, the solver itself (as a compiled library) as well as several interfacing files will become available in your working directory. These files are named according to what you named your solver; in the following we assume “SOLVER_NAME”. Calling the solver from MATLAB or Python is then as simple as:

problem = {}  % a struct of solver parameters
SOLVER_NAME(problem)
import SOLVER_NAME_py  # notice the _py suffix
problem = {}  # a dictionary of solver parameters
SOLVER_NAME_py.SOLVER_NAME_solve(problem)

Note

Don’t give your solver the same name as the script you are calling it from. Doing so will overwrite your calling script with the solver interface. For example, in a script named test_problem.m, choose a name such as test_solver instead of test_problem.

Note

The high-level Python interface provides more convenient access to solvers generated using the high-level interface. This method of calling a solver is only available for solvers generated through the low-level interface, and high-level solvers can only be called from Python through the means described in the high-level interface documentation.

8.15. The QP_FAST algorithm

From FORCESPRO version 6.0.0 a new algorithm was introduced which allows one to generate extremely fast solvers which are especially well-suited for low-level hardware. Generating such a solver is done by specifying the model through the stages object as explained above. The main difference in the solver generation procedure is a tuning step which is explained in details in the following section.

Important

Currently the QP_FAST algorithm is only supported through the MATLAB client of FORCESPRO.

8.15.1. Tuning the QP_FAST algorithm

One of the novel features of the QP_FAST algorithm is that it can be tailored to a specific application by “tuning” it. FORCESPRO supports an automated tuning tool (see Autotuner) for choosing the optimal tuning for a given application. A key step in tuning a fast QP solver is collecting data/problems on which it can be tuned. The way to do this is to first generate a general QP (FORCESPRO) solver which does not require tuning. This is done as follows:

codeoptions.solvemethod = 'PDIP';
generateCode(stages, params, codeoptions, outputs);
# not supported

Now one or more simulations can be performed and all the problems (i.e. the inputs to the solver at run-time) which are solved should be collected and stored in a cell array (here called problems). We refer to problems as the tuning data. Once tuning data has been collected, a fast QP solver can be generated and tuned as follows:

tuningoptions = ForcesAutotuneOptions(problems);
ForcesGenerateQpFastSolver(stages,params,codeoptions,tuningoptions,outputs);
# not supported

The autotuning procedure allows for quite a bit of customization which can be specified through the ForcesAutotuneOptions object (see Autotuner Options).

You can download an example displaying the full workflow by clicking here.

8.15.2. The QP_FAST options

Options specific to the QP_FAST algorithm are specified via the codeoptions.qp_fast field (e.g. codeoptions.qp_fast.warmstart = 0;). The following options can be set:

  • warmstart: Set equal to 1 (default) in order to allow run-time warmstart. Set equal to 0 to not allow run-time warmstart. The run-time warmstart option creates a run-time parameter (problem.warmstart) which if set to 1 will use the solution of the previous call to the solver as initial guess. If set to 0 the solver will use a default initial guess.

  • tol_primal: Set the tolerance for the primal residual (default is 1e-3).

  • tol_dual: Set the tolerance for the dual residual (default is 1e-3).