Explicit Characterisation of Receding Horizon Control Mar a M. - - PowerPoint PPT Presentation

Explicit Characterisation of Receding Horizon Control

Mar´ ıa M. Seron September 2004

Centre for Complex Dynamic Systems and Control

Outline

1

Solution Using Geometry Solution for Arbitrary N Example

2

Implementation of the Explicit Solution Calculation of the Explicit Solution On-line Implementation of the Explicit Solution Suboptimal Solutions and Approximations

Centre for Complex Dynamic Systems and Control

The Receding Horizon Optimal Control Problem

Let the system be given by xk+1 = Axk + Buk,

|uk| ≤ ∆,

(1) where xk ∈ Rn and ∆ > 0 is the input constraint level. Consider the following fixed horizon optimal control problem

P2(x) :

V 

2 (x) min V2({xk}, {uk}),

(2) subject to: xk+1 = Axk + Buk for k = 0, 1, x0 = x, uk ∈ U [−∆, ∆] for k = 0, 1, where the objective function in (2) is V2({xk}, {uk}) 1 2x

2Px2 + 1

2

1

• k=0
• x

kQxk + u kRuk

• .

(3)

Centre for Complex Dynamic Systems and Control

The Receding Horizon Optimal Control Problem

In the objective function we select Q > 0, R > 0 and P satisfying the algebraic Riccati equation P = A PA + Q − K  ¯ RK, (4) where K ¯ R−1BPA,

¯

R R + BPB. (5) Let the control sequence that minimises (3) be

{u

0 , u 1 }.

(6) Then the RHC law is given by the first element of (6) (which depends on the current state x0 = x), that is,

K2(x) = u

0 .

(7)

Centre for Complex Dynamic Systems and Control

Solution of the Associated QP

In the previous lecture we showed that the solution in active region R8 is u(x) =

−Gx − h ∆

• where

G = K + KBKA 1 + (KB)2 h = KB 1 + (KB)2 ∆, if R8 :

       [0 1]u

 (x) ≥ ∆

H1∆[−1 1] ≤ [1 0]Hu

 (x) ≤ H1∆[1 1]

where u

 (x) = −H−1Fx.

Centre for Complex Dynamic Systems and Control

Solution of the Associated QP

The above solution is valid whenever u

 (x) = −H−1Fx belongs to

R8, that is, it satisfies the equations R8 :

       [0 1]u ≥ ∆

H1∆[−1 1] ≤ [1 0]Hu ≤ H1∆[1 1]. Thus, setting u = u

 (x) = −H−1Fx in the above equations we

• btain a characterisation of R8 in the state space:

R8 :

       −[0 1]H−1Fx ≥ ∆

H1∆[−1 1] ≤ −[1 0]Fx ≤ H1∆[1 1].

Centre for Complex Dynamic Systems and Control

Geometric RHC Characterisation for N = 2

If we proceed in a similar way with the remaining regions, we

• btain a characterisation of the RHC law in the form

K2(x) =                          −Kx if x ∈ R, ∆ if x ∈ R1 ∪ R2 ∪ R3, −Gx + h if x ∈ R4, −∆ if x ∈ R5 ∪ R6 ∪ R7, −Gx − h if x ∈ R8,

−8 −6 −4 −2 2 4 6 8 −4 −3 −2 −1 1 2 3 4

x1 k x2 k

R

R1 R2 R3 R4 R5 R6 R7 R8

where the regions R, R1, . . . , R8 in the state space are obtained using geometric arguments as described before. We can see that the above characterisation coincides with the one

• btained previously using Dynamic Programming.

Centre for Complex Dynamic Systems and Control

Geometric RHC Characterisation. More General Cases

For arbitrary horizon N, m ≥ 1 inputs, and more general constraints, as we showed in the previous lecture, the QP solution in a particular active region corresponding to the face

Φℓu = ∆ℓ − Λℓx,

has the form u(x) = H−1/2 ˜

Φ

ℓ[ ˜

Φℓ ˜ Φ

ℓ]−1(∆ℓ − Λℓx)

− H−1/2[I − ˜ Φ

ℓ[ ˜

Φℓ ˜ Φ

ℓ]−1 ˜

Φℓ]H−1/2Fx, Gℓx + hℓ.

Centre for Complex Dynamic Systems and Control

Geometric RHC Characterisation. More General Cases

The above QP solution is valid in the state space region Lℓx ≤ Wℓ, where Lℓ =

• [ ˜

Φℓ ˜ Φ

ℓ]−1(Λℓ − ˜

ΦℓH−1/2F) ˜ Φs ˜ Φ

ℓ[ ˜

Φℓ ˜ Φ

ℓ]−1(Λℓ − ˜

ΦℓH−1/2F) + (Λs − ˜ ΦsH−1/2F)

• Wℓ =
• [ ˜

Φℓ ˜ Φ

ℓ]−1∆ℓ

∆s + ˜ Φs ˜ Φ

ℓ[ ˜

Φℓ ˜ Φ

ℓ]−1∆ℓ

• Centre for Complex Dynamic

Systems and Control

Geometric RHC Characterisation. More General Cases

Thus, the optimal QP solution in an active region (corresponding to active constraints with indices in ℓ) is u(x) Gℓx + hℓ, whenever x satisfies Lℓx ≤ Wℓ. The RHC law KN(x) in the above region is then given by the first m elements of the above solution:

KN(x) = [Im

0](Gℓx + hℓ) if Lℓx ≤ Wℓ. (8)

Centre for Complex Dynamic Systems and Control

Geometric RHC Characterisation. More General Cases

To compute the complete RHC solution in the state space one requires, for example, procedures that enumerate all combinations of active constraints, or recursively partition the state space searching for active regions, or enumerate all parameterised vertices based on double description properties of polyhedra. We will briefly comment on these computational issues later.

Centre for Complex Dynamic Systems and Control

Example

Consider a system with matrices A =

• 0.8955

−0.1897

0.0948 0.9903

• ,

B =

• 0.0948

0.0048

• .

In the objective function we take N = 4, Q =

• 2
• and R = 0.01.

The terminal state weighting matrix P is chosen as the solution of the algebraic Riccati equation P = A PA + Q − K  ¯ RK, where K ¯ R−1BPA and ¯ R R + BPB. We consider input constraints of the form |uk| ≤ 2.

Centre for Complex Dynamic Systems and Control

Example

The state space partition for this case is shown in Figure (a). A “zoom” of this partition is shown in Figure (b). The region denoted by X0 is the projection onto the state space of the constraint polyhedron; in regions X2, X3 and X4 only one constraint is active; in regions X5 and X6 two constraints are active; in region X7 three constraints are active; finally, X1 is the union of all regions where the control is saturated to the value −2.

−2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

x1 x2 X0 X1 X

6

X7

(a) Complete partition.

−1.6 −1.5 −1.4 −1.3 −1.2 −1.1 −1 −0.9 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55

x1 x2 X0 X1 X2 X3 X4 X5 X

6

X7

(b) Zoom.

Centre for Complex Dynamic Systems and Control

Example

The resulting RHC law (8) is

K4(x) = Gix + hi,

if x ∈ Xi, i = 0, . . . , 7, (9) where G0 = −[4.4650 13.5974], h0 = 0, G1 =

[ ],

h1 = −2, G2 = −[5.6901 15.9529], h2 = −0.7894, G3 = −[4.9226 13.8202], h3 = −0.4811, G4 = −[4.5946 13.3346], h4 = −0.2684, G5 = −[6.6778 16.8644], h5 = −1.7057, G6 = −[5.1778 13.4855], h6 = −0.9355, G7 = −[7.4034 16.8111], h7 = −2.6783. Similar expressions hold in the remaining unlabelled regions. These can be obtained by symmetry.

Centre for Complex Dynamic Systems and Control

Example

The figures show the state state partitions for horizons N = 2, N = 3, N = 4 and N = 5.

−3 −2 −1 1 2 3 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

x1 x2 N=2

−3 −2 −1 1 2 3 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

x1 x2 N=3

−3 −2 −1 1 2 3 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

x1 x2 N=4

−3 −2 −1 1 2 3 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

x1 x2 N=5

Centre for Complex Dynamic Systems and Control

Example

We next consider an initial condition x0 = [−1.2 0.53] and simulate the system under the RHC (9). The figure shows the resulting state space trajectory. The trajectory starts in region X4 and moves, successively, into regions X6, X5, X1, X1, X0, and stays in X0 thereafter.

−1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0.1 0.2 0.3 0.4 0.5 0.6

x1 x2

Centre for Complex Dynamic Systems and Control

Recapitulation

The geometric characterisation of QP of the form

˜

u(x) = ˜

Φ

ℓ[ ˜

Φℓ ˜ Φ

ℓ]−1(∆ − Λℓx) − [I − ˜

Φ

ℓ[ ˜

Φℓ ˜ Φ

ℓ]−1 ˜

Φℓ]H−1/2Fx (10) = H1/2(Gℓx + hℓ)

if

• [ ˜

Φℓ ˜ Φ

ℓ]−1(Λℓ − ˜

ΦℓH−1/2F) ˜ Φs ˜ Φ

ℓ[ ˜

Φℓ ˜ Φ

ℓ]−1(Λℓ − ˜

ΦℓH−1/2F) + (Λs − ˜ ΦsH−1/2F)

• Lℓ

x ≤

• [ ˜

Φℓ ˜ Φ

ℓ]−1∆ℓ

∆s + ˜ Φs ˜ Φ

ℓ[ ˜

Φℓ ˜ Φ

ℓ]−1∆ℓ

• Wℓ

,

(11) leads to the RHC characterisation

KN(x) = [Im

0](Gℓx + hℓ) if Lℓx ≤ Wℓ. The same characterisation can be obtained using the KKT

• ptimality conditions, as we showed in the previous lecture.

Centre for Complex Dynamic Systems and Control

Explicit vs Implicit

Recall what we defined as explicit procedures as opposed to implicit, or numerical, procedures: Explicit

x

fp(z) = z2 + 2apz + bp

∂fp ∂z = 0 ⇒ z(p) = −ap

Implicit (numerical) fp(z) = z2 + 2apz + bp

z

z0

p

zk

p

zk+1

p

. . . We have now obtained an explicit characterisation of RHC, which can replace the traditional numerical RHC solution in online implementations.

Centre for Complex Dynamic Systems and Control

Numerical Implementation of RHC

Centre for Complex Dynamic Systems and Control

Explicit Solution

Centre for Complex Dynamic Systems and Control

Implementation of the Explicit Solution

Key issues: Off line How to efficiently calculate KN(x) for all x of interest (region partition and corresponding control law). How to store KN(x) as a function of x (look-up table). On line How to retrieve the value of KN(x) given x.

Centre for Complex Dynamic Systems and Control

Calculation of the Explicit Solution

Naive approach: enumeration of all possible combinations of active constraints. Semi-naive approach: combine the above with dynamic programming. More efficient: use algorithms to implement multiparametric quadratic programming (Dua, Pistikopoulos, Bemporad, T¨

• ndel, Johansen).

Better: combine the above with reachability analysis (Grieder, Borrelli, Torrisi, Morari). New approach: enumeration of all parameterised vertices based on double description properties of polyhedral regions (Olaru and Dumur).

Centre for Complex Dynamic Systems and Control

On-line Implementation of the Explicit Solution

Once the look-up table is precomputed and stored, the on-line computation is reduced to: (i) identifying the “current region”, that is, deciding in which region the current state belongs and (ii) computing the control input using the affine control law corresponding to that region. The first step is the most computationally demanding.

Centre for Complex Dynamic Systems and Control

Searching the Look-up Table

Naive approach: sequential search through the regions of the polyhedral partition. Efficient approach: organise the partition in a binary search tree: each level, or “node”, is associated with one hyperplane

• inequality. When the tree is used “on line” for the current state,

at each node one inequality is evaluated, its sign checked, and the left or right subtree selected based on the sign. Traversing the tree from the “root” to a “leaf” node, the control law corresponding to the current state is found. Tøndel, Johansen and Bemporad (2002) proposed an algorithm to construct a binary search tree with search time logarithmic in the number of polyhedral regions.

Centre for Complex Dynamic Systems and Control

Binary Search Algorithm

Centre for Complex Dynamic Systems and Control

Example: A Simple Binary Search

Consider the double integrator with matrices A =

• 1

1 1

• ,

B =

• 0.5

1

• ,

C =

• 1
• ,

|u(k)| ≤ 1, Q = I2, R = 0.01, and P chosen as the solution of the

algebraic Riccati equation. For N = 5, for example, the state-space partition and RHC law are:

−10 −8 −6 −4 −2 2 4 6 8 10 −5 −4 −3 −2 −1 1 2 3 4 5

x1 x2 X0 X1 X2 X3 X4 X5

K5(x) = Gℓx + hℓ

if x ∈ Xℓ

Xℓ Gℓ hℓ X0

• [0.9653

1.3895] X1

• [0.6154

1.2870]

• 0.4156

X2

• [0.4390

1.2121]

• 0.7982

X3

• [0.3399

1.1665]

• 1.1746

X4

• [0.2771

1.1367]

• 1.5495

X5

• 1

Centre for Complex Dynamic Systems and Control

Example: A Simple Binary Search

The “look-up table” implementation of RHC checks, at each time instant, in which region the current state is and picks the corresponding control law from a stored table. We first consider a simple Matlab code that implements this search in one step. This has only a few lines of code. We will use the double integrator example with N from 2 to 10, in the area −20 ≤ x1 ≤ 20, −8 ≤ x2 ≤ 8, to compare the performance

• f the ‘look-up table” implementation of RHC and Matlab’s QP

algorithm. The total data stored for the look-up table uses 20448 bytes.

Centre for Complex Dynamic Systems and Control

Example: A Simple Binary Search

We grid the state space with 3600 points and, for each point of the grid we evaluate the controller using each method, and for each method we measure computation time we measure number of floating point operations. We then average the latter two measures over the grid. We repeat the same for horizons N = 2, . . . , 10.

Centre for Complex Dynamic Systems and Control

Example: A Simple Binary Search

The left figure is a plot of computation times, and right figure shows the number of floating point operations (flops), for the “look-up table” method (solid line) and QP (dashed-dotted line) as a function of the horizon N.

2 3 4 5 6 7 8 9 10 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 N seconds

Figure: Computation time

2 3 4 5 6 7 8 9 10 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10

5

N flops

Figure: Flops

Centre for Complex Dynamic Systems and Control

Example: A Simple Binary Search

The performance of both algorithms is comparable in relation with computation speed, although the number of flops is much larger for the “look-up table” method. This is expected since, as the horizon increases, there are more regions in the state space partition that have to be checked at each instant. This problem can be easily mitigated by splitting the regions into an orthogonal binary search tree; this increases the total number of regions (as some of the original regions are divided into smaller subregions) but the number of checks at each sample time is smaller and so there are less flops.

Centre for Complex Dynamic Systems and Control

Example: A Simple Binary Search

To accomplish this task we divide the area of interest into the smaller subregions of a grid defined by Nh separating hyperplanes parallel to the axes. The original RHC partition is then intersected with the grid to form a new partition. The maximum number of inequalities to check at each sample time is then Nh plus the maximum number of inequalities of the new RHC regions inside each subregion of the grid.

Centre for Complex Dynamic Systems and Control

Example: A Simple Binary Search

For the double integrator example, the total number of regions

• f the original RHC partition in the area of interest

−20 ≤ x1 ≤ 20, −8 ≤ x2 ≤ 8 is N0 = 165 and the total number

• f hyperplane inequalities to check at each time is 632.

We next grid the area using using five horizontal lines and seven vertical lines (that is, Nh = 12), defining a grid with 48 subregions, as shown in the following figure. The intersection of the original RHC partition with the grid yields a total of 672 regions, but the maximum number of inequalities to check at each sample time is Nh+96=108. The data stored in this case use 164352 bytes (compare with 20448 bytes used before).

Centre for Complex Dynamic Systems and Control

Example: A Simple Binary Search

−20 −15 −10 −5 5 10 15 20 −8 −6 −4 −2 2 4 6 8 x1 x2

Figure: Binary partition using 12 hyperplanes

Centre for Complex Dynamic Systems and Control

Example: A Simple Binary Search

The inclusion of the binary search tree in the algorithm for table look-up does not add too much complexity to the Matlab code, which has approximately 30 lines. To compare the “on-line” performance of the new table look-up algorithm and QP , we proceed as before, averaging computation time and number of floating point operations over a state space grid of 3600 points. In all cases, the binary search tree has Nh = 12 separating hyperplanes.

Centre for Complex Dynamic Systems and Control

Example: A Simple Binary Search

The left figure is a plot of computation times, and the right figure shows the number of floating point operations (flops), for the “binary search look-up table” method (solid line) and QP (dashed-dotted line) as a function of the horizon N.

2 3 4 5 6 7 8 9 10 0.005 0.01 0.015 0.02 0.025 0.03 N seconds

Figure: Computation time

2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 x 10

4

N flops

Figure: Flops

Centre for Complex Dynamic Systems and Control

Example: A Simple Binary Search

We can see from these figures that the table look-up using binary search has much better performance than the standard QP algorithm. Both measures (computation time and number of floating point operations) stay approximately constant for all horizons.

Centre for Complex Dynamic Systems and Control

Suboptimal Solutions and Approximations

For RHC problems of large dimension (large state and input dimensions and long constraint horizons) it may be unrealistic to compute the exact explicit solution, and even if that were possible, its implementation would be impractical due to the large amount of memory that would be required to store a complex region partition. On the other hand, on-line optimisation to solve the associated QP may also be impractical from the computation time perspective.

Centre for Complex Dynamic Systems and Control

Suboptimal Solutions and Approximations

Possible simplifications: Regional solutions. Approximate partition for example, using orthogonal hyperplanes in the tree nodes. Approximate partition using neural networks. Input trajectory parameterisation to reduce degrees of freedom. Allow certain degree of suboptimality and constraint violation. More ... on-going research and growing literature.

Centre for Complex Dynamic Systems and Control