Nonlinear Model Predictive Control – MPC – Tutorial 1: Basic Formulation and Solution in MATLAB – Fusion of Engineering, Control, Coding, Machine Learning, and Science

In this tutorial series, we explain how to formulate and numerically solve different versions of the nonlinear Model Predictive Control (MPC) problem. We implement the solution in MATLAB. The focus of this series is mainly on the practical implementation without going too deep into the theoretical aspects of nonlinear MPC. You are currently reading Tutorial 1. In this tutorial, we explain how to formulate an open-loop nonlinear MPC problem and how to approximate its solution in MATLAB. The YouTube video accompanying this webpage is given below.

Before reading this tutorial and watching the video, we suggest that you go over our previous tutorial on the linear MPC implementation given here. The developed MATLAB code can be found over here (fee needs to be paid). The YouTube tutorial explaining the implementation is given below.

How to Cite this Document and Tutorial:

To cite this document and tutorial:

“Practical Introduction to Nonlinear Model Predictive Control with MATLAB Implementation”. Technical Report, Number 1, Aleksandar Haber, (2024), Publisher: www.aleksandarhaber.com, Link: https://aleksandarhaber.com/nonlinear-model-predictive-control-mpc-tutorial-1-open-loop-formulation-and-solution-in-matlab/.

When citing the report, include the link to this webpage!

Copyright Notices and NOT to Be Used for AI Notice

COPYRIGHT NOTICE: THE TEXT, PHOTOS, AND CODES POSTED ON THIS WEBSITE AND ON THIS WEBPAGE ARE THE OWNERSHIP, INTELLECTUAL PROPERTY, AND COPYRIGHTED BY THE AUTHOR: ALEKSANDAR HABER. THE TEXT AND THE CONTENT OF THIS PAGE SHOULD NOT BE PHYSICALLY COPIED, SHOULD NOT BE COPIED ON OTHER WEBSITES, SHOULD NOT BE REPRINTED, SHOULD NOT BE REPOSTED ON OTHER WEBSITES, SHOULD NOT BE USED AS A LECTURE MATERIAL IN UNIVERSITY COURSES, SHOULD NOT BE USED AS LECTURE MATERIAL IN COURSES ORGANIZED BY AND HOSTED ON ONLINE LEARNING PLATFORMS (such as Udemy, Coursera, etc), AND SHOULD NOT BE USED IN COMMERCIAL SETTING. THE TEXT, PHOTOS, AND CODE PRESENTED ON THIS WEBSITE AND WEBPAGE SHOULD NOT BE INCLUDED IN REPORTS, STUDENT PAPERS, SCIENTIFIC PAPERS, OR IN ANY OTHER PRINTED OR A DIGITAL FORM OR A DOCUMENT WITHOUT WRITTEN CONSENT OF THE AUTHOR. A MONEY FEE MIGHT BE REQUIRED TO REPRINT THE MATERIAL POSTED ON THIS WEBSITE: CONTACT THE AUTHOR: ml.mecheng@gmail.com

CODE COPYRIGHT NOTICE AND LICENSE: THE CODE FILES POSTED ON THIS WEBSITE ARE NOT FREE SOFTWARE AND CODE. IF YOU WANT TO USE THIS CODE IN THE COMMERCIAL SETTING OR ACADEMIC SETTING, THAT IS, IF YOU WORK FOR A COMPANY OR IF YOU ARE AN INDEPENDENT CONSULTANT AND IF YOU WANT TO USE THIS CODE OR IF YOU ARE ACADEMIC RESEARCHER OR STUDENT, THEN WITHOUT MY PERMISSION AND WITHOUT PAYING THE PROPER FEE, YOU ARE NOT ALLOWED TO USE THIS CODE. YOU CAN CONTACT ME AT
ml.mecheng@gmail.com
TO INFORM YOURSELF ABOUT THE LICENSE OPTIONS AND FEES FOR USING THIS CODE. ALSO, IT IS NOT ALLOWED TO (1) MODIFY THIS CODE IN ANY WAY WITHOUT MY PERMISSION. (2) INTEGRATE THIS CODE IN OTHER PROJECTS WITHOUT MY PERMISSION. (3) POST THIS CODE ON ANY PRIVATE OR PUBLIC WEBSITES OR CODE REPOSITORIES. DELIBERATE OR INDELIBERATE VIOLATIONS OF THIS LICENSE WILL INDUCE LEGAL ACTIONS AND LAWSUITS.

NOT TO BE USED FOR AI COPYRIGHT NOTICE: The code and text, as well as all other material on this webpage and the YouTube page, should NOT be used to train an AI algorithm or a large language model, or any type of AI or machine learning algorithm used to recognize, interpret, and draw conclusions from text. Also, it is forbidden to crawl this webpage and to extract information for training an AI algorithm of any sort on the basis of the material presented on this webpage.

Continuous-time Formulation of the Open-Loop Nonlinear Model Predictive Control Problem

Here, for presentation clarity and completeness, we first briefly formulate a continuous-time open-loop MPC problem. We focus on the open-loop control problem. The open-loop means that we predict and compute the complete future control input sequence at the initial time step, without shifting the time step and without observing the state at the shifted time step. That is, we assume that once it is computed at the initial time step, the complete future control sequence is applied at the subsequent time steps. This strategy can easily be modified to develop a closed-loop MPC algorithm (see our previous tutorial given here).

Here it should be kept in mind that the continuous-time formulation is rarely used in practice. The issue is that this problem is infinite-dimensional, and usually, it cannot be solved directly. Instead, this problem is usually discretized, and additional assumptions are introduced before it can be solved. This will be explained in the next section.

We consider a nonlinear dynamical system in the general form (we will introduce examples later on):

(1) $\begin{align*}\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u})\end{align*}$

where

$\mathbf{x} = \mathbf{x}(t) \in \mathbb{R}^{n}$ is the n-dimensional state vector.
$\mathbf{u} = \mathbf{u}(t) \in \mathbb{R}^{m}$ is the m-dimensional control input vector.
$\mathbf{f}$ is a nonlinear state function that maps states and inputs into the n-dimensional derivative vector of states.

The goal of the MPC algorithm is to track a reference (desired) state trajectory. This trajectory should be designed by the user. In typical position control of mechanical systems, this trajectory consists of the user-defined position and velocity profiles. Let the desired reference trajectory be denoted by $\mathbf{x}_{D}(t) \in \mathbb{R}^{n}$ .

Before we formulate the nonlinear MPC problem, we need to introduce the concept of the continuous-time prediction and control horizon denoted by $T$ . The continuous-time prediction and control horizon $T$ denotes the time duration of the time horizon over which we predict the state trajectory as well as formulate and solve the control problem. The prediction and computation are performed at the initial time step (time step 0) of the prediction and control horizon, and the computed control sequence is applied to the actual system during the control horizon.

We assume that the complete state vector at the beginning of the prediction and control horizon is given or directly measured. In the MPC problem formulation given below, we assume that the start time instant of the prediction and control horizon is $0$ . Consequently, we assume that $\mathbf{x}(0)$ is given. The MPC problem is illustrated in Fig. 1 below.

The nonlinear MPC problem for tracking a reference trajectory $\mathbf{x}_{D}(t)$ can be defined as follows.

(2) $\begin{align*}& \min_{u(t)} \int_{0}^{T} \big( \mathbf{x}_{D}(t) - \mathbf{x}(t))^{T} Q(t) \big( \mathbf{x}_{D}(t) - \mathbf{x}(t)) \text{d}t \\& \text{subject to: } \\& \dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u}) \\& \mathbf{u}_{\text{min}} \le \mathbf{u}(t) \le \mathbf{u}_{\text{max}},\;\; 0 \le t \le T\end{align*}$

In (2),

$Q(t) \in\mathbb{R}^{n}$ is the weighting matrix penalizing the state. We can use this matrix to tune the MPC control algorithm.
$\mathbf{u}_{\text{min}}$ and $\mathbf{u}_{\text{max}}$ are the constant and known vectors defining the control input bounds. They encode the physical limitations of the actuators.

Note that the cost function in (2) is an integral over time of the weighted control error (the term $\mathbf{x}_{D}(t) - \mathbf{x}(t)$ is the control error). Also, in this tutorial, for simplicity and brevity, we have only penalize the control inputs and not state trajectories. State constraints can also easily be incorporated into the optimization problem. This will be explained in our future tutorials.

If we would be able to somehow solve the optimization problem (2), then the solution would be a continuous-time control input function

(3) $\begin{align*}\hat{\mathbf{u}}(t,\mathbf{x}(0))\end{align*}$

However, in the most general case, it is practically impossible to compute an analytical form of the control input function given by (3). One of the main difficulties is that the optimization (2) is infinite-dimensional.

The way to solve this problem is to approximate the problem by a finite-dimensional problem. This is explained in the next section.

Discrete-time and Discretized Formulation of the Open-Loop Nonlinear Model Predictive Control Problem

The first step is to introduce a finite-dimensional control input sequence. We assume that instead of searching for a continuous time control input sequence (3), we search for a discrete-time and finite control input sequence

(4) $\begin{align*}\mathbf{u}_{0},\;\; \mathbf{u}_{1}, \;\; \mathbf{u}_{2},\;\; \ldots, \;\; \mathbf{u}_{N-1}\end{align*}$

where

(5) $\begin{align*}\mathbf{u}_{k}=\mathbf{u}(t=kh),\;\; k=0,1,2,3,\ldots, N-1\end{align*}$

and where $h>0$ is a discretization constant, and $k$ is a discrete-time instant. Here, the last time instant $t=(N-1)h$ corresponds to the control horizon $T$ . That is, $t=T=(N-1)h$ . In this way, we split the prediction and control time horizon $[0,T]$ into the $N$ equidistant time instants:

(6) $\begin{align*}t=0, h, 2h, 3h, \ldots, (N-1)h=T\end{align*}$

The positive integer number $N$ is called the discrete-time prediction and control horizon.

Next, we assume that the control inputs are constant between the discrete-time instants $k\cdot h$ and $(k+1)\cdot h$ . That is, we assume the following

(7) $\begin{align*}u(t)=u_{k}=\text{const},\;\; \forall t\in [kh,(k+1)h), \;\; k=0,1,2,\ldots, N-1 \end{align*}$

The control inputs are shown in the figure below.

The next step is to discretize the continuous dynamics (1). Here, we can use explicit or implicit discretization schemes or even the discretization scheme used by MATLAB’s ode45() solver (we actually use the ode45 solver in our MATLAB simulations). For the presentation clarity, we will formally use the discretized dynamics obtained by using an explicit time-forward discretization. By discretizing (1) and under the assumption of the constant input between the time samples, we obtain the discrete-time state equation representing the discretized dynamics:

(8) $\begin{align*}\mathbf{x}_{k+1}=\mathbf{g}(\mathbf{x}_{k},\mathbf{u}_{k})\end{align*}$

where $\mathbf{g}(\cdot)$ is the discretized state-space dynamics, and $\mathbf{x}_{k}$ is the discretized state vector at the time instant $t=k\cdot h$ . That is, $\mathbf{x}_{k}=\mathbf{x}(t=k\cdot h)$ .

Note that the equation (8) is written only for presentation clarity and in order to illustrate the main ideas. In the MATLAB implementation, we will not use this equation, instead, we will use the ode45 solver to obtain a discrete-time state sequence $\mathbf{x}_{k}$ .

Next, we discretize the cost function (2). By approximating the integral by a sum, we obtain the following cost function

(9) $\begin{align*}\min_{\mathbf{u}_{0},\mathbf{u}_{1},\mathbf{u}_{2},\ldots, \mathbf{u}_{N-1} } \sum_{i=1}^{N} h \cdot \big( \mathbf{x}_{D,i} - \mathbf{x}_{i})^{T} Q_{i} \big( \mathbf{x}_{D,i} - \mathbf{x}_{i}) \end{align*}$

In (9):

$\mathbf{x}_{D,i}$ is the discretized reference (desired) trajectory at the time instant $t= i\cdot h$ . That is, $\mathbf{x}_{D,i}=\mathbf{x}_{D}(t=ih)$ .
$Q_{i}$ is the weight function at the time instant $t= i\cdot h$ . That is, $Q_{i}=Q(t=ih)$ .

Finally, we obtain the discrete-time formulation of the nonlinear model predictive control (MPC) problem

(10) $\begin{align*}&\min_{\mathbf{u}_{0},\mathbf{u}_{1},\mathbf{u}_{2},\ldots, \mathbf{u}_{N-1} } \sum_{i=1}^{N} h \cdot \big( \mathbf{x}_{D,i} - \mathbf{x}_{i})^{T} Q_{i} \big( \mathbf{x}_{D,i} - \mathbf{x}_{i}) \\& \text{subject to: } \\& \mathbf{x}_{i+1}=\mathbf{g}(\mathbf{x}_{i},\mathbf{u}_{i}),\;\; 0 \le i \le N-1 \\& \mathbf{u}_{\text{min}} \le \mathbf{u}_{i} \le \mathbf{u}_{\text{max}},\;\; 0 \le i \le N-1\end{align*}$

Let us group all the optimization variables in a single vector

(11) $\begin{align*}\bar{\mathbf{u}}=\begin{bmatrix}\mathbf{u}_{0} \\ \mathbf{u}_{1} \\ \vdots \\ \mathbf{u}_{N-1} \end{bmatrix}\end{align*}$

When implementing the solution of the optimization problem (10), it is important to keep in mind the following

Any optimization solver for the above problem iteratively computes the optimization variable $\bar{\mathbf{u}}$ , starting from some guess of $\bar{\mathbf{u}}$ . When $\bar{\mathbf{u}}$ is fixed, and due to the fact that the initial condition $\mathbf{x}_{0}$ is known, we can simulate the discretized dynamics (8) to obtain the state sequence $\mathbf{x}_{1}, \mathbf{x}_{2},\ldots, \mathbf{x}_{N}$ . That is, in every iteration of the optimization solver, the discrete states in the cost function can be substituted by their numerical values. Since the reference trajectory is known, in every iteration of the optimization solver, we can evaluate the cost function. This is enough for implementing the solution of the optimization problem in MATLAB. We can use an interior point method to solve the problem.
We do not need to know the gradients of the cost function in order to implement the solution in MATLAB. This is explained in the video. However, by computing the gradient, we can significantly speed up the computations. This will be explored in our future tutorials.

Examples

We will test the MPC algorithm by using two examples. The first example is a linear mass-spring-damper system in the state-space form

(12) $\begin{align*}& \dot{x}_{1}=x_{2} \\& \dot{x}_{2} = -\frac{c}{m}x_{1}-\frac{d}{m}x_{2}+\frac{1}{m}u\end{align*}$

where $c$ , $d$ , and $m$ are the model constants.

The second example is a nonlinear system represented by the following state-space model

(13) $\begin{align*}& \dot{x}_{1}=x_{2} \\& \dot{x}_{2} = 2gx_{1}e^{-gx_{1}^{2}}\Big(bx_{1}^{2}+cx_{1}^{3}+dx_{1}^{4} \Big) \\& -e^{-gx_{1}^{2}}\Big(2bx_{1}+3cx_{1}^{2}+4dx_{1}^{3} \Big)-nx_{2}+u \end{align*}$

Control Results – First Example

The first example is defined by the following parameters: $m=1$ , $c=1$ , and $d=1$ .

The figure below shows a simulated phase portrait and a single trajectory.

Fig.3: Phase portrait and a single trajectory of the first example.

The figure below shows the desired and controlled state-space trajectory. The controlled trajectory is controlled by using the solution of the MPC problem.

Fig. 4: Desired and controlled state-space trajectory.

The time representation of the desired value of $x_{1}$ and the controlled value of $x_{1}$ is shown below.

Fig. 5: Desired and controlled $x_{1}$ .

Control Results – Second Example

The second example is defined by the following parameters: $g=1$ , $b=-1$ , $c=-0.1$ , $d=0.5$ , and $n=0.1$ .

The phase portrait is given below and a single trajectory is given in the figure below.

Fig.6: Phase portrait and a single trajectory of the first example.

The figure below shows the desired and controlled state-space trajectory. The controlled trajectory is controlled by using the solution of the MPC problem.

Fig. 7: Desired and controlled state-space trajectory.

The time representation of the desired value of $x_{1}$ and the controlled value of $x_{1}$ is shown below.

How to Cite this Document and Tutorial:

Copyright Notices and NOT to Be Used for AI Notice

Continuous-time Formulation of the Open-Loop Nonlinear Model Predictive Control Problem

Discrete-time and Discretized Formulation of the Open-Loop Nonlinear Model Predictive Control Problem

Examples

Control Results – First Example

Control Results – Second Example

admin

You might also like

Lyapunov Equation: How to Compute the Solution and Origins in Stability Analysis of Dynamical Systems

Compute Fourier Series in Python by Using Symbolic Library and Generate Plots of Approximation Functions

Control Systems Lecture: Laplace Transform