15. Solver Options¶
The default solver options can be loaded when giving a name to the solver with the following command
codeoptions = getOptions('solvername');
stages.codeoptions = forcespro.CodeOptions('solvername') # for the lowlevel interface
codeoptions = forcespro.CodeOptions('solvername') # for the highlevel interface
In the documentation below, we assume that you have created this struct and
named it codeoptions
.
Note
For the lowlevel interface in Python, the codeoptions
struct
has to be an element of the stages
struct.
15.1. General options¶
We will first discuss how to change several options that are valid for all the FORCESPRO interfaces.
15.1.1. Solver name¶
The name of the solver will be used to name variables, functions, but also the MEX file and associated help file. This helps you to use multiple solvers generated by FORCESPRO within the same software project or Simulink model. To set the name of the solver use:
codeoptions.name = 'solvername';
codeoptions.name = 'solvername'
Alternatively, you can directly name the solver when generating the options struct by calling:
codeoptions = getOptions('solvername');
stages.codeoptions = forcespro.CodeOptions('solvername') # for the lowlevel interface
codeoptions = forcespro.CodeOptions('solvername') # for the highlevel interface
15.1.2. Print level¶
To control the amount of information the generated solver prints to the console, set the field printlevel
as outlined in Table 15.1.

Result 
Dependency 


No output will be written. 
(None) 

Summary line after each solve. 


Summary after each iteration of solver. 

Note
For printlevel=0
, the generated solver has no dependency on any system
library. Otherwise, there will be a dependency on <stdio.h>
.
Important
printlevel
should always be set to 0
when recording performance
timings or when deploying the code on an autonomous embedded system.
15.1.3. Maximum number of iterations¶
To set the maximum number of iterations of the generated solver, use:
codeoptions.maxit = 200;
codeoptions.maxit = 200
The default maximum number of iterations for all solvers provided by FORCESPRO
is set to 200
.
15.1.4. Compiler optimization level¶
The compiler optimization level can be varied by changing the field optlevel
from 0
to 3
(default):
codeoptions.optlevel = 0;
codeoptions.optlevel = 0
Important
It is recommended to set optlevel
to 0
during prototyping to evaluate
the functionality of the solver without long compilation times. Then set it
back to 3
when generating code for deployment or timing measurements.
15.1.5. Measure Computation time¶
You can measure the time used for executing the generated code by using:
codeoptions.timing = 1;
codeoptions.timing = 1
By default the execution time is measured. The execution time can be accessed in
the field solvetime
of the information structure returned by the solver. In
addition, the execution time is printed in the console if the flag
printlevel
is greater than 0
.
Important
Setting timing on will introduce a dependency on libraries used for accessing the system clock. Timing should be turned off when deploying the code on an autonomous embedded system.
By default when choosing to generate solvers for target platforms, timing is disabled. You can manually enable timing on embedded platforms by using:
codeoptions.embedded_timing = 1;
codeoptions.embedded_timing = 1
15.1.6. Solver Timeout¶
15.1.6.1. Introduction¶
If you have a critical application which needs to run in a specific timeframe then it’s useful to set a timeout for the solver in order to control its execution time.
The timeout works by checking the execution time of each iteration of the solver and making an estimate for next iterations as:
next_iteration_time = timeout_estimate_coeff * max_iteration_time
where:
max_iteration_time
is the execution time of the currently slowest iterationtimeout_estimate_coeff
is a coefficient used to make the estimate more conservative or forgiving. Its default value is1.20
15.1.6.2. Usage¶
To enable the solver timeout you can use the following codeoption:
% solver_timeout can take values 02
codeoptions.solver_timeout = 1;
# solver_timeout can take values 02
codeoptions.solver_timeout = 1
Setting the option to 1
will enable the timeout and provide the floating point variable solver_timeout
as a runtime parameter.
Setting the option to 2
will additionally provide the floating point variable timeout_estimate_coeff
as a runtime parameter.
Important
For MINLP solvers a timeout is automatically enabled therefore there’s no need to use the above codeoptions.
For more details on how to use it please check section Mixedinteger nonlinear solver
.
Not setting the runtime parameters after enabling them with code generation will result in them taking their default values. The default values for the runtime parameters are:
For
solver_timeout
it’s1.0
which results in timeout being disabledFor
timeout_estimate_coeff
it’s1.20
Important
Since an estimation is required for the timeout, the solvers will always perform
the first iteration (only exception are SQP methods, check the following section
SQP inner QP timeout
).
15.1.6.3. SQP inner QP timeout¶
With the SQP_NLP solve method the QP solved as part of the SQP iteration is also set to timeout based on the remaining time available to the SQP solver. The QP timeout can be useful in cases where the inner QP takes longer time to execute than expected and could otherwise cause the SQP solver to miss the timeout mark (in which case the SQP solver would time out at the start of the next iteration). If the QP times out, the SQP solver will return with the solution from the previous iteration.
If it is deemed more important to solve the whole QP and get a more updated solution rather than having a strict timeout, the inner qp timeout can be disabled with the following codeoption:
% this option is relevant only if codeoptions.solver_timeout is enabled
codeoptions.sqp_nlp.qp_timeout = 0;
# this option is relevant only if codeoptions.solver_timeout is enabled
codeoptions.sqp_nlp.qp_timeout = 0
15.1.6.4. Return Value¶
When solver timeout is enabled, two additional exitflags are available for the user:






The solver timed out and returned the solution found up to the executed iteration 


The timeout provided was too small to even start a single iteration 
If a normal timeout is returned, the outputs of the solver will contain the solution found up to the executed iteration. If an invalid timeout is returned, the outputs of the solver will contain the initialization of the solver (or the previous solution if it exists for SQPs).
15.1.7. Running multiple instances of the same solver¶
Multiple instances of the generated solver can be run in a multithreaded
environment by making use of the expert
C interface. The expert interface
gives the user control over the allocation of the required memory so that a
different block of memory can be preallocated for each thread. This feature is
currently implemented for the algorithms PDIP
and PDIP_NLP
.
The threadsafe expert interface can be requested by setting the option
codeoptions.threadSafeExpert = 1;
codeoptions.threadSafeExpert = 1
and consists of two additional C functions:
<solvername>_solver_mem <solvername>_mem_new()
provides a handle to the memory to be used by a solver.<solvername>_solve_expert(..., <solvername>_solver_mem * mem, ...)
is the solver function (instead of<solvername>_solve(...)
) taking the memory handle as an additional input argument.
In order to ensure threadsafety, each thread must be assigned its own memory handle.
Below an example of how to run multiple solvers in parallel using OpenMP:
// prepare 4 solvers
#define nSolvers = 4
FORCES_NLP_solver_mem mem[nSolvers];
FORCES_NLP_solver_params params[nSolvers];
FORCES_NLP_solver_output output[nSolvers];
FORCES_NLP_solver_info info[nSolvers];
// ... initialize params ...
// initialize memory for each thread
int nThreads = omp_get_num_threads();
int iThread;
for (iThread=0; iThread<nThreads; iThread++)
{
mem[iThread] = FORCES_NLP_solver_mem_new();
}
// run the solvers in parallel
int exitflag;
#pragma omp parallel for
for (iSolver=0; iSolver<nSolvers; iSolver++)
{
int iThread = omp_get_thread_num();
exitflag = FORCES_NLP_solver_solve_expert(¶ms[iSolver], &output[iSolver], &info[iSolver], &mem[iThread], ...);
}
You can find the full code of this example including instructions of how to run
it in the Examples\ThreadSafe
folder that comes with your client.
Alternatively (and for the algorithms not yet supported with threadSafeExpert
= 1
), use the option threadSafeStorage = 1
. This ensures threadsafety
within the default interface by relying on threadlocal storage.
Due to the limited amount of threadlocal storage (depending on the system),
using threadSafeExpert
is recommended over using threadSafeStorage
and threadSafeStorage
will be deprecated in the near future.
Important
When using the codegenerated integrators (see section Codegenerated integrators) with the threadSafeExpert
option enabled, you will have to specify via
the option nlp.max_num_threads
the maximum number of threads on which you wish to run the solver in parallel. For instance, if running the solver on a maximum of 5 threads in parallel one would set
codeoptions.nlp.max_num_threads = 5;
codeoptions.nlp.max_num_threads = 5
Important
The deployment of a Simulink model to a target platform is not supported when the threadSafeExpert option is enabled. To deploy a Simulink model, the threadSafeExpert option must be unset or disabled.
15.1.8. Datatypes¶
The type of variables can be changed by setting the field floattype
as outlined in Table 15.3.
This will effect all floating point variables used inside the solver and the callbacks
generated by the AD tool.

Decimation 
Width (bits) 
Supported algorithms 


64 bit 
Floating point 
PDIP_NLP, SQP_NLP, PDIP, ADMM, DFG, FG 

32 bit 
Floating point 
PDIP, ADMM, DFG, FG 

32 bit 
Fixed point 
ADMM, DFG, FG 

16 bit 
Fixed point 
ADMM, DFG, FG 
Important
Unless running on a resourceconstrained platform, we recommend using double
precision floating point arithmetic to avoid problems in the solver. If
single precision floating point has to be used, reduce the required tolerances
on the solver accordingly by a power of two (i.e. from 1E6
to 1E3
).
When running the solver in double precision arithmetic, it is possible to only use single
precision arithmetic for evaluating the AD tool callbacks. This can be done by setting
the field callback_floattype
; see Table 15.4 and section
Single precision callbacks for details.

Decimation 
Width (bits) 
Supported algorithms 


64 bit 
Floating point 
PDIP_NLP, SQP_NLP 

32 bit 
Floating point 
PDIP_NLP, SQP_NLP 
15.1.9. Code generation server¶
By default, code generation requests are routed to embotech’s server. To send a code generation request to a local server, for example when FORCESPRO is used in an enterprise setting, set the following field to an appropriate value:
codeoptions.server = 'https://yourforcesserver.com:1234';
codeoptions.server = 'https://yourforcesserver.com:1234'
15.1.10. Enforcing solver regeneration¶
In order to avoid unnecessary calls to the codegeneration server, FORCESPRO internally computes a hash of your problem formulation and codeoptions. If this hash is identical to that of an already generated solver, the existing one is reused. In situations where this is not desired, hashing can be disabled as follows:
codeoptions.nohash = 1;
codeoptions.nohash = 1
In that case, the codegen server is always contacted to regenerate a new solver.
15.1.11. Overwriting existing solvers¶
When a new solver is generated with the same name as an existing solver one can
control the overwriting behaviour by setting the field overwrite
as outlined
in Table 15.5.

Result 


Never overwrite. 

Always overwrite. 

Ask to overwrite. 
15.1.12. Skipping the Build of Simulink Sfunction¶
By default, after code generation, the Simulink block is compiled, which may take a very long time for large problems on Windows systems. If you will not use the Simulink block, or want to build it later yourself, you can disable automatic builds by using the following option:
codeoptions.BuildSimulinkBlock = 0;
# does not take effect in Python
15.1.13. Solver info in Simulink block¶
FORCESPRO always generates a Simulink block encapsulating the generated solver. You can add output ports to the Simulink block to obtain the solver exit flag and other solver information (number of iterations, solve time in seconds, value of the objective function) by setting:
codeoptions.showinfo = 1;
codeoptions.showinfo = 1
By default these ports are not present in the Simulink block.
15.1.14. Skipping automatic cleanup¶
FORCESPRO automatically cleans up some of the files that it generates during the code generation, but which are usually not needed any more after building the MEX file. In particular, some intermediate CasADi generated files are deleted. If you would like to prevent any cleanup by FORCES, set the option:
codeoptions.cleanup = 0;
codeoptions.cleanup = 0
The default value is 1
(true).
Important
The library or object files generated by FORCESPRO contain only the solver itself. To retain the CasADi generated files for function evaluations, switch off automatic cleanup as shown above. This is needed if you want to use the solver within another software project, and need to link to it.
15.1.15. MATLAB network communications¶
From version 4.2.0
, a new method is used for the MATLAB client for communicating with the
FORCESPRO codegen server that adapts to MATLAB’s newer implementations.
To revert to the old method, either set:
codeoptions.legacyNetworkConnections = 1;
or disable it by editing the FORCESPRO client.
To do so, please edit the +ForcesWeb/useLegacyVersion.m
function so that it returns true
.
Important
For the client edit change to have an effect codeoptions.legacyNetworkConnections
should not be set.
15.1.16. Python network communications¶
From version 4.3.1
, the Python client supports connections to the FORCESPRO codegen server
through a proxy.
The file forcespro_proxy.py
exists in the FORCESPRO client folder in order to set the configuration for the proxy.
The format of the file is as follows:
# host of the proxy. Can be an IP address ("x.x.x.x") or a DNS record. Set to empty to not use a proxy
host=""
# port number of proxy to connect to. To use default set to 0
port=8888
# Protocol to connect to the proxy (http or https). To use default set to empty
protocol="http"
# Username with which to connect to the proxy. To not use a username set to empty
username="user"
# Password with which to connect to the proxy. To not use a password set to empty
password="pass"
Note
By default the file forcespro_proxy.py
has an empty host entry so that no proxy is used unless set.
15.1.17. Target platform¶
As a default option, FORCESPRO generates code for simulation on the host
platform. To obtain code for deployment on a target embedded platform, set the
field platform
to the appropriate value. The platforms currently supported by
FORCESPRO are given in Table 15.6. In order to acquire
licenses to use a specific platform
, licenses can be requested on the portal
by selecting the platform naming stated in the Portal Selection
.

Description 
Portal Selection 


For the architecture of the host platform. 


For x86_64 based 64bit platforms (detected OS). 


For x86 based 32bit platforms (detected OS). 


For Windows x86_64 based 64bit platforms (supports Microsoft/Intel compiler). 


For Windows x86 based 32bit platforms (supports Microsoft/Intel compiler). 


For Windows x86_64 based 64bit platforms (supports MinGW compiler). 


For Windows x86 based 32bit platforms (supports MinGW compiler). 


For Mac x86_64 based 64bit platforms (supports GCC/Clang compiler). 


For Linux x86_64 based 64bit platforms (supports GCC compiler). 


For Linux x86 based 32bit platforms (supports GCC compiler). 


For Linux x86_64 based 64bit platforms on Docker (supports GCC compiler). 


For Linux x86 based 32bit platforms on Docker (supports GCC compiler). 


For ARM Cortex 32bit processors (Gnueabih machine type). 


For ARM Cortex 64bit processors (Aarch machine type). 


For ARM Cortex 32bit processors using the Integrity toolchain. 


For ARM Cortex 64bit processors using the Integrity toolchain. 


For ARM Cortex M3 32bit processors. 


For the ARM Cortex M4 32bit processors without a floatingpoint unit. 


For the ARM Cortex M4 32bit processors with a floatingpoint unit. 


For the ARM Cortex A7 32bit processors (Gnueabih machine type). 


For the ARM Cortex A8 32bit processors (Gnueabih machine type). 


For the ARM Cortex A9 32bit processors (Gnueabih machine type). 


For the ARM Cortex A15 32bit processors (Gnueabih machine type). 


For the ARM Cortex A53 64bit processors (Gnueabih machine type). 


For the ARM Cortex A72 64bit processors (Gnueabih machine type). 


For the ARM Cortex A15 32bit processors (Gnueabih machine type). 


For the NVIDIA Cortex A57 64bit processors (Aarch machine type). 


For the ARM Cortex A57 64bit processors (Aarch machine type). 


For the ARM Cortex A72 64bit processors (Aarch machine type). 


For 32bit PowerPC based platforms (supports GCC compiler). 


For 64bit PowerPC based platforms (supports GCC compiler). 


For Windows x86 based 32bit platforms (supports MinGW compiler). 


For Windows x86_64 based 64bit platforms (supports MinGW compiler). 


For the dSPACE MicroAutoBox II realtime system (supports Microtec compiler). 


For the dSPACE MicroAutoBox III realtime system (supports Gcc compiler). 


For the dSPACE MicroAutoBox II realtime system (supports Microtec compiler). 


For the dSPACE MicroAutoBox III realtime system (supports Gcc compiler). 


For the dSPACE AutoBox realtime system (supports Gcc compiler). 


For Speedgoat 32bit realtime platforms (supports Microsoft compiler and mainly MATLAB Releases 2018b up to R2020a). 


For Speedgoat 64bit realtime platforms (supports Microsoft compiler and mainly MATLAB Releases 2018b up to R2020a). 


For Speedgoat 64bit realtime platforms (supports MATLAB Releases 2020b onwards). 


For Speedgoat Mobile 32bit realtime platforms (supports Microsoft compiler and Matlab Releases 2018a and earlier). 


For National Instruments compactRIO Linux RTOS platforms (supports NILRT Gcc compiler). 


For Bachmann PLC platforms (supports VxWorks compiler). 

Note
If a solver for another platform is requested, FORCESPRO will still
provide the simulation interfaces for the 'Generic'
host platform to
enable users to run simulations.
15.1.17.1. Cross compilation¶
To generate code for other operating systems different from the host platform, set the appropriate flag from the following list to 1
:
codeoptions.win
codeoptions.mac
codeoptions.gnu
Note that this will only affect the target platform. Interfaces for the host platform will be automatically built.
15.1.17.2. Mac compilation¶
When compiling for mac platforms it’s possible to select the compiler to be used for the web compilation. Select from the available
values gcc
(default) and clang
with the following codeoption:
codeoptions.maccompiler = 'gcc'; % or 'clang'
codeoptions.maccompiler = 'gcc' # or 'clang'
15.1.17.3. SIMD instructions¶
On x86based host platforms, one can enable the sse
field to accelerate the execution of the solver
codeoptions.sse = 1;
codeoptions.sse = 1
On x86based host platforms, one can also add the avx
field to significantly accelerate the
compilation and execution of the convex solver, from version 1.9.0
,
codeoptions.avx = 1;
codeoptions.avx = 1
Note
Currently when options avx
and blckMatrices
are enabled simultaneously, blckMatrices
is automatically disabled.
Note
When sparse parameters are present in the model, the options avx
and neon
are automatically set to zero.
Depending on the host platform, avx
may be automatically enabled. If the machine on which the solver is to be run does not support AVX and the message “Illegal Instruction” is returned at runtime, one can
explicitly disable avx
by setting:
codeoptions.avx = 1;
codeoptions.avx = 1
If the host platform supports AVX, but the user prefers not to have AVX intrinsics in the generated code, one can also keep the default option value of the solver:
codeoptions.avx = 0;
codeoptions.avx = 0
On ‘NVIDIACortexA57’, ‘AARCHCortexA57’ and ‘AARCHCortexA72’ target platforms, one can also enable the field neon
in order to accelerate the execution of the convex solver. From version 1.9.0
, the typical behaviour is that the host platform gets vectorized code based on AVX intrinsics when avx = 1
, and the target platform gets AVX vectorized code if it supports it when avx = 1
and NEON vectorized code if it is one of the above Cortex platforms and neon = 1
.
For single precision, the options are
codeoptions.floattype = 'float';
codeoptions.neon = 1;
codeoptions.floattype = 'float'
codeoptions.neon = 1
For double precision, the options are
codeoptions.floattype = 'double';
codeoptions.neon = 2;
codeoptions.floattype = 'double'
codeoptions.neon = 2
In case one wants to disable NEON intrinsics in the generated target code, the default value of the neon
option is
codeoptions.neon = 0;
codeoptions.neon = 0
If NEON vectorization is being used and there is a mismatch between float precision and the value of the neon option, the solver is automatically generated with the following options:
codeoptions.floattype = 'double';
codeoptions.neon = 2;
codeoptions.floattype = 'double'
codeoptions.neon = 2
and a warning message is raised by the MATLAB client.
Note
From version 1.9.0
, ARMv8A NEON instructions are generated. Hence, target platforms based on ARMv7 and previous versions are currently not supported.
15.1.18. Tips for solving QPs in single precision¶
Solving QPs in single precision can be rather challenging, i.e. nonconverging solves are likely to occur due to the lack of accuracy. In order to mitigate this undesirable behaviour, several options can be tuned to make convergence more robust. They are shown and commented in the code snippet below.
codeoptions.floattype = 'float';
codeoptions.regularize.epsilon = 1e5; % Tolerance on pivot in factorization
codeoptions.regularize.delta = 5e3; % Onthefly regularization coefficient in factorization
codeoptions.regularize.epsilon2 = 1e5; % Tolerance on pivot in factorization
codeoptions.regularize.delta2 = 5e3; % Onthefly regularization coefficient in factorization
codeoptions.accuracy.ineq = 1e4; % infinity norm of residual for inequalities
codeoptions.accuracy.eq = 1e4; % infinity norm of residual for equalities
codeoptions.accuracy.mu = 1e6; % absolute duality gap
codeoptions.accuracy.rdgap = 1e4; % relative duality gap := (pobjdobj)/pobj
codeoptions.init = 1;
codeoptions.floattype = 'float'
codeoptions.regularize.epsilon = 1e5 # Tolerance on pivot in factorization
codeoptions.regularize.delta = 5e3 # Onthefly regularization coefficient in factorization
codeoptions.regularize.epsilon2 = 1e5 # Tolerance on pivot in factorization
codeoptions.regularize.delta2 = 5e3 # Onthefly regularization coefficient in factorization
codeoptions.accuracy.ineq = 1e4 # infinity norm of residual for inequalities
codeoptions.accuracy.eq = 1e4 # infinity norm of residual for equalities
codeoptions.accuracy.mu = 1e6 # absolute duality gap
codeoptions.accuracy.rdgap = 1e4 # relative duality gap := (pobjdobj)/pobj
codeoptions.init = 1;
In general, the rationale behind this tuning is to make the tolerances looser and increase the regularization in the linear algebra. Note that these tips are only applicable to QP solvers. Solving NLPs in single precision is even more challenging and we currently do not offer a set of options to robustify convergence on this type of problems.
15.1.19. MISRA 2012 compliance¶
If your license allows it, add the following field to generate C code that is compliant with the MISRA 2012 rules:
codeoptions.misra2012_check = 1;
codeoptions.misra2012_check = 1
This option makes the generated solver code MISRA compliant. After compilation, the client also downloads a folder
whose name terminates with _misra2012_analysis
. The folder contains one summary of all MISRA violations for the solver source
and header files. Note that the option only produces MISRA compliant code when used with algorithms PDIP
and PDIP_NLP
.
15.1.20. Optimizing code size¶
The size of the solver libraries generated with code option PDIP_NLP
can be reduced by means of the option nlp.compact_code
. By setting
codeoptions.nlp.compact_code = 1;
codeoptions.nlp.compact_code = 1
the user enables the FORCESPRO server to generate smaller code, which results in shorter compilation time and slightly better solve time in some cases. This feature is especially effective on long horizon problems.
Note
The compact_code
option is currently only supported when using the linear systems solver
codeoptions.nlp.linear_solver = 'normal_eqs'
(which is the default choice).
The size of sparse linear algebra routines in the generated code can be reduced by changing the option compactSparse
from 0
to 1
:
codeoptions.compactSparse = 1;
codeoptions.compactSparse = 1
15.1.21. Optimizing Linear Algebra Operations¶
Some linear algebra routines in the generated code have available optimizations that can be enabled by changing the options optimize_<optimization>
from 0
to 1
.
These optimizations change the code in order to make better use of some embedded architectures in which hardware is more limited compared to host PC architectures.
Therefore, these optimizations show better results in embedded platforms such as ARM targets rather than during simulations on host PCs.
The available optimizations are:
Cholesky Division: This option performs the divisions included in the Cholesky factorization more efficiently to reduce its computation time.
Registers: This option attempts to use the architecture’s registers in order to reduce memory operations which can take significant time.
Use Locals: These options (which are separated into
simple
/heavy
/all
in ascending complexity) make better use of data locality in order to reduce memory jumpsOperations Rearrange: This option rearranges operations in order to make more efficient use of data and reduce memory jumps
Loop Unrolling: This option unrolls some of the included loops in order to remove their overhead.
Enable Offset: This option allows the rest of the optimizations to take place in cases where the matrix contains offsets.
codeoptions.optimize_choleskydivision = 1;
codeoptions.optimize_registers = 1;
codeoptions.optimize_uselocalsall = 1;
codeoptions.optimize_uselocalsheavy = 1; % overriden if uselocalsall is enabled
codeoptions.optimize_uselocalssimple = 1; % overriden if uselocalsheavy is enabled
codeoptions.optimize_operationsrearrange = 1;
codeoptions.optimize_loopunrolling = 1;
codeoptions.optimize_enableoffset = 1;
codeoptions.optimize_choleskydivision = 1
codeoptions.optimize_registers = 1
codeoptions.optimize_uselocalsall = 1
codeoptions.optimize_uselocalsheavy = 1 # overriden if uselocalsall is enabled
codeoptions.optimize_uselocalssimple = 1 # overriden if uselocalsheavy is enabled
codeoptions.optimize_operationsrearrange = 1
codeoptions.optimize_loopunrolling = 1
codeoptions.optimize_enableoffset = 1
15.1.22. Dump problem formulation and data¶
The MATLAB client of FORCESPRO provides a builtin tool to dump the problem formulation to reproduce the exact same solver for future reference. This tool is explained in detail in Section 17 and can be turned on by using the setting:
codeoptions.dump_formulation = 1;
Furthermore, you can dump problem structs containing the runtime parameters from C as described in Section 17. This tool is enabled for the host or/and the target platform by setting:
codeoptions.serializeCParamsHost = 1;
codeoptions.serializeCParamsTarget = 1;
codeoptions.serializeCParamsHost = 1
codeoptions.serializeCParamsTarget = 1
15.2. Highlevel interface options¶
The FORCESPRO NLP solver of the highlevel interface implements a nonlinear barrier interiorpoint method. We will now discuss how to change several parameters in the solver.
15.2.1. Integrators¶
When providing the continuous dynamics the user must select a particular
integrator by setting nlp.integrator.type
as outlined in Table 15.7.

Type 
Order 


Explicit Euler Method 
1 

Explicit RungeKutta 
2 

Explicit RungeKutta 
3 

Explicit RungeKutta 
4 

Implicit Euler Method 
1 

Implicit RungeKutta 
2 

Implicit RungeKutta 
4 
The user must also provide the discretization interval (in seconds) and the number of intermediate shooting nodes per interval. For instance:
codeoptions.nlp.integrator.type = 'ERK2';
codeoptions.nlp.integrator.Ts = 0.01;
codeoptions.nlp.integrator.nodes = 10;
codeoptions.nlp.integrator.type = 'ERK2'
codeoptions.nlp.integrator.Ts = 0.01
codeoptions.nlp.integrator.nodes = 10
Tip
Usually an explicit integrator such as RK4 should suffice for most applications. If you have stiff systems, or suspect inaccurate integration to be the cause of convergence failure of the NLP solver, consider using implicit integrators from the table above.
Note
Note that the implicit integrators BackwardEuler
, IRK2
and IRK4
currently
rely on the CasADi AD tool to work.
15.2.1.1. Expert options for implicit integrators¶
The implicit integrators BackwardEuler
, IRK2
and IRK4
do not just evaluate
the differential equation, but actually solve a nonlinear equation to obtain the state
trajectory. This is done by means of Newton iterations, with default values of 10
iterations for BackwardEuler
and 5
iterations for IRK2
and IRK4
. These
default values can be overwritten by using the following option:
codeoptions.nlp.integrator.newtonIter = 3;
codeoptions.nlp.integrator.newtonIter = 3
In order to reduce computational effort, the Jacobian of the nonlinear equation is only computed once by default. If your differential equations are highly nonlinear, it may be worth the effort to recompute it at every Newton iteration. This is achieved by means of the following option:
codeoptions.nlp.integrator.reuseNewtonJacobian = false;
codeoptions.nlp.integrator.reuseNewtonJacobian = false;
Note
In Python, currently the expert options for the implicit integration schemes only work for the codegenerated IRK2
integrator
and not the legacy integrators (see section Codegenerated integrators).
15.2.1.2. Codegenerated integrators¶
From FORCESPRO 4.1.0
, integrators generated on the server are available when using explicit integrators ForwardEuler
, RK2
, RK3
and RK4
, and
the field continuous_dynamics is set in the model structure. From FORCESPRO 4.4.0
the implicit integration scheme IRK2
was added to the list of supported codegenerated integration schemes.
These integrators result in much smaller code size than previously. They also often result in faster run times on embedded targets.
Two different methods are used to compute sensitivities associated to these integrators:
chainrule
, which is the default option, can also be triggered by setting
codeoptions.nlp.integrator.differentiation_method = 'chainrule';
codeoptions.nlp.integrator.differentiation_method = 'chainrule'
vde
, which can be triggered by settings the following option:
codeoptions.nlp.integrator.differentiation_method = 'vde';
codeoptions.nlp.integrator.differentiation_method = 'vde'
When using the vde
option, the following options also need to be set
codeoptions.nlp.sensitivity.Ts = codeoptions.nlp.integrator.Ts;
codeoptions.nlp.sensitivity.nodes = codeoptions.nlp.integrator.nodes / 2;
% When using 'ERK2' or 'ERK4' for the sensitivity computation, the number of nodes for the sensitivity
% needs to be twice the number of nodes for the integrator
codeoptions.nlp.sensitivity.type = 'ERK4'; % Can also be 'ForwardEuler', 'RK2' depending on the application
codeoptions.nlp.sensitivity.Ts = codeoptions.nlp.integrator.Ts
codeoptions.nlp.sensitivity.nodes = codeoptions.nlp.integrator.nodes / 2
# When using 'ERK2' or 'ERK4' for the sensitivity computation, the number of nodes for the sensitivity
# needs to be twice the number of nodes for the integrator
codeoptions.nlp.sensitivity.type = 'ERK4' # Can also be 'ForwardEuler', 'RK2' depending on the application
The vde
option is likely to change the numerical behaviour of the solver but can help for reducing the solve time in some cases, typically by having a looser integration on sensitivity.
Note
The vde
option currently is still in an experimental state and we are working to fully robustify it.
You may give it a try, but be prepared for unexpected behaviour. Also, the RK3
integration method
is currently not supported with the vde
option.
15.2.1.3. Linear subsystem exploitation¶
Often nonlinear optimal control problems contain linear subsystems, meaning part of the differential equation describing the dynamics of the system is linear while another part is nonlinear. By this we mean that the state \(x\) of the system can be split into two parts \(x = (x_1, x_2)\) such that the differential equation \(\dot{x} = c(x,u)\) (\(u\) denoting the control input) governing the dynamics of the system can be written as
Here \(A_1\) and \(B_1\) denote constant matrices and \(c_2\) denotes a nonlinear function. Since FORCESPRO 4.4.0
it is possible to exploit such subsystems for performance by performing parts of the numerical
integration of the system offline. Currently this is supported only for the codegenerated ERK4
integration scheme (see Codegenerated integrators).
FORCESPRO can automatically detect a linear subsystem if it exists. One can activate the detection of linear subsystems by enabling the codeoptions.nlp.integrator.attempt_subsystem_exploitation
option as follows:
codeoptions.nlp.integrator.attempt_subsystem_exploitation = 1;
codeoptions.nlp.integrator.attempt_subsystem_exploitation = 1
Optionally, in combination with setting this option, one can also specify the state indices of the linear subsystem manually. These indices are specified as a numpy array of integers in Python and a vector of indices in Matlab via the
field model.linInIdx
. For example, in a case when \(x\in \mathbb{R}^2\), \(u\in \mathbb{R}\) and the righthandside \(c\) of the ODE describing the dynamics of the system is given by
one would have to specify
model.linInIdx = [1, 2];
model.linInIdx = np.array([0, 1], dtype=np.int)
Note the 1based indexing in Matlab and the 0based indexing in Python. For further details on how to exploit linear subsystems using FORCESPRO, see Controlling a crane using a FORCESPRO NLP solver.
Note
For large systems (more than about 16 states) there might be a considerable overhead in determining the indices of the linear subsystem automatically. In case you encounter such an overhead, you can avoid it by manually specifying
model.linInIdx
as shown above.
15.2.2. Accuracy requirements¶
One can modify the termination criteria by altering the KKT tolerance with respect to stationarity, equality constraints, inequality constraints and complementarity conditions, respectively, using the following fields:
% default tolerances
codeoptions.nlp.TolStat = 1e5; % inf norm tol. on stationarity
codeoptions.nlp.TolEq = 1e6; % tol. on equality constraints
codeoptions.nlp.TolIneq = 1e6; % tol. on inequality constraints
codeoptions.nlp.TolComp = 1e6; % tol. on complementarity
# default tolerances
codeoptions.nlp.TolStat = 1e5 # inf norm tol. on stationarity
codeoptions.nlp.TolEq = 1e6 # tol. on equality constraints
codeoptions.nlp.TolIneq = 1e6 # tol. on inequality constraints
codeoptions.nlp.TolComp = 1e6 # tol. on complementarity
All tolerances are computed using the infinitiy norm \(\lVert \cdot \rVert_\infty\).
15.2.3. Barrier strategy¶
The strategy for updating the barrier parameter is set using the field:
codeoptions.nlp.BarrStrat = 'loqo';
codeoptions.nlp.BarrStrat = 'loqo'
It can be set to 'loqo'
(default) or to 'monotone'
. The default settings
often leads to faster convergence, while 'monotone'
may help convergence for
difficult problems.
15.2.4. Hessian approximation¶
The way the Hessian of the Lagrangian function is computed can be set using the field:
codeoptions.nlp.hessian_approximation = 'bfgs';
codeoptions.nlp.hessian_approximation = 'bfgs'
FORCESPRO currently supports BFGS updates ('bfgs'
) (default) and
GaussNewton approximation ('gaussnewton'
). Exact Hessians will be
supported in a future version. Read the subsequent sections for the
corresponding Hessian approximation method of your choice.
15.2.4.1. BFGS options¶
When the Hessian is approximated using BFGS updates, the initialization of the estimates can play a critical role in the convergence of the method. The default value is the identity matrix, but the user can modify it using e.g.:
codeoptions.nlp.bfgs_init = diag([0.1, 10, 4]);
codeoptions.nlp.bfgs_init = np.diag(np.array([0.1, 10, 4]))
Note that BFGS updates are carried out individually per stage in the FORCESPRO NLP solver, so the size of this matrix is the size of the stage variable. Also note that this matrix must be positive definite. When the cost function is positive definite, it often helps to initialize BFGS with the Hessian of the cost function.
This matrix is also used to restart the BFGS estimates whenever the BFGS updates are skipped several times in a row. The maximum number of updates skipped before the approximation is reinitialized is set using:
codeoptions.nlp.max_update_skip = 2;
codeoptions.nlp.max_update_skip = 2
The default value for max_update_skip
is 2.
In order to set the BFGS initialization through the bfgs_init
codeoption one must first come up with a guess for a good BFGS initialization.
One way to do so is to first run the solver without any userdefined BFGS initialization (i.e. not setting codeoptions.nlp.bfgs_init
)
and using the BFGS matrix reached upon convergence as an inizialization. One can export the BFGS matrix by setting
codeoptions.exportBFGS = 1;
codeoptions.exportBFGS = 1
Istead of specifying the BFGS initialization at codegen one can also specify it at runtime. In order to do this one should set
codeoptions.nlp.parametricBFGSinit = 1;
codeoptions.nlp.parametricBFGSinit = 1
before generating the FORCESPRO solver. Having done this, the generated solver will expect an input problem.BFGSinitLower<stage number>
for every stage.
This is a vector which specifies the BFGS hessian initialization in LOWER TRIANGULAR ROW MAJOR format. Thus, in order to specify e.g. the matrix
for constants \(a_1,a_2,a_3 > 0\) as the BFGS inizialization at stage 6 out of 50 stages in total, one would specify
problem.BFGSinitLower06 = [a_1, 0, a_2, 0, 0, a_3];
problem["BFGSinitLower06"] = np.array([a_1, 0, a_2, 0, 0, a_3])
15.2.4.2. GaussNewton options¶
For problems that have a least squares objective, i.e. the cost function can be expressed by a vectorvalued function \(r_k : \mathbb{R}^n \rightarrow \mathbb{R}^m\) which implicitly defines the objective function as:
the GaussNewton approximation of the Hessian is given by:
and can lead to faster convergence and a more reliable method. When this option
is selected, the functions \(r_k\) have to be provided by the user in the
field LSobjective
. For example if \(r(z)=\sqrt{100} z_1^2 + \sqrt{6}
z_2^2\), i.e. \(f(z) = 50 z_1^2 + 3 z_2^2\), then the following code defines
the leastsquares objective (note that \(r\) is a vectorvalued function):
model.objective = @(z) 0.1* z(1)^2 + 0.01*z(2)^2;
model.LSobjective = @(z) [sqrt(0.2)*z(1); sqrt(0.02)*z(2)];
# not yet implemented
Important
The field LSobjective
will have precedence over objective
, which need
not be defined in this case.
When providing your own function evaluations in C, you must populate the Hessian argument with a positive definite Hessian.
15.2.5. Line search settings¶
The line search first computes the maximum step that can be taken while maintaining the iterates inside the feasible region (with respect to the inequality constraints). The maximum distance is then scaled back using the following setting:
% default fractiontoboundary scaling
codeoptions.nlp.ftbr_scaling = 0.9900;
# default fractiontoboundary scaling
codeoptions.nlp.ftbr_scaling = 0.9900;
15.2.6. Regularization¶
To avoid illconditioned saddle point systems, FORCESPRO employs two different types of regularization, static and dynamic regularization.
15.2.6.1. Static regularization¶
Static regularization of the augmented Hessian by \(\delta_w I\), and of the multipliers corresponding to the equality constraints by \(\delta_c I\) helps avoid problems with rank deficiency. The constants \(\delta_w\) and \(\delta_c\) vary at each iteration according to the following heuristic rule:
where \(\mu\) is the barrier parameter and \(i\) is the number of iterations.
This rule has been chosen to accommodate two goals: First, make the regularization dependent on the progress of the algorithm  the closer we are to the optimum, the smaller the regularization should be in order not to affect the search directions generated close to the solution, promoting superlinear convergence properties. Second, the amount of regularization employed should decrease with the number of iterations to a certain minimum level, at a certain sublinear rate, in order to prevent stalling due to too large regularization. FORCESPRO NLP does not employ an inertiacorrecting linear system solver, and so relies heavily on the parameters of this regularization to be chosen carefully.
You can change these parameters by using the following settings:
% default static regularization parameters
codeoptions.nlp.reg_eta_dw = 1e4;
codeoptions.nlp.reg_beta_dw = 0.8;
codeoptions.nlp.reg_min_dw = 1e9;
codeoptions.nlp.reg_gamma_dw = 1.0/3.0;
codeoptions.nlp.reg_eta_dc = 1e4;
codeoptions.nlp.reg_beta_dc = 0.8;
codeoptions.nlp.reg_min_dc = 1e9;
codeoptions.nlp.reg_gamma_dc = 1.0/3.0;
# default static regularization parameters
codeoptions.nlp.reg_eta_dw = 1e4
codeoptions.nlp.reg_beta_dw = 0.8
codeoptions.nlp.reg_min_dw = 1e9
codeoptions.nlp.reg_gamma_dw = 1.0/3.0
codeoptions.nlp.reg_eta_dc = 1e4
codeoptions.nlp.reg_beta_dc = 0.8
codeoptions.nlp.reg_min_dc = 1e9
codeoptions.nlp.reg_gamma_dc = 1.0/3.0
Note that by choosing \(\delta_w=0\) and \(\delta_c=0\), you can turn off the progress and iteration dependent regularization, and rely on a completely static regularization by \(\delta_{w,\min}\) and \(\delta_{c,\min}\), respectively.
15.2.6.2. Dynamic regularization¶
Dynamic regularization regularizes the matrix onthefly to avoid instabilities due to numerical errors. During the factorization of the saddle point matrix, whenever it encounters a pivot smaller than \(\epsilon\), it is replaced by \(\delta\). There are two parameter pairs: \((\epsilon,\delta)\) affects the augmented Hessian and \((\epsilon_2,\delta_2)\) affects the search direction computation. You can set these parameters by:
% default dynamic regularization parameters
codeoptions.regularize.epsilon = 1e12; % (for Hessian approx.)
codeoptions.regularize.delta = 4e6; % (for Hessian approx.)
codeoptions.regularize.epsilon2 = 1e14; % (for Normal eqs.)
codeoptions.regularize.delta2 = 1e14; % (for Normal eqs.)
# default dynamic regularization parameters
codeoptions.regularize.epsilon = 1e12 # (for Hessian approx.)
codeoptions.regularize.delta = 4e6 # (for Hessian approx.)
codeoptions.regularize.epsilon2 = 1e14 # (for Normal eqs.)
codeoptions.regularize.delta2 = 1e14 # (for Normal eqs.)
15.2.7. Linear system solver¶
The interiorpoint method solves a linear system to find a search direction at every iteration. FORCESPRO NLP offers the following three linear solvers:
'normal_eqs'
(default): Solving the KKT system in normal equations form.'symm_indefinite_fast'
: Solving the KKT system in augmented / symmetric indefinite form, using regularization and positive definite Cholesky factorizations only.'symm_indefinite'
: Solving the KKT system in augmented / symmetric indefinite form, using blockindefinite factorizations.
The linear system solver can be selected by setting the following field:
codeoptions.nlp.linear_solver = 'symm_indefinite';
codeoptions.nlp.linear_solver = 'symm_indefinite'
It is recommended to try different linear solvers when experiencing convergence
problems. The most stable method is 'symm_indefinite'
, while the fastest
solver is 'symm_indefinite_fast'
.
Note
Independent of the linear system solver choice, the generated code is always libraryfree and statically allocated, i.e. it can be embedded anywhere.
The 'normal_eqs'
solver is the standard FORCESPRO linear system solver based on
a full reduction of the KKT system (the socalled normal equations form). It
works well for standard problems, especially convex problems or nonlinear
problems where the BFGS or GaussNewton approximations of the Hessian are
numerically sufficiently well conditioned.
The 'symm_indefinite'
solver is the most robust solver, but still
highspeed. It is based on blockwise factorization of the symmetric indefinite
form of the KKT system (the socalled augmented form). Each block is handled by
symmetric indefinite LDL factorization, with (modified) onthefly
BunchKaufmann permutations leading to boundedness of lower triangular factors
for highest numerical stability. This is our most robust linear system solver,
with only a modest performance penalty (about 30% compared to
'symm_indefinite_fast'
).
The 'symm_indefinite_fast'
solver is robust, but even faster. It is based on
blockwise factorization of the symmetric indefinite KKT matrix, where each
block is handled by a Cholesky factorization. It uses regularization to
increase numerical stability. Currently only used for recedinghorizon/MPClike
problems where dimensions of all stages are equal (minus the first and last
stage, those are handled separately). It is more robust and faster than the
normal equations form. This solver is likely to become the default option in
the future.
15.2.8. Automatic differentiation tool¶
If external functions and derivatives are not provided directly as C code by
the user, FORCESPRO makes use of an automatic differentiation (AD) tool to
generate efficient C code for all the functions (and their derivatives) inside
the problem formulation. Currently, three different AD tools are supported that
can be chosen by means of the setting nlp.ad_tool
as summarized in Table 15.8.

Tool 
URL 


CasADi v3.5.1 


CasADi v2.4.2 


CasADi v3.5.1 


CasADi v3.5.5 


MathWorks Symbolic Math Toolbox 
Note that MathWorks Symbolic Math Toolbox requires an additional license, which is why the default option is set to
codeoptions.nlp.ad_tool = 'casadi';
codeoptions.nlp.ad_tool = 'casadi'
Also note that the use of implicit integrators BackwardEuler
, IRK2
and
IRK4
(see Section 15.2.1) currently still rely on using the
CasADi AD tool.
15.2.9. Reuse of callback code¶
When defining your NLP problem formulation using an AD tool, you may specify objective functions, dynamic equations and constraints separately on each stage. In order to reduce the size of the generated callback code, FORCESPRO will perform a check whether all these callbacks are identical on any two or more stages and if so, only generates the callback code for those stages once. However, checking for exact identity can be tricky and may sometimes lead to false results. By default, FORCESPRO performs a less strict check for identity resulting in less duplicated callback code. If you observe that two stages are wrongly identified as identical, you can enable a more strict check by using the following codeoption:
codeoptions.nlp.strictCheckDistinctStages = 1;
# not yet supported
Note that using this option may be overly conservative and lead to duplicated callback code for different stages that are actually identical.
15.2.10. Safety checks¶
By default, the output of the function evaluations is checked for the presence
of NaN
s or INF
s in order to diagnose potential initialization
problems. In order to speed up the solver one can remove these checks by
setting:
codeoptions.nlp.checkFunctions = 0;
codeoptions.nlp.checkFunctions = 0
15.3. Convex branchandbound options¶
The settings of the FORCESPRO mixedinteger branchandbound convex solver are accessed
through the codeoptions.mip
struct. It is worthwhile to explore different
values for the settings in Table 15.9, as they might have a
severe impact on the performance of the branchandbound procedure.
Note
All the options described below are currently not available with the FORCESPRO nonlinear solver. For mixedinteger nonlinear programs and the available options, please have a look at paragraph Mixedinteger nonlinear solver.
Setting 
Values 
Default 


Any value \(\geq 0\) 


Any value \(\geq 0\) 











Any value \(> 0\) 


Any integer value \(\geq 0\) 

A description of each setting is given below:
mip.timeout
: Timeout in seconds, after which the search is stopped and the best solution found so far is returned.mip.mipgap
: Relative suboptimality after which the search shall be terminated. For example, a value of0.01
will search for a feasible solution that is at most 1%suboptimal. Set to zero if the optimal solution is required.mip.branchon
: Determines which variable to branch on after having solved the relaxed problem. Options are'mostAmbiguous'
(i.e. the variable closest to 0.5) or'leastAmbiguous'
(i.e. the variable closest to 0 or 1).mip.stageinorder
: Stageinorder heuristic: For the branching, determines whether to fix variables in order of the stage number, i.e. first all variables of stage \(i\) will be fixed before fixing any of the variables of stage \(i+1\). This is often helpful in multistage problems, where a timeout is expected to occur, and where it is important to fix the early stages first (for example MPC problems). Options are0
for OFF and1
for ON.mip.explore
: Determines the exploration strategy when selecting pending nodes. Options are'bestFirst'
, which chooses the node with the lowest lower bound from all pending nodes, or'depthFirst'
, which prioritizes nodes with the most number of fixed binaries first to quickly reach a node.mip.inttol
: Integer tolerance for identifying binary solutions of relaxed problems. A solution of a relaxed problem with variable values that are belowinttol
away from binary will be declared to be binary.mip.queuesize
: Maximum number of pending nodes that the branch and bound solver can store. If that number is exceeded during the search, the solver quits with anexitflag
value of2
and returns the best solution found so far.
15.4. Solve methods¶
As a default optimization method the primaldual interiorpoint method is used.
Several other methods are available. To change the solve method set the
solvemethod
field as outlined in Table 15.10.

Method 
Description 


Nonlinear PrimalDual InteriorPoint method 
The Nonlinear PrimalDual InteriorPoint method is the most stable and robust method for most nonlinear problems. 

Sequential Quadratic Programming method 
The Sequential Quadratic Programming method may be more efficient on mildly nonlinear problems. 

PrimalDual InteriorPoint method 
The PrimalDual InteriorPoint method is the most stable and robust method for most convex problems. 

Alternating Direction Methods of Multipliers 
For some problems, ADMM may be faster. The method variant and several algorithm parameters can be tuned in order to improve performance. 

Dual Fast Gradient method 
For some problems with simple constraints, our implementation of the dual fast gradient method can be the fastest option. No parameters need to be tuned in this method. 

Fast Gradient method 
For problems with no equality constraints (only one stage) and simple constraints, the primal fast gradient method can give medium accuracy solutions extremely quickly. The method has several tuning parameters that can significantly affect the performance. 
15.4.1. PrimalDual InteriorPoint Method¶
The PrimalDual InteriorPoint method is the default optimization method for either nonlinear/nonconvex or convex problems. It is a stable and robust method for most of the problems.
15.4.1.1. Solver Initialization¶
The performance of the solver can be influenced by the way the variables are
initialized. The default method (cold start) should work in most cases extremely
reliably, so there should be no need in general to try other methods, unless you
are experiencing problems with the default initialization scheme. To change the
method of initialization in FORCESPRO set the field init
to one of the
values in Table 15.11.

Method 
Initialization method 


Cold start 
Set all primal variables to \(0\), and all dual variables to the square root of the initial complementarity gap \(\mu_0: z_i=0, s_i=\sqrt{\mu_0}, \lambda_i=\sqrt{\mu_0}\). The default value is \(\mu_0=10^6\). 

Centered start 
Set all primal variables to zero, the slacks to the RHS of the corresponding inequality, and the Lagrange multipliers associated with the inequalities such that the pairwise product between slacks and multipliers is equal to the parameter \(\mu_0: z_i=0, s_i=b_{\mathrm{ineq}}\) and \(s_i \lambda_i = \mu_0\). 

Primal warm start 
Set all primal variables as provided by the user. Calculate the residuals and set the slacks to the residuals if they are sufficiently positive (larger than \(10^{4}\)), or to one otherwise. Compute the associated Lagrange multipliers such that \(s_i \lambda_i = \mu_0\). 
15.4.1.2. Initial Complementary Slackness¶
The default value for \(\mu_0\) is \(10^6\). To use a different value, use:
codeoptions.mu0 = 10;
codeoptions.mu0 = 10;
15.4.1.3. Accuracy Requirements¶
The accuracy for which FORCESPRO returns the OPTIMAL flag can be set as follows:
codeoptions.accuracy.ineq = 1e6; % infinity norm of residual for inequalities
codeoptions.accuracy.eq = 1e6; % infinity norm of residual for equalities
codeoptions.accuracy.mu = 1e6; % absolute duality gap
codeoptions.accuracy.rdgap = 1e4; % relative duality gap := (pobjdobj)/pobj
codeoptions.accuracy.ineq = 1e6 # infinity norm of residual for inequalities
codeoptions.accuracy.eq = 1e6 # infinity norm of residual for equalities
codeoptions.accuracy.mu = 1e6 # absolute duality gap
codeoptions.accuracy.rdgap = 1e4 # relative duality gap := (pobjdobj)/pobj
15.4.1.4. Line Search Settings¶
If FORCESPRO experiences convergence difficulties, you can try selecting
different line search parameters. The first two parameters of
codeoptions.linesearch
, factor_aff
and factor_cc
are the backtracking factors
for the line search (if the current step length is infeasible, then it is
reduced by multiplication with these factors) for the affine and combined search
direction, respectively.
codeoptions.linesearch.factor_aff = 0.9;
codeoptions.linesearch.factor_cc = 0.95;
codeoptions.linesearch.factor_aff = 0.9
codeoptions.linesearch.factor_cc = 0.95
The remaining two parameters of the field linesearch
determine the minimum
(minstep
) and maximum step size (maxstep
). Choosing minstep
too high will cause
the generated solver to quit with an exitcode saying that the line search has
failed, i.e. no progress could be made along the computed search direction.
Choosing maxstep
too close to 1 is likely to cause numerical issues, but
choosing it too conservatively (too low) is likely to increase the number of
iterations needed to solve a problem.
codeoptions.linesearch.minstep = 1e8;
codeoptions.linesearch.maxstep = 0.995;
codeoptions.linesearch.minstep = 1e8
codeoptions.linesearch.maxstep = 0.995
15.4.1.5. Regularization¶
During factorization of supposedly positive definite matrices, FORCESPRO applies a regularization to the \(i\)th pivot element if it is smaller than \(\epsilon\). In this case, it is set to \(\delta\), which is the lower bound on the pivot element that FORCESPRO allows to occur.
codeoptions.regularize.epsilon = 1e13; % if pivot element < epsilon …
codeoptions.regularize.delta = 1e8; % then set it to delta
codeoptions.regularize.epsilon = 1e13 # if pivot element < epsilon …
codeoptions.regularize.delta = 1e8 # then set it to delta
15.4.1.6. Multicore parallelization¶
FORCESPRO supports the computation on multiple cores, which is particularly useful for large problems and long horizons (the workload is split along the horizon to multiple cores). This is implemented by the use of OpenMP and can be switched on by using
codeoptions.parallel = 1;
codeoptions.parallel = 1
By default multicore computation is switched off.
When the parallel option is enabled with 1 (codeoptions.parallel = 1
), the maximum number
of threads to be used is set as the maximum number of threads available to OpenMP (max_number_of_threads
).
Additionally, a runtime parameter num_of_threads
is created to control
the number of threads in runtime. The allowed range of values for the runtime
parameter is [1, max_number_of_threads]
. Leaving the parameter unset or
setting a value outside the allowed range will lead in execution with the
maximum number of threads (max_number_of_threads
).
The maximum number of threads can also be set manually during code generation by setting:
% <max_number_of_threads> larger than 1
codeoptions.parallel = <max_number_of_threads>;
# <max_number_of_threads> larger than 1
codeoptions.parallel = <max_number_of_threads>
15.4.2. Alternating Directions Method of Multipliers¶
FORCESPRO implements several optimization methods based on the ADMM framework.
Different variants can handle different types of constraints and FORCESPRO will
automatically choose an ADMM variant that can handle the constraints in a given
problem. To manually choose a specific method in FORCESPRO, use the ADMMvariant
field of codeoptions
:
codeoptions.ADMMvariant = 1; % can be 1 or 2
codeoptions.ADMMvariant = 1 # can be 1 or 2
where variant 1
is as follows:
and variant 2
is as follows:
15.4.2.1. Accuracy requirements¶
The accuracy for which FORCESPRO returns the OPTIMAL flag can be set as follows:
codeoptions.accuracy.consensus = 1e3; % infinity norm of the consensus equality
codeoptions.accuracy.dres = 1e3; % infinity norm of the dual residual
codeoptions.accuracy.consensus = 1e3 # infinity norm of the consensus equality
codeoptions.accuracy.dres = 1e3 # infinity norm of the dual residual
Note that, in contrast to primaldual interiorpoint methods, the required number of ADMM iterations varies very significantly depending on the requested accuracy. ADMM typically requires few iterations to compute medium accuracy solutions, but many more iterations to achive the same accuracy as interiorpoint methods. For feedback applications, medium accuracy solutions are typically sufficient. Also note that the ADMM accuracy requirements have to be changed depending on the problem scaling.
15.4.2.2. Method parameters¶
ADMM uses a regularization parameter \(\rho\), which also acts as the step size in the gradient step. The convergence speed of ADMM is highly variable in the parameter \(\rho\). Its value should satisfy \(\rho > 0\). This parameter can be tuned using the following command:
codeoptions.ADMMrho = 1;
codeoptions.ADMMrho = 1
In some cases it may be possible to let FORCESPRO choose the value \(\rho\) automatically. To enable this feature set:
codeoptions.ADMMautorho = 1;
codeoptions.ADMMautorho = 1
Please note that this does not guarantee that the choice of \(\rho\) will be optimal.
ADMM can also include an ‘overrelaxation’ step that can improve the convergence speed. This step is typically useful for problems where ADMM exhibits very slow convergence and can be tuned using the parameter \(\alpha\). Its value should satisfy \(1 \leq \alpha \leq 2\). This step using the following command:
codeoptions.ADMMalpha = 1;
codeoptions.ADMMalpha = 1
15.4.2.3. Precomputations¶
For problems with timeinvariant data, FORCESPRO can compute full matrix
inverses at code generation time and then implement matrix solves online by
dense matrixvector multiplication. In some cases, especially when the
prediction horizon is long, it may be better to factorize the matrix and
implement matrix solves using forward and backward solves with the precomputed
factors. To manually switch on this option, use the ADMMfactorize
field of
codeoptions
.
When the data is timevarying, or when the prediction horizon is larger than 15 steps, FORCESPRO automatically switches to a factorizationbased method.
codeoptions.ADMMfactorize = 0;
codeoptions.ADMMfactorize = 0
15.4.3. Dual Fast Gradient Method¶
For some problems with simple constraints, our implementation of the dual fast gradient method can be the fastest option. No parameters need to be tuned in this method.
15.4.4. Primal Fast Gradient Method¶
For problems with no equality constraints (only one stage) and simple constraints, the primal fast gradient method can give medium accuracy solutions extremely quickly. The method has several tuning parameters that can significantly affect the performance.
15.4.4.1. Accuracy requirements¶
The accuracy for which FORCESPRO returns the OPTIMAL flag can be set as follows:
codeoptions.accuracy.gmap = 1e5; % infinity norm of the gradient map
codeoptions.accuracy.gmap = 1e5 # infinity norm of the gradient map
The gradient map is related to the difference with respect to the optimal objective value. Just like with other firstorder methods, the required number of FG iterations varies very significantly depending on the requested accuracy. Medium accuracy solutions can typically be computed very quickly, but many iterations are needed to achieve the same accuracy as with interiorpoint methods.
15.4.4.2. Method parameters¶
The user has to determine the step size in the fast gradient method. The convergence speed of FG is highly variable in this parameter, which should typically be set to be one over the maximum eigenvalue of the quadratic cost function. This parameter can be tuned using the following command:
codeoptions.FGstep = 1/1000;
codeoptions.FGstep = 1/1000
In some cases it may be possible to let FORCESPRO choose the step size automatically. To enable this feature set:
codeoptions.FGautostep = 1;
codeoptions.FGautostep = 1
15.4.4.3. Warm starting¶
The performance of the fast gradient method can be greatly influenced by the way the variables are initialized. Unlike with interiorpoint methods, fast gradient methods can be very efficiently warm started with a good guess for the optimal solution. To enable this feature set:
codeoptions.warmstart = 1;
codeoptions.warmstart = 1
When the user turns warm start on, a new parameter z_init_0
is automatically
added. The user should set it to be a good guess for the solution, which is
typically available when solving a sequence of problems.