Skip to content
XeveAbout code and problems

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

April 17, 2020 2 comments Article Uncategorized nspo

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 \vec{x}. The dynamics of the system can then be described with the vector function \vec{\mathrm{f}}(t,x,u) or \vec{\mathrm{ode}}(t,x,u) such that

\dot{\vec{x}}(t) = \vec{\mathrm{ode}}(t,x(t),u(t))

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

\vec{\mathrm{y}}(t) = \vec{\mathrm{h}}(t,x(t),u(t))

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       ]

Tags: CasADi, nlp, nmpc, nonlinear programming, optimal control, python, steady-state

2 comments

  • Akhilesh March 16, 2021 at 13:25 Reply

    Where is the cost function this code is optimizing?

    • nspo April 21, 2021 at 20:05 Reply

      Ha, good question! There’s no cost function in this post. This (and the other two parts of the NMPC post series) actually aren’t doing MPC yet – they are just laying some groundwork. I should definitely get back to writing the next post in the series that’s actually doing NMPC…

Leave a Reply Cancel reply

Your email address will not be published. Required fields are marked *

Calendar

April 2020
M T W T F S S
 12345
6789101112
13141516171819
20212223242526
27282930  
« Jul   May »

Archives

  • January 2021
  • May 2020
  • April 2020
  • July 2018
  • May 2018
  • April 2018
  • March 2018
  • September 2015
  • August 2015
  • June 2015
  • March 2015
  • February 2015
  • September 2014
  • March 2013

Categories

  • Uncategorized

Copyright Xeve 2022 | Theme by ThemeinProgress | Proudly powered by WordPress