H-infinity Norm Calculation via a State Space Formulation

George E Kuster

Thesis submitted to the Faculty of the

Virginia Polytechnic Institute and State University
in partial fulfillment of the requirements for the degree of

Master of Science
in
Mathematics

Lizette Zietsman, Chair

Jeffrey T. Borggaard
Megan Wawro

December 7, 2012
Blacksburg, Virginia

Keywords: Optimization, Systems, Control, H-infinity, Feedback

Copyright 2012, George E Kuster
 H-infinity Norm Calculation via a State Space Formulation

George E Kuster

(ABSTRACT)

There is much interest in the design of feedback controllers for linear systems that minimize
the H∞ norm of a specific closed-loop transfer function. The H∞ optimization problem
initiated by Zames (1981), [23], has received a lot of interest since its formulation. In
H∞ control theory one uses the H∞ norm of a stable transfer function as a performance
measure. One typically uses approaches in either the frequency domain or a state space
formulation to tackle this problem. Frequency domain approaches use operator theory, J-
spectral factorization or polynomial methods while in the state space approach one uses ideas
similar to LQ theory and differential games. One of the key computational issues in the design
of H∞ optimal controllers is the determination of the optimal H∞ norm. That is, determining
the infimum of r for which the H∞ norm of the associated transfer function matrix is less
than r. Doyle et al, [7], presented a state space characterization for the sub-optimal H∞
control problem. This characterization requires that the unique stabilizing solutions to two
Algebraic Riccati Equations are positive semi definite as well as satisfying a spectral radius
coupling condition. In this work, we describe an algorithm by Lin et al, [16], used to calculate
the H∞ norm for the state feedback and output feedback control problems. This algorithm
only relies on standard assumptions and divides the problem into three sub-problems. The
first two sub-problems rely on algorithms for the state feedback problem formulated in the
frequency domain as well as a characterization of the optimal value in terms of the singularity
of the upper-half of a matrix created by the stacked basis vectors of the invariant sub-space
of the associated Hamiltonian matrix. This characterization is verified through a bisection or
secant method. The third sub-problem relies on the geometric nature of the spectral radius of
the product of the two solutions to the Algebraic Riccati Equations associated with the first
two sub-problems. Doyle makes an intuitive argument that the spectral radius condition
will fail before the conditions involving the Algebraic Riccati Equations fail. We present
numerical results where we demonstrate that the Algebraic Riccati Equation conditions fail
before the spectral radius condition fails.
Dedication

To my grandparents.

iii
Acknowledgments

Completing this project would have been an impossible task without the support of numer-
ous others along the way. It would be an gross understatement to say that without the
guidance and enthusiasm of my advisor, Dr. Lizette Zietsman this thesis would not have
come to fruition. Her insight, patience and ability to keep me motivated are several reasons
I remained working steadily.

I would like to thank my fellow graduate students who served as sound boards when I needed
to piece together difficult topics. Specifically Alan Lattimer and Steven Boyce played indis-
pensable roles in helping me synthesize the many topics encountered in this research.

Dr. Vasil Skenderi was the most influential undergraduate mathematics professor I had the
pleasure of taking a course with. Unknowingly, he sparked my interest in research that in-
stilled a desire to attend graduate school.

Above all however, I am grateful to my wife Lianne for her unwavering support, encourage-
ment and understanding throughout this entire process.
George Kuster

iv
Contents

1 Introduction 1

2 Background 3
2.1 State Space Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 State Space Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3 Transfer Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.4 Controllability and Observability . . . . . . . . . . . . . . . . . . . . . . . . 8
2.5 Stabilizability and Detectability . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.6 H∞ norm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3 Algorithms: Frequency Domain 14

3.1 Computing the H∞ norm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2 Bisection Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.3 Two Step Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

4 Algorithms: State space Approach 20

4.1 The Hamiltonian and Associated Riccati Equation . . . . . . . . . . . . . . . 20
4.2 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

5 Conclusions 35
5.1 Conclusion on State Space Formulation . . . . . . . . . . . . . . . . . . . . . 35
5.2 Extending to Mathematics Education . . . . . . . . . . . . . . . . . . . . . . 36

Bibliography 38

v
List of Figures

2.1 Block Diagram: Continuous time state space model. . . . . . . . . . . . . . . 3

2.2 Spring mass damper system. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

3.1 Relationship between singular values and eigenvalues . . . . . . . . . . . . . 17

4.1 Variant secant method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

vi
List of Tables

4.1 Numerical example 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4.2 Numerical example 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.3 Numerical example 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

vii
Chapter 1

Introduction

One of the many challenges in control theory is the design of controllers that result in
sufficient performance of a plant under various types of inputs and disturbances. H∞ control
theory was designed to reduce modeling errors and unknown disturbances in a system, while
providing quantifiable optimization of large scale multivariable problems. In application most
solutions to H∞ control problems are actually sub-optimal controllers. That is, controllers
which achieve a predetermined bound on the H∞ norm of an associated transfer function. In
this formulation of the problem, the solution is robust for a predetermined bound. However,
if one can find the H∞ norm, one can answer the question as to the existence of a controller
for a certain range of perturbations to the system.
George Zames, [23] introduced the formal H-infinity control theory by formulating the
probloem of sensitivity reduction as an optimization problem with an operator norm, in
particular, the H-infinity norm [15]. Where H∞ is the space of all bounded analytic matrix
valued functions in the open right half complex plane. This formulation was based entirely
in the frequency domain. Zames suggested that using the H-infinity norm as a performance
measure would better satisfy the demands of applications in comparison to the popular Lin-
ear Quadratic Gaussian control design [25]. This brought about the design of H-infinity
optimal control, where one aims to find a controller which stabilizes a system while minimiz-
ing the effects of disturbances in the system. This norm is now used to numerically evaluate
the controller’s sensitivity, robustness and performance of the closed loop feedback system.
Shortly afterwards John Doyle developed tools for testing robust stability in an effort to
gauge model uncertainties and achieve the H-infinity control objective of minimizing the
effect of these perturbations on the regulated outputs [18]. In 1984 Doyle, using state space
methods, offered the first solution to a general multivariable H-infinity optimal control prob-
lem [6].

While these frequency domain formulations dealt with the robustness issues, the tools devel-
oped were both notedly complicated and provided limited insight concerning the nature of

1
George Kuster Chapter 1. Introduction 2

the stabilizing solutions. State space formulations were initially developed by Doyle, Francis
and Glover in the mid 1980’s, but were only successful to a degree, that is the solutions to
these early formulations had solutions of high order, see [8] and [11].

In the late 1980’s work by Khargonekar, Petersen, Rotea and Zhou, [13], [14], and [20], pre-
sented the solution to the H-infinity state feedback problem in terms of a constant feedback
gain matrix, and provided a formula for finding this gain matrix in terms of the solution to
a Riccati equation. Furthermore, they established a relationship between linear quadratic
game theory and H-infinity optimal control and gave the solution to the state-feedback
problem by solving an algebraic Riccati equation.
In 1989 Doyle, Glover, Khargonekar and Francis, [7], presented the first general state space
solution to the H-infinity control problem in their now classic paper State space solutions
to standard H2 and H-infinity control problems. Results in this paper, give necessary and
sufficient conditions for the existence of an admissible controller in terms of solutions of
algebraic Riccati equations and a coupling condition involving a spectral radius condition.
Thus formulation resulted in algorithms that uses concave criteria as well as Newton-like
methods, see [9], [10], and [21].
This study focuses on a collection of algorithms used for H∞ -control purposes, in particular
for state feedback, [3], [4], and output feedback problems. The latter though the algorithm
presented by Lin, Wang and Xu [16].
The outline of the remainder is as follows. Chapter 2 provides the necessary background
such as controllability, observability, detectabilty, stabilizability as well as other topics such
as transfer functions and introduces the H∞ norm. Chapter 3 discusses two algorithms based
on frequency domain tools to calculate the H∞ norm of a transfer function associated with
the state feedback case. The first algorithm is the bisection algorithm presented by [3] and
the second algorithm is a two step algorithm originally presented by [4].
The main discussion is found in Chapter 4 where we discuss and implement the algorithm
Lin, Wang and Xu [16] proposed to calculate the H∞ norm of the state and output feedback
control problems. The algorithm uses a state space formulation while appealing to frequency
domain tools to find the infimum of r for which the H∞ norm of the associated transfer
function matrix is less than r. The problem of finding the H∞ norm is broken into three
sub-problems. The first two rely heavily on the algorithms discussed in Chapter 3 and the
third relies on the geometric nature of the spectral radius of the product of the two solutions
to the Algebraic Riccati Equations associated with the first two sub-problems. Concluding
remarks are found in Chapter 5.
Chapter 2

Background

2.1 State Space Formulation

For the initial discussion consider the basic time invariant, continuous time system
ẋ(t) = Ax(t) + Bu(t), x(t0 ) = x0
(2.1.1)
y(t) = Cx(t) + Du(t).

Here A is an n × n matrix, B is n × m, C is r × n and D is r × m. The vector x represents

the state, u the control and y is the measured output. The basic block diagram associated
with this system is given in Figure 1.1.

Figure 2.1: Block Diagram: Continuous time state space model.

Example 2.1.1. To better illustrate the state space formulation, consider the simple example
where two masses are connected by a spring and damper. The control force u is applied to a
mass M1 , which directly effects the position, z of second mass M2 . The input of this system
is the force, and the output is the change of position of the second mass M2 . Here we would
like to control the displacement of M2 .

3
George Kuster Chapter 2. Background 4

Using Newton’s Second Law, we attain the following equations:

M1 ẅ = −b(ẇ − ż) − k(w − z) + u (2.1.2)

M2 z̈ = b(ẇ − ż) + k(w − z). (2.1.3)

Damper  

  ss  
   
!  
     M1              M2  

Spring  

!  
 

Figure 2.2: Spring mass damper system.

In this example, b denotes the damping coefficient, w is the displacement of M1 , z is the

displacement of M2 , and k is the spring constant. This second order system are rewritten as
a system of first order equations:
Let x1 = w, x˙1 = ẇ = x2 , x3 = z, x˙3 = ż = x4 . Then,
b k
x˙2 = − (x2 − x4 ) − (x1 − x3 ) + u, (2.1.4)
m1 m1
b k
x˙4 = (x2 − x4 ) + (x1 − x3 ) + u, (2.1.5)
m2 m1
y = x3 . (2.1.6)

In matrix-vector notation (2.1.4)-(2.1.6) become

       
x˙1 0 1 0 0 x1 0
 x˙2   − k − b k b   x2   1 
 =  m1 m1 m1 m1 ·   m1  · u,
  x3  +  0 (2.1.7)

 x˙3   0 0 0 1 
k b k b
x˙4 m2 m2
− m2
− m2
x4 0
 
x1
   x2 
y = 0 0 1 0 ·  x3  .
 (2.1.8)
x4
George Kuster Chapter 2. Background 5

This spring mass damper system is now of the form (2.1.1)

where
   
0 1 0 0 0
 −k − b k b   1   
m1 m1 m1 m1  B =  m1  C = 0 0 1 0
A=
 0 D = 0.
0 0 1   0 
k b
m2 m2
− m2 − mb2
k
0

2.2 State Space Solutions

If the control u is known, the system (2.1.1) can be solved explicitly. Before we discuss this
in more detail, we review a few basic concepts and follow the discussion in [5].
Definition 2.2.1. Given an n × n matrix A, the matrix exponential eAt is defined as:
∞
At
X (At)k
e = .
k=0
k!

We note two properties of the matrix exponential:

1. eA0 = A0 = I.
∞ ∞ ∞
d At X d Ak tk X Ak tk−1 X Ak tk
e = = = A = AeAt
2. dt k=0
dt k! k=1
(k − 1)! k=0
k! .

It is well known that the solution to the initial value problem

ẋ(t) = Ax(t), x(t0 ) = x0 , (2.2.1)

is given by x(t) = eA(t−t0 ) x0 .

Considering the more general formulation for given u,
ẋ(t) =Ax(t) + Bu(t), x(t0 ) = x0 (2.2.2)
y(t) =Cx(t) + Du(t), (2.2.3)
the solution is obtained through variation of parameters or using the Laplace transform. We
first review the definition of the Laplace transform and the inverse Laplace transform.
Definition 2.2.2. The Laplace Transform of a function f (x), x ∈ [0, ∞) is defined by
the integral Z ∞
ˆ
L(f (x))(s) = f (s) = f (x)e−sx dx,
0
for values of s ∈ C for which this integral is defined.
George Kuster Chapter 2. Background 6

The inverse Laplace transform of fˆ, denoted by L−1 (fˆ), is the function f whose Laplace
transform is fˆ.
To determine the solution to (2.2.2)–(2.2.3), take the Laplace transform of (2.2.2) and solve
for x̂. This results in
x̂(s) = (sI − A)−1 x0 + (sI − A)−1 B û(s). (2.2.4)
Taking the inverse Laplace of (2.2.4) gives

x(t) = L−1 (sI − A)−1 x0 + L−1 (sI − A)−1 Bu(t)

Z t
A(t−t0 )
=e x0 + eA(t−s) Bu(s)ds. (2.2.5)
t0

Rt
Substituting (2.2.5) into (2.2.3) yields y(t) = CeA(t−t0 ) x0 + t0
CeA(t−s) Bu(s)ds + Du(t).

Definition 2.2.3. The matrix eA(t−t0 ) is called the state transition matrix.

In the event that the control u is not predetermined, various control strategies are used to
determine the control that would result in a desired outcome. In this thesis our focus will
be on H∞ control which is described in Section ???

2.3 Transfer Function

In this section we introduce the transfer function associated with a linear time-invariant
system such as (2.1.1) with x(0) = x0 . The transfer function gives a relationship between
the input and output of the system and plays an important role in control theory in general
and provides insight in how disturbances in the system will affect the output.
To derive the transfer function associated with (2.1.1)
ẋ(t) = Ax(t) + Bu(t), x(0) = x0
(2.3.1)
y(t) = Cx(t) + Du(t)

take the Laplace transform as before in (2.2.4) to get

x̂(s) = (sI − A)−1 x0 + (sI − A)−1 B û(s) (2.3.2)

ŷ(s) = C x̂(s) + Dû(s). (2.3.3)

ŷ(s)
Substituting (2.3.2) into (2.3.3), and setting G(s) = Equation (2.3.3) becomes
û(s)
ŷ(s) = G(s)û(s)

or
ŷ(s)
G(s) = . (2.3.4)
û(s)
George Kuster Chapter 2. Background 7

The function G(s) reflects the ratio between the system output ŷ(s) and the system input
û(s) and is called the matrix transfer function or simply the transfer matrix. Entry
(i, j) in the matrix represents the transfer function from the j th input to the ith output.
This transfer function represents properties of the system and is not representative of the
magnitude or nature of the inputs. It also does not serve to provide any information about
the structure of the system.
One main advantage to using transfer functions is they allow for taking a complicated system
in a time domain and represent it as a function of a single variable in frequency domain.
Unfortunately, this approach only works for single input single output systems. This means,
for multiple input-multiple output systems, there are multiple transfer functions used to
represent each input-output relationship. It should be noted that transfer functions are not
unique, that is many systems can share the same transfer function. They are unrelated to
the systems magnitude or the nature of the input.
Definition 2.3.1. The points p at which the transfer function G(p) = ∞ are called the poles
of the system.

If G(∞) is a constant matrix the transfer function is called proper and if G(∞) = 0, the
transfer function is called strictly proper.
Let us use the spring mass damper system from Example 2.1.1 to illustrate how one would
find the transfer function of a given system described by a state space representation.
Example 2.3.1.

From Example 2.1.1 we have

0 1 0 0
   
0
 −k − b k b   1 
 m1 m1 m1 m1
 , B =  m1
 
A=  , C = 0 0 1 0 , D = 0.
  
 0 0 0 1   0 
k b
m2 m2
− m2 − mb2
k
0
Thus
   −1  
s 0 0 0 0 1 0 0 0
0 s 0 0   −4 −2 4 2  1 
G(s) = C(sI−A)−1 B+D =
    
0 0 1 0   −    .
 0 0 s 0   0 0 0 1   0 
0 0 0 s 4 2 −4 −2 0

Simplification leads to
s+2
G(s) = .
s4 + 3s3 + 6s2

Definition 2.3.2. For the continuous time state space representation given by (2.3.1), the
matrix
G(iω) = C(iωI − A)−1 B + D (2.3.5)
is called the frequency response matrix, where ω ∈ R is the frequency.
 George Kuster Chapter 2. Background 8

2.4 Controllability and Observability

Controllability and observability are fundamental ideas of control theory. We follow the
discussion in [5] to introduce these concepts.

Definition 2.4.1. The system described in (2.1.1),

ẋ(t) = Ax(t) + Bu(t),

y(t) = Cx(t) + Du(t)

is called controllable if the system can be driven to any final state x1 = x(t1 ), in finite time
t1 , from any initial state x(0) by selecting the control u(t), 0 ≤ t ≤ t1 accordingly.

The controllability of the system (2.1.1), is referred to as the controllability of the pair
(A, B). The following theorem provides verifiable conditions as to whether or not a system
is controllable.

Theorem 2.4.1. Let A ∈ Rn×n and B ∈ Rn×m (m ≤ n) then the following statements are
equivalent:

(i) The system (2.1.1) is controllable.

(ii) The n × nm controllability matrix CM = [B, AB, A2 B, . . . , An−1 B] has full rank.

(iii) The matrix Z t1

T
WC = eAt BB T eA t dt
0
is nonsingular for any t1 > 0.

(iv) If (λ, x) is an eigenpair of AT , then xT B 6= 0.

(v) Rank(A − λI, B) = n for every eigenvalue λ of A.

(vi) The eigenvalues for A − BK can be arbitrarily assigned by an appropriate choice of K.

Observability and controllability of a system are dual problems and observability is defined
as follows:

Definition 2.4.2. The continuous time system (2.1.1) is said to be observable if there
exists a t1 > 0 such that the initial state x(t0 ) can be uniquely determined from knowledge of
u(t) and y(t) for all t, where 0 ≤ t ≤ t1 .

The observability of the system (2.1.1) is referred to as the observability of the pair (A, C).
The duality means that (A, C) is observable if (A,T C T ) is controllable. Due to the duality of
observability and controllability, observability has similar characterizations as controllability
in Theorem 2.4.1.
 George Kuster Chapter 2. Background 9

Theorem 2.4.2. The following statements are equivalent:

(i) The system (2.1.1) is observable.

T
(ii) The observability matrix OM = [C CA . . . CAn−1 ] has full rank.

(iii) The matrix Z t1

T
WO = eA t C T CeAt dt
0
is nonsingular for any t1 > 0.
 
λI − A
(iv) The matrix has rank n for all eigenvalues of A.
C

(v) If (λ, y) is an eigenpair of A, then Cy 6= 0.

(vi) The eigenvalues for A − LC can be arbitrarily assigned by an appropriate choice of L.

Controllability and observability play important roles in the existence of positive definite
and positive semidefinite solutions to Lyapunov equations.
Definition 2.4.3. Let A be a stable matrix, then the matrix:
Z t1
T
CG = eAt BB T eA t dt (2.4.1)
0

is called the controllability Grammian.

Definition 2.4.4. Let A be a stable matrix, then the matrix:
Z ∞
T
OG = eA t C T CeAt dt (2.4.2)
0

is called the observability Grammian.

Theorem 2.4.3. Let A be a stable matrix. The controllability Grammian CG satisfies the
Lyapunov equation
ACG + CG AT = −BB T (2.4.3)
and is symmetric positive definite if and only if (A, B) is controllable.
Theorem 2.4.4. Let A be a stable matrix. The observability Grammian OG satisfies the
Lyapunov equation
OG A + AT OG = −C T C (2.4.4)
and is symmetric positive definite if and only if (A, C) is observable.

The following is an example illustrating the computation of both the observability and con-
trollability Grammians.
 George Kuster Chapter 2. Background 10

Example 2.4.1. Consider the simple system:

ẋ(t) = Ax(t) + Bu(t),

y(t) = Cx(t) + Du(t),

where    
−4 −8 −12 1  
A =  0 −8 −4  , B =  1  , C = 1 1 1 , D = 0.
0 0 0 1
The matrix A is stable and using the Matlab command lyap to solve the Lyapunov equa-
tions in Theorems 2.4.3 and 2.4.4 we arrive at
 
0.1458 0.0208 0.0208
CG =  0.0208 0.0833 0.0833  = OG .
0.0208 0.0833 0.0833

We note that this matrix is singular, and thus we can conclude that (A, B) is not controllable
and (A, C) is not observable.

2.5 Stabilizability and Detectability

A controller is often designed to stabilize a linear system. This is done by choosing an
appropriate control u such that the closed loop matrix A − BK is stable where K is the
feedback matrix where u(t) = Kx(t). In this section we define stabilizability of a system
and present characterizations of systems that can be stabilized through feedback control.
First consider the uncontrolled system:

ẋ(t) = Ax(t), x(0) = x0 . (2.5.1)

Definition 2.5.1. The equilibrium state of the system (2.5.1), is the vector xe satisfying

Axe = 0 (2.5.2)

If A is a nonsingular matrix, then the only solution is the trivial solution, xe = 0.

Definition 2.5.2. Let xe be an equilibrium state of (2.5.1), then xe is called

1. Stable if for each  > 0 there is a δ > 0 such that ||x(t) − xe || <  for all t ≥ 0 if
||x(t0 ) − xe || < δ.

2. Asymptotically stable if the equilibrium state is stable and δ > 0 exists such that
||x(t0 ) − xe || < δ implies lim x(t) = xe .
t→∞
George Kuster Chapter 2. Background 11

This means that all solutions which start sufficiently close to an asymptotically stable equi-
librium point will not only remain near that point, but will also approach equilibrium as
time increases. Thus if xe = 0 and A is nonsingular, the uncontrolled system above is
asymptotically stable if x(t) → 0 as t → ∞.
The eigenvalues if the matrix A in (2.5.1) characterizes the stability of the system.
Theorem 2.5.1. System (2.5.1) is said to be asymptotically stable if and only if all of the
eigenvalues of the matrix A have strictly negative real parts.
Definition 2.5.3. The matrix A is called a stable matrix if all of the eigenvalues of A
have negative real parts.

As in the case of controllability, the solution to a Lyapunov equation are also characterized
in terms of the stability of (2.5.1).
Theorem 2.5.2. The system:
ẋ(t) = Ax(t),
is asymptotically stable if and only if, for any symmetric positive definite matrix M, there
exists a unique solution to the equation

XA + AT X = −M.

Consider the linear system (2.1.1),

ẋ(t) = Ax(t) + Bu(t),

y(t) = Cx(t) + Du(t).

We consider the problem of finding a controller u which drives this system to a desired state.
That is, find a feedback matrix K such that the feedback matrix (A − BK) is stable. Not
every system can be driven to a specific state, so the importance of this section is that is
presents necessary and sufficient conditions for when a stabilizing feedback matrix exists.
Also suppose we have full information from the state, x(t), that is, at any time t information
about any variable from the state is available. Using linear control, the controller u(t) is of
the form,
u(t) = −Kx(t), (2.5.3)
where K is the feedback matrix.

Substituting u(t) = −Kx(t) into (2.1.1) the closed-loop system is

ẋ(t) = (A − BK)x(t),
(2.5.4)
y(t) = (C − DK)x(t).

Recall that the uncontrolled system (2.5.1) is stable when the matrix A is stable. Thus
stabilizing the system (2.5.4) reduces to finding a matrix K such that (A − BK) is stable. If
George Kuster Chapter 2. Background 12

such a K exists, it is called a stabilizing feedback matrix and we the system is stabilizable.
The pair of matrices (A,B) is known as the stabilizable pair.
The following theorem provides necessary and sufficient conditions for determining if a sta-
bilizing feedback matrix exists for a particular system.

Theorem 2.5.3. Just as in the case of stabilizability, the following statements are equivalent:

1. The matrix pair (A,B) is stabilizable.

2. Rank(A − λI, B) = n for all Re(λ) ≥ 0.

3. For all λ and x 6= 0 such that x∗ A = λx∗ and Re(λ) ≥ 0 one have x∗ B 6= 0.

Corollary 2.5.1. If the pair (A,B) is controllable, then the pair is also stabilizable.

We note that controllability of a system implies stabilizability of a system but stabilizability

does not imply controllability.
Just as observability is the dual of controllability, so is detectability of a system the dual
problem of stabilizability. These system properties plays a important role in control theory
in general and in particular in the results that we will use in H∞ control.

Definition 2.5.4. A system is said to be detectable if there exists a matrix L such that that
given the matrix pair (A, C), A − LC is stable.

Theorem 2.5.4. The following statements are equivalent:

1. (A, C) is detectable.
 
λI − A
2. The matrix has full rank for all Re(λ) ≥ 0.
C

3. For all λ and x 6= 0 such that Ax = λx and Re(λ) ≥ 0, we have Cx 6= 0.

4. (AT , C T ) is stabilizable.

It can be shown that is a system is controllable then the closed loop eigenvalues can be as-
signed arbitrarily (eigenvalue assignment problem). The issue with assigning the eigenvalues
arbitrarily is that there is no methodology to finding an optimal controller such that the
design characteristics are met as efficiently as possible.
This led to the development of numerous other optimal control strategies such as Linear
Quadratic optimal control. The objective of the linear quadratic control is to find a control
u which minimizes some predetermined quadratic cost function. Quadratic optimal control
is beyond the scope of this study. However, note that [7] presents a parallel development of
H∞ control and H2 optimal control which is known as Linear Quadratic Gaussian control.
George Kuster Chapter 2. Background 13

2.6 H∞ norm
First note that the transfer function of any system modeled by linear time invariant ordinary
differential equations is rational and has real valued coefficients as noted by [17] and [1]. Thus,
transfer functions of any such system have poles contained entirely in the left half plane and
are analytic in the right half plane. A space of such functions is called a Hardy space. We
shall denote this space as H∞ . For any multivariable system, the transfer function G(iω) is
a matrix. The singular values of a matrix A, σj (A) are definted as

σj (A) = λj (AAT ),

where λj (E) denotes the jth eigenvalue of the matrix E.

Note the Euclidean norm of a matrix is defined to be

||A|| = max σi (A).

i

Then we have
||G(iω)|| = max σi (G(iω)) = σmax (G(iω)).
i

Definition 2.6.1. The H∞ (induced) norm of the transfer function G(s), denoted as ||G||∞
is defined as
||G||∞ = sup σmax (G(iω)), (2.6.1)
ω∈R

where sup denotes the supremum or least upper bound over all real valued frequencies ω.
ω∈R

This provides another way to characterize the stability of a system, namely if the systems
transfer function G(s) ∈ H∞ then the system is stable. As described in [15] for the case
when G is rational, G ∈ H∞ if and only if G has no poles in the closed right half plane.
Chapter 3

Algorithms: Frequency Domain

3.1 Computing the H∞ norm

In this section we present two algorithms for finding the H∞ -norm of a system of linear
ordinary differential equations. While the algorithms use frequency domain tools, they are
important for understanding the state space algorithm presented later. Approximating the
H∞ -norm of a system requires computing the supremum of the frequency response over all
frequencies, thus the use of iterative process is required. The connection between the H∞ -
norm of a stable transfer function (i.e no poles in the right half plane) and its associated
Hamiltonian matrix plays an important role in the algorithm for the state feedback case
presented later.
Consider the system:
ẋ(t) = Ax(t) + Bu(t),
(3.1.1)
y(t) = Cx(t) + Du(t).
Definition 3.1.1. From the coefficients in the system above, we define Mr , the Hamiltonian
matrix
A + BR−1 DT C BR−1 B T
 
Mr = (3.1.2)
−C T (I + DR−1 DT )C −(A + BR−1 DT C)T )
where R = r2 I − DT D.

Recall the transfer function G(s) defined from coefficient matrices of the system (3.1.1)
G(s) = C(si − A)−1 B + D.
Theorem 3.1.1. Let G(s) be a stable transfer function and let r > 0. Then ||G||∞ < r if
and only if σmax (D) < r and Mr has no purely imaginary eigenvalues.

Theorem 3.1.1 relates the eigenvalues of Mr to the singular values of G(s). This theorem
implies that for r values larger than r∗ = ||G(s)||∞ , Mr will have no purely imaginary
eigenvalues. Further, this implies that if r < r∗ , Mr will have at least one purely imaginary
eigenvalue.

14
George Kuster Chapter 3. Algorithms: Frequency Domain 15

3.2 Bisection Algorithm

This result leads directly to the bisection algorithm presented by Boyd et al [3] where in
each iteration one checks each eigenvalue of Mr to determine the if there exists a purely
imaginary eigenvalue or not. It should be noted that this algorithm converges linearly and
the algorithm is also noted to approximate the H∞ -norm with a relative accuracy of . To
start the algorithm one needs to compute upper and lower bounds of the r-interval. One
could use rlb = 0 and set rub sufficiently large and immediately proceed with the bisection
strategy. However to obtain tighter bounds one could use Hankel singular values of the
system.

Definition 3.2.1. The Hankel singular values are the square roots of the eigenvalues of the
matrix CG OG , where CG and OG are the controllability and observability grammians, defined
by Equations 2.4.3 and 2.4.4 respectively.
Hankel singular values, denoted σHi following convention, are ordered in descending order
such that σH1 is the largest singular value.

The method for calculating rlb and rub from Step 1 (3.2.2), is known as the Enns-Glover
formula. We note that the authors suggest alternative formulas for computing the boundary
values for the bisection algorithm which can be computed with the following formulas:
n p o
rlb = max σmax (D), Trace(OG CG )/n ,
p (3.2.1)
rub = σmax (D) + 2 nTrace(CG OG )

Algorithm 3.2.1. Given system (3.1.1) with coefficient matrices A,B,C,D where A is stable
and error tolerance  > 0.
Step 1: Compute the bounds for the bisection algorithm, rlb and rub , where

P σH1 }
rlb = max{σmax (D),
(3.2.2)
rub = σmax (D) + 2 nj=1 σHi
Step 2: Set r = (rub + rlb )/2
If 2(rub − rlb ) < , end.
Step 3: Compute Mr .
Step 4: Check for purely imaginary eigenvalues of Mr .
If Mr has a purely imaginary eigenvalue, then set rlb = r.
Else set rub = r.

By Theorem 3.1.1 if we find that Mr has no purely imaginary eigenvalues, we can take
r to be the new upper bound for the algorithm. Similarly if we find that Mr does have
purely imaginary eigenvalues, we can take r as the new lower bound. In order to ensure
George Kuster Chapter 3. Algorithms: Frequency Domain 16

the computed value for r∗ is truly an accurate approximation, we stop the iterations when
2(rub − rlb ) < .
It should be pointed out that computing the starting bounds for this algorithm requires
solving two Lyapunov equations (2.4.3) and (2.4.4), multiplying two n × n matrices and
computing the eigenvalues of this product. For large matrices these computations are com-
putationally expense and are impractical. Also, this algorithm relies on very accurate com-
putations of eigenvalues, which requires very specific eigenvalue solvers, see for example [2].

Example 3.2.1. Consider system 3.1.1 where

   
−4 −8 12 1  
A=  0 −8 0  , B = 0  , C = 1 1 1 , D = 0.

0 0 −16 0
Take  = 10−6 .

We must first compute OG and CG . Using the lyap command in Matlab , we attain:
   
0.1250 0 0.1250 0.1250 0 0
OG =  0 0.0625 0  , CG =  0 0 0 .
0.1250 0 0.1250 0 0 0

Using these in Equations 3.2.1 yields:

rlb = 0.0722, rub = 0.4330,

which gives r = 0.2526 for the first iteration and

 
−4.0000 −8.0000 12.0000 15.6735 0 0

 0 −8.0000 0 0 0 0 

 0 0 −16.0000 0 0 0 
Mr =  

 −1.0000 −1.0000 −1.0000 4.0000 0 0 

 −1.0000 −1.0000 −1.0000 8.0000 8.0000 0 
−1.0000 −1.0000 −1.0000 −12.0000 0 16.0000
The eigenvalues of Mr are {8.0000, 16.0000, −0.5714, 0.5714, −8.0000, −16.0000} Since none
of the eigenvalues are purely imaginary we set r = rub and continue the iterations.

After 20 iterations we find that rlb = 0.2499 and rub = 0.2500, thus the error bound is
satisfied and the process is terminated. We arrive at ||G(s)||∞ ≈ .2500.
 George Kuster Chapter 3. Algorithms: Frequency Domain 17

3.3 Two Step Algorithm

Shortly after the Bisection algorithm was presented, a two step algorithm was proposed by
Bruinsma and Steinbuch [4] in which only the lower bound needs to be computed. It should
be noted that a similar algorithm was presented by Boyd and Balakrishnan [3] at around the
same time. It is important to note that the lower bound computation does not require finding
OG or CG , and is thus less expensive to calculate. Bruinsma and Steinbuch’s algorithm also
relies on Theorem 3.1.1 as well as the following theorem.
Theorem 3.3.1. Suppose r > σmax (D) and ω ∈ R. Then the det(Mr − ωiI) = 0 if and only
if for some n, σn (G(ωi)) = r.

As a consequence of this theorem, the imaginary eigenvalues of Mr are exactly the frequen-
cies where some singular value of G(ωi) = r. That is ωi is an eigenvalue of Mr if and only
if r is a singular value of G(ωi) where ω ∈ R. In Figure 3.1 ωi represents the frequency, r
represents the iterations and the mi represent the midpoint between each ωi and ωi+1 . Here
we see Theorem 3.3.1 allows for relating the purely imaginary eigenvalues back to multiple
values of r. Specifically each mi relates back to a value for r. In Figure 3.1 we can see that
the next iteration of r is taken to be the maximum of the singular value of G(s) for each mi .

The algorithm by Bruinsma and Steinbuch builds on the ideas of the bisection algorithm,
but uses Theorem 3.3.1 to search for ||G(s)||∞ using multiple values of ω. The bisection algo-
rithm only searches at one frequency per iteration, where as the two step algorithm searches
multiple frequencies per iteration. This is the main advantage to using this algorithm. The
authors claim this algorithm converges quadratically (they also claim ri < ri+1 , though no
proof of this could be determined).

 
  !!!!!!  
!!

!
!(

  !!
!"

   (!"
)      

!!!!   )    
   
 
 
 
 
     
  !!     !!     !!     !!    
  !!     !!    
!!    
   
Figure 3.1: Relationship between singular values and eigenvalues
George Kuster Chapter 3. Algorithms: Frequency Domain 18

To begin the algorithm let

rlb = max{σmax G(0), σmax G(ωp i), σmax (D)}, (3.3.1)
for finding the lower bound, where ωp = |λi | and λi is a pole of the transfer function G(s)
chosen based on the following criteria:
If G(s) has one or more poles with imaginary parts, then λi is the pole that maximizes the
value  
 Im(λi ) 1 
 Re(λi ) λi  . (3.3.2)
 

If G(s) has poles which are strictly real valued, then λi is the pole which maximizes | λ1i |.
As is the case with the bisection algorithm, this method is for systems where A is stable.
Algorithm 3.3.1. Given system (3.1.1) with coefficient matrices A,B,C,D where A is stable
one can find the H∞ norm using the following steps.
Step 1.1: Compute the transfer function G(s) of the given system.
Step 1.2: Find the poles of G(s).
If the poles are strictly real valued, set ωp = max|λi |.
Else set ωp = λi , where λi maximizes (3.3.2).
Step 1.3: Compute rlb using (3.3.1).
Step 2.1: Compute r = (1 + 2)rlb .
Step 2.2: Compute Mr from (3.1.2).
Sort all of the purely imaginary eigenvalues of Mr , and label them in descending
order ω1 , . . . , ωk .
If Mr has no purely imaginary eigenvalues, set rub = r and proceed to Step 3.
Step 2.3: For i = 1 to k − 1
(i) Compute mi = 21 (ωi + ωi+1 ).
(ii) Compute σmax (G(mi i)) = svdi .
Set rlb = max(svdi ).
Step 3: Compute kGk∞ = 21 (rlb + rub ).

Since the H∞ norm is defined as the maximum singular value along the imaginary axis, we
can sample along the imaginary axis to find the lower bound for the H∞ norm. To ensure
kGk∞ can be calculated using the two step algorithm, one must carefully choose a lower
bound for the algorithm.

After computing the eigenvalues of Mr , one must determine if r is an upper or lower bound,
this is accomplished using Theorem 3.1.1, similar to the bisection algorithm. However, if one
finds purely imaginary eigenvalues ω1 , ω2 , . . . , ωn , by Theorem 3.3.1 we know r is a singular
value for each G(ωi ). This implies there exists frequencies m1 , . . . , mn−1 , each mi between
its respective ωi , ωi+1 such that for some mi , G(ωi) < G(mi ). By computing σmax (G(mi )
for all i, the algorithm searches over a much larger area than the bisection algorithm. After
finding σmax (G(mi ) for all i, we set this value equal to rlb , and start the next iteration.
George Kuster Chapter 3. Algorithms: Frequency Domain 19

Example 3.3.1. Let

 
−1.0000 + 1.0000i −1.0000 −1.0000
A= −1.0000 −2.0000 + 1.0000i −1.0000 .
−1.0000 −1.0000 −2.0000 − 1.0000i

The poles of G(s) are the eigenvalues of A. We find the eigenvalues of A to be {−3.3589 +
0.2858i, −0.3628 + 0.9534i, −1.2784 − 0.2392i}. Thus A is stable.
Next we find that −0.3628 + 0.9534i, maximizes (3.3.2) and using (3.3.1) we
calculate rlb = 0.9907.
Using the above value for rlb we find,
 
−1 + i −1 −1 1.0188 0 0
 −1 −2 + i −1 0 0 0 
 
 −1 −1 −2 − i 0 0 0 
Mr =   −1
.
 −1 −1 1+i 1 1 

 −1 −1 −1 1 2+i 1 
−1 −1 −1 1 1 2−i

Computing the eigenvalues of Mr we get:

−3.1927 + 0.2164i,
3.1927 + 0.2164i,
−1.3224 − 0.1563i,
1.3224 − 0.1563i,
−0.0000 + 0.8597i,
0.0000 + 1.0201i
We find and order the purely imaginary eigenvalues, w1 = −0.0000 + 0.8597i and
w2 = 0.0000 + 1.0201i. Finding the average of these two gives m1 = 0.0000 + 0.9399i.
The maximum singular value of G(m1 i) = 1.0122 = rlb .
We then start the next iteration since there are strictly imaginary eigenvalues of Mr . After
a total of four iterations we arrive at kGk∞ = 1.0121.

Remark: We note one major issue which arrises in implementing this algorithm is the need
for an accurate eigenvalue solver since the depends on the decision if the eigenvalues are
purely imaginary or not.
Chapter 4

Algorithms: State space Approach

Doyle et al [7] presented a state space characterization for the suboptimal H∞ control prob-
lem. This characterization includes the requirement of the unique stabilizing solutions to two
Algebraic Riccati Equations to be positive semidefinite as well as a spectral radius coupling
condition.

4.1 The Hamiltonian and Associated Riccati Equation

We first review the relationship between an Algebraic Riccati Equation and the associated
Hamiltonian matrix and introduce some notation.
The Hamiltonian matrix,  
A −S
H= ,
−Q −AT
where A, Q, S ∈ Rn×n and both Q and S are symmetric has an associated Algebraic Riccati
Equation (ARE),
XA + AT X + Q − XSX = 0.
If H has no eigenvalues on the imaginary axis, H will have n eigenvalues with positive real
parts and n eigenvalues with negative real parts. Denote the n-dimensional stable invariant
subspace of H, by χ− (H), which corresponds to the eigenvalues in the open left half plane
C− . We write  
T1
χ− (H) = Im ,
T2
where T1 and T2 ∈ Rn×n and Im(A) denotes the column space of A.

Theorem 4.1.1. Let (A, B) be stabilizable and (A, Q) be detectable, then the Hamiltonian
matrix H, has n eigenvalues with negative real parts, no eigenvalues on the imaginary axis
and n eigenvalues with positive real parts. Further the associated Riccati equation has a

20
George Kuster Chapter 4. Algorithms: State Space Approach 21

unique stabilizing (symmetric) solution X, and the eigenvalues of A − BK are the stable
eigenvalues of H.

If H satisfies the following two properties,

1. Λ(H) ∩ C0 = ∅ (stability property) where Λ(H) denotes the spectrum of H.

2. The two spaces  

T1
χ− (H) = Im
T2
and  
0
Im
I
are complimentary, in other words T1 is nonsingular;

we can set X = T2 T1−1 . Then X is uniquely determined by H and we say that H ∈ dom(Ric).
If H ∈ dom(Ric) then H has no purely imaginary eigenvalues and the matrix T1 is non-
singular. That is Ric : H → X and Ric(H) = X. One can verify a given H exists in
dom(Ric) with the following lemma from Doyle et al [7].

Lemma 4.1.1. Suppose H ∈ dom(Ric) and X = Ric(H). Then:

1. X is symmetric;

2. X satisfies the ARE XA + AT X + Q − XSX = 0;|||

3. A − SX is stable.

Recall the system,

ẋ(t) = Ax(t) + Bu(t),
(4.1.1)
y(t) = Cx(t) + Du(t),

with the Hamiltonian matrix associated with the system

− r12 BB T
 
A
H(r) = . (4.1.2)
−C T C −AT

Note H(r) is equal to Mr given by (3.1.2), where D = 0.

The Riccati equation associated with the Hamiltonian matrix H(r) is

XA + AT X + Q − XSX = 0,
1
where S = r2
BB T and Q = −C T C.
George Kuster Chapter 4. Algorithms: State Space Approach 22

Recall that the Hamiltonian matrix Mr , is connected to ||G||∞ through Theorem 3.1.1.
Thus, H(r) has no purely imaginary eigenvalues when ||G||∞ < r. Theorem 4.1.1 connects
the Hamiltonian matrix H(r) with its associated ARE and provides conditions under which
H(r) ∈ dom(Ric), or has no purely imaginary eigenvalues. Theorem 4.1.1 further connects
||G||∞ with the existence of a stabilizing solution X. Specifically it lays the groundwork for
determining necessary and sufficient conditions for the existence of an admissible controller
in a state space formulation.
That is if (A, B) is stabilizable and (A, Q) is detectable, then H(r) ∈ dom(Ric), hence there
exists matrices T1 and T2 such that,
 
T1
χ− (H) = Im ,
T2

and X = T2 T1−1 is the stabilizing solution to the associated ARE with X ≥ 0.

Since H(r) ∈ dom(Ric), H(r) has no purely imaginary eigenvalues and by Theorem 3.1.1,
||G||∞ < r. Combining Theorems 3.1.1 and 4.1.1, we have if the stabilizing solution to the
associated ARE exists, then ||G||∞ < r. Thus given system (4.1.1) the problem of finding
the smallest value of r such that H(r) has no purely imaginary eigenvalues becomes finding
the smallest value of r such that a stabilizing solution to the ARE exists where X = T2 T1−1
with X ≥ 0.

Output Feedback

We now consider the output feedback case. A finite-dimensional time-invariant linear continuous-
time dynamical system may be described using the following system of first-order ordinary
differential equations:
ẋ(t) = Ax(t) + B1 w(t) + B2 u(t)
z(t) = C1 x(t) + D12 u(t) (4.1.3)
y(t) = C2 x(t) + D21 w(t).
Here x(t), y(t), z(t) are the state, sensed output and controlled output vectors respectively
and u(t), w(t) are the control and disturbance inputs respectively.
Assume throughout the discussion that follows that the matrices A, B1 , B2 , C1 , C2 , D12 , D21
are real and satisfy the following assumptions:

(A1) (A, B1 ) is stabilizable and (C1 , A) is detectable;

(A2) (A, B2 ) is stabilizable and (C2 , A) is detectable;

T
(A3) D12 [C1 D12 ] = [0 I];
T
(A4) D21 [B1T D21
T
] = [0 I].
George Kuster Chapter 4. Algorithms: State Space Approach 23

The two Hamiltonian matrices associated with system (4.1.3) are

1
 T T

A r 2 B1 B1 − B2 B2
H∞ (r) := , (4.1.4)
−C1T C1 −AT

and its dual,

1
 
AT CT C
r2 1 1
− C2T C2
J∞ (r) := . (4.1.5)
−B1T B1 −A

The two Algebraic Riccati Equations associated with the system are:

1
XA + AT X − X(B2 B2T − B1 B1T )X + C1T C1 = 0, (4.1.6)
r2
and its dual,
1 T
AY + Y T A − Y (C2T C2 − C C1 )Y + B1 B1T = 0. (4.1.7)
r2 1
The four assumptions (A1)—(A4) are standard for control applications. They were originally
formulated for H2 control and since H∞ control was developed in parallel with H2 control
these assumptions followed [7]. Assumption (A1) is necessary for the existence of a stabilizing
controller in the form u = Ky. Assumption (A2) is sufficient for the closed loop system to
be both stabilizable and detectable. If the closed loop system is not stabilizable, then it is
certainly not controllable and therefore the H∞ optimal control problem has no solution.
Assumptions (A3) and (A4) imply there is no overlap between the cost of the control input
and the state. In other words the cost penalty applied to z contains a non-singular and
normalized cost penalty on the control input u. Assumptions (A3) and (A4) can be relaxed
but they are made in an effort to obtain cleaner formulas.
Recall the Hamiltonian matrix H(r) given by (4.1.2), associated with system (4.1.1). It
is important to note the difference between H(r) and H∞ (r) given by (4.1.4). Specifically
the ability to make r12 B1 B1T arbitrarily large or small no longer guarantees a sign definite
matrix in the [1, 2] block of H∞ (r). As a consequence one is not guaranteed the existence
of a stabilizing solution to the ARE (4.1.6) even when H∞ (r) ∈ dom(Ric). In other words
H∞ (r) ∈ dom(Ric) is still a necessary condition for the existence of a stabilizing solution,
but it is now no longer a sufficient condition.

Doyle et al [7] shows that system (4.1.3) which satisfies (A1)—(A4), has an admissible
controller K, such that ||Tzw ||∞ < r if and only if the following necessary and sufficient
conditions are met:

(C1) H∞ (r) ∈ dom(Ric) and X(r) := Ric(H∞ (r)) ≥ 0;

(C2) J∞ (r) ∈ dom(Ric) and Y (r) := Ric(J∞ (r)) ≥ 0;

(C3) ρ(X(r)Y(r)) < r2 .

George Kuster Chapter 4. Algorithms: State Space Approach 24

where Tzw denotes the transfer function from the disturbance w to the controlled output z,
such that
Tzw = C1 (sI − A + B2 K)−1 B1 + D12.

Lin, Wang, Xu [16] present an algorithm to find ||Tzw || using a state space approach. The
algorithm breaks the three conditions (C1)—(C3) into three seperate sub-problems.

1. Finding the smallest r value for which (C1) holds rx , is performed first. In doing so, rx
is distinguished by the singularity of the T1 matrix created by the basis vectors of the
stable invariant subspace of H∞ (r). Using a reformulation (in terms of singularities) of
(C1) a secant method is applied (rather than a Newton type method) directly utilizing
the stable invariant subspace of H∞ (r).
2. The same technique is then applied to J∞ (r), to find ry , the smallest value for r such
that (C2) holds.
3. The geometric nature (strictly decreasing) of ρ(X(r)Y (r)) is then exploited to find r∗ ,
the value for which all three conditions (C1)—(C3), hold.

To find rx the smallest value such that condition (C1) holds, we must first find the interval
such that H∞ (r) has no purely imaginary eigenvalues. It is a well known result that there
is only one such interval, (rlb , ∞) [10], [12], [21]. Furthermore it is a well known result that
rlb = ||GZ ||∞ , where
GZ (s) = B1T (sI − (A − ZC1T C1 )T )−1 C1T ,
with Z = Ric(H2 ) and  
AT −C1 C1T
H2 = .
−B2 B2T −A
Thus to find rlb we can apply either the bisection or two step method to the transfer function
G(iω) = B1T (iωI − (A − ZC1T C1 )T )−1 C1T . (4.1.8)
Alternatively a bisection approach could be used directly on H∞ (r). Upon calculating the
eigenvalues of H∞ (r), one could verify Re(λi ) 6= 0 and im(λi ) = 0 for 1 ≤ i ≤ n. If
Λ(H∞ (r̃)) ∩ C0 = ∅ then r̃ < rlb and if Λ(H∞ (r)) ∩ C0 6= ∅ then r̃ < rlb .

To satisfy the second condition of (C1), X(r) ≥ 0 we note that as a consequence of the
existence of only one interval (rlb , ∞) such that Λ(H∞ (r)) ∩ C0 = ∅ then once rlb is found
H∞ (r) has no purely imaginary eigenvalue for all r > rlb . Then one can find matrices T1 and
T2 which correspond to the stable eigenvalues of H∞ (r) where
 
T1 (r)
χ− (H∞ )(r) = Im ,
T2 (r)

and if T1 (r) is non-singular set X(r) = T2 (r)T1 (r)−1 . When such an X exists and is pos-
itive semi-definite for some r̃ then it positive semi-definite for all r > r̃, and is monotonic
decreasing by the following theorem from [12].
George Kuster Chapter 4. Algorithms: State Space Approach 25

Theorem 4.1.1. If a non-zero X(r̃) = Ric(H∞ (r̃)) exists and is positive semidefinite for
some r̃ > 0, then X(r) = Ric(H∞ (r)) exists and is positive semidefinite for all r ∈ [r̃, ∞)
and obeys the partial order 0 ≤ X(r2 ) ≤ X(r1 ) for r̃ ≤ r1 ≤ r2 < ∞.

Assume 0 < rx < ∞, Theorem 4.1.1 implies the existence of a least r such that X(r) ≥ 0,
where A ≥ 0 denotes the matrix A is positive semi-definite. As a result, one could implement
a bisection algorithm to find such an r noting if X(r) ≥ 0 then r > rx and if X(r)  0 then
r < rx . The following theorem from Lin et al (1998), is the keystone to this, as it provides
the means in which one can implement a bisection or secant method to search for this r. In
fact the following theorem is a reformulation of (C1) from [7].

Theorem 4.1.2. For r > rlb , let X(r) = Ric(H∞ (r)) whenever H∞ (r) ∈ dom(Ric). Suppose
that X(r2 ) ≥ 0 and X(r1 )  0 such that r2 > r1 > rlb . Then the following are true:
(i) There exists at least one r ∈ (r1 , r2 ) such that T1 (r) is singular.
(ii) Let
rs = sup{r|T1 (r) is singular}. (4.1.9)
Then X(r) ≥ 0 for r > rs and X(r)  0 for rlb < r < rs , when X(r) exists.

Proof(outline)
(i) Assume that T1 (r) is nonsingular for all r ∈ (r1 , r2 ).

• Since X(r2 ) ≥ 0 and X(r1 )  0, by continuity of eigenvalues there exists r̃ ∈ (r1 , r2 )

such that X(r̃) is singular and X(r) ≥ 0 for all r > r̃.

• By the monotonicity result in Theorem 4.1.1, we have 0 ≤ X(r) ≤ X(r̃), which implies
the smallest eigenvalue of X(r) is zero and further this applies to all r ≥ r̃.

• However this contradicts the fact that X(r2 ) ≥ 0, thus there exists an r ∈ (r1 , r2 ) such
that T1 (r) is singular.

(ii) Suppose there exists α > rs where X(α)  0, then by Theorem 4.1.1, we have α < r2 .
Thus from (i) there must exist β ∈ (α, r2 ) with X(β) singular. This contradicts the definition
of rs , thus X(r) ≥ 0 for all r > rs .
On the other hand, suppose there exists γ ∈ (rlb , rs ) with X(γ) ≥ 0, then by Theorem 4.1.1
we have X(r) exists for all r ≥ γ. However this is a contradiction because X(rs ) does not
exist. Thus X(r)  0 for r ∈ (rlb , rs ).
Assume that rx < ∞, and rlb < rs .

Then there exists at least one r ∈ (rlb , ∞) such that T1 (r) in,
 
T1 (r)
χ− (H∞ (r)) = Im ,
T2 (r)
George Kuster Chapter 4. Algorithms: State Space Approach 26

is singular and rx = rs , where

rs = sup{r|T1 (r) is singular}.

Thus the problem of finding the smallest value for r such that (C1) holds is now reformulated
into finding rs , the smallest value such that X(r) ≥ 0 for r > rs where rs ≥ rlb > 0.

Note rlb ≤ rs and H∞ (r) ∩ C0 = ∅ for all r > rlb , then by assumptions (A1) and (A2)
and Theorems 4.1.1 and 4.1.2, the matrices T1 and T2 are defined, T1 is non singular and
T2 (r)T1−1 (r) = X(r) ≥ 0 for all r > rs . Thus rs can be approximated with a bisection
strategy, checking the positive semi-definiteness of the matrix T2 (r)T1−1 (r). Explicitly, if
T2 (r)T1−1 (r) ≥ 0 then r > rs and if T2 (r)T1−1 (r)  0 then r < rs .
To increase the rate of convergence,one could apply a variant secant method to search for rs
as well. By Theorem 4.1.2, rs is the largest value such that T1 (r) is singular, then our goal
is to find the largest r such that f (r) = 0, where

f (r) = σmin (T1 (r)).

However singular value plots are non-negative functions and by Theorem 4.1.2 f (r) may
have multiple zeros. For this reason the secant method is not guaranteed to converge using
f (r).

Thus we actually apply the secant method to the new smooth function

−σmin (T1 (r)) if r ≤ rs
f (r) = ,
σmin (T1 (r)) if r ≥ rs

that is, f along with its reflection over the r axis for r < rs . One issue which arrises from
this new function is now one must deduce which case to use, that is f or −f . To make this
determination we appeal once more to Theorem 4.1.2, and check if X(r) = T2 (r)T1−1 (r) ≥ 0
by computing the eigenvalues of the matrix T1T (r)T2 (r) which is congruent to T2 (r)T1−1 (r)
and avoids the inverse calculation for T1 (r).
This leads to an algorithm for the state feedback case. The following is an outline of the
algorithm.
Algorithm 4.1.1. This algorithm computes rx , the smallest value for r such that condition
(C1) holds.
Step 1:
Attain values for r− and r+ , using equations 4.1.10 - 4.1.12.
Compute ΛH∞ (r̃), the eigenvalues of H∞ (r̃), where r̃ = 21 (r− + r+ )
 George Kuster Chapter 4. Algorithms: State Space Approach 27

!(!) = !!"# (!! (!))    

!(!! )  

!(!! )  

r1 rnew     r2    
rs    

−!(!! )  

Figure 4.1: Variant secant method

Step 2:
Determine the next iterations interval [r− , r+ ].

Step 2.1
If the matrix H∞ (r̃) has purely imaginary eigenvalues then set r− = r̃ and move to Step 3.

Step 2.2
If the matrix H∞ (r̃) has no purely imaginary eigenvalues then compute the eigen-matrix
[T1 T2 ]T and set case = 1.
If T2 T1−1 is positive semidefinite, then set r+ = r̃.
If T2 T1−1 is not positive semidefinite, then set r− = r̃.

Step 3: Update r̃.

If r+ − r− ≤ T ol, then use the bisection method.
If r+ − r− > T ol, then use the secant method method, with functions based on the case value
from Step 2.

Step 4: Stop Criteria:

If any of the stop criteria hold, then accept r∗ accordingly
George Kuster Chapter 4. Algorithms: State Space Approach 28

Else return to Step 1.

Step 5: Double check:

Perform Step 1 and Step 2 with r̃ = r∗ .
If T2 T1−1 is positive semi definite find ∆r > 0 such that T1 (r∗ + ∆r) is well conditioned.
If T (r∗ + ∆r)2 T1 (r∗ + ∆r)−1 is positive semi-definite then accept r∗ + ∆r = r∗ .
Else set r− = r∗ + ∆r and return to Step 1.

If T2 T1−1 is not positive semi definite set r+ large and return to Step 1.

We will now explore the complexity of each of the steps above, one at a time.

Step 1
To begin the algorithm one must find an initial interval on which to search for rx , the smallest
value for r for which (C1) holds. In practical cases assuming B2 is not singular, to compute
the upper bound of this interval we can set

σmax (B1 )
r+ = , (4.1.10)
σ0 (B2 )

where σmax is the largest singular value and σmin is the smallest singular value.
Much like the previous algorithms, one could use the Hankel singular values for the upper
and lower bounds as well, where

r− = σH1 , (4.1.11)
X
r+ = 2 σHi . (4.1.12)
i

Note However depending on the system it might be more efficient to set r− = 0 use 4.1.10
for r+ .

Note, if one computes rlb using the two step method or the bisection method on the system
response matrix GZ (iω), 4.1.8 then Step 2.1 below can be avoided since ||GZ ||∞ = rlb , thus
establishing the interval on which ΛH∞ (r) ∩ C0 = ∅ before starting this algorithm.

Step 2
In this step we determine the next upper bound or lower bound of the interval by checking
if the value for r̃ is above or below rx . This involves testing values of r to verify they satisfy
George Kuster Chapter 4. Algorithms: State Space Approach 29

condition (C1), ΛH∞ (r) ∩ C0 6= ∅ and Ric(H) ≥ 0. We must first find rlb , then we seek rs .
We note that rs may or not be the same value as rlb .

Step 2.1, ΛH∞ (r) ∩ C0 6= ∅

Since there exists only one such interval, if ΛH∞ (r) ∩ C0 6= ∅ then r̃ < rlb and we set r̃ = r− .
Doing so requires accurately computing the eigenvalues of a 2n × 2n Hamiltonian matrix.
Error in these calculations could severely effect the choice of r− since it must be determined
if Re(λi ) < T ol and im(λi ) > T ol for all 0 ≤ i ≤ n.

Step 2.2, ΛH∞ (r) ∩ C0 = ∅

Since there can be values for r such that the Hamiltonian matrix H∞ (r) may have no eigen-
values on C0 , but X(r)  0, we must test the positive semi-definite property of T2 T1−1 (r̃).
Here we use the singular values of T1 (r̃) to determine if the value for r̃ is greater than or less
than rx . Theorem 4.1.1 implies that one can impliment a bisection algorithm in the search
for rx and Theorem 4.1.2 provides the means to implement this search, namely searching
over the singular values of T1 (r̃).

There are two major components to this step. First one must compute the 2n×n eigenmatrix
 
T1 (r̃)
,
T2 (r̃)

then it must be determined if T2 (r̃)T1−1 (r̃) ≥ 0.

To do so, one can calculate the n eigenvectors of H∞ (r̃) which correspond to the stable
eigenvalues, order them in a matrix of size 2n × n and divide the matrix such that T1 is the
top half and T2 is the bottom half.
Now we must judge the positive semi-definiteness of T2 T1−1 . Since it is best to avoid the
computation of an inverse matrix, we note that T2 T1−1 andT1T T2 are congruent matrices, so
therefore we check the positive semi-definiteness of T1T T2 instead. To do so we verify that
Re(λi ) ≥ 0 for all 0 ≤ i ≤ n. Again note this requires accurate computations of eigenvalues.
By Theorem 4.1.2 if T2 (r̃)T1−1 (r̃) ≥ 0 then r̃ is an upper bound for rs so we set r̃ = r+ ,
however if it is determined that T2 (r̃)T1−1 (r̃)  0 then r̃ is an a lower bound for rs so then
we set r̃ = r− .

Step 3
In step three of the algorithm we update the next iterate of r̃. Here we use either the bisec-
tion method or secant method. We determine which method to use based on whether or not
George Kuster Chapter 4. Algorithms: State Space Approach 30

(r+ − r̃) is less than some predetermined tolerance, tol. If the difference is less than tol, we
use the secant method and if the difference is larger than tol we use the bisection method to
speed up the iterative process.

The bisection algorithm is standard, that is if (r+ − r− ) > tol set r̃ = 21 (r− + r+ ) and proceed
to Step 4 to check the stop criteria.
The secant method can be applied to both Step 2.1 (case = 0) and Step 2.2 (case =1),
however it is applied to different functions for the for the two different cases since the goal
for each case is different. For the Step 2.1, that is Λ(H∞ (r̃)) ∩ C0 6= ∅, where we let r̃ = r+
such that,
(r+ − r− ) min |λi (H∞ (r− )) − λj (H∞ (r− ))|
i6=j,i,j∈I0
r̃ = . (4.1.13)
min |R(λk (H∞ (r+ )))|+) min |(λi (H∞ (r− )) − λj (H∞ (r− ))|
1≤k≤2n i6=j,i,j∈I0

This is the secant method for the function

(
− min |λi (H∞ (r)) − λj (H∞ (r))| if r ≤ rlb ,
g(r) = i6=j,i,j∈I0 (4.1.14)
mink |R(λk (H∞ (r)))| if r ≥ rlb ,
where λk (H∞ (r)) is the k th eigenvalue of H∞ (r) and
I0 = {i|1 ≤ i ≤ 2n, λ1 (H∞ (r)) ∈ Λ(H∞ (r)) ∩ C0 }.
The main reasons for choosing g(r) as such are, as r → rlb from the right, the real parts of
λi (H∞ (r)) go to zero since some of the eiegnvalues of H∞ (r) are on C0 for r < rlb . Also as
r → rlb from the left the distance between the imaginary eigenvalues of H∞ (r) goes to zero
since the eigenvalues are moving off the imaginary axis.

If case = 1 that is, Λ(H∞ (r̃)) ∩ C0 = ∅ we let r̃ = r− where,

(r+ − r− )σmin (T1 (r− ))

r̃ = r− + . (4.1.15)
σmin (T1 (r+ )) + σmin (T1 (r− ))
We note that this is the secant method applied to the function

−σmin (T1 (r)) if r ≤ rs ,
f (r) = (4.1.16)
σmin (T1 (r)) if r ≥ rs .
It is also important to note that f (r) was chosen so that f is smooth and the secant method
will converge about rs .

Step 4
There are two types of stop criteria used in this algorithm for each case. The first stop
criterion deals with the distance between the previous and the next iteration of r, ∆r̃. The
George Kuster Chapter 4. Algorithms: State Space Approach 31

second stop criterion deals with the magnitude of the functions 4.1.14 and 4.1.15. The two
stop criteria are necessary since ∆r̃ → 0 does not directly imply r → rx . Thus we combine
the convergence of ∆r̃ with the value of equations 4.1.14 or 4.1.15 accordingly.

Suppose the r̂ = r̃ + ∆r̃ is the updated iterate of r̃. For both cases the stop criteria ∆r < tol
was used, were tol = 10−6 .
For the function convergence in case=0, we are concerned with the eigenvalues of the Hamil-
tonian matrix, H∞ (r̃). The secant method uses equation 4.1.14. Since we have g(r) → 0 as
r → rlb then we want |g(r)| < tol and move to Step 5.
For case =1 we are concerned with the singular vales of T1 . For this case we use equation
4.1.15 for the secant method. Since we want to find the smallest value of r such that T1 is
singular we want f (r) = 0, thus we stop when |f (r)| < tol and move to Step 5.

Step 5
It is important to note that the vast majority of the computations in the algorithm rely
on the various matrices being well conditioned. The computation of T2 (r)T1 (r)−1 , is not
trivial and relies on the numerical computations being accurate. This accuracy depends on
the condition of the matrices T2 (r) and T1 (r). Further the secant method applied to 4.1.15
may have converged to an incorrect value of r. Thus we must double check the algorithm
converged to the correct value and also verify the condition of T1 (r) since it must be inverted.
In doing so the actual value for rx may not be attained, but rather a more conservative value
will found.
To ensure convergence of the algorithm to rx = rs we must verify if T2 (r̃)T1−1 (r̃) ≥ 0. If this
holds then we check if T1 (r̃) is well conditioned. If the matrix is well conditioned we stop
and accept r̃ = rx .
On the other hand if T1 (r̃) is ill-conditioned, we must find the smallest perturbation of r̃
such that T1 (r) is well conditioned.
If T2 (r̃)T1−1 (r̃)  0 then we know r̃ < rx thus set r− = r̃ and r+ sufficiently large and return
to Step 1.
It is for this reason that in implementing this algorithm, it would be immensely beneficial to
accurately calculate rlb before hand. The bisection algorithm presented earlier would provide
an approximation for rlb which could then be used as the lower bound in Algorithm 4.1.1.
It should be noted that (C1) and (C2) are dual problems of each other thus the above algo-
rithm is also used to find the value of ry , that is the smallest value of r such that condition
(C2) holds, by replacing A, B1 , B2 and C1 with AT , C1T , C2T and B1T respectively.

Let rm = max{rx , ry }, then by Theorem 4 and the definitions of rx and ry , X(r) :=

Ric(H∞ (r)) and Y (r) := Ric(J∞ (r)) both exist and further X(r), Y (r) ≥ 0 for all r ≥ rm .
George Kuster Chapter 4. Algorithms: State Space Approach 32

From this we see that rm is a lower bound for r∗ , the value for which (C1)—(C3) hold si-
multaneously.

Recall: (C3) ρ(X(r)Y(r)) < r2

Define,
ρ(r) = ρ(X(r)Y(r)). (4.1.17)

To find r∗ we take advantage of the geometrical natures of ρ(r) and r2 . Specifically r2 is and
increasing function of r and ρ(r) is decreasing function of r [12]. One important advantage
of using this geometrical nature is that an upper bound can be easily found.
In implementing a bisection algorithm to approximate the value of r such that (C3) holds,
we note that ρ(r) > r2 for all r < r∗ and ρ(r) < r2 for all r > r∗ . Thus if ρ(rm ) ≤ rm
2
, then
∗ 2
r = rm and we stop. If howeverp ρ(rm ) > rm , by considering the nature of the two curves
we have r∗ > rm and further ρ(rm ) must be an upper bound for r∗ thus we set r+ = rm .
To numerically implement this, we use the standard bisection or secant methods on the
function
f (r) = ρ(r) − r2 .

We stop the iterations when −f (r) < tol, that is when r2 − ρ(r) is positive and sufficiently
small. The following algorithm can be used to compute r∗ , the smallest value of r such that
conditions (C1)—(C3) hold simultaneously.
Algorithm 4.1.2.
Step 1: Find rx and ry using Algorithm 4.1.1.
Step 2: Set r̃ = rm .
Step 3: Compute ρ(r̃).
If ρ(r̃) ≤ r2 , set r∗ = r̃ and stop. p
If ρ(r̃) > r2 , then set r− = r̃ and r+ = ρ(r̃).
Step 4: Update r̃.
If r+ − r− ≥ tol, use the bisection method.
If r+ − r− < tol use the secant method applied to the function ρ(r) − r2 .
Step 5: Stop Criteria
If −f (r̃) < tol holds then set r∗ = r̃ and stop
Else return to Step 4.

4.2 Numerical Examples

The following numerical examples are taken directly from [16]. It should be noted the first
two were originally from [19]. All times were computed using MATLAB’s ’tic toc’ command.
George Kuster Chapter 4. Algorithms: State Space Approach 33

We note that the number of iterations for the examples is much higher than those found
in [16]. This is due to their use much more efficient and accurate numerical solvers, in com-
parison to the use of basic MATLAB commands here. It is noted however that although
the tools used here could be considered inferior, the algorithm still converged to the correct
values.

Example 1
Consider the system,
ẋ(t) = Ax(t) + B1 w(t) + B2 u(t),
z(t) = C1 x(t) + D12 u(t),
y(t) = C2 x(t) + D21 w(t),

Table 4.1: Numerical example 1

Time Iterations Computed value

rlbx 0.707106
rx .029 sec 39 0.707106
rlby 0.707106
ry .314 sec 39 0.707106
r∗ .015 sec 40 0.73205
0.358 sec

This shows the importance of using accurate solvers quite remarkedly. Even though rlb is
approximately equal to rx and ry , in both cases the algorithm required a significant number
of iterations to converge.

Example 2
Consider system 4.1.3 with the following coefficient matrices,

−1 1
 
 
A= , B1 = 0 I , B2 = I,
0 4
   
I I  
C1 = , C2 = I, D12 = , D21 = 0 I .
0 0
Here we see the value of rlbx need not equal rx . This also shows that the r value with
satisfies (C1) and (C2) is a lower bound for r∗ , but need not be close to r∗ .
George Kuster Chapter 4. Algorithms: State Space Approach 34

Table 4.2: Numerical example 2

Time Iterations Computed value

rlbx 0.71837
rx .0297 sec 39 .99̄
rlbx 0.71837
ry .2686 sec 32 .99̄
r∗ .049 sec 88 8.14
.298 sec

The following example was presented as a counter example to the arguement in [7] that (C3)
would fail before (C1) or (C2).

Example 3
Now consider the system with the following coefficient matrices,

− 12 − 21
 1 3 
− 75 7
− 10
 
5 2
A= 3 , B2 = , C2 = 7 ,
0 10
− 32 0 − 12 − 10
   
T 1 0 0 0 T 0 0 1 0
B1 = C1 = , D12 = D21 = .
0 1 0 0 0 0 0 1

Table 4.3: Numerical example 3

Time Iterations Computed value

rlbx 0.69641
rx .277 sec 12 0.6964
rlby 2.1804
ry .432 sec 16 2.767
r∗ .022 sec 29 3.413
.731 sec

Note that when r = 2.25, we have

ρ(2.25) ≈ 4.7916 < 5.0625 = r2 .
However the eigenvalues of Y (r) are .44486 and −6.11457. Thus (C2) fails before (C3). This
is important to note because some algorithms such as that presented by [19] only checks
condition (C3).
Chapter 5

Conclusions

5.1 Conclusion on State Space Formulation

The most intriguing aspect of the paper by Lin, Wang and Xu is their use of frequency
domain tools in the state space formulation. The state space formulation was originally
developed in an effort to move away from the frequency domain approach to H∞ -control.
Yet, Lin Wang and Xu reformulate some of the state space necessary and sufficient conditions
into a statement about singular values.
A major issue which arrises from this reformulation is the accurate calculation of eigenvalues.
It should be noted that for the purposes of this study, MatLab functions were used to compute
all eigenvalues and singular values. This however was not the case in [16] where accurate
eigen-solvers were used.
Accurate computations of eigenvalues is paramount here. The reason for this is the decision
making process inherent in the algorithm presented in Chapter 4. The symmetric properties
of Hamiltonian matrices is the key factor in how the state space formulation operates. Many
eigen-solvers are not capable of preserving this property and as a result loose accuracy in
their computations. In turn because the algorithm here relies on answering yes or no type
questions and removing entire intervals from a search, false negatives or positives could have
devastating results. For instance suppose the algorithm converges to a value such that it
indicates a controller exists for a larger error disturbance than it can actually accommodate.
One may waste valuable resources implementing a control only to find it doesn’t work when
applied.
Future research plans include implementing more accurate eigen-solvers and applying this
algorithm to a large scale model from application type setting. It would be beneficial to see
if this algorithm could be made competitive when compared to other methods currently in
use for estimating the H∞ norm of a system.

35
George Kuster Chapter 5. Conclusions 36

5.2 Extending to Mathematics Education

Since I will be pursuing my Doctorate degree in Mathematics Education, it is pertinent to
view this project as a learning experience not only pertaining to the content of this project
but also how the learning itself occurred. Reflecting back on the process of performing re-
search on more or less my own, there are a few interesting connections noted with various
learning theories. Piaget believed that learning occurred through the complimentary pro-
cesses, assimilation and accommodation. While studying the various topics which encompass
H∞ control, there were numerous occasions which I felt I understood a topic only to find,
upon explaining it to a more knowledgeable other, that the way I assimilated the topic was
not consistent with H∞ control theory. In this instance I was presented with a state of
dis-equilibrium which required accommodation, that is modification to existing schemes.
Often times however it was noted that certain concepts would repeatedly be sources of
confusion, some continuing until late in the research process. The confusion persisted, even
through repeated efforts to clarify the misunderstandings. It was evident during the multiple
discussions about these topics, that the way in which I understood them was not correct,
yet I (perhaps subconsciously) did not accommodate for this new information. There is
a plethora of possibilities which could explain why this did not happen, ranging from the
inability of the more knowledgable sources to relate to my previously established knowledge,
to me not paying close enough attention. However, it might be interesting to apply Tall
and Vinners idea of Concept Image and Concept Definition [22] to this phenomenon. That
is, even though these topics appeared to consistently be sources of confusion, it might be
possible that only the under developed pieces of the overall concept image were being evoked
at any given time.
While it is impossible for a person’s entire concept image to be evoked at any given time,
it has been interesting looking back and reflecting on how I once thought of certain topics
presented in this thesis and comparing this to how I now think about them. For instance
when I first began research on H∞ control, when I read the words “transfer function” I would
simply relate it to a matrix. Now the image is much more developed consisting of concepts
such as linear mappings.
Building on Sfard’s work concerning processes which work on previously established objects,
Zandieh defined process-object pairs [24] in which a series develops. That is roughly, a process
is established as an object which can then be acted on by other processes. While this was
initially developed around topics covered in calculus, it could certainly be extended to what
occurred throughout this research project. Take for instance the transfer function, it could
simply be seen as an object, that is a matrix developed from the coefficients of our system
matrices. It can also be viewed as a ratio, specifically the ratio of the output to the input.
At arguably its deepest level, it could also be seen as a function mapping the inputs to the
outputs. It would be interesting to see how these layers developed during specific processes
conducted in working on this thesis.
Upon reflecting on the research experience there are numerous benefits to completing a Mas-
ters thesis which can be applied to the field Mathematics Education. All too often there
George Kuster Chapter 5. Conclusions 37

is an impression that researching areas in applied mathematics isn’t beneficial for one who
will ultimately study mathematics education. After this experience, I realize this couldn’t
be further from the truth. Not only does completing a Masters Thesis play a key role in
preparing a student for independent research, but further it instills a true appreciation for
what it means to learn and understand a concept. Further, I have gained an appreciation
for my own learning process which is certainly beneficial to understanding how others learn.

Constructivism in Research
Research at the Masters level is a far more independent venture than one typically gets to
experience as an undergraduate. The advisor’s role is to ensure the student stays on track
and their work is appropriate for the task presented. The advisor does not explain or lecture
the student topics concerning the material past the point of redirection or filling gaps in the
student’s understanding. Hopefully, input from the student is met with interest, apprecia-
tion and consideration. That is, the student’s input is valued. The use of more guidance
and less lecture results in the student forming much more personal knowledge concerning
the subject, compared to an assimilation of the lecturer’s knowledge.

Knowing Verus Understanding

One can know that a certain continent is on the other side of the earth, but without ex-
periencing the trip, one don’t have first hand understanding of what this means. Research
provides first hand knowledge of a specific concept, while capturing ones motivation and
building on previous experiences. The learner is actively engaged in constructing their own
knowledge, while simultaneously building confidence in their ability to learn. This increase
in confidence in turn leads to a more actively engaged learner. The important experience
from research is not reading about or getting to know the concepts being researched. The
value of research comes from being actively immersed in the topic of interest. The knowledge
built from participating in this study is no different. The largest amount of knowledge did
not come from reading relevant topics, it came from discussions with those involved in it,
actively implementing algorithms, troubleshooting those algorithms and motivation to seek
connections.

Deeper Understanding
Often times courses are structured in a linear fashion. That is, the course has been designed
to present topics in a predetermined order. This order may or may not be consistent with
the order in which topics would naturally fall and as a consequence may lack motivation.
For illustration purposes only, it is customary to teach integration before infinite series in
many university calculus sequences. Though it may be more natural to learn sequence and
series before integration. While actively participating in research one must learn topics as
they arise, not as they are presented. This provides a significant increase in motivation to
learn these topics
Overall completing this Masters Thesis had a large impact on not only my current under-
standing of control theory, but also on my future work as a mathematics education researcher.
Bibliography

[1] T. Başar and P. Bernhard. H-infinity optimal control and related minimax design prob-
lems: a dynamic game approach. Birkhäuser Boston, 2008.
[2] P. Benner, V. Mehrmann, and H. Xu. A numerically stable, structure preserving method
for computing the eigenvalues of real hamiltonian or symplectic pencils. Numerische
Mathematik, 78(3):329–358, 1998.
[3] S. Boyd, V. Balakrishnan, and P. Kabamba. A bisection method for computing the H∞
norm of a transfer matrix and related problems. Mathematics of Control, Signals, and
Systems (MCSS), 2(3):207–219, 1989.
[4] NA Bruinsma and M. Steinbuch. A fast algorithm to compute the H∞ norm of a transfer
function matrix. Systems & Control Letters, 14(4):287–293, 1990.
[5] B. Datta. Numerical methods for linear control systems. Academic Press, 2003.
[6] J.C. Doyle et al. Advances in multivariable control. In Lecture Notes at ONR/Honeywell
Workshop. Minneapolis, MN, 1984.
[7] J.C. Doyle, K. Glover, P.P. Khargonekar, and B.A. Francis. State-space solutions to
standard H2 and H∞ control problems. Automatic Control, IEEE Transactions on,
34(8):831–847, 1989.
[8] B.A. Francis and J.C. Doyle. Linear control theory with an H∞ optimality criterion.
SIAM Journal on Control and Optimization, 25(4):815–844, 1987.
[9] P. Gahinet. On the game riccati equations arising in H∞ control problem. SIAM Journal
on Control and Optimization, 32(3):635–647, 1994.
[10] P.M. Gahinet and P. Pandey. Fast and numerically robust algorithm for computing the
H∞ optimum. In Decision and Control, 1991., Proceedings of the 30th IEEE Conference
on, pages 200–205. IEEE, 1991.
[11] K. Glover. All optimal Hankel-norm approximations of linear multivariable systems and
their L∞ error bounds. International Journal of Control, 39(6):1115–1193, 1984.
[12] G. Hewer. Existence theorems for positive semidefinite and sign indefinite stabilizing
solutions of H∞ riccati equations. SIAM journal on control and optimization, 31(1):16–
29, 1993.

38
George Kuster Bibliography 39

[13] P.P. Khargonekar, I.R. Petersen, and M.A. Rotea. H∞ -optimal control with state-
feedback. Automatic Control, IEEE Transactions on, 33(8):786–788, 1988.

[14] P.P. Khargonekar, I.R. Petersen, and K. Zhou. Robust stabilization of uncertain linear
systems: quadratic stabilizability and H∞ control theory. Automatic Control, IEEE
Transactions on, 35(3):356–361, 1990.

[15] D.J.N. Limebeer and M. Green. Linear robust control. Prentice-Hall, Englewood Clis,
NJ, 1995.

[16] W.W. Lin, C.S. Wang, and Q.F. Xu. On the computation of the optimal H∞ norms
for two feedback control problems. Linear algebra and its applications, 287(1):223–255,
1999.

[17] K.A. Morris. Introduction to feedback control. Academic Press, Inc., 2000.

[18] A. Packard, M.K.H. Fan, and J. Doyle. A power method for the structured singular
value. In Decision and Control, 1988., Proceedings of the 27th IEEE Conference on,
pages 2132–2137. IEEE, 1988.

[19] P. Pandey, C. Kenney, A. Packard, and AJ Laub. A gradient method for computing the
optimal H∞ -norm. Automatic Control, IEEE Transactions on, 36(7):887–890, 1991.

[20] I. Petersen. Disturbance attenuation and H∞ optimization: a design method based on

the algebraic riccati equation. Automatic Control, IEEE Transactions on, 32(5):427–
429, 1987.

[21] C. Scherer. H∞ -control by state-feedback and fast algorithms for the computation of
optimalH∞ -norms. Automatic Control, IEEE Transactions on, 35(10):1090–1099, 1990.

[22] D. Tall and S. Vinner. Concept image and concept definition in mathematics with partic-
ular reference to limits and continuity. Educational studies in mathematics, 12(2):151–
169, 1981.

[23] G. Zames. Feedback and optimal sensitivity: Model reference transformations, multi-
plicative seminorms, and approximate inverses. Automatic Control, IEEE Transactions
on, 26(2):301–320, 1981.

[24] M. Zandieh. A theoretical framework for analyzing student understanding of the con-
cept of derivative. Research in Collegiate Mathematics Education. IV. CBMS Issues in
Mathematics Education, pages 103–127, 2000.

[25] K. Zhou, J.C. Doyle, K. Glover, et al. Robust and optimal control, volume 40. Prentice
Hall Upper Saddle River, NJ, 1996.