A performance bound for nonlinear control systems

A PERFORMANCE BOUND FOR NONLINEAR CONTROL SYSTEMS

By

RAFAEL J. FANJUL JR.

A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA

1996

Copyright 1996

by

Rafael J. Fanjul Jr.

To my parents, Rafael James and Estela Linares Fanjul,

my grandmother, Margaret Stewart Fanjul,

my aunt, Sheila Stewart,

and

my nephews, Frankie and Rafi Montalvo

ACKNOWLEDGEMENTS

I would like to express my gratitude to my advisor, Professor Hammer, for his

encouragement, guidance and wisdom throughout the course of my studies. Professor

Hammer has always found time to discuss my work and for that I am especially

thankful.

I am grateful to Professor Crisalle for his many hours of assistance in formu-

lating my research. I wish to thank Professor Schwartz for her help during the early

stages of my research. In addition, I would like to thank Professor Anderson and

Professor Principe for serving on my committee.

A number of my fellow students have provided both inspiration and advice.

In particular, I wish to thank Victor Brennan, Patrick Walker, Aiguo Yan, and Kuo-

Huei Yen for their help over the past four-and-a-half years. Finally, I would like to

thank Margaret Fanjul Montalvo and Christopher Riffer for editing the manuscript.

TABLE OF CONTENTS

ACKNOWLEDGEMENTS ..............

ABSTRACT ......................

CHAPTERS

1 INTRODUCTION ...............

1.1 Background ................
1.2 Notation ...... ...........

2 TERMINOLOGY AND BASICS .......

2.1 Right Fraction Representations
2.2 Generalized Right Inverse .

3 THE PERFORMANCE BOUND..

and Coprimeness

3.1 Causality of System s .........................
3.2 Systems with the Lipschitz Norm . . . . .
3.3 Approximate Right Inverse . . . . . .
3.3.1 The Existence Theorem of Best Approximate Right Inverse
3.4 The Measure of Right Singularity . . . . .

4 CALCULATION OF THE PERFORMANCE BOUND ........

4.1 The Estimate of the Performance Bound . . . .
4.2 Practical Implementation Techniques . . . .
4.3 Permanent Magnet Stepper Motor Model Example . .
4.3.1 Simulation ............................
4.3.2 The Estimate of the Performance Bound P() . .
4.4 Aerodynamic Model Example . . . . .
4.4.1 Sim ulation . . . ... . .
4.4.2 The Estimate of the Performance Bound P() . .
4.5 Multivariable Process Control Model Example . . .
4.5.1 Simulation ............................
4.5.2 The Estimate of the Performance Bound P() . .

5 CONCLUSION ..............................

5.1 Sum m ary . . . . . . . .
5.2 Future Directions ..........................

REFERENCES ......................... .......... 77

BIOGRAPHICAL SKETCH .............................. 80

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

A PERFORMANCE BOUND FOR NONLINEAR CONTROL SYSTEMS

By

Rafael J. Fanjul Jr.

December 1996

Chairman: Professor Jacob Hammer
Major Department: Electrical and Computer Engineering

This research focuses on the control of a nonlinear system whose output sub-

ject to an additive disturbance. The main interest is in the investigation of controllers

that reduce the effect of the disturbance on the system output. Usually, it is not possi-

ble to construct a controller that completely eliminates the effects of the disturbance.

It is then of interest to find how well the "best" controller can attenuate the effect of

the disturbance. The main result of this dissertation is a performance bound, that

provides an estimate of the best disturbance attenuation that can be achieved for a

given system, using a causal controller that renders the system internally stable.

An approximate right inverse of a nonlinear system E is introduced to fa-

cilitate the derivation of the performance bound and the development of nonlinear

controllers. An approximate right inverse is constructed to be stable and causal for

implementation purposes. The difference between an approximate right inverse and

a right inverse is that the right inverse may not be both stable and causal. The role

of an approximate right inverse is to approximate a disturbed signal with a signal in

the image of the system E. An approximate right inverse can be constructed for any

system E.

The calculation of the performance bound involves an optimization process of

finding a global maximum of a non-convex function. For several cases, a nonlinear

programming algorithm is developed to handle the optimization. To demonstrate the

application of this performance bound, it is calculated for three practical systems:

1. the voltage control of a permanent magnet stepper motor;

2. the longitudinal control of an aircraft; and

3. the multivariable process control of a regulator that regulates the liquid level

in a pressurized tank.

CHAPTER 1
INTRODUCTION

Over the last 2 decades, there has been considerable interests in the literature

in the derivation of optimal controllers that minimizes the effect of the disturbance on

the output of a control system. This dissertation addresses this question for the case

of nonlinear systems. The main objective is to derive an estimate of the performance

of an optimal controller, by deriving a bound on the effect of the disturbance has on

the system output when the optimal controller is used. Using this bound, we can

then gauge the performance of suboptimal controller, to see how well they compare

to the optimal ones. Suboptimal controllers may be much easier to implement than

their optimal counterparts. Thus, our performance bound can be use to find simple

controllers whose disturbance attenuation properties are close to those of optimal

disturbance attenuating controllers.

The performance bound derived in this dissertation came from the require-

ment to characterize the performance of a controller that reduces the effect of the

disturbance on the output to a minimum. The original work provided is the derivation

of the performance bound and nonlinear programming algorithm in how to calculate

the performance bound. The application of the performance bound for disturbance

attenuation is proposed as future directions.

The basic design problem to which this performance bound relates is the

problem of disturbance attenuation for nonlinear control systems. Specifically, this

dissertation discusses the following configuration. In the configuration of Figure 1.1,

E is the nonlinear system to be controlled. C represents an equivalent controller that

incorporates all the control elements of the loop. The external (or reference) signal

is denoted by v; the disturbance signal is denoted by d; and the output signal is

denoted by z. The closed loop system is required to be internally stable. Internal

stability signifies that a configuration can tolerate small disturbances on its external

and internal ports (including ports within the equivalent controller C) without losing

stability. The equations that describe Figure 1.1 are

z = d +y, y = Eu, u= C(v,z). (1.1)

In Figure 1.1, EC represents the appropriate equivalent system. This can be

expressed in notation

z = c(v, d) (1.2)

where the output signal z is determined by the signals v and d and depends on the

system E as well as the equivalent controller C. We would like to reduce as much

as possible the effects of d on z. Our bound, which is the attainable performance

for control of a nonlinear system whose output is subject to an additive disturbance,

provides an estimate of the minimal effect of d on z. Using the estimates of the min-

imal effect, we can evaluate controllers. The analysis of the dissertation is restricted

to the case of discrete-time systems.

The desired response, with or without the disturbance signal for the configu-

ration in Figure 1.1, is z = Ev. To null out the disturbance signal, the design of the

equivalent controller C would be such that y = Eu = Ev d. Substituting for y in

d

v u ,

Figure 1.1: The block diagram of Ec.

the first equation of (1.1), this yields z = Ev which is our stated goal. For the sake of

argument, assume that the system E has a stable and causal inverse system E-1. The

input signal u to the system E would then be required to be u = E-I(Ev d). The

inverse system E-1 in some cases might not be implementable because it is not sta-

ble or it is not causal. For those cases, approximations of the inverse systems would

then be used for practical implementation. The approximations will be selected in a

manner suited to approximate the signal Ev d with a signal in the image of E.

The composition of the system and its approximate inverse system forms a

nearly identity system. For points inside the image of the system E, the composition

will appear as the identity. For external points, which are points outside the image of

the system E, the composition will produce a point in the image of E that is closest

to the external point. Our performance bound is a norm that gauges the closeness of

the identity system and a nearly identity system formed from the system E and its

approximate inverse.

The performance where Ec can attain by using our bound is determined by the

inherent properties of the system E. The bound is directly related to an approximate

inverse. The role of an approximate inverse is to approximate the signal u with a

signal v in the image of E. Therefore, the approximate invertibility of the system

E provides a measure of the ability to match a set S by a subset of the image of

the system E. A measure of singularity is introduced as an indicator of approximate

invertibility which is then used to derive the performance bound. To demonstrate the

application of this performance bound, it is calculated for three practical systems:

1. the voltage control of a permanent magnet stepper motor;

2. the longitudinal control of an aircraft; and

3. the multivariable process control of a regulator that regulates the liquid level

in a pressurized tank.

The dissertation is organized as follows. Chapter 2 contains the basic notions

of nonlinear systems including the theory of fraction representation. Chapter 3 dis-

cusses an approximate right inverse and the performance bound. Chapter 4 discusses

the main results of the dissertation, including calculations of the performance bound

for applications. Chapter 5 contains a summary and future directions.

1.1 Background

This section contains a qualitative survey of some topics in nonlinear control

that are important to the dissertation. A more technical discussion of these topics is

provided in Chapter 2.

The description of Ec in Figure 1.1 has its origins in the theory of fraction rep-

resentation of nonlinear systems. It can be stated that a right fraction representation

of a nonlinear system E is a factorization of E into a composition form E = PQ-1,

where P and Q are stable systems with Q being invertible. In Hammer [16, 18], tools

for a compact and perceptive statement of results were developed for the theory of

fraction representation of nonlinear systems. The fraction representation E = PQ-1

is said to be coprime when the systems P and Q are right coprime. An attribute of

right coprime fraction representation is that every instability of the inverse system

Q-1 is also an instability of the system E; i.e., there is no cancellation of instabilities

within the composition PQ-1.

A causal (respectively, strictly causal) system E is one where the values of the

output sequence Eu up to and including index i (respectively, i + 1) depend only on

the values of the input sequence u up to index i. A system is bicausal if it is causal

and if it possesses a causal inverse.

In this dissertation, it is assumed that the system E being controlled can be

stabilized, that it is strictly causal, and that it possesses a right coprime fraction

representation of the form E = PQ-1, with Q being bicausal. The stabilization

5

assumption on the system E is necessary because the closed loop system is required to

be internally stable. Strict causality is placed on the system E as a handy assumption

to certify that the closed loop system is well posed. Strict causality is not an essential

condition and it can be replaced with plain causality combined with a well-posedeness

requirement. This dissertation relies on stabilization theory which applies only to

systems possessing right coprime fraction representations.

A result from [22] gives a simple parameterization of the set of all system

responses that can be obtained through internally stable control of a given system.

The control scheme to control E is shown Figure 1.1. The parameterization provides

a clear indication of the effects of the disturbance on the response of the stabilized

closed loop system. This result can be summarized as follows. Let E = PQ-1 be

a right coprime fraction representation (with a bicausal "denominator" Q) of the

system being controlled. Then,

1. For every causal equivalent controller C for which the closed loop system of

Figure 1.1 is internally stable, there exists a stable and causal system q(v, d)

such that

Ec(v, d) = d + Pq(v, d) = [I + PO(v, .)]d (1.3)

where P is the "numerator" of the right coprime fraction representation of E

and I denotes the identity system.

2. Conversely, for every stable and causal system q(v, d), there is an internally sta-

ble control configuration around the system E for which the equivalent system

Ec(v, d) satisfies Ec(v, d) = [I + PO(v, .)]d.

Thus, equation (1.3) provides a complete parameterization of the class of all

responses {Ec} that can be obtained by internally stable control of the system E,

with the stable and causal system 0 serving as the sole parameter. In other words,

for every q4, there is an equivalent controller C that internally stabilizes E and yields

the response (1.3). Conversely, every equivalent controller C that internally stabilizes

E generates an explicit 4.

On a superficial inspection of (1.3), the attainable performance is prescribed

by the "numerator" system P of E because it is the only fixed quantity apart from

the identity. For a linear time-invariant system, it was shown in [50] that the at-

tainable performance was determined by the location of its right half-plane zeros.

For a nonlinear system E, the instabilities of P-1 (the inverse of the "numerator"

system P) would be analogous to the right half-plane zeros of a linear time-invariant

system. Using this analogy it could be inferred that the instabilities of P-1 would

limit performance. The performance bound developed in the dissertation will further

confirm this hypothesis.

1.2 Notation

The following notational convention will apply unless otherwise stated.

?m

S(Rm)
Dx

f(E) > A

lul

p(u)

IIL |
B

Lip(), S(RP)

II1 IILip
M, (U, S2)

oe(S1, S2)

P(E)

The set of m-dimensional real vectors.

The set of m-dimensional real vector sequences.

The A-step shift operator.
The system E has latency of at least A.

The -norm of the sequence u.

The weighted -norm of the sequence u.

The Lipschitz semi-norm of the system E.

The set of all stable, causal and recursive systems E: D S(R)
where D is a bounded subset of S(3?m).

The subset of B with ||EI| < 0c.

The Lipschitz norm of the system E.
The distance from any u E S1 and the set S2.

The maximal distance from any u E Si and the set S2-

The right singularity of measure of the system E.

CHAPTER 2
TERMINOLOGY AND BASICS

This chapter contains a summary of the principal mathematical results which

are needed in the dissertation. The presentation is for discrete-time time-invariant

nonlinear systems. Almost all the necessary mathematics are contained in [9, 14, 17,

18, 22, 27].

The set of real numbers is denoted by R. The set of m-dimensional real

vectors is denoted by Rm. The set of all sequences u is denoted by S(Rm) where

u = {uo, 1, U2, ...} of m-dimensional real vectors ui E RM", i = 0, 1,2,.... Given

the sequence u E S(Rm), the ith element is denoted by ui. The set of elements

{ui, ui+i,..., Uj where j > i > 0 is denoted by ui. A system E, from an input/output
perspective, is a map E: S(Rm) --+ S(RP) which transforms input sequences of m-

dimensional real vectors into output sequences of p-dimensional real vectors. The

image of a subset S C S("m) through the system E is denoted by E[S]. The entire

image of system the E is denoted by Im where ImE d- E[S(Rm)]. Given a system

E: S(Rm) S(RP) and an input sequence u E S(Rm), we denote by Eu]i = y, the

ith element of the output sequence y = Eu, and by y = Eu]j the set of elements

{yi, yi+1,... j} where j > i > 0 are integers.
There are two kinds of binary operations, addition and composition. For a pair

of systems E, E2: S(Rm) -- S(RP), the sum is defined, as usual, by (iE + E2)u =

Elu + E2u for all sequences u E S(Rm); the right side of the last formula is the
usual elementwise addition of sequences of real vectors. Composition is the usual

composition of maps.

Definition 2.0.1 A system E: S(Rm) -- S(RP) is causal (respectively, strictly caus-

al) if it satisfies the following condition. For every integer i > 0 and for every pair

of input sequences u,v E S(Rm) satisfying uo = V, the output sequences satisfy

Eu]o = Ev]b (respectively, Eu]o' = Ev]'+).

A system M: S(Rm) -+ S(,m) is bicausal if it is causal and if it possesses a causal

inverse.

A system E: S(Rm) -- S(RP) is called a recursive system if there is a pair of

integers r, yi > 0 and a function f: (,),+71 X (Rm)P+1 --+ RP such that, for every input

sequence u E S(Rm), the corresponding output sequence y = Eu can be computed

recursively in the form

Yk+q+l = f(Yk,..., Yk+,, Uk,..., Uk+p) (2.1)

for all integers k > 0. The initial conditions yo,..., y, must of course, be specified

and fixed. The function f is called a recursion function of E. For causality (strict

causality), the condition q +1 > p (q > y_) is placed on system E. The class of strictly

causal systems include every system E: S(Rm) -+ S(RP) that can be represented in

the form

Xk+1 f(xk, Uk)

yk = h(xk), k = 1,2,.... (2.2)

Here, u E S(Rm) is the input sequence; y E S(RP) is the output sequence; and

x E S(Rn) is an intermediate sequence of "states." In the case that the maps f: Rn x
Rm -,* and h: R" -+ R are continuous, then the system (2.2) constitutes a

continuous realization of the system E.

For a real number 0 > 0, the set of all vectors in m with components in the

closed interval [-0, ] is denoted by [-0, ]m. The set of all sequences u E S(~m)

with elements ui belonging to [-0,0]m for all integers i > 0 is denoted by S(Om).

Thus, S(Om) consists of all sequences bounded by 0. It follows then that a system

E: S(Rm) -- S(RP) is BIBO (Bounded-Input Bounded-Output)-stable if for every real

number 0 > 0, there exists a real number M > 0 such that E[S(0m)] C S(MP). A

sequence u E S(R"m) is said to be bounded if there is a real number 0 > 0 such that

u e S(Om).

The basic notion of stability that is used in this dissertation is related to con-

tinuity with respect to a metric. Two norms are particularly useful in this context to

derive a metric: the e-norm and the weighted ~-norm. The t"-norm is denoted by

- ; for a vector a = (al, a2,... ,am) Rm, it is simply lal d= max{ ai|, la21,. ., *amm}.
For a sequence u E S(Rm), the e"-norm is given by

Iul df sup IuIl. (2.3)
i>O
The weighted t"-norm is denoted by p, and is given by

p,(u) f sup (1 + e)-i u, > 0 (2.4)
i>O
for a sequence u e S(Rm). For purposes of this dissertation different values of e will

not affect our results; hence, the subscript will be dropped and the weighted e-norm

is simply denoted by p. The use of the weighted oo-norm simplifies mathematical

arguments over the -norm because the bounded set of sequences S(0m) is compact

with respect to p. The norm p induces a metric p on a given S(Rm), for every pair

of elements u, v S(Rm), by p(u, v) =f p(u v). Formally, the notion of stability

employed in this dissertation is as follows.

Definition 2.0.2 A system E: S(Rm) --- S(RP) is stable with respect to the metric

p if it is BIBO-stable, and if the restriction E: S(am) -- S(RP) is continuous with

respect to p for every real number a > 0.

Definition 2.0.2 is usually referred to as input/output stability. The following

concept, which describes a weak form of uniform continuity with respect to the "-

norm, plays a fundamental role in stabilization theory (see [17]).

Definition 2.0.3 A stable system E: S(3m) -- S(RP) is differentially bounded if
there is a pair of real numbers e, 0 > 0 such that, for every pair of sequences u E

S(Rm) and v E S(em), one has IE(u + v) E(u)Il 8.

So far, only stability properties of individual systems have been mentioned.
When several individual systems are combined into a composite system, a stronger

notion of stability is required, and it is usually referred to as internal stability. Internal

stability guarantees desirable stability properties of the composition, and takes into

account the effects of various disturbances and noises that may affect the component

systems. Consider a composite system E(s) that consists of s individual systems,

labeled El,.. ., E where SE: S(Rmm()) -- S(RP()), i = 1,..., s. Individual entries in

the list E1,..., ES may represent summers, multipliers, etc. Let u E S(Rm) be the

external input sequence of the composite system, and let y E S(Rp) be its output

sequence. Let uj E S(Rm(j)) be the input sequence of the system Ej within the

configuration, and let y3 E S(Rp(j)) be its output sequence. The interconnections

among the subsystems are then characterized by a set of qualities ui = yj(), which

determine to which output each input is connected. The external signal u is now

augmented by s new input signals qi E S(Rm(i)), i = 1,..., s, and set u'i yj(i) + .

For each i, the 7i acts as an additive disturbance on the input port of the system E.

The disturbances are all assumed to be bounded by a real number 6 > 0, so that in

fact i'i E S(bm(i)), i= 1,...,s.

Let E*s*:S(Rm) x S(Rm(i)) x ... x S(GRm()) -_ S(RP) x S(Rp(x)) x ... x

S(R(S)): (u, 71, ..., .) -> E*S*(u, 1,..., 9) denote the system induced by the in-
terconnected system E(S) and the disturbances, having the input signals u, 1i/,... ,7/

and the output signals y, y1,..., ys, respectively.

Definition 2.0.4 The composite system E(s) is internally stable if the system E*`* is

stable in the sense of Definition 2.0.2. The composite system E() is strictly internally

stable if, besides being stable, the system E*`* is also differentially bounded.

Definition 2.0.5 A system E: S(Rm) --+ S(RP) is entirely stabilizable if there is a
strictly internally stable control configuration that stabilizes E over the entire input

space S(Rm).

2.1 Right Fraction Representations and Coprimeness

A right fraction representation of a system E: S(Rm) -+ S(RP) is determined
by three quantities: a subset S C S(Rq9), q > 0, called the factorization space; and
two stable systems P: S -+ S(RP) and Q: S -- S(Rm), where Q is invertible, such

that E = PQ-1. A right fraction representation E = PQ-1 is coprime whenever

the stable systems P and Q are right coprime according to the following definition

([16, 18]). (Let G: S1 -- S2 be a map, where S1 C S(Rm) and S2 C S(Rn) are subsets.
For a subset S C S(R"), we denote G*[S] the inverse image of S through G, i.e., the
set of all sequences u E S1 satisfying Gu E S.)

Definition 2.1.1 Let S C S(Rq() be a subset. Two stable systems P: S -* S(RP) and

Q: S --+ S(Rm) are right coprime whenever the following conditions hold.

1. For every real number r > 0 there is a real number 0 > 0 such that

P*[S(TP)] n Q*[S(Tm)] c S(0q).

2. For every real number 7 > 0 the set S n S(rq) is a closed subset of S(rq) (with

respect to the topology induced by p).

The concept of a homogeneous system is of key importance to the theory of
right coprime fraction representations of nonlinear systems. A homogeneous system

has the property of being a continuous map whenever its outputs are bounded. The

precise definition is as follows.

Definition 2.1.2 A system E: S(Rm) -- S(RP) is a homogeneous system if the fol-

lowing holds for every real number a > 0: for every subset S C S(am) for which

there exists a real number r > 0 satisfying E[S] C S(TP), the restriction of E to the
closure 3S of S in S(am) is a continuous map E: S -- S(7r).

The importance of homogeneous systems to the dissertation is stated in the
next two theorems.

Theorem 2.1.3 [17] An injective system E: S(Rm) -* S(RP) has a right coprime

fraction representation if and only if it is a homogeneous system.

Theorem 2.1.4 [17] Let E: S(Rm) S(RP) be a recursive system. If E has a recur-

sive representation yk+7+1, = f(Yk, *, Yk+n,, Uki,..., k+,) with a continuous recursion

function f, then E is a homogeneous system.

The following two theorems are important for stabilization theory. Theo-

rem 2.1.5 defines the technical details in the stabilizing system E in Figure 2.1 with
r = B-1 and

ization, which is the denominator system contains the exact information about the

instabilities of the system, of right coprime fraction representation for stabilization.

Theorem 2.1.5 [19] Let E: S(am) -+ S(RP) be a system with a bounded input space

S(am), a > 0, and assume it has a right coprime fraction representation E = PQ-',
where P: S -+ S(RP) and Q: S -+ S(am), and where S C S(Rq) for some integer q >

s +I I I I Y

Figure 2.1: The block diagram of E(,).

0. Then, for every stable system M: S -- S with a stable inverse system M-: S S,
there exists a pair of stable systems A: S(RP) -+ S(Rq) and B: S(am) -+ S(Rq) such
that AP + BQ = M.

Theorem 2.1.6 [22] Let E: S(?m) -+ S(RP) be a system having a right coprime
fraction representation E = PQ-1, where P: S S(P), Q: S -- S( ?m), S C

S("R), and where Q is bicausal. Let D: S(R') -- S(Rm) be any stable and causal
system for which the composition ED: S(R~ ) -+ S(RP) is stable. Then, there is a
stable and causal system b: (R') -- S such that D = Q0.

In Figure 1.1, E: S(Rm) S(RP) is a strictly causal system that needs to be
controlled. It will be convenient to regard the external input sequence v as a fixed
"parameter" while regarding the disturbance d as an external input, i.e., to consider

appropriate partial functions. Since no restriction will be placed on v and d, this will
have no effect on the validity of the final result. In the same spirit, we shall use the
notation

I(v)z C(v,z)

in which v can be intuitively viewed as a parameter of the system 1(v), while z is
its input. Control is achieved by a causal nonlinear dynamic controller C: S(Rm) x
S(RP) -, S(Rm): (v, z) C(v, z). For Figure 2.2, T(v)z f C(v, z) was defined with
the system T(v): S(RP) -- S(Rm): z J(v)z = u. A result of [22] is now stated

d z

Figure 2.2: The block diagram of Eq().

formally in Theorem 2.1.7. which furnishes a parameterization of the class of all

systems that can be obtained from a given system E by internally stable control.

Theorem 2.1.7 [22] Let E: S(Rm) -- S(3RP) be a strictly causal system having

a right coprime fraction representation E = PQ-1 with a bicausal denominator

Q: S(Rm) -- S(Rm). Assume that E can be strictly internally stabilized by a con-

troller that admits representation Ci(s,z) = G-(z)[s + Tz] = u. Then, referring to

(1.2), the following is true. The class of all input&disturbance/output responses Ec
that can be achieved through internally stable control of E is given by

{Ec(v, d) = [I + Pq(v)]d, 0(.): S(Rm) x S(WP) -+ S(Rm): (v, d) '-4 (v)d is a stable
and causal system}.

A special case of Theorem 2.1.7 is when the system E is stable. In that case,

the right coprime factorization of E = PQ-1 can be taken as P = E and Q = I, and

the following is obtained.

Corollary 2.1.8 [22] Let : S(Rm) --+ S(P) be a strictly causal, stable and differ-

entially bounded system. Then, referring to (1.2), the following is true. The class

of all input&disturbance/output responses Ec that can be achieved through internally

stable control of E is given by

{Ec(v, d) = [I + EO(v)]d, 0(.): S(Rm) x S(RP) -- S(Rm): (v, d) -* 0(v)d is a stable
and causal system}.

2.2 Generalized Right Inverse

For the nonlinear recursive system E: S(Rm) -+ S(,RP) to possess a right in-

verse, the system E is required to be surjective. By restricting the range of E to the

image of E, a surjective system Er: S(Rm) Im E is obtained with its right inverse

system E*: Im E --+ S(?m). Let Eg: S(Re) -- S(m) be any extension of E* from

Im E to the whole space S(RP). Then, for every element y Im E, it is evident that

EE-y = EE*y = y. The system E3 is a generalized right inverse of the system E. The

generalized right inverse of system E is non-unique when system E is not bijective.

The next theorem and its corollaries from [14] state that a recursive system has a

recursive generalized right inverse.

Theorem 2.2.1 [14] A recursive system E: S(Rm) S(RP) has a recursive gener-

alized right inverse Eg: S(RP) -- S(m).

A couple of notes from the proof of Theorem 2.2.1. Let the system E*: Im E -+ S(R~m)

be the recursive system represented by

Uk+- = g(Uk,.. Uk+--l1, Yk ... I Yk+p+n+1). (2.5)

In order to extend the domain of E* from Im E to all S(RP), let ge: ('m), x (' ),++2
-- m be any extension of the function g. Then, the recursive system Eg: S(Rep) -*

S(Rm) represented by

Uk+p = ge(uk,..., Uk+- 1, Yk, Yk+,p+t+1) (2.6)

is a generalized right inverse of E.

Corollary 2.2.2 [14] Let E: S(Rm) -- S(PF) be a recursive isomorphism. Then, its

inverse E-1: S(RP) -- S(Rm) is also a recursive isomorphism.

Theorem 2.2.1 cannot be generalized to the case of an arbitrary domain D.

Nevertheless, Theorem 2.2.1 can be generalized to the case when the domain D is

in the following particular form. A subset > C S(Rm) is recursive if there exists

an integer ( and a function a assigning each point (ao,...,ac ) E (sm)4+1 a subset

a(ac) C m such that D consists of exactly all sequences u E S(Rm) satisfying

Uk+4+1 E a(uk+ ) for all integers k. The function a is called the generating function
of D. For instance, if II: S(Rq) -- S(Rm) is a recursive system, then ImII is a

recursive subset of S(Rm).

17

Corollary 2.2.3 [14] A recursive system E: D -- S(WP), where D is a recursive

subset of S("m), has a recursive right inverse E*: Im E -- D.

CHAPTER 3
THE PERFORMANCE BOUND

This chapter starts by discussing causality of systems and Lipschitz norms

which will be useful in establishing the existence theorem of a best approximate

right inverse. The performance bound will later be defined as a lower bound for the

performance index of the approximate right inverse optimization problem.

3.1 Causality of Systems

Recall a causal system E is one for which the values of the output sequence

Eu up to and including index i depend only on the values of the input sequence u up

to index i. This was precisely defined in Definition 2.0.1.

For an integer A, the A-step shift operator is denoted by DX, defined, on any

sequence u by

D u]k ef Uk-A (3.1)

for all integer k for which Uk-A exists.

Definition 3.1.1 Let : S1 -+ S2 be a system, where Si C S(~m) and S2 C S(RP).

The system E has latency of at least A if there is an integer A such that, for every

pair of input sequences u,v E S1 and for every integer k > 0, the equality Uo = Vk

implies Eu]k+A = -v]k+

We write C(E) > A if the system E has latency of at least A. Intuitively,

the latency represents a 'time delay' incurred in the propagation of changes from the

input of E to the output of E. It is a simple consequence of the definition that a

system E is causal if and only if L(E) > 0, and it is strictly causal if and only if

(E) > 1. From [17], a few simple properties of latency are listed.

Theorem 3.1.2 [17] Let E: S1 -- S2 and E2:S2 -* S3, where S1 C S(Rm), S2 C
D I[S(RP)] and S3 C DA"+A2 [S(Rq)], be systems, each with a well-defined latency,
respectively (E1) > A1 and (E2) > A2. Under this hypothesis the composition

E d= E2E has a well-defined latency, and (E) > A1 +A2. In particular if A1 +A2 > 0
then E is a causal system.

Theorem 3.1.3 [17] Let E: S1 -- S2, where S1 C S(Rm) and S2 C S(RP), be a re-

cursive system with a recursive representation yk+n+1 = f(Yk, Yk+,, uk,..., Uk+,)

(and some fixed initial conditions). Under this hypothesis E has well-defined latency,

and (C) > 7 + 1- i.

Theorem 3.1.4 [17] A system E:S1 -- S2, where Si C S(Rm) and S2 C S(RP),

has well-defined latency if and only if there exists an integer A such that the system

D-'E: Si -+ D-A[S2] is a causal system.

The system E: S1 -* S2, where S1 C S(Rm) and S2 C S(RP), is a recursive

system with a recursive representation yk+,+l = f(yk,. Yk+,, Uk,... Uk+,). By
restricting the range of E to the image of E, we obtain the map Er: S1 E[S1]

which is evidently surjective and possesses a right inverse E*: E[S1] -- Si such that

EE*y = y for y E [S1]. From Theorem 2.2.1, the system E* is recursive with a

recursive representation given by

Uk+, = g(Uk, ,Uk+p-1, Yk, ,, Yk+g+,4+1). (2.5)

Examining (2.5), 7 + 1 shifts are required on the output of a right inverse system E*
to guarantee causality, i.e., L(D(n+I)E*) > 0.

Lemma 3.1.5 Let E: S1 -- S2, where S1 C S(sRm) and S2 C S(RP), be a recursive

system with a recursive representation yk++1i = f(Yk,* Yk+, k,,.. Uk+,). Let

E*: E[S1] -+ Si be a right inverse of E by restricting the range of E to the image of

E with a recursive representation uk+g = g(uk,... Uk+~ yk,..* yk+p+r+1). Under
this hypothesis the system E has a latency of at most y + 1, i.e., (D(7+l)rE*) > 0.

3.2 Systems with the Lipschitz Norm

In the previous chapter we discussed two norms: the e-norm and the weight-
ed ~-norm. We will need to define a new norm for a set of systems that has
the convergence property which is that every Cauchy sequence in a compact metric
space (X, p) converges to some point of X. The Lipschitz norm for a set of systems
is introduced here because it has the necessary convergence property. The Lipschitz
norm will be derived from a semi-norm. The semi-norm definition from [41] is given
below.

Definition 3.2.1 A semi-norm on a vector space X is a real-valued function p on
X for all x and y in X and all scalars a such that

1. p(x + y) 5 p(x) + p(y)

2. p(ax) = Ialp(x)

3. p(x) 0 if x $ 0.

In this section, we let (S(Rm), p) and (S(RP), p) be two metric spaces over the
set of all sequences of m-dimensional and p-dimensional real vectors. Let B be the
set of all stable, causal, and recursive systems E: D -- S(RP) where D is a bounded
subset of S(Rm). Introduce an operator I I: B -+ W+ defined by

| de_ sf sup p((u) (2)) (3.2)
|S| = sup (3.2)
ul,uz26D,uiU2 P (u U2)

The number ||jE| has the following properties.

Lemma 3.2.2 If E, T E B and a E 3, then

1. I|E1 = 0 if and only if E is a constant system on TD;

2. IlaE| = lal II~EI; and

3. |ii + T 1| I< | |11 + I11| |.

The number ||lE| is called the Lipschitz semi-norm ([7]) of the system E on D.
A semi-norm has the property that I|E|l = 0 does not necessarily imply that E = 0.
In fact, it can be easily seen that flE|| = 0 if and only if E is a constant system (need
not be zero) that maps all sequences from D to the same sequence in S(RP). The
Lipschitz semi-norm is like a derivative bound or gain for all systems in B. Let's
define Lip(D, S(RP)) as the subset of B satisfying IIEIl < 00.
It is clear that an element E of B is in Lip(D, S(RP)) if and only if there is a
number L > 0 such that

p (E(ul) E(u2)) < Lp (u u2)

for all ui, u2 E ). Moreover, ISEll is the "least" such number L. It is also evident
that a system with the Lipschitz semi-norm is both bounded and continuous on its
domain. The Lipschitz semi-norm will be used to define the boundary of a set from
which an approximate right inverse can be chosen. The semi-norm I | can be made
into a norm as seen in the following theorem.

Theorem 3.2.3 Let uo be an element of D and let

def
IIIILE P p( (uo)) + I|II

Sp (E(uo))+ sup (E 2), (3.3)
u1,U2 V,U1 U2 P (Ul U2)
then the number IIE|ILip defines a norm for all E E Lip(D, S(RP)).

IIEI Lip will be called the Lipschitz norm ([7]) of the system E defined by
uo E D. A convenient choice of u0 is of course u0 = 0 if 0 E D, where note that E(0)
is not zero in general. To prove the theorem, it amounts to showing that |IIIIL; = 0
implies E = 0, the zero system. This, however, is an immediate consequence of
part 1. of Lemma 3.2.2.

Theorem 3.2.4 [7] The family Lip(D, S(RP)) is a complete metric space under the
Lipschitz norm 1l| lLip.

A complete metric space (Lip(D, S(RP)), IILip,) has the property that every

Cauchy sequence of systems {Ek} converges.

3.3 Approximate Right Inverse

In [50], the concept of an approximate inverse was formulated for linear sys-

tems. The problem of disturbance attenuation was addressed with a configuration

very similar to Figure 1.1. As a result of [50], the approximate invertibility of the

system was shown to be a necessary and sufficient condition for disturbance attenu-

ation. The optimization problem from [50] used the t' -norm II |1o and a weighted

semi-norm || 1w with the property II| w < I |loo. The weighted semi-norm has

a weighing filter where W will denote a stable and causal system of unit norm such

that IIdllw = IIWdlloo.

The definition from [50] of an approximate inverse is as follows. The set of all

stable, causal, and linear systems is denoted by Cs. For any stable and strictly causal

system P, an approximate right inverse of P is any stable and causal system TP for

which III PklIIw < Il11iw. The right singularity measure of P, denoted by /(P), is

y(P) inf II P|PIw. (3.4)
TECs

In general, P(P) is a number in the interval 0 < P(P) < 1.

For nonlinear systems we cannot define a weighted semi-norm with the with

the property II I1w < II* |oo. As a result, we cannot use (3.4) to derive a right

singularity measure of E between 0 and 1. Nevertheless we can define an approximate

inverse optimization problem. Let's first review some properties of right inverses for

nonlinear systems. In Section 2.2 the generalized right inverse was introduced for a

class of recursive systems in the form of (2.1). For systems E that are not bijective,

the generalized right inverse of systems E are non-unique. For the systems that are

bijective, the generalized right inverse is stable but not necessarily causal.

By assumption, the system E is strictly causal. This implies (E) > 1. The

design integer A, where A > 0, is used to characterize the class of approximate right

inverse systems that takes in account the latency of E. The design integer A is

selected before the optimization process. In the three examples of Chapter 4, the

design integer was selected at the minimum latency of E, i.e., when (E) > Ax, the

design integer A = A1.

The approximate right inverse optimization problem of the stable, strictly

causal and recursive system E: D -- S(RP) with D, a compact subset of S(Rm), is

shown in Figure 3.1 with e = [I(u) D-^EI(u)] where e E S(RP). The performance

index is

inf sup p(e) (3.5)
9'EA uES(aP)
where A is the set of all stable and causal systems T: S(ap) -- S(Rm). The system T

is a stable approximate right inverse of E while the best approximant T* is not guar-

anteed to be stable. To guarantee stability of the best approximant T*, the systems

E and IF are chosen from a class of systems Lip(D, S(RP)) and Lip(S(aP), S(Rm)),

respectively. The set A has to be compact subset of Lip(S(aP), S(Rm)). The next two

Figure 3.1: The approximate right inverse optimization problem of E: inf sup p(e).

paragraphs will explain why we need the set A to be compact to guarantee stability

of the best approximant a*.

Approximation theory is needed to solve the problem of finding a best ap-

proximant I*. In solving the best approximation problem, additional hypotheses are

placed to guarantee the existence and uniqueness of the best approximation. Our

practical problem is that we need approximate inverses that are stable, causal, and

recursive but for now we will take approximate inverses that are only stable and

causal. We shall be therefore be interested in a basic existence Theorem 3.3.1, which

gives sufficient conditions to guarantee existence of closest point.

Theorem 3.3.1 [6] Let K denote a compact set in a metric space. To each point p

of the space there corresponds a point in K of minimum distance from p.

In the previous section, the class of systems with the Lipschitz norm was

introduced so that we can create a compact set of systems in a metric space. Using

the hypothesis of Theorem 3.3.1 as a starting point, we can then prove the existence

of best approximate right inverse.

3.3.1 The Existence Theorem of Best Approximate Right Inverse

For clarity, the following sets are defined for this subsection only. Let a > 0 be

a real number. Denote by Lip(S(aP), S(R )) the set of all stable and casual systems

:. S(aP) -+ S(Rm) satisfying I1 11 < o0. Let D be a compact subset of S(Rm) and

let Lip(D, S(RP)) be the set of all stable and causal systems E: D -- S(RP) satisfying

IIEII < oo. Let G > 0 be a real number. We denote by AG the compact subset

Lip(S(aP), S(Rm)) given by

AG = {T1T E Lip(S(aP),S(Rm)) and |I11 < G}. (3.6)

Theorem 3.3.2 Let E E Lip(D, S(RP)) be a strictly causal and recursive system for

a given real number G > 0 and for a design integer A > 0. There is a stable and

causal system V* E AG such that

sup p (I(u) D- (u)) < sup p ((u) D (u)) (3.7)
uES(aP) u6S(aP)

for all E AG-

Proof:The right hand side of (3.7) is a continuous function of T. If we take T = 0,

then supES(OC) p(I(u) D-^EI(u)) = a. There exists an infimum of the right hand

side of (3.7) since it is continuous with respect to T and it is bounded between 0 and

a. Let

r = inf sup p (I(u) D-AE(u)) (3.8)
IEAG uES(aP)
and suppose that

sup p(I(u)-D- D Ek(u)) =rk, k=0,1,2,..., (3.9)
uES(aP)

with rk monotonically decreasing to r. Then, for k sufficiently large

sup p (I(u)- D-AEk(u)) < r + 1 < M. (3.10)
uES(aP)

Since the sequences of systems {I D- Ek} are bounded and monotonic, it follows

that the sequence of systems {I D-A Ek) converges. This implies that for every

e > 0 there exists an integer N such that I, n > N implies

sup ([1(u) -D-C (u)]- [I(u) D- i(u)]) < e (3.11)
uES(aP)

or

sup p (D~ E (u) D- A J(u)) < e. (3.12)
uES(aP)
The Lip(S(aP), S(RP)) systems set is a subset of the stable and causal sys-

tems F: S(a) -+ S(RP) satisfying IrF|l < oo. The composition of two systems of

bounded Lipschitz norm is also a system of bounded Lipschitz norm; hence, for

all E Lip(S(aP),S(Rm)) and a given EE Lip(D, S(RP)), the composite system

ET C Lip(S(aP), S(RP)). Now we define fk Lip(S(aP), S(RP)) such that fk dW Ek
and rewriting (3.12) we have

IID-Xf, D- fmIILip < (3.13)

This establishes limk-oo D-Afk = D-Af* where the system D- f* is the limit of the
sequence of systems {D-fk}. To translate D-Xf* into a Lip(S(aP), S(RP)) system
we get DA(D-Af*) = f*]}. The system f*]o has its first A outputs yo, yl,...,y,-i
set to zero because they were not used in the convergence of the sequence {D-Xfk}. It
follows by the completeness of Lip(S(aP), S(RP)) that the system f*] E Lip(S(aP),

s(RP).
The sequence of systems {D- fk} converges to the limit system D-Af*. There-
fore, the sequence of systems {D-AXEY k} converges also to the limit system D-Xf*.
In other words, for every e > 0 and for all u E S(am), there exists an integer N such
that 1, n > N, this implies

sup p (D-AE,(u) D- APEII(u)) < (3.14)
uES(aP)
From the properties of a system of bounded Lipschitz norm we can we write

sup p (D-XEV,(u) D-AEXl(u))
uES(aP)
< IID-'IILip sup p('I(u)- ''(u)). (3.15)
ueS(aP)
We will assume IID-^EIILip 0, which implies that E is not the zero system. The
sequence of systems {Ik} is a continuous mapping from a compact metric space

(S(aP), p) into a compact metric space (D, p). Now the sequence of systems { k} is
a sequence in a compact metric space (AG, I1 IILip). Therefore some subsequence of

{ k} converges to a system in AG. Thus we have limjoo, kj = T* where I* AG
for the subsequence of systems { k, }.
Now we can define a subsequence of {D- fk} by {D- fk,} such that
D-Afk, df D-^ qk,. Every subsequence of {D-Afk} converges to D-'f* so

limj-oo D-^fk, = D- f*; therefore

lim D-^ EI = D- f*. (3.16)

The system T* E AG is best approximate right inverse of the system E with optimal

sequence of systems { kj } E AG that converges to T*. U

3.4 The Measure of Right Singularity

As has been previously noted, for every pair of elements u, v E S(Rm), the

metric p, is defined

p,(u,v) d sup (1 + e)-i ui vil, > 0. (3.17)
i>o
The distance from any u E S("Rm) and the set S2 C S(Rm) is defined by the number

M, (u, S2) d' inf p,(u,v). (3.18)
vES2
We can now define the maximal distance from any u E S1 C S(Rm) and the set S2

by the number

e($Sl, S2) = sup eM,(u,S2) = sup inf p,(u,v). (3.19)
uES1 uES1 VES2
From Theorem 3.3.2, there exists a best approximate right inverse system

T* E AG for a design integer A > 0 such that the inequality

sup p, (I(u) D-A^E (u)) < sup p, (I(u) D-AEY (u)) (3.20)
uES(aP) uES(aP)

Figure 3.2: The performance bound optimization problem of E: supues(ap) pc(e).

holds for all T E AG and for all c > 0. Using T* for I with e = [I(u) D- E *(u)]
in Figure 3.1, Figure 3.2 represents the performance bound optimization problem of
the causal and recursive system E E Lip(D, S(RP)), where the performance index is

sup p,(e) (3.21)
uES(0aP)
where e E S(RP). Now we want to put a lower bound on the term supuEs(aP) p,(e)

supUEs(aP) p~(I(u) D-^AE*(u)). Using (3.18) and (3.19) we can write

sup p, (I(u)- D-AIE*(u)) = sup eM: (u,D-AEJ*[S(aP)])
uES(caP) uES(aP)
= Of (S(aP), D-AE*[S(aP)])

> (S(aP), D-ImE). (3.22)

Note that even if Im E is unbounded, e (S(aP), D-D^ImE) is bounded since the first
operand of (-, .), S(co), is bounded.

Definition 3.4.1 For any system E e Lip(D, S(~P)) that is a strictly causal recur-
sive of the system (2.1), the right singularity measure, denoted by P(E), for a real

a > 0 and a design integer A > 0 is

P(E) def O, (S(aP), D- Im E) (3.23)

The performance bound is finally defined as the right singularity measure P(E).
It provides a measure of the ability to match a set by subset of the image of a system
as seen in (3.23). The performance bound came from the requirement to provide an
estimate of the minimal effect of the disturbance signal d on the output signal z of
Figure 1.1. By assumption the system E: D -+ S(RP) is stable where D is a compact
subset of S(Rm). In that case, the right coprime factorization of E = PQ-1 can be

taken as P = E and Q = I where I: D -- D, the identity system, and P: D -- E[D].

When the system E: S(am) -+ S(RP) is unstable, the right coprime factoriza-
tion of E = PQ-1 is taken as P: S E[S(am")] and Q: S -- S(am) where S C S(Rm)

and Q is bicausal. From Theorem 2.1.5, any stable and bicausal system M: S -- S

with a stable inverse system M-l: S S can be selected, so there exists a pair of sta-

ble systems A: S(RP) -- S(Rm) and B: S(am) -- S(Rm) such that AP + BQ = M.

The stabilization of E can be seen in Figure 2.1 with r = B-1 and p = A. The

closed loop response of the stabilization is PM-1. Now M can be selected as the

identity system I: S --+ S. The factorization space S can be taken as S(am) and [21]

describes the construction of the stabilizing controllers that yield the factorization

space of S(am). As a result the closed loop response is P: S(am) --+ E[S(am)].

The preceding paragraphs illustrate that the "numerator" system P directly

drives the performance bound by the image of its system E. The image of E for

discrete-time nonlinear systems is analogous to the right half-plane zeros of linear

continuous-time systems. Both of these quantities are fixed and cannot be altered by

feedback. Now we can use the performance bound for systems of practical origin.

CHAPTER 4
CALCULATION OF THE PERFORMANCE BOUND

4.1 The Estimate of the Performance Bound

It turns out that the calculation of the right singularity measure P(E) for
nonlinear systems is rather laborious. We develop in this section an estimate P(E)
of P(S) which is easier to calculate. This estimate satisfies

P(C) < Pc).

(4.1)

From Theorem 3.3.2, there exists a best approximate right inverse system T* E AG
such that the inequality

sup P. (I(u) D-~E*(u)) < sup p, (I(u)- D-A^ (u)) (3.20)
uES(ceP) ueS(aP)

holds for all E AG and for a design integer A > 0. Theorem 3.3.2 demonstrates
the existence of a best approximate right inverse system but does not provide the
construction of a best approximate right inverse system.

Figure 4.1: The estimate of the performance bound optimization problem of the
system E with an approximate right inverse system A: 'P() = supUES(,) pM(e).

Using a causal approximate right inverse system I E Lip(S(ac), S(R3m)) for a
given design integer A > 0 from Figure 3.1, Figure 4.1, with e = [I(u) D-X ET(u)],
represents the estimate of the performance bound optimization problem of the causal
and recursive system E E Lip(D, S(RP)), where the performance index is

sup pe(e) (4.2)
ueS(ceP)
where e S(RP). Now we want to express the term supues(Op) Pe(e) = supUEs(op) PI

(I(u) D-^ET(u)) using the operator g,(, .). Using (3.18) and (3.19) we can write

sup p, (I(u) D- EI (u)) = sup eM (u, D-'~EI[S(aP)])
uES(aP) uES(aP)
= g (S(aP), D- A [S()aP]). (4.3)

Definition 4.1.1 For any system E, E Lip(D, S(RP)) that is strictly causal and re-
cursive, an estimate of the performance bound, denoted by P(E), given real a > 0
and design integer A > 0, is given by

P(E) de (S(a), D-E[S(P)]), (4.4)

where T E Lip(S(aP), S(Rm)) is a causal approximate right inverse system.

4.2 Practical Implementation Techniques

For the applications presented here, we will be interested in the response of
the systems over a finite interval of time, say for a discrete set [0, T], where T > 0

is a fixed integer. Let CT(R") be the set of all functions h: [0, T] --+ ". Denote by

hi| supiE[o,T] Ih(i) the e"-norm on CTR(") and by pe(h) f supie[o,T](1+e)-ih(i)I,
for e > 0, the norm p, on CT(Rn). Let h CT(Rn) be any function. Then

sup (1 + e)-~ih(i)l sup Ih(i)l
iE[O,T] iE[O,T]
and

(1 + e)T sup (1 + e)-'h(i)l = sup (1 + e)T- lh(i)l > sup Ih(i)|,
iE[O,T] ie[O,T] ie[O,T]

so that

p,(h) < Ihl < (1 + e)Tp,(h) and (1 + e)-T IhI p,(h) < Ihl (4.5)

on CT(Rn). These inequalities imply the equivalence of the norms I I and p,(.) on
CT(R"). We will use the f"-norm for its simplicity in the calculation of the estimate
of the performance bound P(E() (so the norms I I and p,(.) are interchangeable over

finite intervals of time).

For the mathematical analysis presented in Chapter 2 and Chapter 3, the

signals of interest are always bounded. For implementation purposes, the maxi-

mum peak deviation of the input signal, maxi>o Auij where Au, = ui+l u, is

limited. To limit the maximum peak deviation, the signals are transformed by using

a discrete-time lowpass filter. A new family of sequences, denoted by Snp,(a), is

formed by using a second order discrete-time lowpass filter with a cutoff frequency

QLP in rad/sec on S(a). The transformation between S(a) and SnLp(ao) is bijective
since the lowpass filter is one-to-one and the image of the lowpass filter is SnOL(a).

Thus, the transformed space SoLp(a) is isomorphic to the input space S(a).

The alternate definitions of P(E) are given in (4.3) and using the -norm

(for e = 0) we get

P(E) = sup I(u) D-'~E(u) = sup eMo (u, D-AEk[S(a)]) (4.6)
uES(oP) uES(aP)

II +

LPF --

I C -D-

S F 4 T Fuct----- ----S --- ---- --- -
u E S(e P)T e R+
Figure 4.2: The Function f: S(ac) *+: u H f(u).

Define a function f: S(aP) R+ that incorporates the lowpass filter, as displayed in

Figure 4.2, such that

P(E)= sup f(u), (4.7)
UES(ceP)
where the function f is defined as follows

f(u) J I (LPF(u)) D- EW (LPF(u)) (4.8)

The parameters of the function f are the system E, an approximate inverse system

1P, the design integer A, the optimization set S(aP), and the lowpass filter cutoff

frequency QLP. The function f must also satisfy f(u) = gMo(u, D-'EII[SOLP(cP)])

where QM (', ) is defined in (3.18). It then follows that f(.) is a convex function when

its input or optimization set S( P) and the image set ETS[SOnLP(a)] are convex since

eM0 is a distance function from a point in the optimization set to the image set. The
distance function peM by its construction is a convex function. In practice, the image

set is not convex; however, studying the optimization problem of finding a global

maximum of a convex function will give us insight to the non-convex optimization

problem.

Unlike the minimization of a convex function, once a local maximum has

been found there is, more or less by definition, no local information to tell you how

to proceed to a higher local maximum. In particular, there is no local criterion

for deciding whether a given local maximum is really the global maximum. Hence,

maximizing a convex function is usually a much harder task than minimizing a convex

function. Theorem 4.2.1 shows that a convex function achieves a maximum over a

compact polyhedral set at an extreme point.

Theorem 4.2.1 [3] Let g: S(RP) -, R be a convex function, and let S be a nonempty

compact polyhedral set in S(3P). For the problem of maximizing g(u) subject to u E S,

there is an optimal solution, u, where u is an extreme point ofS such that g(u) > g(u)

for all u e S.

Using Theorem 4.2.1, he global maximum of f(u) maybe found by evaluating

all the extreme points of S(ap) and finding the maximum on this set. For example,

if we have a one-dimensional set of sequences which is a finite length of 250 elements,

the number of extreme points, and hence the number of evaluations of f, would be

2250 = 1.8903 x 1075. This procedure is laborious as far as computation is concerned.

Although the image set S is non-convex, the property that the optimization set is

convex is enough to suspect that global maximum would be found on an extreme

point of the optimization set. If we take an interior point ui from the optimization

set which is a nonempty compact polyhedral set and form a neighborhood N,,(ui)

of radius ri, then we can find a Ui+1 E Nr,(ui) such that gMo(ui+i, S) > eMo(ui, S)

and ui+l is closer to the boundary than ui. The radius r, of a neighborhood Nr(ut)

is selected so that we do not accumulate in an interior local maximum. Since the

optimization set is convex and eMo is a distance function, the sequence {ui} will

converge to a point on the boundary. Therefore, the global maximum will be on

the boundary. For simulation purposes, the global maximum is assumed to be an

extreme point.

The three systems of practical origin that have been selected for investigation

are all fourth order nonlinear systems. Each of the systems has a distinct memoryless

nonlinear subsystem. These subsystems represent, respectively, dead-zones, position

constraints, and rate constraints for the problems: permanent magnet stepper motor

model example, multivariable process control model example, and aerodynamic model

example. Nonlinear programming is used to calculate the estimate of the performance

bound. The algorithms that are developed take advantage of the traits of each of

the nonlinear systems and of the finite bandwidth of the input signal space. A useful

trait of a system is a property that would help reduce the region of optimization

in the search for a global maximum. It will be shown in the next sections that the

computation is not as laborious as previously conjectured.

4.3 Permanent Magnet Stepper Motor Model Example

This section is concerned with voltage control of a permanent magnet (PM)

stepper motor. This example comes from [5], wherein the authors developed a model-

based control law using exact linearization implemented on an industrial setup.

The PM stepper motor model consists of a slotted stator with two phases

and a permanent magnet rotor. One side of the rotor is a north pole, and the other

side is a south pole. The teeth on each side of the rotor are out of alignment by a

tooth-width. The equations of motion for the PM stepper motor model are given by:

di 1
d- -[va Ri, + Kmw sin(Nr)] (4.9)
dt L
dib 1
S[b [vb-Rib Kmwcos(NrO)] (4.10)
dt L

[-Kmia sin(NO) + Kmib cos(N,) Bw] (4.11)
dt
d w (4.12)

where Va, Vb and ia, ib are the voltages and currents in phases A and B, respectively.

Further, w is the rotor (angular) speed, 0 is the rotor (angular) position, B is the

viscous friction coefficient, J is the combined inertia of the motor and the translation

table, Km is the motor torque constant, R is the resistance of the phase winding, L is

the inductance of the phase winding, and N, is the number of rotor teeth. Table 4.1

shows the values of the stepper motor parameters used in [5].

Table 4.1: The PM Stepper Motor Model Parameters
Parameter Value
L 1.5 mH
R 0.55 Q
J 4.5 x 10-5 kg-m2
Km 0.19 N-m/A
N, 50
B 8.0 x 10- N-m-s/rad

For the PM stepper motor model, there is an appropriate nonlinear coordinate

transformation which is known as the direct-quadrature (DQ) transformation. The

DQ transformation for the phase voltages and currents is defined as follows:

Vd cosdef CO(NrO) sin(Nr ) Va (4.13)
Vq -sin(NO) cos(NO) Vb
id idef cos(N ) Sin(Nr ) a(4.
q] -sin(Nr,) cos(NO,) j ib

The direct current id corresponds to the component of the stator magnetic field

along the axis of the rotor magnetic field, while the quadrature current iq corresponds

to the orthogonal component. The application of the DQ transformation to the

original system equations (4.9) through (4.12) yields:

did 1
= [Vd Rid+ NrLiq] (4.15)
dt L
di 1
= [v,- Ri,- NrLwid Kmw] (4.16)
dw 1
d =- [Kiq Bw] (4.17)
dt J
= (4.18)
dt

where Vd is the direct voltage, v, is the quadrature voltage, id is the direct current,

iq is the quadrature current, w is the angular velocity, and 0 is the angular position.

The quadrature component iq of the current produces torque while the direct

component id does not produce any torque. However, in order to attain higher speeds,

it is necessary to apply a negative direct current in order to cancel the effect of the

back-emf of the motor which results from the "back-emf" term Kmw in (4.16). The

direct current id will be forced to be negative by the correct choice direct voltage Vd.

The optimal lead-angle approach is used in [5] to properly choose vd which becomes

a nonlinear function of the states. The equations of motion are discretized using first

order Euler integration. Let Xk = [id iq k Ok k]'. Equations (4.15) through (4.18)

take the following form:

Xk+1 = f(xk) + g(xk) vq (4.19)

Wk = [0 0 1 0] k, (4.20)

where xk E R4, and f, g are vector fields on R4.
The analytical expression for the input dead-zone is

VCqk 1 if VCk > 1
Vqk =DZ(vcqk)= 0 if -1 < cq, < 1 (4.21)
vcqk + 1 if vcqk < -1.
The system E: D -- S(R) : vcq '- w, where V C S(R), is the composition of
the equations of motion with the dead-zone for which vcq is the quadrature voltage
command. Figure 4.3 shows the block diagram of the system E.

4.3.1 Simulation

The signal of interest is the angular velocity signal w. It is assumed that all
signals of interest start at rest (i.e., w, = 0 for all n < no). In one of the applications
presented in [5], the desired speed trajectory required the motor to go from to zero
to 1350 RPM in 10 ms with a sample rate of 10 kHz. The bandwidth of w was set to
50 Hz. The rise time was about 10.6 ms and the absolute bound of Iwl < 1350 RPM
was imposed.

Dead-Zone

V qkv k Xk+1 =f() + g(k) vq k k+1

ck = [0 0 10] X k

Figure 4.3: The block diagram of the system E: D -+ S(R) : vcq w for the PM
Stepper Motor Model Example.

The block diagram of an approximate right inverse system T: S2,.so(a) --

S(R) : we t Vqinv is shown in Figure 4.4. The first block calculates the raw quadra-

ture voltage vk for a given angular velocity command wek. The analytical expression

for the continuous dead-zone inverse (CDI) is given by

V1 + 1 if v' > 1
qk qk G-
vqinvy -= CDI(v) = G- v if <, < (4.22)
Vqnv qk qk G-1 qk G-i
-1 if V < 1
qk qk G-1

Here G = IICDI(.-)I represents the Lipschitz semi-norm of the CDI system. The

range of G is 1 < G < oo. In the CDI block, the voltage quadrature command

Vqinvk is determined by using (4.22) on the raw quadrature voltage v' The voltage

quadrature command Vqinvk is then used to update the states xinvk of the system T.

Two integration of Vq are required to reach w. This results in the discretized

model of at least a two step delay. Therefore, (E) > 2. From Figure 4.4 it follows

that (T) > 0. Theorem 3.1.2 applied to the composite system CE: S2,.5o(a) -+

S(R) : we '- w yields C(EI) > 2. For the simulation, g(xk) = [0 !At 0 0]' # 0 in

(4.19) for all vcq E D and a given x0. The design integer is taken as A = 2 for an

approximate right inverse system F.

For the simulation of 40 milli-seconds, a sampling period of 100 micro-seconds

yields a sequence of 400 points. In the previous section, the global maximum was

-.----.--.--.---- .----------- ---- -.-__ -_ _- __---------- -- _--- - ----- -----, -_ --_ -_- _-_--- -

k X= f(k-l) + g( k-1)V k Continuous qinvk
q 'k Dead-Zone Inverse
k = [0 0 1 0]
kb

Xk- Dead-Zone

z Xk kx ) k-.)) vq k vq k /

L-.---- -.-.- .- ,----- - -----------_. ----_ -. ---------__ ---. -_----. ._- .-- -- ..- .

Figure 4.4: The block diagram of an approximate right inverse system il: S2,.so(a) --*
S(R) : we -H vqinv for the PM Stepper Motor Model Example.

shown to be on the boundary. For simulation purposes, the global maximum is

assumed to be an extreme point. The pre-filtered binary random sequences have

only the values 1350 RPM, the extreme points of S(1350). Table 4.2 shows the

Monte Carlo results for 4 different system gains G of the CDI system: 2, 10, 100

and 1000, where the error Iwe D-2 E(wc)I is evaluated for each of the sequences.

As can be seen from Table 4.2, as G oo, the error wec D-2'E(wc)l 0 in

a linear fashion. For the limiting case as G -- oo, an approximate right inverse

system T is not continuous. However, the performance bound P(E) appears to be

directly proportional to the system gain G. This observation will be proven true in

the following subsection.

From Table 4.2, the maximum value of wc D-2 E(wc)l, for G = 10 rad/sec,

is 0.02686. This corresponds to the 18th sequence in the Monte Carlo run. In the

following subsection we will calculate the estimate of the performance bound P(E)

over the optimization set S21.50(1350). A first estimate is P(E) > 0.02686. In

Figure 4.5, the 18th sequence for we, the filtered square wave, is plotted with the

response sequence EC(wc), without the D-2 shift since the errors are less than 0.01%

of the magnitude of the signal. The absolute angular velocity error sequence is plotted

in Figure 4.6. The quadrature voltage vq and the absolute quadrature voltage error,

Iv, -Vqk, are plotted in Figure 4.7 and Figure 4.8, respectively. Note, that Figure 4.6
and Figure 4.8 appear to be the same figure but of different scales.

4.3.2 The Estimate of the Performance Bound P(S)

In this section, we will find an estimate of the performance bound P(S) for the

voltage control problem of a PM stepper motor by using the optimization definition

as follows:

P() = sup wc D-'2 (wc) (4.23)
C ES2,.5o(1350)

Table 4.2: Family of Sequences, S2,.50(1350), Comparison
Approximate Right Inverse Error Iwc D-2 (wc)l
Sequence Gain G = 2 Gain G = 10 Gain G = 100 Gain G = 1000
1 0.1339 0.02628 0.001151 0.0002244
2 0.1341 0.02588 0.001021 0.0000000
3 0.1342 0.02510 0.001979 0.0001157
4 0.1343 0.02546 0.001346 0.0000000
5 0.1332 0.02623 0.001882 0.0000000
6 0.1336 0.02631 0.000332 0.0000000
7 0.1343 0.02652 0.002509 0.0000000
8 0.1340 0.02608 0.000279 0.0000000
9 0.1341 0.02659 0.001415 0.0000000
10 0.1326 0.02628 0.001764 0.0000000
11 0.1343 0.02641 0.001312 0.0002664
12 0.1339 0.02650 0.002210 0.0000549
13 0.1344 0.02617 0.000710 0.0000000
14 0.1332 0.02537 0.002150 0.0000000
15 0.1334 0.02670 0.002611 0.0000000
16 0.1328 0.02683 0.002629 0.0000000
17 0.1334 0.02662 0.002360 0.0001910
18 0.1336 0.02686 0.001415 0.0000000
19 0.1341 0.02619 0.002322 0.0000000
20 0.1341 0.02636 0.002085 0.0000224
21 0.1337 0.02680 0.002059 0.0000000
22 0.1342 0.02615 0.002300 0.0002322
23 0.1337 0.02464 0.000000 0.0000000
24 0.1343 0.02566 0.001393 0.0000000
25 0.1342 0.02677 0.002103 0.0002603
max 0.1344 0.02686 0.002629 0.0002664

Figure 4.5:
S27.50(1350)
Carlo run.

0 5 10 15 20 25 30 35 40
Time (milli-seconds)

The angular velocity response w = -TI(wc) and the command wc C
with system gain G = 10 are shown for the 18th sequence in the Monte

0.031

0.025 F

0.02 1

0.015k

0.005 ....

n

0 5 10 15 20 25 30
Time (milli-seconds) [dt=0.1 msec]

35 40

Figure 4.6: The absolute angular velocity error, IWCk D-2 E'(W k)l is shown for the
18th sequence in the Monte Carlo run of qc e S2,.50(1350) with system gain G = 10
where Iwe D-2T(wc)| = 0.02686.

15 20 2!
Time (milli-seconds)

Figure 4.7: The quadrature voltage
Carlo run of we E S2-.5o(1350) with

Vq is shown for the 18th sequence
system gain G = 10.

40

in the Monte

0.09 H

-0.08

S0.07.
w
(D
0 0.06 ..

0 0 .0 5 ..... . ...... .
0 0

6 .
S0.02

0.01
0.02 ;;

0.01

0 5 10 15 20 25 30
Time (milli-seconds) [dt=0.1 msec]

35 40

Figure 4.8: The absolute quadrature voltage error, Iv' Vqk is shown for the 18th
sequence in the Monte Carlo run of we E S2..5o(1350) with system gain G = 10 where
1v; v = 0.0993.

I 1 I 1

where the system E and an approximate right inverse system I are illustrated in

Figure 4.3 and Figure 4.4, respectively.

The first block of Figure 4.4 calculates the raw quadrature voltage v' for a

given angular velocity command wc,. In the CDI block, the voltage quadrature

command Vqinv is determined by using (4.22) on the raw quadrature voltage v' In

Figure 4.3 (vcqk = Vqinvk), the quadrature voltage vqk is the output of the DZ block

by using (4.21). The composition of (4.22) and (4.21) yields

V' if v' > 1
qk qk G-1
Gv' -1 if 1< <' 1
Gqk G qk G-1
v = DZ CDI v )) = 0 if -1 < 1 (4.24)
G v' +1 if < V <-
qk G-1 qk G
v' if v' < -
qk qk G-1

It follows then that for v' = -
1k G'

sup v Vq = sup V' DZ (CDIv)I = 1 (4.25)
v'qs(E) v'ES()

For G = 10, = 0.1, and this is consistent with Figure 4.8 where |v' Vq,

0.09993 < 1
-- G
The maximum quadrature voltage error occurs when the magnitude of the raw

quadrature voltage v' k = Since there is a feedback loop inside an approximate

inverse system T, the quadrature voltage errors are not accumulative with respect to

the angular velocity errors. Evaluating the right hand side of (4.23) with Iv'k = -

will give us the estimate of the performance bound which requires two integration

of the errors as follows:

1 1
at IV
SAvk I G a I G
I' Zqk = G at l =

60 Km At2
IAq2+ = 2xG- L J

or
60 Km
(E) = m At2. (4.26)
27rG-L J

Using the optimization set we E S2r.50(1350), and for G = 10, this yields P(E) =

0.02688. The estimate is 0.00002 higher than the value from Table 4.2 of 0.02686.

The estimate of the performance bound P(E) is directly proportional to the system

gain G. Figure 4.6 is a scaling of Figure 4.8 by a factor of L.j-_At2. In this example,

we have found a closed form solution for an estimate of the performance bound

given by P(E) which was verified by simulation. In the next 2 examples, nonlinear

programming algorithms will be used to calculate P(E).

4.4 Aerodynamic Model Example

This section will provide an example of the performance bound calculation for

the longitudinal control of an aircraft. Control of pitch attitude of an aircraft can

be achieved by deflecting all or portion of either a forward or aft tail surface. The

derivation of the equations of motion of the aircraft with the assumption that the

roll rate, yaw rate and y-velocity component in the body-axes frame are all zero is

found in [4]. These equations are:

1
i = [X (u, w, q, 6) + T mg sin(O)] wq (4.27)
m
1
v = [Z (u, w, q, 6) + mg cos(0)] + uq (4.28)
m
q = M(u,w,q,0,6)/Iy (4.29)

0 = q (4.30)

where u and w are the body-axis x-velocity component and the body-axis z-velocity

component, respectively. The body-axis pitch rate and the Euler pitch angle are

given by q and 0, respectively. The longitudinal control input or elevator deflection

is given by b. X and Z are the aerodynamic forces and M is aerodynamic moment.

The following are constants: T, thrust; m, mass; g, gravity; and I,,, moment of

inertia. The aerodynamic forces and moments are represented by the means of the

aerodynamic stability coefficients as seen below:

a = arctan(w/u)

a = (cos(a))2 (u wu)

Q = p(u2+w2)

L = (CL + CL, + CLqq+ CL,6) QS

D = (CD CDa) QS

X = L sin(a) D cos(a) (4.31)

Z = -Lcos(a) Dsin(a) (4.32)

M = (Cm + Cma + Cm,, + Cmqq + Cm6) QSc. (4.33)

Here a and i are the angle of attack and the angle of attack rate, respectively.

Q is the dynamic pressure. The aerodynamic moment, M, and the aerodynamic

forces, X and Z, are directly proportional to t he dynamic pressure. L and D are the

aerodynamic forces in the wind axes which are known as lift and drag The following

are constants: S, the wing reference area and c, the wing mean aerodynamic cord.

Table 4.3: The Aerodynamic Model Parameters
Parameter Value Parameter Value
m 1247.390 kg CD 0.05
T 1600 N CD, 0.33
g 9.81 m-s-2 CL 0.41
Iy, 4066 kg-m2 CL, 4.44
p 1.225 kg.m-3 CLq 3.8.c/(2uo)
S 17.094 m2 CL6 0.355
c 1.737 m Cm0 -0.683
uo 53.8931 m-s-1 Cm& -4.36-c/(2uo)
vo -0.0931634 m-s-1 Cmq -9.96-c/(2uo)
qo 0 rad-s-1 Cm6 -0.923
0o 0 rad
-i1 0.00128510 rad

The reference geometry, mass characteristics, and aerodynamic characteris-

tics of the general aviation airplane: NAVION were taken from Appendix B of the

reference [33]. The values of the parameters are shown in Table 4.3. The equa-
tions of motion are discretized using second order Runge-Kutta integration. Let

Xk = [Uk Vk qk Ok]'. Equations (4.27) through (4.30) take the following form:

Xk+1 = f(xk) + g(xk) 6k (4.34)

qk = [0 0 1 0] k (4.35)

where xk E R4 and f, g are vector fields on R4.
The actuator has a first-order response modelled by G(s) = 10/(s + 10) which
has rate constraints of 15 deg/sec and position constraints of 15 deg. The system
E: I -- S(R) : bc i-9 q, where VD C S(R), is the composition of the equations of
motion with the actuator dynamics for which bc is the elevator deflection command.
Figure 4.9 shows the block diagram of the system E.

4.4.1 Simulation

The signal of interest is the pitch rate signal q. It is assumed that all signals
of interest start at rest (i.e., q, = 0 for all n < no). The absolute bound of jq\ < 1.15
deg/sec was imposed so that after 10 seconds of flight, the aircraft will always be
within its parameters of good flight conditions. Typical pitch rate commands are
generated by the longitudinal autopilot. The altitude hold mode and the mach hold
mode have time constants of approximately 100 milli-seconds ([4]). This corresponds
to a restriction on bandwidth of q to be approximately 10 rad/sec.

Figure 4.9: The block diagram of the system E: D -- S(R) : 6c q for the Aerody-
namic Model Example.

The block diagram of an approximate right inverse system : Sn,,L(a) --

S(R) : qc -- 5inv is shown in Figure 4.10. The first block calculates the raw elevator

deflection 6' for a given pitch rate command qc,. In the limiting logic block, the

elevator deflection command 6inv is determined by first using the constraints on the

raw elevator deflection 6' and then using the inverse actuator dynamics. The elevator

deflection command 6inv is then used to update the states xi,, of the system i.

One step delay of 6c is required to reach q. Therefore, C(E) > 1. From Fig-

ure 4.10 it follows that (C ) > 0. Theorem 3.1.2 applied to the composite system

ED: Sp,,(a) -- S(R) : qc q yields (ET) > 1. For the considered flight condi-

tions, g(xk) 7 0 in (4.34) for all bc E D and a given x0. The design integer is taken

as A = 1 for an approximate right inverse system '.

For a simulation of 10 seconds, a sampling period of 40 milli-seconds yields

a sequence of 250 points. The pre-filtered binary random sequences have only the

values 1.15 deg/sec because we are only interested in the extreme points of S(1.15).

Table 4.4 shows the Monte Carlo results for 3 different cutoff frequencies fLP: 10,

15, and 20 rad/sec, where the error Iqc D-1El'(qc)| is evaluated for each of

the sequences. For practical purposes, as can be seen from Table 4.4, the system

': Sio(1.15) -* S(R) is a right inverse of the system E.

----------------------"------------------------------------------------- --
QCk = (x )+ g(x k1) 8k Limiting k
q = [O0 1 0] x 1 1
b
Xk-I Position Rate
i Constraints Constraints

X,= f(x,,)+ g(x ,) k- k Z I a,

Figure 4.10: The block diagram of an approximate right inverse system (P: S L, (a) ->
S(R) : qc 6inv for the Aerodynamic Model Example.

Table 4.4: Family of Sequences, SoL,(1.15), Comparison
Approximate Right Inverse Error Iqc D-1 (qc)
Sequence QLP = 10 rad/sec QLP = 15 rad/sec RLp = 20 rad/sec
1 4.996E-16 0.9075 1.353
2 4.163E-16 1.7650 2.265
3 3.331E-16 0.4226 1.553
4 2.776E-16 0.3620 1.121
5 3.608E-16 0.5654 1.592
6 3.331E-16 0.4446 1.356
7 3.331E-16 1.4340 2.153
8 2.498E-16 0.9444 1.788
9 3.331E-16 0.6539 1.896
10 4.441E-16 1.2530 1.931
11 2.220E-16 1.2270 2.121
12 5.551E-16 0.9812 1.371
13 4.996E-16 0.7821 1.471
14 5.551E-16 0.7177 1.488
15 4.441E-16 0.6425 1.351
16 4.441E-16 0.7448 1.602
17 5.551E-16 0.4841 1.454
18 4.996E-16 0.4835 1.586
19 4.441E-16 0.6658 1.585
20 3.331E-16 0.7861 1.675
21 3.608E-16 0.5702 1.544
22 3.331E-16 0.6071 1.809
23 3.886E-16 0.6047 1.635
24 3.053E-16 0.4456 1.797
25 4.441E-16 0.9205 1.413
max 5.551E-16 1.7650 2.265

1.5 ...\ .. . .... ... \ I- q
1.5 .qc

S-1
0.5-

ca
0-

-0.5-
o

1-1.5 -

-21----
0 1 2 3 4 5
Time (second,

Figure 4.11: The shifted pitch rate response q
qc E S15(1.15) are shown for the 2nd sequence in

1.81 1 1 1--

1.6

1.4

L.

Ii.,

0 1 2 3 4 5 6 7
Time (seconds) [dt=0.04 sec]

6 7 8 9 10
S)

= D-1EY(qc) and the command
the Monte Carlo run.

n A

Figure 4.12: The absolute pitch rate error, lqc, D-1E (qck)I, is shown for the 2nd
sequence in the Monte Carlo run of qc E S15(1.15) where Iqc D-1E'(qc)[ = 1.765.

I I

0.4-

0.2

Go n~nl i

(JU......... 1....hM ILA J..4. 4

8 9 10

a 0.
o
0

0-
cO)

0

(U

-1.

Figure 4.13:
Carlo run of

0 1 2 3 4 5 6 7 8 9 10
Time (seconds)

The elevator deflection 6 is shown for the 2nd sequence in the Monte
qc E S15(1.15).

4 5 (
Time (seconds)

Figure 4.14: The elevator deflection rate 6 is shown for the the 2nd sequence in the
Monte Carlo run of qc E S15(1.15).

From Table 4.4, the maximum value of Iqc D-1'E(qc)I for QLP = 15

rad/sec is 1.765. This corresponds to the 2nd sequence in the Monte Carlo run. In

the following subsection we will calculate the estimate of the performance bound P(E)

over the optimization set S15(1.15). A first estimate is P(S) > 1.765. In Figure 4.11,

the 2nd sequence of qc, the filtered square wave, is plotted with the corresponding

shifted sequence EI(qc) which is shifted one step to the left by 40 milli-seconds. The

absolute pitch rate error sequence is plotted in Figure 4.12. The elevator deflection 6

and elevator deflection rate 6 are plotted in Figure 4.13 and Figure 4.14, respectively.

As can be seen from Figure 4.12, the maximum error occurs at 8.0 seconds. Examining

the plots between 7 and 8 seconds reveals that the pitch rate response is not keeping

up with the pitch rate command. The elevator deflection is never position constrained

but is rate constrained from 7.44 to 8.12 seconds for a total of 680 milli-seconds. As

seen from this example, when we increase the bandwidth of the pitch rate command

from 10 to 15 rad/sec, the rate constraints become effective and thus the measure of

right singularity is increased.

4.4.2 The Estimate of the Performance Bound P(S)

In this section, we will find an estimate of the performance bound P(E) for

the longitudinal control problem of an aircraft by using the optimization definition

as follows:

P(E)= sup qc D-'ET(qc) (4.36)
qcES15(1.15)
where the system E and an approximate right inverse system T are illustrated in

Figure 4.9 and Figure 4.10, respectively. The discussion following Theorem 4.2.1

in Section 4.2 suggests that the maximum of the right hand side of (4.36) is found

by evaluating all extreme points and then finding the maximum. In Section 4.4.1,

25 random extreme points out of possible 2250 (1.8903 x 1075) extreme points were

evaluated with the conclusion P(E) > 1.765.

A
1r---Ni --- Nu 2

qc = 1.15

0

qc= -1.15

4Np 4---Np 2

Figure 4.15: Signal Construction's Nomenclature.

It is time to take advantage of the traits of the system. The first block of

Figure 4.10 calculates the raw elevator deflection 6' for a given pitch rate command

qck. In mathematical terms it is solving for 5' in the following equations:

S= (qck k-1) (4.37)
Mk = At

Mk = (Cm + Cmak-1 +Cmqqk-1 +

Cmba( (Uk-1, Wk-1,) + Cm6 6)Qk-1Sc. (4.38)

As can be seen from (4.38), the raw aerodynamic moment Mk' is affine with respect

to the dynamic pressure Qk-1. It turns out that it is also affine with respect to the

raw elevator deflection 65. Thus, S~ is inversely proportional to Qk-1 and decreasing

dynamic pressure causes the elevator deflection to increase. This increases the likeli-

hood for the elevator deflection to be rate constrained. Therefore, the optimization

region of (4.36) can be reduced to a subset of S15(1.15) which contains a lower than

initial dynamic pressure after a period of flight. A pitch-up command, say a com-

mand qc = 1.15 deg/sec for at least 2 seconds after initialization, will decrease the

dynamic pressure.

Two problems of maximization of convex functions are:

1. Once a local maximum has been found. There is no local information to tell

you how to proceed to a higher local maximum.

2. There is no local criterion for deciding whether a given local maximum is really

the global maximum.

These features are also applicable to the maximization of non-convex function since its

maximum for simulation purposes will be an extreme point. The signal construction

and the grid search method will be techniques used to find a higher local maximum

to get around the first problems. These two methods are explained in the following

paragraphs. After an extensive search for the signal that maximizes (4.36) and an

appropriate exiting condition, we can say that we have qualitatively found the global

maximum. This alleviates the second problem.

The signal construction will always be done in the pre-filtered space S(1.15)

with out any loss of generality since the pre-filtered space S(1.15) is isomorphic to

the filtered space S15(1.15). The baseline signal will be qc = 1.15 deg/sec for the

entire simulation of 10 seconds. The construction of the signals begins with a pulse

qc = -1.15 for Npl steps and starting at Nx steps. An optimization takes place

over N1 and Npl. The signal construction's nomenclature is shown in Figure 4.15. A

second pulse is constructed with Nu2 steps at qc = 1.15 and Np2 steps at qc = -1.15.

The optimization takes place over N2, Nu2, and Np2 with Npl fixed. This procedures

continues until the (i + 1)-th pulse optimization does not produce higher maxima.

A signal with i pulses has 2i parameters, (Ni, Npl,..., Npi, Nu2,..., Nui). In

the search grid method the starting parameter Ni (see Figure 4.15) of the pulse train

is still optimized over its range. The other 2i 1 parameters (Npl,..., Npi, Nu2,...,

Nu;) are varied over a selected range. The selected range is based on the available

o,

0
E
E
o

CL
(D

C

ca
0
c -1
(E-1

4 5
Time (seconds)

Figure 4.16: The shifted pitch rate response q = D-1E'(qc) and the command

qc S15(1.15) are shown for the first pulse optimization (N1 = 242, Npi = 4).

0 2

o2
W
4Ri 1.5

-o

0

4 5 6
Time (seconds) [dt=0.04 sec]

Figure 4.17: The absolute pitch rate error, Iqck D-'E (qc,)l, is shown for the first

pulse optimization (N1 = 242, Npl = 4) of qc S15(1.15) where \qc D-1E(qc) =

1.646.

0 1 2 3 4 5 6 7 8 9 10
Time (seconds)

Figure 4.18: The shifted pitch rate response q = D-1' (qc) and the command
qc e S15(1.15) are shown for the second pulse optimization (N2 = 222, Nu2 =7,
Np2 = 7).

3|5----i-:---i--:- i- i- - ------------ i-

2.5-

ia,

S1.5 -
C.

0.5

0

0 1 2 3 4 5 6 7 8 9 10
Time (seconds) [dt=0.04 sec]

Figure 4.19: The absolute pitch rate error, Iqck D-',!F(qc,)I, is shown for the
second pulse optimization (N2 = 222, Nu2 = 7, Np2 = 7) of qc E S15(1.15) where
Iqc D-'1T(qc)I = 2.713.

t 5 (
Time (seconds)

Figure 4.20: The shifted pitch rate response
qc E S15(1.15) are shown for the third pulse
Np3 = 7).

q = D-'ET(qc)
optimization (N3

and the command
= 195, Nu3 = 8,

0 1 2 3 4 5 6 7
Time (seconds) [dt=0.04 sec]

8 9 10

Figure 4.21: The absolute pitch rate error, Iqc, D-1E,(qck)I, is shown for the
third pulse optimization (N3 = 195, Nua = 8, Np3 = 7) of qc E S15(1.15) where
Iqc D-1'E(qc)j = 2.856.

2.

0.

.5

1-

5

2.5

2
(S

"o
Ti 1
C

0.5
E
E
o
0

8 -0.5
C
0
. -1

n" -1.5

a _O

-2.51' 1 1 1
0 1 2 3 4 5 6 7 8 9 10
Time (seconds)

Figure 4.22: The shifted pitch rate response q = D-1'ET(qc) and the command qc e
S15(1.15) are shown for the search grid method optimization (N3 = 215, Npl = 2,
Nu2 = 7, Np2 = 5, Nu3 = 7, Np3 = 8).

0 1 2 3 4 5 6
Time (seconds) [dt=0.04 sec]

7 8 9 10

Figure 4.23: The absolute pitch rate error, Iqck D-1EC(qck)l, is shown for the
search grid method optimization (N3 = 215, Np1 = 2, NU2 = 7, Np2 = 5, Nu3 = 7,
Np3 = 8) of qc E S15(1.15) where P(E) = 2.908.

ci)
c -0.5
o

12

S-1.5

0 1 2 3 4 5 6 7 8 9 10
Time (seconds)

Figure 4.24: The elevator deflection 6 is shown for the search grid method optimiza-
tion (N3 = 215, Npl = 2, Nu2 = 7, Np2 = 5, Nu3 = 7, Np3 = 8) of qc E Si5(1.15).

-15L i i I ( e o
0 1 2 3 4 5 6
Time (seconds)

7 8 9 10

Figure 4.25: The
optimization (N3
qc E S15(1.15).

elevator deflection
= 215, Npl = 2,

rate S is shown for the the
Nu2 = 7, Np2 = 5, Nu3

search grid method
= 7, Np3 = 8) of

computing power. All points in the selected range are evaluated and then the max-

imum is found. The maximum of the search grid method is qualitatively the global

maximum provided that the selected range is suitable.

The results of the nonlinear programming algorithm, which is the integration

of the signal construction and the grid search method, are shown in Table 4.5 and

Figures 4.16-4.25. Figure 4.16 and Figure 4.17 represent the first pulse optimization

with plots of the sequences qc, D-1'E(qc) and Iqc, D-1EI(qck)j. Figure 4.18

and Figure 4.19 represent the second pulse optimization with plots of the sequences

qc, D-1E'(qc) and Iqck D-1YE(qc,)l. Figure 4.20 and Figure 4.21 represent
the third pulse optimization with plots of the sequences qc, D-1E' (qc) and Iqck -

D-lE(qc,)I.

The third pulse optimization parameters are (Npi, Nu2, Np2, Nu3, Np3)

(4, 7, 7,8, 7). The search grid was taken to be (2:6,5:9,5:9,6: 10,5:9) which yields

3125 grid points that have to be optimized over the variable N3. The result of the

search grid method and the estimate of the performance bound are:

(N3, Np, N2, Np2,Nu3, Np3) = (215, 2, 7, 5, 7, 8)

P(E) = 2.908.

The search grid sequence of qc E S15(1.15) is examined. In Figure 4.22, the se-

quence qc, the filtered square wave, is plotted with the corresponding shifted se-

quence E((qc) which is shifted one step to the left by 40 milli-seconds. The absolute

pitch rate error sequence is plotted in Figure 4.23. The elevator deflection 6 and

elevator deflection rate 6 are plotted in Figure 4.24 and Figure 4.25, respectively. As

can be seen from Figure 4.23, the maximum error occurs at 9.88 seconds. Examin-

ing the plots between 8.6 and 10 seconds reveals that the pitch rate response does

not keep up with the pitch rate command. The elevator deflection is never position

constrained but is rate constrained for a total 1.2 seconds out of the last 1.4 seconds

of simulation.

Table 4.5: The Estimate of the Performance Bound P(E) for the Aerodynamic Model
Example
Optimization Method Parameters Iqc D-'a (c)j
(Ni, Np, N2 NP2, Nu3, Np3) for qc E S15(1.15)
First Pulse (242,4) 1.646
Second Pulse (222, 4,7,7) 2.713
Third Pulse (195,4,7,7,8,7) 2.856
Search Grid (215, 2, 7, 5, 7, 8) 2.908

-P(E) = 2.908

4.5 Multivariable Process Control Model Example

This section is concerned with multivariable control of a plant consisting of a

cylindrical tank containing a water column pressurized by air. This example comes

from [10, 11] wherein the authors developed a multivariable predictive control strategy

with a variable receding horizon.

The process control model consists of a cylindrical tank which contains a water

column pressurized by air. Figure 4.26 is a simplified process diagram. A centrifugal

pump feeds water into the tank through the pneumatic valve 1, and a compressor

feeds air through the pneumatic valve 2. Liquid can drain from the tank through

the manual valve 3 while air can exit through the relief valve 4. Measurements of

the process output-the height of the liquid column ya and the gauge pressure of

the enclosed air y2-are obtained by means of pressure transducers. Valves 1 and 2

are manipulated through electro-pneumatic transducers by voltage signals ul and u2.

The height of the tank is 105 cm and the diameter is 7.5 cm.

The dynamics of the process control model can be described in terms of four

states. Let xz denote the liquid-level height [m], x2 the air pressure inside the column

valve 2 2
Air
Air

U2 ..................Vr ........ y

valve 1 valve 3

Figure 4.26: The Multivariable Process Control Model with inputs: valve signals ul
and U2; with outputs: liquid level yl and the air pressure Y2*

[bar gauge], x3 the opening of water valve 1 (lift) dimensionlesss], and x4 the opening

of air valve 2 (lift) dimensionlesss]. The equations of motion of the process control

model are given by:

X1 = (qf qid) (4.39)
1 + 2
x2 -= S( ) [(qif qid) + (qgf qgd)] (4.40)
S(L xi)
1
3 = (glu x3) (4.41)
TI
1
4 I- (ggU2 X4) (4.42)

where L is the height of the column and S is the cross-sectional area of the column.

The volumetric flow rates of the liquid feed and liquid discharge streams are repre-

sented by the variables qlf and qid, respectively; and the gas feed and gas discharge

rates by the variables qgf and qgd, respectively.

The liquid and air flow rates are modeled using a steady-state relationships

appropriate for the square-root and equal-percentage characteristics of the valves

used in the process control model:

qif = 3Ki /plf x2- piG(xi + AL') (4.43)

qid = xidKldV2 + piG(xi + AL") (4.44)

qgf = x4Kgf(Pgf 2)0.65 (4.45)

qgd = KgdVx (4.46)

where the constants Kij represent the valve coefficients, G is the acceleration of

gravity, pi is density of the liquid, pgf is the inlet pressure of the air, plf is the inlet

pressure of the liquid, AL' is the height difference between valve 1 and the bottom of

the column, and AL" is the height difference between valve 3 and the bottom of the

column. Specific values and definitions of all the multivariable process control model

parameters are given in Table 4.6. The states and inputs must satisfy the constraints
0 0 aa

0 < X < X < [ (4.47)
0 X3 1 0 u2 u2
S0 X4 1

where u1 = u2 = 10[V] (the maximum signal voltage) and x1 and 22 are bounds

established by the user. The bounds of the output vector are identical to those of

the states x1 and x2 because yl = xz and y2 = x2.

Table 4.6: The Multivariable Process Control Model Parameters
Parameter Value Parameter Value
S 0.0043 m2 Kid 4.0 x 10-4 m3s-1 bar
L 1.05 m Kif 6.9 x 10-4 m3s-1 bar
AL' 0.16 m Kgd 9.8 x 10-4 m3.s-1 .bar
AL" 0.40 m Kgf 3.1 x 10-4 m.s-1 bar
pif 0.5 bar T1 1.5 s
pi 1.0 x 103 kg.m-3 g1 0.1 V-1
pgf 0.4 bar 7g 0.3 s
g9 0.1 V-1

The equations of motion are discretized using first order Euler integration.

Let Xk = [xlk X2k k x4k]'. Equations (4.39) through (4.42) take the following form:

xk+ = f(xk) + [gI(xk) 92(xk) [U (4.48)

ylk 1 0 00]
Y2 1 0 0 Zk (4.49)
Y2,k 0 1 0 0

where xk E R4, and f, gl, g2 are vector fields on R4. The system E: D --+ S(R) :

uc -4 y of interest, where V C S(R2), is the composition of the equations of motion

with the input/states constraints for which uc is the voltage command signals ul and

u2. Figure 4.27 shows the block diagram of the system E.

Input State
Constraints Constraints
xcL xC Yk+i

Uk k+1 Z-1

Xk+1 =f(xLk) +g(xL k) k
L --------------------------------------------------------.-----------

Figure 4.27: The block diagram of the system E: V -+ S(2R) : uc y for the
Multivariable Process Control Model Example.

4.5.1 Simulation

The signals of interest are the liquid-level height yl and the gauge pressure

y2 of the enclosed air. The bounds of (4.47) have been set to X1 = 0.75 m and

2 = 0.30 bar. It is assumed that all signals of interest start at their midpoints (i.e.,

yin = 0.375 m and y2n = 0.15 bar for all n < no). From [11], the time constant of

the liquid-level yl (under zero gauge-pressure of air) varies from 50-70 s and the time

constant of the air pressure y2 is 1.3 s. The bandwidth of yl was set at 0.3 rad/sec

to stress an approximate right inverse system j while the bandwidth of y2 was set

at 1/1.3 rad/sec.

The set Ds C S(R2) is defined with initial conditions yli = 0.375 and y2, =

0.15 for all i < 0, and with jlyi 0.375[ < 0.357 and \y2, 0.151 < 0.1428, for all

i > 0. The domain V~ C S('2) of an approximate right inverse system q takes the

sequences of Ds's components 1 and 2 that are passed through a second order discrete-

time low pass filter with cutoff frequencies 0.3 and 1/1.3 rad/sec, respectively. The

outputs of the filters 1 and 2 are then clipped by [0.0, 0.75] and [0.0, 0.30], respectively.

Input '
Constraints

yc with design integer A = 2 for the Multivariable Process Control
Yck U.
Xk-1 Optimization k
State
Constraints
---- L--k----I----^ -- -- _
z1 "- Xk k = f(xL k-)+ g(xL k.) U. XLk
I_______--__----- I I--k-1x-
-------------------------------------------------------------------------------------

Figure 4.28: The block diagram of an approximate right inverse system T: DP -+
S(2) : yc -+ uinv with design integer A = 2 for the Multivariable Process Control
Model Example.

Input
Constraints
YCk Uinvk-1
1 2 -u
k-2 diag Optimization r i
Sw d i State
Constraints

z1 Xk-l = f(xL k-2)+ g(xL k-2) U Ik-I xk-1
I--------------------------------------------------.z-------------

Figure 4.29: The block diagram of an approximate right inverse system I: D> -D+
S(R2) : yc uinv with design integer A = 3 for the Multivariable Process Control
Model Example.

The block diagrams of approximate right inverse systems T.: D -+ S(R')

: ye ui,,n with design integers A = 2 and A = 3 are shown in Figure 4.28 and

Figure 4.29, respectively. Two integration of ul are required to reach yl and sim-

ilarly for u2 and Y2. Two integration of ul are also required to reach y2 but three

integration of u2 are required to reach yl.

For the design integer A = 2 case, the system 9 is causal. The ux optimization

block in Figure 4.28 optimizes ulk to give a minimum value of lyc, D-2 E (yc k) each

step. For the design integer A = 3 case, the system I is strictly causal. The ul and u2

optimization block in Figure 4.29 optimizes Ulk and us2 to give a minimum value of

IYck -D-3E1(yck) each step. The procedure involves undertaking a U2, optimization
if the ul, optimization fails so as to reduce the error IYc, D-3 1E(yck)I. The signal

u2k is optimized to give a minimum value of lyCk+e D-3YE(yc,+i ) provided the
result does not augment the error IYCk D-3~ A (yCk,).

For the simulation of 60 seconds, a sampling period of 150 milli-seconds yields

a two-dimensional sequence of 400 points. The pre-filtered binary random sequences

{ylck} only have values of 0.018, 0.0732 m for k > 1, and the pre-filtered binary
random sequences {Y2Ck only have values of 0.0072, 0.2928 bar for k > 1 because we

are only interested in extreme points of Ds. Table 4.7 shows the Monte Carlo results

for 2 different design integers A of the system T: A = 2 and A = 3 where the error

lyc D-x'^(yc)j is evaluated for each of the sequences. Also included in Table 4.7
is the average number of floating point operations (flops) per Monte Carlo sequence

run for the ul optimization block and for the ul and u2 optimization block. On

the average the ul and u2 optimization subroutine requires 3x greater the number

of flops as the ul optimization subroutine, but a single execution of the ul and u2

subroutine can require 4+ x greater the number of flops than the ul subroutine.

4.5.2 The Estimate of the Performance Bound P((E)

In this section, we will find an estimate of the performance bound P(E) for

the multivariable process control problem of a regulator that regulates the liquid level

Table 4.7: Family of Sequences, DT, Comparison
Error yce D-'ET(yc)l Average Number of Flops
Sequence A = 2 A=3 A =2 A = 3
1 0.02481 0.02190 137.0 367.8
2 0.03742 0.03138 140.3 449.3
3 0.01584 0.01506 142.5 396.0
4 0.05440 0.04699 150.0 432.7
5 0.02408 0.02255 135.5 447.5
6 0.02064 0.01732 137.3 441.0
7 0.02147 0.01917 136.1 379.4
8 0.02613 0.02389 140.9 426.5
9 0.01185 0.01124 135.5 401.5
10 0.00601 0.00589 139.7 374.2
11 0.02080 0.02089 144.9 459.0
12 0.03080 0.02840 137.0 439.1
13 0.01270 0.01113 139.1 404.9
14 0.02746 0.02446 135.5 456.5
15 0.03962 0.03634 144.0 513.7
16 0.02189 0.01942 140.3 413.9
17 0.01852 0.01713 134.6 454.2
18 0.01745 0.01819 136.1 467.9
19 0.01757 0.01527 138.5 441.9
20 0.01637 0.01393 143.1 399.8
21 0.03920 0.03406 144.6 435.2
22 0.03529 0.02838 138.2 418.2
23 0.01935 0.01764 140.9 464.1
24 0.01943 0.01733 138.5 420.2
25 0.01594 0.01466 139.4 393.5
max 0.0544 0.04699 150.0 513.7
min 134.6 367.8
mean 139.6 427.9
st dev 3.703 33.90

67

0.8 I I ,

yl
0.7- / y1C

'0.6 r /
E
E
S0.5 -

/ : \

Q \
0.4 -
0

0.3 -
.0
2 0.2 -

0.1 -

0 1
ol.- ---- i ----- i---- ^~-i--^- ------

0 10 20 30 40 50 60
Time (seconds)

Figure 4.30: For A = 2, the liquid height response yl = D-2'IE(ylc) and the com-
mand Ylc E Dp are shown for the square wave optimization (N1 = 104, N12 = 103).

0.8

7 y1
0.7- \ y -C

S0.6- /
E

00.5- I

U I
S0.4 .
o

S0.3 -

o 0.2 -

0 .1 ..... ... .

0 I \t I
0 10 20 30 40 50 60
Time (seconds)

Figure 4.31: For A = 3, the liquid height response yi = D-3'EI(ylc) and the com-
mand ylc E D D are shown for the square wave optimization (N1i = 78, N2 = 73).

10 20 30 40 50
Time (seconds)

Figure 4.32: For A = 2, the air pressure response y2 = D-2'E(y2c) and the command

y2C E DpS are shown for the square wave optimization (Ni, = 104, NI1 = 103).

0 10 20 30 40 50 60
Time (seconds)

Figure 4.33: For A = 3, the air pressure response Y2 = D-3'T(y2c) and the command

y2e E 'DT are shown for the square wave optimization (N, = 78, 12 = 73).

,0.25

S0.2
o

( 0.15
a
0

! 0.1

Co
0

0.5 F

C
0.4
(c
16

20.3
W
LU
_3
<0
.0 0.2

0.1 F

U
0 10 20 30 40 50
Time (seconds) [dt=0.15 sec]

Figure 4.34: For A = 2, the absolute errors, IYlCk D-2EJ((ylCk)I and lY2Ck -

D-2E(y2ck) are shown for the square wave optimization (N, = 104, N12 = 103)

of yc e DEZ where P(S) = 0.5983 and Iy2c D-2ET(y2c) = 0.0059.

0 10 20 30 40 50
Time (seconds) [dt=0.15 sec]

Figure 4.35: For A = 3, the absolute errors, \ylck D-3'E(ylck)I and |Y2ck -

D-3 J(y2ck)I, are shown for the square wave optimization (N1, = 78, N1 = 73) of

yc E DR where !P(E) = 0.4544 and \y2c D-3EX(y2c)j = 0.2495.

/

I N

rN
I N

rN
I N

IN
I N
I N
I N

\I

> 6

E 5
E
o
0
14

3

U
0 10 20 30 40 50 1
Time (seconds)

Figure 4.36: For A = 2, the valve commands ul and u2 are shown for the square wave
optimization (N1i = 104, N12 = 103) of yc E DP,.

0 10 20 30 40 50
Time (seconds)

Figure 4.37: For A = 3, the valve commands ul and uz are
optimization (N11 = 78, N1, = 73) of yc Dv.

shown for the square wave

-ul
- u2

-- I

4

- - - -

in a pressurized tank by using the optimization definition as follows:

P(E)= sup yc D-AE (yc)l (4.50)
ycEV",

where the system E and an approximate right inverse system T1 with design integers

A: 2 and 3 are illustrated in Figure 4.27, Figure 4.28, and Figure 4.29, respectively.

In Section 4.4.2, a nonlinear programming algorithm was developed to calcu-

late the estimate of the performance bound P(E). It exploits the traits of the system.

Due to the sluggishness of the liquid height response, only one transition was needed

from high to low in both components of the signal yc. The transition is at N1i steps

for ylc and N1, steps for Y2c. The results of the nonlinear programming algorithm

are shown in Table 4.8 and Figures 4.30-4.37.

Table 4.8: The Estimate of the Performance Bound P(E) for the Multivariable Pro-
cess Control Model Example
Parameters lye D-'EZ (yc)|
Design Integer (N,, NA12) for yc E Dp>

A = 2 (104, 103) P(E) = 0.5983

A = 3 (78,73) P(E) = 0.4544

Figure 4.30 and Figure 4.31 represent the optimal liquid height response and

command for A = 2 and A = 3, respectively. The ul and u2 subroutine improves

the liquid height response for A = 3, at the expense of the air pressure response.

Figure 4.32 and Figure 4.33 represent the optimal air pressure response and command

for A = 2 and A = 3, respectively. The ul and u2 subroutine increases the air pressure

to help drain the tank while the ul subroutine is operating with near zero gauge

pressure on yl which has a time constant between 50-70 s.

The optimal absolute error sequences are plotted in Figure 4.34 and Fig-

ure 4.35 for A = 2 and A = 3, respectively. In the ul & u2 subroutine, the signal u2,

72

is optimized to give a minimum value of IYck+1 D-3ET(yck+ )I provided the error

Iyck -D3-'ITy(yck)I does not augment. In the process, the error IY2Ck -D-3EJ!(y2ck)I
is augmented but never above the value of the error IylCk D-3 E(ylC,)I which is

evident in Figure 4.35. The optimal valve commands ul and u2 are plotted in Fig-

ure 4.36 and Figure 4.37 for A = 2 and A = 3, respectively. For A = 2, the signal ul

is constrained throughout most of the simulation while for A = 3 the signal u2, which

influences yl in three steps, is optimally used to reduce the error ylcc,-D-3' 1(ylCk)j.

CHAPTER 5
CONCLUSION

5.1 Summary

This dissertation presents new results for the problem of disturbance atten-

uation for nonlinear control systems. Of primary significance is the derivation of a

performance bound that provides a measure of the ability to match a set by subset

of the image of a system. Our bound arose from the desire of providing an estimate

of the minimal effect of the disturbance signal d on the output signal z of Figure 1.1.

It turns out that the calculation of the performance bound is rather labori-

ous because its definition (Definition 3.4.1) uses Im E. In practice, we will find an

estimate of the performance bound that uses an approximate right inverse system

in its calculation instead of Im E. The estimate of the performance bound involves

an optimization problem which relates to finding a global maximum of a non-convex

function. Nonlinear programming is used to calculate the estimate of the performance

bound.

Three systems of practical origin were selected for the analysis of the per-

formance bound. In the PM stepper motor example, a closed form estimate of the

performance bound was found which was verified by simulation. In the aerodynamic

model example, a nonlinear programming algorithm was determined to compute the

estimate of the performance bound. The multivariable process model example pro-

vided us insight into how different design integers can influence the estimate of the

performance bound.

5.2 Future Directions

The results of this dissertation serve to enhance many of the tools used to

obtain the optimal controller for nonlinear control systems. In the process of com-

pleting one particular objective, other promising areas frequently appear as directions

for additional research. One of the more intriguing of such areas emerging from this

work is the derivation of a bound that handles constant bias disturbances.

Referring to Figure 1.1, and using the parameterization (1.3) with the assump-

tion that the system E is stable (i.e., P = E and Q = I), we get

Ec(v, d) = d + EZ(v, d). (5.1)

The stable and causal system q is selected such that

0(v, d) = T (Ev d) (5.2)

where I is an approximate right inverse system. All that is needed to define a

performance bound is an optimization set. In this case, the optimization set is the

image of [Ev d] where the domain of [Ev d] is a bounded set. The performance

bound is given by

P(E) = ~ (Im [Ev d], D-AIm) (5.3)

The bound of (5.3) is a theoretical estimate of the minimal effect of the disturbance

d on the output z.

Let's use A = 0 initially in the bound optimization problem. Note that there

is another A associated with the approximate inverse optimization problem. Both are

equal for most applications. In this application the design integer A will be associated

with an approximate inverse and A > 0. Now we need to modify our definition of the

norm p, for the bound. For a sequence u E S(Rm) and for a j > 0, the norm p,j is

defined

p,j (u) d sup (1 + e)-i i, C > 0. (5.4)
i>j

The class of disturbances is a constant bias; i.e., the magnitude of the disturbance is
unknown, but it is constant (i.e., a step disturbance). Let us define our performance
bound for A = 0, denoted by R(E), such that

R(E) d_ inf sup p,,x (d + EC(v, d) Ev) (5.5)

where Ev is our desired response. The supremum is over the input signal space and
the infimum is over the stable and causal systems q. The stable and causal system
Scan be chosen as

q(v,d) = I (D-Ev d) (5.6)

where the design integer A > 0 of an approximate right inverse system I is the least
latency of the system E.
It will be shown that the bound R(E) is well suited for a constant bias dis-
turbance d. If the system I were a right inverse of the system E, then EXI = D'I.
Let us substitute (5.6) into (5.5), we get

R(E) = sup pe,A (d + E (D- v- d) Ev)

= sup p,,A (d + D (I (D- Ev d) Ev)

= sup pe, (d + (Ev D d) Ev)

= sup p, (d D"d)

R(E) = 0

where the supremum is taken over the input signal space.
The bound R(E) can be calculated using techniques of Chapter 4. The op-
timization set is the image of [D-AEv d] where the domain of [D-^Ev d] is a
bounded set. For the PM stepper motor example, the estimate of the performance
bound P(E) of (4.26) turns out to be the performance bound R(E) because the
performance bound for this problem is independent of the optimization set. Thus

( = m At2 (5.7)
27rG.L- J

76

for the PM stepper motor example.

Throughout this dissertation, no new results were found for the design of the

equivalent controller C of Figure 1.1 because the results from [22] were sufficient. In

this emerging area of of deriving a performance bound 7(E) to handle a constant

bias, attention should be given to the controllers that implement the new control laws.

There will be other performance requirements in addition to the bound requirement,

and so the performance index will be augmented to include additional variables.

REFERENCES

[1] M. A. Armstrong, Basic Topology. New York: Springer-Verlag, 1983.

[2] K. J. Astrom and B. J. Wittenmark, Adaptive Control, Second Edition. New
York: Addison-Wesley Publishing Company, 1995.

[3] M. S. Bazaraa, H. D. Sherali and C. M. Shetty, Nonlinear Programming: Theory
and Algorithms. New York: John Wiley & Sons, Inc., 1993.
[4] J. H. Blakelock, Automatic Control of Aircraft and Missiles, Second Edition.
New York: John Wiley & Sons, Inc., 1991.
[5] M. Bodson, J. N. Chiasson, R. T. Novotnak and R. B. Rekowski, "High-
performance nonlinear feedback control of a permanent magnet stepper motor",
IEEE Transactions on Control Systems Technology, Vol. 1, No. 1, pp. 5-14, 1993.
[6] E. W. Cheney, Introduction to Approximation Theory, Second Edition. New
York: Chelsea Publishing Company, 1982.
[7] R. J. P. de Figueiredo and G. Chen, Nonlinear Feedback Control Systems: an
Operator Theory Approach. New York: Academic Press, Inc., 1993.
[8] C. A. Desoer and W. S. Chan, "The feedback interconnection of lumped linear
time-invariant system", Journal of The Franklin Institute, Vol. 300, No. 5 & 6,
pp. 335-351, 1975.
[9] D. S. Dummit and R. M. Foote, Abstract Algebra. Englewood Cliffs, NJ: Prentice
Hall, 1991.
[10] M. A. Eggiman, O. D. Crisalle and R. Longchamp, "A linear-programming
predictive controller with variable horizon" in Proc. 1992 Amer. Contr. Conf.,
Chicago, IL, pp. 1568-1575, 1992.
[11] M. A. Eggiman, O. D. Crisalle and R. Longchamp, "Design of a multivari-
able linear-programming predictive controller with variable receding horizon"
submitted for publication in Journal of Process Control.
[12] B. Etkin, Dynamics of Flight: Stability and Control. New York: John Wiley &
Sons, Inc., 1982.
[13] M. W. Hirsch and S. Smale, Differential Equations, Dynamical Systems and
Linear Alegebra. New York: Academic Press, Inc., 1974.
[14] J. Hammer, "Nonlinear systems: stability and rationality", International Jour-
nal of Control, Vol. 40, No. 1, pp. 1-35, 1984.

78

[15] J. Hammer, "On nonlinear systems, additive feedback, and rationality", Inter-
national Journal of Control, Vol. 40, No. 5, pp. 953-969, 1984.
[16] J. Hammer, "Nonlinear systems, stabilization and coprimeness", International
Journal of Control, Vol. 42, No. 1, pp. 1-20, 1985.
[17] J. Hammer, "Stabilization of nonlinear systems", International Journal of Con-
trol, Vol. 44, No. 5, pp. 1349-1381, 1986.
[18] J. Hammer, "Fraction representation of nonlinear systems: a simplified ap-
proach", International Journal of Control, Vol. 46, No. 2, pp. 455-472, 1987.
[19] J. Hammer, "Assignment of dynamics for nonlinear recursive feedback systems",
International Journal of Control, Vol. 48, No. 3, pp. 1183-1212, 1988.
[20] J. Hammer, "Robust stabilization of nonlinear systems", International Journal
of Control, Vol. 49, No. 2, pp. 629-653, 1989.
[21] J. Hammer, "Fraction representations of nonlinear systems and non-additive
state feedback", International Journal of Control, Vol. 50, No. 5, pp. 1981-
1990, 1989.
[22] J. Hammer, "Internally stable nonlinear systems with disturbance: A parameter-
ization", IEEE Transactions on Automatic Control, Vol. 39, No. 2, pp. 300-314,
1994.
[23] A. Isidori, Nonlinear Control Systems, Third Edition. New York: Springer-
Verlag, 1995.
[24] T. Kenjo, Stepping motors and their microprocessor control. New York: Oxford
University Press, 1984.
[25] H. K. Khalil, Nonlinear Systems, Second Edition. New York: Macmillan Pub-
lishing Company, 1996.
[26] A. N. Kolmogorov and S. V. Fomin, Introductory Real Analysis. New York:
Dover Publications, Inc., 1975.
[27] K. Kuratowski, Introduction to Set Theory and Topology. New York: Pergamon
Press, 1972.
[28] S. Mac Lane and G. Birkhoff, Algebra, Third Edition. New York: Chelsea
Publishing Company, 1993.
[29] 0. L. Mangasarian, Nonlinear Programming. Philadelphia: Society for Industrial
and Applied Mathematics, 1994.
[30] D. McLean, Automatic Flight Control Systems. Englewood Cliffs, NJ: Prentice
Hall, 1990.
[31] D. McRuer, I. Ashkenas and D. Graham, Aircraft Dynamics and Automatic
Control. Princeton, NJ: Princeton University Press, 1973.
[32] J. R. Munkres, Topology: A First Course. Englewood Cliffs, NJ: Prentice Hall,
1975.

79

[33] R. C. Nelson, Flight Stability and Automatic Control. New York: McGraw-Hill,
Inc., 1989.
[34] H. Nijmeijer and A. J. van der Schaft, Nonlinear Dynamical Control Systems.
New York: Springer-Verlag, 1990.
[35] R. J. Norris, Analysis of multivariable control systems in the presence of struc-
tured uncertainties. Ph.D. thesis, University of Florida, Gainesville, 1990.
[36] A. D. B. Paice and J. B. Moore "Robust stabilization of nonlinear plants via left
coprime factorizations", Systems & Control Letters, Vol. 15, No. 2, pp. 125-135,
1990.
[37] T. J. Rivlin, An Introduction to the Approximation of Functions. New York:
Dover Publications, Inc., 1969.
[38] R. T. Rockafellar, Convex Analysis. Princeton, NJ: Princeton University Press,
1970.
[39] H. H. Rosenbrock, State-space and Multivariable Theory. New York: John Wiley
& Sons, Inc., 1970.
[40] H. L. Royden, Real Analysis. New York: The Macmillan Company, 1963.
[41] W. Rudin, Functional Analysis. New York: McGraw-Hill, Inc., 1973.
[42] W. Rudin, Principles of Mathematical Analysis. New York: McGraw-Hill, Inc.,
1976.
[43] G. Tao and P. V. Kokotovid, "Adaptive control of plants with unknown dead-
zones", IEEE Transactions on Automatic Control, Vol. 39, No. 1, pp. 59-68,
1994.
[44] G. Tao and P. V. Kokotovic, Adaptive Control of Systems with Acutator and
Sensor Nonlinearity. New York: John Wiley & Sons, Inc., 1996.
[45] J. van Tiel, Convex Analysis: An Introductory Text. New York: John Wiley &
Sons, Inc., 1984.
[46] M. S. Verma and L. R. Hunt, "Right coprime factorizations and stabilization
for nonlinear systems", IEEE Transactions on Automatic Control, Vol. 38, No.
2, pp. 222-231, 1993.
[47] M. Vidyasagar, Nonlinear Systems Analysis, Second Edition. Englewood Cliffs,
NJ: Prentice Hall, 1993.
[48] G. Zames, "On the input-output stability of time-varying nonlinear feedback
systems, Part I", IEEE Transactions on Automatic Control, Vol. AC-11, No. 2,
pp. 228-238, 1966.
[49] G. Zames, "On the input-output stability of time-varying nonlinear feedback
systems, Part II", IEEE Transactions on Automatic Control, Vol. AC-11, No.
3, pp. 465-476, 1966.
[50] G. Zames, "Feedback and optimal sensitivity: Model reference transformations,
multiplicative seminorms, and approximate inverses", IEEE Transactions on
Automatic Control, Vol. AC-26, No. 2, pp. 301-320, 1981.

BIOGRAPHICAL SKETCH

Rafael J. Fanjul Jr. was born in West Palm Beach, Florida, on May 28, 1963.

He received his bachelor's degree in mechanical engineering from Georgia Institute

of Technology in 1985. He then entered graduate school at Georgia Institute of

Technology and received a master's degree in electrical engineering in 1986. From

1986 to 1992 he was employed by Martin Marietta, Orlando, Florida, as a guidance

and control engineer. During his tenure at Martin Marietta, he attended night school

at the University of Central Florida and received a master's degree in mathematical

science in 1992. He has been a Professional Engineer in the State of Florida since

September 1990. In August 1992, he was granted an educational leave of absence from

Martin Marietta to pursue a doctoral degree. Since August 1992, he has been a Ph.D.

student and pre-doctoral fellow at the University of Florida under the supervision of

Professor J. Hammer.

I certify that I have read this stupid and that iIn my opinion it conforms to
acceptable standards of scholarly presentation and is fully adequate. in scope ',and
quality, as a dissertation for the degree of Doctor of Philosophy.

,^-- t$\A
Jacolb Haainer. Chairman
Professor of Electrical and Com'piiuter
Engineering

I certify that 1 have read this study and that in my opinion it conforms to
acceptable standards of scholarly presentation and is fully adequate. in scope and
quality, as a dissertation for the degree of Doctor of Philosophy.

1John M. MI. Anderson
Assistant Professor of Electrical and
computer r EngineerineI

I certify that I have read this sludy and that in my opinion it conforms to
acceptable standards of scholarly presentation and is fully adequate. in scope and
quality, as a dissert at ion for the degree of Doctor of Philosophy.

Oscar D. Crisalle
Associate Professor of chemical l
Engineering

I certify that I have read this study and that in my opinion it conforms to
acceptable standards of scholarly presentation and is fully adequate. in scope and
quality, as a dissertation for the degree of Doctor of Philosophy.

:al and Computer

Engineering

I certify that I have read this study and that in my opinion it conforms to
acceptable standards of scholarly presentation and is fully adequate, in scope and
quality, as a dissertation for the degree of Doctor of Philosophy.

Carla A. Schwartz
Assistant Professor of Electrical and
Computer Engineering

This dissertation was submitted to the Graduate Faculty of the Department
of Electrical and Computer Engineering in the College of Engineering and to the
Graduate School and was accepted as partial fulfillment of the requirements for the
degree of Doctor of Philosophy.

December 1996 -_______
Winfred M. Phillips
Dean, College of Engineering

Karen A. Holbrook
Dean, Graduate School

LD
178C
1996

UNIVERSITY OF FLORIDA
S111111111111 111 1111 1 l l I I 1111111
3 1262 08554 4434

	Title Page
	Dedication
	Acknowledgement
	Table of Contents
	Abstract
	Introduction
	Terminology and basics
	The performance bound
	Calculation of the performance...
	Conclusion
	References
	Biographical sketch

MILLISECOND	CLASS.METHOD	MESSAGE
0	sobekcm_page_globals.constructor
0	sobekcm_page_globals.constructor	Application State validated or built
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.constructor	Navigation Object created from URI query string
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.display_item	Retrieving item or group information
0	sobekcm_page_globals.get_entire_collection_hierarchy	Retrieving hierarchy information
0	sobekcm_assistant.get_entire_collection_hierarchy
0	cached_data_manager.retrieve_item_aggregation
0	cached_data_manager.retrieve_item_aggregation	Found item aggregation on local cache
0	item_aggregation_builder.get_item_aggregation	Found 'all' item aggregation in cache
0	system.web.ui.page.page_load (ufdc.page_load)
0	sobekcm_page_globals.constructor.on_page_load
0	html_echo_mainwriter.add_style_references	Adding style references to HTML
0	html_echo_mainwriter.add_text_to_page	Reading the text from the file and echoing back to the output stream
33	html_echo_mainwriter.add_text_to_page	Finished reading and writing the file