‘w‘
w

x».

 

 

 

 

 

.‘ .
‘1‘. . 1137
.{t .‘

.i ' ~"

I
~v

 

{ ‘
'uw
I“.

'J

'{u
‘14 .1
~

- - 1
1 \ I 4 M
I
0‘ :421'13‘0‘!‘ -
4 (24:13”, 3‘ . ‘.
. 1'! i"

. w?
, ~ ‘3’,“-
'. (“‘73 V ‘ . , .
. , 5i!“ I" (.‘\ I , '. ‘ ' r»!{‘.{ ‘ - I I " h‘ " ‘ I
v I ‘i ' ’ I . ' . . ‘ ) I ‘ ‘ "1 . 4 >
. J 3““. 3 . I - ,3: u; n I w ,. v v . q ' i- ' - n ‘ V. ' ‘i ‘ I ' "f ' '. “fa":
. \ n ' '4' h ‘5" “5""
. . 1‘. “.1 9'
. .

‘1
V‘ r:
hm“.
. (L. ‘
.; .I .d‘
\. n I l '

.. ‘1}? [x a“ .
"\fx'fi'sreff '-‘

:‘t
1" ‘33:!“
. , ..
1!.“ J!

L.“ '
~?-:.rf-).'
.L a: ;
. .« - ..
r397)" ”2—4:;
w

A
1 "‘g
I.“ 1 1‘
”1:15.39:
2 “5H

,1

 

 

 

 

 

 

 

 

 

gmofltflj

MICHIGAN STATE UNIVERSITY LIBRARIES

III/I l/l/I/llll ///I//////I///II//////I Ill/II/II/I/lllllllllll

,fi 7 7 3 1293 00569 1823

" LIBRARY *
Michigan State
University 4—1

 

 

 

This is to certify that the

dissertation entitled

MODELING OF SINGULARLY PERTURBED SYSTEMS
WITH APPLICATION TO MARKOV CHAINS

presented by
Rabah Wasel Aldhaheri

has been accepted towards fulfillment
of the requirements for

Ph.D. degree in Electrical
Engineering

 

 

Major professor
Date W I 98’ 8'

MSU is an Affirmative Action/Equal Opportunity Institution 042771

 

 

 

MSU

LIBRARIES

v

 

 

RETURNING MATERIALS:
Place in book drop to
remove this checkout from
your record. FINES will
be charged if book is
returned after the date
stamped below.

 

 

 

 

 

 

 

MODELING OF SINGULARLY PERTURBED SYSTEMS
WITH APPLICATION TO MARKOV CHAINS

By

Rabah Wasel Aldhaheri

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Electrical Engineering and
Systems Science

1988

A

U

ABSTRACT

MODELING OF SINGULARLY PERTURBED SYSTEMS
WITH APPLICATION TO MARKOV CHAINS

By

Rabah Wasel Aldhaheri

A new method for modeling a two—time scale system in the singularly perturbed
form is presented. The method uses an "ordered" real Schur form decomposition
which can be efficiently computed using standard subroutines from EISPACK. Three
results are given in Chapter 2. First, it is shown that any two-time scale system can
be modeled in the singularly perturbed form via a transformation into an "ordered"
real Schur form, followed by balancing. Second under some conditions, the fast or
the slow state variables can be chosen among the physical variables. Third, a pro-
cedure is given to achieve modeling by permuting the original state variables. The

conditions under which this latter procedure works are necessary and sufficient.

In the third chapter of this dissertation, a nearly completely decomposable Mar-
kov chain that has many applications in queueing networks and inventory control is
considered. This class of system is characterized by high dimensionality and ill-
conditioning, which motivates researchers to look for efficient techniques to over-
come these problems. We propose a general transformation to decompose this large
and ill-conditioned Markov chain system into a reduced-order and well-conditioned

aggregated system. It is shown that the present transformations known in the

Rabah Aldhaheri

literature are a subclass of our proposed transformation. A block diagonal transforma-
tion which is a subclass of the general one is also proposed to simplify and reduce

the amount of computation involved in computing the aggregated matrix.

In the fourth chapter, an algorithm is developed to compute the optimal policies
that attain the minimum average cost per stage for controlled Markov chains. This
algorithm overcomes the large dimensionality and the ill-conditioning problem associ-
ated with many controlled Markov chain problems. Suboptimal policies can be com-

puted with less computations via a modification of the above algorithm.

 

ACKNOWLEDGMENTS

I would like to thank my major advisor, Dr. Hassan K. Khalil for his patience,
guidance and encouragement during the course of this research. Thanks to my gui-
dance committee: Dr. R. Schlueter, Dr. R. O. Barr, and Dr. D. Yen for their comments

and suggestions.

Finally, I wish to express my sincere thanks and appreciation to my parents for
their continued prayers and encouragement. Special thank to my wife, Fadila, and my

children for their patience, help and understanding throughout my studies.

iv

Table of Contents

LIST OF TABLES .....................................................................................................
LIST OF FIGURES ...................................................................................................
INTRODUCTION ......................................................................................................

2 A REAL SCHUR FORM METHOD FOR MODELING
SINGULARLY PERTURBED SYSTEMS ...........................................................

2.1 Introduction ...............................................................................................
2.2 Real Schur Form .......................................................................................
2.3 General Modeling Method .......................................................................
2.4 Modeling with Physical Fast or Slow Variables .....................................
2.5 Modeling by Permutation .........................................................................
2.6 Summary ...................................................................................................

3 EXACT AGGREGATION OF NEARLY-COMPLETELY
DECOMPOSABLE MARKOV CHAINS (NCDMC) ..........................................

2.1 Background ................................................................................................
3.2 Two—time Scale Decomposition of NCDMC ..........................................
3.3 Properties of Exact Aggregated System ..................................................
3.4 First Order Approximation of the Aggregated System ...........................

3.5 Higher Order Approximations of the Perron-Frobenius

Eigenvector of NCD Systems ..................................................................

3.6 Block Diagonal Decoupling Transformation ...........................................

vii

viii

l

5

5

10

l3

17

26

3O

33

33

40

43

48

53

55

 

 

 

 

3.7 Computational Considerations ................................................................. 57

4 OPTIMAL AND NEAR OPTIMAL CONTROL FOR NEARLY
COMPLETELY DECOMPOS ABLE CONTROLLED

MARKOV CHAINS .............................................................................................. 64
4.1 Problem Statement .................................................................................... 64

4.2 Algorithms for Computing the Optimal Control Policy ......................... 67

4.3 Decomposition of the Average Cost for a Fixed Policy ......................... 74

4.4 Optimal Policies for the Average Cost Problem ..................................... 78

4.5 Near Optimal Policies for the Average Cost Problem ............................ 81

4.6 Example: Minimizing the Average Cost Per Stage ................................ 87

5 CONCLUSIONS AND FURTHER RESEARCH ................................................ 94
APPENDIX ................................................................................................................ 97
LIST OF REFERENCES ........................................................................................... 99

vi

2.1
4.1

4.2

4.3

LIST OF TABLES

The IL"... of the ten rows of the matrix Q ................................................................ 29
The optimal policies and the average cost for

different values of 6 evaluated by the algorithm 4.4 ............................................... 90
The near optimal policy and the corresponding

average cost ................................................................................................................ 91
The optimal policy and the minimum average cost

evaluated by he policy iteration method [41]

applied to the full order system ................................................................................ 92

vii

LIST OF FIGURES

2.1 RC circuit with R3 much larger than R1 and R2 ........................................................ 21
2.2 Modeling flow chart ................................................................................................. 32

viii

CHAPTER 1

INTRODUCTION

Analysis and modeling of large scale systems is one of the challenging tasks to
system analysts and control engineers. Difficulties stem from the fact that most large
scale systems suffer from the so called " curse of dimensionality." Of course, the
problem will be worse if the large dimension of the system model is accompanied by
ill-conditioning problems. So, this situation creates an urgent need for efficient model

reduction and decomposition methods.

The presence of multiple time structures in many large scale systems made the
singular perturbation technique one of the most effective methods for decomposing a
large system into smaller subsystems of lower dimensions. The analysis and design
techniques for singularly perturbed systems have been well documented in [1, 2]. The
first task we face in applying singular perturbation techniques is the modeling task,
that is, how to put the multi-time scale system in the singularly perturbed form. A
few approaches are available in the control theory literature to handle this task. Some
of those approaches rely on the underlying physical properties of the system. They
require a very deep insight into the physical nature of the system in order to achieve
the modeling step. Other approaches start from a given mathematical model, without
further insight into the physical nature of the system; yet, they require computing the
eigenvectors of the system [5, 6, 9] and in this case it is assumed that the system has

distinct eigenvalues.

In Chapter 2 of this dissertation, a new method for modeling singularly per-

turbed systems is presented. This method makes use of the real Schur form

decomposition [16, 17]. It is numerically stable and very efficient, and it uses stan-
dard subroutines from EISPACK [18]. Three results are given in this chapter, first, a
systematic procedure is given to put any two-time-scale system in the singularly per-
turbed form by transforming the system into an ordered real Schur form followed by
balancing its elements. This procedure does not preserve the physical state variables
of the system. It introduces new variables composed of linear combinations of the
physical ones. Second, we give conditions under which the system can be modeled in
the singularly perturbed form with all the fast or the slow state variables chosen from
the original physical variables. Third, we give necessary and sufficient conditions
under which all the state variables can be preserved as physical ones and the singu-

larly perturbed form is achieved by permutation only.

Finite state Markov chains is one of the applications where large scale systems
techniques have been applied. This area is very rich theoretically and it is popular to
so many researchers. However, the practical usefulness of this area has been severely
limited due to the extremely large dimensions of most Markov chains. The computa—
tional burden of these problems has disscourged researchers and engineers from using
Markov chains for modeling purposes. In the past few years, some applications in
queueing networks [25] and hydraulic scheduling [30, 44] considered a class of Mar-
kov chains where the states can be clustered into a small number of groups such that
there is a strong interaction within each group while the interaction among the groups
is weak compared to the interaction within the groups. This class is known as weakly
coupled Markov chains or nearly completely decomposable Markov chains. In addi-
tion to the high dimensionality, this class is ill-conditioned. The strong and the weak
transition probabilities between the states can be viewed as a two-time-scale property.

Thus singular perturbation techniques can be applied to this class of systems.

In Chapter 3, a general transformation is given to decompose the nearly com-
pletely decomposable Markov chain into slow and fast parts and as a result of this, a
high order and ill—conditioned system is reduced to a small order and well-
conditioned system. This transformation enabled us to compute the exact as well as
0(8") approximations of the steady state probability distribution of a Markov chain
which is usually encountered in queueing network problems [25, 46]. It is shown that
this transformation is more general than the one considered in [11, 30]. In fact, the
transformation in [11, 30] is shown to be a subclass of this transformation. Moreover,
it is shown that all the transformations that satisfy the conditions of Section 3.2 pro-
duce the same 0(a) approximation of the aggregated matrix. A special class of "block
diagonal transformations" is also given to simplify and reduce the amount of compu-
tations required to form the aggregated matrix. This transformation depends only on

the dimensions of the subsystems. It is independent of the system parameters.

In Chapter 4, we consider the controlled Markov chain of the same class of sys-
tems treated in Chapter 3. The control problem which we consider in this chapter is:
What are the optimal policies that minimize the average cost per stage over the
infinite horizon? [40]. Again this problem is very ill-conditioned and has very large
dimension. So, the conventional algorithms are not practical to handle this type of
problems. In this chapter, we specialize the policy iteration algorithm of [41] to
computing the optimal policies of nearly completely decomposable Markov chains.
The proposed algorithm has two steps: a value determination step where the average
cost and certain dual variables are computed for a fixed policy, and a policy improve-
ment step where the dual variables from the previous step are used to compute a new
policy by minimizing a Hamiltonian function. This new policy is passed to the value
determination step and the cycle is repeated till convergence occurs. The new feature
about this algorithm is that in the value determination step, the algorithm computes

the average cost per stage by using the reduced order aggregated matrix which we

developed in Chapter 3 while in the policy improvement step, the algorithm performs
the minimization on the full system. The reason for doing it this way is that in con-
trolled Markov chains, each row of the probability transition matrix depends on one
particular control; therefore, the minimization step can be performed point wise, i.e.,
each row can be minimized separately. This proposed algorithm would not be possi-
ble if the transformation of [11,30] is used because that transformation depends on

the control policy.

Suboptimal policies can be computed with less computational effort. A modified
version of the above algorithm is given to compute the suboptimal policies and the
associated average cost per stage. This modified version is independent of e. The

convergence of this algorithm is shown in Section 4.5.

Finally, in Chapter 5 we draw conclusions and point to a number of research

directions to which this dissertation leads.

CHAPTER 2

A REAL SCHUR FORM METHOD FOR MODELING
SINGULARLY PERTURBED SYSTEMS

2.1 Introduction

Singular perturbation methods have been used in a wide range of control prob—
lems, see for example [1,2]. The interest in singular perturbation methods stems from
their ability to exploit the multiple time scale structure of dynamic models to decom-
pose analysis and design problems into simpler problems which are solved using
simpler models. The first task we face in applying singular perturbation techniques is
a modeling task, that is, how to mathematically describe a multiple time scale system
in the singularly perturbed form. This modeling task is usually achieved with a great
deal of physical insight into the system, see for example [1; chapter 1] and [2;
chapters 4 and 5]. Attempts to approach this modeling task starting from a given
mathematical model, and without additional information about the underlying physi-
cal system, have mainly concentrated on linear time-invariant systems, which is the
class of systems considered in this chapter; In this case the standard singularly per-

turbed form (in a fast time scale) is described by

221.1(8) mire] [x1]
(2. 1)

’5 = '4" = [7121(2) 1122(2) x;

where e is a small positive parameter, x1 2 R" and x2 6 R’”. The matrix A22 is non-
singular uniformly in e, i.e., 326422) 2 K > 0, where g(.) denotes the minimum singular

value and K is independent of e. The matrices Ail- are 0(1), i.e., their elements are

 

bounded from above by constants independent of e. 1' The singularly perturbed sys-
tem (2.1) is a two-time-scale system in the sense that its eigenvalues are clustered

into two groups; 11 slow eigenvalues which satisfy
I x, l s K,- 5, i=1, ....... ,n (2.2)
and m fast eigenvalues which satisfy
l k,- l 2 Ki, i= n+1, ....... ,n+m . (2.3)

Any matrix A whose eigenvalues satisfy the conditions (2.2) and (2.3) will be
referred to as a two-time-scale matrix. In general, a two-time-scale matrix does not
have to be in the singularly perturbed form (2.1), since this form is dependent on the

choice of state variables.

The modeling problem for linear time-invariant systems can be stated as fol-
lows: given a two-time-scale matrix A, find a similarity transformation 8 such that
S“AS is in the singularly perturbed form. From a theoretical viewpoint, this problem
can always be solved by transforming A into its real Jordan form. This solution,
however, is usually not acceptable from a practical viewpoint for three reasons. First,
the computation of the real Jordan form presents severe practical difficulties when A
is defective or close to a defective matrix [3]. This is particularly harmful when the
real Jordan form computation is done for the purpose of modeling a two-timescale
system in the singularly perturbed form. A major advantage of singular perturbation
techniques is a time scale decomposition based on the clustering of the eigenvalues
into slow and fast groups. Such decomposition will be well-conditioned even when

there are multiple or very close eigenvalues within each group [1]. This advantage

 

T Notice that 0(8) does not imply a lower bound. Numerically a quantity is considered 0(8)
if its magnitude is sufficiently smaller than %. For example, if 8 = 0.01, then 10 is 0(1); also
0.01 is 0(1) since it is 0(8) and any 0(8) quantity is 0(1).

will be lost if we have to compute the real Jordan form because multiplicity or
almost multiplicity of eigenvalues within each group will result in ill-conditioned
computations. Second, the computation of the real Jordan form could be burden,
especially when the dimension of A is high. Third, the transformation into real Jor-
dan form will, in general, change all state variables. This is a disadvantage when the
original state variables are physical ones. In such a case it is preferable in solving the
modeling problem to try to retain as many original state variables as possible. The
optimum case, of course, is to be able to model the system by only permuting the

original state variables. However, this in not always possible.

Previous work on the modeling problem is given in [4] — [15]. In [4] — [8] the

main element of the modeling procedure is finding a solution of a matrix Riccati-type

I, 0
equation, L e R’” x " say, such that the similarity transformation [L l ] transforms the

31*

system into a block triangular form [0 B] with the eigenvalues of 81 being the

2 .
slow ones and the eigenvalues of B; being the fast ones. The methods of [4] - [6]
start by grouping the state variables into two groups and then check a norm condition
that is sufficient for the existence of L. The three methods differ in the way they
achieve the state grouping. It is done by comparing row norms of the matrix A in
[4], by comparing row norms of the eigenvectors corresponding to slow eigenvalues
in [5] and by comparing row norms of all eigenvectors in [6]. In [7] the matrix L is
determined from the eigenvectors corresponding to the slow eigenvalues. In [8] it is
determined without calculating the eigenvalues and the eigenvectors of A, but with
manipulation of the characteristic polynomial of A that requires decomposing this

polynomial as a product of two polynomials corresponding to the slow and fast

eigenvalues.

 

Another approach is taken in [9], where left and right eigenvectors are computed
to determine the participation of each mode in every state variables. This information
is used to classify state variables as slow and fast ones. The method is useful only
when the modeling problem can be solved by permutation. In all above methods,
whenever eigenvectors are computed it is assumed that the corresponding eigenvalues
are distinct.

A special case that received a great deal of attention is the case when the matrix

A can be represented as
A=Ao+8A1. (2.4)

When A0 is nonsingular, the system has only one time scale as 8—) 0. When A0 is
singular with semisimple null structure, the system has at least two time scales and
can be transformed into the singularly perturbed form by a transformation that is
dependent only on A0. One such transformation is given in [10]. If the rank of the

matrix A0 is assumed to be m and its dimension to be (u + m), then the transformation
P1 . n x (Mm) °
rs chosen such that the rows of P, e R span the left null space of A0, i.e.,
P2

PIAO = 0, while the rows of P2 6 R"x(”*"') span the left range space of A0. Applying

this transformation to (2.1) gives us

5’1 [8141105) 8A12(€)] [Yr]
(2.5)

5'2 = 5 A21(8) A2203) Y2

where, y, = Plx and y2 = sz. This case and its multiple time scale extensions have
been studied by several investigators, e.g., [101-[15]. It should be noted that the
system (2.5) is a special case of the singularly perturbed system (2.1), where the
two-one block is 0(a). This special form is a result of the choice of P2 described

above. It is shown in [13] that one does not have to choose P2 such that its rows

span the left range space of A0. Rather, it is enough to choose P2 as any matrix that

PI

will result in a nonsingular [ . such a choice will still bring the system into the

P2

 

singularly perturbed form (2.1), but A2105) will be 0(1), in general. This observation
on the choice of P2 is behind the modeling with physical fast variables idea of sec-
tion 2.4. It is also behind the modeling of Markov chains in Chapter 3.

In this chapter a new method for modeling singularly perturbed systems is
presented. The method makes use of the real Schur decomposition [l6] - [17]. In
Section 2.2 we provide some background material on the computation of the real
Schur decomposition. The computation is numerically stable and efficient. It uses
standard subroutines from EISPACK [18] and from [20,21]. In fact, the real Schur
decomposition is the standard routine to compute eigenvalues in EISPACK. There-
fore, the computational effort needed in‘ our method is basically the computational
effort needed to compute eigenvalues using EISPACK. The method can handle
matrices of dimension of the order of a few hundreds. The stable Schur decomposi-
tion can almost always be used in lieu of the Jordan decomposition [16]. It is gaining
popularity in solving control problems, as for example, in solving Riccati equation

[19].

In Section 2.3 we give a procedure for modeling any two-time-scale system in
the singularly perturbed form by transforming the system into an "ordered" real Schur
form, and then balancing the elements of the Schur form. As a byproduct, the singu-
larly perturbed system will be block triangular. In Section 2.4 we give conditions
under which the system can be modeled in the singularly perturbed form with all fast
or all slow variables chosen fiom the original state variables. It is shown that those
conditions hold whenever A takes the form (2.4) and A0 is singular with semisimple
null structure. Finally, in Section 2.5 we give necessary and sufficient conditions for

the existence of a permutation that transforms a two-time-scale system into the

 

10

singularly perturbed form.

2.2 Real Schur Form

For any real matrix A e N“ there exists an orthogonal matrix Q e R'x',i.e.,

QTQ = 1,, such that

 

 

311 $12 . . . SIN-1
0 $22 . . . Sm
. 0 s33 . . .
QTAQ=T= . (2.6)
p . o . . sNN‘

where each S,-,- is either 1 x 1 matrix corresponding to a real eigenvalue or 2 x 2
matrix corresponding to a pair of complex conjugate eigenvalues of A. The transfor-
mation (2.6) is called the real Schur decomposition, and T is called the real Schur

form (RSF) of A. Moreover, suppose T is partitioned as

T11 T12] (2 7)

_ T _

T' Q AQ ‘ [0 7'22
where the blocks T11 e R” "P and T22 e 119"?”< ('7’) have no common eigenvalues i.e.,
MT“) n MTn) = Q, then the first p orthogonal columns of Q span the invariant sub-

space associated with the eigenvalues of T11 [16].

The most common general purpose algorithm used to compute (2.6) is the
double-Francis QR algorithm . Before the QR process is applied, A is initially
reduced to upper Hessenberg matrix H e R' x ' and this is accomplished by a finite
sequence of similarity transformations known as elementary reflector matrices. This
algorithm is very efficient and numerically stable and it is the most preferred method
for calculating the eigenvalues of a real matrix. For further details on the real Schur

decomposition and its computation, the reader may consult any standard text book on

 

11

numerical linear algebra, e.g. [16,17].

The above algorithm does not guarantee any special ordering of the eigenvalues
along the diagonal of T. In our work with two-time-scale matrices we are interested
in a real Schur form in which the eigenvalues are clustered into two separate groups:
slow and fast eigenvalues, i.e., if T is partitioned as in (2.7), then the eigenvalues of
T11 will be the slow (fast) ones and the eigenvalues of T22 will be the fast (slow)
ones. Fortunately, the ordering of eigenvalues along the diagonal of T can be
changed by interchanging adjacent pairs of eigenvalues using orthogonal transforma-
tions [16,pp. 240 - 242]. Stewart [20] introduced an algorithm written as a Fortran
subroutine [HQR3] to compute a RSF with the diagonal blocks S1,- ordered so that the
eigenvalues appear in descending order of absolute value along the diagonal of T.
Subroutine (HQR3) achieves the ordering of eigenvalues we need in (2.7). Actually,
it goes beyond what we need because in our case the eigenvalues of T11 and T2 are
not required to be of any particular ordering. However, for convenience, we will use
this subroutine to compute the RSF (2.7) where in this case T“ and T22 have the fast
and the slow eigenvalues respectively. To fix notation we will proceed with this par-

ticular ordering.

Definition (ORSF): If the diagonal blocks of the real Schur form (2.6) are ordered

such that it can be partitioned as in (2.7) with

T11 6 RM , MT! 1) = {1;}le = the fast ( largest in magnitude) eigenvalues
and

T22 e R""" , MT”) = {79]}: m1 = the slow ( smallest in magnitude) eigenvalues
where n=r-m, then T is called ordered real Schur form (ORSF).

We will always deal with a two-time-scale matrix A with m fast eigenvalues
and 11 slow ones. These will be the same integers used in the ORSF of A. In this

case the first m orthogonal columns of Q form a basis for the invariant subspace of A

12

associated with the fast eigenvalues, MT“).
The following sequence of subroutines is used to compute an ordered real Schur

form (ORSF):

1. ORTHES to reduce the two-time scale matrix A to upper Hessenberg form
using orthogonal similarity transformation.

2. ORTRAN to accumulate the orthogonal similarity transformation used in reduc-
tion of the matrix A to upper Hessenberg form

3. HQR3 to compute the real Schur form in descending order of absolute value
along its diagonal

The subroutines ORTHES and ORTRAN are available in Eispack[18]

 

13

2.3 General Modeling Method

In this section we show that any two-time—scale system can be put into the
singularly perturbed form by an orthogonal transformation followed by scaling. The
first step is to transform the matrix A into an ORSF using an orthogonal matrix Q,

i.e.,

Tu T12] (2 8)

T = QT A Q = [0 T22
where T11, T12 and T22 are m x m, m x n and n x n, respectively, and MT“) and MT”)
are the fast and the slow eigenvalues of the system. We discussed in section 2 how
to perform this step. For this system to be in the singularly perturbed form, the ele-
ments of T“ and T12 should be 0(IX(T11)I) and elements of T22 should be 0(IMT2QI).
That is, T11 and T12 should be 0(1) while T22 should be 0(8). In general the elements
of Ti,- do not satisfy these requirements and balancing may be needed to reduce the
order of magnitude of the elements of T11 and T12 to the order of magnitude of MT“),
and the order of magnitude of elements of T22 to the order of magnitude of MT”).
The idea of matrix balancing is well known in numerical analysis, e.g. [22,23]. The
Eispack subroutine BALANC [18] balances a real matrix by sealing. This subroutine,
however, balances only the irreducible part of the matrix, i.e., that part of the matrix
which can not be reduced to a block triangular structure by permutation. With the
quasi-triangular structure of the RSF, it is clear that subroutine BALANC is not use-
ful in our case. Therefore, we devised our own algorithm to balance the elements of
T. Similar to BALAN C, the diagonal scaling elements are chosen to be exact powers
of the radix base used to implement the algorithm on a digital computer. In this case
no rounding errors will occur if the the algorithm is implemented in floating point

arithmetic.

14

To describe the balancing algorithm, consider the RSF (2.6) and suppose that
the diagonal blocks S,-,- are ordered as in the ORSF. The goal of the balancing algo-

rithm is to satisfy the condition

agents-9|, i=1, ....... ,N (2.9)

where L,- is the largest (in magnitude) element of (S;,-, Sign ........... , SW) and B is the

radix base employed, usually 2, 10 or 16 . If condition (2.9) is satisfied, the system

i
l3

the on-diagonal matrices 55;. Of course 1 x 1 blocks need no balancing. For 2 x 2

will be in the singularly perturbed form, provided that e < . We start by balancing

matrices, the S,-,- elements could be much larger than the magnitude of the complex

eigenvalues. For example, the matrix

-(or+l) -(0t+1)

a + J—
ct-t-l
has eigenvalues :51- :t j 123-, irrespective of the value of or. Hence its elements can be

made arbitrarily larger than the eigenvalues. To balance a 2 x 2 matrix, consider

an 012
[021 022] with complex conjugate eigenvalues I: out: j y, y: 0. Our goal is to

achieve the condition lay- I s B I A. I. We do this in two steps. First we achieve the

conditions

IaiilSBlll and lalzafllsBlllz (2.10)
and then balance the outer diagonal elements by scaling. Notice that the condition
lama” I 3 Bl A I2 guarantees that lag]- l S B l M can be achieved by sealing. Notice
also that (2.10) will always hold when a“ and an have the same sign, because in this

case lag] < 2 IN S. B W and Ian a21l s W2 s B W2. If (2.10) does not hold we can

achieve it by using a Givens rotation, see [16] problem P7 .4—7. Let 6 be given by

15

012 ‘I‘ 021
‘1 (art " “27)2 + (012 + “21)2

csTarr 012 cs 0‘ 812
‘5 C “21 022 "S C = 821 or (2'12)

where c = cos(9), s = sin(8) and 812621 = — 72.

cos(20) = (2.11)

 

 

It can be verified that

Now, if one of the outer diagonal elements, an say, violates the condition

Ialzl s BIN, it is multiplied by B'“, where o is a positive integer determined from

i012I
0 0+1
[3 s T: s a (2.13)

The other outer diagonal element is multiplied by B“. The second condition of (2.10)
guarantees that B" Iaml S B IN. Notice that during balancing the block Sta the blocks
Si], for i >j and Sfl, forj < i , will change.

After balancing the diagonal blocks Sir, we proceed to balance the off diagonal

 

. Li
blocks Silo. For that we compute the ratio W S 50' . If
B“ — S B“. c > 0. (2.14)

S IMSQI

the ith block row is multiplied by B‘° and the ith block column is multiplied by B“.
This procedure is performed starting from the block SNN and moving up, so that mul-

tiplication of the columns by a scaling factor does not alter previously balanced rows.

The previous procedure puts the system in a singularly perturbed form with the

first m state variables as the fast ones and the last n variables as the slow ones. The

0 1,,
block permutation [I 0} will put the system in the standard singularly perturbed

form (2.1).

16

Example 2.1: Consider the following system

115.810 -4.548 —13.181 -13.862 —0.104
13.768 3.258 13.677 13.213 0.402
)2: 2.745 3.844 -0.703 —0.736 1.093 x
5.316 -0.938 6.153 6.955 -1.022
1-2-447 —2.981 —5.547 —1.183 -2.979.

 

The eigenvalues of this system are :

-5.032. -1.97:l: j 0.1431 and — 0.154 :t j 0.149,

 

which can be clustered as three fast eigenvalues and two slow eigenvalues. A reason-

able choice of e is e = 0.1.

The ORSF of this matrix, computed using the procedure of section 2.2, is

 

 

15.032 —0.133 0.747 32.530 11.377.

0 -l.736 -3.5025 —11.350 0.4454

r: 0 0.0215 —2.204 4.450 -3.550
0 0 0 - 0.2431 0.0555

_ 0 o 0 —0.5421 -0.0646‘

Taking B = 2, it is clear that T is not balanced; so we should apply the balancing
algorithm we discussed in this section. The 2 x 2 diagonal blocks have diagonal ele-
ments of the same sign, therefore no rotation is needed in this example and the diag-

onal scaling matrix that put the system in singularly perturbed form is

D = diag [ 2, 2, 2, -;—, 1 ]. For this D, the ORSF becomes

q

 

 

15.032 -0.133 0.747 8.1325 5.6885

0 -1.736 —3.5025 --2.8375 0.2227

D'lTD= 0 0.0215 -2.204 1.1125 —1.775
0 0 0 —0.2431 0.111
_ 0 0 0 —0.2711 -0.06464

Finally, block permutation yields the following singularly perturbed form

 

 

‘- 02431 0.111 0 0 0
—0.2711 —0.0646 0 0 0
8.1325 5.6885 -5.032 -0.133 0.747
-2.8375 0.2227 0 —l.736 -3.5025
_1.1125 —1.775 0 0.0215 -2.204‘

17

2.4 Modeling With Physical Fast or Slow Variables

When the state variables of a given two-time-scale system are physical vari-
ables, it is usually preferable in transforming the system into the singularly perturbed
form to retain as many physical variables as possible. The general modeling method
of section 2.3 does not , in general, retain any of the original states variables. In this
section we show that if some additional conditions are satisfied, then it is possible to
obtain a singularly perturbed model in which all the fast or all the slow variables are

chosen from the original physical state variables.
Theorem 2.1
Let A be a two-time-scale matrix, and (2.8) be an ordered real Schur decompo-

sition. Suppose that T11 and T12 are 0(1), T22 is 0(8) and there exists a permutation

matrix matrix P such that

A W11 W12
1’ Q = P [Q1 Q2] = W= W21 W22 (2-15)
where the m x m matrix W2? is 0(1). Then the transformation matrix
1,, —W11W2'11- e E
U‘1P= 0 I P (2.16)

where E is 0(1), puts A into the singularly perturbed form.

Remarks

1. The permutation of the rows of Q is equivalent to the permutation of the origi-
nal state variables. It does not affect the ORSF. In particular, T = QTAQ can be
written as T = (PQ)T P A PT (PQ).

2. The in columns of Q1 are orthogonal. Hence it is always possible to permute its
rows such that W21 is nonsingular. The assumption that W2]1 is 0(1) is an addi-

tional requirement on the linear independence of the rows of W21. It is

18

equivalent to saying that the minimum singular value of W2, is bounded from

below by a constant independent of 8.
3. The assumption that W231 is 0(1) guarantees that both U and U '1 are 0(1).

4. The term 8 E in (2.16) is included to allow dropping 0(8) terms of Wu W211,

which simplifies the transformation.
Proof
From the fact that WTW = I, , the following expressions can be verified :
' W127 = W12 - W11 W511 W22
‘ Wiz W11 Wzil + Wiz = 0
° Wzil = Wit + Wit W11 W511

The first of these expressions shows that W13 is 0(1). Using the above expressions,

we obtain
" 5 5 W21 WET - 5 5 W22 T11 T12 Wit W2i1+€ WirE 311 312
U‘1PAPTU = r r 3'
W21 W22 0 T22 W12 8 W125 321 322
where

311 = W12T722Wi2 " EE [W21T11Wir + Werrzwiz + szTzzWizl g 8 411
3,, = -e [szrTuWEi - Wa’TzzWizsr — e2 E
x [werrrwir + Wertzwiz + W22722W1T2] E 3 e 412
321 = WerrrWir + W21T12Wi2 + W22722Wi2 2 A21
322 = W21711W211+ E [W21T11Wi1 + Werrzwiz + W22T22Wi2] E g 422
It is clear that 4,,- are 0( 1); hence the system is in the singularly perturbed form.

Theorem 2.1 shows a procedure to model a two-time-scale system in the singularly

perturbed form with physical fast variables. The significance of this theorem depends

19

on how general its conditions are. The following theorem shows that the conditions
of Theorem 2.1 hold for the case A =A0 + 8A1(8), where A0 is singular and has a
semisimple null structure. This case, and its multiple-time scale extensions, have been
studied by many researchers, e.g. [10] - [15].
Theorem 2.2

Let A(8) = A0 + 8 A1(8), where A0 is singular and has a semisimple null structure,
A1(8) is 0(1), and 8 > O is sufficiently small. Then any ordered real Schur form
decomposition of A has the properties
(i) T11 and T12 are 0(1)
(ii) T22 is 0(8)

W11 W12
(iii) There exists a permutation P such that PQ = W = [W21 W22], where W3} is 0(1).

Proof

The semisimple null structure of A0 implies that there exists an orthogonal

matrix Q0 such that

1 R11 R12
QOAOQO" 0 0

where R11 is nonsingular. Moreover, there exists a permutation of the rows of Q0 such

that

PQ W Wir W92
0 ‘ ° ' we. we.
where W1 is nonsingular. Since W0 is independent of 8, (W‘2’1)‘1 is 0(1).

Therefore A(8) can be represented as

Dr 1(8) 012(8)

A(8) = PTWo 0(a) W3 P = PT W0 8 021(8) 8 022(8)] WgP (2.17)

 

20

where Di,- (8) are 0(1). Using (2.17) in any ordered real Schur decomposition (2.8) of

A(8) yields

T11 =Lir M11 +5551 M21 (218)

T12 = Lit M12 + 8 Lit M22 (2-19)

0 = Ll, Mll + e. 152 M21 (2.20)

722 = LB 14,, + 8 LI; M22 (2.21)

where
L11(8) [42(8)
L“) = W31) 903) = W5 W = [131(6) law]

and

011(8) 012(8) [3411(8) M12(8)]
M(8) =

021(2) 022(2) “8) = Mate) Mate)

 

 

The fact that Li,- and Di,- are 0(1) implies that T“ and T12 are 0(1). Moreover, since
I 7).le 1)| are bounded from below by a constant independent of 8, it follows that, for
sufficiently small 8, Tfil is 0(1). Using this in (2.18) shows that Lfi‘ and Mfil are
0(1), which implies from (2.20) that L12 is 0(8). Finally (2.21) shows that T22 is
0(8). Now, the fact that L12 is 0(8) together with LTL = 1, imply that L21 is 0(8).
Using this in
W21: W31 L11 + W022 [21: W91 L11 + 0(8)
proves that W231 is 0(1), for sufficiently small 8.

Theorem 2.2 shows that there is a large class of interesting physical problems for
which the conditions of Theorem 2.1 are satisfied. It is important, however, to note
that Theorem 2.1 does not require that the matrix A(8) be modeled in the form

A =A0 + 8A1, where A0 is singular with a semisimple null structure. This modeling

 

21

step in itself may not be easy and may require physical insight into the problem. For
example, in [10] this form was obtained by recognizing a network structure with
weak coupling and 8 was taken as a measure of the weak coupling. In another exam-
ple [l], 8 was introduced as a reciprocal of a high-gain parameter in a high gain feed-
back system. Without such insight it may be difficult to model A as A = A0 + 8A1.
Theorem 2.1 alleviates the need for this modeling problem, because its conditions are
checked on an ordered real Schur decomposition of the matrix A itself, i.e., we do

not need to know the matrix A0.
We illustrate Theorems 2.1 and 2.2 by the following two examples.

Example 2.2 Consider the following RC circuit
R1 R2 R3
5"" Ct ’_ C; T C3

Fig. 2.1 RC circuit with R3 much larger than R1 and R2

 

 

 

 

 

 

Choose the capacitance voltages V1, V2 and V3 as state variables and let C1 = C3 = C,
C2 = %, R1=R2= r, R3 =R and i = 8. Here a small value of 8 represents a weak

0 i O t t 0
connection between the capacrtors. Let t= #2:, t = —% and 8 = i, where rd rs the
r

physical time variable.

22

The state equation for the above circuit is given by

—1 l 0 V1
8%=-:—:=A(8)V= 2 —4 2 V, where v: V2
0 1 -1-e V3

which takes the form A(8) =Ao+ 841, where A0 is singular with a semisimple null

structure.

Let 8 = 0.1. In this case the eigenvalues of A are -5.01021, -1.05184 and - 0.03795,
which shows the two—time-scale property of A. The ORSF of A, computed as in sec-

tion 2.2, is given by

—5.01021 — 0.06573 1.41267
0 —1.05184 - 0.00601
0 0 — 0.03795

It is clear that T12 is 0(1) and T22 is 0(8). We can consider Tu to be 0(1), by view-
ing 5.01021 as an 0(1) quantity. This is a close call because :15— = 10. Such a situation

arises because 8 = 0.1 is not very small. For smaller values of 8 the 0(1) elements

will be clearly far from in

The matrix Q, partitioned as in (2.15), is given by

- 0.23483 0.67639 0.69810
Q = 0.94173 - 0.01963 0.33580
- 0.24084 — 0.73628 0.63237

_1 1.05468 - 0.02812 . .
Here Q21 = _ 0.34499 —1.34898 = 0(1). Hence all the conditions of Theorem 2.1

are satisfied with P = 13. Therefore W = Q and

—W11W2‘11=[O.481 0.906]=[0.5 1.0]—[0.019 0.094]

To simplify the transformation we may drop the 0(8) part [0.019 0.094]. This can

be done by choosing - E = [ 0.19 0.94] such that -— W11 W231 — 8 E = [ 0.5 1.0 ].

23

Therefore

1.0 0.5 1.0
U“: 0.0 1.0 0.0
0.0 0.0 1.0

and the singularly perturbed system is given by

 

Ari

dt 0 0

dV2 ‘8 y
_ = 2 .5 0 V2
d1 0 1 —1-8
dv3 V3
dt

 

The eigenvalues of the slow and the fast subsystems are
MA“ — A12 4:3 421) = - 0.036364 and 10122) = — 1.1 and —5 which are C(93) and 0(8)
close to the exact eigenvalues, respectively.

In this example we introduced one state variable y as a slow state and retained V2

and V3 as fast states. 0n the other hand, if we choose the transformation

= -1.0 1.0 0.0 ,

[151] 0.4 0.2 0.4
0.0 1.0 -1.0

the singularly perturbed form (2.5) can be written as

 

id Y1
d t
- 0.48 — 0.168 0.248 Yr
4 )2
T:- = 0 - 3 “-2 yz
8 —2 + 0.48 -3 - 0.68
d YB Y3
d t

 

By using this transformation, none of the state variables are kept as physical ones.

24

Example 2.3 To illustrate that it is not necessary to model the matrix A(8) in the

form A = A0 + 8A1, where A0 is singular, let us consider the following system

—1 1 0
A(8): 1—2 1
0 1 -2-8

It is clear that A0 is nonsingular. This suggests that A(8) has only one time scale.
However, for 8 = 0.1 the eigenvalues for A are - 0.2083 , - 1.6085 and - 3.2832,
which shows a two-time-scale property with one slow and two fast eigenvalues. This
does not represent a contradiction. For, even though A(O) is nonsingular, its
minimum singular value is 0.2083, i.e. , it is 0(8). We apply Theorem 2.1, computing

the ORSF as in section 2.2, to obtain

-3.2832 0 0 - 0.3172 — 0.5869 — 0.7450
T = 0 -1.6085 0 and Q = 0.7243 0.3571 - 0.5898
0 0 - 0.2083 - 0.6122 0.7267 - 0.3118

Notice that the ORSF is diagonal because A(8) is symmetric [16]. All the conditions

of Theorem 2.1 are satisfied with P = 13 and E is chosen to yield

1.0 0.8 0.4
U“= 0.0 1.0 0.0
0.0 0.0 1.0

The singularly perturbed system is given by

- 0.2 — 0.04 0.04
U'bw: 1.0 —2.80 0.60
0 1.0 - 2.10
The eigenvalues of the slow subsystem is — 0.2083, compared with the exact eigen-

value — 0.2083, while the eigenvalues of the fast subsystem are -1.6 and —3.3, com-

pared with the exact eigenvalues —1.6085 and -3.2832.

25

Note that if the conditions of Theorem 2.1 hold on AT, where AT is the transpose
of the two-time-scale matrix A, then all the slow variables can be chosen from the
original physical state variables while the fast states are introduced as a linear combi-
nation of all or some of the original state variables. Such a result can be easily seen

if we apply the transformation matrix (2.16) on AT to yield

U “IPATPTU = (2.22)

 

8311 8312
321 322

where 8,,- are 0(1). If we take the transpose of the two sides of (2.22), we obtain

ePit Bit
UTPAPTU ‘7 = a 31-2 3% (2.23)
1,, 0
Using the scaling matrix S'1 = I , the system (2.23) becomes
0 _"£.
8
88B 88%
s-IUTPAPTU—Ts = (2.24)

Biz Biz

This can be illustrated by the following example

Example 2.4 Consider Example 2.2 again, but this time compute the ORSF of AT(8)

to obtain
—5.01021 0.05682 1.41303 0.40581 — 0.68641 0.60345
T: 0 -l.05184 - 0.01038 and Q = - 0.81370 0.02932 0.58055
0 0 - 0.03795 0.41619 0.72662 0.54663

All the conditions of Theorem 2.1 are satisfied with P = 13. Now, choose B such that

1.0 0.0 0.0
UTP= —1.0 1.0 0.0
—1.0 0.0 1.0

II.

 

26

The singularly perturbed system, (2.24) is given by

0.0 8 0.0
S’IUTAU “’5 = 0.0 - 5.0 2.0
-1.0 0.0 —1 — e

The eigenvalues of the slow and the fast subsystems are - 0.0364 and { -1.1 , — 5},
respectively. Again, these eigenvalues are C(82) and 0(8) close to the exact slow and

fast eigenvalues of the system. In this example the slow variable is taken as V1, while

the fast variables are taken as :15- (V2 — V1) and :1:- (V3 - V1).

2.5 Modeling by Permutation

When can a two-time-scale system be modeled in the singularly perturbed form
by reordering its state variables? The following theorem gives necessary and

sufficient conditions in terms of the ordered real Schur decomposition.
Theorem 2.3

Let A be a two-time-scale matrix, and (2.8) be an ordered real Schur decompo-

sition. There exists a permutation matrix P such that PAPT is in the singularly per-

turbed form if and only if
(i) T11 and T12 are 0(1)
(ii) T22 is 0(8)

on. W11 W12 I o
(111) PQ = W = W21 W22 , where the n x m matrix W11 rs 0(8)

proof:
Sufficiency: From the orthogonality of W we have
Wirwzr = [m — Wirwtr
92(W21) = 1 — 620V“) (2.25)

where g(.) and at) denote the minimum and maximum singular values respectively.

27

Since W“ is 0(8),
6204'”) s K 22 (2.26)

where K is independent of 8. Using (2.26) in (2.25) yields

92(W21)21- K 62

which implies that
6'2 I s ——1
(W21) 1_ K82

Hence, for sufficiently small 8, W5} is 0(1). Thus, all the conditions of Theorem 2.1
are satisfied and the system can be put in the singularly perturbed form via the
transformation (2.16). But in the current case Wqu} is 0(8). Therefore, choosing
8E = - WHWQ}, the matrix U'l becomes the identity matrix, and the transformation

(2.16) reduces to a permutation.

Necessity; This follows as a special case of Theorem 2.2. Suppose there is a permu-

tation matrix P such that PAPT is in the singularly perturbed form, i.e.,

EA11 ISA12] (2 27)

T =

P A P [ A21 A22
Then, all the conditions of Theorem 2.2 are satisfied which implies conditions (i) and
(ii). Moreover, the orthogonal matrix Q, in the proof of Theorem 2.2 can be chosen

0 1,,
as Q0 = P7 [I 0]. Furthermore, the matrix W0 in the same proof can be chosen as

0 1,,
W0=PQ0= 1m 0

Due to this special form of WO, W is given by

w = w0 1. (2.28)

[[21 10.2]
= L11 L12

28

It was shown that L21 is 0(8); hence W11 is 0(8), which completes the proof of the

theorem.

Example 2.5 : Consider the lOth-order system described by

L0.05 0.05 0 0 0 0 0 0 0 q
0.05 -0.46 0.05 o 0.36 0 0 0 0 0
0.05 -0.46 0.05 0 0.36 0 0 0 0
0 0.05 —0.41 0 0 0.36 o 0 0
0.2 0 0 -0.25 0.05 0 0 0 0
A = 0 0.2 0 0.05 -0.66 0.05 0.36 o 0
0 0.2 0 0.05 -0.61 0 0.36 0

0.2 0 —0.25 0.05 o
0 0.2 0.05 -0.61 0.36
0 0 0 0.2 —0.2j

 

 

'oooooooo

COCO
GOO
GOO
GOO

The eigenvalues of this model are

— 1.0121, — 0.9065, - 0.6421, - 0.6063, — 0.3792, — 0.2156, - 0.1077, - 0.0654, —- 0.02516
and 0.

For 4 slow and 6 fast states, take 8 = 0.5. Calculating the ORSF of this model yields

— 1.0121 0.0136 0.0097 0.1550 0.0645 0.00811
0 - 0.9065 0.0169 -0.0578 0.1569 0.0346
0 0 -0.6421 -0.0001 -0.0117 —-0.0083
T“: o 0 0 -0.6063 0.0271 -0.1686 '
0 0 0 0 —0.3792 0.0073
0 0 0 0 0 -0.2156

 

 

q

—0.0071 -0.0054 0.0061 0.0399
-0.0245 -0.0383 -0.0276 0.0021
0.1230 -0.0174 -0.0954 0.0314
712 = -0.0185 —0.0248 —0.0292 —0.0659
0.0650 0.0946 0.0569 —0.0344
—0.0003 -0.0099 -0.0316 —0.1021

 

 

and

29

-0.1077 -0.0160 -0.0051‘ —0.0086
0 -0.0654 -0.0161 —0.0118
T22“ 0 0 -0.0252 -0.0123
0 0 0 0 ]

 

 

This shows that T11 and T12 are 0(1) and T22 is 0(8). In fact in this example T12 hap-
pens to be 0(8). Therefore conditions (i) and (ii) of Theorem 2.3 are satisfied. Condi-
tion (iii) is satisfied if the 10 x 6 matrix Q has four rows which are 0(8). The ll.ll,, of

each row of Q is given in the following table

Table 2.1 The IL", of the ten rows of the matrix Q1

 

 

ll.ll,, 0.0744 0.8782 0.7544 0.6460 0.4307

 

 

ll.ll,, 0.6358 0.6247 0.3746 0.6727 0.3332

 

 

 

 

 

 

 

 

which shows that rows 1, 5, 8 and 10 are 0(8). Permuting the rows of Q accordingly,

we verify that

i0.084 0.003 0.001 —0.003 0.001 0.001
1 0.509 0.094 —0.003 0.109 0.002 0.0
I: 0.002 0.265 0.077 0.507 0.155 0.180 = 0(8)
_ 0 -0.002 0.113 -0.001 0.228 0.478‘

 

 

Thus the system can be modeled in the singularly perturbed form by regrouping the

states 1, 5, 8 and 10 as the slow states and 2, 3, 4, 6, 7 and 9 as the fast states.

30

2.6 Summary

A reliable and numerically stable algorithm is developed to transform any two-
time-scale system into a singularly perturbed form. This is accomplished in two
steps: first, transform the matrix A into an ORSF; second, use the balancing algo-
rithm developed in this chapter to balance the ORSF such that T11 and T12 are 0(1)

and T22 is 0(8).

If we are interested in the physical state variables of the system, sufficient con-
ditions are derived to put the system in the singularly perturbed form whereas the fast
or the slow physical states are retained The conditions hold when the matrix A can
be modeled in the form A =Ao + 841, where A0 is singular with a semisimple null
structure. The significance of these conditions over previous results in this area is that
we do not require that the matrix A be modeled in the above form. Such modeling
step may require a priori knowledge of the physical nature of the system. An addi-
tional assumption on the orthogonal transformation matrix Q yields necessary and
sufficient conditions for a two-time-scale matrix to be modeled in the singularly per-

turbed form while retaining all the states as physical ones.

The modeling procedure for the special cases of Theorems 2.1 and 2.3 can be
combined with the general modeling method of section 2.3 in one modeling pro-
cedure presented in the modeling flow chart of Fig. 2.2. According to this flow
chart, we first compute the ORSF of the matrix A and check whether T11 and T12 are
0(1) and T22 is 0(8). If they are not satisfied, then we use the general method to put
the system in the singularly perturbed form. If they are satisfied, we check whether
there are 11 rows of the matrix Q1 which are of order 0(8). If this is the case, then we
can obtain the singularly perturbed form by permutation; otherwise we check the

existence of a permutation of the rows of the matrix Q such that W} is 0(1). If this

31

is true, then a singularly perturbed form with fast variables chosen from the physical
ones is achieved; otherwise, the general method is used.

32

 

compute
ORSF

 

 

 

T11.T12 = 0(1)

no
and

 

T22 = 0(8)

yes

Are there

es
y 11 rows of

 

 

 

 

Apply the
permutation of
Theorem 2.3

 

 

Qt =0(€)

110

Is there
no

 

a permutation
.1. W211 = 0(1

yes

 

Apply the
transformation
of Theorem 2.1

 

 

 

Fig. 2.2 Modeling flow chart

 

 

 

General
Method

 

 

CHAPTER3

EXACT AGGREGATION OF NEARLY-COMPLETELY
DECOMPOSABLE MARKOV CHAINS (NCDMC)

3.1 Background

Finite Markov chains have many applications in biological, physical and social
sciences, as well as in economics and engineering. Queueing networks which are
used in computer system modeling are one of the areas which received a great deal

of attention in the past few years [25,27 - 28].

A fundamental problem in these applications is to compute the unique stationary

probability distribution vector which satisfies

1tP=rt, 1t,->0, rtl,,=1 (3.1)

where P e R" x " is a transition probability matrix, 1,, is a column vector consisting
of n ones and it is the left eigenvector corresponding to the unit eigenvalue of the
matrix P.

In many applications, the number of the states of a Markov chain model is so
large that a direct solution method, such as an LU factorization, is not possible. Also,

the number of iterations required by using an iterative method, such as the power

method or the

Fortunately, in many large systems [24] the states can be clustered into a small
number of groups such that there is a strong interaction within each group while the
interaction among the groups is weak compared to the interaction within the groups.
This class of systems is known in the literature as nearly completely decomposable

Markov Chain.

33

34

The structure of this class of Markov chains can be represented by the following

transition probability matrix

 

 

P=l,,+A+eB=l,,+Q (3.2)
where
Q=A+83. (3.3)
34, 0 . . 0.
0 A2 0 .
A= 0 0 . (3.4)
L0 . 0 . AN
N
andAiiSannixnimatrixJ:l,2,.....,Nand2n,-=n.
ill

The matrices P and (In, + A9 , i= 1, 2, N are stochastic. Hence, the row sums of B

and A,- are zero. It is assumed that the Markov chain has a single ergodic class, i.e.,
all states communicate with each other. In this case the Markov chain and the matrix
P are called irreducible. This assumption is equivalent to the assumption that equa-
tion (3.1) has a unique solution, The probabilities 1:,- for all i are positive [35, pp.
181-182], and the vector 1: is called the Perron-Frobenius eigenvector, or simply, P-F
eigenvector. Furthermore, the matrix P has a simple unity eigenvalue, and the matrix

Q has a unique zero eigenvalue. It is also assumed that each of the matrices (1,}, + A9

is irreducible. Hence A,- has a unique zero eigenvalue.

The maximum degree of coupling 8 is very small and is defined by

8=|lP—l,,-All,,<<1 (3.5)

The problem of interest is to solve (3.1) or, equivalently, to solve

rtQ=0, 1:1,,=1 (3.6)

35

Simon and Ando [24] were the first to propose this class of systems which is
known as nearly-completely decomposable systems. They gave examples from
economics to illustrate this class of systems. Their analysis is based on the following
argument:

1. Over a relatively short period of time a local equilibrium is reached among each
group separately and in this period we can analyze the system as if there was no

interaction among the groups.

2. In the long run, the interaction among groups can not be neglected and the
whole system moves towards a global equilibrium which is defined by the
steady-state probability. The equilibrium maintained in the short run will remain
approximately the same because the interaction among the groups will not

influence the relations between the states within each group.

The approximate Perron-Frobenius eigenvector was obtained by approximate
aggregation as it will be shown in this section. Courtois [25] developed the theory of
the NCDMC and applied his method to several computer system models [25 - 27].

According to his method the approximate P-F eigenvector is computed as follow:

i. Compute the P—F eigenvector of each block, i.e.,

viAi=Ov V; l~=1, for all i=1,2, ....... ,N (3.7)

ii. Form the approximate aggregated matrix P1 as

"v10 00‘ 1"1 0 0 0

0 V2 . . 0 182 ' A
P1=IN+8. 0 . 0 B . 0 0 =1N+ev13W1=IN+er (3.8)

 

 

 

 

 

Iiii

36

where 1718 RN“ and W1 8 RH” are defined by

 

 

 

 

'v1000' 14000

.. 0 v2 . . 0 11:2 - ~

V1: . 0 . 0 and W1: . 0 . 0 (3'9)
_0 . . v”. 0 . . 1,,”

iii. Compute the P—F eigenvector of P1 that satisfies
X P1 = X, X 1,, = 1
or

XQO=01 X 1N=1 (3.10)

iv. The approximate P-F eigenvector of the whole system is

 

 

'v, 0 0 0l
0 v2 .

1t=X . 0 . 0 (3.11)
_0 . . vN‘

Courtois [26] showed that it is an 0(8) approximation of 7:.

Vantilborgh [29] carried this approximation one more step to come up with an C(82)

approximation of the P-F eigenvector it. He started by partitioning Q = A + 8 B as

a 0.
9‘ 0.0.

where Q 8 RIM") X """N’ and 0,, e 11"" x "N

He then formed a matrix

in = Q: - Q]. Q“ QC (3.12)

a.

37

This partition is not based on decomposing the zero and nonzero eigenvalues.
So, he was confronted with a near singular matrix Q. To overcome this problem, he

made use of the NCD structure of Q to expand it in 8. Thus he obtained

 

 

 

 

Ill 0 0 .VI 0 0
0 1,.2 . . 0 V2 . .
8-1”“ . 0 . . 90-1 . 0 .. (3°13)
lflN-r J . . VN_1-l

where Q0 is obtained by partitioning the matrix Q, of equation (3.10) as

Q0 Q0.
Q0: Q01 40» (3.14)

“’th6 qu iS 1X1 311d Q0 6 RN_1 XN-l.

Again Q0 is a near singular matrix; therefore, another expansion in 8 should be pur-

sued. Hence, equation (3.12) becomes

 

 

 

 

 

 

1"! 0 ' 0 FVI O . 0
l' A 0 In: - ' 1 0 V2 . . A l
QN=QN-QL . 0 . . Q0_ . 0 . . QC=QN‘4Q0 P (3-15)
I’m-1‘ _. . . VN_1j
where q e R"” x(N-l) and p e R(N—1)x "" are given by
1‘1 0 ' 0 -V1 0 . 0
A 0 1,,2 . . 0 v2 . . A
q=QL . 0 . . and p=. 0 .. Qc (3.16)
L lnN-l i vN‘l

 

 

38
The dominant part of Q5‘ is given by [29]
lIII-l( X‘p lap-1X“

where

.. _ 1
l—XN

 

N—l
2X1

i=1

and X is defined by (3.10). Substitute (3.17) and (3.18) into (3.16) to obtain

4 1N—1 2‘ P + 0(82)

 

Q'=Q - ..
N N Xpl,

N

In general, the matrices Q; for l = 1, ....... , N are defined by

1
X?) 1.

I

 

Q; = Q, - q 1~_, 2‘12 + 0022)

(3.17)

(3.18)

(3.19)

(3.20)

where q, p and I? are defined by the same procedure as in block N, i.e., by (3.16) to

(3.18).

Now, this method can be summarized as follow :

1. Compute v, for all i, Q0 and X by using (3.7), (3.8) and (3.10)

2. Compute q and p that correspond to the Ith block as defined in (3.16)

3. Form the matrix Q}, I=1, ..... , N as in (3.20)

4. Repeat steps i - iv of Courtois’ 0(8) approximation method with the new

matrices Q; replacing A,- to obtain an C(82) approximation of rt.

This method requires more than double the computations needed for the 0(8) approx-

imation.

Phillips and Kokotovic [l 1] gave a singular perturbation interpretation to Cour-

tois’ aggregation of NCD systems. They defined the continuous time Markov process

by

39

£2792 = 1‘0) ( 381 + B) , rt(t) 1,, = 1 (321)

where 1t(t) is n—dimensional probability distribution vector at time t

Consider the similarity transformation

[, ]=..[w, a] m. [e n] M (3.22.
2
where the matrices 17, and W, are chosen as in (3.9); hence
V,A=0, AW,=0 and v,W,=IN, (3.23)

while the matrices 172 and W; are chosen to satisfy
172W, = O , 91W; = 0 Vsz = [rt-N

which ensures that

Substituting (3.22) into (3.21) gives

#3th = gm 1?, B W, + MI) 17, B W, (3.24)
£3552 = e §(:)17,BW2 + 11(1) 172(A + 8 B) W2 . (3.25)

Equations (3.24) and (3.25) are in the singularly perturbed form. The 0(8) approxi-
mate stationary probability distribution vector is obtained, as t —-) oo, by setting 8 = 0.

This yields
£0 (1?, B W,) = 0 and no = 0 (3.26)
and the approximate 1: is given by substituting £0 and 110 in (3.22). Therefore,

a: = to 17, + 0(a) . (3.27)

40

Equations (3.26) and (3.27) are equivalent to (3.10) and (3.11), respectively. Notice
that this transformation is similar to the one discussed in chapter 2 where A has a
semisimple null structure, 17, spans the left null space of A and 172 is perpendicular to
the null-space of A. This is the reason for having a weak coupling between 11 and i
in equation (3.25). As mentioned in chapter 2, this transformation is a special case.
There is a wide class of transformations which can put the system in the singularly
perturbed form, but the coupling between the 110) and E,(t) will not be necessarily
weak.

In section 3.2, we will propose a more general transformation to put the system
in the singularly perturbed form and will show later in section 3.6 that the transfor-

mation proposed in [11] is a subclass of this transformation.

3.2 Two-time Scale Decomposition of NCDMC

The nearly-completely decomposable Markov chain matrix Q is very ill condi-
tioned because it has one zero eigenvalue and N - 1 eigenvalues which are very close
to zero. Therefore, a direct method to solve the steady-state problem will suffer from
ill-conditioning. However, this matrix exhibits a two-time scale property which can
be used to decompose the system into slow and fast parts and, as a result of this, the

high order system is reduced to a low order one.

The transformation proposed in this section is more general than the transforma-
tion considered in [11]. This transformation does not put any restriction on choosing

the matrix V1. The transformation can be described as follows:

W, is chosen as in (3.9); hence A W1: 0; W2 can be any matrix such that the

transformation

r = [w, w,] (328)

.l
:1."-

41

is nonsingular. The matrix F‘1 is written as

_ V
r 1= [,4] (3.29)

Therefore, V, W, = IN, V, W2 = 1H,, V, W2 = 0 and V2 W, = 0. The dimensions of

V,, V,, W, and Wzarean,(n—N)xn,an and nx(n-N), respectively.

Since we are interested in the stationary probability disuibution, we apply this
transformation to equation (3.1). Let §= 1: W, and n = 1: W2 then, equation (3.1)

becomes

rt[W, W2] [,2] P = 1:. (3.30)

Multiply (3.30) from the right by I‘ to get

1N + 6 V13 W1 V,(A + e 3) W2
[8 n] 8 V23 W, 1,.” + 11204 + e B) W2 = [5 ‘1] « (331)

Because A has N zero eigenvalues and

I 0 V,A W2

it is seen that V2A W2 is nonsingular. This implies that for 8 sufficiently small,

V2(A + 8 B) W, is nonsingular. From equation (3.31) we have
8 (IN + e V,B W,) + e 11 V23 W, = g (3.33)
and
g V,(A + e B) W2 + 11(1..~ + V2(A + e B) W2) = n . (3.34)
Therefore,

= -§ V,(A + e B) W2[V2(A + as) Wfli—l . (3.35)

.m

42

Substitution of (3.35) into (3.33) yields

a [1,, + 8[V,B W, - V,(A + 23) W2[ V2(A + 23) W2]-1 V23 W,]] = g (3.36)

Let
Q, = [V,B W, - V,(A + e3) W2[ V2(A + as) W2]_l V23 W,] (3.37)
and
P, = [N + 8 Q, . (3.38)
Then
€n=§ 63%
or simply
t Q. = 0 . (3.40)

Again from (3.6), it 1,, = 1, therefore,
1c[W, W2] [2] 1,,=1 . (3.41)

Noting that 1,I = W, 1”, (3.41) can be written as

I
[e nmy]
Hence
t 1,, = 1 . (3.42)

Thus, the solution of the full order system (3.6) is completely characterized by

1‘ =& V1 + I] V2 (3.43)

43

where 5:, is defined by (3.40) and (3.42), and n is defined by (3.35)

The existence of a unique solution of (3.40) and (3.42) is discussed in section

3.3. In sections 3.4 and 3.5, first and higher order approximations of the solution of

equations (3.40), (3.42) are given, respectively.

3.3 Properties of the Exact Aggregated System

The aggregated matrix Q, is a reduced form of the original matrix Q and it

inherits some of its properties. In this section some properties of Q, are given in the

form of lemmas or theorems.
Lemma 3.1

The row sums of the matrix Q, are equal to zero.
Proof :

From (3.37), Q, can be rewritten as
-1
Q, = [V, — V,(A + as) W2[V2(A + 83) W2] V2]B W, = V, B W,
where

V, = [V, — V,(A + 88) W2[ V2(A + 83) W44 V2] .

Recall from section 3.1, that each row sum of B is equal to zero.

Let the ith row sum of B W, = C, then

"1 "r “‘2
C: 2b,.» 212,4 .............. + z b,,=zb,,-=0.
j =n,+l jan—nml j=l

(3.44)

(3.45)

(3.46)

Thus, the ith row sum of B W, = 0. Now we need to prove that the ith row sum of

V, B W, = 0. Denote the ith row of V, by V,,- and the jth column of B W, by Ci. Then

the ith row sum of V, B W, is

44
N . N _

2V..-C’= .rZC’=0 (3.47)
i=1 i=1

and this completes the proof.

As a result of this lemma, the matrix Q, has a zero eigenvalue. Another useful pro-
perty of Q, which will be used in this chapter and in chapter 4 is given in the follow-

ing lemma.
Lemma 3.2

Let the matrix Q, be defined as in (3.37) and ‘11,,- be defined as the matrix Q,

with its ith column replaced by the vector 1”. Then

i. The matrix Q, has a unique zero eigenvalue.

ii. The matrix ‘11,,- is nonsingular, for i = 1, 2, N;

iii. Equations (3.40) and (3.42) have a unique solution and g,- > 0 , i = l, 2, ..... , N.
Proof :

i. From the ergodicity assumption, the matrix Q has a unique zero eigenvalue and

equation (3.6) has a unique solution; therefore,
Rank( Q , 1,, ) = n . (3.48)

The rank of the matrix is invariant to pre and post multiplication by nonsingular

matrices; therefore,

1‘ 0
Rank(Q, 1,)=Rank[l"1(Q, 1,,)[O 1]]=Pank[r"or, r4 1,]
8 V18 W1 V1(A + 8 3) W2 1N
=R‘"‘" 8 V,B W, V,(A + as) W, 0 = "- (3'49)

Since V,(A + 88) W, is nonsingular, an elementary row operation brings the above

matrix into the form

45

8 Q, 0 1N
8 V28 W1 V2(A '1' 88) W2 0

(3.50)
where the elementary row operation matrix is
1,, — V,(A + e B) W,( V,(A + e B) W,)‘1
0 1,.”
Hence
8 Q: 0 1N
““49: 10‘1““ e V,B W, V,(A+eB)W, 0
=n-N+Rank(Q,, 1N)=n.
Therefore,
Rank(Q,, 1N)=N. (3.51)

Thus, Rank( Q, ) = N - 1. Moreover, since 1,, spans the null space of the matrix Q,,

(3.51) implies
R(Q.) e me.) e R”. (3.52)
Hence, Q, has a semi-simple null structure. This implies that Q, has a unique zero

eigenvalue.

ii. To prove that ‘11,,- is nonsingular, let us prove first that any N - 1 columns of Q,
are linearly independent. Any column of Q, can be replaced by a zero column

via an elementary column operation (ECO), which can be represented by

Q, E,- = [q}, 43, ..... , qi‘l, 0 , qf,+1,....,qf] (3.53)

46

where q; is the ith column of the matrix Q, and E, e R” x” is a matrix which can be

obtained by performing the same ECO on IN. It takes one of the following forms:

 

 

 

 

 

 

"10000“ "11000" 10001'
1100. 0100. 01001
El— . O I , E2: 1 I and EN: . 0 1 . . (3.54)
10001, 91001, 00001,
Since E,- is nonsingular,
Rank(Q,E,-)=Rank(Q,)=N—1. (3.55)

Hence, any N — 1 columns of Q, are linearly independent.

Now, the range space of Q, is spanned by any N - 1 columns of Q, and the vector
11,, spans the null-space of Q,. Using (3.52), we conclude that ‘1’,- is a nonsingular
matrix.

iii. The proof of this part follows from parts (i) and (ii). Since equations (3.40),

(3.42) can be rewritten as

i [(13, q}, ...... .qf,“ , 1N,q§+‘, ....... 4:1]: [0. 0, ...... ,0, 1 ,0, ....... ,0]
01'
g‘Ps-i = 6,:

The nonsingularity of ‘1’,- irnplies that the equation has a unique solution. Moreover,

the relationship

§=75 W1

implies that g,- > 0, since 1:,- > O.

47

Since Q, preserves some of the properties of Q, like zero row sums, the unique-
ness of the zero eigenvalue and the property that the left eigenvector corresponding
to the unit eigenvalue is positive and its sum equals 1, a natural question that comes
to mind is whether or not P, = IN + 8 Q, is a stochastic matrix. The answer is no, in

general.

The following example is a counter example that shows that the matrix P, may not

be a stochastic

Example 3.1: Consider the irreducible transition probability matrix

’ 1

 

 

From (3.2) to (3.4) the matrices A and B are

 

 

 

 

—.5.50000 -1-10011
.5-.50000 0—10001
00-5500 11—1—100
A: 0 0 5—.5 0 03nd3=010—100
0 0 0 0—.5 .5 .5 .5 .5 .5—1—1
0 0 0 0 .5—.5 0.5 0.5 0—14
LetI‘beequalto

100000

100100

010000

l"=010010

,001000

001001

 

 

48

From (3.37)

—1.5 - 1.58 - .58 1.5 + 28
Q,= 1.5+28 —1.5—8 —8
0.75 + .58 0.75 + 8 -1.5 - 1.58

It is clear that the mauix P, = IN + 8 Q, is not a stochastic matrix because of the
negative sign of the (12) and (23) elements. Notice, however, that, the matrix

P, = IN + 8 Q,, is stochastic.

3.4 First Order Approximation of the Aggregated System
An 0(8) approximation of the Perron-Frobenius eigenvector can be obtained by

setting 8 = 0 in equation (3.44) to obtain

Q0 = Q,.L 0:- [V1 — V1 A W2( V2 A Wfl-l V2]B W1 = V0 B W1 (3.56)

where
Vo = [V, - V, A W,( V, A W,)'1 V,] (3.57)
and equation (3.38) becomes
P, = [N + 8 Q,, (3.58)

This approximation reduces the amount of computations involved in computing the
matrix V,. In this section the properties of Q0, V0 and P, are presented and a closed

form expression for the P-F eigenvector of the stochastic matrix is given.
Theorem 3.1

Let A be defined as in equation (3.4) and V,, as in (3.57). Then
i. V,, W, =1” (3.59)

ii. V0.4 W,=0 (3.60)

49
iii. Vo A = 0 (3.61)

iv. V0 is a block (row) diagonal matrix and it is nonnegative.

v. The ith diagonal block of the matrix V0 is the P-F eigenvector of the matrix
(Ai + In, )-

Proof :
Parts (i) and (ii) follow by multiplying equation (3.57) from the right by W,

and A W, respectively.

iii. The matrix 1‘ of (3.28) is nonsingular, therefore, the matrices A and A l“ have
the same range space. But,

AF=A[W, W,] = [0 A W,]. (3.62)

Thus, the range space of A is spanned by the (n - N) columns of the matrix
A W,. Therefore, V,, A W, = 0 implies that V,, A = 0.

iv. To prove this part, let

 

 

V11 V12 . . Vm
V21 V22 . . V2”

V,,: . . . . . (3.63)
LVN1 VNz - - VNN

where V,,, for all i and j, are vectors matched with the partitions of the matrix A.
From (i) and (iii)
V,,-A,- = 0, v, 1, = 1 for 1‘: 1,2, ....... ,N (3.64)
By (3.64),
V,,- (Apt-1,.) = V,, and v-- 1 ,=1,t'=1, ...... , N. (3.65)

art,

But (A,- + 1,) is a stochastic matrix with a single ergodic class. Therefore, the vector

50

v,,- is the PF eigenvector which is unique and positive.
Now, to prove that the matrix V0 is block diagonal, we should prove that V,,- A,- = 0 for
i a: j implies that V,,- = 0.
From (iii), V,,- A,- = 0, i at j, or
V,,- (A,- + I,),) = V,,- . (3.66)
Since (A,- + 1,!) is a stochastic matrix with a single ergodic class, v,- must be equal to
K V,,, where K is any constant. From part (i), V,,- 1,, = 1 and V,,- 1,): 0. This shows that
V,,- 1n,- = vaj 1,} = 0 , (3.67)
which implies that K = 0. Hence, V,,- = 0 for all i at j.
This completes the proof of (iv).
v. This part follows from (i), (iii) and (iv)
The last part of the above theorem shows that V,, = V,, where V, is given by
(3.9). Hence, the matrix P, of (3.58) is the same as P, of (3.8). Therefore, P, is
independent of the choice of W,. In other words, all transformations of the form

(3.28) produce the same matrix P,. Furthermore, the expression (3.57) for V0 gives a

closed form expression for the P-F eigenvectors of the stochastic matrices (1,, + A,).

Such closed form expression is important especially when the stochastic matrix is a
function of a certain parameter as in the controlled Markov chain problem of Chapter

4.
Theorem 3.2

The matrix P, of (3.58) is a stochastic matrix with a single ergodic class.
Proof :

From (3.59) and (3.61), V,, W, = 1N and V,, A = 0. Therefore P, = 1N + 8 Q,, can

be written as

51

P1=Vo(IN+A+€B) W1=V0PW1

(3.68)

where P is defined by (3.1). Similar to the proof of lemma 3.1, it can be shown that

the row sums of V,, P W, are equal to one. Moreover, since V0 , P and W, are non-

negative matrices, P, is also nonnegative. Hence, P, is nonnegative and its row sums

are equal to l, which implies that P, is a stochastic matrix. The ergodicity follows

from the ergodicity of P; see, for example, [25].

As a result of this theorem, the matrix P, has a unique left eigenvector

corresponding to the unit eigenvalue, i. e.,
50 P1 = 50 , go 1N = 1
or simply,
50 Q0 = 0
Theorem 3.3
i. P, is an C(82) approximation of P,
ii. 8,, is an 0(8) approximation of §
iii. 1:0, given by
uo=tovr+nov2=§ Vi—to m M v24 W2)“ V2=§ovo
is an 0(8) approximation of 1:.
Proof :

i. From (3.38) and (3.44)
fl=m+8KBWL
Expand V, in 8, to obtain

P, = 1N + a (V0 + 0(8)) 3 W, = P, + 00:2)

(3.69)

(3.70)

(3.71)

52
From lemma 3.2
i W..- = e.- (3.72)

where e, is an N-dimensional row vector with the ith element equal to 1 and the

rest of the elements equal to zero, i.e.,

§[q},q§, ...... ,qu-l,1N,q§+1, ....... ,qfi’]=[0,0, ...... ,0,1.0, ....... ,0]. (3.73)
Similarly,

50 5’01 = 61 (3.74)

80[q,1,,q,2,, ...... ,qgl,rN,q‘+1, ....... fly]: [0,0, ...... ,o,1,0, ....... ,0] (3.75)

where <10 is the ith column of the matrix Q,, and the right hand vector is the same as

the one defined in (3.72). As in Lemma 3.2, the matrices ‘I’,,- and ‘Po, are nonsingular.

If Q, is expanded in 8, we obtain

Q, = Q,, + 0(8) ; therefore, ‘I’,,- = ‘P0, + 0(8) . (3.76)

From equations (3.72) — (3.76)

( E - §o ) ‘Pol = 0(6) (3.77)

Since ‘PO, is nonsingular,

E = to + 0(a) . (3.78)
Substituting (3.35) in (3.43) and using (3.45), we obtain
1c=§Vs=(§o+0(€))(Vo+0(i-Z))

= éo Vo + 0(8) = to + 0(8). (3.79)

53

3.5 Higher Order Approximations of the Perron-Frobenius Eigenvector of NCD

Systems

In some of the applications, an 0(8) approximation may not be good enough to

meet certain specifications; so, a higher order approximation is needed. In [26] and

[29], a procedure was given to get an C(82) approximation of the steady state proba-

bility. In this section we develop a straight forward procedure to obtain an 0(8")

approximation. The only computational effort we need in this procedure is some

additional matrix multiplications which can be done easily in the presence of today’s

VLSI technology by using, for example, an array processor.

It is shown in the appendix that V, can be written as
V,= V,,—e V03 W,(V,(A+8B)W,)‘1V,.
Substituting (3.80) into (3.44), we get
Q,=VOB W, —e VOBW,(V,(A+8B)W,)'1V,B W, .
If we expand ( V,(A + 8 B) W,)'1 in 8, we obtain

( V,(A + e B) W,) ‘1 = ( V,A WJ‘ f: e"[- ( V,B W,)( V,A W,)“]k .
b=0

Substitution of (3.82) into (3.81) yields

Q, = Q,, + 8Q, + ...... + rat-10,, + 08*)
where

Q,=-VOBW,(V,AW ‘1 V,BW,,

j-I
Q,- = - V0 8 W3 v24 W2)" [— < v28 sz V24 W04] v23 W,. j 2 2.

Now, let the (k+1)th order approximation of P, be P, then

Pk=l+8Qo+82Q1+ ....... +8ka_1.

(3.80)

(3.81)

(3.82)

(3.83)

(3.84)

(3.85)

(3.86)

54

Theorem 3.4

If P, is defined by (3.86), then é, is an 0(8‘) approximation of g and 1C), is an

O(8*) approximation of rt, where
gkpk=§kr §k1N=1v (337)
“k = SI: V1 + 1'11: V2 (3-33)

and n, is an 0(8") approximation of 1] obtained by expanding the right hand side of
(3.35) and dropping the 8"“ and higher order terms.

Proof :

From (3.87)
ng-k-‘O, Sh 1N=1

Q,, = Q0 + a Q, + ....... + s“ Q,., . (3.89)

From Lemma 3.2 and Theorem 3.3 it follows that
gt ‘Pk,’ = e,- (3.90)

where e,- is defined in (3.72) and ‘1’,- is obtained from Q, by replacing its ith column

by IN. Now, if we expand Q, in 8, we get
Q, = Q“, + 003") . Hence, ‘1’,- = ‘11,,- + 0(5") . (3.91)
It follows from equations (3.72), (3.90) and (3.91) that
(t - t.) ‘1'...- = 00:"). (3.92)
From (3.91) it is clear that for 8 sufficiently small, ‘1’,- is nonsingular. Thus
§= it + 00:"). (3.93) >
From (3.93) and from the fact that n is a function of i, it can be shown that

rt = i, V, + n, V, + 0(8‘) = 1:, + 0(8") (3.94)

55

3.6 Block Diagonal Decoupling Transformation

From sections 3.3 and 3.4 it is clear that there is a wide class of transformations

which give the same aggregated (first order) matrix, (3.58). Moreover, any order of

approximation of the P-F eigenvector it can be achieved via an expansion of the

exact reduced order system, (3.44). In this section a subclass of the general transfor-

mation matrix discussed previously is given. This transformation has a block diagonal

structure. The idea of the transformation is to specialize the choice of W, in (3.28).

Since W, is block diagonal, we may choose W, to be block diagonal as well. Such a

choice will result in block di

agonal matrices

V, , V, , W, and W, defined in section 3.2 become

r

 

 

fl

Vino 000

0V?)00

 

 

V, and V,. 80 the matrices

V,— 0 0 V?) . , (3.95)
L0 0 0 0 Vs”) _
iv?) 0 0 0 o
0 V3” 0 0
V,— 0 0 V331 . , (3.96)
0 0 0 0 V3”)
P ,1 0 0 0 0 q
1,2 0 0
W, .. 0 1,3 . . and (3.97)

 

 

56

w,= o 0 1V9). . , (3.98)

 

 

where
V1” W1°=1. V10 WS"=0. V5” W1°=0 and V3" WS)=I.,.1. (3.99)

In this transformation, computation of V0 is easier and more efficient. Recall from

(3.57) that

V,,: V,- V,AW,(V,AW,)‘1V,.

In this case V, , V, , ( V, A W7) and (V,A W,) are block diagonal. Each diago-
nal block of the matrix V, A W, is equal to v? A,- we 6 11'1"" "1". This simplifies
the matrix multiplications and inversion involved in computing V0

The transformation given in [11] belongs to this class of transformations. As we dis-
cussed in 3.1, in that transformation V?) is chosen to be the P-F eigenvector of
(A,+1,,,) , i=1,, ...... , N, i.e, 0A,: 0 , V?) 1,,,=1 , while the choice of W?) is the
same as our choice. The matrices V8) and WE” are chosen to satisfy the conditions of

(3.99).

We now propose a specific choice of W,” that renders a particular simple

transformation. This choice proceeds as follows:

ForA,-,i=1,2, ....... ,N, let Moe R’” and W,e R"""""’1 be

WP: 1,,, and V1130: [,0 J (3.100)

,-1

 

hull.

57

For this choice, V?) e R1x "‘ and V50 8 19".“< "" ’1 can be computed to be
Vi” = [1, 0, 0, ..... , 0] (3,101)

and

V,.): '1 ‘ I (3.102)
' ; ni—l

Our choice of transformation (3.100) — (3.102) is much simpler than the one
proposed in [11] because it is independent of the system parameters. This is very
important when we deal with optimal control problem as we will see in chapter 4.

Moreover, the transformation F1 A F can be obtained by inspection.

One more feature of our transformation, (3.100) — (3.102) is that we can relate
it to what we have done in chapter 2, where we keep the fast variables as physical
ones while the slow states are taken as linear combinations of the other states. From
(3.30), E: 1: W, represents the slow states and it is a linear combination of some of

the states, where

8: [8,, 8,, ....... . 8”] and (3.103)
"1 "1+": 1:
§1=§1t,', §2=. 2116i, ........ ,§N=. .2 17C, (3.104)

The " fast " states are equal to n = 1: W,, where

111: 1:2 , 112 = 763 , ...., “"1: “a, +1 , .......... , "M = It" (3.105)

3.7 Computational Considerations

So far we have outlined a method for computing the stationary probability dis-

tribution of an NCDMC, section 3.3, as well as a procedure for obtaining

58

approximations to any desired degree of accuracy, section 3.5. We have also intro-
duced a block diagonal transformation, section 3.6, to be used with this method In
this section we discuss the computational advantages of the proposed method. To put
our discussions into perspective, we discuss also Courtois’ and Vantilborgh’s

methods.
The common motivation for all these methods can be summarized in two points:

1. To reduce the higher order n-dimensional system (3.1) to a small order N-

dimensional system, where N is the number of groups in the system.
2. To overcome the ill-conditioning of the problem

For an 0(a) approximation of the stationary probability distribution, Courtois
[25] succeeded in achieving these two goals. His method requires computing the P-F
eigenvectors of the groups (1,,-t- A9, 1': 1, 2, N. The computation of the P-F eigen-
vectors can be performed either by a direct method or by an iterative method, such as
the methods given in [31 — 34]. From the eigenvectors the matrix $71 of (3.9) and the
approximate aggregated matrix of (3.8) are computed. These two equations are used
to compute the approximate stationary probability of the whole system. Vantilborgh
[29] tried to extend Courtois’ work in order to obtain higher order approximations.
We summarized his procedure for obtaining C(82) approximations in section 3.1.
From that summary we can see that he needs to repeat the first step of the 0(a)
approximation by computing the P-F eigenvectors of the new mauices
(1,,; + Q;),i = 1, 2,....., N. 80, if it is required to compute the O(e*) approximation, the
P-F eigenvectors computation must be done It times. Also in the course of forming
the new matrices QE, so many approximations are done to get rid of near singular
matrices, e.g., (3.14) and (3.17). This process of forming Q; is very involved and the
complexity will grow tremendously when higher than C(82) approximation is

required.

59

We claim that our proposed method, outlined in sections 3.3 to 3.6, is a natural
and a straight-forward extension of Courtois’s work. The 0(8) approximation of sec-
tion 3.4 is equivalent to Courtois’ method. The matrix V0 can be computed exactly as
in Courtois’ work, using the fact V0 = VI. Alternatively, it can be computed using the
closed form expression (3.57), which requires inversion of V2 A W2. When the block
diagonal transformation is used, this matrix is block diagonal. Hence its inversion
reduces to the inversion of N diagonal blocks, of dimension ("i - 1) each. The two
methods for computing V0 are comparable, since computing the P-F eigenvector of a
matrix of dimension n,- is roughly equivalent to the inversion of a matrix of dimen-
sion n,- - 1. The computation of ( V2 A W2)_1 is preferable as we move towards higher

order approximations, as we will see shortly.
To compute the exact stationary probability distribution using the method of
section 3.3, we need to compute V,. Recalling that

-1
Vs = V0 '- 8 V0 B W2 [V2( A + 8 8) W2] V2 (3.80)

we see that the main step, beyond computing V0, is the inversion of

[V2( A + e B) W2]. As we mentioned before, when the block diagonal transformation

is used, the matrix V2 A W; is block diagonal. But [V2( A + e B) W2] is not a block

diagonal matrix even though it is 0(a) close to a block diagonal one. This suggests

that one should exploit the closeness to a block diagonal matrix to simplify the inver-

sion of [V2( A + e B) W2]. This was done in section 3.5 by expanding

-1
[V2( A + e B) W2] as a power series in 8. Although through this expansion pro-

cedure we need only to compute the inverse of the block diagonal matrix V2 A W2,
the amount of associated matrix multiplications grows exponentially as we get to

higher order terms in the expansion. An alternative idea to exploit the closeness of

6O

—1
[V2( A + e B) W2] to a block diagonal matrix is to compute [V2( A + e B) W2] itera-
tively. This idea is developed below.

Letxé VzA W2 and Y2 v23 W2. Then
Z é [V2( A + e B) W2]-1 = [X + eY]" (3.106)
where X is a block diagonal nonsingular matrix.
Multiply (3.106) from the right by (X + 81’) to obtain
Z [X+ eY] = I .
Hence,
Z=X"-eZYX‘13f(Z) (3.107)

where f maps Z into Z.

The solution of this equation can be obtained via successive approximations if f (.) is

a contraction mapping. We have
mo 4%) = a [22-21] L
where L $- YX“ and 21, z, e 2. Thus

llf(Zl)-f(Zy)llSe||Zz—ZlllIILll.

Since A, B, V2 and W; are 0(1), IlLl|=O(1). This implies that ellLll<1 for
sufficiently small 8. Therefore, f is a contraction mapping. Thus equation (3.107) has

a unique solution which gives f (Z) = Z.

The successive approximation for solving the above equation is described as fol-

lows:

o Let20=0

61

o Fork=1,2, ......

2*” =X’1—esz (3.108)
0 As k -> oo , Z" —) [V2( A + a B) W2]_l. Moreover,

uz"-zu=0(e"). (3.109)
To show (3.109), subtract (3.108) from (3.107) to get
z—z"=e [zk-z] L.
Now, fork=0,z—z1 =—eZL=0(e) and for k=1
Z—Zz=—e(Zl—Z)L=EZZL2=0(62)
and by induction we can show that (4.109) is true.

Notice that the only matrix inversion required in this process is computing the
inverse of the block diagonal matrix V2 A W2 which is already computed before

when we compute the matrix V0.
We conclude this section by the following example

Example 3.2 : Consider the example given in [26,29]. Its transition Probability

 

 

matrix is

10.85 0 0.149 0.0009 0 0.00005 0 0.00005‘
0.1 0.65 0.249 0 0.0009 0.00005 0 0.00005
0.1 0.8 0.0996 0.0003 0 0 0.0001 0
0 0.0004 0 0.7 0.2995 0 0.0001 0

P = 0.0005 0 0.0004 0.399 0.6 0.0001 0 0
0 0.0005 0 0 0.00005 0.6 0.2499 0.15
0.00003 0 0.00003 0.00004 0 0.1 0.8 0.0999
_0 0.00005 0 0 0.00005 0.1999 0.25 0.55

62

The matrices A,- for 1' = 1,2 and 3 are chosen as

-0.15 0 0.15 [413 0.3 :l -0.4 0.25 0.15
4
0.1 0.8 -0.9 0.2 0.25 - 0.45

From (3.3) and (3.5), e = 0.001 and the matrix B is

 

'0 0 — 1 0.9 0 0.05 0 0.05 '
0 0 — 1 0 0.9 0.05 0 0.05
0 0 - 0.4 0.3 0 0 0.1 o
0 0.4 0 0 - 0.5 0 0.1 0
3 = 0.5 0 0.4 — 1 0 0.1 0 0
0 0.05 0 0 0.05 0 — 0.1 0
0.03 0 0.03 0.04 0 0 0 - 0.1
_0 0.05 0 0 0.05 - 0.1 0 0

 

We first compute the 0(8) approximation. Then the exact steady state probability is

computed. First, from (3.5 8) the approximate aggregated matrix is

0.9991 10 0.000790 0.000100
P1 = 0.000614 0.999286 0.000100
0.000056 0.000044 0.999900

The steady state probability of (3.69) is £0 = [0.22258 0.27742 0.50000] and
from (3.71),

no = [0.089032, 0.092903, 0.040645, 0.158526, 0.118894, 0.12037, 0.277778, 0.101852] .

The steady state probability computed by Courtois’ method is

1:0 = [0.089029, 0.0929, 0.040644, 0.15853, 0.118897, 0.12037, 0.277777, 0.101852] .

The two methods are equivalent in the case of the first order approximation. The

slight difference in no computed by the two methods is due to numerical errors.

III!

63

The exact aggregated matrix computed as in (3.38) is

0.999109 0.000791 0.000100
P, = 0.000614 0.999286 0.(X)OIOO
0.000056 0.000044 0.999900

The steady state state probability for this matrix is equal to
§= [0.222529 0.277471 0.500000]
and the exact steady state probability as defined in (3.43) equals

it = [0.089283, 0.092758. 0.040488. 0.158533. 0.118938, 0.120416. 0.277795. 0.101789] .
If we compute it by using (3.1), we obtain

1: = [0.089069, 0.092959, 0.040509. 0.158529. 0.118935. 0.120385, 0.277795, 0.101819] .

Notice here in this example that the difference between the exact steady state
probability evaluated by the proposed method and by using the direct method, i.e., by
using (3.1) is very small. That is because the order of this example is eight only.
However, the condition number of the matrix A +8 B, with one columns replaced by
the vector 1,, is 16554 compared with 15.2345, the condition number of Q, with one
column replaced by the vector 1”. This shows that the condition number of the

aggregated matrix is much better than of the full order matrix.

It.

CHAPTER 4

OPTIMAL AND NEAR OPTIMAL CONTROL FOR
NEARLY CONIPLETELY DECOMPOSABLE
CONTROLLED MARKOV CHAINS

4.1 Problem Statement

For the finite state Markov chain, the transition probability is defined by
Pi](u)=prob{Xk+l =j/X= i, u,-=u}, i,j= 1, 2, ..... , n and u 6 U0).

So, P004) is the probability that the next state is j given that the current state is z’ and

the control applied is u e U(r), where the set U0) is a compact set for all i. If we

. . . . T
assume that the control u rs statronary, then the statronary polrcy u = [1:1, uz, ..... , un]

can take any value in U 6 U0) x U(2) x ...... x U(n). The transition probability matrix

P(u) is given by

.1’110‘1) P12(u1) - - - P111011)-
1021042) P2200.) - - « 1’24“?)

P(u)= . . . . . . . (4.1)
PM“) P .1204") - - - PM“)

 

 

Notice here that the matrix has the same properties as the uncontrolled Markov chain,
i.e., for each u, the row sum is 1 and its elements are nonnegative. Note also that the
ith row of P(u) depends only on “i-

In this chapter we assume that P(u) takes the nearly completely decomposable

form, i.e., P(u) can be written as

P(u) = 1,, + A(u) + 8 B(u) = 1,, + Q(u) (4.2)

65

where A(u), B(u), e and Q(u) are defined by (3.3) — (3.5). We assume that for all
u,- 6 U0), the matrix P(u) has a single ergodic class. Moreover, for each u associated

with the groups, i= 1, 2, ...... , N, the stochastic matrix 1.6+ A,(u) has a single ergodic
class. This assumption is equivalent to the following:

For every u e U there are unique stationary probability vectors rt(u) and v,(u) such

that
1c(u) Q(u) = 0, 1t(u) 1,, = 1 (4.3)
and
v,(u) A,(u) = 0 , v,- 1,”: 1, i= 1,...., N . (4.4)
The problem of interest which we are going to deal with in this section is: what
is the policy u e U that minimizes the average cost per stage

1
M+l

 

J(u) = lim
M-D on

M
E 2 f (x1. . 1101)) (4.5)
e-o

where f (i, u,-) is the instantaneous cost. This cost is assumed to be stationary and con-

tinuous in u,- e U(i) for all 1'.

Under the ergodicity assumption [36, 40],

 

_. 1 M " _. _ .
1(u)_Ml-)lm~M+lEb§o§Pr[Xk—l]f(x.’ul)
1

 

= lim E i p00.) P"(u) flu) (4.6)

M-)°°M+l ‘30

where po(u) is the initial probability distribution and

F'f(1’ “1).
f (2. “9
f (u) =

 

 

{<4 a».

66

But,

lim
M—Hn M+1

 

M
E 2; p00.) P"(u) = p004) 1,,1r(u) = 1:01).
bao

Substituting this into (4.6) we obtain
J ( u ) = Mu) f (u) . (47)

Notice here that J( u) is independent of the initial probability po(u); so, the optimal
cost per stage is equal for all states. Now, the minimum average cost per stage turns

out to be:

1’04) = m{ «(10f (14)} (4.8)

where 1c(u) is defined by (4.3). The existence of u e U that minimizes (4.8) is

guaranteed because of the compactness of U [36].

The minimum cost of (4.8) can be obtained by solving the two equations (4.3) and
(4.7) for all u e U; then choosing the policy that gives the minimum average cost.
This procedure is possible if the order of the matrix Q(u) is small and the number of
policies among which we choose is finite and small. Unfortunately, this is not the
case for most practical problems. Equation (4.3) usually represents a very high order
system. In addition to that, it is very ill-conditioned as we discussed in chapter 3.
These difficulties enforced the researchers [36 — 43] to look for an alternative
approach to tackle this problem. In section 4.2, four numerical methods are presented
to solve the minimum average cost problem in its general form, and the fifth method

deals specifically with the class of nearly-completely decomposable systems.

In section 4.3 we will make use of the structure of Q(u) to overcome the high
dimension and the ill-conditioning problems by using the two-time scale property to
decompose the system into small and well-conditioned subsystems. This will help us

solve the optimal control problem of (4.8) which is the topic of section 4.4.

67

In section 4.5, we consider the solution of (4.8) as 8-) 0. This results in a
reduced problem for which a near-optimal policy is obtained. Finally, in section 4.6,

an example is given to illusu’ate the optimal and the near-optimal solutions.

4.2 Algorithms for Computing the Optimal Control Policy
For c e R" and l. e R, equation (4.7) can be written as [36, 39]

2. 1,, = Q(u) c + f (u) (4.9)
where we have n equations in n + 1 unknowns. The following lemma is recalled from
[36, 39].

Lemma 4.1 Under the above assumption about Q(u) and f (u) and for u e U,
i. If (2,, c) is a solution to (4.9), then A = 1(a).

ii. If (71. c) is a solution to (4.9), then (2., c + 5 1,,) is also a solution for any scalar

5.

iii. A solution always exists.

In the literature there are several algorithms to compute the stationary policy

that minimizes the cost (4.9). These methods are described briefly in this section.

1. Policy Iteration Method [41]

This method is composed of two steps: the value determination and the policy

improvement steps.
a. Value Determination Operation

For a given admissible stationary policy u" at the kth iteration, the average cost
according to Lemma 4.1 is equal to 1", i.e., J(u") = 2", where the superscript denotes

the iteration number. The corresponding unknown dual variables 1:" e R" are given

68
by
xk=f(i,u§)+f;q,j(u§)c"(j) ,i=1,2, ..... ,n
j-l

or in the matrix form

71"l,,=f(u")+Q(u")c". (4.10)

Recall from Lemma 4.1 that if (7).", c") is a solution to (4.10), then so is (7.", c" + 8 1,,)
for every constant 8. This is equivalent to setting one of the n elements of c" to zero.

Therefore, equation (4.10) with
c"(s)=0, (4.11)
where s is any state, has a unique solution; see [40, 41].

b. Policy Improvement

The next stationary policy u"+1 is determined by minimizing
H,- (um, e ) = min {f(i, u.) + 2; gym) c"(j)},i=1,2, ..... , n. (4.12)
"t5 ”(0 j..1

In the matrix form this minimization can be written as

”04"“, e ) == mirtr]{f( u ) + Q(u) 0"}. (4.13)

Because in the controlled Markov chain every element in a row is dependent on one
particular policy 14,, the minimization (4.12) is a point wise minimization. So we can
find the minimum policy for each row separately. Then the iteration will proceed
with this new policy to the value determination step and the cycle is repeated. If u"
gives the minimum cost in (4.12), we choose a?“ = u" even if there are other controls

giving the minimum beside uf‘ [40]. If the policies on two successive iterations are

equal, then the optimal policy is reached, i.e., f (u) = 3." and the iteration is stopped.

69

In this algorithm, it is assumed that the system is a complete ergodic Markov
processes, i.e., through the iterations, we assume that every stationary policy gives a
Markov chain with a single ergodic class [41]. It was proved in [41] that in this
algorithm the policy evaluated at each iteration gives a lower cost than the previous
one. This guarantees that the algorithm converges. Moreover, if the number of poli-

cies is finite, the iteration converge in a finite number of steps.

The implementation of the algorithm is based on being able to solve the linear

equation which is defined by (4.10) and (4.11).

The following methods do not require a matrix inversion of order n in each

iteration; but they solve the optimization problem iteratively.

2. Successive Approximations Method [42]
This algorithm is summarized as follows:

Pick any initial value for co; for k = 1, 2, ....... , perform the minimization
do) = min {f(i, 14,-) + 2 4,, (11,) c"‘1(/) + 64(1)} ,i = 1, 2, ..... , n. (4.14)
“t5 (10') F1

Set 2." = C" (s) and c" (z) = C" (r) — 2." for 1': 1, 2, ...... n. The iteration is repeated until

the algorithm converges.
It is shown in [42] that J(u‘) = klim 71". This method does not require solving n linear

equations, however, the rate of convergence is low [37].

3. Gauss—Seidel Method [43]

This method is an iterative method as the successive approximation one except
that at every iteration the new C" (r) is used in computing C" (1) for j > i, i.e., equation

(4.14) is replaced by

70
i-l n
0‘0) = min {ma u.) + 2; q,,-(u,-)C*'1(,) + c"‘1(i)+ z q,,(u,-) $10)}, i= l,...,n . (4.15)
“15 (1(1) 1‘31 jag

It has been found in [40] that the convergence of this method is better than the suc-
cessive approximations method. The convergence of the Gauss-Seidel method is, in

general, guaranteed if the matrix Q(u) is diagonally dominant [45].

4. Varaiya’s Algorithm [36, 37]
This algorithm can be described in the following steps:
0 Given the tolerance 8 > 0, choose co.

9 Fork=1,2, ...... ,h(c")eR"isdefinedby
h(¢5=m{f(u)+ Q(u) 9“}- (4.16)

o For the minimum of (4.16) determine

m0 = max h.- 0") and 110*) = min h.- (.4) (4.17)
where i is the ith element of the vector )1 (c"), and

1:05:71!- fimc‘) and K(c*)=r(c*)_r(c5 1,,. (4.18)

5:1
. If Re") — a (c") < 8, stop and set u‘ = u" and 10.") = 11(6).
Else
0"” = c" + A K(c") (4.19)

where the step size A is a small positive constant.

0 The procedure continues till convergence is achieved.

71

5. Nested Optimization Method

The nearly completely decomposable Markov chain of (4.2) was considered in
[30, 38] where a near optimal average cost is computed. This method employed the

transformation (3.22) to transform equation (4.9) into the aggregated form

J 1,,, = [V,(u) — e F (u) 192] B(u)W1c‘1 + f (u) (4.20)
where
F (u) = 1?, (u) B(u)W2(u) [9g A(u) + e B(u)) 0720.)] ‘1 (4.21)
and
f“ (u) = [171 (u) - e F 00172] f (u) . (4.22)

If we try to apply Varaiya’s algorithm to (4.20) to obtain the minimum average cost,

we obtain
1: (6’0 = .f‘éi‘b{[‘71(“) - e 1700172] MW? #00} (4.23)

where in this case as k—) .., h (610—) J 1,,. But, equation (4.23) is nonlinear and the
transformation of (3.22) depends on u. This makes the minimization of (4.23) very

complicated task whether we use Varaiya’s algorithm or any other known algorithm.

Now, if we let 6 equal zero, (4.23) becomes

’2 (a’f) =mig {17104) B(u)W1 51‘3“ V104) f (14)} =mig {9104) [B(u)W1 61+)" 00]}. (424)

Notice that each row of B(u)W1(’:"f and f (u) is a function of only one control

u,- ,i = 1, 2,....,n. Now, partition u and f (u) as follows

72

 

 

 

 

"u“ "f(u‘)"
u2 f (uz)

u = . and f (u) = . (4.25)
1“". .f (”N )1

where u" e R"‘ e Ui is the control associated with the fast subsystem A), i = 1, N

and f (u‘) is the corresponding instantaneous cost. Now, let g"(u) be equal to

 

 

' '1
g" (u‘)
g" (uz)
g’ru) = 3(a) W. 61‘ +f(u) = . . (4.26)
.8" 04”)
Hence equation (4.24) becomes
h (fl) = grim» g" (1‘)} (4.27)

where g" (u) can be considered as the instantaneous cost for (4.27).

Recall that 17104) is block diagonal; therefore, equation (4.24) can be written as
It, (61‘) = min{v,(u) g" (18)}. 1': 1, 2...... N. (4.28)
n" e U‘

Equation (4.28) shows that for a fixed dual variable 61‘, the minimization problem of
(4.27) can be separated into N different optimization problems for the fast subsystems

A,(u) , i = 1, 2,...., N. Varaiya’s algorithm can now be applied as follows:
For a fixed at and for each i= 1, 2,..., n, the minimum average cost can be computed

by
h,(é’f,di)= mirlrf{A,(u‘)df+g*(u‘)},i=1,2, ..... ,N. (4.29)
u‘e

If the algorithm converges, then as j —9 oo , h,- ( é’f , dj) -> h,- ( c’f ) for i = 1, 2...... N.

 

73

Notice here that the minimum cost for (4.27) is computed at 61‘. Now é’f is updated by
the same procedure as described in method 4 to obtain 6’1"“. This updated dual vari-

able is used to update the cost
g"+1 (u) = 800W, 6’1‘“+f(u).

The minimization of (4.28) is repeated to obtain h (é’fH ) from (4.27) and the cycle
continues till the algorithm converges. Notice that in this algorithm the minimization
of the aggregated system (4.27) is replaced by the minimization of N fast subsystems
defined by (4.28) and (4.29). Thus, the optimization is nested. So, this algorithm
decomposes the optimization problem into (N + 1) subproblems: one slow time-scale

problem and N fast time-scale problems.

Let us finally note that, because of the complexity of (4.23), the optimal solution
of the average cost is very difficult. This makes seeking a near-optimal rather than

optimal solution a must if we decide to choose this method.

In section 4.3 we will see that the optimal solution for the average cost per
stage is possible by using a transformation which does not depend on the control u

as the one proposed in section 3.6.

A comparison between the first, second and fourth methods is given in [37],
where the rate of convergence, number of operations and the CPU time are compared
for the three methods. This comparison shows that the first method has an advantage
over the other two methods in the CPU time and the rate of convergence. However,
this method does not work well when the matrix Q(u) is near-singular. This led the

authors of [37] to favor the fourth method.

In the following section we will treat the near-singularity and the high order

problems in order to get benefit of the advantages of the first method.

it

74

4.3 Decomposition of the Average Cost for a Fixed Policy
We have seen in the previous section that the cost J(u) = 7. is defined by the fol-
lowing two equations
1 l. = Q04) 0 +f (u) (4.9)
and
c(s)=0. (4.30)
Due to the high dimensions of Q(u) and the ill-conditioning problem, solving the
above two equations is impractical.

In this section the cost 1(a) is decomposed into fast and aggregate components
by using the two time scale property of nearly decomposable Markov chains and the
transformation developed in Chapter 3. As discussed in Chapter 3, the transformation
can be chosen to independent of the system parameters. A specific example is given
in section 3.6. Throughout this chapter, whenever we refer to the transformation of
Chapter 3, it is implied that the transformation is independent of the system parame-
ters; hence, independent of u. The decomposition can be done as follows: From

equation (4.9)
I“). 1,,: I“1 Q(u)I‘I‘" c+r-‘f(u). (4.31)
substitute (3.28) and (3.29) into (4.31) to obtain

A 1N 8 V1 800 W1 V1(A(u) +8800) w, Vlc v1 f(u)
0 = eVzB(u) Wl V2(A(u) +eB(u)) W2 vzc + V2f(u) (4.32)

VIC = _:' Cl and V2C = C2 . (4.33)

Therefore (4.32) becomes

75

A. 1N = V1 B(U) W1 C1 + V1( 1404) + 8 800) W2 C2 + V1f(U) (4.34)
0 = V2 B(u) W1 cl + V2( A(u) + 8 B(u)) W2 (:2 + V2 f (u) . (4.35)
It was shown in Chapter 3 that V2( A(u) + e B(u)) W2 is nonsingular; thus,
-1
C2 = " [V2( A(u) + 8 B(u)) W2] [V2 B(U) W1 Cl 4‘ V2f(ll)] . (4.36)

Substituting this into (4.34) yields

7» 1N = V. 00 B(u) W1 ct + V.- (u)f (u) = (2.-(1061 +13 04) (4.37)
where
v, (u) = Vll -- V,( A(u) + e B(u)) W2[ V2( A(u) + e B(u)» W2] '1 V2 , (4.38)
Q. (u) = V..- (u) B(u) W1 (4.39)
and
f. (u) = V. (u)f (u) . (4.40)

Now, the linear system of n equations (4.9) is reduced to a smaller order system of
order N, where N is defined in Chapter 3. Similarly, if we apply this transformation

to equation (4.7), we obtain

1(a) = “(101‘ F'1f(u)= [5 (u) 11 (14)] [:‘;:::J (441)
2
where g (u) = it (u) W1 and n (u) = 1: (u) W2. Using (3.35), we obtain
11 (u) = — t 00 V.( 400 + e 800) W2 [Va 400 + e 300) will (4.42)
substitution of (4.42) into (4.41), yields
1(a) = § (14)}:- (u) (4.43)

where i (u) is defined by

76

g 0‘) Q: (14) = 0 (444)
Notice here that the properties of V, and Q, which are defined in Chapter 3 still hold
for V, (u) and Q, (u).

From Chapter 3 we know that Q, (u) has some of the properties of Q(u) such
as the row sum of Q, (u) equals zero and the uniqueness of the zero eigenvalue, but
( 1N + s Q, (u)) may not be a stochastic matrix. Therefore, the existence of solution
of (4.37) does not follow from standard properties of Markov chains. The following

two theorems show the existence and uniqueness of the solution of (4.37).
Theorem 4.1

For a fixed policy u, consider the N linear equations with (N + 1) unknowns
A e R and c, e R" of equation (4.37). Then

i. If (A, Cl) is a solution to (4.37), then 2. = 1(a).

ii. If 0., Cl) is a solution to (4.37), then (2., C1 + 5 1N) is also a solution for every

5.
iii. A solution always exists.
Proof :

i. Multiply (4.37) on the left by § (u) to obtain
7» = 5 (14) Q.- (u) 01 + 5 (10f. (u) = § (10f. (u) = 1(a)

ii. The proof of this part comes from the fact that the row sum of Q, (u) equals
zero; therefore 5 Q, (u) IN = 0. That is to say 6 IN 6 N ( Q, (u) ), where

N ( Q, (u) ) is the null space of Q, (u).

iii. Recall that A x = d has a solution if d e R (A ), where, R (A ) is the range space

of A.

 

FE

 

77
Now, from equations (4.43) and (4,45)
a (u) [100 1N -f. 00] = o .

So, i (u) is orthogonal to [1(a) 1N - f, (u)]. This implies that
N (23.00) .L [1001. 1.00] . (4.46)

But
N ($00) = R ( Q. (“Dc (4-47)
where R ( Q, (u))‘ is the orthogonal complement of the range space R ( Q, (14)). Sub-

stituting (4.47) into (4.46) yields

R ( Q. 04»: -I. [1001. -f. 00] . (4.48)

Therefore, the columns of [1(a) 11v - f, (14)] are in the range space of Q, (u) which

means that there exist a vector c, e R" such that
Q. (u) or = 1(a) In -f. (u) .
Thus, Cl is a solution.

The above theorem proves the existence of a solution. The following theorem

proves the uniqueness of the solution under certain conditions
Theorem 4.2
For c1, f, (u) e R”, let the ith element of c, be equal to zero, then
A 1N = Q. (u) 01 +f. (u)
with
cl (1') = 0 (4.49)

has a unique solution.

78

Proof :

To prove that (4.37) and (4.49) have a unique solution, it is enough to prove

that
‘I‘,,- (u) = [q] (u), q} (u), ...... ,qu (u) , — 1,.,,q;'.+1 (u), ....... .qfi" 0.)]

where, qf, (u) is the ith column of Q, (u), is nonsingular. The proof of this is a direct
result of Lemma 3.2.

So, by using the two-time scale property of the full order system, we could reduce

the order of the system from n to N and the system which we need to solve now has

a unique solution and is defined by (4.37) and (4.49).

At this stage we have overcome the drawback of the policy iteration method
which requires solving n, near-singular linear equations at each iteration. Instead, we

solve N, well-conditioned linear equations.

In the next section, the policy iteration method is employed to solve the optimal

control problem and in the following section it will be used for the near-optimal case.

4.4 Optimal Policies for the Average Cost Problem

In this section the policy iteration method [41] is used to minimize the cost of

the aggregated system, i.e.,
A IN = mirlrj {Q, (u) c, + f, (u)} (4.50)
“E

The value determination method for aggregated form is used to compute 2." and ch
for a fixed policy u" at the kth iteration. This can be done by solving N linear equa-

tions in N unknowns as follows

3" 1~= Q.( u") of +12.( u") (4.51)

79

c’,‘( t) = 0 (4.52)
where t is any state.

For a reason which will be discussed shortly, in addition to of, we need to com-

pute c". This can be done by computing c’z‘ first; by using (4.36)
-1
c’i=— [V2000 ”30!» W2] [V.B(u*) W. .§+ mm]; (4.53)
using equation (4.33), the dual variable c" can then be written as

Ck: [W1 W2] V161 = -1- W1 C’f '1' W2 C5 (4.54)
V2Ck 8

Once we compute 7." and c" at the kth iteration, we can go to the policy
improvement step to compute the next policy. Now, for the full order system, it was
shown in section 4.2 that the minimization step is a point wise minimization, i.e.,
each row can be minimized separately. But, if we try to minimize the aggregated sys-
tem (4.50), we will loose this property. Moreover, we need to disaggregate the
reduced policies to the original ones in order to choose the minimum policies among
the given physical ones. Thus, this step will be applied to the full order system;

therefore, the minimization step (4.13) become
1104"“, e) = grim.) + ( A(u) + e B(u)) c"}. (4.55)
Substitute (4.54) into (4.55) and from the fact that A(u) W, = 0, we obtain
H (u"+1, e) = 321% f (u) + B(u) W, of + ( A(u) + e B(u)) W2 c5}. (4.56)

Since W, and W2 are independent of u, we still have each row in (4.56) dependent
on one control u,- e U(t). So, the point wise minimization can still be performed in

this step.

80

It is important to notice that because the transformation proposed in Chapter 3 is
independent on u, this step could be done without violating the point wise minimiza-

tion. This step would not be possible if W2(u) of [11] is used instead of W2.
We conclude this section by summarizing the above algorithm.
1. Value Determination Step

0 For an admissible stationary policy u" at the kth iteration, compute the aggre-
gated matrix Q,(u") and the aggregated instantaneous cost flu") as defined in

(4.39) and (4.40), respectively.

0 Solve for 7." and of which are uniquely determined by (4.51) and (4.52). Set
J(u") = 2." and by (4.53), compute c5. Now, pass the dual variables of and c5

to the policy improvement step.
2. Policy Improvement Step

To compute the next policy, i.e., 14"“, use equation (4.56) with the new of and
c’fi computed by the first step. This minimization, as we discussed previously, is per-
formed one row at a time. The new policy u"+1 is sent to the value determination step
to obtain new values for 2"” , of” and 0’5“ and the iteration continues until the the

algorithm converges; for which 14"“ = u" and J(u‘) = 8.".

The convergence of this algorithm is a by-product of the theorem given in [41]
because the similarity transformation we applied to the system will not change its

convergence.

We conclude this section by noting that with this algorithm, instead of solving n
ill-conditioned linear equations in each iteration, we only solve N well-conditioned

linear equations, where N < n.

fir;

81

4.5 Near Optimal Policies for the Average Cost Problem

The most expensive part in the algorithm proposed in section 4.4 is computing
V,( u) which is defined by (4.38). In this section we wish to take advantage of the
singularly perturbed form of the nearly decomposable Markov chain to obtain a
near-optimal average cost per stage by approximating V,( u ). By letting a = 0 in

(4.34) and (4.35) we obtain
7‘0 1N = V1 B(u) W16 10 + V1 A(u) W2¢20 + Vrf (u) (457)
0 = V2 B(u) W1C10 + V2 1404) W2C20 + V2f(u) . (4.58)
Since V2 A(u) W2 is nonsingular,
-r
020 = — [V2400 W.] [V. B(u) W. 010 + w 00] (4.59)

Substituting (4.59) into (4.57) yields

7\0 1N = [V1 " V1 A(u) W2( V2 404) W2)"1 V2] [B(u) W1 0 10 +f 01)]

g V0 (u) B(u) W1 0 10 + V0 (u)f (u) 3 Q0 (u) C 10 +f0 (u) (4.60)
where
-1
V0 (1‘) = V1 - V1 :40!) W2[ V2 A(U) W2] V2 , (4.61)
Q0 (14) = Vo (u) 800 W1 (4.62)
and
fo (u) = Vo (u)f (u) . (4.63)

Equations (4.60) — (4.63) are the unperturbed form of equations (4.37) — (4.40).

By the same procedure, equations (4.43) — (4.45) can be written as

1004) = 50 (u) f0 (u) (4-64)

82

where 5,0 (u) is defined by
£0 (14) Q0 0!) = 0 . in (u) > 0 (4.65)

From Theorem 3.2 of Chapter 3 we notice that the results of Theorems 4.1 and 4.2

hold for 2.0 , 010 and that
)‘0 IN = Q6 (u) 010 +fo (u) (4.67)
with
(:10( i ) = 0 , (4.68)
where i is any state, has a unique solution.

In the following theorem we want to show in what sense is the optimal solution

of (4.64) a near optimal solution of the original problem (4.43).
Theorem 4.3

Let f (u) be a uniformly bounded function of u for all u e U. If u; is a policy

that minimizes (4.64), i.e.,

Jo(u6) = flit}! 1004) = gig/{Eu (u)fo (14)} (4.69)
and if u‘ is a policy that minimizes (4.43), i.e.,

104') = girl), 1(a) = git}, {5 (11)}: (14)} . (4.70)

Then

10.5 ) = 104‘) + 0( e ) (4.71)

 

83

Proof :

From Theorem 3.3 and equations (3.80) and (3.83)

E (u) = §o (u) + 0( 8 ) . (4.72)
V, (u) = Vo (u) + 0( e) (4.73)

and
Q. (u) = Q0 (10+ 0( 8 ) . (4-74)

Substitute (4.72) —— (4.74) into (4.43) to obtain

J(u) = [Q(u) + 0( 8)]

x [Vo(u)+0(e)]f(u)=§o(u)fo(u)+0(e). (4.75)

Therefore,

J(ut') ) =J(u6) + 0( e) = mgr}, {16(14)}+ 0( a)

= 31% {£000 Vo (u)f(U) }+ 0( 8)

-ueU

min{(§(u)+0(e))(V,(u)+0(e))f(u)}+0(e).
Because If (u) I s K for all u, we have

1045 ) = flit}, {§ (14) V. (u)f(u)}+ 0( 8)

= min] {J(u) }+ 0( e ) = J(u') + 0( e ) (4.76)
are

From the previous discussion we found that letting e = 0 gives us a near

optimal solution. If the near optimal solution satisfies our needs, then the computation

 

84

of the value determination step will be much easier because instead of computing
Q, (u), we need to compute Q0 (u). As it was shown in section 3.7, we can use the
block diagonal transformation to simplify the computation of Q0 (14). If we need
higher order approximation of Q, (u), we can use the technique used in sections (3.5)
and (3.7). The higher order approximation of Q, (u) gives a higher order approxima-
tion to J(u'). Computing the near optimal average cost per stage 1004) is performed
by the same algorithm presented in the previous section. In this case the two steps of

the algorithm are modified as follows
1. Value determination Step

For the kth iteration and for a given stationary policy u", compute Q0 (u) and
f0 (u) as defined in (4.62) and (4.63). Solve for 25 and cfo by using equations (4.67)

and (4.68) and use equation (4.59) to compute c’fo .
2. Policy Improvement Step

Use the dual variables cfo and c5, computed from the value determination step

to compute the next policy 14"“ as follows
H (um, 0) = min] { f (u) + B(u) w, c’fo + A(u) W2 .50 }. (4.77)
“6

Notice here that the minimization step is independent of e. The new policy 11"“ is
passed to the value determination step and the cycle continues until the convergence

is achieved.

In the following theorem we want to prove that using this modified version, the
algorithm gives a policy that attains the minimum average cost per stage for the

problem described by (4.67) and (4.68).

85

Theorem 4.4

The algorithm described above produces a monotonically decreasing cost and
converges to a policy that gives the minimum cost per stage to equations (4.67) and

(4.68).
proof :

Suppose a policy f is chosen over the policy u by the policy Improvement rou-
tine (4.77). Let X0 , 510 and 6'20 be the corresponding average cost and dual variables
computed by (4.67), (4.68) and (4.59). Then (4.57) can be written as

X0=V1 3(3) W1 Z"10+V14(F) W250+V1f(i) (4.78)
Now, if it u, we want to prove that X0 5 2.0.
Because the policy improvement step (4.77) chose '17 over u,
f( F) + A( 7) W2 6'20 + B (7) W1 010 Sf(u) + A(u) W2 C20 + B(u)W1010o
This can be written as
Y=f(u)-f(17)+(A(u)- A( 17)) W2 620 + ( B(u) - 3( 7)) W1 010 Z 0 (4.79)
where y e R".

The value determination step for i and u is given respectively, by (4.78) and (4.57).
Subtract (4.78) from (4.57) to obtain

(7‘0 ‘ X0) 1N = V1 B(u) W1 6‘10 - V13( 17) W1 E10 + V1 A(u) W2 6‘20
—V1A(17)W2?20+V1f(u)—V,f(17).(4.80)
Substitution of (4.79) into (4.80) yields

(40"):0) 1N=V1[Y+A(F)W2020+B(17)W1Clo]

—v,3(a)w,r,o—V,A(17)W,r,o.

86

This equation can be written as
0.0—2,) 1,,: V1 7+ VI A( 17) W2( c20- 6‘20 ) + Vl B( E) W,( cm- 5,0). (4.81)
Let 23 = 20 — 20, c‘fo = 010 — 3'10 and 0‘20 = cm — 'C-zo ; therefore equation (4.81) becomes
23 1,,, = V,y + V, A( F) ”V2020 + V, B( 7) W140. . (4.82)
If we apply the same procedure to (4.58), we obtain
0 = V2 B( F) ch‘fo + V2 A( Ti) ch‘zo + V2 7 (4.83)
V2 A( F) W2 is nonsingular; hence,
‘30 = - [V2 AC 3') W2]-l [V2 3(17)W1€io + V2 7] - (4.84)
Substituting this in (4.82) yields
71-3 1N = Q0( 17) Cio 4' V0( 17)7. (485)
Notice here that equation (4.85) takes the same form of (4.60). From (4.64) and
Theorem 4.1, the solution of (4.60) satisfies
20 = 104) = 50 ((4)18 04) = in 04) V0 (u)f(u)
where (:0 (u) is defined by equations (4.65) and (4.66). Similarly, 26 can be written
as
26=§<E>Vo(17)1(. (4.86)

Because Q0 ( E) has the single ergodic class property, the elements of E ( F) are
positive. From Theorem 3.1, part v and equation (4.79), we see that V,, ( F) and y

are nonnegative, therefore,

23>0 => 2022,).

Therefore, the average cost is monotonically decreasing. Since it is bounded from

below, it is convergent.

87

The second step of the proof is to show that it converges to a policy which
gives the lowest average cost. Suppose it does not, i.e., let us assume that 20 < 20,
i.e., 26 > 0 and suppose that the policy improvement step (4.77) converges to u rather
than 17. If this the case, then equation (4.79) yields 7 S 0 which means by (4.86) that
23 s 0. But this contradicts the assumption that 2,, < 20; therefore, the policy will

converge to “12'. This completes the proof.

4.6 Example : Minimizing the Average Cost Per Stage

In this section, we consider the example given in [30]. The algorithm presented
in sections 4.4 and 4.5 are applied to this example to compute the optimal and near-
optimal policies and the average cost per stage corresponding to these policies. Also
this example is solved by the policy iteration method [41] applied to the full order
system. The control problem in this example can be visualized as one of maintenance
scheduling. Similar to this problerrr, hydro scheduling with multiple time scales, was

considered in [44]. The probability transition matrix for this example is given by

.45 .45 0 .05 .05 0

 

 

0 0 0
.27 .36 .27 .03 .04 .03 0 0 0
.0 .72 0 0 .08 .02 0 0 0
.5u .5u .3u .45 - .5u .45 - .5u 0 .05 .05 0
.3u .4u .2u .27 - .3u .36 — .4u 0.27 - .3u .03 .04 .03
0 .8u 0 O .72 - .8u .18 - .2u 0 .08 .02
0 O 0 .5u .5u 0 .5 - .5u .5 - .5u 0
O O 0 .3u .4u .3u .3 - .3u .4 — .4u .3 - .3u
_ 0 0 0 O .8u .2u 0 .8 - .8u .2 - .Zud

The state variables are defined by the variable G and D as follow:

88

 

 

state G D
x, 2 0
x2 2 1
x3 2 2
x4 1 0
x5 1 1
x6 1 2
x7 0 0
x3 0 1
x9 0 2

 

 

 

 

 

 

where G is the number of power generating units available and D is the demand in
terms of generating units needed. The control variable u e [ 0.02 , 0.2 ]. The problem
is to find the policy u (G , D) that minimizes the average cost per stage (4.5) if the

instantaneous cost is

2
f(G. D. u(G. D» = [( D — or] + K(u(G. D»2

where K = 30 and (D — G)" = max ((6 - D), 0).

From (4.2), for 8 = 0.2

 

 

 

— 0.5 0.5 0 0 0 0 0 0 0
0.3 - 0.6 0.3 0 0 0 0 0 0
O 0.8 - 8 0 0 0 O 0 0
0 0 0 -— 0.5 0.5 0 O 0 0
A(u) = 0 0 0 0.3 - 0.6 0 3 O 0 0
0 0 0 0 0.8 - 8 0 0 0
0 0 O 0 O 0 - 0.5 0.5 0
0 0 0 0 0 0 0.3 — 0.6 0.3
L 0 0 0 0 0 0 0 0.8 - .8‘
and B(u) equals
- .25 — .25 0 .25 .25 0 0 0 0
— .15 - .2 — .15 .15 .2 .15 0 O 0
0 — .4 - .1 0 .4 .1 0 O 0
2.514 2.511 0 — .25 — 2.5a — .25 — 2.5u 0 .25 .25 0
1514 2a 1.5u - .15 —1.5u - .2 — 2u — .15 - 1.5u .15 .2 .15
0 4a u 0 -.4-4u -.1—u 0 .4 .1
O 0 0 2.514 2.5a O - 2.5a —2.5u O
0 0 0 1.5u 2u 1.5u - 1.5a -2u —1.5u
0 0 O

 

0 4a u 0 —4u -u

The optimal policy u" and the minimum average cost per stage are computed for
different values of e by using the algorithm developed in section 4.4. These results

are given in table 4.1.

90

Table 4.1 The optimal policies and the average cost for different values of 8
evaluated by the algorithm of section 4.4

 

8 The minimum policy a. J (14.)

 

0.2 .02 .02 .02 .10513 .11178 .11213 .12740 .14478 .15023 0.6709471

 

0.1 .02 .02 .02 .10761 .11087 .11102 .13376 .14287 .14590 0.6714597

 

0.05 .02 .02 .02 .10880 .11041 .11048 .13716 .14183 .14343 0.6716016

 

10"2 .02 .02 .02 .10972 .11004 .11005 .14000 .14096 .14129 0.6716502

 

10’3 .02 .02 .02 .10992 .10995 .10995 .14066 .14076 .14079 0.6716522

 

10—4 .02 .02 .02 .10994 .10994 .10994 .14073 .14074 .14074 0.6716520

 

10—5 .02 .02 .02 .10994 .10994 .10994 .14073 .14073 .14073 0.6716521

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

For all the values of e the algorithm converges in 5 to 7 iterations. If we let 8 = 0,
the near optimal policy 14,", and the average cost per stage are computed by the

method given in section 4.5. The following table gives u; and J ( us)

 

91

Table 4.2 The near optimal policy and the corresponding average cost

 

0.02

0.02

The 0.02
minimum 0.10994
policy 0.10994
0. 10994
0. 14073
0. 14073
0. 1407 3

5.

 

J( .43) 0.6716522

 

 

 

 

Notice here that u; and J (143) are independent of 6. Also as e —-) 0, u" —-> us.

Finally, this example is solved by the policy iteration method [41] which is

described in section 4.2 and the optimal policies for different 8 is given in table 4.3

92

Table 4.3 The optimal policy and the minimum average cost evaluated by the

policy iteration method [41] applied to the full order system.

 

8 The minimum policy u" J (ui)

 

0.2 .02 .02 .02 .10513 .11178 .11213 .12740 .14478 .1502?) 0.6709469

 

0.1 .02 .02 .02 .10761 .11087 .11102 .13376 .14287 .14590 0.6714593

 

0.05 .02 .02 .02 .10880 .11041 .11048 .13716 .14183 .14343 0.6716016

 

10-2 .02 .02 .02 .10972 .11004 .11005 .14000 .14096 .14129 0.6716472

 

10.3 .02 .02 .02 .10993 .10996 .10996 .14066 .14076 .14079 0.6716321

 

10"4 .02 .02 .02 .10992 .10992 .10992 .14072 .14073 .14073 0.671725

 

10.5 .02 .02 .02 .10996 .10996 .10996 .14091 .14091 .14091 0.6692431

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

When 6 gets smaller, the convergence of this algorithm is not predicted because
of the near singularity of the matrix Q(u), for example, in this example at e = 10‘3,
the algorithm does not converge and the minimum average cost fluctuates between

.6716139, .6716321 and .6716230 all the way up to 100 iterations.

93

9 At 8 = 10"4 the program takes only 6 iterations, but it does not converge to the

right one.

0 At a = 10'5 the cost fluctuates between .6701525 and .6692431 up to 100 itera-

tions.

None of these problems happened when the optimal policy is evaluated by using
the algorithm of section 4.4, see Table 4.1.

CHAPTER 5

CONCLUSIONS AND FURTHER RESEARCH

The real Schur form decomposition is used to develop a reliable and numeri-
cally stable algorithm that transforms any two-time scale system into the singularly
perturbed form. The algorithm comprises two step: first, transform the matrix A into
an ordered real Schur form; second, balance the elements of the ordered real Schur
form such that T1, and Ta are O( 1) and T22 is O(e). If we are interested in the physi-
cal state variables of the system, sufficient conditions are derived to put the two-time
scale system in the singularly perturbed form with all the fast or all the slow state
variables chosen from the original physical variables. These sufficient conditions hold
when the matrix A is given by A =Ao + 6A,, where A0 is a singular matrix with a
semisimple null structure. This case has been extensively studied in the literature.
The significance of these conditions over the previous results in this area is that we
do not require that the matrix A be modeled in the form A = A0 + e A,. Such modeling
form, in general, requires a priori knowledge about the physical nature of the system.
Also, necessary and sufficient conditions are given to retain all the physical states

and the modeling step is achieved by permutation only.

In Chapter 3, stochastic systems which can be modeled as large finite-state
nearly completely decomposable Markov chains are considered. A general transfor-
mation is proposed to decompose the large and the ill-conditioned Markov chain sys-
tem into a reduced order and well-conditioned aggregated matrix (3.37). This
transformation enabled us to compute the exact steady state probability distribution of

the NCDMC (3.43), a problem which is usually encountered in queueing network. It

94

 

95

is shown that the transformations available in the literature to handle this class of
problems are a subclass of our proposed transformation. It is also shown that all the
transformations that satisfy the conditions of Section 3.2 give the same first order
approximation of the aggregated matrix (3.58). A block diagonal transformation,
which is a subclass of the general one, is proposed to simplify and reduce the amount

of computations required to form the aggregated matrix.

Finally, in Chapter 4, the controlled Markov chain of the same class as the one
treated in Chapter 3 is considered. An algorithm composed of two steps, value deter-
mination and policy improvement steps, is proposed. This algorithm computes
optimal policies that attain the minimum average cost per stage [40]. This algorithm

overcomes the ill-conditioning problem associated with the NCDMC.

Suboptimal policies can be computed with less computational effort. In Section
4.5, a modified version of the above algorithm is given to compute the suboptimal
policies and the average cost per stage when 8 = 0. The convergence of this algorithm
is shown. The above algorithm and its modified version are applied to an example to
compute the optimal and suboptimal policies and the average cost per stage

corresponding to these policies.

The possibility of further research and applications of the subject of this disser-
tation is very high. In Chapter 2, for example, the block diagonalization of equation
(2.8) is needed to decompose the singularly perturbed form into two separate subsys-
tems, slow and fast. Since the algorithm proposed in this chapter balances the ele-
ments of the ith off diagonal block to be O( 2 ( S,,-) ), the solution of the Sylvester
equation which block diagonalizes equation (2.8) is expected to be well conditioned
and 0(1). To confirm this expectation it is necessary to show that the separation

between T1, and T22 is O(1).

96

The transformation of Theorem 2.1 is applied to two-time scale systems so that
the state variables of one scale are preserved to be physical ones. A multi-time scale
transformation which preserves the state variables of more than one time scale as

physical variables is a point that needs further investigation.

In Chapter 4, we consider the optimal solution of the average cost per stage
(4.5). An alternative optimal control formulation is to minimize a discounted cost,

defined by
J(u): m%ZB"f(xk,u(xQ) (5.1)
“5 i=0

where B is the discount rate. The algorithm considered in Sections 4.4 and 4.5 could
be extended to this problem to compute optimal and suboptimal policies that minim-

ize the discounted cost of equation (5.1).

In equation (4.43), we have closed-form expressions for i (u) and f, (u). The
minimization of this equation is possible; but disaggregating the optimal reduced pol-
icies into the original ones is an open question and it has to be the subject of further

consideration.

In the algorithm we proposed in Section 4.4, a symbolic inversion of the matrix
V2 (A (u) + e B (u) ) W2 which is required in the value determination step, will reduce

the amount of computations drastically.

I think the theoretical richness of this area combined with the numerous applica-
tions which could be modeled as Markov chains should result in many more contri-

butions to Markov chain modeling and applications in the near future.

APPENDIX

APPENDIX

To show that

-l
V,.: V0-8V03W2 [V2(A+8B)W2] V2,

let us expand [V,( A + e B) 012]" in e and substitute into (3.45) to obtain
°° k
V, = V,— V,(A+ 83) W2 [( V2.1 W2)“: e"[— ( V23 W2)( V2A W '1] ] V2
b=0

-1 —l
=V1—V1(A+EB)W2[V2AW2] {I-EV2BW2 [V2AW2]

2 -12 3 --13
+8 V2BW2 [V2AW2] -8 V2BW2[V2AW2] + ..... }V2.

This can be written as

-1 -l
VS=V1—(V1AW2)[V2AW2] V2—8{ VIBW2 [V2AW2] V2
-1 —l
-V,AW2[V2AW2] V219W2 [VZAWZ] V2}
-1 -1
+22{ V,BW2[V2AW2] (VZBW2)[V2AW2] V2

-1 -l 2
—(V,A W7) [Vzsz] [VZBWZ [V2.4 W2] ] V2} + e3....
Substitute for V0 to obtain

-1 -l
VS=V0’8[V1—(V1AW2)[V2AW2] V2]BW2[V2AW2] V2

-1
+e2 [V,—(V,AW,) [V271 W2] V2]BW2

97

98

—l -l
X[V2AW2] (V2BW2) [V2AW2] V2+ ......

The terms in the square brackets equal to V0; therefore,

-1 -l
VS=Vo-EV03W2{ [V2AW2] -8[V2AW2]
-1 2 —l
[V2BW2] [V2AW2] +8 [V2AW2]

[V2 B W?) [V2 A W2]—l]2 — } V2

-1
=V0—8V03W2[V2(A+SB)W2] V2.

LIST OF REFERENCES

10.

ll.

12.

LIST OF REFERENCES

P. V. Kokotovic, H. K. Khalil and J. O’Reilly, Singular Perturbation Methods in
Control: Analysis and Design. New York: Academic Press, 1986.

P. Kokotovic and H. Khalil, Ed., Singular Perturbations in Systems and Control,
IEEE press, 1986.

G. H. Golub and J. H. Wilkinson, "Ill-conditioned Eigensystems and the Com-
putation of the Jordan Canonical Form," SIAM Rev., vol. 18, pp. 578-619,
1976.

J. H. Chow and P. Kokotovic, "Eigenvalue Placement in Two-Time-Scale Sys-
tems," Proc. of IFAC Symposium on large scale systems, 1976.

J. Anderson, "Decoupling of Two-Time-Scale Linear Systems," Proc. of Joint
Automatic Control, 1978.

G. P. Syrcos and P. Sannuti, "Singular Perturbation Modeling of Continuous and
Discrete Physical System," Int. J. Control, vol. 37, pp. 1007-1022, 1983.

B. Avramovic, P. V. Kokotovic, J. R. Winkelrnan and J. H. Chow, "Area
Decomposition for Electromechanical Models of Power Systems," Automatica,
vol. 16, pp. 637-648, 1980.

A. J. Fossard, M. Berthelot and J. F. Magni, "On Coherency-based Decomposi-
tion Algoritrns," Automatica, vol. 19, pp. 247-253, 1983.

I. J. Perez-Arriaga, G. C. Verghese and F. C. Schweppe, "Selective Modal
Analysis with Application to Electric Power Systems, Part 1, IEEE Trans. power
Appar. and Sys., vol. PAS-101, 1982

J. H. Chow, Ed., Time-Scale Modeling of Dynamic Networks with Application
to Power System," Lecture Notes In Control Information Science, vol. 46, New
York: Springer-Verlag, 1982.

R. G. Phillips and P. V. Kokotovic, "A Singular Perturbation Approach to
Modeling and Control of Markov Chains," IEEE Trans. Automat. Contr., vol.,
AC-26. pp. 1087-1094, 1981.

M. Coderch, A. S. Willsky, S. S. Sastry and D. A. Castanon, "Hierarchical
Aggregation of Linear Systems with Multiple Time Scales," IEEE Trans.

99

HH-

13.

14.

15.

16.

17.

18.

19.

20.

21.

22.

23.

24.

25.

100

Automat. Contr., AC-28, pp. 1017-1030, 1983.

H. K. Khalil, "Time Scale Decomposition of Linear Implicit Singularly Per-
turbed Systems," IEEE Trans. Automat. Contr., vol. AC-29, pp. 1054-1056,
1984.

X.-C Lou, G. C. Verghese, A. S. Willsky and M. Vidyasagar, "An Algebraic
Approach to Analysis and Control of Time-Scales," Proc. 1984 ACC, June
1984, pp. 1368-1372.

X.-C Lou, "An Algebraic Approach to the Analysis and Control of Time-
Scales," PH.D. Dissertation, Dept. Elec. Eng. Comput. Sci. M.I.T., 1985.

G. H. Golub and C. F. Van Loan, Matrix Computations, Baltimore, The Johns
Hopkins University Press, 1983.

G. W. Stewart, Introduction to Matrix Computations, New York, Academic
Press, 1973.

B. T. Smith et al., Matrix Eigensystems Routines-EISPACK Guide, 2nd ed.
(Lecture Notes in Computer Science), vol. 6, New York: Springer-Verlag, 197 6.

A. J. Laub, "A Schur Method for Solving Algebraic Riccati Equations," IEEE
Trans. Automat. Contr., vol. AC-24, pp. 913-921, 1979.

G. W. Stewart, "Algorithm 506: HQR3 and EXCHG: Fortran Subroutines for
Calculating and Ordering the Eigenvalues of a Real Upper Hessenberg Matrix,"
ACM Trans. Math. Software, vol. 2, pp. 275-280, 1976.

D. S. Flamm and R. W. Walker, "Remark On Algorithm 506," ACM Trans
Math. Software, vol. 8, pp. 219-220, 1982.

B. N. Parlett and C. Reinch, "Balancing a Matrix for Calculation of Eigenvalues
and Eigenvectors," Numer. Math., vol. 13, pp. 293-304, 1969.

E. E. Osborne, "On Pre-Conditioning of Matrices," Jour. ACM 7, pp. 338-345,
1960.

H. A. Simon and A. Ando, "Aggregation of Variables in Dynamic Systems,"
Econometrica 29, pp. 111-138, 1961.

P. J. Courtois, Decomposability Queueing and Computer System Applications.
New York, Acadirrric Press, 1977.

 

26.

27.

28.

29.

30.

31.

32.

33.

34.

35.

36.

37.

38.

39.

 

101

P. J. Courtois, "Error Analysis in Nearly-Completely Decomposable Stochastic
Systems," Econometrica, vol. 43, pp. 691-709, 1975.

P. J. Courtois, "Decomposability, Instabilities, and Saturation in Multiprogram-
ming Systems," Commun. ACM 18, pp. 371-377, 1975.

A. Brandwajn, "A Model of a Time Sharing Virtual Memory System Solved
Using Equivalence and Decomposition Methods," Acta Informatica 4, 11-47,
1974.

H. Vantilborgh, "Aggregation with an Error of C(82)," J. ACM vol. 32, pp.
162-190, 1985.

R. G. Phillips, "Near Optimal Policies for Large Scale Markovian Decision
Processes," Proc. Int. Large Scale Sys. Symp., pp. 507-512, 1982.

W. J. Stewart, "A Comparison of Numerical Techniques in Markov Modeling,"
Commun. ACM 21, pp. 144-152, 1978.

J. R. Koury, D. F. McAllister and W. J. Stewart, "Iterative Methods For Com-
putin g Stationary Distributions of Nearly Completely Decomposable Markov
Chains," SIAM J. Alg. Disc. Math. 5. PP. 164-186, 1984.

D. F. McAllister, G. W. Stewart and W. J. Stewart, "On a Rayleigh-Ritz
Refinement Technique for Nearly Uncoupled Stochastic Matrices," Linear A1 g.
Applications 60, pp. 1-25, 1984.

W-L Cao and W. J. Stewart, "Iterative Aggregation / Disaggregation Techniques
for Nearly Uncoupled Markov Chains," J. ACM 32, pp. 702-719, 1985.

J. L. Doob, Stochastic Processes. New York: Wiley Press, 1953.

P. Varaiya, "Optimal and Suboptimal Stationary Controls for Markov Chains,"
IEEE Trans. on Automat. Contr., vol. AC-23, pp. 388-394, 1978.

J. L. Popyack, R. L. Brown, and C. C. White, 111, "Discrete Versions of an
Algorithm Due Varaiya," IEEE Trans. on Automat. Contr., vol. AC-24, June
1979. ‘

S. H. Javid, "Nested Optimization of Weakly Coupled Markov Chains," Proc. of
the 18th Annual Allerton Conf. on Communication, Control and Computing,
Univ. of Illinois, pp. 881 - 890, 1980.

D. Teneketzis, S. H. Javid, and B. Sridhar, "Control of Weakly-Coupled Markov
Chains," Proc. 1980 CDC, pp. 137-142, 1980.

40.

41.

42.

43.

45.

46.

102

D. P. Bertsekas, Dynamic Programming: Deterministic and Stochastic Models.
Prentice-Hall, Inc., 1987.

R. A. Howard, Dynamic Programming and Markov Processes. Cambridge, MA:
MIT Press, 1960.

D. J. White, "Dynamic Programming, Markov Chains, and the Method of Suc-
cessive Approximations," J. Math. Anal. Appl., vol. 6, pp. 373-376, 1963.

H. Kushner, Introduction to Stochastic Control. New York: Holt, Rinehart and
Winston, 1971.

A. Cohen et.al., "Research and Development of a Unified Approach to Opera-
tions Scheduling for Electric Power Under Uncertainty," Final Report, USDOE
Contract No. DEACOl-79ET‘29243, Oct., 1982. of Energ. Sys., Final Report,
1982.

J. Stoer and R. Bulirsch, Introduction to Numerical Analysis. New York:
Springer-Verlag, 1980.

L. Kleirrrock, Queueing Systems, Vol. I and II, New York: Wiley, 197 6.