17. Solver Options

FORCESPRO expects a set of solver options that help to customize the generated solver code. Section Default options describes how to obtain a set of default options, while further customizations are described in the remaining sections.

In the documentation below, we assume that you have created this struct and named it codeoptions.

Note

For the low-level interface in Python, the codeoptions struct has to be an element of the stages struct.

17.1. Default options

Default solver options can be loaded by giving a name to the solver with the following command

codeoptions = getOptions('solvername');

This function is very generic and does not tailor the default options to the algorithm used. It is therefore recommended to use the following function instead (which has been introduce with release 6.0.0):

codeoptions = ForcesGetDefaultOptions( 'solvername','algorithm','floattype' );

In addition to the solver name, this function expects two more input arguments:

  • a string specifying the algorithm, which can take any value listed in Table 17.14 or 'SQP_NLP_fast' or 'QP_FAST' to denote the SQP Fast (see Different SQP variants) or QP Fast algorithm (see The QP_FAST algorithm), respectively;

  • a string 'floattype', which currently can be set to either 'double' or 'float' (see Table 17.5).

17.2. General options

We will first discuss how to change several options that are valid for all the FORCESPRO interfaces.

17.2.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';

Alternatively, you can directly name the solver when generating the options struct by calling:

codeoptions = getOptions('solvername');

17.2.3. Maximum number of iterations

To set the maximum number of iterations of the generated solver, use:

codeoptions.maxit = 200;

The default maximum number of iterations for all solvers provided by FORCESPRO is set to 200.

17.2.4. Parametric number of iterations

The maximum number of iterations of the generated solver can be made parametric with the following codeoption:

codeoptions.parametric_iterations = 1;

This will generate the runtime parameter maxit which can be provided as input to the solver. The runtime parameter needs to be initialized with a value in the range (0, codeoptions.maxit] (i.e. the runtime iterations parameter needs to be set lower or equal to the hardcoded maximum number of iterations).

Note

The parametric number of iterations feature is currently not available for the SQP solvers.

17.2.5. Compiler optimization level

The compiler optimization level can be varied by changing the field optlevel from 0 to 3 (default):

codeoptions.optlevel = 0;

Important

It is recommended to set optlevel to 0 or 1 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.

Note

In very rare cases, when the solver has a large amount of parameters or the problem is relatively large, compiling with codeoptions.optlevel = 0; has produced a stack overflow in the past. If you encounter such a situation, please try to increase the value of codeoptions.optlevel and submit a bug report to support@embotech.com.

17.2.6. Measure Computation time

You can measure the time used for executing the generated code by using:

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;

17.2.7. Solver Timeout

17.2.7.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 iteration

  • timeout_estimate_coeff is a coefficient used to make the estimate more conservative or forgiving. Its default value is 1.20

17.2.7.2. Usage

To enable the solver timeout you can use the following codeoption:

% solver_timeout can take values 0-2
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 Mixed-integer 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’s -1.0 which results in timeout being disabled

  • For timeout_estimate_coeff it’s 1.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).

17.2.7.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;

17.2.7.4. Return Value

When solver timeout is enabled, two additional exitflags are available for the user:

Table 17.2 Timeout exitflags

Exitflag Name

Value

Description

TIMEOUT_<SOLVERNAME>

2

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

INVALID_TIMEOUT_<SOLVERNAME>

-12

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).

17.2.8. Early-terminate solver

You can terminate a solver during execution from C based on an external event by setting the option

codeoptions.solver_exit_external = 1;

Setting the option to 1 will enable the external termination by providing a runtime parameter solver_exit_external which can take 2 values:

Table 17.3 solver_exit_external parameter

Value of solver_exit_external

Description

0

solver runs normally

1

solver terminates at the end of the current iteration

When a solver is early-terminated, it always finishes the current iteration and returns its solution.

A solver that was early-terminated before convergence returns a dedicated exitflag:

Table 17.4 Early-termination exitflag

Exitflag Name

Value

Description

EXIT_EXTERNAL_<SOLVERNAME>

3

The solver was terminated because solver_exit_external = 1 was set externally

Note

To enable external and asynchronous update of the variable solver_exit_external, it is declared as volatile (ensuring visibility) and of type sig_atomic_t (ensuring atomicity). If you use this parameter within a multithreaded environment to terminate a solver on another thread, please note that this assumes a cache coherent architecture to work reliably as the volatile qualifier itself does not guarantee thread safety.

17.2.10. Datatypes

The type of variables can be changed by setting the field floattype as outlined in Table 17.5. This will effect all floating point variables used inside the solver and the callbacks generated by the AD tool.

Table 17.5 Data type options

floattype

Decimation

Width (bits)

Supported algorithms

'double' (default)

64 bit

Floating point

PDIP_NLP, SQP_NLP, PDIP, QP_FAST, ADMM, DFG, FG

'float'

32 bit

Floating point

SQP_NLP, PDIP, QP_FAST, ADMM, DFG, FG

'int'

32 bit

Fixed point

ADMM, DFG, FG

'short'

16 bit

Fixed point

ADMM, DFG, FG

Important

Unless running on a resource-constrained 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 1E-6 to 1E-3).

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 17.6 and section Single precision callbacks for details.

Table 17.6 Callback data type options

floattype

Decimation

Width (bits)

Supported algorithms

'double' (default)

64 bit

Floating point

PDIP_NLP, SQP_NLP

'float'

32 bit

Floating point

PDIP_NLP, SQP_NLP

17.2.11. Code generation server

By default, code generation requests are routed to embotech’s main FORCESPRO server (https://forces.embotech.com) which always provides the most up-to-date release of FORCESPRO. To send a code generation request to a custom 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';

17.2.12. Enforcing solver regeneration

In order to avoid unnecessary calls to the code-generation 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;

In that case, the codegen server is always contacted to re-generate a new solver.

17.2.13. Overwriting existing solvers

When a new solver is generated with the same name as an existing solver one can control the overwriting behavior by setting the field overwrite as outlined in Table 17.7.

Table 17.7 Overwrite existing solver options

overwrite

Result

0

Never overwrite.

1

Always overwrite.

2 (default)

Ask to overwrite.

17.2.16. 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;

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.

17.2.17. MATLAB network communications

From version 5.0.0, the MATLAB client will perform connections to a REST interface for communicating with the FORCESPRO codegen server.

To revert to an old method, either set:

% WSDL connection
codeoptions.server_connection = ForcesWeb.ServerConnections.WSDL;
% WSDL legacy connection
codeoptions.server_connection = ForcesWeb.ServerConnections.WSDL_legacy;

or change it by editing the FORCESPRO client. To do so, please edit the +ForcesWeb/defaultServerConnection.m function so that it returns the selected ForcesWeb.ServerConnections value.

Important

Setting the codeoptions.server_connection option will override the value in +ForcesWeb/defaultServerConnection.m

The connections are performed over HTTPS. For the server verification, a check via SSL certificates is performed. To perform this verification, a client certificate is needed on the client side. This certificate, if missing, can be acquired by accessing the FORCESPRO server via a browser (the browser will have the option to download the client certificate – the certificate with the entire chain of certificates should be selected). MATLAB stores the client certificates to be used for this verification in an installation specific file. In case this file cannot be accessed/changed due to admin rights, it is possible to manually select the files that will be checked for the certificates for the FORCESPRO connections. This can be done by editing the file: +ForcesWeb/SSLCertificateFiles.m in the FORCESPRO client. This file contains the entries:

% path to certificates file that contains client certificate to authenticate for codegen server. To use default set to empty
codegen = '';
% path to certificates file that contains client certificate to authenticate for storage server. To use default set to empty
storage = '';

The user can download the certificates with the entire chain from the codegen server (by default https://forces.embotech.com but also depends on the server the user has selected) and from the storage server (https://forcesblob.embotech.com) and then assign the path to those certificate files to the file above. This way, MATLAB will use those files for authentication instead of the default file.

Note

In most systems, the default certificate file contains the most common certificates. Therefore, for most users these certificate changes are not required. The above file is provided for cases of strict IT configurations in which certificates need to be specifically enabled for a successful verification.

17.2.18. Python network communications

From version 5.0.0, the Python client will perform connections to a REST interface for communicating with the FORCESPRO codegen server.

To revert to the old method, either set:

# WSDL connection
codeoptions.server_connection = WSDL

or change it by editing the FORCESPRO client. To do so, please edit the default_forcespro_connection.py function so that it returns the selected server_connections value.

Important

Setting the codeoptions.server_connection option will override the value in default_forcespro_connection.py

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.

The connections are performed over HTTPS. For the server verification, a check via SSL certificates is performed. To perform this verification, a client certificate is needed on the client side. This certificate, if missing, can be acquired by accessing the FORCESPRO server via a browser (the browser will have the option to download the client certificate – the certificate with the entire chain of certificates should be selected). Python reads a system specific file to select the client certificates to be used for this verification. In case this file cannot be accessed/changed due to admin rights, it is possible to manually select the files that will be checked for the certificates for the FORCESPRO connections. This can be done by editing the file: forcespro_certificates.py in the FORCESPRO client. This file contains the entries:

# path to certificates file that contains client certificate to authenticate for codegen server (path separator on Windows is "\\"). To use default set to empty
codegen=""
# path to certificates file that contains client certificate to authenticate for storage (path separator on Windows is "\\"). To use default set to empty
storage=""

The user can download the certificates with the entire chain from the codegen server (by default https://forces.embotech.com but also depends on the server the user has selected) and from the storage server (https://forcesblob.embotech.com) and then assign the path to those certificate files to the file above. This way, Python will use those files for authentication instead of the default file.

Note

In most systems, the default certificate file contains the most common certificates. Therefore, for most users these certificate changes are not required. The above file is provided for cases of strict IT configurations in which certificates need to be specifically enabled for a successful verification.

17.2.19. 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 17.8. 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.

Table 17.8 Target platforms supported by FORCESPRO

platform

Description

Portal Selection

'Generic' (default)

For the architecture of the host platform.

'x86_64' (Engineering Node)

'x86_64'

For x86_64 based 64-bit platforms (detected OS).

'x86_64'

'x86'

For x86 based 32-bit platforms (detected OS).

'x86'

'Win-x86_64'

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

'x86_64'

'Win-x86'

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

'x86'

'Win-MinGW-x86_64'

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

'x86_64'

'Win-MinGW-x86'

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

'x86'

'Mac-x86_64'

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

'x86_64'

'Gnu-x86_64'

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

'x86_64'

'Gnu-x86'

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

'x86'

'Docker-Gnu-x86_64'

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

'Docker-Gnu-x86_64'

'Docker-Gnu-x86'

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

'Docker-Gnu-x86'

'ARM-Generic'

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

'ARM-Generic-Gnu'

'ARM-Generic64'

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

'ARM-Generic64-Gnu'

'Integrity-ARM32'

For ARM Cortex 32-bit processors using the Integrity toolchain.

'Integrity-ARM32'

'Integrity-ARM64'

For ARM Cortex 64-bit processors using the Integrity toolchain.

'Integrity-ARM64'

'ARM Cortex-M3'

For ARM Cortex M3 32-bit processors.

'ARM-Cortex-M3'

'ARM-Cortex-M4-NOFPU'

For the ARM Cortex M4 32-bit processors without a floating-point unit.

'ARM-Cortex-M4'

'ARM-Cortex-M4'

For the ARM Cortex M4 32-bit processors with a floating-point unit.

'ARM-Cortex-M4'

'ARM-Cortex-A7'

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

'ARM-Cortex-A7'

'ARM-Cortex-A8'

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

'ARM-Cortex-A8'

'ARM-Cortex-A9'

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

'ARM-Cortex-A9'

'ARM-Cortex-A15'

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

'ARM-Cortex-A15'

'ARM-Cortex-A53'

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

'ARM-Cortex-A53'

'ARM-Cortex-A72'

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

'ARM-Cortex-A72'

'TI-Cortex-A15'

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

'TI-Cortex-A15'

'NVIDIA-Cortex-A57'

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

'NVIDIA-Cortex-A57'

'AARCH-Cortex-A53'

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

'AARCH-Cortex-A53'

'AARCH-Cortex-A57'

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

'AARCH-Cortex-A57'

'AARCH-Cortex-A72'

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

'AARCH-Cortex-A72'

'PowerPC'

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

'PowerPC-Gnu'

'PowerPC64'

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

'PowerPC64-Gnu'

'MinGW32'

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

'x86'

'MinGW64'

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

'x86_64'

'dSPACE-MABII'

For the dSPACE MicroAutoBox II real-time system (supports Microtec compiler).

'dSPACE-MABII-Microtec'

'dSPACE-MABIII'

For the dSPACE MicroAutoBox III real-time system (supports Gcc compiler).

'dSPACE-MABIII-Gcc'

'dSPACE-MABXII'

For the dSPACE MicroAutoBox II real-time system (supports Microtec compiler).

'dSPACE-MABII-Microtec'

'dSPACE-MABXIII'

For the dSPACE MicroAutoBox III real-time system (supports Gcc compiler).

'dSPACE-MABIII-Gcc'

'dSPACE-AutoBox'

For the dSPACE AutoBox real-time system (supports Gcc compiler).

'dSPACE-AutoBox-Gcc'

'dSPACE-MicroLabBox'

For the dSPACE MicroLabBox real-time system (supports Gcc compiler).

'dSPACE-MicroLabBox-Gcc'

'dSPACE-SCALEXIO'

For the dSPACE SCALEXIO real-time system (supports Gcc compiler).

'dSPACE-SCALEXIO-Gcc'

'Speedgoat-x86'

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

'Speedgoat-x86'

'Speedgoat-x64'

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

'Speedgoat-x64'

'Speedgoat-QNX'

For Speedgoat 64-bit real-time platforms (supports MATLAB Releases 2020b onwards).

'Speedgoat-QNX'

'Speedgoat-Legacy-x86'

For Speedgoat Mobile 32-bit real-time platforms (supports Microsoft compiler and MATLAB Releases 2018a and earlier).

'Speedgoat-x86'

'NI-cRIO'

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

'NI-cRIO'

'AURIX'

For Infineon AURIX platforms.

on special request

'IAtomE680_Bachmann'

For Bachmann PLC platforms (supports VxWorks compiler).

'IAtomE680-VxWorks'

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.

17.2.19.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.

17.2.19.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'

17.2.19.3. SIMD instructions

On x86-based host platforms, one can enable the sse field to accelerate the execution of the solver

codeoptions.sse = 1;

On x86-based 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;

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 run-time, one can explicitly disable avx by setting:

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;

On ‘NVIDIA-Cortex-A57’, ‘AARCH-Cortex-A53’, ‘AARCH-Cortex-A57’ and ‘AARCH-Cortex-A72’ 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 behavior 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;

For double precision, the options are

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;

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;

and a warning message is raised by the MATLAB client.

Note

From version 1.9.0, ARMv8-A NEON instructions are generated. Hence, target platforms based on ARMv7 and previous versions are currently not supported.

17.2.19.4. C Standard

FORCESPRO generates solver code in C, which is intended to adhere to the C99 standard. If you are using a compiler that is meant to compile C90 code and fails to compile your FORCESPRO solver, try enabling the following option:

codeoptions.c90 = 1;

This does not guarantee the solver code to stick to the C90 standard, but enables a few adjustments that may render compilation with those legacy compilers possible.

17.2.20. Tips for solving QPs in single precision

Solving QPs in single precision can be rather challenging, i.e. non-converging solves are likely to occur due to the lack of accuracy. In order to mitigate this undesirable behavior, 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 = 1e-5;  % Tolerance on pivot in factorization
codeoptions.regularize.delta = 5e-3;    % On-the-fly regularization coefficient in factorization
codeoptions.regularize.epsilon2 = 1e-5; % Tolerance on pivot in factorization
codeoptions.regularize.delta2 = 5e-3;   % On-the-fly regularization coefficient in factorization

codeoptions.accuracy.ineq = 1e-4;       % infinity norm of residual for inequalities
codeoptions.accuracy.eq = 1e-4;         % infinity norm of residual for equalities
codeoptions.accuracy.mu = 1e-6;         % absolute duality gap
codeoptions.accuracy.rdgap = 1e-4;      % relative duality gap := (pobj-dobj)/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.

17.2.21. 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;

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.

17.2.22. 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;

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;

17.2.23. 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 jumps

  • Operations 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; % overridden if uselocalsall is enabled
codeoptions.optimize_uselocalssimple = 1; % overridden if uselocalsheavy is enabled
codeoptions.optimize_operationsrearrange = 1;
codeoptions.optimize_loopunrolling = 1;
codeoptions.optimize_enableoffset = 1;

17.2.24. Optimizations for small problems

For small convex problems with a short control horizon and/or a small number of controls, the sparse FORCESPRO multistage formulation might be disadvantageous compared to a standard QCQP form. For that situation, FORCESPRO offers the option

codeoptions.condense = 1;

which automatically maps the sparse problem to a dense problem with all dynamic states eliminated. Please consider the requirements described in Condensing (automatic state elimination) for being able to use this option with your formulation.

17.2.25. Dump problem formulation and data

The FORCESPRO client provides a built-in tool to dump the problem formulation to reproduce the exact same solver for future reference. This tool is explained in detail in Section 20. In case you want to use the so-called Legacy dumps (which are only supported by the MATLAB client), you can turn on those dumps the setting:

codeoptions.dump_formulation = 1;

Furthermore, you can dump problem structs containing the runtime parameters from C as described in Section 20. This tool is enabled for the host or/and the target platform by setting:

codeoptions.serializeCParamsHost = 1;
codeoptions.serializeCParamsTarget = 1;

Additionally, you can ask FORCESPRO to generate a MEX function for performing the above-mentioned serialization of runtime parameters by setting:

codeoptions.serializeCParamsMex = 1;

17.2.26. Identifying FORCESPRO solver

In order to be able to identity and distinguish the different FORCESPRO solvers, each solver is assigned a unique GUID (32 hex characters) called Solver Id. The Solver Id is available for all solvers in the header file as a C comment in the format:

/* Solver Id: xxxxxxxx-xxxx-xxxx-xxxx-xxxxxxxxxxxx */

This Solver Id is also available during runtime in the different FORCESPRO interfaces. It is returned from the solver as an unsigned integer array of length 8 each element of which represents 4 hex characters of the full GUID. These elements can be used to generate the GUID as:

solver_id[0] /* hex: a1a2a3a4 */
solver_id[1] /* hex: b1b2b3b4 */
solver_id[2] /* hex: c1c2c3c4 */
solver_id[3] /* hex: d1d2d3d4 */
solver_id[4] /* hex: e1e2e3e4 */
solver_id[5] /* hex: g1g2g3g4 */
solver_id[6] /* hex: h1h2h3h4 */
solver_id[7] /* hex: i1i2i3i4 */

/* Solver Id: a1a2a3a4 b1b2b3b4 - c1c2c3c4 - d1d2d3d4 - e1e2e3e4 - g1g2g3g4 h1h2h3h4 i1i2i3i4 */

The Solver Id is available during runtime through the info struct when using the following interfaces: * C Interface * MATLAB MEX Interface * Python Interface * MPC Toolbox Plugin MATLAB Interface (with the exception of the dense QPs)

The Solver Id is available during runtime as a Simulink Block output port when using the following interfaces: * Simulink Interface (when solver is generated with codeoptions.showinfo = 1) * MPC Toolbox Plugin Simulink Interface (the output port for the Solver Id can be enabled via the Block Parameters of the MPC Toolbox Plugin FORCESPRO Simulink Block – with the exception of the dense QPs)

17.3. High-level interface options

The FORCESPRO NLP solver of the high-level interface implements a nonlinear barrier interior-point method. We will now discuss how to change several parameters in the solver.

17.3.1. Integrators

When providing the continuous dynamics the user must select a particular integrator by setting nlp.integrator.type as outlined in Table 17.9.

Table 17.9 Integrators options

nlp.integrator.type

Type

Order

'ForwardEuler'

Explicit Euler Method

1

'ERK2'

Explicit Runge-Kutta

2

'ERK3'

Explicit Runge-Kutta

3

'ERK4' (default)

Explicit Runge-Kutta

4

'BackwardEuler'

Implicit Euler Method

1

'IRK2'

Implicit Runge-Kutta

2

'IRK4'

Implicit Runge-Kutta

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;

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 and their legacy variant cannot be used in combination with CasADi MX expressions (see Automatic differentiation expression class).

17.3.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;

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 = 0;

17.3.1.2. Code-generated 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';

  • vde, which can be triggered by settings the following option:

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

The vde option is likely to change the numerical behavior 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 behavior. Also, the RK3 integration method is currently not supported with the vde option.

17.3.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

\begin{align} \dot{x}_1 &= A_1 x_1 + B_1 u \\ \dot{x}_2 &= c_2(x_1,x_2,u). \end{align}

Here \(A_1\) and \(B_1\) denote constant matrices and \(c_2\) denotes a 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 Code-generated 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;

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 right-hand-side \(c\) of the ODE describing the dynamics of the system is given by

\begin{equation*} c(x,u) = \begin{pmatrix} x_2 \\ x_1 \\ \cos(x_1) + \sin(x_2) + x_3 + u \end{pmatrix}, \end{equation*}

one would have to specify

model.linInIdx = [1, 2];

Note the 1-based indexing in MATLAB and the 0-based 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.

Moreover, note that linear subsystems currently cannot be exploited when using CasADi MX expressions (see Automatic differentiation expression class).

17.3.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 = 1e-5; % inf norm tol. on stationarity
codeoptions.nlp.TolEq = 1e-6;   % tol. on equality constraints
codeoptions.nlp.TolIneq = 1e-6; % tol. on inequality constraints
codeoptions.nlp.TolComp = 1e-6; % tol. on complementarity

All tolerances are computed using the infinity norm \(\lVert \cdot \rVert_\infty\).

17.3.3. Barrier strategy

The strategy for updating the barrier parameter is set using the field:

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.

17.3.4. Hessian approximation

The way the Hessian of the Lagrangian function is computed can be set using the field:

codeoptions.nlp.hessian_approximation = 'bfgs';

FORCESPRO currently supports BFGS updates ('bfgs') (default) and Gauss-Newton approximation ('gauss-newton'). Exact Hessians will be supported in a future version. Read the subsequent sections for the corresponding Hessian approximation method of your choice.

17.3.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.:

model.bfgs_init = diag([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. It’s also possible to specify the initialization of the BFGS per stage (especially when the dimensions of the problem vary per stage):

model.bfgs_init = cell(1, model.N);
for i = 1:model.N
  model.bfgs_init{i} = diag([0.1, 10, 4]);
end

Finally, it’s possible to specify a initialization separately for the final stage:

model.bfgs_initN = diag([0.1, 10, 4]); % overrides model.bfgs_init{model.N}

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 re-initialized is set using:

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 user-defined BFGS initialization (i.e. not setting model.bfgs_init) and using the BFGS matrix reached upon convergence as an initialization. One can export the BFGS matrix by setting

% diagonal of BFGS
codeoptions.exportBFGS = 1;
% lower triangular of BFGS
codeoptions.exportBFGS = 2;

Instead of specifying the BFGS initialization at codegen one can also specify it at run-time. In order to do this one should set

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

\[\begin{split}\begin{pmatrix} a_1 & 0 & 0 \\ 0 & a_2 & 0 \\ 0 & 0 & a_3 \end{pmatrix}\end{split}\]

for constants \(a_1,a_2,a_3 > 0\) as the BFGS initialization at stage 6 out of 50 stages in total, one would specify

problem.BFGSinitLower06 = [a_1, 0, a_2, 0, 0, a_3];

17.3.4.2. Gauss-Newton options

For problems that have a least squares objective, i.e. the cost function can be expressed by a vector-valued function \(r_k : \mathbb{R}^n \rightarrow \mathbb{R}^m\) which implicitly defines the objective function as:

\[f_k(z_k,p_k) = \frac{1}{2} \lVert r_k(z_k,p_k) \rVert_2^2 \,,\]

the Gauss-Newton approximation of the Hessian is given by:

\[\nabla_{xx}^2 L_k \approx \nabla r_k(z_k,p_k) \nabla r_k(z_k,p_k)^\top\]

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 least-squares objective (note that \(r\) is a vector-valued 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)];

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.

17.3.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 fraction-to-boundary scaling
codeoptions.nlp.ftbr_scaling = 0.9900;

17.3.6. Regularization

To avoid ill-conditioned saddle point systems, FORCESPRO employs two different types of regularization, static and dynamic regularization.

17.3.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:

\[\begin{split}\delta_w & = \eta_w \min(\mu,\lVert c(x) \rVert))^{\beta_w} \cdot (i+1)^{-\gamma_w} + \delta_{w,\min} \\ \delta_c & = \eta_c \min(\mu,\lVert c(x) \rVert))^{\beta_c} \cdot (i+1)^{-\gamma_c} + \delta_{c,\min} \\\end{split}\]

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 inertia-correcting 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 = 1e-4;
codeoptions.nlp.reg_beta_dw = 0.8;
codeoptions.nlp.reg_min_dw = 1e-9;
codeoptions.nlp.reg_gamma_dw = 1.0/3.0;

codeoptions.nlp.reg_eta_dc = 1e-4;
codeoptions.nlp.reg_beta_dc = 0.8;
codeoptions.nlp.reg_min_dc = 1e-9;
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.

17.3.6.2. Dynamic regularization

Dynamic regularization regularizes the matrix on-the-fly 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 = 1e-12; % (for Hessian approx.)
codeoptions.regularize.delta = 4e-6; % (for Hessian approx.)
codeoptions.regularize.epsilon2 = 1e-14; % (for Normal eqs.)
codeoptions.regularize.delta2 = 1e-14; % (for Normal eqs.)

17.3.7. Linear system solver

The interior-point method solves a linear system to find a search direction at every iteration. FORCESPRO NLP offers the following four linear solvers:

  • 'normal_eqs' (default): Solving the KKT system in normal equations form.

  • 'symm_indefinite': improved variant of 'symm_indefinite_legacy' introduced in FORCESPRO version 5.0.0; roughly as efficient as normal_eqs but more robust.

  • 'symm_indefinite_fast': Solving the KKT system in augmented / symmetric indefinite form, using regularization and positive definite Cholesky factorizations only. This is often the fastest solver but may be less numerical stable than symm_indefinite.

  • 'symm_indefinite_legacy': Solving the KKT system in augmented / symmetric indefinite form; may be removed in a future release

The linear system solver can be selected by setting the following field:

codeoptions.nlp.linear_solver = 'symm_indefinite';

It is recommended to try different linear solvers as the robustness and speed of the solvers are problem-dependent. The overall most robust method is symm_indefinite, which is also very efficient. For certain problems normal_eqs and 'symm_indefinite_fast' may be slightly faster than symm_indefinite but possibly also slightly less numerically stable.

Note

Independent of the linear system solver choice, the generated code is always library-free and statically allocated, i.e. it can be embedded anywhere.

Note

From FORCESPRO version 5.0.0 onwards, the option symm_indefinite refers to an improved version; use symm_indefinite_legacy to restore the previous default.

The 'normal_eqs' solver is the standard FORCESPRO linear system solver based on a full reduction of the KKT system (the so-called normal equations form). It works well for standard problems, especially convex problems or nonlinear problems where the BFGS or Gauss-Newton approximations of the Hessian are numerically sufficiently well conditioned.

The 'symm_indefinite' solver is numerically more robust than 'normal_eqs' and symm_indefinite_fast and typically similarly efficient. It is an improved variant of the 'symm_indefinite_legacy'. Furthermore, it implements iterative refinement which further improves numerical stability (see Iterative refinement).

The 'symm_indefinite_fast' solver is typically the fastest solver. Currently only used for receding-horizon/MPC-like problems where dimensions of all stages are equal (except for the first and last stage, those are handled separately).

The 'symm_indefinite_legacy' solver is the most robust one, but currently replaced by an at least equally robust improved variant.

17.3.7.1. Iterative refinement

The linear solver 'symm_indefinite' supports iterative refinement to further improve numerical stability. Iterative refinement is recommended for problems that don’t converge due to numerical issues but can be safely disabled (default) for most problems. In order to enable iterative refinement, set codeoptions.nlp.refinement_steps to the desired number of steps > 0. Two types of iterative refinement are implemented which can be selected by setting codeoptions.nlp.refinement_type as outlined in Table 17.10.

Table 17.10 Options for setting iterative refinement type

nlp.refinement_type

Description

0

Includes additional modification (default)

1

Strictly based on the original linear system

17.3.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, two different AD tools (or four different AD tool versions) are supported that can be chosen by means of the setting nlp.ad_tool as summarized in Table 17.11.

Table 17.11 Automatic differentiation tool options

nlp.ad_tool

Tool

URL

'casadi'

CasADi (as in path or v3.5.5)

CasADi

'casadi-2.4.2'

CasADi v2.4.2

CasADi

'casadi-3.5.1'

CasADi v3.5.1

CasADi

'casadi-3.5.5'

CasADi v3.5.5

CasADi

'symbolic-math-tbx'

MathWorks Symbolic Math Toolbox

MathWorks

Note that MathWorks Symbolic Math Toolbox requires an additional license, which is why the default option is set to

codeoptions.nlp.ad_tool = 'casadi';

Also note that the use of implicit integrators BackwardEuler, IRK2 and IRK4 (see Section 17.3.1) currently still rely on using the CasADi AD tool.

17.3.9. Automatic differentiation expression class

The AD tool CasADi supports two different types of symbolic expressions, called SX and MX. While SX expressions tend to be leaner and offer best performance for mathematical “operations that are naturally written as a sequence of scalar operations” (see CasADi’s symbolic framework)), MX expressions offer advantages for matrix operations.

FORCESPRO only supported SX expressions prior to version 6.1.0 and still uses them by default. However, starting with version 6.1.0, FORCESPRO offers to make use of MX expressions by setting the following codeoption:

codeoptions.nlp.ad_expression_class = 'MX';

All supported values of that codeoption nlp.ad_expression_class are summarized in Table 17.12. Note that this option is only supported when using CasADi as AD Tool and that MX expressions are only supported for CasADi versions 3.5.1 and 3.5.5.

Table 17.12 Automatic differentiation expression classes

nlp.ad_expression_class

Supported AD Tools

'SX'

CasADi v2.4.2, v3.5.1, v3.5.5

'MX'

CasADi v3.5.1, v3.5.5

Note

Use CasADi MX expressions with care. In particular, when switching your FORCESPRO formulation between using SX and MX expressions, minor numerical differences may arise that can lead to (usually minor) differences in number of iterations or optimal solutions.

Note

Use of CasADi MX expressions is currently not supported in combination with any of the following FORCESPRO features: legacy and symbolic dump tool, implicit legacy integrators, linear subsystem exploitation, MINLP formulations and inside the MathWorks MPC Plugins.

17.3.10. Re-use 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;

Note that using this option may be overly conservative and lead to duplicated callback code for different stages that are actually identical.

17.3.11. Safety checks

By default, the output of the function evaluations is checked for the presence of NaNs or INFs 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;

17.4. Convex branch-and-bound options

The settings of the FORCESPRO mixed-integer branch-and-bound convex solver are accessed through the codeoptions.mip struct. It is worthwhile to explore different values for the settings in Table 17.13, as they might have a severe impact on the performance of the branch-and-bound procedure.

Note

All the options described below are currently not available with the FORCESPRO nonlinear solver. For mixed-integer nonlinear programs and the available options, please have a look at paragraph Mixed-integer nonlinear solver.

Table 17.13 Branch-and-bound options

Setting

Values

Default

mip.timeout

Any value \(\geq 0\)

31536000 (1 year)

mip.mipgap

Any value \(\geq 0\)

0

mip.branchon

'mostAmbiguous', 'leastAmbiguous'

'mostAmbiguous'

mip.stageinorder

0 (OFF), 1 (ON)

1 (ON)

mip.explore

'bestFirst', 'depthFirst'

'bestFirst'

mip.inttol

Any value \(> 0\)

1E-5

mip.queuesize

Any integer value \(\geq 0\)

1000

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 sub-optimality after which the search shall be terminated. For example, a value of 0.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: Stage-in-order 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 are 0 for OFF and 1 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 below inttol 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 an exitflag value of -2 and returns the best solution found so far.

17.5. Solve methods

As a default optimization method the primal-dual interior-point method is used. Several other methods are available. To change the solve method set the solvemethod field as outlined in Table 17.14.

Table 17.14 Solve methods

solvemethod

Method

Description

'PDIP_NLP'

Nonlinear Primal-Dual Interior-Point method

The Nonlinear Primal-Dual Interior-Point method is the most stable and robust method for most nonlinear problems.

'SQP_NLP'

Sequential Quadratic Programming method

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

'PDIP'

Primal-Dual Interior-Point method

The Primal-Dual Interior-Point method is the most stable and robust method for most convex problems.

'ADMM'

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.

'DFG'

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.

'FG'

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.

17.5.1. Primal-Dual Interior-Point Method

The Primal-Dual Interior-Point 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.

17.5.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 17.15.

Table 17.15 PDIP solver initialization

init

Method

Initialization method

0 (default)

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\).

1

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\).

2

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\).

17.5.1.2. Initial Complementary Slackness

The default value for \(\mu_0\) is \(10^6\). To use a different value, use:

codeoptions.mu0 = 10;

17.5.1.3. Accuracy Requirements

The accuracy for which FORCESPRO returns the OPTIMAL flag can be set as follows:

codeoptions.accuracy.ineq = 1e-6;  % infinity norm of residual for inequalities
codeoptions.accuracy.eq = 1e-6;    % infinity norm of residual for equalities
codeoptions.accuracy.mu = 1e-6;    % absolute duality gap
codeoptions.accuracy.rdgap = 1e-4; % relative duality gap := (pobj-dobj)/pobj

17.5.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;

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 = 1e-8;
codeoptions.linesearch.maxstep = 0.995;

17.5.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 = 1e-13; % if pivot element < epsilon …
codeoptions.regularize.delta = 1e-8;    % then set it to delta

17.5.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

where variant 1 is as follows:

\begin{align*} \text{minimize} \quad & \frac{1}{2} y^\top H y + f^\top y \\ \text{subject to} \quad & Dy=c \\ & \underline{z} \leq z \leq \bar{z} \\ & y = z \end{align*}

and variant 2 is as follows:

\begin{align*} \text{minimize} \quad & \frac{1}{2} y^\top H y + f^\top y \\ \text{subject to} \quad & Dy=c \\ & A y = z \\ & z \leq b \end{align*}

17.5.2.1. Accuracy requirements

The accuracy for which FORCESPRO returns the OPTIMAL flag can be set as follows:

codeoptions.accuracy.consensus = 1e-3;  % infinity norm of the consensus equality
codeoptions.accuracy.dres = 1e-3;    % infinity norm of the dual residual

Note that, in contrast to primal-dual interior-point 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 achieve the same accuracy as interior-point 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.

17.5.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;

In some cases it may be possible to let FORCESPRO choose the value \(\rho\) automatically. To enable this feature set:

codeoptions.ADMMautorho = 1;

Please note that this does not guarantee that the choice of \(\rho\) will be optimal.

ADMM can also include an ‘over-relaxation’ 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;

17.5.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.

17.5.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.

17.5.4.1. Accuracy requirements

The accuracy for which FORCESPRO returns the OPTIMAL flag can be set as follows:

codeoptions.accuracy.gmap = 1e-5;  % 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 first-order 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 interior-point methods.

17.5.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;

In some cases it may be possible to let FORCESPRO choose the step size automatically. To enable this feature set:

codeoptions.FGautostep = 1;

17.5.4.3. Warm starting

The performance of the fast gradient method can be greatly influenced by the way the variables are initialized. Unlike with interior-point 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;

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.