2) **Gradient step: While** *j* ≤ *N* **Do**

– Integrate backward in time

$$
\tilde{\boldsymbol{\lambda}}\_{k}^{(j)}(\tau) = -H\_{\mathbf{x}}(\tilde{\mathbf{x}}\_{k}^{(j)}(\tau), \tilde{\boldsymbol{\lambda}}\_{k}^{(j)}(\tau), \tilde{\boldsymbol{u}}\_{k}^{(j)}(\tau)), \quad \tilde{\boldsymbol{\lambda}}\_{k}^{(j)}(T) = V\_{\mathbf{x}}(\tilde{\mathbf{x}}\_{k}^{(j)}(T))\tag{26}
$$

– Compute the search direction

$$\bar{s}\_{\dot{k}}^{(j)}(\tau) = -H\_{\rm u}(\bar{\mathbf{x}}\_{\dot{k}}^{(j)}(\tau), \bar{\boldsymbol{\lambda}}\_{\dot{k}}^{(j)}(\tau), \bar{u}\_{\dot{k}}^{(j)}(\tau)), \quad \tau \in [0, T] \tag{27}$$

– Compute the step size *<sup>α</sup>*(*j*) *<sup>k</sup>* by (approximately) solving the line search problem

$$\mathfrak{a}\_k^{(j)} = \arg\min\_{a>0} J\left(\mathfrak{x}\_{k'}\psi(\mathfrak{a}\_k^{(j)} + a\mathfrak{s}\_k^{(j)})\right) \tag{28}$$

– Compute the new control trajectory

$$
\pi\_k^{(j+1)}(\tau) = \psi\left(\vec{u}\_k^{(j)}(\tau) + \alpha\_k^{(j)}\vec{s}\_k^{(j)}(\tau)\right) \tag{29}
$$

– Integrate forward in time

$$\dot{\mathfrak{x}}\_{k}^{(j+1)}(\tau) = f(\mathfrak{x}\_{k}^{(j+1)}(\tau), \mathfrak{u}\_{k}^{(j+1)}(\tau)) , \quad \mathfrak{x}\_{k}^{(j+1)}(0) = \mathfrak{x}\_{k} \tag{30}$$

$$\text{--}\quad\text{Quit if }|f(\mathbf{x}\_{k'}\boldsymbol{\pi}\_{k}^{(j+1)}) - f(\mathbf{x}\_{k'}\boldsymbol{\pi}\_{k}^{(j)})| \le \varepsilon\_{j}.\text{ Otherwise set }j \gets j+1 \text{ and return to }2\text{).}$$

Table 1. Gradient projection method for solving OCP (3).

Nonlinear Model Predictive Control 11

A Real-Time Gradient Method for Nonlinear Model Predictive Control 19

*<sup>α</sup>*<sup>ˆ</sup> <sup>=</sup> <sup>−</sup> *<sup>c</sup>*<sup>1</sup> 2*c*<sup>2</sup>

case, the interval [*α*1, *α*3] can be adapted by a scaling factor to track the minimum point of the line search problem (28) over the single gradient iterations. Table 2 summarizes the overall

In general, the gradient method in Table 1 is stopped if the convergence criterion is fulfilled for some tolerance *ε<sup>J</sup>* > 0. In practice this can lead to a large number of iterations that moreover varies from one MPC iteration to the next. In order to ensure a real-time feasible MPC implementation, the gradient algorithm is stopped after *N* iterations and the re-initialization

. (36)

*<sup>k</sup>* approximately solves the line

(*j*)

<sup>3</sup> , *<sup>κ</sup>*<sup>+</sup> <sup>=</sup> <sup>3</sup>

*<sup>α</sup>* = 0.1, *<sup>ε</sup>*<sup>+</sup>

2 )

*α*<sup>1</sup> if *J*<sup>1</sup> + *ε<sup>g</sup>* ≤ min{*J*2, *J*3} *α*<sup>3</sup> if *J*<sup>3</sup> + *ε<sup>g</sup>* ≤ min{*J*1, *J*2}

*α*<sup>2</sup> ←

*α*<sup>1</sup> + *α*<sup>3</sup>

<sup>2</sup> (39)

*<sup>α</sup>* = 0.9)

*<sup>k</sup>* )) at the sample points (32)

(37)

(38)

*<sup>k</sup>* is set to one of the interval bounds *α*<sup>1</sup> or *α*3. In this

If *c*<sup>2</sup> > 0, then the polynomial *g*(*α*) has a minimum at the point

If in addition *α*ˆ lies inside the interval [*α*1, *α*3], then *α*ˆ = *α*

of the algorithm is done as outlined in Section 3.1.

1) **Initialization: Default values and tolerances**

– Set interval adaptation tolerances *ε*−

– Compute the cost values *Ji* :<sup>=</sup> *<sup>J</sup>*(*xk*, *<sup>ψ</sup>*(*u*(*j*)

– Compute the approximate step size *α*

2) **Approximate line search**

[*α*1, *α*3] ←

⎧ ⎪⎨

⎪⎩

– Set interval adaptation factors *κ*−, *κ*<sup>+</sup> (e.g. *κ*<sup>−</sup> = <sup>2</sup>

– Set polynomial tolerances *εc*, *ε<sup>g</sup>* (e.g. *ε<sup>c</sup>* = 10−5, *ε<sup>g</sup>* = 10−6) – Set initial line search interval (32) (e.g. *α*<sup>1</sup> = 10−2, *α*<sup>3</sup> = 10−1)

if *c*<sup>2</sup> > *ε<sup>c</sup>* : *α*

else (*c*<sup>2</sup> ≤ *εc*) : *α*

*<sup>κ</sup>*+[*α*1, *<sup>α</sup>*3] if *<sup>α</sup>*<sup>ˆ</sup> <sup>≥</sup> *<sup>α</sup>*<sup>1</sup> <sup>+</sup> *<sup>ε</sup>*<sup>+</sup>

*κ*−[*α*1, *α*3] if *α*ˆ ≤ *α*<sup>1</sup> + *ε*<sup>−</sup>

Table 2. Adaptive line search for the gradient algorithm in Table 1.

[*α*1, *α*3] else

*<sup>α</sup>* , *<sup>ε</sup>*<sup>+</sup>

(*j*)

(*j*) *<sup>k</sup>* =

(*j*) *<sup>k</sup>* =

– Adapt the line search interval [*α*1, *α*3] for the next gradient iteration according to

– Set interval adaptation limits *α*min, *α*max (e.g. *α*min = 10−5, *α*max = 1.0)

– Compute the polynomial coefficients (35) and the candidate point (36)

*<sup>α</sup>* (e.g. *ε*<sup>−</sup>

*<sup>k</sup>* + *αis*

*<sup>k</sup>* according to

⎧ ⎪⎨

⎪⎩

⎧ ⎪⎨

⎪⎩

*α*<sup>2</sup> else

*<sup>α</sup>* (*α*<sup>3</sup> − *α*1) and *α*<sup>3</sup> ≤ *α*max

*<sup>α</sup>* (*α*<sup>3</sup> − *α*1) and *α*<sup>1</sup> ≥ *α*min ,

(*j*)

*α*<sup>1</sup> if *α*ˆ < *α*<sup>1</sup> *α*<sup>3</sup> if *α*ˆ > *α*<sup>3</sup> *α*ˆ else

(*j*)

algorithm for the approximate line search and the interval adaptation.

search problem (28). Otherwise, *α*

*ψ* = [*ψ*1,..., *ψm*] <sup>T</sup> in (28) represents a projection function of the form

$$\psi\_i(u\_i) = \begin{cases} u\_i^- & \text{if } u\_i < u\_i^- \\ u\_i^+ & \text{if } u\_i > u\_i^+ \\ u\_i & \text{else} \end{cases}, \quad i = 1, \dots, m \tag{31}$$

which guarantess the adherence of the input constraints [*u*−, *u*+]. For the real-time implementation within a suboptimal MPC scheme, the line search problem (28) can be solved in an approximate manner (see Section 4.2). Finally, the control trajectory *u*¯ (*j*+1) *<sup>k</sup>* (*τ*), *τ* ∈ [0, *T*] follows from evaluating (29) with *s*¯ (*j*) *<sup>k</sup>* (*τ*) and the step size *α* (*j*) *k* .

The convergence properties of the gradient (projection) method are investigated, for instance, in Dunn (1996); Leese (1977); Nikol'skii (2007). In particular, Dunn (1996) proved under certain convexity and regularity assumptions that the gradient method exhibits a linear rate of convergence of the form (15).

#### **4.2 Adaptive line search**

The line search (28) represents a scalar optimization problem that is often solved approximately. The most straightforward way is to use a fixed step size *α* throughout all gradient iterations. This, however, usually leads to a slow rate of convergence.

An attractive alternative to a constant step size is to use a polynomial approximation with an underlying interval adaptation. To this end, the cost functional *J* � *xk*, *ψ*(*u*¯ (*j*) *<sup>k</sup>* + *αs*¯ (*j*) *<sup>k</sup>* ) � in the line search problem (28) is evaluated at three sample points

$$
\mathfrak{a}\_1 < \mathfrak{a}\_2 < \mathfrak{a}\_3 \quad \text{with} \quad \mathfrak{a}\_2 = (\mathfrak{a}\_1 + \mathfrak{a}\_3)/2 \tag{32}
$$

that are used to construct a quadratic polynomial approximation *g*(*α*) of the form

$$J\left(\mathbf{x}\_{k'}\psi(\vec{u}\_k^{(j)} + a\vec{s}\_k^{(j)})\right) \approx \mathfrak{g}(\mathfrak{a}) := c\_0 + c\_1\mathfrak{a} + c\_2\mathfrak{a}^2. \tag{33}$$

The coefficients *c*0, *c*1, *c*<sup>2</sup> are obtained by solving the set of equations

$$J\left(\mathbf{x}\_{k\prime}\boldsymbol{\Psi}(\bar{\boldsymbol{u}}\_{k}^{(j)} + \boldsymbol{u}\_{l}\bar{\mathbf{s}}\_{k}^{(j)})\right) =: J\_{l} = \mathbf{g}\left(\boldsymbol{u}\_{l}\right), \quad i = 1, 2, 3\tag{34}$$

with the explicit solution

$$\begin{aligned} c\_{0} &= \frac{\mathfrak{a}\_{1} \left(\mathfrak{a}\_{1} - \mathfrak{a}\_{2}\right) \mathfrak{a}\_{2} \mathfrak{f}\_{3} + \mathfrak{a}\_{2} \mathfrak{a}\_{3} \left(\mathfrak{a}\_{2} - \mathfrak{a}\_{3}\right) \mathfrak{f}\_{1} + \mathfrak{a}\_{1} \mathfrak{a}\_{3} \left(\mathfrak{a}\_{3} - \mathfrak{a}\_{1}\right) \mathfrak{f}\_{2}}{\left(\mathfrak{a}\_{1} - \mathfrak{a}\_{2}\right) \left(\mathfrak{a}\_{1} - \mathfrak{a}\_{3}\right) \left(\mathfrak{a}\_{2} - \mathfrak{a}\_{3}\right)} \\ c\_{1} &= \frac{\left(\mathfrak{a}\_{2}^{2} - \mathfrak{a}\_{1}^{2}\right) \mathfrak{f}\_{3} + \left(\mathfrak{a}\_{1}^{2} - \mathfrak{a}\_{3}^{2}\right) \mathfrak{f}\_{2} + \left(\mathfrak{a}\_{3}^{2} - \mathfrak{a}\_{2}^{2}\right) \mathfrak{f}\_{1}}{\left(\mathfrak{a}\_{1} - \mathfrak{a}\_{2}\right) \left(\mathfrak{a}\_{1} - \mathfrak{a}\_{3}\right) \left(\mathfrak{a}\_{2} - \mathfrak{a}\_{3}\right)} \\ c\_{2} &= \frac{\left(\mathfrak{a}\_{1} - \mathfrak{a}\_{2}\right) \mathfrak{f}\_{3} + \left(\mathfrak{a}\_{2} - \mathfrak{a}\_{3}\right) \mathfrak{f}\_{1} + \left(\mathfrak{a}\_{3} - \mathfrak{a}\_{1}\right) \mathfrak{f}\_{2}}{\left(\mathfrak{a}\_{1} - \mathfrak{a}\_{2}\right) \left(\mathfrak{a}\_{1} - \mathfrak{a}\_{3}\right) \left(\mathfrak{a}\_{2} - \mathfrak{a}\_{3}\right)}}. \end{aligned} \tag{35}$$

If *c*<sup>2</sup> > 0, then the polynomial *g*(*α*) has a minimum at the point

$$
\hat{\mathfrak{a}} = -\frac{\mathfrak{c}\_1}{2\mathfrak{c}\_2} \,. \tag{36}
$$

If in addition *α*ˆ lies inside the interval [*α*1, *α*3], then *α*ˆ = *α* (*j*) *<sup>k</sup>* approximately solves the line search problem (28). Otherwise, *α* (*j*) *<sup>k</sup>* is set to one of the interval bounds *α*<sup>1</sup> or *α*3. In this case, the interval [*α*1, *α*3] can be adapted by a scaling factor to track the minimum point of the line search problem (28) over the single gradient iterations. Table 2 summarizes the overall algorithm for the approximate line search and the interval adaptation.

In general, the gradient method in Table 1 is stopped if the convergence criterion is fulfilled for some tolerance *ε<sup>J</sup>* > 0. In practice this can lead to a large number of iterations that moreover varies from one MPC iteration to the next. In order to ensure a real-time feasible MPC implementation, the gradient algorithm is stopped after *N* iterations and the re-initialization of the algorithm is done as outlined in Section 3.1.

#### 1) **Initialization: Default values and tolerances**


10 Will-be-set-by-IN-TECH

<sup>T</sup> in (28) represents a projection function of the form

*<sup>i</sup>* if *ui* < *u*<sup>−</sup>

*<sup>i</sup>* if *ui* <sup>&</sup>gt; *<sup>u</sup>*<sup>+</sup>

which guarantess the adherence of the input constraints [*u*−, *u*+]. For the real-time implementation within a suboptimal MPC scheme, the line search problem (28) can be solved

*<sup>k</sup>* (*τ*) and the step size *α*

The convergence properties of the gradient (projection) method are investigated, for instance, in Dunn (1996); Leese (1977); Nikol'skii (2007). In particular, Dunn (1996) proved under certain convexity and regularity assumptions that the gradient method exhibits a linear rate

The line search (28) represents a scalar optimization problem that is often solved approximately. The most straightforward way is to use a fixed step size *α* throughout all

An attractive alternative to a constant step size is to use a polynomial approximation with an

*i*

*i* ,

(*j*) *k* .

> � *xk*, *ψ*(*u*¯ (*j*) *<sup>k</sup>* + *αs*¯

<sup>≈</sup> *<sup>g</sup>*(*α*) :<sup>=</sup> *<sup>c</sup>*<sup>0</sup> <sup>+</sup> *<sup>c</sup>*1*<sup>α</sup>* <sup>+</sup> *<sup>c</sup>*2*α*<sup>2</sup> . (33)

=: *Ji* = *g*(*αi*), *i* = 1, 2, 3 (34)

*α*<sup>1</sup> < *α*<sup>2</sup> < *α*<sup>3</sup> with *α*<sup>2</sup> = (*α*<sup>1</sup> + *α*3)/2 (32)

*i* = 1, . . . , *m* (31)

(*j*+1)

*<sup>k</sup>* (*τ*), *τ* ∈ [0, *T*]

(*j*) *<sup>k</sup>* )

� in the

*ψi*(*ui*) =

⎧ ⎪⎨

*u*−

*u*+

in an approximate manner (see Section 4.2). Finally, the control trajectory *u*¯

gradient iterations. This, however, usually leads to a slow rate of convergence.

that are used to construct a quadratic polynomial approximation *g*(*α*) of the form

(*j*) *<sup>k</sup>* ) �

> (*j*) *<sup>k</sup>* ) �

*<sup>c</sup>*<sup>2</sup> <sup>=</sup> (*α*<sup>1</sup> <sup>−</sup> *<sup>α</sup>*2) *<sup>J</sup>*<sup>3</sup> <sup>+</sup> (*α*<sup>2</sup> <sup>−</sup> *<sup>α</sup>*3) *<sup>J</sup>*<sup>1</sup> <sup>+</sup> (*α*<sup>3</sup> <sup>−</sup> *<sup>α</sup>*1) *<sup>J</sup>*<sup>2</sup>

*<sup>c</sup>*<sup>0</sup> <sup>=</sup> *<sup>α</sup>*<sup>1</sup> (*α*<sup>1</sup> <sup>−</sup> *<sup>α</sup>*2) *<sup>α</sup>*<sup>2</sup> *<sup>J</sup>*<sup>3</sup> <sup>+</sup> *<sup>α</sup>*2*α*<sup>3</sup> (*α*<sup>2</sup> <sup>−</sup> *<sup>α</sup>*3) *<sup>J</sup>*<sup>1</sup> <sup>+</sup> *<sup>α</sup>*1*α*<sup>3</sup> (*α*<sup>3</sup> <sup>−</sup> *<sup>α</sup>*1) *<sup>J</sup>*<sup>2</sup> (*α*<sup>1</sup> − *α*2) (*α*<sup>1</sup> − *α*3) (*α*<sup>2</sup> − *α*3)

(*α*<sup>1</sup> <sup>−</sup> *<sup>α</sup>*2) (*α*<sup>1</sup> <sup>−</sup> *<sup>α</sup>*3) (*α*<sup>2</sup> <sup>−</sup> *<sup>α</sup>*3) .

underlying interval adaptation. To this end, the cost functional *J*

The coefficients *c*0, *c*1, *c*<sup>2</sup> are obtained by solving the set of equations

line search problem (28) is evaluated at three sample points

*J* � *xk*, *ψ*(*u*¯ (*j*) *<sup>k</sup>* + *αs*¯

*J* � *xk*, *ψ*(*u*¯ (*j*) *<sup>k</sup>* + *αis*¯

*c*<sup>1</sup> = � *α*2 <sup>2</sup> <sup>−</sup> *<sup>α</sup>*<sup>2</sup> 1 � *J*<sup>3</sup> + � *α*2 <sup>1</sup> <sup>−</sup> *<sup>α</sup>*<sup>2</sup> 3 � *J*<sup>2</sup> + � *α*2 <sup>3</sup> <sup>−</sup> *<sup>α</sup>*<sup>2</sup> 2 � *J*1 (*α*<sup>1</sup> <sup>−</sup> *<sup>α</sup>*2) (*α*<sup>1</sup> <sup>−</sup> *<sup>α</sup>*3) (*α*<sup>2</sup> <sup>−</sup> *<sup>α</sup>*3) (35)

(*j*)

*ui* else

⎪⎩

*ψ* = [*ψ*1,..., *ψm*]

follows from evaluating (29) with *s*¯

of convergence of the form (15).

**4.2 Adaptive line search**

with the explicit solution


$$\begin{array}{rcl} \text{if } c\_2 > c\_c: & \quad \quad a\_k^{(j)} = \begin{cases} \alpha\_1 & \text{if } \land < \alpha\_1\\ \alpha\_3 & \text{if } \land > \alpha\_3\\ \lor & \text{else} \end{cases} \end{array} \tag{37}$$

$$\text{else } (c\_2 \le \varepsilon\_\varepsilon): \quad a\_k^{(j)} = \begin{cases} \alpha\_1 & \text{if } f\_1 + \varepsilon\_\mathcal{S} \le \min\{f\_2, f\_3\} \\ \alpha\_3 & \text{if } f\_3 + \varepsilon\_\mathcal{S} \le \min\{f\_1, f\_2\} \\ \alpha\_2 & \text{else} \end{cases} \tag{38}$$

– Adapt the line search interval [*α*1, *α*3] for the next gradient iteration according to

$$\begin{array}{ll} \{a\_1, a\_3\} \leftarrow \begin{cases} \kappa^+[a\_1, a\_3] & \text{if } \Upk \ge a\_1 + \varepsilon\_{\mathfrak{a}}^+(a\_3 - a\_1) \text{ and } a\_3 \le a\_{\max} \\ \kappa^-[a\_1, a\_3] & \text{if } \Upk \le a\_1 + \varepsilon\_{\mathfrak{a}}^-(a\_3 - a\_1) \text{ and } a\_1 \ge a\_{\min}, \quad a\_2 \leftarrow \frac{a\_1 + a\_3}{2} \\ \ [a\_1, a\_3] & \text{else} \end{cases} \end{array} \tag{39}$$

Table 2. Adaptive line search for the gradient algorithm in Table 1.

Nonlinear Model Predictive Control 13

A Real-Time Gradient Method for Nonlinear Model Predictive Control 21

Pendulum link inner outer

length *li* [m] 0.323 0.480

mass *mi* [kg] 0.853 0.510 moment of inertia *Ji* [N m s2] 0.013 0.019 friction constant *di* [N m s] 0.005 0.005

distance to center of gravity *ai*

1, *φ*2, *φ*˙ 2]

> 0 Δ*x*¯

with Δ*x*¯ = *x*¯ − *xSP* and Δ*u*¯ = *u*¯ − *uSP* is used, which penalizes the distance to a desired setpoint (*xSP*, *uSP*), i.e. 0 = *f*(*xSP*, *uSP*). The symmetric and positive definite weighting

The CLF condition in Assumption 3 is approximately satisfied by solving the Riccati equation

to account for the fast dynamics of the double pendulum and the highly unstable behavior in

The suboptimal MPC scheme together with the gradient method were implemented as Cmex functions under MATLAB. The functions that are required in the gradient method are

<sup>4</sup> For the linearized (stabilizable) system Δ*x*˙ = *A*Δ*x* + *b*Δ*u*, it can be shown that the CLF inequality (9) is exactly fulfilled (in fact, (9) turns into an equality) for the terminal cost *V*(*x*) = Δ*x*T*P*Δ*x* and the linear (unconstrained) feedback law *<sup>q</sup>*(*x*) = <sup>−</sup>*R*−1*b*T*P*Δ*<sup>x</sup>* with *<sup>P</sup>* following from the Riccati equation (47).

the setpoint (*xSP*, *uSP*). <sup>4</sup> The sampling time Δ*t* and the prediction horizon *T* are set to

Table 3. Mechanical parameters of the double pendulum in Figure 2.

<sup>T</sup>(*T*)*P*Δ*x*¯(*T*) + *<sup>T</sup>*

For the MPC formulation, a quadratic cost functional (3a)

matrices *Q*, *R* in the integral part of (45) are chosen as

*xSP*,*uSP* and *<sup>b</sup>* <sup>=</sup> *<sup>∂</sup> <sup>f</sup>*

*∂u* 

[m]

With the state vector *x* = [*xc*, *x*˙*c*, *φ*1, *φ*˙

*J*(*xk*, *u*¯) = Δ*x*¯

the general nonlinear system

where *A* = *<sup>∂</sup> <sup>f</sup>*

*∂x* 

the inverted position.

**5.2 Simulation results**

*i* = 1 *i* = 2

0.215 0.223

<sup>T</sup>, the second-order ODEs (42) can be written as

<sup>T</sup>(*τ*)*R*Δ*u*¯(*τ*) d*τ* , (45)

*x*˙ = *f*(*x*, *u*), *x*(0) = *x*<sup>0</sup> . (44)

<sup>T</sup>(*τ*)*Q*Δ*x*¯(*τ*) + Δ*u*¯

*Q* = diag(10, 0.1, 1, 0.1, 1, 0.1), *R* = 0.001 . (46)

*PA* <sup>+</sup> *<sup>A</sup>*T*<sup>P</sup>* <sup>−</sup> *PbR*−1*b*T*<sup>P</sup>* <sup>+</sup> *<sup>Q</sup>* <sup>=</sup> 0 , (47)

*xSP*,*uSP* describe the linearization of the system (44) around

Δ*t* = 1 ms , *T* = 0.3 s (48)

Fig. 2. Inverted double pendulum on a cart.

#### **5. Example – Inverted double pendulum**

The inverted double pendulum on a cart is a benchmark problem in control theory due to its highly nonlinear and nonminimum-phase dynamics and its instability in the upward (inverted) position. The double pendulum in Figure 2 consists of two links with the lengths *li* and the angles *φi*, *i* = 1, 2 to the vertical direction. The displacement of the cart is given by *xc*. The mechanical parameters are listed in Table 3 together with their corresponding values (Graichen et al., 2007). The double pendulum is used in this section as benchmark example for the suboptimal MPC scheme and the gradient algorithm in order to show its performance for a real-time MPC implementation.

#### **5.1 Equations of motion and MPC formulation**

Applying the Lagrangian formalism to the double pendulum leads to the equations of motion (Graichen et al., 2007)

$$M(\phi)\ddot{\phi} + c(\phi, \dot{\phi}, \ddot{\mathfrak{x}}\_{\varepsilon}) = 0 \tag{40}$$

with the generalized coordinates *φ* = [*φ*1, *φ*2] <sup>T</sup> and the functions

$$M(\phi) = \begin{bmatrix} J\_1 + a\_1^2 m\_1 + l\_1^2 m\_2 & a\_2 l\_1 m\_2 \cos(\phi\_1 - \phi\_2) \\\\ a\_2 l\_1 m\_2 \cos(\phi\_1 - \phi\_2) & J\_2 + a\_2^2 m\_2 \end{bmatrix} \tag{41a}$$

$$\mathcal{L}(\boldsymbol{\phi}, \dot{\boldsymbol{\phi}}) = \begin{bmatrix} d\_1 \dot{\phi}\_1 + d\_2 (\dot{\phi}\_1 - \dot{\phi}\_2) + l\_1 m\_2 a\_2 \sin(\phi\_1 - \phi\_2) \dot{\phi}\_2^2 - (a\_1 m\_1 + l\_1 m\_2) \left[ \mathcal{g} \sin(\phi\_1) + \cos(\phi\_2) \ddot{x}\_c \right] \\\ d\_2 (\phi\_2 - \phi\_1) - a\_2 m\_2 \left[ \mathcal{g} \sin(\phi\_2) + l\_1 \sin(\phi\_1 - \phi\_2) \dot{\phi}\_1^2 + \cos(\phi\_2) \ddot{x}\_c \right] \end{bmatrix}. \tag{41b}$$

The acceleration of the cart *x*¨*c* serves as control input *u*. Thus, the overall model of the double pendulum can be written as the second-order ordinary differential equations (ODE)

$$\begin{aligned} \ddot{x}\_{\mathbb{C}} &= u \\ \ddot{\phi} &= -M^{-1}(\phi)c(\phi, \dot{\phi}, u) \end{aligned} \tag{42}$$

The acceleration of the cart is limited by the constraints

$$
\mu \in \left[ -6, +6 \right] \text{ m/s}^2 \text{ .} \tag{43}
$$

12 Will-be-set-by-IN-TECH

φ1(t)

φ2(t)

a1

0

<sup>1</sup>*m*<sup>2</sup> *a*2*l*1*m*<sup>2</sup> cos(*φ*<sup>1</sup> − *φ*2)

<sup>1</sup> <sup>−</sup> *<sup>φ</sup>*˙ <sup>2</sup>) + *<sup>l</sup>*1*m*2*a*<sup>2</sup> sin(*φ*<sup>1</sup> <sup>−</sup> *<sup>φ</sup>*2)*φ*˙ <sup>2</sup>

�

*x*¨*<sup>c</sup>* = *u*

<sup>2</sup>*m*<sup>2</sup>

pendulum can be written as the second-order ordinary differential equations (ODE)

m1, l1, J<sup>1</sup>

a2

m2, l2, J<sup>2</sup>

Fig. 2. Inverted double pendulum on a cart.

**5. Example – Inverted double pendulum**

performance for a real-time MPC implementation.

**5.1 Equations of motion and MPC formulation**

with the generalized coordinates *φ* = [*φ*1, *φ*2]

*<sup>a</sup>*2*l*1*m*<sup>2</sup> cos(*φ*<sup>1</sup> <sup>−</sup> *<sup>φ</sup>*2) *<sup>J</sup>*<sup>2</sup> <sup>+</sup> *<sup>a</sup>*<sup>2</sup>

*d*2(*φ*<sup>2</sup> − *φ*1) − *a*2*m*<sup>2</sup>

The acceleration of the cart is limited by the constraints

<sup>1</sup>*m*<sup>1</sup> + *l* 2

motion (Graichen et al., 2007)

*J*<sup>1</sup> + *a*<sup>2</sup>

<sup>1</sup> + *d*2(*φ*˙

*M*(*φ*) =

*c*(*φ*, *φ*˙) =

⎡ ⎣

� *d*1*φ*˙ xc

The inverted double pendulum on a cart is a benchmark problem in control theory due to its highly nonlinear and nonminimum-phase dynamics and its instability in the upward (inverted) position. The double pendulum in Figure 2 consists of two links with the lengths *li* and the angles *φi*, *i* = 1, 2 to the vertical direction. The displacement of the cart is given by *xc*. The mechanical parameters are listed in Table 3 together with their corresponding values (Graichen et al., 2007). The double pendulum is used in this section as benchmark example for the suboptimal MPC scheme and the gradient algorithm in order to show its

Applying the Lagrangian formalism to the double pendulum leads to the equations of

⎤

The acceleration of the cart *x*¨*c* serves as control input *u*. Thus, the overall model of the double

*<sup>g</sup>* sin(*φ*2) + *<sup>l</sup>*<sup>1</sup> sin(*φ*<sup>1</sup> <sup>−</sup> *<sup>φ</sup>*2)*φ*˙ <sup>2</sup>

<sup>T</sup> and the functions

u(t) = xc(t)

*M*(*φ*)*φ*¨ + *c*(*φ*, *φ*˙, *x*¨*<sup>c</sup>* ) = 0 (40)

⎦ (41a)

<sup>2</sup> − (*a*1*m*<sup>1</sup> + *l*1*m*2)[*g* sin(*φ*1) + cos(*φ*2)*x*¨*c*]

*<sup>φ</sup>*¨ <sup>=</sup> <sup>−</sup>*M*−1(*φ*)*c*(*φ*, *<sup>φ</sup>*˙, *<sup>u</sup>*). (42)

*<sup>u</sup>* <sup>∈</sup> [−6, <sup>+</sup>6] m/s<sup>2</sup> . (43)

<sup>1</sup> + cos(*φ*2)*x*¨*<sup>c</sup>*

�

� . (41b)


Table 3. Mechanical parameters of the double pendulum in Figure 2.

With the state vector *x* = [*xc*, *x*˙*c*, *φ*1, *φ*˙ 1, *φ*2, *φ*˙ 2] <sup>T</sup>, the second-order ODEs (42) can be written as the general nonlinear system

$$
\dot{\mathbf{x}} = f(\mathbf{x}, \boldsymbol{\mu}) \,, \quad \mathbf{x}(0) = \mathbf{x}\_0 \,. \tag{44}
$$

For the MPC formulation, a quadratic cost functional (3a)

$$J(\mathbf{x}\_k, \tilde{\mathbf{u}}) = \Delta \tilde{\mathbf{x}}^{\mathsf{T}}(T) P \Delta \tilde{\mathbf{x}}(T) + \int\_0^T \Delta \tilde{\mathbf{x}}^{\mathsf{T}}(\tau) Q \Delta \tilde{\mathbf{x}}(\tau) + \Delta \tilde{\mathbf{u}}^{\mathsf{T}}(\tau) R \Delta \tilde{\mathbf{u}}(\tau) \, d\tau,\tag{45}$$

with Δ*x*¯ = *x*¯ − *xSP* and Δ*u*¯ = *u*¯ − *uSP* is used, which penalizes the distance to a desired setpoint (*xSP*, *uSP*), i.e. 0 = *f*(*xSP*, *uSP*). The symmetric and positive definite weighting matrices *Q*, *R* in the integral part of (45) are chosen as

$$Q = \text{diag}(10, 0.1, 1, 0.1, 1, 0.1), \quad R = 0.001. \tag{46}$$

The CLF condition in Assumption 3 is approximately satisfied by solving the Riccati equation

*PA* <sup>+</sup> *<sup>A</sup>*T*<sup>P</sup>* <sup>−</sup> *PbR*−1*b*T*<sup>P</sup>* <sup>+</sup> *<sup>Q</sup>* <sup>=</sup> 0 , (47)

where *A* = *<sup>∂</sup> <sup>f</sup> ∂x xSP*,*uSP* and *<sup>b</sup>* <sup>=</sup> *<sup>∂</sup> <sup>f</sup> ∂u xSP*,*uSP* describe the linearization of the system (44) around the setpoint (*xSP*, *uSP*). <sup>4</sup> The sampling time Δ*t* and the prediction horizon *T* are set to

$$
\Delta t = 1 \text{ ms} \,, \quad T = 0.3 \text{ s} \tag{48}
$$

to account for the fast dynamics of the double pendulum and the highly unstable behavior in the inverted position.

#### **5.2 Simulation results**

The suboptimal MPC scheme together with the gradient method were implemented as Cmex functions under MATLAB. The functions that are required in the gradient method are

<sup>4</sup> For the linearized (stabilizable) system Δ*x*˙ = *A*Δ*x* + *b*Δ*u*, it can be shown that the CLF inequality (9) is exactly fulfilled (in fact, (9) turns into an equality) for the terminal cost *V*(*x*) = Δ*x*T*P*Δ*x* and the linear (unconstrained) feedback law *<sup>q</sup>*(*x*) = <sup>−</sup>*R*−1*b*T*P*Δ*<sup>x</sup>* with *<sup>P</sup>* following from the Riccati equation (47).

Nonlinear Model Predictive Control 15

A Real-Time Gradient Method for Nonlinear Model Predictive Control 23

MPC iter. *N* / CPU time [ms] / mean cost sampling step sampling step value [–]

Table 4. CPU time consumption of the real-time MPC scheme for different numbers of

between *N* = 1 and *N* = 10 iterations per sampling step are comparatively small compared

To investigate this point more precisely, Table 4 lists the required CPU time for different MPC settings. The computations were performed on a computer with an Intel i7 CPU (M620, 2.67 GHz) 5, 4 GB of memory, and the operating system MS Windows 7 (64 bit). The overall MPC scheme compiled as Cmex function under MATLAB 2010b (64 bit). All evaluated tests in Table 4 show that the required CPU times are well below the actual sampling time of Δ*t* = 1 ms. The CPU times are particularly remarkable in view of the high complexity of the nonlinear pendulum model (40)-(42), which illustrates the real-time feasibility of the suboptimal MPC scheme. The last column in Table 4 shows the average cost value that is obtained by integrating the cost profiles in Figure 3 and dividing through the simulation time of 5 s. This index indicates that the MPC scheme increases in terms of control performance for

From these numbers and the simulation profiles in Figure 3, the conclusion can be drawn that *N* = 2 gradient iterations per MPC step represents a good compromise between control performance and the low computational demand of approximately 100 *μ*s per MPC step.

Suboptimal solution strategies are efficient means to reduce the computational load for a real-time MPC implementation. The suboptimal solution from the previous MPC step is used for a warm-start of the optimization algorithm in the next run with the objective to reduce the suboptimality over the single MPC steps. Section 3 provides theoretical justifications for a

A suitable optimization algorithm is the gradient method in optimal control, which allows for a time and memory efficient calculation of the single MPC iterations and makes the overall MPC scheme suitable for very fast or high dimensional dynamical systems. The control performance and computational efficiency of the gradient method is illustrated in Section 5 for a highly nonlinear and complex model of a double pendulum on a cart. The suboptimal MPC scheme based on a real-time implementation of the gradient method was

suboptimal MPC scheme with a fixed number of iterations per sampling step.

<sup>5</sup> Only one core of the i7 CPU was used for the computations.

gradient iterations *N* per sampling step.

to the increase of numerical load.

larger numbers of *N*.

**6. Conclusions**

 0.053 0.0709 0.095 0.0641 0.133 0.0632 0.212 0.0610 0.405 0.0590

Fig. 3. MPC results for the double pendulum on a cart.

computed under the computer algebra system MATHEMATICA and are exported to MATLAB as optimized C code. The numerical integrations of the canonical equations (25)-(30) are performed by discretizing the time interval [0, *T*] with a fixed number of 30 equidistant points and using a second order Runge-Kutta method. The nonlinear model (44), respectively (41), is used within the MPC scheme as well as for the simulation of the pendulum.

The considered simulation scenario consists of an initial error around the origin (*xSP* = 0, *uSP* = 0) and a subsequent setpoint step of 1 m in the cart position at time *t* = 2s(*xSP* = [1 m, 0, 0, 0, 0, 0] <sup>T</sup>, *uSP* = 0). Figure 3 shows the simulation results for a two-stage scenario (initial error and setpoint change at *t* = 2 s). Already the case of one gradient iteration per sampling step (*N* = 1) leads to a good control performance and a robust stabilization of the double pendulum. Increasing *N* results in a more aggressive control behavior and a better exploitation of the control constraints (43).

The lower plots in Figure 3 show the (discrete-time) profiles of the cost value *J*(*xk*, *u*¯ (*N*) *<sup>k</sup>* ) and of the optimization error Δ*J*(*N*)(*xk*) = *J*(*xk*, *u*¯ (*N*) *<sup>k</sup>* ) <sup>−</sup> *<sup>J</sup>*∗(*xk*). In order to determine <sup>Δ</sup>*J*(*N*)(*xk* ), the optimal cost *J*∗(*xk*) was computed in each step *xk* by solving the OCP (3) for the double pendulum with a collocation-based optimization software. It is apparent from the respective plots in Figure 3 that the cost as well as the optimization error rapidly converge to zero which illustrates the exponential stability of the double pendulum in closed-loop and the incremental improvement of the algorithm. It is also seen in these plots that the performance improvement



Table 4. CPU time consumption of the real-time MPC scheme for different numbers of gradient iterations *N* per sampling step.

between *N* = 1 and *N* = 10 iterations per sampling step are comparatively small compared to the increase of numerical load.

To investigate this point more precisely, Table 4 lists the required CPU time for different MPC settings. The computations were performed on a computer with an Intel i7 CPU (M620, 2.67 GHz) 5, 4 GB of memory, and the operating system MS Windows 7 (64 bit). The overall MPC scheme compiled as Cmex function under MATLAB 2010b (64 bit). All evaluated tests in Table 4 show that the required CPU times are well below the actual sampling time of Δ*t* = 1 ms. The CPU times are particularly remarkable in view of the high complexity of the nonlinear pendulum model (40)-(42), which illustrates the real-time feasibility of the suboptimal MPC scheme. The last column in Table 4 shows the average cost value that is obtained by integrating the cost profiles in Figure 3 and dividing through the simulation time of 5 s. This index indicates that the MPC scheme increases in terms of control performance for larger numbers of *N*.

From these numbers and the simulation profiles in Figure 3, the conclusion can be drawn that *N* = 2 gradient iterations per MPC step represents a good compromise between control performance and the low computational demand of approximately 100 *μ*s per MPC step.

## **6. Conclusions**

−10 −5 0 5

−10

0

0.02

<sup>T</sup>, *uSP* = 0). Figure 3 shows the simulation results for a two-stage scenario

opt. error Δ

computed under the computer algebra system MATHEMATICA and are exported to MATLAB as optimized C code. The numerical integrations of the canonical equations (25)-(30) are performed by discretizing the time interval [0, *T*] with a fixed number of 30 equidistant points and using a second order Runge-Kutta method. The nonlinear model (44), respectively (41),

The considered simulation scenario consists of an initial error around the origin (*xSP* = 0, *uSP* = 0) and a subsequent setpoint step of 1 m in the cart position at time *t* = 2s(*xSP* =

(initial error and setpoint change at *t* = 2 s). Already the case of one gradient iteration per sampling step (*N* = 1) leads to a good control performance and a robust stabilization of the double pendulum. Increasing *N* results in a more aggressive control behavior and a better

(*N*)

the optimal cost *J*∗(*xk*) was computed in each step *xk* by solving the OCP (3) for the double pendulum with a collocation-based optimization software. It is apparent from the respective plots in Figure 3 that the cost as well as the optimization error rapidly converge to zero which illustrates the exponential stability of the double pendulum in closed-loop and the incremental improvement of the algorithm. It is also seen in these plots that the performance improvement

The lower plots in Figure 3 show the (discrete-time) profiles of the cost value *J*(*xk*, *u*¯

is used within the MPC scheme as well as for the simulation of the pendulum.

J(N) [-]

0.04

angle φ1 [deg]

0

10

angle φ2 [deg]

0 1 2 3 4

0 1 2 3 4

0 1 2 3 4

*<sup>k</sup>* ) <sup>−</sup> *<sup>J</sup>*∗(*xk*). In order to determine <sup>Δ</sup>*J*(*N*)(*xk* ),

time [s]

(*N*) *<sup>k</sup>* ) and

time [s]

time [s]

0 1 2 3 4

N = 1 N = 2 N = 5 N = 10

0 1 2 3 4

0 1 2 3 4

exploitation of the control constraints (43).

of the optimization error Δ*J*(*N*)(*xk*) = *J*(*xk*, *u*¯

time [s]

Fig. 3. MPC results for the double pendulum on a cart.

time [s]

time [s]

0

−5

0

[1 m, 0, 0, 0, 0, 0]

0.5

cost value

J [-] 1

control u [m/s2]

0

5

0.5

cart xc [m] 1

Suboptimal solution strategies are efficient means to reduce the computational load for a real-time MPC implementation. The suboptimal solution from the previous MPC step is used for a warm-start of the optimization algorithm in the next run with the objective to reduce the suboptimality over the single MPC steps. Section 3 provides theoretical justifications for a suboptimal MPC scheme with a fixed number of iterations per sampling step.

A suitable optimization algorithm is the gradient method in optimal control, which allows for a time and memory efficient calculation of the single MPC iterations and makes the overall MPC scheme suitable for very fast or high dimensional dynamical systems. The control performance and computational efficiency of the gradient method is illustrated in Section 5 for a highly nonlinear and complex model of a double pendulum on a cart. The suboptimal MPC scheme based on a real-time implementation of the gradient method was

<sup>5</sup> Only one core of the i7 CPU was used for the computations.

Nonlinear Model Predictive Control 17

A Real-Time Gradient Method for Nonlinear Model Predictive Control 25

*<sup>k</sup>* (*T*)) + *<sup>T</sup>*

0 *l*(*x*¯ ∗ *<sup>k</sup>* (*τ*), *u*¯<sup>∗</sup>

prove that Γ*<sup>α</sup>* contains the CLF region *Sβ*, consider *xk* ∈ *S<sup>β</sup>* and the bound on the optimal cost

with the CLF trajectories *x*¯*q*(*τ*), *x*¯*q*(0) = *xk*, and *u*¯*q*(*τ*) = *q*(*x*¯*q*(*τ*)). Similar to the proof of

Hence, (51)-(52) and definition (10) show that *J*∗(*xk*) ≤ *V*(*xk* ) ≤ *β* < *α* for all *xk* ∈ *Sβ*, which

0 *l*(*x*¯ *<sup>q</sup>*(*τ*), *u*¯

 *T* 0 *l*(*x*¯ *<sup>q</sup>*(*τ*), *u*¯

0

with the quadratic cost functional (53), the linear dynamics (54) and some initial state *<sup>x</sup>*<sup>0</sup> <sup>∈</sup> **<sup>R</sup>***n*. The admissible input set U[0,*T*] is assumed to be convex and the weighting matrices *P*, *Q*, *R* are symmetric and positiv definite. A useful property of the linear-quadratic problem is the strong

<sup>2</sup> [0,*T*] <sup>≤</sup> *<sup>J</sup>*(*u*) + *<sup>J</sup>*(*v*) <sup>−</sup> <sup>2</sup>*J*( <sup>1</sup>

for some constant *C* > 0 and all control functions *u*, *v* ∈ U[0,*T*]. To show this, first consider the control term of the cost functional (53) and the right-hand side of (55), which can be written

(*<sup>u</sup>* <sup>+</sup> *<sup>v</sup>*)T*R*(*<sup>u</sup>* <sup>+</sup> *<sup>v</sup>*) <sup>d</sup>*<sup>t</sup>* <sup>=</sup> <sup>1</sup>

subj. to *<sup>x</sup>*˙ <sup>=</sup> *Ax* <sup>+</sup> *Bu x*(0) = *<sup>x</sup>*<sup>0</sup> , *<sup>x</sup>* <sup>∈</sup> **<sup>R</sup>***n*, *<sup>u</sup>* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* (54)

*<sup>q</sup>*(*T*)) + *<sup>T</sup>*

*<sup>q</sup>*(0)) <sup>−</sup>

**8. Appendix B – Verification of Assumption 6 for linear-quadratic OCPs**

The following lines show that Assumption 6 is fulfilled for OCPs of the form

*<sup>k</sup>* (*t*)) d*τ*

= *α* . (50)

*<sup>k</sup>* (*T*) ∈/ *Sβ*. This, however, is a

*<sup>k</sup>* (*T*) ∈ *S<sup>β</sup>* for all *xk* ∈ Γ*α*. To

*<sup>q</sup>*(*τ*)) d*t* (51)

*<sup>q</sup>*(*τ*))d*τ* . (52)

*x*T(*t*)*Qx*(*t*) + *u*T(*t*)*Ru*(*t*) d*t* , (53)

(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*)T*R*(*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*) <sup>d</sup>*<sup>t</sup>* .

<sup>2</sup> *v*) (55)

<sup>2</sup> *<sup>u</sup>* <sup>+</sup> <sup>1</sup>

2 T 0

This allows one to derive a lower bound on the optimal cost

<sup>∗</sup>(*xk*) = *V*(*x*¯

≥ *β* +

= *β* 1 + *ml MV T* 

From this last line it can be concluded that *xk* ∈/ Γ*<sup>α</sup>* for all *x*¯<sup>∗</sup>

contradiction to the previous assumption and implies that *x*¯∗

<sup>∗</sup>(*xk*) ≤ *V*(*x*¯

*<sup>q</sup>*(*T*)) <sup>≤</sup> *<sup>V</sup>*(*x*¯

min*u*∈U[0,*T*] *<sup>J</sup>*(*u*) = *<sup>x</sup>*T(*T*)*Px*(*T*) + <sup>T</sup>

*<sup>C</sup>*||*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*||<sup>2</sup>

2

*<sup>u</sup>*T*Ru* <sup>+</sup> *<sup>v</sup>*T*Rv* <sup>−</sup> <sup>1</sup>

*Lm*

∗

 *T* 0 *ml β MV* d*τ*

*J*

*J*

*V*(*x*¯

Lemma 1, the CLF inequality (9) implies that

proves that Γ*<sup>α</sup>* contains *Sβ*.

convexity property (Allaire, 2007)

in the form

 T 0

also experimentally validated for a laboratory crane (Graichen et al., 2010) and for a helicopter with three degrees-of-freedom (Graichen et al., 2009), both experiments with sampling times of 1-2 milliseconds. The applicability of the gradient-based MPC scheme to high dimensional systems is demonstrated in (Steinböck et al., 2011) for a reheating furnace in steel industry.

#### **7. Appendix A – Reachability of CLF region (Theorem 1)**

This appendix proves the statements 1 and 2 in Theorem 1 concerning the reachability of the CLF region *S<sup>β</sup>* by the MPC formulation (3) without terminal constraints. The discrete-time case was investigated in Limon et al. (2006). The following two lemmas generalize these results to continuous-time systems as considered in this chapter. Lemma 1 represents an intermediate statement that is required to derive the actual result in Lemma 2.

**Lemma 1.** *Suppose that Assumptions 1-3 are satisfied. If x*¯∗ *<sup>k</sup>* (*T*) <sup>∈</sup>/ *<sup>S</sup><sup>β</sup> for any xk* <sup>∈</sup> **<sup>R</sup>***n, then x*¯∗ *<sup>k</sup>* (*τ*) ∈/ *S<sup>β</sup> for all times τ* ∈ [0, *T*]*.*

*Proof.* The proof is accomplished by contradiction. Assume that *x*¯∗ *<sup>k</sup>* (*T*) ∈/ *S<sup>β</sup>* and that there exists a time *τ*ˆ ∈ [0, *T*) such that *x*¯<sup>∗</sup> *<sup>k</sup>* (*τ*ˆ) ∈ *Sβ*. Starting at this point *x*¯<sup>∗</sup> *<sup>k</sup>* (*τ*ˆ), consider the residual problem

$$\hat{f}^\*(\mathfrak{x}\_k^\*(\mathfrak{t})) = \min\_{\mathfrak{u} \in \mathcal{U}\_{[T-\mathfrak{t}]}} \left\{ V\Big(\mathfrak{x}(T-\mathfrak{t}; \mathfrak{x}\_k^\*(\mathfrak{t}), \mathfrak{u}) \Big) + \int\_0^{T-\overline{\mathfrak{t}}} l\Big(\mathfrak{x}(\tau; \mathfrak{x}\_k^\*(\mathfrak{t}), \mathfrak{u}), \mathfrak{u}(\tau) \Big) \, \mathrm{d}\tau \right\}$$

subject to the dynamics (3b), for which the optimal trajectories are *u*ˆ∗(*τ*) = *u*¯∗ *<sup>k</sup>* (*τ*ˆ + *τ*) and *x*ˆ∗(*τ*) = *x*¯∗ *<sup>k</sup>* (*τ*ˆ + *τ*), *τ* ∈ [0, *T* − *τ*ˆ] by the principle of optimality. Since *x*¯<sup>∗</sup> *<sup>k</sup>* (*τ*ˆ) ∈ *S<sup>β</sup>* by assumption, the CLF inequality (9) with *x*¯*q*(0) = *x*¯<sup>∗</sup> *<sup>k</sup>* (*τ*ˆ) leads to the lower bound

$$\begin{aligned} V(\mathfrak{x}\_k^\*(\mathfrak{t})) &\geq V(\mathfrak{x}^q(T-\mathfrak{t})) + \int\_0^{T-\mathfrak{t}} l(\mathfrak{x}^q(\tau), \mathfrak{u}^q(\tau)) \mathrm{d}\tau \\ &\geq \widehat{f}^\*(\mathfrak{x}\_k^\*(\widehat{\tau})) \\ &\geq V(\widehat{\mathfrak{x}}^\*(T-\widehat{\tau})) = V(\widetilde{\mathfrak{x}}\_k^\*(T)) > \beta. \end{aligned}$$

The last line, however, implies that *x*¯∗ *<sup>k</sup>* (*τ*ˆ) ∈/ *Sβ*, which contradicts the previous assumption and thus proves the lemma.

**Lemma 2.** *Suppose that Assumptions 1-3 are satisfied and consider the compact set* Γ*<sup>α</sup> defined by* (10)*. Then, for all xk* ∈ Γ*α, the endpoint of the optimal state trajectory satisfies x*¯<sup>∗</sup> *<sup>k</sup>* (*T*) ∈ *Sβ. Moreover, S<sup>β</sup>* ⊆ Γ*α.*

*Proof.* We will again prove the lemma by contradiction. Assume that there exists a *xk* ∈ Γ*<sup>α</sup>* such that *x*¯∗ *<sup>k</sup>* (*T*) ∈/ *Sβ*, i.e. *V*(*x*¯<sup>∗</sup> *<sup>k</sup>* (*T*)) > *β*. Then, Lemma 1 states that *x*¯<sup>∗</sup> *<sup>k</sup>* (*τ*) ∈/ *S<sup>β</sup>* for all *τ* ∈ [0, *T*], or using (5),

$$\left|\left|\vec{x}\_{k}^{\*}(\tau)\right|\right|^{2} > \frac{\beta}{M\_{V}} \quad \forall \, \tau \in \left[0, T\right]. \tag{49}$$

16 Will-be-set-by-IN-TECH

also experimentally validated for a laboratory crane (Graichen et al., 2010) and for a helicopter with three degrees-of-freedom (Graichen et al., 2009), both experiments with sampling times of 1-2 milliseconds. The applicability of the gradient-based MPC scheme to high dimensional systems is demonstrated in (Steinböck et al., 2011) for a reheating furnace in steel industry.

This appendix proves the statements 1 and 2 in Theorem 1 concerning the reachability of the CLF region *S<sup>β</sup>* by the MPC formulation (3) without terminal constraints. The discrete-time case was investigated in Limon et al. (2006). The following two lemmas generalize these results to continuous-time systems as considered in this chapter. Lemma 1 represents an

*<sup>k</sup>* (*τ*ˆ) ∈ *Sβ*. Starting at this point *x*¯<sup>∗</sup>

 *T*−*τ*ˆ 0

∗

*<sup>k</sup>* (*T*)) > *β*. Then, Lemma 1 states that *x*¯<sup>∗</sup>

*<sup>k</sup>* (*T*)) > *β* .

*l x*¯(*τ*; *x*¯ ∗

*<sup>k</sup>* (*τ*ˆ) leads to the lower bound

*l*(*x*¯*q*(*τ*), *u*¯*q*(*τ*))d*τ*

*<sup>k</sup>* (*τ*ˆ) ∈/ *Sβ*, which contradicts the previous assumption

∀ *τ* ∈ [0, *T*] . (49)

*<sup>k</sup>* (*T*) <sup>∈</sup>/ *<sup>S</sup><sup>β</sup> for any xk* <sup>∈</sup> **<sup>R</sup>***n, then*

*<sup>k</sup>* (*τ*ˆ), *u*¯), *u*¯(*τ*)

*<sup>k</sup>* (*T*) ∈/ *S<sup>β</sup>* and that there

*<sup>k</sup>* (*τ*ˆ), consider the residual

 d*τ* 

*<sup>k</sup>* (*τ*ˆ + *τ*) and

*<sup>k</sup>* (*τ*ˆ) ∈ *S<sup>β</sup>* by

*<sup>k</sup>* (*T*) ∈ *Sβ. Moreover,*

*<sup>k</sup>* (*τ*) ∈/ *S<sup>β</sup>* for all

intermediate statement that is required to derive the actual result in Lemma 2.

**7. Appendix A – Reachability of CLF region (Theorem 1)**

**Lemma 1.** *Suppose that Assumptions 1-3 are satisfied. If x*¯∗

*Proof.* The proof is accomplished by contradiction. Assume that *x*¯∗

*x*¯(*T* − *τ*ˆ; *x*¯

∗ *<sup>k</sup>* (*τ*ˆ), *u*¯) + *T*−*τ*ˆ 0

*<sup>k</sup>* (*τ*ˆ + *τ*), *τ* ∈ [0, *T* − *τ*ˆ] by the principle of optimality. Since *x*¯<sup>∗</sup>

<sup>∗</sup>(*T* − *τ*ˆ)) = *V*(*x*¯

**Lemma 2.** *Suppose that Assumptions 1-3 are satisfied and consider the compact set* Γ*<sup>α</sup> defined by*

*Proof.* We will again prove the lemma by contradiction. Assume that there exists a *xk* ∈ Γ*<sup>α</sup>*

*MV*

*<sup>q</sup>*(*<sup>T</sup>* <sup>−</sup> *<sup>τ</sup>*ˆ)) +

subject to the dynamics (3b), for which the optimal trajectories are *u*ˆ∗(*τ*) = *u*¯∗

*<sup>k</sup>* (*τ*) ∈/ *S<sup>β</sup> for all times τ* ∈ [0, *T*]*.*

exists a time *τ*ˆ ∈ [0, *T*) such that *x*¯<sup>∗</sup>

*<sup>k</sup>* (*τ*ˆ)) = min

*u*¯∈U[*T*−*τ*ˆ]

*V*(*x*¯ ∗

*<sup>k</sup>* (*T*) ∈/ *Sβ*, i.e. *V*(*x*¯<sup>∗</sup>

The last line, however, implies that *x*¯∗

and thus proves the lemma.

 *V* 

assumption, the CLF inequality (9) with *x*¯*q*(0) = *x*¯<sup>∗</sup>

*<sup>k</sup>* (*τ*ˆ)) ≥ *V*(*x*¯

<sup>≥</sup> <sup>ˆ</sup>*<sup>J</sup>* ∗(*x*¯ ∗ *<sup>k</sup>* (*τ*ˆ))


≥ *V*(*x*ˆ

(10)*. Then, for all xk* ∈ Γ*α, the endpoint of the optimal state trajectory satisfies x*¯<sup>∗</sup>

*<sup>k</sup>* (*τ*)||<sup>2</sup> <sup>&</sup>gt; *<sup>β</sup>*

*x*¯∗

problem

ˆ*J* ∗(*x*¯ ∗

*x*ˆ∗(*τ*) = *x*¯∗

*S<sup>β</sup>* ⊆ Γ*α.*

such that *x*¯∗

*τ* ∈ [0, *T*], or using (5),

This allows one to derive a lower bound on the optimal cost

$$\begin{split} f^\*(\mathbf{x}\_k) &= V(\mathbf{\bar{x}}\_k^\*(T)) + \int\_0^T l(\mathbf{\bar{x}}\_k^\*(\tau), \bar{\mathbf{u}}\_k^\*(t)) \, \mathrm{d}\tau \\ &\ge \beta + \int\_0^T m\_l \frac{\beta}{M\_V} \mathrm{d}\tau \\ &= \beta \left( 1 + \frac{m\_l}{M\_V} T \right) \\ &= \alpha \,. \end{split} \tag{50}$$

From this last line it can be concluded that *xk* ∈/ Γ*<sup>α</sup>* for all *x*¯<sup>∗</sup> *<sup>k</sup>* (*T*) ∈/ *Sβ*. This, however, is a contradiction to the previous assumption and implies that *x*¯∗ *<sup>k</sup>* (*T*) ∈ *S<sup>β</sup>* for all *xk* ∈ Γ*α*. To prove that Γ*<sup>α</sup>* contains the CLF region *Sβ*, consider *xk* ∈ *S<sup>β</sup>* and the bound on the optimal cost

$$J^\*(\mathfrak{x}\_k) \le V(\mathfrak{x}^q(T)) + \int\_0^T l(\mathfrak{x}^q(\tau), \mathfrak{u}^q(\tau)) \, \mathrm{d}t \tag{51}$$

with the CLF trajectories *x*¯*q*(*τ*), *x*¯*q*(0) = *xk*, and *u*¯*q*(*τ*) = *q*(*x*¯*q*(*τ*)). Similar to the proof of Lemma 1, the CLF inequality (9) implies that

$$V(\mathfrak{x}^q(T)) \le V(\mathfrak{x}^q(0)) - \int\_0^T l(\mathfrak{x}^q(\tau), \mathfrak{u}^q(\tau)) d\tau. \tag{52}$$

Hence, (51)-(52) and definition (10) show that *J*∗(*xk*) ≤ *V*(*xk* ) ≤ *β* < *α* for all *xk* ∈ *Sβ*, which proves that Γ*<sup>α</sup>* contains *Sβ*.

#### **8. Appendix B – Verification of Assumption 6 for linear-quadratic OCPs**

The following lines show that Assumption 6 is fulfilled for OCPs of the form

$$\min\_{\mathbf{u}\in\mathcal{U}\_{[0,T]}} J(u) = \mathbf{x}^{\mathsf{T}}(T)\mathbf{P}\mathbf{x}(T) + \int\_{0}^{\mathsf{T}} \mathbf{x}^{\mathsf{T}}(t)\mathbf{Q}\mathbf{x}(t) + \mathbf{u}^{\mathsf{T}}(t)\mathbf{R}u(t)\,\mathrm{d}t\,,\tag{53}$$

$$\text{subj. to } \qquad \dot{\mathbf{x}} = A\mathbf{x} + Bu \quad \mathbf{x}(0) = \mathbf{x}\_0 \quad \mathbf{x} \in \mathbb{R}^n, \quad u \in \mathbb{R}^m \tag{54}$$

with the quadratic cost functional (53), the linear dynamics (54) and some initial state *<sup>x</sup>*<sup>0</sup> <sup>∈</sup> **<sup>R</sup>***n*. The admissible input set U[0,*T*] is assumed to be convex and the weighting matrices *P*, *Q*, *R* are symmetric and positiv definite. A useful property of the linear-quadratic problem is the strong convexity property (Allaire, 2007)

$$\left|\mathbb{C}||\mu - v||\_{L\_2^u[0,T]}^2 \le f(\mu) + f(v) - 2f(\frac{1}{2}\mu + \frac{1}{2}v) \tag{55}$$

for some constant *C* > 0 and all control functions *u*, *v* ∈ U[0,*T*]. To show this, first consider the control term of the cost functional (53) and the right-hand side of (55), which can be written in the form

$$\int\_0^\mathsf{T} \boldsymbol{u}^\mathsf{T} \boldsymbol{R} \boldsymbol{u} + \boldsymbol{v}^\mathsf{T} \boldsymbol{R} \boldsymbol{v} - \frac{1}{2} (\boldsymbol{u} + \boldsymbol{v})^\mathsf{T} \boldsymbol{R} (\boldsymbol{u} + \boldsymbol{v}) \, \mathrm{d}t = \frac{1}{2} \int\_0^\mathsf{T} (\boldsymbol{u} - \boldsymbol{v})^\mathsf{T} \boldsymbol{R} (\boldsymbol{u} - \boldsymbol{v}) \, \mathrm{d}t \, \boldsymbol{v} \, \mathrm{d}t$$

Nonlinear Model Predictive Control 19

A Real-Time Gradient Method for Nonlinear Model Predictive Control 27

Findeisen, R. (2006). *Nonlinear Model Predictive Control: A Sampled-Data Feedback Perspective*,

Fontes, F. (2001). A general framework to design stabilizing nonlinear model predictive

Fontes, F., Magni, L. & Gyurkovics, E. (2007). Sampled-data model predictive control

Graichen, K., Egretzberger, M. & Kugi, A. (2010). Suboptimal model predictive control of a

Graichen, K., Kiefer, T. & Kugi, A. (2009). Real-time trajectory optimization under input

Graichen, K. & Kugi, A. (2010). Stability and incremental improvement of suboptimal

Graichen, K., Treuer, M. & Zeitz, M. (2007). Swing-up of the double pendulum on a

Grüne, L. & Pannek, J. (2011). *Nonlinear Model Predictive Control: Theory and Algorithms*,

Hsu, J. & Meyer, A. (1968). *Modern Control Principles and Applications*, McGraw-Hill, New York. Ito, K. & Kunisch, K. (2002). Receding horizon optimal control for infinite dimensional systems, *ESAIM: Control, Optimisation and Calculus of Variations* 8: 741–760. Jadbabaie, A., Yu, J. & Hauser, J. (2001). Unconstrained receding horizon control of nonlinear

Kothare, S. D. O. & Morari, M. (2000). Contractive model predictive control for constrained nonlinear systems, *IEEE Transactions on Automatic Control* 45(6): 1053–1071.

Lee, Y., Kouvaritakis, B. & Cannon, M. (2002). Constrained receding horizon predictive control

Leese, S. (1977). Convergence of gradient methods for optimal control problems, *Journal of*

Limon, D., Alamo, T., Salas, F. & Camacho, E. (2006). On the stability of constrained MPC without terminal constraint, *IEEE Transactions on Automatic Control* 51(5): 832–836. Mayne, D., Rawlings, J., Rao, C. & Scokaert, P. (2000). Constrained model predictive control:

Michalska, H. & Mayne, D. (1993). Robust receding horizon control of constrained nonlinear

Nikol'skii, M. (2007). Convergence of the gradient projection method in optimal control

Ohtsuka, T. (2004). A continuation/GMRES method for fast computation of nonlinear

Parisini, T. & Zoppoli, R. (1995). A receding-horizon regulator for nonlinear systems and a

systems, *IEEE Transactions on Automatic Control* 38(11): 1623–1633.

problems, *Computational Mathematics and Modeling* 18(2): 148–156.

systems, *IEEE Transactions on Automatic Control* 46(5): 776–783.

Lee, E. & Markus, L. (1967). *Foundations of Optimal Control*, Wiley, New York.

for nonlinear systems, *Automatica* 38(12): 2093–2102.

*Optimization Theory and Applications* 21(3): 329–337.

stability and optimality, *Automatica* 36(6): 789–814.

receding horizon control, *Automatica* 40(4): 563–574.

neural approximation, *Automatica* 31(10): 1443–1451.

for nonlinear time-varying systems: Stability and robustness, *in* R. Findeisen, F. Allgöwer & L. Biegler (eds), *Assessment and Future Directions of Nonlinear Model*

laboratory crane, *8th IFAC Symposium on Nonlinear Control Systems (NOLCOS 2010)*,

constraints for a flatness-controlled laboratory helicopter, *European Control Conference*

MPC without terminal constraints, *IEEE Transactions on Automatic Control*

cart by feedforward and feedback control with experimental validation, *Automatica*

Vol. 1087, Fortschritt-Berichte VDI Reihe 8.

Budapest (Hungary).

55(11): 2576–2580.

Springer, London.

43(1): 63–71.

controllers, *Systems & Control Letters* 42(2): 127–143.

*(ECC) 2009*, Budapest (Hungary), pp. 2061–2066.

*Predictive Control*, LNCIS 358, Springer, Berlin, pp. 115–129.

The same simplifications can be used for the state-dependent terms in (53) since the linear dynamics (54) ensures that the superposition of two input signals *w*(*t*) = <sup>1</sup> <sup>2</sup> *<sup>u</sup>*(*t*) + <sup>1</sup> <sup>2</sup> *v*(*t*) yield a corresponding superposed state response *xw*(*t*) = <sup>1</sup> <sup>2</sup> *xu*(*t*) + <sup>1</sup> <sup>2</sup> *xv*(*t*) with *xw*(0) = *x*0. Hence, the right-hand side of (55) can be written as

$$\begin{aligned} f(u) + f(v) - 2f(\frac{1}{2}u + \frac{1}{2}v) &= \frac{1}{2} \Delta x^{\mathsf{T}}(T) P \Delta x(T) + \frac{1}{2} \int\_{0}^{\mathsf{T}} \Delta x^{\mathsf{T}}(t) Q \Delta x(t) \, \mathrm{d}t \\ &+ \frac{1}{2} \int\_{0}^{\mathsf{T}} \left( u(t) - v(t) \right)^{\mathsf{T}} R \left( u(t) - v(t) \right) \, \mathrm{d}t \\ &\ge C ||u - v||\_{L\_{2}^{\mathsf{u}}[0,T]}^{2} \end{aligned}$$

with Δ*x*(*t*) = *xu*(*t*) − *xv*(*t*) and the constant *C* = *λ*min(*R*)/2. Since *J*(*u*) is strongly (and therefore also strictly) convex on the convex set U[0,*T*] , it follows from standard arguments (Allaire, 2007) that there exists a global and unique minimum point *u*<sup>∗</sup> ∈ U[0,*T*]. Moreover, since <sup>U</sup>[0,*T*] is convex, <sup>1</sup> <sup>2</sup> (*<sup>u</sup>* <sup>+</sup> *<sup>u</sup>*∗) ∈ U[0,*T*] for all *<sup>u</sup>* ∈ U[0,*T*] such that *<sup>J</sup>*( <sup>1</sup> <sup>2</sup> *<sup>u</sup>* <sup>+</sup> <sup>1</sup> <sup>2</sup> *u*∗) ≥ *J*(*u*∗). Hence, the strong convexity inequality (55) can be turned into the quadratic growth property

$$\mathbb{C} \| |\boldsymbol{\mu} - \boldsymbol{\mu}^\*| \|\_{L^{\boldsymbol{w}}\_2[0,T]}^2 \le \boldsymbol{f}(\boldsymbol{\mu}) + \boldsymbol{f}(\boldsymbol{\mu}^\*) - \boldsymbol{2} \boldsymbol{f}(\frac{1}{2}\boldsymbol{\mu} + \frac{1}{2}\boldsymbol{\mu}^\*) \le \boldsymbol{f}(\boldsymbol{\mu}) - \boldsymbol{f}(\boldsymbol{\mu}^\*) \quad \forall \, \boldsymbol{\mu} \in \mathcal{U}\_{[0,T]}.$$

This shows that Assumption 6 is indeed satisfied for linear-quadratic OCPs of the form (53).

#### **9. Acknowledgements**

This work was supported by the Austrian Science Fund under project no. P21253-N22.

#### **10. References**

Allaire, G. (2007). *Numerical Analysis and Optimization*, Oxford University Press, New York. Berkovitz, L. (1974). *Optimal Control Theory*, Springer, New York.


18 Will-be-set-by-IN-TECH

The same simplifications can be used for the state-dependent terms in (53) since the linear

<sup>Δ</sup>*x*T(*T*)*P*Δ*x*(*T*) + <sup>1</sup>

*L<sup>m</sup>* <sup>2</sup> [0,*T*]

*u*(*t*) − *v*(*t*)

<sup>2</sup> (*<sup>u</sup>* <sup>+</sup> *<sup>u</sup>*∗) ∈ U[0,*T*] for all *<sup>u</sup>* ∈ U[0,*T*] such that *<sup>J</sup>*( <sup>1</sup>

<sup>2</sup> *xu*(*t*) + <sup>1</sup>

2 T 0

> <sup>T</sup> *R*

<sup>2</sup> *<sup>u</sup>*(*t*) + <sup>1</sup>

<sup>2</sup> *xv*(*t*) with *xw*(0) = *x*0. Hence,

 d*t*

, it follows from standard

<sup>2</sup> *<sup>u</sup>* <sup>+</sup> <sup>1</sup>

Δ*x*T(*t*)*Q*Δ*x*(*t*) d*t*

*u*(*t*) − *v*(*t*)

<sup>2</sup> *u*∗) ≤ *J*(*u*) − *J*(*u*∗) ∀ *u* ∈ U[0,*T*] .

<sup>2</sup> *v*(*t*) yield

.

<sup>2</sup> *u*∗) ≥

dynamics (54) ensures that the superposition of two input signals *w*(*t*) = <sup>1</sup>

<sup>2</sup> *<sup>v</sup>*) = <sup>1</sup> 2

> + 1 2 T 0

<sup>≥</sup> *<sup>C</sup>*||*<sup>u</sup>* <sup>−</sup> *<sup>v</sup>*||<sup>2</sup>

with Δ*x*(*t*) = *xu*(*t*) − *xv*(*t*) and the constant *C* = *λ*min(*R*)/2. Since *J*(*u*) is strongly

arguments (Allaire, 2007) that there exists a global and unique minimum point *u*<sup>∗</sup> ∈ U[0,*T*]

*J*(*u*∗). Hence, the strong convexity inequality (55) can be turned into the quadratic growth

This shows that Assumption 6 is indeed satisfied for linear-quadratic OCPs of the form (53).

This work was supported by the Austrian Science Fund under project no. P21253-N22.

Allaire, G. (2007). *Numerical Analysis and Optimization*, Oxford University Press, New York.

scheme with guaranteed stability, *Automatica* 34(10): 1205–1217.

Cannon, M. & Kouvaritakis, B. (2002). Efficient constrained model predictive control with asymptotic optimality, *SIAM Journal on Control and Optimization* 41(1): 60–82. Chen, H. & Allgöwer, F. (1998). A quasi-infinite horizon nonlinear model predictive control

DeHaan, D. & Guay, M. (2007). A real-time framework for model-predictive control

Diehl, M., Ferreau, H. & Haverbeke, N. (2009). Efficient numerical methods for nonlinear

Diehl, M., Findeisen, R., Allgöwer, F., Bock, H. & Schlöder, J. (2005). Nominal stability of

problems, *SIAM Journal on Control and Optimization* 34(4): 1270–1290.

of continuous-time nonlinear systems, *IEEE Transactions on Automatic Control*

MPC and moving horizon estimation, *in* L. Magni, D. Raimondo & F. Allgöwer (eds), *Nonlinear Model Predictive Control – Towards New Challenging Applications*, pp. 391–417.

real-time iteration scheme for nonlinear model predictive control, *IEE Proceedings*

<sup>2</sup> conditions amd the gradient projection method for optimal control

<sup>2</sup> *<sup>u</sup>* <sup>+</sup> <sup>1</sup>

a corresponding superposed state response *xw*(*t*) = <sup>1</sup>

<sup>2</sup> *<sup>u</sup>* <sup>+</sup> <sup>1</sup>

(and therefore also strictly) convex on the convex set U[0,*T*]

<sup>2</sup> [0,*T*] <sup>≤</sup> *<sup>J</sup>*(*u*) + *<sup>J</sup>*(*u*∗) <sup>−</sup> <sup>2</sup>*J*( <sup>1</sup>

Berkovitz, L. (1974). *Optimal Control Theory*, Springer, New York.

*Control Theory and Applications* 152(3): 296–308.

the right-hand side of (55) can be written as

*<sup>J</sup>*(*u*) + *<sup>J</sup>*(*v*) <sup>−</sup> <sup>2</sup>*J*( <sup>1</sup>

Moreover, since <sup>U</sup>[0,*T*] is convex, <sup>1</sup>

*L<sup>m</sup>*

52(11): 2047–2057.

Dunn, J. (1996). On *l*

*<sup>C</sup>*||*<sup>u</sup>* <sup>−</sup> *<sup>u</sup>*∗||<sup>2</sup>

**9. Acknowledgements**

**10. References**

property


**2** 

*Poland* 

Joanna Zietkiewicz

*Poznan University of Technology,* 

*Department of Control and Robotics,* 

**Feedback Linearization and LQ Based** 

Feedback linearization is a powerful technique that allows to obtain linear model with exact dynamics (Isidori,1985), (Slotine & Li, 1991). Linear quadratic control is well known optimal control method and with its dynamic programming properties can be also easily calculated (Anderson & Moore, 1990). The combination of feedback linearization and LQ control has been used in many algorithms in Model Predictive Control applications for many years and it is used also in the current papers (He De-Feng et al.,2011), (Margellos & Lygeros, 2010). Another problem apart from finding the optimal solution on a given horizon (finite or infinite) is the constrained control. A method which uses the advantages of feedback linearization, LQ control and applying signals constraints was proposed in (Poulsen et al., 2001b). In every step it is based on interpolation between the LQ optimal control and a feasible solution – the solution that fulfils given constraints. A feasible solution is obtained by taking calculated from LQ method optimal gain for a perturbed reference signal. The compromise between the feasible and optimal solution is calculating by minimization of one

variable – the number of degrees of freedom in prediction is reduced to one variable.

whole algorithm is operating in a discretized system.

Feedback linearization relies on choosing new state input and variables and then compensating nonlinearities in state equations by nonlinear feedback. The signals from nonlinear system are constrained, they are accessible from linear model through nonlinear equations. Therefore in the interpolation a nonlinear numerical method has to be used. The

There are several problems while using the method. One of them is that signals from nonlinear system can change its values within given one discrete time interval, while we assume that variables of linear model are unchanged. Those values should be considered as constrained. Another problem is finding the basic feasible perturbed reference signal which will provide well control performance. Method proposed in (Poulsen et. al, 2001b) gives good results if the weight matrices in cost function and the sampling interval are well chosen. Often it is difficult to choose these parameters and in general the solution may provide not only unfeasible signals (violating constraints), but also signals which violate assumption for system equations (like assumption of nonzero values in a denominator of a

**1. Introduction** 

fraction).

**Constrained Predictive Control** 

*Institute of Control and Information Engineering,* 

