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