NMPC with CasADi and Python – Part 1: ODE and steady state

CasADi is a powerful open-source tool for nonlinear optimization. It can be used with MATLAB/Octave, Python, or C++, with the bulk of the available resources referencing the former two options. This post series is intended to show a possible method of developing a simulation for an example system controlled by Nonlinear Model Predictive Control (NMPC) using CasADi and Python.

In this post, a file describing the system equations and a script to determine a steady-state setpoint will be developed. This older post contains similar code for CasADi inside MATLAB.

Before using it, the latest version of CasADi for Python should be installed using the official installation instructions. Afterwards, it can be loaded in Python with

import casadi

Tutorials often use

from casadi import *

which leads to shorter code but in some cases might possibly create problems due to names conflicting with other modules.

We first want to create a file containing the system equations. It is usually a good idea to transform the system equations, which may be known as a set of ordinary differential equations or may yet have to be derived, into a state-space representation such that the state of the system is completely described by the state vector . The dynamics of the system can then be described with the vector function or such that

This is called the state equation. In this series, we ignore the output equation

and assume that all states of the system are known at all times.

We now create a file ode.py with a function ode(t, x, u) that contains our state equation. The content of ode.py obviously depends on the specific system, but an example set of equations for a system with four states shall be provided here:

def ode(t, x, u):
    dx0 = casadi.tan(x[3]) * x[1] + 2*u[1] - 1
    dx1 = x[0] - casadi.sin(x[1])
    dx2 = 13 * x[3] + u[0] + 1
    dx3 = 2*x[0]

On each content line, the derivative of a specific state is defined. The return value of the function must be a vector and cannot simply be a Python list like in

    return [dx0, dx1, dx2, dx3]

Instead, the vertcat function can be used for vertical concatenation:

    return casadi.vertcat(dx0, dx1, dx2, dx3)

The content of ode.py should now be similar to the following:

import casadi

def ode(t, x, u):
    dx0 = casadi.tan(x[3]) * x[1] + 2*u[1] - 1
    dx1 = x[0] - casadi.sin(x[1])
    dx2 = 13 * x[3] + u[0] + 1
    dx3 = 2*x[0]
    
    return casadi.vertcat(dx0, dx1, dx2, dx3)

If we now were to call the function ode from the python console with arbitrary arguments, we would find that it returns a CasADi-specific datatype:

>>> from ode import ode
>>> t=0
>>> x=[0.1, 0.2, 0.3, 0.4]
>>> u=[0.5, 0.6]
>>> ode(t, x, u)
DM([0.284559, -0.0986693, 6.7, 0.2])

Now, let’s try and create a file steadystate.py to determine a point where our system is in a steady state. First, we import CasADi and the ode function:

import casadi
from ode import ode

Then, we go on to create an instance of the Opti class. Opti is a helper framework provided by CasADi that simplifies a lot of things when you are working with nonlinear programming.

ocp = casadi.Opti()

Next, we define symbolic variables for our states and control inputs. We have four states and two control inputs:

# define decision variables
nx = 4
nu = 2
X = ocp.variable(nx,1)
U = ocp.variable(nu,1)

Note that X and U do not actually contain numeric values but serve more as a placeholder for defining equations. Next, we define the constraints of our problem: we want to find a steady-state solution, therefore all four state derivatives should zero. This can be expressed pretty elegantly:

# define steady−state problem
ocp.subject_to( ode(0,X,U) == 0)

That’s it for the problem definition! Now all that’s left is telling CasADi which solver it should use (ipopt seems to be a pretty good default choice if you don’t know any better), telling it to solve the problem, and extracting the values we want to know:

# solve
ocp.solver('ipopt')
SS_sol = ocp.solve()

x_SS = SS_sol.value(X)
u_SS = SS_sol.value(U)

print("Steady state solution: \n x={} \n u={}".format(x_SS, u_SS))

The complete file steadystate.py should now look like this:

import casadi
from ode import ode

ocp = casadi.Opti()

# define decision variables
nx = 4
nu = 2
X = ocp.variable(nx,1)
U = ocp.variable(nu,1)

# define steady−state problem
ocp.subject_to( ode(0,X,U) == 0)

# solve
ocp.solver('ipopt')
SS_sol = ocp.solve()

x_SS = SS_sol.value(X)
u_SS = SS_sol.value(U)

print("Steady state solution: \n x={} \n u={}".format(x_SS, u_SS))

When you run the script, it will hopefully tell you lots of debug information and

EXIT: Optimal Solution Found.

If CasADi is not able to find a solution, it could for example mean that there is an error in the state equations, that CasADi needs your help by providing an initial guess, that there is some numeric problem, or that it is just unable to solve the problem at all. In our case, though, it provides the following result:

Steady state solution: 
 x=[ 0.00000000e+00  4.91924738e-33  0.00000000e+00 -7.64705882e-02] 
 u=[-0.00588235  0.5       ]