Structure of the Hessian Graham C. Goodwin September 2004 Centre - PowerPoint PPT Presentation

Structure of the Hessian Graham C. Goodwin September 2004 Centre for Complex Dynamic Systems and Control

1. Introduction We saw in earlier lectures that a core ingredient in quadratic constrained optimisation problems is the Hessian matrix H . So far we have simply given an “in principle” approach to the evaluation of this matrix. That is, for a system with state equations x k + 1 = Ax k + Bu k , we can compute the Hessian as H = Γ  Q Γ + R , where Γ , Q , R . Centre for Complex Dynamic Systems and Control

This will be satisfactory for simple problems. However, for more complex problems (for example, high order systems or problems having mixed stable and unstable modes) this “brute force” approach may fail. A hint as to the source of the difficulties is that the direct way of computing the Hessian depends on powers of the system matrix A . Clearly, if the system has unstable modes, then some entries of Γ will diverge as N increases. Centre for Complex Dynamic Systems and Control

We will show in this lecture that this problem can be resolved by focusing attention on the stable and unstable parts of the system separately. Centre for Complex Dynamic Systems and Control

2. Separating Stable and Unstable Modes Our eventual goal in this lecture is to gain a better understanding of the structure of the Hessian matrix particularly for large optimisation horizons. However, as mentioned above, the straightforward approach to evaluating the Hessian will often meet difficulties for open loop unstable plants due to exponential divergence of the system impulse response. Centre for Complex Dynamic Systems and Control

One way of addressing this problem is to recognise that there is an intimate connection between “stability” and “causality.” In particular, a system having all modes unstable becomes stable if viewed in reverse time, that is, as an anti-causal system. This line of reasoning leads to an alternative viewpoint in which unstable modes are treated differently. We show below that this leads to an equivalent problem formulation with a different Hessian having different properties. Centre for Complex Dynamic Systems and Control

Consider a discrete time linear system and suppose that it has no eigenvalues on the unit circle. We can then partition the state vector x k ∈ R n as � � x s k x k = , x u k where the states x s k and x u k are associated with the stable and unstable modes, respectively. Centre for Complex Dynamic Systems and Control

We can then factor the state equations into stable and unstable parts as follows: � � � � � � � � x s x s A s 0 B s k + 1 k = + u k , (1) x u x u 0 A u B u k + 1 k � � � x s � k y k = C x k = C s C u , x u k where u k ∈ R m and y k ∈ R p ( p ≥ m ) . The eigenvalues of A s have moduli less than one, and the eigenvalues of A u have moduli greater than one. Centre for Complex Dynamic Systems and Control

We can then express the solution of the system equations as k − 1 � A k − 1 − j x s k = A k s x s 0 + B s u j for k = 1 , . . . , N , (2) s j = 0 N − 1 � k = A − ( N − k ) A k − 1 − j x u µ − B u u j for k = 0 , . . . , N − 1 . (3) u u j = k x u N � µ, Centre for Complex Dynamic Systems and Control

We then need an equality constraint that both µ and the sequence of control signals { u 0 , . . . , u N − 1 } need to satisfy in order to bring the unstable states back to their correct initial values, that is, N − 1 � A − j − 1 A − N B u u j = x u u µ − 0 . (4) u j = 0 Centre for Complex Dynamic Systems and Control

We are thus led to the following equivalent statement of the optimisation problem. P N ( x ) : V  N ( x ) � min V N ( { x k } , { u k } , µ ) , (5) subject to: � � x s k x k = for k = 0 , . . . , N , x u k k − 1 � A k − 1 − j x s k = A k s x s 0 + B s u j for k = 1 , . . . , N , s j = 0 N − 1 � x u k = A − ( N − k ) A k − 1 − j µ − B u u j for k = 0 , . . . , N − 1 , u u j = k (6) Centre for Complex Dynamic Systems and Control

� � x s x 0 = 0 = x , x u 0 u k ∈ U for k = 0 , . . . , N − 1 , x k ∈ X for k = 0 , . . . , N − 1 , � � x s N x N = ∈ X f ⊂ X , µ where N − 1 � �  � � x s x s V N ( { x k } , { u k } , µ ) � 1 + 1 � N N ( x  k Qx k + u  P k Ru k ) . (7) 2 µ µ 2 k = 0 Centre for Complex Dynamic Systems and Control

The above formulation of the problem, at the expense of the introduction of the additional optimisation variable µ , avoids exponentially diverging terms in the computation of the Hessian matrix. Thus, at least intuitively, it would seem to be more apposite for studying the structure of the problem for large horizons. Centre for Complex Dynamic Systems and Control

3. Numerical Issues in the Computation of the Hessian We represent the time evolution of the system output using the usual vector notation. Thus, let �  , � y = y  y  . . . y  1 2 N (8) �  . � u = u  u  u  . . . 0 1 N − 1 Centre for Complex Dynamic Systems and Control

y = ( Γ s + Γ u ) u + Ω s x s 0 + Ω u µ, (9) � ����� �� ����� � ¯ Γ where     C s A s C s B s 0 · · · 0             C s A 2    C s A s B s C s B s 0  · · ·      s        Ω s =   , Γ s =   , . . . .    ...      . . . .      .   . . .                  C s A N C s A N − 1 C s A N − 2 B s B s · · · C s B s s s s Centre for Complex Dynamic Systems and Control

 C u A − ( N − 1 )   C u A − ( N − 1 )  C u A − 1 C u A − 2 0 u B u u B u . . . B u     u u          C u A − ( N − 2 )   C u A − ( N − 2 )  C u A − 1     0 0 u B u B u . . .         u u          .   . . .  ... ...     Ω u = . , Γ u = − . . .         . . . .                 C u A − 1 C u A − 1     0 0 · · · 0 u B u         u             C u 0 0 0 0 · · · Centre for Complex Dynamic Systems and Control

The matrix ¯ Γ � Γ s + Γ u has the form  ¯ ¯ ¯ ¯  h 0 h − 1 h − 2 h − ( N − 1 ) . . . . . .       ¯ ¯ ¯ ¯   h 1 h 0 h − 1 . . . . . . h − ( N − 2 )          . . .  ... ...   . . .    . . .        ¯ Γ � Γ s + Γ u =   , (10)  . . .  ... ...   . . .    . . .          . .    . .  ¯ ¯   h 0 h − 1 . .          ¯ ¯ ¯ ¯  h N − 1 h N − 2 . . . . . . h 1 h 0 Centre for Complex Dynamic Systems and Control

where { ¯ h k : k = − ( N − 1 ) , . . . , N − 1 } , is a finite subsequence of the infinite sequence { ¯ h k : k = −∞ , . . . , ∞} defined by { ¯ h k : k = 0 , . . . , ∞} � { C s B s , C s A s B s , C s A 2 s B s , . . . } , (11) { ¯ h k : k = − 1 , . . . , −∞} � {− C u A − 1 u B u , − C u A − 2 u B u , − C u A − 3 u B u , . . . } . (12) Centre for Complex Dynamic Systems and Control

Consider, P = Q and Q = C  C and R = ρ I m > 0 . V N as follows 1 : V N ( x , u , y , µ ) = 1 2 ( x  Qx + y  y + ρ u  u ) . 1 We keep the function V N but change its arguments as appropriate. Centre for Complex Dynamic Systems and Control

V N ( x , u , µ ) = 1 � � ¯ �  � ¯ � � Γ u + Ω s x s Γ u + Ω s x s + ρ u  u x  Qx + 0 + Ω u µ 0 + Ω u µ 2 = 1 Γ  � � Γ  ¯ Γ + ρ I ] u + u  ¯ 2 u  [¯ Ω s x s 0 + Ω u µ + 1 � �  � � Ω s x s Ω s x s 0 + Ω u µ 0 + Ω u µ . (13) 2 Centre for Complex Dynamic Systems and Control

With respect to the new optimisation variables ( u , µ ) , the modified Hessian of the quadratic objective function is � ¯ Γ  ¯ � ¯ Γ + ρ I Γ  Ω u H ′ � . (14) u ¯ Ω  Γ Ω  u Ω u Centre for Complex Dynamic Systems and Control

For future use, we extract the left-upper submatrix, which we call the “regularised sub-Hessian”: Γ  ¯ H N � ¯ ¯ Γ + ρ I . (15) Centre for Complex Dynamic Systems and Control

The above problem formulation has the ability to ameliorate the numerical difficulties encountered when dealing with unstable plants. Indeed, we see that all terms depend only on exponentially decaying quantities. Centre for Complex Dynamic Systems and Control