Welcome to CasADi’s documentation!

1. Introduction

CasADi is an open-source software tool for numerical optimization in general and optimal control (i.e. optimization involving differential equations) in particular. The project was started by Joel Andersson and Joris Gillis while PhD students at the Optimization in Engineering Center (OPTEC) of the KU Leuven under supervision of Moritz Diehl.

This document aims at giving a condensed introduction to CasADi. After reading it, you should be able to formulate and manipulate expressions in CasADi’s symbolic framework, generate derivative information efficiently using algorithmic differentiation, to set up, solve and perform forward and adjoint sensitivity analysis for systems of ordinary differential equations (ODE) or differential-algebraic equations (DAE) as well as to formulate and solve nonlinear programs (NLP) problems and optimal control problems (OCP).

CasADi is available for C++, Python and MATLAB/Octave with little or no difference in performance. In general, the Python API is the best documented and is slightly more stable than the MATLAB API. The C++ API is stable, but is not ideal for getting started with CasADi since there is limited documentation and since it lacks the interactivity of interpreted languages like MATLAB and Python. The MATLAB module has been tested successfully for Octave (version 4.0.2 or later).

1.1. What CasADi is and what it is not

CasADi started out as a tool for algorithmic differentiation (AD) using a syntax borrowed from computer algebra systems (CAS), which explains its name. While AD still forms one of the core functionalities of the tool, the scope of the tool has since been considerably broadened, with the addition of support for ODE/DAE integration and sensitivity analysis, nonlinear programming and interfaces to other numerical tools. In its current form, it is a general-purpose tool for gradient-based numerical optimization – with a strong focus on optimal control – and CasADi is just a name without any particular meaning.

It is important to point out that CasADi is not a conventional AD tool, that can be used to calculate derivative information from existing user code with little to no modification. If you have an existing model written in C++, Python or MATLAB/Octave, you need to be prepared to reimplement the model using CasADi syntax.

Secondly, CasADi is not a computer algebra system. While the symbolic core does include an increasing set of tools for manipulate symbolic expressions, these capabilities are very limited compared to a proper CAS tool.

Finally, CasADi is not an “optimal control problem solver”, that allows the user to enter an OCP and then gives the solution back. Instead, it tries to provide the user with a set of “building blocks” that can be used to implement general-purpose or specific-purpose OCP solvers efficiently with a modest programming effort.

1.2. Help and support

If you find simple bugs or lack some feature that you think would be relatively easy for us to add, the simplest thing is simply to write to the forum, located at http://forum.casadi.org/. We check the forum regularly and try to respond as quickly as possible. The only thing we expect for this kind of support is that you cite us, cf. Section 1.3, whenever you use CasADi in scientific work.

If you want more help, we are always open for academic or industrial cooperation. An academic cooperation usually take the form of a co-authorship of a peer reviewed paper, and an industrial cooperation involves a negotiated consulting contract. Please contact us directly if you are interested in this.

1.3. Citing CasADi

If you use CasADi in published scientific work, please cite the following:

@Article{Andersson2018,
  Author = {Joel A E Andersson and Joris Gillis and Greg Horn
            and James B Rawlings and Moritz Diehl},
  Title = {{CasADi} -- {A} software framework for nonlinear optimization
           and optimal control},
  Journal = {Mathematical Programming Computation},
  Year = {In Press, 2018},
}

1.4. Reading this document

The goal of this document is to make the reader familiar with the syntax of CasADi and provide easily available building blocks to build numerical optimization and dynamic optimization software. Our explanation is mostly program code driven and provides little mathematical background knowledge. We assume that the reader already has a fair knowledge of theory of optimization theory, solution of initial-value problems in differential equations and the programming language in question (C++, Python or MATLAB/Octave).

We will use Python and MATLAB/Octave syntax side-by-side in this guide, noting that the Python interface is more stable and better documented. Unless otherwise noted, the MATLAB/Octave syntax also applies to Octave. We try to point out the instances where has a diverging syntax. To facilitate switching between the programming languages, we also list the major differences in Chapter 10.

2. Obtaining and installing

CasADi is an open-source tool, available under LGPL license, which is a permissive license that allows the tool to be used royalty-free also in commercial closed-source applications. The main restriction of LGPL is that if you decide to modify CasADi’s source code as opposed to just using the tool for your application, these changes (a “derivative-work” of CasADi) must be released under LGPL as well.

The source code is hosted on GitHub and has a core written in self-contained C++ code, relying on nothing but the C++ Standard Library. Its front-ends to Python and MATLAB/Octave are full-featured and auto-generated using the tool SWIG. These front-ends are unlikely to result in noticeable loss of efficiency. CasADi can be used on Linux, OS X and Windows.

For up-to-date installation instructions, visit CasADi’s install section: http://install.casadi.org/

3. Symbolic framework

At the core of CasADi is a self-contained symbolic framework that allows the user to construct symbolic expressions using a MATLAB inspired everything-is-a-matrix syntax, i.e. vectors are treated as n-by-1 matrices and scalars as 1-by-1 matrices. All matrices are sparse and use a general sparse format – compressed column storage (CCS) – to store matrices. In the following, we introduce the most fundamental classes of this framework.

3.1. The SX symbolics

The SX data type is used to represent matrices whose elements consist of symbolic expressions made up by a sequence of unary and binary operations. To see how it works in practice, start an interactive Python shell (e.g. by typing ipython from a Linux terminal or inside a integrated development environment such as Spyder) or launch MATLAB’s or Octave’s graphical user interface. Assuming CasADi has been installed correctly, you can import the symbols into the workspace as follows:

from casadi import *
import casadi.*

Now create a variable x using the syntax:

x = MX.sym("x")
x
x = MX.sym('x')
x =

x

This creates a 1-by-1 matrix, i.e. a scalar containing a symbolic primitive called x. This is just the display name, not the identifier. Multiple variables can have the same name, but still be different. The identifier is the return value. You can also create vector- or matrix-valued symbolic variables by supplying additional arguments to SX.sym:

y = SX.sym('y',5)
Z = SX.sym('Z',4,2)
[y_0, y_1, y_2, y_3, y_4] 
[[Z_0, Z_4], 
 [Z_1, Z_5], 
 [Z_2, Z_6], 
 [Z_3, Z_7]]
y = SX.sym('y',5)
Z = SX.sym('Z',4,2)
y =

[y_0, y_1, y_2, y_3, y_4]

Z =


[[Z_0, Z_4], 
 [Z_1, Z_5], 
 [Z_2, Z_6], 
 [Z_3, Z_7]]

which creates a 5-by-1 matrix, i.e. a vector, and a 4-by-2 matrix with symbolic primitives, respectively.

SX.sym is a (static) function which returns an SX instance. When variables have been declared, expressions can now be formed in an intuitive way:

f = x**2 + 10
f = sqrt(f)
print(f)
sqrt((10+sq(x)))
f = x^2 + 10;
f = sqrt(f);
display(f)
f =

sqrt((10+sq(x)))

You can also create constant SX instances without any symbolic primitives:

  • B1 = SX.zeros(4,5): A dense 4-by-5 empty matrix with all zeros
  • B2 = SX(4,5): A sparse 4-by-5 empty matrix with all zeros
  • B4 = SX.eye(4): A sparse 4-by-4 matrix with ones on the diagonal
B1: @1=0, 
[[@1, @1, @1, @1, @1], 
 [@1, @1, @1, @1, @1], 
 [@1, @1, @1, @1, @1], 
 [@1, @1, @1, @1, @1]]
B2: 
[[00, 00, 00, 00, 00], 
 [00, 00, 00, 00, 00], 
 [00, 00, 00, 00, 00], 
 [00, 00, 00, 00, 00]]
B4: @1=1, 
[[@1, 00, 00, 00], 
 [00, @1, 00, 00], 
 [00, 00, @1, 00], 
 [00, 00, 00, @1]]

Note the difference between a sparse matrix with structural zeros and a dense matrix with actual zeros. When printing an expression with structural zeros, these will be represented as 00 to distinguish them from actual zeros 0.

The following list summarizes the most commonly used ways of constructing new SX expressions:

  • SX.sym(name,n,m): Create an \(n\)-by-\(m\) symbolic primitive
  • SX.zeros(n,m): Create an \(n\)-by-\(m\) dense matrix with all zeros
  • SX(n,m): Create an \(n\)-by-\(m\) sparse matrix with all structural zeros
  • SX.ones(n,m): Create an \(n\)-by-\(m\) dense matrix with all ones
  • SX.eye(n): Create an \(n\)-by-\(n\) diagonal matrix with ones on the diagonal and structural zeros elsewhere.
  • SX(scalar_type): Create a scalar (1-by-1 matrix) with value given by the argument. This method can be used explicitly, e.g. SX(9), or implicitly, e.g. 9 * SX.ones(2,2).
  • SX(matrix_type): Create a matrix given a numerical matrix given as a NumPy or SciPy matrix (in Python) or as a dense or sparse matrix (in MATLAB/Octave). In MATLAB/Octave e.g. SX([1,2,3,4]) for a row vector, SX([1;2;3;4]) for a column vector and SX([1,2;3,4]) for a 2-by-2 matrix. This method can be used explicitly or implicitly.
  • repmat(v,n,m): Repeat expression \(v\) \(n\) times vertically and \(m\) times horizontally. repmat(SX(3),2,1) will create a 2-by-1 matrix with all elements 3.
  • (Python only) SX(list): Create a column vector (\(n\)-by-1 matrix) with the elements in the list, e.g. SX([1,2,3,4]) (note the difference between Python lists and MATLAB/Octave horizontal concatenation, which both uses square bracket syntax)
  • (Python only) SX(list of list): Create a dense matrix with the elements in the lists, e.g. SX([[1,2],[3,4]]) or a row vector (1-by-\(n\) matrix) using SX([[1,2,3,4]]).

3.1.1. Note on namespaces

In MATLAB, if the import casadi.* command is omitted, you can still use CasADi by prefixing all the symbols with the package name, e.g. casadi.SX instead of SX, provided the casadi package is in the path. We will not do this in the following for typographical reasons, but note that it is often preferable in user code. In Python, this usage corresponds to issuing import casadi instead of from casadi import *.

Unfortunately, Octave (version 4.0.3) does not implement MATLAB’s import command. To work around this issue, we provide a simple function import.m that can be placed in Octave’s path enabling the compact syntax used in this guide.

3.1.2. Note for C++ users

In C++, all public symbols are defined in the casadi namespace and require the inclusion of the casadi/casadi.hpp header file. The commands above would be equivalent to:

#include <casadi/casadi.hpp>
using namespace casadi;
int main() {
  SX x = SX::sym("x");
  SX y = SX::sym("y",5);
  SX Z = SX::sym("Z",4,2)
  SX f = pow(x,2) + 10;
  f = sqrt(f);
  std::cout << "f: " << f << std::endl;
  return 0;
}

3.2. DM

DM is very similar to SX, but with the difference that the nonzero elements are numerical values and not symbolic expressions. The syntax is also the same, except for functions such as SX.sym, which have no equivalents.

DM is mainly used for storing matrices in CasADi and as inputs and outputs of functions. It is not intended to be used for computationally intensive calculations. For this purpose, use the builtin dense or sparse data types in MATLAB, NumPy or SciPy matrices in Python or an expression template based library such as eigen, ublas or MTL in C++. Conversion between the types is usually straightforward:

C = DM(2,3)

C_dense = C.full()
from numpy import array
C_dense = array(C) # equivalent

C_sparse = C.sparse()
from scipy.sparse import csc_matrix
C_sparse = csc_matrix(C) # equivalent
C = DM(2,3);
C_dense = full(C);
C_sparse = sparse(C);

More usage examples for SX can be found in the example pack at http://install.casadi.org/. For documentation of particular functions of this class (and others), find the “C++ API” on http://docs.casadi.org/ and search for information about casadi::Matrix.

3.3. The MX symbolics

Let us perform a simple operation using the SX above:

x = SX.sym('x',2,2)
y = SX.sym('y')
f = 3*x + y
print(f)
print(f.shape)
@1=3, 
[[((@1*x_0)+y), ((@1*x_2)+y)], 
 [((@1*x_1)+y), ((@1*x_3)+y)]]
(2, 2)
x = SX.sym('x',2,2);
y = SX.sym('y');
f = 3*x + y;
disp(f)
disp(size(f))
@1=3, 
[[((@1*x_0)+y), ((@1*x_2)+y)], 
 [((@1*x_1)+y), ((@1*x_3)+y)]]
   2   2

As you can see, the output of this operation is a 2-by-2 matrix. Note how the multiplication and the addition were performed element-wise and new expressions (of type SX) were created for each entry of the result matrix.

We shall now introduce a second, more general matrix expression type MX. The MX type allows, like SX, to build up expressions consisting of a sequence of elementary operations. But unlike SX, these elementary operations are not restricted to be scalar unary or binary operations (\(\mathbb{R} \rightarrow \mathbb{R}\) or \(\mathbb{R} \times \mathbb{R} \rightarrow \mathbb{R}\)). Instead, the elementary operations that are used to form MX expressions are allowed to be general multiple sparse-matrix valued input, multiple sparse-matrix valued output functions: \(\mathbb{R}^{n_1 \times m_1} \times \ldots \times \mathbb{R}^{n_N \times m_N} \rightarrow \mathbb{R}^{p_1 \times q_1} \times \ldots \times \mathbb{R}^{p_M \times q_M}\).

The syntax of MX mirrors that of SX:

x = MX.sym('x',2,2)
y = MX.sym('y')
f = 3*x + y
print(f)
print(f.shape)
((3*x)+y)
(2, 2)
x = MX.sym('x',2,2);
y = MX.sym('y');
f = 3*x + y;
disp(f)
disp(size(f))
((3*x)+y)
   2   2

Note how the result consists of only two operations (one multiplication and one addition) using MX symbolics, whereas the SX equivalent has eight (two for each element of the resulting matrix). As a consequence, MX can be more economical when working with operations that are naturally vector or matrix valued with many elements. As we shall see in Chapter 4, it is also much more general since we allow calls to arbitrary functions that cannot be expanded in terms of elementary operations.

MX supports getting and setting elements, using the same syntax as SX, but the way it is implemented is very different. Test, for example, to print the element in the upper-left corner of a 2-by-2 symbolic variable:

x = MX.sym('x',2,2)
print(x[0,0])
x[0]
x = MX.sym('x',2,2);
x(1,1)
ans =

x[0]

The output should be understood as an expression that is equal to the first (i.e. index 0 in C++) structurally non-zero element of x, unlike x_0 in the SX case above, which is the name of a symbolic primitive in the first (index 0) location of the matrix.

Similar results can be expected when trying to set elements:

x = MX.sym('x',2)
A = MX(2,2)
A[0,0] = x[0]
A[1,1] = x[0]+x[1]
print('A:', A)
A: (project((zeros(2x2,1nz)[0] = x[0]))[1] = (x[0]+x[1]))
x = MX.sym('x',2);
A = MX(2,2);
A(1,1) = x(1);
A(2,2) = x(1)+x(2);
display(A)
A =

(project((zeros(2x2,1nz)[0] = x[0]))[1] = (x[0]+x[1]))

The interpretation of the (admittedly cryptic) output is that starting with an all zero sparse matrix, an element is assigned to x_0. It is then projected to a matrix of different sparsity and an another element is assigned to x_0+x_1.

Element access and assignment, of the type you have just seen, are examples of operations that can be used to construct expressions. Other examples of operations are matrix multiplications, transposes, concatenations, resizings, reshapings and function calls.

3.4. Mixing SX and MX

You can not multiply an SX object with an MX object, or perform any other operation to mix the two in the same expression graph. You can, however, in an MX graph include calls to a Function defined by SX expressions. This will be demonstrated in Chapter 4. Mixing SX and MX is often a good idea since functions defined by SX expressions have a much lower overhead per operation making it much faster for operations that are naturally written as a sequence of scalar operations. The SX expressions are thus intended to be used for low level operations (for example the DAE right hand side in Section 4.4), whereas the MX expressions act as a glue and enables the formulation of e.g. the constraint function of an NLP (which might contain calls to ODE/DAE integrators, or might simply be too large to expand as one big expression).

3.5. The Sparsity class

As mentioned above, matrices in CasADi are stored using the compressed column storage (CCS) format. This is a standard format for sparse matrices that allows linear algebra operations such as element-wise operations, matrix multiplication and transposes to be performed efficiently. In the CCS format, the sparsity pattern is decoded using the dimensions – the number of rows and number of columns – and two vectors. The first vector contains the index of the first structurally nonzero element of each column and the second vector contains the row index for every nonzero element. For more details on the CCS format, see e.g. Templates for the Solution of Linear Systems on Netlib. Note that CasADi uses the CCS format for sparse as well as dense matrices.

Sparsity patterns in CasADi are stored as instances of the Sparsity class, which is reference-counted, meaning that multiple matrices can share the same sparsity pattern, including MX expression graphs and instances of SX and DM. The Sparsity class is also cached, meaning that the creation of multiple instances of the same sparsity patterns is always avoided.

The following list summarizes the most commonly used ways of constructing new sparsity patterns:

  • Sparsity.dense(n,m): Create a dense \(n\)-by-\(m\) sparsity pattern
  • Sparsity(n,m): Create a sparse \(n\)-by-\(m\) sparsity pattern
  • Sparsity.diag(n): Create a diagonal \(n\)-by-\(n\) sparsity pattern
  • Sparsity.upper(n): Create an upper triangular \(n\)-by-\(n\) sparsity pattern
  • Sparsity.lower(n): Create a lower triangular \(n\)-by-\(n\) sparsity pattern

The Sparsity class can be used to create non-standard matrices, e.g.

print(SX.sym('x',Sparsity.lower(3)))

[[x_0, 00, 00], 
 [x_1, x_3, 00], 
 [x_2, x_4, x_5]]
disp(SX.sym('x',Sparsity.lower(3)))

[[x_0, 00, 00], 
 [x_1, x_3, 00], 
 [x_2, x_4, x_5]]

3.5.1. Getting and setting elements in matrices

To get or set an element or a set of elements in CasADi’s matrix types (SX, MX and DM), we use square brackets in Python and round brackets in C++ and MATLAB. As is conventional in these languages, indexing starts from zero in C++ and Python but from one in MATLAB. In Python and C++, we allow negative indices to specify an index counted from the end. In MATLAB, use the end keyword for indexing from the end.

Indexing can be done with one index or two indices. With two indices, you reference a particular row (or set or rows) and a particular column (or set of columns). With one index, you reference an element (or set of elements) starting from the upper left corner and column-wise to the lower right corner. All elements are counted regardless of whether they are structurally zero or not.

M = SX([[3,7],[4,5]])
print(M[0,:])
M[0,:] = 1
print(M)
[[3, 7]]
@1=1, 
[[@1, @1], 
 [4, 5]]
M = SX([3,7;4,5]);
disp(M(1,:))
M(1,:) = 1;
disp(M)
[[3, 7]]
@1=1, 
[[@1, @1], 
 [4, 5]]

Unlike Python’s NumPy, CasADi slices are not views into the data of the left hand side; rather, a slice access copies the data. As a result, the matrix \(m\) is not changed at all in the following example:

M = SX([[3,7],[4,5]])
M[0,:][0,0] = 1
print(M)

[[3, 7], 
 [4, 5]]

The getting and setting matrix elements is elaborated in the following. The discussion applies to all of CasADi’s matrix types.

Single element access is getting or setting by providing a row-column pair or its flattened index (column-wise starting in the upper left corner of the matrix):

M = diag(SX([3,4,5,6]))
print(M)

[[3, 00, 00, 00], 
 [00, 4, 00, 00], 
 [00, 00, 5, 00], 
 [00, 00, 00, 6]]
M = diag(SX([3,4,5,6]));
disp(M)

[[3, 00, 00, 00], 
 [00, 4, 00, 00], 
 [00, 00, 5, 00], 
 [00, 00, 00, 6]]
print(M[0,0])
print(M[1,0])
print(M[-1,-1])
3
00
6
disp(M(1,1))
disp(M(2,1))
disp(M(end,end))
3
00
6

Slice access means setting multiple elements at once. This is significantly more efficient than setting the elements one at a time. You get or set a slice by providing a (start , stop , step) triple. In Python and MATLAB, CasADi uses standard syntax:

print(M[:,1])
print(M[1:,1:4:2])
[00, 4, 00, 00]

[[4, 00], 
 [00, 00], 
 [00, 6]]
disp(M(:,2))
disp(M(2:end,2:2:4))
[00, 4, 00, 00]

[[4, 00], 
 [00, 00], 
 [00, 6]]

In C++, CasADi’s Slice helper class can be used. For the example above, this means M(Slice(),1) and M(Slice(1,-1),Slice(1,4,2)), respectively.

List access is similar to (but potentially less efficient than) slice access:

M = SX([[3,7,8,9],[4,5,6,1]])
print(M)
print(M[0,[0,3]], M[[5,-6]])

[[3, 7, 8, 9], 
 [4, 5, 6, 1]]
[[3, 9]] [6, 7]
M = SX([3 7 8 9; 4 5 6 1]);
disp(M)
disp(M(1,[1,4]))
disp(M([6,numel(M)-5]))

[[3, 7, 8, 9], 
 [4, 5, 6, 1]]
[[3, 9]]
[[6, 7]]

3.6. Arithmetic operations

CasADi supports most standard arithmetic operations such as addition, multiplications, powers, trigonometric functions etc:

x = SX.sym('x')
y = SX.sym('y',2,2)
print(sin(y)-x)

[[(sin(y_0)-x), (sin(y_2)-x)], 
 [(sin(y_1)-x), (sin(y_3)-x)]]
x = SX.sym('x');
y = SX.sym('y',2,2);
disp(sin(y)-x)

[[(sin(y_0)-x), (sin(y_2)-x)], 
 [(sin(y_1)-x), (sin(y_3)-x)]]

In C++ and Python (but not in MATLAB), the standard multiplication operation (using *) is reserved for element-wise multiplication (in MATLAB .*). For matrix multiplication, use A @ B or (mtimes(A,B) in Python 3.4+):

print(y*y, y@y)

[[sq(y_0), sq(y_2)], 
 [sq(y_1), sq(y_3)]] 
[[(sq(y_0)+(y_2*y_1)), ((y_0*y_2)+(y_2*y_3))], 
 [((y_1*y_0)+(y_3*y_1)), ((y_1*y_2)+sq(y_3))]]
disp(y.*y)
disp(y*y)

[[sq(y_0), sq(y_2)], 
 [sq(y_1), sq(y_3)]]

[[(sq(y_0)+(y_2*y_1)), ((y_0*y_2)+(y_2*y_3))], 
 [((y_1*y_0)+(y_3*y_1)), ((y_1*y_2)+sq(y_3))]]

As is customary in MATLAB, multiplication using * and .* are equivalent when either of the arguments is a scalar.

Transposes are formed using the syntax A.T in Python, A.T() in C++ and with A or A.' in MATLAB:

print(y)
print(y.T)

[[y_0, y_2], 
 [y_1, y_3]]

[[y_0, y_1], 
 [y_2, y_3]]
disp(y)
disp(y')

[[y_0, y_2], 
 [y_1, y_3]]

[[y_0, y_1], 
 [y_2, y_3]]

Reshaping means changing the number of rows and columns but retaining the number of elements and the relative location of the nonzeros. This is a computationally very cheap operation which is performed using the syntax:

x = SX.eye(4)
print(reshape(x,2,8))
@1=1, 
[[@1, 00, 00, 00, 00, @1, 00, 00], 
 [00, 00, @1, 00, 00, 00, 00, @1]]
x = SX.eye(4);
reshape(x,2,8)
ans =

@1=1, 
[[@1, 00, 00, 00, 00, @1, 00, 00], 
 [00, 00, @1, 00, 00, 00, 00, @1]]

Concatenation means stacking matrices horizontally or vertically. Due to the column-major way of storing elements in CasADi, it is most efficient to stack matrices horizontally. Matrices that are in fact column vectors (i.e. consisting of a single column), can also be stacked efficiently vertically. Vertical and horizontal concatenation is performed using the functions vertcat and horzcat (that take a variable amount of input arguments) in Python and C++ and with square brackets in MATLAB:

x = SX.sym('x',5)
y = SX.sym('y',5)
print(vertcat(x,y))
[x_0, x_1, x_2, x_3, x_4, y_0, y_1, y_2, y_3, y_4]
x = SX.sym('x',5);
y = SX.sym('y',5);
disp([x;y])
[x_0, x_1, x_2, x_3, x_4, y_0, y_1, y_2, y_3, y_4]
print(horzcat(x,y))

[[x_0, y_0], 
 [x_1, y_1], 
 [x_2, y_2], 
 [x_3, y_3], 
 [x_4, y_4]]
disp([x,y])

[[x_0, y_0], 
 [x_1, y_1], 
 [x_2, y_2], 
 [x_3, y_3], 
 [x_4, y_4]]

There are also variants of these functions that take a list (in Python) or a cell array (in Matlab) as inputs:

L = [x,y]
print(hcat(L))

[[x_0, y_0], 
 [x_1, y_1], 
 [x_2, y_2], 
 [x_3, y_3], 
 [x_4, y_4]]
L = {x,y};
disp([L{:}])

[[x_0, y_0], 
 [x_1, y_1], 
 [x_2, y_2], 
 [x_3, y_3], 
 [x_4, y_4]]

Horizontal and vertical split are the inverse operations of the above introduced horizontal and vertical concatenation. To split up an expression horizontally into \(n\) smaller expressions, you need to provide, in addition to the expression being split, a vector offset of length \(n+1\). The first element of the offset vector must be 0 and the last element must be the number of columns. Remaining elements must follow in a non-decreasing order. The output \(i\) of the split operation then contains the columns \(c\) with \(\textit{offset}[i] \le c < \textit{offset}[i+1]\). The following demonstrates the syntax:

x = SX.sym('x',5,2)
w = horzsplit(x,[0,1,2])
print(w[0], w[1])
[x_0, x_1, x_2, x_3, x_4] [x_5, x_6, x_7, x_8, x_9]
x = SX.sym('x',5,2);
w = horzsplit(x,[0,1,2]);
disp(w(1)), disp(w(2))

{
  [1,1] =

  [x_0, x_1, x_2, x_3, x_4]

}

{
  [1,1] =

  [x_5, x_6, x_7, x_8, x_9]

}

The vertsplit operation works analogously, but with the offset vector referring to rows:

w = vertsplit(x,[0,3,5])
print(w[0], w[1])

[[x_0, x_5], 
 [x_1, x_6], 
 [x_2, x_7]] 
[[x_3, x_8], 
 [x_4, x_9]]
w = vertsplit(x,[0,3,5]);
disp(w{1}), disp(w{2})

[[x_0, x_5], 
 [x_1, x_6], 
 [x_2, x_7]]

[[x_3, x_8], 
 [x_4, x_9]]

Note that it is always possible to use slice element access instead of horizontal and vertical split, for the above vertical split:

w = [x[0:3,:], x[3:5,:]]
print(w[0], w[1])

[[x_0, x_5], 
 [x_1, x_6], 
 [x_2, x_7]] 
[[x_3, x_8], 
 [x_4, x_9]]
w = {x(1:3,:), x(4:5,:)};
disp(w{1}), disp(w{2})

[[x_0, x_5], 
 [x_1, x_6], 
 [x_2, x_7]]

[[x_3, x_8], 
 [x_4, x_9]]

For SX graphs, this alternative way is completely equivalent, but for MX graphs using horzsplit/vertsplit is significantly more efficient when all the split expressions are needed.

Inner product, defined as \(<A,B> := \text{tr}(A \, B) = \sum_{i,j} \, A_{i,j} \, B_{i,j}\) are created as follows:

x = SX.sym('x',2,2)
print(dot(x,x))
(((sq(x_0)+sq(x_1))+sq(x_2))+sq(x_3))
x = SX.sym('x',2,2);
disp(dot(x,x))
(((sq(x_0)+sq(x_1))+sq(x_2))+sq(x_3))

Many of the above operations are also defined for the Sparsity class (Section 3.5), e.g. vertcat, horzsplit, transposing, addition (which returns the union of two sparsity patterns) and multiplication (which returns the intersection of two sparsity patterns).

3.7. Querying properties

You can check if a matrix or sparsity pattern has a certain property by calling an appropriate member function. e.g.

y = SX.sym('y',10,1)
print(y.shape)
(10, 1)
y = SX.sym('y',10,1);
size(y)
ans =

   10    1

Note that in MATLAB, obj.myfcn(arg) and myfcn(obj, arg) are both valid ways of calling a member function myfcn. The latter variant is probably preferable from a style viewpoint.

Some commonly used properties for a matrix A are:

  • A.size1() The number of rows
  • A.size2() The number of columns
  • A.shape (in MATLAB “size”) The shape, i.e. the pair (nrow,*ncol*)
  • A.numel() The number of elements, i.e \(\textit{nrow} * \textit{ncol}\)
  • A.nnz() The number of structurally nonzero elements, equal to A.numel() if dense.
  • A.sparsity() Retrieve a reference to the sparsity pattern
  • A.is_dense() Is a matrix dense, i.e. having no structural zeros
  • A.is_scalar() Is the matrix a scalar, i.e. having dimensions 1-by-1?
  • A.is_column() Is the matrix a vector, i.e. having dimensions \(n\)-by-1?
  • A.is_square() Is the matrix square?
  • A.is_triu() Is the matrix upper triangular?
  • A.is_constant() Are the matrix entries all constant?
  • A.is_integer() Are the matrix entries all integer-valued?

The last queries are examples of queries for which false negative returns are allowed. A matrix for which A.is_constant() is true is guaranteed to be constant, but is not guaranteed to be non-constant if A.is_constant() is false. We recommend you to check the API documentation for a particular function before using it for the first time.

3.8. Linear algebra

CasADi supports a limited number of linear algebra operations, e.g. for solution of linear systems of equations:

A = MX.sym('A',3,3)
b = MX.sym('b',3)
print(solve(A,b))
(A\b)
A = MX.sym('A',3,3);
b = MX.sym('b',3);
disp(A\b)
(A\b)

3.9. Calculus – algorithmic differentiation

The single most central functionality of CasADi is algorithmic (or automatic) differentiation (AD). For a function \(f: \mathbb{R}^N \rightarrow \mathbb{R}^M\):

\[y = f(x),\]

forward mode directional derivatives can be used to calculate Jacobian-times-vector products:

\[\hat{y} = \frac{\partial f}{\partial x} \, \hat{x}.\]

Similarly, reverse mode directional derivatives can be used to calculate Jacobian-transposed-times-vector products:

\[\bar{x} = \left(\frac{\partial f}{\partial x}\right)^{\text{T}} \, \bar{y}.\]

Both forward and reverse mode directional derivatives are calculated at a cost proportional to evaluating \(f(x)\), regardless of the dimension of \(x\).

CasADi is also able to generate complete, sparse Jacobians efficiently. The algorithm for this is very complex, but essentially consists of the following steps:

  • Automatically detect the sparsity pattern of the Jacobian
  • Use graph coloring techniques to find a few forward and/or directional derivatives needed to construct the complete Jacobian
  • Calculate the directional derivatives numerically or symbolically
  • Assemble the complete Jacobian

Hessians are calculated by first calculating the gradient and then performing the same steps as above to calculate the Jacobian of the gradient in the same way as above, while exploiting symmetry.

3.9.1. Syntax

An expression for a Jacobian is obtained using the syntax:

A = SX.sym('A',3,2)
x = SX.sym('x',2)
print(jacobian(A@x,x))

[[A_0, A_3], 
 [A_1, A_4], 
 [A_2, A_5]]
A = SX.sym('A',3,2);
x = SX.sym('x',2);
jacobian(A*x,x)
ans =


[[A_0, A_3], 
 [A_1, A_4], 
 [A_2, A_5]]

When the differentiated expression is a scalar, you can also calculate the gradient in the matrix sense:

print(gradient(dot(A,A),A))

[[(A_0+A_0), (A_3+A_3)], 
 [(A_1+A_1), (A_4+A_4)], 
 [(A_2+A_2), (A_5+A_5)]]
gradient(dot(A,A),A)
ans =


[[(A_0+A_0), (A_3+A_3)], 
 [(A_1+A_1), (A_4+A_4)], 
 [(A_2+A_2), (A_5+A_5)]]

Note that, unlike jacobian, gradient always returns a dense vector.

Hessians, and as a by-product gradients, are obtained as follows:

[H,g] = hessian(dot(x,x),x)
print('H:', H)
H: @1=2, 
[[@1, 00], 
 [00, @1]]
[H,g] = hessian(dot(x,x),x);
display(H)
H =

@1=2, 
[[@1, 00], 
 [00, @1]]

For calculating a Jacobian-times-vector product, the jtimes function – performing forward mode AD – is often more efficient than creating the full Jacobian and performing a matrix-vector multiplication:

A = DM([[1,3],[4,7],[2,8]])
x = SX.sym('x',2)
v = SX.sym('v',2)
f = mtimes(A,x)
print(jtimes(f,x,v))
[(v_0+(3*v_1)), ((4*v_0)+(7*v_1)), ((2*v_0)+(8*v_1))]
A = [1 3;4 7;2 8];
x = SX.sym('x',2);
v = SX.sym('v',2);
f = A*x;
jtimes(f,x,v)
ans =

[(v_0+(3*v_1)), ((4*v_0)+(7*v_1)), ((2*v_0)+(8*v_1))]

The jtimes function optionally calculates the transposed-Jacobian-times-vector product, i.e. reverse mode AD:

w = SX.sym('w',3)
f = mtimes(A,x)
print(jtimes(f,x,w,True))
[(((2*w_2)+(4*w_1))+w_0), (((8*w_2)+(7*w_1))+(3*w_0))]
w = SX.sym('w',3);
f = A*x;
jtimes(f,x,w,true)
ans =

[(((2*w_2)+(4*w_1))+w_0), (((8*w_2)+(7*w_1))+(3*w_0))]

4. Function objects

CasADi allows the user to create function objects, in C++ terminology often referred to as functors. This includes functions that are defined by a symbolic expression, ODE/DAE integrators, QP solvers, NLP solvers etc.

Function objects are typically created with the syntax:

f = functionname(name, arguments, ..., [options])

The name is mainly a display name that will show up in e.g. error messages or as comments in generated C code. This is followed by a set of arguments, which is class dependent. Finally, the user can pass an options structure for customizing the behavior of the class. The options structure is a dictionary type in Python, a struct in MATLAB or CasADi’s Dict type in C++.

A Function can be constructed by passing a list of input expressions and a list of output expressions:

x = SX.sym('x',2)
y = SX.sym('y')
f = Function('f',[x,y],\
           [x,sin(y)*x])
print(f)
f:(i0[2],i1)->(o0[2],o1[2]) SXFunction
x = SX.sym('x',2);
y = SX.sym('y');
f = Function('f',{x,y},...
           {x,sin(y)*x});
disp(f)
f:(i0[2],i1)->(o0[2],o1[2]) SXFunction

which defines a function \(f : \mathbb{R}^{2} \times \mathbb{R} \rightarrow \mathbb{R}^{2} \times \mathbb{R}^{2}, \quad (x,y) \mapsto (x,\sin(y) x)\). Note that all function objects in CasADi, including the above, are multiple matrix-valued input, multiple, matrix-valued output.

MX expression graphs work the same way:

x = MX.sym('x',2)
y = MX.sym('y')
f = Function('f',[x,y],\
           [x,sin(y)*x])
print(f)
f:(i0[2],i1)->(o0[2],o1[2]) MXFunction
x = MX.sym('x',2);
y = MX.sym('y');
f = Function('f',{x,y},...
           {x,sin(y)*x});
disp(f)
f:(i0[2],i1)->(o0[2],o1[2]) MXFunction

When creating a Function from expressions like that, it is always advisory to name the inputs and outputs as follows:

f = Function('f',[x,y],\
      [x,sin(y)*x],\
      ['x','y'],['r','q'])
print(f)
f:(x[2],y)->(r[2],q[2]) MXFunction
f = Function('f',{x,y},...
      {x,sin(y)*x},...
      {'x','y'},{'r','q'});
disp(f)
f:(x[2],y)->(r[2],q[2]) MXFunction

Naming inputs and outputs is preferred for a number of reasons:

  • No need to remember the number or order of arguments
  • Inputs or outputs that are absent can be left unset
  • More readable and less error prone syntax. E.g. f.jacobian('x','q') instead of f.jacobian(0,1).

For Function instances – to be encountered later – that are not created directly from expressions, the inputs and outputs are named automatically.

4.1. Calling function objects

MX expressions may contain calls to Function-derived functions. Calling a function object is both done for the numerical evaluation and, by passing symbolic arguments, for embedding a call to the function object into an expression graph (cf. also Section 4.4).

To call a function object, you either pass the argument in the correct order:

r0, q0 = f(1.1,3.3)
print('r0:',r0)
print('q0:',q0)
r0: [1.1, 1.1]
q0: [-0.17352, -0.17352]
[r0, q0] = f(1.1,3.3);
disp(r0)
disp(q0)
[1.1, 1.1]
[-0.17352, -0.17352]

or the arguments and their names as follows, which will result in a dictionary (dict in Python, struct in MATLAB and std::map<std::string, MatrixType> in C++):

res = f(x=1.1, y=3.3)
print('res:', res)
res: {'q': DM([-0.17352, -0.17352]), 'r': DM([1.1, 1.1])}
res = f('x',1.1,'y',3.3);
disp(res)

  scalar structure containing the fields:

    q =

    [-0.17352, -0.17352]

    r =

    [1.1, 1.1]

When calling a function object, the dimensions (but not necessarily the sparsity patterns) of the evaluation arguments have to match those of the function inputs, with two exceptions:

  • A row vector can be passed instead of a column vector and vice versa.
  • A scalar argument can always be passed, regardless of the input dimension. This has the meaning of setting all elements of the input matrix to that value.

When the number of inputs to a function object is large or changing, an alternative syntax to the above is to use the call function which takes a Python list / MATLAB cell array or, alternatively, a Python dict / MATLAB struct. The return value will have the same type:

arg = [1.1,3.3]
res = f.call(arg)
print('res:', res)
arg = {'x':1.1,'y':3.3}
res = f.call(arg)
print('res:', res)
res: [DM([1.1, 1.1]), DM([-0.17352, -0.17352])]
res: {'q': DM([-0.17352, -0.17352]), 'r': DM([1.1, 1.1])}
arg = {1.1,3.3};
res = f.call(arg);
disp(res)
arg = struct('x',1.1,'y',3.3);
res = f.call(arg);
disp(res)

{
  [1,1] =

  [1.1, 1.1]

  [1,2] =

  [-0.17352, -0.17352]

}

  scalar structure containing the fields:

    q =

    [-0.17352, -0.17352]

    r =

    [1.1, 1.1]

4.2. Converting MX to SX

A function object defined by an MX graph that only contains built-in operations (e.g. element-wise operations such as addition, square root, matrix multiplications and calls to SX functions, can be converted into a function defined purely by an SX graph using the syntax:

sx_function = mx_function.expand()

This might speed up the calculations significantly, but might also cause extra memory overhead.

4.3. Nonlinear root-finding problems

Consider the following system of equations:

\[\begin{split}\begin{aligned} &g_0(z, x_1, x_2, \ldots, x_n) &&= 0 \\ &g_1(z, x_1, x_2, \ldots, x_n) &&= y_1 \\ &g_2(z, x_1, x_2, \ldots, x_n) &&= y_2 \\ &\qquad \vdots \qquad &&\qquad \\ &g_m(z, x_1, x_2, \ldots, x_n) &&= y_m, \end{aligned}\end{split}\]

where the first equation uniquely defines \(z\) as a function of \(x_1\), ldots, \(x_n\) by the implicit function theorem and the remaining equations define the auxiliary outputs \(y_1\), ldots, \(y_m\).

Given a function \(g\) for evaluating \(g_0\), ldots, \(g_m\), we can use CasADi to automatically formulate a function \(G: \{z_{\text{guess}}, x_1, x_2, \ldots, x_n\} \rightarrow \{z, y_1, y_2, \ldots, y_m\}\). This function includes a guess for \(z\) to handle the case when the solution is non-unique. The syntax for this, assuming \(n=m=1\) for simplicity, is:

z = SX.sym('x',nz)
x = SX.sym('x',nx)
g0 = sin(x+z)
g1 = cos(x-z)
g = Function('g',[z,x],[g0,g1])
G = rootfinder('G','newton',g)
print(G)
G:(i0,i1)->(o0,o1) Newton
z = SX.sym('x',nz);
x = SX.sym('x',nx);
g0 = sin(x+z);
g1 = cos(x-z);
g = Function('g',{z,x},{g0,g1});
G = rootfinder('G','newton',g);
disp(G)
G:(i0,i1)->(o0,o1) Newton

where the rootfinder function expects a display name, the name of a solver plugin (here a simple full-step Newton method) and the residual function.

Rootfinding objects in CasADi are differential objects and derivatives can be calculated exactly to arbitrary order.

4.4. Initial-value problems and sensitivity analysis

CasADi can be used to solve initial-value problems in ODE or DAE. The problem formulation used is a DAE of semi-explicit form with quadratures:

\[\begin{split}\begin{aligned} \dot{x} &= f_{\text{ode}}(t,x,z,p), \qquad x(0) = x_0 \\ 0 &= f_{\text{alg}}(t,x,z,p) \\ \dot{q} &= f_{\text{quad}}(t,x,z,p), \qquad q(0) = 0 \end{aligned}\end{split}\]

For solvers of ordinary differential equations, the second equation and the algebraic variables \(z\) must be absent.

An integrator in CasADi is a function that takes the state at the initial time x0, a set of parameters p, and a guess for the algebraic variables (only for DAEs) z0 and returns the state vector xf, algebraic variables zf and the quadrature state qf, all at the final time.

The freely available SUNDIALS suite (distributed along with CasADi) contains the two popular integrators CVodes and IDAS for ODEs and DAEs respectively. These integrators have support for forward and adjoint sensitivity analysis and when used via CasADi’s Sundials interface, CasADi will automatically formulate the Jacobian information, which is needed by the backward differentiation formula (BDF) that CVodes and IDAS use. Also automatically formulated will be the forward and adjoint sensitivity equations.

4.4.1. Creating integrators

Integrators are created using CasADi’s integrator function. Different integrators schemes and interfaces are implemented as plugins, essentially shared libraries that are loaded at runtime.

Consider for example the DAE:

\[\begin{split}\begin{aligned} \dot{x} &= z+p, \\ 0 &= z \, \cos(z)-x \end{aligned}\end{split}\]

An integrator, using the “idas” plugin, can be created using the syntax:

x = SX.sym('x'); z = SX.sym('z'); p = SX.sym('p')
dae = {'x':x, 'z':z, 'p':p, 'ode':z+p, 'alg':z*cos(z)-x}
F = integrator('F', 'idas', dae)
print(F)
F:(x0,p,z0,rx0[0],rp[0],rz0[0])->(xf,qf[0],zf,rxf[0],rqf[0],rzf[0]) IdasInterface
x = SX.sym('x'); z = SX.sym('z'); p = SX.sym('p');
dae = struct('x',x,'z',z,'p',p,'ode',z+p,'alg',z*cos(z)-x);
F = integrator('F', 'idas', dae);
disp(F)
F:(x0,p,z0,rx0[0],rp[0],rz0[0])->(xf,qf[0],zf,rxf[0],rqf[0],rzf[0]) IdasInterface

Integrating this DAE from 0 to 1 with \(x(0)=0\), \(p=0.1\) and using the guess \(z(0)=0\), can be done by evaluating the created function object:

r = F(x0=0, z0=0, p=0.1)
print(r['xf'])
0.1724
r = F('x0',0,'z0',0,'p',0.1);
disp(r.xf)
0.1724

The time horizon is assumed to be fixed [1] and can be changed from its default [0, 1] by setting the options “t0” and “tf”.

4.4.2. Sensitivity analysis

From a usage point of view, an integrator behaves just like the function objects created from expressions earlier in the chapter. You can use member functions in the Function class to generate new function objects corresponding to directional derivatives (forward or reverse mode) or complete Jacobians. Then evaluate these function objects numerically to obtain sensitivity information. The documented example “sensitivity_analysis” (available in CasADi’s example collection for Python, MATLAB and C++) demonstrate how CasADi can be used to calculate first and second order derivative information (forward-over-forward, forward-over-adjoint, adjoint-over-adjoint) for a simple DAE.

4.5. Nonlinear programming

Note

This section assumes familiarity with much of the material that comes above. There is also a higher-level interface in Chapter 9. That interface can be learned stand-alone.

The NLP solvers distributed with or interfaced to CasADi solves parametric NLPs of the following form:

(4.5.1)\[\begin{split}\begin{array}{cc} \begin{array}{c} \text{minimize:} \\ x \end{array} & f(x,p) \\ \begin{array}{c} \text{subject to:} \end{array} & \begin{array}{rcl} x_{\textrm{lb}} \le & x & \le x_{\textrm{ub}} \\ g_{\textrm{lb}} \le &g(x,p)& \le g_{\textrm{ub}} \end{array} \end{array}\end{split}\]

where \(x \in \mathbb{R}^{nx}\) is the decision variable and \(p \in \mathbb{R}^{np}\) is a known parameter vector.

An NLP solver in CasADi is a function that takes the parameter value (p), the bounds (lbx, ubx, lbg, ubg) and a guess for the primal-dual solution (x0, lam_x0, lam_g0) and returns the optimal solution. Unlike integrator objects, NLP solver functions are currently not differentiable functions in CasADi.

There are several NLP solvers interfaced with CasADi. The most popular one is IPOPT, an open-source primal-dual interior point method which is included in CasADi installations. Others, that require the installation of third-party software, include SNOPT, WORHP and KNITRO. Whatever the NLP solver used, the interface will automatically generate the information that it needs to solve the NLP, which may be solver and option dependent. Typically an NLP solver will need a function that gives the Jacobian of the constraint function and a Hessian of the Lagrangian function (\(L(x,\lambda) = f(x) + \lambda^{\text{T}} \, g(x))\) with respect to \(x\).

4.5.1. Creating NLP solvers

NLP solvers are created using CasADi’s nlpsol function. Different solvers and interfaces are implemented as plugins. Consider the following form of the so-called Rosenbrock problem:

\[\begin{split}\begin{array}{cc} \begin{array}{c} \text{minimize:} \\ x,y,z \end{array} & x^2 + 100 \, z^2 \\ \begin{array}{c} \text{subject to:} \end{array} & z+(1-x)^2-y = 0 \end{array}\end{split}\]

A solver for this problem, using the “ipopt” plugin, can be created using the syntax:

x = SX.sym('x'); y = SX.sym('y'); z = SX.sym('z')
nlp = {'x':vertcat(x,y,z), 'f':x**2+100*z**2, 'g':z+(1-x)**2-y}
S = nlpsol('S', 'ipopt', nlp)
print(S)
S:(x0[3],p[],lbx[3],ubx[3],lbg,ubg,lam_x0[3],lam_g0)->(x[3],f,g,lam_x[3],lam_g,lam_p[]) IpoptInterface
x = SX.sym('x'); y = SX.sym('y'); z = SX.sym('z');
nlp = struct('x',[x;y;z], 'f',x^2+100*z^2, 'g',z+(1-x)^2-y);
S = nlpsol('S', 'ipopt', nlp);
disp(S)
S:(x0[3],p[],lbx[3],ubx[3],lbg,ubg,lam_x0[3],lam_g0)->(x[3],f,g,lam_x[3],lam_g,lam_p[]) IpoptInterface

Once the solver has been created, we can solve the NLP, using [2.5,3.0,0.75] as an initial guess, by evaluating the function S:

r = S(x0=[2.5,3.0,0.75],\
      lbg=0, ubg=0)
x_opt = r['x']
print('x_opt: ', x_opt)

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.3, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:        3
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        2

Total number of variables............................:        3
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        1
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  6.2500000e+01 0.00e+00 9.00e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.8457621e+01 1.48e-02 4.10e+01  -1.0 4.10e-01   2.0 1.00e+00 1.00e+00f  1
   2  7.8031530e+00 3.84e-03 8.76e+00  -1.0 2.63e-01   1.5 1.00e+00 1.00e+00f  1
   3  7.1678278e+00 9.42e-05 1.04e+00  -1.0 9.32e-02   1.0 1.00e+00 1.00e+00f  1
   4  6.7419924e+00 6.18e-03 9.95e-01  -1.0 2.69e-01   0.6 1.00e+00 1.00e+00f  1
   5  5.4363330e+00 7.03e-02 1.04e+00  -1.7 8.40e-01   0.1 1.00e+00 1.00e+00f  1
   6  1.2144815e+00 1.52e+00 1.32e+00  -1.7 3.21e+00  -0.4 1.00e+00 1.00e+00f  1
   7  1.0214718e+00 2.51e-01 1.17e+01  -1.7 1.33e+00   0.9 1.00e+00 1.00e+00h  1
   8  3.1864085e-01 1.04e-03 7.53e-01  -1.7 3.58e-01    -  1.00e+00 1.00e+00f  1
   9  0.0000000e+00 3.19e-01 0.00e+00  -1.7 5.64e-01    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  0.0000000e+00 0.00e+00 0.00e+00  -1.7 3.19e-01    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 10

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   0.0000000000000000e+00    0.0000000000000000e+00


Number of objective function evaluations             = 11
Number of objective gradient evaluations             = 11
Number of equality constraint evaluations            = 11
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 11
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 10
Total CPU secs in IPOPT (w/o function evaluations)   =      0.003
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
               t_proc [s]   t_wall [s]    n_eval
           S      0.00467      0.00467         1
       nlp_f      1.2e-05     1.19e-05        11
       nlp_g      2.7e-05     2.43e-05        11
  nlp_grad_f        2e-05     1.86e-05        12
  nlp_hess_l      1.2e-05     1.39e-05        10
   nlp_jac_g      1.4e-05     1.39e-05        12
x_opt:  [0, 1, 0]
r = S('x0',[2.5,3.0,0.75],...
      'lbg',0,'ubg',0);
x_opt = r.x;
disp(x_opt)

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.3, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:        3
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        2

Total number of variables............................:        3
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        1
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  6.2500000e+01 0.00e+00 9.00e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.8457621e+01 1.48e-02 4.10e+01  -1.0 4.10e-01   2.0 1.00e+00 1.00e+00f  1
   2  7.8031530e+00 3.84e-03 8.76e+00  -1.0 2.63e-01   1.5 1.00e+00 1.00e+00f  1
   3  7.1678278e+00 9.42e-05 1.04e+00  -1.0 9.32e-02   1.0 1.00e+00 1.00e+00f  1
   4  6.7419924e+00 6.18e-03 9.95e-01  -1.0 2.69e-01   0.6 1.00e+00 1.00e+00f  1
   5  5.4363330e+00 7.03e-02 1.04e+00  -1.7 8.40e-01   0.1 1.00e+00 1.00e+00f  1
   6  1.2144815e+00 1.52e+00 1.32e+00  -1.7 3.21e+00  -0.4 1.00e+00 1.00e+00f  1
   7  1.0214718e+00 2.51e-01 1.17e+01  -1.7 1.33e+00   0.9 1.00e+00 1.00e+00h  1
   8  3.1864085e-01 1.04e-03 7.53e-01  -1.7 3.58e-01    -  1.00e+00 1.00e+00f  1
   9  0.0000000e+00 3.19e-01 0.00e+00  -1.7 5.64e-01    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  0.0000000e+00 0.00e+00 0.00e+00  -1.7 3.19e-01    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 10

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   0.0000000000000000e+00    0.0000000000000000e+00


Number of objective function evaluations             = 11
Number of objective gradient evaluations             = 11
Number of equality constraint evaluations            = 11
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 11
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 10
Total CPU secs in IPOPT (w/o function evaluations)   =      0.005
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
               t_proc [s]   t_wall [s]    n_eval
           S      0.00972      0.00975         1
       nlp_f      1.8e-05      1.8e-05        11
       nlp_g      7.5e-05      6.8e-05        11
  nlp_grad_f      3.8e-05     3.76e-05        12
  nlp_hess_l      2.7e-05      2.6e-05        10
   nlp_jac_g      2.3e-05     2.32e-05        12
[0, 1, 0]

4.6. Quadratic programming

CasADi provides interfaces to solve quadratic programs (QPs). Supported solvers are the open-source solvers qpOASES (distributed with CasADi) and OOQP as well as the commercial solvers CPLEX and GUROBI.

There are two different ways to solve QPs in CasADi, using a high-level interface and a low-level interface. They are described in the following.

4.6.1. High-level interface

The high-level interface for quadratic programming mirrors that of nonlinear programming, i.e. expects a problem of the form (4.5.1), with the restriction that objective function \(f(x,p)\) must be a convex quadratic function in \(x\) and the constraint function \(g(x,p)\) must be linear in \(x\). If the functions are not quadratic and linear, respectively, the solution is done at the current linearization point, given by the “initial guess” for \(x\).

If the objective function is not convex, the solver may or may not fail to find a solution or the solution may not be unique.

To illustrate the syntax, we consider the following convex QP:

(4.6.1)\[\begin{split} \begin{array}{cc} \begin{array}{c} \text{minimize:} \\ x,y \end{array} & x^2 + y^2 \\ \begin{array}{c} \text{subject to:} \end{array} & x+y-10 \ge 0 \end{array}\end{split}\]

To solve this problem with the high-level interface, we simply replace nlpsol with qpsol and use a QP solver plugin such as the with CasADi distributed qpOASES:

x = SX.sym('x'); y = SX.sym('y')
qp = {'x':vertcat(x,y), 'f':x**2+y**2, 'g':x+y-10}
S = qpsol('S', 'qpoases', qp)
print(S)
S:(x0[2],p[],lbx[2],ubx[2],lbg,ubg,lam_x0[2],lam_g0)->(x[2],f,g,lam_x[2],lam_g,lam_p[]) MXFunction
x = SX.sym('x'); y = SX.sym('y');
qp = struct('x',[x;y], 'f',x^2+y^2, 'g',x+y-10);
S = qpsol('S', 'qpoases', qp);
disp(S)
S:(x0[2],p[],lbx[2],ubx[2],lbg,ubg,lam_x0[2],lam_g0)->(x[2],f,g,lam_x[2],lam_g,lam_p[]) MXFunction

The created solver object S will have the same input and output signature as the solver objects created with nlpsol. Since the solution is unique, it is less important to provide an initial guess:

r = S(lbg=0)
x_opt = r['x']
print('x_opt: ', x_opt)


####################   qpOASES  --  QP NO.   1   #####################

    Iter   |    StepLength    |       Info       |   nFX   |   nAC    
 ----------+------------------+------------------+---------+--------- 
       0   |   1.866661e-07   |   ADD CON    0   |     1   |     1   
       1   |   8.333622e-10   |   REM BND    1   |     0   |     1   
       2   |   1.000000e+00   |    QP SOLVED     |     0   |     1   
x_opt:  [5, 5]
r = S('lbg',0);
x_opt = r.x;
disp(x_opt)


####################   qpOASES  --  QP NO.   1   #####################

    Iter   |    StepLength    |       Info       |   nFX   |   nAC    
 ----------+------------------+------------------+---------+--------- 
       0   |   1.866661e-07   |   ADD CON    0   |     1   |     1   
       1   |   8.333622e-10   |   REM BND    1   |     0   |     1   
       2   |   1.000000e+00   |    QP SOLVED     |     0   |     1   
[5, 5]

4.6.2. Low-level interface

The low-level interface, on the other hand, solves QPs of the following form:

(4.6.2)\[\begin{split} \begin{array}{cc} \begin{array}{c} \text{minimize:} \\ x \end{array} & \frac{1}{2} x^T \, H \, x + g^T \, x \\ \begin{array}{c} \text{subject to:} \end{array} & \begin{array}{rcl} x_{\textrm{lb}} \le & x & \le x_{\textrm{ub}} \\ a_{\textrm{lb}} \le & A \, x& \le a_{\textrm{ub}} \end{array} \end{array}\end{split}\]

Encoding problem (4.6.1) in this form, omitting bounds that are infinite, is straightforward:

H = 2*DM.eye(2)
A = DM.ones(1,2)
g = DM.zeros(2)
lba = 10.
H = 2*DM.eye(2);
A = DM.ones(1,2);
g = DM.zeros(2);
lba = 10;

To create a solver instance, instead of passing symbolic expressions for the QP, we now pass the sparsity patterns of the matrices \(H\) and \(A\). Since we used CasADi’s DM-type above, we can simply query the sparsity patterns:

qp = {}
qp['h'] = H.sparsity()
qp['a'] = A.sparsity()
S = conic('S','qpoases',qp)
print(S)
S:(h[2x2,2nz],g[2],a[1x2],lba,uba,lbx[2],ubx[2],x0[2],lam_x0[2],lam_a0,q[],p[])->(x[2],cost,lam_a,lam_x[2]) QpoasesInterface
qp = struct;
qp.h = H.sparsity();
qp.a = A.sparsity();
S = conic('S','qpoases',qp);
disp(S)
S:(h[2x2,2nz],g[2],a[1x2],lba,uba,lbx[2],ubx[2],x0[2],lam_x0[2],lam_a0,q[],p[])->(x[2],cost,lam_a,lam_x[2]) QpoasesInterface

The returned Function instance will have a different input/output signature compared to the high-level interface, one that includes the matrices \(H\) and \(A\):

r = S(h=H, g=g, \
      a=A, lba=lba)
x_opt = r['x']
print('x_opt: ', x_opt)


####################   qpOASES  --  QP NO.   1   #####################

    Iter   |    StepLength    |       Info       |   nFX   |   nAC    
 ----------+------------------+------------------+---------+--------- 
       0   |   1.866661e-07   |   ADD CON    0   |     1   |     1   
       1   |   8.333622e-10   |   REM BND    1   |     0   |     1   
       2   |   1.000000e+00   |    QP SOLVED     |     0   |     1   
x_opt:  [5, 5]
r = S('h', H, 'g', g,...
      'a', A, 'lba', lba);
x_opt = r.x;
disp(x_opt)


####################   qpOASES  --  QP NO.   1   #####################

    Iter   |    StepLength    |       Info       |   nFX   |   nAC    
 ----------+------------------+------------------+---------+--------- 
       0   |   1.866661e-07   |   ADD CON    0   |     1   |     1   
       1   |   8.333622e-10   |   REM BND    1   |     0   |     1   
       2   |   1.000000e+00   |    QP SOLVED     |     0   |     1   
[5, 5]

4.7. For-loop equivalents

When modeling using expression graphs in CasADi, it is a common pattern to use of for-loop constructs of the host language (C++/Python/Matlab).

The graph size will grow linearly with the loop size \(n\), and so will the construction time of the expression graph and the initialization time of functions using that expression.

We offer some special constructs that improve on this complexity.

4.7.1. Map

Suppose you are interested in computing a function \(f : \mathbb{R}^{n} \rightarrow \mathbb{R}^{m}\) repeatedly on all columns of a matrix \(X \in \mathbb{R}^{n \times N}\), and aggregating all results in a result matrix \(Y \in \mathbb{R}^{m \times N}\):

N = 4
X = MX.sym("X",1,N)

print(f)

ys = []
for i in range(N):
  ys.append(f(X[:,i]))

Y = hcat(ys)
F = Function('F',[X],[Y])
print(F)
f:(i0)->(o0) SXFunction
F:(i0[1x4])->(o0[1x4]) MXFunction
N = 4;
X = MX.sym('X',1,N);

disp(f)

ys = {};
for i=1:N
  ys{end+1} = f(X(:,i));
end
Y = [ys{:}];
F = Function('F',{X},{Y});
disp(F)
f:(i0)->(o0) SXFunction
F:(i0[1x4])->(o0[1x4]) MXFunction

The aggregate function \(F : \mathbb{R}^{n \times N} \rightarrow \mathbb{R}^{m \times N}\) can be obtained with the map construct:

F = f.map(N)
print(F)
map4_f:(i0[1x4])->(o0[1x4]) Map
F = f.map(N);
disp(F)
map4_f:(i0[1x4])->(o0[1x4]) Map

CasADi can be instructed to parallelize when \(F\) gets evaluated. In the following example, we dedicate 2 threads for the map task.

F = f.map(N,"thread",2)
print(F)
threadmap2_map2_f:(i0[1x4])->(o0[1x4]) ThreadMap
F = f.map(N,'thread',2);
disp(F)
threadmap2_map2_f:(i0[1x4])->(o0[1x4]) ThreadMap

The map operation supports primitive functions \(f\) with multiple inputs/outputs which can also be matrices. Aggregation will always happen horizontally.

The map operation exhibits constant graph size and initialization time.

4.7.2. Fold

In case each for-loop iteration depends on the result from the previous iteration, the fold construct applies. In the following, the x variable acts as an accumulater that is initialized by \(x_0 \in \mathbb{R}^{n}\):

x = x0
for i in range(N):
  x = f(x)

F = Function('F',[x0],[x])
print(F)
F:(i0)->(o0) MXFunction
x = x0;
for i=1:N
  x = f(x);
end

F = Function('F',{x0},{x});
disp(F)
F:(i0)->(o0) MXFunction

For a given function \(f : \mathbb{R}^{n} \rightarrow \mathbb{R}^{n}\), the result function \(F : \mathbb{R}^{n} \rightarrow \mathbb{R}^{n}\) can be obtained with the fold construct:

F = f.fold(N)
print(F)
fold_f:(i0)->(o0) MXFunction
F = f.fold(N);
disp(F)
fold_f:(i0)->(o0) MXFunction

In case intermediate accumulator values are desired as output (\(\mathbb{R}^{n} \rightarrow \mathbb{R}^{n \times N}\)), use mapaccum instead of fold.

The fold/mapaccum operation supports primitive functions \(f\) with multiple inputs/outputs which can also be matrices. The first input and output are used for accumulating, while the remainder inputs are read column-wise over the iterations.

The map/mapaccum operation exhibits a graph size and initialization time that scales logarithmically with \(n\).

Footnotes

[1]for problems with free end time, you can always scale time by introducing an extra parameter and substitute \(t\) for a dimensionless time variable that goes from 0 to 1

5. Generating C-code

The numerical evaluation of function objects in CasADi normally takes place in virtual machines, implemented as part of CasADi’s symbolic framework. But CasADi also supports the generation of self-contained C-code for a large subset of function objects.

C-code generation is interesting for a number of reasons:

  • Speeding up the evaluation time. As a rule of thumb, the numerical evaluation of autogenerated code, compiled with code optimization flags, can be between 4 and 10 times faster than the same code executed in CasADi’s virtual machines.
  • Allowing code to be compiled on a system where CasADi is not installed, such as an embedded system. All that is needed to compile the generated code is a C compiler.
  • Debugging and profiling functions. The generated code is essentially a mirror of the evaluation that takes place in the virtual machines and if a particular operation is slow, this is likely to show up when analyzing the generated code with a profiling tool such as gprof. By looking at the code, it is also possible to detect what is potentially done in a suboptimal way. If the code is very long and takes a long time to compile, it is an indication that some functions need to be broken up in smaller, but nested functions.

5.1. Syntax for generating code

Generated C code can be as simple as calling the generate member function for a Function instance.

x = MX.sym('x',2)
y = MX.sym('y')
f = Function('f',[x,y],\
      [x,sin(y)*x],\
      ['x','y'],['r','q'])
f.generate('gen.c')
print(open('gen.c','r').read())
/* This file was automatically generated by CasADi.
   The CasADi copyright holders make no ownership claim of its contents. */
#ifdef __cplusplus
extern "C" {
#endif

/* How to prefix internal symbols */
#ifdef CODEGEN_PREFIX
  #define NAMESPACE_CONCAT(NS, ID) _NAMESPACE_CONCAT(NS, ID)
  #define _NAMESPACE_CONCAT(NS, ID) NS ## ID
  #define CASADI_PREFIX(ID) NAMESPACE_CONCAT(CODEGEN_PREFIX, ID)
#else
  #define CASADI_PREFIX(ID) gen_ ## ID
#endif

#include <math.h>

#ifndef casadi_real
#define casadi_real double
#endif

#ifndef casadi_int
#define casadi_int long long int
#endif

/* Add prefix to internal symbols */
#define casadi_copy CASADI_PREFIX(copy)
#define casadi_f0 CASADI_PREFIX(f0)
#define casadi_s0 CASADI_PREFIX(s0)
#define casadi_s1 CASADI_PREFIX(s1)

/* Symbol visibility in DLLs */
#ifndef CASADI_SYMBOL_EXPORT
  #if defined(_WIN32) || defined(__WIN32__) || defined(__CYGWIN__)
    #if defined(STATIC_LINKED)
      #define CASADI_SYMBOL_EXPORT
    #else
      #define CASADI_SYMBOL_EXPORT __declspec(dllexport)
    #endif
  #elif defined(__GNUC__) && defined(GCC_HASCLASSVISIBILITY)
    #define CASADI_SYMBOL_EXPORT __attribute__ ((visibility ("default")))
  #else
    #define CASADI_SYMBOL_EXPORT
  #endif
#endif

static const casadi_int casadi_s0[6] = {2, 1, 0, 2, 0, 1};
static const casadi_int casadi_s1[5] = {1, 1, 0, 1, 0};

void casadi_copy(const casadi_real* x, casadi_int n, casadi_real* y) {
  casadi_int i;
  if (y) {
    if (x) {
      for (i=0; i<n; ++i) *y++ = *x++;
    } else {
      for (i=0; i<n; ++i) *y++ = 0.;
    }
  }
}

/* f:(x[2],y)->(r[2],q[2]) */
static int casadi_f0(const casadi_real** arg, casadi_real** res, casadi_int* iw, casadi_real* w, void* mem) {
  casadi_int i;
  casadi_real *rr;
  const casadi_real *cs;
  casadi_real *w0=w+0, w1;
  /* #0: @0 = input[0][0] */
  casadi_copy(arg[0], 2, w0);
  /* #1: output[0][0] = @0 */
  casadi_copy(w0, 2, res[0]);
  /* #2: @1 = input[1][0] */
  w1 = arg[1] ? arg[1][0] : 0;
  /* #3: @1 = sin(@1) */
  w1 = sin( w1 );
  /* #4: @0 = (@1*@0) */
  for (i=0, rr=w0, cs=w0; i<2; ++i) (*rr++)  = (w1*(*cs++));
  /* #5: output[1][0] = @0 */
  casadi_copy(w0, 2, res[1]);
  return 0;
}

CASADI_SYMBOL_EXPORT int f(const casadi_real** arg, casadi_real** res, casadi_int* iw, casadi_real* w, void* mem){
  return casadi_f0(arg, res, iw, w, mem);
}

CASADI_SYMBOL_EXPORT void f_incref(void) {
}

CASADI_SYMBOL_EXPORT void f_decref(void) {
}

CASADI_SYMBOL_EXPORT casadi_int f_n_in(void) { return 2;}

CASADI_SYMBOL_EXPORT casadi_int f_n_out(void) { return 2;}

CASADI_SYMBOL_EXPORT const char* f_name_in(casadi_int i){
  switch (i) {
    case 0: return "x";
    case 1: return "y";
    default: return 0;
  }
}

CASADI_SYMBOL_EXPORT const char* f_name_out(casadi_int i){
  switch (i) {
    case 0: return "r";
    case 1: return "q";
    default: return 0;
  }
}

CASADI_SYMBOL_EXPORT const casadi_int* f_sparsity_in(casadi_int i) {
  switch (i) {
    case 0: return casadi_s0;
    case 1: return casadi_s1;
    default: return 0;
  }
}

CASADI_SYMBOL_EXPORT const casadi_int* f_sparsity_out(casadi_int i) {
  switch (i) {
    case 0: return casadi_s0;
    case 1: return casadi_s0;
    default: return 0;
  }
}

CASADI_SYMBOL_EXPORT int f_work(casadi_int *sz_arg, casadi_int* sz_res, casadi_int *sz_iw, casadi_int *sz_w) {
  if (sz_arg) *sz_arg = 4;
  if (sz_res) *sz_res = 3;
  if (sz_iw) *sz_iw = 0;
  if (sz_w) *sz_w = 3;
  return 0;
}


#ifdef __cplusplus
} /* extern "C" */
#endif

x = MX.sym('x',2);
y = MX.sym('y');
f = Function('f',{x,y},...
      {x,sin(y)*x},...
      {'x','y'},{'r','q'});
f.generate('gen.c');
type('gen.c')
/* This file was automatically generated by CasADi.
   The CasADi copyright holders make no ownership claim of its contents. */
#ifdef __cplusplus
extern "C" {
#endif

/* How to prefix internal symbols */
#ifdef CODEGEN_PREFIX
  #define NAMESPACE_CONCAT(NS, ID) _NAMESPACE_CONCAT(NS, ID)
  #define _NAMESPACE_CONCAT(NS, ID) NS ## ID
  #define CASADI_PREFIX(ID) NAMESPACE_CONCAT(CODEGEN_PREFIX, ID)
#else
  #define CASADI_PREFIX(ID) gen_ ## ID
#endif

#include <math.h>

#ifndef casadi_real
#define casadi_real double
#endif

#ifndef casadi_int
#define casadi_int long long int
#endif

/* Add prefix to internal symbols */
#define casadi_copy CASADI_PREFIX(copy)
#define casadi_f0 CASADI_PREFIX(f0)
#define casadi_s0 CASADI_PREFIX(s0)
#define casadi_s1 CASADI_PREFIX(s1)

/* Symbol visibility in DLLs */
#ifndef CASADI_SYMBOL_EXPORT
  #if defined(_WIN32) || defined(__WIN32__) || defined(__CYGWIN__)
    #if defined(STATIC_LINKED)
      #define CASADI_SYMBOL_EXPORT
    #else
      #define CASADI_SYMBOL_EXPORT __declspec(dllexport)
    #endif
  #elif defined(__GNUC__) && defined(GCC_HASCLASSVISIBILITY)
    #define CASADI_SYMBOL_EXPORT __attribute__ ((visibility ("default")))
  #else
    #define CASADI_SYMBOL_EXPORT
  #endif
#endif

static const casadi_int casadi_s0[6] = {2, 1, 0, 2, 0, 1};
static const casadi_int casadi_s1[5] = {1, 1, 0, 1, 0};

void casadi_copy(const casadi_real* x, casadi_int n, casadi_real* y) {
  casadi_int i;
  if (y) {
    if (x) {
      for (i=0; i<n; ++i) *y++ = *x++;
    } else {
      for (i=0; i<n; ++i) *y++ = 0.;
    }
  }
}

/* f:(x[2],y)->(r[2],q[2]) */
static int casadi_f0(const casadi_real** arg, casadi_real** res, casadi_int* iw, casadi_real* w, void* mem) {
  casadi_int i;
  casadi_real *rr;
  const casadi_real *cs;
  casadi_real *w0=w+0, w1;
  /* #0: @0 = input[0][0] */
  casadi_copy(arg[0], 2, w0);
  /* #1: output[0][0] = @0 */
  casadi_copy(w0, 2, res[0]);
  /* #2: @1 = input[1][0] */
  w1 = arg[1] ? arg[1][0] : 0;
  /* #3: @1 = sin(@1) */
  w1 = sin( w1 );
  /* #4: @0 = (@1*@0) */
  for (i=0, rr=w0, cs=w0; i<2; ++i) (*rr++)  = (w1*(*cs++));
  /* #5: output[1][0] = @0 */
  casadi_copy(w0, 2, res[1]);
  return 0;
}

CASADI_SYMBOL_EXPORT int f(const casadi_real** arg, casadi_real** res, casadi_int* iw, casadi_real* w, void* mem){
  return casadi_f0(arg, res, iw, w, mem);
}

CASADI_SYMBOL_EXPORT void f_incref(void) {
}

CASADI_SYMBOL_EXPORT void f_decref(void) {
}

CASADI_SYMBOL_EXPORT casadi_int f_n_in(void) { return 2;}

CASADI_SYMBOL_EXPORT casadi_int f_n_out(void) { return 2;}

CASADI_SYMBOL_EXPORT const char* f_name_in(casadi_int i){
  switch (i) {
    case 0: return "x";
    case 1: return "y";
    default: return 0;
  }
}

CASADI_SYMBOL_EXPORT const char* f_name_out(casadi_int i){
  switch (i) {
    case 0: return "r";
    case 1: return "q";
    default: return 0;
  }
}

CASADI_SYMBOL_EXPORT const casadi_int* f_sparsity_in(casadi_int i) {
  switch (i) {
    case 0: return casadi_s0;
    case 1: return casadi_s1;
    default: return 0;
  }
}

CASADI_SYMBOL_EXPORT const casadi_int* f_sparsity_out(casadi_int i) {
  switch (i) {
    case 0: return casadi_s0;
    case 1: return casadi_s0;
    default: return 0;
  }
}

CASADI_SYMBOL_EXPORT int f_work(casadi_int *sz_arg, casadi_int* sz_res, casadi_int *sz_iw, casadi_int *sz_w) {
  if (sz_arg) *sz_arg = 4;
  if (sz_res) *sz_res = 3;
  if (sz_iw) *sz_iw = 0;
  if (sz_w) *sz_w = 3;
  return 0;
}


#ifdef __cplusplus
} /* extern "C" */
#endif

This will create a C file gen.c containing the function f and all its dependencies and required helper functions. We will return to how this file can be used in Section 5.1 and the structure of the generated code is described in Section 5.3 below.

You can generate a C file containing multiple CasADi functions by working with CasADi’s CodeGenerator class:

f = Function('f',[x],[sin(x)])
g = Function('g',[x],[cos(x)])
C = CodeGenerator('gen.c')
C.add(f)
C.add(g)
C.generate()
f = Function('f',{x},{sin(x)});
g = Function('g',{x},{cos(x)});
C = CodeGenerator('gen.c');
C.add(f);
C.add(g);
C.generate();

Both the generate function and the CodeGenerator constructor take an optional options dictionary as an argument, allowing customization of the code generation. Two useful options are main, which generates a main entry point, and mex, which generates a mexFunction entry point:

f = Function('f',[x],[sin(x)])
opts = dict(main=True, \
            mex=True)
f.generate('gen.c',opts)
f = Function('f',{x},{sin(x)});
opts = struct('main', true,...
              'mex', true);
f.generate('gen.c',opts);

This enables executing the function from the command line and MATLAB, respectively, as described in Section 5.2 below.

If you plan to link directly against the generated code in some C/C++ application, a useful option is with_header, which controls the creation of a header file containing declarations of the functions with external linkage, i.e. the API of the generated code, described in Section 5.3 below.

Here is a list of available options for the CodeGenerator class:

Option Default value Description
verbose true Include comments in generated code
mex false Generate an MATLAB/Octave MEX entry point
cpp false Generated code is C++ instead of C
main false Generate a command line interface
casadi_real double Floating point type
casadi_int long long int Integer type
with_header false Generate a header file
with_mem false Generate a simplified C API
indent 2 Number of spaces per indentation level

5.2. Using the generated code

The generated C code can be used in a number of different ways:

  • The code can be compiled into a dynamically linked library (DLL), from which a Function instance can be created using CasADi’s external function. Optionally, the user can rely on CasADi to carry out the compilation just-in-time.
  • The generated code can be compiled into MEX function and executed from MATLAB.
  • The generated code can be executed from the command line.
  • The user can link, statically or dynamically, the generated code to his or her C/C++ application, accessing the C API of the generated code.
  • The code can be compiled into a dynamically linked library and the user can then manually access the C API using dlopen on Linux/OS X or LoadLibrary on Windows.

This is elaborated in the following.

5.2.1. CasADi’s external function

The external command allows the user to create a Function instance from a dynamically linked library with the entry points described by the C API described in Section 5.3. Since the autogenerated files are self-contained [1], the compilation – on Linux/OSX – can be as easy as issuing:

gcc -fPIC -shared gen.c -o gen.so

from the command line. Or, equivalently using MATLAB’s system command or Python’s os.system command. Assuming gen.c was created as described in the previous section, we can then create a Function f as follows:

f = external('f', './gen.so')
print(f(3.14))
0.00159265
f = external('f', './gen.so');
disp(f(3.14))
0.00159265

We can also rely on CasADi performing the compilation just-in-time using CasADi’s Importer class. This is a plugin class, which at the time of writing had two supported plugins, namely 'clang', which invokes the LLVM/Clang compiler framework (distributed with CasADi), and 'shell', which invokes the system compiler via the command line:

C = Importer('gen.c','clang')
f = external('f',C)
print(f(3.14))
0.00159265
C = Importer('gen.c','clang');
f = external('f',C);
disp(f(3.14))
0.00159265

We will return to the external function in Section 6.3.

5.2.2. Calling generated code from MATLAB

An alternative way of executing generated code is to compile the code into a MATLAB MEX function and call from MATLAB. This assumes that the mex option was set to “true” during the code generation, cf. Section 5.1. The generated MEX function takes the function name as its first argument, followed by the function inputs:

%mex gen.c -largeArrayDims  % Matlab
mex gen.c -DMATLAB_MEX_FILE % Octave

disp(gen('f', 3.14))
 0.0015927

Note that the result of the execution is always a MATLAB sparse matrix by default. Compiler flags -DCASASI_MEX_ALWAYS_DENSE and -DCASASI_MEX_ALLOW_DENSE may be set to influence this behaviour.

5.2.3. Calling generated code from the command line

Another option is to execute the generated code from the Linux/OSX command line. This is possible if the main option was set to “true” during the code generation, cf. Section 5.1. This is useful if you e.g. want to profile the generated with a tool such as gprof.

When executing the generated code, the function name is passed as a command line argument. The nonzero entries of all the inputs need to be passed via standard input and the function will return the output nonzeros for all the outputs via standard output:

# Command line
echo 3.14 3.14 > gen_in.txt
gcc gen.c -o gen
./gen f < gen_in.txt > gen_out.txt
cat gen_out.txt # returns 0.00159265 0.00159265

5.2.4. Linking against generated code from a C/C++ application

The generated code is written so that it can be linked with directly from a C/C++ application. If the with_header option was set to “true” during the code generation, a header file with declarations of all the exposed entry points of the file. Using this header file requires an understanding of CasADi’s codegen API, as described in Section 5.3 below. Symbols that are not exposed are prefixed with a file-specific prefix, allowing an application to link against multiple generated functions without risking symbol conflicts.

5.2.5. Dynamically loading generated code from a C/C++ application

A variant of above is to compile the generated code into a shared library, but directly accessing the exposed symbols rather than relying on CasADi’s external function. This also requires an understanding of the structure of the generated code.

In CasADi’s example collection, codegen_usage.cpp demonstrates how this can be done.

5.3. API of the generated code

The API of the generated code consists of a number of functions with external linkage. In addition to the actual execution, there are functions for memory management as well as meta information about the inputs and outputs. These functions are described in the following. Below, assume that the name of function we want to access is fname. To see what these functions actually look like in code and when they are called, we refer to the codegen_usage.cpp example.

5.3.1. Reference counting

void fname_incref(void);
void fname_decref(void);

A generated function may need to e.g. read in some data or initialize some data structures before first call. This is typically not needed for functions generated from CasADi expressions, but may be required e.g. when the generated code contains calls to external functions. Similarly, memory might need to be deallocated after usage.

To keep track of the ownership, the generated code contains two functions for increasing and decreasing a reference counter. They are named fname_incref and fname_decref, respectively. These functions have no input argument and return void.

Typically, some initialization may take place upon the first call to fname_incref and subsequent calls will only increase some internal counter. The fname_decref, on the other hand, decreases the internal counter and when the counter hits zero, a deallocation – if any – takes place.

5.3.2. Number of inputs and outputs

casadi_int fname_n_in(void);
casadi_int fname_n_out(void);

The number of function inputs and outputs can be obtained by calling the fname_n_in and fname_n_out functions, respectively. These functions take no inputs and return the number of input or outputs (casadi_int is an alias for long long int).

5.3.3. Names of inputs and outputs

const char* fname_name_in(casadi_int ind);
const char* fname_name_out(casadi_int ind);

The functions fname_name_in and fname_name_out return the name of a particular input or output. They take the index of the input or output, starting with index 0, and return a const char* with the name as a null-terminated C string. Upon failure, these functions will return a null pointer.

5.3.4. Sparsity patterns of inputs and outputs

const casadi_int* fname_sparsity_in(casadi_int ind);
const casadi_int* fname_sparsity_out(casadi_int ind);

The sparsity pattern for a given input or output is obtained by calling fname_sparsity_in and fname_sparsity_out, respectively. These functions take the input or output index and return a pointer to a field of constant integers (const casadi_int*). This is a compact representation of the compressed column storage (CCS) format that CasADi uses, cf. Section 3.5. The integer field pointed to is structured as follows:

  • The first two entries are the number of rows and columns, respectively. In the following referred to as nrow and ncol.
  • If the third entry is 1, the pattern is dense and any remaining entries are discarded.
  • If the third entry is 0, that entry plus subsequent ncol entries form the nonzero offsets for each column, colind in the following. E.g. column \(i\) will consist of the nonzero indices ranging from colind[i] to colind[i+1]. The last entry, colind[ncol], will be equal to the number of nonzeros, nnz.
  • Finally, if the sparsity pattern is not dense, i.e. if nnz \(\ne\) nrow * ncol, then the last nnz entries will contain the row indices.

Upon failure, these functions will return a null pointer.

5.3.5. Memory objects

A function may contain some mutable memory, e.g. for caching the latest factorization or keeping track of evaluation statistics. When multiple functions need to call the same function without conflicting, they each need to work with a different memory object. This is especially important for evaluation in parallel on a shared memory architecture, in which case each thread should access a different memory object.

void* fname_alloc_mem(void);

Allocates a memory object which will be passed to the numerical evaluation.

int fname_init_mem(void* mem);

(Re)initializes a memory object. Returns 0 upon successful return;

int fname_free_mem(void* mem);

Frees a memory object. Returns 0 upon successful return;

5.3.6. Work vectors

int fname_work(casadi_int* sz_arg, casadi_int* sz_res, casadi_int* sz_iw, casadi_int* sz_w);

To allow the evaluation to be performed efficiently with a small memory footprint, the user is expected to pass four work arrays. The function fname_work returns the length of these arrays, which have entries of type const double*, double*, casadi_int and double, respectively.

The return value of the function is nonzero upon failure.

5.3.7. Numerical evaluation

int fname(const double** arg, double** res,
          casadi_int* iw, double* w, void* mem);

Finally, the function fname, performs the actual evaluation. It takes as input arguments the four work vectors and a memory object created using fname_alloc_mem (or NULL if absent). The length of the work vectors must be at least the lengths provided by the fname_work command and the index of the memory object must be strictly smaller than the value returned by fname_n_mem.

The nonzeros of the function inputs are pointed to by the first entries of the arg work vector and are unchanged by the evaluation. Similarly, the output nonzeros are pointed to by the first entries of the res work vector and are also unchanged (i.e. the pointers are unchanged, not the actual values).

The return value of the function is nonzero upon failure.

Footnotes

[1]An exception is when code is generated for a function that in turn contains calls to external functions.

6. User-defined function objects

There are situations when rewriting user-functions using CasADi symbolics is not possible or practical. To tackle this, CasADi provides a number of ways to embed a call to a “black box” function defined in the language CasADi is being used from (C++, MATLAB or Python) or in C. That being said, the recommendation is always to try to avoid this when possible, even if it means investing a lot of time reimplementing existing code. Functions defined using CasADi symbolics are almost always more efficient, especially when derivative calculation is involved, since a lot more structure can typically be exploited.

Depending on the circumstances, the user can implement custom Function objects in a number of different ways, which will be elaborated on in the following sections:

6.1. Subclassing FunctionInternal

All function objects presented in Chapter 4 are implemented in CasADi as C++ classes inheriting from the FunctionInternal abstract base class. In principle, a user with familiarity with C++ programming, can implement a class inheriting from FunctionInternal, overloading the virtual methods of this class. The best reference for doing so is the C++ API documentation, choosing “switch to internal” to expose the internal API.

Since FunctionInternal is not considered part of the stable, public API, we advice against this in general, unless the plan is to contribution to CasADi’s source.

6.2. Subclassing Callback

The Callback class provides a public API to FunctionInternal and inheriting from this class has the same effect as inheriting directly from FunctionInternal. Thanks to cross-language polymorphism, it is possible to implement the exposed methods of Callback from either Python, MATLAB/Octave or C++.

The derived class consists of the following parts:

  • A constructor or a static function replacing the constructor
  • A number of virtual functions, all optional, that can be overloaded in order to get the desired behavior. This includes the number of inputs and outputs using get_n_in and get_n_out, their names using get_name_in and get_name_out and their sparsity patterns get_sparsity_in and get_sparsity_out.
  • An optional init function called duing the object construction.
  • A function for numerical evaluation.
  • Optional functions for derivatives. You can choose to supply a full Jacobian (has_jacobian, get_jacobian), or choose to supply forward/reverse sensitivities (has_forward, get_forward, has_reverse, get_reverse).

For a complete list of functions, see the C++ API documentation for Callback.

The usage from the different languages are described in the following.

6.2.1. Python

class MyCallback(Callback):
  def __init__(self, name, d, opts={}):
    Callback.__init__(self)
    self.d = d
    self.construct(name, opts)

  # Number of inputs and outputs
  def get_n_in(self): return 1
  def get_n_out(self): return 1

  # Initialize the object
  def init(self):
     print('initializing object')

  # Evaluate numerically
  def eval(self, arg):
    x = arg[0]
    f = sin(self.d*x)
    return [f]

The implementation should include a constructor, which should begin with a call to the base class constructor using Callback.__init__(self) and end with a call to initialize object construction using self.construct(name, opts).

This function can be used as any built-in CasADi function with the important caveat that when embedded in graphs, the ownership of the class will not be shared between all references. So it is important that the user does not allow the Python class to go out of scope while it is still needed in calculations.

# Use the function
f = MyCallback('f', 0.5)
print(f(2))

x = MX.sym("x")
print(f(x))
initializing object
0.841471
f(x){0}

6.2.2. MATLAB

In MATLAB, a custom function class can be defined as follows, in a file MyCallback.m:

classdef MyCallback < casadi.Callback
  properties
    d
  end
  methods
    function self = MyCallback(name, d)
      self@casadi.Callback();
      self.d = d;
      construct(self, name);
    end

    % Number of inputs and outputs
    function v=get_n_in(self)
      v=1;
    end
    function v=get_n_out(self)
      v=1;
    end

    % Initialize the object
    function init(self)
      disp('initializing object')
    end

    % Evaluate numerically
    function arg = eval(self, arg)
      x = arg{1};
      f = sin(self.d * x);
      arg = {f};
    end
  end
end

This function can be used as any built-in CasADi function, but as for Python, the ownership of the class will not be shared between all references. So the user must not allow a class instance to get deleted while it is still in use, e.g. by making it persistent.

% Use the function
f = MyCallback('f', 0.5);
res = f(2);
disp(res)

x = MX.sym('x');
disp(f(x))

6.2.3. C++

In C++, the syntax is as follows:

#include "casadi/casadi.hpp"
using namespace casadi;
class MyCallback : public Callback {
  // Data members
  double d;
public:
  // Constructor
  MyCallback(const std::string& name, double d,
             const Dict& opts=Dict()) : d(d) {
    construct(name, opts);
  }

  // Destructor
  ~MyCallback() override {}

  // Number of inputs and outputs
  casadi_int get_n_in() override { return 1;}
  casadi_int get_n_out() override { return 1;}

  // Initialize the object
  void init override() {
    std::cout << "initializing object" << std::endl;
  }

  // Evaluate numerically
  std::vector<DM> eval(const std::vector<DM>& arg) const override {
    DM x = arg.at(0);
    DM f = sin(d*x);
    return {f};
  }
};

A class created this way can be used as any other Function instance, but with the important difference that the user is responsible to managing the memory of this class.

int main() {
  MyCallback f("f", 0.5);
  std::vector<DM> arg = {2};
  std::vector<DM> res = f(arg);
  std::cout << res << std::endl;
  return 0;
}

6.3. Importing a function with external

The basic usage of CasADi’s external function was demonstrated in Section 5.2 in the context of using autogenerated code. The same function can also be used for importing a user-defined function, as long as it also uses the C API described in Section 5.3.

The following sections expands on this.

6.3.1. Default functions

It is usually not necessary to define all the functions defined in Section 5.3. If fname_incref and fname_decref are absent, it is assumed that no memory management is needed. If no names of inputs and outputs are provided, they will be given default names. Sparsity patterns are in general assumed to be scalar by default, unless the function corresponds to a derivative of another function (see below), in which case they are assumed to be dense and of the correct dimension.

Furthermore, work vectors are assumed not to be needed if fname_work has not been implemented.

6.3.2. Meta information as comments

If you rely on CasADi’s just-in-time compiler, you can provide meta information as a comment in the C code instead of implementing the actual callback function.

The structure of such meta information should be as follows:

/*CASADIMETA
:fname_N_IN 1
:fname_N_OUT 2
:fname_NAME_IN[0] x
:fname_NAME_OUT[0] r
:fname_NAME_OUT[1] s
:fname_SPARSITY_IN[0] 2 1 0 2
*/

6.3.3. Derivatives

The external function can be made differentiable by providing functions for calculating derivatives. During derivative calculations, CasADi will look for symbols in the same file/shared library that follows a certain naming convention. For example, you can specify a Jacobian for all the outputs with respect to all inputs for a function named fname by implementing a function named jac_fname. Similary, you can specify a function for calculating one forward directional derivative by providing a function named fwd1_fname, where 1 can be replaced by 2, 4, 8, 16, 32 or 64 for calculating multiple forward directional derivatives at once. For reverse mode directional derivatives, replace fwd with adj.

This is an experimental feature.

6.4. Just-in-time compile a C language string

In the previous section we showed how to specify a C file with functions for numerical evaluation and meta information. As was shown, this file can be just-in-time compiled by CasADi’s interface to Clang. There exists a shorthand for this approach, where the user simply specifies the source code as a C language string.

body =\
'r[0] = x[0];'+\
'while (r[0]<s[0]) {'+\
' r[0] *= r[0];'+\
'}'

f = Function.jit('f',body,\
      ['x','s'],['r'])
print(f)
f:(x,s)->(r) JitFunction
body =[...
'r[0] = x[0];',...
'while (r[0]<s[0]) {',...
' r[0] *= r[0];',...
'}'];

f = Function.jit('f',body,...
      {'x','s'},{'r'});
disp(f)
f:(x,s)->(r) JitFunction

These four arguments of Function.jit are mandatory: The function name, the C source as a string and the names of inputs and outputs. In the C source, the input/output names correspond to arrays of type casadi_real_t containing the nonzero elements of the function inputs and outputs. By default, all inputs and outputs are scalars (i.e. 1-by-1 and dense). To specify a different sparsity pattern, provide two additional function arguments containing vectors/lists of the sparsity patterns:

sp = Sparsity.scalar()
f = Function.jit('f',body,\
     ['x','s'], ['r'],\
     [sp,sp], [sp])
print(f)
f:(x,s)->(r) JitFunction
sp = Sparsity.scalar();
f = Function.jit('f',body,...
     {'x','s'}, {'r'},...
     {sp,sp}, {sp});
disp(f)
f:(x,s)->(r) JitFunction

Both variants accept an optional 5th (or 7th) argument in the form of an options dictionary.

6.5. Using lookup-tables

Lookup-tables can be created using CasADi’s interpolant function. Different interpolating schemes are implemented as plugins, similar to nlpsol or integrator objects. In addition to the identifier name and plugin, the interpolant function expects a set of grid points with the corresponding numerical values.

The result of an interpolant call is a CasADi Function object that is differentiable, and can be embedded into CasADi computational graphs by calling with MX arguments. Furthermore, C code generation is fully supported for such graphs.

Currently, two plugins exist for interpolant: 'linear' and 'bspline'. They are intended to behave simiarly to MATLAB/Octave’s interpn with the method set to 'linear' or 'spline' – corresponding to a multilinear interpolation and a (by default cubic) spline interpolation with not-a-knot boundary conditions.

In the case of bspline, coefficients will be sought at construction time that fit the provided data. Alternatively, you may also use the more low-level Function.bspline to supply the coefficients yourself. The default degree of the bspline is 3 in each dimension. You may deviate from this default by passing a degree option.

We will walk through the syntax of interpolant for the 1D and 2D versions, but the syntax in fact generalizes to an arbitrary number of dimensions.

6.5.1. 1D lookup tables

A 1D spline fit can be done in CasADi/Python as follows, compared with the corresponding method in SciPy:

import casadi as ca
import numpy as np
xgrid = np.linspace(1,6,6)
V = [-1,-1,-2,-3,0,2]
lut = ca.interpolant('LUT','bspline',[xgrid],V)
print(lut(2.5))
# Using SciPy
import scipy.interpolate as ip
interp = ip.InterpolatedUnivariateSpline(xgrid, V)
print(interp(2.5))
-1.35
-1.3500000000000005

In MATLAB/Octave, the corresponding code reads:

xgrid = 1:6;
V = [-1 -1 -2 -3 0 2];
lut = casadi.interpolant('LUT','bspline',{xgrid},V);
lut(2.5)
% Using MATLAB/Octave builtin
interp1(xgrid,V,2.5,'spline')
ans =

-1.35

ans = -1.3500

Note in particular that the grid and values arguments to interpolant must be numerical in nature.

6.5.2. 2D lookup tables

In two dimensions, we get the following in Python, also compared to SciPy for reference:

xgrid = np.linspace(-5,5,11)
ygrid = np.linspace(-4,4,9)
X,Y = np.meshgrid(xgrid,ygrid,indexing='ij')
R = np.sqrt(5*X**2 + Y**2)+ 1
data = np.sin(R)/R
data_flat = data.ravel(order='F')
lut = ca.interpolant('name','bspline',[xgrid,ygrid],data_flat)
print(lut([0.5,1]))
# Using Scipy

interp = ip.RectBivariateSpline(xgrid, ygrid, data)
print(interp.ev(0.5,1))
0.245507
0.24550661674668917

or, in MATLAB/Octave compared to the built-in functions:

xgrid = -5:1:5;
ygrid = -4:1:4;
[X,Y] = ndgrid(xgrid, ygrid);
R = sqrt(5*X.^2 + Y.^2)+ 1;
V = sin(R)./R;
lut = interpolant('LUT','bspline',{xgrid, ygrid},V(:));
lut([0.5 1])
% Using Matlab builtin
interpn(X,Y,V,0.5,1,'spline')
ans =

0.245507

ans =  0.24551

In particular note how the values argument had to be flatten to a one-dimensional array.

6.6. Derivative calculation using finite differences

CasADi 3.3 introduced support for finite difference calculation for all function objects, in particular including external functions defined as outlined in Section 6.2, Section 6.3 or Section 6.4 (for lookup tables, Section 6.5, analytical derivatives are available).

Finite difference derivative are disabled by default, with the exception of Function.jit, and to enable it, you must set the option 'enable_fd' to True/true:

f = external('f', './gen.so',\
   dict(enable_fd=True))

e = jacobian(f(x),x)
D = Function('d',[x],[e])
print(D(0))
1
f = external('f', './gen.so',...
    struct('enable_fd',true));

e = jacobian(f(x),x);
D = Function('d',{x},{e});
disp(D(0))
1

cf. Section 5.1.

The 'enable_fd' options enables CasADi to use finite differences, if analytical derivatives are not available. To force CasADi to use finite differences, you can set 'enable_forward', 'enable_reverse' and 'enable_jacobian' to False/false, corresponding to the three types of analytical derivative information that CasADi works with.

The default method is central differences with a step size determined by estimates of roundoff errors and truncation errors of the function. You can change the method by setting the option 'fd_method' to 'forward' (corresponding to first order forward differences), 'backward' (corresponding to first order backward differences) and 'smoothing' for a second-order accurate discontinuity avoiding scheme, suitable when derivatives need to be calculated at the edges of a domain. Additional algorithmic options for the finite differences are available by setting 'fd_options' option.

7. The DaeBuilder class

The DaeBuilder class in CasADi is an auxiliary class intended to facilitate the modeling complex dynamical systems for later use with optimal control algorithms. This class can be seen as a low-level alternative to a physical modeling language such as Modelica (cf. Section 7.3), while still being higher level than working directly with CasADi symbolic expressions. Another important usage it to provide an interface to physical modeling languages and software and be a building blocks for developing domain specific modeling environments.

Using the DaeBuilder class consists of the following steps:

  • Step-by-step constructing a structured system of differential-algebraic equations (DAE) or, alternatively, importing an existing model from Modelica
  • Symbolically reformulate the DAE
  • Generate a chosen set of CasADi functions to be used for e.g. optimal control or C code generation

In the following sections, we describe the mathematical formulation of the class and its intended usage.

7.1. Mathematical formulation

The DaeBuilder class uses a relatively rich problem formulation that consists of a set of input expressions and a set of output expressions, each defined by a string identifier. The choice of expressions was inspired by the functional mock-up interface (FMI) version 2.0 [1]

7.1.1. Input expressions

‘t’
Time \(t\)
‘c’
Named constants \(c\)
‘p’
Independent parameters \(p\)
‘d’
Dependent parameters \(d\), depends only on \(p\) and \(c\) and, acyclically, on other \(d\)
‘x’
Differential state \(x\), defined by an explicit ODE
‘s’
Differential state \(s\), defined by an implicit ODE
‘sdot’
Time derivatives implicitly defined differential state \(\dot{s}\)
‘z’
Algebraic variable, defined by an algebraic equation
‘q’
Quadrature state \(q\). A differential state that may not appear in the right-hand-side and hence can be calculated by quadrature formulas.
‘w’
Local variables \(w\). Calculated from time and time dependent variables. They may also depend, acyclically, on other \(w\).
‘y’
Output variables \(y\)

7.1.2. Output expressions

The above input expressions are used to define the following output expressions:

‘ddef’
Explicit expression for calculating \(d\)
‘wdef’
Explicit expression for calculating \(w\)
‘ode’
The explicit ODE right-hand-side:
\(\dot{x} = \text{ode}(t,w,x,s,z,u,p,d)\)
‘dae’
The implicit ODE right-hand-side: \(\text{dae}(t,w,x,s,z,u,p,d,\dot{s}) =0\)
‘alg’
The algebraic equations:
\(\text{alg}(t,w,x,s,z,u,p,d) = 0\)
‘quad’
The quadrature equations: \(\dot{q} = \text{quad}(t,w,x,s,z,u,p,d)\)
‘ydef’
Explicit expressions for calculating \(y\)

7.2. Constructing a DaeBuilder instance

Consider the following simple DAE corresponding to a controlled rocket subject to quadratic air friction term and gravity, which loses mass as it uses up fuel:

\[\begin{split}\begin{aligned} \dot{h} &= v, \qquad &h(0) = 0 \\ \dot{v} &= (u - a \, v^2)/m - g, \qquad &v(0) = 0 \\ \dot{m} &= -b \, u^2, \qquad &m(0) = 1 \end{aligned}\end{split}\]

where the three states correspond to height, velocity and mass, respectively. \(u\) is the thrust of the rocket and \((a,b)\) are parameters.

To construct a DAE formulation for this problem, start with an empty DaeBuilder instance and add the input and output expressions step-by-step as follows.

dae = DaeBuilder()
# Add input expressions
a = dae.add_p('a')
b = dae.add_p('b')
u = dae.add_u('u')
h = dae.add_x('h')
v = dae.add_x('v')
m = dae.add_x('m')
# Add output expressions
hdot = v
vdot = (u-a*v**2)/m-g
mdot = -b*u**2
dae.add_ode('hdot', hdot)
dae.add_ode('vdot', vdot)
dae.add_ode('mdot', mdot)
# Specify initial conditions
dae.set_start('h', 0)
dae.set_start('v', 0)
dae.set_start('m', 1)
# Add meta information
dae.set_unit('h','m')
dae.set_unit('v','m/s')
dae.set_unit('m','kg')
dae = DaeBuilder;
% Add input expressions
a = dae.add_p('a');
b = dae.add_p('b');
u = dae.add_u('u');
h = dae.add_x('h');
v = dae.add_x('v');
m = dae.add_x('m');
% Add output expressions
hdot = v;
vdot = (u-a*v^2)/m-g;
mdot = -b*u^2;
dae.add_ode('hdot', hdot);
dae.add_ode('vdot', vdot);
dae.add_ode('mdot', mdot);
% Specify initial conditions
dae.set_start('h', 0);
dae.set_start('v', 0);
dae.set_start('m', 1);
% Add meta information
dae.set_unit('h','m');
dae.set_unit('v','m/s');
dae.set_unit('m','kg');

Other input and output expressions can be added in an analogous way. For a full list of functions, see the C++ API documentation for DaeBuilder.

7.3. Import of OCPs from Modelica

An alternative to model directly in CasADi, as above, is to use an advanced physical modeling language such as Modelica to specify the model. For this, CasADi offers interoperability with the open-source JModelica.org compiler, which is written specifically with optimal control in mind. Model import from JModelica.org is possible in two different ways; using the JModelica.org’s CasadiInterface or via DaeBuilder’s parse_fmi command.

We recommend the former approach, since it is being actively maintained and refer to JModelica.org’s user guide for details on how to extract CasADi expressions.

In the following, we will outline the legacy approach, using parse_fmi.

7.3.1. Legacy import of a modelDescription.xml file

To see how to use the Modelica import, look at thermodynamics_example.py in CasADi’s example collection.

Assuming that the Modelica/Optimica model ModelicaClass.ModelicaModel is defined in the files file1.mo and file2.mop, the Python compile command is:

from pymodelica import compile_jmu
jmu_name=compile_jmu('ModelicaClass.ModelicaModel', \
  ['file1.mo','file2.mop'],'auto','ipopt',\
  {'generate_xml_equations':True, 'generate_fmi_me_xml':False})

This will generate a jmu-file, which is essentially a zip file containing, among other things, the file modelDescription.xml. This XML-file contains a symbolic representation of the optimal control problem and can be inspected in a standard XML editor.

from zipfile import ZipFile
sfile = ZipFile(jmu_name','r')
mfile = sfile.extract('modelDescription.xml','.')

Once a modelDescription.xml file is available, it can be imported using the parse_fmi command:

dae = DaeBuilder()
ocp.parse_fmi('modelDescription.xml')

7.4. Symbolic reformulation

One of the original purposes of the DaeBuilder class was to reformulate a fully-implicit DAE, typically coming from Modelica, to a semi-explicit DAE that can be used more readily in optimal control algorithms.

This can be done by the python{make_implicit} command:

ocp.make_explicit()
ocp.make_explicit();

Other useful commands available for an instance ocp of DaeBuilder include:

  • print ocp Print the optimal optimal control problem to screen
  • ocp.scale_variables() Scale all variables using the nominal attribute for each variable
  • ocp.eliminate_d() Eliminate all independent parameters from the symbolic expressions

For a more detailed description of this class and its functionalities, we again refer to the API documentation.

7.5. Function factory

Once a DaeBuilder has been formulated and possibly reformulated to a satisfactory form, we can generate CasADi functions corresponding to the input and output expressions outlined in Section 7.1. For example, to create a function for the ODE right-hand-side for the rocket model in Section 7.2, simply provide a display name of the function being created, a list of input expressions and a list of output expressions:

f = dae.create('f',\
     ['x','u','p'],['ode'])
f = dae.create('f',...
     {'x','u','p'},{'ode'});

Using a naming convention, we can also create Jacobians, e.g. for the ‘ode’ output with respect to ‘x’:

f = dae.create('f',\
     ['x','u','p'],\
     ['jac_ode_x'])
f = dae.create('f',...
     {'x','u','p'},
     {'jac_ode_x'});

Functions with second order information can be extracted by first creating a named linear combination of the output expressions using add_lc and then requesting its Hessian:

dae.add_lc('gamma',['ode'])
hes = dae.create('hes',\
  ['x','u','p','lam_ode'],\
  ['hes_gamma_x_x'])
dae.add_lc('gamma',{'ode'});
hes = dae.create('hes',...
  {'x','u','p','lam_ode'},...
  {'hes_gamma_x_x'});

It is also possible to simply extract the symbolic expressions from the DaeBuilder instance and manually create CasADi functions. For example, dae.x contains all the expressions corresponding to ‘x’, dae.ode contains the expressions corresponding to ‘ode’, etc.

Footnotes

[1]FMI development group. Functional Mock-up Interface for Model Exchange and Co-Simulation. https://www.fmi-standard.org/, July 2014. Specification, FMI 2.0. Section 3.1, pp. 71–72

8. Optimal control with CasADi

CasADi can be used to solve optimal control problems (OCP) using a variety of methods, including direct (a.k.a. discretize-then-optimize) and indirect (a.k.a. optimize-then-discretize) methods, all-at-once (e.g. collocation) methods and shooting-methods requiring embedded solvers of initial value problems in ODE or DAE. As a user, you are in general expected to write your own OCP solver and CasADi aims as making this as easy as possible by providing powerful high-level building blocks. Since you are writing the solver yourself (rather than calling an existing “black-box” solver), a basic understanding of how to solve OCPs is indispensable. Good, self-contained introductions to numerical optimal control can be found in the recent textbooks by Biegler [1] or Betts [2] or Moritz Diehl’s lecture notes on numerical optimal control.

8.1. A simple test problem

To illustrate some of the methods, we will consider the following test problem, namely driving a Van der Pol oscillator to the origin, while trying to minimize a quadratic cost:

(8.1.1)\[\begin{split}\begin{array}{lc} \begin{array}{l} \text{minimize:} \\ x(\cdot) \in \mathbb{R}^2, \, u(\cdot) \in \mathbb{R} \end{array} \quad \displaystyle \int_{t=0}^{T}{(x_0^2 + x_1^2 + u^2) \, dt} \\ \\ \text{subject to:} \\ \\ \begin{array}{ll} \left\{ \begin{array}{l} \dot{x}_0 = (1-x_1^2) \, x_0 - x_1 + u \\ \dot{x}_1 = x_0 \\ -1.0 \le u \le 1.0, \quad x_1 \ge -0.25 \end{array} \right. & \text{for} \, 0 \le t \le T \\ x_0(0)=0, \quad x_1(0)=1, \end{array} \end{array}\end{split}\]

with \(T=10\).

In CasADi’s examples collection [3], you find codes for solving optimal control problems using a variety of different methods.

In the following, we will discuss three of the most important methods, namely direct single shooting, direct multiple shooting and direct collocation.

8.1.1. Direct single-shooting

In the direct single shooting method, the control trajectory is parameterized using some piecewise smooth approximation, typically piecewise constant.

Using an explicit expression for the controls, we can then eliminate the whole state trajectory from the optimization problem, ending up with an NLP in only the discretized controls.

In CasADi’s examples collection, you will find the codes direct_single_shooting.py and direct_single_shooting.m for Python and MATLAB/Octave, respectively. These codes implement the direct single shooting method and solves it with IPOPT, relying on CasADi to calculate derivatives. To obtain the discrete time dynamics from the continuous time dynamics, a simple fixed-step Runge-Kutta 4 (RK4) integrator is implemented using CasADi symbolics. Simple integrator codes like these are often useful in the context of optimal control, but care must be taken so that they accurately solve the initial-value problem.

The code also shows how the RK4 scheme can be replaced by a more advanced integrator, namely the CVODES integrator from the SUNDIALS suite, which implements a variable stepsize, variable order backward differentiation formula (BDF) scheme. An advanced integrator like this is useful for larger systems, systems with stiff dynamics, for DAEs and for checking a simpler scheme for consistency.

8.1.2. Direct multiple-shooting

The direct_multiple_shooting.py and direct_multiple_shooting.m codes, also in CasADi’s examples collection, implement the direct multiple shooting method. This is very similar to the direct single shooting method, but includes the state at certain shooting nodes as decision variables in the NLP and includes equality constraints to ensure continuity of the trajectory.

The direct multiple shooting method is often superior to the direct single shooting method, since “lifting” the problem to a higher dimension is known to often improve convergence. The user is also able to initialize with a known guess for the state trajectory.

The drawback is that the NLP solved gets much larger, although this is often compensated by the fact that it is also much sparser.

8.1.3. Direct collocation

Finally, the direct_collocation.py and direct_collocation.m codes implement the direct collocation method. In this case, a parameterization of the entire state trajectory, as piecewise low-order polynomials, are included as decision variables in the NLP. This removes the need for the formulation of the discrete time dynamics completely.

The NLP in direct collocation is even larger than that in direct multiple shooting, but is also even sparser.

Footnotes

[1]Lorenz T. Biegler, Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical Processes , SIAM 2010
[2]John T. Betts, Practical Methods for Optimal Control Using Nonlinear Programming , SIAM 2001
[3]You can obtain this collection as an archive named examples_pack.zip in CasADi’s download area

9. Opti stack

The Opti stack is a collection of CasADi helper classes that provides a close correspondence between mathematical NLP notation, e.g.

(9.1)\[\begin{split} \begin{array}{cc} \begin{array}{c} \text{minimize} \\ x,y \end{array} & (y-x^2)^2 \\ \begin{array}{c} \text{subject to} \end{array} & x^2+y^2=1 \\ & x+y \ge 1, \end{array}\end{split}\]

and computer code:

opti = casadi.Opti()

x = opti.variable()
y = opti.variable()

opti.minimize(  (y-x**2)**2   )
opti.subject_to( x**2+y**2==1 )
opti.subject_to(       x+y>=1 )

opti.solver('ipopt')


sol = opti.solve()

print(sol.value(x))
print(sol.value(y))

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

0.7861513776531158
0.6180339888825889
opti = casadi.Opti();

x = opti.variable();
y = opti.variable();

opti.minimize(  (y-x^2)^2   );
opti.subject_to( x^2+y^2==1 );
opti.subject_to(     x+y>=1 );

opti.solver('ipopt');

sol = opti.solve();

sol.value(x)
sol.value(y)

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

ans =  0.78615
ans =  0.61803

The main characteristics of the Opti stack are:

  • Allows natural syntax for constraints.
  • Indexing/bookkeeping of decision variables is hidden.
  • Closer mapping of numerical data-type to the host language: no encounter with DM.

9.1. Problem specification

Variables: Declare any amount of decision variables:

  • x = opti.variable(): scalar
  • x = opti.variable(5): column vector
  • x = opti.variable(5,3): matrix
  • x = opti.variable(5,5,'symmetric'): symmetric matrix

The order in which you declare the variables is respected by the solver. Note that the variables are in fact plain MX symbols. You may perform any CasADi MX operations on them, e.g. embedding integrator calls.

Parameters: Declare any amount of parameters. You must fix them to a specific numerical value before solving, and you may overwrite this value at any time.

p = opti.parameter()
opti.set_value(p, 3)
p = opti.parameter();
opti.set_value(p, 3)

Objective: Declare an objective using an expression that may involve all variables or parameters. Calling the command again discards the old objective.

opti.minimize(   sin(x*(y-p))   )
opti.minimize(   sin(x*(y-p))   )

Constraints: Declare any amount of equality/inequality constraints:

  • opti.subject_to( sqrt(x+y) >= 1): inequality
  • opti.subject_to( sqrt(x+y) > 1): same as above
  • opti.subject_to( 1<= sqrt(x+y) ): same as above
  • opti.subject_to( 5*x+y==1 ): equality

You may also throw in several constraints at once:

opti.subject_to([x*y>=1,x==3])
opti.subject_to({x*y>=1,x==3});

You may declare double inequalities:

opti.subject_to( opti.bounded(0,x,1) )
opti.subject_to( 0<=x<=1 );

When the bounds of the double inequalities are free of variables, the constraint will be passed on efficiently to solvers that support them (notably IPOPT).

You may make element-wise (in)equalities with vectors:

x = opti.variable(5,1)

opti.subject_to( x*p<=3 )
x = opti.variable(5,1);

opti.subject_to( x*p<=3 )

Elementwise (in)equalities for matrices are not supported with a natural syntax, since there is an ambiguity with semi-definiteness constraints. The workaround is to vectorize first:

A = opti.variable(5,5)
opti.subject_to( vec(A)<=3 )
A = opti.variable(5,5);
opti.subject_to( A(:)<=3 );

Each subject_to command adds to the set of constraints in the problem specification. Use subject_to() to empty this set and start over.

Solver: You must always declare the solver (numerical back-end). An optional dictionary of CasADi plugin options can be given as second argument. An optional dictionary of solver options can be given as third argument.

p_opts = {"expand":True}
s_opts = {"max_iter": 100}
opti.solver("ipopt",p_opts,
                    s_opts)
p_opts = struct('expand',true);
s_opts = struct('max_iter',100);
opti.solver('ipopt',p_opts,
                    s_opts);

Initial guess: You may provide initial guesses for decision variables (or simple mappings of decision variables). When no initial guess is provided, numerical zero is assumed.

opti.set_initial(x, 2)
opti.set_initial(10*x[0], 2)
opti.set_initial(x, 2);
opti.set_initial(10*x(1), 2)

9.2. Problem solving and retrieving

9.2.1. Solving

After setting up the problem, you may call the solve method, which constructs a CasADi nlpsol and calls it.

sol = opti.solve()

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.3, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:        2
Number of nonzeros in inequality constraint Jacobian.:        2
Number of nonzeros in Lagrangian Hessian.............:        3

Total number of variables............................:        2
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        1
Total number of inequality constraints...............:        1
        inequality constraints with only lower bounds:        1
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 1.00e+00 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.5180703e+00 2.32e-01 3.95e+08  -1.0 1.11e+00    -  9.90e-01 1.00e+00h  1
   2  9.1720733e-01 1.38e-02 8.65e+08  -1.0 1.05e-01  10.0 1.00e+00 1.00e+00f  1
   3  8.9187743e-01 4.67e-05 9.29e+06  -1.0 7.19e-03    -  1.00e+00 1.00e+00h  1
   4  8.9179090e-01 5.46e-10 2.17e+02  -1.0 2.42e-05    -  1.00e+00 1.00e+00h  1
   5  8.5961118e-01 2.43e-04 3.39e+00  -1.0 1.56e-02    -  1.00e+00 1.00e+00f  1
   6  5.0768568e-01 3.58e-02 2.18e+00  -1.0 1.89e-01    -  3.86e-01 1.00e+00f  1
   7  4.4751644e-01 4.12e-04 6.30e+07  -1.0 1.96e-02   9.5 1.00e+00 1.00e+00h  1
   8  4.4707538e-01 4.25e-08 2.27e+04  -1.0 2.53e-04    -  1.00e+00 1.00e+00h  1
   9  4.4688522e-01 9.34e-09 2.20e+00  -1.0 9.32e-05    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  2.9995484e-02 3.46e-01 1.19e+00  -1.0 5.68e-01    -  5.45e-01 1.00e+00f  1
  11  3.9964883e-04 3.56e-02 5.81e-02  -1.0 2.09e-01    -  1.00e+00 1.00e+00h  1
  12  3.7581923e-06 5.63e-04 5.70e-03  -2.5 2.69e-02    -  1.00e+00 1.00e+00h  1
  13  7.6573980e-10 1.48e-06 3.13e-05  -3.8 1.11e-03    -  1.00e+00 1.00e+00h  1
  14  4.8778322e-14 2.53e-10 7.52e-09  -5.7 1.29e-05    -  1.00e+00 1.00e+00h  1
  15  8.8029044e-20 1.60e-14 6.47e-13  -8.6 9.88e-08    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 15

                                   (scaled)                 (unscaled)
Objective...............:   8.8029044107981477e-20    8.8029044107981477e-20
Dual infeasibility......:   6.4749359007033628e-13    6.4749359007033628e-13
Constraint violation....:   1.5987211554602254e-14    1.5987211554602254e-14
Complementarity.........:   2.5060006300474834e-09    2.5060006300474834e-09
Overall NLP error.......:   2.5060006300474834e-09    2.5060006300474834e-09


Number of objective function evaluations             = 16
Number of objective gradient evaluations             = 16
Number of equality constraint evaluations            = 16
Number of inequality constraint evaluations          = 16
Number of equality constraint Jacobian evaluations   = 16
Number of inequality constraint Jacobian evaluations = 16
Number of Lagrangian Hessian evaluations             = 15
Total CPU secs in IPOPT (w/o function evaluations)   =      0.005
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
               t_proc [s]   t_wall [s]    n_eval
       nlp_f      2.6e-05      2.3e-05        16
       nlp_g      5.3e-05     4.74e-05        16
  nlp_grad_f      3.7e-05     3.64e-05        17
  nlp_hess_l      6.2e-05     5.95e-05        15
   nlp_jac_g        5e-05      4.9e-05        17
      solver      0.00608      0.00608         1
sol = opti.solve();

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.3, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:        2
Number of nonzeros in inequality constraint Jacobian.:        2
Number of nonzeros in Lagrangian Hessian.............:        3

Total number of variables............................:        2
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        1
Total number of inequality constraints...............:        1
        inequality constraints with only lower bounds:        1
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 1.00e+00 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.5180703e+00 2.32e-01 3.95e+08  -1.0 1.11e+00    -  9.90e-01 1.00e+00h  1
   2  9.1720733e-01 1.38e-02 8.65e+08  -1.0 1.05e-01  10.0 1.00e+00 1.00e+00f  1
   3  8.9187743e-01 4.67e-05 9.29e+06  -1.0 7.19e-03    -  1.00e+00 1.00e+00h  1
   4  8.9179090e-01 5.46e-10 2.17e+02  -1.0 2.42e-05    -  1.00e+00 1.00e+00h  1
   5  8.5961118e-01 2.43e-04 3.39e+00  -1.0 1.56e-02    -  1.00e+00 1.00e+00f  1
   6  5.0768568e-01 3.58e-02 2.18e+00  -1.0 1.89e-01    -  3.86e-01 1.00e+00f  1
   7  4.4751644e-01 4.12e-04 6.30e+07  -1.0 1.96e-02   9.5 1.00e+00 1.00e+00h  1
   8  4.4707538e-01 4.25e-08 2.27e+04  -1.0 2.53e-04    -  1.00e+00 1.00e+00h  1
   9  4.4688522e-01 9.34e-09 2.20e+00  -1.0 9.32e-05    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  2.9995484e-02 3.46e-01 1.19e+00  -1.0 5.68e-01    -  5.45e-01 1.00e+00f  1
  11  3.9964883e-04 3.56e-02 5.81e-02  -1.0 2.09e-01    -  1.00e+00 1.00e+00h  1
  12  3.7581923e-06 5.63e-04 5.70e-03  -2.5 2.69e-02    -  1.00e+00 1.00e+00h  1
  13  7.6573980e-10 1.48e-06 3.13e-05  -3.8 1.11e-03    -  1.00e+00 1.00e+00h  1
  14  4.8778322e-14 2.53e-10 7.52e-09  -5.7 1.29e-05    -  1.00e+00 1.00e+00h  1
  15  8.8029044e-20 1.60e-14 6.47e-13  -8.6 9.88e-08    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 15

                                   (scaled)                 (unscaled)
Objective...............:   8.8029044107981477e-20    8.8029044107981477e-20
Dual infeasibility......:   6.4749359007033628e-13    6.4749359007033628e-13
Constraint violation....:   1.5987211554602254e-14    1.5987211554602254e-14
Complementarity.........:   2.5060006300474834e-09    2.5060006300474834e-09
Overall NLP error.......:   2.5060006300474834e-09    2.5060006300474834e-09


Number of objective function evaluations             = 16
Number of objective gradient evaluations             = 16
Number of equality constraint evaluations            = 16
Number of inequality constraint evaluations          = 16
Number of equality constraint Jacobian evaluations   = 16
Number of inequality constraint Jacobian evaluations = 16
Number of Lagrangian Hessian evaluations             = 15
Total CPU secs in IPOPT (w/o function evaluations)   =      0.005
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
               t_proc [s]   t_wall [s]    n_eval
       nlp_f      2.3e-05     2.12e-05        16
       nlp_g      5.9e-05     5.44e-05        16
  nlp_grad_f      4.2e-05     3.99e-05        17
  nlp_hess_l      6.7e-05     6.59e-05        15
   nlp_jac_g      5.5e-05     5.68e-05        17
      solver      0.00718      0.00718         1

The call will fail with an error if the solver fails to convergence. You may still inspect the non-converged solution (see Section 9.3).

You may call solve any number of times. You will always get an immutable copy of the problem specification and its solution. Consecutively calling solve will not help the convergence of the problem.

To warm start a solver, you need to explicitly transfer the solution of one problem to the initial value of the next.

sol1 = opti.solve()
print(sol1.stats()["iter_count"])

# Solving again makes no difference
sol1 = opti.solve()
print(sol1.stats()["iter_count"])

# Passing initial makes a difference
opti.set_initial(sol1.value_variables())
sol2 = opti.solve()
print(sol2.stats()["iter_count"])

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

15
15
4
sol1 = opti.solve();
sol1.stats.iter_count

% Solving again makes no difference
sol1 = opti.solve();
sol1.stats.iter_count

% Passing initial makes a difference
opti.set_initial(sol1.value_variables());
sol2 = opti.solve();
sol2.stats.iter_count

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

ans =  15
ans =  15
ans =  4

In order to initialize the dual variables, e.g. when solving a set of similar optimization problems, you can use the following syntax:

sol = opti.solve()
lam_g0 = sol.value(opti.lam_g)
opti.set_initial(opti.lam_g, lam_g0)

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

sol = opti.solve();
lam_g0 = sol.value(opti.lam_g);
opti.set_initial(opti.lam_g, lam_g0);

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

9.2.2. Numerical value at the solution

Afterwards, you may retrieve the numerical values of variables (or expressions of those variables) at the solution:

  • sol.value(x): value of a decision variable
  • sol.value(p): value of a parameter
  • sol.value(sin(x+p)): value of an expression
  • sol.value(jacobian(opti.g,opti.x)): value of constraint jacobian

Note that the return type of value is sparse when applicable.

9.2.3. Numerical value at other points

You may pass a list of overruling assignment expressions to value. In the following code, we are asking for the value of the objective, using all optimal values at the solution, except for y, which we set equal to 2. Note that such statement does not modify the actual optimal value of y in a permanent way.

obj = (y-x**2)**2

opti.minimize(obj)


print(sol.value(obj,[y==2]))

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

1.909830056703819
obj = (y-x^2)^2;

opti.minimize(obj);



sol.value(obj,{y==2})

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

ans =  1.9098

A related usage pattern is to evaluate an expression at the initial guess:

print(sol.value(x**2+y,opti.initial()))

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

0.0
sol.value(x**2+y,opti.initial())

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

ans = 0

9.2.4. Dual variables

In order to obtain dual variables (Lagrange multipliers) of constraints, make sure you save the constraint expression first:

con = sin(x+y)>=1
opti.subject_to(con)
sol = opti.solve()

print(sol.value(opti.dual(con)))

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

0.258679501736367
con = sin(x+y)>=1;
opti.subject_to(con);
sol = opti.solve();

sol.value(opti.dual(con))

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

ans =  0.25868

9.3. Extras

It may well happen that the solver does not find an optimal solution. In such cases, you may still access the non-converged solution through debug mode:

opti.debug.value(x)

Related, you may inspect the value of an expression, at the initial guess that you supplied to the solver:

opti.debug.value(x,opti.initial())

In case the solver stops due to problem infeasibility, you may identify the problematic constraints with:

opti.debug.show_infeasibilities()

In case the solver reports NaN/Inf at a certain location, you may find out which constraint or variable is to blame by looking at its description:

opti.debug.x_describe(index)
opti.debug.g_describe(index)

You may specify a callback function; it will be called at each iteration of the solver, with the current iteration number as argument. To plot the progress of the solver, you may access the non-converged solution through debug mode:

opti.callback(lambda i: plot(opti.debug.value(x)))
opti.callback(@(i) plot(opti.debug.value(x)))

The callback may be cleared from the Opti stack by calling the Callback function without arguments.

10. Difference in usage from different languages

10.1. General usage

Example:

Table 10.1.1 General usage
Task Python C++ MATLAB/Octave
Starting CasADi
from casadi import *
#include <casadi/casadi.hpp>
using namespace casadi;
import casadi.*
Printing object
print(A)
std::cout << A;
disp(A)
Printing with type information A <ENTER> (interactive), print(repr(A))
std::cout << repr(A);
A <ENTER> (interactive), disp(repr(A))
Get (extended) representation, more=false by default A.str(more) str(A, more) str(A, more)
Calling a class function
SX.zeros(3,4)
SX::zeros(3,4)
SX.zeros(3,4)
Creating a dictionary (e.g. for options)
d = dict(opt1=opt1)
#or
d = {'opt1':opt1}
#or
d = {}; a['opt1'] = opt1
Dict a;
a["opt1"] = opt1;
a = struct;
a.opt1 = opt1;

#or
a = struct('opt1',opt1);
Creating a symbol
MX.sym("x",2,2)
MX::sym("x",2,2)
MX.sym('x',2,2)
Creating a function
Function("f",[x,y],[x+y])
Function("f",{x,y},{x+y})
Function('f',{x,y},{x+y})
Calling a function (form 1)
z=f(x,y)
z = f({x,y})
z=f(x,y)
Calling a function (form 2)
res = f(x=x,y=y)
res["z"]
auto res =\
f({{"x",x},{"y",y}});
res["z"]
res=f('x',x,'y',y)
res.z
Create an NLP solver
nlp = {"x":x,"f":f}
nlpsol("S","ipopt",nlp)
MXDict nlp =\
{{"x",x},{"f",f}};
nlpsol("S","ipopt",nlp);
nlp=struct('x',x,'f',f);
nlpsol('S','ipopt',nlp);

10.2. List of operations

The following is a list of the most important operations. Operations that differ between the different languages are marked with a star (*). This list is neither complete, nor does it show all the variants of each operation. Further information is available in the API documentation.

Table 10.2.1 General usage
Task Python C++ MATLAB/Octave
Addition, subtraction
x+y, x-y, -x
x+y, x-y, -x
x+y, x-y, -x
Elementwise multiplication, division
x*y, x/y
x*y, x/y
x.*y, x./y
Matrix multiplication
x @ y
# or before Py3.4:
mtimes(x,y)
mtimes(x,y)
x*y
Linear system solve
solve(A,b)
solve(A,b)
A\b
Natural exponential function and logarithm
exp(x), log(x)
exp(x), log(x)
exp(x), log(x)
Exponentiation
x**y
pow(x,y)
x^y, x.^y
Square root
sqrt(x)
sqrt(x)
sqrt(x)
Trigonometric functions
sin(x), cos(x), tan(x)
sin(x), cos(x), tan(x)
sin(x), cos(x), tan(x)
Inverse trigonometric
asin(x), acos(x)
asin(x), acos(x)
asin(x), acos(x)
Two argument arctangent
atan2(x, y)
atan2(x, y)
atan2(x, y)
Hyperbolic functions
sinh(x), cosh(x), tanh(x)
sinh(x), cosh(x), tanh(x)
sinh(x), cosh(x), tanh(x)
Inverse hyperbolic
asinh(x), acosh(x)
asinh(x), acosh(x)
asinh(x), acosh(x)
Inequalities
a<b, a<=b, a>b, a>=b
a<b, a<=b, a>b, a>=b
a<b, a<=b, a>b, a>=b
(Not) equal to
a==b, a!=b
a==b, a!=b
a==b, a~=b
Logical and
logic_and(a, b)
a && b
a & b
Logical or
logic_or(a, b)
a || b
a | b
Logical not
logic_not(a)
!a
~a
Round to integer
floor(x), ceil(x)
floor(x), ceil(x)
floor(x), ceil(x)
Modulus after division
fmod(x, y)
fmod(x, y)
mod(x, y)
Modulus after division
fabs(x)
fabs(x)
abs(x)
Norm
norm_1(x)
norm_1(x)
norm(x,1)
Sign function
sign(x)
sign(x)
sign(x)
(Inverse) error function
erf(x), erfinv(x)
erf(x), erfinv(x)
erf(x), erfinv(x)
Elementwise min and max
fmin(x, y), fmax(x, y)
fmin(x, y), fmax(x, y)
min(x, y), max(x, y)
Global min and max
mmin(x), mmax(x)
mmin(x), mmax(x)
mmin(x), mmax(x)
Index of first nonzero
find(x)
find(x)
find(x)
If-then-else
if_else(c, x, y)
if_else(c, x, y)
if_else(c, x, y)
Transpose
A.T
A.T()
A',A.'
Inner product
dot(x, y)
dot(x, y)
dot(x, y)
Horizontal/vertical concatenation
horzcat(x, y)
vertcat(x, y)
hcat([x, y])
vcat([x, y])
horzcat(x, y)
vertcat(x, y)
MX::horzcat({x,y})
MX::vertcat({x,y})
[x,y]
[x; y]
c = {x,y};
[c{:}]
vertcat(c{:})
Horizontal/vertical split (inverse of concatenation)
vertsplit(x)
horzsplit(x)
vertsplit(x)
horzsplit(x)
vertsplit(x)
horzsplit(x)
Element access
# 0-based
A[i,j]
A[i]
// 0-based
A(i,j)
A(i)
%  1-based
A(i,j)
A(i)
Element assignment
# 0-based
A[i,j] = b
A[i] = b
// 0-based
A(i,j) = b
A(i) = b
%  1-based
A(i,j) = b
A(i) = b
Nonzero access
# 0-based
A.nz[k]
// 0-based
A.nz(k)
(currently unsupported)
Nonzero assignment
# 0-based
A.nz[k] = b
// 0-based
A.nz(k) = b
(currently unsupported)
Project to a different sparsity
project(x, s)
project(x, s)
project(x, s)