5. Lowlevel interface¶
FORCES PRO supports designing solvers and controllers via MATLAB and Python scripts. Regardless of the language used, a Simulink block is always created such that you can plug your advanced formulation directly into your simulation models, or download it to a realtime target platform.
The lowlevel interface gives advanced optimization users the full flexibility when designing custom optimization solvers and MPC controllers based on nonstandard formulations.
5.1. Supported problem class¶
The FORCES PRO lowlevel interface supports the class of convex multistage quadratically constrained programs (QCQPs) of the form
for \(i=1,...,N\) and \(k=1,...,M\). To obtain a solver for this optimization program using the FORCES PRO 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
FORCES PRO supports all problem data to be parametric, i.e. to be unknown at code generation time. Read Section 9 to learn how to use parameters correctly.
5.2. Expressing the optimization problem in MATLAB or Python¶
In the following, we describe how to model a problem of the above form with FORCES PRO. First make sure that the FORCES PRO client is on the MATLAB/Python path. See Section 2 for more details.
Currently, Python versions 2.7, 3.4 and 3.6 are supported. To use the FORCES PRO client, Python needs the following packages:
suds (Python2) / sudsjurko (Python3)
scipy (also numpy if not installed with scipy)
matplotlib (needed for the examples)
Python users first have to import the FORCES PRO module.
import sys
if sys.version_info[0] == 3:
if sys.version_info[1] == 6:
from forcespro36 import *
elif sys.version_info[1] == 4:
from forcespro34 import *
else:
from forcespro3 import *
else:
from forcespro import *
5.3. 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);
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.
5.4. 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
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
5.5. 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
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).
5.6. Equality constraints¶
The equality constraints for each stage, which are given by the matrices \(C_i\), \(D_i\) and the vector \(c_i\), have to be provided in the following form:
stages(i).eq.C = ...;
stages(i).eq.c = ...;
stages(i).eq.D = ...;
stages.eq[ i ]['C'] = ...
stages.eq[ i ]['c'] = ...
stages.eq[ i ]['D'] = ...
5.7. 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:
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
stages.ineq[ i ]['b']['lb'] = ... # lower bounds
stages.ineq[ i ]['b']['ubidx'] = ... # index vector for upper bounds
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’].
5.8. 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
stages.ineq[ i ]['A'] = ... # Jacobian of linear inequality
stages.ineq[ i ]['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.
5.9. 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
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.
5.9.1. Example¶
To express the two quadratic constraints
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[31]['q']['idx'] = np.zeros((2,), dtype=object) % index vectors
stages.ineq[31]['q']['idx'][0] = np.array([3,5])
stages.ineq[31]['q']['idx'][1] = np.array([1])
stages.ineq[31]['q']['Q'] = np.zeros((2,), dtype=object) % Hessians
stages.ineq[31]['q']['Q'][0] = np.array([1.0 0],[0 2.0])
stages.ineq[31]['q']['Q'][1] = np.array([5])
stages.ineq[31]['q']['l'] = np.zeros((2,), dtype=object) % linear terms
stages.ineq[31]['q']['l'][0] = np.array([0], [1])
stages.ineq[31]['q']['l'][1] = np.array([0])
stages.ineq[31]['q']['r'] = np.array([3],[1]) % RHSs
5.10. 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[1] = 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, FORCES PRO will add a branchandbound procedure to the standard convex solver it generates.
5.11. Declaring Solver Outputs¶
FORCES PRO 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 in MATLAB, use the function
output = newOutput(name, maps2stage, idxWithinStage)
To define a standard output as a slice of the primal solution vector in Python, use the function
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 in MATLAB, use the function
output = newOutput(name, maps2stage, idxWithinStage, maps2const)
To define an output as a slice of certain Lagrange multipliers in Python, use the function
stages.newOutput(name, maps2stage, idxWithinStage, maps2const)
where the remaining argument maps2const
specifies the constraint associated with the Lagrange multipliers being requested.
maps2const 
Constraint 


Equalities 

Upper bounds 

Lower bounds 

Polytopic bounds 
5.11.1. Example¶
To define an output to be the first two elements of the primal solution vector in MATLAB, use the following command:
output1 = newOutput('u0', 1, 1:2)
To perform the same action in Python, use the following command:
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 in MATLAB, use the following command:
output2 = newOutput('dual_eq0', 2, [1 3], 'r')
To perform the same action in Python, use the following command:
stages.newOutput('dual_eq0', 2, list(1,3), 'r')