ECE 602 Fall 2010

State Equation Solution

by

Stanislaw H. Żak

September 14, 2010

Linear systems can be modeled using convolutions and in the case of lumped systems, we can use state-space equations. Analytical computation of convolutions can be quite involved. An efficient method of obtaining a solution to a model of a lumped system is to convert the transfer function into a state-space model and then solve the state equation. We begin with a time-invariant case and then proceed with the time-varying case.

1 State Equation Solution of Uncontrolled System

Given a dynamical system’s time-invariant model of the form

ẋ(t) = Ax(t)  (1)

subject to an initial condition

x(0) = x_{0},  (2)

where A ∈ ℝ^{n×n}. The solution to (1) subject to (2) is

x(t) = e^{At}x_{0}  (3)

The above is a generalization of the situation when A is a scalar. The exponential matrix e^{At} is defined as

e^{At} = I_{n} + tA + \frac{t^{2}}{2!}A^{2} + \frac{t^{3}}{3!}A^{3} + ⋯

To verify that (3) solves (1) subject to (2), we use the fact that

\frac{d}{dt}e^{At} = Ae^{At} = e^{At}A.

Suppose now that the initial condition x(0) = x_{0} is replaced with a more general condition

x(t_{0}) = x_{0},

where t_{0} is a given initial time instant. Then, the solution to ẋ(t) = Ax(t) subject to x(t_{0}) = x_{0} is

x(t) = e^{A^{(t−t_{0})}}x_{0}

The matrix e^{A^{(t−t_{0})}} is often written as

e^{A^{(t−t_{0})}} = Φ (t, t_{0})

and is called the state transition matrix, because it relates the state at any instant of time t_{0} to the state at any other time t.

2 Computing State Transition Matrix of LTI Systems

For an overview of methods to calculate the state transition matrix for a linear, timeinvariant (LTI) system, we refer the reader to the classic paper on the subject by Moler and Van Loan [2], where the numerical issues related to the computation of the matrix exponential are discussed. In what follows, we present a few simple methods to compute the matrix exponential.

2.1 The Inverse Laplace Transform Method

In the case when the matrix A is of low dimensions, we can use the formula that results from applying the Laplace transform to the time-invariant system ẋ = Ax, x(0) = x_{0}. The Laplace transform of the above matrix differential equation is

sX(s) − x_{0} = AX(s).

Hence

X(s) = (sI_{n} − A)^{−1} x_{0},

and the inverse transform yields

e^{At} = L^{−1} {(sI_{n} − A)^{−1}}  (4)

where L^{−1} denotes the inverse Laplace transform operator.

Example 1 We use the inverse Laplace method to compute e^{At} for

A = [ \begin{matrix}0 & 3 \\ 2 & 1 \end{matrix} ] .

We have

e^{At} = L^{−1} ([sI_{2} − A]^{−1})
= L^{−1} ([ \begin{matrix}s & −3 \\ −2 & s − 1 \end{matrix} ]^{−1})
= L^{−1} (\frac{1}{(s − 3)(s + 2)} [ \begin{matrix}s − 1 & 3 \\ 2 & s \end{matrix} ])
= [ \begin{matrix}0.4e^{3t} + 0.6e^{−2t} & 0.6 (e^{3t} − e^{−2t}) \\ 0.4 (e^{3t} − e^{−2t}) & 0.6e^{3t} + 0.4e^{−2t} \end{matrix} ] .  (5)

2.2 Diagonalization Method

A convenient way to calculate e^{At} is when the matrix A has distinct eigenvalues. In such a case the n-by-n matrix A has n linearly independent eigenvectors as per the following theorem.

Theorem 1 Suppose the characteristic equation Δ(λ) = 0 has n distinct roots, λ_{1}, λ_{2}, …, λ_{n}. Then, there exist n linearly independent vectors v_{1}, v_{2}, …, v_{n} such that

Av_{i} = λ_{i}v_{i}, i = 1, 2, …, n.

Proof See Exercise 3.

□

Let

T = [ v_{1} v_{2} ⋯ v_{n} ] .

The matrix T is invertible because the eigenvectors v_{i}, i = 1, 2, …, n are linearly independent. Then, because T^{−1}T = I_{n}, we have

T^{−1}AT = T^{−1}A [ v_{1} v_{2} ⋯ v_{n} ]
= T^{−1} [ Av_{1} Av_{2} ⋯ Av_{n} ]
= T^{−1} [ λ_{1}v_{1} λ_{2}v_{2} ⋯ λ_{n}v_{n} ]

= T^{−1}T [ \begin{matrix}λ_{1} & 0 & 0 & ⋯ & 0 \\ 0 & λ_{2} & 0 & ⋯ & 0 \\ ⋮ & ⋮ & & ... & ⋮ \\ 0 & 0 & 0 & ⋯ & λ_{n} \end{matrix} ]

= Λ.  (6)

We use the above result to compute the transition matrix e^{At}. Recall that

e^{At} = I_{n} + At + \frac{A^{2}t^{2}}{2!} + \frac{A^{3}t^{3}}{3!} + ⋯

Note that

Λ^{k} = [ \begin{matrix}λ_{1}^{k} & 0 & 0 & ⋯ & 0 \\ 0 & λ_{2}^{k} & 0 & ⋯ & 0 \\ ⋮ & ⋮ & & ... & ⋮ \\ 0 & 0 & 0 & ⋯ & λ_{n}^{k} \end{matrix} ] , k ≥ 1

Hence

e^{Λt} = [ \begin{matrix}e^{λ_{1}t} & 0 & 0 & ⋯ & 0 \\ 0 & e^{λ_{2}t} & 0 & ⋯ & 0 \\ ⋮ & ⋮ & & ... & ⋮ \\ 0 & 0 & 0 & ⋯ & e^{λ_{n}t} \end{matrix} ] .

From (6), we obtain

A = T ΛT^{−1}.

Therefore,

A^{k} = (T ΛT^{−1})^{k}
= (T ΛT^{−1}) (T ΛT^{−1}) × ⋯ × (T ΛT^{−1})
= T Λ^{k}T^{−1},

and hence

e^{At} = e^{T ΛT^{−1}}
= I_{n} + T ΛT^{−1}t + \frac{(T ΛT^{−1})^{2} t^{2}}{2!} + \frac{(T ΛT^{−1})^{3} t^{3}}{3!} + ⋯
= T (I_{n} + Λt + \frac{Λ^{2}t^{2}}{2!} + \frac{Λ^{3}t^{3}}{3!} ) T^{−1}
= T e^{Λt}T^{−1}.  (7)

Example 2 We use the diagonalization method to compute e^{At} for

A = [ \begin{matrix}0 & 3 \\ 2 & 1 \end{matrix} ] .

Note that this is the same matrix A as in Example 1. The eigenvalues of A are distinct so the diagonalization method is applicable here. We first compute an eigenvector, v_{1}, corresponding to the eigenvalue λ_{1} = 3 by solving the eigenvalue-eigenvector equation,

(A − λ_{1}I_{2}) v_{1} = 0.

That is,

(A − λ_{1}I_{2}) v_{1} = [ \begin{matrix}−3 & 3 \\ 2 & −2 \end{matrix} ] [ \begin{matrix}v_{11} \\ v_{21} \end{matrix} ] = [ \begin{matrix}0 \\ 0 \end{matrix} ] .

We take

v_{1} = [ \begin{matrix}1 \\ 1 \end{matrix} ] .

We next compute an eigenvector corresponding to the eigenvalue λ_{2} = −2 by solving the equation

(A − λ_{2}I_{2}) v_{2} = [ \begin{matrix}2 & 3 \\ 2 & 3 \end{matrix} ] [ \begin{matrix}v_{12} \\ v_{22} \end{matrix} ] = [ \begin{matrix}0 \\ 0 \end{matrix} ] .

We take

v_{2} = [ \begin{matrix}−3 \\ 2 \end{matrix} ] .

We use the above eigenvectors to form the transformation matrix

T = [ v_{1} v_{2} ] = [ \begin{matrix}1 & −3 \\ 1 & 2 \end{matrix} ] .

Its inverse is

T^{−1} = \frac{1}{5} [ \begin{matrix}2 & 3 \\ −1 & 1 \end{matrix} ] .

We are now ready to compute e^{At} using (7),

e^{At} = T e^{Λt}T^{−1}
= \frac{1}{5} [ \begin{matrix}1 & −3 \\ 1 & 2 \end{matrix} ] [ \begin{matrix}e^{3t} & 0 \\ 0 & e^{−2t} \end{matrix} ] [ \begin{matrix}2 & 3 \\ −1 & 1 \end{matrix} ]
= \frac{1}{5} [ \begin{matrix}2e^{3t} + 3e^{−2t} & 3 (e^{3t} − e^{−2t}) \\ 2 (e^{3t} − e^{−2t}) & 3e^{3t} + 2e^{−2t} \end{matrix} ] .

The result is the same as when we used the inverse Laplace transform method.

The above example illustrates that the transformation T producing the diagonal form of A is a powerful computational tool. The existence of a diagonalizing T is not always guaranteed. The inability to transform A to diagonal form is connected with the existence of repeated eigenvalues that we discuss next.

3 Computing Functions of a Square Matrix

In this section, we discuss methods of computing functions of a square matrix that is not necessarily diagonalizable. The methods discussed allow us to easily compute the matrix exponential and hence the state transition matrix of LTI systems. A major tool used by the methods discussed is the so called Jordan canonical form that we discuss next.

3.1 The Jordan Canonical Form

If a square n-by-n matrix A has repeated eigenvalues, then it may not be diagonalizable. Although some matrices possessing repeated eigenvalues can be transformed into diagonal form, in many cases it is only possible to obtain the so-called “Jordan canonical form” or just the “Jordan from,” for short. To illustrate the point, we consider the case when the characteristic equation, Δ(λ) = (λ−λ_{1})^{n}, that is, the matrix A has just one eigenvalue with multiplicity n. For simplicity, we analyze the case when n = 4 and suppose that

rank [ A − λ_{1}I_{4} ] = 3.

That is, the nullity of the matrix [ A − λ_{1}I_{4} ] is one, which means that there is only one eigenvector corresponding to the eigenvalue λ_{1}. Thus such a matrix A cannot be diagonalized. We need three more linearly independent vectors to represent A in a new basis, where it has a particularly simple form that is as close as possible to being diagonal.

To proceed, we need a notion of a generalized eigenvector that we define now.

Definition 1 A vector v is called a generalized eigenvector of grade k if

(A − λI_{n})^{k} v = 0

and

(A − λI_{n})^{k−1} v ≠ 0.
◇

For k = 1, we have

(A − λI_{n}) v = 0,

which means that v is an ordinary eigenvector corresponding to the eigenvalue λ.

Suppose now that A ∈ ℝ^{4×4} and v is a generalized eigenvector of grade n = 4 corresponding to the repeated eigenvalue λ_{1} whose multiplicity is also 4. Then by Definition 1,

(A − λ_{1}I_{4})^{4} v = 0 but (A − λ_{1}I_{4})^{3} v ≠ 0.

Let

v_{1} = (A − λ_{1}I_{4})^{3} v.

Then,

(A − λ_{1}I_{4}) v_{1} = (A − λ_{1}I_{4})^{4} v
= 0.

Hence

Av_{1} = λ_{1}v_{1}.  (8)

We also have to have

(A − λ_{1}I_{4})^{2} v ≠ 0,

otherwise we would have (A − λ_{1}I_{4})^{3} v = 0 as well, which would mean that v is not a generalized eigenvector of grade 4.

Let

v_{2} = (A − λ_{1}I_{4})^{2} v.

Note that

(A − λ_{1}I_{4}) v_{2} = (A − λ_{1}I_{4})^{3} v
= v_{1}.

Therefore,

Av_{2} = v_{1} + λ_{1}v_{2}.  (9)

We proceed in a similar way by letting

v_{3} = (A − λ_{1}I_{4}) v.

We observer that

(A − λ_{1}I_{4}) v_{3} = (A − λ_{1}I_{4})^{2} v
= v_{2}.

Hence,

Av_{3} = v_{2} + λ_{1}v_{3}.  (10)

Finally, we let

v_{4} = v

and we note that

(A − λ_{1}I_{4}) v_{4} = (A − λ_{1}I_{4}) v
= v_{3}.

Therefore,

Av_{4} = v_{3} + λ_{1}v_{4}.  (11)

Combining now (8)–(11), we obtain

A [ v_{1} v_{2} v_{3} v_{4} ] = [ λ_{1}v_{1} v_{1} + λ_{1}v_{2} v_{2} + λ_{1}v_{3} v_{3} + λ_{1}v_{4} ]

= [ v_{1} v_{2} v_{3} v_{4} ] [ \begin{matrix}λ_{1} & 1 & 0 & 0 \\ 0 & λ_{1} & 1 & 0 \\ 0 & 0 & λ_{1} & 1 \\ 0 & 0 & 0 & λ_{1} \end{matrix} ] .  (12)

Let

T = [ v_{1} v_{2} v_{3} v_{4} ] .

The matrix T is invertible by construction, that is, the set of vectors {v_{1}, v_{2}, v_{3}, v_{4}} forms a basis for ℝ^{4}. Therefore we can represent (12) as

T^{−1}AT = [ \begin{matrix}λ_{1} & 1 & 0 & 0 \\ 0 & λ_{1} & 1 & 0 \\ 0 & 0 & λ_{1} & 1 \\ 0 & 0 & 0 & λ_{1} \end{matrix} ]

= J.  (13)

The representation of A in the basis {v_{1}, v_{2}, v_{3}, v_{4}} has the form of a Jordan block of order 4.

Example 3 Consider the matrix

A = [ \begin{matrix}5 & 2 & 4 \\ 0 & 5 & 3 \\ 0 & 0 & 5 \end{matrix} ] .

The characteristic polynomial is Δ(λ) = (λ − 5)^{3} and

rank [ A − λI_{3} ]|_{λ=5} = 2,

that is, the nullity of the characteristic matrix is 1. We proceed to find v_{1} such that

[A − λI_{3}]|_{λ=5} v_{1} = [ \begin{matrix}0 & 2 & 4 \\ 0 & 0 & 3 \\ 0 & 0 & 0 \end{matrix} ] [ \begin{matrix}v_{11} \\ v_{21} \\ v_{31} \end{matrix} ]
= 0.

We take

v_{1} = [ \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ] .

Then we proceed to find v_{2} such that

[A − λI_{3}]|_{λ=5} v_{2} = v_{1},

that is, we solve for v_{2} satisfying

[ \begin{matrix}0 & 2 & 4 \\ 0 & 0 & 3 \\ 0 & 0 & 0 \end{matrix} ] [ \begin{matrix}v_{12} \\ v_{22} \\ v_{32} \end{matrix} ] = [ \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ] .

We can take

v_{2} = [ \begin{matrix}0 \\ 1/2 \\ 0 \end{matrix} ] .

Finally we find v_{3} such that

[A − λI_{3}]|_{λ=5} v_{3} = v_{2},

that is, we solve for v_{3} satisfying

[ \begin{matrix}0 & 2 & 4 \\ 0 & 0 & 3 \\ 0 & 0 & 0 \end{matrix} ] [ \begin{matrix}v_{13} \\ v_{23} \\ v_{33} \end{matrix} ] = [ \begin{matrix}0 \\ 1/2 \\ 0 \end{matrix} ] .

We can take

v_{3} = [ \begin{matrix}0 \\ −2/6 \\ 1/6 \end{matrix} ] .

The transformation matrix T has the form

T = [ v_{1} v_{2} v_{3} ] = [ \begin{matrix}1 & 0 & 0 \\ 0 & 1/2 & −2/6 \\ 0 & 0 & 1/6 \end{matrix} ] ,

and

T^{−1} = [ \begin{matrix}1 & 0 & 0 \\ 0 & 2 & 4 \\ 0 & 0 & 6 \end{matrix} ] .

The matrix A in the new basis has the form

T^{−1}AT = [ \begin{matrix}1 & 0 & 0 \\ 0 & 2 & 4 \\ 0 & 0 & 6 \end{matrix} ] [ \begin{matrix}5 & 2 & 4 \\ 0 & 5 & 3 \\ 0 & 0 & 5 \end{matrix} ] [ \begin{matrix}1 & 0 & 0 \\ 0 & 1/2 & −2/6 \\ 0 & 0 & 1/6 \end{matrix} ]

= [ \begin{matrix}5 & 1 & 0 \\ 0 & 5 & 1 \\ 0 & 0 & 5 \end{matrix} ] ,

which is a third order Jordan block.

If the rank of the characteristic matrix evaluated at λ_{1} is n − 2, then the equation

[A − λ_{1}I_{n}] v = 0,

has two linearly independent solutions. Therefore, A has two linearly independent eigenvectors. Let u_{1} and v_{1} be linearly independent eigenvectors associated with the eigenvalue λ_{i}. In order to transform A into its Jordan form, we need to find n − 2 generalized eigenvectors. In this case we have two chains of generalized eigenvectors, {u_{1}, u_{2}, …, u_{k}} and {v_{1}, u_{2}, …, v_{l}}, where k + l = n. In the basis

{u_{1}, u_{2}, …, u_{k}, v_{1}, u_{2}, …, v_{l}}

the matrix A will have a block diagonal form

T^{−1}AT = [ \begin{matrix}J_{1} & O \\ O & J_{2} \end{matrix} ] ,

where J_{1} is a Jordan block of order k and J_{2} is a Jordan block of order l.

Example 4 Let

A = [ \begin{matrix}3 & 1 & −1 \\ 0 & 6 & −3 \\ 0 & 3 & 0 \end{matrix} ] .

The characteristic polynomial

Δ(λ) = (λ − 3)^{3}

and

rank [A − λI_{3}]|_{λ=3} = [ \begin{matrix}0 & 1 & −1 \\ 0 & 3 & −3 \\ 0 & 3 & −3 \end{matrix} ]
= 1

As two linearly independent solutions to

[A − λI_{3}]|_{λ=3} v = 0,

we can take

u_{1} = [ \begin{matrix}1 \\ 3 \\ 3 \end{matrix} ] and v_{1} = [ \begin{matrix}0 \\ 1 \\ 1 \end{matrix} ] .

Now we have to find one generalized eigenvector that together with u_{1} and v_{1} would form a basis for ℝ^{3} and enable us to construct a transformation that brings A into its Jordan form. We try first u_{1} to construct u_{2} such that

[A − λI_{3}]|_{λ=3} u_{2} = u_{1},

that is, we solve for u_{2} satisfying

[ \begin{matrix}0 & 1 & −1 \\ 0 & 3 & −3 \\ 0 & 3 & −3 \end{matrix} ] [ \begin{matrix}u_{12} \\ u_{22} \\ u_{32} \end{matrix} ] = [ \begin{matrix}1 \\ 3 \\ 3 \end{matrix} ] .

We can take

u_{2} = [ \begin{matrix}1 \\ 1 \\ 0 \end{matrix} ] .

We form the transformation matrix

T = [ u_{1} u_{2} v_{1} ] = [ \begin{matrix}1 & 1 & 0 \\ 3 & 1 & 1 \\ 3 & 0 & 1 \end{matrix} ] .

Its inverse is

T^{−1} = [ \begin{matrix}1 & −1 & 1 \\ 0 & 1 & −1 \\ −3 & 3 & −2 \end{matrix} ] .

It is then easy to check that

T^{−1}AT = [ \begin{matrix}J_{1} & O \\ O & J_{2} \end{matrix} ] = [ \begin{matrix}3 & 1 & 0 \\ 0 & 3 & 0 \\ 0 & 0 & 3 \end{matrix} ] ,

where

J_{1} = [ \begin{matrix}3 & 1 \\ 0 & 3 \end{matrix} ] and J_{2} = 3.

In general, the Jordan canonical form is a block diagonal matrix, where each diagonal block is a Jordan block corresponding to a possibly different eigenvalue, that is,

T^{−1}AT = [ \begin{matrix}λ_{1} & 1 & 0 & ⋯ & 0 \\ 0 & λ_{1} & 1 & ⋯ & 0 \\ 0 & 0 & λ_{1} & ⋯ & 0 \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & 0 & ⋯ & λ_{1} \\ & & & & & . & . & . \\ & & & & & & & & λ_{r} & 1 & 0 & ⋯ & 0 \\ & & & & & & & & 0 & λ_{r} & 1 & ⋯ & 0 \\ & & & & & & & & 0 & 0 & λ_{r} & ⋯ & 0 \\ & & & & & & & & ⋮ & ⋮ & ⋮ & ... & ⋮ \\ & & & & & & & & 0 & 0 & 0 & ⋯ & λ_{r} \end{matrix} ]

= [ \begin{matrix}J_{1} & O & ⋯ & O \\ O & J_{2} & & O \\ ⋮ & & ... & ⋮ \\ O & O & ⋯ & J_{q} \end{matrix} ] ,  (14)

where each r_{i}×r_{i} Jordan block has the form J_{i} = λ_{i} if r_{i} = 1, or, for r_{i} > 1 the corresponding Jordan block has a non-diagonal structure. Note that

n = r_{1} + r_{2} + ⋯ + r_{q}

and q ≥ r, which means that there may be more than one Jordan blocks associated with one eigenvalue as was the case in Example 4.

Example 5 To illustrate the general case depicted in (14), we consider, as in Chen [1, p. 60], a 5-by-5 matrix A with repeated eigenvalue λ_{1} of multiplicity 4 and a simple eigenvalue λ_{2}, that is, the characteristic polynomial of A is

Δ(λ) = (λ − λ_{1})^{4}(λ − λ_{2}).

Depending on the nullity of the matrix [A − λ_{1}I_{5}], the Jordan form of the matrix A make take one of the following forms:

• If the nullity of [A − λ_{1}I_{5}] is 1, then there is only one eigenvector, say u_{1} associated with λ_{1} and we need to compute three generalized eigenvectors to form a basis {u_{1}, u_{2}, u_{3}, u_{4}, v_{1}}, where v_{1} is an eigenvector associated with λ_{2}. The matrix A in this basis will have the form:

Â_{1} = T _{1}^{−1}AT_{1} = [ \begin{matrix}λ_{1} & 1 & 0 & 0 & ⋮ & 0 \\ 0 & λ_{1} & 1 & 0 & ⋮ & 0 \\ 0 & 0 & λ_{1} & 1 & ⋮ & 0 \\ 0 & 0 & 0 & λ_{1} & ⋮ & 0 \\ ⋯ & ⋯ & ⋯ & ⋯ & ⋯ & ⋯ \\ 0 & 0 & 0 & 0 & ⋮ & λ_{2} \end{matrix} ] .

Note that we can re-order the basis vectors as {v_{1}, u_{1}, u_{1}, u_{3}, u_{4}}. In this basis the matrix A representation has the form

Â_{2} = T _{2}^{−1}AT_{2} = [ \begin{matrix}λ_{2} & ⋮ & 0 & 0 & 0 & 0 \\ ⋯ & ⋯ & ⋯ & ⋯ & ⋯ & ⋯ \\ 0 & ⋮ & λ_{1} & 1 & 0 & 0 \\ 0 & ⋮ & 0 & λ_{1} & 1 & 0 \\ 0 & ⋮ & 0 & 0 & λ_{1} & 1 \\ 0 & ⋮ & 0 & 0 & 0 & λ_{1} \end{matrix} ] .

• If the nullity of [A − λ_{1}I_{5}] is 2, then there are two eigenvectors associated with λ_{1}, that is, there are two Jordan blocks associated with λ_{1}. The matrix A may take one of the following forms:

Â_{3} = T _{3}^{−1}AT_{3} = [ \begin{matrix}λ_{1} & 1 & 0 & ⋮ & 0 & & 0 \\ 0 & λ_{1} & 1 & ⋮ & 0 & & 0 \\ 0 & 0 & λ_{1} & ⋮ & 0 & & 0 \\ ⋯ & ⋯ & ⋯ & ⋯ & ⋯ & ⋯ \\ 0 & 0 & 0 & ⋮ & λ_{1} & ⋮ & 0 \\ & & & ⋯ & ⋯ & ⋯ & ⋯ \\ 0 & 0 & 0 & & 0 & ⋮ & λ_{2} \end{matrix} ]

or

Â_{4} = T _{4}^{−1}AT_{4} = [ \begin{matrix}λ_{1} & 1 & ⋮ & 0 & 0 & & 0 \\ 0 & λ_{1} & ⋮ & 0 & 0 & & 0 \\ ⋯ & ⋯ & ⋯ & ⋯ & ⋯ & ⋯ \\ 0 & 0 & ⋮ & λ_{1} & 1 & ⋮ & 0 \\ 0 & 0 & ⋮ & 0 & λ_{1} & ⋮ & 0 \\ & & ⋯ & ⋯ & ⋯ & ⋯ & ⋯ \\ 0 & 0 & & 0 & 0 & ⋮ & λ_{2} \end{matrix} ] .

• If the nullity of [A − λ_{1}I_{5}] is 3, then there are three eigenvectors associated with λ_{1}, that is, there are three Jordan blocks associated with λ_{1}. The matrix A may take the following form:

Â_{5} = T _{5}^{−1}AT_{5} = [ \begin{matrix}λ_{1} & 1 & ⋮ & 0 & & 0 & & 0 \\ 0 & λ_{1} & ⋮ & 0 & & 0 & & 0 \\ ⋯ & ⋯ & ⋯ & ⋯ & ⋯ \\ 0 & 0 & ⋮ & λ_{1} & ⋮ & 0 & & 0 \\ & & ⋯ & ⋯ & ⋯ & ⋯ & ⋯ \\ 0 & 0 & & 0 & ⋮ & λ_{1} & ⋮ & 0 \\ & & & & ⋯ & ⋯ & ⋯ & ⋯ \\ 0 & 0 & & 0 & & 0 & ⋮ & λ_{2} \end{matrix} ] .

• Finally, if the nullity of [A − λ_{1}I_{5}] is 4, then there are four eigenvectors associated with λ_{1}, that is, there are four Jordan blocks associated with λ_{1}. The Jordan form of the matrix matrix A is a diagonal matrix,

Â_{6} = T _{6}^{−1}AT_{6} = [ \begin{matrix}λ_{1} & 0 & 0 & 0 & 0 \\ 0 & λ_{1} & 0 & 0 & 0 \\ 0 & 0 & λ_{1} & 0 & 0 \\ 0 & 0 & 0 & λ_{1} & 0 \\ 0 & 0 & 0 & 0 & λ_{2} \end{matrix} ] .

To summarize our discussion, we need the following definitions.

Definition 2 For a given eigenvalue λ_{i} of an n-by-n matrix A, its algebraic multiplicity is the multiplicity of it in the characteristics polynomial Δ(λ), that is, the algebraic multiplicity of λ_{i} is the sum of orders of all Jordan blocks associated with λ_{i}.

◇

Definition 3 For a given eigenvalue λ_{i} of an n-by-n matrix A, its geometric multiplicity is the nullity of the associated characteristic matrix [ A − λ_{i}I_{n} ], that is, the geometric multiplicity of λ_{i} equals the number of Jordan blocks associated with it.

◇

Then we have the following theorem:

Theorem 2 A given n-by-n matrix A is diagonalizable if and only if geometric multiplicities of its eigenvalues are the same as their algebraic multiplicities, that is, if and only if A has n linearly independent eigenvectors.

We close our discussion on the Jordan canonical form with an observation by Owens [3, p. 34] that in almost all situations of practical relevance, the matrix A has distinct eigenvalues. So in a sense the cases described in this section can be considered “pathological” from the applications point of view. In the literature, see for example [1, p. 61], if A cannot be diagonalized then it is called to be defective. Calculation of the Jordan canonical is numerically ill conditioned. One can try MATLAB’s function jordan. However, the matrix A must be known exactly, which in MATLAB’s language means that its elements must be integers or ratios of small integers. Finally, to diagonalize a square matrix using MATLAB, one can call the function eig. If, however, A is not diagonalizable, then the function eig will give incorrect solution.

3.2 Computing Polynomial Functions of a Square Matrix

If A is an n-by-n square matrix then its k-th power is A^{k} = \underbrace{AA × ⋯ × A}_{k times},

where

A^{0} = I_{n}.

Suppose now that we have a polynomial

f(λ) = λ^{n} + α_{n−1}λ^{n−1} + ⋯ + α_{1}λ + α_{0}.

Then we define f(A) as

f(A) = A^{n} + α_{n−1}A^{n−1} + ⋯ + α_{1}A + α_{0}I_{n}.

For example, if

f(λ) = λ^{2} + 2λ − 3 = (λ + 3)(λ − 1),

then

f(A) = A^{2} + 2A − 3I_{n} = (A + 3I_{n})(A − I_{n}).

Now if the matrix A is block diagonal, for example, such as

A = [ \begin{matrix}A_{1} & O & O \\ O & A_{2} & O \\ O & O & A_{3} \end{matrix} ] ,

where A_{i} are square matrices and O denotes a matrix of zeros of appropriate dimensions, then we have

A^{k} = [ \begin{matrix}A_{1}^{k} & O & O \\ O & A_{2}^{k} & O \\ O & O & A_{3}^{k} \end{matrix} ]

and

f(A) = [ \begin{matrix}f(A_{1}) & O & O \\ O & f(A_{2}) & O \\ O & O & f(A_{3}) \end{matrix} ] .

Consider a similar matrix Â of the given matrix A, that is,

Â = T^{−1}AT ,

where T ∈ ℝ^{n×n} is an invertible matrix. The above can be re-written as

A = T ÂT^{−1},

and we have

A^{k} = (T ÂT^{−1}) (T ÂT^{−1}) × ⋯ × (T ÂT^{−1})
= T Â^{k}T^{−1},

or

Â^{k} = (T^{−1}AT ) (T^{−1}AT ) × ⋯ × (T^{−1}AT )
= T^{−1}A^{k}T .

Thus, we have

f (A) = T f(Â)T^{−1}  (15)

or

f(Â) = T^{−1}f (A) T .  (16)

Now let λ_{i} be the i-th eigenvalue of A and suppose that the algebraic multiplicity of λ_{i} is n_{i}. Then the characteristic polynomial of A can be written as,

Δ(λ) = det(λI_{n} − A) = ∏_{i=1}^{r} (λ − λ_{i})^{n_{i}} .  (17)

We now discuss a very useful property of Jordan blocks that we will utilize in the proof of the very imporatnt theorem of linear algebra, the theorem of Cayley and Hamilton. Suppose that we have a fourth order Jordan block,

J = [ \begin{matrix}λ_{i} & 1 & 0 & 0 \\ 0 & λ_{i} & 1 & 0 \\ 0 & 0 & λ_{i} & 1 \\ 0 & 0 & 0 & λ_{i} \end{matrix} ] .

Then we have

(J − λ_{i}I) = [ \begin{matrix}0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \end{matrix} ] .

Let us now compute successive powers of (J − λI). We obtain

(J − λ_{i}I)^{2} = [ \begin{matrix}0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{matrix} ] , (J − λ_{i}I)^{3} = [ \begin{matrix}0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{matrix} ] ,

and

(J − λ_{i}I)^{k} = O for k ≥ 4.

A matrix with above property is called a nilpotent matrix.

In general, if we have a Jordan block of order r_{i} then

(J − λ_{i}I)^{k} = O for k ≥ r_{i}.  (18)

We have all we need to state and prove the theorem of Cayley and Hamilton.

Theorem 3 (Cayley-Hamilton theorem) A square matrix satisfies its characteristic polynomial equation, that is,

Δ(A) = O.

Proof Suppose that Â is the Jordan canonical form of A. Thus there exists a non-singular matrix T such that

Â = T^{−1}AT = [ \begin{matrix}J_{1} & O & ⋯ & O \\ O & J_{2} & & O \\ ⋮ & & ... & ⋮ \\ O & O & ⋯ & J_{q} \end{matrix} ] .

Let us now consider the matrix polynomial Δ(A). We use (15) to obtain

Δ(A) = T Δ(Â)T^{−1}
= T ∏_{i=1}^{r} (Â − λ_{i}I_{n})^{n_{i}} T^{−1}

= T ∏_{i=1}^{r} [ \begin{matrix}J_{1} − λ_{i}I_{r_{1}} & O & ⋯ & O \\ O & J_{2} − λ_{i}I_{r_{2}} & & O \\ ⋮ & & ... & ⋮ \\ O & O & ⋯ & J_{q} − λ_{i}I_{r_{q}} \end{matrix} ]^{n_{i}} T^{−1}

= T ∏_{i=1}^{r} [ \begin{matrix}(J_{1} − λ_{i}I_{r_{1}})^{n_{i}} & O & ⋯ & O \\ O & (J_{2} − λ_{i}I_{r_{2}})^{n_{i}} & & O \\ ⋮ & & ... & ⋮ \\ O & O & ⋯ & (J_{q} − λ_{i}I_{r_{q}})^{n_{i}} \end{matrix} ] T^{−1}.

The order of any Jordan block associated with λ_{i} cannot exceed the corresponding algebraic multiplicity. Taking into account the nilpotent property given by (18), we obtain

Δ(A) = T [ \begin{matrix}O & O & ⋯ & O \\ O & (J_{2} − λ_{1}I_{r_{2}})^{n_{1}} & & O \\ ⋮ & & ... & ⋮ \\ O & O & ⋯ & (J_{q} − λ_{1}I_{r_{q}})^{n_{1}} \end{matrix} ] ×

⋯ × [ \begin{matrix}(J_{1} − λ_{r}I_{r_{1}})^{n_{r}} & O & ⋯ & O \\ O & (J_{2} − λ_{r}I_{r_{2}})^{n_{r}} & & O \\ ⋮ & & ... & ⋮ \\ O & O & ⋯ & O \end{matrix} ] T^{−1}

= T [ \begin{matrix}O & O & ⋯ & O \\ O & O & & O \\ ⋮ & & ... & ⋮ \\ O & O & ⋯ & O \end{matrix} ] T^{−1}

= O,  (19)

which completes the proof of the theorem.

□

It follows from the Cayley-Hamilton theorem that Δ(A) = A^{n} + α_{n−1}A^{n−1} + ⋯ + α_{1}A + α_{0}I_{n} = O, that is,

A^{n} = −α_{n−1}A^{n−1} − ⋯ − α_{1}A − α_{0}I_{n}.  (20)

Thus A^{n} can be expressed as a linear combination of the matrices

{I_{n}, A, …, A^{n−1}}.

If we multiply both sides of (20) by A, then we obtain

A^{n+1} = −α_{n−1}A^{n} − ⋯ − α_{1}A^{2} − α_{0}A.  (21)

Thus, A^{n+1} can be expressed as a linear combination of

{A, A^{2}, …, A^{n}}.

Using (20), we conclude that A^{n+1} can be expressed as a linear combination of {I_{n}, A, …, A^{n−1}}. Using the same argument, we conclude that any A^{k}, k ≥ 0, can be expressed as a linear combination of {I_{n}, A, …, A^{n−1}}. Therefore, any matrix polynomial f(A) can also be expressed as

f(A) = β_{0}I_{n} + β_{1}A + ⋯ + β_{n−1}A^{n−1}

for some β_{i}, i = 0, 1, …, n − 1. The above property is used to efficiently compute matrix polynomials of A, which can then be used to compute the matrix exponential e^{At} as well as the state transition matrix for discrete-time systems. In what follows, we present such a method and illustrate it with examples. We follow the presentation of Chen [1, pp. 63–67].

Suppose we are given a polynomial f(λ) and we are interested in calculating f(A). We begin with dividing f(λ) by Δ(λ) to obtain

f(λ) = q(λ)Δ(λ) + h(λ),  (22)

where q(λ) is the quotient, h(λ) is the remainder, and

deg_{λ} h(λ) < n.

We now evaluate f(A) using (22) and the Cayley-Hamilton theorem to obtain

f(A) = q(A)Δ(A) + h(A)
= q(A)O + h(A)
= h(A).

Thus instead of evaluating the matrix polynomial f(A) that may contain very high powers of A, we can find f(A) by computing h(A), where the highest power of A is less than n. In addition, rather than performing the division (22), we can find h(λ) by assuming that it has the form

h(λ) = β_{0} + β_{1}λ + ⋯ + β_{n−1}λ^{n−1}

and then solve n equations for β_{i}, i = 0, 1, …, n − 1. Consider first the case when A has distinct eigenvalues. To obtain equations involving the coefficients β_{i}, we note that

f(λ_{i}) = q(λ_{i})Δ(λ_{i}) + h(λ_{i}) = h(λ_{i}), i = 1, 2, …, n.

Example 6 For a given

f(λ) = λ^{7},

we compute

f(A) = A^{7},

where

A = [ \begin{matrix}0 & 2 \\ 1 & 1 \end{matrix} ] .  (23)

The characteristic polynomial of A is

Δ(λ) = λ^{2} − λ − 2 = (λ + 1)(λ − 2).

Let

h(λ) = β_{0} + β_{1}λ.

Then we have two equations:

f(λ_{i}) = h(λ_{i}), i = 1, 2,

that is,

f(−1) = (−1)^{7} = −1 = β_{0} − β_{1}

and

f(2) = 2^{7} = 128 = β_{0} + 2β_{1}.

We represent the above linear equations for β_{0} and β_{1} in the matrix format as

[ \begin{matrix}1 & −1 \\ 1 & 2 \end{matrix} ] [ \begin{matrix}β_{0} \\ β_{1} \end{matrix} ] = [ \begin{matrix}−1 \\ 128 \end{matrix} ] .

Solving, we obtain β_{0} = 42 and β_{1} = 43. Hence

f(A) = A^{7} = 42I_{2} + 43A = [ \begin{matrix}42 & 86 \\ 43 & 85 \end{matrix} ] .

Example 7 We can use the above method to evaluate e^{At}, where A is given by (23). In other words, let f(λ) = e^{λt}; we wish to find

f(A) = e^{At}.

Let h(λ) = h(λ) = β_{0} + β_{1}λ. Then

f(−1) = h(−1)

yields

e^{−t} = β_{0} − β_{1},

and

f(2) = h(2)

yields

e^{2t} = β_{0} + 2β_{1}.

We thus have two equations for β_{0} and β_{1},

[ \begin{matrix}1 & −1 \\ 1 & 2 \end{matrix} ] [ \begin{matrix}β_{0} \\ β_{1} \end{matrix} ] = [ \begin{matrix}e^{−t} \\ e^{2t} \end{matrix} ] .

Solving the above gives

[ \begin{matrix}β_{0} \\ β_{1} \end{matrix} ] = \frac{1}{3} [ \begin{matrix}2 & 1 \\ −1 & 1 \end{matrix} ] [ \begin{matrix}e^{−t} \\ e^{2t} \end{matrix} ] = \frac{1}{3} [ \begin{matrix}2e^{−t} + e^{2t} \\ −e^{−t} + e^{2t} \end{matrix} ] .

Hence

e^{At} = β_{0}I_{2} + β_{1}A = \frac{1}{3} [ \begin{matrix}2e^{−t} + e^{2t} & −2e^{−t} + 2e^{2t} \\ −e^{−t} + e^{2t} & e^{−t} + 2e^{2t} \end{matrix} ] .

When the characteristic polynomial has repeated zeros, we use the derivatives of f(λ) and h(λ) to obtain additional equations to determine the coefficients of h(λ) as per the following theorem.

Theorem 4 Chen [1, p. 64] We are given f(λ) and an n × n matrix A with characteristic polynomial

Δ(λ) = ∏_{i=1}^{m}(λ − λ_{i})^{n_{i}},

where n = ∑_{i=1}^{m} n_{i}. Let

h(λ) = β_{0} + β_{1}λ + ⋯ + β_{n−1}λ^{n−1}.

It is a polynomial of degree n − 1 with n unknown coefficients. These n unknowns are to be solved from the following set of n equations:

f^{(l)}(λ_{i}) = h^{(l)}(λ_{i}) for l = 0, 1, …, n_{i} − 1 and i = 1, 2, …, m

where

f^{(l)}(λ_{i}) = \frac{d^{l}f(λ)}{dλ^{l}} |_{λ=λ_{i}}

and h^{(l)}(λ_{i}) is similarly defined. Then we have

f(A) = h(A)

and h(λ) is said to be equal f(λ) on the spectrum of A.

Example 8 We compute e^{At} for

A = [ \begin{matrix}3 & 0 & 1 \\ −1 & 1 & 2 \\ −2 & 0 & 0 \end{matrix} ] .

The characteristic polynomial is

Δ(λ) = (λ − 1)^{2}(λ − 2).

Let f(λ) = e^{λt} and let h(λ) = β_{0} + β_{1}λ + β_{2}λ^{2}. Then using Theorem 4, we obtain

f(1) = h(1), e^{t} = β_{0} + β_{1} + β_{2},
f′(1) = h′(1), te^{t} = β_{1} + 2β_{2},
f(2) = h(2), e^{2t} = β_{0} + 2β_{1} + 4β_{2}.

Thus we have a set of equations

[ \begin{matrix}1 & 1 & 1 \\ 0 & 1 & 2 \\ 1 & 2 & 4 \end{matrix} ] [ \begin{matrix}β_{0} \\ β_{1} \\ β_{2} \end{matrix} ] = [ \begin{matrix}e^{t} \\ te^{t} \\ e^{2t} \end{matrix} ] .

Solving the above gives

[ \begin{matrix}β_{0} \\ β_{1} \\ β_{2} \end{matrix} ] = [ \begin{matrix}0 & −2 & 1 \\ 2 & 3 & −2 \\ −1 & −1 & 1 \end{matrix} ] [ \begin{matrix}e^{t} \\ te^{t} \\ e^{2t} \end{matrix} ] = [ \begin{matrix}−2te^{t} + e^{2t} \\ 2e^{t} + 3te^{t} − 2e^{2t} \\ −e^{t} − te^{t} + e^{2t} \end{matrix} ] .

Hence

e^{At} = h(A)
= β_{0}I_{3} + β_{1}A + β_{2}A^{2}

= (−2te^{t} + e^{2t}) [ \begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix} ] + (2e^{t} + 3te^{t} − 2e^{2t}) [ \begin{matrix}3 & 0 & 1 \\ −1 & 1 & 2 \\ −2 & 0 & 0 \end{matrix} ]

+ (−e^{t} − te^{t} + e^{2t}) [ \begin{matrix}7 & 0 & 3 \\ −8 & 1 & 1 \\ −6 & 0 & −2 \end{matrix} ]

= [ \begin{matrix}−e^{t} + 2e^{2t} & 0 & −e^{t} + e^{2t} \\ 6e^{t} + 5te^{t} − 6e^{2t} & e^{t} & 3e^{t} + 5te^{t} − 3e^{2t} \\ 2e^{t} − 2e^{2t} & 0 & 2e^{t} − e^{2t} \end{matrix} ] .

Example 9 Suppose we are given the Jordan block of order n_{i},

J = [ \begin{matrix}λ_{i} & 1 & 0 & ⋯ & 0 & 0 \\ 0 & λ_{i} & 1 & ⋯ & 0 & 0 \\ 0 & 0 & λ_{i} & ⋯ & 0 & 0 \\ ⋮ & ⋮ & & ... & & ⋮ \\ 0 & 0 & 0 & ⋯ & λ_{i} & 1 \\ 0 & 0 & 0 & ⋯ & 0 & λ_{i} \end{matrix} ] .

The characteristic polynomial of J is Δ(λ) = (λ − λ_{i})^{n_{i}}. For computational simplicity we select the (n_{i} − 1)-st order polynomial h(λ) of the form

h(λ) = β_{0} + β_{1}(λ − λ_{i}) + β_{2}(λ − λ_{i})^{2} + ⋯ + β_{n_{i}−1}(λ − λ_{i})^{n_{i}−1}.

The condition f(λ) = h(λ) on the spectrum of J yields

β_{0} = f(λ_{i}), β_{1} = f′(λ_{i}), β_{2} = \frac{f′′(λ_{i})}{2!}, …, β_{n_{i}−1} = \frac{f^{(n_{i}−1)}(λ_{i})}{(n_{i} − 1)!}.

Hence

f(J) = f(λ_{i})I_{n_{i}} + \frac{f′(λ_{i})}{1!}(J − λ_{i}I_{n_{i}}) + \frac{f′′(λ_{i})}{2!}(J − λ_{i}I_{n_{i}})^{2} + ⋯

⋯ + \frac{f^{(n_{i}−1)}(λ_{i})}{(n_{i} − 1)!}(J − λ_{i}I_{n_{i}})^{n_{i}−1}.

Using the nilpotent property of the matrix (J − λ_{i}I_{n_{i}}), we obtain

f(J) = [ \begin{matrix}f(λ_{i}) & f′(λ_{i})/1! & f′′(λ_{i})/2! & ⋯ & f^{(n_{i}−1)}(λ_{i})/(n_{i} − 1)! \\ 0 & f(λ_{i}) & f′(λ_{i})/1! & ⋯ & f^{(n_{i}−2)}(λ_{i})/(n_{i} − 2)! \\ 0 & 0 & f(λ_{i}) & ⋯ & f^{(n_{i}−3)}(λ_{i})/(n_{i} − 3)! \\ ⋮ & ⋮ & ⋮ & & ⋮ \\ 0 & 0 & 0 & ... & f′(λ_{i})/1! \\ 0 & 0 & 0 & ⋯ & f(λ_{i}) \end{matrix} ] .

In particular, if f(λ) = e^{λt}, then we have

e^{J t} = [ \begin{matrix}e^{λ_{i}t} & te^{λ_{i}t} & t^{2}e^{λ_{i}t}/2! & ⋯ & t^{n_{i}−1}e^{λ_{i}t}/(n_{i} − 1)! \\ 0 & e^{λ_{i}t} & te^{λ_{i}t} & ⋯ & t^{n_{i}−2}e^{λ_{i}t}/(n_{i} − 2)! \\ 0 & 0 & e^{λ_{i}t} & ⋯ & t^{n_{i}−3}e^{λ_{i}t}/(n_{i} − 3)! \\ ⋮ & ⋮ & ⋮ & & ⋮ \\ 0 & 0 & 0 & ... & te^{λ_{i}t} \\ 0 & 0 & 0 & ⋯ & e^{λ_{i}t} \end{matrix} ] .

4 State Equation Solution of Controlled System

Consider a model of a controlled dynamic system,

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

subject to a given initial condition x(0) = x_{0}, where B is an n × m matrix, and m < n as is usually in practice. Premultiplying both sides of (24) by the integrating factor, e^{−At}, yields

e^{−At}ẋ(t) = e^{−At}Ax(t) + e^{−At}Bu(t).

We then rearrange the above equation as

e^{−At}ẋ(t) − e^{−At}Ax(t) = e^{−At}Bu(t).  (25)

Noting that

\frac{d}{dt} (e^{−At}x(t)) = −Ae^{−At}x(t) + e^{−At}ẋ(t),

we can represent (25) in the form

\frac{d}{dt} (e^{−At}x(t)) = e^{−At}Bu(t).

Integrating both sides of the resulting equation gives

e^{−At}x(t) − x(0) = ∫_{0}^{t} e^{−Aτ}Bu(τ)dτ.

Hence

x(t) = e^{At}x(0) + ∫_{0}^{t} e^{A(t−τ)}Bu(τ)dτ.

If the initial condition is x(t_{0}) = x_{0}, then integration from t_{0} to t yields

x(t) = e^{A(t−t_{0})}x(t_{0}) + ∫_{t_{0}}^{t} e^{A(t−τ)}Bu(τ)dτ  (26)

As in the case of the uncontrolled system, we can also use the Laplace transform technique to obtain the solution to (24) subject to the initial condition x(0) = x_{0}. Indeed, taking the Laplace transform of (24) and performing some manipulations yields

X(s) = (sI_{n} − A)^{−1} x(0) + (sI_{n} − A)^{−1} BU(s).

Hence

x(t) = L^{−1} {(sI_{n} − A)^{−1} x(0)} + L^{−1} {(sI_{n} − A)^{−1} BU(s)} .

If x_{0} = x(t_{0}), then

x(t) = L^{−1} {e^{−t_{0}s} (sI_{n} − A)^{−1} x(t_{0})} + L^{−1} {e^{−t_{0}s} (sI_{n} − A)^{−1} BU(s)} .

5 From Continuous to Discrete Models

Our objective here is to present a method for constructing linear discrete-time models from given linear continuous-time models.

We use a sample-and-hold device that transforms a continuous signal, f(t), into the staircase signal,

f(kh), kh ≤ t < (k + 1)h,

where h is the sampling period. The sample and zero-order hold (ZOH) operation is illustrated in Figure 1.

Figure 1: Sample and zero-order hold (ZOH) element operating on a continuous function.

Suppose that we are given a continuous-time model,

ẋ(t) = Ax(t) + Bu(t), x_{0} = x(t_{0})  (27)
y(t) = Cx(t).  (28)

The solution to the state equation is

x(t) = e^{A(t−t_{0})}x(t_{0}) + ∫_{t_{0}}^{t} e^{A(t−τ)}Bu(τ)dτ.  (29)

We assume that the input to the system is generated by a sample and zero-order device and has the form,

u(t) = u(k), kh ≤ t < (k + 1)h.

Let t_{0} = kh and t = (k + 1)h and let us use shorthand notation,

x[k] = x(kh)

Then taking into account that u(t) is constant on the interval [kh, (k+1)h) we represent (29) as

x[k + 1] = e^{Ah}x[k] + ∫_{kh}^{(k+1)h} e^{A(kh+h−τ)}Bu[k]dτ
= e^{Ah}x[k] + ∫_{kh}^{(k+1)h} e^{A(kh+h−τ)}Bdτ u[k].  (30)

Consider now the second term on the right-hand side of the above equation. Let

η = kh + h − τ.

Then we can represent (30) as

x[k + 1] = e^{Ah}x[k] + ∫_{kh}^{(k+1)h} e^{A(kh+h−τ)}Bdτ u[k]
= e^{Ah}x[k] + ∫_{0}^{h} e^{Aη}Bdη u[k]
= A_{d}x[k] + B_{d}u[k],  (31)

where

A_{d} = e^{Ah} and B_{d} = ∫_{0}^{h} e^{Aη}Bdη.

The discrete output equation has the form

y[k] = Cx[k].

We point out that the above discrete model, employing the zero-order hold, is exact in the sense that the state, input, and output of this model are exactly the same at the sampling times as the state, input, and output of the underlying continuous-time model.

Example 10 We construct an exact linear discrete state-space model for the continuous system model,

ẋ(t) = [ \begin{matrix}0 & 1 \\ 0 & 0 \end{matrix} ] x(k) + [ \begin{matrix}2 \\ 1 \end{matrix} ] u(t)
y(t) = [ 1 0 ] x(t).

The sampling period is h = 1 sec.

We first evaluate A_{d} = e^{Ah}. Note that A^{2} = O. Hence

A_{d} = I_{2} + Ah
= [ \begin{matrix}1 & 1 \\ 0 & 1 \end{matrix} ] .

Next, we evaluate

B_{d} = ∫_{0}^{h} e^{Aτ}Bdτ.

We obtain

B_{d} = ∫_{0}^{h} e^{Aτ}Bdτ
= ∫_{0}^{h} [ \begin{matrix}1 & τ \\ 0 & 1 \end{matrix} ] [ \begin{matrix}2 \\ 1 \end{matrix} ] dτ
= ∫_{0}^{1} [ \begin{matrix}2 + τ \\ 1 \end{matrix} ] dτ
= [ \begin{matrix}2τ + \frac{1}{2}τ^{2} \\ τ \end{matrix} ]|_{0}^{1}
= [ \begin{matrix}2.5 \\ 1 \end{matrix} ] .

The discrete model then has the form

x[k + 1] = [ \begin{matrix}1 & 1 \\ 0 & 1 \end{matrix} ] x[k] + [ \begin{matrix}2.5 \\ 1 \end{matrix} ]
y[k] = [ 1 0 ] x[k].

6 Summary and Notes

7 Exercises

Exercise 1 Write a MATLAB program that animates the motion of the double pendulum system discussed in class. Obtain a plot of the motion of the pendulum tip versus time on the time interval [0, 5 sec]. Use the initial condition,

[ θ_{1}(0) θ̇_{1}(0) θ_{2}(0) θ̇_{2}(0) ] = [ 135^{∘} 0 180^{∘} 0 ] .

Recommended step-size for solving differential equations is dt = 0.0025.

Exercise 2 Suppose that x = [ 1 2 3 ]^{⊤} is the representation of the given vector in the natural basis {e_{1}, e_{2}, e_{3}}. Find the representation of this vector in the basis {e′_{1}, e′_{2}, e′_{3}}, where

e′_{1} = e_{1} + e_{3}
e′_{2} = e_{2} − e_{3}
e′_{3} = e_{3}.

Exercise 3 Prove Theorem 1.

Exercise 4 Find a similarity transformation that transforms

A = [ \begin{matrix}λ & 0 & 0 & 0 \\ 0 & λ & 0 & 1 \\ 1 & 0 & λ & 0 \\ 0 & 0 & 1 & λ \end{matrix} ]

into its Jordan form. Do not fix λ; construct a transformation that works for any λ. Verify your calculations using MATLAB’s Symbolic Toolbox.

Exercise 5 Calculate the solution for the system,

ẋ = [ \begin{matrix}0 & −1 \\ 1 & −1 \end{matrix} ] x + [ \begin{matrix}1 \\ 0 \end{matrix} ] u, x(0) = [ \begin{matrix}x_{1}(0) \\ x_{2}(0) \end{matrix} ] ,

where u(t) is the unit step, that is, u(t) = 1(t). Plot the system trajectories on the time interval [0, 2] for x(0) = 0.

Exercise 6 Construct an exact linear discrete state-space model for the system,

ẋ = [ \begin{matrix}0 & −1 \\ 1 & −1 \end{matrix} ] x + [ \begin{matrix}1 \\ 0 \end{matrix} ] u
y = [ 0 1 ] x

The sampling period is h = 0.2 sec. Verify you calculations using MATLAB’s function c2d. Find the transfer function of the discrete systems and its step response. Plot the step response of the discrete system on the time interval [0, 2 sec].

Exercise 7 Use the Kronecker product to solve the following Sylvester’s equation,

AX + XB = C,

where

A = [ \begin{matrix}2.7694 & 0.7254 & −0.2050 \\ −1.3499 & −0.0631 & −0.1241 \\ 3.0349 & 0.7147 & 1.4897 \end{matrix} ] , B = [ \begin{matrix}1.4090 & 0.7172 & 0.7269 & 0.8884 \\ 1.4172 & 1.6302 & −0.3034 & −1.1471 \\ 0.6715 & 0.4889 & 0.2939 & −1.0689 \\ −1.2075 & 1.0347 & −0.7873 & −0.8095 \end{matrix} ] ,

and

C = [ \begin{matrix}−2.9443 & −0.7549 & −0.1022 & 0.3129 \\ 1.4384 & 1.3703 & −0.2414 & −0.8649 \\ 0.3252 & −1.7115 & 0.3192 & −0.0301 \end{matrix} ] .

Verify your solution using MATLAB’s function lyap.

Exercise 8 (i) Does the rank of the matrix

Q = [ \begin{matrix}α^{4} & α^{3} & α^{2} \\ α^{3} & α^{2} & α \\ α^{2} & α & 1 \end{matrix} ] , α ∈ ℝ,

depend on the value of α? Is Q positive definite, positive semidefinite, negative definite, negative semidefinite, or indefinite? Justify your answer.

(ii) Find the range of the parameter γ for which the following quadratic form is positive

definite;

f(x_{1}, x_{2}) = γx_{1}^{2} + x_{2}^{2} + 4x_{1}x_{2}.

(iii) Find the range of values of the parameter ξ for which the function

f(x_{1}, x_{2}, x_{3}) = 2x_{1}x_{3} − x_{1}^{2} − x_{2}^{2} − 5x_{3}^{2} − 2ξx_{1}x_{2} − 4x_{2}x_{3},

is negative semi-definite.

References

[1] C.-T. Chen. Linear System Theory and Design. Oxford University Press, New York, third edition, 1999.

[2] C. Moler and C. Van Loan. Nineteen dubious ways to compute the exponential of a matrix. SIAM Review, 20(4):801–836, October 1978.

[3] D. H. Owens. Multivariable and Optimal Systems. Academic Press, London, 1981.