Adaptive Horizon Model Predictive Control and Al'brekht's Method

TL;DR

This paper introduces Adaptive Horizon Model Predictive Control (AHMPC), which dynamically adjusts the prediction horizon to stabilize nonlinear control systems efficiently, leveraging Al'brekht's method for terminal cost and feedback.

Contribution

It proposes a novel adaptive horizon scheme for MPC that reduces computational complexity while ensuring stability, based on Al'brekht's method for terminal conditions.

Findings

AHMPC achieves stabilization with shorter horizons.

The method improves computational efficiency for nonlinear systems.

It extends MPC applicability to faster and more complex dynamics.

Abstract

A standard way of finding a feedback law that stabilizes a control system to an operating point is to recast the problem as an infinite horizon optimal control problem. If the optimal cost and the optmal feedback can be found on a large domain around the operating point then a Lyapunov argument can be used to verify the asymptotic stability of the closed loop dynamics. The problem with this approach is that is usually very difficult to find the optimal cost and the optmal feedback on a large domain for nonlinear problems with or without constraints. Hence the increasing interest in Model Predictive Control (MPC). In standard MPC a finite horizon optimal control problem is solved in real time but just at the current state, the first control action is implimented, the system evolves one time step and the process is repeated. A terminal cost and terminal feedback found by Al'brekht's methoddefined in a neighborhood of the operating point is used to shorten the horizon and thereby make the nonlinear programs easier to solve because they have less decision variables. Adaptive Horizon Model Predictive Control (AHMPC) is a scheme for varying the horizon length of Model Predictive Control (MPC) as needed. Its goal is to achieve stabilization with horizons as small as possible so that MPC methods can be used on faster and/or more complicated dynamic processes.

Equations

$$x^{+}$$

$$x$$

$$u$$

$$y=h(x,u)$$

$$u\in \mathbb{U},$$

$$\kappa (x)\in \mathbb{U},$$

$$V(x^0)=\underset{\mathbf{u}}{\min }\sum _{k=0}^{\infty }l(x(k),u(k))$$

$$\begin{array}[]{rllllllllllllllllllllllllllllllll}V(x)&=&\mbox{min}_{u}\left\{V(f(x,u))+l(x,u)\right\}\\ \kappa(x)&=&\mbox{argmin}_{u}\left\{V(f(x,u))+l(x,u)\right\}\\ \end{array}$$

$$\kappa (x)\in \mathbb{U},$$

$$\sum _{k=0}^{N-1}l(x(k),u(k))+V_f(x(N))$$

$${\mathbf{u}}_N=(u_N(0),u_N(1),\dots ,u_N(N-1))$$

$$V_N(x^0)$$

$${\mathbf{u}}_N^{\ast }$$

$${\mathbf{x}}_N^{\ast }$$

$$V_f(x)>V_f(f(x,{\kappa }_f(x)))\ge 0$$

$$x_N^{\ast }(k+1))$$

$${\kappa }_f(x_N^{\ast }(k))$$

$$f(x_N^{\ast }(k),{\kappa }_f(x_N^{\ast }(k)))$$

$$h(x_N^{\ast }(k),{\kappa }_f(x_N^{\ast }(k)))$$

$$V_f(x_N^{\ast }(k))$$

$$V_f(x_N^{\ast }(k))-V_f(x_N^{\ast }(k+1))$$

$$l(x,u)$$

$$f(x,u)$$

$$P$$

$$K$$

$$x^{+}$$

$$V(x)$$

$$0$$

$$f(x,u)$$

$$l(x,u)$$

$$V(x)$$

$$\kappa (x)$$

$$V^{[3]}(x)-V^{[3]}((F+GK)x)$$

$${\left({\kappa }^{[2]}(x)\right)}^{\prime }\left(R+G^{\prime }PG\right)$$

Full text

Adaptive Horizon Model Predictive Control (AHMPC) is a scheme for varying the horizon length of Model Predictive Control (MPC) as needed. Its goal is to achieve stabilization with horizons as small as possible so that MPC methods can be used on faster and/or more complicated dynamic processes.

1 Introduction

Model Predictive Control (MPC) is a way to steer a discrete time control system to a desired operating point. We will present an extension of MPC that we call Adaptive Horizon Model Predictive Control (AHMPC) which adjusts the length of the horizon in MPC while confirming in real time that stabilization is occuring for a nonlinear system.

We are not the first to consider adaptively changing the horizon length in MPC, see [11], [14], [15], [16]. In these papers the horizon is changed so that a terminal constraint is satisfied by the predicted state at the end of horizon. In [5] the horizon length is adaptively changed to ensure that the infinite horizon cost of using the finite horizon MPC scheme is not much more than the cost of the corresponding infinite horizon optimal control problem.

Adaptive horizon tracking is discussed in [12] and [4]. In [10] an adaptive parameter estimation algorithm suitable for MPC was proposed, which uses the available input and output signals to estimate the unknown system parameters. In [7] is a detailed analysis of the impact of the optimization horizon and the time varying control horizon on stability and performance of the closed loop is given.

2 Review of Model Predictive Control

We briefly describe MPC following the definitive treatise of [17]. We largely follow their notation.

We are given a controlled, nonlinear dynamics in discrete time

[TABLE]

where the state $x\in{I\!\!R}^{n\times 1}$, the control $u\in{I\!\!R}^{m\times 1}$ and $x^{+}(k)=x(k+1)$. Typically this a discretization of a controlled, nonlinear dynamics in continuous time. The goal is to find a feedback law $u(k)=\kappa(x(k))$ that drives the state of the system to some desired operating point. A pair $(x^{e},u^{e})$ is an operating point if $f(x^{e},u^{e})=x^{e}$. We conveniently assume that, after state and control coordinate translations, the operating point of interest is $(x^{e},u^{e})=(0,0)$.

The controlled dynamics may be subject to constraints such as

[TABLE]

and possibly constraints involving both the state and control

[TABLE]

A control $u$ is said to be feasible at $x\in\mathbb{X}$ if

[TABLE]

Of course the stabilizing feedback $\kappa(x)$ that we seek needs to be feasible, that is, for every $x\in\mathbb{X}$,

[TABLE]

An ideal way to find a stabilizing feedback is to choose a Lagrangian $l(x,u)$ (aka running cost) that is nonnegative definite in $x,u$, positive definite in $u$ and then to solve the infinite horizon optimal control problem of minimzing the quantity $\sum_{k=0}^{\infty}l(x(k),u(k))$ over all choices of infinite control sequences

${\bf u}=(u(0),u(1),\ldots)$ subject to the dynamics (1), the constraints (2, 3, 4) and the initial condition $x(0)=x^{0}$. Assuming the minimum exists for each $x^{0}\in\mathbb{X}$, we define the optimal cost function

[TABLE]

Let ${\bf u}^{*}=(u^{*}(0),u^{*}(1),\ldots)$ be a minimzing control sequence with corresponding state sequence ${\bf x}^{*}=(x^{*}(0)=x^{0},x^{*}(1),\ldots)$. Minimizing control and state sequences need not be unique but we shall generally ignore this problem because we are using optimization as a path to stabilization. The key question is whether the possibly nonunique solution is stabilizing to the desired operating point. As we shall see AHMPC nearly verifies stabilization in real time.

If a pair $V(x)\in{I\!\!R},\kappa(x)\in{I\!\!R}^{m\times 1}$ of functions satisfy the infinite horizon Bellman Dynamic Program Equations (BDP)

[TABLE]

and the constraints

[TABLE]

for all $x\in\mathbb{X}$ then it is not hard to show that $V(x)$ is the optimal cost and $\kappa(x)$ is an optimal feedback on $\mathbb{X}$. Under suitable conditions a Lyapunov argument can be used to show that the feedback $\kappa(x)$ is stabilizing.

The difficulty with this approach is that it is generally impossible to solve the BDP equations on a large domain $\mathbb{X}$ if the state dimension $n$ is greater than $2$ or $3$. So both theorists and practicioners have turned to Model Predictive Control (MPC).

In MPC one chooses a Lagrangian $l(x,u)$, a horizon length $N$, a terminal domain $\mathbb{X}_{f}\subset\mathbb{X}$ containing $x=0$ and a terminal cost $V_{f}(x)$ defined and positive definite on $\mathbb{X}_{f}$. Then one considers the problem of minimizing

[TABLE]

by choice of feasible

[TABLE]

subject to the dynamics (1), the constraints (8), the final condition $x(N)\in\mathbb{X}_{f}$ and the initial condition $x(0)=x^{0}$. Assuming this problem is solvable, let $V_{N}(x^{0})$ denote the optimal cost,

[TABLE]

and let

[TABLE]

denote optimal control and state sequences when the horizon length is $N$. We then define the MPC feedback law $\kappa_{N}(x^{0})=u_{N}^{*}(0)$.

The terminal set $\mathbb{X}_{f}$ is controlled invariant (aka viable) if for each $x\in\mathbb{X}_{f}$ there exists a $u\in\mathbb{U}$ such that $f(x,u)\in\mathbb{X}_{f}$ and the constraint (4) is satisfied. If $V_{f}(x)$ is a control Lyapunov function on $\mathbb{X}_{f}$ then, under suitable conditions, a Lyapunov argument can be used to show that the feedback $\kappa_{N}(x)$ is stabilizing on $\mathbb{X}_{f}$. See [17] for more details.

AHMPC requires a little more, the existence of terminal feedback $u=\kappa_{f}(x)$ defined on the terminal set $\mathbb{X}_{f}$ that leaves it positively invariant, if $x\in\mathbb{X}_{f}$ then $f(x,\kappa_{f}(x))\in\mathbb{X}_{f}$, and which makes $V_{f}(x)$ a strict Lyapunov function on $\mathbb{X}_{f}$ for the closed loop dynamics, if $x\in\mathbb{X}_{f}$ and $x\neq 0$ then

[TABLE]

If $x\in\mathbb{X}_{f}$ then AHMPC does not need to solve (9) to get $u$, it just takes $u=\kappa_{f}(x)$. A similar scheme has been called dual mode control in [11].

The advantage of solving the finite horizon optimal control problem (9) over solving the infinite horizon problem (6) is that it may be possible to solve the former on-line as the process evolves. If it is known that the terminal set $\mathbb{X}_{f}$ can be reached from the current state $x$ in $N$ or fewer steps then the finite horizon $N$ optimal control problem is a feasible nonlinear program with finite dimenionsal decision variable ${\bf u}_{N}\in{I\!\!R}^{m\times N}$. If the time step is long enough, if $f,h,l$ are reasonably simple and if $N$ is small enough then this nonlinear program possibly can be solved in a fraction of one time step for ${\bf u}^{*}_{N}$. Then the first element of this sequence $u_{N}^{*}(0)$ is used as the control at the current time. The system evolves one time step and the process is repeated at the next time. Conceptually MPC computes a feedback law $\kappa_{N}(x)=u^{*}_{N}(0)$ but only at values of $x$ when and where it is needed.

Some authors like [6], [8] do away with the terminal cost $V_{f}(x)$ but there is a theoretical reason and a practical reason to use one. The theoretical reason is that a control Lyapunov terminal cost facilitates a proof of asymptotic stability via a simple Lyapunov argument, see [17]. But this is not a binding reason because under suitable assumptions, asymptotic stability can be shown even when there is no terminal cost provided the horizon is sufficiently long. The practical reason is more important, when there is a terminal cost one can usually use a shorter horizon $N$. A shorter horizon reduces the dimension $mN$ of the decision variables in the nonlinear programs that need to be solved on-line. Therefore MPC with a suitable terminal cost can be used for faster and more complicated systems.

An ideal terminal cost $V_{f}(x)$ is $V(x)$ of the corresponding infinite horizon optimal control provided that the latter can be accurately computed off-line on a reasonably large terminal set $\mathbb{X}_{f}$. For then the infinite horizon cost (6) and (9) will be the same. One should not make too much of this fact as stabilization is our goal, the optimal control problems are just a means to accomplish this. This in contrast to Economic MPC where the cost and the associated Lagrangian are chosen to model real world costs.

3 Adaptive Horizon Model Predictive Control

For AHMPC we assume that we have the following.

• •

A discrete time dynamics $f(x,u)$ with operating point $x=0,u=0$.

• •

A Lagrangian $l(x,u)$, nonegative definite in $(x,u)$ and positive definite in $u$.

• •

State constraints $x\in\mathbb{X}$ where $\mathbb{X}$ is a neighborhood of $x=0$.

• •

Control constraints $u\in\mathbb{U}$ where $\mathbb{U}$ is a neighborhood of $u=0$.

• •

Mixed constraints $h(x,u)\in\mathbb{Y}$ which are not active at the operating point $x=0,u=0$.

• •

The dynamics is recursively feasible on $\mathbb{X}$, that is, for every $x\in\mathbb{X}$ there is a $u\in\mathbb{U}$ satisfying $h(x,u)\in\mathbb{Y}$ and $f(x,u)\in\mathbb{X}$

• •

A terminal cost $V_{f}(x)$ defined and nonnegative definite on some neighborhood $\mathbb{X}_{f}$ of the operating point $x=0,u=0$. The neighborhood $\mathbb{X}_{f}$ need not be known explicitly.

• •

A terminal feedback $u=\kappa_{f}(x)$ defined on $\mathbb{X}_{f}$ such that the terminal cost is a valid Lyapunov function on $\mathbb{X}_{f}$ for the closed loop dynamics using the terminal feedback $u=\kappa_{f}(x)$.

One way of obtaining a terminal pair $V_{f}(x),\ \kappa_{f}(x)$ is to approximately solve the infinite horizon dynamic program equations (BDP) on some neighborhood of the origin. For example if the linear part of the dynamics and the quadratic part of the Lagrangian constitute a LQR problem satisfying the standard conditions then one can let $V_{f}(x)$ be the quadratic optimal cost and $\kappa_{f}(x)$ be the linear optimal feedback of the LQR problem. Of course the problem with such terminal pairs $V_{f}(x),\ \kappa_{f}(x)$ is that generally there is no way to estimate the terminal set $\mathbb{X}_{f}$ on which the feasibility and Lyapunov conditions are satisfied. It is reasonable to expect that they are satisfied on some terminal set but the extent of this terminal set is very difficult to estimate.

In the next section we show how higher degree Taylor polynomials for the optimal cost and optimal feedback can be computed by the discrete time version of [3] found in [1] because this can lead to a larger terminal set $\mathbb{X}_{f}$ on which the feasibility and the Lyapunov conditions are satisfied. It would be very difficult to determine what this terminal set is but we do not need to do this.

AHMPC mitigates this last difficulty just as MPC mitigates the problem of solving the infinite horizon Bellman Dynamic Programming equations BDP (7). MPC does not try to compute the optimal cost and optimal feedback everywhere, instead it computes them just when and where they are needed. AHMPC does not try to compute the set $\mathbb{X}_{f}$ on which $\kappa_{f}(x)$ is feasible and stabilizing, it just tries to determine if the end state $x_{N}^{*}(N)$ of the currently computed optimal trajectory is in a terminal set $\mathbb{X}_{f}$ where the feasibility and Lyapunov conditions are satisfied.

Suppose the current state is $x$ and we have solved the horizon $N$ optimal control problem for ${\bf u}^{*}_{N}=(u^{*}_{N}(0),\ldots,u^{*}_{N}(N-1))$, ${\bf x}^{*}_{N}=(x^{*}_{N}(0)=x,\ldots,x^{*}_{N}(N))$ The terminal feedback $u=\kappa_{f}(x)$ is used to compute $M$ additional steps of this state trajectory

[TABLE]

for $k=N,\ldots,N+M-1$.

Then one checks that the feasibility and Lyapunov conditions hold for the extended part of the state sequence,

[TABLE]

for $k=N,\ldots,N+M-1$ and some Class K function $\alpha(s)$. For more on Class K functions we refer the reader to [9].

If (11-15) hold for all for $k=N,\ldots,N+M-1$. then we presume that $x^{*}_{N}(N)\in\mathbb{X}_{f}$, a set where $\kappa_{f}(x)$ is stabilizing and we use the control $u_{N}^{*}(0)$ to move one time step forward to $x^{+}=f(x,u^{*}_{N}(0))$. At this next state $x^{+}$ we solve the horizon $N-1$ optimal control problem and check that the extension of the new optimal trajectory satisfies (11-15).

If (11-15) does not hold for all for $k=N,\ldots,N+M-1$. then we presume that $x^{*}_{N}(N)\notin\mathbb{X}_{f}$. We extend current horizon by $L\geq 1$ and if time permits we solve the horizon $N+L$ optimal control problem at the current state $x$ and then check the feasibility and Lyapunov conditions again. We keep increasing $N$ by $L$ until these conditions are satisfied on the extension of the trajectory. If we run out of time before they are satisfied then we use the last computed $u_{N}^{*}(0)$ and move one time step forward to $x^{+}=f(x,u^{*}_{N}(0))$. At $x^{+}$ we solve the horizon $N+L$ optimal control problem.

How does one choose the extended horizon $M$ and the class ${\it K}$ function $\alpha(\cdot)$? If the extended part of the state sequence is actually in the region where the terminal cost $V_{f}(x)$ and the terminal feedback $\kappa_{f}(x)$ well approximate the solution to infinite horizon optimal control problem then the dynamic programing equations (7) should approximately hold. In other words

[TABLE]

If this does not hold throughout the extended trajectory we should increase both the horizon $N$. We can also increase the extended horizon $M$ but this may not be necessary. If the Lyapunov and feasibilty conditions are going to fail somewhere on the extension it is most likely this will happen at the beginning of the extension. Also we should choose $\alpha(\cdot)$ so that $\alpha(|x|)<|l(x,\kappa_{f}(x))|/2$.

The nonlinear programming problems generated by employing MPC on a nonlinear system are generally nonconvex so the solver might return local rather than global minimizers. In which case there is no guarantee that an MPC approach is actually stabilizing. AHMPC mitigates this difficulty by checking that stabilization is occuring. If (11-15) don’t hold even after the horizon $N$ has been increased substantially then this is a strong indication that the solver is returning locally rather than globally minimizing solutions and these local solutions are not stabilizing. To change this behaviour one needs to start the solver at a substantially different initial guess. Just how one does this is an open research question. It is essentially the same question as which initial guess should one pass to the solver at the first step of MPC.

The actual computation of the $M$ additional steps (10) can be done very quickly because the closed loop dynamics function $f(x,\kappa_{f}(x))$ can be computed and compiled before hand. Similarly the feasibility and Lyapunov conditions (11-15) can be computed and compiled before hand. The number $M$ of additional time steps is a design parameter. One choice is to take $M$ a fraction of the current $N$.

4 Choosing a Terminal Cost and a Terminal Feedback

A standard way of obtaining a terminal cost $V_{f}(x)$ and a terminal feedback $\kappa_{f}(x)$ is to solve the Linear Quadratic Regulator (LQR) using the quadratic part of the Lagrangian and the linear part of dynamics around the operating point $(x^{e},u^{e})=(0,0)$. Suppose

[TABLE]

Then the LQR problem is to find $P,K$ such that

[TABLE]

Under mild assumptions, the stabilizability of $F,G$, the detectability of $Q^{1/2},F$, the nonnegative definiteness of $[Q,S;S^{\prime},R]$ and the positive definiteness of $R$, there exist an unique nonnegative definite $P$ satisfying the first equation (LABEL:P) which is called the discrete time algebraic Riccati equation. Then $K$ given by (17) puts all the poles of the closed loop linear dynamics

[TABLE]

inside of the open unit disk. See [2] for details. But $l(x,u)$ is a design parameter so if we choose $[Q,S;S^{\prime},R]$ to be positive definite and then $P$ will be positive definite.

If we define the terminal cost to be $V_{f}(x)={1\over 2}x^{\prime}Px$ then we know that it is positive definite for all $x\in{I\!\!R}^{n\times 1}$. If we define the terminal feedback to be $\kappa_{f}(x)=Kx$ then we know by a Lyapunov argument that the nonlinear closed loop dynamics

[TABLE]

is locally asymptotically stable around $x^{e}=0$. The problem is that we don’t know what is the neighborhood $\mathbb{X}_{f}$ of asymptotic stability and computing it off-line can be very difficult in state dimensions higher that two or three.

There are other possible choices for the terminal cost and terminal feedback. Al’brekht [3] showed how the Taylor polynomials of the optimal cost and the optimal feedback could be computed for some smooth, infinite horizon optimal control problems in continuous time. Aguilar and Krener [1] extended this to some smooth, infinite horizon optimal control problems in discrete time. The discrete time Taylor polynomials of the optimal cost and the optimal feedback may be used as the terminal cost and terminal feedback in an AHMPC scheme so we briefly review [1].

Since we assumed that $f(x,u)$ and $l(x,u)$ are smooth and the constraints are not active at the origin, we can simplify the BDP equations. The simplified Bellman Dynamic Programming equations (sBDP) are obtained by setting the derivative with respect to $u$ of the quantity to be minimized in (7) to zero. The result is

[TABLE]

If the quantity to be minimized is strictly convex in $u$ then the BDP equations and sBDP equations are equivalent. But if not then BDP implies sBDP but not vice versa.

Suppose the discrete time dynamics and Lagrangian have Taylor polynomials around the operating point $x=0,u=0$ of the form

[TABLE]

for some integer $d\geq 1$ and where [j] indicates homogeneous polynomial terms of degree $j$.

Also suppose the infinite horizon optimal cost and optimal feedback have similar Taylor polynomials

[TABLE]

We plug these polynomials into sBDP and collect terms of lowest degree. The lowest degree in (18) is two while in (19) it is one. The result is the discrete time Ricatti equations (LABEL:P, 17).

At the next degrees we obtain the equations

[TABLE]

Notice these are linear equations in the unknowns $V^{[3]}(x)$ and $\kappa^{[2]}(x)$ and the right sides of these equations involve only known quantities. Moreover $\kappa^{[2]}(x)$ does not appear in the first equaion. The eigenvalues of the linear operator

[TABLE]

are of the form of $1$ minus the products of three eigenvalues of $F+GK$. Since the eigenvalues of $F+GK$ are all inside the open unit disk this operator (22) is invertible. Having solved (20) for $V^{[3]}(x)$ we can readily solve (21) for $\kappa^{[2]}(x)$ since we have assumed $R$ is positive definite

The higher degree equations are similar, at degrees $d+1,d$ they take the form

[TABLE]

The ”Known Quantities” involve the terms of the Taylor polynomials of $f,l$ and previously computed $V^{[i+1]},\kappa^{[i]}$ for $1\leq i<j$. Again the equations are linear in the unknowns $V^{[k+1]},\kappa^{[k]}$ and the first equation does not involve $\kappa^{[k]}$. The eigenvalues of the linear operator

[TABLE]

are of the form of $1$ minus the products of $k+1$ eigenvalues of $F+GK$. We have written MATLAB code to solve these equations to any degree and in any dimensions. The code is quite fast. It found the Taylor polynomials of the optimal cost to degree 6 and the optimal feedback to degree 5 in 0.12 sec. on a laptop using 3.1 GHz Intel Core i5.

5 Completing the Squares

The infinite horizon optimal cost is certainly nonnegative definite and if we choose $Q>0$ then it is positive definite. That implies that its quadratic part $V^{[2]}(x)={1\over 2}x^{\prime}Px$ is positive definite.

But its Taylor polynomial of degree $d+1$,

[TABLE]

need not be positive definite for $d>1$. This can lead to problems if we define this Taylor polynomial to be our terminal cost $V_{f}(x)$ because then the nonlinear program solver might return a negative cost $V_{N}(x)$ (9). The way around this difficulty is to ”complete the squares”.

Theorem Suppose a polynomial $V(x)$ is of degrees two through $d+1$ in $n$ variables $x_{1},\ldots,x_{n}$. If the quadratic part of $V(x)$ is positive definite then there exist a nonnegative definite polynomial $W(x)$ of degrees two through $2d$ such that the part of $W(x)$ that is of degrees two through $d+1$ equals $V(x)$. Moreover, we know that $W(x)$ is nonegative definite because it is the sum of $n$ squares.

Proof: We start with the quadratic part of $V(x)$, because it is positive definite it must be of the form ${1\over 2}x^{\prime}Px$ where $P$ is a positive definite $n\times n$ matrix. We know that there is an orthogonal matrix $T$ that diagonalizes $P$

[TABLE]

where $\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{n}>0$. We make the linear change of coordinates $x=Tz$. We shall show that $V(z)=V(Tz)$ can be extended with higher degree terms to a polynomial $W(z)$ of degrees two through $2d$ which is a sum of $n$ squares. We do this degree by degree. We have already showed that the degree two part of $V(z)$ is a sum of $n$ squares,

[TABLE]

The degrees two and three parts of $V(z)$ are of the form

[TABLE]

Consider the expression

[TABLE]

then

[TABLE]

Let $\delta_{1,i_{2},i_{3}}={\gamma_{1,i_{2}.i_{3}}\over\lambda_{1}}$ then the degrees two and three parts of

[TABLE]

have no terms involving $z_{1}$.

Next consider the expression

[TABLE]

then

[TABLE]

Let $\delta_{2,i_{2},i_{3}}={\gamma_{2,i_{2},i_{3}}\over\lambda_{2}}$ then the degrees two and three parts of

[TABLE]

have no terms involving either $z_{1}$ or $z_{2}$.

We continue on in this fashion defining $\delta_{3}(z),\ldots,\delta_{n}(z)$ such that

[TABLE]

has no terms of degrees either two or three.

We redefine

[TABLE]

This does not change the degree two and three terms of ${\lambda_{1}\over 2}(\delta_{1}(z))^{2}$ and its degree four terms are of the form

[TABLE]

If we let $\delta_{1,i_{2},i_{3},i_{4}}={\gamma_{1,i_{2},i_{3},i_{4}}\over\lambda_{1}}$ then we cancel the degree four terms involving $z_{1}$ in

[TABLE]

Next we redefine

[TABLE]

Again this does not change the degree two and three terms of ${\lambda_{2}\over 2}(\delta_{2}(z))^{2}$ and its degree four terms are of the form

[TABLE]

If we let $\delta_{2,i_{2},i_{3},i_{4}}={\gamma_{2,i_{2},i_{3},i_{4}}\over\lambda_{2}}$ then we cancel the degree four terms involving $z_{2}$ in

[TABLE]

We continue on in this fashion. The result is a sum of squares whose degree two through four terms equal $V(z)$.

Eventually we define

[TABLE]

and

[TABLE]

At degree $d=3$ we solved a linear equation from the quadratic coefficients of $\delta_{1}(z),\ldots,\delta_{n}(z)$ to the cubic coefficients of $V(z)$. We restricted the domain of this mapping by requiring that the quadratic part of $\delta_{i}(z)$ does not depend on $z_{1},\ldots,z_{i-1}$. This made the restricted mapping square, the dimensions of the domain and the range of the linear mapping are the same. We showed that the restricted mapping has an unique solution. If we drop this restriction then the overall dimension of the domain is

[TABLE]

while the dimension of the range is

[TABLE]

So the unrestricted mapping has more unkowns than equations and hence there are multiple solutions.

But the restricted solution that we constructed is a least squares solution to the unrestricted equations because $\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{n}>0$. To see this consider the coefficient $\gamma_{1,1,2}$ of $z_{1}^{2}z_{2}$. If we allow $\gamma_{2}(z)$ to have a term of the form $\delta_{2,1,1}z_{1}^{2}$ we can also cancel $\gamma_{1,1,2}$ by choosing $\delta_{1,1,2}$ and $\delta_{2,1,1}$ so that

[TABLE]

Because $\lambda_{1}\geq\lambda_{2}$ a least squares solution to this equation is $\delta_{1,1,2}={\gamma_{1,1,2}\over\lambda_{1}}$ and $\delta_{2,1,1}=0$. Because $T$ is orthogonal $W(x)=W(T^{\prime}z)$ is also a least squares solution.

6 Example

Consider a double pendulum that we wish to stabilize to straight up. The first two states are the the angles between the two links and straight up measured in radians counterclockwise. The other two states are their angular velocities. The controls are the torques applied at the base of the lower link and at the joint between the links. The links are assumed to massless, the base link is one meter long and the other link is two meters long. There is a mass of two kilograms at the joint between the links and a mass of one kilogram at the tip of the upper link. There is linear damping at the base and the joint with both coefficients equal to $0.5$ sec*-1*. The resulting continuous time dynamics is discretized using Euler’s method with time step $0.1$ sec.

We simulated AHMPC with two different terminal costs and terminal feedbacks. In both cases the Lagrangian was

[TABLE]

The first pair $V^{2}_{f}(x),\kappa^{1}_{f}(x)$ was found by solving the infinite horizon LQR problem obtained by taking the linear part of the dynamics around the operating point $x=0$ and the quadratic Lagrangian (26). Then $V^{2}_{f}(x)$ is a positive definite quadratic and $\kappa^{1}_{f}(x)$ is linear.

The second pair $V^{6}_{f}(x),\kappa^{5}_{f}(x)$ was found using the discrete time version of Al’brekht’s method. Then $V^{6}_{f}(x)$ is the Taylor polynomial of the optimal cost to degree $6$ and $\kappa^{5}_{f}(x)$ is the Taylor polynomial of the optimal feedback to degree $5$. But $V^{6}_{f}(x)$ is not positive definite so we completed the squares to get $V^{10}_{f}(x)$ which is positive definite.

In all the simulations we imposed the control constraint $|u|_{\infty}\leq 5$ and started at $x(0)=(0.9\pi,0.9\pi,0,0)$ with an initial horizon of $N=50$ time steps. The extended horizon was kept constant at $M=5$. The class K function was taken to be $\alpha(s)=s^{2}/10$.

If the Lyapunov and feasibility conditions were violated the horizon $N$ was increased by $5$ and the finite horizon nonlinear program was solved again without advancing the system. If after three tries the Lyapunov and feasibility conditions were still not satisfied then the first value of the control sequence was used, the simulation was adavnced one time step and the horizon was increased by $5$.

If the Lyapunov and feasibility conditions were satisfied over the extended horizon then the simulation was advanced one time step and the horizon $N$ was decreased by $1$.

The simulations were first run with no noise and the results are shown in the following figures. Both methods stabilized the links to straight up in about $t=80$ times steps ($8$ seconds). The degree $2d=10$ terminal cost and the degree $d=5$ terminal feedback seems to do it a little more smoothly and with shorter maximum horizon $N=65$ versus $N=75$ for LQR ($d=1$).

Then we added noise to the simulations. At each advancing step a Gaussian random vector with mean zero and covariance $0.0004$ times the identity was added to the state. The next figures show the results using the degree $10$ terminal cost and degree $5$ terminal feedback. Notice that the horizon converges to zero after $64$ time steps and stays at zero for awhile where the terminal feedback is used. But then the noise causes the horizon to jump to $15$ before it settles back to zero. The LQR terminal cost and feedback failed to stabilize the pendula.

We also considered $d=3$ so that after completing the squares the terminal cost is degree $2d=6$ and the terminal feedback is degree $d=3$. It stabilized the noiseless simulation with a maximum horizon of $N=80$ which is greater than the maximum horizons for both $d=1$ and $d=5$. But it did not stabilize the noisy simulation. Perhaps the reason is revealed by Taylor polynomial approximations to $\sin x$ as shown in Figure 10. The linear approximation in green overestimates the magnitude of $\sin x$ so the linear feedback is stronger than it needs to be to overcome gravity. The cubic approximation in blue underestimates the magnitude of $\sin x$ so the cubic feedback is weaker than it needs to be to overcome gravity. The quintic approximation in orange overestimates the magnitude of $\sin x$ so the quintic feedback is also stronger than it needs to be but by a lesser margin than the linear feedback to overcome gravity. This may explain why the degree $5$ feedback stabilizes the noise free pendula in a smoother fashion than the linear feedback, compare Figures 1 and 2 and Figures 3 and 4.

7 Conclusion

Adaptive Horizon Model Predictive Control is a scheme for varying the horizon length in Model Predictive Control as the stabilization process evolves. It adapts the horizon in real time by testing Lyapunov and feasibility conditions on extensions of optimal trajectories returned by the nonlinear program solver. In this way it seeks the shortest horizons consistent with stabilization.

AHMPC requires a terminal cost and terminal feedback that stabilizes the plant in some neighborhood of the operating point but that neighborhood need not be known explictly. Higher degree Taylor polynomial approximations to the optimal cost and the optimal feedback of the coresponding infinite horizon optimal control problems can be found by an extension of Al’brekht’s method [1]. The higher degree Taylor polynomial approximations to optimal cost need not be positive definite but they can be extended to nonnegative definite polynomials by completing the squares. These nonnegative definite extensions and the Taylor polynomial approximations to the optimal feedback can be used as terminal costs and terminal feedbacks in AHMPC. We have shown by an example that a higher degree terminal cost and feedback can outperform using LQR to define a degree two terminal cost and a degree one terminal feedback.