Model Predictive Control Basics

A hands-on tutorial with Python and CasADi

The MPC loop. Image by author.

Quick Summary

In this article we will:

  • Cover the basic ideas.
  • Code up a solver in Python.
  • Play with a simple linear system: the double integrator.
  • Get all the code here: GitHub Repository

1. Introduction

Model predictive control (MPC) is a popular feedback control methodology where a finite-horizon optimal control problem (OCP) is iteratively solved with an updated measured state on each iteration.

The maths behind it is relatively simple and intuitive (especially when compared to things like robust control) and it is easy to code up an MPC controller. Other pros are that it can effectively handle hard and soft constraints on the state and control and it can generally be used on nonlinear systems with nonconvex constraints.

1.2 Running Example

Throughout the article I will consider the double integrator with a zero-order hold control as the running example in the code. The continuous time system reads:

[ \dot{x}_1(t) = x_2(t), \quad \dot{x}_2(t) = u(t), ]

with ( t \in \mathbb{R} ) denoting time. Here ( x_1(t) \in \mathbb{R} ) is the position whereas ( x_2(t) \in \mathbb{R} ) is the velocity. You can think of this system as a 1kg block sliding on a frictionless table, with ( u(t) ) the applied force.

Running example: the double integrator. Image by author.

If we constrain the control to be piecewise constant over intervals of length 0.1 seconds, we get the discrete-time system:

[ x_{k+1} = A x_k + B u_k, ]

with ( k \in \mathbb{Z} ), where,

[ A = \begin{pmatrix} 1 & 0.1 \ 0 & 1 \end{pmatrix}, \quad B = \begin{pmatrix} 0 \ 0.1 \end{pmatrix} ]

and ( x_k \in \mathbb{R}^2, \ u_k \in \mathbb{R}. )

You can use the scipy package's cont2discrete function to get this discrete time system, as follows:

import numpy as np
from scipy.signal import cont2discrete

A = np.array([[0, 1],[0, 0]])
B = np.array([[0],[1]])
C = np.array([[1, 0],[0, 1]])
D = np.array([[0, 0],[0, 0]])
dt = 0.1  # in seconds
discrete_system = cont2discrete((A, B, C, D), dt, method='zoh')
A_discrete, B_discrete, *_ = discrete_system

2. Optimal Control Problem

We’ll consider the following discrete-time optimal control problem (OCP):

[ \text{OCP}(\bar{x}): \begin{cases} \min_{u,x} \sum_{k=0}^{K-1}(x_k^T Q x_k + u_k^T R u_k) + x_K^T Q_K x_K \ \text{s.t. } x_{k+1} = A x_k + B u_k, \ x_0 = \bar{x}, \ \text{for } k \in [0:K-1], \end{cases}]

where,

  • K denotes the finite horizon over which we solve the OCP,
  • x denotes the state at step k,
  • u denotes the control at step k,
  • Q, R, and Q_K are matrices that specify the cost function.

Therefore, the optimal control problem is to find a control and state sequence that minimises the cost function subject to the dynamics, as well as constraints on the state and control.

2.1 Coding an OCP solver

CasADi’s opti stack makes it really easy to set up and solve the OCP.

First, some preliminaries:

from casadi import *

n = 2  # state dimension
m = 1  # control dimension
K = 100  # prediction horizon
x_bar = np.array([[0.5],[0.5]])  # 2 x 1 vector
Q = np.array([[1. , 0], [0. , 1.]])
R = np.array([[1]])
Q_K = Q

Now we define the problem’s decision variables:

opti = Opti()
x_tot = opti.variable(n, K+1)  # State trajectory
u_tot = opti.variable(m, K)    # Control trajectory

Next, we impose the dynamic constraints and set up the cost function:

opti.subject_to(x_tot[:, 0] == x_bar)
cost = 0
for k in range(K):
    x_tot_next = get_x_next_linear(x_tot[:, k], u_tot[:, k])
    opti.subject_to(x_tot[:, k+1] == x_tot_next)
    cost += mtimes([x_tot[:,k].T, Q, x_tot[:,k]]) + \
               mtimes([u_tot[:,k].T, R, u_tot[:,k]])

cost += mtimes([x_tot[:,K].T, Q_K, x_tot[:,K]])

def get_x_next_linear(x, u):
    # Linear system
    A = np.array([[1. , 0.1],[0. , 1.]])
    B = np.array([[0.005],[0.1]])
    return mtimes(A, x) + mtimes(B, u)

Now let's add the control and state constraints:

opti.subject_to(opti.bounded(-u_max, u_tot, u_max))
opti.subject_to(opti.bounded(x_1_min, x_tot[0,:], x_1_max))

And solve:

opts = {