Welcome to CasADi’s documentation!¶

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 = {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.

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

#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]]

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]]

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

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

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

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]]

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

[x_0, x_1, x_2, x_3, x_4]
[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\):

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

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

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.

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]]

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

[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]]

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

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

      casadi.DM object with properties:


    r =

      casadi.DM object with properties:

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] =

    casadi.DM object with properties:


  [1,2] =

    casadi.DM object with properties:


}
  scalar structure containing the fields:

    q =

      casadi.DM object with properties:


    r =

      casadi.DM object with properties:

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:

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:

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 controls u, 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 at a number of output times. The control vector u is assumed to be piecewise constant and has the same grid discretization as the output grid.

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.

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,z0,p,u[0],adj_xf[],adj_zf[],adj_qf[])->(xf,zf,qf[0],adj_x0[],adj_z0[],adj_p[],adj_u[]) 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,z0,p,u[0],adj_xf[],adj_zf[],adj_qf[])->(xf,zf,qf[0],adj_x0[],adj_z0[],adj_p[],adj_u[]) IdasInterface

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

4.5. Nonlinear programming¶

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

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

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

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

This is Ipopt version 3.14.11, running with linear solver MUMPS 5.4.1.

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  3.7092062e-66 3.19e-01 2.57e-32  -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
Variable bound 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 seconds in IPOPT                               = 0.007

EXIT: Optimal Solution Found.
           S  :   t_proc      (avg)   t_wall      (avg)    n_eval
       nlp_f  |  59.00us (  5.36us)  14.26us (  1.30us)        11
       nlp_g  | 132.00us ( 12.00us)  29.65us (  2.70us)        11
  nlp_grad_f  |  80.00us (  6.67us)  19.74us (  1.64us)        12
  nlp_hess_l  |  60.00us (  6.00us)  13.88us (  1.39us)        10
   nlp_jac_g  |  67.00us (  5.58us)  16.22us (  1.35us)        12
       total  |  31.09ms ( 31.09ms)   7.77ms (  7.77ms)         1
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 is Ipopt version 3.14.11, running with linear solver MUMPS 5.4.1.

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  3.7092062e-66 3.19e-01 2.57e-32  -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
Variable bound 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 seconds in IPOPT                               = 0.008

EXIT: Optimal Solution Found.
           S  :   t_proc      (avg)   t_wall      (avg)    n_eval
       nlp_f  |  63.00us (  5.73us)  15.03us (  1.37us)        11
       nlp_g  | 134.00us ( 12.18us)  30.13us (  2.74us)        11
  nlp_grad_f  |  91.00us (  7.58us)  21.51us (  1.79us)        12
  nlp_hess_l  |  62.00us (  6.20us)  15.23us (  1.52us)        10
   nlp_jac_g  |  68.00us (  5.67us)  17.38us (  1.45us)        12
       total  |  34.91ms ( 34.91ms)   8.80ms (  8.80ms)         1
[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), OOQP, OSQP and PROXQP 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.

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

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]

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;

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

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.

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

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

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

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

F = f.fold(N)
print(F)

fold_f:(i0)->(o0) MXFunction

F = f.fold(N);
disp(F)

fold_f:(i0)->(o0) MXFunction

x = SX.sym("x")
f0 = Function("f0",[x],[sin(x)])
f1 = Function("f1",[x],[cos(x)])
f2 = Function("f2",[x],[tan(x)])
f_cond = Function.conditional('f_cond', [f0, f1], f2)
print(f_cond)

f_cond:(i0,i1)->(o0) Switch

x = SX.sym('x');
f0 = Function('f0',{x},{sin(x)});
f1 = Function('f1',{x},{cos(x)});
f2 = Function('f2',{x},{tan(x)});
f_cond = Function.conditional('f_cond', {f0, f1}, f2);
disp(f_cond);

f_cond:(i0,i1)->(o0) Switch

x = SX.sym("x")
f_true = Function("f_true",[x],[cos(x)])
f_false = Function("f_false",[x],[sin(x)])
f_cond = Function.if_else('f_cond', f_true, f_false)
print(f_cond)

f_cond:(i0,i1)->(o0) Switch

x = SX.sym('x');
f_true = Function('f_true',{x},{cos(x)});
f_false = Function('f_false',{x},{sin(x)});
f_cond = Function.if_else('f_cond', f_true, f_false);
disp(f_cond);

f_cond:(i0,i1)->(o0) Switch

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 3.6.4.
 *  It consists of: 
 *   1) content generated by CasADi runtime: not copyrighted
 *   2) template code copied from CasADi source: permissively licensed (MIT-0)
 *   3) user code: owned by the user
 *
 */
#ifdef __cplusplus
extern "C" {
#endif

/* How to prefix internal symbols */
#ifdef CASADI_CODEGEN_PREFIX
  #define CASADI_NAMESPACE_CONCAT(NS, ID) _CASADI_NAMESPACE_CONCAT(NS, ID)
  #define _CASADI_NAMESPACE_CONCAT(NS, ID) NS ## ID
  #define CASADI_PREFIX(ID) CASADI_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

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

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};

/* 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, int 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, int mem){
  return casadi_f0(arg, res, iw, w, mem);
}

CASADI_SYMBOL_EXPORT int f_alloc_mem(void) {
  return 0;
}

CASADI_SYMBOL_EXPORT int f_init_mem(int mem) {
  return 0;
}

CASADI_SYMBOL_EXPORT void f_free_mem(int mem) {
}

CASADI_SYMBOL_EXPORT int f_checkout(void) {
  return 0;
}

CASADI_SYMBOL_EXPORT void f_release(int 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 casadi_real f_default_in(casadi_int i) {
  switch (i) {
    default: return 0;
  }
}

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 3.6.4.
 *  It consists of: 
 *   1) content generated by CasADi runtime: not copyrighted
 *   2) template code copied from CasADi source: permissively licensed (MIT-0)
 *   3) user code: owned by the user
 *
 */
#ifdef __cplusplus
extern "C" {
#endif

/* How to prefix internal symbols */
#ifdef CASADI_CODEGEN_PREFIX
  #define CASADI_NAMESPACE_CONCAT(NS, ID) _CASADI_NAMESPACE_CONCAT(NS, ID)
  #define _CASADI_NAMESPACE_CONCAT(NS, ID) NS ## ID
  #define CASADI_PREFIX(ID) CASADI_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

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

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};

/* 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, int 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, int mem){
  return casadi_f0(arg, res, iw, w, mem);
}

CASADI_SYMBOL_EXPORT int f_alloc_mem(void) {
  return 0;
}

CASADI_SYMBOL_EXPORT int f_init_mem(int mem) {
  return 0;
}

CASADI_SYMBOL_EXPORT void f_free_mem(int mem) {
}

CASADI_SYMBOL_EXPORT int f_checkout(void) {
  return 0;
}

CASADI_SYMBOL_EXPORT void f_release(int 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 casadi_real f_default_in(casadi_int i) {
  switch (i) {
    default: return 0;
  }
}

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:

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.

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

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

0.00159265

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

0.00159265

C = Importer('gen.c','shell')
f = external('f',C)
print(f(3.14))

0.00159265

C = Importer('gen.c','shell');
f = external('f',C);
disp(f(3.14))

0.00159265
CasADi - 2023-11-07 01:51:05 WARNING("Failed to remove ./tmp_casadi_compiler_shellP0YrAk.so") [.../casadi/solvers/shell_compiler.cpp:66]
CasADi - 2023-11-07 01:51:05 WARNING("Failed to remove ./tmp_casadi_compiler_shellP0YrAk.o") [.../casadi/solvers/shell_compiler.cpp:67]

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

disp(gen('f', 3.14))

1.5927e-03

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

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

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

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

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

void* fname_alloc_mem(void);

int fname_init_mem(void* mem);

int fname_free_mem(void* mem);

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

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

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 work with finite differences (enable_fd), 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. Also see the callback.py example.

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

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]

# 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}

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

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

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

#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};
  }
};

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.

/*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.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'],\
      {"compiler":"shell"})
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'},...
      struct('compiler','shell'));
disp(f)

f:(x,s)->(r) JitFunction
CasADi - 2023-11-07 01:50:52 WARNING("Failed to remove jit_tmpuFVG1B.c") [.../casadi/core/function_internal.cpp:128]
CasADi - 2023-11-07 01:50:52 WARNING("Failed to remove ./tmp_casadi_compiler_shellCEFXwL.so") [.../casadi/solvers/shell_compiler.cpp:66]
CasADi - 2023-11-07 01:50:52 WARNING("Failed to remove ./tmp_casadi_compiler_shellCEFXwL.o") [.../casadi/solvers/shell_compiler.cpp:67]

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],\
     {"compiler":"shell"})
print(f)

f:(x,s)->(r) JitFunction

sp = Sparsity.scalar();
f = Function.jit('f',body,...
     {'x','s'}, {'r'},...
     {sp,sp}, {sp},...
     struct('compiler','shell'));
disp(f)

f:(x,s)->(r) JitFunction
CasADi - 2023-11-07 01:50:54 WARNING("Failed to remove jit_tmp9M1NLU.c") [.../casadi/core/function_internal.cpp:128]
CasADi - 2023-11-07 01:50:54 WARNING("Failed to remove ./tmp_casadi_compiler_shell7aMPMD.so") [.../casadi/solvers/shell_compiler.cpp:66]
CasADi - 2023-11-07 01:50:54 WARNING("Failed to remove ./tmp_casadi_compiler_shell7aMPMD.o") [.../casadi/solvers/shell_compiler.cpp:67]

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.

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

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

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

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