Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Robust control of periodically time-varying systems
(USC Thesis Other)
Robust control of periodically time-varying systems
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
ROBUST CONTROL OF PERIODICALLY TIME-VARYING
SYSTEMS
by
Serkan Kalender
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(MECHANICAL ENGINEERING)
December 2007
Copyright 2007 Serkan Kalender
ii
ACKNOWLEDGEMENTS
I would like to thank my Ph.D. advisor Dr. Henryk Flashner. This
dissertation could not have materialized without his guidance and encouragement. I
would also like to thank Dr. Bingen Yang and Dr. Petros Ioannou for serving in my
dissertation committee for their helpful recommendations. Additionally, I would like
to thank Dr. Paul K. Newton, Dr. Firdaus E. Udwadia for serving in my guidance
committee.
I am grateful to the University of Southern California Aerospace and
Mechanical Engineering Department for granting me scholarships throughout my
stay here.
Finally, I wish to thank my parents and my brother Volkan, for their support
throughout my research.
iii
TABLE OF CONTENTS
Acknowledgements ii
List of Tables v
List of Figures viii
Abstract xi
Chapter 1: Introduction 1
Chapter 2: Characteristics of Periodic Mechanical Control Systems 7
2.1 Hamiltonian Control Systems 7
2.1.1 Canonical Transformations 10
2.1.2 Lagrange and Poisson Brackets 12
2.1.3 Invariants in Hamiltonian Mechanics 14
2.1.3.1 Integral Invariants of Poincarè 14
2.1.4 Symplecticity of Poincarè Maps 15
2.1.5 Periodic Hamiltonian Systems 17
2.2 Properties of Symplectic Matrices 19
2.3 Stability of Equilibrium Points of Nonlinear Periodically Time
Varying Dynamical Systems 19
Chapter 3: Discrete-Time Control of Linear Time-Periodic Systems 23
3.1 Discrete-Time Formulation of Linear Time-Periodic
Control Systems 24
3.1.1 Examples 29
3.1.1.1 Undamped and Unforced Mathieu Equation 29
3.1.1.2 Underdamped and Unforced Mathieu Equation 30
3.2 Control System Formulation 31
3.2.1 Fast Sampling Control System Formulation 32
3.3 Controller Design 34
3.3.1 Full State Feedback Optimal Regulator Design 35
3.3.2 Output Feedback Regulator Design 37
3.3.2.1 Synchronous Sampling Observer and Control
Law Update 38
3.3.2.2 Fast Sampling Observer 39
3.4 Examples 40
3.4.1 Dead-beat Full State Regulator for Mathieu
Equation 40
3.4.2 Dead-beat Output Feedback Regulator for Mathieu
Equation 47
iv
3.4.2.1 Dead-beat Regulator with Synchronous
Sampling Observer 47
3.4.2.2 Dead-beat Regulator with Fast
Sampling Observer 50
3.4.3 Multivariable Feedback Design Example:
Control Design for Periodically Excited Inverted Double
Pendulum 53
3.4.3.1 Full State Feedback Regulator Designs 55
3.4.3.1.1 Dead-beat Controller Design 56
3.4.3.1.2 Optimal Controller Design 56
3.4.3.2 Fast-Sampled Full State Feedback Regulator
Designs 60
3.4.3.1.1 Dead-beat Controller Design 61
3.4.3.1.2 Optimal Controller Design 62
3.4.3.3 Output Feedback Regulator Designs 65
3.4.3.1.1 Synchronous Sampling Observer 65
3.4.3.1.2 Fast Sampling Observer 66
Chapter 4: Robust Stability Analysis 70
4.1 Discrete-Time Representation of Periodic Systems with
Real Parametric Uncertainties 71
4.2 Examples 73
4.2.1 Active Control of Helicopter Rotor Blade Flapping
Motion in Forward Flight 73
4.2.2 Double Inverted Pendulum with Load Direction
Uncertainty 82
4.3 Error Analysis for Truncated Point-Mappings 84
4.3.1 Uncertain Underdamped Mathieu Equation 90
4.4 Robust Stability Analysis by Structured Singular Value
Theory 93
4.4.1 Computation of Structured Singular Value 96
4.4.2 Structured Singular Value Analysis Examples 97
. 4.4.2.1 Lightly Damped Mathieu Equation 97
4.4.2.2 Robust Stability Analysis of Multivariable
Feedback Design Example with Respect to
Load Parameter Uncertainties 101
Chapter 5: Engineering Application Example
Spacecraft Attitude and Momentum Unloading Control
System Design 111
5.1 Spacecraft Momentum Unloading Problem 112
5.1.1 Spacecraft Attitude Dynamics 113
5.1.2 Attitude Kinematics 113
5.1.1 Gravity Gradient Torque 114
v
5.2 Magnetic Momentum Unloading 115
5.2.1 Earth Magnetic Field Models 115
5.2.2 Straight and Tilted Dipole Model Approximations 116
5.3 Linearized Model 118
5.3.1 Attitude Control Loop 119
5.3.2 Momentum Unloading Control Loop 119
5.4 Design Example 120
5.4.1 Momentum Unloading Control System Design Using
Straight Dipole Model 121
5.4.2 Momentum Unloading Control System Design Using
Tilted Dipole Model 124
Chapter 6: Conclusion 132
References 134
Appendix-A: Derivation of Point-Mapping Algorithm 139
Appendix-B: Derivation of Truncated Point-Mapping Algorithm 142
Appendix-C: Helicopter Rotor Dynamics 147
Appendix-D: Structured Singular Value Theory 157
vi
LIST OF TABLES
3.1 Effect of step size on accuracy of the computational method 29
3.2 Pulse, linear, and quadratic function parameterized control
input matrices for fast sampling observer 50
3.3 State transition matrices for fast-sampled state feedback design 60
3.4 Input transition matrices for fast-sampled state feedback design 61
3.5 Dead-beat feedback gains for fast-sampled full-state feedback design 61
3.6 Steady-state optimal feedback gains for fast-sampled full-state feedback
design 62
4.1 Comparison of state and input matrices—direct integration versus linear
approximation 78
4.2 Comparison of state and input matrices-direct integration versus linear
approximation with assuming quadratic terms in advance ratio as an
independent perturbation 79
4.3 Spectral radius for the rotor control system--exact versus linear
approximation 80
4.4 State transition matrices with load direction uncertainty-- linear versus
second order approximation 83
4.5 Error analysis for linear and second order approximation 84
4.6 Computed values of the nominal state transition matrix and state
transition matrices due to the uncertainty 91
4.7 Computed norm bounds on error in state transition matrix due to
truncation operation 91
4.8 Exact norm bounds on error in state transition matrix due to
truncation operation 92
4.9 Stability bounds by µ analysis 100
vii
4.10 Comparison of state transition matrices obtained by first, second, and
third order truncated approximations 105
4.11 Computed norm bounds on state transition matrix error 106
4.12 Computed norm bounds on error for linear and quadratic function
parameterizations 107
4.13 Lower and upper bounds on µ for optimal linear and quadratic function
parametrized controllers 108
4.14 Lower and upper bounds on µ for dead-beat linear and quadratic function
parametrized controllers 108
5.1 Input transition matrices 122
5.2 Nominal and uncertainty matrices due to the uncertainties δ
1
, and δ
2
127
5.3 Feedback gains 129
5.4 Computed µ-upper bounds 129
viii
LIST OF FIGURES
3.1 Comparison of trajectories: point-mapping versus numerical
integration 30
3.2 Block diagram representation of T -periodic sampled data system 31
3.3 Block diagram for the discrete-time state variable system 31
3.4 Discrete-time formulation of the closed loop system 32
3.5 Discrete-time formulation of the closed loop system with fast sampling 33
3.6 Multi-rate sampling 34
3.7 Controller-observer configuration 38
3.8 Full state feedback regulator block diagram 40
3.9 Closed loop response, and control effort for dead-beat regulator with
pulse function parameterized controller for Mathieu equation 42
3.10 Closed loop response, and control effort for dead-beat regulator with
linear function parameterized controller for Mathieu equation 44
3.11 Closed loop response, and control effort for dead-beat regulator with
quadratic function parameterized controller for Mathieu equation 46
3.12 State and estimated state histories for Mathieu equation with estimated
state feedback control (Pulse-function parameterization) 48
3.13 State and estimated state histories for Mathieu equation with estimated
state feedback control (Linear function parameterization) 48
3.14 State and estimated state histories for Mathieu equation with estimated
state feedback control (Quadratic function parameterization) 49
3.15 State and estimated state histories for Mathieu equation with estimated
state feedback using fast sampling observer
(Pulse function parameterization) 51
ix
3.16 State and estimated state histories for Mathieu equation with estimated
state feedback using fast sampling observer
(Linear function parameterization) 52
3.17 State and estimated state histories for Mathieu equation with estimated
state feedback using fast sampling observer
(Quadratic function parameterization) 52
3.18 Inverted double pendulum subjected to periodic Forcing 53
3.19 Closed loop state trajectories and control forces for the dead-beat linear
function parametrized controller 58
3.20 Closed loop state trajectories and control forces for the dead-beat
quadratic, steady state optimal linear, and quadratic function
parametrized controllers 59
3.21 Closed loop state trajectories for fast-sampling controllers 63
3.22 Control forces for fast-sampling configuration 64
3.23 Closed loop state trajectories for output feedback controllers 68
3.24 Output feedback control forces 69
4.1 Full state feedback pulse function parameterized controller response
for rigid rotor flapping model 81
4.2 General control configuration 94
4.3 M-∆ Connection 98
4.4 Closed loop state trajectories with the quadratic function parametrized
dead-beat controller for different parameter values 109
4.5 Closed loop state trajectories with the quadratic function parametrized
steady-state optimal controller for different parameter values 110
5.1 Spacecraft attitude and momentum unloading control system 114
5.2 Simulated time history of the attitude angles for the combined attitude
and momentum control system design (using straight dipole model) 123
x
5.3 Simulated time history of the momentum wheel angular momentum
for the combined attitude and momentum control system design
(using straight dipole model) 123
5.4 M-∆ Connection 128
5.5 Simulated time history of the attitude angles for the combined attitude
and momentum control system design (using tilted dipole model) 130
5.6 Simulated time history of the momentum wheel angular momentum
for the combined attitude and momentum control system design
(using tilted dipole model) 131
C.1 Single main rotor/ Tail rotor configuration 147
C.2 Rotor in forward flight 149
C.3 Rotor blade motion 150
C.4 Helicopter rotor control with swash plate mechanism 150
C.5 Aerodynamic forces on the blade section 151
C.6 Velocity components 152
C.7 Forces acting on the blade section 155
D.1 Canonical robust control configuration 158
D.2 N∆-structure 159
D.3 M∆-structure 160
xi
ABSTRACT
In this dissertation, a practical control design method for linear periodically
time varying systems is presented. The method is based on parametrization of
control, and obtaining an equivalent discrete-time control system via a point-
mapping computational algorithm. The design method simplifies the design
procedure, and guarantees asymptotic stability of the closed loop system. In addition,
the designed controllers can directly be implemented digitally. Moreover, it is shown
that the method can easily be applied to various control configurations, as well as
multi-rate configurations. The design method is then extended for analysis of
robustness of the control design with respect to plant parametric uncertainties. This
is achieved by computation of approximate discrete time dynamics of the perturbed
system by truncated point-mappings. By computing an upper norm bound on the
error due truncated approximations, the robustness analysis of the system due to
parametric uncertainties is formulated as a discrete time structured singular value
problem. Simulation studies with mathematical models as well as engineering
application models, show good performance of the control system design method
and demonstrate the ability of the extension of the method to assess the robustness of
the closed-loop system with respect to parametric uncertainties.
1
Chapter 1
Introduction
Nonlinear periodically time varying systems appear in various engineering
applications such as structures under periodic loading (see Pandiyan and Sinha [47],
Bolotin [8]), helicopter rotor control applications (see Arcara and Bittanti [2], Calico
and Wiesel [11,59], Mc Killip [41]), high-speed rotating machinery (see Tondl,
Ruijgrok, Verhulst and Nabergoj [57]) and robotics. Models of this class of dynamic
systems consist of a set of coupled usually nonlinear ordinary differential equations
with periodically time varying coefficients (NLPTV). The investigation of stability
of a particular periodic solution of a nonlinear autonomous system also results
NLPTV equations. Because of their importance to many areas of engineering, this
type of problems has attracted a great amount of theoretical research (see Hale [27],
Urabe [58]).
The two main categories for analysis of this class of problems are the
perturbation analysis and numerical simulation (see Bogoliobov and Mitropolsky [7],
Nayfeh and Mook [43], Thompson and Stewart [56]). Though an analytical treatment
of the problem is possible with the perturbation methods, they are valid only locally.
In addition perturbation methods require a great deal of challenge for higher order
2
systems and systems with large number of parameters. Though being well suited for
investigating the system dynamics in state space, stability analysis and global
domain of attraction for the fixed points of the system cannot be addressed
theoretically by numerical simulation studies. Moreover, for parametric studies such
as bifurcation analysis requires a great amount of computational time.
Because of the difficulties mentioned above with the numerical simulation
studies, combined numerical and analytical methods have been developed. A way to
analyze nonlinear periodically time varying systems is based on the idea of Poincare
Map, is to formulate the equivalent discrete time dynamics by defining a point
mapping (see Bernussou [6]) .The point mapping techniques have been widely used
for understanding the global dynamics of the NLPTV systems and date back to
Poincare (see Guckenheimer and Holmes [25]). An algorithm for approximating the
mapping in analytical form has been introduced by Flashner [19], and Flashner and
Hsu [21]. For a modification of the algorithm see Flashner and Guttalu [20].
An alternate quantitative approach uses Floquet-Lyapunov theorem (see
Yakubovich and Starzhinskii [63]) that transform the quasilinear periodic systems
into a new set of equation with time invariant linear parts. The stability and
bifurcation characteristic of the original systems is preserved by this transformation.
To obtain the transformation one has to compute the state transition matrix of the
linear part as function of time. Sinha and coworkers developed a computational
procedure by which the L-F transformation is obtained in terms of Chebyshev
3
polynomials (see Sinha and Wu [53]), and they also a proposed an algorithm to
obtain the state transition matrix symbolically (see Sinha and Butcher [54]).
For control design for periodic systems, the full nonlinear differential
equations with periodic coefficients governing the dynamics of the system should be
analyzed for the stability analysis and controller synthesis. However, in many
engineering applications the local behavior about a known solution and stabilization
of this solution of interest is sufficient for associated stability analysis and control
system design.
Upon linearization, the above class of NLPTV systems yield a system
described by a set of linear ordinary differential equations with periodically time
varying coefficients (LPTV). In control system design formulation, the above
approach of linearization results in a LPTV with uncertain parameters. The stability
analysis for such systems is usually performed by the Floquet-Lyapunov theory with
the assumption that the parametric uncertainties are negligible. Another approach for
stability analysis of uncertain LPTV systems is based on the concept of “Integral
Quadratic Constraint (IQC),” first proposed by Yakubovich, and further developed
by Megretski and his colleagues [33, 34, 40]. The other approach based on the
concept of “lifting” transforms the system into an equivalent time invariant discrete-
time system, and is mainly used for stability analysis and control design for LPTV
discrete time systems (see Bamieh and Pearson [3], Chen and Francis [14], Dellerud
[16]).
4
A conservative approach for robust stability analysis for LPTV systems
neglects the periodic nature of the underlying model --i.e. assuming a nominal LTI
plant and representing the periodic coefficients as time varying uncertainties in
addition to the dynamic and/or parametric in the original model --, and uses operator
theoretic methods such as small gain and passivity theorems (see Harris and Valenca
[28], Khalil [36], Narendra and Taylor [42], Sundareshan and Thathachar [55]).
One approach to reduce the conservativeness introduced by allowing a larger
class of perturbations in the above is to introduce the bound on the derivative of the
periodic terms and to use passivity theory (see Willems [62], Zames and Kallman
[66], Fu, Souza and Dasgupta [24], Megretski [40]).
The other main approaches for decreasing the level of conservativeness can
be grouped as characteristic polynomial and linear fractional transformation
frameworks. The first approach applies to the cases where the uncertainties in the
model are parametric. This approach leads to Kharitinov-type results such as
“Kharitinov Theorem”--restricted to the cases where each coefficient of the
characteristic polynomial is only a function of one parameter--, “Edge Theorem” --
restricted to the case where each coefficient of the characteristic polynomial is
dependent upon the parameters affinely, and “Mapping Theorem” -- the parametric
uncertainties enter to the characteristic polynomial in a multi-linear fashion--, etc.
(see Barmish for a complete discussion on Kharitinov -type theorems [4], and
Mansour, Balemi and Truöl [39]). Similar results for systems with parametric
5
uncertainty in state space form can be found in the literature using Lyapunov type
arguments in which the robust stability problem is to construct a common quadratic
Lyapunov function for all possible parametric uncertainties. This condition reduces
to check whether the vertice matrices satisfy Lyapunov inequalities when the entries
of the state matrix depend linearly on the parametric uncertainties (see Horisberger
and Belanger [30], Boyd and Yang [9]).
The second approach applies both parametric and dynamic uncertainties, and
is based on the linear fractional transformation (LFT) framework. This approach is
termed as the structured singular value theory (see Doyle [17] for μ -analysis, see
Safonov [50] multivariable stability margin k
m
–analysis).
In this dissertation, a new methodology for designing control systems for
linear periodically time varying systems is presented. The proposed methodology
formulates the control design problem for linear periodic systems via control vector
parameterizations, in discrete-time domain by applying the point-mapping technique.
The approach yields an equivalent discrete-time linear time invariant representation
of the original problem that allows the known LTI control design techniques be
applied. The formulation allows stabilization of trivial equilibrium as well as higher
periodic solutions without a need for finding the steady state solution of the plant.
Another advantage of the design methodology used is that the designed controller
can be implemented directly, thus rendering the need for analysis effort required in
discretization of the control law.
6
The method is then extended for analysis of robustness of the control design
with respect to plant parametric uncertainties. This is achieved by computation of
approximate discrete time dynamics of the perturbed system by truncated point-
mappings. By computing an upper norm bound on the error due truncated
approximations, the robustness analysis of the system due to parametric uncertainties
is then formulated as a discrete time structured singular value problem.
This dissertation is organized as follows. In chapter 2, characteristics of
Hamiltonian control systems and specifically periodic Hamiltonian systems are
provided. In chapter 3, discrete-time representation of the system dynamics, and the
formulation of a sampled-data control problem are introduced. In chapter 4, the
extension the methodology for analysis of robustness with respect to parametric
uncertainties is presented. Application of the design methodology and its extension
for robust stability analysis to spacecraft attitude and momentum unloading control
design problem is then presented in chapter 5. Finally concluding remarks are given
in the final chapter.
7
Chapter 2
Characteristics of Periodic Mechanical Control Systems
A large class of aerospace and mechanical systems can be formulated using
the Euler-Lagrangian, or equivalently Hamiltonian equations of motion and they will
be referred to as Hamiltonian control systems. In this chapter, some of the qualitative
features of the above class of control systems are presented, (see Goldstein[26]), and
Nijmeijer and van der Schaft [45]), the definition of symplectic matrices and their
properties are given, and finally stability results for general nonlinear periodically
time varying systems are presented.
2.1 Hamiltonian Control Systems
A mechanical control system with n-degrees of freedom is described by n
Lagrange equations of motion:
e
i
i i
d L L
F
dt q q
∂ ∂
− =
∂ ∂
ɺ
(2.1)
where
1
( ,..., )
n
q q q = , and
1
( ,..., )
n
q q q = ɺ ɺ ɺ are the generalized coordinates and
generalized velocities respectively. L is the Lagrangian and is defined as the
difference between the kinetic and potential energy of the system
(i.e. ( , ) ( ) L T q q V q = − ɺ ).
e
i
F are the generalized external forces acting on the system.
8
Disregarding dissipative forces and interpreting the external forces in
e
i
F
(2.1) as control variables (inputs), one obtains a control system. If the system is
under-actuated, a more general control system can also be defined as:
, 1, 2,...
0 , 1, 2,...
i
i i
u i m
d L L
i m m m dt q q
= ∂ ∂
− =
= + + ∂ ∂
ɺ
(2.2)
where
i
u are the control inputs.
In general, for an n-degree of freedom mechanical system with a Lagrangian
dependent on the control input u , and in the absence of other external forces, is
called a Lagrangian control system and the equations of motion are given as follows:
( , , , ) ( , , , )
0
i i
d L q q u t L q q u t
dt q q
∂ ∂
− =
∂ ∂
ɺ ɺ
ɺ
1,2,... i n = (2.3)
Defining the generalized momenta as:
( , , , )
i
i
L q q u t
p
q
∂
=
∂
ɺ
ɺ
1,2,... i n = (2.4)
The n-second order equations of motion in (2.3)- in terms of the generalized
coordinates and generalized velocities-, can be described as set of 2n first order
equations in terms of the generalized coordinates and generalized momenta by
defining the Hamiltonian as the Legendre transform of the Lagrangian such that:
1
( , , , ) ( , , , )
n
i i
i
H q p u t p q L q q u t
=
= −
∑
ɺ (2.5)
9
By taking the differential of the Hamiltonian defined above one obtains:
1 1
n n
i i i
i i
i i i
H H H H
dH dq dp du
q p u t
= =
∂ ∂ ∂ ∂
= + + +
∂ ∂ ∂ ∂
∑ ∑
(2.6)
Substituting the definition of the Hamiltonian in equation (2.5), in equation
(2.6) yields:
1 1 1 1 1
n n n n n
i i i i i i i
i i i i i
i i i
L L L L
dH q dp p dq dq dq du dt
q q u t
= = = = =
∂ ∂ ∂ ∂
= + − − − −
∂ ∂ ∂ ∂
∑ ∑ ∑ ∑ ∑
ɺ ɺ ɺ
ɺ
(2.7)
By the definition of generalized momenta in equation (2.5) it follows that:
1 1
0
n n
i i i
i i
i
L
p dq dq
q
= =
∂
− =
∂
∑ ∑
ɺ ɺ
ɺ
(2.8)
And by the Lagrange’s equation it follows that:
i
i
L
p
q
∂
=
∂
ɺ (2.9)
Substituting equations (2.8) and (2.9) into equation (2.7), reduces the
equation (2.7) to the form below:
1 1 1
n n n
i i i i
i i i
i i
L L L
dH q dp dq du dt
p u t
= = =
∂ ∂ ∂
= − − −
∂ ∂ ∂
∑ ∑ ∑
ɺ
ɺ
(2.10)
Moreover note that form equation (2.5) it follows that:
1 1
n n
i i
i i
i i
L H
du du
u u
= =
∂ ∂
=
∂ ∂
∑ ∑
(2.11)
Comparing equation (2.6) with equation (2.10), the following set of 2n+1
relation is obtained:
10
( , , , )
( , , , )
i
i
i
i
H
q q p u t
p
H
p q p u t
q
∂
=
∂
∂
=−
∂
ɺ
ɺ
( , , , ) L H q p u t
t t
∂ ∂
=−
∂ ∂
(2.12)
The equations in (2.12) without the control input u are called as the
canonical Hamiltonian equations, and the system with governing equations of
motions as above is called a Hamiltonian control system (see Nijmeijer and van der
Schaft [45]).
In addition one can easily show, by taking the total time derivative of the
Hamiltonian, that:
dH H
dt t
∂
=
∂
(2.13)
Thus, formally H − , and t can also be considered as a pair of canonical
variables.
2.1.1 Canonical Transformations
In this section the concept of canonical transformations in classical
mechanics is revisited for the control free Hamiltonian system described by the set of
equations in (2.12).
11
Definition 2.1 (Canonical Transformation): Consider the transformation of
coordinates from an independent set of coordinates -
i
q and momenta-
i
p , to a new
set of coordinates
i
q′ and
i
p′ defined as:
( , , )
( , , )
i i i
i i
q q q p t
p p q p t
′ ′ =
′ ′ =
(2.14)
Equation (2.14) states that the new coordinates are defined as function of
both old coordinates and old momenta.
If there exists a function ( , , ) F q p t ′ ′ such that the equations of motion are in
the Hamiltonian form:
i
i
i
i
F
q
p
F
p
q
∂
′=
′ ∂
∂
′=−
′ ∂
ɺ
ɺ
(2.15)
Then the transformation is said to be a canonical transformation. And the
function F is the new Hamiltonian in the new coordinate set. In other words, the
canonical transformations have the property of preserving the Hamiltonian structure.
In addition the new coordinates q′ and
i
p′ should satisfy the Hamilton’s least action
principle in the new and the old coordinates such that:
( )
2
1
( , , ) 0
t
i i
t
p q F q p t dt δ ′ ′ ′ ′ − =
∑
∫
ɺ (2.16a)
( )
2
1
( , , ) 0
t
i i
t
p q H q p t dt δ ′ ′ ′ ′ − =
∑
∫
ɺ (2.16b)
12
Thus the integrands above can differ by a total time derivative of an arbitrary
function K . Since the variation vanishes at the end points, the variation of the
integral of difference term is zero for any function K . This function K is called as
the generating function of the canonical transformation.
2.1.2 Lagrange and Poisson Brackets
The Lagrange bracket for any pair of dynamical quantity f and g is defined
as:
{ }
1
,
n
i i i i
i
q p p q
f g
f g f g
=
∂ ∂ ∂ ∂
= −
∂ ∂ ∂ ∂
∑
(2.17)
The Lagrange brackets defined above are invariant under canonical transformations
(i.e. Lagrange brackets are satisfied in any set canonical coordinates used). Choosing
the functions as the individual canonical coordinates results relations of the form:
{ }
, 0
i j
q q = (2.18a)
{ }
, 0
i j
p p = (2.18b)
{ }
,
i j ij
q p δ = (2.18c)
Equations (2.18) above are referred as fundamental Lagrange brackets (see
Goldstein [26])
Similarly the Poisson brackets are invariant under canonical transformations
and are defined as:
[ ]
1
,
n
i
i i i i
f g f g
f g
p q q p
=
∂ ∂ ∂ ∂
= −
∂ ∂ ∂ ∂
∑
(2.19)
13
And the following relations hold:
, 0
i j
q q =
(2.20a)
, 0
i j
p p =
(2.20b)
,
i j ij
q p δ =
(2.20c)
By choosing the function g as the Hamiltonian, f as the generalized
coordinates
i
q and the generalized momenta
i
p , one gets the canonical equations of
motion in terms of Poisson brackets below:
[ ] ,
i i
i
H
q H q
p
∂
= =
∂
ɺ (2.21a)
[ ] ,
i i
i
H
p H p
q
∂
=− =
∂
ɺ (2.21b)
The relations above are special cases for the total time derivative of a
function ( , , ) f q p t ,i.e.
1
n
i i
i
i i
df f f f
q p
dt q p t
=
∂ ∂ ∂
= + +
∂ ∂ ∂
∑
ɺ ɺ (2.22)
Equation (2.22) can be written in terms of the Hamiltonian as:
[ ]
1
,
n
i
i i i i
df f H f H f f
f H
dt q p p q t t
=
∂ ∂ ∂ ∂ ∂ ∂
= − + = +
∂ ∂ ∂ ∂ ∂ ∂
∑
(2.23)
If the function is chosen to be the Hamiltonian itself, from (2.23) we obtain
dH H
dt t
∂
=
∂
(2.24)
14
Thus the equations of motion (2.11) and (2.12), can be written with the
Poisson brackets such that:
[ ] ,
i i
q q H = ɺ (2.25a)
[ ] ,
i i
p p H =− ɺ (2.25b)
[ ] ,
dH H
H H
dt t
∂
= +
∂
(2.25c)
2.1.3 Invariants in Hamiltonian Mechanics
In section 2.1.1, canonical transformations are defined as transformations that
preserve the Hamiltonian structure. It is also of interest whether other expressions are
invariant under canonical transformations. One such set is known as Poincarè
integral invariants are presented in this section. In the following the presentation is
carried through using a time-independent Hamiltonian, although same results will
hold for the time-dependent Hamiltonian, by simply extending the phase space by
including time as an independent variable, with its conjugate variable H − .
2.1.3.1 Integral Invariants of Poincarè
The integral
i i
i
A
dp dq
∑
∫∫
is invariant under canonical transformations, where
the integral is evaluated over any arbitrary two-dimensional surface A in the phase
space. Similarly the integral
i i k k
i
s
dp dq dp dq
∑
∫∫∫∫
is an invariant under a canonical
15
transformation, where the integral is evaluated over any arbitrary four-dimensional
surface S in the phase space.
Integral invariants can be extended up to the phase space dimension 2n ,
leading to the invariant
1 1
... ... ...
n n
dp dp dq dq
∫ ∫
. The invariance of this last integral
under a canonical transformation, is just a statement of Liouville’s theorem, i.e.
volume of the phase space is preserved under a canonical transformation.
In modern treatments of the subject differential forms are used for
representation of the integral invariants of Poincarè stated above. (See Arnold [1])
2.1.4 Symplecticity of Poincarè Maps
From the invariance properties above, and using Stoke’s theorem, the first
Poincarè integral invariant for a Hamiltonian system can be written as:
( , )
i i i i
i i
C C
p dq p q
′
=
∑ ∑
∫ ∫
(2.26)
where C is a closed curve enclosing the area A .
For a tube of trajectories in phase space encircled by some closed curve C ,
we have
i i i i
i i
C C
p dq p dq
′
=
∑ ∑
∫ ∫
(2.27)
where C′ is the curve enclosing the tube of trajectories after it has evolved under a
canonical transformation.
16
The integrals in equation (2.27) above, corresponds to the sum of areas
projected onto the phase space planes ( , )
i i
p q . This is why the preservation is called
symplectic rather than area preserving. Actually, this shows that the Hamiltonian
flow itself is a canonical transformation itself.
Note that if the curves are defined on fixed two-dimensional surfaces in
phase space then the Poincarè Map (or surface of section) will be a 2n-2 dimensional
surface. Consider as an example, the surface is defined as
1
q =constant then (2.27)
becomes
1 1 1 1
2 2
n n
i i i i
i i
C C
p dq p dq p dq p dq
= =
′
+ = +
∑ ∑
∫ ∫
(2.28)
And the
1 1
p dq terms about both curves will be zero, thus (2.28) reduces to:
2 2
n n
i i i i
i i
C C
p dq p dq
= =
′
=
∑ ∑
∫ ∫
(2.29)
The above line integrals in (2.29) correspond to projected areas on the phase
space planes ( , )
i i
p q such that:
2 2
i
i
n n
i i i i
i i
A
A
dp dq dp dq
= =
′
=
∑ ∑
∫∫ ∫∫
(2.30)
where
i
A , and
i
A
′
are the projected areas on the ( , )
i i
p q planes of C and C′
respectively. Thus the Poincarè Map (or surface of section) is a symplectic mapping.
17
2.1.4 Periodic Hamiltonian Systems
For an explicitly time-dependent T -periodic Hamiltonian system, i.e.
( ) ( ) , , , , H p q t T H p q t + = (2.31)
The equations of motion for the system with the above Hamiltonian are
written as:
( , , )
i i
i
H
q Q p q t
p
∂
= =
∂
ɺ
( , , )
i i
i
H
p P p q t
q
∂
=− =
∂
ɺ 1,2,..., i N = (2.32)
Note that the functions
i
Q and
i
P are T -periodic, i.e. ( ( , , ) ( , , )
i i
Q p q t T Q p q t + = ,
and ( , , ) ( , , )
i i
P p q t T P p q t + = )
A point mapping that maps a phase point ( ) ( ), ( ) q t nT p t nT + + is to a phase
point ( ) ( ( 1) ), ( ( 1) ) q t n T p t n T + + + + is:
( 1) ( , )
( 1) ( , )
n n
n n
q n Q p q
p n P p q
+ =
+ =
1,2,..., i N = (2.33)
Defining a state vector ( ) ( ) ( ) ( )
T
x n q n p n = , Eq. (2.33) can be written in a
more compact form as:
( ) ( 1) ( ) x n G x n + = (2.34)
The mapping ( ) ( ) G x n is symplectic, and satisfies the Poisson brackets:
18
, 0
, 0
,
i j
i N j N
i j N ij
G G
G G
G G δ
+ +
+
=
=
=
1, 2,..., i n = ; 1, 2,..., j n = (2.35)
And the Liouville theorem can be restated for the mapping as:
det 1
G
x
∂
=
∂
(2.36)
Let H be the linear part of the mapping (2.34) linearized around a fixed point
of the mapping, i.e.
( )
x G x
∗ ∗
= given as:
x x
G
H
x ∗
=
∂
=
∂
(2.37)
The linear part of the mapping is also symplectic, and the matrix H is
symplectic, i.e.
0
T
H JH J − = (2.38)
where J is a 2 2 N N × skew-symmetric matrix defined as:
0
0
N
N
I
J
I
−
=
Note that the above is valid if the mapping is expressed in canonical
coordinates. However, it can be shown that for a linear system defined in some other
coordinate system ( 1) ( ) z n Az n + = is a Hamiltonian system if and only if there exists
an invertible skew-symmetric matrixW such that:
0
T
A WA W − = (2.39)
19
2.2 Properties of Symplectic Matrices
Symplectic matrices arise for linear Hamiltonian systems as presented in the
sections before, as well as in control theory for discrete time systems such as solution
of linear quadratic optimal control problems (see Chen [13]), the solution of
H
∞
problems, etc. In this section the properties of symplectic matrices are presented.
Definition 2.2: A 2 2 n n × matrix M is symplectic if and only if
T
M JM J = where
0
0
n
n
I
J
I
−
=
Theorem 2.2 (spectrum of symplectic matrices) If ( , ) x λ are an eigenpair of a
symplectic matrix M , then
1
( , ) Jx λ
−
are an eigenpair of
T
M . And the spectrum of
the matrix ( ) M Λ is symmetric with respect to the unit circle.
In addition to the theorem for the spectrum of a symplectic matrix, the
product of two symplectic matrices is also symplectic. And the determinant of a
symplectic matrix is unity, which is restatement of Liouville theorem.
2.3 Stability of Equilibrium Points of Nonlinear Periodically Time
Varying Dynamical Systems
Consider a non-autonomous nonlinear time periodic dynamical system,
which is described by a set of n equations:
( ) ( ) ( ), x t F x t t = ɺ (2.40)
where
n
x R ∈ is the state vector, and :
n n
F R R R
+
× → is a vector valued function
periodic in time with period T .
20
By taking snapshots of state ( ) x t at discrete time instances t nT = ,
0,1,2,... n= the surface of section, “Poincarè Mapping” of (2.1) is obtained as:
( ) ( ) ( 1) ( ) x n T P x nT + = (2.41)
An equilibrium state ( )
e
x t of (2.1) satisfying ( ) ( ), 0
e
F x t t = t ∀ , is obviously
a fixed-point of the mapping (2.41) and satisfies
( )
x P x
∗ ∗
= (2.42)
And the local stability property of the equilibrium state is determined by the
eigenvalues of the linear part of the mapping (2.41) given as:
x x
P
H
x ∗
=
∂
=
∂
(2.43)
The equilibrium state of (2.40) is then
i) Locally asymptotically stable if all the eigenvalues of H are inside the
unit disk
ii) Unstable if at least one of the eigenvalues of H is outside the unit disk
And if some of the eigenvalues of H are on the unit disk while the remaining
eigenvalues are inside the unit disk, the linear analysis is inconclusive. (See Jury
[29])
Similarly, a KT periodic solution of the continuous time system (2.40), is a
K -periodic solution of the map (1), (2),..., ( 1) x x x K
∗ ∗ ∗
− satisfying
( 1) ( ( )) x i P x i
∗ ∗
+ = 1, 2,.., 1 i K = − , and
( )
(1) ( (1))
K
x P x
∗ ∗
= (2.44)
21
And the local stability of the K -periodic solution is determined by the
eigenvalues of the matrix
( ) ( 1)... (1) H H K H K H = − (2.45)
where
( )
( )
x x i
P
H i
x ∗
=
∂
=
∂
The above results are straightforward and can be found in any textbook on
dynamical systems. However, in general the mapping in (2.41) cannot be found
exactly, thus one is forced to find an approximation to the mapping. This can be
achieved using the point-mapping method proposed by Flashner [19]. Equivalently,
(2.1) can be linearized around the equilibrium point yielding a linear periodic system
of the form:
( ) ( ) ( ) x t A t x t δ δ = ɺ (2.46)
where
( )
x x
F
A t
x ∗
=
∂
=
∂
, and ( ) ( ) A t A t T = + .
And the mapping for the linearized system above can be found using the
linear version of the point-mapping method such that (for details see Section 3.1):
( 1) ( ) x n Hx n + = (2.47)
22
By applying the Liouville theorem the determinant of the mapping is given
as:
( )
0
det( ) exp trace ( )
T
H A d τ τ
=
∫
(2.48)
The equation (2.48) is known as Liouville’s Formula.
23
Chapter 3
Discrete-Time Control of Linear Time-Periodic Systems
In this chapter, the proposed methodology for control design for the LPTV
systems is explained. The underlying idea is to reformulate the control problem by
parametrizing the control input, and applying point-mapping technique to obtain an
equivalent LTI discrete-time representation. Then, the known control design
techniques for LTI systems can be applied. Another advantage of the design
methodology used is that the designed controller can be implemented directly, thus
rendering the need for analysis effort required in discretization of the control law.
The organization of this chapter is as follows: In section 3.1 discrete time
representation of the LPTV control system is developed. In section 3.2, and section
3.3 general applicability of the method to with different control configurations are
shown following with a set of examples in section 3.4.
24
3.1 Discrete-Time Formulation of Linear Time-Periodic Control
Systems
Consider the problem of finding a control law ( ) t u for the dynamical system
( ) ( ) ( ) ( ) ( )
( ) ( ) ( )
t t t t t
t t t
x = A x + B u
y = C x
ɺ
(3.1)
where t R
+
∈ denotes time,
n
R ∈ x is the state vector,
m
R ∈ u is the control input
vector. The matrices ( )
n n
t R
×
∈ A , ( )
n m
t R
×
∈ B and ( )
n p
t R
×
∈ C are periodic with a
period T such that ( ) ( ) ( ) ( ) ( ) ( ) t +T t t +T t t +T t A = A , B = B , C = C .
Let the components of the control vector
i
u ( t ) be parameterized by a set of
predetermined functions
1 2
( ) ( ) ... ( )
i
i i il
t t t t ψ ψ ψ =
(i)
φ ( ) and a vector of
constants
1 2
...
i
T
i i il
a a a =
(i)
α as follows:
( )
1
( ) ( ) ( ) , 1,
i
l
i
i i ij ij
j
u t t t a i m φ α ψ
=
= = =
∑
… (3.2)
where ( ), 1, , 1,
ij
t i m j l ψ = = … … are piecewise functions defined over the period T .
This control parameterization allows for implementation of different
interpolators (data holds) used in discrete-time control systems. Let
( i )
i1 i
( t ) 1, l 1 φ ψ = = = (3.3)
i.e. a pulse function over the period T , then the parameterization models a zero-
order hold. It can also be interpreted as selecting the magnitude of the control over
the entire period at the beginning of each period.
25
By choosing a linear interpolating function with continuity constraint at the
end of the period such that
( i )
i1 i
t
( t ) 1 , l 1
T
= = − = φ ψ (3.4)
we get a first order hold. From control design point of view, this parameterization
(interpolator) can be viewed as determination of the (constant) rate of change of the
control effort over the period in the beginning of the period. In this case the control
magnitude is determined by the previous period. When the following set of functions
is chosen
[ ]
T
( i )
i1 i2 i
i1 i2
( t ) , l 2
t
1, 1
T
= =
= = −
φ ψ ψ
ψ ψ
(3.5)
both magnitude and rate of change of the control effort are determined in the
beginning of each period. Similarly, the parameterization of Eq. (3.2) allows
representation of higher order polynomial interpolators (holds) [23]. In particular one
can use basis functions such as cardinal basis for n
th
order Hermite polynomials and
B-splines [15].
Another possible choice of interpolating functions is a set of harmonics given
by:
[ ]
( )
( )
( )
T
( i )
i0 i2m
i0
i( 2 j 1 ) i 2 j
( t )
1,
cos j t , sin j t , j 1, ,m,
2
T
φ ψ ψ
ψ
ψ ω ψ ω
π
ω
−
=
=
= = =
=
⋯
…
(3.6)
26
In this case the interpolation functions can be interpreted as a multi-blade
transformation used developing rotorcraft control design [31].
Note that usually in discrete-time control systems identical digital to analog
converters (interpolators) are implemented in each control loop:
[ ]
( i ) ( k )
k j
2 l
( t ) ( t ) , l =l =l, i,k 1, ,m φ φ φ
φ ψ ψ ψ
1
= = =
=
…
…
(3.7)
The parametrization in Eq. (3.2) can be also used to model the variation of
the control effort within the periodT . Divide the period T into
s
m intervals of
duration
s
T T / m Δ = . Assume that within each interval T Δ the control effort is
parameterized as:
( ) ( ) ( ) ( ) ( )
1
( ) ( ) ( ) ( )
1 2
( ) ( ) ( ) , 1,
( ) ( ) ( ) ( )
r
r
l
r r T r r r
i i j ij s
j
r r r r
l
u t t t a r m
t t t t
φ α ψ
φ ψ ψ ψ
=
= = =
=
∑
…
⋯
(3.8)
Then the control over the entire period T can be written as in Eq. (3.2):
( ) (1) (2)
(1)
(2)
( )
( ) ( ) ( ) ( ) ( ) , 1,
:
s
s
m
i i i
i
i
i
m
i
u t t t t t i m φ α φ φ φ α
α
α
α
α
= = =
=
⋯ …
(3.9)
Note that one can impose continuity conditions on the control effort by
imposing continuity at the limits of the subintervals. In this case, choosing splines as
basis functions is convenient. It should also be noted that one is not restricted to
27
using polynomials as basis functions, and other classes of functions such as
harmonics can also be employed.
Eq. (3.2) can be rewritten in a more compact form as:
(1)
(2)
( )
( ) ( ) ,
( )
( )
( )
( )
m
t t
t
t
t
t
=
=
k
u Ξ a
φ 0 0
0 φ
Ξ
0
0 0 φ
⋯
⋱ ⋮
⋮ ⋮ ⋱
⋯
(3.10)
where
u
n
R ∈
k
a is a constant vector over the period, and ( )
u
m n
t R
×
∈ Ξ is a known
matrix.
Then the state equations of dynamical system are given by:
( ) ( ) ( ) ( )
k
t t t t = + x A x F a ɺ (3.11)
where ( ) ( ) ( ) t t t F = B Ξ , and satisfies ( ) ( ) t t +T F = F .
Note that this approach can be easily extended to systems with different
sampling periods in each control loop yielding in multi-rate systems.
Consider a mapping of the state ( ) t x from t 0 = to t T = , called also point-
mapping, then the system in Eq. (3.11) is formulated as a discrete-time system
represented by a set of difference equations:
for 0,1, 2...
k+1 k k
k = x = Hx + Ga (3.12)
where
1 k+
x , and
k
x denote the state of the system at t=t
0
+(k+1)T, and t= t
0
+kT
respectively.
28
Note that (3.12) is a linear time invariant discrete time system. Formally, the
solution of (3.1) is:
0 0
0
( ) ( , ) ( , ) ( ) ( )
t
x t t t x t B u d τ τ τ τ =Φ + Φ
∫
(3.13)
where
0 0
( ) x t x = -initial state,
0
( , ) t t Φ - the transition matrix
Yielding
0 0
( , ) t T t =Φ + H and
0
0
( , ) ( )
t T
t
t F d τ τ τ
+
= Φ
∫
G (3.14)
Since an analytical solution for H and G is possible only in special cases, a
computational method proposed by Flashner and Hsu [19, 21], is modified to
compute the matrices H and G .The constant matrices H and G are computed
using a computational algorithm based on Runge-Kutta method (see Appendix-A)
yielding
1 1
1 1
1
...
( )
N N
i N
N j i
i j N
−
+ −
= =
=
= +
∑ ∏
H Φ Φ Φ
G Ψ Φ Ψ
(3.15)
where
4
( )
0 0
1
( , ( 1) )
j
m
j m
m
t jh t j h I h d k
=
= + + − = +
∑
Φ Φ
,
( )
( ) ( ) ( 1)
j j j
m m m
m
k A I c A
−
= +
4
( )
1
j
m
j m
m
h d k
=
=
∑
Ψ ,
( ) ( ) ( ) ( 1)
j j j j
m m m m
m
k F c A F
−
= +
And the constants
m
c and
m
d are Runge-Kutta method constants
29
3.1.1 Examples
To demonstrate the effectiveness of the discretization algorithm, the
following examples are carried out.
3.1.1.1 Undamped and Unforced Mathieu Equation
Consider π-periodic Mathieu equation
0 1
1 0.2cos(2 ) 0
x x
t
=
− −
ɺ (3.16)
Since (3.16) is a Hamiltonian system, the point mapping should be an area-
preserving more precisely a symplectic map. Thus the determinant of the map should
be unity. The state transition matrix computed by point-mapping algorithm H , and
determinant of the mapping for various number of integration step sizes are shown in
table 3.1.
Table 3.1: Effect of step size on accuracy of the computational method
Number of
Integration
Steps
H
det( ) H
N=50 -1.02317556994653 0.15462570195382
0.15953546115221 -1.00145885759592
0.99999995473152
N=100
-1.01776449012161 0.15495669065366
0.15992389247579 -1.00689431166621
0.99999999858481
N=200
-1.01505072743651 0.15503957501474
0.16002099694562 -1.00961415978075
0.99999999995577
N=400
-1.01369243851069 0.15506030470879
0.16004527286985 -1.01097396985882
0.99999999999862
N=800
-1.01301302304190 0.15506548767087
0.16005134183947 -1.01165376561233
0.99999999999995
30
3.1.1.2 Underdamped and Unforced Mathieu Equation
Consider π-periodic under damped Mathieu equation
0 1
1 0.2cos(2 ) 0.2
x x
t
=
− − −
ɺ ,
1
(0)
1
x
=
(3.17)
The trajectories of the states obtained by numerical integration and the states
at discrete time instances (t kπ = , 1,2,... k = ) obtained by point-mapping algorithm
is shown in Fig. 3.1 below:
Figure 3.1: Comparison of trajectories: point-mapping versus numerical integration
31
3.2 Control System Formulation
The system in Eq. (3.11) represents a continuous plant controlled by a control
k
a updated every period T . The control law interfaces with the system by a
generalized hold (interpolator) represented by the matrix ( )
m l
t R
×
∈ Ξ as shown in
Fig. 3.2.
Figure 3.2: Block diagram representation of T -periodic sampled data system
The output vector at the sampling instant is given by:
k k k
= y C x (3.18)
An equivalent discrete-time representation of the system is shown in Fig. 3.3.
Figure 3.3: Block diagram for the discrete-time state variable system
( ) t Ξ
( ) B t
∫
( ) A t
k
a
( ) u t
( ) x t
k
x
+
+
T T
a
k
+
+
k
x
Delay
C
H
G
y
k
1 k+
x
32
Assuming that the output variables are measured at time instances t kT = (for
k=1,2,3 …), the discretization procedure in Eq. (3.15), transforms the LPTV control
problem in Eq. (3.1) into equivalent feedback discrete-time control design problem
shown in Fig 3.4. The objective is to find a control law ( )
k k
K e = a that satisfies the
continuous- time performance specifications.
Figure 3.4: Discrete-time formulation of the closed loop system
3.2.1 Fast Sampling Control System Formulation
The proposed formulation is also applicable to systems with sampling faster
than the system’s periodicity. Assume output sampling period, /
s s
T T m = , and
introduce the variables
s
n m
k R
×
∈
_
x , 1
s
n m
k R
×
+ ∈
_
x ,
s
l m
R
×
∈
_
a ,
s
r m
R
×
∈
_
y as:
1
0 0 0
0 0 0
{ ( ) ( / ) ... ( ( 1) / )}
{ ( / ) ( 2 / ) ... ( ( 1) }
k
k
T
s s s
T
s s
t kT t kT T m t kT m T m
t kT T m t kT T m t k T
+
= + + + + + −
= + + + + + +
_
_
x x x x
x x x x
0 0 0
{ ( ) ( / ) ... ( ( 1) / )}
T
k
s s s
t kT t kT T m t kT m T m = + + + + + −
_
a a a a
0 0 0
{ ( ) ( / ) ... ( ( 1) / )}
T
s s s k
t kT t kT T m t kT m T m = + + + + + −
_
y y y y (3.19)
-
r(t)
y(t)
K(.) ( ) t Ξ
System
T
33
The discrete time formulation of system dynamics given is:
1 k k k +
= +
_ _ _
x H x G a (3.20)
where { }
i
diag H = H , and { }
i
diag G = G
0 0
( , ( 1) )
i
s s
T T
H t i t i
m m
=Φ + + − ,
0
0
( 1)
( , ) ( )
s
s
T
t i
m
i
T
t i
m
G t F d τ τ τ
+
+ −
= Φ
∫
1, 2,...,
s
i m =
The computation for
i
H ,and
i
G are then carried out using the point-mapping
algorithm, and the output vector is given as:
k
k
=
_ _
y Cx (3.21)
where { }
i
diag C = C and ( ) ( )
0
1 /
i s
C C t i T m = + − for i ∀ .
After computation of the equivalent system matrices defined above, the
control design problem for the periodic system in Eq. (3.1) can then similarly be
formulated as a feedback sampled-data control design problem as shown in Fig. 3.5.
In the rest of this dissertation, a system with sampling faster than system periodicity
will be referred to as fast-sampled system.
Figure 3.5: Discrete-time formulation of the closed loop system with fast sampling
-
r(t)
y(t)
K(.) ( ) t Ξ
System
T/m
s
+
34
A more general framework for the control problem will arise if the measured
physical signals are sampled with different sampling rates. This is usually
encountered in complex multivariable control systems for several reasons. First it
may be impossible to sample all the physical signals at one rate. Second it may be
more desirable to sample some signals with different bandwidth with different rates
to have a better trade-off between performance and implementation cost. The multi-
rate sampling shown below correspond to sampling r-channels of y periodically with
periods
i
T
m
, 1, 2,..., i r = respectively and
i
m ’s are integers. By choosing appropriate
new state variables, the equivalent discrete time LTI system can be found using the
point-mapping algorithm.
Figure 3.6: Multi-rate sampling
3.3 Controller Design
The discretization algorithm to formulate the equivalent LTI control system
design problem for LPTV control systems allows us to use the well known design
techniques such as pole placement, loop shaping, linear quadratic regulator, H
2
and
:
:
T/ m
1
T/ m
r
y
1
(t)
y
1
(k T/ m
1
)
y
r
(k T/ m
r
)
( ) ( ) ( ) ( ) ( )
( ) ( )
x t A t x t B t u t
y t Cx t
•
= +
=
35
H
∞
for both formulations described above. To demonstrate the versatility of the
approach as well as its ease in designing control laws for the LPTV systems, we will
consider full state and output feedback regulator design problems in the following
section.
3.3.1 Full State Feedback Optimal Regulator Design
As discussed before, known discrete-time methods, such as pole placement or
LQR, can be used to design a controller for the discrete-time representation of the
LPTV system. However, using this approach the inter-sample behavior of the control
system is not included in the design. To address this issue, an optimal control design
approach is modified to include the inter-sampling behavior by employing the point-
mapping discretization method introduced before. This approach is presented in the
following.
Inter-sample behavior of the control system and evaluation of performance
with respect to a class of control vector parameterizations can be handled by the
following performance index:
( )
0
lim ( ) ( ) ( ) ( ) ( ) ( )
f
T T T
t
J t t t t t t dt
∞
→∞
= +
∫
e Q e u R u (3.22)
where the error is defined as ( ) ( ) ( ) t t t = − e r y
The optimal controller will be mathematically defined as an optimization as
follows:
( )
min
opt
u t
u J = (3.23)
36
Subjected to the plant dynamics given as:
( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( )
t t t t t t
t t t
= + +
=
x A x B u r
y C x
ɺ
(3.24)
where ( )
n n
t
×
∈ A ℝ , ( )
n m
t
×
∈ B ℝ , ( )
p n
t
×
∈ C ℝ are periodic in time
To simplify derivation, consider a finite-time optimal regulator problem, i.e.
assume ( ) 0 t = r , given by
{ }
0
1
( )
min ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
f
t
T T T
f f
u t
t
J t t t t t t dt t P t = + +
∫
y Q y u R u x x (3.25)
such that (3.24) holds. With t
0
= 0, and t
f
→∞, the performance criteria stated in Eq.
(3.22) is recovered.
Using the discretization approach presented before, it can be shown that
assuming output sampling period T, where T is system periodicity, the continuous
optimal regulator problem is equivalent to the discrete-time problem with
performance index
{ }
1
1
0
2 ( ) ( )
f
f f
k k
T T T T T
k k k k k k k f f k
k
J Q Z R C t PC t
= −
=
′ ′ = + + +
∑
x x x a a a x x (3.26)
Subject to the discrete-time dynamics given in Eq. (3.11), where
( 1)
( ) ( ) ( ) ( ) ( )
k T
T T
kT
Q t t t t t dt
+
′=
∫
Θ C Q C Θ
( )
( 1)
( ) ( ) ( ) ( ) ( )
k T
T T
kT
Z t t t t t dt
+
=
∫
Θ C Q C Γ
37
( )
( 1)
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
k T
T T T
kT
R t t t t t t t t dt
+
′= +
∫
Γ C Q C Γ Ξ R Ξ
( ) ( , )
k k
t t t t Θ = + Φ
( ) ( , ) ( ) ( )
k
k
t t
k k
t
t t t d τ τ τ τ
+
Γ = +
∫
Φ B Ξ
The integrals for performance weights above are computed using a numerical
quadrature by employing the transition matrix and the matrices H, and G computed
during system discretization by the proposed point-mapping algorithm at time
instances t
n
= nT/N. The solution to the optimal regulator problem
k k k
K =− a x
(k=1,2, … , k
f
)
is then found invoking dynamic programming [5] as:
( ) ( )
1
1 1
T T
k k k
K P R P Z
−
+ +
′ ′ = + + G G G H (3.27)
where
k
P is computed is computed using Riccati equation is given by:
( )( ) ( )
1
1
T T T T T
k k k k k k
P P P P Z P R P Z Q
−
+
′ ′ = − − + + + + H H H G G G G H (3.28)
with the terminal condition
1
f
k
P P =
It should be noted that the above design methodology can be applied to the
case when the sampling period is faster than the periodicity shown in Fig. 3.5.
3.3.2 Output Feedback Regulator Design
An output feedback control system, calls for the design of compensator that
includes design an estimator (observer), and a feedback gain that uses the estimated
states to generate the control input sequence [13]. In the following, two approaches
38
for the controller-estimator configuration differ in sampling and control law update
rates will be presented.
3.3.2.1 Synchronous Sampling Observer and Control Law Update
Consider the discrete-time representation of the LPTV plant with a controller
update period equal to system periodicity T, see Eq. (3.11). Then the state vector can
be estimated using an observer given by:
1
ˆ ˆ ˆ ( )
k k k k k
L
+
= + − + x Hx y Cx Ga (3.29)
where
n r
L R
×
∈ is an observer gain to be determined.
The feedback control law is given by:
ˆ
k k
K =− a x (3.30)
The above closed-loop control configuration of periodic system is shown in
Fig. 3.7.
Figure 3.7: Controller-observer configuration
K
( ) t Ξ
Time-Periodic
System
Observer
Equation (25)
k
y
k
a
ˆ
k
x
( ) t y
T
T
39
The separation principle ensures that the design of state feedback and the
observer can be carried out independently.
3.3.2.2 Fast Sampling Observer
The output sampling rate can be increased to improve state estimation
accuracy, thus enhancing the closed-loop system performance. Let the output
sampling time be given by /
s s
T T m = , i.e. m
s
times faster than system’s periodicity,
and define additional state variables as follows:
0 0 0
1
0 0 0
ˆ ˆ ˆ ˆ { ( ) ( / ) ... ( ( 1) / )}
ˆ ˆ ˆ ˆ { ( / ) ( 2 / ) ... ( )}
T
k
s s s
T
k
s s
t kT t kT T m t kT m T m
t kT T m t kT T m t kT T +
= + + + + + −
= + + + + + +
_
_
x x x x
x x x x
Then the observer’s state equation is:
1
ˆ ˆ ˆ { }
k k k k k
L
+
= + − +
_ _ _ _
x H x y C x G a 1, 2,...,
s
i m = (3.31)
where { }
i
diag H = H , { }
i
diag G = G , { }
i
diag C = C ,
0 0
( , ( 1) )
i
s s
T T
H t i t i
m m
=Φ + + − ,
0
0
( 1)
( , ) ( )
s
s
T
t i
m
i
T
t i
m
G t F d τ τ τ
+
+ −
= Φ
∫
( ) ( )
0
1 /
i s
C C t i T m = + − 1, 2,...,
s
i m =
It is easy to show that by augmenting the observer equations above with the
control law in equation (3.30), and applying a similarity transformation that the
separation principle holds for the controller-fast observer connection.
40
In the formulation of both synchronous and fast sampling observers, it was
assumed that the controller update period equals the LPTV periodicity. Both designs
can be easily modified to allow controller update period faster than system’s
periodicity.
3.4 Examples
3.4.1 Dead-beat Full State Regulator for Mathieu Equation
Consider π-periodic Mathieu equation
0 1
1 0.2cos(2 ) 0
x x
t
=
− −
ɺ
0
( )
1
u t
+
(3.32)
It is assumed that the entire state vector is are available at the start of each
period
k
t kT = . Deadbeat controller designs for different choices of basis functions--
( )
ij
t φ ’s—(see Equation (3.2)) are presented in the following.
Figure 3.8: Full state feedback regulator block diagram
∫
( ) Bt
( ) A t
( ) x t ( ) x t
•
+
( ) u t
K
( ) t Ξ
( ) 0 r t =
Compensator D (z)
+
+ -
41
i) Pulse Function Parameterization
Here the control input is parameterized as follows:
( ) ( ) ( )
k k
u t t t φ =Ξ = a a (3.33)
where
1 ( 1)
( )
0
k T t kT
t
else
φ
− ≤ <
=
(3.34)
The equivalent discrete time dynamics computed by point-mapping algorithm
is:
1
-1.0137 0.1551 2.1345
0.1600 -1.0110 -0.1675
k k k +
= +
x x a (3.35)
The state feedback gains ( )
1 2
K k k = for the dead-beat discrete-time
controller
k k
K =− a x for the discrete-time model given in Eq. (3.35) are the gains
that place the closed loop poles at the origin in the complex plane, and is computed
such that:
( )
2
det zI H GK z − + = (3.36)
Equation (3.36) yields two equations with two unknown gains to be
determined as follows:
1 2
1 2
2.1345 0.1675 2.0247 0
2.1320 0.1717 1 0
k k
k k
− + =
+ + =
(3.37)
Then the feedback gains are the unique solution to the set of linear equations
above, and are calculated as:
( ) -0.7120 3.0144 K =
42
The closed loop system response and the control effort are shown in Fig. 3.9
below.
Figure 3.9: Closed loop response, and control effort for dead-beat regulator with
pulse function parameterized controller for Mathieu equation
43
ii) Linear Parameterization
Here the control input is parameterized as:
1 2
( ) ( ) [ ( ) ( )]
k k
u t t t t ψ ψ =Ξ = a a (3.38)
where
1
1 ( 1)
( )
0
k T t kT
t
else
ψ
− ≤ <
=
,
2
( 1)
( )
0
t k T t kT
t
else
ψ
− ≤ <
=
By point-mapping the equivalent multi-input discrete system is:
1
1.0137 0.1551
0.1600 1.011
k k +
−
=
−
x x +
2.1345 3.1952
-0.1675 1.7730
k
a (3.39)
The dead-beat controller gain for the multi-input case is not unique. The
additional degrees of freedom to select the control gains can be used to impose
additional conditions. From a practical point of view, one of the advantageous
choices is to determine gain matrices with small gains to achieve less control effort.
On the other hand, one of the important aspects of selection of the feedback gain
matrix is the sensitivity of the closed loop poles to the perturbations in the state or
the input matrices. The latter approach is known as robust pole assignment problem
and several algorithms are available. Here we use the algorithm proposed by Kautsky
and Nichols [37] that optimizes the choice of eigenvectors for a robust solution.
The dead-beat controller gain for the system in Eq. (3.39) is given by:
-0.5345 0.8115
0.0398 -0.4935
K
=
44
The closed loop system response and the control effort are shown in Fig. 3.10
below.
Figure 3.10: Closed loop response, and control effort for dead-beat regulator with
linear function parameterized controller for Mathieu equation
45
iii) Quadratic Parameterization
If the control input is parameterized as:
[ ]
1 2 3
( ) ( ) ( ) ( ) ( )
k k
u t t t t t ψ ψ ψ =Ξ = a a (3.40)
where
1
1 ( 1)
( )
0
k T t kT
t
else
ψ
− ≤ <
=
2
( 1)
( )
0
t k T t kT
t
else
ψ
− ≤ <
=
2
3
( 1)
( )
0
t k T t kT
t
else
ψ
− ≤ <
=
By point-mapping the equivalent multi-input discrete system is:
1
1.0137 0.1551
0.1600 1.011
k k +
−
=
−
x x +
2.1345 3.1952 5.8774
-0.1675 1.7730 5.8972
k
a (3.41)
The feedback gain matrix is then computed as:
-0.3794 0.3952
-0.2101 0.1773
0.0795 -0.2135
K
=
The closed loop system response and the control effort are shown in Fig.
3.11.
46
Figure 3.11: Closed loop response, and control effort for dead-beat regulator with
quadratic function parameterized controller for Mathieu equation
47
From the simulations, first, it is seen that all of the three controllers stabilize
the open loop unstable plant. The first design using pulse function parameterized
control achieves approximately dead-beat response in two steps showing the
accuracy of the discretization algorithm. The second and the third designs with linear
and quadratic parameterized control vector results dead-beat response in one step,
and the tracking error has been reduced considerably for the second and the third
design. The required control power, (
2
( ) u t dt
∫
), has been reduced approximately by
43 %, and 34% with the second and the third design respectively.
3.4.2 Dead-beat Output Feedback Regulator for Mathieu Equation
3.4.2.1 Dead-beat Regulator with Synchronous Sampling Observer
Consider again the example discussed in section 3.4.1. Assume that only the
position is measured, i.e. ( ) (1 0) ( ) y t x t =
Provided that the pair ( , ) H C is detectable, the dead-beat observer gains for
the three different choices of basis functions used in the examples are:
( ) -2.0247 6.7515 L=
΄
The dead-beat feedback gains for the pulse function, linear function, and
quadratic function parameterized control inputs are same as before. The simulated
responses of the system with estimated state feedback for the three cases are shown
in Fig. 3.12, Fig. 3.13, and Fig. 3.14.
48
Figure 3.12: State and estimated state histories for Mathieu equation with estimated
state feedback control (Pulse-function Parameterization)
Figure 3.13: State and estimated state histories for Mathieu equation with estimated
state feedback control (Linear parameterization)
49
Figure 3.14: State and estimated state histories for Mathieu equation with estimated
state feedback control (Quadratic Parameterization)
The above simulation shows that the estimation error for both states goes to
zero in two periods. As the theory suggests the controller-estimator configuration
from all the three cases have dead-beat response. The dead-beat response is achieved
in nine periods with the pulse function parameterized controller compared to two
periods with the full state feedback case. The dead -beat response is achieved in 3
periods for the linear and quadratic function parameterized controllers compared to
one period with the full state feedback.
50
3.4.2.2 Deadbeat Regulator with Fast Sampling Observer
If the position output for the system is measured twice a period, then using
the observer equation (see Eq. (3.31)) fast sampling dead-beat observer gains
1
L and
2
L can be calculated as:
( )
1
-0.0035 -1.0011 L =
T
, ( )
2
-0.0004 -0.9996 L =
T
where
1
H ,
2
H are same for all the three different basis functions and are computed
by the point-mapping algorithm as:
1
-0.0803 0.9910
-1.0070 0.0767
H
=
,
2
0.0784 1.0004
-1.0058 0.0788
H
=
1
G ,and
2
G for the pulse-function, linear function, and quadratic function
parameterized control input are computed by the point mapping algorithm and are
shown in the table below:
Table 3.2: Pulse, linear, and quadratic function parameterized control
input matrices for fast sampling observer
Pulse-Function
Parameterization 1
1.0116
1.0164
G
=
2
0.9884
0.9339
G
=
Linear Function
Parameterization 1
1.0116 0.5771
1.0664 1.0326
G
=
2
0.9884 2.1171
0.9339 2.4348
G
=
Quadratic
Function
Parameterization
1
1.0116 0.5771 0.4717
1.0664 1.0326 1.1648
G
=
2
0.9884 2.1171 4.6753
0.9339 2.4348 6.4634
G
=
51
The simulated response of the system with estimated state feedback from fast
sampling observer, with the three different choices of basis functions as control
vector parameterization before are shown in Fig. 3.15, Fig. 3.16, and Fig. 3.17.
As can be seen from the simulations, the estimation error with fast sampling
reduces to zero in one period. And the dead-beat response is achieved in four periods
for the linear function parameterized controller with fast sampling observer. The
dead-beat response is achieved in three periods with the linear and quadratic function
parameterized controllers. It is clear that the transient response is greatly improved
with fast sampling observer.
Figure 3.15: State and estimated state histories for Mathieu equation with estimated
state feedback using fast sampling observer (Pulse-function parameterization)
52
Figure 3.16: State and estimated state histories for Mathieu equation with estimated
state feedback using fast sampling observer (Linear function parameterization)
Figure 3.17: State and estimated state histories for Mathieu equation with estimated
state feedback using fast sampling observer (Quadratic function parameterization)
53
3.4.3 Multivariable Feedback Design Example: Control Design for
Periodically Excited Inverted Double Pendulum
The equations of motion for the inverted double pendulum shown in Fig.
3.18 (see Bolotin [8]) acted upon by periodic loading P(t) as shown are:
( )
2 2
1 2 1
2 2
2 1
2
2
2
2 2 1 1 2 2 1 1
2
2
1 2 1 1 2 2 2
3 cos( )
cos( )
sin( ) 2 sin( )
sin( ) sin ( 1)
ml ml
ml ml
ml k k Pl u
ml k k Pl u
θ θ θ
θ θ
θ
θ θ θ θ θ αθ θ
θ θ θ θ θ α θ
−
−
− − + − − +
=
− + − − − +
ii
ii
i
i
(3.42)
where α is the load-direction parameter, k is the spring stiffness,
1
u
, and
2
u are
control torques applied to the bottom and top masses, respectively.
Figure 3.18: Inverted double pendulum subjected to periodic forcing
54
Linearizing the equations of motion around the equilibrium position
1 2
0 θ θ = = yields:
1 1 2
0.5 ( 3) 0.5 (2 ) k p k p θ θ θ = − + −
ɺɺ
2 2
1 2
/ 2 / 2 u ml u ml + −
( ) ( ) { }
2 2
2 1 2 1 2
0.5 (5 ) 1.5 2 /2 3 /2 k p p k u ml u ml θ θ α θ = − + − − − +
ɺɺ
(3.43)
where
2
/ k k ml = ,and ( )
1 2
cos( ) / p P P t l k ω = +
Define the state vector as
( )
1 2 1 2
( ) ( ) ( ) ( )
T
t t t t θ θ θ θ = z(t)
ɺ ɺ
. Then the
dynamic equations of motion in state space form are:
( ) ( )
1
2
2
0 0 1 0
0 0
0 0 0 1 ( ) 0 0
1
( ) ( )
0.5 ( 3) 0.5 (2 ) 0 0 ( ) 1/ 2 1/ 2
0.5 (5 ) 1.5 2 0 0
1/ 2 3/ 2
u t
t t
k p k p u t ml
k p p k α
= +
− −
−
− − −
−
z z ɺ
(3.44)
The nominal values of the parameters are k =2,
1
/ Pl k =2,
2
/ P l k =0.7, α =1,
ω =2. Using the point-mapping algorithm the nominal state transition matrix defined
in Eq. (3.11) (T π = ) is:
-0.1170 -0.2154 0.0598 0.3238
-0.7841 -1.3004 2.3811 -0.6260
-0.1877 0.5808 -1.6615 0.2994
-2.4894 0.8357 -1.2989 0.2442
=
H
The eigenvalues of H (Floquet multipliers) are: λ
1,2
= 0.0118±0.9999i, λ
3
=
2.45 , λ
4
= 0.4082, yielding an unstable open loop system.
55
The two sets of basis functions considered here are:
i) Linear Function Parameterization
1
2
1 ( 1)
( )
0
( 1)
( )
0
i
i
k T t kT
t
else
t k T t kT
t
else
ψ
ψ
− ≤ <
=
− ≤ <
=
1,2 i= (3.45)
ii) Quadratic Function Parameterization
1
( )
i
t ψ , and
2
( )
i
t ψ same as in Eq. (3.45) and
2
3
( 1)
( )
0
i
t k T t kT
t
else
ψ
− ≤ <
=
1,2 i= (3.46)
3.4.3.1 Full State Feedback Regulator Designs
Assuming that the entire state vector is available at times t
k
=kT, the input
transition matrix G for the linear function and quadratic function parameterized
control vectors were computed by the point-mapping algorithm (see Eq. (3.11)) as:
0.6662 1.3790 -0.2881 -1.2017
1.5422 0.6932 -0.5834 1.4911
=
-0.1965 1.1811 0.0610 -1.6409
0.8269 1.0274 -0.4089 0.6604
linear
G ,
0.6662 1.3790 2.8463 -0.2881 -1.2017 -2.7982
1.5422 0.6932 -0.0672 -0.5834 1.4911 4.7126
-0.1965 1.1811 3.8057 0.0610 -1.6409 -4.8021
0.8269 1.0274 0.6518 -0.4089 0.6604 4.6552
quadratic
=
G
56
3.4.3.1.1 Dead-beat Controller Design
Dead-beat controller gain matrices for the discrete-time plant with the linear
and quadratic function parametrized controllers were computed as:
6.2619 -6.0002 8.9030 -2.8206
-1.5546 0.7514 -1.9477 -0.3140
12.7795 -11.5194 14.8272 -7.6650
-1.2795 0.4774 -0.9045 -0.3555
linear
K
=
0.8124 -1.2778 2.2869 -0.2359
-0.3436 0.2046 -0.2434 0.3391
-0.3855 0.1740 -0.5425 -0.2079
=
-0.1759 0.3203 -0.6165 -0.0239
0.8922 -1.2079 1.8238 -0.9051
-0.6913 0.5364 -0.8684 0.1749
quadratic
K
3.4.3.1.2 Optimal Controller Design
Selecting performance weights as Q= I
(4x4)
, R= I
(2x2)
, and using the
trapezoidal approximation for calculating the weights for the equivalent discrete-time
regulator problem given in Eq. (3.26), the steady-state optimal control gains for the
linear and quadratic function parameterized controllers are obtained as:
0.8606 -1.2543 2.1368 -0.0779
-0.7839 0.6246 -1.2051 0.1265
1.5872 -0.5709 0.6421 0.3705
-0.5896 0.2007 -0.1553 -0.1604
linear
K
=
57
1.2985 -1.2624 2.3715 0.0080
-2.0137 0.9349 -2.0961 0.1376
0.5340 0.9349 -2.0961 -0.0618
0.7839 0.1189 -0.3171 0.8482
0.6813 -0.9598 1.4069 -1.0057
-0.3794 0.3519 -0.4682 0.2636
quadratic
K
=
Simulation studies using the linearized model show that all four controllers
defined above stabilized the plant, and a dead-beat response is achieved after one
period for both dead-beat controllers.
The simulation studies of the controlled system using the non-linear model,
however, show a qualitative difference between linear and quadratic function
parametrized dead-beat controllers. The quadratic function parameterized dead-beat
controller stabilizes the non-linear plant for a significantly larger range of initial
conditions. The simulated response of the linear function parametrized controller
with a small initial condition
1
0.05 θ = rad and
2
0.1 θ = rad is shown in Fig. 3.19.
58
0 1 2 3 4
-0.1
-0.05
0
0.05
0.1
time (x T sec)
θ
1
(ra d )
0 1 2 3 4
0
0.2
0.4
time (x T sec)
θ
2
(ra d /s)
0 1 2 3
0
0.1
0.2
0.3
time (x T sec)
u
1
0 1 2 3 4
0
0.2
0.4
0.6
time (x T sec)
u
2
Figure 3.19: Closed loop state trajectories and control forces for the dead-beat linear
function parameterized controller
For comparison purposes, the simulated responses, and the control forces
with quadratic function parameterized dead-beat controller, linear and quadratic
function parameterized controllers with initial conditions
1
0.3 θ = rad and
2
0.25 θ =
rad are shown Fig. 3.20.
59
0 1 2 3
-0.1
0
0.1
0.2
0.3
time (x T sec)
θ
1
(rad)
0 1 2 3
-0.2
0
0.2
0.4
time (x T sec)
θ
2
(rad/s)
Quadratic Dead-beat
Quadratic Optimal
Linear Optimal
0 1 2 3
0
0.5
1
time (x T sec)
u
1
0 1 2 3
0
0.5
1
time (x T sec)
u
2
Quadratic Dead-beat
Quadratic Optimal
Linear Optimal
Figure 3.20: Closed loop state trajectories and control forces for the dead-beat
quadratic, steady state optimal linear, and quadratic function
parametrized controllers
60
3.4.3.2 Fast Sampled Full State Feedback Regulator Designs
Assuming the sampling period is half of the base period of the system (i.e.
/ 2
s
T T = ) , the discrete time representation of the linear periodic system given as
in Eq. (3.44) is:
1 k k k +
= +
_ _ _
x H x G a (3.48)
where
{ }
{ }
1
( ) ( / 2)
( / 2) ( )
k
k
T
T
kT kT T
kT T kT T
+
= +
= + +
_
_
x z z
x z z
{ } ( ) ( / 2)
T
k kT kT T = +
_
a a a
1 2
{ , } diag H H = H , and
1 2
{ , } diag G G = G
The computed matrices
1 2
, H H and
1 2
, G G (see Eq. (3.20)) for linear and
quadratic function parametrized controllers are tabulated in tables 3.3 and 3.4
respectively.
Table 3.3: State transition matrices for fast-sampled state feedback design
H
1
H
2
0.3486 -0.2927 1.0518 -0.0151
1.8275 -0.6908 1.3364 0.5016
-0.6098 -0.0740 -0.0149 0.1941
0.6519 -1.7176 1.9714 -0.9209
-0.3677 0.3117 0.9397 0.0222
1.8538 -0.5681 1.2991 0.6137
-1.5839 0.2507 -0.0256 -0.1680
0.3272 -0.7435 1.7028 -0.3166
61
Table 3.4: Input transition matrices for fast-sampled state feedback design
Linear Function
Parametrized Controller
Quadratic Function Parametrized
Controller
G
1
0.4721 0.2671 -0.4375 -0.2469
-0.0684 -0.1336 0.8275 0.6177
0.3419 0.4104 -0.1496 -0.3107
0.5328 -0.0253 -0.1587 0.7386
0.4721 0.2671 0.2204 -0.4375 -0.2469 -0.2063
-0.0684 -0.1336 -0.1480 0.8275 0.6177 0.5655
0.3419 0.4104 0.5025 -0.1496 -0.3107 -0.4292
0.5328 -0.0253 -0.2431 -0.1587 0.7386 1.1859
G
2
0.5280 1.1338 -0.5629 -1.2096
-0.1429 -0.3954 0.9895 2.2501
0.6666 1.6438 -0.8699 -2.0708
0.2081 0.1339 0.5540 1.9633
0.5280 1.1338 2.5068 -0.5629 -1.2096 -2.6734
-0.1429 -0.3954 -1.0634 0.9895 2.2501 5.2461
0.6666 1.6438 4.1640 -0.8699 -2.0708 -5.0824
0.2081 0.1339 -0.4628 0.5540 1.9633 6.2494
3.4.3.2.1 Dead-beat Controller Design
The dead-beat controller gain matrix
1
K for the first half of the period and
2
K for the second half of the period (i.e.
1 2
( ) ( / 2) k k K kT K kT T =− =− − +
_ _
a K x z z ),
for both linear and quadratic function parametrized controllers are tabulated below.
Table 3.5: Dead-beat feedback gains for fast-sampled full-state feedback design
Linear Function
Parametrized Dead-Beat
Controller Gains
Quadratic Function Parametrized
Dead-Beat Controller Gains
K
1
8.5514 -1.0182 7.8997 1.8969
-10.4075 -0.4543 -7.5119 -2.0319
4.6928 0.1783 3.0227 2.1015
-4.6339 -1.5682 -2.6370 -2.2332
6.3056 -1.1072 6.2353 1.3940
-2.0184 -0.3108 -1.1709 -0.3729
-5.3397 -0.0756 -4.0463 -1.0377
3.5436 -0.1434 2.3516 1.5236
-0.2551 -0.5255 0.0397 -0.2434
-2.7877 -0.6649 -1.7034 -1.2680
K
2
19.0532 5.5369 18.6277 6.0953
-8.2482 -2.6136 -7.1025 -2.4154
11.9012 2.7121 7.2053 5.5732
-4.6492 -1.5528 -2.6565 -2.2155
4.5248 1.2079 4.8899 1.5569
4.5346 1.1960 4.9909 1.5750
-2.7148 -0.8085 -2.5636 -0.8496
3.0329 0.4550 1.9284 1.3941
3.1494 0.4363 2.0188 1.4437
-1.6548 -0.4221 -0.9923 -0.7763
62
3.4.3.2.2 Optimal Controller Design
The steady-state optimal controller gain matrix
1
K for the first half of the
period and
2
K for the second half of the period are obtained using the recursive
relations given before (i.e.
1 2
( ) ( / 2) k k K kT K kT T =− =− − +
_ _
a Kx z z ). Steady state
optimal controller gains for linear and quadratic function parametrizations are
tabulated below.
Table 3.6: Steady-state optimal feedback gains for fast-sampled full-state feedback
design
Linear Function Parametrized
Steady-state Optimal Controller
Gains
Quadratic Function
Parametrized Steady-state
Optimal Controller Gains
K
1
1.1842 -1.2620 2.2265 -0.0138
-1.3522 0.9248 -1.6195 0.1022
0.6674 0.3047 -0.3807 0.8430
0.6087 -0.9232 1.1407 -0.7309
1.5698 -1.4910 2.5509 -0.0084
-3.2434 2.1129 -3.2779 0.1202
1.3214 -0.8285 1.1479 -0.0339
-0.1744 1.1314 -1.3179 0.9541
3.7687 -3.9241 4.5321 -1.2293
-2.1407 1.9810 -2.2444 0.3511
K
2
0.0284 0.1994 -0.8296 -0.1158
-0.1186 -0.2649 1.5170 0.0837
-0.2763 0.2347 -0.9275 0.1239
0.7700 -0.4157 1.1679 0.1011
0.1091 -0.2501 -0.0168 0.1133
-0.2672 0.7924 -0.5510 -0.3902
0.0799 -0.3940 0.7306 0.2070
-0.6895 0.1469 0.2419 -0.4215
1.4873 -0.1769 -1.0871 1.0440
-0.2091 -0.1112 0.7579 -0.2684
Simulation studies with the nonlinear model show that all of the fast
sampling regulators do stabilize the plant for a wide range of initial conditions.
System response and the control forces as function of time for designs with fast
sampling regulator are shown Figs. 3.21 and 3.22, respectively. As can be seen in
the figures, dead-beat response is achieved for both linear and quadratic function
parametrized dead-beat designs in half of the system period, however the required
63
control effort is much more larger than the optimal designs. Fast sampling steady
state optimal regulators with both the quadratic and linear function parameterization
stabilized the trivial equilibrium point within approximately two periods compared to
approximately three periods for the steady state optimal designs before (see Fig.
3.20). The computed performance index for the steady state optimal linear function
parameterized controller with fast sampling was 0.4058 whereas the performance
index of the steady state optimal linear function parameterized controller with the
sampling period equal to the periodicity of the systems was 0.5781. Similarly, the
computed performance index for the steady state optimal quadratic function
parameterized controller with fast sampling is 0.3747 compared to a performance
index of 0.4817 achieved before.
0 1 2 3
0
0.1
0.2
0.3
time (x T sec)
θ
1
(rad)
Quadratic Dead-beat
Quadratic Optimal
Linear Dead-beat
Linear Optimal
0 1 2 3
-0.1
0
0.1
0.2
0.3
time (x T sec)
θ
2
(rad)
Quadratic Dead-beat
Quadratic Optimal
Linear Dead-beat
Linear Optimal
Figure 3.21: Closed loop state trajectories for fast-sampling controllers (T
s
=T/2)
64
0 1 2 3
-2
0
2
4
6
8
time (x T sec)
u
1
0 1 2 3
-0.2
-0.1
0
0.1
0.2
time (x T sec)
u
1
Quadratic Dead-beat
Linear Dead-beat
Quadratic Optimal
Linear Optimal
0 1 2 3
-2
0
2
4
time (x T sec)
u
2
Quadratic Dead-beat
Linear Dead-beat
0 1 2 3
-0.2
-0.1
0
0.1
time (x T sec)
u
2
Quadratic Optimal
Linear Optimal
Figure 3.22: Control forces for fast-sampling (T
s
=T/2) configuration
65
3.4.3.3 Output Feedback Regulator Designs
In general, an infinite number of output feedback controllers, depending upon
the choice of controller update period, observer update period and control
parametrization, can easily be designed with the proposed method. For illustration
purposes we will consider quadratic function parametrized output feedback
controllers with synchronous (T
s
=T), and fast sampling observers (T
s
=T/2), where
the controller update period is T.
Since the separation principle holds for both synchronous and fast-sampling
observers, the steady-state optimal and dead-beat feedback gains are same as given
before.
3.4.3.3.1 Synchronous Sampling Observer
The selection of the observer gain in Eq. (3.29) can easily be done by pole-
placement considering the estimation error equation. Assuming both angles are
measured, the observer gain placing the poles of the estimation error equation at
[0.001, 0.002, 0, 0] is:
-0.2359 -0.4308
-1.5676 -2.6017
0.3991 1.5628
-2.0335 1.5979
L
=
66
3.4.3.3.2 Fast Sampling Observer
Assuming the sampling period is half of the base period of the system (i.e.
observer update period / 2
s
T T = ) , the fast sampling observer equation is given as:
1
ˆ ˆ ˆ { }
k k k k k
L
+
= + − +
_ _ _ _
x H x y C x G a (3.49)
where
{ }
{ } 1
ˆ ˆ ˆ ( ) ( / 2)
ˆ ˆ ˆ ( / 2) ( )
T
k
T
k
kT kT T
kT T kT T +
= +
= + +
_
_
x x x
x x x
{ } ( ) ( / 2)
T
k kT kT T = +
_
a a a
1 2
{ , } diag H H = H ,
1 2
{ , } diag G G = G ,
1 2
{ , } diag C C = C , and
1 2
{ , } L diag L L =
Assuming again both angles are measured, the observer gains L
1
, and L
2
,
then can be obtained
by considering the estimation error equation which leads to
placing the roots of the following polynomials:
1 1 1
2 2 2
det( ) 0
det( ) 0
zI H L C
zI H L C
− + =
− + =
(3.50)
The observer gains which place the roots of both of the following
polynomials at [0.001, 0.002, 0,0] are as follows:
1
-0.2264 0.1251
3.2097 -1.0546
0.1850 -0.4130
-4.0468 0.5987
L
=
,
2
0.0381 0.0281
3.9675 -1.3191
-2.0407 0.3615
0.1120 -1.0377
L
=
67
Simulation studies of the closed-loop system with the linear model showed,
as expected by linear system theory, that estimation errors of both synchronous and
fast-sampling observers reach to zero independently of the control force, and that the
convergence is faster for fast-sampling observer. Because its faster convergence rate,
the compensators with the fast sampling observer were able to stabilize the non-
linear plant for a wider range of initial conditions. Moreover, they have a smoother
response compared to compensators with synchronous sampling observer.
Time histories of the angular positions and control forces for dead-beat
controller with synchronous and fast sampling observers are shown Fig. 3.23.
Similarly, angular positions, and the control forces for the steady-state optimal
controller with synchronous fast sampling observer are shown Fig. 3.24.
68
0 1 2 3 4 5 6 7 8
-0.2
-0.1
0
0.1
time (x T sec)
θ
1
(rad)
Dead-beat output feedback T
observ er
=T
Dead-beat output feedback T
observ er
=T/2
0 1 2 3 4 5 6 7 8
-0.1
0
0.1
0.2
time (x T sec)
θ
1
(rad)
Optimal output feedback T
observ er
=T
Optimal output feedback T
observ er
=T/2
0 1 2 3 4 5 6 7 8
-0.5
0
0.5
time (x T sec)
θ
2
(rad)
0 1 2 3 4 5 6 7 8
-0.4
-0.2
0
0.2
0.4
time (x T sec)
θ
2
(rad)
Dead-beat output feedback T
observ er
=T
Dead-beat output feedback T
observ er
=T/2
Optimal output feedback T
observ er
=T
Optimal output feedback T
observ er
=T/2
Figure 3.23: Closed loop state trajectories for output feedback controllers
69
0 1 2 3 4 5 6 7 8
-1
0
1
time (x T sec)
u
1
0 1 2 3 4 5 6 7 8
-1
-0.5
0
0.5
time (x T sec)
u
1
Dead-beat output feedback T
observ er
=T
Dead-beat output feedback T
observ er
=T/2
Optimal output feedback T
observ er
=T
Optimal output feedback T
observ er
=T/2
0 1 2 3 4 5 6 7 8
-0.5
0
0.5
1
1.5
time (x T sec)
u
2
0 1 2 3 4 5 6 7 8
-0.4
-0.2
0
0.2
time (x T sec)
u
2
Dead-beat output feedback T
observ er
=T
Dead-beat output feedback T
observ er
=T/2
Optimal output feedback T
observ er
=T
Optimal output feedback T
observ er
=T/2
Figure 3.24: Output Feedback Control Forces
70
Chapter 4
Robust Stability Analysis
In the preceding chapter, the proposed methodology for stability analysis and
controller synthesis problem for linear periodically time varying systems was
presented. The underlying assumption was that all the parameters in the state space
model of the physical system were known exactly. However, in practice, the plant
model is an imprecise representation of the actual plant due to parameter uncertainty,
unmodeled dynamics, operating point variations, sensor noise, and unknown
disturbance inputs.
In this chapter the linear periodically time varying system with real
parametric uncertainties of the form below will be considered.
( ) ( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( )
i i i i
i i
t t t t t t t
t t t
δ δ
=
∑ ∑
= A + A x + B + B u x
y C x
ɺ
(4.1)
where
( ) t A , ( ) t B , ( ) t C - Nominal system matrices satisfying
( ) ( ) ( ) ( ) ( ) ( ) t +T t t +T t t +T t A = A , B = B , C = C
( )
i
t A , ( )
i
t B - Perturbation matrices satisfying
( ) ( )
i i
t T t + = A A , and ( ) ( )
i i
t T t + = B B
71
1
i
δ ≤ i =1,2,..,p - Parametric uncertainties (normalized)
In the following section an extension of the discretization method given in
section 3.1 will be applied to the perturbed case.
4.1 Discrete-Time Representation of Periodic Systems with Real
Parametric Uncertainties
Consider (4.1) with control vector parameterized as:
( ) ( )
k
t t =Ξ u a for
0 0
( 1) t kT t t k T + ≤ < + + (4.2)
Applying the discretization algorithm described in chapter 3, the discrete time
representation of the system in Eq. (4.1) is:
( ) ( )
1
( ) ( )
k k k +
= + + + x H H δ x G G δ a (4.3)
where
0
( )
k
t kT = + x x , and
( )
1 2
, ,...
p
δ δ δ = δ
The system matrices above are multinomial expressions in the uncertain
parameters such that:
( ) ...
i i ij i j ijk i j k
i i j i j k
δ δδ δδ δ = + + +
∑ ∑∑ ∑∑∑
H δ H H H (4.4)
( ) ...
i i ij i j ijk i j k
i i j i j k
δ δδ δδ δ = + + +
∑ ∑∑ ∑∑∑
G δ G G G (4.5)
Truncating the series yields an r
th
order truncated approximation of the
discrete time representation of the system as follows:
( ) ( )
( ) ( )
1
( ) ( )
r r
k k k +
= + + + x H h δ x G g δ a (4.6)
where H, and G are the nominal system matrices,
( )
( )
r
h δ ,
( )
( )
r
g δ are the r
th
-order
truncated approximations of the system matrices due to the parametric uncertainties.
72
The main steps for the computational algorithm to obtain the truncated
approximations to the system matrices due to uncertainties given in Eq. (4.6) (i.e.
computation of H , G ,
i
H ,
i
G ,
ij
H
ij
G ,
iij
H ,
iij
G , etc) are presented in the
following.
1. The nominal system matrices from the initial time t
0
to the j
th
integration step
t
0
+jh (i.e.
, 1 0 0
ˆ
( , )
j j
t t jh
−
= + Φ Φ , and
, 1 0 0
ˆ
( , )
j j
t t jh
−
= + Ψ Ψ ) are computed
using the recursions below:
, 1 1, 2
ˆ ˆ
j j j j j − − −
= Φ Φ Φ ,
1,0
ˆ
n n
I
×
= Φ (4.7)
, 1 1, 2
ˆ ˆ
j j j j j j − − −
= + Ψ Φ Ψ Ψ ,
1,0
ˆ
n l ×
= Ψ 0 (4.8)
Note for j=N, the nominal system matrices H , and G are obtained as in Eq.
(3.15)
2. Denote the r
th
order truncation operation as
r
f , the uncertainty matrices at
each integration step are give as:
ˆ ˆ ˆ ˆ
j j
r
r r r r r
j j j-1 j, j-1 j-1
( ) ( ) ( ) ( ) ( )
h (δ) = Φ h (δ) + h (δ)Φ + h (δ)h (δ) (4.9)
ˆ
ˆ ˆ ˆ
j j-1 j j-1
r
r r r r r r
j j, j-1 j j
( ) ( ) ( ) ( ) ( ) ( )
g (δ) = h (δ)Ψ + Φ g (δ) + g (δ) + h (δ)g (δ) (4.10)
where the truncated approximations due to uncertainties in each integration step j are
given as :
{ }
( ) ( ) ( ) ( ) ( )
1 1, , 2 2, , 3 3, , 4 4, ,
( ) ( ) ( ) ( ) ( )
j
r r r r r
j H j H j H j H
h d d d d = + + + h δ k δ k δ k δ k δ (4.11)
{ }
( ) ( ) ( ) ( ) ( )
1 1, , 2 2, , 3 3, , 4 4, ,
( ) ( ) ( ) ( ) ( )
j
r r r r r
j G j G j G j G
h d d d d = + + + g δ k δ k δ k δ k δ (4.12)
73
In the above equations,
( )
, ,
( )
r
i j H
k δ , and
( )
, ,
( )
r
i j G
k δ are the r
th
order truncated
approximations of the Runge-Kutta parameters multiplying the state and the control
input respectively (see the detailed description of the computational method in the
appendix B).
Thus, to obtain the truncated approximations to system matrices due to the
parametric uncertainties
( )
( )
r
h δ , and
( )
( )
r
g δ defined in Eq.(4.6):
i. The truncated approximations due to uncertainties in each integration
step j are computed using Eq. (4.11) and Eq. (4.12)
ii. The operations defined in Eq. (4.9) and Eq. (4.10) are used
recursively for N integration steps
4.2 Examples
4.2.1 Active Control of Helicopter Rotor Blade Flapping Motion in Forward
Flight
The mathematical model for the linearized flapping motion of a single rigid
rotor blade around the steady state periodic equilibrium is given by (see Johnson
[32]):
2
2
2
4 4
1 sin 1 cos sin 2
8 3 8 3
d d
d d
β γ β γ
μ ψ μ ψ μ ψ β
ψ ψ
+ + + + +
=
2
8
1 sin 2 sin 2 ( )
8 3
γ
μ ψ μ ψ θ ψ
+ +
4
2 sin
8 3
γ
μ ψ λ
− +
(4.13)
where
β - is the flap angle
74
θ - is the pitch angle
t ψ =Ω - is the non-dimensional time, which is the angle of rotation of the blade
around flapping angle, and is called the azimuth angle
Ω - is the angular speed of the rotor shaft
μ - is the advance ratio (ratio of the velocity of the helicopter with respect to the
surrounding air to the velocity of the blade tip)
γ - is the ratio of the aerodynamic forces and the inertial forces on the blade and is
termed the Lock number
λ =is the air inflow ratio.
In the following λ is assumed to be a bounded disturbance acting on the
system. In more accurate models, the inflow ratio is considered as a state having its
governing equations (see Appendix-B for details).
Without an active control the input command by the pilot is transmitted to the
blade with a conventional swash plate mechanism, thus the pilot command is
composed of collective pitch and cyclic pitch such that:
1 1
( ) cos( ) sin( )
p o c s
θ ψ θ θ ψ θ ψ = + + (4.14)
And the stability of the trivial solution can be improved by an active control
input in addition to the pilot control.
( ) ( ) ( )
p c
θ ψ θ ψ θ ψ = + (4.15)
The active control input can either be imposed on the conventional swash
plate where sensors and actuators are in the non-rotating frame, or the control input
75
is directly applied to each blade by means of sensors and actuators attached to
individual blades. The first approach is called Higher Harmonic Control, and the
latter approach is called Individual Blade Control (see Calise, Wasikowski, and
Schrage [12]).
Defining the states as
T
d
x
d
β
β
ψ
=
ɺ and the control vector ( )
c
u θ ψ = ; Eq.
(4.13) is rewritten in state-space form as following:
( ) ( ) x A x B u ψ ψ = + ɺ
( ) y C x ψ = (4.16)
where
2
0 1
( )
4 4
1 cos sin 2 1 sin
8 3 8 3
Aψ
γ γ
μ ψ μ ψ μ ψ
=
− + + − +
(4.17)
2
0
( )
8
1 sin 2 sin 2
8 3
Bψ
γ
μ ψ μ ψ
=
+ +
(4.18)
1 0
( )
0 1
C ψ
=
(4.19)
In the following it is assumed that the feedback control is implemented
through a flap torque actuator located at the blade root and the controller has an
access to both flap position and velocity by means of sensors located on the blade.
For the following parameter values 8 γ = , and for the forward flight velocities
corresponding to advance ratio values 0 0.5 μ ≤ ≤ , the perturbed flapping equation
76
(4.16) is written in the form (4.1) with nominal advance ratio 0.25
n
μ = , and
normalized perturbation matrices
1
A and
1
B with respect to the maximum
perturbation 0.25
p
μ = as:
( ) ( )
1 1 1 1
( ) ( ) ( ) ( ) ( ) ( ) ( )
n n
x A A x B B u t ψ ψ δ ψ ψ ψ δ ψ = + + + ɺ (4.20)
where
2
0 1
( )
4 4
1 cos sin 2 1 sin
3 3
n
n n n
A ψ
μ ψ μ ψ μ ψ
=
− + + − +
(4.21)
2
0
( )
8
1 sin 2 sin 2
3
n
n n
B ψ
μ ψ μ ψ
=
+ +
(4.22)
1
0 0
( )
1 1
cos sin
3 3
A ψ
ψ ψ
=
− −
(4.23)
1
0
( )
2
sin
3
B ψ
ψ
=
(4.24)
With a pulse function parameterized control vector, the flapping equation can
be discretized as a first order approximation as:
( ) ( )
1 1 1 1 k k k
δ δ
+
≅ + + +
1
x H H x G G a (4.25)
where the nominal and the first order approximations for the state and the control
input vector are obtained by the point-mapping algorithm as follows:
77
0.7453 -0.0437
0.0465 0.6090
=
H and
-1.7850
-1.1290
=
G
1
-0.0006 -0.0690
0.0938 -0.0093
=
H and
1
-2.2558
-1.5005
=
G (4.26)
Note that the above choice of control vector parameterization is the most
limited case. Even when a conventional swash plate mechanism is used as an
actuator, the control vector parameterized as
1 1
( ) cos( ) sin( )
p o c s
θ ψ θ θ ψ θ ψ = + + can
be implemented without introducing additional hardware.
The table below shows the comparison between the state and input matrices
for the linear approximation and direct integration. As can be seen from the table the
approximation is close to the true values in the advance ratio values between 0.125 to
0.375, and the control input matrix is highly sensitive to the variation in the advance
ratio μ. The results show that a low gain feedback can stabilize the flapping motion
of the blades in the forward flight envelope even in this extreme flight conditions,
where the uncertainty in the advance ratio is one hundred percent. Though in this
case where only one parameter is uncertain, i.e. the worst-case variation occurs either
in the maximum or the minimum of the parameter value, one may check the robust
stability by checking the stability of extreme plants. Thus, the linear approximation
can be used in preliminary control design. Another advantage of the linear
approximation is that it gives an insight to the worst-case parameter variation
directions.
78
Table 4.1: Comparison of state and input matrices—direct integration versus linear
approximation
Advance
Ratio
H G H
(Linear
Approximation)
G
(Linear
Approximation)
μ = 0 0.6747 -0.0083
0.0083 0.6757
0.3253
-0.0083
0.7459 0.0253
-0.0473 0.6183
0.4708
0.3715
Μ = 0.125 0.6918 -0.0174
0.0171 0.6587
-0.6307
-0.4801
0.7456 -0.0092
-0.0004 0.6137
-0.6571
-0.3788
Μ =0.375 0.8416 -0.0825
0.1051 0.5314
-3.2112
-2.0457
0.7450 -0.0782
0.0934 0.6043
-2.9129
-1.8792
μ =0.5 0.9911 -0.1272
0.2037 0.4339
-4.9964
-3.3282
0.7447 -0.1127
0.1403 0.5997
-4.0408
-2.6295
The main reason for the magnitude of difference in the linear approximation
comes from neglecting the second order terms in μ. One way to have a better
approximation is to retain higher order values of μ in the analysis. An easier way of
analyzing the effects of parameter perturbations in this specific problem is to assume
the perturbations in μ
2
terms as an independent perturbation; this renders the need for
higher order approximation algorithm though being a conservative approach by
allowing perturbation directions that are not possible in the original problem,
Let
2
1 δ ≤ , then the perturbed flapping equation is written as follows:
( ) ( )
1 1 2 2 1 1 2 2
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
n n
x A A A x B B B u t ψ ψ δ ψ δ ψ ψ ψ δ ψ δ ψ = + + + + + ɺ (4.27)
where the perturbation matrices
2
A ,and
2
B are:
2
0 0
( )
1
sin( ) 0
12
A ψ
ψ
=
−
,
2
0
( )
1
sin
6
B ψ
ψ
=
79
With a pulse function parameterized control vector, the flapping equation can
be discretized as a first order approximation as:
( ) ( )
1 1 1 2 2 1 1 2 2 k k n k
δ δ δ δ
+
≅ + + + + + x H H H x G G G a (4.28)
where
2
H and
2
G are computed by the algorithm are:
2
-0.2040
-0.0090
=
G
2
0.0985 0.0014
-0.0061 -0.0800
=
H
For the extreme values of the advance ratio the above approach yields the
state and control input matrices as shown in Table 4.2. With the latter approach the
linear approximation is greatly improved for the maximum value of the advance
ratio. For the minimum value of advance ratio corresponding to the hovering flight,
the approximation becomes worse since the above approach allows perturbation in
negative direction for the μ
2
terms which is not physically possible.
Table 4.2: Comparison of state and input matrices-direct integration versus linear
approximation with assuming quadratic terms in advance ratio as an independent
perturbation
Advance
Ratio
H G H
(Linear
Approximation)
G
(Linear
Approximation)
μ = 0 0.6747 -0.0083
0.0083 0.6757
0.3253
-0.0083
0.6474 0.0239
-0.0413 0.6984
0.6749
0.3805
μ =0.5 0.9911 -0.1272
0.2037 0.4339
-4.9964
-3.3282
0.8432 -0.1113
0.1342 0.5197
-4.2449
-2.6385
80
Control Design
By checking the spectral radius for the extreme plants using the linear
approximation as a guideline, the control law is chosen such that the nominal plant
closed loop poles are located at 0.5 and 0.6.
k
= -[-0.2378 0.1507]
k
a x (4.29)
The response of the closed loop system for the nominal, maximum and
minimum advance ratios in Fig. 4.1 shows that the following error for the flap states
are reduced to virtually zero by the state feedback controller in equation (4.29) in
less than 3 rotor revolutions for all possible forward flight conditions.
The table below shows the spectral radius for the controlled system for the
worst –case parameter variations namely for hover flight (μ=0), and maximum
forward flight speed (μ=0.5).
Table 4.3: Spectral radius for rotor control system--exact versus linear
approximation
Advance
Ratio
Spectral Radius
ρ
s
Exact
Spectral Radius
ρ
s
Linear
Approximation
Approach I
Spectral Radius
ρ
s
Linear
Approximation
Approach II
μ = 0 0.7469 0.8514 0.7859
μ =0.5 0.4283 0.7456 0.5569
81
Figure 4.1: Full state feedback pulse function parameterized controller response for
rigid rotor flapping model
82
4.2.2 Double Inverted Pendulum with Load Direction Uncertainty
The linearized equation of motion for the uncontrolled double inverted
pendulum around the trivial equilibrium position θ
1
= θ
2
=0, is given in state space
(see chapter 3 section 3.4.3) as:
( ) ( ) { }
0 0 1 0
0 0 0 1
( ) ( )
0.5 ( 3) 0.5 (2 ) 0 0
0.5 (5 ) 1.5 2 0 0
t t
k p k p
k p p k α
=
− −
− − −
z z ɺ (4.30)
where
( )
1 2 2 2
( ) ( ) ( ) ( ) ( )
T
t t t t t θ θ θ θ = z
ɺ ɺ
Consider the load direction parameter is uncertain such that:
0.8 1.2 α ≤ ≤ (4.31)
In state space form, the uncertainty in the load direction parameter is then
represented as:
( )
1 1
( ) ( ) ( ) ( ) ( ) ( )
n n
t A t A t t B t t δ = + + z z u ɺ (4.32)
where the nominal plant and control input matrices are same as before, the
normalized perturbation
1
( ) /
n n
= − δ α α α ,and the perturbation matrix due to the
uncertainty is:
1 1
0 0 0 0
0 0 0 0
( ) ( )
0 0 0 0
0 0.2 0 0
A t A t
kp
π
= + =
−
(4.33)
By point mapping, the discretized nominal state and nominal control input
matrices are same as before.
83
And for this specific problem, for all values of load direction parameter the
point mapping should be symplectic, i.e. the mapping is area preserving, and the
eigenvalues come in pairs such that 1/
i j
λ λ = . Thus the error in approximation -due
to the truncation of higher order terms in the uncertain parameter-, can be found by
checking how closely the state transition matrix H, satisfies symplectic matrix
relations. First and second order approximation for the state transition matrices for
the extreme values of load direction parameter are shown below in Table 4.4. The
eigenvalues and the symplecticity of the eigenpairs of the state transition matrix by
the first order and second order approximation are shown in Table 4.5.
Table 4.4: State transition matrices with load direction uncertainty-- linear versus
second order approximation
Load
Direction
Parameter
Value
Linear Approximation
( ) H α =
1 1 n
H H δ +
Second order Approximation
( ) H α =
2
1 1 1 11 n
H H H δ δ + +
α = 0.8
0.0062 0.0745 0.0559 0.4865
0.1494 -2.3466 3.9457 -0.6781
-0.4032 0.8849 -1.9984 0.3308
-2.6138 -0.2880 0.5110 -0.3421
0.0011 0.1334 0.0478 0.4969
0.5786 -2.3364 4.2693 -0.4818
-0.4924 0.8772 -2.0580 0.2853
-2.0039 -0.8189 1.2286 -0.2774
α = 1.2
-0.2401 -0.5054 0.0636 0.1612
-1.7176 -0.2541 0.8165 -0.5739
0.0277 0.2767 -1.3246 0.2680
-2.3650 1.9593 -3.1088 0.8305
-0.2451 -0.4464 0.0554 0.1717
-1.2884 -0.2439 1.1402 -0.3776
-0.0615 0.2690 -1.3842 0.2225
-1.7551 1.4284 -2.3912 0.8952
84
Table 4.5: Error analysis for linear and second order approximation
4.3 Error Analysis for Truncated Point-Mappings
In general, the closeness of the r-th order truncated approximations to the
exact values for the transition matrices in the range of uncertain parameters can be
checked for some special cases, such as for Hamiltonian systems (as shown in the
inverted double pendulum example), and for linear time invariant systems. However,
the robust stability of the system for a given controller and/or the uncontrolled
system with the truncated approximation yields only approximate answers. In order
have a quantitative measure for the closeness of the approximations a special
algorithm is developed to obtain a (consistent) norm bound on the error introduced
by the truncation operation for the state and input transition matrices. The
development of the algorithm is explained in detail as follows.
Load
Direction
Parameter
Value
Eigenvalues of the
State Transition
Matrix
(First Order
Approximation)
Product of
eigenvalue
pairs
i j
λλ
Det(H) Eigenvalues of the
State Transition
Matrix
(Second Order
Approximation)
Product of
eigenvalue
pairs
i j
λλ
Det(H)
α = 0.8 -0.1178 + 1.1268i
-0.1178 - 1.1268i
-0.2998
-4.1455
1.2836
1.2428
1.5952
-0.0315 + 1.0473i
-0.0315 - 1.0473i
-0.2491
-4.3585
1.0978
1.0857
1.1919
α = 1.2 0.5299 + 1.1721i
0.5299 - 1.1721i
-1.3148
-0.7332
1.6546
0.9640
1.5952
0.5745 + 0.8146i
0.5745 - 0.8146i
-0.7074
-1.4196
0.9936
1.0042
0.9979
85
The discrete time representation of the system with the uncertainty structure
given in Eq. (4.1) is formally:
( ) ( )
1
( ) ( )
k k k +
= + x H δ x G δ a (4.34)
where
( )
0
( ) , ,
o
t T t =Φ + H δ δ (4.35)
( ) ( )
0
( )
( ) , ,
o
t T
t
t F d τ τ τ
+
= Φ
∫
G δ δ (4.36)
A r
th
order of approximation of the discrete time representation of the system
by the point-mapping algorithm is given in Eq. (4.2) .Let the errors introduced by the
r
th
order approximation of the state and input transition matrices due to the
uncertainties by the point mapping algorithm be defined as:
( )
( ) ( )
r r
−
( )
h H δ h δ
ɶ
≜ (4.37)
( )
( ) ( )
r r
−
( )
g G δ g δ ɶ ≜ (4.38)
Note that the error in both state and input transition matrices are due to
truncation operation and due to the numerical method used to integrate Eq. (4.35),
and Eq. (4.36). Since the norm bound on the error introduced by a consistent
numerical method of order p, is in the order of h
p
where h is the integration step size,
the approximation error due to integration method we used can be made arbitrarily
small by choosing a sufficiently small integration step. Thus, in the following
analysis, the errors introduced by the truncation operation are considered only.
86
For each integration step, the errors due to the uncertainties with r
th
order
truncation are explicitly written as:
{ }
( ) ( ) ( ) ( ) ( )
1 1, , 2 2, , 3 3, , 4 4, ,
j
r r r r r
j H j H j H j H
h d d d d = + + + h k k k k
ɶ ɶ ɶ ɶ ɶ
(4.39)
{ }
( ) ( ) ( ) ( ) ( )
1 1, , 2 2, , 3 3, , 4 4, ,
j
r r r r r
j G j G j G j G
h d d d d = + + + g k k k k
ɶ ɶ ɶ ɶ
ɶ (4.40)
By denoting ( )
( )
, 0
m
n j n m
A A t jh c h = + + , ( )
( )
, 0
m
i j i m
A A t jh c h = + + ,
( )
( )
, 0
m
n j n m
F F t jh c h = + + , ( )
( )
, 0
m
i j i m
F F t jh c h = + + ,t he error terms defined in Eq.
(4.37), and Eq. (4.38) are given as:
( )
1, ,
r
j H n n ×
= k 0
ɶ
for r >0 (4.41)
( )
1, ,
r
j G n l ×
= k 0
ɶ
for r >0 (4.42)
{ }
(2) (1)
2 , ,
( )
2, ,
for r=1
for r>1
i i j i i j
r
i i
j H
n n
c h A A δ δ
×
=
∑ ∑
k
0
ɶ
(4.43)
{ }
(2) (1)
2 , ,
( )
2, ,
for r=1
for r>1
i i j i i j
r
i i
j G
n l
c h A F δ δ
×
=
∑ ∑
k
0
ɶ
(4.44)
(3) (3) (1)
3 , , 2, ,
(3) (2) (1) (1) (2) (1)
, 3 , 2 , , , 2 ,
( )
3, ,
(3) (1)
3 , 2, ,
for r=1
n j i i j j H
i
i i j i i j n j i i j n j i i j
i i i i
r
j H
i i j j H
i
c h A A k
A c h A I c h A A A c h A
c h A k
δ
δ δ δ δ
δ
+
+ + + +
=
∑
∑ ∑ ∑ ∑
∑
k
ɶ
ɶ
ɶ
{ }
for r=2
for r>2
n n ×
0
(4.45)
87
(3) (3) (1) (3) (2) (2) (1)
3 , , 2, , 3 , , 2 , ,
( )
(3) (1) 3, ,
3 , 3, ,
for r=1
n j i i j j G i i j i i j n j i i j
i i i i
r
j G
i i j j G
i
c h A A k c h A F c hA F
c h A k
δ δ δ δ
δ
+ + +
=
∑ ∑ ∑ ∑
∑
k
ɶ
ɶ
ɶ
{ }
for r=2
for r>2
n l ×
0
(4.46)
( )
(4) (4) (1)
4 , , 3, ,
(4) (3) (2) (1) (2) (1)
, 4 3 2 , , , , 2 ,
(4) (1) ( )
4 , 3, , 4, ,
for r=1
n j i i j j H
i
i i j n j n j i i j i i j n j
i i i
r
i i j j H j H
i
c h A A k
A c h c hc hA A A A I c hA
c h A k
δ
δ δ δ
δ
+
+ + +
=
∑
∑ ∑ ∑
∑
k
ɶ
ɶ ɶ
{ }
(4) (2)
4 , 3, ,
for r=2
for r=3
i i j j H
i
n n
c h A k δ
×
∑
0
ɶ
for r>3
(4.47)
{ }
(4) (4) (1)
4 , , 3, ,
(4) (3) (2) (2) (1)
4 , 3 , , 2 , ,
(4) (1) ( )
, 4 3, , 4, ,
for r=1
n j i i j j G
i
i i j i i j n j n j n j
i i
r
i i j j G j G
i
c h A A k
c h A c h A F c hA F
A c hk
δ
δ δ
δ
+
+ +
=
∑
∑ ∑
∑
k
ɶ
ɶ
ɶ
{ }
(4) (2)
, 4 3, ,
for r=2
for r=3
i i j j G
i
n l
A c hk δ
×
∑
0
ɶ
for r>3
(4.48)
88
To explicitly define the errors for the state and input transition matrix
( ) r
h
ɶ
,
and
( ) r
g ɶ respectively, N integration steps should be carried out. To simplify the
notation the functional dependency of the r
th
order truncated state and input state
transition matrices due to the uncertainties as well as error terms are omitted in the
following. Note that the input for the whole period is a constant, i.e.
k
= a a .
Beginning with the first step of integration,
( ) ( ) ( ) ( )
1 1 1 1
( ) ( ) ( ) ( )
0 1 0 1
r r r r
x t h x t + = + + + + + Φ h h Ψ g g a
ɶ
ɶ (4.49)
( ) ( ) ( ) ( )
2 2 2 2
( ) ( ) ( ) ( )
0 2 0 2
2
r r r r
x t h x t h + = + + + + + + Φ h h Ψ g g a
ɶ
ɶ (4.50)
Substituting Eq.(4.49) in Eq.(4.50), and rearranging the terms yield:
( ) ( ) ( ) { }
1 2 2 1 2 1 2 1 1
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
0 2 ,1 2 1 2 1 0
2 ( )
r r r r r r r r r
n
x t h x t + = + + + + + + + + Φ H Φ h h Φ h h Φ h h h Φ h h
ɶ ɶ ɶ
( ) {
1 2 2 2 1
( ) ( ) ( ) ( ) ( )
2 1 2 2 1
r r r r r
+ + + + + + Φ Ψ Ψ Φ g h Ψ g h g
( ) ( ) }
2 1 2 1 1 2
( ) ( ) ( ) ( ) ( ) ( )
2 1
r r r r r r
+ + + + + + Φ h g h Ψ g g g a
ɶ
ɶ ɶ ɶ (4.51)
With the recursive operations defined in Eq. (4.7), Eq (4.8), and denoting the
r
th
order truncation operation, and complementary truncation operation as
r
f , and
r
f
−
respectively, Eq. (4.51) is equivalent to:
( )
{ 1 2 2 1 2 1
( ) ( ) ( ) ( ) ( ) ( )
0 2,1 2 1
ˆ
2
r r
r r r r r r
x t h
−
+ = + + + + Φ Φ h h Φ h h h h
( ) ( )}
2 1 2 1 1
( ) ( ) ( ) ( ) ( )
2 1 0
( )
r r r r r
x t + + + + + Φ h h h Φ h h
ɶ ɶ ɶ
+
{ 1 2 2 2 1 2 1
( ) ( ) ( ) ( ) ( ) ( ) ( )
2,1 2 1
ˆ
r r
r r r r r r r
−
+ + + + + Ψ Φ g h Ψ g h g h g
( ) ( ) }
2 1 2 1 1 2
( ) ( ) ( ) ( ) ( ) ( )
2 1
r r r r r r
+ + + + + + Φ h g h Ψ g g g a
ɶ
ɶ ɶ ɶ (4.52)
89
Thus the error terms in the second integrations step are given explicitly as:
( ) ( )
2 1 2 1 2 1 1
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
2, 2 1
r
r r r r r r r r
H
E
−
= + + + + + h h Φ h h h Φ h h
ɶ ɶ ɶ
(4.53)
( ) ( )
2 1 2 1 2 1 1 2
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
2, 2 1
r
r r r r r r r r r
G
E
−
= + + + + + + h g Φ h g h Ψ g g g
ɶ
ɶ ɶ ɶ (4.54)
Similarly, the states in the third integration step are:
( ) ( ) ( ) ( )
3 3 3 3
( ) ( ) ( ) ( )
0 3 0 3
3 2
r r r r
x t h x t h + = + + + + + + Φ h h Ψ g g a
ɶ
ɶ
( ) ( )
{ 3 3 1 2 2 1
( ) ( ) ( ) ( ) ( ) ( ) ( )
3 2,1 2 1 2, 0
ˆ
r
r r r r r r r
H
E x t
= + + + + + +
Φ h h Φ Φ h h Φ h h
ɶ
+
} 1 2 2 2 1
( ) ( ) ( ) ( ) ( ) ( )
2,1 2 1 2,
ˆ
r
r r r r r r
G
E
+ + + + +
Ψ Φ g h Ψ g h g a
( )
3 3
( ) ( )
3
r r
+ + + Ψ g g a ɶ (4.55)
Using the operations defined in Eq. (4.9) and Eq. (4.10), the error terms in the
third integration step are explicitly given as:
( ) ( )
3 2 3 3
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
3, 3 2, 2,1 2 2,
ˆ ˆ
r
r r r r r r r r
H H H
E E E
−
= + + + + + h h Φ h h Φ h
ɶ
(4.56)
( ) ( )
3 1 3 3 2 3
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
3, 3 2, 2,1 2,
ˆ
ˆ ˆ
r
r r r r r r r r r
G G G
E E E
−
= + + + + + + h g Φ h h Ψ g g
ɶ
ɶ (4.57)
Thus, by inspection the errors in the j
th
integration step are as follows:
( ) ( )
1
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
, 1, , 1 1 1,
ˆ ˆ ˆ
j j j j
r
r r r r r r r r
j H j j H j j j j H
E E E
−
−
− − − −
= + + + + + h h Φ h h Φ h
ɶ
(4.58)
( ) ( )
1 1
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
, 1, 1, 2 1,
ˆ
ˆ ˆ
j j j j j j
r
r r r r r r r r r
j G j j G j j j G
E E E
− −
−
− − − −
= + + + + + + h g Φ h h Ψ g g
ɶ
ɶ (4.59)
To compute the norm bounds on the state and input transition matrices due to
r
th
order truncated approximations, first the norms bounds on errors for each
integration step given in Eq. (4.39), and Eq. (4.40) should be computed. Any
90
consistent matrix norm can be used. Then the norm bound on the errors can be
calculated using the recursive equations (4.58) and (4.59). Based on the above, a
special algorithm is developed to compute an upper bound on the errors introduced
by truncation operation ,using maximum singular value as the norm, is developed.
In the next section, a numerical example is given to illustrate the algorithm.
4.3.1 Uncertain Underdamped Mathieu Equation
Consider the lightly damped Mathieu equation, with the amplitude of the
periodic forcing uncertain, such that:
2
2
cos(2 )
d x dx
c kx a t x
dt dt
+ + = (4.60)
where c = 0.2, k =1, 0 0.5 a ≤ ≤
With [ ]
T
x x = x ɺ , the system is given in state space as:
{ }
1 1 n
A A δ = + x x ɺ
(4.61)
where
1
1 δ ≤
0 1
1 0.25cos(2 ) 0.2
n
A
t
=
− + −
(4.62)
1
0 0
0.25cos(2 ) 0
A
t
=
(4.63)
By point-mapping algorithm, the first, second, and third order truncated
approximation in uncertainties are:
( )
1 1 1 k k
δ
+
≈ + x H H x (4.64)
( )
2
1 1 1 1 11 k k
δ δ
+
≈ + + x H H H x (4.65)
91
( )
2 3
1 1 1 1 11 1 111 k k
δ δ δ
+
≈ + + + x H H H H x (4.66)
Where the computed values for the nominal state transition matrix, and the
state transition matrices due to the uncertainty are tabulated below:
Table 4.6: Computed values of the nominal state transition matrix
and the state transition matrices due to the uncertainty
H
-0.72915618643517 0.01155943281998
-0.01155943281998 -0.73146807299917
H
1
0.04322550318409 0.14879784664539
0.14152305855385 0.01346593385501
H
11
-0.01435172045961 -0.00204257368744
-0.00201449774116 -0.01394320572213
H
111
-0.00003283971841 -0.00007446062498
0.00133862266599 -0.00001794759342
And the approximation errors due to the truncation in uncertainties in are
computed as: (using maximum singular value as the norm)
Table 4.7: Computed norm bounds on error in state transition
matrix due to truncation operation
Truncation Order
(r)
Computed Upper Norm Bound on Error
( )
( )
max
r
h σ
ɶ
1 0.04870025599542
2 0.00261925601691
3 3.852179951222627 x 10
-5
92
Note that when
1
1 δ = the amplitude of the periodic forcing is zero, thus the
system becomes a LTI system. For this case, the state transition matrix can be
computed analytically, and the error introduced by the truncation operation can be
calculated. The table below shows the errors with respect to truncation order.
Table 4.8: Exact norm bounds on error in state transition matrix due to truncation
operation
Truncation
Order
(r)
( ) r
h
ɶ
Norm Bound on Actual
Error
( )
( )
max
r
σ h
ɶ
1
-0.01437812367759 -0.00211613776281
-0.00065254374310 -0.01395489612503
0.01558581972813
2
-0.00002640321798 -0.00007356407537
0.00136195399806 -0.00001169040291
0.00136224867341
3
0.00000643650043 0.00000089654961
0.00002333133207 0.00000625719051
2.500272919958065 x 10
-5
The example clearly shows that the computational algorithm enables us to
compute an upper bound on the errors due to truncation operation. Computation of
the upper norm bound on the error can also be used to the select the truncation order
for a given set of uncertainty range.
93
4.4 Robust Stability Analysis by Structured Singular Value
Theory
The discrete time representation of the perturbed systems with an r
th
order
approximation is given as:
( ) ( )
( ) ( ) ( ) ( )
1
( ) ( )
r r r r
k k k +
= + + + + + g x H h δ h x G g δ a ɶ
ɶ
k k
= y Cx (4.67)
where
H
, G
are the nominal state and input transition matrices respectively
( )
( )
r
h δ , and
( )
( )
r
g δ are the state and input transition matrices due to the uncertainties
( ) r
h
ɶ
, and
( ) r
g ɶ are the error terms introduced in Eq. (4.37) and Eq. (4.38) respectively
From the singular value properties we know that the maximum norm of a
matrix is less than the maximum singular value. Let the compute norm bounds for
the error matrices for the state and input matrices be c
1
and c
2
respectively. Then the
family of plants described by the above equation is a subset of the following family
of plants described by:
( ) ( )
( ) ( )
1 1 1 1 2 2 1
( ) ( )
r r
k p p k p p k
c c δ δ
+ + + + +
= + + + + + x H h δ H x G g δ G a
k k
= y Cx (4.68)
Where
1
1
p
δ
+
≤ ,
2
1
p
δ
+
≤
H
p+1
is a n by n and G
p+1
is a m by l matrices with their all elements equal to one.
94
For robust stability analysis of the perturbed plant, the discrete time system
given in Eq. (4.67) with output feedback of the form
k k
K =− a y , then can be
transformed to general control configuration as shown below in Fig. 4.2, by pulling
out the uncertainties.
Figure 4.2: General Control Configuration
The 2 by 2 block system matrix P maps the discrete signals u
∆
, and
k
a
to the
signals y and y
∆
and partitioned conformally with the signals such that:
11 12
21 22 k k
P P
P P
Δ Δ
=
y u
y a
(4.69)
And the transfer function matrix from u
∆
to y
∆
is related to the P and K by a
lower linear fractional transformation (LFT) as:
( )
1
11 12 22 21
( , )
l
M F P K P P K I P K P
−
= + − ≜ (4.70)
P
K
∆
y
k
a
k
y
∆
u
∆
95
Similarly the transfer function from a
k
to y
k
is related to P and ∆ by an upper
LFT such that:
( )
1
22 21 11 12
( , )
k k
y u u
T F P P P I P P
−
= Δ + Δ − Δ ≜ (4.71)
By the small gain theorem, a sufficient condition for the robust stability of
the closed loop system is that ||∆||
∞
y u
T
Δ Δ
∞
<1. Since the uncertainties are normalized
such that ||∆||
∞
≤ 1, a sufficient condition for robust stability is given as:
1
y u
T
Δ Δ
∞
< (or ||M||
∞
< 1) (4.72)
The above conditions can be arbitrarily conservative in the presence of
structured uncertainties, and they are only necessary and sufficient when the
uncertainty has no structure.
Definition (Structured Singular Value)[17] Let M be a given complex matrix and
let ∆=diag{∆
i
} denote a set of uncertainty matrices with ||∆||
∞
<1, and with a given
block diagonal structure (i.e. some of the blocks may be repeated and some may be
restricted to be real). The structured singular value is defined by:
max
1
( )
inf{ ( ) | det( ) 0}
M
I M
μ
σ − =
Δ
Δ Δ
≜ (4.73)
And if no such structured Δ exists then ( ) 0 M μ =
Δ
For discrete-time system, the necessary and sufficient condition for robust
stability is given as:
( )
( ) 1
j
M e
θ
μ
Δ
< , [0, ] θ π ∀ ∈ (4.74)
96
4.4.1 Computation of Structured Singular Value
As defined above, the structured singular value directly provides robustness
conditions for general types of uncertainties. The problem is the computation of the
structured singular value. In fact, it has been shown that the calculation of the exact
value of μ of an uncertain model with a general uncertainty structure is a NP-hard
problem [18]. However, algorithms for computation of upper and lower bounds on
the structured singular value even in the case of general type (mixed
parametric/dynamic) uncertainty structure has been implemented in Matlab Robust
Control Toolbox [65].
In addition to computation of μ, obtaining M-∆ structure has been automated
by the new version of robust control toolbox which in general requires tedious
calculations for a complicated uncertainty structure.
The above formulation also allows us to formulate the robust control problem
as discrete μ -synthesis problem. The main advantage of the above LFT formulation
is that performance robustness can be studied easily incorporating the performance
specification as a fictitious uncertainty block to the general control formulation
shown in Fig. 4.2.
97
4.4.2 Structured Singular Value Analysis Examples
4.4.2.1 Lightly Damped Mathieu Equation
In this example, we revisit the example given in section 4.2.1. The main idea
of using this example is to investigate the conservativeness of our approach
compared to the conventional ones. And the results obtained by the approach can
readily be compared to exact (numerically obtained) values. We reformulate the
problem as finding the critical forcing amplitude for which the system becomes
unstable for the lightly damped Mathieu equation given in equation (4.60).
The simplest and the most conservative approach to find the critical value of
the forcing amplitude is by using small gain theorem. The lightly damped Mathieu
equation can be written as a feedback interconnection of the linear time invariant
dynamics of the plant with the periodic uncertain forcing terms such that:
( ) ( ) ( ) x t Ax t Bu t = + ɺ
( ) ( ) y t Cx t =
( ) ( ) z t t y =−Δ (4.75)
where
0 1
1 0.2
A
=
− −
,
0
1
B
=
( ) 1 0 C = , ( ) t a Δ ≤
98
The system described by Eq. (4.75) is then put in M −Δ connection shown
in the figure below where M (s) =C ( sI – A)
-1
B =
2
1
0.2 1 s s + +
Figure 4.3: M −Δ connection
Then by using the small-gain theorem, the system is robustly stable if
( ) ( )
max max
( ) ( ) 1 M j j σ ω σ ω Δ < (4.76)
Using Eq. (4.76), the critical value for the perturbation is found to be 0.2,
which is a restatement of the circle criterion. Note that with the above formulation of
the problem the periodicity of the perturbation term is completely ignored and is
assumed to be nonlinear time varying with sector bounds.
With the first order truncated point mapping approximation (i.e. r=1), the
equivalent discrete-time system is given by:
( )
(1)
1 1 1 k k
δ
+
= + + x H H h x
ɶ
(4.77)
where the norm bound using maximum singular value as the norm has been
computed as
( )
(1)
max
0.04870025599542 h σ ≤
ɶ
M
∆
z
y
99
From the singular value properties we know that
(1) (1)
max max
( ) ( ) h h σ ≤
ɶ ɶ
thus,
the family of uncertain matrices given in Eq. (4.77) will be a subset of the one
below:
( )
1 1 1 2 2 k k
δ δ
+
= + + x H H H x (4.78)
where
(1)
2 max
1 1
( )
1 1
σ
=
H h
ɶ
,
2
1 δ ≤
In the above simple case the system can be written as:
( )
1 k l r k
E E
+
= + Δ x H x (4.79)
where ( )
1 2 2 2
,
x
diag I δ δ Δ= ,
(1)
max
(1)
max
1 0 ( )
0 1 ( )
l
E
σ
σ
=
h
h
ɶ
ɶ
,
1
(1) (1)
max max
( ) ( )
r
E
σ σ
=
H
h h
ɶ ɶ
With the above factorization, we can represent the uncertain system in
M −Δ connection shown in Fig. 4.3 as:
1
2 2
0
k r k
k l x k
E
E
+
=
x H x
z w
(4.80)
By small μ-theorem the system is robustly stable if ( ) ( ) 1
zw
T z μ < .
where the transfer function from w to z is given as
1
( ) ( )
zw l n r
T z E zI H E
−
= − .
100
In general it is well known that computation of μ is non-polynomial, in other
words the computational complexity increases combinatorically with the number of
parameters. However practical algorithms are developed, and have been
implemented in MATLAB Robust Control Toolbox. (see Appendix-C for a brief
review of μ-theory, and computational algorithms) . Moreover, the above steps of
getting the M −Δ connection have also been mechanized with the new Robust
Control Toolbox even with systems with more general uncertainty structure.
Note that the computed upper norm bounds on the error matrices for different
truncation orders are conservative estimates of the exact norm bounds. And the
sufficient condition for robust stability for the uncertain system is given as:
( ) ( ) 1
zw
T z μ < (4.81)
In the table below, the stability bounds obtained by μ –analysis with and
without the norm bounds on error by different levels of truncation order are shown.
Table 4.9: Stability bounds by μ –analysis
Truncation
order-r
Stability bound obtained
with the addition of norm
bound on error (sufficient
condition)
Stability bound obtained
without considering the error
introduced by p-mapping
algorithm
1 0.35242220909411 0.40950465174423
2 0.39706775192473 0.40103602801025
3 0.40081126280686 0.40086974138558
101
All of the discrete time representation of the plants with different truncation
orders have been found to be not robustly stable in the range of the parameter
1
δ .
Thus, the computed norm bounds on the transition matrix error due to the truncated
approximations in the uncertainty are still an upper norm bounds on the exact error
matrices. Therefore, the system will be stable if the forcing amplitude is smaller than
the critical values as shown in the second column in the table. For comparison
purposes the critical values computed by neglecting the error terms are also shown in
the third column.
The value for the forcing amplitude for which the system remains stable is
computed via point-mapping algorithm after several iterations is found to be
0.40084. As can be seen from the table, the critical values obtained by μ –analysis
becomes less and less conservative estimate to the exact value when the truncation
order is increased.
4.4.2.2. Robust Stability Analysis of Multivariable Feedback Design
Example with Respect to Load Parameter Uncertainties
In this example, the robust stability of the steady-state optimal and dead-beat
feedback designs with linear and quadratic function parametrized controllers for the
inverted double pendulum example (see chapter 4 section 4.3) with respect to
parametric uncertainties in the load parameters will be considered.
102
The linearized equation of motion for the double inverted pendulum example
around the trivial equilibrium position
1 2
0 θ θ = = is given as:
1 1 2
0.5 ( 3) 0.5 (2 ) k p k p θ θ θ = − + −
ɺɺ
2 2
1 2
/ 2 / 2 u ml u ml + −
( ) ( ) { }
2 2
2 1 2 1 2
0.5 (5 ) 1.5 2 /2 3 /2 k p p k u ml u ml θ θ α θ = − + − − − +
ɺɺ
(4.82)
where
2
/ k k ml = ,and ( )
1 2
cos( ) / p P P t l k ω = +
By defining the state vector as
( )
1 2 1 2
( ) ( ) ( ) ( )
T
t t t t θ θ θ θ = z(t)
ɺ ɺ
, the
equations of motion for the double pendulum can be rewritten in state space form as:
( ) ( )
1
2
2
0 0 1 0
0 0
0 0 0 1 ( ) 0 0
1
( ) ( )
0.5 ( 3) 0.5 (2 ) 0 0 ( ) 1/ 2 1/ 2
0.5 (5 ) 1.5 2 0 0
1/ 2 3/ 2
u t
t t
k p k p u t ml
k p p k α
= +
− −
−
− − −
−
z z ɺ
(4.83)
Consider the constant and the periodic loading are uncertain such that the
load parameters
1 1
/ p Pl k = , and
2 2
/ p P l k = are such that:
1
1.8 2.2 p ≤ ≤
2
0.6 0.8 p ≤ ≤ (4.84)
The parametric uncertainties are then represented in state space form as
follows:
( )
1 2
( ) ( ) ( ) ( ) ( ) ( ) ( )
n 1 2 n
t t t t t t t δ δ = + + + z A A A z B u ɺ (4.85)
103
where the nominal plant and control input matrices are same as before,
1
δ , and
2
δ
are the normalized perturbation parameters with respect to constant loading p
1
and
periodic loading p
2
.
The perturbation matrices due to the uncertainties are:
0 0 0 0
0 0 0 0
0.5 0.5 0 0
0.5 0.5 0 0
t
k k
k k
=
−
−
1
A ( )
0 0 0 0
0 0 0 0
0.5 cos(2 ) 0.5 cos(2 ) 0 0
0.5 cos(2 ) 0.5 cos(2 ) 0 0
t
k t k t
k t k t
=
−
−
2
A ( )
First order approximation (see Eq. (4.86)), and second order approximation
(see Eq. (4.87)), and third order approximation (see Eq. Eq. (4.88)) for the state
transition matrices are given below:
(1)
1 2
( )
n 1 2
δ δ = + + h H H H (4.86)
(2) (1) 2 2
1 2 1 2 12 11 22
δδ δ δ = + + + h h H H H (4.87)
(3) (2) 2 2
1 2 1 2 112 221
δ δ δδ = + + h h H H
3 3
1 2 111 222
δ δ + + H H (4.88)
As before in the load direction parameter uncertainty analysis, for all values
of the load parameters, the point mapping should be symplectic, thus the error in
approximation can be found by checking how closely the state transition matrix H,
satisfies symplectic matrix relations.
104
The perturbation matrices due to uncertainties are computed by point-
mapping algorithm are shown below:
1
-0.3579 0.2240 -0.2706 -0.0064
0.2102 -0.5517 0.3546 -0.3910
-0.3548 0.5663 -0.5811 0.2984
-0.5645 0.0222 0.1358 -0.3285
=
H ,
2
0.0973 -0.0119 0.0069 0.0428
0.1870 -0.1394 0.0438 0.0065
-0.0412 0.1097 -0.0944 0.0520
0.1748 -0.0629 0.1231 0.0523
=
H
11
-0.0509 0.0192 -0.0097 -0.0165
0.0749 -0.0824 0.0352 -0.0269
-0.1290 0.1151 -0.0655 0.0241
0.0650 -0.1124 0.0700 -0.0679
=
H ,
12
0.0049 0.0168 -0.0079 0.0186
0.0192 -0.0199 0.0024 -0.0025
-0.0294 0.0516 -0.0221 0.0258
0.0447 -0.0271 0.0102 0.0071
=
H
22
-0.0038 0.0023 -0.0001 -0.0006
-0.0019 0.0011 -0.0006 -0.0001
-0.0027 0.0040 0.0001 0.0010
-0.0059 0.0006 -0.0006 -0.0028
=
H ,
112
-0.0006 0.0028 -0.0013 0.0022
0.0008 -0.0015 0.0004 -0.0007
-0.0038 0.0066 -0.0024 0.0033
0.0033 -0.0027 0.0002 0.0002
=
H
111
0.0000 -0.0034 0.0017 -0.0034
0.0045 -0.0031 0.0005 0.0004
-0.0094 0.0038 -0.0006 -0.0032
0.0129 -0.0124 0.0043 -0.0026
=
H ,
221
-0.0009 0.0008 -0.0001 0.0000
0.0000 -0.0001 0.0000 -0.0001
=
-0.0010 0.0015 -0.0003 0.0005
-0.0007 -0.0003 0.0002 -0.0007
H
222
0.0567 -0.0668 0.0180 -0.0215
0.0339 -0.0128 0.0027 0.0099
1.0e-003
-0.0731 -0.0006 0.0163 -0.0533
0.2379 -0.1526 0.0205 0.0276
= ×
H
The eigenvalues and the symplecticity of the eigenpairs of the state transition
matrix by the first order, second order and third order approximation for the extreme
values of the load parameters are shown in Table 4.10.
105
Table 4.10: Comparison of state transition matrices obtained by first, second, and
third order truncated approximations
Load
Parameter
Values
Product of
Eigenvalue
pairs
r =1
λ
i
λ
j
Det(H) Product of
Eigenvalue
pairs
r =2
λ
i
λ
j
Det(H) Product of
eigenvalue
pairs
r =3
λ
i
λ
j
Det(H)
p
1
= 1.8
p
2
= 0.6
1.2270
1.0924
1.3403
0.9788
0.9924
0.9713
1.0003
0.9995
0.9998
p
1
= 1.8
p
2
= 0.8
1.1788
0.9817
1.1572
0.9842
1.0114
0.9954
1.0002
0.9984
0.9987
p
1
= 2.2
p
2
= 0.6
1.1370
1.0178
1.1572
1.0273
0.9856
1.0125
1.0021
0.9971
0.9992
p
1
= 2.2
p
2
= 0.8
1.0716
1.2508
1.3403
1.0117
1.0294
1.0415
1.0007
0.9997
1.0004
p
1
= 1.8
p
2
= 0.7
1.1971
1.0326
1.2361
0.9831
1.0043
0.9873
1.0000
0.9986
0.9986
p
1
= 2.2
p
2
= 0.7
1.0984
1.1253
1.2361
1.0185
1.0042
1.0228
1.0013
0.9978
0.9991
p
1
= 2
p
2
= 0.6
1.0011
1.0118
1.0128
0.9998
1.0001
0.9999
1.0000
1.0000
1.0000
p
1
= 2
p
2
= 0.8
1.0009
1.0120
1.0128
1.0002
0.9999
1.0001
1.0000
1.0000
1.0000
106
The table shows that the accuracy of the mappings increases dramatically by
retaining the higher order terms as well as the linear terms in the uncertainties. With
the third order approximation, the error in approximation for the product of the
eigenvalue pairs, and the error in approximation in the determinant of the state
transition matrices are less than 0.2 percent in all the extreme cases considered
above. In addition it can be seen that the effect of change in periodic loading
parameter p
2
, has a less effect on the eigenvalues of the state transition matrix,
though the uncertainty in the periodic loading parameter is more than that of the
constant loading parameter p
1
percentage wise.
The observation above is verified with the computation of upper bounds of
the error matrix defined in Eq. (4.42) (see Table 4.11) , where the norm bound on
error matrix is drastically reduced when the truncation order is incresed.
Table 4.11: Computed Norm Bounds on State Transition Matrix Error
Truncation
Order
Computed Upper Norm Bound on Error
( )
( )
max
r
h σ
ɶ
1 14.6991
2 0.7160
3 0.0275
Similar to the truncated approximations for the state transition matrix, the
approximations for the input matrices for the linear and quadratic function
107
parameterized control vectors due to the uncertainties in the load parameters are
given as:
(1)
1 1 2 2
( )
n
δ δ = + + g G G G (4.89)
(2) (1) 2 2
1 2 12 1 11 2 22
δδ δ δ = + + + g g G G G (4.90)
(3) (2) 2 2
1 2 112 1 2 221
δ δ δδ = + + g g G G
3 3
1 111 2 222
δ δ + + G G (4.91)
Similarly, the error matrices for linear and quadratic function parameterized
controllers are shown in Table 4.12 below.
Table 4.12: Computed norm bounds on error for linear and quadratic function
parameterizations
Truncation
Order
Linear Function
Parameterized Controller
( )
( )
max
r
g σ ɶ
Quadratic Function
Parameterized Controller
( )
( )
max
r
g σ ɶ
1 2.8185 3.0257
2 0.1381 0.1443
3 0.0061 0.0063
As can be seen form the norm bounds on error, the third order approximation
is quite accurate. The upper and lower bounds on the structured singular value for the
closed loop system with nominal optimal controllers and dead-beat controllers
designed in the previous part of the paper are tabulated in the tables 4.13 and 4.14.
108
Table 4.13: Lower and Upper bounds on µ for optimal linear and quadratic function
parameterized controllers
Controller µ-Upper
Bound
µ- Lower Bound
Linear Function Parameterized Optimal
Controller
1.2517 0.8260
Quadratic Function Parameterized Optimal
Controller
0.6157 0.5587
Table 4.14: Lower and Upper bounds on µ for dead-beat linear and quadratic
function parameterized controllers
Controller µ-Upper
Bound
µ- Lower Bound
Linear Function Parameterized
Dead-beat Controller
5.43774 5.2493
Quadratic Function Parameterized Dead-beat
Controller
0.9378 0.8563
The tables 4.13, and 4.14 show that linear function parameterized steady-state
optimal controller and dead-beat controllers are not possibly robustly stable with
respect to the plant parametric uncertainties. The computation of the closed loop
eigenvalues at
1
δ =-1, and
2
δ =-1 for both of the controllers verified that the system is
not robustly stable indeed. On the other hand, the quadratic function parameterized
steady-state optimal controller and the dead-beat controller robustly stabilize the
trivial equilibrium point of the linear model for the double inverted pendulum.
The full blown nonlinear simulations asserted that the closed loop system
with the quadratic function parameterized controllers robustly stabilizes for a wide
range of initial conditions. The simulated responses of both designs for the nonlinear
109
model with an initial condition of x=[0.3 rad 0.25 rad 0 0]
T
, for different load
parameter values are shown in the figures 4.4 and 4.5. The simulations above shows
that the nominal optimal controller design is more robustly stable than the dead-beat
controller design, and this verifies the results obtained by structured singular value
analysis.
0 10 20 30 40 50
-0.2
-0.1
0
0.1
0.2
0.3
0.4
time(sec)
x 1 (ra d )
0 10 20 30 40 50
-0.4
-0.2
0
0.2
0.4
0.6
time(sec)
x 2 (ra d )
p1=2.2, p2=0.8
p1=2.2, p2=0.6
p1=1.8, p2=0.8
p1=1.8, p2=0.6
0 10 20 30 40 50
-0.4
-0.2
0
0.2
time(sec)
x 3 (ra d / s )
p1=2.2, p2=0.8
p1=2.2, p2=0.6
p1=1.8, p2=0.8
p1=1.8, p2=0.6
0 10 20 30 40 50
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
time(sec)
x 4 (ra d / 4 )
closed loop response of --quadratic function parametrized-- controller
p1=2.2, p2=0.8
p1=2.2, p2=0.6
p1=1.8, p2=0.8
p1=1.8, p2=0.6
p1=2.2, p2=0.8
p1=2.2, p2=0.6
p1=1.8, p2=0.8
p1=1.8, p2=0.6
Figure 4.4: Closed loop state trajectories with the quadratic function parameterized
dead-beat controller for different parameter values
110
0 5 10 15 20 25 30 35 40 45 50
-0.1
0
0.1
0.2
0.3
time(sec)
x 1 ( r a d )
p1=2.2, p2=0.6
p1=2.2, p2=0.8
p1=1.8, p2=0.8
p1=1.8, p2=0.6
0 5 10 15 20 25 30 35 40 45 50
-0.1
0
0.1
0.2
0.3
time(sec)
x 2 ( r a d )
p1=2.2, p2=0.6
p1=2.2, p2=0.8
p1=1.8, p2=0.8
p1=1.8, p2=0.6
0 5 10 15 20 25 30 35 40 45 50
-0.2
-0.1
0
0.1
time(sec)
x 3 ( r a d / s )
p1=2.2, p2=0.6
p1=2.2, p2=0.8
p1=1.8, p2=0.8
p1=1.8, p2=0.6
0 5 10 15 20 25 30 35 40 45 50
-0.2
-0.1
0
0.1
time(sec)
x 4 ( r a d / 4 )
p1=2.2, p2=0.6
p1=2.2, p2=0.8
p1=1.8, p2=0.8
p1=1.8, p2=0.6
Figure 4.5: Closed loop state trajectories with the quadratic function parameterized
steady-state optimal controller for different parameter values
111
Chapter 5
Engineering Application Example
Spacecraft Attitude and Momentum Unloading Control
System Design
In this chapter, the control design approach explained in detail in Chapter 3,
and the extension of the design approach for robust stability analysis described in
Chapter 4 are applied to the momentum-unloading problem of a three-axis stabilized
spacecraft where the attitude control is achieved via momentum exchange devices,
and the momentum unloading is achieved via magnetic torquers onboard. This
chapter is organized as follows. In section 5.1 momentum unloading problem and the
attitude dynamics of an earth pointing spacecraft is presented. The concept of
magnetic momentum unloading is briefly described in section 5.2. The linearized
equations for the combined attitude control loop and magnetic momentum unloading
control loop is presented in section 5.3. In the last section of the chapter, the control
design method for linear periodically time varying systems is applied to the magnetic
momentum unloading control design with earth’s magnetic field modeled as a
straight dipole model. In the second part, the design using the extended approach is
given with the magnetic field modeled as a tilted dipole.
112
5.1 Spacecraft Momentum Unloading Problem
Spacecraft attitude control while preventing the actuators (i.e. Control
Moment Gyros CMG’s or momentum wheels) saturate is called momentum-
unloading problem. The angular momentum of the actuators varies almost
periodically having a constant and an almost periodic part. The momentum
unloading is necessary to prevent momentum buildup in the actuators over the
mission.
The momentum unloading control system and the attitude control system can
be considered as two interacting loops with the attitude control torques are
disturbances to the momentum unloading control system and vice versa. The
difference in bandwidths of the attitude control system and the momentum unloading
control system justifies the separate design procedures for both control systems. It is
advantageous to generate the unloading via the satellites interaction with its orbital
environment, which does not require reaction jets for unloading. The three sources
commonly used for momentum unloading are earth’s magnetic field, satellite gravity
gradient, and atmospheric drag. In the following, the proposed methodology will be
applied to magnetic unloading, which is the most commonly used unloading scheme.
In the following, attitude dynamics of an earth pointing spacecraft and the concept of
magnetic momentum unloading is briefly described.
113
5.1.1 Spacecraft Attitude Dynamics
The combined attitude control system and momentum unloading system with
typical magnetic torquers is shown in Fig. 5.1. The nonlinear equations describing
the attitude control, and momentum unloading system is given as:
s s c m gg d
ω + × = + + + H H T T T T
ɺ
w w c
ω + × =− H H T
ɺ
(5.1)
where
s
H , and
w
H are the angular momentum vector of the spacecraft and
momentum wheels respectively,
c
T is the attitude control torque supplied by the
momentum exchange devices, and
m
T is the magnetic unloading torque, and ω is the
absolute angular velocity of the spacecraft expressed in body reference frame,
gg
T is
the gravity gradient torque, and
d
T is the disturbance torque
5.1.2 Attitude Kinematics
For the Euler angle rotation sequence of φ θ ψ → → , The kinematic relation
between the rates of change of Euler angles to the absolute angular velocity of the
body axes for a spacecraft in a circular orbit is:
( )
( )
0
0
0
1 tan sin tan cos cos sin
0 cos sin cos cos sin sin sin
0 sin / cos cos / cos sin cos cos sin sin
x
y
z
φ θ φ θ φ ω ω θ φ
θ φ φ ω ω φ ψ φ θ ψ
ψ φ θ φ θ ω ω φ ψ φ θ ψ
+
= − + +
− + − +
ɺ
ɺ
ɺ
(5.2)
where
0
ω is the orbital rate of the spacecraft.
114
The attitude dynamics of the earth pointing spacecraft is completely
described by Eq. (5.1) together with the kinematic relation given in Eq. (5.2).
Figure 5.1: Spacecraft attitude and momentum unloading control system
5.1.3 Gravity Gradient Torque
For the (3-2-1) Euler angle rotation sequence, the gravity gradient torque in
body reference coordinate frame is expressed as:
( )
2
0
0
3 0
0
z y x
z x y
y x z
a a a
a a I a
a a a
ω
−
= −
−
gg
T (5.3)
where I is the spacecrafts inertia tensor, sin
x
a θ =− , sin cos
y
a φ θ = ,and
cos cos
z
a φ θ =
Attitude
Control Law
CMGs
or
MWs
Spacecraft
Dynamics
Attitude
Sensors
Momentum
Unloading
Control Law
Magnetic
torquers
T
m
T
d
Disturbance
Torques
Magnetic
Control
Torques
θ
r
θ
a
+
-
+ +
Earth’s
Magnetic
Field
H Angular Momentum Vector
115
5.2 Magnetic Momentum Unloading
A magnetic moment generated by the electromagnets (torque coils, torque
rods) or gimbaled permanent magnets onboard the spacecraft causes a magnetic
torque, which can be applied to unload the accumulated angular momentum on the
momentum wheels or control moment gyros. The magnetic torque as expressed in
spacecraft body reference frame is:
m
T = M ×B
(5.4)
where B is the earth’s Magnetic field.
5.2.1 Earth Magnetic Field Models
The earth Magnetic field is generally defined as [60]:
V =−∇ B (5.5)
where V is earth’s geomagnetic potential and it is given as the solution of the
Laplace equation:
2
0 V ∇ = (5.6)
A solution of equation (5.6) is given as:
( )
1
1 0
( , , ) cos( ) sin( ) ( )
n
k n
m m m
n n n
n m
R
V r R a m b m L
r
θ φ φ φ θ
+
= =
= +
∑ ∑
(5.7)
where
R - Equatorial radius of earth
( )
m
n
L θ - Schmidt normalized Legendre functions
116
m
n
a ,
m
n
b - Gaussian Coefficients
, , r θ φ - Spherical Coordinates of Magnetic Coordinate Frame
Then the components of the magnetic field expressed in Magnetic
Coordinate Frame are:
( )
2
1 0
( 1) cos( ) sin( ) ( )
n
k n
m m m
r n n n
n m
R
B n a m b m L
r
φ φ θ
+
= =
= + +
∑ ∑
( ) ( )
2
1 0
cos( ) sin( ) ( )
n
k n
m m m
n n n
n m
R d
B a m b m L
r d
θ
φ φ θ
θ
+
= =
=− +
∑ ∑
(5.8)
( )
2
1 0
1
cos( ) sin( ) ( )
sin( )
n
k n
m m m
n n n
n m
R
B m a m b m L
r
φ
φ φ θ
θ
+
= =
=− +
∑ ∑
5.2.2 Straight and Tilted Dipole Approximations
The straight and tilted dipole model approximations are obtained by retaining
only the first and the first and the second terms of the Eq. (5.8) respectively. Note
that, the magnetic field components above expressed above are in magnetic
coordinate frame thus complicates the simulations necessary for attitude and
momentum unloading control system. Wheeler [61] developed a simpler approach
where the magnetic field components are expressed in orbit reference frame, and
arrived at the following model for tilted dipole model:
( ) ( ) ( ) ( ) ( ) ( ) { }
3
/ cos sin cos 1/2 sin 1 cos cos 1 cos cos
R
x e
B M r i i i ε α ε α μ α μ = − + − − − +
(5.9a)
( )[ ]
3
/ cos cos sin sin cos
R
y e
B M r i i ε ε μ = − − (5.9b)
117
( ) ( ) ( ) ( ) ( ) { }
3
/ 2cos sin sin sin 1 cos sin 1 cos sin
R
z e
B M r i i i ε α ε α μ α μ = − + + − − − +
(5.9c)
where
e
M = earth dipole moment
r = the radial distance of the spacecraft from the earth center
ε = tilt angle
i = orbit inclination
α = true anomaly + argument of perigee
e
μ α = +Δ−Ω
e
α = angle form vernal equinox east to Greenwich meridian
Δ = magnetic dipole east longitude
Ω = orbit right ascension
By setting the tilt angle to zero, we arrive at the straight dipole model in orbit
reference coordinate frame:
3
cos( )sin( )
e
Rx
M
B i
r
α
=
(5.10a)
3
cos( )
e
Ry
M
B i
r
−
=
(5.10b)
3
2sin( )sin( )
e
Rz
M
B i
r
α
=
(5.10c)
where
i - Orbit Inclination
α - True Anomaly + Argument of Perigee
118
5.3 Linearized Model
Since the goal of the attitude control system for an earth-pointing
configuration is to align the body axes with the orbit reference frame - assuming the
body axes coincide with the principal axes-, the linearized equations of motion are:
( ) ( )
2
0 0
4
x x x
x y z y z x c d m
I I I I I I T T T φ ω φ ω ψ + − + − − = + +
ɺɺ
ɺ
(5.11a)
( )
2
0
3
y y y
y x z c d m
I I I T T T θ ω θ + − = + +
ɺɺ
(5.11b)
( ) ( )
2
0 0
z z z
z y x y z x c d m
I I I I I I T T T ψ ω ψ ω φ + − + − − = + +
ɺ
ɺɺ
(5.11c)
0
x z x
w w c
h h T ω + =−
ɺ
(5.11d)
y y
w c
h T =−
ɺ
(5.11e)
0
x x z
w w c
h h T ω − =−
ɺ
(5.11f)
Note that the design of the attitude control and momentum unloading control
law is in principal a multi-variable time- varying control problem where the time
variation is due to the magnetic field. However, the momentum unloading control
system and the attitude control system can be considered as two interacting loops
with the attitude control torques are disturbances to the momentum unloading control
system and vice-versa. Then the design of the attitude control system and the
momentum unloading system can be designed separately. This separate design
approach is justified due to the difference in bandwidth of both control systems. The
119
state space representation of the linearized attitude control loop and linearized
momentum unloading loop are given below:
5.3.1 Attitude Control Loop
Defining the state as [ ]
T
z φ θ ψ φ θ ψ =
ɺ ɺ
ɺ , the attitude control loop is given
as:
2
0 0
2
0
2
0 0
0 0 0 1 0 0 0 0 0
0 0 0 0 1 0 0 0 0
0 0 0 0 0 1 0 0 0
4 / 0 0 0 0 / 1/ 0 0
0 3 / 0 0 0 0 0 1/ 0
0 0 / / 0 0 0 0 1/
x x x
y y
z z z
a I d I I
b I I
c I d I I
ω ω
ω
ω ω
= +
− −
−
− −
z z u ɺ
(5.12)
where
( )
y z
a I I = −
,
( )
x z
b I I = −
,
( )
y x
c I I = −
,
( )
y z x
d I I I = − −
5.3.2 Momentum Unloading Control Loop
Defining the state as [ ]
x x z
T
w w w
x h h h = , the momentum unloading control loop
is given as:
0
0
0 0
0 0 0
0 0
m
ω
ω
= +
−
x x T ɺ
(5.13)
Provided that the attitude control system is designed properly, the body
reference frame will be coincident with the orbit reference frame, thus allowing
120
replacement of the magnetic field expressed in body reference frame by the magnetic
field expressed in orbit reference frame. With the control vector composed of
magnetic moments created by the electro magnets onboard, and defined as
[ ]
T
x y z
M M M = u , the momentum unloading system is described in state-space as:
0
0
0
0 0
0 0 0 0
0 0
0
r r
r r
r r
z y
z x
y x
B B
B B
B B
ω
ω
−
= +
−
−
x x u ɺ
(5.14)
The control design for the attitude control loop, do not possess any challenges
since the system is linear time invariant. On the other hand the momentum unloading
system is time-varying due to the time varying nature of the magnetic field
complicating the design process.
5.4 Design Example
For the validation of the new momentum unloading approach a small
spacecraft in a circular orbit with an altitude of 500 kilometers, and an orbit
inclination of 30 degrees is considered. The principle moments of inertia of the
spacecraft are: I
x
= 100 kg-m
2
, I
y
= 150 kg-m
2
I
z
= 80 kg-m
2
. The attitude control
system is designed with a proportional /rate controller of the form
i i i
c p i d i
T K K ϕ ϕ =− − ɺ , using the linearized equations for the attitude control loop given
in Eq. (5.11). In the following, the magnetic momentum unloading control design
with earth’s magnetic field modeled as a straight dipole model is given first. And in
121
the second part, the design using the extended approach is given with the magnetic
field modeled as a tilted dipole.
5.4.1 Momentum Unloading Control System Design Using Straight Dipole
Model
With sampling four times the orbit period, and each control vector
component parameterized as:
(1) (2) (3) (4) 1 2 3 4
( ) ( ) ( ) ( ) ( ) ( )
T
i
i i i i i
u t t t t t t φ φ φ φ α α α α = =
φ α
,
1,2,3 i=
(5.15)
where
( )
1 m-1 / 4 m / 4
( )
0 else
m
T t T
t φ
≤ <
=
1, 2,3,4 m=
The equivalent discrete time dynamics for the momentum-unloading loop
described by Eq. (5.14) is:
1 1 1
, 1 , k k k
H G
+
= +
1 1
x x a
2 2 2
2, 1 2, k k k
H G
+
= + x x a
3 3 3
3, 1 3, k k k
H G
+
= + x x a
4 4 4
4, 1 4, k k k
H G
+
= + x x a
(5.16)
where
( ) ( ) ( )
,
[ ( 1) / 4 ( 1) / 4 ( 1) / 4 ]
x y z
T
m k w w w
h kT m T h kT m T h kT m T = + − + − + − x
( ) ( ) ( )
, 1
[ / 4 / 4 / 4 ]
x y z
T
m k w w w
h kT mT h kT mT h kT mT
+
= + + + x
1 1 1
[ ]
m m m m T
k
α α α = a 1, 2,3,4 m=
122
The state transition matrices, and the input transition matrices computed by
the point-mapping algorithm, are given below:
1 2 3 4
0 0 1
0 1 0
1 0 0
H H H H
= = = =
−
Table 5.1: Input transition matrices
1
G
-0.2074 -0.2821 0.2074
0.2395 0 0.1197
-0.2074 0.0599 -0.2074
2
G
-0.2074 -0.0599 0.2074
0.2395 0 - 0.1197
-0.2074 0.2821 - 0.2074
3
G
-0.2074 0.2821 0.2074
-0.2395 0 -0.1197
-0.2074 -0.0599 -0.2074
4
G
-0.2074 0.0599 0.2074
-0.2395 0 0.1197
-0.2074 -0.2821 -0.2074
With the transition matrices computed, we can design a discrete-time
controller of the form
m m
k k
K =− a x , (m=1,2,3,4) using pole placement technique.
The simulation shown below in Fig. 5.2 and Fig. 5.3 are for the combined attitude
control and momentum unloading control system response where the feedback gains
for the unloading loop are chosen as dead-beat gains. In the simulations the
disturbance torques acting on the spacecraft are assumed to be
5
1 10 [1 1 1]
−
= ×
d
T
N-m and the simulations are carried out with the full nonlinear model, and straight
dipole model is used for modeling earth’s magnetic field.
123
0 1 2 3 4 5 6 7
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
0.06
0.08
0.1
Tim e(x Orbital period) Tim e(x Orbital period) Tim e(x Orbital period) Tim e(x Orbital period)
ti m e e v o l u ti o n o f e u l e r a n g l e s ( d e g r e e s ) ti m e e v o l u ti o n o f e u l e r a n g l e s ( d e g r e e s ) ti m e e v o l u ti o n o f e u l e r a n g l e s ( d e g r e e s ) ti m e e v o l u ti o n o f e u l e r a n g l e s ( d e g r e e s )
φ φ φ φ
θ θ θ θ
ψ ψ ψ ψ
Figure 5.2: Simulated time history of the attitude angles for the combined attitude
and momentum control system design (using straight dipole model)
0 1 2 3 4 5 6 7
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
Time(x Orbital period)
Ang. momentum of RW in pitch,yaw,roll directions(Nmsec)
h
w x
h
w y
h
w z
Figure 5.3: Simulated time history of the momentum wheel angular momentum for
the combined attitude and momentum control system design (using straight dipole
model)
124
The simulations show that, the designed control system achieves an accuracy
of 0.005 degrees in attitude angles, and the angular momentum of the wheel is tightly
bounded.
5.4.2 Momentum Unloading Control System Design Using Tilted Dipole
Model
With the magnetic field modeled as a straight dipole, the magnetic
momentum unloading system is periodically time varying with the orbital period.
This allowed us to use the methodology described before. However, the magnetic
field components are not periodic with the orbital period for the tilted dipole model
due to the rotation of earth and regression of the ascending node. Thus the design
approach cannot be applied. However, the tilted dipole model can be put in a form
consisting of a nominal periodic part and a time varying part as:
1 2
( ) ( ) ( )cos ( )sin
n
B t B t B t B t μ μ = + +
(5.17)
where ( )
n
B t is the nominal input transition matrix,
1
( ) B t and
2
( ) B t are the
perturbation matrices due to the time varying terms cosμ and sinμ .
For a circular orbit
0
t α ω =
( )
0 0 e e
t t μ α α ω ω
Ω
= +Δ−Ω= + +Δ− Ω +
0
α =initial angle from vernal equinox
e
ω = earth rotation rate
125
0
Ω = initial orbit right ascension
ω
Ω
=orbit right ascension rate
And the nominal and perturbation matrices are given as:
( )( )
3
0 2sin sin cos
( ) / cos 2sin sin 0 sin cos
cos sin cos 0
n e
i i
B t M r i i
i i
α
ε α α
α
−
=
− −
( )( )
3
1
0 2cos sin sin
( ) / sin 2cos sin 0 cos cos
sin cos cos 0
e
i i
B t M r i i
i i
α
ε α α
α
−
= −
−
( )( )
3
2
0 2cos 0
( ) / sin 2cos 0 sin
0 sin 0
e
B t M r
α
ε α α
α
= −
−
Note that the perturbation matrices are also periodic in time with the orbital
period (i.e. ( ) ( )
i i
B t B t T = + ), and the perturbation parameters are bounded such that
sin 1 μ ≤ , and cos 1 μ ≤ . Denoting the perturbation parameters
1
cos δ μ = ,
2
sin δ μ = , the momentum unloading control system is then can be put in the
following form:
( )
1 1 2 2
( ) ( ) ( ) ( ) ( ) ( )
n n
t A t t B t B t B t δ δ = + + + x x u ɺ
(5.18)
The above is in the same form as the discussed form in Chapter 4 with null
plant perturbation matrices. This formulation allows us to use the discrete time
126
formulation for the perturbed system. Again with sampling four times over the orbit
period, and pulse function parameterized control vector components as described in
Eq. (5.15), the equivalent discrete time dynamics for the perturbed momentum-
unloading system is:
( ) ( )
1 1 1 1 1
, 1 1 2 1, 1 2
( , ) ( , )
k n k n k
H H G G δ δ δ δ
+
= + + +
1
x x a
( ) ( )
2 2 2 2 2
2, 1 1 2 2, 1 2
( , ) ( , )
k n k n k
H H G G δ δ δ δ
+
= + + + x x a
( ) ( )
3 3 3 3 3
3, 1 1 2 3, 1 2
( , ) ( , )
k n k n k
H H G G δ δ δ δ
+
= + + + x x a
( ) ( )
4 4 4 4 4
4, 1 1 2 4, 1 2
( , ) ( , )
k n k n k
H H G G δ δ δ δ
+
= + + + x x a
(5.19)
where
( ) ( ) ( )
,
[ ( 1) / 4 ( 1) / 4 ( 1) / 4 ]
x y z
T
m k w w w
h kT m T h kT m T h kT m T = + − + − + − x
( ) ( ) ( )
, 1
[ / 4 / 4 / 4 ]
x y z
T
m k w w w
h kT mT h kT mT h kT mT
+
= + + + x
1 1 1
[ ]
m m m m T
k
α α α = a 1, 2,3,4 m=
In this special case where there is no perturbation to the plant matrix, the
input uncertainty matrices due to the perturbations are null, and the input uncertainty
matrices are only due to terms that are linear in the uncertainties. (i.e.
1 2 1 1 2 2
( , )
i i i
G G G δ δ δ δ = + ) . The Table 5.2 shows the nominal and uncertainty
matrices due to the perturbations computed by the extended point-mapping
algorithm.
127
Table 5.2: Nominal and uncertainty matrices due to the uncertainties δ
1
and δ
2
Nominal Input
Transition Matrices
G
n
i
Perturbation matrices
due to δ
1
G
1
i
Perturbation matrices
due to δ
2
G
2
i
i =1 -0.2036 -0.2769 0.2036
0.2351 0 0.1175
-0.2036 0.0588 -0.2036
-0.0228 0.0932 0.0228
-0.0791 0 -0.0396
-0.0228 -0.0198 -0.0228
0 0.0228 0
-0.0914 0 0.0457
0 -0.1077 0
i =2 -0.2036 -0.0588 0.2036
0.2351 0 -0.1175
-0.2036 0.2769 -0.2036
-0.0228 0.0198 0.0228
-0.0791 0 0.0396
-0.0228 -0.0932 -0.0228
0 -0.1077 0
0.0914 0 0.0457
0 0.0228 0
i =3 -0.2036 0.2769 0.2036
-0.2351 0 -0.1175
-0.2036 -0.0588 -0.2036
-0.0228 -0.0932 0.0228
0.0791 0 0.0396
-0.0228 0.0198 -0.0228
0 -0.0228 0
0.0914 0 -0.0457
0 0.1077 0
i =4
-0.2036 0.0588 0.2036
-0.2351 0 0.1175
-0.2036 -0.2769 -0.2036
-0.0228 -0.0198 0.0228
0.0791 0 -0.0396
-0.0228 0.0932 -0.0228
0 0.1077 0
-0.0914 0 -0.0457
0 -0.0228 0
Then the design procedure for momentum unloading control law of the form
m m m
k k
K =− a x will be as follows:
i. Design a controller
m m m
k k
K =− a x
ii. Check the Schur stability of the family of closed loop matrices
( )
1 1 2 2
1,2,3,4
m m m m m m
c n n
H H G G G K m δ δ = − + + =
(5.20)
iii. Repeat the process until the family of plants described by Eq. (5.20) is
Schur stable.
The crucial step in the design procedure is the second step, and this is done
by using structured singular value theory explained in detail in the following.
128
The simplest factorization of the closed loop system for is given as:
( )
, 1 ,
m m m
m k n n l r m k
G K E E
+
= − + Δ x H x (5.21)
where ( )
1 3 3 2 3 3
,
x x
diag I I δ δ Δ= , [ ]
3 3 3 3
|
m
l
E I I
× ×
= ,
1
2
m m
m
r
m m
G K
E
G K
−
=
−
With the above factorization, we can then represent the uncertain systems
(for m=1,2,3,4) in M −Δ connection (see Fig. 5.4) as:
, 1 ,
, ,
2 2
0
m m
m k m k
n n r
m
m k m k
l x
H G K E
E
+
−
=
x x
z w
(5.22)
Figure 5.4: M −Δ connection
Note that, the system is robustly stable if
( )
( ) 1
m
zw
T z μ < (for m=1,2,3,4) by
small μ-theorem, where the transfer function from w to z is given as:
1
( ) ( )
m m m m m
zw l n n r
T z E zI H G K E
−
= − +
M
∆
w
k
z
k
129
Following the design procedure, the feedback gains are chosen as given in
Table 5.3, and the computed µ-bounds with these gains are tabulated in Table 5.4.
Table 5.3: Feedback Gains
K
1
0.9838 2.1883 -1.2249
-0.2406 -0.4873 -1.3921
1.9637 -0.3736 1.2574
K
2
1.4756 1.5809 -1.3200
-1.3833 1.5748 0.6926
1.6626 -1.2894 1.4395
K
3
0.9838 -2.1883 -1.2249
0.2406 -0.4873 1.3921
1.9637 0.3736 1.2574
K
4
1.4756 -1.5809 -1.3200
1.3833 1.5748 -0.6926
1.6626 1.2894 1.4395
Table 5.4: Computed µ-upper bounds
µ µ µ µ
1
µ µ µ µ
2
µ µ µ µ
3
µ µ µ µ
4
0.4139
0.4327
0.4139
0.4327
The simulation shown below in Fig. 5.5, and Fig. 5.6 are for the combined
attitude control and momentum unloading control system response where the
130
designed feedback gains for the unloading loop are given in Table 5.3. In the
simulations the disturbance torques acting on the spacecraft are again assumed to be
5
1 10 [1 1 1]
−
= ×
d
T N-m and the simulations are carried out with the full nonlinear
model, and tilted dipole model is used for modeling earth’s magnetic field.
The simulations show that, the designed control system is robustly stable due
to the variation of the magnetic field due to the regression of the ascending node and
earth’s rotation. The performance of the control system is slightly reduced, however
an accuracy of 0.01 degrees in attitude angles is achieved and the angular momentum
of the wheels are bounded between -0.3 and 0.3 N-ms.
0 2 4 6 8 10 12 14 16
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
0.025
0.03
time(x Orbital Period)
time evolution of euler angles(degrees)
φ
θ
ψ
Figure 5.5: Simulated time history of the attitude angles for the combined attitude
and momentum control system design (using tilted dipole model)
131
0 2 4 6 8 10 12 14 16
-0.5
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
time(x Orbital Period)
Ang. momentum of RW in pitch,yaw,roll directions(Nmsec)
h
w x
h
w y
h
w z
Figure 5.6: Simulated time history of the momentum wheel angular momentum for
the combined attitude and momentum control system design (using tilted dipole
model)
132
Chapter 6
Conclusion
In this dissertation, a practical approach for designing control systems for
linear periodically time varying systems is presented. The formulation is based on
parametrization of control, and obtaining an equivalent discrete-time control system
via a point-mapping computational algorithm.
The design method for linear periodically time-varying systems simplifies the
design procedure, guarantees asymptotic stability of the closed loop system. In
addition to the simplicity of the design approach, the designed controllers can
directly be implemented digitally. Moreover, the method can easily be applied to
various control configurations, as well as multi-rate configurations.
The extension of the algorithm to formulate the discrete time representation
of the perturbed system also allowed us to analyze robustness of the control design
for linear periodically time-varying systems with respect to parametric uncertainties.
A discrete-time dynamic model for a time-periodic system with real parametric
uncertainties is obtained using a truncated point-mapping algorithm. The error
bounds due to truncated point-mapping approximation are computed and a
robustness analysis problem of the system due to parametric uncertainties is
133
formulated using discrete-time structured singular value theory. This extension of the
approach also can be used for parametric stability analysis for time-periodic systems.
Simulation studies with mathematical models as well as engineering
application models, show good performance of the control system and demonstrate
the ability of the approach to assess the robustness of the closed-loop system to
parameter uncertainties.
134
References
[1] Arnold V.I., Mathematical Methods of Classical Mechanics, Springer-Verlag,
New York, 1978
[2] Arcara P., Bittanti S., “Periodic Control of Helicopter Rotors for Attenuation of
Vibrations in Forward Flight”, IEEE Transactions on Control Systems Technology,
Vol. 8, No. 6, pp. 883-894, November 2000
[3] Bamieh B, Pearson J.B., “ A general Framework for Linear Periodic Systems with
Application to H
∞
Sampled-Data Control Syatems”, IEEE Transactions on Automatic
Control, Vol. 37, pp.418-435, 1992
[4] Barmish B.R., New Tools for Robustness of Liner Systems, Macmillan, New York,
1994
[5] Bellman R.E., Dynamic Programming, Princeton Univ. Press, N.J., 1957
[6] Bernussou, Point Mapping Stability, Oxford, Pergamon, 1977
[7] Bogoliobov N.N., Mitropolsky Y.A., Asymptotic Methods in the Theory of
Nonlinear Oscillations, Hindustan Publishing Company, Delhi, 1961
[8] Bolotin V.V., The dynamic Stability of Elastic Systems, Holden Day, San
Francisco, 1964
[9] Boyd S., Yang Q., “ Structured and Simultaneous Lyapunov Functions for System
Stability Problems”, International Journal of Control, Vol.49, pp. 2215-2240
[10] Braatz R.D., Young P.M., Doyle J.C., Morari M., “Computational Complexity of
µ Calculation”, IEEE Transactions on Automatic Control, Vol. 39, pp. 1000-1002,
1994
[11] Calico R.A., Wiesel W.E., “Stabilization if Helicopter Blade Flapping”, Journal
of the American Helicopter Society, pp.59-64, October 1986
[12] Calise A.J., Wasikowski M.E., Schrage D.P., “Optimal Output Feedback for
Linear Time-Periodic Systems”, Journal of Guidance, Control, and Dynamics, Vol.
15, No. 2,pp. 416-423, March-April 1992
[13] Chen C., Linear Systems Theory, Oxford University Press, New York, 1999
[14] Chen T., Francis B., Optimal Sampled-Data Control Systems, Springer-Verlag,
London, 1995
135
[15] De Boor C., A Practical Guide to Splines, Springer-Verlag, New York, 2001
[16] Dellerud G.E., Control of Uncertain Sampled Data Systems, Birkhäuser, Boston,
1996
[17] Doyle J.C., “Analysis of Feedback Systems with Structured Uncertainty”, IEEE
Proceedings, Pt. D, 129:pp.242-250, 1982
[18] Fan M.K.H, Tits A.L., Doyle J.C., “ Robustness in the Presence of Mixed
Parametric Uncertainty and Unmodeled Dynamics”, IEEE Transaction of Automatic
Control, Vol. 36, No. 1, pp. 25-38, 1991
[19] Flashner H., A point Mapping Study of Dynamical Systems, Ph.D. Dissertation,
University of California, Berkeley, 1979
[20] Flashner H., Guttalu R.S., “Analysis of Non-linear Non-autonomous Systems by
Truncated Point Mappings”, International Journal of Non-linear Mechanics, Vol. 24,
No. 4, pp. 327-344,1989
[21] Flashner H., Hsu C.S., “ A study of Nonlinear Periodic Systems via The Point
Mapping Method”, International Journal for Numerical Methods in Engineering”,
Vol. 19, pp. 185-215, 1983
[22] Francis B.A., Khargonekar P.P., Robust Control Theory, Springer-Verlag, New
York, 1995
[23] Franklin G.F., Powell J.D., Workman M.L., Digital Control of Dynamic Systems,
Addison-Wesley, Reading MA.1998
[24] Fu M., Souza C.E., Dasgupta S., “Robustness Analysis of Linear Systems with
Time-varying Parameters”, Proceedings of 38
th
Conference on Decision and Control
Phoenix, Arizona, December 1989
[25] Guckenheimer J., Holmes P., Nonlinear Oscillations, Dynamics Systems, and
Bifurcation of Vector Fields, Springer-Verlag, New York, 1983
[26] Goldstein, H., Classical Mechanics, Addison-Wesley, Reading, 1957
[27] Hale J.K., Oscillations in Nonlinear Systems, McGraw-Hill, New York, 1963
[28] Harris C.J., Valenca J.M.E., The Stability of Input-Output Dynamical Systems,
Academic Press, London, 1983
136
[29] Hinrichsen D., Son N.K., “ The Complex Stability Radius of Discrete-Time
Systems and Symplectic Pencils”, Proceedings of 28
th
Conference on Decision and
Control Tampa, Florida, December 1989
[30] Horisberger, H.P, Belanger P.R., “Regulators for Linear Time Invariant Plants
with Uncertain Parameters”, IEEE Transactions on Automatic Control, Vol. AC-21,
pp 705-708,1976
[31] Hohenemser, K.H., Yin S.K., “Some Applications of the Method of Multiblade
Coordinates”, Journal of American Helicopter Society, Vol.17,No.3 , pp. 3-12, 1972
[32] Johnson W., Helicopter Theory, Dover, 1980
[33] Jönsson U.T., Kao C., Megretski A., “Analysis of Periodically Forced Uncertain
Feedback Systems”, IEEE Transactions on Circuits and Systems—I: Fundamental
Theory and Applications, Vol. 50, No.2, pp. 244-258,February 2003
[34] Jönsson U.T., Rantzer A.,” Systems with Uncertain Parameters-Time-Variations
with Bounded Derivatives”, Proceedings of the 33
rd
Conference on Decision and
Control Lake Buena Vista, Florida, 1994
[35] Jury E.I., Inners and Stability of Dynamic Systems, John Wiley and Sons, New
York, 1974
[36] Khalil H.K., Nonlinear Systems, Prentice Hall, Upper Saddle River, New Jersey,
2002
[37] Kautsky, J. and N.K. Nichols, “Robust Pole Assignment in Linear State
Feedback”, International Journal of Control, Vol. 41, pp. 1129-1155, 1985
[38] Kwakernaak, H., Sivan R., Linear Optimal Control Systems, Wiley Interscience,
New York, 1972,
[39] Mansour M., Balemi S., Truöl W., Robustness of Dynamic Systems with
Parameter Uncertainties, Birkhäuser, Berlin, 1992
[40] Megretski A., “Frequency-Domain Criteria of Robust Stability for Slowly Time-
Varying Systems”, IEEE Transactions on Automatic Control, Vol. 40, No.1, pp. 153-
155, January 1995
[41] McKillip R.M., “Periodic Model-Following for the Control-Configured
Helicopter”, Journal of the American Helicopter Society, pp. 4-12, July 1991
137
[42] Narendra K.S., Taylor J.H., Frequency Domain Criteria for Absolute Stability,
Academic Press, New York, 1973
[43] Nayfeh A.H., Mook D.T., Nonlinear Oscillations, John Wiley and Sons, New
York, 1979
[44] Newman S., The Foundations of Helicopter Flight, John Wiley and Sons, New
York, 1994
[45] Nijmeijer H., Van Der Schaft A.J., Nonlinear Dynamical Control Systems,
Springer-Verlag, New York, 1990
[46]Qui L., Davison E.J., “The Stability Robustness Determination of State Space
Models with Real Unstructured Perturbations”, Mathematical Control Signals and
Systems, Vol. 4, pp. 247-267, 1991
[47] Pandiyan R., Sinha S.C., “Time Varying Controller Synthesis for Nonlinear
Systems Subjected to Periodic Parametric Loadings”, Journal of Vibration and
Control, Vol.7, pp. 73-90, 2001
[48] Sanchez-Pena R.S., Sznaier M., Robust Systems Theory And Applications, John
Wiley and Sons, New York, 1998
[49] Skogestad S., Postlethwaite I., Multivariable Feedback Control, John Wiley and
Sons, Chichester, 1996
[50] Safonov M.G., “Stability Margins for Diagonally Perturbed Multivariable
Feedback Systems”, IEEE Proceedings, Pt. D, 129:251-256, 1982
[51] Safonov M.G., Limebeer D.J.N., Chiang R.Y., “simplifying the H ∞ Theory via
Loop-shifting, Matrix-pencil and Descriptor Concepts”, International Journal of
Control, Vol. 50 , No. 6, pp. 2467-2488, 1989
[52] Sinha S.C., Joseph P., “Control of General Dynamic Systems with Periodically
Varying Parameters via Lyapunov Floquet Transformation”, Journal of Dynamic
Systems, Measurement, and Control, Vol. 116, pp. 650-658, 1994
[53] Sinha S.C., Wu D.H., “An efficient Computational Scheme for the Analysis of
Periodic Systems”, Journal of Sound and Vibration, Vol. 151, pp.91-117, 1991
[54] Sinha S.C., Butcher E.A., “Symbolic Computation of Fundamental Solution
Matrices for Time-Periodic Dynamical Systems”, Journal of Sound and Vibration,
Vol. 206(1), pp.61-85, 1997
138
[55] Sundareshan M.K., Thathachar M.A.L, “ L
2
Stability of Linear Time Varying
Systems – Conditions Involving Noncausal Multipliers”, IEEE Transactions on
Automatic Control, Vol. AC-17, No.4, pp.504-510, August 1972
[56] Thompson J.M.T., Stewart H.B., Nonlinear Dynamics and Chaos, John Wiley
and Sons, Chichester, 1986
[57] Tondl A., Ruijgrok T., Verhulst F., Nabergoj R., Autoparametric Resonance in
Mechanical Systems, Cambridge University Press, 2000
[58] Urabe M, Nonlinear Autonomous Oscillations, Academic Press, New York, 1967
[59] Webb S.G., Calico R.A., Wiesel W.E., “Time-Periodic Control of a Multiblade
Helicopter”, Journal of Guidance, Vol. 14, No.6, pp. 1301-1308, November-
December 1991
[60] Wertz J.R., Spacecraft Attitude Determination and Control, D. Reidel Publishing
Company, Dordrecht, 1978
[61] Wheeler P.C., Magnetic Attitude Control of Rigid Axially Symmetric Spinning
Satellites in Circular Earth Orbits, SUDAER No.224, Department of Aeronautics and
Astronautics, Stanford University, April 1965
[62] Willems J.C., “ On the Asymptotic Stability of the Null Solution of Linear
Differential Equations with Periodic Coefficients”, IEEE Transactions on Automatic
Control, Vol. AC-13, No.1, pp. 65-72, February 1968
[63] Yakubovich V.A., Starzhinski V.M., Linear Differential Equations with Periodic
Coefficients, John Wiley and Sons, New York, 1975
[64] Young P.M., Doyle J.C., “Computation of µ with Real and Complex
Uncertainties”, 29
th
IEEE Conference in Decision and Control, pp. 1230-1235, 1990
[65] Young P.M., Newlin M., Doyle J.C., “ Practical Computation of the Mixed µ
Problem”, Proceedings of the American Control Conference, Chicago, pp. 2190-
2194
[66] Zames G., Kallman R.R., “ On Spectral Mappings, Higher Order Circle Criteria,
and Periodically Varying Systems”, IEEE Transactions on Automatic Control, pp.
649-652, December 1970
139
APPENDIX-A
Derivation for Point Mapping Algorithm
Consider the linear periodically time varying system
( ) ( ) ( ) ( ) ( ) t t t t t = + x A x B u ɺ (A.1)
With the parameterized control vector ( ) ( )
k
u t t a =Ξ the discrete time
representation of the system by the point-mapping algorithm is given by:
1 k k k +
= + x Hx Ga (A.2)
Integrating Eq. (A.1) with classical R-K method with a fixed integration step
size / h T N = , the state in each step j is given by:
4
0 0
1
( ) ( ( 1) )
j
m m
m
t jh t j h d k
=
+ = + − +
∑
x x (A.3)
The Runge-Kutta vectors are given as follows:
( ) ( ) ( )
1, 0 0 0 j
t jh t jh t jh = + + + + k A x F a (A.4)
where ( ) ( ) ( ) t t t = F B Ξ , and a is constant over the integration period.
( ) ( ) ( )
2, 0 2 0 2 1, 0 2 j j
t jh c h t jh c h t jh c h = + + + + + + +
k A x k F a (A.5)
Denoting ( ) ( )
( ) ( )
0 0
,
m m
j m j m
A t jh c h t jh c h = + + = + + A F F ,and substituting Eq.
(A.5) in Eq. (A.4), and collecting the terms in ( )
0
t jh + x & a yields:
{ }
( )
{ }
(2) (1) (2) (1) (2)
2, 2 0 2 j j j j j j
I c h t jh c h = + + + +
k A A x A F F a (A.6)
140
Similarly,
( )
(3) (3)
3, 0 3 2, j j j j
t jh c h = + + +
k A x k F a (A.7)
( )
(4) (4)
4, 0 4 3, j j j j
t jh c h = + + +
k A x k F a (A.8)
By inspection of the above Runge-Kutta coefficients, one can compute the
state transition matrix and in each step as:
4
( )
1
j
m
j m
m
I h d k
=
= +
∑
Φ (A.9)
where
( )
( ) ( ) ( 1)
j j j
m m m
m
k I c
−
= + A A
Similarly, the nominal input matrices in each step can be expressed as:
4
( )
1
j
m
j m
m
h d k
=
=
∑
Ψ (A.10)
where
( ) ( ) ( ) ( 1)
j j j j
m m m m
m
k c
−
= + F A F
To obtain the point-mapping, N integration steps should be carried out.
Beginning with first step of integration,
( ) ( ) ( ) ( )
0 1 0 1
t h t + = + x Φ x Ψ a (A.11)
( ) ( ) ( ) ( )
0 2 0 2
2 t h t h + = + + x Φ x Ψ a (A.12)
Substituting (A.11) in (A.12) yields
( ) ( ) ( ) ( ) ( )
0 2 0 2 2
2 t h t h + = + + + x Φ x Ψ a Ψ a ( ) ( )
2 1 0
t = Φ Φ x ( )
2 1 2
+ + Φ Ψ Ψ a (A.13)
141
From the above expression, the recursive nature of the expression can be
seen. Let us define the following recursive operations:
, 1 1, 2
ˆ ˆ
j j j j j − − −
= Φ Φ Φ (A.14)
, 1 1, 2
ˆ ˆ
j j j j j j − − −
= + Ψ Φ Ψ Ψ (A.15)
The operations above calculate the nominal system matrices from t
0
to t
0
+jh.
(i.e.
, 1 0 0
ˆ
( , )
j j
t t jh
−
= + Φ Φ , and
, 1 0 0
ˆ
( , )
j j
t t jh
−
= + Ψ Ψ ) with the initial conditions
1,0
ˆ
n n
I
×
= Φ , and
, 1
ˆ
j j n l − ×
= Ψ 0 .
Then the state transition matrix H, and input transition matrix G are
computed using the recursive operations defined in Eq. (A.14) and Eq. (A.15)
respectively with j=N, yielding
1 1
1 1
1
...
( )
N N
i N
N j i
i j N
−
+ −
= =
=
= +
∑ ∏
H Φ Φ Φ
G Ψ Φ Ψ
(A.16)
142
APPENDIX-B
Derivation for Truncated Point Mapping Algorithm
Let the structure of the linear periodically time varying system with
parametric uncertainties has a form:
( ) ( ) ( ) ( ) ( ) ( ) ( )
i i i i
i i
t t t t t t t δ δ
= + + +
∑ ∑
x A A x B B u ɺ (B.1)
With the parameterized control vector ( ) ( )
k
u t t a =Ξ , a r
th
truncated
approximation of the discrete time representation of the system by the point-mapping
algorithm is given by:
( ) ( )
( ) ( )
1
( ) ( )
r r
k k k +
= + + + x H h δ x G g δ a (B.2)
where
( )
( )
r
h δ ,
( )
( )
r
g δ are the r
th
order truncated approximations of the state and input
transition matrices uncertainties
Integrating Eq. (B.1) with classical R-K method with a fixed integration step
size / h T N = , the state in each step j is given by:
4
0 0
1
( ) ( ( 1) )
j
m m
m
t jh t j h d k
=
+ = + − +
∑
x x (B.3)
The Runge-Kutta vectors are given as follows:
( ) ( ) ( ) ( ) ( )
1, 0 0 0 0 0 j i i i i
i i
t jh t jh t jh t jh t jh δ δ
= + + + + + + + +
∑ ∑
k A A x F F a (B.4)
where ( ) ( ) ( ) t t t = F B Ξ , and ( ) ( ) ( )
i i
t t t = F B Ξ
143
( ) ( ) ( )
( ) ( )
2, 0 2 0 2 0 2 1,
0 2 0 2
j i i j
i
i i
i
t jh c h t jh c h t jh c h
t jh c h t jh c h
δ
δ
= + + + + + + +
+ + + + + +
∑
∑
k A A x k
F F a
(B.5)
Denoting ( ) ( ) ( )
( ) ( ) ( )
, 0 , 0 , 0
, ,
m m m
n j m i j i m n j m
A t jh c h t jh c h t jh c h = + + = + + = + + A A A F F
( )
( )
, 0
m
i j i m
t jh c h = + + F F , and substituting Eq. (B.5) in Eq. (B.4), and collecting the
terms in ( )
0
t jh + x & a yields:
( )
(2) (2) (1) (1)
2, , , 2 , , 0
(2) (2) (1) (1) (2) (2)
, , 2 , , , ,
j n j i i j n j i i j
i i
n j i i j n j i i j n j i i j
i i i
I c h t jh
c h
δ δ
δ δ δ
= + + + +
+ + + + +
∑ ∑
∑ ∑ ∑
k A A A A x
A A F F F F a
(B.6)
Similarly,
( )
(3) (3) (3) (3)
3, , , 0 3 2, , , j n j i i j j n j i i j
i i
t jh c h δ δ
= + + + + +
∑ ∑
k A A x k F F a (B.7)
( )
(4) (4) (4) (4)
4, , , 0 4 3, , , j n j i i j j n j i i j
i i
t jh c h δ δ
= + + + + +
∑ ∑
k A A x k F F a (B.8)
By inspection of the above Runge-Kutta coefficients, one can compute the
nominal state transition matrix in each step as: and the nominal input matrices
recursively using:
,
4
( )
1
n j
m
j m
m
I h d k
=
= +
∑
Φ (B.9)
where
( )
, , ,
( ) ( ) ( 1)
n j n j n j
m m m
m
k I c
−
= + A A
144
Similarly, the nominal input matrices in each step can be computed
recursively using:
,
4
( )
1
n j
m
j m
m
h d k
=
=
∑
Ψ (B.10)
where
, , , ,
( ) ( ) ( ) ( 1)
n j n j n j n j
m m m m
m
k c
−
= + F A F
To compute the truncated approximations to the state and input transition
matrices due to uncertainties in each integration step, multiplication and truncation
operations on Eqs (B.5-B.8) should be carried out. Denoting the r
th
–order truncated
approximations in the uncertain parameters multiplying the state and the control
input as
( )
, ,
r
i j H
k , and
( )
, ,
r
i j G
k respectively, the state and the input state transition matrices
due to the uncertainties in each integration step are computed using
{ }
( ) ( ) ( ) ( ) ( )
1 1, , 2 2, , 3 3, , 4 4, ,
( ) ( ) ( ) ( ) ( )
j
r r r r r
j H j H j H j H
h d d d d = + + + h δ k δ k δ k δ k δ (B.11)
{ }
( ) ( ) ( ) ( ) ( )
1 1, , 2 2, , 3 3, , 4 4, ,
( ) ( ) ( ) ( ) ( )
j
r r r r r
j G j G j G j G
h d d d d = + + + g δ k δ k δ k δ k δ (B.12)
To obtain the point-mapping, N integration steps should be carried out.
Beginning with first step of integration,
( ) ( ) ( ) ( )
1 1
( ) ( )
0 1 0 1
r r
t h t + = + + + x Φ h x Ψ g a (B13)
( ) ( ) ( ) ( )
2
( ) ( )
0 2 2 0 2
2
r r
t h t h + = + + + + x Φ h x Ψ g a (B14)
145
Substituting (B.14) in (B.13) yields
( ) ( ) ( ) ( ) ( )
2 2
( ) ( ) ( )
0 2 2 0 2 2
2
r r r
t h t h + = + + + + + + x Φ h x Ψ g a Ψ g a
( ) ( )
1 2 2 1
( ) ( ) ( ) ( )
2 1 2 1 0
r r r r
t = + + + Φ Φ Φ h h Φ h h x
( )
2 1 2 2 1
( ) ( ) ( ) ( ) ( )
2 1 2 1 2
r r r r r
+ + + + + + Φ Ψ Ψ h Ψ Φ g g h g a (B.15)
From the above expression, the recursive nature of the expression can be
seen. Let us define the following recursive operations
, 1 1, 2
ˆ ˆ
j j j j j − − −
= Φ Φ Φ (B.16)
, 1 1, 2
ˆ ˆ
j j j j j j − − −
= + Ψ Φ Ψ Ψ (B.17)
The operations above calculate the nominal system matrices from t
0
to t
0
+jh.
(i.e.
, 1 0 0
ˆ
( , )
j j
t t jh
−
= + Φ Φ , and
, 1 0 0
ˆ
( , )
j j
t t jh
−
= + Ψ Ψ ) with the initial conditions
1,0
ˆ
n n
I
×
= Φ , and
, 1
ˆ
j j n l − ×
= Ψ 0 .
Denoting the r
th
order truncation operation as
r
f , we can similarly express
the uncertainty matrices at each integration step as a recursive operation as shown
below:
( ) ( ) ( ) ( ) ( )
1 , 1 1
ˆ ˆ ˆ ˆ
( ) ( ) ( ) ( ) ( )
j j
r
r r r r r
j j j j j j − − −
= + + h δ Φ h δ h δ Φ h δ h δ (B.18)
1 1
( ) ( ) ( ) ( ) ( ) ( )
, 1
ˆ
ˆ ˆ ˆ ( ) ( ) ( ) ( ) ( ) ( )
j j j j
r
r r r r r r
j j j j j
− −
−
= + + + g δ h δ Ψ Φ g δ g δ h δ g δ (B.19)
Then the nominal and state transition matrices given in Eq. (B.2) are
computed using the recursive operations defined in Eq. (B.9) and Eq. (B.10)
146
respectively with j=N. And the r
th
order truncated mappings for the state and the
input transition matrices given in Eq. (B.2) are computed using the recursive
operations defined in Eq. (B.18) and Eq. (B.19) respectively with j=N.
147
APPENDIX-C
Helicopter Rotor Dynamics
We consider a rotorcraft configuration that consists of single main rotor and
single tail rotor. In this configuration the torque is balanced by the moment generated
around the helicopter’s center of gravity by the tail rotor thrust. The helicopter’s
vertical, lateral and longitudinal motion is controlled by changing the main rotor
thrust magnitude and direction, and the orientation of the helicopter is controlled by
the pitching and rolling moments around the center of gravity generated by the main
rotor thrust/direction control and yawing moment generated by the tail rotor.
Figure C.1: Single main rotor/ Tail rotor configuration (courtesy of Westland
Helicopters)
148
A mathematical model of the helicopter dynamics is a complex interaction of
aerodynamic forces and structural dynamics of the fuselage and rotors. The
periodicity of the model stems from the periodic aerodynamic loads on the rotor
blades. This can be seen by the airflow velocity experienced by the blades around the
azimuth in forward flight as shown in Fig. C.2. The rotor disc plane is divided into
two the “advancing side”(0
o
≤Ψ≤180
o
) and the “retreating side” (180
o
≤Ψ≤360
o
). On
the advancing side the airflow velocity experienced by the blades is greater than the
airflow velocity experienced by the blade on the retreating side. Thus, if the angle of
attack for the blades is kept constant, the lift generated on the advancing side will be
greater than on the retreating side creating a rolling moment around the center of
gravity of the helicopter. The rolling torque imbalance is solved by connecting the
blades to the rotor hub with hinges, allowing the blades to move freely out of the
plane of rotation, called “flapping”. Although, this hinged connection relieves the
rolling effect on the fuselage in forward flight (see Fig. C.3 for the terminology used
for rotate blade orientation angles), introduces a backward rotor plane tilt, thus
producing a thrust force component opposite to the flight direction. The tilt also
causes a nose-up pitching moment on the fuselage. To oppose the tendency of the
backward disk tilt, the pitch angles of the blades are altered periodically so as to
drive the blades flap up at the rear of the rotor disk and flap down at the front of the
rotor disk. This is achieved by introducing a “cyclic pitch”, namely by increasing the
blade pitch at the retreating side and decreasing it at the advancing side. In addition
to the cyclic pitch, constant pitch commands to the blades, termed “collective pitch”
149
are introduced to change the thrust force magnitude. The cyclic and collective pitch
controls are applied to the blades by a swash plate mechanism in most of the
helicopters (see Fig. C.4)
Figure C.2: Rotor in forward flight
The flapping motion also introduces loading in tangential (in-plane) direction
due to the Coriolis forces. To avoid loading in the tangential direction the blades are
connected to the rotor hub with a lag hinge. The lag hinges are usually accompanied
with mechanical lag dampers that mitigate the ground resonance instability.
Ψ
0
o
90
o
180
o
270
o
Ω
Flight
Direction
150
Figure C.3 shows the blade degrees of freedom and Fig. C.4 illustrates the
basic arrangement of fully articulated rotors with flap, lag and pitch hinges. There
are also new designs without lag and flap hinges (“hingeless rotors”) where the blade
motion occurs through bending at the blade root and teetering (“gimbaled”) rotors
without a lag hinge.
Figure C.3: Rotor blade motion (θ – blade pitch angle; β – blade flap angle; ζ –
blade lag angle)
Figure C.4: Helicopter rotor control with swash plate mechanism
β
ζ
θ
151
Simplified Model for Blade Flapping Motion in Forward Flight
The following assumptions are made (see Johnson [32]):
• A rotor with no hinge offset, no hinge spring,
• There is no coupling between pitch and flap motions (the only motion is rigid
flapping),
• The effects of reverse flow, blade tip loss and root cut out are neglected,
• The blade weight is neglected.
Aerodynamic Forces
Under the assumptions above, the aerodynamic forces on the blade sections
can be found by applying Blade Element Theory. The theory assumes that each blade
section is a two dimensional airfoil and the effect of the rotor wake is considered to
be purely an induced velocity at the blade section. Consider the air velocity seen by
the blade section as shown in Fig. C.5.The components of the velocity of the air with
respect to the blade are V
x
–tangential to the disk plane; V
y
–perpendicular to rotor
disk plane, and V
r
–radial to the disk plane. Note there is also a radial component of
aerodynamic force (dR) in radial direction.
Figure C.5: Aerodynamic forces on the blade section
dL
dD
φ
θ
dT
dH
α V
y
V
x
152
Figure C. 6: Velocity components
From the above figures the velocity components are given as:
sin
x
V r V =Ω + Ψ
cos
r
V V = Ψ (C.1)
tan cos sin
y
d
V V v r V
dt
β
α β = + + + Ψ
where
ν is the velocity of air downward (downwash) through the disk
r is the radial distance of the blade section to the rotor hub
R is the blade radius
Ψ
Ω
V
r
V
x
x
y
Forward
Velocity
V
r
β
z
Ω
V
r
R
153
The usual convention is to non-dimensionalize the velocity components
above by the blade tip speed R Ω . Thus, defining the advance ratio and downwash
ratio denoted by μ and
i
λ such that
V
R
μ =
Ω
,
i
R
ν
λ =
Ω
The above velocity components can be written in a non-dimensional form:
x
u sin
r
R
μ = + Ψ
r
u = cos μ Ψ (C.2)
y
u tan cos sin
i
r d
R dt
β
μ α λ μ β = + + + Ψ
Ω
The aerodynamic lift and drag forces can be expressed in terms of the lift and
the drag coefficients as:
2
1
2
L
dL V cdrC ρ =
2
1
2
D
dD V cdrC ρ = (C.3)
where c is the section blade chord and
L
C and
D
C are the lift and the drag coefficient
respectively.
Then the normal, in-plane and radial forces are:
cos sin dT dL dD ϕ ϕ = −
sin cos dH dL dD ϕ ϕ = + (C.4)
sin
radial
dR dT dD β =− +
where
radial
dD is the radial drag force due to the radial flow along the blade.
154
Blade Flapping Motion
The flapping motion of the blade now can be derived. Since the weight of the
blade is neglected, the forces acting on the blade are only the aerodynamic forces and
the centrifugal forces (see Fig. C.7). Taking moments around the flap hinge yields:
2
0 0
( sin )
R R
b
I r r mdr rdTdr β β =− Ω +
∫ ∫
ii
(C.5)
where
2
R
b
o
I mr dr =
∫
is moment of inertia of the blade around the flap hinge
m is mass per unit length of the blade
Using dimensionless quantities and dimensionless time variable t ψ =Ω , Eq.
(C.5) can be written as:
1 2
2
0
( )
d dT
r dr
d ac
β
β γ + =
Ψ
∫
(C.6)
where γ , the Lock number, is defined as:
4
b
acR
I
ρ
γ =
155
Assuming
• The angles β, θ, and φ are small
• A constant lift curve slope,
L
C aα =
• Constant chord length
• Uniform downwash ratio over the rotor disk
The aerodynamic force perpendicular to the blade section ( dT ) can be
approximated as:
2 2
( )
2 2
y
x x
x
u
ac ac
dT u u
u
α θ ≅ = − (C.7)
Using (C.7), the integral on the right hand side of the equation (C.6) –the
aerodynamic flap moment—can be evaluated as:
1
0
dT
r dr
ac
∫
= ( )
2
2
1 1
sin sin sin
8 3 4 6 4
i
μ μ μ
ψ ψ θ μα λ ψ
+ + − + +
+
1 1
sin cos sin
6 6 6 4
d
d
μ β μ
ψ μ ψ ψ β
+ − +
Ψ
(C.8)
β
Ω
Aerodynamic
Forcing ( dT )
Centrifugal
Force r
Figure C.7: Forces acting on the blade section
156
Thus, the equation of motion for rigid flapping is:
2
2
1 1
sin 1 cos sin
6 6 6 4
d d
d d
β μ β μ
γ ψ γ μ ψ ψ β
+ + + + + =
Ψ Ψ
2
2
1
sin sin
8 3 4
μ μ
ψ ψ θ
+ +
( )
1
sin
6 4
i
μ
ψ μα λ
− + +
(C.9)
where
0 1 1
cos sin
c s
θ θ θ ψ θ ψ = + + is the collective pitch and cyclic pitch control.
157
APPENDIX-D
Structured Singular Value Theory
In this appendix, a short review of canonical robust control problem leading
to closely related H
∞,
and structured singular value theory is given. The results for
both continuous and discrete-time systems are presented.
Canonical Robust Control Problem
For systems of general form as shown in the figure d.1 below, the canonical
robust control problem is mathematically described as follows. Given a multi-input
multi-output (MIMO) plant transfer function matrix P, find a stabilizing controller K
such that the closed loop transfer function matrix remains stable for all n diagonal
uncertainty matrices ∆
i
with ||∆
i
||
∞
≤ 1. Each ∆
i
for i=1,…,n-1 represents a specific
type of uncertainty either complex or real. Note that the n
th
uncertainty block is not
an actual uncertain gain. The fictitious uncertainty of size ||∆
n
||
∞
≤ 1 is wrapped as a
purely imaginary loop form the disturbance inputs d to regulated error output signals
e to represent the performance specification. This form as shown in Fig. D.1 is useful
for controller synthesis.
158
Figure D.1. Canonical robust control configuration
If the controller is given a priori, and we want to analyze the performance of
the uncertain system then the N∆-structure shown in Fig. D.2 is used.
And with the generalized plant P partitioned as:
1 11 12 1
2 21 22 2
y P P u
y P P u
=
(D.1)
The transfer function matrix from u
1
to y
1
- is related to the P and K by a
lower linear fractional transformation (LFT) as:
( )
1
11 12 22 21
( , )
l
N F P K P P K I P K P
−
= + − ≜ (D.2)
P
1
2
1
0 ... 0
: :
: : . :
0 0 ...
n−
Δ
Δ
Δ
∆
n
K
e d
y
2
u
2
u
1
y
1
159
≡
1
2
1
0 ... 0
: :
diag{ }
: : . :
0 0 ...
i
n−
Δ
Δ
Δ= Δ =
Δ
1
2
0 ... 0
: :
ˆ
: : . :
0 0 ...
n
Δ
Δ
Δ=
Δ
Figure D.2: N∆-structure
Similarly, the transfer function matrix N can be partitioned as:
11 12
21 22
N N y u
N N d e
Δ Δ
=
(D.3)
Then the closed loop transfer function from d to e is related to N and ∆ by an
upper LFT such that:
( )
1
22 21 11 12
( , )
ed u
T F N N N I N N
−
= Δ + Δ − Δ ≜ (D.4)
For robust stability analysis, the system is rearranged to M∆- structure as
shown in Fig. D.3, where M=N
11.
N
∆
∆
n
N
ˆ
Δ
160
Figure D.3: M∆-structure
By the small gain theorem, the feedback system is stable if its loop gain is
less than one. H
∞
norm is the relevant measure of the gain for linear time-invariant
systems.
Thus, the sufficient condition for the robust stability of the closed loop
system is that ||∆||
∞
y u
T
Δ Δ
∞
<1. Since the uncertainties are normalized such that
||∆||
∞
≤ 1, a sufficient condition for robust stability is given as:
1
y u
T
Δ Δ
∞
< (or ||M||
∞
< 1) (D.5)
Similarly a sufficient condition for robust performance is by small gain
theorem
ˆ
1
y u
T
Δ Δ
∞ ∞
Δ < , and is equivalent to:
1 1
1
y u
T
∞
< (||N||
∞
< 1) (D.6)
The above conditions can be arbitrarily conservative in the presence of
structured uncertainties as shown in Fig. D.1, and they are only necessary and
sufficient when he uncertainty has no structure (i.e. full block complex uncertainty).
M
∆
161
Definition (Structured Singular Value) [17]: Let M be a given complex matrix
and let ∆=diag{∆
i
} denote a set of uncertainty matrices with ||∆||
∞
<1, and with a
given block diagonal structure (i.e. some of the blocks may be repeated and some
may be restricted to be real). The structured singular value is defined by:
max
1
( )
inf{ ( ) | det( ) 0}
M
I M
μ
σ − =
Δ
Δ Δ
≜ (D.7)
And if no such structured Δ exists then ( ) 0 M μ =
Δ
From the definition of the structured singular value, a necessary and
sufficient condition for robust stability is:
( ) ( ) 1 M j μ ω
Δ
< , ω ∀ ∈ℝ (D.8)
Due to the equivalence of robust performance, and robust stability as
depicted in Fig. D.2, the necessary and sufficient condition for robust performance is
given as:
( )
ˆ
( ) 1 N j μ ω <
Δ
, ω ∀ ∈ℝ (D.9)
Similarly, if the underlying linear time invariant system is a discrete-time
system, the necessary and sufficient conditions for robust stability and robust
performance are given respectively as:
( )
( ) 1
j
M e
θ
μ
Δ
< , [0, ] θ π ∀ ∈ (D.10)
( )
ˆ
( ) 1
j
N e
θ
μ
Δ
< , [0, ] θ π ∀ ∈ (D.11)
162
Computation of Structured Singular Value
As defined above, the structured singular value directly provides robustness
conditions for general types of uncertainties. The problem is the computation of the
structured singular value. In fact, it has been shown that the calculation of the exact
value of μ of an uncertain model with a general parametric uncertainty structure is a
NP-hard problem [10]. Because of this constraint on the computation of μ, research
was directed in two directions. The first research direction looks for an exact or
approximate with guaranteed bounds polynomial time analysis of parametric
uncertainty structures, and the second by studying approximate methods that can
bound in polynomial time. The former first research direction the formulation is
based on in terms of the characteristic polynomial (see Barmish [4]), whereas the
latter is based on the structured singular value in a LFT setup [14]. The following
section describes the research efforts aimed to construct computable bounds on the
structured singular value in the case of general type (mixed parametric/dynamic)
uncertainty structure.
Mixed (Real/complex) μ
For a general uncertainty structure which has repeated scalar real blocks,
repeated scalar complex blocks, and full complex blocks such that:
1 1
1 1 1 1
1
{ ,..., , ,..., , ,..., }
r r p r r c c n c c
p p n n
r k k r k k c k k c k k m
diag I I I I δ δ δ δ
× × × ×
= Δ Δ Δ (D.12)
163
where the repeated real and complex parameters are
i
r
δ ∈ℝ i =1,…,p (
i
r
k -times
repeated),
i
c
δ ∈ℂ i =1,…,n (
i
c
k -times repeated) and the dynamic uncertainties
are full block complex uncertainties such that
d d
i i
k k
i
×
Δ ∈ℂ
The set of matrices Q , D ,G ,
D ,
G are defined next as follows:
[-1 , 1] , 1,..,
1, 1,...,
, 1,...,
i
j j
k k
i p
j p p n
I k p n p n m
δ
δ δ
∗
∗
∈ =
Δ∈ = = + +
Δ Δ = = + + + +
Q Δ ≜ (D.13)
1 1
diag{ ,..., , ,..., }
d d d d
m m m m
p n k k m k k
D D d I d I
+ × ×
D≜ (D.14)
where
0
i i
D D
∗
= > ,
r r
i i
k k
i
D
×
∈ℂ for 1,.., i p =
0
i i
D D
∗
= > ,
c c
i i
k k
i
D
×
∈ℂ for 1,..., i p p n = + +
0
j
d > ,
j
d ∈ℝ for 1,..., j m =
1 1 1 1
1
diag{ ,..., ,0 ,...0 ,0 ,...,0 }
c c c c d d d d
n n m m
p k k k k k k k k
G G
× × × ×
G≜ (D.15)
where
r r
i i
k k
i i
G G
×
∗
= ∈ℂ for 1,.., i p =
1 1
diag{ ,..., , ,..., }
d d d d
m m m m
p n k k m k k
D D d I d I
+ × ×
D≜ (D.16)
where
det( ) 0
i
D ≠ ,
r r
i i
k k
i
D
×
∈ℂ for 1,.., i p =
164
det( ) 0
i
D ≠ ,
c c
i i
k k
i
D
×
∈ℂ for 1,..., i p p n = + +
0
j
d ≠ ,
j
d ∈ℝ for 1,..., j m =
1 1 1 1
1
diag{g ,..., ,0 ,...0 ,0 ,...,0 }
c c c c d d d d
n n m m
P k k k k k k k k
g
× × × ×
G≜ (D.17)
where
g
i
∈ℝ for 1,.., i p =
1
i
i p
r
i
P k
=
=
=
∑
In addition to the defined matrices above, define the real spectral radius as
( ) max{ , is a real eigenvalue of }
R i i
A A ρ λ λ ≜ , the following bounds are obtained:
( ) ( ) ( )
1
sup inf
R
D
Q
MQ M DMD ρ μ
−
∈
∈
≤ ≤
Δ
D
Q
(D.18)
Actually the lower bound is exact (i.e. ( ) ( ) sup
R
Q
MQ M ρ μ
∈
=
Δ
Q
) as proven by
Young and Doyle [64]. However the lower bound computation is a nonconvex
optimization problem even with pure dynamic uncertainty.
For computation of the upper bound the following theorem due to Fan, Tits,
Doyle [18] is shown below.
165
Theorem D.1: Let
( ) ( ) { }
inf inf | 0
D
G
M DM j GM M G D α α α
∗ ∗
∈
∈
= ∈ + − − ≤
D
G
ℝ
then if 0 α ≤ , ( ) 0 M μ =
Δ
, otherwise ( ) M μ α ≤
Δ
.
The upper bound computation is a convex optimization problem since it is a
linear matrix inequality in both D , and G .
For computation purposes the following upper bound characterization is
given by Young [65].
Theorem D.2: If there is β > 0, and there exists matrices
D∈ D , and
G∈G such
that
( )
( )
1/ 4 1/ 4
1 2 2
max
1
1 I G DM D jG I G σ
β
− −
−
+ − + ≤
, then ( ) M μ β ≤
Δ
.
From the above theorem, the minimization problem for computing μ,
alternatively can be carried out as computation of the minimal β such that the above
inequality holds. The advantage of this alternate form is that the minimization of the
maximum singular value norm provides numerical advantages. Solving the above
minimization problem provides an initial guess for the upper bound computation
given in theorem D.1. This has been implemented in Robust Control Toolbox.
166
Robust Controller Synthesis
By the definition of μ, it follows that a necessary and sufficient condition for
the LTI continuous system to be stable is:
1 1
max
( ) 1
y u
T j σ μ ω
Δ
Δ <
for ω ∀ . (D.19)
Informally, the above condition is known as ‘”small μ theorem”. From the
condition above, if the controller K(s) stabilizes the general plant P, and the
resultant closed loop system transfer matrix satisfies the condition
1 1
( ) 1
y u
T j μ ω
Δ
<
for all frequencies then the canonical robust control problem is solved.
Mathematically, the above leads to the μ-synthesis optimization:
1 1 ,
( ) stabilizing
inf sup ( )
y u
optimum
K s
T j
ω
μ μ ω
Δ Δ
≜ (D.20)
As explained before, the computation of μ is a nonconvex optimization and is
generally a NP-hard problem. The solution to μ- synthesis problem defined above is
generally even harder and there exists no exact solution. However, there exists
algorithms such as DK -iteration which iteratively solves scaled H
∞
synthesis
problems, and they yield good sub-optimal solutions to the μ- synthesis problem. In
the order this appendix to serve as a complete reference, H
∞
optimal control problem,
and its solution in terms of H
∞
Riccati equations is given in the following section.
Then the
DK -iteration algorithm is explained briefly.
167
H
∞
-Controller Synthesis
Since solving μ- synthesis optimization is generally difficult, one can start
with a conservative singular value upper bound on μ. This leads to the H
∞
optimal
control problem which is mathematically defined as:
1 1 1 1
max, max
( ) stabilizing ( ) stabilizing
inf sup ( ) inf ( )
optimum y u y u
K s K s
T j T j
ω
σ σ ω ω
∞
=
≜ (D.21)
This problem is also closely related to the standard H
∞
control problem which
is to find a stabilizing controller such that
1 1
( ) 1
y u
T jω
∞
< .
H
∞
Controller Riccati Equations
Suppose the plant given in figure D.1. has a state space realization as:
1 2
1 1 11 12 1
2 2 21 22 2
x A B B x
y C D D u
y C D D u
=
ɺ
(D.22)
To simplify the exposure suppose further that
11
D =0,
22
D =0,
12 12
T
D D I = ,
12 1
0
T
D C = ,
21 21
T
D D I = ,
1 21
0
T
B D =
The assumptions above, by a suitable choice of variables can be made to hold
in nearly all cases as shown in [51].
168
The H
∞
controller which solves the standard H
∞
control problem is given in
observer reconstructed state-feedback form as:
ɵ
1
2 2
( )
T
u B X I YX x
−
=− −
ɵ
( )
ɵ ɵ
( ) 1 1 2 2 2 2 2
T T
x A YC C x B u YC y C x = + + + −
i
(D.23)
where X ,and Y are positive definite the solutions to the following algebraic Riccati
equations:
( )
2 2 1 1 1 1
0
T T T T
XA A X X B B B B X C C + − − + = (D.24)
( )
2 2 1 1 1 1
0
T T T T
AY YA Y C C C C Y B B + − − + = (D.25)
The solution to the H
∞
optimal control problem stated in Eq. (D.21) is
computed by the technique called γ-iteration. The idea is to scale the plant matrix P
by multiplying its first row by a positive scalar γ, and then to test whether a solution
to the standard H
∞
control problem exist. The greatest value of γ for which the
solution exists yields the minimal value
max,
1
optimum
σ
γ
= .
μ- Synthesis with DK -iteration
The DK -iteration algorithm is a combination of H
∞
-synthesis and μ-
analysis. The algorithm proceeds as follows:
• Initialize with D(s)=I, K(s)=0, and
1 1
11
( ) ( )
y u
T s P s =
169
Step 1: Compute the solution to the H
∞
control problem K(s) for the scaled
plant
1
11 12
21 22
0 0
( )
0 0
P P D D
P s
P P I I
−
′ =
Step 2: Compute the diagonal scaling ( ) D jω which minimizes
( )
1 1
1
max y u
DT D σ
−
Step 3: Find a polynomial transfer function D(s) that has a frequency
response close to the frequency response of ( ) D jω computed in step
2, using a curve fitting method, and go to step 1.
The above algorithm works well in many practical cases where all the
uncertainties are complex, albeit the fact that there is no theoretical justification that
the algorithm converges to a joint minimum in K(s) and D(s). There are also
algorithms as often called as DGK-iteration which handles real uncertainties with
less conservatism.
Abstract (if available)
Abstract
In this dissertation, a practical control design method for linear periodically time varying systems is presented. The method is based on parametrization of control, and obtaining an equivalent discrete-time control system via a point-mapping computational algorithm. The design method simplifies the design procedure, and guarantees asymptotic stability of the closed loop system. In addition, the designed controllers can directly be implemented digitally. Moreover, it is shown that the method can easily be applied to various control configurations, as well as multi-rate configurations. The design method is then extended for analysis of robustness of the control design with respect to plant parametric uncertainties. This is achieved by computation of approximate discrete time dynamics of the perturbed system by truncated point-mappings. By computing an upper norm bound on the error due truncated approximations, the robustness analysis of the system due to parametric uncertainties is formulated as a discrete time structured singular value problem. Simulation studies with mathematical models as well as engineering application models, show good performance of the control system design method and demonstrate the ability of the extension of the method to assess the robustness of the closed-loop system with respect to parametric uncertainties.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Periodic solutions of flexible systems under discontinuous control
PDF
High-accuracy adaptive vibrational control for uncertain systems
PDF
Adaptive control with aerospace applications
PDF
Control of spacecraft with flexible structures using pulse-modulated thrusters
PDF
Distributed adaptive control with application to heating, ventilation and air-conditioning systems
PDF
Nonlinear control of flexible rotating system with varying velocity
PDF
Data-driven H∞ loop-shaping controller design and stability of switched nonlinear feedback systems with average time-variation rate
PDF
Modeling and dynamic analysis of coupled structure-moving subsystem problem
PDF
Dynamic analysis and control of one-dimensional distributed parameter systems
PDF
Dynamic modeling and simulations of rigid-flexible coupled systems using quaternion dynamics
PDF
On the synthesis of dynamics and control for complex multibody systems
PDF
On the synthesis of controls for general nonlinear constrained mechanical systems
PDF
Practical adaptive control for systems with flexible modes, disturbances, and time delays
PDF
Verification, learning and control in cyber-physical systems
PDF
Time-varying closed-loop modeling of circulatory control in sleep-disordered breathing
PDF
Medium and high-frequency vibration analysis of flexible distributed parameter systems
PDF
Integrated control of traffic flow
PDF
Optimal design, nonlinear analysis and shape control of deployable mesh reflectors
PDF
Form finding and shape control of space deployable truss structures
PDF
Coordinated freeway and arterial traffic flow control
Asset Metadata
Creator
Kalender, Serkan
(author)
Core Title
Robust control of periodically time-varying systems
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Mechanical Engineering
Publication Date
09/17/2007
Defense Date
04/09/2007
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
control vector parametrization,OAI-PMH Harvest,parametric uncertainties,point-mapping,robust control,sampled-data Systems,structured singular value theory,time-periodic systems
Language
English
Advisor
Flashner, Henryk (
committee chair
), Ioannou, Petros A. (
committee member
), Yang, Bingen (
committee member
)
Creator Email
kalender@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m822
Unique identifier
UC1210950
Identifier
etd-Kalender-20070917 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-554308 (legacy record id),usctheses-m822 (legacy record id)
Legacy Identifier
etd-Kalender-20070917.pdf
Dmrecord
554308
Document Type
Dissertation
Rights
Kalender, Serkan
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
control vector parametrization
parametric uncertainties
point-mapping
robust control
sampled-data Systems
structured singular value theory
time-periodic systems