Tutorial:
How to solve Non Linear Programming (NLP) problems using CasADi and Interactive Python
Vladimiros Minasidis
Outline
$$\begin{aligned}
& \underset{x\in R^n}{\text{minimize}}& &f(x;p) \\
& \text{subject to :}& & c_i(x;p) \leq 0, \, i \in \mathcal{I}\\
& & & c_i(x;p) = 0, \, i \in \mathcal{E}
\end{aligned}$$
IPython: Productive Interactive Computing
IPython: Productive Interactive Computing
- IPython is a very powerfull tool that incorporates the common shell-like usage patterns, while keeping 100% syntactic compatibility with the Python language itself.
- It's designed in a language-agnostic way to facilitate interactive computing in any language.
- IPython kernel supports multi-language integration, letting you for example mix Python code with Cython, R, Octave, and scripting in Bash, Perl or Ruby.
IPython: Productive Interactive Computing
IPython provides a rich architecture for interactive computing with:
- Powerful interactive shells (terminal and Qt-based).
- A browser-based notebook with support for code, text, mathematical expressions, inline plots and other rich media.
- Support for interactive data visualization
- Flexible, embeddable interpreters to load into your own projects.
more information on http://ipython.org/
IPython: Productive Interactive Computing
IPython provides a rich architecture for interactive computing with:
- Powerful interactive shells (terminal and Qt-based).
e.g. to run any command at the system shell, simply prefix it with ! :
!ping www.ntnu.no
A very nice cheatsheet with IPython commands can be found on http://damontallen.github.io/IPython-quick-ref-sheets/
IPython: Productive Interactive Computing
IPython provides a rich architecture for interactive computing with:
- A browser-based notebook with builtin support for code, text, LaTeX, mathematical expressions, inline plots ...
Cell containing Python code:
%pylab inline
from IPython.core.display import HTML
# Generate a HTML5 compatible wav player
def wavPlayer(filepath):
src = """
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<title>Simple Test</title>
</head>
<body>
<audio controls="controls" style="width:600px" >
<source src="%s" type="audio/wav" />
Your browser does not support the audio element.
</audio>
</body>
"""%(filepath)
display(HTML(src))
# Play the file
wavPlayer("files/arribba.wav")
IPython: Productive Interactive Computing
IPython provides a rich architecture for interactive computing with:
- A browser-based notebook with builtin support for code, text, LaTeX, mathematical expressions, inline plots ...
Cell containing text and mathematical expressions:
Simple spectral analysis
An illustration of the Discrete Fourier Transform using windowing, to reveal the frequency content of a sound signal.
$$ X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2 \pi i}{N} k n} \quad \quad k = 0, \dots, N-1 $$
We begin by loading a datafile using SciPy's audio file support:
IPython: Productive Interactive Computing
IPython provides a rich architecture for interactive computing with:
- A browser-based notebook with builtin support for code, text, LaTeX, mathematical expressions, inline plots ...
Let's plot the spectrogramm of the arribba.wav file:
%matplotlib inline
from scipy.io import wavfile
rate, x = wavfile.read('files/arribba.wav')
from matplotlib import pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ax1.plot(x); ax1.set_title('Raw audio signal')
ax2.specgram(x); ax2.set_title('Spectrogram');
IPython: Productive Interactive Computing
IPython provides a rich architecture for interactive computing with:
- Support for interactive data visualization
import plotly.plotly as py
from plotly.graph_objs import *
py.sign_in("BananaJoe","flye8ny912")
trace1 = Scatter(
x=[1,2,3],
y=[5,6,7],
mode='markers',
marker=Marker(
color=['blue','red','yellow'],
size=[10,20,30]
)
)
data = Data([trace1])
# Make Figure object
fig = Figure(data=data)
IPython: Productive Interactive Computing
IPython provides a rich architecture for interactive computing with:
- Support for interactive data visualization
# Send to Plotly and show in notebook
py.iplot(fig, filename='random_stuff')
IPython: Productive Interactive Computing
You can even combine some old legacy code with IPython:
# Importing the neccessary libraries
from numpy import arange
from matplotlib.pyplot import plot
from IPython.html.widgets import interact
# Some IPython magic
%matplotlib inline
%load_ext fortranmagic
Fortran code cell
%%fortran
subroutine f1(x, z)
real, intent(in) :: x
real, intent(out) :: z
z = sin(x)
end subroutine f1
IPython: Productive Interactive Computing
You can even combine some old legacy code with IPython:
If you have a F90 compiler you might see some IPython magic interaction:) )
@interact
def ifort(w=5):
xx = arange(0,10,0.1)
yy = [f1(w*x) for x in xx]
plot(xx,yy)
IPython: Productive Interactive Computing
IPython provides a rich architecture for interactive computing with:
- A lot of amazing free open-source extensions but they could need some dirty fixes to work as promised.
Dont get too confident with those or they might bite :)
from IPython.display import YouTubeVideo
YouTubeVideo("ZbspVpA-Ezk")
- On the positive side: The community is really active and willing to help even with smallest issues.
IPython: Productive Interactive Computing
One of the best extensions written for IPython:
NBconvert -- a tool which allows you to convert an .ipynb notebook document file into various static formats.
The currently supported export formats are:
HTML- renders a full static HTML page.
LaTeX- generates.texLatex file. Can automatically generate PDF files from the.tex.
Python- convert a notebook to an executable Python scripts
And the best one is saved for the last
Reveal.js- generates a Reveal.js HTML slideshow which must be served by an HTTP server.
IPython: Productive Interactive Computing
So again: What is an Interactive Python ?
IPython is a suite of tools whose goal it is to cover the whole scientific workflow, from interactive data analysis to publication.
It’s main focus is Python, but people behind the project are ambitious about including other languages. IPython can already magically whizz data back and forth between Python, R, Octave and Julia.
CasADi
CasADi
Symbolic framework for algorithmic differentiation and numeric optimization
A low-level open-source tool, written in self-contained C++ code, for efficient implementation of algorithms for nonlinear numerical optimization
CasADi efficiently calculates the relevant derivative or ODE/DAE sensitivity information as needed by the NLP solver
Full-featured Python front end, as well as back ends to state-of-the-art codes such as Sundials (CVODES, IDAS and KINSOL), IPOPT, WORHP, SNOPT and KNITRO
Developed by Joel Andersson and Joris Gillis at the Optimization in Engineering Center, OPTEC of the K.U. Leuven under supervision of Moritz Diehl
CasADi
Symbolic framework for algorithmic differentiation and numeric optimization
CasADi: can evaluate exact (up to machine precision) derivatives very efficiently using Algorithmic Differentiation
Say that we want to obtain the $\nabla_x{F(x)}$ where $F: R^n \rightarrow R^m$
The idea is quite simple:
$F(x)$ is decomposed into sequence of elementary operations (EO). Typical EOs in computer software are $x+y$, $x*y$, $sin(x)$ and $x^y$ .
Create a computational graph (aka information flow diagramm) for those elementary operations
Automatic Differentiation tools then follow the graphs and apply
the chain ruleto these elementary operations and evaluate the derivatives.
CasADi
Symbolic framework for algorithmic differentiation and numeric optimization
For example :
Let's say we want to calculate $\frac{\partial{f}}{\partial{x_1}}$ where $f(x_1,x_2) = x_1x_2 + sin(x_1)$. $f$ can be decomposed into the following elementary operations and those can be differentiated using chain rule:
| Elem. operations (EO) | Differentiated EO |
|---|---|
| $w_1 = x_1$ | $w'_1$ |
| $w_2 = x_2$ | $w'_2$ |
| $w_3 = w_1 w_2$ | $w'_3 = w'_1 w_2 + w_1 w'_2$ |
| $w_4 = \sin(w_1)$ | $w'_4 = \cos(w_1)w'_1$ |
| $w_5 = w_3 + w_4$ | $w'_5 = w'_3 + w'_4 $ |
This is how the computational graph looks like:
source http://en.wikipedia.org/wiki/Automatic_differentiation
CasADi
Symbolic framework for algorithmic differentiation and numeric optimization
Since we are interested only in $\frac{\partial{f}}{\partial{x_1}}$ we can start the evaluation of the differentiated elementary operations by setting $w'_1 = 1$ and $w'_2 = 0$
| Elem.operations (EO) | Differentiated EO |
|---|---|
| $w_1 = x_1$ | $w'_1 = 1$ |
| $w_2 = x_2$ | $w'_2 = 0$ |
| $w_3 = w_1 w_2$ | $w'_3 = w'_1 w_2 + w_1 w'_2 = 1 x_2 + x_1 0 =x_2$ |
| $w_4 = \sin(w_1)$ | $w'_4 = \cos(w_1)w'_1 = \cos(x_1) 1$ |
| $w_5 = w_3 + w_4$ | $w'_5 = w'_3 + w'_4 = x_2 + \cos(x_1)$ |
CasADi
Symbolic framework for algorithmic differentiation and numeric optimization
The Wikipedia page doesn't have the computational graph for the $\nabla_x{f(x)}$ but these can be generated also from CasADi.
# Load pregenerated images since we haven't seen any CasADi code yet
display(Image(filename='files/fgraph.png'))
display(Image(filename='files/jgraph.png'))
CasADi
Symbolic framework for algorithmic differentiation and numeric optimization
Two different methods for Algorithmic Differentiation:
- Forward accumulation (aka Forward Sensitivity) - which was used in the example. This method traverses the graph from the start and adds all contributions from the the differentiated elementary functions. Very cheap for calculating $\nabla_x{f(x)}^Tx$ product.
- Reverse accumulation (aka Adjoint Sensitivity)- which traverses the graph backwards. Stores all the structural information and its very cheap for calculating $x^T\nabla_x{f(x)}$ but requires extra storage.
Both methods have approximately the same computational cost as the one for evaluating $f$
CasADi
Symbolic framework for algorithmic differentiation and numeric optimization
CasADi allows you to create symbolic expressions using syntax similar to e.g.
Symbolic Math Toolbox for MATLAB or SymPy
Let's look at an example:
# Import CasADi and numpy libraries
import casadi.casadi as cd
import numpy as np
# Declare a symbolic NxM x variable
N = 3
M = 2
x = cd.SX.sym("x",N,M)
print x
# Create a symbolic function y = f(x)
y = x**2 - 2*x
print y
CasADi
Symbolic framework for algorithmic differentiation and numeric optimization
We can get the jacobian of an expression quite easily:
# Get the Jacobian of y w.r.t x
jac = cd.jacobian(y,x)
print jac
CasADi
Symbolic framework for algorithmic differentiation and numeric optimization
We can perform various scalar operations between symbolic variables and numerical varibles:
# Declare a new symbolic variable and operate on y
z = cd.SX.sym("z",N,M)
print y*cd.sin(z)
# Add a number to the expression y
B3 = cd.DMatrix.ones(N,M)*4 # a dense N-by-M matrix with all numerical fours
print B3
print y + B3
# We can also repeat an expression NxM number of times
y_rep = cd.SX.repmat(y + z,2,1)
print y_rep
CasADi
Symbolic framework for algorithmic differentiation and numeric optimization
The symbolic instance SX has a small overhead and it's ideal for scalar operations but it's also limited to just scalar operations:
In order to do matrix operation, need to change the symbolic class SX to MX:
e.g. $\sum{x_i^2+x_i}$ will be declared using MX symbolics
n = 6s
m = 1
# MX symbolics
X = cd.MX.sym("X",n,m)
Y = cd.sumAll(X**2 + X)
print Y
MX symbolics have a large overhead but are more generic. SX and MX can be mixed to construct functions that are both fast and generic.
CasADi
Symbolic framework for algorithmic differentiation and numeric optimization
In order to access the full power of CasADi symbolic functions, these need to be declared using "SXFunction" or "MXFunction classes" e.g. :
# Declare MX symbolics
X = cd.MX.sym("X",n,m)
Y = cd.sumAll(X**2 + X)
print "Y = ", Y
# Declare an MX function using MXFunction class
F = cd.MXFunction(cd.nlpIn(x=X),cd.nlpOut(f=Y))
F.init()
# Define the evaluating point
X0 = 2*cd.DMatrix.ones(n,m)
# Evaluate MX function
print "F[{0}] = ".format(X0), F([X0])
# Get the jacobian of y wrt x using MX Function object
J = F.jacobian("x","f")
J.init()
# Evaluate the jacobian at x0
J_num = J([X0])
print "J_num = ", J_num
# Calculate the Hessian
H = F.hessian("x","f")
H.init()
print "H = ", H([X0])
CasADi
Symbolic framework for algorithmic differentiation and numeric optimization
CasADi features, among others, a backend to a very popular large scale NPL solver IPOPT.
Parametric NLPs in a CasADi are formulated as follows:
$$\begin{aligned} & \underset{x}{\text{minimize}}& &f(x;p) \\ & \text{subject to :}& & x_{lb} \leq x \leq x_{ub}\\ & & & g_{lb}[i] \leq g_i(x;p) \leq g_{ub}[i], \, \forall i \in \{1,...,k\} \end{aligned}$$ where:
- $x\in R^n$ are decision variables
- $p\in R^m$ are parameters
- Equality constraints can be formulated as $g_{lb}[i] = g_{ub}[i]$ for some $i$
CasADi automatically generates exact Jacobian and exact Hessian information and provides them to the NLP solve if required information.
How to solve an NLP using CasADi in Python
Lets look at a toy example on how to use CasADi and IPOPT:
Mathematical problem formulation:
$$\begin{aligned} & \underset{x}{\text{min}}& &f(x;p) = x_1^2 + x_2^2 + x_3^2 \\ & \text{subject to :}& & x_1, x_2, x_3 \geq 0\\ & & & 6x_1 + 3x_2 + 2x_3 - p_1 = 0\\ & & & p_2x_1 + x_2 - x_3 - 1 = 0 \end{aligned}$$
IPython code:
# Import libraries
import casadi as cd
import numpy as np
import matplotlib.pyplot as plt
# Decision variables
x = cd.SX.sym("x",3)
# Parameters
p = [5.00,1.00]
# Objective function
f = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]
# Concatenate nonlinear constraints
g = cd.vertcat(( \
6*x[0] + 3*x[1] + 2*x[2] - p[0],
p[1]*x[0] + x[1] - x[2] - 1))
# Nonlinear bounds
lbg = [0.00, 0.00]
ubg = [0.00, 0.00]
# Input bounds for the optimization variables
lbx = [0.00, 0.00, 0.00]
ubx = [ cd.inf, cd.inf, cd.inf]
# Initial guess for the decision variables
x0 = [0.15, 0.15, 0.00]
# Create NLP solver
nlp = cd.SXFunction(cd.nlpIn(x=x),cd.nlpOut(f=f, g=g))
solver = cd.NlpSolver("ipopt", nlp)
# Initialize solver
solver.init()
# Pass the bounds and the initial values
solver.setInput( x0, "x0")
solver.setInput(lbx, "lbx")
solver.setInput(ubx, "ubx")
solver.setInput(lbg, "lbg")
solver.setInput(ubg, "ubg")
# Solve NLP
solver.evaluate()
The problem is solved ! Lets retract the information:
# Print the solution
print "----"
print "Minimal cost " , solver.getOutput("f")
print "----"
print "Optimal solution"
print "x = " , solver.output("x").data()
print "----"
Summary
- IPython Notebook is a very powerfull tool. It keeps all your results, figures, and outputs in a single file, in a plain-text based format, which you can put under version-control, email, edit, and view on any platform. After your analysis is done you can run your notebook through a simple tool which will produce a publishable document in a myriad of formats.
- CasADi CasADi efficiently calculates the relevant derivative or ODE/DAE sensitivity information as needed by the NLP solver. Using the full-featured Python front end the time required to implementback ends to state-of-the-art codes such as Sundials (CVODES, IDAS and KINSOL)and IPOPT can be dramatically reduced
Optimal Control problem formulation
$$\begin{equation} min_u\,J = E\,(\,\textbf{x}(t_0),t_0,\textbf{x}(t_f),t_f\,) + \int_{t_0}^{t_f} \mathcal{L}\,(\,\textbf{x}(t),\textbf{u}(t),t\,)\operatorname{dt} \end{equation}$$ $\textit{subject to }:$ $$\begin{align} \dot{\textbf{x}}(t) = &\, \textbf{f}\,(\textbf{x}(t),\textbf{u}(t),t) &\textit{Dynamic constraints}\\ &\textbf{g}(\textbf{x}(t),\textbf{u}(t),t) \geq \textbf{0} &\textit{Algerbaic path constraints}\\ &\boldsymbol{\phi}\,[\,\textbf{x}(t_0),t_0,\textbf{x}(t_f),t_f\,] = 0& \textit{Boundary conditions} \end{align}$$
In order to access the full power of CasADi symbolic functions, these need to be declared using "SXFunction" or "MXFunction classes" e.g. :