No items found.

Consider a linear time-invariant system given by:

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

where $x_k \in \mathbb{R}^n$ is the state, $u_k \in \mathbb{R}^m$ is the control input, and $A \in \mathbb{R}^{n \times n}$ and $B \in \mathbb{R}^{n \times m}$ are constant matrices. The goal of linear MPC, given a reference trajectory $x_{ref,k}$ and $u_{ref,k}$, is to find a sequence of control inputs $u_k$ that will drive the system to the reference trajectory.

Given a horizon of $N$ steps, at every time $t_i$ when the state (or its estimate) is $x_i$, t

\begin{align}

\text{minimize}_{u_{i+k}} \quad & (x_{i+N} - x_{ref,{i+N}})^T Q_N (x_{i+k} - x_{ref,{i+N}}) + \\

& \sum_{k=0}^{N-1} (x_{i+k} - x_{ref,{i+k}})^T Q (x_{i+k} - x_{ref,{i+k}}) + \\

& \sum_{k=0}^{N-1} (u_{i+k} - u_{ref,{i+k}}) R (u_{i+k} - u_{ref,{i+k}}) \\

\text {subject to} \quad & x_{i+k+1} = A x_{i+k} + B u_{i+k}, \quad k = 0, \cdots, N-1 \\

\quad & x_{i} = x_{i}^{0}, \\

\quad & x_{lb} < x_{i+k} < x_{ub}, \\

\quad & u_{lb} < u_{i+k} < u_{ub},

\end{align}

where $Q_N \in \mathbb{R}^{n \times n}$, $Q \in \mathbb{R}^{n \times n}$, and $R \in \mathbb{R}^{m \times m}$ are positive definite weighting matrices, $(\cdot)_{lb}$ and $(\cdot)_{ub}$ are lower and upper bounds. The first equality constraint is the system dynamics, the second equality constraint is the initial condition, and the last two inequalities are the state and control input constraints, representing admissible values for the state and control input, respectively. The objective function to be minimized is quadratic. The first term represents the terminal cost for deviations of the terminal state from reference terminal state. Similarly, the second and third terms represent the stage cost for deviations of the state and control input from their respective references.

While the solution for $u_{i+k}$ is obtained for $k=[0,1,2,\cdots, N-1]$, only the first control input $u_i$ is applied to the system. The optimization problem is then solved again at the next time step, and the first control input of the new solution is applied to the system. This process is repeated at every time step as the plant's state evolves in time.

For **nonlinear systems** to be controlled by linear MPC, there are two approaches. First, the system may be linearized around an equilibrium point to obtain a linear time-invariant system. Second, the nonlinear system bay be linearized around the time-varying reference trajectory. Presently, only the first approach is implemented in Collimator. Other limitations of the [.code]LinearDiscreteTime[.code] block in Collimator are:

1. $Q_N$ = $Q$ is assumed.

2. bounds on $x$ are not implemented.

3. $u_{lb}$ and $u_{ub}$ are assumed to be the same for all the $m$ components of the control vector.

4. $x_{ref}$ is assumed to be a constant (so the block is presently only suitable for regulation of state as opposed to tracking)

The block expects an [.code]LTISystem[.code] continous plant, and internally converts the plant into a discrete-time representation with Euler discretization.

****

A continuous-time linearized model of a Cessna aircraft at a speed of 128.2 m/s and an altitude of 5,000m is given by the following [1,2]:

\begin{align}

\dot{x} = \begin{bmatrix}-1.2822 & 0 & 0.98 & 0 \\

0 & 0 & 1 & 0 \\

-5.4293 & 0 & -1.8366 & 0 \\

-128.2 & 128.2 & 0 & 0\end{bmatrix}x + \begin{bmatrix}-0.3 \\

0 \\

- 17 \\

0 \end{bmatrix}u

\end{align}

The components of the state vector $x$ are:

$x_1$: the angle of attack (specified in radians): This is the angle between the airplane's longitudinal axis and the aircraft's velocity vector.

$x_2$: the pitch angle (specified in radians). This is the angle between the airplane's longitudinal axis and the horizontal plane.

$x_3$: the pitch rate, representing the rate of change of the pitch angle $x_2$ (specified in radians/sec).

$x_4$: the aircrafts's altitude (specified in meters) relative to the reference altitude of 5,000 m, at which the linearized model is available.

The control variable $u$ is the angle (specified in radians) of the elevator, a controllable surface on the aircraft's tail. The elevator can be changed to control the aircraft's pitch.

Our goal with MPC is to design a controller that regulates the linearized system's altitude to zero (i.e., the aircraft's altitude to 5,000m), subject to the constraint of the control input magnitude not exceeding 0.262 radians.

**References**:

[1] F. Borrelli∗, C. Jones†, M. Morari, Model Predictive Control Part III – Feasibility and Stability, Course notes. (Available online)

[2] L. Lessard, Model Predictive Control, Course notes. (Available online)

Let's first create the matrices for the linearized Cessna model and specify the MPC parameters. We implement the LTI system with full observation so that its output can be used by the MPC.

Let's first observe the MPC controller when the constraint on $u$ is not specified.

While the aircraft's altitute is regulated to zero, we see that the the magnitudes of the applied control input are too large. This is natural as we did not specify any constraints to the MPC. Let's specify them now.

We follow the same procedure as above, except when creating the MPC block, we additionally pass the [.code]lbu[.code] and [.code]ubu[.code] arguments.

The contrainsts are well respected now. The minor violations in the initial two steps can potentially be eliminated by tighter tolerances.

Here, we demonsrate how a linearized version of a nonlinear plant can be used in MPC.

We use a damped pendulum model, with the following dynamics:

\begin{equation}

\dot{x} = \begin{bmatrix}\dot{\theta} \\ \dot{\omega}\end{bmatrix} = \begin{bmatrix}\omega \\

-\dfrac{g}{L} \sin(\theta) - \dfrac{b}{m L^2} \omega + \dfrac{u}{m L^2}\end{bmatrix}.

\end{equation}

The state vector is $x = [\theta, \omega]^T$ and the control input is the torque $u$. The pendulum parameters are mass $m$, length $L$, damping coefficients $b$, and acceleration due to gravity $g$.

We observe that the pendulum is driven to the zero vector state, while not violating the constraint that torque magnitude should not exceed 5 N.m.