\documentclass[12pt]{article}
\usepackage{latexsym}
\renewcommand{\baselinestretch}{1.2}
\setlength{\oddsidemargin}{0in}
\setlength{\evensidemargin}{0in}
\setlength{\textwidth}{6.5in}
\setlength{\topmargin}{0in}
\setlength{\textheight}{8.5in}
%\input setbmp
%\input{psfig.sty}
\usepackage{epsfig}
\newtheorem{defn}{Definition}
\newtheorem{thm}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{example}{Example}
\newtheorem{corollary}{Corollary}
\newtheorem{proposition}{Proposition}
\newtheorem{algorithm}{Algorithm}
\newtheorem{exercise}{Exercise}
%\usepackage[centertags]{amsmath}
%\usepackage{amsfonts}
\newcommand{\markex}{\nopagebreak\begin{center}\rule{5in}{.01in}\end{center}}
\newcommand{\markdefn}{\nopagebreak\hfill$\Diamond$}
\newcommand{\mb}{\mbox}
\newcommand{\bm}{\boldmath}
\def\b#1{\mb{\bm$#1$}}
\newfont{\Bb}{msbm10 scaled\magstep1}
\newfont{\Fr}{eufm10 scaled\magstep1}
\newfont{\Sc}{eusm10 scaled\magstep1}

%-------------
%\usepackage{latexsym}
%\renewcommand{\baselinestretch}{1.55}
%\setlength{\topmargin}{-0.25in}
%\setlength{\textheight}{8.5in}
%\setlength{\oddsidemargin}{0in}
%\setlength{\textwidth}{6.5in}
%\newtheorem{example}{Example}
%\newtheorem{thm}{Theorem}
%\newcommand{\markex}{\nopagebreak\begin{center}\rule{5in}{.02in}\end{center}}
%\newcommand{\mb}{\mbox}
%\newcommand{\bm}{\boldmath}
%\def\b#1{\mb{\bm$#1$}}
%\newfont{\Bb}{msbm10 scaled\magstep1}
%\pagestyle{empty}
\begin{document}
\noindent{\Large ECE 680}\hfill{\Large Fall 2011}
\\
\begin{center}
%{\LARGE \bf Model-based Predictive Control (MPC)\footnote{handout prepared by S.~H.~\.{Z}ak}}
{\LARGE \bf Model-based Predictive Control (MPC)}\\

by\\

{\large Stanislaw H. \.{Z}ak}
\end{center}

\section{Introduction}
The model-based predictive control~(MPC) methodology is also referred to as the moving horizon control or the receding horizon control. The idea behind this approach can be explained using an example of driving a car. The driver looks at the road ahead of him and taking into account the present state and the previous action predicts his action up to some distance ahead, which we refer to as the prediction horizon. Based on the prediction, the driver adjusts the driving direction. The MPC main idea is illustrated in Figure~\ref{fig:MPC_car}.
\begin{figure}[h]
\centering \epsfig{file=MPC_car_2009.eps,width=11cm}
\caption{The driver predicts future travel direction based on the current state of the car and the current position of the steering wheel.} \label{fig:MPC_car}
\end{figure}

The MPC is constructed using control and optimization tools. The objective of this write-up is to introduce the reader to the linear MPC which refers to the family of MPC schemes in which linear models of the controlled objects are used in the control law synthesis.

\noindent\rule[.1in]{3.in}{.01in}\\{\copyright 2009 by Stanislaw H. \.{Z}ak}\\

In the MPC approach, the current control action is computed on-line rather than using a pre-computed, off-line, control law.

A model predictive controller uses, at each sampling instant, the plant's current input and output measurements, the plant's current state, and the plant's model to
\begin{itemize}
\item calculate, over a finite horizon, a future control sequence that optimizes a given performance index and satisfies constraints on the control action;
\item use the first control in the sequence as the plant's input.
\end{itemize}
The MPC strategy is illustrated in Figure~\ref{fig:discrete_MPC}, where $N_p$ is the prediction horizon, $u(t+k|t)$ is the predicted control action at $t+k$ given $u(t)$. Similarly, $y(t+k|t)$ is the predicted output at $t+k$ given $y(t)$.
\begin{figure}[h]
\centering \epsfig{file=discrete_MPC_2009b.eps,width=8cm}
\caption{Controller action construction using model-based predictive control~(MPC) approach.} \label{fig:discrete_MPC}
\end{figure}

\section{Basic Structure of MPC}
In Figure~\ref{fig:MPC_state_f}, we show a basic structure of an MPC-controlled plant, where we assume that the plant's state is available to us.
\begin{figure}[h]
\centering \epsfig{file=MPC_closed_2009.eps,width=12cm}
\caption{State feedback model predictive controller.} \label{fig:MPC_state_f}
\end{figure}



\section{From Continuous to Discrete Models}
Our objective here is to present a method for constructing linear discrete-time models from given linear continuous-time models. The obtained discrete models will be used to perform computations to generate control commands.

We use a sample-and-hold device that transforms a continuous signal, $f(t)$, into the staircase signal,
\[
f(kh),\quad kh\le t< (k+1)h,
\]
where $h$ is the sampling period. The sample and zero-order hold~(ZOH) operation is illustrated in Figure~\ref{fig:sample_hold}
\begin{figure}[h]
\centering \epsfig{file=sample_zoh_2009.eps,width=12cm}
\caption{Sample and zero-order hold~(ZOH) element operating on a continuous function.}
\label{fig:sample_hold}
\end{figure}

Suppose that we are given a continuous-time model,
\begin{eqnarray*}
\dot{\b x}(t) &=& \b A\b x(t) + \b B\b u(t),\quad \b x_0=\b x(t_0)\\
\b y(t) &=& \b C\b x(t).
\end{eqnarray*}
The solution to the state equation is
\begin{equation}
\b x(t)=e^{\b A(t-t_0)} \b x(t_0) +\int_{t_0}^t e^{\b A(t-\tau)} \b B \b u(\tau) d \tau.
\label{eq:contin_model}
\end{equation}
We assume that the input to the system is generated by a sample-and-hold device and has the form,
\[
\b u(t)=\b u(k), \quad kh\le t< (k+1)h.
\]
Let $t_0=kh$ and $t=(k+1)h$ and let us use shorthand notation,
\[
\b x(kh)=\b x(k)
\]
Then taking into account that $\b u(k)$ is constant on the interval $[kh, (k+1)h)$ we represent~(\ref{eq:contin_model}) as
\begin{eqnarray}
\b x(k+1) &=& e^{\b Ah}\b x(k)+\int_{kh}^{(k+1)h} e^{\b A(kh+h-\tau)}\b B\b u(k) d\tau\nonumber\\
&=& e^{\b Ah}\b x(k)+\int_{kh}^{(k+1)h} e^{\b A(kh+h-\tau)}\b B d\tau\, \b u(k).
\label{eq:pre_chang_lim}
\end{eqnarray}
Consider now the second term on the right-hand side of the above equation. Let
\[
\eta=kh+h-\tau.
\]
Then we can represent~(\ref{eq:pre_chang_lim}) as
\begin{eqnarray}
\b x(k+1) &=& e^{\b Ah}\b x(k)+\int_{kh}^{(k+1)h} e^{\b A(kh+h-\tau)}\b B d\tau\, \b u(k)\nonumber\\
&=& e^{\b Ah}\b x(k)+\int_{0}^{h} e^{\b A\eta}\b B d\eta\, \b u(k)\nonumber\\
&=& \b{\Phi} \b x(k) +\b{\Gamma}\b u(k),
\label{eq:final_d_time}
\end{eqnarray}
where
\[
\b{\Phi}=e^{\b Ah}\quad\mb{and}\quad \b{\Gamma}= \int_{0}^{h} e^{\b A\eta}\b B d\eta.
\]
The discrete output equation has the form
\[
\b y(k) = \b C\b x(k).
\]


\section{Simple Discrete-Time MPC}
The development a discrete MPC controller in this section follows that of~\cite[Chapter~1]{wang09}. We consider a discretized model of a dynamic system of the form,
\begin{eqnarray}
\b x(k+1) &=& \b{\Phi} \b x(k) +\b{\Gamma}\b u(k)\label{eq:d_state_mo}\\
\b y(k) &=& \b C\b x(k),\label{d_output_mo}
\end{eqnarray}
where $\b{\Phi}\in \mb{\Bb R}^{n\times n}$, $\b{\Gamma}\in \mb{\Bb R}^{n\times m}$, and $\b{C}\in \mb{\Bb R}^{p\times n}$.

Applying the backward difference operator, $\Delta \b x(k+1)=\b x(k+1)-\b x(k)$, to~(\ref{eq:d_state_mo}) gives
\begin{equation}
\Delta \b x(k+1) = \b{\Phi} \Delta \b x(k) +\b{\Gamma}\Delta \b u(k),\label{eq:diff_state_mo}
\end{equation}
where $\Delta \b u(k+1)=\b u(k+1)-\b u(k)$.

We now apply the backward difference operator to~(\ref{d_output_mo}) to obtain
\begin{eqnarray*}
\Delta \b y(k+1) &=& \b y(k+1)-\b y(k)\\
&=& \b C\b x(k+1) - \b C\b x(k)\\
&=& \b C\Delta \b x(k+1).
\end{eqnarray*}
Substituting  into the above~(\ref{eq:diff_state_mo}) yields
\[
\Delta \b y(k+1) = \b C\b{\Phi} \Delta \b x(k) +\b C\b{\Gamma}\Delta \b u(k).
\]
Hence,
\begin{equation}
\b y(k+1) = \b y(k) + \b C\b{\Phi} \Delta \b x(k) +\b C\b{\Gamma}\Delta \b u(k).
\label{eq:diff_o_output}
\end{equation}
We combine~(\ref{eq:diff_state_mo}) and~(\ref{eq:diff_o_output}) into one equation to obtain
\begin{equation}
\left[\begin{array}{c}
\Delta\b x(k+1)\\
\b y(k+1)
\end{array}\right]=\left[\begin{array}{cc}
\b{\Phi} & \b O\\
\b C\b{\Phi} & \b I_p
\end{array}\right] \left[\begin{array}{c}
\Delta\b x(k)\\
\b y(k)
\end{array}\right] + \left[\begin{array}{c}
\b{\Gamma}\\
\b C\b{\Gamma}
\end{array}\right] \Delta\b u(k).
\label{eq:aug_d_state}
\end{equation}
We represent~(\ref{d_output_mo}) as
\begin{equation}
\b y(k) = \left[\begin{array}{cc}
\b O & \b I_p
\end{array}\right]\left[\begin{array}{c}
\Delta\b x(k)\\
\b y(k)
\end{array}\right].
\label{eq:out_d_aug}
\end{equation}
We now define the augmented state vector,
\begin{equation}
\b x_a(k)=\left[\begin{array}{c}
\Delta\b x(k)\\
\b y(k)
\end{array}\right].
\label{eq:aug_st_vec}
\end{equation}
Let
\begin{equation}
\b{\Phi}_a=\left[\begin{array}{cc}
\b{\Phi} & \b O\\
\b C\b{\Phi} & \b I_p
\end{array}\right] ,\quad \b{\Gamma}_a=\left[\begin{array}{c}
\b{\Gamma}\\
\b C\b{\Gamma}
\end{array}\right],\quad\mb{and}\quad \b C_a=\left[\begin{array}{cc}
\b O & \b I_p
\end{array}\right].
\label{eq:aug_mo_mats}
\end{equation}
Using the above notation, we represent~(\ref{eq:aug_d_state}) and~(\ref{eq:out_d_aug}) in a compact format as
\begin{eqnarray}
\b x_a(k+1) &=& \b{\Phi}_a \b x_a(k) +\b{\Gamma}_a\Delta\b u(k)\label{eq:d_aug_state_mo}\\
\b y(k) &=& \b C_a\b x_a(k),\label{d_aug_output_mo}
\end{eqnarray}
where $\b{\Phi}_a\in \mb{\Bb R}^{(n+p)\times (n+p)}$, $\b{\Gamma}_a\in \mb{\Bb R}^{(n+p)\times m}$, and $\b{C}_a\in \mb{\Bb R}^{p\times (n+p)}$.


Suppose now that the state vector $\b x_a$ at each sampling time, $k$, is available to us. Our control objective is to construct a control sequence,
\begin{equation}
\Delta \b u(k), \Delta \b u(k+1),\ldots , \Delta \b u(k+N_p-1),
\label{eq:d_pred_c_seq}
\end{equation}
where $N_p$ is the prediction horizon, such that a given cost function and constraints are satisfied. The above control sequence will result in a predicted sequence of the state vectors,
\[
\b x_a(k+1|k), \b x_a(k+2|k),\ldots , \b x_a(k+N_p|k),
\]
which can then be used to compute predicted sequence of the plant's outputs,
\begin{equation}
\b y(k+1|k), \b y(k+2|k),\ldots , \b y(k+N_p|k).
\label{eq:pred_outs}
\end{equation}
Using the above information, we can compute the control sequence~(\ref{eq:d_pred_c_seq}) and then apply $\b u(k)$ to the plant to generate $\b x(k+1)$. We repeat the process again, using $\b x(k+1)$ as an initial condition to compute $\b u(k+1)$, until we reach the boundary of the control horizon, $N_c$.

We now present an approach to construct $\b u(k)$ given $\b x(k)$. Using the plant model parameters and the measurement of $\b x_a(k)$ we evaluate the augmented states over the prediction horizon successively applying the recursion formula~(\ref{eq:d_aug_state_mo}) to obtain,
\begin{eqnarray*}
\b x_a(k+1|k) &=& \b{\Phi}_a \b x_a(k) +\b{\Gamma}_a\Delta\b u(k)\\
\b x_a(k+2|k) &=& \b{\Phi}_a \b x_a(k+1|k) +\b{\Gamma}_a\Delta\b u(k+1)\\
&{}& \quad =\b{\Phi}_a^2 \b x_a(k) + \b{\Phi}_a \b{\Gamma}_a\Delta\b u(k) + \b{\Gamma}_a\Delta\b u(k+1)\\
{} &\vdots& {}\\
\b x_a(k+N_p|k) &=& \b{\Phi}_a^{N_p} \b x_a(k) + \b{\Phi}_a^{N_p - 1} \b{\Gamma}_a\Delta\b u(k) + \cdots + \b{\Gamma}_a\Delta\b u(k+N_p-1)
\end{eqnarray*}
We represent the above set of equations in the form,
\begin{equation}
\left[\begin{array}{c}
\b x_a(k+1|k)\\
\b x_a(k+2|k)\\
\vdots\\
\b x_a(k+N_p|k)
\end{array}\right]=\left[\begin{array}{c}
\b{\Phi}_a\\
\b{\Phi}_a^2\\
\vdots\\
\b{\Phi}_a^{N_p}
\end{array}\right]\b x_a(k)+\left[\begin{array}{cccc}
\b{\Gamma}_a & {}& {}& {}\\
\b{\Phi}_a \b{\Gamma}_a & \b{\Gamma}_a &{}& {}\\
\vdots &{}& \ddots &{}\\
\b{\Phi}_a^{N_p - 1} \b{\Gamma}_a &{}& \cdots& \b{\Gamma}_a
\end{array}\right]\left[\begin{array}{c}
\Delta\b u(k)\\
\Delta\b u(k+1)\\
\vdots\\
\Delta\b u(k+N_p-1)
\end{array}\right].
\label{eq:many_eqs_in_one}
\end{equation}
We wish to design a controller that would force the plant output, $\b y$, to track a given reference signal, $\b r$. Using~(\ref{d_aug_output_mo}) and~(\ref{eq:many_eqs_in_one}), we compute the sequence of predicted outputs~(\ref{eq:pred_outs}),
\begin{eqnarray}
&{}&\left[\begin{array}{c}
\b y(k+1|k)\\
\b y(k+2|k)\\
\vdots\\
\b y(k+N_p|k)
\end{array}\right]=\left[\begin{array}{c}
\b C_a\b x_a(k+1|k)\\
\b C_a\b x_a(k+2|k)\\
\vdots\\
\b C_a\b x_a(k+N_p|k)
\end{array}\right]\nonumber\\
&{}&\quad =\left[\begin{array}{c}
\b C_a\b{\Phi}_a\\
\b C_a\b{\Phi}_a^2\\
\vdots\\
\b C_a\b{\Phi}_a^{N_p}
\end{array}\right]\b x_a(k)\\
&{}&\qquad +\left[\begin{array}{cccc}
\b C_a\b{\Gamma}_a & {}& {}& {}\\
\b C_a\b{\Phi}_a \b{\Gamma}_a & \b C_a\b{\Gamma}_a &{}& {}\\
\vdots &{}& \ddots &{}\\
\b C_a\b{\Phi}_a^{N_p - 1} \b{\Gamma}_a &{}& \cdots& \b C_a\b{\Gamma}_a
\end{array}\right]\left[\begin{array}{c}
\Delta\b u(k)\\
\Delta\b u(k+1)\\
\vdots\\
\Delta\b u(k+N_p-1)
\end{array}\right].
\label{eq:out_eqntocom_in}
\end{eqnarray}
We write the above compactly as
\begin{equation}
\b Y=\b W\b x_a(k)+\b Z\Delta \b U,
\label{eq:compact_Y}
\end{equation}
where
\[
\b Y=\left[\begin{array}{c}
\b y(k+1|k)\\
\b y(k+2|k)\\
\vdots\\
\b y(k+N_p|k)
\end{array}\right],\quad \Delta\b U=\left[\begin{array}{c}
\Delta\b u(k)\\
\Delta\b u(k+1)\\
\vdots\\
\Delta\b u(k+N_p-1)
\end{array}\right],
\]
and
\[
\b W=\left[\begin{array}{c}
\b C_a\b{\Phi}_a\\
\b C_a\b{\Phi}_a^2\\
\vdots\\
\b C_a\b{\Phi}_a^{N_p}
\end{array}\right],\quad\mb{and}\quad \b Z=\left[\begin{array}{cccc}
\b C_a\b{\Gamma}_a & {}& {}& {}\\
\b C_a\b{\Phi}_a \b{\Gamma}_a & \b C_a\b{\Gamma}_a &{}& {}\\
\vdots &{}& \ddots &{}\\
\b C_a\b{\Phi}_a^{N_p - 1} \b{\Gamma}_a &{}& \cdots& \b C_a\b{\Gamma}_a
\end{array}\right].
\]

Suppose now that we wish to construct a control sequence, $\Delta \b u(k), \ldots , \Delta \b u(k+N_p-1)$, that would minimize the cost function
\begin{equation}
J(\Delta \b U)=\frac{1}{2}\left(\b r_p-\b Y  \right)^{\top}\b Q \left(\b r_p-\b Y  \right) + \frac{1}{2}\Delta \b U^{\top}\b R\Delta \b U,
\label{eq:cost_function}
\end{equation}
where $\b Q=\b Q^{\top}>0$ and $\b R=\b R^{\top}>0$ are real symmetric positive semi-definite and positive-definite weight matrices, respectively. The multiplying scalar, $1/2$, is just to make subsequent manipulations cleaner. Finally, the vector $\b r_p$ consists of the values of the command signal at sampling times, $k, k+1,\ldots , k+N_p-1$. The selection of the weight matrices, $\b Q$ and $\b R$ reflects our control objective to keep the tracking error $\|\b r_p-\b Y\|$ ``small'' using the control actions that are ``not too large.''

We first apply the first-order necessary condition~(FONC) test to $J(\Delta \b U)$,
\[
\frac{\partial J}{\partial \Delta \b U}=\b 0^{\top}.
\]
Then, we solve the above equation for $\Delta \b U=\Delta \b U^*$, where
\begin{eqnarray*}
\frac{\partial J}{\partial \Delta \b U}&=& -\left( \b r_p -\b W\b x_a - \b Z\Delta \b U \right)^{\top}\b Q\b Z + \Delta \b U^{\top}\b R\\
&=&\b 0^{\top}.
\end{eqnarray*}
Performing simple manipulations yields
\[
-\b r_p^{\top}\b Q\b Z + \b x_a^{\top}\b W^{\top}\b Q\b Z +\Delta \b U^{\top}\b Z^{\top}\b Q \b Z  + \Delta \b U^{\top}\b R=\b 0^{\top}.
\]
Applying the transposition operation to both sides of the above equation and rearranging terms, we obtain
\[
\left(\b R+ \b Z^{\top}\b Q\b Z  \right)\Delta\b U =\b Z^{\top}\b Q \left(\b r_p - \b W\b x_a  \right).
\]
Note that the matrix $\left(\b R+ \b Z^{\top}\b Q\b Z  \right)$ is invertible, and in fact, positive definite because $\b R=\b R^{\top} > 0$ and $\b Z^{\top}\b Q\b Z$ is also symmetric and at least positive semi-definite. Hence, $\Delta \b U$ that satisfies the FONC is
\begin{equation}
\Delta \b U^* = \left(\b R+ \b Z^{\top}\b Q\b Z  \right)^{-1}\b Z^{\top}\b Q \left(\b r_p - \b W\b x_a  \right).
\label{eq:big_delta_U}
\end{equation}
Now, applying the second derivative test to $J(\Delta \b U)$, which we refer to as the second-order sufficiency condition~(SONC), we obtain
\begin{eqnarray*}
\frac{\partial^2 J}{\partial \Delta \b U^2} &=& \b R + \b Z^{\top}\b Q\b Z\\
&>& 0,
\end{eqnarray*}
which implies that $\Delta \b U^*$ is a strict minimizer of $J$.

Using~(\ref{eq:big_delta_U}), we compute $\Delta \b u(k)$,
\begin{eqnarray}
\Delta \b u(k) &=& \overbrace{\left[\begin{array}{cccc}
\b I_m & \b O & \cdots & \b O
\end{array}\right]}^{N_p \mb{\footnotesize{  block matrices}}}\left(\b R+ \b Z^{\top}\b Q\b Z  \right)^{-1}\b Z^{\top}\b Q \left(\b r_p - \b W\b x_a  \right)\nonumber\\
&=& \b K_{\b r}\b r_p -\b K_{\b x}\Delta\b x(k) -\b K_{\b y}\b y(k),
\label{eq:lin_sys_dMPC}
\end{eqnarray}
where
\[
\b K_{\b r}=\left[\begin{array}{cccc}
\b I_m & \b O & \cdots & \b O
\end{array}\right]\left(\b R+ \b Z^{\top}\b Q\b Z  \right)^{-1}\b Z^{\top}\b Q, \quad \b K_{\b x}=\b K_{\b r}\b W\left[\begin{array}{c}
\b I_n\\
\b O
\end{array}\right],
\]
and
\[
\b K_{\b y}=\b K_{\b r}\b W\left[\begin{array}{c}
\b O\\
\b I_p
\end{array}\right].
\]
An implementation of the above controller, using a discrete-time integrator, is shown in Figure~\ref{fig:lin_sys_dMPC}.
\begin{figure}
\centering \epsfig{file=lin_sys_dMPC_2009.eps,width=14cm}
\caption{Discrete-time MPC for linear time-invariant systems.} \label{fig:lin_sys_dMPC}
\end{figure}





\section{MPC With Constraints}
An attractive feature of the model-based predictive control approach is that a control engineer can incorporate different types of constraints on the control action. We consider three types of such constraints.
\subsection{Constraints on the Rate of Change of the Control Action}
Hard constraints on the rate of change of the control signal can be expressed as
\begin{equation}
\Delta u_i^{\min} \le \Delta u_i(k) \le \Delta u_i^{\max}, \quad i=1,2,\ldots, m
\label{eq:const_rate_of_ch}
\end{equation}
Let
\[
\Delta \b u^{\min} = \left[\begin{array}{ccc}
\Delta u_1^{\min}& \cdots & \Delta u_m^{\min}
\end{array}\right]^{\top}\quad\mb{and}\quad \Delta \b u^{\max} = \left[\begin{array}{cccc}
\Delta u_1^{\max}& \cdots & \Delta u_m^{\max}
\end{array}\right]^{\top}.
\]
Then, we can express~(\ref{eq:const_rate_of_ch}) as
\begin{equation}
\Delta \b u^{\min} \le \Delta \b u(k) \le \Delta \b u^{\max}.
\label{eq:vec_const_rate_of_ch}
\end{equation}
The above can be equivalently represented as
\begin{equation}
\left[\begin{array}{c}
-\b I_m\\
\b I_m
\end{array}\right]\Delta \b u(k) \le \left[\begin{array}{c}
-\Delta \b u^{\min}\\
\Delta \b u^{\max}
\end{array}\right].
\label{eq:con_in_delU}
\end{equation}
Using the above approach we can represent constraints on the rate of change of the control action over the whole prediction horizon, $N_p$, in terms of $\Delta \b U(k)$, by augmenting the above inequality to incorporate constraints for the remaining sampling times. That is, if the above constraints are  imposed on the rate of change of the control action for all sampling times within the prediction horizon, then this can be expressed in terms of $\Delta \b U(k)$ as
\begin{equation}
\left[\begin{array}{ccccc}
-\b I_m & \b O & \cdots & \b O & \b O\\
\b I_m & \b O & \cdots & \b O & \b O\\
\b O & -\b I_m & \cdots & \b O & \b O\\
\b O & \b I_m & \cdots & \b O & \b O\\
\vdots& {}& {} & {} & \vdots\\
\b O & \b O & \cdots & -\b I_m & \b O\\
\b O & \b O & \cdots & \b I_m & \b O\\
\b O & \b O & \cdots & \b O & -\b I_m\\
\b O & \b O & \cdots & \b O & \b I_m
\end{array}\right]\left[\begin{array}{c}
\Delta\b u(k)\\
\Delta\b u(k+1)\\
\vdots\\
\Delta\b u(k+N_p-2)\\
\Delta\b u(k+N_p-1)
\end{array}\right]\le
\left[\begin{array}{c}
-\Delta \b u^{\min}\\
\Delta \b u^{\max}\\
-\Delta \b u^{\min}\\
\Delta \b u^{\max}\\
\vdots\\
-\Delta \b u^{\min}\\
\Delta \b u^{\max}\\
-\Delta \b u^{\min}\\
\Delta \b u^{\max}
\end{array}\right].
\label{eqn:full_rate_co}
\end{equation}
On the other hand, if the constraints on the rate of change of the control are imposed only on the first component of $\Delta \b U(k)$, then we express this in terms of $\Delta \b U(k)$ as
\[
\left[\begin{array}{ccccc}
-\b I_m & \b O & \cdots & \b O & \b O\\
\b I_m & \b O & \cdots & \b O & \b O
\end{array}\right]\Delta \b U(k)\le \left[\begin{array}{c}
-\Delta \b u^{\min}\\
\Delta \b u^{\max}
\end{array}\right].
\]
%\begin{equation}
%\left[\begin{array}{c}
%-\b I_{N_p}\\
%\b I_m
%\end{array}\right]\Delta \b u \le \left[\begin{array}{c}
%-\Delta \b u^{\min}\\
%\Delta \b u^{\max}
%\end{array}\right].
%\label{eq:Pred_ho_con_in_delU}
%\end{equation}
\subsection{Constraints on the Control Action Magnitude}
Hard constraints on the control action magnitude at the time sampling $k$ have the form
\begin{equation}
u_i^{\min} \le u_i(k) \le u_i^{\max},\quad i=1,2\ldots, m
\label{eq:con_on mag}
\end{equation}
We now express constraints on the control action magnitude over the whole prediction horizon in term of $\Delta \b U(k)$. To accomplish our goal, we first note that
\begin{eqnarray*}
\b u(k)&=&\b u(k-1)+\Delta \b u(k)\\
&=& \b u(k-1) + \left[\begin{array}{cccc}
\b I_m & \b O & \cdots & \b O
\end{array}\right]\Delta \b U(k),
\end{eqnarray*}
where, recall that $\Delta \b U(k)=\left[\begin{array}{cccc}
\Delta \b u(k) & \Delta \b u(k+1) & \cdots & \Delta \b u(k+N_p-1)
\end{array}\right]^{\top}$.

Similarly,
\begin{eqnarray*}
\b u(k+1)&=&\b u(k)+\Delta \b u(k+1)\\
&=& \b u(k-1) + \left[\begin{array}{cccc}
\b I_m & \b O & \cdots & \b O
\end{array}\right]\Delta \b U(k) +\Delta \b u(k+1)\\
&=& \b u(k-1) + \left[\begin{array}{cccc}
\b I_m & \b I_m & \cdots & \b O
\end{array}\right]\Delta \b U(k).
\end{eqnarray*}
Continuing in this manner, we obtain
\begin{equation}
\left[\begin{array}{c}
\b u(k)\\
\b u(k+1)\\
\vdots\\
\b u(k+N_p-1)
\end{array}\right]=\left[\begin{array}{c}
\b I_m\\
\b I_m\\
\vdots\\
\b I_m
\end{array}\right]\b u(k-1)+ \left[\begin{array}{cccc}
\b I_m & \b O &  \cdots & \b O\\
\b I_m & \b I_m & \cdots & \b O\\
\vdots & {}&  \ddots& \vdots\\
\b I_m & \b I_m &  \cdots & \b I_m\\
\end{array}\right]\left[\begin{array}{c}
\Delta\b u(k)\\
\Delta\b u(k+1)\\
\vdots\\
\Delta\b u(k+N_p-1)
\end{array}\right].
\label{eq:detailed_con_mag}
\end{equation}
Let
\[
\b E=\left[\begin{array}{c}
\b I_m\\
\b I_m\\
\vdots\\
\b I_m
\end{array}\right]\quad\mb{and}\quad \b H=\left[\begin{array}{cccc}
\b I_m & \b O &  \cdots & \b O\\
\b I_m & \b I_m & \cdots & \b O\\
\vdots & {}&  \ddots& \vdots\\
\b I_m & \b I_m &  \cdots & \b I_m\\
\end{array}\right].
\]
Then, we can represent~(\ref{eq:detailed_con_mag}) as
\begin{equation}
\b U(k) = \b E\b u(k-1) + \b H\Delta \b U(k).
\label{eq:simp_rep_mag_con}
\end{equation}
Suppose now that we are faced with constructing a control action subject to the following constraints,
\[
\b U^{\min} \le \b U(k) \le \b U^{\max}.
\]
The above constraints can be equivalently represented as
\[
\left[\begin{array}{c}
-\b U(k)\\
\b U(k)
\end{array}\right] \le \left[\begin{array}{c}
-\b U^{\min}\\
\b U^{\max}
\end{array}\right],
\]
Taking into account~(\ref{eq:simp_rep_mag_con}), we write the above as
\begin{equation}
\left[\begin{array}{c}
-\left(\b E\b u(k-1) + \b H\Delta \b U(k)  \right)\\
\b E\b u(k-1) + \b H\Delta \b U(k)
\end{array}\right]\le\left[\begin{array}{c}
-\b U^{\min}\\
\b U^{\max}
\end{array}\right].
\label{eq:more_det_rep}
\end{equation}
The above, in turn, can be represented as
\begin{equation}
\left[\begin{array}{c}
-\b H \\
\b H
\end{array}\right]\Delta\b U(k)\le\left[\begin{array}{c}
-\b U^{\min}+\b E\b u(k-1)\\
\b U^{\max}-\b E\b u(k-1)
\end{array}\right].
\label{eq:contro_in_co}
\end{equation}
In a special case, when a control designer elects to impose constraints only on the first component of $\Delta \b U(k)$, that is, on $\Delta \b u(k)$ only, then this scenario is expressed in terms of $\Delta \b U(k)$ as
\begin{equation}
\left[\begin{array}{cccc}
-\b I_m &\b O & \cdots & \b O\\
\b I_m &\b O & \cdots & \b O
\end{array}\right]\Delta\b U(k)\le\left[\begin{array}{c}
-\b U^{\min}+\b E\b u(k-1)\\
\b U^{\max}-\b E\b u(k-1)
\end{array}\right].
\label{eq:mag_co_first}
\end{equation}

\subsection{Constraints on the Plant Output}
Recall from~(\ref{eq:compact_Y}) that the predicted plant output is, $\b Y(k)=\b W\b x_a(k)+\b Z\Delta \b U(k)$. Suppose now that the following constraints are imposed on the predicted plant's output,
\[
\b Y^{\min} \le \b Y(k) \le \b Y^{\max}.
\]
We represent the above as
\[
\left[\begin{array}{c}
-\b Y(k)\\
\b Y(k)
\end{array}\right] \le \left[\begin{array}{c}
-\b Y^{\min}\\
\b Y^{\max}
\end{array}\right].
\]
Taking into the account the expression for $\b Y$ given by~(\ref{eq:compact_Y}), we represent the above as
\begin{equation}
\left[\begin{array}{c}
-\b Z\\
\b Z
\end{array}\right]\Delta \b U(k) \le \left[\begin{array}{c}
-\b Y^{\min} + \b W \b x_a(k)\\
\b Y^{\max} - \b W\b x_a(k)
\end{array}\right].
\label{eq:plan_out_cons_det}
\end{equation}

The above discussion clearly demonstrates that we need an effective method of minimizing a function of many variables, $J(\Delta \b U)$, subject to inequality constraints such as given by~(\ref{eq:con_in_delU}), (\ref{eq:more_det_rep}), and~(\ref{eq:plan_out_cons_det}), or their combination. In the following, we present a powerful method for solving such optimization problems. This is the method that we shall use to implement our MPCs.


\section{An Optimizer for Solving Constrained Optimization Problems}
If follows from the discussion in the previous section that at each sampling time the MPC calls for a solution to a constrained optimization problem of the form,
\begin{eqnarray*}
\mb{minimize}&{}& J(\Delta \b U)\\
\mb{subject to}&{}& \b g(\Delta \b U)\le \b 0,
\end{eqnarray*}
where $\b g(\Delta \b U)\le \b 0$ contains inequality constrains given by~(\ref{eq:con_in_delU}), (\ref{eq:more_det_rep}), and~(\ref{eq:plan_out_cons_det}), or their combination. In this section, we present an iterative methods for solving the above optimization problems. To proceed, we first present the descent gradient method, followed by the Newton's method, for solving unconstrained optimization problems of the form,
\[
\mb{minimize } J(\Delta \b U)
\]

\subsection{Gradient Descent Method}
For the sake of simplicity, we denote the argument of a function of many variables as $\b x$, where $\b x\in \mb{\Bb R}^N$, and the function that we will be minimizing will be denoted as $f$, where $f: \mb{\Bb R}^N \to \mb{\Bb R}$.

The method of the gradient descent is based on the following property of the gradient of a differentiable function, $f$, on $\mb{\Bb R}^N$:
\begin{thm}
At a given point $\b x^{(0)}$, the vector
\[
\b v=-\nabla f\left( \b x^{(0)} \right)
\]
points in the direction of most rapid decrease of $f$ and the rate of increase of $f$ at $\b x^{(0)}$ in the direction $\b v$ is $-\left\Arrowvert\nabla f\left( \b x^{(0)} \right)\right\Arrowvert$, equivalently, the rate of decrease of $f$ at $\b x^{(0)}$ in the direction $\b v$ is $\left\Arrowvert\nabla f\left( \b x^{(0)} \right)\right\Arrowvert$.
\end{thm}

Thus, if we wish to minimize a differentiable function, $f$, then moving in the direction of the negative gradient is a good direction. The gradient descent algorithm rests on the above observation and has the form,
\[
\b x^{(k+1)}=\b x^{(k)}-\alpha \nabla f\left( \b x^{(k)}\right),
\]
where $\alpha >0$ is a step size.

\subsection{Second-Order Necessary Conditions for a Minimum}
Suppose we are given a function $f$ of one variable $x$. Recall a well-known theorem referred to as the second-order Taylor's formula or the extended law of mean~\cite[p.~1]{peressini88}:
\begin{thm}
Suppose that $f(x)$, $f'(x)$, $f''(x)$ exist on the closed interval $[a, b]=\{x\in \mb{\Bb R}: a\le x \le b\}$. If $x^*$, $x$ are any two different points of $[a, b]$, then there exists a point $z$ strictly between $x^*$ and $x$ such that
\[
f(x)=f\left(x^* \right) + f'\left(x^* \right)\left(x-x^* \right) + \frac{f''(z)}{2}\left(x-x^* \right)^2.
\]
\end{thm}
Using the above formula we observe that if
\begin{itemize}
\item $f'\left(x^* \right)=0$, and
\item $f''\left(x^* \right)>0$,
\end{itemize}
then
\[
f(x)=f\left(x^* \right)+\mb{ a positive number}
\]
for all $x$ ``close'' to $x^*$. Indeed, if $f''(x)$ is continuous at $x^*$ and $f''\left(x^* \right)>0$, then $f''\left(x \right)>0$ for all $x$ in some neighborhood of $x^*$. Therefore,
\[
f(x) > f\left(x^* \right)\quad \mb{ for all } x \mb{ close to } x^*,
\]
which means that $x^*$ is a strict local minimizer of $f$.

We now extend the above result for functions of many variables. To proceed, we need the second-order Taylor's formula for such functions that can be found in~\cite[p.~11]{peressini88}.
\begin{thm}
Suppose that $\b x^*$, $\b x$ are points in $\mb{\Bb R}^N$ and that $f$ is a function of $N$ variables with continuous first and second partial derivatives on some open set containing the line segment
\[
\left[\b x^*, \b x  \right]=\left\{\b w \in \mb{\Bb R}^N: \b w=\b x^* + t\left(\b x-\b x^*  \right); 0\le t\le 1  \right\}
\]
joining $\b x^*$ and $\b x$. Then there exists a $\b z\in \left[\b x^*, \b x  \right]$ such that
\begin{equation}
f(\b x)=f\left(\b x^* \right) + \nabla f\left(\b x^* \right)^{\top} \left(\b x-\b x^* \right) + \frac{1}{2} \left(\b x-\b x^* \right)^{\top}\b F\left( \b z \right)\left(\b x-\b x^* \right),
\label{eq:many_var_Taylor}
\end{equation}
where $\b F\left( \cdot \right)$ is the Hessian of $f$, that is, the second derivative of $f$.
\end{thm}
Hence if
\begin{itemize}
\item $\b x^*$ is a critical point, that is, $\nabla f\left(\b x^* \right)=\b 0$, and
\item $\b F\left( \b x^* \right) >0$,
\end{itemize}
then using~(\ref{eq:many_var_Taylor}), we conclude that
\[
f(\b x)=f\left(\b x^* \right) + 0 + \mb{a positive number}
\]
for all $\b x$ in a neighborhood of $\b x^*$. Therefore for all $\b x\ne \b x^*$ in some neighborhood of $\b x^*$, we have
\[
f(\b x)>f\left(\b x^* \right),
\]
which implies that $\b x^*$ is a strict local minimizer.

\subsection{Newton's Method}
The idea behind Newton's method for function minimization is to minimize the quadratic approximation rather than the function itself as illustrated in Figure~\ref{fig:newton_many}.
 \begin{figure}[h]
\centering \epsfig{file=newton.eps,width=10cm}
\caption{Newton's method minimizes the quadratic approximation of the objective function that utilizes first and second derivatives of the optimized function.}
\label{fig:newton_many}
\end{figure}
Newton's method seeks a critical point, $\b x^*$, of a given function. If at this critical point we have $\b F\left( \b x^* \right) >0$, then $\b x^*$ is a strict local minimizer of $f$.

We can obtain a quadratic approximation $q$ of $f$ at $\b x^*$ from the second-order Taylor series expansion of $f$ about $\b x^*$,
\[
q(\b x)= f\left(\b x^* \right) + \nabla f\left(\b x^* \right)^{\top} \left(\b x-\b x^* \right) + \frac{1}{2} \left(\b x-\b x^* \right)^{\top}\b F\left( \b x^* \right)\left(\b x-\b x^* \right).
\]
Note that
\[
q\left( \b x^*\right)=f\left( \b x^* \right),\quad \nabla q\left( \b x^*\right)=\nabla f\left( \b x^* \right),
\]
as well as their Hessians, that is, their second derivatives evaluated at $\b x^*$ are equal.

A critical point of $q$ can be obtained by solving the algebraic equation,
\[
\nabla q (\b x)=\b 0,
\]
that is, by solving the equation
\[
\nabla q(\b x)= \nabla f\left(\b x^* \right) + \b F\left( \b x^* \right)\left(\b x-\b x^* \right)=\b 0.
\]

Suppose now that we have a quadratic approximation of $f$ at a point $\b x^{(k)}$, that is,
\[
q(\b x)= f\left(\b x^{(k)}\right) + \nabla f\left(\b x^{(k)} \right)^{\top} \left(\b x-\b x^{(k)} \right) + \frac{1}{2} \left(\b x-\b x^{(k)} \right)^{\top}\b F\left( \b x^{(k)} \right)\left(\b x-\b x^{(k)} \right).
\]
We assume that $\det \b F\left( \b x^{(k)} \right) \ne 0$. Denoting the solution to the above equation as $\b x^{(k+1)}$, we obtain
\[
\b x^{(k+1)}=\b x^{(k)} - \b F\left( \b x^{(k)}\right)^{-1}\nabla f\left(\b x^{(k)} \right)
\]
The above is known as the Newton's method for minimizing a function of many variables $f$. Note that $\b x^{(k+1)}$ is a critical point of the quadratic function $q$ that approximates $f$ at $\b x^{(k)}$.

A computationally efficient representation of Newton's algorithm has the form,
\[
\b x^{(k+1)} =\b x^{(k)}-\Delta \b x^{(k)},
\]
where $\Delta \b x^{(k)}$ is obtained by solving the equation,
\[
\b F\left( \b x^{(k)}\right)\Delta \b x^{(k)} = \nabla f\left(\b x^{(k)} \right).
\]

\section{Minimization Subject to Equality Constraints}
We now discuss the problem of finding a point $\b x\in \mb{\Bb R}^N$ that minimizes $f(\b x)$ subject to equality constraints,
\[
\left.\begin{array}{ccc}
h_1(\b x) &=& 0\\
h_2(\b x) &=& 0\\
{}&\vdots &{}\\
h_M(\b x) &=& 0
\end{array}\right\}
\]
where $M\le N$. We write the above equality constraints in a compact form as
\[
\b h(\b x)=\b 0,
\]
where $\b h:\mb{\Bb R}^N \to \mb{\Bb R}^M$. We refer to the set of points satisfying the above constraints as the surface.

We now introduce the notion of the tangent plane to the surface $S=\{\b x:\b h(\b x)=\b 0\}$ at a point $\b x^*\in S$ . First, we define a curve on the surface $S$ as a family of points $\b x(t)\in S$ continuously parameterized by $t$ for $t\in [a, b]$. The curve is differentiable if $\dot{\b x}(t)=d\b x(t)/dt$ exists. A curve is said to pass through the point $\b x^*\in S$ if $\b x^*=\b x\left(t^* \right)$ for some $t^* \in [a, b]$. The tangent plane to the surface $S=\{\b x:\b h(\b x)=\b 0\}$ at $\b x^*$ is the collection of the derivatives at $\b x^*$ of all differentiable curves on $S$ that pass through $\b x^*$.

To proceed, we need the notion of a regular point of the constraints. We say that $\b x^*$ satisfying the constraints, that is, $\b h\left(\b x^*  \right)=\b 0$, is a regular point of the constraints if the gradient vectors,
\[
\nabla h_1\left(\b x^*  \right), \ldots , \nabla h_M\left(\b x^*  \right)
\]
are linearly independent.

One can show, see, for example~\cite[p.~298]{luenberger84}, that the tangent space at a regular point $\b x^*$, denoted $T\left(\b x^*  \right)$, to the surface $\{\b x: \b h\left(\b x  \right)=\b 0\}$ at the regular point $\b x^*$ is
\begin{equation}
T\left(\b x^*  \right)=\left\{\b y: \left[\begin{array}{c}
\nabla h_1\left(\b x^*  \right)^{\top}\\
  \vdots\\
  \nabla h_M\left(\b x^*  \right)^{\top}
  \end{array}
  \right]\b y=\b 0  \right\}.
\label{eq:tangent_plane}
\end{equation}
With the above notions in place, we are ready to state and prove the following result which is known as the first-order necessary condition~(FONC) for function minimization subject to equality constraints.
\begin{thm}
Let $\b x^*$ be a local minimizer (or maximizer) of $f$ subject to the constraints $\b h(\b x)=\b 0$ and let $\b x^*$ be a regular point of the constraints. Then there exists a vector $\b{\lambda}^*$ such that
\[
\nabla f\left(\b x^*  \right) + \left[\begin{array}{ccc}
\nabla h_1\left(\b x^*  \right)&\cdots &\nabla h_M\left(\b x^*  \right)
\end{array}\right]\b{\lambda}^* = \b 0.
\]
\label{thm:lagrange_multipliers}
\end{thm}
{\bf Proof}\quad Let $\b x(t)$ be a differentiable curve passing through $\b x^*$ on the surface $S=\{\b x:\b h(\b x)=\b 0\}$ such that $\dot{\b x}\left(t^*\right)=\b y$ where $t^*\in [a, b]$. Note that the tangent plane is given by~(\ref{eq:tangent_plane}). Because $\b x^*$ is a local minimizer of $f$ on $S$, we have
\[
\left.\frac{d}{dt}f(\b x(t))\right\arrowvert_{t=t^*}=0.
\]
Applying the chain rule to the above gives
\[
\nabla f\left(\b x^*  \right)^{\top}\b y = 0.
\]
Thus $\nabla f\left(\b x^*  \right)$ is orthogonal to the tangent space $T\left(\b x^*  \right)$. That is, $\nabla f\left(\b x^*  \right)$ is a linear combination of the gradients $\nabla h_1\left(\b x^*  \right), \ldots , \nabla h_M\left(\b x^*  \right)$. This fact can be expressed as
\[
\nabla f\left(\b x^*  \right) + \left[\begin{array}{ccc}
\nabla h_1\left(\b x^*  \right)&\cdots &\nabla h_M\left(\b x^*  \right)
\end{array}\right]\b{\lambda}^* = \b 0
\]
for some constant vector $\b{\lambda}^* \in \mb{\Bb R}^M$.\\
\hspace*{5.5in}$\Box$

The vector $\b{\lambda}^*$ is called the vector of Lagrange multipliers.

We now introduce the Lagrangian associated with the constrained optimization problem,
\[
l(\b x,\b{\lambda})=f(\b x)+\b{\lambda}^{\top}\b h(\b x).
\]
Then the FONC can be expressed as
\begin{eqnarray*}
\nabla_{\b x} l(\b x,\b{\lambda}) &=& \b 0\\
\nabla_{\b{\lambda}} l(\b x,\b{\lambda}) &=& \b 0.
\end{eqnarray*}
Note that the second of the above condition is equivalent to $\b h(\b x)=\b 0$, that is,
\[
\nabla_{\b{\lambda}} l(\b x,\b{\lambda}) = \b h(\b x)=\b 0.
\]
Equivalently the FONC can be written as
\begin{eqnarray*}
\nabla_{\b x} l(\b x,\b{\lambda}) &=& \b 0\\
\b h(\b x) &=& \b 0.
\end{eqnarray*}
We now apply Newton's method to solve the above system of equations iteratively,
\[
\left[\begin{array}{c}
\b x^{(k+1)}\\
\b{\lambda}^{(k+1)}
\end{array}\right]=\left[\begin{array}{c}
\b x^{(k)}\\
\b{\lambda}^{(k)}
\end{array}\right]+\left[\begin{array}{c}
\b d^{(k)}\\
\b y^{(k)}
\end{array}\right],
\]
where $\b d^{(k)}$ and $\b y^{(k)}$ are obtained by solving the matrix equation,
\[
\left[\begin{array}{cc}
\b L\left(\b x^{(k)},\b{\lambda}^{(k)}\right) & D\b h\left(\b x^{(k)}\right)^{\top}\\
D\b h\left(\b x^{(k)}\right) & \b O
\end{array}\right]\left[\begin{array}{c}
\b d^{(k)}\\
\b y^{(k)}
\end{array}\right]=\left[\begin{array}{c}
-\nabla_{\b x} l(\b x,\b{\lambda})\\
-\b h\left(\b x^{(k)}\right)
\end{array}\right],
\]
where $\b L\left(\b x,\b{\lambda}\right)$ is the Hessian of $l(\b x,\b{\lambda})$ with respect to $\b x$, and $D\b h(\b x)$ is the Jacobian matrix of $\b h(\b x)$, that is,
\[
D\b h(\b x)=\left[\begin{array}{c}
\nabla h_1\left(\b x  \right)^{\top}\\
  \vdots\\
  \nabla h_M\left(\b x  \right)^{\top}
  \end{array}
  \right].
\]

The above algorithm is also referred to as in the literature as sequential quadratic programming~(SQP); see, for example~\cite[Section~15.5]{nash96}.
\section{A Lagrangian Algorithm for Equality Constraints}

The first-order Lagrangian algorithm for the optimization problem involving minimizing $f$ subject to the equality constraints, $\b h(\b x)=\b 0$, has the form,
\begin{eqnarray*}
\b x^{(k+1)} &=& \b x^{(k)}-\alpha_k \left(\nabla f\left(  \b x^{(k)}\right)+D\b h\left(\b x^{(k)}  \right)^{\top}\b{\lambda}^{(k)}  \right)\\
\b{\lambda}^{(k+1)}  &=& \b{\lambda}^{(k)} +\beta_k \b h\left(\b x^{(k)}  \right),
\end{eqnarray*}
where $\alpha_k$ and $\beta_k$ are positive constants. Note that the update for $\b x^{(k)}$ is a descent gradient for minimizing the Lagrangian with respect to $\b x$, while the update for $\b{\lambda}^{(k)}$ is a gradient ascent for maximizing the Lagrangian with respect to $\b{\lambda}$. For an analysis of the above algorithm, we recommend~\cite[pp.~522--524]{chong08}

\section{Minimization Subject to Inequality Constraints}
We now discuss the problem of finding a point $\b x\in \mb{\Bb R}^N$ that minimizes $f(\b x)$ subject to inequality constraints,
\[
\left.\begin{array}{ccc}
g_1(\b x) &\le& 0\\
g_2(\b x) &\le& 0\\
{}&\vdots &{}\\
g_P(\b x) &\le& 0
\end{array}\right\}
\]
We write the above inequality constraints in a compact form as
\[
\b g(\b x)\le\b 0,
\]
where $\b h:\mb{\Bb R}^N \to \mb{\Bb R}^P$.

Let $\b x^*$ be a point satisfying the constraints, that is, $\b g\left(\b x^*  \right)\le\b 0$, and let $J$ be the set of indices $j$ for which $g_j\left(\b x^* \right)=0$. Then $\b x^*$ is said to be a regular point of the constraints $\b g(\b x)\le\b 0$ if the gradient vectors,
\[
\nabla g_j\left(\b x^*  \right), \quad j\in J,
\]
are linearly independent.

We say that a constraint $g_j(\b x)\le 0$ is active at $\b x^*$ if $g_j\left(\b x^* \right)=0$. Thus the index set $J$, defined above, contains indices of active constraints.

\begin{thm}
Let $\b x^*$ be a regular point and a local minimizer of $f$ subject to $\b g(\b x)\le\b 0$. Then there exists a vector $\b{\mu}^*\in \mb{\Bb R}^P$ such that
\begin{enumerate}
\item $\b{\mu}^* \ge 0$,
\item $\nabla f\left(\b x^*  \right) + \left[\begin{array}{ccc}
\nabla g_1\left(\b x^*  \right)&\cdots &\nabla g_P\left(\b x^*  \right)
\end{array}\right]\b{\mu}^* = \b 0$,
\item $\b{\mu}^{*\top} \b g\left(\b x^*  \right)=0$.
\end{enumerate}
\label{thm:kkt_multipliers}
\end{thm}
{\bf Proof}\quad First note that because $\b x^*$ is a relative minimizer over the constraint set $\{\b x: \b g(\b x)\le \b 0\}$, it is also a minimizer over a subset of the constraint set obtained by setting the active constraints to zero. Therefore, for the resulting equality constrained problem, by Theorem~\ref{thm:lagrange_multipliers}, we have
\begin{equation}
\nabla f\left(\b x^*  \right) + \left[\begin{array}{ccc}
\nabla g_1\left(\b x^*  \right)&\cdots &\nabla g_P\left(\b x^*  \right)
\end{array}\right]\b{\mu}^* = \b 0,
\label{eq:app_lagrang_ine}
\end{equation}
where $\mu_j^*=0$ if $g_j\left(\b x^*  \right)<0$.
This means that $\b{\mu}^{*\top} \b g\left( \b x^*\right)=0$. Thus $\mu_i$ may be non-zero only if the corresponding constraint is active, that is, $g_i\left(  \b x^* \right)=0$.

It thus remains to show that $\mu^*\ge \b 0$. We prove this by contraposition. We suppose that for some $k\in J$, we have $\mu_k^*<0$. Consider next a surface formed by all other active constraints, that is, the surface
\begin{equation}
\{\b x: g_j(\b x)=0, j\in J, j\ne k\}.
\label{eq:surface_minus_one}
\end{equation}
Also consider a tangent space to the above surface at $\b x^*$. By assumption $\b x^*$ is a regular point of the active constraints. Therefore, there exists a vector $\b y$ such that
\[
\nabla g_k\left(\b y \right)^{\top}\b y<0.
\]
Let now $\b x(t)$, $t\in [-a, a]$, $a>0$, be a curve on the surface~(\ref{eq:surface_minus_one}) such that
\[
\dot{\b x}(0)=\b y.
\]
Note that for small $t$, the curve $\b x(t)$ is feasible. We apply the transposition operator to the both sides of~(\ref{eq:app_lagrang_ine}) to obtain
\[
\nabla f\left(\b x^*  \right)^{\top} + \sum_{i=1}^P \mu_i^*
\nabla g_i\left(\b x^*  \right)^{\top} = \b 0^{\top}.
\]
Post-multiplying the above by $\b y$ an taking into account the fact that $\b y$ belongs to the tangent space to the surface~(\ref{eq:surface_minus_one}) gives
\begin{eqnarray*}
\nabla f\left(\b x^*  \right)^{\top}\b y &=& -\mu_k^*\nabla g_k\left(\b x^*  \right)^{\top}\b y\\
&<& 0.
\end{eqnarray*}
Suppose, without loss of generality, that $\|\b y\|_2=1$. Then, we have that
\begin{eqnarray*}
\frac{d f(\b x(t))}{dt}&=&\nabla f\left( \b x^* \right)^{\top}\b y\\
&<& 0,
\label{eq:rate_increase}
\end{eqnarray*}
that is, the rate of increase of $f$ at $\b x^*$ in the direction $\b y$ is negative. This would mean that we could decrease the value of $f$ moving just slightly away from $\b x^*$ along $\b y$ while, at the same time, preserving feasibility. But this contradicts the minimality of $f\left(\b x^* \right)$. In summary, if $\b x^*$ is a relative minimizer then we also have the components of $\b{\mu}^*$ all non-negative.\\
\hspace*{5.5in}$\Box$


The vector $\b{\mu}^*$ is called the vector of the Karush-Kuhn-Tucker~(KKT) multipliers.


\section{A Lagrangian Algorithm for Inequality Constraints}
\label{sec:lag_inequality}
We now present a first-order Lagrangian algorithm for the optimization problem involving inequality constraints,
\begin{eqnarray*}
\mb{minimize}&{}& f(\b x)\\
\mb{subject to}&{}& \b g(\b x)\le \b 0,
\end{eqnarray*}
where $\b g:\mb{\Bb R}^N \to \b g:\mb{\Bb R}^P$. The Lagrangian function is
\[
l(\b x,\b{\mu})=f(\b x)+\b{\mu}^{\top}\b g(\b x).
\]
The first-order Lagrangian algorithm for the above optimization problem involving minimizing $f$ subject to the inequality constraints, $\b g(\b x)\le\b 0$, has the form,
\begin{eqnarray*}
\b x^{(k+1)} &=& \b x^{(k)}-\alpha_k \left(\nabla f\left(  \b x^{(k)}\right)+D\b g\left(\b x^{(k)}  \right)^{\top}\b{\mu}^{(k)}  \right)\\
\b{\mu}^{(k+1)}  &=& \left[\b{\mu}^{(k)} + \beta_k \b g\left(\b x^{(k)}  \right)\right]_+,
\end{eqnarray*}
where the operation $[\cdot]_+=\max(\cdot,0)$ is applied component-wise.

For an analysis of the above algorithm, we recommend~\cite[pp.~524--528]{chong08}

\begin{example}\label{ex1:DMPC}\rm
We illustrate the design of the discrete MPC subject to the constraints on the control input on a plant modeled by the following continuous state-space model,
\begin{equation}
\left.\begin{array}{ccl}
\dot{\b x}&=&\left[\begin{array}{cc}
-0.1 & -3.0\\
1 & 0
\end{array}\right]\b x + \left[\begin{array}{c}
1\\
0
\end{array}\right]u\\
y&=& \left[\begin{array}{cc}
0 & 10
\end{array}\right]\b x.
\end{array}\right\}
\label{eq:plant_ss_mo}
\end{equation}
The system with the transfer function as the above system was used by Wang~\cite[p.~70]{wang09} to test her discrete MPC designs. We use different optimizer than Wang.  The reader can compare the results of two different MPC synthesis methods. The MPC controller's objective is to force the plant's output to track the unit step.

We proceed to generate the discrete state-space model employing the MATLAB's command \verb!c2dm!. We use the sampling time $h=0.1$. The command input is the unit step. The weight matrices in the cost function~(\ref{eq:cost_function}) are $\b R=0.01 \b I_3$ and $\b Q=\b I_3$. The prediction horizon is $N_p=3$. Our objective is to design a discrete-MPC subject to the constraints on the control of the form,
\[
-0.3\le u(k)\le 0.5
\]
First, we consider the case when the constraints are imposed only on the first component of $\b U(k)$. The constraints take the form of~(\ref{eq:mag_co_first}). Thus, in our example, the constraints can be written as,
\[
\left[\begin{array}{ccc}
-1 & 0 & 0\\
1 & 0 & 0
\end{array}\right]\left[\begin{array}{c}
\Delta u(k)\\
\Delta u(k+1)\\
\Delta u(k+2)
\end{array}\right]\le \left[\begin{array}{c}
0.3 + u(k-1)\\
0.5 - u(k-1)
\end{array}\right].
\]
We employed the first-order Lagrangian algorithm for inequality constraints presented in Section~\ref{sec:lag_inequality}, where $\alpha_k=\beta_k=0.005$. We set zero initial conditions. Plots of the control effort and the plant's output are shown in Figure~\ref{fig:DMPC_ex1}. We wish to emphasize that the simulations were performed for the discretized plant. Thus, the control and output were evaluated at the sampling instances and marked on the plots using circles. For the purpose of visualization the control action was presented in the form of a staircase function while the plant output values were connected by straight lines.
\begin{figure}[h]
\centering \epsfig{file=DMPC_ex1.eps,width=11cm}
\caption{Plots of control effort and the plant's output versus time of the closed-loop system of Example~\ref{ex1:DMPC} for the case when the constraint was imposed only on the first component of $\b U(k)$.} \label{fig:DMPC_ex1}
\end{figure}

Next we tested the discrete MPC where the constraints, $-0.3\le u(k)\le 0.5$, were imposed on all the elements of $\b U(k)$. We use~(\ref{eq:contro_in_co}) to express the constraints, where
\[
\left[\begin{array}{ccc}
-1 & 0 & 0\\
-1 & -1 & 0\\
-1 & -1 & -1\\
1 & 0 & 0\\
1 & 1 & 0\\
1 & 1 & 1
\end{array}\right]\left[\begin{array}{c}
\Delta u(k)\\
\Delta u(k+1)\\
\Delta u(k+2)
\end{array}\right]\le \left[\begin{array}{c}
0.3 + u(k-1)\\
0.3 + u(k-1)\\
0.3 + u(k-1)\\
0.5 - u(k-1)\\
0.5 - u(k-1)\\
0.5 - u(k-1)
\end{array}\right].
\]
\begin{figure}[h]
\centering \epsfig{file=DMPC_ex2.eps,width=11cm}
\caption{Plots of control effort and the plant's output versus time of the closed-loop system of Example~\ref{ex1:DMPC} for the case when the constraints were imposed on all the components of $\b U(k)$.} \label{fig:DMPC_ex2}
\end{figure}
\end{example}
\begin{example}\label{ex1:CMPC}\rm
In this example, we apply the discrete MPC controller from Example~\ref{ex1:DMPC} to the continuous plant modeled by~(\ref{eq:plant_ss_mo}). The MPC controller's objective is to force the plant's output to track the unit step. The zero-order hold is applied to the discrete MPC controller's output sequence resulting in a piece-wise constant input to the continuous plant model. In Figure~\ref{fig:CMPC_ex1}, we show plots of the plant output as well as the control effort versus time. As in the previous example, the prediction horizon is $N_p=3$ and the constraints on the control are
\[
-0.3\le u(k)\le 0.5
\]
\begin{figure}[h]
\centering \epsfig{file=CMPC_ex1.eps,width=11cm}
\caption{Plots of control effort and the plant's output versus time of the closed-loop system of Example~\ref{ex1:CMPC} for the case when the constraints were imposed only on the first component of $\b U(k)$.} \label{fig:CMPC_ex1}
\end{figure}
\end{example}
\section{Observer-Based MPC}

If the plant state is not accessible, we need to use an observer in the controller implementation as shown in Figure~\ref{fig:MPC_observer}.
\begin{figure}
\centering \epsfig{file=MPC_observer_2009.eps,width=12cm}
\caption{MPC implementation when the plant's state is not available.} \label{fig:MPC_observer}
\end{figure}
\section{Notes}
The key reason for huge popularity of the MPC approach is its ability to systematically take into account constraints thus allowing processes to operate at the limits of achievable performance~\cite{kouvaritakis01}. For some impressive industrial applications of MPCs see~\cite{kouvaritakis01,camacho04}. Nonlinear model predictive controllers are presented in~\cite{mayne90, michalska95}.

\bibliography{Root2}
\bibliographystyle{plain}
\end{document}
