Model Predictive Control Basics

Quick

In this text we’ll:

1. Introduction

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

Within the OCP one uses a of the plant to search out the optimal open loop control over the finite horizon considered. Since the model cannot capture the true plant’s behaviour 100%, and since in the actual world the system is subjected to noise and disturbances, one only applies the primary portion of the optimal open loop control, and frequently the state to re-solve the OCP. This closes the loop and creates a .

The maths behind it is comparatively easy and intuitive (especially in comparison to things like robust control) and it is straightforward to code up an MPC controller. Other pros are that it might effectively handle hard and soft constraints on the state and control (hard constraints are strict, whereas soft constraints are enforced via penalisation in the fee function) and it might generally be used on nonlinear systems with nonconvex constraints (depending on how nasty these are after all!). The one con is that you’ll want to solve a optimization problems “online” in real time, which generally is a problem in case you’re controlling a quick system or have limited computing resources.

1.2 Running Example

Throughout the article I’ll consider the double integrator with a zero-order hold control (i.e. a piecewise constant control) because the running example within the code. The continual time system reads:

[
begin{align*}
dot x _1(t) & = x_2(t),
dot x _2(t) & = u(t),
end{align*}
]

with (t in mathbb{R} ) denoting time. Here ( x_1(t) in mathbb{R} ) is the whereas ( x_2(t) in mathbb{R} ) is the You’ll be able to consider this method as a 1kg block sliding on a frictionless table, with (u(t)) the applied force.

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

[
begin{align*}
x_{k+1} &= A x_k + Bu_k,
end{align*}
]

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

[
begin{align*}
A :=
left(
begin{array}{cc}
1 & 0.1
0 & 1
end{array}
right), ,,
B :=
left(
begin{array}{c}
0.05
1
end{array}
right)
end{align*}
]

and (x_k in mathbb{R}^2 ), (u_k in mathbb{R} ). (Note that I take advantage of (x_k) and (u_k) to confer with the discrete-time system’s state and control to make the notation neater in the remaining of the article.)

You need to use the scipy packages’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 next discrete-time optimal control problem (OCP):

[
begin{equation}
{mathrm{OCP}(bar x):}begin{cases}
minlimits_{mathbf{u},mathbf{x}} quad & sum_{k=0}^{K-1} left(x_k^{top}Qx_k + u_k^top R u_k right)+ x_{K}^top Q_K x_{K}
mathrm{subject to:}quad & x_{k+1} = Ax_k + Bu_k, &k in[0:K-1], & dots(1)
quad & x_0 = bar x , & & dots(2)
quad & x_kin [-1,1]times(-infty,infty),& k in[1:K], & dots(3)
quad & u_kin [-1,1],& k in[0:K-1], & dots(4)
end{cases}
end{equation}
]

where,

• (Kinmathbb{R}_{geq 0}) denotes the over which we solve the OCP,
• (kinmathbb{Z}) denotes the discrete time step,
• ([p:q]), with (p,qinmathbb{Z}), denotes the set of integers ({p,p+1,…,q}),
• (bar x in mathbb{R}^2 ) denotes the of the dynamical system,
• (x_kinmathbb{R}^2 ) denotes the at step (k),
• (uinmathbb{R}) denotes the at step (k),
• (Qinmathbb{R}^{2times 2}), (Rinmathbb{R}) and (Q_K in mathbb{R}^{2times 2}) are square positive definite matrices that specify the fee function ((R) is scalar here because (u) is scalar).

Furthermore, we’ll let,

• (mathbf{u}:= (u_0,u_1,…,u_{K−1})inmathbb{R}^K ) denote the
• (mathbf{x}:= (x_0,x_1,…,x_{K})inmathbb{R}^{2(K+1)} ) denote the

To be rigorous, we’ll say that the pair ((mathbf{u}^{*}, mathbf{x}^{*})inmathbb{R}^K times mathbb{R}^{2(K+1)}) is a to ( mathrm{OCP}(bar{x})) provided that it minimises the fee function over all admissible pairs, that’s,

[
begin{equation*}
J(mathbf{u}^{*}, mathbf{x}^{*}) leq J(mathbf{u}, mathbf{x}),quad forall (mathbf{u},mathbf{x})in Omega
end{equation*}
]

where (J:mathbb{R}^K times mathbb{R}^{2(K+1)}rightarrow mathbb{R}_{geq 0} ) is,

[
begin{equation*}
J(mathbf{u},mathbf{x}) :=left( sum_{k=0}^{K-1 }x_k^top Q x_k + u_k^top R u_k right) + x_K^top Q_K x_K
end{equation*}
]

and (Omega ) denotes all admissible pairs,

[
Omega :={(mathbf{u},mathbf{x})in mathbb{R}^{K}times mathbb{R}^{2(K+1)} : (1)-(4),, mathrm{hold}}.
]

Thus, the optimal control problem is to search out a control and state sequence, (( mathbf{u}^{*},mathbf{x}^{*})), that minimises the fee function, subject to the dynamics, (f), in addition to constraints on the state and control, (x_kin[-1,1]times(-infty,infty)), (u_k in [-1,1] ), for all (k). The associated fee function is important to the controller’s performance. Not only within the sense of ensuring the controller behaves well (for instance, to forestall erratic signals) but additionally in specifying the the closed loop state will settle at. More on this in Section 4.

Note how (mathrm{OCP}(bar x) ) is parametrised with respect to the initial state, (bar x ). This comes from the elemental idea behind MPC: that the optimal control problem is iteratively solved with an updated measured state.

2.1 Coding an OCP solver

CasADi’s stack makes it very easy to establish and solve the OCP.

First, some preliminaries:

from casadi import *

n = 2 # state dimension
m = 1 # control dimension
K = 100 # prediction horizon

# an arbitrary initial state
x_bar = np.array([[0.5],[0.5]]) # 2 x 1 vector

# Linear cost matrices (we'll just use identities)
Q = np.array([[1. , 0],
            [0. , 1. ]])
R = np.array([[1]])
Q_K = Q

# Constraints for all k
u_max = 1
x_1_max = 1
x_1_min = -1

Now we define the issue’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 arrange the fee function:

# Specify the initial condition
opti.subject_to(x_tot[:, 0] == x_bar)

cost = 0
for k in range(K):
    # add dynamic constraints
    x_tot_next = get_x_next_linear(x_tot[:, k], u_tot[:, k])
    opti.subject_to(x_tot[:, k+1] == x_tot_next)

    # add to the fee
    cost += mtimes([x_tot[:,k].T, Q, x_tot[:,k]]) +      
                           mtimes([u_tot[:,k].T, R, u_tot[:,k]])

# terminal cost
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)

The code mtimes([x_tot[:,k].T, Q, x_tot[:,k]]) implements matrix multiplication, (x_k^{top} Q x_k ).

We now add the control and state constraints,

# constrain the control
opti.subject_to(opti.bounded(-u_max, u_tot, u_max))

# constrain the position only
opti.subject_to(opti.bounded(x_1_min, x_tot[0,:], x_1_max))

and solve:

# Say we wish to minimise the fee and specify the solver (ipopt)
opts = {"ipopt.print_level": 0, "print_time": 0}
opti.minimize(cost)
opti.solver("ipopt", opts)
    
solution = opti.solve()

# Get solution
x_opt = solution.value(x_tot)
u_opt = solution.value(u_tot)

We are able to plot the answer with the repo’s plot_solution() function.

from MPC_tutorial import plot_solution

plot_solution(x_opt, u_opt.reshape(1,-1)) # must reshape u_opt to (1,K)

3. Model Predictive Control

The answer to ( mathrm{OCP}(bar x) ), ( (mathbf{x}^{*},mathbf{u}^{*}) ), for a given (bar x), provides us with an control, ( mathbf{u}^{*} ). We now by iteratively solving ( mathrm{OCP}(bar x) ) and updating the initial state (that is the MPC algorithm).

[
begin{aligned}
&textbf{Input:} quad mathbf{x}^{mathrm{init}}inmathbb{R}^2
&quad bar x gets mathbf{x}^{mathrm{init}}
&textbf{for } k in [0:infty) textbf{:}
&quad (mathbf{x}^{*}, mathbf{u}^{*}) gets argmin mathrm{OCP}(bar x)
&quad mathrm{apply} u_0^{*} mathrm{ to the system}
&quad bar x gets mathrm{measured state at } k+1
&textbf{end for}
end{aligned}
]

3.1 Coding the MPC algorithm

The remainder is kind of straightforward. First, I’ll put all our previous code in a function:

def solve_OCP(x_bar, K):
    ....

    return x_opt, u_opt

Note that I’ve parametrised the function with the initial state, (bar x), in addition to the prediction horizon, (K). The MPC loop is then implemented with:

x_init = np.array([[0.5],[0.5]]) # 2 x 1 vector
K = 10
number_of_iterations = 150 # must after all be finite!

# matrices of zeros with the proper sizes to store the closed loop
u_cl = np.zeros((1, number_of_iterations))
x_cl = np.zeros((2, number_of_iterations + 1))

x_cl[:, 0] = x_init[:, 0]

x_bar = x_init
for i in range(number_of_iterations):
    _, u_opt = solve_OCP(x_bar, K)
    u_opt_first_element = u_opt[0]

    # save closed loop x and u
    u_cl[:, i] = u_opt_first_element
    x_cl[:, i+1] = np.squeeze(get_x_next_linear(x_bar,   
                                                u_opt_first_element))

    # update initial state
    x_bar = get_x_next_linear(x_bar, u_opt_first_element)

Again, we will plot the closed loop solution.

plot_solution(x_cl, u_cl)

Note that I’ve “measured” the plant’s state by utilizing the function get_x_next_linear(). In other words, I’ve assumed that our model is 100% correct.

Here’s a plot of the closed loop from a bunch of initial states.

4. Further topics

4.1 Stability and Recursive Feasibility

Two of a very powerful elements of an MPC controller are of the iteratively invoked OCP and of the closed loop. In other words, if I actually have solved the OCP at time (k), will there exist an answer to the OCP at time (k+1)? If there exists an answer to the OCP at each time step, will the closed-loop state asymptotically settle at an equilibrium point (i.e. will or not it’s stable)?

Ensuring that an MPC controller exhibits these two properties involves fastidiously designing the fee function and constraints, and selecting a protracted enough prediction horizon. Going back to our example, recall that the matrices in the fee function were simply chosen to be:

[
Q = left( begin{array}{cc}
1 & 0
0 & 1
end{array}
right),, Q_K = left( begin{array}{cc}
1 & 0
0 & 1
end{array}
right),, R = 1.
]

In other words, the OCP penalises the space of the state to the origin and thus drives it there. As you may probably guess, if the prediction horizon, (K), may be very short and the initial state is positioned very near the constraints at (x_1=pm 1), the OCP will find solutions with insufficient “foresight” and the issue shall be infeasible at some future iteration of the MPC loop. (You may also experiment with this by making (K) small within the code.)

4.2 Some Other Topics

MPC is an lively field of research and there are a lot of interesting topics to explore.

What if the total state can’t be measured? This pertains to and . What if I don’t care about asymptotic stability? This (often) has to do with How do I make the controller to noise and disturbances? There are just a few ways to take care of this, with probably being the best-known.

Future articles might give attention to a few of these topics.

5. Further reading

Listed here are some standard and excellent textbooks on MPC.

[1] Grüne, L., & Pannek, J. (2016).

[2] Rawlings, J. B., Mayne, D. Q., & Diehl, M. (2020).

[3] Kouvaritakis, B., & Cannon, M. (2016).