Citation
An investigation of controlling large space structures using adaptive independent modal control techniques

Material Information

Title:
An investigation of controlling large space structures using adaptive independent modal control techniques
Creator:
Hughes, Brian David
Place of Publication:
Denver, CO
Publisher:
University of Colorado Denver
Publication Date:
Language:
English
Physical Description:
179 leaves : illustrations ; 29 cm

Subjects

Subjects / Keywords:
Large space structures (Astronautics) ( lcsh )
Large space structures (Astronautics) ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaf 102).
Thesis:
Submitted in partial fulfillment of the requirements for the degree, Master of Science, Department of Electrical Engineering, Department of Computer Science and Engineering.
Statement of Responsibility:
by Brian David Hughes.

Record Information

Source Institution:
|University of Colorado Denver
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
20868362 ( OCLC )
ocm20868362
Classification:
LD1190.E54 1989m .H83 ( lcc )

Downloads

This item has the following downloads:


Full Text
AN INVESTIGATION OP CONTROLLING LARGE SPACE STRUCTURES
USING ADAPTIVE INDEPENDENT MODAL CONTROL TECHNIQUES
by
Brian David Hughes
B.GeoE., University of Minnesota, 1982
B.EE., Auburn University, 1985
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillment
of the requirements for the degree of
Masters of Science
Department of Electrical Engineering and Computer Science
1989


This thesis for the Master of Science degree by-
Brian David Hughes
as been approved for the
Department of
Electrical Engineering and Computer Science
by
Sanaa A. Azim
Date_


Hughes, Brian David (M.S., Electrical Engineering)
An Investigation of Controlling Large Space Structures Using Adaptive
Independent Modal Techniques
Thesis directed by Assistant Professor Sanaa Abdel-Azim
and reviewed by Professor Edward T. Wall
In this paper the basic principles in Large Space Structure
(LSS) modeling, estimation, and control are introduced. Several
different control techniques (both adaptive and non-adaptive) are
investigated and simulated on a computer. The concept of modal
control is explored and its advantages and disadvantages are found.
New ideas for using adaptive control on both a deterministic and
stochastic systems with modal state feedback are presented.
The form and content of this abstract are approved. I recommend its
publication.
Signed
Sanaa Abdel-Azim
Signed


CONTENTS
Chapter
I. INTRODUCTION ......................-...................... 1
II. SYSTEM MODELING ............................................ 3
2.1 Modal Theory ........................................ 10
2.2 Modal Superpostion Form ............................. 15
2.3 Computer Simulation ................................. 16
III. DIGITAL FILTERING ......................................... 20
IV. MODAL CONTROL ............................................. 29
4.1 Effects of Modeling Errors, Residual Modes, and
Spillover ........................................... 37
4.2 Distributed Actuator Control ........................ 39
4*5 Passive Control ..................................... 40
4*4 Using The Classical Techniques of Root Locus ........ 41
45 Findings from Direct Velocity Feedback
Simulations and Independent Modal Velocity
Feedback Simulations ................................ 56
V. INTRODUCTION TO ADAPTIVE PREDICTION AND CONTROL ........... 61
5.1 DARMA Modeling ...................................... 64
5.2 Adaptive Prediction ................................- 65
VI. D-STEP-AHEAD CONTROL ...................................... 72
VII. POLE-PLACEMENT ADAPTIVE CONTROL ........................... 75
7.1 Compensator Design for Pole Placement ..... 77
7.2 Compensator Design with a Self Tuning Adaptive
Controller .......................................... 85


VIII. AN ADAPTIVE CONTROL SCHEME POR INDEPENDENT MODAL CONTROL 88
8.1 Adaptive Control With Stochastic Noise .............. 91
8.2 Kalman Filtering .................................. 92
IX. CONCLUSIONS ....................................... 100
REFERENCES ....................................................... 102
APPENDIX
A. COMPUTER PROGRAMS USED FOR DESIGN ........................ 103
B. COMPUTER SIMULATION PROGRAMS ............................ 131
C. EXAMPLE OF A RESTRICTED COMPLEXITY ESTIMATOR ............. 178


TABLES
Table
1 Z-Plane Modal Magnitudes Versus Gain


FIGURES
Figure
1. A Large Space Structure ..................................... 4
2. A Spring and Damper Model ................................... 4
3* Model for a LSS ............................................. 6
4. Equivalent Spring Diagram ................................... 7
5. Formulas for Lumped Systems .............-.................. 9
6. Computer Simulation in the Physical Matrix Form ............ 18
7* Computer Simulation in the Modal Matrix Form ............... 19
8. Low Pass Filter Simulation ................................. 23
9* Band Pass Filter Simulation ................................ 26
10. Comparison of Filter Frequency Response .................... 28
11. 1st Modal Damped System Simulation ......................... 31
12. 1st and 2nd Modal Damped System Simulation ................. 32
13* Proportional Tip Velocity Feedback ........................ 41
14. Root-Locus Plot Using Direct Proportional
Velocity Feedback .......................................... 44
15* Tip Velocity Feedback with Proportional Gain ............... 46
16. Modal Tip Velocity Feedback ............................... 47
17. Root-Locus with 1st Mode Filtered ........................ 49
18. Simulation with 1st Mode Dampened ........................ 51
19* Independent Modal Control for the 1st and 2nd Modes ........ 52
20. Diagram for the 2nd Mode ........................ 53
21. Simulation with 1st and 2nd Modes Dampened ............... 55


22. Control Input for Direct Velocity Feedback ....... 57
25 Control Input for Modal Velocity Feedback .................. 58
24- Adaptive Predictor Simulation .............................. 70
25* Pole Placement Control Plant ............................... 76
26. Independent Modal Control Using a Compensator .............. 79
21. Simulation of Compensator Modal Control
With Gain Equal to 4.5 ..................................... 80
28. Simulation of Compensator Modal Control
With Gain Equal to 10 .................................... 81
29. Compensator Independent Modal Control for Modes 1 and 2 ... 83
30. Simulation for Independent Modal Control
for the 1st and 2nd Modes .................................. 84
51 Adaptive Self Tuning Regulator Compensator Control Scheme 85
32. Adaptive STR Control Scheme ................................ 88
33- Adaptive STR Control with a Kalman Filter ........ 94
34. Noise Probability Density Function ......................... 96
33- Self Tuning Parameter Estimator with Stochastic Noise .... 96
36. Simulation Results of an Adaptive Predictor
with Kalman Filtering....................................... 98


CHAPTER I
INTRODUCTION
This Thesis can he divided into two parts. The first part
(Chapters 2, 3 and 4) investigates the basics of modal control and
Large Space Structure (LSS) modeling. The second part (Chapters 5,
6, 7, and 8) explores adaptive estimation and control strategies for
both deterministic and stochastic systems.
In the first part of the thesis some of the principles of
modal control are discussed. This is not an easy task; modal control
theory covers a plethora of areas: system mathematical modeling,
vibration mechanics, structural mechanics, dynamics, modern control,
numerical analysis, estimation and optimization theory, filtering,
etc [l].
Therefore, several different practices and principles need to
be explored. First, the basics of system modeling are briefly
examined. Next, some of the necessary mathematical manipulations on
matrices are reviewed. Then some of the background in digital
filtering is given. This paper will use digital control; digital
filtering will be used as a way to separate the different modes.
A technique to optimize the gains for different modes when
using independent modal feedback control is presented. This
technique uses the classical Root Locus method combined with modern
numerical analysis. It is shown that this technique works quite well
for finding a suitable or optimal gain.


2
Finally, a number of computer simulations are made to verify
the system model and prove the control system results. Both the real
physical model (i.e., real within the framework of modeling theory)
and an ideal modal model are simulated. Much of the work in this
thesis was in developing computer programs. The programs were used
to simulate the system, to make numerical calculations (such as
finding roots of polynomials), and to design the control system.
The second part of the thesis explores different techniques
to adaptively self tune and control a large space structure. First,
the fundamentals of prediction and adaptive prediction are reviewed.
Then these fundamentals are carried forward and used to generate self
tuning adaptive control systems. A couple of different adaptive
control schemes are presented.
The last part of the thesis deals with stochastic control
using adaptive techniques. In a real space system, it would be
expected that there are many areas where noise becomes a factor in
controlling a large space structure. Therefore, this thesis will
investigate ways to reduce the effects of noise disturbances in
estimating and controlling the vibrational dynamics of the system.


CHAPTER II
SYSTEM MODELING
In modeling theory, a physical system can be represented by
using a mathematical model. The mathematical model is the set of
algebraic, vector, and differential equations which describe the
behavior of the system when various inputs and design parameters are
varied. In creating the model for this thesis, the principle of
"lumping" is used. This principle will model a distributed mass
system into individual nodal points. Each node contains a fictional
"point mass" which is approximately equal to one half of the mass per
unit length times the length between the different nodes.
The elastic properties of the physical system can be modeled
into equivalent spring and dampers (i.e., the distance between the
nodes can be modeled with massless springs and dampers). The problem
is then reduced to that of a simple 2nd order mass-spring-damper
model. In modeling theory, just about any linear physical system can
be represented by a system of equivalent mass-spring-dampers. Often
times with very accurate results.
In lumped system models, like finite element models, the
basic theory is to represent the physical system by dividing it into
a number of small segments. Each segment is mathematically modeled
by a single lumped mass and the spring and damper equivalent. A
basic rule is the more segments (or nodes) the model is divided into,
the more closely the model represents the true system. There are


4
limitations however: each node modeled into the system creates
another order in the model. The computations become excessive with
higher order models (especially if the computer used is not a main
frame). For this thesis, the system model is lumped into 4 nodes.
Figure 1 illustrates a such physical system.
Figure 1 A LSS
Assume a series of trusses are attached to a large space
structure (LSS); the LSS.is assumed as rigid for modeling purposes.
For this model the requirement that only the very tip of the truss
will contain the sensors and control actuators is imposed. Some of
the reasons for this requirement could be weight constraints,
manufacturing constraints (i.e., the truss is assembled in space unit
by unit from the end), etc. A very simplified model of the structure
using springs, masses and dampers is illustrated in Figure 2.


5
A large space structure will be of lightweight material such
as aluminum. Aluminum has a modules of elasticity E = 2.5 x 1011 pa
(N/m^). Aluminum structures typically have very low inherent damping
(around 0.2$) so the damping coefficients ... C4 can be assumed to
be equal to zero for simplification purposes. Note that in making
the assumption that Cn is zero, we are actually modeling the worst
case scenario where no inherent damping is present. Therefore, a
control design with this assumption will be slightly conservative.
By assuming this undamped case, the model is thus reduced
into a simple mass-spring equation which is easily transformed into
its mathematical equivalent using basic 2nd order differential
equations. These differential equations can be arranged into the
matrix form below.
[m]X(t) + [k]X(t) = [b]f(t) (2.1)
Where [m] is the mass matrix, [k] is the stiffness matrix,
[b] is the force distribution matrix, is the displacement at
time t, and f(t) is the force at time t. Equation (2.1 ) will be
referred to throughout this report as the physical form of the model.
By defining the state vector V = [ X X ]^, equation (2.1)
can be transformed into its state-space equivalent form
V = [A]V + [b]u (2.2)
where u is the input vector of forces, [a] is the state variable
matrix, and [b] is the input distribution matrix.


6
Figure 3 illustrates the axial vibration model for the
undamped (free body) case.
Figure 3- Example of a model for an undamped LSS
Transforming the axial vibration model into a set of second
order differential equations,
mix] = K^xi + K2(x2 - x^)
m2X2 = - K2(x2 - x-|) + K^(x3 - x2)
m3x*3 = - K3U3 - x2) + K4_(x4 - X3) (2.3)
m4*x*4 = - K4U4 - X3) f(t)
Equation (2.3) can be placed into the state matrix equivalent
system of equation (2.2) as follows:
m-| 0 0 0 (Ki + K2) -k2 0 0 "0
0 m2 0 0 -k2 (k2 + K3) , k3 n 0 0
0 0 m3 0 *t + 0 -K3 (K3 + k4) -k4 *t = 0
0 0 0 1114 0 0 -k4 k4 -1
given [m]x(t) + [K]x(t) = [b]f(t) (2.4)
set V = [ x x ]T (2.5)
so
0
I
-nr1K
V
0
nT^b
0
f(t)
T
(2.6)


7
To find the values for the stiffness matrix [k], the
equivalent spring values, for each node, need to be determined. The
equivalent spring for a uniform bar in elastic tension-compression is
shown below.
Figure 4- Equivalent spring diagram
P A*E
k = = ---- (2.7)
£ L
Where A is the cross sectional area
E is the modules of Elasticity
L is the length of the bar
Note, if complete analysis on the model is required, the
equivalent mass and springs constants in the bending and torsional
directions would also need to be analyzed. However, for this thesis
just the axial modes (modes in the tension-compression direction)
will be controlled. The bending and torsional vibration modes can be
modeled into their equivalent mass-spring-damper system in an
analogous fashion. Therefore, any of the principles applied to the
axial vibration model can likewise be applied for torsional and
bending modes.


8
For this model the stiffness values for K are derived as
follows.
E*A
ki = ------------- (2.8)
( ^i ^i-1)
Where kj_ is the equivalent spring constant for the ith node of the
model.
is the equivalent area for the ith node,
li lj__-| is the length between the ith node and the i-1 node
To model the individual k values, a cross sectional area of
5 x 105 m2 was chosen, and each nodal segment was selected to be
10 meters in length. The value for E was chosen as 2.5 x 10 N/m^.
This is 1000 times less than the true value for aluminum. This was
done for computational purposes only (i.e., to achieve lower
frequency modes than would be expected for a real physical system).
The IBM-PC computer used had some problems in finding the eigenvalues
for the higher frequency modes due to limits in the machine. The
resulting values for the k's for each node or segment was 750 N/m2.
Figure 5 on the following page contains a summary of the
formulas used for calculating the lumped mass into nodes [2].


9
- concentrated mi
m4 diitributed mm
m m tumped matt of equivalent tyMem
Figure 5. Formulas for lumped systems
Source: Shigley, J.E., Mechanical Engineering Design,
McGraw-Hill, New Jersey, 1977.
The approximate distributed mass of each section is
md = P*A*L = 1 Kg; where p is the density of around 3000 Kg/m^.
The formulas above were used with many liberal approximations in
order to create the model. Specifics are not considered, since this
thesis concentrates on controlling a space system, not modeling it.
Thus, the lumped mass "m" for the equivalent modeled system was
selected as 1 Kg for the first three nodes. The mass at the 4th node
(the tip) will have a lumped mass of m = 1.5 Kg. This represents the
addition of a concentrated mass "mc" of the control actuator and
sensor as well as the distributed mass.
The physical state-space matrix system in the form of
equation (2.2) with the numerical values for all k's and m's is as
follows.


10
v(t) =
0 0 0 0 -1500 750 0 0 0
0 0 0 0 750 -1500 750 0 0
0 0 0 0 0 750 -1500 750 0
0 0 0 0 0 0 500 -500 0
1 0 0 0 0 0 0 0 V(t) + .66667
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 0 0
V = [ X X ]T and V = C X X ]T.
f(t)
In order to simulate this system on a computer, it is
converted into its discrete equivalent form. This was done by
running a version of the Standford Program "CD.BAS" which was
modified to make it more user friendly; it is called "CD2.BAS". This
program uses the series expansion method to convert the matrices into
their discrete equivalent. The sampling rate was chosen as .05
seconds (20 Hz). As will be shown in later sections, this is a good
choice for sampling rate since it is greater than twice the highest
modeled frequency of the structure. Note, by Shannons sampling
theorem, for a mode to be controllable it must be sampled at a rate
equal or greater than twice the mode's frequency. See Appendix A for
the listing of the program.
2.1 MODAL THEORY
Given the second order physical system model of a linear
structure
[m]x + [c]X + [k]x = [b]u (2.1.1)
The structural eigenvalue problem associated, with the X = (^e^
solutions, for the undamped free vibration case of
[m]x + [k]X = 0
(2.1.2)


11
is det[ m\2 _ k ] = 0, eigenvalues: {Xi^ = ~C0i^> i = 1,2,...n}
(2.1.3)
and [ mX^2 k ]

(2.1.4)
For this case, all eigenvalues (Xq = iC0, i^ = -1 ) are
purely imaginary and eigenvectors are real [3]. These eigenvectors
make up the modal matrix and determine the mode shapes. The
corresponding values of CO which solve (2.1.3) are the modal
frequencies.
The coordinate transform between the structure's physical
coordinates (x) and the modal (sometimes called principal)
coordinates (7"p is given by
X = T7} (2.1.5)
where T is the modal matrix made up of the eigenvectors
T [ Cp 1 (p 2> >
Now, by substitution of equation (2.1.5) into the state-space
representation of (2.1.2) the following system is derived
[ -T-1m-1cT ] ' [-T1m1kT]
[I] [0]
~[-T1m-1b]~
B =
L [0] J
(2.1.6)
(2.1.7)
For this thesis, an undamped system is modeled. For the
undamped case the damping matrix [c] is zero. Therefore, the upper
left term in equation (2.1.6) will likewise be zero. The resulting
state-space representation for dynamic matrix [a] is as follows
(matrix [b] remains unchanged):


A
[O] ' [-T-1m-1kT]
_ [I] [o]
12
(2.1.8)
To convert the physical system state-space form into the
modal state-space form, a satisfactory solution for the eigenvector
matrix [t] must first he found. To find [t], a version of the
Stanford Program "MODALCOMP" called "MC2.BAS' was used which again
was modified to make more user friendly. The program uses the
numerical analysis technique of Leverrier, Bairstow, and Gauss to
calculate the eigenvalues, the eigenvector matrix [t], the inverse
eigenvector matrix [t]-1, and the modal dynamics matrix [t]-1[f][t].
Where [p] is simply the dynamics matrix of the state space form:
x = [f]x + [g]u
Note, to make this program work, [p] must he set as nf^k; then the
problem will he in its equivalent state space form (i.e., the
eigenvalues found solve the equation det[\2 m^k] = 0 ). See
Appendix A for a program listing and a computer printout of the
results.
Prom the results of the program, the Modal state-space matrix
system can he constructed. The program gives the results for the
upper right partition of dynamics matrix [a] in equation (2.1.8).
The Modal state-space input distribution matrix [b] is calculated as
follows:
[b] = [T1m-1b]
-
0 .08833
0 = .18642
0 -.27711
.666667 .36828
(2.1.9)


13
The modal state-space matrix system for the model is shown
below in equation (2.1.10).
f(t) -
(2.1.10)
0 0 0 0 -2610 06 0 0 0 .08833
0 0 0 0 0 -1658. 20 0 0 .18642
0 0 0 0 0 0 -657 63 0 -.27711
0 0 0 0 0 0 0 -74 .111 .36828
1 0 0 0 0 0 0 0 <(t) * 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 0 0
l(t)
Note : X = T7J
X = TT}
i
(2.1.11)
As before, the state-space matrix was converted into its
discrete time equivalent using "CD2.BAS". The results can be found
in Appendix A.
The eigenvector matrix found by the modified Standford
program is not unique. In fact there is an infinite number of
eigenvector matrices that can exist and still be correct; any
eigenvector may be scaled by some constant:
ur = cr <£r (2.1 .12)
Where cr is a scaling constant whose units are such that cj} r^m(pT has
the dimension of mass, and Ur is the scaled equivalent eigenvector.
It is well known that the eigenvectors are orthogonal with
respect to m and k. In other words
r^m<^s = S^r^S^s = 0 where s i r. (2.1.13)
These eigenvectors can be normalized by using the orthogonality
conditions (sometime called orthonormal relations) to scale the
eigenvector matrix into its normalized form. This is quite simply


U
done by scaling the rth mode so that the generalized mass, or modal
mass, defined by
Mr = UrTmUr
Mr = cr(pr^m(pror
has a specified value of 1. In other words, the modal mass
[m] = [i]. It isn't hard by substituting equation (2.1.12) into
(2.1.14) to derive the following equation for the scaling constant
for each eigenvector:
(2.1.14)
(2.1.15)
c
(2.1.16)
For this system the scaling constants were calculated to be
ci =.680713, c2 = .672263, C3 = .612563, and C4 = .606862.
By multiplying each eigenvector by its scaling constant the
normalized eigenvector (modal) matrix is calculated as
T-
n
-.459917
.680713
-.547586
.129758
.672263
-.141807
-.642347
.277302
.545389
.612563
.142617
-.452384
.19771
.375884
.516913
.606862
In the normalized modal matrix form, the following properties
are true.
[Tnh1 = [Tn]T[m]
[l] = [Tn]T[m][Tn]
[U)2] = [TnMk][Tn]
From equation (2.1.17),
Tn"1
-.459917
.672262
545390
.197712
.680712
-.141805
.612563
375882
-.547588
-.642346
.142617
.516913
.194636
415953
-.678576
.606862
(2.1.17)
(2.1.18)
(2.1.19)


15
Using the normalized equations, the state-space modal equations are
as follows.
<(t) -
-[Tn]T[c][Tn] -[Tn]T[k][Tn] > + [TnMbI
[l] 0 f(t) 0
u
(t)
or
f (t)
-[TnMc][Tn] -[0)2] + ' OnKb]'
. Cl] (t)
0 0
u(t)
since there is no inherent damping [c] in the model, the equations
can he simplified to
? (t)
0 -[0)2]" ~[Tn]T[-b]
_ [I] 0 ?(t) 0
u(t)
2.2 MODE SUPERPOSITION FORM
Another technique for evaluating the modal vibrations in a
structure is called the modal superposition form. This technique
uses mode shapes that are superimposed on each other to define the
deflections in the structure. The deflections for the structure is
given by
y(t,x)
(2.2.1)
where t is time, x is position location, 1 is upper limit of modes (4
for the model system), ^(x) is the kth mode shape at position x,
and Ui^t) is the kth modal amplitude as a function of time. If
there are point actuators located at positions x X4 then the
modal amplitudes each satisfy are in equation (2.2.2)


16
(2.2.2)
VkCt) + 2£kfLkUk(t) +sTlkUk(t) / ui(t)v//k(xi), k = 1.................4
i^1
where J tk and ^ k are the kth modal frequency and damping
coefficients respectively; u-^(t) is the ith actuator force on the
structure [3].
2.3 COMPUTER SIMULATION
To simulate the mathematical models, programs were written in
BASIC and Pascal (see Appendix B for the listing of the simulation
programs). Both the modal state matrix form and the physical state
matrix form of the physical system are simulated. The data gathered
is the velocity and displacement of the tip of the structure with
respect to time. To generate the vibrations in the structure, a unit
impulse input "u" of 10 newtons applied at the tip is simulated.
This impulse force could be representative of a force caused by
firing of retro rockets on the space shuttle, docking of another
space vehicle, etc.
As mentioned previously, the mathematical models are first
converted into their discrete equivalent form. This was done by
using the program "CD2.BAS", which is a modified version of the
Standford Program called "Condistran".
Figures 6 and 7 illustrate the graphical representation of
tip velocity versus time, and tip displacement versus time, for the


17
physical state matrix form and the modal state matrix form of the
system model respectively. As was expected, the results for the two
models are virtually identical.


,3
TIP VELOCITY vs TIKE
Figure 6. Computer Simulation in the Physical Matrix Form.


,3
IIP VELOCITY vs TIKE
Figure 7. Computer Simulation in the Modal Matrix Form.


CHAPTER III
DIGITAL FILTERING
In order to use modal control, it is necessary to separate
the system dynamic vibration response into its separate modes. To do
this, causal digital filtering techniques are used. Causal digital
techniques can he broken down into two basic forms: nonrecursive and
recursive. Butterworth recursive filters are used for this thesis.
Recursive filtering has the advantage of being faster to implement in
real time. The basic transfer function equation for the recursive
filter is shown below.
Y(z) b0 + b^ z1 + + bgZ-cl
HipU) ------------------:--------------- (3-D
X(z) 1 + a^z1 + apZP
Where Hqp(z) designates a low pass filter. A digital filter equation
such as (5.1) is sometimes modified into the following equivalent
form:
B(1 + b'jz-^ + + bqZ_cl)
H(z) = -------------------------- (3-2)
1 + a^z_1 + + apZP
For the physical system, the structure will be controlled by
implementing sensor and controllers (actuators) on the end or tip of
the structure. Therefore, only the velocity value at the tip needs
to be filtered.
The digital filter can be implemented on a computer by the
following simple algorithm formula,


21
Y(n) =
y~Bk*X(n-k)
k=0
P
J^Ak*Y(n-k)
k=1
(3-5)
where Y(n) is the output at time n, X(n) is the input state at time
n> Bk is the transfer functions numerator coefficients, and Ak is the
denominator coefficients. The order of the numerator is q and the
denominator order is p. In this case, Y(n) will he the modal
velocity of the desired mode at the tip. X(n) will be the actual
velocity of the tip.
The first mode (lowest frequency mode) can be filtered by
using a low pass filter (LPF). All other modes must use a band pass
filter (BPF).
In order to determine the filter transfer function equation,
a program was written that calculates the poles of the Butterworth in
both the Z and S domain, as well as the required gain. See Appendix
A for a program listing. The zeros of a Butterworth filter are all
at (z = -1) by definition. For the LPF design, a bandwidth of around
2.25 Hz is desired. This is slightly greater than the 1st modal
frequency as will be shown in the next chapter. The program was run
and it gave the results for the gain and denominator coefficients of
the digital filter. The numerator coefficients of the filter were
found by applying the following simple equation:
coefficient bk is
8
£
k=0
(n)<
(n k)! k!
(n is the order = 8)
(3.4)
The results are summarized as follows:


22
numerator coefficients denominator coefficients
bo = 1
bi = 8 a1 = -4.385951
b2 = 28 a2 = 8.920154
b3 = 56 a3 =-10.821606
b4 = 70 a4 = 8.489805
b5 = 56 a5 = -4.385521
^6 ss 28 a6 = 1.450788
b? = 8 a? = -.280147
b8 = 1 a8 = .024117
Gain is B = 5-328643 x 10"5
To compare the results of the filtered modes to that of the
ideal modes, the modal equation form of the ideal output of the
sensors can he found by using the equation below.
N N .
4i(t) - Sqij - (3.5)
i=i
Note: <^) is the eigenvector modal matrix
^ is the modal velocity vector
"s" subscript denotes the mode.
i" subscript denotes the node (nodal position of the sensor).
h
Thus for a chosen mode s qs(t) = $ is ^s (3-6)
is the predicted "ideal" output of the sensor at position "i" of the
system.
Programs were written that simulate both the ideal sensor
output using the modal state matrix system along with equation (3-6),
and the output of the physical sensor found by digital filtering the
physical matrix system.
Figure 8 shows the response of the digital filter, the ideal
filter response, and the difference (error) between them
respectively.


23
Figure 8. Low Pass Filter Simulation.


24
To create a BPF, a LPF Butterworth is designed with the
desired bandwidth, then the following transformation property is used
to shift the frequency response.
Hbp(z) = Hlp(e-jk>cz) + Hip(ejO;cZ) (3.7)
Thus, the order of the BPF will be twice that of the LPF.
The basic equation for a BPF is shown below.
bo + biz*^ + + bp+qZ-P-(l
Hbp(z) = ------------------------------ (3.8)
1 + a^z-^ + + a2pa2p
For the BPF design, a LPF with the appropriate bandwidth was
first designed by using the same filter program as before. Then the
BPF was calculated by utilizing equation (3*7) and putting it into
the form below,
(3.9)
Hbp(z) = B*-
nk=i (1" 2kejwc z-1)
B*-
H k=1 0 ~ Z]ze~3U)c Z-1)
^m=1 (1 Pme^C z-1) JT m=1 (1 pme"jWc z1)
where zk and pm are the LPF's zeros and poles respectively andGi)cis
the center frequency for the BPF. Using equation (3*9) to calculate
the transfer function, the resulting BPF is as follows.
numerator coefficients denominator coefficients
b0 = 1
bi = .3456895 a1 = -2.728852
b2 = -49.36748 a2 = 6.753185
b3 = -11.23136 a3 -10.36797
b4 = 286.1777 a4 = 14.35491
b5 = 39.3709 a5 = -14.79166
b6 = -420.0358 a6 = 13-79967
b7 = -30.828 a 7 = -9.915372
b8 = 190.56 a8 = 6.447651
b9 = 5.39381 a9 = -3-104981
b10= -21.87708 a10= 1.348126
b11 = -.102038 *11 -.3582316
bi2= .2951724 a12= .0871268
Gain is B = 2 x 8.576548 x 106 = 1.71531 x 10-5
I


25
As before, computer simulations were made of the band pass
filter output, the ideal modal filter output, and the difference
(error) between them. Figure 9 on the following page shows the
results.


Figure 9. Band Pass Filter Simulation .


27
By examining both the 1st mode's LPF and the 2nd mode's BPF
difference plots, it is obvious that the digital filter output
differs quite a bit from the ideal filter output at the start of the
simulation. In the simulations, the output of the digital filter was
different than the output of the ideal filter for about one second,
after this the outputs of the two filters match (relatively
speaking). It appears as though the filter undergoes a "settling
time" in which the filter becomes equivalent to the ideal filter.
Part of the reason for this is the filter coefficients do not contain
an initial value. Therefore, the filter will need as many sample
cycles as the order of the filter before it gives truly useful data.
The intuitive way to reduce this settling time is to increase
the sampling rate. This would lessen the time it takes for all of
the coefficients of the filter to contain an initial value and for
the filter to produce legitimate output information. Unfortunately,
this method is not without penalties. If the sampling rate is
increased, the spacing between the modes is decreased in the digital
frequency plot. As the spacing between the modes is decreased, the
digital filter design will need to be of higher order to prevent
overlapping (spillover) of the modes. A higher order filter will
need more samples cycles for all its coefficients to contain initial
values; thus, it is a no win situation. Figure 10 illustrates this
point.


28
MAGNITUDE VERSUS FREQUENCY
fs = 20 Hz
MAGNITUDE VERSUS FREQUENCY
fs = 40 Hz
Figure 10. Comparison of Filter Frequency Response


CHAPTER IV
MODAL CONTROL
As was shown in the section on Modal Theory, the modal state
matrix form is a diagonal dominant matrix system. The upper left
partition of matrix [a] contains the damping terms. The damping form
can he put into its 2nd order equivalent in the usual way:
damping term = -2 C kWk (4.1)
where C k and OJ ^ is the damping coefficient and frequency for the
kth order mode. Eor the 1st mode, C\) i was shown to he 8.609 rad/sec.
If the damping coefficient is chosen as ^] = .7, then the damping
term is calculated by using equations (4.1) as -12.0526. The modal
domain state mat rix system is now
o 0 0 0 -2610.06 0 0 0
0 0 0 0 0 -1658.2 0 0
0 0 0 0 0 0 -657.6 0
0 0 0 -12.0526 0 0 0 -74.1
Ac = 1 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0
0 0 1 0 0 0 0 o
0 0 0 1 0 0 0 0 J
(4.2)
The modal form input distribution matrix B is the same as in the
Modal Theory section.
If the input is set as u = -GC^ where G is the gain matrix
and C is the output distribution matrix, then the closed loop form of
the modal state space equation can be written as follows:
f(t) = £c]?(t) + B*fd(t)
Y(t) 'cf(t)
(4.5)


30
where fd(t) is the disturbance force,
Ac = [ A GC ],
and
G =
0
0
0
12.0526
0
0
0
0
C=[0001 0000 ].
(4.4)
(4.5).
The same process can be used to simulate damping of the 1st
and 2nd modes together. The 2nd mode frequency is CJ 2 = 25.644
rad/sec. Prom equation (4.1) C*2 = *7, thus -2^>2^2 = -35*9016.
The modal state space form for the closed loop matrix Ac is now as
follows (note, again matrix B does not change):
0 0 0 0 2610.06 0 0 0
0 0 0 0 0 -1658.2 0 0
0 0 -35.902 0 0 0 -657.6 0
0 0 0 - 12.0526 0 0 0 -74.1
1 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0
0 0 : 1 0 0 0 0 0
0 0 0 1 0 0 0 0
The modally damped systems are now simulated. As before the
systems are first transformed into their discrete equivalent forms
using the Standford program Condistran. The following pages show the
plot of the "Ideal Modal Response" at the tip of the system with the
1st mode dampening coefficient of .70; and also with both the 1st and
2nd modes dampened with a coefficient of .70. See Appendix B for a
listing of the simulation programs.


.3
IIP VELOCITY vs TIME
IDEAL MODAL SOLUTION HIIH 1st MODE DAMPENED ? .7
Figure 11. 1st Modal Damped System Simulation.
LO


,3
TIP VELOCITY vs TIME
IDEAL MODAL SOLUTION
1st AND 2nd MODES DAMPENED 0 ,7
-.3
TIME (seconds)
Figure 12. 1st and 2nd Modal Damped System Simulation.
LO
N3


53
Unfortunately, even though modal control is fairly easy to
understand and very easy to show mathematically, it is very difficult
to implement physically. This is because it requires distributed
controllers, and the best that can be done is to simulate the
distributed controller by a series of discrete point actuators.
For the model, suppose it is desired that the physical model
is dampened so-the response is the same as the modal model with the
1st mode dampened at 0.7* The calculations are as follows: First,
the diagonal elements of the modal matrix are set equal to the
required modal gain matrix.
Km = [ -T1nr1cT ] = diagonal { -2^Cl)i ). i = 1, , 4
Then the physical gain matrix is found by use of the following
equation.
Kp = [ -T*2£c0iT-1 ] = -TKmT-1
Note: Km = T"1KpT
Kp TKmT-1
(4.7)
(4.8)
(4.9)
For the first mode
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 -12.053
Thus,
-.47114 -.8957 -1.232 -2.169
-.89568 -1.703 -2.3419 -4.124
-1 .2317 -2.342 -3.220 -5-671
-1.446 -2.749 -3.781 -6.658


34
The same basic equations can be carried out in the normalized
eigenvalue form. For the normalized form the modal gain is
Kmnorm = -[Tn]T[c][Tn] = diagonal { -2^0J }, i = 1,...,4 (4.10)
so Kmnorm = Km
Thus the normalized physical gain matrix is
Kpnorm = (4-11)
Note: Kpnorm = [Tn]^[Km][Tn]^ (normalized physical gain matrix)
Km = [Tn]T[Kpnornl][Tn] (modal gain matrix)
By taking advantages of some of the properties of normalized
eigenvector matrices, the following equalities can be found.
Kpnorm = [Tn][m][Km][Tn]1 (4.12)
= [Tn][m][Km][Tn]T[m]
The normalized physical gain matrix is easily calculated.
Kpnorm
-.471137 -.8957 -1.252 -2.1692
-.89571 -1.7029 -2.5418 -4.1240
-1.25177 -2.5418 -5.2205 -5.6713
-2.16918 -4.1240 -5.6713 -9.9872
Note that the normalized physical gain matrix is equal to the
previously found gain matrix except for the multiplication of m.
Kpnorm = Kp[m] (4.13)
Unfortunately, to implement these physical gain matrices
(Kp and Kpnorm) the physical system would need to have an actuator
and a sensor on each node. As previously stated, for this physical
structure model only the tip of the structure contains the sensor and
actuators. Now the question is, is this system modal-state
controllable by the use of a single actuator location with a direct


35
proportional feedback gain matrix Kp. That is, could the physical
system response be the same as the ideal modal response with the
single tip actuator using proportional velocity feedback. The answer
is quite obviously no. The physical form of the gain matrix
(Kp or Kpnorm) could not possibly be implemented with a single
actuator located only at the tip (node 4) This is also obvious by
looking at the feedback fb(k) = B*C*gain. Since the input
distribution matrix B = [ 0 0 0 1 0000 ]Tf ana the observer
distribution matrix C = [ 0 0 0 1 0000] for the first mode, only
the gain matrix coefficient k^^ can contain a non-zero value as shown
below.
0 0 0 0 I
0 0 0 0 1
B*C*gain = 0 0 0 0 1 0 (4.14)
0 0 0 k44 1

0 1 0
So the question is, how can k^^ be chosen such that the
physical system will most closely match the modal controlled system.
The PACOSS project designers used the following formula to generate
the feedback gain to achieve the desired active damping in the
specific modes [4].
K = T-T[ 2 foc OJoc]T-1 (4.15)
na x nc 00 00
where Toc is the normalized modal matrix of observed and controlled
modes, and ^*oc is the desired active damping ratio of the controlled
modes. Note that the formulation of equation (4*15) assumes the
number of actuators, na, is equal to the number of controlled modes,
nc. If na is greater than nc, the gain matrix is overspecified and a


36
least squares (minimum norm) solution may be used. If the na is less
than nc, the gain matrix is underspecified and a pseudoinverse
minimum norm approach would have to be-used.
The system presented in this thesis has more controlled modes
than actuators since actuator control is limited to the tip of the
structure; so this "underspecified" system must use a pseudoinverse
approach. The problem with the pseudoinverse approach is that the
system will show a loss of closed loop damping ratios from the design
values when this approach is used. In a later section it will be
shown that this is certainly true for this system.
Even though the system cannot be modal-state controllable by
using a single point actuator, there are some possible advantages of
trying to dampen the structure by independent modes. For a large
space structure certain modes may not be important and need not be
considered. If they are not important, no control energy should be
wasted by controlling them. Therefore, some of the possible
independent modal solutions methods mentioned in this section can
still hold merit for this single point actuator problem.
In section 4*4 another method will be considered that is
different than equation (4.15) This method will use the classical
approach of the root locus technique to achieve the best possible
results in terms of transient responses. It will be shown that this
method has some benefits over equation (4-15) for an underspecified
case, since the transient response can easily be optimized.


37
4.1 EFFECTS OF MODELING ERRORS, RESIDUAL MODES, AND SPILLOVER
When trying to use modal control on a structure, the
uncertainties of modelling errors, residual modes and spillover
effects must be considered. For this system a 4th order model was
used to mathematically represent the physical system. Thus, there
are only 4 modes in the model. In an actual structure there are a
large number of modes. The effect of these unmodeled modes, called
residual modes, caused by this modeling truncation must be examined.
When the sensor outputs are contaminated by the residual
modes through the term C*V(k), where V(k) is the velocity output of
the sensor, this is called observation spillover. When the feedback
control excites the residual modes through the term B*f(k), where
f(k) is the feedback control force, this is called control spillover.
The effects of the spillover can be summarized as follows. When only
control spillover is present, the effect is to excite the residual
modes at their natural frequencies. When both observer and control
spillover is present, the results are much more complicated and can
possibly lead to instabilities in the system [5].
For the system used in this thesis, only the 1st and 2nd
modes will be controlled. The 3rd and 4th modes can be thought of as
the residual modes.
So far in the discussions the actuators and sensors have been
considered collocated. The question of the effect of this constraint
on the control system design must be considered. For a real physical
system collocated sensor and actuator pairs are more desirable than
non-collocated sensors and actuators due to their robustness. If it


38
is not possible to collocate the sensors and actuators due to some
physical design constraint, then there is a greater potential for
generating difficult control problems. For example, if there is
modeling error or model truncation there is a possibility that the
control scheme for the non-collocated system could produce an
unstable system. On the other hand, a collocated sensor/actuator
system (assuming ideal sensors and actuators) will always produce a
positive definite system, regardless of modeling errors [4]
Therefore, the collocated system will be asymptotically stable by the
definition of Liapunov.
For this paper the false assumption that the actuators (and
sensors) are ideal and have unity gain is intentionally made. In
reality, a physical actuator transfer function is usually non-linear
but can be approximated with a linear 1st, 2nd, or higher order
transfer function (TF). A simple 1st order TF for an actuator can be
as follows: .
TF(s) = ( tcons-{;*s + 1)^ (4.1.1)
Note that this TF tends to create a phase lag into the system.
Therefore, it should be kept in mind that the assumptions of ideal
actuators and sensors are not totally valid in the real sense.
Nevertheless, these assumptions will be used in this thesis for
analytical purposes.


39
4.2 DISTRIBUTED ACTUATOR CONTROL
There has been work done recently in developing a distributed
actuator for controlling flexible beams. The current design is to
use a piezoelectric polymer film that can be polarized, or made
piezo-electrically active [6]. Voltages or fields applied across the
film results in a longitudinal strain. The field can be varied
spatially; likewise the strain will vary spatially. It is possible
that a control system could be designed where the applied field could
be spatially related to the mode shape ^(x) so the feedback will be
applied in the form of a distributed mode shape. If there is an
accurate knowledge of the mode shapes for a given system, then this
system should work. Plus there would be several advantages:
theoretically, only one sensor is needed to control the modes (each
modal frequency would need to be filtered), the piezoelectric film
should be much more reliable than the modern day actuator, and the
control system would be inherently robust due to the basic physics of
the system (no moving mechanical parts for the controller).
. Thus, in the future it is quite possible that a true
distributed control load could be applied, and modal control could be
implemented directly. This may only be feasible for simple linear
beam type structures without complex mode shapes and would work best
for the bending modes; however, theoretically it is possible that it
could work for the axial or rotational modes also.


40
4-3 PASSIVE CONTROL
So far this thesis has only mentioned the principles of
actively controlling large flexible structures. In reality, as much
design time, if not more, will be spent on considering ways to use
passive damping in the LSS design. This is due to the much greater
reliability of a passive dampened structure versus an actively
dampened structure.
Every structure will possess some inherent damping. For an
aluminum structure, an inherent dampening ratio of 0.2$ is often
assumed. There is a lot of research going on in developing new
materials in the form of composites that can be used to increase this
damping ratio. In addition, there are several studies and articles
written on how to find the optimal design of a structure to maximize
the passive damping as well as find the actuator locations on the
structure that will produce maximum active dampening effects.
By using passive dampening together with active control,
several advantages can be achieved: a reduction in the energy
requirements needed to control the structure, significant reductions
in maximum actuator torque requirements, less complex active control
systems involving fewer actuators and sensors, and weight savings due
to smaller and fewer actuators needed. All these benefits lead to
more robust and reliable systems which will more efficiently
accomplish their mission goals [4]*


41
4-4 USING THE CLASSICAL TECHNIQUES OF ROOT LOCUS
In Chapter 3> it was shown how the vibration dynamics of a
physical system could he modeled by digital filtering each of the
dominant modes. It may be desireable to control the system by
independently damping the dominant modes. This requires finding a
suitable gain for each of the modes.
One simple way to find the suitable modal gains is to use the
classical Root-Locus technique. The Bairstow numerical algorithm can
be used in conjunction with the Root-Locus to numerically find the
exact poles at certain gains. To illustrate this technique, the
direct velocity feedback method with a proportional feedback gain is
first considered. The feedback will be generated by velocity
readings taken from the sensor on the tip of the structure. The
readings are then fed to the control actuators. The control f
actuators are also located on the tip of the structure. This is
often referred to as collocated sensor-actuator control. The
illustration below shows the problem.
Figure 13* Proportional Tip Velocity Feedback
Where R is the disturbance input, V@4 is the velocity at the tip
(node 4), and K is the gain.


42
To find a suitable gain, the basic principal of Root-Locus
was used. Instead of using the graphical techniques, the computer's
ability to mathematically find the poles using the numerical analysis
method called Bairstow's algorithm is utilized. Using this
algorithm, a program was written that will find the zeros (closed
loop poles) of the function
[ D(z) + K*N(z) ] = 0 (4.4.1)
which is the denominator of the closed loop transfer function. That
is G(z) = N(z) (4.4.2)
where D(z) is the denominator coefficients and N(z) is the numerator
coefficients of the open loop system transfer function G(z)
respectively; K is the proportional feedback gain. By using the
program, the Root-Locus plot can be directly created by entering
different values for K and receiving the pole locations for the modes
at that K. The system transfer function G(z) was obtained by using
the modified version of the Stanford Program "Condistran". The
program uses the recursive Leverrier algorithm for finding the
transfer function denominator and numerator polynomial coefficients.
To determine the best (optimum) value for gain K that gives
the fastest transient settling time, the formulas listed in equations
(4.4.3) through (4.4.6) can be used. Note that in the Z-domain all
points inside the unit circle are stable. The closer to the origin
these points lie, the faster the settling time for the system. If
the poles (roots of the equations) are plotted in the Z-plane, the
polar results are as follows:


43
the polar form N 11 CD (4.4.3)
£> ln r
the damping ratio t = I (4.4.4)
"Vin2r 02
the natural frequency n = 1/tVln^r + 02 (4.4.5)
the time constant T = 1/tu) = t/ln r (4.4.6)
where t is the sampling period.
Obviously, from equation (4-4.6), to get the best transient
response the time constant, and thus the magnitude of r, must be
minimized.
Figure 14 shows the root locus plot as well as the numerical
values obtained for different gains in the z-domain. The real and
imaginary coordinates of the poles of the characteristic equation are
printed for the different values of the gain K. These coordinates
can easily be changed into the polar form of equation (4.4.3), and
the time constant for each gain K can be calculated.


44
Z-llane
Figure 14. Root-Locus Plot Using Direct Proportional
Velocity Feedback.


45
From Figure 14 it is evident that a proportional gain of
K= 32 is a very good choice to dampen the vibration modes (especially
for modes 1 and 2). Unfortunately, due to physical constraints, a
gain of 32 may he more than a real actuator could produce. This is
one of the reasons why modal control will he important for future
large space structures. With modal control, only certain modes will
need to he dampened (usually the lowest order modes). No effort will
he wasted on damping modes that are not dominant and do not have a
large effect on the overall structure. As will he shown later, when
modal control is used, the gains will he noticeably smaller (albeit,
the transient vibration damping will not he quite as good).
To simulate the system using direct proportional velocity
feedback, the discrete form of the state space matrix system is used
where X(k+1) = [Ad]*l(k) + [Bd]*R(k) Kp*[Bd]*V@4(k) (4.4.7)
so u(k) = [ R(k) Kp*V@4(k) ],
[Ad] and [Bd] are the discrete forms for the state space matrices,
Kp is the proportional feedback gain, and V@4(k) is the sensor
velocity readings at the tip of the structure. Note that V@4(k) is
the 4th element of vector X(k)
(i.e., X = [ X! X2 X3 V@4 X5 X6 X7 X8 ]* ).
Figure 15 shows the transient vibration behavior of the
proportional tip velocity feedback problem. The system was again
initially agitated by applying a 10 newton impulse on the tip of the
structure. The top graph shows the response using the close to
optimal gain value of 32. The bottom graph shows the response of a
gain value of 10 for comparison purposes.


,3
46
IIP VELOCITY vs TIME
,3 TIP VELOCITY vs TIME
T DIRECT VELOCITY FEEDBACK
WITH GAIN K : 10
i
-.3 _1_ TIME (seconds)
Figure 15. Tip Velocity Feedback with Proportional Gain.


47
For modal control, the individual modal velocity can be
found by filtering. Then the mode can be independently controlled by
multiplying the modal velocity by a suitable feedback gain. For the
first mode, a low pass filter can be used. Again the tip of the
structure is used and the velocity sensor and control actuators are
collocated. The illustration below shows the system.
Figure 16. Modal Tip Velocity Feedback
Here, designates a gain for the first mode. This time the
characteristic equation to be solved is
[ D(z)*DF(z) + K-|*N(z)*NF(z) ] = 0. (4-4.8)
Where DF(z) and NF(z) are the low pass filter transfer function
denominator and numerator coefficients respectively.
It is expected that since only the first mode is being
controlled, the other three modes will not vary a lot as the gain
value is changed. By looking at the actual plot on the following


48
page, this appears to he the case. Also note that a gain of around
Ki = 4*5 gives the best transient time response.
Figure 17 on the following page show the results of the
Z-plane poles for different K-| values on the Root-Locus plot.


49
Z-Plane
Figure 17.
Root-Locus with the
1st Mode Filtered.


50
To simulate this system, where the 1st mode is dampened using
actuator control at the tip of the structure, the discrete state
space matrix system is shown below.
X(k+1) = [Ad]*X(k) + [Bd]*E(k) K1*[Bd]*V1(k) (4-4.8)
note: u(k) = [ R(k) K-j^V^k) ]
K1 is the 1st modal feedback gain (chosen as 4.5), and Vi (k) is the
1st mode's velocity (the modal velocity) at time pulse (k) at the tip
of the structure. As before, the system is agitated with a 10 newton
impulse force.
Figure 18 contains two plots: the top plot shows the tip
velocity versus time; the bottom plot shows the tip velocity for the
first mode only. See the Appendix B for listings of the simulation
programs used to develop the data for the plots.


3
TIP VELOCITY vs TIHE
FEEDBACK GAIN Hi = 4,5
51
,3 TIP MODAL VELOCITY vs TIME


52
After examining the top plot of Figure 18, it is obvious that
the system vibration response will not be suitable with only one mode
controlled. Therefore, the next mode of highest importance must be
looked at. In most cases, this will be the next mode (mode 2). As
stated previously, this mode and all other higher modes will need to
be filtered with a band pass filter. The control scheme will now
include an additional BPF and another feedback loop. Figure .19
illustrates the system.
Figure 19- Independent Modal Control For The 1st and 2nd Modes
Since the root locus technique only works when one gain is
changed at a time, only the root locus plot using the feedback loop
for the second mode is generated. As was shown earlier, when
changing the gain for a particular mode, the other modes are not
effected to a great extent. They tend to be mutually exclusive.
Therefore, with narrow band filtering it should be'possible to
generate a root locus plot for one mode, find the suitable gain.


53
then use another independent root locus plot for a different mode;
and find a suitable gain for this mode.
Figure 20 illustrates the system used to generate the root
locus plot for the second mode.
Figure 20. Diagram for the 2nd Mode
Table 1 shows the numerical results for the Root-Locus plot.
The table shows how the magnitude of the z-plane poles change as the
gain K.2 was changed. As predicted, by changing the gain K2, only
mode 2 was effected; that is, the other modes were not effected in a
relative sense.
Gain K r for mode 1 r for mode 2 r for mode 5
0 1 1 1
1 1 .99594 1
2 .99999 .99161 1
3 .99999 .98708 1
4 1 .98262 1
5 1 .97945 1
5-7 99999 97874 1
6 .99999 .97889 1
7 .99999 .98085 1
10 99999 .99210 1
15 .99999 1.01002 (unstable!) 1
Table 1. Z-plane Modal Magnitude versus Gain K2


54
Figure 21 show the transient vibration behavior of the system
using modal control to dampen the 1st and 2nd mode at the tip of the
structure. The first modal feedback gain Ki is chosen as 4*5, the
second modal feedback gain K2 is chosen as 5*70. As before, the
system is agitated with a 10 newton impulse force.
To simulate this system, the discrete state space matrix
system is shown below:
(4.4.9)
X(k+1) = [Ad]*X(k) + [Bd]*R(k) K^fBdJ^V! (k) K2*[Bd]*V2(k)
where K-j and K2 are the 1st and 2nd modal feedback gains, and V^(k)
and V2(k) are the first and second modes' velocities (modal
velocities) at time pulse (k) at the tip of the structure.
Note: u(k) = [ R(k) K^V^k) K2*V2(k) ].
Two plots are shown in Figure 21. The top plot shows the tip
velocity versus time; the bottom plot shows the tip velocity for the
second mode only. The tip velocity plot for the first mode is
virtually identical to what was found in Figure 18. See Appendix B
for listings of the simulation programs used to develop the data for
the plots.


Figure 21. Simulation with 1st and 2nd Modes Dampened.


56
4*5 FINDINGS FROM DIRECT VELOCITY FEEDBACK SIMULATIONS AND
INDEPENDENT MODAL VELOCITY FEEDBACK SIMULATIONS
So far this thesis has looked into the theories and
simulated results of using independent modal control and direct
velocity feedback. Before proceeding into the subject of adaptive
control, it might be useful to take a further look into what can be
found from the results of the simulations. Figures 22 and 23
contain plots of the control input required for the direct velocity
methods and the control input required for the independent modal
control technique.
The top plot of Figure 22 shows the control input for direct
velocity feedback simulation with a gain of 10. The bottom plot
shows the control input for direct velocity feedback simulation with
a gain of 52; note that the scale for this particular graph is
8 to -8 Newtons, all the other plots of Figures 22 and 23 have a
scale of 4 to -4 Newtons. The plots of Figure 23 show the control
input for independent modal control. Both the first and second
modal control inputs are plotted.
After reviewing the control input plots as well as the
system response plots from previous sections, some obvious findings
and conclusions can be made about the usefulness of independent
modal control over direct velocity feedback.


4-
57
CONTROL INPUT vs TIHE
Figure 22. Control Input for Direct Velocity Feedback.


4
CONTROL INPUT vs TIME.
HODAL VELOCITY FEEDBACK 2nd HODE
HUH GAIN K2 : 5,?
58
4 __1_
TIHF. (seconds)
Figure 23. Control Input for Modal Velocity Feedback.


59
Some of the advantages of Independent Modal Velocity Feedback
over that of direct velocity feedback are listed below:
1 No control energy is wasted on unimportant modes. That is, the
modes that do not have a large effect on overall structure stability
will not be included in the control loop.
2. The maximum control force requirement is much less.
3. The rate of change in control force required is not as large and
is smooth in nature. In other words, the rate of change in
acceleration (sometimes called a jerk) required of the control
actuator is not as overwhelming as it is for direct velocity feedback
control. This will allow for a less rigid (beefed up) actuator
requirement; in a real system, the limitations of the physical
actuator's frequency response will prohibit large "jerks" anyway.
Some of the disadvantages of Independent Modal Velocity
Feedback are listed below:
1. It is slightly more complicated to implement (i.e., an addition
of a filter is needed in the control loop for each controlled mode).
2. The simulated transient response for independent modal velocity
feedback is not quite as good as direct velocity feedback. However,
as pointed out above in advantage number 3 it is very difficult for
a physical actuator to produce a control input into the structure
when large changes in acceleration are required. Therefore, the
simulated results for direct velocity feedback control (where large
jerks are present) may not be a good representation of a physical
system's actual response.


60
From these findings it can be assumed that the independent
modal control technique as a minimum holds a degree of merit and
deserves further investigation. The remaining chapters of this
thesis will investigate adaptive control techniques to further
improve on the system response when using independent modal control
combined with an adaptive controller.


CHAPTER V
INTRODUCTION OP ADAPTIVE PREDICTION AND CONTROL
The benefits of using Independent Modal Control and Direct
Velocity Feedback Control have been investigated in chapter
However, there has been some very liberal assumptions made in this
investigation. First, It was assumed that the math model was ideal
(no modeling errors of any kind are present). Second, the model was
constant with time. Third, the physical system behaved in a linear
fashion. Finally, it was assumed that no noise or stochastic
disturbances were present.
These assumptions raise some very obvious questions: what
happens if the math model does not exactly match the physical system
(which it never will in a real system); what happens to the control
law if the system changes with time (due to thermal expansion,
structure material aging, additional constructions, etc.); will the
control law remain stable with nonlinearities; and how will the
control system handle stochastic noise (external physical forces,
sensor/actuator noise, etc.).
Since it is impossible to build and test the control design
on a large space structure on the ground due to the presence of
gravity, all these questions raise a degree of uncertainty about any
control scheme a designer may choose. It would be very beneficial to
have a control system that can adjust itself to account for the


62
unknowns and to optimize the performance once unknowns become known
(or at least become predictable). This brings up the concept of
adaptive prediction and adaptive control.
Although the concept of adaptive prediction and control has
been around for several decades, it has only been in the last decade
or so with the introduction of the microchip that the idea has been
practical to physically implement. The fundamental idea with
adaptive prediction is to take the given input and output data up to
the present time and adjust the predictor model parameters (in real
time) so that past predictions closely match the observed data. Then
these new parameters are used to generate future predictions.
Adaptive control can be broken up into two different
classes: indirect or explicit, and direct or implicit. The
fundamental idea with indirect adaptive control is to adjust the
controller based on the adaptively predicted model; that is, the
parameters of the unknown plant are estimated from input/output data
and these estimates are used to generate a feedback control function.
The fundamental concept with the direct adaptive controller
is to parameterize the system directly in terms of the control law
parameters; no explicit plant identification is involved prior to
generating the feedback control signal [7].
The indirect adaptive control problem can be defined as a
two tier process: prediction and control. First, a prediction scheme
must be implemented so that the system's unknown parameters are found
and can be used to predict the system's future response. This can be
done with a variety of techniques:


62
unknowns and to optimize the performance once unknowns become known
(or at least become predictable). This brings up the concept of
adaptive prediction and adaptive control.
Although the concept of adaptive prediction and control has
been around for several decades, it has only been in the last decade
or so with the introduction of the microchip that the idea has been
practical to physically implement. The fundamental idea with
adaptive prediction is to take the given input and output data up to
the present time and adjust the predictor model parameters (in real
time) so that past predictions closely match the observed data. Then
these new parameters are used to generate future predictions.
Adaptive control can be broken up into two different
classes: indirect or explicit, and direct or implicit. The
fundamental idea with indirect adaptive control is to adjust the
controller based on the adaptively predicted model; that is, the
parameters of the unknown plant are estimated from input/output data
and these estimates are used to generate a feedback control function.
The fundamental concept with the direct adaptive controller
is to parameterize the system directly in terms of the control law
parameters; no explicit plant identification is involved prior to
generating the feedback control signal [7].
The indirect adaptive control problem can be defined as a
two tier process: prediction and control. First, a prediction scheme
must be implemented so that the unknown system parameters are found
and can be used to predict the future response of the system. This
can be done with a variety of techniques:


63
1) By using a conventional linear approach; such as the linear
projection algorithm shown Below.
y(t+d) =o:(z1 )y(t) +^(z1)u(t) (5.1)
Where CX(z^ ) and ^(z-') are derived from the systems transfer
function as will be shown in a later section.
2) Designing a restricted complexity predictor:
See Appendix C for an illustration for one possible restricted
complexity predictor algorithm that can be used for modal estimation.
3) Using an adaptive predictor--using either the direct or
indirect approach:
One type of prediction algorithm called a d-step-ahead
projection adaptive prediction is shown below.
y(t+d) = 0(t,z-1)y(t) + j@(t,z-1)u(t) (5.2)
Note how equation (5*2) differs from (5*1) CX(z-^) and (z"are
now functions of time "t".
The second step of the adaptive process is to find a
suitable control scheme using adaptive control techniques. Some of
the popular schemes are
1) d-step-ahead adaptive control.
Where the control input is determined at each point in time so
the system output at a future time instant is brought to a desired
value.
2) Adaptive algorithms for closed loop pole assignments.
Where the poles of the closed loop system are adjusted to suit
a desired system performance.


64
5.1 PARMA MODELING
When working with adaptive prediction and control, a very
common type of system modeling is the deterministic autoregressive
moving average (DARMA) model. The model is expressed as
A(z-1)y(t) =B(z-1)u(t); t >_ 0 (5-1.1)
where
A(z"1) = Ao + A-|z1 + ... + Anz_n; Ao is nonsingular (5-1.2)
B(z"1) = (Bo + B^"1 + + Bmz-ra)z-d (5-1.5)
Note, if u(t) is a random "white noise" process then the model is
called an ARMA model (the deterministic part is taken out).
After comparing equations (5-1.1), (5.1.2), and (5.1*3) it is
clear that the DARMA model for a single input single output (SISO)
system is simply the system transfer function.
To verify the DARMA model for the system used in this thesis,
the transfer function found from the Standford Program "Conditran" is
used to simulate the system using a DARMA routine. The results were
identical to the results from the physical state space and modal
state space forms. The first and second modes were also simulated by
multiplying the system transfer function with the digital filter
transfer function and then again simulating with DARMA. As expected,
the results are identical with previous findings. See Appendix B for
a program listing.


65
5-2 ADAPTIVE PREDICTION
Before getting right into adaptive prediction, some of the
basics of prediction methods are first examined. Starting with the
DARMA model
A(z-1)y(k) = B(z-1)u(k) (5.2.1)
the SISO case can be shown as
A(z-1) = 1 + k\z~^ + k-2z~2 + + Anzn (5*2.2)
B(z-^ ) = z-d(Bo + Biz"^ + + B^z--*-) (5.2.5)
= z-dB'(z-1) (5.2.4)
which is simply the system transfer function.
A conventional linear predictor for the system output at time
(k+d) has the form
y(k+d) =CX(z-1)y(k) + ^(z1)u(k) (5-2.5)
where CX(z-^) = G(z"1) (5.2.6)
^(z-1 ) = F(z-1)B,(z-') (5.2.7)
F(z-1) and G(z-1) are unique polynomials satisfying
1 = F(z-1)A(z"1) + zdG(z-1) (5.2.8)
F(z-^ ) = 1 + F^z1 + + F(j_-|zd+' (5.2.9)
G(z-') = Go + G-|z' + + Gn_izn+1 (5.2.10)
The coefficients for F(z-1) and G(z-1) can be derived as follows:
Fo = 1 (5.2.11)
i-1
pi = -S i = 1.......d 1 (5.2.12)
j=0
d-1
Gj_ = FjAj_+(j_j; i = 0, .., n 1 (5*2.15)
j=0
Note, for MIMO systems F(z1) and G(z1 ) are matrices, however the


66
same basic principles apply, Only the "1" in equations (5-2.8),
(5-2-9)> and (5.2.11) needs to be changed to the identity matrix [i].
For the one-step-ahead predictor (where d = 1), the process
of finding the coefficient for F(z-1) and G(z1) is quite trivial.
The results turn out to be F(z-^ ) = 1, and
G(z~^ ) = zA(z1) (i.e., Gj_ = A-j_+-| for i = 0, ..., n-1)
Thus 0<(z1 ) = zA(z~1) = A) + A2Z-^ + + Anzn+,' (5.2.14)
^3(z1) = z-1B'(z1) = Boz-1 + B1Z-2 + + Bqz"1-1 (5-2.15)
The conventional linear predictor of equation (5.2.5) is
often written in the following form:
y(k+d) =0oTC^(k) (5.2.16)
where (£>(k)^ = [ y(k), ..., y(k-n+l), u(k), ..., U(k-l-d+l)] (5.2.17)
n is the order of the denominator of the transfer function and 1 is
the numerator order.
Q0^ contains the coefficients of C*(z1) and jS(zT^).
0oT = [ C*(z-1) : ^(z_1)] (5-2.18)
Equation (5.2.16) fully describes the conventional linear
predictor. This form works quite well assuming the transfer function
is mathematically modeled accurately and the transfer function is
time invariant. The next logical place to investigate is how to
improve on the predictor by adaptively changing the coefficients of
the predictor "on line" to better match the resulting data when the
assumption of a perfect time invariant model is not valid.
There are two fundamental types of adaptive predictors:
direct and indirect. In the direct method the parameters of the
model are estimated directly. No calculations are required to


67
determine the predictor from the estimated model. For the indirect
method the parameters are first estimated in an arbitrary model for
the system. Then the model is transformed into the required
predictor format. The discussion below should help clarify this
concept.
For the direct adaptive predictor, the output of the
predictor can be described as
y(t + d, Q) = Cjb(t)T Q (t) (5-2.19)
The predictor is really only a special model for the system. A
widely used form to estimate the parameters of the system |0(t)} is
as follows: A * A
0(t) = 0(t-O + M(t-1) (t-d)e(t) (5.2.20)
where M(t-1) is an algorithm gain, (£>(t-d) is the input and output
state variables, and e(t) is the error between the predicted output
and the actual system output.
Taking equation (5.2.20), an adaptive d-step-ahead predictor
based on the conventional linear projection scheme can be written as
follows:
(5.2.21)
6(t) = § (t D c + d)T^(t d)[y(t) - d)T0 (t 1)l
where c is a small constant that is added to prevent the possibility
of division by zero.
The d-step-ahead prediction of y(t) is given by
i T A
y(t) = (p(t d) 0(t d)
(5.2.22)


68
Another more sophisticated type of direct adaptive predictor
is the adaptive d-step-ahead predictor based on least squares. The
algorithm for this predictor is (5-2.23)
0(t) = 0 (t - 1) + P(t - d -1 )<£>(t d)[y(t) -Cjb(t - d) 0 (t - 1)]
1 + Cp(t d)TP(t d - l)<£>(t - d)
P(t - d) = P(t - d - 1) - P(t - d l)c£)(t d)(£) (t - d)TP(t - d - 1)
1 + cjbCt d)TP(t- d - 1)cf)(t - d)
Note that this predictor is analogous to the Kalman filter predictor
applied on a non-stochastic system.
In contrast to the direct predictor, the indirect adaptive
predictor first fits a model to the data and then manipulates the
estimated model into the form of a d-step-ahead predictor.
Starting with0(t) from equations (5-2.22) or (5.2.23),
A A
A(t,z1) and B(t, z-1) are found as follows:
A(t,z"^) = 1 + a-|(t)z' + + an(t)zn
B'(t,z 1) = h(t) + + bm(t)z
m
(5-2.24)
(5.2.25)
Next, P(t,z1) and G(t,z1) is found by solving the following
equation
A .A
1 = P(t, z1 )A(t, z~1 ) + z-<*G(t,z-1)
Lastly, C*(t,z-1) and A(t,z-1) are found by
A *
0(t,z-^) = G(t,z-')
A A
p(t,z1) = p(t,z-1 )B'(t,z_1 )
Thus, the adaptive d-step-ahead predictor is given by
A
y(k+d) = C*(t,z1)y(k) + yS(t,z1)u(k) (5.2.29)
Prom the derivations above, it is easy to see that the direct
(5.2.26)
(5.2.27)
(5.2.28)
method has the advantage of fewer mathematical calculations required.


69
However, an indirect approach has more flexibility and is often
superior to the direct approach in terms of accuracy. For this
thesis, the direct approach is used.
To test the direct adaptive predictor algorithms, simulations
were run using a one-step-ahead predictor (i.e., d=l). An
intentional error was entered into the estimator's initial model
parameters to see if the predictor would converge toward the true
model. The true values for the estimator can be found by using the
following function:
H(z) = C(zl A)"1B (5.2.30)
First, a standard was established by running a predictor with no
adaptive estimation. Then an adaptive predictor based on linear
projection was run. The results were very promising since the error
did indeed converge toward zero. The results were much better than
when using a non-adaptive predictor/estimator.
Next, an adaptive predictor based on least squares was
simulated. This time the results were even better. The adaptive
estimator's error converged toward zero at a faster rate than with
the projection predictor. A least squares estimator is generally
better than a estimator based on projection. Also, it has the
ability to be successfully used on non-deterministic (stochastic)'
systems, whereas the projection estimator tend to be very sensitive
to state variable errors and may not be convergent on stochastic
systems [8]. See Figure 24 for the simulated results. Appendix B
contains the program listings.


Output for Error Output for Error Output for Error
ffiSOLUTC VftlUE fCR M (BROR 70
to idiptive predictioa
feSOLUTE VftlUE FR THE EftftCfi
Miptive Prediction hsed on Projection
2 3 < 5 i 1 ft
Time (Seconds)
Figure 24. Adaptive P-redictor Simulation.


71
A large space structure control designer should be very-
encouraged by the results of the adaptive predictor. The adaptive
predictor actually refined the parameters of the predictor in time
using only the data found from the tip velocity sensor and the known
input data. Therefore, it is not necessary to be overly concerned
about not having a perfect model when generating a control design;
the adaptive predictor will converge towards that actual system in
time (i.e., will minimize the error between the model and the
physical system). It is only necessary to assure that the control
scheme will take advantage of the updated model in time (to be an
adaptive controller).
In the future it may be possible to create an adaptive
controller that will first find an accurate model for itself (either
starting completely from scratch or starting from a rough estimate)
and then self generate (or self tune) a control system that will
optimize the desired system response. This would really be the
"optimal" control scheme for a space structures since the danger of
uncertainties and modeling errors ruining the control system design
would be eliminated (or at least reduced).


CHAPTER VI
D-STEP-AHEAD CONTROL
In researching ways to improve on transient results of the
Direct Velocity Feedback and the Independent Modal Velocity Feedback
controllers, the possibility of using a d-step-ahead controller was
investigated. The logic for choosing the d-step-ahead controller is,
since the system vibrates in a predictable fashion, it may be
possible to control the system with a series of forces timed to
cancel the vibrations. To simplify the computations, the popular
one-step-ahead controller (i.e., d = 1) was chosen.
The control law for a one-step-ahead controller is rather
simple. Starting with the predictor found in the previous chapter
y*(t + d) = C*(z-1)y(t) + j6(z1)u(t); t >_ 0 (6.1)
where y (t + d) is now the desired future output at time (t + d). It
is trivial to derive the following equation for u(t):
u(t) = y*(t + d) CX(z-1)y The d-step-ahead controller is somewhat analogous to a
deadbeat controller. That is, like in deadbeat control, the system
is brought to its nominal condition (zero in this case) by applying
the correct inputs at the given time. One of the problems with
deadbeat control and one-step-ahead control is that the control force
required to bring the system to its desired output y (t + 1) may be
excessive. To reduce the chance of excessively large control input,


73
an alternate form of equation (6.2) called the weighted form is shown
in equation (6.3).
u(t) = /3o[y*(t + d) C*(z~1)y(t) fi'(z~1)u(t -1)] (6.3)
^ X
Where ^ is a chosen weighting constant. It is also possible to
further impose a penalty on large control inputs as necessary to
prevent over stressing the controller (system actuator) by varying
X between a defined set of boundaries.
To test the possibility of using one-step-ahead control on
the chosen system model, a simulation program was written that uses a
DARMA model and a projection estimator to find y(t + 1). The desired
output y*(t + 1) is selected by multiplying y(t + 1) by a chosen
constant. The feedback control input is found from equation (6.2).
The results of this simulation proved somewhat promising.
The system output could be chosen to approach zero at just about any
rate desired, since the control law requires y(t) to track the
j v . v
desired y (t). However, as predicted above, if y (,t + 1 ) is chosen
to approach the nominal zero value too fast, the control input force
required becomes very large. See Appendix B for the simulation
program listings.
Next, one-step-ahead control was applied to control the first
mode independently. To do this the estimator was first multiplied by
the low pass filter transfer function. Then the desired output
y (t + d) is actually the desired value for the first mode of the
system. A DARMA model was used for the system model. This time, the
results did not show as much promise. The overall system was not
stable even though the first mode was dampened.


74
Even though the results above has some limited success, it is
logical to wonder whether the DARMA model represent a true system
model. To answer this question the simulations were again made using
a state space representation of the physical system (as opposed to a
DARMA model). The control input was again found from the one-step-
ahead control law. The results show the system does not reach its
desired nominal steady state value.
Overall, the results of one-step-ahead proved disappointing.
The first problem is the control system is very sensitive to
inaccuracies in system modeling. This is very undesirable for LSS
control. Secondly, even though the control scheme has the desireable
effect of driving the system to a nominal value in a short period in
time, the penalty is in the large control efforts required. The
control efforts required are even more severe than those of direct
velocity feedback control discussed earlier. These large control
efforts can actually excite the system modes.
Based on the results of the computer simulations, as well as
findings from other readings, it is clear that step-ahead control
will not work well on lightly damped systems (such as large space
structures). The control input required did not diminish with time.
Instead, new control effort is required just to dampen the modes
excited from previous control inputs. Also, it is not feasible to
use d-step-ahead control to control the modes independently. Even
though it may be possible to drive a particular mode to its nominal
value, in doing so the other modes are actually excited. Therefore,
the overall system stability is unsatisfactory.


CHAPTER VII
POLE PLACEMENT ADAPTIVE CONTROL
Another popular control scheme is Pole Placement. The basic
idea is to assign the closed loop poles to a desireable location on
the Z-plane. To very briefly go over the theory, start by defining
the vectors L(z"^ ) and P(z-^) (of order n-1) to be
L(z-1) = (z-1) = F(z-1)B'(z-1)} P(z-1) = (z-1) = G(z-1) (7.1)
Notice that these vectors are the same as those used in the
D-Step-Ahead section. If the system is described by a DARMA model
A(z_1)y(t) =B(z"1)u(t) (7-2)
then (neglecting to include all the intermediate level derivations)
the following equation can be obtained:
A(z-1)L(z~1) + B(z-1)P(z1) = A*(z-1) (7.3)
where A*(z1) is the desired closed-loop polynomial and has an order
of (2n-l). Note, A(z-^ ) and B(z^) must be relatively prime
polynomials (their element matrix is nonsingular).
The implication of equation (7*3) is that given A(z~') and
B(z"1 ) a fairly simple algebraic procedure can be used to solve the
set of linear equations to calculate L(z-^ ) and P(z1) that will
result in assigning the closed loop poles to any A*(z-') locations
desired. One special case is when A*(z1) is chosen to be 1 (all
closed loop poles are at the origin); this is the deadbeat case.


76
Figure 25 shows the control scheme.
Figure 25* Pole Placement Control Plant
In equation (7.3), the variables in A(z1) and B(z1) are
stationary (do not change with time). As stated before, a LSS will
have changing parameters, so the constraint of equation (7*3) may not
be practical in the long run. If the an adaptive predictor is added
to the control system to estimate the system parameters, then the
system becomes a self tuning adaptive controller. Equation (7*3) is
now written as
A A A A
A(t,z1)L(t,z-1) + B(t,z-1)p(t,z-1) =A*(z1) (7.4)
A . A
where A(t,z1) and B(t,z"') are time varying parameters found by the
estimator.
The adaptive control scheme is as follows. First, readings
are taken from the output y(t) and the known input u(t). These
variables are used by the adaptive predictor to calculate predicted
output as well as update the estimator model of A(t,z-^) and
A
B(t,z-1). Then, an on-line numerical analysis program (such as a
Gaussian elimination algorithm) will solve the set of linear
equations that will satisfy equation (7.4).


77
7.1 COMPENSATOR DESIGN FOR POLE PLACEMENT CONTROL
As a variation from the strict pole placement method,
consider using a compensator in the control system. Using a
compensator technique is somewhat analogous to the pole placement
method. However, there is a distinct difference. Whereas the pole
placement method will give a direct analytical location for pole
assignment, the compensator design techniques are more of a trial and
error art where the compensator tries to "steer" the root locus plot
onto a suitable pole location.
Staying with the Root Locus techniques, the compensator
design method is as follows. If it is desired to independently
control the 1st mode, than it is necessary that the effects of the
controller on the other dominant modes are eliminated (for this case
modes 1,2, and 3 are considered dominant, mode 4 can be considered as
a residual mode). To do this the compensator zeros are placed very
near the physical system poles (e.g., directly beneath them) for
modes 2 and 3 Likewise the compensator poles should be placed
directly under the 2nd and 3rd modes' zeros. The idea behind this
scheme is that the 2nd and 3rd modes will remain stationary as the
gains are changed since the poles will approach the zeros in the same
relative location. Note, that the system poles and zeros are all on
the unit circle since the physical system has no inherent damping.
Placing the compensator poles and zeros a short distance beneath the
exact system zeros and pole locations, as opposed to directly on top
of them, should eliminate the possibility creating system instability


78
by accidentally placing the compensator zeros a small amount outside
the unit circle.
Now the gain must be found that will give a suitable
transient response. As done previously in the Root Locus section of
this paper, the numerical technique of Bairstow can be used to find
the minimum (or desirable) value for the "r" in the z-plane. In a
sense, the compensator can be thought of as the digital filter used
previously. As before the following characteristic equation is
solved by the numerical Bairstow algorithm.
D(z)*CD(z) + K1*n(z)*CN(z) = 0 (7.1.1)
Where D(z) and N(z) are the physical system denominator and numerator
transfer function parameters respectively; likewise, CD(z) and CN(z)
are the compensator denominator and numerator transfer function
parameters.
To test this theory the compensator for the 1st mode was
created by placing the compensator zeros directly beneath the 2nd and
3rd modes' poles and placing the compensator poles directly beneath
the 2nd and 3rd modes' zeros. The following Compensator was found:
(7.1.2)
C(z) = 1 + .32814z~1 + 1.489163z~2 + .3280576z~5 + .9999z~4
1 - .1946001z-1 + 1.2710072-^ .1944628z^ + .99979z"4
Next, the numerical root locus technique of equation (7.1.1)
was used to find a suitable gain. The results of the program show
that a gain of K-|= 24.6 will produce an optimal response (r = .573)
however smaller gain can be used and still produce quite a good
response. Now the system can be simulated by adding the compensator


79
into the control loop. The system is simulated with a gain of both
4-5 and 10. The diagram below illustrates the control system.
R .£r\ h ^ Physical System v@
i r
K') ^ Uompprisator .
Figure 26. Independent Modal Control Using a Compensator
The results of the simulations were extremely good. By
examining the plots of Figures 27 and 28 on the following pages it is
easy to see that not only did the effects of 1st mode quickly
deteriorate, but the control input was in the form of a sinusoidal
force equation with a frequency equal to the first mode. This is
quite desirable for a couple of reasons: first it shows that the
controller is only controlling the 1st mode as designed; and
secondly, the physical actuator will be more able to create such a
force input requirement (i.e., there are no large jerks required from
the control actuator).


Figure 27. Simulation of Compensator Modal
Control With Gain Equal to 4.5.


TIP VELOCITY vs TIHE
INDEPENDENT MODAL CONTROL USING ft COMPENSATOR
Feedback Gain is 10
81
Figure 28. Simulation of Compensator Modal
Control With Gain Equal to 10.


82
Next a compensator was designed to control the 2nd mode by
canceling the control effects of the 1st and 3rd modes. The
resulting compensator is
C(z) = 1 .92024z-1 + .368992z~2 .9202875z~5 + .999952z~4 (7.1.3)
1 .23808z-1 + .23791 Hz-* .999831z"^
Note that the order of the denominator is less than the numerator.
This is because the 1st mode has only one zero (at z = 1; another
zero is assumed to be at negative infinity). The gain for the
compensator was found as before, and the system was simulated. The
results of the system where not nearly as good for mode 1. This is
because as the gain for the second mode compensator was varied, the
value of "r" for the second mode changed very little. The root locus
plot tended to follow an arc just inside the unit circle. Thus, when
simulated, the 2nd mode's controller had little impact on dampening
the second mode.
To try to fix this situation, a slight modification to the
control algorithm scheme was made. Instead of cancelling the zeros
for modes 1 and 3, cancel the zeros for modes 2 and 3 instead.
Therefore, the root locus plot will act as though the 1st mode's
zeros are the 2nd mode's zeros. The new compensator was calculated
as: (7.1.4)
C(z) = 1 .920259z-1 + .369026z~2 .9202919z"3 + .9999234z4
1 .194600z-1 + 1.27100z2 .1944626z^ + .9997986z4
The characteristic equation was numerically solved as before to find
a suitable gain. Any gain from 4 to 34 appears to give a very good
settling time for this compensator (with 34 being the best). Again,
the system was simulated. The gain was chosen as 8 to reduce the


83
large control force needed if a higher gain was used. The results
of the simulations are very promising. The 2nd mode was indeed
independently controlled with this compensator controller.
Next, "both modes 1 and 2 where independently controlled at
the same time using the compensator independent modal control scheme.
As before, the gain for mode 1 (K1) was chosen as 4.5 and K2 was
chosen as 8. The figure below illustrates the control system.
Figure 29. Compensator Independent Modal Control For Modes 1 and 2
The results of the simulations found in Figure 30 are
excellent. The system control response was almost comparable with
the ideal modal control response shown back in Figure 12. By
properly selecting a good compensator, independent modal control
appears to work extremely well.


Figure 30. Simulation for Independent Modal
Control for the 1st and 2nd Modes.


85
7.2 COMPENSATOR DESIGN WITH A SELF TUNING ADAPTIVE CONTROLLER
Although the compensator design worked quite well for this
known system, could it be used on an system with uncertainties in the
parameters? If an adaptive estimator is built into the system as
shown in Figure 31, then it is proposed that the control system can
self tune itself to account for modeling uncertainties.
Figure 31 Adaptive STR Compensator Control Scheme
The control system works as follows. First the output from
the tip sensor and the known input from the control actuators are fed
into the adaptive estimator. The adaptive estimator in turn
generates the latest estimate for the system Transfer Function (TF)
coefficients using the method of least squares.
Next, the transfer function coefficients are fed into the
digital compensator function generator (DCFG). The DCFG will use a
numerical analysis technique to calculate the poles and zeros of the


86
system transfer function. Then by using the pole and zero cancelling
process described in section 7.1, the DCFG will assign the zeros and
poles to each compensator, calculate the compensator coefficients,
and adaptively update the compensators in the control loop.
Finally, the coefficients of the system transfer function as
well as the digital compensator coefficients are fed into the Optimal
Gain Generator (OGG). The OGG will first multiply the system
transfer function with the compensator transfer functions. The
resulting function can then be transformed into its characteristic
equations as shown in equation (7.1.1). A numerical analysis
iteration algorithm (such as the Newton method) combined with the
root finding method would then be used to find the proportional gain
Kn for each mode n that results in a desired value for "r". Where
"r" is the distance from the origin to the nth root (pole) location
on the z-plane. The desired value for "r" may be its minimum value,
or some other value that will produce a suitable settling time. The
new gain value is in turn fed to the control actuators thus
completing the control loop. It should be found that once the
physical system is adequately known, small variations in the
parameters (due to the inherent changing nature) will not require
large changes in the gain. Therefore, some computer processing time
could be saved by only updating the control gains at certain
intervals, as opposed to at after each sample cycle.
In the following chapter this same self tuning principle will
be used in another variation of a similar control scheme. In
addition, the effects of stochastic noise on the system and how it


87
can be reduced will be reviewed. The findings from the next section
for reducing stochastic noise is directly applicable to the adaptive
compensator control scheme just described.


CHAPTER VIII
AN ADAPTIVE CONTROL SCHEME FOR INDEPENDENT MODAL CONTROL
Another possible adaptive control scheme that takes into
account all the findings and theories described in the first part of
this thesis is shown in Figure 32.
Figure 32. Adaptive STR Control Scheme
The control system works as follows. First the output from
the tip sensor and the known input from the control actuators are fed
into the adaptive estimator. The adaptive estimator in turn


89
generates the latest estimate for the system Transfer Function (TF)
using the method of least squares.
Next, the transfer function coefficients are fed into the
digital filter generator (DFG). The DFG will use a numerical
analysis technique (such as Bairstow's) to calculate the poles and
zeros of the system's transfer function. The digital filter for each
control mode is now created by centering the filter response to the
system poles (modal frequency). A fairly simple computer algorithm
process will generate the digital filter coefficients of the transfer
function. The filter could be a Butterworth, Chebyshev, or any other
recursive filter. Thus an adaptive digital filter is created for
each controlled mode.
Finally, the coefficients of the system transfer function as
well as the digital filters coefficients are fed into the optimal
gain generator (OGG). The OGG will first multiply the system
transfer function with the digital filters transfer functions. The
resulting function can then be transformed into its characteristic
equation. A numerical analysis iteration algorithm (such as the
newton method) combined with the root finding method would then be
used to find the proportional gain Kn for each mode n that results in
the minimum value for "r". Where "r" is the distance from the origin
to the nth root (pole) location on the z-plane. The new gain value
is in turn fed to the control actuators to complete the control loop.
The control system just described will adaptively optimize
the state variable feedback system when proportional feedback is
used. That is, for a proportional feedback system, this technique


90
will produce the optimal transient response in terms of damping or
settling time. Listed below are some additional advantages in using
this control scheme:
1) It is adaptive so the control system will change with time as
the parameters of the physical system change (due to thermal
expansion/contraction, aging, mass changes due to additions on the
LSS, etc.).
2) It uses the simplicity of proportional feedback lawno complex
algorithm is needed to generate the control law.
3) The nature of proportional feedback makes the control system
rather robust.
On the other hand, this control scheme does have some
drawbacks. Mainly, the system is considered deterministic (no
stochastic disturbances are considered in the design). Also, the
transient control response is limited by the nature of proportional
feedback methods. In other words, generally it is not possible too
select the desired pole locations that give the desired transient
response and then find the proportional gain needed. This control
scheme is limited to finding the best possible solution with the
inherent limitations of proportional feedback. Finally, a relatively
large amount of numerical data processing is required. This could
have the undesirable effect of overloading the capabilities of the
onboard computer. It is possible to reduce the data processing
burden by only updating the adaptive filter and control gains at
certain intervals rather than at each sample period.


91
8.1 ADAPTIVE CONTROL WITH STOCHASTIC NOISE
The next logical step is to consider the effect stochastic
noises. The state-space model for the stochastic system can be
written as follows:
x(k+l) = A*x(k) + B*u(k) + v-j (k) (-8.1.1)
Y(k) = C(k) + v 2(k)
Where A is considered to be relatively constant with time and v-| (k)
and V2(k) are the disturbance noise and the measurement noise
respectively. The disturbance noise for a Large Space Structure can
be caused by a number of factors: unmodeled (residual) modes,
nonlinear structure behavior, onboard disturbances (motors, pumps,
etc), and environmental conditions (thermal gradients, solar
pressure, gravity gradients, atmospheric effects, etc) [l]. The
measurement noise can be caused by the limitations in the sensor
design itself.
By looking at equation (8.1.1) it is easy to see that even if
the parameters of the physical system were defined exactly, the
estimator would not realize this and try to adjust the parameters
such that it matched the output results. That is, the parameter
estimator would assume that the noise was caused by a deterministic
inputs and would begin to chase these "ghosts" that just are not
there. Therefore, it is desireable to create a filter that could
filter out the noise and leave only the actual deterministic output
such that the parameter estimator can give the optimal estimate for
the system parameters. To do this, the Kalman filter will be
considered.


Full Text

PAGE 1

AN INVESTIGATION OF CONTROLLING LARGE SPACE STRUCTURES USING ADAPTIVE INDEPENDENT MODAL CONTROL TECHNIQUES by Brian David Hughes B.GeoE., University of Minnesota, 1982 B.EE., Auburn University, 1985 A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Masters of Science Department of Electrical Engineering and Computer Science 1989

PAGE 2

This thesis for the Master of Science degree by Brian David as been approved for the Department of Electrical Engineering and Computer Science by Sanaa A. Azim Edward T. Wall Jan T. Bialasiewicz Date

PAGE 3

Hughes' Brian David (M.S., Electrical Engineering) An of Controlling Large Space Structures Using Adaptive Independent Modal Techniques Thesis directed by Assistant Professor Sanaa Abdel-Azim and reviewed by Professor Edward T. Wall In this paper the basic principles in Large Space Structure (LSS) modeling, estimation, and control are introduced. Several different control techniques (both adaptive and non-adaptive) are investigated and simulated on a computer. The concept of modal control is explored and its advantages and disadvantages are found. New ideas for using adaptive control on both a deterministic and stochastic systems with modal state feedback are presented. The form and content of this abstract are approved. I recommend its publication. Signed Sanaa Abdel-Azim Signed Edwarn. Wall

PAGE 4

CONTENTS Chapter I. INTRODUCTION .. -.. II. SYSTEM MODELING 3 2.1 Modal Theory . . . . . . . . . . . . . . . . . 1 0 2.2 Modal Superpostion Form 15 2.3 Computer Simulation 16 III. DIGITAL FILTERING 20 IV. v. VI. VII. MODAL CONTROL .......................................... 29 / 4.1 Effects of Modeling Errors, Residual Modes, and Spillover . . . . . . . . . . . . . . . . . . . . 37 4.2 Distributed Actuator Control 39 4.3 Passive Control 40 4.4 Using The Classical Techniques of Root Locus 41 4.5. Findings from Direct Velocity Feedback Simulations and Independent Modal Velocity Feedback Simulations 56 INTRODUCTION TO ADAPTIVE PREDICTION AND CONTROL 61 DARMA Modeling 64 5.2 Adaptive Prediction . 65 D-STEP-AHEAD CONTROL ................................... 72 POLE-PLACEMENT ADAPTIVE CONTROL ...................... 75 7.1 Compensator Design for Pole Placement 77 7.2 Compensator Design with a Self Tuning Adaptive Controller . . . . . . . . . . . . . . . . . . . 85

PAGE 5

VIII. AN ADAPTIVE CONTROL SCHEME FOR INDEPENDENT MODAL CONTROL 88 8.1 Adaptive Control With Stochastic Noise 91 8. 2 Kalman Fi 1 te ring .............. -. . .. . . . . . . . . 92 IX. CONCLUSIONS 1 00 REFERENCES 1 02 APPENDIX A. COMPUTER PROGRAMS USED FOR DESIGN 103 B. COMPUTER SIMULATION PROGRAMS ......................... 131 C. EXAMPLE OF A RESTRICTED COMPLEXITY ESTIMATOR .. 178

PAGE 6

TABLES Table 1. Z-Plane Modal Magnitudes Versus Gain 53

PAGE 7

FIGURES Figure A Space Structure 4 2. A Spring .and Damper Model 4 Model for a LSS 6 4. Equivalent Spring Diagram 7 5. Formulas for Lumped Systems 9 6. Computer Simulation in the Physical Matrix Form ,,,,,,,,,, 18 7. Computer Simulation in the Modal Matrix Form . 19 8. Low Pass Filter Simulation 23 g. Band Pass Filter Simulation 26 10. Comparison of Filter Frequency Response .... 28 11. 1st Modal Damped System Simulation ... 31 12. 1st and 2nd Modal Damped System Simulation . 32 13. Proportional Tip Velocity Feedback 41 14. Root-Locus Plot Using Direct Proportional Velocity Feedback .......................................... 44 15. Tip Velocity Feedback with Proportional Gain . 46 16. Modal Tip Velocity Feedback 47 17. Root-Locus with 1st Mode Filtered 49 18. Simulation with 1st Mode Dampened . 51 19. Independent Modal Control for the 1st and 2nd Modes . 52 20. Diagram for the 2nd Mode .......... 53 21. Simulation with 1st and 2nd Modes Dampened ....... 55

PAGE 8

22. Control Input for Direct Velocity Feedback . 57 23. Control Input for Modal Velocity Feedback ...... 58 24. Adaptive Predictor Simulation ..... 70 25. Pole Placement Control Plant 76 26. Independent Modal Control Using a Compensator 79 27. Simulation of Compensator Modal Control With Gain Equal to 4.5 .. 80 28. Simulation of Compensator Modal Control With Gain Equal to 10 ..... 81 29. Compensator Independent Modal Control for Modes 1 and 2 .. 83 30. Simulation for Independent Modal Control for the 1st and 2nd Modes 84 31. Adaptive Self Tuning Regulator Compensator Control Scheme 85 32. Adaptive STR Control Scheme ... 88 33.. Adaptive STR Control with a Kalman Filter . 94 34. Noise Probability Density Function .... 96 35-Self Tuning Parameter Estimator with Stochastic Noise 96 36. Simulation Results of an Adaptive Predictor with Kalman Filtering ..................................... 98

PAGE 9

CHAPTER I INTRODUCTION This Thesis can be divided into two parts. The first part (Chapters 2, 3 and 4) investigates the basics of modal control and Large Space Structure (LSS) modeling. The second part (Chapters 5, 6, 7, and 8) explores adaptive estimation and control strategies for both deterministic and stochastic systems. In the first part of the thesis some of the principles of modal control are discussed. This is not an easy task; modal control theory covers a plethora of areas: system mathematical modeling, vibration mechanics, structural mechanics, dynamics, modern control, numerical analysis, estimation and optimization theory, filtering, etc [1]. Therefore, several different practices and principles need to be explored. First, the basics of system modeling are briefly examined. Next, some of the necessary mathematical manipulations on matrices are reviewed. Then some of the background in digital filtering is given. This paper will use digital control; digital filtering will be used as a way to separate the different modes. A technique to optimize the gains for different modes when using independent modal feedback control is presented. This technique uses the classical Root Locus method combined with modern numerical analysis. It is shown that this technique works quite well for finding a suitable or optimal gain.

PAGE 10

2 Finally, a number of computer simulations are made to verify the system model and prove the control system results. Both the real physical model (i.e., real within the framework of modeling theory) and an ideal modal model are simulated. Much of the work in this thesis was in developing computer programs. The programs were used to simulate the system, to make numerical calculations (such as finding roots of polynomials), and to design the control system. The second part of the thesis explores different techniques to adaptively self tune and control a large space structure. First, the fundamentals of prediction and adaptive prediction are reviewed. Then these fundamentals are carried forward and used to generate self tuning adaptive control systems. A couple of different adaptive control schemes are The last part of the thesis deals with stochastic control using adaptive techniques. In a real space system, it would be expected that there are many areas where noise becomes a factor in controlling a large space structure. Therefore, this thesis will investigate ways to reduce the effects of noise disturbances in estimating and controlling the vibrational dynamics of the system.

PAGE 11

CHAPTER II SYSTEM MODELING In modeling theory, a physical system can be represented by using a mathematical model. The mathematical model is the set of algebraic, vector, and differential equations which describe the behavior of the system when various inputs and design parameters are varied. In creating the model for this thesis, the principle of "lumping" is used. This principle will model a distributed mass system into individual nodal points. Each node contains a fictional "point mass" which is approximately equal to one half of the mass per unit length times the length between the different nodes. The elastic properties of the physical system can be modeled into equivalent spring and dampers (i.e., the distance between the nodes can be modeled with massless springs and dampers). The problem is then reduced to that of a simple 2nd order mass-spring-damper model. In modeling theory, just about any linear physical system can be represented by a system of equivalent mass-spring-dampers. Often times with very accurate results. In lumped system models, like finite element models, the basic theory is to represent the physical system by dividing it into a number of small segments. Each segment is mathematically modeled by a single lumped mass and the spring and damper equivalent. A basic rule is the more segments (or nodes) the model is divided into, the more closely the model represents the true system. There are

PAGE 12

4 limitations however: each node modeled into the system creates another order in the model. The computations become excessive with higher order models (especially if the computer used is not a main frame). For this thesis, the system model is lumped into 4 nodes. Figure 1 illustrates a such physical system. Figure 1. A LSS Assume a series of trusses are attached to a large space structure (LSS); the LSS is assumed as rigid for modeling purposes. For this model the requirement that only the very tip of the truss will contain the sensors and control actuators is imposed. Some of the reasons for this requirement could be weight constraints, manufacturing constraints (i.e., the truss is assembled in space unit by unit from the end), etc. A very simplified model of the structure using springs, masses and dampers is illustrated in Figure 2. Figure 2. Spring and damper model

PAGE 13

5 A large space structure will be of lightweight material such as aluminum. Aluminum has a modules of elasticity E = 2.5 x 1011 PA (N/m2). Aluminum structures typically have very low inherent damping (around 0.2%) so the damping coefficients C1 C4 can be assumed to be equal to zero for simplification purposes. Note that in making the assumption that Cn is zero, we are actually modeling the worst case scenario where no inherent damping is present. Therefore, a control design with this assumption will be slightly conservative. By assuming this undamped case, the model is thus reduced into a simple mass-spring equation which is easily transformed into its mathematical equivalent using basic 2nd order differential equations. These differential equations can be arranged into the matrix form below. [m]X(t) + [k]X(t) [b]f(t) ( 2. 1 ) Where [m] is the mass matrix, [k] is the stiffness matrix, [b] is the force distribution matrix, X(t) is the displacement at time t, and f(t) is the force at time t. Equation (2.1) will be referred to throughout this report as the physical form of the model . By defining the state vector V = [X X ]T, equation (2.1) can be transformed into its state-space equivalent form V = [A]V + [B]u where u is the input vector of forces, [A] is the state variable and [B] is the input distribution matrix. (2.2)

PAGE 14

6 Figure 3 illustrates the axial vibration model for the undamped (free body) case. K2 K3 K4 +f(t) m2 m3 m41 -f( t) 11 X1 12 X2 13 X3 14 X4 Figure 3. Example of a model for an undamped LSS Transforming the axial vibration model into a set of second order differential equations. K2(x2 x1) + K3(x3 x2) K3(x3 x2) + K4(x4 x3) (2.3) Equation (2.3) can be placed into system of equation (2.2) as follows: 0 m2 0 0 given set so 0 g }t lK1 + K2) -K2 0 -K2 (K2 + K3) m3 + 0 -K3 0 m4 0 0 [m]x(t) + [K]x(t) V = [ ;; x ]T v G [b]f(t) K4(x4 x3) -f(t) the state matrix equivalent 0 -K3 (K3 + K4) -K4 =[g}t K4 -1 (2.4) (2.5) (2.6)

PAGE 15

7 To find the values for the stiffness matrix [K], the equivalent spring values, for each node, need to be determined. The equivalent spring for a uniform bar in elastic tension-compression is shown below. Figure 4. Equivalent spring diagram p A*E k (2.7) L Where A is the cross sectional area E is the modules of Elasticity L is the length of the bar Note, if complete analysis on the model is required, the equivalent mass and springs constants in the bending and torsional directions would also need to be analyzed. However, for this thesis just the axial modes (modes in the tension-compression direction) will be controlled. The bending and torsional vibration modes can be modeled into their equivalent mass-spring-damper system in an analogous fashion. Therefore, any of the principles applied to the axial vibration model can likewise be applied for torsional and bending modes.

PAGE 16

For this model the stiffness values for K are derived as follows. E*A 8 (2.8) Where ki is the equivalent spring constant for the ith node of the model. Ai is the equivalent area for the ith node. li -li-1 is the length between the ith node and the i-1 node To model the individual k values, a cross sectional area of 3 x 1o-5 m2 was chosen, and each nodal segment was selected to be 10 meters in length. The value for E was chosen as 2.5 x 108 N/m2. This is 1000 times less than the true value for aluminum. This was done for computational purposes only (i.e., to achieve lower frequency modes than would be expected for a real physical system). The IBM-PC computer used had some problems in finding the eigenvalues for the higher frequency modes due to limits in the machine. The resulting values for the k's for each node or segment was 750 N/m2. Figure 5 on the following page contains a summary of the formulas used for calculating the lumped mass into nodes [2].

PAGE 17

"' man "' dilctibutcd n1111 "' lumrcd mau ur cquioaknl rcm ,,.; .. c c ... cilc-r r'l system ""' "'J "'t t .l'l Equivalent "' sysrem '" "' + "' "' + 1'" Simrlt t .... .., Fiac4 t...m 1'1 Ill "'c f t r:-'1 1'1 1/2 f "' t "' "' "' + 111 ,;. '" + 0.31 Figure 5. Formulas for lumped systems Source: Shigley, J.E., Mechanical Engineering Design, McGraw-Hill, New Jersey, 1977. The approximate distributed mass of each section is md = p *A*L = 1 Kg; where p is the density of around 3000 Kg/m3. The formulas above were used with many liberal approximations in 9 order to create the model. Specifics are not considered, since this thesis concentrates on controlling a space system, not modeling it. Thus, the lumped mass "m" for the equivalent modeled system was selected as 1 Kg for the first three nodes. The mass at the 4th node (the tip) will have a lumped mass of m = 1.5 Kg. This represents the addition of a concentrated mass "m0 of the control actuator and sensor as well as the distributed mass. The physical state-space matrix system in the form of equation (2.2) with the numerical values for all k's and m's is as follows.

PAGE 18

0 0 0 0 0 0 0 0 V(t) = 1 0 0 0 0 1 0 0 . V = [ X where 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 -1500 750 0 0 0 0 0 0 750 -1500 750 0 0 0 o 0 0 750 -1500 500 0 0 0 0 JT [ .. x ]T. X and V = X 0 0 750 -500 0 0 0 0 V(t) + 0 0 0 0 .66667 0 0 0 In order to simulate this system on a computer, it is converted into its discrete equivalent form. This was done by running a version of the Standford Program "CD.BAS" which was 10 f(t) modified to make it more user friendly; it is called "CD2.BAS". This program uses the series expansion method to convert the matrices into their discrete equivalent. The sampling rate was chosen as .05 seconds (20Hz). As will be shown in later sections, this is a good choice for sampling rate since it is greater than twice the highest modeled frequency of the structure. Note, by Shannon's sampling theorem, for a mode to be controllable it must be sampled at a rate equal or greater than twice the mode's frequency. See Appendix A for the listing of the program. 2. 1 MODAL THEORY Given the second order physical system model of a linear structure [m]X + [c]X + [k]X = [b]u (2.1 .1) The structural eigenvalue problem associated. with the X solutions, for the undamped free vibration case of [m]x + [k]x = o (2.1.2)

PAGE 19

11 is det[ mA2 -k] = 0, eigenvalues: {A_i2 -u.Ji2, i = 1,2, n} (2.1.3) and k J
PAGE 20

12 A [ [o] [r] (2.1.8) To convert the physical system state-space form into the modal state-space form, a satisfactory solution for the eigenvector matrix [T] must first be found. To find [T], a version of the Stanford Program "MODALCOMP" called "MC2.BAS" was used which again was modified to make more user friendly. The program uses the numerical analysis technique of Leverrier, Bairstow, and Gauss to calculate the eigenvalues, the eigenvector matrix [T], the inverse eigenvector matrix [T]-1, and the modal dynamics matrix [T]-1[F][T]. Where [F] is simply the dynamics matrix of the state space form: x = [F]x + [G]u Note, to make this program work, [F] must be set as m-1k; then the problem will be in its equivalent state space form (i.e., the eigenvalues found solve the equation -m-1k] = 0 ). See Appendix A for a program listing and a computer printout of the results. From the results of the program, the Modal state-space matrix system can be constructed. The program gives the results for the upper right partition of dynamics matrix [A] in equation (2.1.8). The Modal state-space input distribution [B] is calculated as follows: f .08833] .18642 -.27711 .36828 (2.1.9)

PAGE 21

13 The modal state-space matrix system for the model is shown below in equation (2.1 .10). (2.1.10) 0 0 0 0 -2610.06 0 0 0 .08833 0 0 0 0 0 -1658.20 0 0 .18642 0 0 0 0 0 0 -657.63 0 -.27711 0 0 0 0 0 0 0 -74.111 .36828 t) = 1 0 0 0 0 0 0 0 (t) + 0 U(t) 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 0 0 Note X TTJ [;] and{= (2.1.11) X TTJ As before, the state-space matrix was converted into its discrete time equivalent using "CD2.BAS". The results can be found in Appendix A. The eigenvector matrix found by the modified Standford program is not unique. In fact there is an infinite number of eigenvector matrices that can exist and still be correct; any eigenvector may be scaled by some constant: Ur = Cr c;f; r (2.1.12) Where Cr is a scaling constant whose units are such that cprTmgbr has the dimension of mass, and Ur is the scaled equivalent eigenvector. It is well known that the eigenvectors are orthogonal with respect to m and k. In other words where s f r. (2.1.13) These eigenvectors can be normalized by using the orthogonality conditions (sometime called orthonormal relations) to scale the eigenvector matrix into its normalized form. This is quite simply

PAGE 22

14 done by scaling the rth mode so that the generalized mass, or modal mass, defined by Mr = UrTmur Mr crcfJrTmcf;>rcr has a specified value of 1. In other words, the modal mass (2.1.14) (2.1.15) [M] =[I]. It isn't hard by substituting equation (2.1.12) into (2.1 .14) to derive the following equation for the scaling constant for each eigenvector: ci2 = m11 cpa 2 + m22cp2i 2 + + mnncp ni 2 (2.1.16) For this system the scaling constants were calculated to be c1 =.680713, c2 = .672263, c3 = .612563, and c4 = .606862. By multiplying each eigenvector by its scaling constant the normalized eigenvector (modal) matrix is calculated as [-.459917 .672263 -545389 .19171 J .680713 -.141807 .612563 .375884 Tn = --547586 -.642347 .142617 .516913 .129758 .277302 --452384 .606862 In the normalized modal matrix form, the following properties are true. [Tn]T[m] [Tn]T[m][Tn] [Tn]T[k][Tn] From equation (2.1 .17), T -1 n [-.459917 .672262 .545390 .197712 .680712 -.141805 .612563 .375882 -.547588 -.642346 .142617 .516913 .194636] .415953 -.678576 .606862 (2.1.17) (2.1.18) (2.L19)

PAGE 23

15 Using the normalized equations, the state-space modal equations are as follows. (t) or ( t) = [-[Tn]T[c][Tn] [I] = [-[Tn]T[c][Tn] [I] -[ wo2J J (t) + [[Tn]T[bj u(t) 0 + [ [Tn]T[b]J u(t) 0 since there is no inherent damping [c] in the model, the equations can be simplified to [ 0 -[ w 2)] + [ [Tn)T[b)J { ( t) u(t) [I] 0 0 2.2 MODE SUPERPOSITION FORM Another technique for evaluating the modal vibrations in a structure is called the modal superposition form. This technique uses mode shapes that are superimposed on each other to define the deflections in the structure. The deflections for the structure is given by 1 y(t,x) = L k = 0 (2.2.1) where t is time, x is position location, 1 is upper limit of modes (4 for the model system), is the kth mode shape at position x, and lJk(t) is the kth modal amplitude as a function of time. If there are point actuators located at positions x1 ,x4 then the modal amplitudes each satisfy are in equation (2.2.2)

PAGE 24

16 Vk( t) + 2 C kfhvk( t) + Sl'-k Vk( t) )4 : u1 ( t)\h(x1), (2.2.2) k 1 4 = 1 where Jlk and t k are the kth modal frequency and damping coefficients respectively; ui(t) is the ith actuator force on the structure [3]. 2.3 COMPUTER SIMULATION To simulate the mathematical models, programs were written in BASIC and Pascal (see Appendix B for the listing of the simulation programs). Both the modal state matrix form and the physical state matrix form of the physical system are simulated. The data gathered is the velocity and displacement of the tip of the structure with respect to time. To generate the vibrations in the structure, a unit impulse input "u" of 10 newtons applied at the tip is simulated. This impulse force could be representative of a force caused by firing of retro rockets on the space shuttle, docking of another space vehicle, etc. As mentioned previously, the mathematical models are first converted into their discrete equivalent form. This was done by using the program "CD2.BAS", which is a modified version of the Standford Program called "Condistran". Figures 6 and 7 illustrate the graphical representation of tip velocity versus time, and tip displacement versus time, for the

PAGE 25

17 physical state matrix form and the modal state matrix form of the system model respectively. As was expected, the results for the two models are virtually identical.

PAGE 26

,3 -.. -.3 .B25 T v -.925 _.._ I VJ I v TIP QLOCITY vs TIH Jj "-1 TIMJ Cseco fl.P DISPLACEXOO 11s TIHE A \ I v v rtKE Cseconas> Figure 6. Comp_uter Sirn_u1atiop in the Physical Matrix Form. 18 \ I

PAGE 27

19 .3 TIP UELOCITY vs 11Kt r--\ I I \ tl /j N II \ \( t \ I,J v -.3 "-TIKl (seco as> .e2s TIP vs !IKE _.. Nl 1\ \ \ 'I . . . \J v v v -.e2S LTIKE (seconas> Figure 7. Simulation in the Modal Matrix Form.

PAGE 28

CHAPTER III DIGITAL FILTERING In order to use modal control, it is necessary to separate the system dynamic vibration response into its separate modes. To do this, causal digital filtering techniques are used. Causal digital techniques can be broken down into two basic forms: nonrecursive and recursive. Butterworth recursive filters are used for this thesis. Recursive filtering has the advantage of being faster to implement in real time. The basic transfer function equation for the recursive filter is shown below. Y(z) bo + b1z-1 + ... + bqz-q Hlp(z) = ( 3.1 ) X(z) + a1z-1 + ... apz-P Where Hlp(z) designates a low pass filter. A digital filter equation such as (3.1) is sometimes modified into the following equivalent form: H(z) B ( 1 + b 1 z -1 + + bq z-q ) 1 + a1z-1 + + apz-P (3.2) For the physical system, the structure will be controlled by implementing sensor and controllers (actuators) on the end or tip of the structure. Therefore, only the velocity value at the tip needs to be filtered. The digital filter can be implemented on a computer by the following simple algorithm formula,

PAGE 29

q Y(n) = k=O p L Ak*Y(n-k) k=1 21 (3.3) where Y(n) is the output at time n, X(n) is the input state at time n, Bk is the transfer functions numerator coefficients, and Ak is the denominator coefficients. The order of the numerator is q and the denominator order is p. In this case, Y(n) will be the modal velocity of the desired mode at the tip. X(n) will be the actual velocity of the tip. The first mode (lowest frequency mode) can be filtered by using a low pass filter (LPF). All other modes must use a band pass filter (BPF). In order to determine the filter transfer function equation, a program was written that calculates the poles of the Butterworth in both the Z and S domain, as well as the required gain. See Appendix A for a program listing. The zeros of a Butterworth filter are all at (z -1) by definition. For the LPF design, a bandwidth of around 2.25 Hz is desired. This is slightly greater than the 1st modal frequency as will be shown in the next chapter. The program was run and it gave the results for the gain and denominator coefficients of the digital filter. The numerator coefficients of the filter were found by applying the following simple equation: 8 coefficient bk is L k=O (n)! (n -k) !k! The results are summarized as follows: (n is the order 8) (3.4)

PAGE 30

22 numerator coefficients denominator coefficients bo = 1 b1 = 8 a1 = -4.383951 b2 = 28 a2 8.920154 b3 = 56 a3 =-10.821606 b4 70 a4 8.489805 b5 = 56 a5 -4.385521 b6 = 28 a6 1 .450788 b7 = 8 a7 -.280147 bs = 1 as .024117 Gain is B = 5.328643 x 1o-5 To compare the results of the filtered modes to that of the ideal modes, the modal equation form of the ideal output of the sensors can be found by using the equation below. Note: {is "s'' "i" Thus for a N qi(t) = j=1 = N j=1 the eigenvector modal matrix the modal velocity vector subscript denotes the mode. subscript denotes the node (nodal position .. is s chosen mode s qs(t) = (3.5) of the sensor). (3.6) is the predicted "ideal" output of the sensor at position "i" of the system. Programs were written that simulate both the ideal sensor output using the modal state matrix system along with equation (3.6), and the output of the physical sensor found by digital filtering the physical matrix system. Figure 8 shows the response of the digital filter, the ideal filter response, and the difference (error) between them respectively.

PAGE 31

m r.ocAL urtocltY vs tiME msr PoOH I DUL FILTtR f\ r I I I l v v v v v v v v .3 TIME .3 .,....,..-TIP UtLOCitY vs TIKt 1st (fiLTtRtD UltH 8th I I \) I I v v v v v v v v v -.3 -'-TIX . 2 TIP VELOCITY tRSOR Yl 1st MODE I ABSIUtiH Ut u I .2 TIKI lstcondsl Figure B. Low Pass Filter Simulation. 23

PAGE 32

24 To create a BPF, a LPF Butterworth is designed with the desired bandwidth, then the following transformation property is used to shift the frequency res.ponse. (3.7) Thus, the order of the BPF will be twice that of the LPF. The basic equation for a BPF is shown below. + (3.8) + For the BPF design, a LPF with the appropriate bandwidth was first designed by using the same filter program as before. Then the BPF was calculated by utilizing equation (3.7) and putting it into the form below, (3.9) q IT k=1 (1 zkejWcz-1) q T1 k=1 (1 zke-jGUc z-1) B*---------------------+ B IT!=1 (1 -pmejWcz-1) n P m=1 ( 1 Pme-jW<:: z-1) where Zk and Pm are the LPF's zeros and poles respectively and CVcis the center frequency for the BPF. Using equation (3.9) to calculate the transfer function, the resulting BPF is as follows. numerator coefficients denominator coefficients bo 1 b1 .3456895 a1 = -2.728852 b2 -49.36748 a2 = 6.753185 b3 -11.23136 a3 = -10.36797 b4 286.1777 a4 = 14.35491 b5 39.3709 a5 = -14.79166 b6 -420.0358 a6 = 13.79967 b7 = -30.828 a7 = -9.915372 bs = 190.56 as = 6.447651 b9 = 5.39381 a9 = -3.104981 b10= -21.87708 a10= 1.348126 b11= -.102038 a11= -.3582316 b12= .2951724 a12= .0871268 Gain is B = 2 X 8.576548 X 10-6 1.71531 X 10-5

PAGE 33

As before, computer simulations were made of the band pass filter output, the ideal modal filter output, and the difference (error) between them. Figure 9 on the following page shows the results. 25

PAGE 34

,3 .3 TIXE (seconds! .3 ::. """":"'" .. TIP KODAL utLOCITY vs liKE Znd KODE FILTRD IIITH A 12 ORDR BUTTRIIORTII Brr liKE Cmondsl ,) TIP KODnL VELOCitY ERROR vs liKE lsi KODt I ABS(Vti It Vidrall I .3 tl HE C uconds I Figure 9. Band Pass Filter Simulation 26

PAGE 35

27 By examining both the 1st mode's LPF and the 2nd mode's BPF difference plots, it is obvious that the digital filter output differs quite a bit from the ideal filter at the start of the simulation. In the simulations, the output of the digital filter was different than the output of the ideal filter for about one second, after this the outputs of the two filters match (relatively speaking). It appears as though the filter undergoes a "settling time" in which the filter becomes equivalent to the ideal filter. Part of the reason for this is the filter coefficients do not contain an initial value. Therefore, the filter will need as many sample cycles as the order of the filter before it gives truly useful data. The intuitive way to reduce this settling time is to increase the sampling rate. This would lessen the time it takes for all of the coefficients of the filter to contain an initial value and for the filter to produce legitimate output information. Unfortunately, this method is not without penalties. If the sampling rate is increased, the spacing between the modes is decreased in the digital frequency plot. As the spacing between the modes is decreased, the digital filter design will need to be of higher order to prevent overlapping (spillover) of the modes. A higher order filter will need more samples cycles for all its coefficients to contain initial values; thus, it is a no win situation. Figure 10 illustrates this point.

PAGE 36

0 0 MAGNITUDE VERSUS FREQUENCY fs = 20 Hz Order Order \ \ Mode 1 Mode 2 Mode 3 frequency MAGNITUDE VERSUS FREQUENCY fs = 40 Hz 3 Mode 4 Mode t R frequency 4 Figure 10. Comparison of Filter Frequency Response 28 Pi Pi

PAGE 37

CHAPTER IV MODAL CONTROL As was shown in the section on Modal Theory, the modal state matrix form is a diagonal dominant matrix system. The upper left partition of matrix [A] contains the damping terms. The damping form can be put into its 2nd order equivalent in the usual way: damping term = -2 t k W k ( 4. 1 ) where t k and W k is the damping coefficient and frequency for the kth order mode. For the 1st mode, GU1 was shown to be 8.609 rad/sec. If the damping coefficient is chosen as = .7, then the damping term is calculated by using equations (4.1) as -12.0526. The modal domain state matrix system is now 0 0 0 0 -2610.06 0 0 0 0 0 0 0 0 -1658.2 0 0 0 0 0 0 0 0 -657.6 0 0 0 0 -12.0526 0 0 0 -74.1 Ac 1 0 0 0 0 0 0 0 (4.2) 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 The modal form input distribution matrix B is the same as in the Modal Theory section. If the input is set as u = -GC{ where G is the gain matrix and C is the output distribution matrix, then the closed loop form of the modal state space equation can be written as follows: = .+ B*fd(t) Y(t) = C{(t) (4.3)

PAGE 38

where fd(t) is the disturbance force, Ac = [ A-GC ], and G 0 0 0 12.0526 0 0 0 0 c 30 (4.4) [ 0 0 0 1 0 0 0 0 ]. (4.5). The same process can be used to simulate damping of the 1st and 2nd modes together. The 2nd mode frequency is UU2 = 25.644 rad/sec. From equation (4.1) 2 = .1, thus = -35.9016. The modal state space form for the closed loop matrix Ac is now as follows (note, again matrix B does not change): 0 0 0 0 -2610.06 0 0 0 0 0 0 0 0 -1658.2 0 0 0 0 -35.902 0 0 0 -657.6 0 0 0 0 -12.0526 0 0 0 -74.1 Ac 1 0 0 0 0 0 0 0 (4.6) 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 The modally damped systems are now simulated. As before the systems are first transformed into their discrete equivalent forms using the Standford program Condistran. The following pages show the plot of the "Ideal Modal Response" at the tip of the system with the 1st mode dampening coefficient of .70; and also with both the 1st and 2nd modes dampened with a coefficient of .10. See Appendix B for a listing of the simulation programs.

PAGE 39

.3 IDEAL J\ ,, I } v -.3 H TIP vs TIME r 11 l II I I v v TJMr f,./1\,.,.. ... J_\ Figure 11. 1st Modal Damped System Simulation. f I w

PAGE 40

.3 -.3 TIP UELOCITY vs TIME IDEAL MODAL SOLUTION 1st AND 2nd HODES DAMPENED@ .? TIME (seconds) Figure 12. 1st and 2nd Modal Damped System Simulation. w N

PAGE 41

:;:; Unfortunately, even though modal control is fairly easy to understand and very easy to show mathematically, it is very difficult to implement physically. This is because it requires distributed controllers, and the best that can be done is to simulate the distributed controller by a series of discrete point actuators. For the model, suppose it is desired that the physical model is dampened sothe response is the same as the modal model with the 1st mode dampened at 0.7. The calculations are as follows: First, the diagonal elements of the modal matrix are set equal to the required modal gain matrix. diagonal { L i = 1, I 4 (4.7) Then the physical gain matrix is found by use of the following

PAGE 42

34 The same basic equations can be carried out in the normalized eigenvalue form. For the normalized form the modal gain is Kmnorm = 1 '4 (4.1 0) so Kmnorm Km Thus the normalized physical gain matrix is Kpnorm = (4.11) Note: Kpnorm = [Tn]-T[Km][Tn]-1 (normalized physical gain matrix) (modal gain matrix) By taking advantages of some of the properties of normalized eigenvector matrices, the following equalities can be found. Kpnorm = [Tn][m][Km][Tn]-1 [Tn][m][Km][Tn]T[m] The normalized physical gain matrix is easily calculated. -.471137 -.8957 -1.232 -2.1692 -.89571 -1.7029 -2.3418 -4.1240 Kpnorm -1.23177 -2.3418 -3.2205 -5.6713 -2.16918 -4.1240 -5.6713 -9.9872 Note that the normalized physical gain matrix is equal to the (4.12) previously found gain matrix except for the multiplication of m. (4.13) Unfortunately, to implement these physical gain matrices (Kp and Kpnorm) the physical system would need to have an actuator and a sensor on each node. As previously stated, for this physical structure model only the tip of the structure contains the sensor and actuators. Now the question is, is this system modal-state controllable by the use of a single actuator location with a direct

PAGE 43

35 proportional feedback gain matrix Kp. That is, could the physical system response be the same as the ideal modal response with the single tip actuator using proportional velocity feedback. The answer is quite obviously no. The physical form of the gain matrix (Kp or Kpnorm) could not possibly be implemented with a single actuator located only at the tip (node 4). This is also obvious by looking at the feedback fb(k) = B*C*gain. Since the input distribution matrix B [ 0 0 0 1 0 0 0 0 ]T, and the observer distribution matrix C [ 0 0 0 1 0 0 0 0 ] for the first mode, only the gain matrix coefficient k44 can contain a non-zero value as shown below. 0 0 0 0 0 0 0 0 B*C*gain 0 0 0 0 0 (4.14) 0 0 0 k44 0 0 So the question is, how can k44 be chosen such that the physical system will most closely match the modal controlled system. The PACOSS project designers used the following formula to generate the feedback gain to achieve the desired active damping in the specific modes [4]. K na x nc (4.15) where Toe is the normalized modal matrix of observed and controlled modes, and is the desired active damping ratio of the controlled modes. Note that the formulation of equation (4.15) assumes the number of actuators, na, is equal to the number of controlled modes, nc. If na is greater than nc, the gain matrix is overspecified and a

PAGE 44

36 least squares (minimum norm) solution may be used. If the na is less than nc, the gain matrix is underspecified and a pseudoinverse minimum norm approach would have to beused. The system presented in this thesis has more controlled modes than actuators since actuator control is limited to the of the structure; so this "underspecified" system must use a pseudoinverse approach. The problem with the pseudoinverse approach is that the system will show a loss of closed loop damping ratios from the design values when this approach is used. In a later section it will be shown that this is certainly true for this system. Even though the system cannot be modal-state controllable by using a single point actuator, there are some possible advantages of trying to dampen the structure by independent modes. For a large space structure certain modes may not be important and need not be considered. If they are not important, no control energy should be wasted by controlling them. Therefore, some of the possible independent modal solutions methods mentioned in this section can still hold merit for this single point actuator problem. In section 4.4 another method will be considered that is different than equation (4.15). This method will use the classical approach of the root locus technique to achieve the best possible results in terms of transient responses. It will be shown that this method has some benefits over equation (4.15) for an underspecified case, since the transient response can easily be optimized.

PAGE 45

37 4.1 EFFECTS OF MODELING ERRORS, RESIDUAL MODES, AND SPILLOVER When trying to use modal control on a structure, the uncertainties of modelling errors, residual modes and spillover effects must be considered. For this system a 4th order model was used to mathematically represent the physical system. Thus, there are only 4 modes in the model. In an actual structure there are a large number of modes. The effect of these unmodeled modes, called residual modes, caused by this modeling truncation must be examined. When the sensor outputs are contaminated by the residual modes through the term C*V(k), where V(k) is the velocity output of the sensor, this is called observation spillover. When the feedback control excites the residual modes through the term B*f(k), where f(k) is the feedback control force, this is called control The effects of the spillover can be summarized as follows. When only control spillover is present, the effect is to excite the residual modes at their natural frequencies. When both observer and control spillover is present, the results are much more complicated and can possibly lead to instabilities in the system [5]. For the system used in this thesis, only the 1st and 2nd modes will be controlled. The 3rd and 4th modes can be thought of as the residual modes. So far in the discussions the actuators and sensors have been considered collocated. The question of the effect of this constraint on the control system design must be considered. For a real physical system collocated sensor and actuator pairs are more desirable than non-collocated sensors and actuators due to their robustness. If it

PAGE 46

38 is not possible to collocate the sensors and actuators due to some physical design constraint, then there is a greater potential for generating difficult control problems. For example, if there is modeling error or model truncation there is a possibility that the control scheme for the non-collocated system could produce an unstable system. On the other hand, a collocated sensor/actuator system (assuming ideal sensors and actuators) will always produce a positive definite system, regardless of modeling errors [4]. Therefore, the collocated system will be asymptotically stable by the definition of Liapunov. For this paper the false assumption that the actuators (and sensors) are ideal and have unity gain is intentionally made. In reality, a physical actuator transfer function is usually non-linear but can be approximated with a linear 1st, 2nd, or higher order transfer function (TF). A simple 1st order TF for an actuator can be as follows: TF(s) = ( tconst*s + 1)-1 (4.1.1) Note that this TF tends to create a phase lag into the system. Therefore, it should be kept in mind that the assumptions of ideal actuators and sensors are not totally valid in the real sense. Nevertheless, these assumptions will be used in this thesis for analytical purposes.

PAGE 47

39 4.2 DISTRIBUTED ACTUATOR CONTROL There has been work done recently in developing a distributed actuator for controlling flexible beams. The current design is to use a piezoelectric polymer film that can be polarized, or made piezo-electrically active [6]. Voltages or fields applied across the film results in a longitudinal strain. The field can be varied spatially; likewise the strain will vary spatially. It is possible that a control system could be designed where the applied field could be spatially related to the mode shape so.the feedback will be applied in the form of a distributed mode shape. If there is an accurate knowledge of the mode shapeg for a given system, then this system should work. Plus there would be several advantages: theoretically, only one sensor is needed to control the modes (each modal frequency would need to be filtered), the piezoelectric film should be much more reliable than the modern day actuator, and the control system would be inherently robust due to the basic physics of the system (no moving mechanical parts for the controller) Thus, in the future it is quite possible that a true distributed control load could be applied, and modal control could be implemented directly. This may only be feasible for simple linear beam type structures without complex mode shapes and would work best for the bending modes; however, theoretically it is possible that it could work for the axial or rotational modes also.

PAGE 48

40 4.3 PASSIVE CONTROL So far this thesis has only mentioned the principles of actively controlling large flexible structures. In reality, as much design time, if not more, will be spent on considering ways to use passive damping in the LSS design. This is due to the much greater reliability of a passive dampened structure versus an actively dampened structure. Every structure will possess some inherent damping. For an aluminum structure, an inherent dampening ratio of 0.2% is often assumed. There is a lot of research going on in developing new materials in the form of composites that can be used to increase this damping ratio. In addition, there are several studies and articles written on how to find the optimal design of a structure to maximize the passive damping as well as find the actuator locations on the structure that will produce maximum active dampening effects. By using passive dampening together with active control, several advantages can be achieved: a reduction in the energy requirements needed to control the structure, significant reductions in maximum actuator torque requirements, less complex active control systems involving fewer actuators and sensors, and weight savings due to smaller and fewer actuators needed. All these benefits lead to more robust and reliable systems which will more efficiently accomplish their mission goals [4].

PAGE 49

4.4 USING THE CLASSICAL TECHNIQUES OF ROOT LOCUS In Chapter 3, it was shown how the vibration dynamics of a physical system could be modeled by digital filtering each of the dominant modes. It may be desireable to control the system by independently damping the dominant modes. This requires finding a suitable gain for each of the modes. 41 One simple way to find the suitable modal gains is to use the classical Root-Locus technique. The Bairstow numerical algorithm can be used in conjunction with the Root-Locus to numerically find the exact poles at certain gains. To illustrate this technique, the direct velocity feedback method with a proportional feedback gain is first considered. The feedback will be generated by velocity readings taken from the sensor on the tip of the structure. The readings are then fed to the control actuators. The control actuators are also located on the tip of the structure. This is often referred to as collocated sensor-actuator control. The illustration below shows the problem. R + u Physical system v Figure 13. Proportional Tip Velocity Feedback Where R is the disturbance input, V@4 is the velocity at the tip (node 4), and K is the gain.

PAGE 50

42 To find a suitable gain, the basic principal of Root-Locus was used. Instead of using the graphical techniques, the computer's ability to mathematically find the poles using the numerical analysis method called Bairstow's algorithm is utilized. Using this algorithm, a program was written that will find the zeros (closed loop poles) of the function [ D(z) + K*N(z) ] = 0 (4.4.1) which is the denominator of the closed loop transfer function. That is G(z) N(z) DTZT where D(z) is the denominator coefficients and N(z) is the numerator coefficients of the open loop system transfer function G(z) respectively; K is the proportional feedback gain. By using the program, the Root-Locus plot can be directly created by entering different values for K and receiving the pole locations for the modes at that K. The system transfer function G(z) was obtained by using the modified version of the Stanford Program "Condistran". The program uses the recursive Leverrier algorithm for finding the transfer function denominator and numerator polynomial coefficients. To determine the best (optimum) value for gain K that the fastest transient settling time, the formulas listed in equations (4.4.3) through (4.4.6) can be used. Note that in the Z-domain all points inside the unit circle are stable. The closer to the origin these points lie, the faster the settling time for the system. If the poles (roots of the equations) are plotted in the Z-plane, the polar results are as follows:

PAGE 51

the polar form z = the damping ratio t the natural frequency n = the time constant T = where t is the sampling period. -ln r + 82 1 /fV1n2 r + 82 1/Cw = -t/ln r 43 (4.4-3) (4.4.4) (4.4.5) (4.4.6) Obviously, from equation (4.4.6), to get the best transient response the time constant, and thus the magnitude of r, must be minimized. Figure 14 shows the root locus plot as well as the numerical values obtained for different gains in the z-domain. The real and imaginary coordinates of the poles of the characteristic equation are printed for the different values of the gain K. These coordinates can easily be changed into the polar form of equation (4.4.3), and the time constant for each gain K can be calculated.

PAGE 52

I o o I I I lea oo IOo I .... I........ I 4. :::::!:! ..... 0 t 0 0 mode 4 .' ; . . .. : :: : : : :: : .. ... "0 ...... ,, ::::: ::: ... --+i ;_; :.;_i : :-: -: :-+-; ;--; :-: ,;...;: :::::::::I .. . . .. ...... ::L: \ 0. 0 ------1--. i .. : : : : : : ;.; : : : : : : ':; : : : ... :: '!:: : : : :. : 0 0 ........ ... I :::,: ......... . -..... . ...... .. 0 0 -:-: r.. 0 ,-:-: _;.. +. II o oiO o II II 00 I 0111 I o I 0 o o I I o e o I 1 o o 0 .. . . ... ........ 0 .. . . ... 0 . . .. .. . . . ... . . .. .. :: :!::::: ... :: :; ......... II Ill oae 0 I: o I oo eo ... . ... ....... .., : 0 :: : :: : : :.: : : : : : :: .;.:.:.:-.f.:.' ::.;":..::.:,;..:..:. .. :::J,:... :..:..::..:..:..::..:..:.! -----. ... ...... ... . .. ...... 0 : :: : : :.:! : : : : : :: 0 0 ... . . . . ...... ....... 0 ......... . . .. . . . . .. .. .. .... . . .. . 0 o I o 0 I o ...... --1--_ .. _. ....;..j...:. .;._ . 0 .. . ..... . . . . .... 0 . . .. . . .. ... .. . .. . ... . ... 0 Figure 14. Root-Locus Plot Us.ing Direct .Proportional Velocity .Feedback. 44 K = 10 K = 32.

PAGE 53

45 From Figure 14 it is evident that a proportional gain of K= 32 is a very good choice to dampen the vibration modes (especially for modes 1 and 2). Unfortunately, due to physical constraints, a gain of 32 may be more than a real actuator could produce. This is one of the reasons why modal control will be important for future large space structures. With modal control, only certain modes will need to be dampened (usually the lowest order modes). No effort will be wasted on damping modes that are not dominant and do not have a large effect on the overall structure. As will be shown later, when modal control is used, the gains will be noticeably smaller (albeit, the transient vibration damping will not be quite as good). To simulate the system using direct proportional velocity feedback, the discrete form of the state space matrix system is used where so X(k+1) u(k) [Ad]*X(k) + [Bd]*R(k) Kp*[Bd]*V(k) [ R(k) Kp*V(k) ]1 (4.4-7) [Ad] and [Bd] are the discrete forms for the state space matrices, Kp is the proportional feedback gain, and V@4(k) is the sensor velocity readings at the tip of the structure. Note that V@4(k) is the 4th element of vector X(k) (i.e., x = [ x1 x2 x 3 V X5 x6 X7 Xa ]T ). Figure 15 shows the transient vibration behavior of the proportional tip velocity feedback problem. The system was again initially agitated by applying a 10 newton impulse on the tip of the structure. The top graph shows the response using the close to optimal gain value of 32. The bottom graph shows the response of a gain value of 10 for comparison purposes.

PAGE 54

.3 -.3 .3 -.3 Figure 15. TIP UELOClTY vs TIHE DIRECT FEEDBACX WITH GAIN H : 32 TIME (seconds) TIP UELOCITY TIME DIRECT UELOCITY FEEDBACK WITH GAIN K : 19 TIME (seconds) Tip Velocity Feedback with Proportional Gain 46

PAGE 55

47 For modal control, the individual modal velocity can be found by filtering. Then the mode can be independently controlled by multiplying the modal velocity by a suitable feedback gain. For the first mode, a low pass filter can be used. Again the tip of the structure is used and the velocity sensor and control actuators are collocated. The illustration below shows the system. u Physical system Figure 16. Modal Tip Velocity Feedback Here, K1 designates a gain for the first mode. This time the characteristic equation to be solved is [ D(z)*DF(z) + K1*N(z)*NF(z) ] = O. (4.4.8) Where DF(z) and NF(z) are the low pass filter transfer function denominator and numerator coefficients respectively. It is expected that since only the first mode is being controlled, the other three modes will not vary a lot as the gain value is changed. By looking at the actual plot on the following

PAGE 56

48 page, this appears to be the case. Also note that a gain of around K1 = 4.5 gives the best transient time response. Figure 17 on the following page show the results of the Z-plane poles for different K1 values on the Root-Locus plot.

PAGE 57

t. ... ......... ...... .. . . . . .... .. Z-Plane :::: ::!:: It I ::::::::1 t I' .......... ......... .. . ... . . .... . . ......... t . ..... . ...... . . .... . .. I I I ell . . . . . . . . ........ . . . I a a ' I ...... :;::::t:: '''' ... ...... -. I a a I I a ....... . . . .. .. ::::::::: i::::f:;: ... ... . . . .. :::::. !.: .: .. : .:.1::::: ... ... \ ...... . . . . . . 1 1 : 1 : : ; 1 1 : : 1 : I : :: : : li : . : : : : = : : : : : : : : : : = . . . . . ............. .... :::::: ;_:,.:..;.:-:.:. .L I .. :.:..:. :::::::: ...... I I I I . . . . . ... ...... .. . ... . . .. ..... ::::::::: :::.::::: ............ . . . . ........ ..:. .. ,. ................. : : : : : : : :. : : : : :: : ........ ,.. . . . . ...... J .: : : : : : a a ... a a a I I I I I a I I I
PAGE 58

50 To simulate this system, where the 1st mode is dampened using actuator control at the tip of the structure, the discrete state space matrix system is shown below. note: X(k+1) = [Ad]*X(k) + [Bd]*R(k)-K1*[Bd]*V1(k) u(k) = [ R(k) K1*V1 (k) ] (4.4.8) K1 is the 1st modal feedback gain (chosen as 4.5), and V1(k) is the 1st mode's velocity (the modal velocity) at time pulse (k) at the tip of the structure. As before, the system is agitated with a 10 newton impulse force. Figure 18 contains two plots: the top plot shows the tip velocity versus time; the bottom plot shows the tip velocity for the first mode only. See the Appendix B for listings of the simulation programs used to develop the data for the plots.

PAGE 59

.3 TIP VELOCITY vs TIHE 51 -r--FEEDBACH CAIH H1 : 4.5 I r I ,. II lj I( I -.3 I-TIME (seconds) .3 TIP MODAL UELOCIIVvs TIME FEEDBACK GAIN Kl : 4.5 -.3 TIME (seconds) Figure 18. Simulation with 1st Mode Dampened.

PAGE 60

52 After examining the top plot of Figure 18, it is obvious that the system vibration response will not be suitable with only one mode controlled. Therefore, the next mode of highest importance must be looked at. In most cases, this will be the next mode (mode 2). As stated previously, this mode and all other higher modes will need to be filtered with a band pass filter. The control scheme will now include an additional BPF and another feedback loop. Figure .19 illustrates the system. u Physical system v Figure 19. Independent Modal Control For The 1st and 2nd Modes Since the root locus technique only works when one gain is changed at a time, only the root locus plot using the feedback loop for the second mode is generated. As was shown earlier, when changing the gain for a particular mode, the other modes are not effected to a great extent. They tend to be mutually exclusive. Therefore, with narrow band filtering it should be'possible to generate a root locus plot for one mode, find the suitable gain.

PAGE 61

53 then use another independent root locus plot for a different mode; and find a suitable gain for this mode. Figure 20 illustrates the system used to generate the root locus plot for the second mode. R+ u Physical system ... -, r K2 -BPFL -, I ... Figure 20. Diagram for the 2nd Mode Table 1 shows the numerical results for the Root-Locus plot. The table shows how the magnitude of the z-plane poles change as the gain K2 was changed. As predicted, by changing the gain K2, only mode 2 was effected; that is, the other modes were not effected in a relative sense. Gain K 0 1 2 3 4 5 5.7 6 7 10 15 r for mode 1 1 .99999 .99999 1 1 .99999 .99999 .99999 .99999 .99999 r for mode 2 1 .99594 99161 .98708 .98262 .97945 .97874 .97889 .98085 .99210 1.01002 (unstable!) r for mode 3 Table 1. Z-plane Modal Magnitude versus Gain K2

PAGE 62

54 Figure 21 show the transient vibration behavior of the system using modal control to dampen the 1st and 2nd mode at the tip of the structure. The first modal feedback gain K1 is chosen as 4.5, the second modal feedback gain K2 is chosen as 5.70. As before, the system is agitated with a 10 newton impulse force. To simulate this system, the discrete state space matrix system is shown below: (4.4.9) X(k+1) = [Ad]*X(k) + [Bd]*R(k)-K1*(Bd]*V1(k)-K2*[Bd]*V2(k) where K1 and K2 are the 1st and 2nd modal feedback gains, and V1(k) and V2(k) are the first and second modes' velocities (modal velocities) at time pulse (k) at the tip of the structure. Note: u(k) = [ R(k) K1*V1(k) K2*V2(k) ]. Two plots are shown in Figure 21. The top plot shows the tip velocity versus time; the bottom plot shows the tip velocity for the second only. The tip velocity plot for the first mode is virtually identical to what was found in Figure 18. See Appendix B for listings of the simulation programs used to develop the data for the plots.

PAGE 63

.3 -.3 .3 -,3 Figure 21. TIP vs TIME FEEDBACK GAINS Kl: 4.5 I K2: 5.7 TIME (seconds) TIP MODAL vs TIME SECOND MODE FEEDBACK CAINS K1: 4.5 & X2: 5.7 TIME (seconds) Simulation with 1st and 2nd Modes Dampened. 55

PAGE 64

4.5 FINDINGS FROM DIRECT VELOCITY FEEDBACK SIMULATIONS AND INDEPENDENT MODAL VELOCITY FEEDBACK SIMULATIONS So far this thesis has looked into the theories and simulated results of using independent modal control and direct velocity feedback. Before proceeding into the subject of adaptive 56 control, it might be useful to take a further look into what can be found from the results of the simulations. Figures 22 and 23 contain plots of the control input required for the direct velocity methods and the control input required for the independent modal control technique. The top plot of Figure 22 shows the control input for direct velocity feedback simulation with a gain of 10 The bottom plot shows the control input for direct velocity feedback simulation with a gain of 32; note that the scale for this particular graph is 8 to -8 Newtons, all the other plots of Figures 22 and 23 have a scale of 4 to -4 Newtons. The plots of Figure 23 show the control input for independent modal control. Both the first and second modal control inputs are plotted. After reviewing the control input plots as well as the system response plots from previous sections, some obvious findings and conclusions can be made about the usefulness of independent modal control over direct velocity feedback.

PAGE 65

CONTROL INPUT vs liME DIRECT VELOCITY FEEDBACK GAlli K : lB -4 TIME (seconds) 8 COHTROL INPUT vs TIME DIRECT FEEDBACK CAIU K : 32 -e TIME
PAGE 66

4 -4 4 -4 Figure CONTROL INPUT vs TIHE HODAL UELOCI tY F-EEDBACX 2na HODE CAIN K2 : 5.7 TIME (seconas) CONTROL INPUT vs TIME HODAL VELOCITY FEEDBACK 1st MODE WITH CAIH Kl : 4,5 TIHE (seconcts) Control Iriput' for Modal Velocity Feedback. 58

PAGE 67

59 Some of the advantages of Independent Modal Velocity Feedback over that of direct velocity feedback are listed below: 1. No control energy is wasted on unimportant modes. That is, the modes that do not have a large effect on overall structure stability will not be included in the control loop. 2. The maximum control force requirement is much less. 3. The rate of change in control force required is not as large and is smooth in nature. In other words, the rate of change in acceleration (sometimes called a jerk) required of the control actuator is not as overwhelming as it is for direct velocity feedback control. This will allow for a less rigid (beefed up) actuator requirement; in a real system, the limitations of the physical actuator's frequency response will prohibit large "jerks" anyway. Some of the disadvantages of Independent Modal Velocity Feedback are listed below: 1. It is slightly more complicated to implement (i.e., an addition of a filter is needed in the control loop for each controlled mode). 2. The simulated transient response for independent modal velocity feedback is not quite as good as direct velocity feedback. However, as pointed out above in advantage number 3, it is very difficult for a physical actuator to produce a control input into the structure when large changes in acceleration are required. Therefore, the simulated results for direct velocity feedback control (where large jerks are present) may not be a good representation of a physical system's actual response.

PAGE 68

60 From these findings it can be assumed that the independent modal control technique as a minimum holds a degree of merit and deserves further investigation. The remaining chapters of this thesis will investigate adaptive control techniques to further improve on the system response when using independent modal control combined with an adaptive controller.

PAGE 69

CHAPTER V INTRODUCTION OF ADAPTIVE PREDICTION AND CONTROL The benefits of using Independent Modal Control and Direct Velocity Feedback Control have been investigated in chapter 4. However, there has been some very liberal assumptions made in this investigation. First, It was assumed that the math model was ideal (no modeling errors of any kind are present). Second, the model was constant with time. Third, the physical system behaved in a linear fashion. Finally, it was assumed no noise or stochastic disturbances were present. These assumptions raise some very obvious questions: what happens if the math model does not exactly match the physical system (which it never will in a real system); what happens to the control law if the system changes with time (due to thermal expansion, structure material aging, additional constructions, etc.); will the control law remain stable with nonlinearities; and how will the control system handle stochastic noise (external physical forces, sensor/actuator noise, etc.). Since it is impossible to build and test the control design on a large space structure on the ground due to the presence of gravity, all these questions raise a degree of uncertainty about any control scheme a designer may choose. It would be very beneficial to have a control system that can adjust itself to account for the

PAGE 70

62 unknowns and to optimize the performance once unknowns become known (or at least become predictable). This brings up the concept of adaptive prediction and adaptive control. Although the concept of adaptive prediction and control has been around for several decades, it has only been in the last decade or so with the introduction of the microchip that the idea has been practical to physically implement. The fundamental idea with adaptive prediction is to take the given input and output data up to the present time and adjust the predictor model parameters (in real time) so that past predictions closely match the observed data. Then these new parameters are used to generate future predictions. Adaptive control can be broken up into two different classes: indirect or explicit, and direct or implicit. The fundamental idea with indirect adaptive control is to adjust the controller based on the adaptively predicted model; that is, the parameters of the unknown plant are estimated from input/output data and these estimates are used to generate a feedback control function. The fundamental concept with the direct adaptive controller is to parameterize the system directly in terms of the control law parameters; no explicit plant identification is involved prior to generating the feedback control signal [7]. The indirect adaptive control problem can be defined as a two tier process: prediction and control. First, a prediction scheme must be implemented so that the system's unknown parameters are found and can be. used to predict the system's future response. This can be done with a variety of techniques:

PAGE 71

62 unknowns and to optimize the performance once unknowns become known (or at least become predictable). This brings up the concept of adaptive prediction and adaptive control. Although the concept of adaptive prediction and control has been around for several decades, it has only been in the last decade or so with the introduction of the microchip that the idea has been practical to physically implement. The fundamental idea with adaptive prediction is to take the given input and output data up to the present time and adjust the predictor model parameters (in real time) so that past predictions closely match the observed data. Then these new parameters are used to generate future predictions. Adaptive control can be broken up into two different classes: indirect or explicit, and direct or implicit. The fundamental idea with indirect adaptive control is to adjust the controller based on the adaptively predicted model; that is, the parameters of the unknown plant are estimated from input/output data and these estimates are used to generate a feedback control function. The fundamental concept with the direct adaptive controller is to parameterize the system directly in terms of the control law parameters; no explicit plant identification is involved prior to generating the feedback control signal [7]. The indirect adaptive control problem can be defined as a two tier process: prediction and control. First, a prediction scheme must be implemented so that the unknown system parameters are found and can be used to predict the future response of the system. This can be done wi th a variety of techniques:

PAGE 72

63 1) By using a conventional linear approach; such as the linear projection algorithm shown below. y(t+d) =Ql(z-1)y(t) + {3(z-1)u(t) Where CX(z-1) and {J(z-1) are derived from the systems transfer function as will be shown in a later section. 2) Designing a restricted complexity predictor: (5 .1 ) See Appendix C for an illustration for one possible restricted complexity predictor algorithm that can be used for modal estimation. 3) Using an adaptive predictor--using either the direct or indirect approach: One type of prediction algorithm called a d-step-ahead projection adaptive prediction is shown below. y(t+d) = CX(t,z-1)y(t) + {3(t,z-1)u(t) (5.2) Note how equation (5.2) differs from (5.1 )-C:X(z-1) and {3 (z-1) are now functions of time "t". The second step of the adaptive process is to find a suitable control scheme using adaptive control techniques. Some of the popular schemes are 1) d-step-ahead adaptive control. Where the control input is determined at each point in time so the system output at a future time instant is brought to a desired value. 2) Adaptive algorithms for closed loop pole assignments. Where the poles of the closed loop system are adjusted to suit a desired system performance.

PAGE 73

5.1 DARMA MODELING When working with adaptive prediction and control, a very common type of system modeling is the deterministic autoregressive moving average (DARMA) model. The model is expressed as 64 A(z-1 )y(t) = B(z-1)u(t); t > 0 (5.1.1) where A(z-1) B(z-1 ) Ao + A1z-1 + + Anz-n; Ao is nonsingular (Bo + B1z-1 + + Bmz-m)z-d (5.1.2) (5.1.3) Note, if u(t) is a random "white noise" process then the model is called an ARMA model (the deterministic part is taken out). After comparing equations (5.1.1), (5.1.2), and (5.1.3) it is clear that the DARMA model for a single input single output (SISO) system is simply the system transfer function. To verify the DARMA model for the system used in this thesis, the transfer function found from the Standford Program "Conditran" is used to simulate the system using a DARMA routine. The results were identical to the results from the physical state space and modal state space forms. The first and second modes were also simulated by multiplying the system transfer function with the digital filter transfer function and then again simulating with DARMA. As expected, the results are identical with previous findings. See Appendix B for a program listing.

PAGE 74

65 5.2 ADAPTIVE PREDICTION Before getting right into adaptive prediction, some of the basics of prediction methods are first examined. Starting with the DARMA model A(z-1)y(k) B(z-1)u(k) the SISO case can be shown as A(z-1) B( z-1 ) z-d(Bo + B1z-1 + + Blz-1) z-dB' (z-1) which is simply the system transfer function. (5.2.1) (5.2.2) (5.2.3) (5.2.4) A conventional linear predictor for the system output at time (k+d) has the form y(k+d) = CX(z-1)y(k) + {3 (z-1 )u(k) where CX(z-1) G(z-1) {3(z-1) F(z-1)B'(z-1) F(z-1) and G(z-1) are unique polynomials satisfying F(z-1) G( z-1 ) 1 + F1z-1 + + Fd-1z-d+1 Go+ G1z-1 + + Gn-1z-n+1 (5.2.5) (5.2.6) (5.2.7) (5.2.8) (5.2.9) (5.2.10) The coefficients for F(z-1) and G(z-1) can be derived as follows: Fo = (5.2.11) i-1 Fi -2: FjAi-j; i = 1 ... d -1 j=O (5.2.12) d-1 Gi = -L; FjAi+d-j; i = 0, ... n j=O (5.2.13) Note, for MIMO systems F(z-1) and G(z-1) are matrices, however the

PAGE 75

66 same basic principles apply, Only the "1" in equations (5.2.8), (5.2.9), and (5.2.11) needs to be changed to the identity matrix [I]. For the one-step-ahead predictor (where d = 1), the process of finding the coefficient for F(z-1) and G(z-1) is quite trivial. The results turn out to be F(z-1) = 1, and G( z-1 ) = zA(z-1) (i.e., Gi = Ai+1 for i = 0, ... n-1) Thus CX(z-1) zA(z-1) = A1 + A2z-1 + ... + A z-n+1 (5.2.14) n (3 (z-1) z-1B'(z-1) = Boz-1 + B1z-2 + ... + Blz-1-1 (5.2.15) The conventional linear predictor of equation (5.2.5) is often written in the following form: y(k+d) = BoTcp (k) (5.2.16) where = [ y(k), y(k-n+1), u(k), U(k-l-d+1)] (5.2.17) n is the order of the denominator of the transfer function and 1 is the numerator order. tjoT contains the coefficients of CX(z-1) and J3Cz-1). tj0T = [ O<(z-1) : {3Cz-1) J (5.2.18) Equation (5.2.16) fully describes the conventional linear predictor. This form works quite well assuming the transfer function is mathematically modeled accurately and the transfer function is time invariant. The next logical place to investigate is how to improve on the predictor by adaptively changing the coefficients of the predictor "on line" to better match the resulting data when the assumption of a perfect time invariant model is not valid. There are two fundamental types of adaptive predictors: direct and indirect. In the direct method the parameters of the model are estimated directly. No calculations are required to

PAGE 76

67 determine the predictor from the estimated model. For the indirect method the parameters are first estimated in an arbitrary model for the system. Then the model is transformed into the required predictor format. The discussion below should help clarify this concept. For the direct adaptive predictor, the output of the predictor can be described as ,. y(t + d, 9) = cpCt)T e (t) (5.2.19) The predictor is really only a special model for the system. A widely used form to estimate the parameters of the system {fj(t)} is as follows: ,. ,. ,. e c t) 8 (t-1) + M(t-1) cp (t-d)e(t) (5.2.20) where M(t-1) is an algorithm gain, is the input and output state variables, and e(t) is the error between the predicted output and the actual system output. Taking equation (5.2.20), an adaptive d-step-ahead predictor based on the conventional linear projection scheme can be written as follows: (5.2.21) $(t -d) T ect) = e ct-1) c +(t d)Tcpct d)[y(t) -cpctd) 8 ct-1)] where c is a small constant that is added to prevent the possibility of division by zero. The d-step-ahead prediction of y(t) is given by ,. T ,. y(t)=cpCt-d) BCt-d) (5.2.22)

PAGE 77

68 Another more sophisticated type of direct adaptive predictor is the adaptive d-step-ahead predictor based on least squares. The algorithm for this predictor is (5.2.23) Note that this predictor is analogous to the Kalman filter predictor applied on a non-stochastic system. In contrast to the direct predictor, the indirect adaptive predictor first fits a model to the data and then the estimated model into the form of a d-step-ahead predictor. Starting withB(t) from equations (5.2.22) or (5.2.23), A A(t,z-1) and B(t,z-1) are found as follows: A A(t,z-1) = 1 + a1 (t)z-1 + + an(t)z-n A A B'(t,z-1) = bO(t) + + bm(t)z-m Next, F(t,z-1) and G(t,z-1) is found by solving the following equation "' 1 = F(t,z-1)A(t,z-1) + z-da(t,z-1) Lastly,CX(t,z-1) and f3Ct,z-1) are found by A CX(t,z-1) = G(t,z-1) .. {:3 ( t, z-1 ) = F ( t, z-1 ) B' ( t, z-1 ) Thus, the adaptive d-step-ahead predictor is given by y(k+d) = CX(t,z-1)y(k) + )9Ct,z-1)u(k) (5.2.24) (5.2.25) (5.2.26) (5.2.27) (5.2.28) (5.2.29) From the derivations above; it is easy to see that the direct method has the advantage of fewer mathematical calculations required.

PAGE 78

However, an indirect approach has more flexibility and is often superior to the direct approach in terms of accuracy. For this thesis, the direct approach is used. 69 To test the direct adaptive predictor algorithms, simulations were run using a one-step-ahead predictor (i.e., d = 1). An intentional error was entered into the estimator's initial model parameters to see if the predictor would converge toward the true model. The true values for the estimator can be found by using the following function: H(z) = C(zi -A)-1B First, a standard was established by running a predictor with no adaptive estimation. Then an adaptive predictor based on linear projection was run. The results were very promising since the error did indeed converge toward zero. The results were much better than when using a non-adaptive predictor/estimator. Next, an adaptive predictor based on least squares was simulated. This time the results were even better. The adaptive estimator's error converged toward zero at a faster rate than with the projection predictor. A least squares estimator is generally better than a estimator based on projection. Also, it has the ability to be successfully used on non-deterministic (stochastic) systems, whereas the projection estimator tend to be very sensitive to state variable errors and may not be convergent on stochastic systems [s]. See Figure 24 for the simulated results. Appendix B contains the program listings.

PAGE 79

U25 1-1 1-1 e.o2 1-1 9.015 I 5. Ul 'l ;::l o U:JS T 'I' fe.SOlUH YAlU[ ICR Ttl[ [RRCI! tdapUV! prtdicllcn lai "I' 'I' 3 e "I' "I"' "I' 5 1 a Time (Seconds) 1'1![ M1pUve PredicliCI'I 8md en 2 3 5 ' P!rolutr YRlU[ 1'1![ lust Sq.rlres fl:l1pth1 PrtdicUcn Time (Seconds) 9.03r------------------, 0 1-1 1-1 O.e2 1-1 0 'H 0.815 ;::l p. ;::l 0 Ul i 'c r 9 .. ......... 1 3 4 56 '8 Time (Seconds) Figure 24. Adaptive P-redictor Simulation. 70

PAGE 80

71 A large space structure control designer should be very encouraged by the results of the adaptive predictor. The adaptive predictor actually refined the parameters of the predictor in time using only the data found from the tip velocity sensor and the known input data. Therefore, it is not necessary to be overly concerned about not having a perfect model when generating a control design; the adaptive predictor will converge towards that actual system in time (i.e., will minimize the error between the model and the physical system). It is only necessary to assure that the control scheme will take advantage of the updated model in time (to be an adaptive controller). In the future it may be possible to create an adaptive controller that will first find an accurate model for itself (either starting completely from scratch or starting from a rough estimate) and then self generate (or self tune) a control system that will optimize the desired system response. This would really be the "optimal" control scheme for a space structures since the danger of uncertainties and modeling errors ruining the control system design would be eliminated (or at least reduced).

PAGE 81

CHAPTER VI D-STEP-AHEAD CONTROL In researching ways to improve on transient results of the Direct Velocity Feedback and the Independent Modal Velocity Feedback controllers, the possibility of using a d-step-ahead controller was investigated. The logic for choosing the d-step-ahead controller is, since the system vibrates in a predictable fashion, it may be possible to control the system with a series of forces timed to cancel the vibrations. To simplify the computations, the popular one-step-ahead controller (i.e., d = 1) was chosen. The control law for a one-step-ahead controller is rather simple. Starting with the predictor found in the previous chapter y*(t +d)= + J3Cz-1)u(t); t > 0 ( 6 .1 ) where y (t +d) is now the desired future at time (t +d). It is trivial to derive the following equation for u(t): u(t) y*(t +d)CX(z-1)y-(t) @'(z-1)u(t -1) (6.2) {3o The d-step-ahead controller is somewhat analogous to a deadbeat controller. That is, like in deadbeat control, the system is brought to its nominal condition (zero in this case) by applying the correct inputs at the given time. One of the problems with deadbeat control and one-step-ahead control is that the control force required to bring the system to its desired output y*(t + 1) may be excessive. To reduce the chance of excessively large control input,

PAGE 82

73 an alternate form of equation (6.2) called the weighted form is shown in equation (6.3). u(t) = @o[y*(t + d)
PAGE 83

74 Even though the results above has some limited success, it is logical to wonder whether the DARMA model represent a true system model. To answer this question the simulations were again made using a state space representation of the physical system (as opposed to a DARMA model). The control input was again found from the one-stepahead control law. The results show the system does not reach its desired nominal steady state value. Overall, the results of one-step-ahead proved disappointing. The first problem is the control system is very sensitive to inaccuracies in system modeling. This is very undesirable for LSS control. Secondly, even though the control scheme has the desireable effect of driving the system to a nominal value in a short period in time, the penalty is in the large control efforts required. The control efforts required are even more severe than those of direct velocity feedback control discussed earlier. These large control efforts can actually excite the system modes. Based on the results of the computer simulations, as well as findings from other readings, it is clear that step-ahead control will not work well on lightly damped systems (such as large space structures). The control input required did not diminish with time. Instead, new control effort is required just to dampen the modes excited from previous control inputs. Also, it is not feasible to use d-step-ahead control to control the modes independently. Even though it may be possible to drive a particular mode to its nominal value, in doing so the other modes are actually excited. Therefore, the overall system stability is unsatisfactory.

PAGE 84

CHAPTER VII POLE PLACEMENT ADAPTIVE CONTROL Another popular control scheme is Pole Placement. The basic idea is to assign the closed loop poles to a desireable location on the Z-plane. To very briefly go over the theory, start by defining the vectors L(z-1) and P(z-1) (of order n-1) to be L(z-1) (z-1) = F(z-1)B'(z-1); G(z-1) (7. 1 ) Notice that these vectors are the same as those used in the D-Step-Ahead section. If the system is described by a DARMA model A(z-1)y(t) B(z-1)u(t) (7.2) then (neglecting to include all the intermediate level derivations) the following equation can be obtained: (7-3) where A*(z-1) is the desired closed-loop polynomial and has an order of (2n-1). Note, A(z-1) and B(z-1) must be relatively prime polynomials (their element matrix is nonsingular). The implication of equation (7.3) is that given A(z-1) and B(z-1) a fairly simple algebraic procedure can be used to solve the set of linear equations to calculate L(z-1) and P(z-1) that will result in assigning the closed loop poles to any A*(z-1) locations desired. One special case is when A*(z-1) is chosen to be 1 (all closed loop poles are at the origin); this is the deadbeat case.

PAGE 85

76 Figure 25 shows the control scheme. PLANT P(z-1) 1 {u(t)} ... 1 B(z-1) y L(z-1) A(z-1) r (t) I Figure 25. Pole Placement Control Plant In equation (7.3), the variables in A(z-1) and B(z-1) are stationary (do not change with time). As stated before, a LSS will have changing parameters, so the constraint of equation (7.3) may not be practical in the long run. If the an adaptive predictor is added to the control system to estimate the system parameters, then the system becomes a self tuning adaptive controller. Equation (7.3) is now written as A( t, z-1 )L( t, z-1) + B( t, z-1 )P( t, z-1) = A*(z-1) (7.4) ,... ,... where A(t,z-1 ) and B(t,z-1) are time varying parameters found by the estimator. The adaptive control scheme is as follows. First, readings are taken from the output y(t) and the known input u(t). These variables are used by the adaptive predictor to calculate predicted output as well as update the estimator model of A(t,z-1) and ,.... 1 B( t, z-) Then, an on-line numerical analysis program (such as a Gaussian elimination algorithm) will solve the set of linear equations that will satisfy equation (7.4).

PAGE 86

77 7.1 COMPENSATOR DESIGN FOR POLE PLACEMENT CONTROL As a variation from the strict pole placement method, consider using a compensator in the control system. Using a compensator technique is somewhat analogous to the pole placement method. However, there is a distinct difference. Whereas the pole placement method will give a direct analytical location for pole assignment, the compensator design techniques are more of a trial and error art where the compensator tries to "steer" the root locus plot onto a suitable pole location. Staying with the Root Locus techniques, the compensator design method is as follows. If it is desired to independently control the 1st mode, than it is necessary that the effects of the controller on the other dominant modes are eliminated (for this case modes 1 ,2, and 3 are considered dominant, mode 4 can be considered as a residual mode). To do this the compensator zeros are placed very near the physical system poles (e.g., directly beneath them) for modes 2 and 3. Likewise the compensator poles should be placed directly under the 2nd and 3rd modes' zeros. The idea behind this scheme is that the 2nd and 3rd modes will remain stationary as the gains are changed since the poles will approach the zeros in the same relative location. Note, that the system poles and zeros are all on the unit circle since the physical system has no inherent .damping. Placing the compensator poles and zeros a short distance beneath the exact system zeros and pole locations, as opposed to directly on top of them, should eliminate the possibility creating system instability

PAGE 87

7.8 by accidentally placing the compensator zeros a small amount outside the unit circle. Now the gain must be found that will give a suitable transient response. As done previously in the Root Locus section of this paper, the numerical technique of Bairstow can be used to find the minimum (or desirable) value for the "r" in the z-plane. In a sense, the compensator can be thought of as the digital filter used previously. As before the following characteristic equation is solved by the numerical Bairstow algorithm. D(z)*CD(z) + K1*N(z)*CN(z) = 0 (7.1.1) Where D(z) and N(z) are the physical system denominator and numerator transfer function parameters respectively; likewise, CD(z) and CN(z) are the compensator denominator and numerator transfer function parameters. To test this theory the compensator for the 1st mode was created by placing the compensator zeros directly beneath the 2nd and 3rd modes' poles and placing the compensator poles directly beneath the 2nd and 3rd modes' zeros. The following Compensator was found: (7.1.2) c ( z) = 3"'2""8'-::<1 4-;,.zT---1 r+_ ..::..9=-'99=-:9:-=z:-::---4 """"ll"' 1 .1946001z-l + 1.271007z-2.1944628z-3 + .99979z-4 Next, the numerical root locus technique of equation (7.1.1) was used to find a suitable gain. The results of the program show that a gain of K1= 24.6 will produce an optimal response (r = .573) however smaller gain can be used and still produce quite a good response. Now the system can be simulated by adding the compensator

PAGE 88

79 into the control loop. The system is simulated with a gain of both 4.5 and 10. The diagram below illustrates the control system. R + u Physical System v@ 4 y y K1 Compensator J -Figure 26. Independent Modal Control Using a Compensator The results of the simulations were extremely good. By examining the plots of Figures 27 and 28 on the following pages it is easy to see that not only did the effects of 1st mode quickly deteriorate, but the control input was in the form of a sinusoidal force equation with a frequency equal to the first mode. This is quite desirable for a couple of reasons: first it shows that the controller is only controlling the 1st mode as designed; and secondly, the physical actuator will be more able to create such a force input requirement (i.e., there are no large jerks required from the control actuator).

PAGE 89

.3 -r-,3 _.._ 4 -4 TIP UELOCITV vs TIME INDEPENDENT MODAL CONTROL USING A COMPENSATOR Feedback Gain is 4.5 TIME (seconds) CONTROL INPUT vs TIME INDEPENDENT MODAL CONTROL USING A COMPENSATOR Feedback Gain is 4.5 TIME (seconds) Figure .27. Simulation of Compensator Modal Control With Gain Equal to 4.5. 80

PAGE 90

.3 -.3 4 -4 TIP UELOClTY vs tlKE IHDEPEHDEHT HODAL CONTROL USING A COMPENSATOR Feedhack Cain is 19 TIME (seconds) CONTROL INPUT. vs TIME. INDEPENDENT MODAL CONTROL USING A COMPENSATOR Feedback Gain is 19 TIHE (seconds) Figure 28. Simulation of Compensator Modal Control With Gain Equal to 10. 81

PAGE 91

82 Next a compensator was designed to control the 2nd mode by canceling the control effects of the 1st and 3rd modes. The resulting compensator is C(z) = 1 .92024z-1 + .368992z-2-.9202875z-3 + .999932z-4 (7.1.3) 1 -.23808z-l + .2379111z-2-.999831z-3 Note that the order of the denominator is less than the numerator. This is because the 1st mode has only one zero (at z = 1; another zero is assumed to be at negative infinity). The gain for the compensator was found as before, and the system was simulated. The results of the system where not nearly as good for mode 1. This is because as the gain for the second mode compensator was varied, the value of "r" for the second mode changed very little. The root locus plot tended to follow an arc just inside the unit circle. Thus, when simulated, the 2nd mode's controller had little impact on dampening the second mode. To try to fix this situation, a slight modification to the control algorithm scheme was made. Instead of cancelling the zeros for modes 1 and 3, cancel the zeros for modes 2 and 3 instead. Therefore, the root locus plot will act as though the 1st mode's zeros are the-2nd mode's zeros. The new compensator was calculated as: (7.1.4) C(z) = .920259z-1 + .369026z-2 .9202919z-3 + .9999234z-4 .194600z-l + 1.27100z-2.1944626z-3 + .9997986z-4 The characteristic equation was numerically solved as before to find a suitable gain. Any gain from 4 to 34 appears to give a very good settling time for this compensator (with 34 being the best). Again, the system was simulated. The gain-was chosen as 8 to reduce the

PAGE 92

83 large control force needed if a higher gain was used. The results of the simulations are very promising. The 2nd mode was indeed independently controlled with this compensator controller. Next, both modes 1 and 2 where independently controlled at the same time using the compensator independent modal control scheme. As before, the gain for mode 1 (K1) was chosen as 4.5 and K2 was chosen as 8. The figure below illustrates the control system. R + u Physical System v 1st Mode's Compensator 2nd Mode's Compensator 4 Figure 29. Compensator Independent Modal Control For Modes 1 and 2 The results of the simulations found in Figure 30 are excellent. The system control response was almost comparable with the ideal modal control response shown back in Figure 12. By properly selecting a good compensator, independent modal control appears to work extremely well.

PAGE 93

.3 .3 TIP V[L((IfY vs mE IHDIPDIDDit ftODAL COHtROL USIHC A COIIPDISATOR ftodal Fudhck Cains Kl : 4.5 and X2 : 8 till (m onlsl IHP!It vs JIKI lmEPOOOO MODAL US HC A COIIPIMSAIOR lsi Modal Control lnPIII with Cdn Kl : 4,5 tiKI lsmndsl COHtROL IHPUT YS IIJU: IHDJ:POOD!I I!ODAL COI!tROL US HC A COIIPDISATOR 2nd Modd Control lnPIII v th Cain )'2 : 8 Tift lmondsl Figure 30. Simulation for Independent Modal Control for the 1st and 2nd Modes. 84

PAGE 94

85 7.2 COMPENSATOR DESIGN WITH A SELF TUNING ADAPTIVE CONTROLLER Although the compensator design worked well for this known system, could it be used on an system with uncertainties in the parameters? If an adaptive estimator is built into the system as shown in Figure 31, then it is proposed that the control system can self tune itself to account for modeling uncertainties. R I Physical System 1st Mode s Compensator ft 2nd Mode's Compensator I I Digital Compensator Function Generator Figure 31 Adaptive STR Compensator Control Scheme The control system works as follows. First the output from the tip sensor and the known input from the control actuators are fed into the adaptive estimator. The adaptive estimator in turn generates the latest estimate for the system Transfer Function (TF) coefficients using the method of least Next, the transfer function coefficients are fed into the digital compensator function generator (DCFG). The DCFG will use a numerical analysis to calculate the poles and zeros of the

PAGE 95

86 system transfer function. Then by using the pole and zero cancelling process described in section 7.1, the DCFG will assign the zeros and poles to each compensator, calculate the compensator coefficients, and adaptively update the compensators in the control loop. Finally, the coefficients of the system transfer function as well as the digital compensator coefficients are fed into the Optimal Gain Generator (OGG). The OGG will first multiply the system transfer function with the compensator transfer functions. The resulting function can then be transformed into its characteristic equations as shown in equation (7.1 .1). A numerical analysis iteration algorithm (such as the Newton method) combined with the root finding method would then be used to find the proportional gain Kn for each mode n that results in a desired value for "r". Where "r" is the distance from the origin to the nth root (pole) location on the z-plane. The desired value for "r" may be its minimum value, or some other value that will produce a suitable settling time. The new gain value is in turn fed to the control actuators thus completing the control loop. It should be found that once the physical system is adequately known, small variations in the parameters (due to the inherent changing nature) will not require large changes in the gain. Therefore, some computer processing time could be saved by only updating the control gains at certain intervals, as opposed to at after each sample cycle. In the following chapter this same self tuning principle will be used in another variation of a similar control scheme. In addition, the effects of stochastic noise on the system and how it

PAGE 96

87 can be reduced will be reviewed. The findings from the next section for reducing stochastic noise is directly applicable to the adaptive compensator control scheme just described.

PAGE 97

CHAPTER VIII AN ADAPTIVE CONTROL SCHEME FOR INDEPENDENT MODAL CONTROL Another possible adaptive control scheme that takes into account all the findings and theories described in the first part of this thesis is shown in Figure 32. System / 11 I 1 TF coefficients Figure 32. Adaptive STR Control Scheme The control system works as follows. First the output from the tip sensor and the known input from the control actuators are fed into the adaptive estimator. The adaptive estimator in turn

PAGE 98

89 generates the latest estimate for the system Transfer Function (TF) using the method of least squares. Next, the transfer function coefficients are fed into the digital filter generator (DFG). The DFG will use a numerical analysis technique (such as Bairstow's) to calculate the poles and zeros of the system's transfer function. The digital filter for each control mode is now created by centering the filter response to the system poles (modal frequency). A fairly simple computer algorithm process will generate the digital filter coefficients of the transfer function. The filter could be a Butterworth, Chebyshev, or any other recursive filter. Thus an adaptive digital filter is created for each controlled mode. Finally, the coefficients of the system transfer function as well as the digital filters coefficients are fed into the optimal gain generator (OGG). The OGG will first multiply the system transfer function with the digital filters transfer functions. The resulting function can then be transformed into its characteristic equation. A numerical analysis iteration algorithm (such as the newton method) combined with the root finding method would then be used to find the proportional gain Kn for each mode n that results in the minimum value for "r". Where "r" is the distance from the origin to the nth root (pole) location on the z-plane. The new gain value is in turn fed to the control actuators to complete the control loop. The control system just described will adaptively optimize the state variable feedback system when proportional feedback is used. That is, for a proportional feedback system, this technique

PAGE 99

90 will produce the optimal transient response in terms of damping or settling time. Listed below are some additional advantages in using this control scheme: 1) It is adaptive so the control system will change with time as the parameters of the physical system change (due to thermal expansion/contraction, aging, mass changes due to additions on the LSS, etc.). 2) It uses the simplicity of proportional feedback law--no complex algorithm is needed to generate the control law. 3) The nature of proportional feedback makes the control system rather robust. On the other hand, this control scheme does have some drawbacks. Mainly, the system is considered deterministic (no stochastic disturbances are considered in the design). Also, the transient control response is limited by the nature of proportional feedback methods. In other words, generally it is not possible too select the desired pole locations that give the desired transient response and then find the proportional gain needed. This control scheme is limited to finding the best possible solution with the inherent limitations of proportional feedback. Finally, a relatively large amount of numerical data processing is required. This could have the undesirable effect of overloading the capabilities of the onboard computer. It is possible to reduce the data processing burden by only updating the adaptive filter and control gains at certain intervals rather than at each sample period.

PAGE 100

8.1 ADAPTIVE CONTROL WITH STOCHASTIC NOISE The next logical step is to consider the effect stochastic noises. The state-space model for the stochastic system can be written as follows: 91 x(k+1) Y(k) A*x(k) + B*u(k) + v1(k) C(k) + v2(k) (8.1.1) Where A is considered to be relatively constant with time and v1(k) and v2(k) are the disturbance noise and the measurement noise respectively. The disturbance noise for a Large Space Structure can be caused by a number of factors: unmodeled (residual) modes, nonlinear structure behavior, onboard disturbances (motors, pumps, etc), and environmental conditions (thermal gradients, solar pressure, gravity gradients, atmospheric effects, etc) [1]. The measurement noise can be caused by the limitations in the sensor design itself. By looking at equation (8.1.1) it is easy to see that even if the parameters of the physical system were defined exactly, the estimator would not realize this and try to adjust the parameters such that it matched the output results. That is, the parameter estimator would assume that the noise was caused by a deterministic inputs and would begin to chase these "ghosts" that just are not there. Therefore, it is desireable to create a filter that could filter out the noise and leave only the actual deterministic output such that the parameter estimator can give the optimal estimate for the system parameters. To do this, the Kalman filter will be considered.

PAGE 101

92 8.2 KALMAN FILTERING Kalman filtering theory is a topic that has been quite thoroughly reviewed over the last 25 to 30 years. Therefore, the basic theory fundamentals and the associated discussion are left out of the paper. Instead, just a brief review of the algorithm for implementing the Kalman filter on line is given. Assuming the control designer has a reasonably good idea of the noise covariance matrices Q, R, and S where Q = E{v1(t)v1(s)} (8.1 .2) R = E{v2(t)v2(s)} (8.1.3) and s E{v1(t)v2(s)}, (8.1 .4) if it is further assumed that v1(t) and v2(t) are uncorrelated (i.e., the covariance matrix S = 0), then implementing the discrete Kalman Filter is fairly simple. First, an initial assumption for the error covariance matrix Po must be made where Po E{x(o) xo )( x( 0) xo) T) (8.1.5) Usually, any reasonable positive definite choice for Po should be adequate. Often, Po is chosen as the identity matrix [I]. The Kalman filter gain itself K is usually taken as zero initially. The following algorithm shows the steps that can be taken to implement the Kalman filter on line.

PAGE 102

93 (8.2.1) x(-k) = A(k-1)*x(k-1+) + B*u(k) (State estimate) 2. P(-k) = A(k-1)*P(k-1+)A(k-1)T + Q(k-1) (Error covariance) x(k+) P(k+) x(-k) + K(k)[Y(k) -C(k)x(-k)] [I -K(k)C(k)]P(-k) (State estimate) (Error covariance) 5. K(k+) = P(-k)C(k)T[C(k)P(-k)C(k)T + R(k)]-1 (Kalman Gain Matrix) where the sign (-) means just prior to and (+) means just after the discrete time increment. It is quite simple to implement these five steps on a computer. Notice how the state variable matrices A, B, and C from equation (8.2.1) are now considered time variant. Recall that the adaptive least squares estimator and predictor will produce the latest estimates of the DARMA model parameters. To implement the output of the parameter estimator into the state variable matrices A, B, and C directly, the matrices could be transformed into their observable form. In the observable form, the parameters of the matrices are as follows. lg 0 0 -an-1 bn-1 (8.2.2). 1 0 A = ... -a4 B = b3 c = [ 0 0 0 1 0 -a3 b2 0 1 0 -a1 b1 0 0 1 -ao bo The parameters for {a0 an-1, b0 bn-11 are found directly from the least squares adaptive estimator output 8 (k). When the Kalman filter is put into the observable state-space model form, the states x(k) do not have the same physical meaning as the physical J state-space form (i.e., they no .longer are the velocities at the node

PAGE 103

94 locations). Now they are the past values for V However, the output and input Y(k) and U(k) are still the same as before. Notice that the observable form of (8.2.2) is the exact equivalent of the DARMA model. If the Kalman filter is entered into the proposed control scheme, the control system will look as shown in Figure 33. / System I I I I I TF coefficients y Figure 33. Adaptive STR Control System with Kalman Filtering

PAGE 104

95 The control system is similar to the previous scheme. However, now the output from the sensor is fed into the Kalman filter instead of into the adaptive parameter estimator. Also, the known input of the system u(k) is fed into the the Kalman filter in addition to being fed into the adaptive estimator. The output of the Kalman filter is the estimated (noise filtered) output. This output is fed into the adaptive parameter estimator and is used to generate the parameters for which are the transfer function coefficients. As before these updated coefficients are fed into the DFG and OGG; in addition they are fed back into the Kalman filter to update the filter's state matrices. In other words, the Kalman filter estimates the states and the adaptive estimator estimates the parameters. One of the great features of the Kalman filter is that if the noise is Gaussian, the filter gives the minimum variance estimate of the state. If the gaussian assumption is removed, the filter will still give the linear minimum variance estimate of the state. Therefore, it is a quite attractive scheme and is well used in practice. The adaptive control scheme with the Kalman filter just described could be directly implemented on the compensator independent modal control scheme from the previous chapter with very little alterations. To test how well the Kalman filter combined with a least squares parameter estimator can find the actual parameters of the system in a noisy environment, another-pair of computer simulations were performed. An intentional error is made for the initial

PAGE 105

96 parameters of the estimator model. Sensor noise is added to the actual physical system model. To make the noise approach a Gaussian distribution, the computer's random function generator is convoluted three times. Figure 34 illustrates the noise probability density functions. E{x} -.1 1 -.1 1 -. 1 Figure 34. Noise Probability Density Functions For the computer simulations, the noise probability density function range in value from -0.1 to +0.1. This represents a 3-sigma (standard deviation) noise dispersion of 33% of the actual measurement values (i.e., the measurement readings expected values range from .-0.3 to +0.3). The simulated system is shown in Figure 35. v2 u 1_-+ ___ Y __ ----t""..,l Kalman I P' 1 y A e ... I Least Squares Estimator 1 .... Figure 35. Self Tuning Parameter Estimator With Stochastic Noise

PAGE 106

97 The parameters found from the Least Squares Estimator are now compared to the actual system parameters. Two simulations were made: First, the resulting error of the least squares estimator without the Kalman filter in a noisy environment was found. The error was defined as the difference between the actual output found by using the actual deterministic system parameters, and the predicted values found by using the parameters of the least squares estimator. The second simulation found the resulting error when the Kalman filter is added as was shown in Figure 35. Figure 36 shows the results of the two computer simulations.

PAGE 107

... .. e ... .. 0 ... .. e .. .... 98 ABSOLUTE VALUE FOR THE ERROR With NOISE 0.026 . 0.022 o.o: O,QUi 0,0,4 o.o1z 0.01 o.ooe o.OOG 0.004 o.ooz -.. 11-I i .... .. 0 ... ... r ''I' 0 z -' 4 G 7 ,."'" (Seconds) ABSOLUTE VALUE FOR THE ERROR With NOISE ba:rt with t
PAGE 108

99 By analyzing Figure 36, it is easy to see that the Kalman filter does improve parameter estimator for this noisy system. The improvements were not as spectacular as one might hope. However, if the covariance matrices are chosen more carefully, it can be expect that even further improvements in the parameter estimation are possible. It is also interesting to note that despite the addition of noise, the least squares estimator alone (without the Kalman filter) still converges the parameter estimates toward their actual values; albeit, at a very slow rate. In this discussion it is assumed that the properties of the stochastic noise is fairly well known (covariance matrices Q, R, and S). This assumption is probably not reasonable for a LSS with many unknowns and untestable properties. In Goodwin and Sin's book "Adaptive Filtering Prediction and Control" an algorithm called "Improved Extended Kalman Filter for Linear Systems with Unknown Parameters" is given [s]. The algorithm takes into consideration the incorrect specifications of noise covariances and incorrect parameter estimates. This algorithm may better suit the stochastic control for a LSS with all its unknowns. Unfortunately, the algorithm contains many partial derivatives and is quite difficult to implement. However, it may be of interest for future studies.

PAGE 109

CHAPTER IX CONCLUSION This thesis presents some of the considerations that are needed for LSS modal control design. It was shown that a linear system can be modeled with a series of 2nd order differential equations. Later it was demonstrated that the model parameters can be further refined with the adaptive estimation technique of Least Squares. Several different control schemes were presented. If the performance index is based on the control input force required, then independent modal control is much better than other control designs such as direct. velocity feedback or step ahead control. Two different modal control methods were presented: filtering.the modes independently and then applying proportional feedback, or using a compensator to independently control the modes. Both schemes work fairly well. However, the compensator control method could produce a better transient response; the independently filtered control method was inherently more robust. The non-modal control method of direct velocity feedback worked well in the simulations, howev-er the control input required was unsatisfactory due to the large jerks needed. The d-step-ahead control method was overall unsatisfactory due to the control force requirements and the control sensitivity of a lightly damped system.

PAGE 110

1 01 It was shown that either one of the independent modal control I techniques could be combined with an adaptive estimator to form a iself tuning regulator adaptive control system. In a noisy environment, which is the inherent environment for a LSS, adding a Kalman filter into the control loop to estimate the states can further improve the control. In this case the adaptive estimator will estimate the system parameters and self tune the control system; the Kalman filter will estimate the system states.

PAGE 111

102 REFERENCES 1. BALAS, M.J., "Trends in Large Space Structure Control Theory: Fondest Hopes, Wildest Dreams." IEEE Transactions on Automatic Control, VOL AC-27, NO. 3. June 1982. 2. Shigley, J.E., Mechanical Engineering Design, McGraw-Hill, New York, 1977 3. Atluri, S.N., and A.K. Amos, Large Space Structures: Dynamics and Control. Springer-Verlag, Berlin, 1988. 4. Gehling, R.N., "Effect of Passive Damping on Active Control Energy Requirements." PACOSS, F33615-82-C-3222, AFWAL, March 1986. 5. BALAS, M.J., "Feedback Control of Flexible Systems." 6. IEEE Transactions on Automatic Control, VOL. AC-23, NO. 4. August 1978. Burke, S.E. and J.E. Hubbard Jr (1988). Control Design for Flexible Beams." pp 619-627. "Distributed Actuator Automatica, Sep 1988, 1. Chalam, V.V., Adaptive Control Systems Techniques and Applications. Marcel Dekker, Inc., New York, 1987. 8. Goodwin, G.C., and K.S. Sin, Adaptive Filtering Prediction and Control, Prentice-Hall, New Jersey, 1984.

PAGE 112

APPENDIX A COMPUTER PROGRAMS USED FOR DESIGN The program code and run listings contained in Appendix A were used to design and develop the system model as well as the different control schemes for this thesis.

PAGE 113

100 REft t 110 REn I fiNDS DISCRETE VERSION 120 REft CONTINUOUS SYSTEft, I,, 130 REft I X{t+1}:fDX{t}+GDtU{t}, REft GIVEN THE CONTINUOUS VERSION 150. REft I XDOT=fX+GU 1b0 REn AND THE SAftPLE TinE T. 170 REft I 180 REft I THEN fiNDS TRANSfER fUNCTION REn 1 Y{Z}/U{Zl. 200 REft I WHERE Y:HX, GIVEN H. 210 REft lltllllltllllllllllllllltttlll 220 REft 230 INPUT "Enter the order of the system :",N Dift f{N,N},fD{N,N},Ul{N,N},E{N,N},Y{N,N},Z{N,N},S{N,N},Sl{N,N} 250 INPUT "Enter the sampling time Ts :",T 2b0 INPUT "Enter the number of expansion terms desired :,n 270 If {ftl3 AND ftC25} THEN 300 280 PRINT"Choose another value. Normally a nlllter between 8 and 15 gives a sat is factory'" PRINT"result.": GOTO 2b0 300 PRINT: PRINT"ENTER THE VALUES fOR CONTINUOUS TiftE DYNAftiCS ftATRIX (f)" 310 PRINT: fOR I : l TO N: fOR J : l TO N 320 PRINT"f{";I;",";J;"}: ";: INPUT f{I,d} 330 NEXT J,l : PRINT PRINT"ENTER THE VALUES fOR CONTINUOUS TlftE CONTROL DISTRIBUTION VECTOR lGl 350 PRINT: fOR I : 1 TO N: PRINT '"'i;.J;i: ";: INPUT UU: NEXT : PRINT 3b0 PRINT"ENTER THE VALUES fOR CONTINUOUS TlftE OUTPUT DISTRIBUTION VECTOR lHJ'" 370 PRINT: fOR I : 1 TO N: PRINT "Hi;I;"} : ";: INPUT H{I}! NEXT 380 DEf fN R{S} : INT {100000! 1 S + .5} I 100000! REft REn 1 fiNDS fD AND &D 1 REft PRINT: PRINT TAB{ lS}"NUftBER Of TERftS TO GO:" PRINT TAB{ 10}; fOR I : 1 TO N: fOR J : l TO N:Ul{I,J} : O: NEXT J,I: fOR I : 1 TO N:Ul{I,I} : l: NEXT fOR I : 1 TO N: fOR J : l TO N:fD{I,J} : Ul{I,J}: NEXT J,I fOR t : ft TO 1 STEP 1: If t I l THEN fOR I : 1 TO N: fOR J : l TO N:E{I,J}: T 1 fD{I,J}: NEXT J,I fOR I : 1 TO N: fOR J : l TO N:Y{I,J} : {T I t} I f{I,J}: NEXT J,I 500 fOR I : 1 TO N: fOR J : l TO N:Z{I,J} : 0: NEXT J,I: fOR I : 1 TO N: 510 fOR J : l TO N: fOR l : l TO N 520 Z{I,J} : Z{I,J} t Y{I,L} I fD{L,J}: NEXT L,J,I 530 fOR I : 1 TO N: fOR J : l TO N:fD{I,J} : Z{I,J} + Ul{I,J}: NEXT J,I: PRINT t:: NEXT t 550 PRINT : FOR I : l TO N: FOR t : 1 TO N:GD{I} : GD{I} t HI,O 1 G{[}: 5b0 NEXT C,l 570 PRINT: PRINT TAB{ lO}"TRANSITION ftATRIX fD IS:": PRINT 560 fOR I: l TON: PRINT TAB{ 10};: fOR J : l TO N:.PRINT fN R{fD{I,J}}:: 590 NEXT : PRINT : NEXT bOO PRINT: PRINT TAB{ lO}"CONTROL DISTIB. VECTER 'D IS:": PRINT blO fOR I : l TO N: PRINT TAB{ 10}: fN UU}:, NEXT : PRINT 104

PAGE 114

b20 REn 105 b30 REM 1 FINDS COEff. DENOM.& NUn. POLYNOMIALS USING LEVERRIER LAGORITHM REM b50 PRINT: Of DENOftiNATOR & NUMERATOR POLYNOMIALS IN THE fOL LOVING fORft:":PRINT bbO PRINT TAB{12} "ncoeff{1}1Z"{n-1} t nqoeff{2}1Z.{n-2} 1 1 ncoeff{n}" b70 PRINT"T{Z} : ------------------------------------------------------------b71 PRINT TAB{9}"Zn t dcoeff{l}IZ"{n-1} t t t dcoeff{n}" b85 PRINT: PRINT TAB{?} "coeft # DENOft. COEff NUn. COEff-b90 FOR I : l TO N:S{I,I} : 1: NEXT : fOR L : 1 TO N: FOR I : l TO N: 700 G1{1} : 0: NEXT : fOR I : l TO N 710 fOR J : l TO N:G1{I} : G1{1} t S{I,J} 1 GD{J}: NEXT J,I:N{L} : 0: 720 fOR I : 1 TO N 730 N{L} : N{l} I H{l} I NEXT : fOR I : l TO N: fOR J : l To N: S1{l,J} : O: NEXT J,I 750 fOR I : l TO N: fOR J : 1 TO N: fOR ( : l TO N: 7b0 S1{I,J} : S1{I,J} t S{I,C} I fD{C,J}: NEXT (,J,I 770 D : O: fOR I : 1 TO H:D : D Sl{I,I} I L: NEXT :D{l} : D: 780 PRINT TAB{ 1D}:L," "D," "N{L} 790 fOR I : l TO N:Sl{I,I} : Sl{I,I} t D: NEXT : fOR I : l TO N: 800 fOR J : 1 TO N:S{I,J} : Sl{I,J} 810 NEXT J,I,L 820 REM 830 REft 1 fiNDS POLES Of TRANSfER fUNCTION 1 REn 850 PRINT: PRINT TAB{ 22}"POLES ARE:": PRINT TAB{ 12}" REAL.," IftAG.": GOSUB 1000 l!bO REn 870 REn 1 fiNDS ZEROS Of TRANSfER fUNCTION & THE GAIN I 860 REn 890 L1 : 1: fOR I : 1 TO N: If ABS {N{I}} > .00001 THEN 910 900 L1 : 1 t I: NEXT 910 fOR I : l TO {N L1}:D{I} : N{l t Ll} I N{L1}: NEXT :N : N L1: 920. If N : 0 THEN 930 PRINT: PRINT TAB{ 22}"ZEROS ARE:": PRINT TAB{ 12}" REAL"," InAG.": '0 SUB 1000: GOTO 9b0 950 PRINT TAB{ 10}"THERE ARE NO ZEROS." 9b0 PRINT: PRINT TAB{ 10}"GAIN IS:" fN R{N{ll}}: END 970 REM 960 REn a SUBROUTINE: fiNDS ROOTS Of POLYNOMIAL USING BAIRSTOY ALGORITHM 990 REM 1000 J : l:n : O:E : .OOODl:RO : 2:50 : l lOlO l : D:Rc : RO:Sc : SO:ftl : N 2 1 n: If Ml : 2 THEN 1020 If ftl > 2 THEN 1050 1030 R{J} : D{1}:l{J} : O: GOTO 1220 R2 : D{l}:Sc : D{2}: GOSUB 12b0: GOTO 1220 1050 B{O} : 1:8{1} : D{l} Rc: FOR [ : 2 TO nl: lObO B{(} : D{[} R2 I 8{( l} S2 I B{( 2}: NEXT

PAGE 115

1070 P{O} : O:P{1} : 1: fOR ( : 2 TO ft1: 10&0 P{(} : 8{( 1} R2 I P{( -1} S2 I P{[ 2}: NEXT 1090 Q{O} : O:Q{1} : 0: fOR ( : 2 TO ft1: 1100 Q{(} : B{[ -2} -R2 I Q{( -1} -S2 I Q{[ -2}: NEXT 1110 R3 : B{ft1 1}:S3 : B{ft1} f R2 I B{ftl : P{ft1 -l}: 1120 = Pn1l R2 1 Pn1 1} Bnl 1} 1130 R5 : Q{ft1 1}:S5 : Q{ftl} t R2 I Q{ft1 l}: Db = 1 ss -RS 1 sq:Rb = s3 1 Rs -ss R3J 1 Db 11so Sb = {R3 1 1 s3J 1 Db:R2 = R2 Rb:s2 = s2 sb: 11b0 If ABS {Rb} E OR ABS {Sb} E THEN 1160 1170 GOSUB l2b0: GOTO 1200 1160 If L ) N 1 50 THEN PRINT "POLYNOftiAL SOLVER DID NOT CONVERGE": END 1190 L : l t 1: GOTO 1050 1200 ft : ft t 1:J : J t : N 2 1 ft: FOR [ : 1 TO N3:D{[} : 8{[}: NEXT : 1210 GOTO 1010 1220 fOR [ : l TO N: PRINT TAB{ 10}: FN l{R{[}}," ftfN R{I{[}}: NEXT : RETU RH 1230 REn REft 1 SUBROUTINE: SOLVES QUADRATIC EQUATION 1250 REft lebo P9 = R2 1 R2 -q 1 s2: If P9 o THEN 129o i270 R1 : SQR { P9}:RJ} : R2 I 2:RJ t 1} : R2 I 2:I{J} : Rl I 2: 1280 I{J t 1} = -R1 I 2: RETURN 1290 R1 : SQR {P9}:1{J} : O:I{J t l} : O:R{J} = {R1 R2} I 2: 1300 R{J t l} : { -Rl R2} I RETURN

PAGE 116

RUN Enter the order of the system :8 Enter the ing time Ts :.05 Enter the number of expansion terms desired :15 ENTER THE VALUES fOR CONTINUOUS TinE DYNAftiCS ftATRIX lfJ f{ 1 1 } : 0 f{ l } : 0 f{ 1 3 } : 0 f{ 1 1 } : 0 f{ 1 1 5 } : -1500 f{ 1 I L } : 750 f{ l 1 1 } : 0 rt 1 I a l = o f{ 2 l } : 0 f{ 2 1 2 } : 0 f{ 2 1 3 } : 0 f{ l 1 } : 0 f{ 2 1 s } : :750 f{ 2 1 b } : -1500 f{ i! 1 1 } : 750 f{ .j! 1 6 } : f 0 f{ 3 1 l } : 0 f{ 3 1 2 } : 0 f{ 3 1 3 } : 0 f{ 3 1 } : 0 f{ 3 1 5 } : 0 f{ 3 b } : 750 f{ 3 1 1 } : -1500 f{ 3 8 } : 750 f{ 1 1 } : 0 1 i!}: 0 f{ 1 3 } : 0 f{ 1 } : 0 f{ 1 5 } : 0 f{ b } : 0 f{ 1 } : 500 f{ 6 } : -500 f{ 5 1 } : 1 f{ 5 2 } : 0 f{ 5 3 } : 0 f{ 5 } : 0 f{ 5 5 } : 0 f{ 5 1 b } : 0 f{ 5 1 } : 0 F< s a l = o 107

PAGE 117

f{ b l } : 0 f{ b } : l r< L 3 l = r o f{ b } : 0 f{ b s } : 0 f{ b b } : 0 f{ b 1 } : 0 f{ L 6 } : t 0 f{ 1 ; l } : 0 f{ 7 } : 0 f{ 1 3 } : f l f{ 7 } : 0 f{ 7 5 } : 0 f{ 1 L } : D f{ 7 1 } : 0 f{ 1 6 } : 0 f{ 6 l } : t 0 f{ 6 } : t 0 f{ 6 3 } : 0 f{ 6 } : 1 f{ 6 s } : t 0 f{ 8 L } : f 0 F{ a 1 } : o f{ 6 6 } : f 0 ENTER THE VALUES fOR CONTINU9US TinE CONTROL USTRIBUTION VECTOR (G) H l } : f 0 G{3}:f0 G{ } : .bbbbbbb7 H S } : f 0 G{b}:fO G{ 7 } : t 0 G{ a J = t o ENTER THE VALUES FOR CONTINUOUS TIME OUTPUT DISTRIBUTION VECTOR lHl H{ l } : t 0 H{ 2 } : t 0 H{ 3 } : 0 H{ } : t 0 H{ 5 } : 0 H{ b } : t 0 H{ 7 } : f 0 H{ a } = t 1 108

PAGE 118

TRANSITION ftATRIX fD IS: -.25866 .46602 .09909 .0073 -29.68679 6.30802 .60676 -.l59b .10882 3.65011 -23.57877 7.36637 .09909 .18614 .60668 b.30802 q,38929 -2q,9b2Sq l3,q2902 .ooqa7 .07255 .40579 .53918 6.952679 .02523 .D106 .D0111 .DODDS .258b8 ,q6602 .099D9 .0073 .0106 .02b3q .OlD6q .OD119 . 46602 .1596 ,q7089 .10882 .00111 .02b02 .01217 .099D9 .186lq .60868 .00079 .00811 .07255 .51bb2 CONTROL DISTIB. VECTER GD IS: 3.6b0134E-05 8.ll32b9E-03 2.737l34E-02 2.370317-07 6.962379-06 1.104155E-Oq 7.Sq92EOq 109

PAGE 119

110 COEffiCIENTS Of l NUftERlTOR POLYHOftiALS IN THE fOLLOWING fORft: -ncoeff{UZ{n-U t ncoeff{c}r{n-2} t ... t ncoeff{n} T{l} : ----------------,-----zn dcoeff{1lZ"{n} dcoeff{2}Z{nc} t dcoeff{n} coeff # C6Eff HUM. COEff J. 2 l.758U5E-03 3 .36&3131 2.31.,31 c.SO,lLLE-03 5 .36&3132 2.50,1bbED3 b 2.3b,30,ED3 1 .11556,3 a ,,,,,,L POLES ARE: REAL InAG. ,,0676 ,,0676 .26qs, .26qs, -. ,56L5 6,371 -,6,371 -.63252 .ssq .63252 -.ssq ZEROS ARE: REAL !MAG. .50057 .aLSb, .50057 -.AbSb, .saoq2 -.a1q32 .saoq2 -.3507 ,]bq, -.3507 -. 0 GUN IS: ,00075

PAGE 120

250 INPUT "ENTER THE ORDER Of THE SYSTEn : ", N 260 PRINT: PRINT "ENTER THE TERnS Of THE [f) nATRIX BY ROVtt": PRINT 270 Din f{N,N},fl{N,N},A{N,N},T{N,N},AI{N,N},U{N,N},D{N},R{N},I{N},Rl{N} 280 Din Il{N},E{N},B{N},P{N}: fOR I : l TO N: fOR J : 1 TO N PRINT "f{":I:",";J;"} : ": : INPUT f{I,J} : NEXT J,I 300 PRINT: INPUT "ENTER THE NUnBER Of CONTROLS : ", Nl : PRINT HO PRINT "ENTER THE TERns Of THE CONTROL DISTRIBUTION nATRIX IGJ BY COLUnN" 320 PRINT: Din G{N,Nl},C{Nl,N}: fOR J : 1 TO Nl: fOR I : 1 TO N 330 PRINT "G{":I:",":J:"} : ": : INPUT G{l,J} : NEXT I.J: PRINT PRINT "EtHER THE TERnS Of fHE GAIN ftATRIX [(] uBY ROWu": PRINT 350 fOR I : l TO Nl: fOR J : l TO N: PRINT "C{":I:",":J:"} : ": 360 INPUT C{I,J} : NEXT J,I : PRINT 370 INPUT "ENTER THE NUMBER Of MEASUREMENTS : ",N2 : PRINT 380 PRINT "ENTER THE TERMS Of THE OUTPUT DISTRIBUTION nATRIX IHI uBY ROWn" PRINT: Din H{N2,N},l{N,N2}: fOR I : 1 TO N2: fOR J : l TO N qoo PRINT "H{":I:",":J:"} : ": : INPUT H{I,J} : NEXT J,I: PRINT PRINT"ENTER THE TERMS Of THE ESTIMATOR GAIN MATRIX Ill BY COLUMN": PRINT fOR J : l TO N2: fOR I : 1 TO N: PRINT "l{":I:",":J:"} : ": INPUT l{I,J} : NEXT l,J : PRINT DEf fN R{S} : INT {100000! r S t ,5} I 100000! REM REM I CALCULATES f-GC-LH I no REn -fOR I : l TO N: FOR J : 1 TO N: fOR { : l TO Nl: f{l,J} : f{I,J} G{l,(} I {{[,J}: NEXT (,J,l 500 fOR I : l TO N: fOR J = l TO N: fOR t : l TO N2: 510 f{I,J} : f{!,J} l{I,(} I H{(,J}: NEXi (,J,I 520 REM 530 REft REn t fiNDS COEffiCIENTS Of CHARACTERISTIC EQUATION Of f-GC-LH 550 REM 560 PRINT" COEffiCIENTS Of CHARACTERISTIC EQUATION IN THE fOllOWING fORM" 570 PRINT: PRINT" sn t COEff{l}S"{n-1} t COEff{2}S"{n-2} 1 1 580 PRINT: PRINT TAB{10} "COEff NO.","COEff VALUE" fOR I : l TO N:S{l,I} : 1: NEXT : l = 1 TO N: fOR I : l TO N: bOO fOR J : l TO N:A{I,J} : 0: NEXT J,l blO fOR I : l TO N: fOR J : l TO N: fOR ( : 1 TO N: 620 A{l,J} : A{l,J} t 1 f{(,J}: NEXT t,J,l b3D D : 0: FOR I : l TO N:D : D -A{l,I} I L: NEXT :D{l}: D: PRINT TAB{lO}" "D: FOR I : l TO N 1;50 AU,U : A{l,U t D: NEXT : fOR I = 1 TO N: fOR J = 1 TO N: tbO S{I,J} : A{I,J}: b70 NEXT J,I,l 660 REn REM fiNDS ROOTS Of CHARACTERISTIC EQUATION 700 REn 710 J : 1:n : O:E : .00001:RO : .S:SO : 0 720 R2 : RO:S2 : SO:l : O:ftl : N -2 n: If Ml : 2 THEN 750 730 If nl 1 2 THEN 760 111

PAGE 121

7q0 R{J} : : O: GOTO 750 R2 : D{l}:S2 : D{2}: GOSUB ,bO: GOTO 7b0 8{0} : 1:8{1} : D{l} R2: fOR ( : 2 TO ft1: 770 B{t} : D{t} R2 I B{K l} S2 I B{t 2}: NEXT 760 P{O} : O!:P{l} = 1: FORt = 2 TO nl: 790 P{t} : B{t 1} R2 I P{t 1} S2 I P{t 2}: NEXT 600 Q{O} : O:Q{l} : O: FOR K = 2 TO ft1: 610 Q{t} : 8{( -2} R2 I Q{K 1} S2 I Q{t 2}: NEXT 620 R3 : B{ftl l}:S3 : 8{ft1} t R2 B{ftl l}:Rq : P{ftl 1}: 630 sq = Pn1l R2 1 Pnl -1} Bn1 -lJ aqo RS = Qn1 u:ss = Qnn R2 1 Q{ft1 n:DL = Rq 1 S5 RS sq: 650 Rb : {Sl I RS ss I Rl} I Db 6b0 Sb : {R3 I sq Rq 1 Sl} I Db:R2 : R2 t RL:S2 : S2 t Sb: 670 If ABS _{Rb} J E OR ABS{Eb}T E THEN 680 ,bO: GOTO 920 6,0 If l > N 1 50 THEN PRINT TAB{ lO}ftPOLYNOftiAl SOLVER DIDNOT CONVERGrw: If l> N1SO THEN END l : l t 1: GOTO 7b0 ,20 ft : ft t 1:J : J t 2: fOR [ : 1 TO {N 2 I ft}:D{[} : B{K}: NEXT : GOTO 720 930 REft 9q0 REft I SUBROUTINE: FINDS ROOTS Of QUADRATIC EQUATION I REft 9b0 P, : R2 1 R2 q 1 If P, > : 0 THEN ,,0 ,70 Rl : SQR { P9}:R{J} : R2 I 2:R{J t 1} : R2 I 2: 960 I{J} : Rl I 2:I{J t l} : R1 I 2: RETURN 990 R1 : SQR {P9}:I{J} : O:I{J t 1} : O:R{J} : {R1 -R2} I 2: 1000 R{J t l} : { -R1 R2} I 2: RETURN 1010 REft 1020 REft 1 ORDERS EIGENVALUES BY DECREASING ftAGNITUDE I 1030 REft 10qo fOR J : l TO {N 1}:ftft : {R{J}} 2: fOR I : {J t 1} TO N: 1050 ftA{I} : {R{I}} 2 t {I {I}} 2 lObO If ftA{l} c : ftft THEN 1090 1070 ftft : ftA{I}:RR : R{I}:II : I{I}:L : I:R{l} : R{J}:l{L} : I{J}: 1040 R{J} : RR:I{J} : II NEXT I,J UOO I : l 1110 If I{l} : 0 THEN 11q0 1120 If I{I} > 0 THEN 11b0 1130 GOTO 1160 11q0 If I : N THEN 1200 1150 I : I t l: GOTO 1110 llbO If {I t l}_ : N THEN 1200 1170 I : I + 2: GOTO 1110 1160 !{I} : -I{I}:I{I + 1} : -I{I}: If {I t l} : H THEN 1200 1190 I : I + 2: GOTO 1110 1200 PRINT: PRINT TA8{13}ftEIGENYALUES ARE:ft:PRINT TAB{10}ft REAL InAGINARY 1210 fOR t : 1 TO N: PRINT TAB{ 10} fN R{R{[}}:ft ft fN R{I{K}}: NEXT 112

PAGE 122

REn 1 conPUTES EIGENVECTOR nATRIX. T 1 113 REM 1250 fOR Il : 1 TO N: If I{Il} l 0 THEN 1300 l2b0 IF I{I1} : 0 THEN 1290 1270 A : R{I1}:B : I{I1}: GOSUB 13b0 12BO 1300 A : R{I1}: GOSUB 1300 NEXT Il 1310 PRINT: PRINT TAB{ 10}"EIGENVECTOR ftATRIX T, IS:": PRINT 1320 fOR I : 1 TO N: PRINT TAB{ 10}:: FOR J : 1 TON: 1330 PRINT FN R{T{l,J}},: 1340 NEXT : PRINT : NEXT : GOTO 1700 1350 REft 13b0 REn 1 fiNDS REALIIftAG. PARTS Of COftPlEX EIGENVECTOR t 1HO REft 1360 fOR I : l TO N: fOR J : 1 TO N:fl{I,J} : -f{I,J}: NEXT J,I: fOR I : 1 TO N:f1{I,I} : f1{I,I} t A NEXT : fOR I : 1 TO N: fOR J : 1 TO N:A{I,J} : 0: NEXT J,I fOR I : 1 TO N: fOR J : 1 TO N: fOR : 1 TO N: 1420 A{I,J} : A{I,J} t fl{l,[} I fl{(,J}: NEXT 1430 fOR I : 1 TO N:A{I,l} : A{I,I} t B 1 B: NEXT :n2 : 2: GOSUB 2340 1440 fOR I : 1 TO N:R1{I} : X{I}:I1{I} : 0: NEXT : fOR I : 1 TO N: 1450 fOR J : 1 To N 14b0 I1{I} : I1{I} t f1{I,J} I R1{J}: NEXT J,I: fOR I : 1 TO N: 1470 IUI} : Il{I} I B: NEXT. 1460 REn 1 NORnALIZES conPLEX conPONENT={1,0l 1 REft 1500 fOR I = 1 TO N:ft{I} : {R1{I}} 2 t {I1{I}} NEXT 1510 nn : ft{U:l : 1: fOR I : 2 TO N: If ft{l} I nn THEN L : I 1520 If ft{I} I nn THEN nn : ft{I} 1530 NEXT : fOR I = 1 To N:RR{Il = {R1{1l 1 R1{l} 11{1} Il{lll 1 nn 1540 II{I} = {I1{I} Rl{l} Rl{I} Il{l}} 1 nn: NEXT 1550 fOR I : 1 TO N:T{I,Il} : RR{I}:T{I,I1 t 1} : II{I}: NEXT : RETURN 15b0 REn 1570 REn 1 fiNDS COftPONENTS Of REAL EIGENVECTOR 15M REn fOR I : 1 TO N: fOR J : 1 TO N:A{l,J} : -f{I,J}: NEXT J,I 1b00 fOR I : 1 TO N:A{l,I} : A{I,l} t A: NEXT :n2 : 1: GOSUB 2340 1bl0 REft 1b20 REft 1 NORftAliZES lARGETS REAL COftPONENT:l I 1b30 REn 1b40 fOR I : 1 TO N:ft{I} : {X{!}} 2: NEXT :nn : ft{l}:L : l lb50 FOR I : 2 TO N: IF ft{l} I nn THEN L : 1: If ft{I} I nn THEN nn : ft{l} 1bb0 NEXT : fOR I : 1 TO N:T{l,l1} : X{!} I X{L}: NEXT : RETURN 1b70 REM lbBO REn 1 conPUTES INV{Tl REn 1700 fOR I : 1 TO N: fOR J : 1 TO N:A{I,J} : T{l,J}: NEXT J,I: GOSUB 2090 1710 PRINT: PRINT TAB{lO}"THE INVERSE Of THE ftATRIX, INV{T}, IS:" 1720 PRINT

PAGE 123

1730 fOR I : 1 TO N: PRINT TAB{ 10}:: fOR J : 1 TO N: PRINT fN R{AI{I,J}},: NEXT : PRINT : NEXT 1750 REft 17bO REn 11 conPuTES INV{Tl 1 f 1 T 11 1770 REft !760 fOR I = 1 TO N: FOR J : 1 TO N:FUI,J} : 0: NEXT J,I 1790 fOR I : 1 TO N: fOR J : 1 TO N: FOR t : 1 TO N: 1600 F1{l,J} : f1{I,J} + f{I,t} I T{(,J}: NEXT t,J,I 1610 fOR I : 1 TO N: fOR J : 1 TO N:A{I,J} : 0: NEXT J,I fOR I : 1 TO N: fOR J : 1 TO N: fOR t : 1 TO N: 1630 1{I,J} : A{I,J} + AI{I,t} 1 f1{t,J}: NEXT t,J,I 16qO PRINT: PRINT TAB{10}"THE ftODAL DYNAftiCS ftATRIX, INV{T}IftT, IS:": PRINT 1650 fOR I : 1 TON: PRINT TAB{ 10}:: fOR J = 1 TON: PRINT fN R{A{I,J}},: 16b0 NEXT : PRINT : NEXT 1&70 GOTO 1'UO 1&60 REft 1&90 REft I COftPUTES Lft=INV{T}IL I 1'100 REft 1'110 fOR I l TO N: fOR J : 1 TO Nl: fOR [ : l TO N: Gl{I,J} : 'l{I,J} t AI{!,[} I L{t,J}: NEXT [,J,I 1'130 PRINT: PRINT TAB{lO}"ftODAL INPUT ftATRIX, lft:INV{T}Il, IS: ": PRINT 19qo fOR I = l TO N 1950 PRINT TAB{ 10}:: fOR J : l TO N1: PRINT fN R{G1{I,J}},: NEXT : PRINT : 1'1b0 NEXT : GOTO 1970 REft 1960 REft I {OftPUTES (ft:(IT I 19'10 REft fOR I : l TO fOR J : l TO N: fOR t : l TO N: Hl{l,J} : H1{I,J} + ({!,[} 1 T{t,J}: NEXT [,J,I PRINT: PRINT TAB{ lO}"ftODAL OUTPUT ftATRIX, Cft:C1T, IS: ": PRINT fOR I : l PRINT TAB{ 10}:: fOR J : l TON: PRINT fN R{Hl{l,J}},: 20qo NEXT : PRINT : NEXT : END REft REft 2070 REft 1 SUBROUTINE: fiNDS AI=INV{A} 1 REft 20'10 fOR I = l TO N: fOR J : l TO N:U{I,J} = 0: NEXT J,l: 2100 fOR I = l TO N:U{I,I} = l: NEXT 2110 fOR J = l TO {N 1}:ftft = ABS {l{J,J}}:l : J: fOR I : {J + 1} TO N: If ABS {A{I,J}} I nn THEN l : I If ABS {A{I,J}} I nn THEN nn : ABS {l{I,J}} clqo NEXT : If L = J THEN 2150 FOR I : J TO N:Al{I} : l{J,I}: NEXT : fOR P : l TO N:Cl{P} : U{J,P}: NEXT 114

PAGE 124

21b0 fOR I : J TO N:A{J,I} : A{l,I} I A{l,J}: NEXT 2170 fOR P : 1 TO N:U{J,P} : U{L,P} I A{l,J}: NEXT 2180 fOR I : 1 TO N:A(l,I} : A1(I}: NEXT 21,0 fOR P : 1 TO N:U{l,P} : C1{P}: NEXT : GOTO 2220 2200 fOR P : 1 TO N:U{J,P} : U{J,P} I A{J,J}: NEXT 2210 fOR I : N TO J STEP 1:A{J,I} : A{J,I} I A{J,J}: NEXT 2220 FOR p : 1 TO'N: fORt : {J t 1} TO N:U{t,P} : U{t,P} A{t,J} I U{J,P}: 2230 NEXT [,P fOR [ : {J t 1} TO N: fOR n : N TO J STEP 1: 2250 A{[,ft} : A{[,ft} A{t,J} I A{J,ft}: NEXT n.-t,J 22b0 fOR P : 1 TO N:U{N,P} : U{N,P} I A{N,N}: NEXT : fOR P : 1 TO N: 2270 fOR I : {N 1} TO 1 STEP 1 2260 AI{I,P} : U{l,P}: FOR t : {I t U TO N: 22,0 AI{I,P} : AI{l,P} AU,O I AI{t,P}: NEXT t,I,P 2300 RETURN 2310 REft 2320 REft 1 SUBROUTINE: fiNDS EIGENVECTORS 'IVEN ftATRIX l RANt DEfiCENCY ftc 1 2330 REn fOR I : l TO N:Z{I} : I: NEXT : FOR J : l TO {N ftc} 2350 REft 1 fiNDS LARGEST ElEftENT IN PIVOT ROW I CLUftN I TO RI,HT l BELOV 1 23LO nn : ABS {A{J,J}}:R : J:C : J: fOR [ : J TO N: fOR I : J TO N: 2370 If ABS {l{I,t}} c : nn THEN 23,0 2380 R : I:C : [:ftn : ABS {l{R,C}} NEXT I,r: If R : J THEN REn I INTERCHANGES ROWS I fOR I : J TO N:A1{I} : A{J,I}: NEXT : fOR I : J TO N:A{J,I} : A{R,I}: NEXT fOR I : J TO N:A{R,I} : Al{I}: NEXT IF C : J THEN 2500 REn I INTECHAN,ES COLUftNS I fOR I : 1 TO N:A1{I} : A{I,J}: NEXT : fOR I : 1 TO N:A{I,J} : A{I,C}: NEXT A2 : Z{J}:Z{J} : Z{C}:Z{C} : A2: fOR I : l TO N:A{I,C} : A1{I}: NEXT REft I DIVIDES ROV BY PIVOT I ZEROS COLUftN ABOVE I BELOI PIVOT 1 2500 FOR n : N TO J STEP -1:A{J,ft} : A{J,ft} I A{J,J}: NEXT : 2510 fOR [ : {J t 1} TO N 2520 FOR n : N TO J STEP -l:A{[,ft} : A{t,ft} -A{[,J} I A{J,ft}: NEXT n,t: 2530 If J : 1 THEM 2Sb0 fOR [ : 1 TO {J 1}: fOR ft : N TO J 1: 2550 A{[,ft} : A{[,ft} -l{[,J} I A{J,ft}: NEXT ft,[ 25b0 NEXT J:X1{N} : 1: FOR I : 1 TO {N 1}:X1{1} : A{I,N}: NEXT 2570 fOR I : 1 TO N:X{Z{I}} : X1{I}: NEXT : RETURN 115

PAGE 125

RUN ENTER THE ORDER Of THE SYSTEn : q ENTER THE TERnS Of THE If! ftATRil ttBY ROWtt f{ 1 l } : -1500 f{ l 2 } : 750 f{ 1 3 } : 0 f{ 1 } : 0 f{ 2 1 } : 750 f{ 2 2 } : -1500 f{ 2 3 } : 750 f{ 2 } : 0 f{ 3 l } : 0 f{ 3 2 } : ?50 f{ 3 3 } : -1500 f{ ] } : 750 f{ 1 l } : 0 f{ 2 } : 0 f{ 3 } : 500 f{ } : -500 ENTER THE NUnBER Of CONTROLS : 0 ENTER THE NUMBER Of MEASUREftENTS : 0 ll6

PAGE 126

COEffiCIENTS Of CHARACTERISTIC EQUATION IN THE fOLLOWING fORn S"n t COEff{1}S"{n-1} t COEff{2}tS"{n-2} t t COEff{n} COEff NO. 1 2 3 q EIGENVALUES ARE: COEFF VllUE 5000 7500000 3.375Et0'! 2.10'1375Etll REAL InAGINARV bl0.057 0 -1b58.20q 0 -b57.b285 0 -7q,Ul45 0 EI,ENYECTOR nATRIX T, IS: 1 1 1 -.60443 -.'1555 .32579 .lJ'I3'1 .23262 .&5176 -.13651 l THE INVERSE Of THE EIGENVECTOR MATRIX, INV{T}, IS: -.31307 .46337 . 372.75 -.0'15B -.Hir82 .27'163 .3340'1 .37523 .06136 .226U .3137 .55242 THE MODAL DYNAMICS MATRIX, INV{T}fT, IS: -2610.056 .0006 -.00006 .00003 .0015 -1658.205 0 .00021. .000'12 b57.b263 .00002 .OOU -.DODH .00027 -74.11141. 117

PAGE 127

TRANSITION ftATRIX fD IS: 0 D D D 0 D 0 0 0 0 D 0 D 0 0 0 0 0 0 D 0 .90878 D D 0 -3.59218 0 0 0 0 D 0 0 .Ocl95 0 0 0 D D 0 0 ,03736 D 0 0 0 o o o o o o .90678 CONTROL DISTIB. VECTER GD IS: 09H77E -03 l. 7650b3E -02 l.b26bl7E-oq q,s326bE-oq 118

PAGE 128

119 COEffiCIENTS Of IENOftiNATOR l NUftERATOR POLYNOftiALS IN THE fOLLOWING fORft: ncoeff{l}l"{nl} t ncoeff{2}al"{n} t t ncoeff{n} T{Z} : zn t dcoeff{l}aZ"{nl} t dcoeff{2}Z"{n} t t dcoeff{n} coeff # DENOft. COEff NUft. COEff l .l75588L 2 l.35L755E 3 -.3883127 5 -.3883l2L b 7 .1755887 8 POLES ARE: REAL IftAG. .'10878 .'10878 .2M5'J 'JS8b5 -. 'J58b5 .8'1371 ,8,371 -.83252 .553'1, -.83252 -.553'1, ZEROS ARE: REAL lftAG. .8,371 . -.8'1371 .'J58b5 .2845'1 ,58b5 -.83252 .55399 -.83252 .553'1, l 0 UIN IS:

PAGE 129

lrDD REM rr fiNDS BUTTERWORTH POLES IN BOTH Z AND S PLANE u UO REM 120 PI = 130 INPUT THE ORDER Of THE FILTER ";p INPUT "ENTER THE SAMPLE RATE fs ";FS 150 INPUT "ENTER THE PASS BAND tl ";Fl lbO INPUT "ENTER THE STOP BAND f2 ";f2 170 BV = ((f1 + f2l/21/IFS/21 t PI DIM PSRC2tPI,PSII21Pl,PZIC2tPI,PZRC2tPl 190 N = 1 200 REM 210 REM II POLES IN S -PLANE II 220 REM 230 PRINT : LPRINT POLES Of BUTTERWORTH IN S -PLANE : LPRINT FOR K = l TO 2rP 150 PSR = COSC21K1Pl/C21Pl + Pl/{21Pll1BW 2b0 PSI = SINC2lK&PI/l21PI + PI/C2&Pil1BW 270 If PSR < O _THEN : PSRCNI = PSR : PSICNI =PSI : N = N + 1 NEXT K 290 fOR I = 1 TO N -1 300 lPRINT" P";I;" = J";PSICII 310 NEXT I 320 REM 330 REM II PLOTS POLES IN S PLANE II REM 350 INPUT" NTER RETURN TO CONTINUE";A 3b0 CLS: SCREEN 2 : VINDOV 1-3,21-Cl,-2) 370 LINE C -3,01 Cl,OI : LINE <0,21 C0,-21 LINE c-.o5,li-C.o5,ll : LINE C-.05,-li-C.OS,-11 390 LINE C-2,-.0SI C-2,.051 : LINE C-J.,-.QSI C-J.,.D5l LOCATE 1,25 : PRINT "BUTTERWORTH fiLTER POLES" qJ.o LOCATE 2,30 : PRINT "! PUN" q20 fOR I = l TO P : CIRCL[(PSRUJ ,pSJ
PAGE 130

bOO B = 8/2-P blD FOR I = l TO P b20 LPRINT" Pz";I;" = ";PZR(!I;" j";PZI!II b30 NEXT I LPRINT : LPRINT n B = ";B b50 INPUT "PRESS RETURN TO CONTINUE ";A bbO CLS : SCREEN 2 : WINDOW 1-2,1.51 (2,-1.51 : LINE ( -2,01 12,01 b7D LINE (Q,l.SI (Q,-l.Sl fOR I = D TO 3bD STEP b PSET (CU I 700 NEXT I 710 LOCATE 1,30 : PRINT "BUTTERWORTH POLES" 720 LOCATE : PRINT "Z PLANE" 730 fOR I= l TOP : CIRCLE :NEXT I GOTO 121

PAGE 131

ENTER THE ORDER Of THE FILTER 6 ENTER THE SAftPLE RATE fs 20 ENTER THE PASS BAND f1 1.5 ENTER THE STOP BAND f2 3 POLES Of BUTTERWORTH IN S PLANE p 1 : P 2 : p 3 : -.5677313 p : p 5 : p b : -.5677313 p 1 : p 6 : j j .5677312 j j .l3790U j-.1379013 j-.b9327b3 POLES IN THE Z PLANE Pz 1 : j Pz 2 : .5567875 j 39b6l6b Pz 3 : j Pz : ,qbq5207 j 7.73'l'l7q[-02 Pz S : j -1. 739977[ -02 Pz b ,q937713 Pz 7 : .5567876 j-.3'lb6l87 Pz 6 J-.5bS3q 8 : 122

PAGE 132

60 INPUT RENTER NUMBER Of INCREMENTS ft;NWC 61 PI = ,0 Din 8Z(12l, 8Pil2l ,1 Bn = 100 BZ
PI THEN END 3bD GOTO 220 123

PAGE 133

...:r C'
PAGE 134

N = 9 lOS INPUT "ENTER K ";KK 110 DEfDBl D,s,R,Q,B 120 DIM NlNJ,DlN>.B 155 NEXT I 170 U = 1: fOR I = l TO N: If ABS !Nm I > -00001 THEN 190 160 l1 = 1 + I: NEXT fOR I = 1 TO fN Ul :D(ll = NU LU I NfLU: NEXT : N = N -u: 200 If N = 0 THEN 250 210 PRINT TABt lOI" ROOTS ARE :": 220 PRINT TAB< 101" REAl. !MAG" GOSUB 310: If Q = l THEN PRINT TAB< 101" 0 0" GOTO 270 250 If Q = 0 THEN PRINT TAB< lOI"NO ROOTS." 2b0 PRINT TAB< 10I"ONE ROOT AT (0,01." 270 ENJ 260 REn REM I SUBROUTIN: fiNDS ROOTS Of POLYNOMIAL USING BIRSTOV'S AlGORTIHM t 300 REn 310 J = 1:n = o:E = .oooo1:Ro = .s:so = o 320 l = O:R2 = Ro:s2 = SO:nl = N 2 I n: If Ml = 2 THEN 350 330 If M1 > 2 THEN 3b0 340 RlJl = Dl1l:IlJl = o: GOTO 530 350 R2 = Dlll :s2 = D !21: GOSUB 550: GOTO 530 3b0 BlOl = 1!Blll = Dlll -R2: fOR ( = 2 TO M1! 370 8(Kl = P [ OR ABS E THEN 460 GOSUB SSO: GOTO 510 490 If l > N I SO THEN PRINT TAB< lOI"POLYNOMAIL SOVER, biDK'T END 500 L = l + 1: GOTO 3b0 SlO n = n + l!J = J + 2:N3 = N -2 r n: fORK= l TO N3:b
PAGE 135

530 DEf fN RCSl = INT ClOOOOO!IS + .511100000! FORK= l TON: PRINT lOl; fH fN RCllKll: NEXT K: RETURN 550 P, = R2 I R2 q I S2: If P, > = 0 THEN 560 SbO Rl = SQR C P,l:RtJl = R2 1 2:RtJ + ll = -R2 I 2:ICJl = Rl I 2: 570 HJ + ll = Rl I 2: RETURN Rl = SQR CP,l:ICJl = O:ItJ + ll = O:R(Jl = lRlR2l I 2: 5,0 RCJ + ll = ( -Rl R2l I 2: RETURN 126

PAGE 136

100 N ; 110 INPUT "ENTER K 120 DEfDBL D,S,R,Q,B 130 DIM NINJ,DINJ,BINl,PINl,QINJ,RINl,I(Nl,NCINJ,DClNl BS Ba .= S-326bq3E -OS 1qo N((21 = 2 737l3q[ -02 : NCOI = -23B6bb : N((qJ = 150 NCISI = : NC(bl = 2.,qs2q : NC(7) = 3lbD7b : NCI6l = 2-105112 lbO = 727036 : NCUOI = -NCI,I : NCIU) = -NCI61 : NCI121 ; -NCI7l 170 NCn3l = -NCibl : NCml = NCI5l : NCU51 ; -N((q) : NCUbl ; -RCI31 160 NCI171 = -NC(2l : NCill = 0 : DCI1l = l 1,0 DC<21 = -q.20A3b2: DCI3l = 6-Sb3265 : DCm = -U.qS37, : DCI51 = l2.3qbb2 200 DC 2 THEN q]O q10 R(JI = Dlll:l(Jl = o: GOTO bOD q20 R2 = Dlli:S2 = Dl2l: GOSUB b20: 'OTO bOO q30 BIOJ = l:Bill = D(ll -R2: FORK= 2 TO Ml: qqo BIKI = DIKl -R2 I 811: -ll S2 I BIK-21: NEXT qso PIOl = Q:Pill = 1: fOR K = 2 TO Ml: qbO P IKI = -8 IK 11 -R2 rP IK ll S2 I P (K 21: NEXT q70 Q(OI = O:Q(ll = 0: FOR K = 2 TO M1: q&Q Q(Kl = B + R2 1 QIMl-l>:Db = Rq 1 ISRs r sq: 520 Rb = IS3 l R5 SS I R31 I Db 530 Sb = IR3 I sq Rq I S3l I Db:R2 = R2 + Rb:S2 = S2 Sb: sqo If ABS E OR ABS E THEN SbO 550 GOSUB b20: COTO 560 127

PAGE 137

5b0 If l > N I 50 THEN PRINT TABt lOl"POLYNOMAIL iOVER DIDN'T CONVER,E": END 570 L = L + 1: GOTO 560 n = n + l:J = J + 2:N3 = N2 In: fORt= 1 TO = BlKJ: HEXT : 5,0 'OTO 3,0 bOD DEF fN RtSJ = INT UDODDDIIi + .SJ/100000_! L10 fORK= 1 TON: PRINT TAB< lOJ; fN R = 0 THEN LSD L30 R1 = SQR (P,I:RtJl = R2 I 2:RtJ + 11 = R2 I 2:ItJl = Rl I 2: LqQ ItJ + 11 = R1 I 2: RETURN LSD Rl = SQR
PAGE 138

N = 21 129 INPUT RENTER K n;KK 120 DEFDBL J,s,R,Q,B 130 DIM N = Dnl Rc: FORK = 2 TO n1: qqo BCKI = DCKI -R2 l BCK-11 S2 I BCK-21: NEXT qso PCOl = Q:PCll = 1: fOR K = 2 TO n1: qbo PCKI = BCK-11 -R2 I P E OR ABS E THEN SbD SSO GOSUB b2D: GOTO 560

PAGE 139

5b0 If l ) N r 50 THEN PRINT TAB( lDl"POLYNOMAIL SOVER DIDN'T CONVERGE": END 570 l = l + l: GOTO q30 n = ft + l:J = J + 2:N3 = N2 In: FORK= l TO H3:D(KI = BlKI: NEXT : GOTO bOD DEF fN R(SI = INT llDDDOD!IS + .SIIlDDDOD! blO FORK= l TON: PRINT TAB( lOJ; FN R(R(KJJ;" "FN Rli(KJI: NEXT K: RETURN b20 = R2 I R2 q I S2: IF P, ) = 0 THEN LSD b30 Rl = SQR (P,I:R(JI = R2 I 2:R(J + 11 = R2 I 2:I(JI = Rl I 2: bqQ I(J + ll = Rl I 2: RETURN bSD. Rl = SQR
PAGE 140

APPENDIX B COMPUTER SIMULATION PROGRAMS The programs contained in Appendix B are written in both BASIC and Pascal and were'used to simulate and test the different control methods for this thesis.

PAGE 141

10 INPUT ftTime of simulation desired ft 20 OPEN wow,tl,ftSiftDATw 100 Din ,6h 8 (6), C lbl X 161 XNEW (f. I 110 All,11 = -.2Sbb& : Atl,21 = : All,31 = : = -.0073 120 All,Sl = : All,bl = 3-65011 : A(l,71 = b : All,6l = .&0&7b 130 At2,1l = : A(2,21 = -.lS% : H2,3l = : = .10M2 Al2,5l = 3 : A<2,bl = -23-57677 : A(2,7l = q.36929 : A12,6l = 7.3&b37 150 Al3,ll = : Al3,2l = : At3,3l = -.l6blq : = .b06b6 1b0 At3,Sl = b : At3,bl = q.]6929 : At3,71 = : At3,&l = 170 A(q,}J = .ooq67 : A(q,2J = -07255 : Alq,JJ = .qos79 : A(q,qJ = .Slbb2 160 A (q,sJ = -53916 : A(q,bJ = q.12q25 : A(q,?J = 6-952b79 : Atq,6l = -}q.qq]5b 190 A(5,1l = : Al5,2l = -OlOb : AC5,3l = .QQlll : A1S,qJ = .QOODS 200 AI5,Sl =-.256b6 : A = 2bD U6,51 = .ooq&7 : AC&,bl = .07255 : AC&,7l = .qQS79 : Al6,&l = .Slbb2 270 REn 260 Bill = .ooooq : B<2> = .ooo11 : BC31 = .oo&1l : Blq> = -02737 210 8(5) = o : scb> = .oooo1 : am = .noon : B<6J = .ooo75 300 REft 310 fOR I = l TO 6 320 xm = o 330 NEXT I 3qo REM 350 T = 0 3b0 ft: = 0 370 If T = 0 THEN ft: = -10 360 fOR I = 1 TO 6 390 XNEWIII = 0 qoo FOR J = 1 TO 6 q}Q lNEWIII = XNEWlll + AtJ,JltX!Jl q20 NEXT J q3o XNEW
PAGE 142

100 INPUT "Time of simulation desired ;TFIN 110 OPEN "O",a1,"MODESin" Un A<8,81, 8(6) 1 ((61 1 XC61, XIIOH6), 130 fOR I = l TO 6 lqO fOR J = 1 TO 6 150 Asx sqo NEXT I 133

PAGE 143

550 REft 5b0 REft II IS TIP IS TIP VELOCITY II 570 REn SAD PRINT REn bOD FOR I = l TO A blD X(ll = XNW(II b20 NEXT I b30 T = T + .os If T < TfiN THEN bSD CLOSE bbD END 134

PAGE 144

100 INPUT "Time of simulation desired" ;TFIN 110 OPEN 120 DIM Al6,6l, 8(6), CC6l XC61, XNEWt6l, 130 FOR I = l TO 6 fOR J = l TO 6 150 A
PAGE 145

100 INPUT RTirne of simulation desired n ;TfiN 110 OPEN Roft,a1,ftSIMfiLTn 120 DIM ACA ,AI, 8 !AI, C !AI X !AI XNEW CAl, Xf (6), Yf!Al 130 ACl,ll = -.25Ab6 : All,2l = .4bb02 : All,3l = : All,qJ = -.0073 AC1,51 = : Al1,bl = 3.65011 : All,7l = b.30602 : All,61 = .6067b 150 Al2,ll = : A<2,2l = : Al2,3l = : Al2,ql = .10662 1b0 AC2,5l = 3-65011 : AC2,Ll = -23.57677 : A(2,7l = : Al2,6l = 7-36b37 170 Al3,1l = : AC3,21 = : AC3,Jl = -.l6b14 : ACJ,q) = .b06b6 lAD A<3,5l = b-30602 : A<3,bl = : A(3,7l = : AB,Al = A(q,ll = : A(q,2l = -07255 : = : A
PAGE 146

hlO NEXT J beD = XNW(Il + BCilrfK b30 NEXT I LqQ REn LSD REn II RECURSIVE BUTTERWORTH fiLTER II bbO REn b70 N = INT(TI20l b60 If H > A THN H = 6 b,O XfCNI = XCq) 700 y = 0 710 fOR [ = 0 TO N 720 Y = Y + BZCKIIXfCN-Kl 730 NEXT K 7qofOR K = l TON 750 Y = Y BPCKIIYfCN-KI 7b0 NEXT [ 770 REn 760 REn II Y IS nODAL TIP VELOCITYfOUND BY THE fiLTER II 7,0 REn 600 PRINT ll,T,Y AlD REn 620 REn II SET NEW VALUES II 630 REn 6qo fOR I = l TO A 650 X(!) = XNEV(l) AbO NEXT I A70 YHNI = Y A60 If H < 6 THEN ,30 6,0 FOR I = l TO N ,DO YfCI-ll = YfCII ,10 XFCI-ll = Xf(Il ,20 NEXT I ,30 T = T + .os ,qo If T < TfiN THEN 550 %0 CLOSE ,bO END 137

PAGE 147

100 INPUT of simulation desired 110 OPEN o,a1,sinFilT2. 120 Dift A<6,6l, B<6J, ({6), XC6J, XHEWl6l, Xfll2l, Yf(l2J, BZ
PAGE 148

blD XNEW 12 THEN N = 12 720 Xf(NJ = X(q) 730 y = 0 7qo fOR K = 0 TO N 7SO Y = Y + BZ(KliXf(H-KJ 7LO NEXT K 170 fOR K = l TO N 760 Y & Y BP
PAGE 149

10 INPUT time of simulation desired 15 INPUT Enter gain K ;KK 20 OPEN o,al,VEtfBK10 100 Din ACA,&J, BCAJ, CC4J XC6J, XNEWC4J 110 ACl,1J = -.25Ab& : A = 0 qoo fOR J = 1 TO 6 XNEV
PAGE 150

100 INPUT time of simulation desired ;TFIK llD OPEN ft0ft,a1,DAnPlft 120 Dlft A<6,6J, 8(6), Cl6l Xl6), XNEWI6), Tft(q,qJ 130 fOR I = 1 TO 6 1q0 fOR J = l TO 6 150 A = .aoqa, : 8<31 = -.a1D3b : 8
PAGE 151

SSD REn SbD REn II xq IS TIP IISPLACEnHT, yq IS TIP VELOCITY, YnOE IS TH II Sbl REft II ftOIAl VELOCITY Of THE ftAJOR ftOOE AT THE TIP II 570 REn PRINT al,r,xq,yq,vnonE REn bOO fOR I R l TO blD Xll) = XNEWli) L20 NEXT I L3D T = T + DS If T < TfiH THEN LSD CLOSE bbD END 142

PAGE 152

100 INPUT "Time of simulation desired :TfiN 110 OPEN "0",#1,"DAftP12" 120 Dlft A{A,a}, 8{8}, ({6} X{6}, XNEV{6}, 130 fOR I : l TO 6 fOR J : 1 TO 6 150 A{I,J} : 0 1LO NEXT 170 NEXT I lAO A{l,l} : -.63252 : A{l,S} : A{2,2} : -,qqabS : A{2,b} : 200 A{3,3} : -.Ob&q7 : A{3,7} : l1.b0528 210 : : : 220 A{S,l} : .OlQ&q : A{S,S} : .&3252 230 A{b,2} : : A{b,b} : A{7,3}: .017bS : A{7,7}: .Sb512 250 A{6,q}: .03bql : A{A,A}: 2b0 REft 270 8{1} : : 8{2} : : 8{3} : : 8{q} : .013q1 260 8{5} : .OOOOL : 8{b} : .0001L : 8{7} : .00016 : 8{6} : .00037 REft 300 T,1} : .b7Sbq : T,2} : l : T,3} : : T{l,-q} : 310 1,1} = 1 : t,2} = -.llo1q : 1,3} = 1 : = 320 T,1} = -.aoqq3: 1{3,2} = : T,3} = .23262 : = .65178 330 .110L2: T{q,2}: ,q12q,: T{q,3}: -.73651 : 1 REn 350 fOR I : 1 TO 6 3b0 X{I} : 0 370 NEXT I 360 REft T : 0 f[ : 0 q10 If T : 0 THEN f[ : q20 fOR I : 1 TO a q30 XNEV{I} : 0 fOR J : 1 TO a q50 XNEV{I} : XNEV{I} t A{I,J}tX{J} qto NEXT J q70 XNEV{l} : XNEV{I} t B{I}rf( NEXT I REn 5oo yq = o : xq = o 510 fOR I : 1 TO 520 yq : yq t T{q,I}rX{I} 530 xq : xq t T{q,IJX{Itq} 5qo NEXT I 143

PAGE 153

550 REn 5b0 REft ir IS TIP DISPliCEftENT, IS TIP VELOCITY 11 570 REft 560 PRINT #l,T,Xq,yq,YftODE 5,0 REn LOO fOR I : l TO 6 LlO X{l} : XNEV{I} NEXT I b30 T : T t .OS If T c TfiH THEN qoo LSD CLOSE LbO END 144

PAGE 154

10 INPUT "Time of simulation desired :TfiN 145 15 INPUT "Enter gain ( ":tt 20 OPEN "0",#1,"VElf8t10" lOU Din l{8,8}, 8{8}, C{8} X{8}, 110 l,1} = -.256b6 : A,2} : ,qbb02 : A,3} : : : -.0073 120 1{1,5} : 1{l,b}: 3.65011 : l{1,7}: b.30602 : l{1,6}: .6067b 130 A,1} : : l{2,2} : -.15,b : A{2,3} : : : .10882 l{2,5} : 3.65011 : l{2,b} : -23.57877 : 1{2,7} : : l{2,8} : 7.38637. 15U 1{3,1} : .0,,0, : l{3,2} : : A{3,3} : : : .60868 1b0 l{3,5} : 6.30802: l{3,b} : l{3,7} : : l{3,6} : 170 : : : .07255 : : : : .51b62 180 : .53,16: : q,,2425 : : 8.,5267, : : 1,0 l{5,1} : .02523 : A{5,2} : .OlOL : A{5,3} : .00111 : A{S,q} : .DODOS 200 l{S,S} :-.25668 : l{S,b} : ,qb602 : l{5,7} : .0,,0, : A{5,8} : .0073 210 l{6,1} : .01Db : l{b,2} : .02b3q : A{b,l} : : l{b,q} : .0011, 220 1{6,5} : ,qbb02 : l{b,b} : -.15,L : l{b,7} : ,q708, : l{b,8} : .10882 230 1{7,1} : .00111 : l{7,2} : .01064 : l,3} : .02602 : l{7,q} : .01217 240 1{7,5} : .0,,0, : A{7,b} : : l,7} : -.16b1q : 1{7,8} : .b08b& 250 l{8,1} : : A,2} : .0007, : l,3} : .00811 : : ,Oq10b 2LO A,5} = .ooq57 : A,bl : .07255 : A,7} = : A,8} = .51662 270 REn 280 8{1} : 3.bb013E-05: 8{2} : .0007,21: 8{3} : 8.113301E-03: : .0273713 2,0 8{5} : 2.37E-07 : 8{b}: 8{7} : 8{8}: 300 REn 310 fOR I : 1 TO 8 320 X{I} : 0 330 NEXT I 340 REft 350 T : 0 3b0 f[ : 0 370 If T : 0 THEN f( : -10 380 fOR I : 1 TO 6 3,0 : 0 qoo fOR J : 1 TO 8 q10 XNEW{I} : XNEW{I} t l{I,JJX{J} q20 NEXT J q30 XHEV{I} : XNEW{I} t B{I}tf[ NEXT I REft qq5 REft 11 X{8} IS TIP DISPLACEftENT, X{4} IS TIP VELOCITY 11 REft 450 PRINT #l,T,X{6},X{4} 455 REn 4b0 fOR I : l TO 8 470 X{I} : XNEW{l} 480 NEXT I 4,U T : T t .OS SOU If T < TfiN THEN 3LO 505 ClOSE 510 END

PAGE 155

100 INPUT "Time of simulation desired :TfiN 110 OPEN Din A{6,8}, 8{8}, ({8} X{8}, XNEV{8}, Xf{6}, Yf{6}, [{6} 130 A{l,l} : .25666 : A{1,2} : : A{l,3} : .0,909 : : -.0073 A{l,S} : .88679: A{l,b} : 3.85011 : A,7} : b.30802 : A{l,8} : .6067b 150 A{2,1} : : A{2,2} : -.159b : A{2,3} : : : .10882 lbO A{2,S} : 3.65011 : A{2,L} : .57677 : A{2,7} : : t{2,8} : 7.38b37 170 A{3,l} : .09909 : A{3,2} : : A{3,3} : : : .b08b8 180 A{3,5} : b.30602 : A{3,b}: : A{3,7} : A{3,6}: 190 : : : .07255 : : : : .516b2 200 .53918: 6.952679: 210 l{S,l} : .D2523 : l{5,2} : .01Db : A{S,3} : .D0111 : : .DDD05 A{5,5} :-,258L8 : A{S,L} : : l{S,7} : .0990,: A{5,8}: .OD73 230 A{L,l} : .OIOL : A{L,2} : : A{b,3} : : : .00119 A{L,S} : i{b,L} : -.1596: A{b,7} : A{b,8} : .10882 2SD l{7,1} = .00111 : A,2} : : A,3} : .02602 : : .01217 2b0 A{7,S}: .09909: A{7,b} : ,q7Q89 : A{7,1} = : A,8} = .L08b8 27D A{8,l} : : A{6,2} : .00079 : l{6,3} : .00611 : : 280 A{li,S} : : A{6,b} : .07255 : A,1} : : A{8,6} : .51bb2 290 REft 300 REft II INPUT VECTOR (8] II 31D REn 320 8{1} : 3.LbOBqE-05 : 8{2} : : 8{3} : 8.1B2b9E 330 8{5} : 2.370377E-07 : B{L} : b.9b2379EOb 8{7} : : 8{8}.: 350 REn 3b0 REft II BUTTERWORTH filTER ZEROS II 370 REn 360 8ft : 390 Bl{O} : 1 : 8l{l} : 8 : 81{2} : 211 : BZ{3} : Sb : 8l{n : 70 'IDD BZ{5} : 81{3} : BZ{L} : 8Z{2} : BZ{7} : 8Z{l} : 81{8} : BZ{O} fOR I : 0 TO 6 : BZ{I} : BZ{I}rBn : NEXT I '120 REn REII n BUTTERWORTH FILTER POlES u REn BP{l}: '1.363951: BP{2}: : 8P{3}: .82lb0b# qbQ.BP{q} : a.q89605: BP{5} : : BP{b} : 1,q50?66 BP} : : BP} : .02q111 REft n = q,s 500 REft 510 REft II VAlUES II 520 REft 53D FOR I : 1 TO 6 X{I} : 0 550 NEXT I 5b0 REn 570 T : 0 560 REft 590 R[ft II BEGIN SIMUlATION II bOO REn 146

PAGE 156

blO n: = o b20 If T : 0 THEN ft : -LO b30 fOR I : l TO 6 XNEW{I} : 0 bSO fOR J : 1 TO 6 bbO XNEW{I} : XNEW{I} 1 A{I,J}rX{J} b10 NEXT J b60 XNEV{I} : XNEW{I} 1 B{l}rf( (liB{l}IY L'IO NEXT I 100 REII 110 REII II RECURSIVE BUTTERWORTH fiLTER II 120 REII 730 N : INT{TI20} If N t 6 THEN N : 6 150 Xf{H} : 1b0 y : 0 170 fOR C : 0 TO H 160 y : y I BZ{t}IXf{N[} 7'10 NEXT C 600 fOR t : 1 TO H 610 Y : Y BP{C}IYf{N-C} 620 NEXT t 630 REft REft 11 Y IS ftODAL TIP VELOCITY fOUND BY THE fiLTER, II 650 REII u IS THE TIP YELO_(ITY lND X{6} IS THE TIP DISPLlCEftENT n 6b0 REII 610 PRINT #1, T, Y 660 REII 6'10 REII II SET NEV VALUES II '100 REII ,10 fOR I : 1 TO 6 ,20 X{I} : XHEV{I} ,30 NEXT I ,qo Yf{N} : Y ,SO If H < 6 THEN 1000 ,bO fOR I : 1 TO N ,70 Yf{I-1} : Yf{l} ,60 Xf{l-1} : Xf{l} ,'10 NEXT I 1000 T : T I .05 1010 If T < TfiN THEN blO 1020 CLOSE 1030 END 147

PAGE 157

100 INPUT "Time of sllllllation desired" ;TfiN 110 OPEN "O",al,"fllfDBK-1&2" 120 Din A(6,6l, 8(6), ((6J, Xt6l, XHW(6J, XflAl, Yf(6), BZl&l, BP
PAGE 158

blD BZ2(q) = 2Ab.l777 : BZ2C5l = : BZ2Cbl = -q2A.0356 b20 BZ2(7) = : BZ2W = l9D.Sb : BZ2m. = 5Al b3D BZ2(l0l = -2l-A7707 : BZ2Clll = -.1020379 : BZ2ll2l = q bqQ fOR I = 0 TO 12 : BZ2liJ = BZ2Cil121BM2 : NEXT I b50 REft bbD REM II 1st ftODAl GAIN IS K1, THE 2nd MODAL 'AIN IS K2 II L70 REft LAD Kl = q.s: K2 = 5.7 REft 700 REft II BEGIN SIMULATION n 710 REM 720 fOR I = 1 TO A 730 xm = o 7qo NEXT I 750 REM 7LDT=O 770 fK = D ?AD If T = D THEN fK = -10 790 FOR I = 1 TO a ADD XNEW(I) = 0 610 fOR J = l TO a 620 XNEWCII = XNEWCII + ACI,JliX(J) a3D NEXT J aqo XNEV a THN N = a 910 Xf(NJ = XHEW(q) 920 Y = D ,30 fOR K = D TO N ,qo Y = Y + BZCKJIXFCN-KJ ,50 NEXT K 9b0 fOR K = 1 TO N Y = Y BPCKJaYf(H-KJ ,60 NEXT K REM 1000 REM n RECURSIVE BUTTERWORTH BAND PASS fiLTER FOR THE 2nd MODE n 1010 REM 1020 N2 = INTCTI201 1030 If N2 > 12 THEN N2 = 12 Joqo Xf2(N21 = XHEW(q) 1050 V2 = D lObO FOR K = n TO N2 1070 V2 = V2 + BZ2lKJzXF2CN2-KI lOaD NEXT K 1090 FOR K = l TO N2 1100 V2 = Y2BP2(Kl1YF2lN2-Kl UlO NEXT K 149

PAGE 159

U20 REM 1130 REM II Y IS nODAl TIP VELOCITY Of THE fiRST nODE fOUND BY THE fiLTER II llqo REM It V2 IS THE MODAL TIP VELOCITY fOR THE SECOND MODE. X(ql IS It 1150 REM II THE TIP VELOCITY AND XC6J IS THE TIP DISPLACfnENT II ULO REM U70 PRINT 11, T, y, V2, X(q), XC6J U60 REM 11,0 REM II SET NEW VALUES fOR. THES STATE VARIABLES II 1200 REM 1210 fOR I = l TO 6 1220 X
PAGE 160

PROGRAM URMA: CONST td : 0.05 : YlR ( : integer: t,tf,yd,Uin : real: A,B,U,Y : array[l 6] of real: UTA: text: BEGIN {1 INITIAliZE PARAMETERS r} Afll :: 0.1755611 AI2J Af31 :: -0.3863131 A(q) :: 0.371q5q, A(5J :: -0.3663132 AlLI :: Al71 :: 0.17556,3 Al&l :: l fOR { :: l TO 6 DO BEGIN U((J :: 0 : Y(t)::O: END: t :: O: : 8(11 :: : : 8121 := : : 8131 := 1.351161E-c : : :: -b.b,06,1E-3 : : 8[51 :: b.b,061,E-3 : : Blbl :: -1.35,,61E-c : : 8[71 :: -l.cqlS&SE-2 : : 8161 := : VRITE{'ENTER THE DESIRED TinE Of THE SIMULATION:'}: READLN{tf}: ASSIGN{DATA, 'ARMASift'}: REVRIT[{DATA}: REPEAT {1 Begin Simulation Routine 1} yd :: 0 : If t : 0 THEN Uin :: -10 ElSE Uln :: O: fORt:: l TO 6 DO yd :: yd B!t)IU[,() l[(JtY{,(J: FOR ( :: 1 TO 1 DO {1 Reset the values for U{k} and Y{k} 1} BEGIN U!tl :: U(CtlJ: Y!tl :: Y!ttll: END: Y(6) :: yd: U(6J := Uin: ',yd}: t :: t I td; UNTil tf c t : END. 151

PAGE 161

PROGRAft DARftAlstmode: CONST td = 0. OS : gain : 5.328bq]-5: VAR [ : integer: t,tf,yd,Uin : real: A,B,U,Y : array(l lbl of real: UTA: text: BEGIN {1 INITIALIZE PARAftETERS 1} lUI := -q,C!o&3bc : BUl . UcJ -6.5b3265 : 8[2) 0.2H6bb -. U31 :: :am o.a?cmu . : sm -.USJ :: : 8[5) Ubi 13.61107 : Blbl 3.0lb07b .-.l[1) :: -lq,25103 : 8(1) 2.105112 U8J : 8[8) -0.127038 .-. :: : :: -8[6) ulOI : = 12.Mb0. 6 : BUOI : = -8(1) Allll :: : 8[111 :: -B[bJ All2J :: : BU2J :: -8151 um := : Bl13J := -am Ami:= 1.qusss : sm1 : = -sm ll151 :: -0.275"1123 : 8(151 :: -8{21 A[lb] :: 0.024U7 : BUbl :: -8(1) fOR [ :: 1 TO lb DO BEGIN 8([) :: gainBitl: U([J :: D : Y([) :: 0 : END: t :: 0: VRITE{'ENTER THE DESIRED TinE Of THE SiftULATION:'}: READLN{tf}: ASSIGN{DATA, 'ARftASift.ftDl'}: REVRITE{DATU: REPEAT { Begin Simulation Routine 1} yd :: 0 : If t : 0 THEN Uin := -10 ELSE Uin :: D: FOR [ :: l TO lb DO yd := yd t BlrlUIU-Cl l([JYU7-Cl: 152

PAGE 162

FOR [ :: l TO 15 DO {t Keget the valueg for U{k} and Y{k} t} BEGIN Ultl :: Ul(tll: Y(t) :: Y[[tlJ: END: YU!.J :: yd! UU!.l :: Utn: ',yd}: t :: t I td: UNTIL tf c t : END. 153

PAGE 163

PROGRAM ADAPT_EST_USING_PROJECTION: CONST td:O.OS: epsilon : 0.0001 : c : O.OOL : :const : O.S : TYPE vector : array[1 8) of real: matrix : array[l 6,1 81 of real: VAR J, t : integer: t, tf, yd, Ydplusl, fd, Uin, Error, Denom, Ydsys : real: AS : matrix: A, 8, U, Y, BS, X, Xnew: vector: DATA: text: BEGIN {1 INITIALIZE PARlftETERS 1} All) :: 0.1 u21 := o.q l(3] :: -0.3 :: O.l l(5J :: -0.3 Ut.l :: A[1] :: 0.1 A(6J :: lo : 8[11 :: 2.?E-2 {1 : 8[21 :: 1.2E-2: {1 : 8(3] :: 1.3E-2 : {1 : em := -t..bE-1 : : 8[51 :: la.bE-3 : : B[bl :: -lo.lE-2 : : 8{71 := -1.2E-2 : : 8[8) :: -2.7E-2 : DARftA ftODEl 1} AND THE ESTiftATOR 1} ftODEl ARE IDENTICAL 1} fOR THE SISO CASE 1} AS{lo,1):: -0.256b8: AS[1,2J:: o,qbb02: lSll,lJ:: AS{1,q):: -0.0073: AS[1,5):: -2,.861.7,: AS[1,b):: 3.85011: AS[1,7]:: b.308D2: AS[1,8]:: 0.60876: AS[2,1J:: o,qbb02: AS[2,2]:: -0.15,b: AS[2,l]:: AS{2,q):: 0.108e2: ASt2,5J:: 3.65011: lS[2,bJ:: -23.5767?: AS[2,?):: AS[2,8J:: ?.36b37: AS[3,1J:: 0,0,,0,: AS(],2J:: o.q706't: lSH,lJ:: -0.16b.J.q: AS(],q):: O.b0flb6: AS{3,SJ:: b.30802: AS[3,b):: AS[3,7J:: -zq.'tb2Sq: AS(3,8):: l],qc,Oc: o.ooqa7: AS[q,zJ:: 0.07255: AS(q,]J:: Q,qOS7't: AS[q,q):: O.Slbbc: 0.53,18: AS[q,LJ:: AS(q,?J:: AS[q,8):: -.J,q,qqJSb: AS[5,l):: 0.02523: AS[5,2J:: O.OlOb: AS[5,3J:: O.OOLH: AS-(S,q):: 0.00005: AS[5,SJ:: -0.256b6: AS[S,bJ:: O,qbb02: AS[5,7J:: AS[5,6l:= 0.0073: AS[b,1J:: D.010b: AS[b,2J:: 0.02b3q: AS[b,lJ:: O.OlObq: ASib,q):: 0.0011,: AS{b,SJ:: o,qbb02: AS[b,bJ:: -0.15,b: AS[b,?J:: o.q7061: AS[b,6):: 0.10682: AS[7,1J:: 0.00111: AS(7,2):: O.OlObq: AS[7,3):: 0.02b02: AS[7,q):: 0.01217: AS[7,SJ:: O.O't,O,: AS(7,b]:: o,q706't: AS(1,7J:: -O.l8b1q: AS[7,6):: O b06b8: ASt6,1J:= o.ooooq: ASt6,2J:= o.ooo1,: asta.lJ:= o.ooa11: o.oq.J.ob: Asta,sJ:= o.ooqa?: o.o7255: ASl6,7l:= o.5lbb2: BS[l):: 3.bb013qE-5: BS(2J:: 7.'t2l027E-q: BSnJ:: 8.U32b,E-3: BS(q):: 2.73713qE-2: BSlSJ:: 2.370377E-7: BS[b):: 8S[7J:: BSl8J:: 154

PAGE 164

fOR [ :: a DO BEGIN u ([] :: 0 ; Y([) :: 0 : xm := o END: t :: 0: U In :: O: WRITE{'ENTER THE DESIRED TIME Of THE SIMULATION:'}: READLN{tt}: ASSIGN{DATA,'ADAPTEST.PRJ"}: REWRITHDAT U: REPEAT { Begin Sirrulation Routine 1} {1 Sirrulate the System using State ftatrix Equations 1} If t = 0 THEN fd := ELSE fd :: O: fOR J DO BEGIN XNEV(JJ :: 0; fOR { :: DO XNEV(JJ :: XNEV(JJ t ASlJ,[)IX[tJ: XHEVfJJ :: XNEV(JJ t BSlJIfd t BSJJ1Uin: END: yd :: {1 Update Theta --the estimator's parameters 1} Ydsys := 0 : fOR [ :: Ydsys :: Ydsys t B[C)tU(,-[) A{[)IY['Hl Error :: yd Ydsys: Denom :: c : fOR a DO Denom :: Denom t U([J U(t::J t Y(C)IY[C) : fORt:: := 1 TO 6 DO BEGIN Uti : = A.m aconstYl,-tlError : 8[t::) : = BUI aconstUI,-t::l Error/Denom : END: {1 System DARMA Step Ahead Estimator Routine 1} fdplusl :: ydAlll: fOR[ :: 2 DO :: B!t::JUllO-tl A(t::)IY(lO-t::l: 155

PAGE 165

fOR ( :: l TO 7 DO {I Reset the values for andY{(} 1} BEGIN Ultl :: U([tlJ: Y[[J :: Y([tll: END: Yl61 :: yd: Ul81 :: fd Uln: fOR [ :: l TO 8 DO {I Reset the values for I} Xltl :: XNEVltl: RHElN{DATA, t: q:2,' ',yd: ', Ydplusl: ',Error}: t :: t f td: UNTIL {tf epsilon < t}:
PAGE 166

PROGRAM UAPT _EST _nodi f ied_Least_Squares: CONST td : 0.05 : epsilon: 0.0001 : TYPE vector : array(l 6) of real: vector2: array(1 1b) of real: matrix : array(l 6,l 6) of real: matrix2 : array(1 lb,l lbl of real: YAR J, ( : integer: t, tf, yd, Ydplusl, fd, Uin, Error, Denom, Ydsys : real: AS : matrix: P : matrix2: A, B, U, Y, BS, X, Xnew SigP, Psig, Theta, Sigma DATA BEGIN {1 INITIALIZE PARAMETERS 1} : vector: : vector2: : ted: AUJ :: 0.1 ll21 :: U3J :: -0.3 A(q) :: 0.3 A[5) :: -0.3 Albl : = o.q A[7) :: 0.1 A(6) :: l : Blll := 2.7E {1 DARftA ftODEL 1} : 8[21 := l.2E-2 : {1 which is the transfer 1} : 8(3] :: l.3E-2 : {1 function of the system 1} : B(q) :: -b.bE-3 : : 8{51 :: b.bE-3 : : 8[b) :: -l.3E-2 : : 8(71 :: -l.2E: : 8[6) :: -2.7E-2: AS[l,1J:: ASU,2J:: o,qbb02: ASU,3J:: AS[l,qJ:: -0.0073: AS{l,S]:: AS(l,b):: 3.65011: AS[l,7):: b.30602: AS[1,8):: 0.6087b: AS[2,1):: o,qbb02: AS[2,2]:: lSt2,3J:: AS(c,q):: 0.10882: AS{2,5]:: 3.65011: AS[2,b):: -23.57877: AS{2,7):: AS[2,8):: 7.38b37: ASt3,l):: AS(3,2):: AS(3,3J:: -0.18blq: AS(3,q):: O.b08b8: lS(3,5):: b.30802: AS(3,b]:: AS(3,7):: AS(3,6):: AS{q,l):: O.ooqa7: AS[q,2):: 0.07255: AS(q,3):: AS(q,q):: O.Slbb2: lS[q,S):: lS[q,b):: q,12q2s: AS[q,7):: lS(q,a):: -1q,qqJSb: ASI5,1):: 0.02523: AS(5,2):: O.OlOb: lS(5,3):: 0.00111: AS[S,q):: 0.00005: AS{S,SJ:: ASIS,LJ:: Q,qbb02: AS[5,7J:: 0.01101: AS(5,8J:: 0.0073: AS[b,1J:: 0.010L: ASIL,2J:: 0.02b3q: AS!b,3J:: 0.010bq: AS(b,q);: AS[b,SJ:: AS(b,bl:= AS{b,7J:: o.q7081: ASib,8J:: 0.10882: ASl7,1J:: 0.00111: AS[7,2J:: 0.010bq: ASI7,3J:: 0.02b02: ASl7,qJ:: 0.01217: lS,5):: ASI7,bJ:: AS[7,7J:: -0.18blq: lS[7,6J:: O.b08b8: ASI6,l):: O.OOOOq: AS[6,2J:: ASI6,3J:: 0.00811: AS[&,qJ:: O,Qq10b: AS{8,5]:: O.ooq67: lS{6,b]:: 0.07255: AS[8,7):: o,qos71: AS[8,8J:: O.Slbb2: 8SI1J:: 3.bb013qE-S: BSl2J:: BSI3J:: 8.ll32b1E: 8S(q):: 2.737l3qE-2: BS(SJ:: 2.370377: BS(b]:: BSI7J:: 1.lOqlS5E-q: BSI8J:: 7,5q12E-q: 157

PAGE 167

fOR ( :: l TO lb DO fOR J :: l TO lb DO P[(,JJ :: O.Ol; fOR [ :: l TO lb DO P(t,!J :: l: fOR [ :: l TO 6 DO BEGIN u ([J :: 0 ; Y([J :: 0 : X[[J :: 0 : Sigma[[) :: Y(,[J; Slgma([t&l :: U(,[J; theta(() :: l([J: Theta([t6] :: 8([): END: t :: O: Uin := O: frror :: O: \ VRITH'ENTER THE DESIRED TIME Of THE SIMUUTION: '}: RUDLN{tf}: ASSIGN{DATl, 'ADAPTEST.LSQ'}: iEVRITEUATA}: REPEAT {1 Begin Simulation Routine 1} {I Simulate the System using State ftatrlx Equations 1} If t : 0 THEN fd :: ELSE fd :: O: fOR J :: l TO 8 DO BEGIN XHEV[J] :: O: FOR [ :: 8 DO XNEK(JI :: XNEV[JI t AS(J,[)IXi[J: XNEV(J) :: XNEV[J] t BS[J)Ifd t BS(J]IUin: END: yd := xm: If t I 0 THEN BEGIN {1 Update Theta's Parameters-the estimator's 1} {1 Parmeters Using the ftethod of least Squares 1} fOR J :: l TO lb DO {1 Generate denominator : l Sigma' x P x Sigma 1} BEGIN SigP(JJ :: O: {1 SigP : Sigma' K p 1} fOR [ :: l TO lb DO SlgP[JJ :: SlgP[J) t Slgma([)IP[(,JJ: END: Denom :: l: fOR [ :: l TO DO Denom :: Denom t SlgP([)ISigma((J: 158

PAGE 168

fOR J :: l TO lb DO BEGIN Psig(JJ :: 0: {1 Psig : P{t -d} x Signa 1} fOR t :: l TO lb DO Psig(JJ :: Psig[JJ 1 PlJ,[)ISigma[tJ: END: Ydsys :: 0 : {1 Generate Error : y -Sigma' x Theta 1} fort:: l TO lb DO Ydsys := Ydsys 1 Sigma[[)IThetaltl: Error :: yd Ydsys : fORt :: l TO lb DO {1 Update Theta 1} Theta([) := Theta[[) 1 Psig([)IError/Denom: fOR J :: l TO lb DO {I Update P{t -d} 1} FORt :: 1 TO lb DO P[J,[J :: P(J,tl Psig[J)tSigP([J/Denom: END: {1 System DARMA Step Ahead Estimator Routine 1} Ydplusl :: UiniBUJ yd1Alll: fOR [ :: 2 TO 6 DO Ydplusl :: Ydplusl I 8[[)1U(l0-[) A([)IY(lO-tJ: fOR [ := l TO 1 DO {1 Reset the values for 1} BEGIN Uf[J := U(tlll: Yltl := Y[tll): END: Y(6) :: yd: U[6) :: fd 1 Uln: fOR [ :: 1 TO 6 DO {I Reset the values for X{k}, A{k}, 1} BEGIN {I B{k}, and 1} Sigmaltl :: Sigma([16) :: Xltl :: XNEVltl: Altl ::Theta([); Bltl := Theta([16): END: VRITElN{DATA,t:q:2,' ,yd,. ', Ydplusl,' ',Error}: t : t I td; UNTil {tf I epsilon < t}: ClOSE{DUA}: END. 159

PAGE 169

PROGRAM DARftA_VITH_fEEDBAU: CONST td : 0.05 : epsilon: 0.0001: TYPE vector : array[l 8) of real: VAR [ : integer: t, tf, yd, Ydesired, Ydplusl, fd, Uin : real: A,8,U,Y : vector: DATA: BEGIN {r INITIALIZE PARAMETERS r} All) :: 0.17558,1 :: Al31 :: .3883131 Am := A!SJ :: AILJ :: Alll :: 0.17558,3 A(8] :: 1 fOR [ :: 1 TO 8 DO BEGIN um := o : vm := o : END: t :: D: Uin :: 0: : 8{1] :: : {1 : :: : {r : 8[3] :: 1.35,,8lE-2 : {1 : 8(q) :: -b.b,08,,E-3 : {r : 8[51 :: L.L,08,,E-3 : : 8(LJ :: -1.35,,81E-2 : : 8[7) :: -1.24l56SE-2 : : 8(8) :: DARftl ftODEl r} AND THE ESTiftlTOR r} ftODEl ARE IDENTICAl r} fOR THE SISO CASE r} WRITE{'ENTER THE DESIRED TiftE Of THE SiftULATION:'}: READLN{tf}: 'ARftAfDBt'}: REVRITE{DATAl: 160

PAGE 170

REPEAT {1 Begin Simulation Routine 1} yd := 0 : {1 Simulate the System using DARftA 1} If t : 0 THEN fd ;: -10 ELSE fd :: O: fOR ( :: l TO 8 DO yd :: yd t 8[() A[(J {1 System Estimator Routine 1} Ydplusl := UintB(lJ ydlt[ll: fOR ( :: 2 TO 8 DO Ydplusl :: Ydplusl t 8([)tU(l0[) t(()1Y(10-(J: {1 .Control Law Routine 1} Ydesired :: Ydplusl10,15: Uin :: Ydesired t ydiA(lJ: fOR [ :: 2 TO 8 DO Uln :: Uln {8[[JtU[l0[) -A([)IY(lQ[J}: Uln :: Uin/B[lJ: fOR [ :: 1 TO 1 DO {1 Reset the values for and 1} BEGIN U[[) :: U[[tll: Y[[J :: Y([tll: END: Y{8) :: yd: U!81 :: fd t Uin: VRITELNUATA, t: 2,1 I ,yd, I 1 Ydplusl, 1 ',U in}: t :: t t td: UNTIL {tf t epsilon c t}: CLOSHDATA}: END. 161

PAGE 171

PKOGKAM CONST td : 0.05 : 0.0001 : TYPE vector : array[l 6) af real: matrix: array[l 6,1 6) of real: VAK J, 1: : integer: t, tf, yd, Ydesired, dplus1, fd, Uin, Uin_old : real: AS : matrix: A, 8, U, Y, BS, X, Xnew: vector: DATA: text: BEGIN Alll :: 0.17556,1 At21 := Alll :: -0.3663131 :: Al51 :: -0.3663132 Albl :: Al7J :: 0.1755613 A(6) :: 1 {1 INITUUZE PARAftETERS 1} : 8(1} :: 2.7]7l]qE : {I DAKMA MODEL 1} : 8121 :: : {1 : 8[31 :: 1.35,161E-2 : {1 : := -b.b10611E-3 : {1 : Bl5J :: b.b10&11E-3 : : BlbJ := -1.351161E-2 : : am : = -2 : : 8[6) :: : AND THE ESTIMATOR 1} nODEL ARE IDENTICAL 1} fOR THE SISO CASE 1} ASU,1J:: -0.25666: ASU,cl:= O,qbb02: ASl1,3J:: 0.01101: -0.0073: AS[1,5J:: -21.66b7,: AS[1,bl:= 3.65011: ASl1,71:= b.30602: !S[l,6):: 0.60676: ASl2,1J:: o.qbb02: AStc,cJ:: -0.1516: AS(2,3J:: o,q70&1: AS[2,qJ:: 0.10662: ASl2,5J:: 3.65011: AS[2,bl:: -23.57677: AS(c,71:= AS{c,61:= 7.36637: AS(3,1J:: 0.09,01: AS[3,2J:: AS[3,3J:: AS(3,q):: 0.60666: AS[3,5J:: 6.30602: AS[l,bJ:: AS(3,7J:: lSl3,6J:: 0.07255: AS(q,3):: 0.51662: AS(q,5):: 0.53116: AS(q,b):: 6.'152671: AS!5,11:: 0.02523: AS[5,2):: O.OlOb: AS[5,3J:: 0.00111: AS{S,q):: 0.00005: ASI5,5):: -0.25666: AS[5,b):: o.qbb02: AS[5,7J:: 0.0110'1: AS[5,6):: 0.0073: AS[b,1):: 0.0106: AS[b,2):: 0.02b3q: AS{b,3J:: 0.010bq: AS[b,q):: 0.00111: ASlb,51:= ASib,bJ:: -0.15,6: ASib,7):: AS[b,6):: 0.10662: AS[7,1J:: 0.00111: AS[7,2):: O.OlObq: AS[7,3):: 0.02602: 0.01217: AS[7,5J:: 0.01'10,: AS!7,bJ:: AS(7,7):: -O.l8blq: AS(7,8):: O.b08b6: AS[6,l):: AS[6,2J:: 0.00079: 0.006.1: AS[8,q):: AS[6,5):: AS(6,b):: 0.07255: AS(6,7]:: AS(8,6]:: O.Slbb2: BSl1J:: BS(2J:: 7.121027-q: BSl3J:: 6.1132b1E-3: BS(5J:: 2.370377E-7: BS(b):: b.9b2379E-b: BS!7J:: BS!8):: FOR [ :: 1 TO 6 DO BEGIN U[[) :: 0 : Y([J :: 0 : Xl[J:=O: END: 162

PAGE 172

t :: D: Uin :: O: Uln_old :: O: URITE{'ENTER THE DESIRED TinE Of THE SIMULATION: '}: READLN{tf}: ASSIGN {DATA,' ARftAfBL Sift'}: REWRITE {D lT A}: REPEAT {1 Begin Simulation Routine 1} {1 Simulate the System using State ftatrix Equations 1} If t : 0 THEN fd :: -10 ELSE fd :: 0: fOR J :: l TO DO BEGIN XNEWIJJ :: O: fOR ( :: 1 TO DO XNEW(J) :: XNEV(JJ t XNEW(JJ :: XNEW(JJ t BS(JJfd t BS{J) Uin: END: yd :: {1 System DARftA step ahead Estimator Routine 1} Ydplusl ;: UinBill -ydrA[lJ: fORt :: 2 DO Ydplusl :: Ydplusl + BltJUilD-(J -AltlYilO-tl: {1 One-Step-Ahead feedback {ontrol law Routine 1} Ydesired :: YdpluslrO.a: Uin :: Ydesired t ydiA(lJ: fOR C :: 2 TO a DO Uln :: Uln {B(CJrU(lO-tJ A(tJIY(lD-tJ}: Uln :: Uin/8(11: fOR [ :: 1 TO 1 DO {1 Reset the values for U{l:} and W:} 1} BEGIN U([J :: U1Kt1J: Y{[) :: Y([ t l): END: YlaJ := yd: :: fd t Uin_old: Uin_old :: Uin: fOR [ ;: l TO DO {1 Reset the values for X{l:} 1} Xltl :: XNEVI'tJ: ',yd,' ',Ydplusl,' ',Uin}: t :: t f td: UNTIL {tf t epsilon< t}; CLOSHDAT Al: END. 163

PAGE 173

PROGRAn Compensator_Pole_Piacement_lst_nODE_FEEDBAC(: CONST td : 0.05 : epsilon: 0.0001: gain : q,5 : TYPE vector = array(1 8) of real: matriK : array(1 8,1 . 6) of real: VAR J, (, N : integer: t, tf, yd, ydfilt, fd, Uin: real: AS : matr iK: BS, X, Xnew : vector: Xf, Yf, Lz, Pz : array[O 41 of real: DATA : teKt: BEGlN {1 INITIALIZE PARAMETERS 1} AS(l,ll:= -0.256b6: ASll,2J:: 0.4bb02: lS[l,3):: ASI1,4J:: -0.0073: ASll,5J:: -29.66b79: lS(1,LJ:: 3.65011: ASI1,?J:: b.30602: lS(1,6J:: 0.6067b: lS(.2,lJ:: 0.4bb02: lSl2,2J:: -D.l5'lb: AS(2,3J:: 0.47069: lS[2,q):: 0.10662: ASI2,5J:: 3.65011: AS12,bJ:: -23.57617: AS(2,7]:: q,WI29: AS(2,6J:: 7.36b37: AS(3,lJ:: 0.09909: AS(3,2J:: O.q7069: AS(3,3J:: -0.16blq: AS[3,4J:: O.b06be: AS(3,5J:: b.30602: lS(3,b):: q,36929: AS(3,7J:: -2q,9b25q: AS(3,6):: lS(q,1):: 0.00467: I.S{4,2):: 0.07255: AS(q,3J:: o,qo57'1: AS.(q,q):: 0.51bb2: AS(q,5J:= 0.53916: lS(q,b):: q,92q25: AS(4,7J:: AS(q,a):: -14,q43Sb: lS[5,1J:: 0.02523: ASIS,2J:: O.DlOb: l:S(S,3):: 0.00111: AS(5,q):: 0.00005: lS(5,5J:: -0.256b6: ASlS,bJ:: O,qbb02: AS15,7):: AS(5,6):: 0.0073: AS(b,1J:: 0.010b: AS(b,2J:: 0.02b34: AS(b,3):: O.Ol0b4: AS(b,q);: AS{b,5):: O,qbbD2: AS(b,b):: AS(b,7):: AS(b,6):: 0.10662: AS(1,1J:: 0.00111: AS(1,.2):: 0.010bq: AS(7,3):: 0.02b02: AS(7,q):: 0.01217: AS(7,5):: ASI7,b):: O.q706'1: AS(7,7):: -0.16bl4: AS(7,6J:: O.b08b6: AS(6,1):: 0.00004: AS{6,2):: AS(8,3):: AS(6,4):: 0.04lOb: lSI6,SJ:: 0.00467: lSl6,LJ:: 0.07255: AS(.6,7):: 0.40579: lS(6,6J:: 0.5lbb2: BSllJ:: 3.bb013qE-5: BSI2J:: 7.'12l027E-4: BS(3J:: BS(q):: 2.73713qE-2: BSISJ:: 2.370377E-7: BSfbJ:: 8Sl7J:: BS[6J:: lz(ll :: -0.1'1qb001 : Lil21 :: 1.271007 : Lz131 :: Lz(q) := 0.'1'1'17'1'1 : PzlOI := l : PzBJ := 0.32814 : PzW :: l,q8'l1L3 : Pzl31 :: 0.326057b: Pz(q) := : fOR ( := 1 TO 6 DO Xl(J :: O: t :: O: Uin :: 0: N :: O: 164

PAGE 174

WRITE{'ENTER THE DESIRED TIM Of THE SIMULATION: '}: READLN{tf}: ASSIGN{DATA, REWRITE{DATAJ: REPEAT {1 Begin Simulation Routine 1} {1 Simulate the System using State ftatrix Equations 1} If t : 0 THEN fd :: -10 ELSE fd :: O: fOR J :: l TO 6 DO BEGIN XNEV(JJ :: 0: fOR [ :: l TO 8 DO XNEVfJJ :: XNEV(JJ t AS(J,[JIX([J: XNEVlJJ :: KNEVlJJ t BSlJJ1fd t BS{J)IUin: END: yd := xm: {1 Simulate the Digital Compensator 1} If N t q THEN N :: Xf(NJ :: -yd: Ydfi. lt :: O: fORt :: 0 TO H DO Ydfilt :: Ydfllt t Pz((]IXf(N-tl: fORt :: l TO N DO Ydfilt :: Y$ffllt -lz([)IYflN-[J: Uln :: galnYdfllt: fOR [ :: l TO 6 DO {I Reset the values for X{l:.} 1} Xltl :: XNEW(tl: Yf(NJ :: Ydfilt: If N : q THEN {1 Reset the filter values t} fOR [ :: l TO DO BEGIN Xf([-lJ :: Xf(tJ: Yf(t-11 :: Yf(tJ: END: WRITELN{DATA,t:q:2,. ',yd,' ',Uin}: t::tttd: INC{N}: UNTil {tf t epsilon ( t}: ClOSE {DATA}: END. 165

PAGE 175

991 :o :: N :o =: u1n :o :: 1 :o =: mx oo 'l 01 t =: l : :: : bfb202b'O:: IE)Zd : b"I2Db"IE'O :: (2)Zd: bb5202b'O:: (UZJ: f :: IO)Zd : "l'lblbbb'O :: : :: ([)Zl : L00tl2't :: (2)Zl : :: ({)Zl :!(V)SB =!ILI5B :"l-3bLE29b'"l ::("1)59 :l-3LLEOLE'2 :!(5)59 !E-3b"I2Ett''l ::( =!12159 ::({)59 :2"1"1t5'D :!(l'9lSY :ss2LD'D =!I"J'9JSY =!15'9lSY =!IE'IlJSY :bLOOO'O =!I2'9JSY ::(t'91SY =!l!I'LlSY :!(L'llSY =!l"l'llSY :bObbO'O =:ts'LlSY :l{2t0'0 :20"120'0 =!IE'LlSY =!l2'llSY :tttOD'D ::(t'llSY :!(!I'"IJSY =!ll'"IJSY :"lb5t'O-:!["l'"llSY =![5'"11SY :bttOO'O =!I2'"11SY :"IOtD'O ::(t'"IISY :ELOO'O =![9'5JSY :bDbbO'O =![L'51SY :!["J'5JSY :9"1952'0-=!(5'5JSY :sooooo :tttooo ==tE'sJsY :"Jotoo ==t2'51SY :E252oo ==tt'SJSY :bL"I25b''l :!ltb'0 :s52LO'O =!I!I'EJSY ::tl'BSY :![9'EJSY :20110E'"I ::ts'ElSY :qcJ90-'i'O :!(E'EJSY :!f2'ElSY :bObbO'O ::(t'EJSY :H99E'L ::(9'2lSY =!IL'21SY :u!IL5'E2-=!l9'2JSY :no59'E =![5'2l5Y :299ot'O :: :: IE'2l5Y :"lbsto: !(2'21SY :: lt'21SY :9L!IO!I'D :![9'tJSY :2090E''i :!(L'tlSY :tt059'E =!l"l'tlSY :bl"l99'b2-:!(5'tl5Y :ElDO'O-:bObbD'O ::(E'tlSY :![2't1SY :99952'0-::(t't)SY {1 3ZI1YI1INI 1} ::Jxa:J : YlYG : Jl!aJ JO 'Q)Al!JJI! : Zd 'Zl 'JA '.:IX :JO:JOaA : MaUX 'X '58 :x1 J:Jew : SY :1eaJ: u!n ':JI!lPA pA 'H ':J :Ja6alU! : H '1 r :1eaJ JO (!l''t'll''t)AeJJe : X!Jll!W :1eaJ 10 [9''t)AeJJe: J010aA JdAl : 9 : u1e6 : toooo = uoJ!sda : so o = Pl lSNOJ

PAGE 176

WRITE{'ENTER THE DESIRED TIME Of THE SIMULATION: '}: READLN{tf}: ASSIGN{DATA, 'COMPf8(2.A6'}: REWRITHDATH: REPEAT {1 Begin Simulation Routine 1} {1 Simulate the System using State Matrix Equations 1} If t : 0 THEN fd := ELSE fd :: O: fOR J :: l TO 5 DO BEGIN XNEW(J] :: O: fOR ( :: l TO 5 DO XNEW(JJ :: XNEW(JJ t AS(J,()IXftl: XNEK[JJ :: XNEW(JJ t BS(J)Ifd t BS[JJIUin: END: yd := xm: {1 Simulate the Digital Compensator 1} If H THEN H :: XflNI := yd; Ydfilt := 0; fOR ( :: 0 TO H DO Ydfllt :: Ydfllt t Pz(()IXf[NtJ: fOR [ := l TO H DO Ydfilt :: ldfilt tzt(JIYf(N-tJ: Uln :: -galniYdfilt: fOR ( :: l TO 6 DO {1 Reset the values for 1} XltJ :: XNEV(!J: Yf!Nl :: Ydfilt; If N: {1 Reset the filter values 1} fORt :: l TO DO BE, IN Xf[t-11 :: Xf{tJ: Yflt-11 :: Yf(tJ: END: ',yd,' ',Uin}: t :: t td: INC {N}: UNTIL {tf t epsilon l t}: CLOS[{DATH: END. 167

PAGE 177

PROGRAn Compensator_Pole_Piacement_lst_and_2nd_ftODE_fEEDBAC[_alternate_method: CONST td : 0.05 : epsilon: 0.0001 : gain : : gain2 : 8 : TYPE vector : 8) of real: matrix : arrayU 811 11) of real: VAR J1 [I N : integer: t1 tf, yd, ydfilt1 fd, Uin1 Ul, U2 : real: AS : matrix: BS1 X, Xnew Xf 1 Yf 1 Yf2, lz, Pz, lz21 Pz2 DATA : vector: : array(O of real: : text: BEGIN {1 INITIALIZE PARAftETERS 1} AS(l,1J:: -0.2581.8: AS[1,3):: 0.0'1'10'1: -0.0073: ASU,51:= -2'1.881.7'1: 3.85011: AS(l,?J:: ld0802: ASU,8J:: 0.8087b: AS(2,1J:: ASI2.2J:: -0.15'11.: AS(2,3J:: 0.10882: 3.85011: AS(2,bJ:: ;5?877: AS12,7J:: AS[2,8J:: 7.381.37: AS[31 1):: 0.0'1'10'1: AS(3,2J:: AS(l,3J:: O.b08b6: AS[3,5J:: 1..30802: AS[3,b):: AS[3,7J:: AS.[3,8J:: 0.07255: 0.51bb2: 0.53'118: 8.'1521.7'1: AS[S,l]:: 0.02523: AS,2):: 0.0101.: AS,3):: O.OOlU: 0.00005: AS[S,SJ:: -0.2581.8: AS,1.):: o,qbi.D2: AS!S,7):: AS[5,8J:: 0.0073: D. OlD!.: AS(b,31:: AS O.DDU'I: AS(b,5):: o.qbi.02: AS(b,b):: -D.lS'Ib: AS(b,7):: 0.10882: .AS(7,lJ:: 0.00111: AS(712):: 0.010bq: 0.021.02: 0.01217: ASI7,5}:: 0 .0'1'10'1: ASI711.]:: O.b08b8: AS(8,1J:: AS(81 2):: 0.0007,: AS(8,3):: 0.00811: O,OqlOb: AS[81 5):: AS(8,bJ:: 0.07255: AS,81:= 0.51bb2: BS(lJ:: BS):: BS(3J:: 8.ll32b'IE-3: 2.737l3qE-2: BS(5):: 2.370377E-7: BS(b):: b.'lb237'1E-b: BS(7J:: 1.lOqlSSE-q: BS[8):: 7,Sq'I2E-q: lz[11 :: -0.1'1qb001 : lzl21 :: 1.271007 : lz{31 :: -0.1'1qqb28: lz(q) :: 0.'1'1'17'1'1 : Pz(OJ ::. 1 : Pzlll :: : Pz(21 :: l.q8'llbl : Pz(3J :: 0.]280571. : Pz(q) :: 0.'1'1'1'1 : lz2111 :: -O.l'lqb001: lz2121 :: 1.271007 : lz2131 := :: O.'l'l'l7'18b : Pz2101 :: 1 : Pz2UI :: -0.'12025'1'1 : P-z2l21 :: 0.3b'I02b'l: Pz2(3) :: -0.'1202'11'1 : :: : fOR ( :: l TO 8 DO X(() :: 0: 168

PAGE 178

t :: 0: Uin :: 0: N :: 0: WRITH'ENTER THE DESIRED TinE Of THE SinUUTION: 'J: READLN{tf}: ASSIGN{DATA, 'CnPfB(l2.A'}: REWRITE{DATA}: REPEAT {1 Begin Simulation Routine 1} {1 Simulate the System using State ftatrix Equations 1} If t: 0 THEN fd ;: -10 ELSEfd ;: O: FOR J :: 1 TOft DO BEGIN XNEUJJ :: O: FORt ;: 1 TO & DO XNEW[JJ ;: XNEW(JJ AS[J,[)IX{(J: XNEW[Jl :: XNEW(JJ BS(JJfd BS{J)IUin: END: yd := xm: {1 Simulate the Digital Compensators 1} If N t THEN N :: Xf[NJ :: -yd: Ydfllt ;: O: fORt:: 0 TON DO Ydfilt := Ydfilt I Pz(()IXf[N-tl: fOR ( :: l TO N DO Ydfllt :: Ydfllt Lz{[)IYf[H-t): Ol :: gainiYdfllt: Ydfllt2 :: O: fOR ( :: 0 TON DO Ydfllt2 :: Ydfilt2 Pz2(()1Xf(N-tl: FOR ( :: 1 TO N DO Ydfllt2 :: Ydfllt2-lz2[()1Yf2[H-(J: 02 :: -galn21Ydfilt2: Uln :: Ul 02: FOR ( :: 1 TO.ft DO {1 Reset the values for X{t} I} X((J ;: XNEW[tJ:. Yf[NJ ;: Ydfllt: Yf2[NI :: Ydf i I t2: If N {t Reset the filter values 1} FOR [ := 110 DO BEGIN Xf([ll :: Xf(t): Yf([-11 :: Yfltl: Yf2(t-11 :: Yf2[tl: END: ',yd,' ',U1,' ',U2,' ',Uin}: t ;: t I td; INC{NJ: UNTIL {tf t epsilon< t}: CLOSHDATAJ: END. 169

PAGE 179

PROGRAn ADAPT_EST_least_Squares_with_Noise: CONST td : 0.05 : epsilon = 0.0001 : TYPE vector : array(1 6) of real: vector( : arrayU .. 11.) of real: matriK : array(1 6,1 6) of real: matrhcc : arrayU .. 11.,1 .. 11.1 of real: VAR J, [ : integer: t, tf,.yd, Ydplus1, fd, Uin, Error, Denom, ldsys, ldnoise : real: : matrix: P : matrlx2: A, 8, U, l . BS, X, Xnew SlgP, Pslg, Theta, Sigma DATA BEGIN {1 INITIALIZE PARAftETERS 1} : vector: : vector2: : text: AUI :: O.l U21 := am := -o.3 :: 0.3 A(5) :: -0.3 Ubi :: A :: 0.1 U&l :: .L : BUI :: 2;7-2:' {t URnA ftODEL 1} : 8[2) :: 1.2E-c : {1 which is the transfer 1} : 8(3] :: 1.3E-2 : {1 function of the system 1} : :: -b.I.E-3 : : 8(51 :: l..bE-3 : : 8(1.1 :: -1.3-2 : : B := -1.2-2 : : 8[6) :: -2.7-2: ASI.L,1J:: -0.2561.6: AS(1,2J:: AS(1,3):: -0.0073: ASI.L,51: = AS11,bl:: 3.650U; ASU, 71:: 1..30602: ASU,&l:: 0.60871.; AS{2,1J:= AS(2,2J:: AS(c,3J:: 0.10662; AS(2,SJ: = 3.6SO.ll; ASIC,bl:: -23.51617: ASI2, 11:: AS(2,&1:: 7.361.31: AS,11:= AS,31:= 0.1.061.6: AS(3,5J:: 1..30602; AS(3,1.J:: AS(3,1J:: AS(3,6]:: 0.07255: 0.511.1.2: AS(q,I.J:: AS(q,7J:: AS(q,6):: -1q,qq]SI.: AS,1):: 0.02523: AS(5,2J:: 0.0101.: AS(5,3):: 0.00111: AS(5,q):: 0.00005: AS(5,5J:: -0.2561.6: AS(5,1.J:: o.qbi.02: AS(5,7):: AS(5,6):: O.OOH: 0.0101.; ASib,2J:: AS.,31:: O.OlOI.q; AS(I.,q):: AS(b,5J:: AS(L,bJ:: -0.15'11.: AS(1.,7J:: AS(b,6):: 0.10662; AS(7,lJ:: 0.00111; AS(1,2J:: lS(7,3J:: 0.021.02: 0.01217: AS17,5):: 0.0'1'10'1; AS(7,1.J:: AS[7,7J:: -O .L&blq: AS,6):: O.b06b6; AS[6,1):: O.ooooq: AS(II,c):: AS(l,3J:: 0.00611: AS[II,q):: O,Oq10b: AS(6,5):: O.ooq67: lS(6,1.):: 0.07255: AS[6,7):: O,q057'1: AS[II,6):: BSllJ:: BS[cJ:: BS[JJ:: 2.7371JqE-2: BS:: 2.370377-7: BS[b]:: BS[7J:: 1 .LOql55E-q: BSIIIJ:: 7.5q'I2E-q: 170

PAGE 180

fOR [ ::. l TO lb DO FOR J ;: l TO lb DO PI[,JJ ;: 0; FOR C :: l TO lb DO P[t:,C) :: 1: FOR [ :: l DO BEGIN U([J :: 0 : Y[() :: 0 : X[t:J :: 0 : Sigma[tl :: -Y[,-tl: Slgma[(t6] :: U(,-[J: Theta1tJ :: Altl: Theta([t6J :: Bltl: END: t :: O: U in :: O: Error :: O: WRITH'ENTER THE DESIRED TinE Of THE SiftULATION: '}: READLN{tf}: ASSIGN{DATA. 'ADAPEST.STC'l: REWRITHDATA}: REPEAT {a Begin Simulation Routine 1} {1 Simulate the System using State ftatrix Equations 1} If t : 0 THEN fd :: -10 ELSE fd :: O: fOR J :: l TO 8 DO BEGIN XNEVIJI :: D: fORt :: l TO 8 DO XHEV(J) :: XNEV(J) t AS[J,[JXICI: XHU[J) :: XNEV[JJ t BS[JJFd t BS(JJtUin: END: yd :: X(qJ: {I Adding sensor noise to the system I} Randomize: Ydnoise :: yd t O.ll{random-O.S}I{random-O.S}I{random-0.5}: If t I 0 THEN BE IN {1 Update Theta's Parametersthe estimator's 1} {I Parmetets Using the ftethod of least Squares t} fOR J :: 1 TO lb DO {1 Generate denominator : l t Sigma' x P x Sigma } BEGIN SigP(J) :: O: {1 SigP : Sigma' X p I} fORt :: 1 TO lb DO SigP(J) ;: SigP(J) t Sigma(t:JP(t,JJ: END: 171

PAGE 181

Denom :: 1: fOR ( :: 1 TO 1b DO Denom :: Denom t SigP(()ISigma(CJ: fOR J :: 1 TO 11. DO {I Update P{t -d} 1} fOR C :: 1 TO 1b DO P(J,[J :: P(J,CI Psig(J)ISigP((]/Denom: Ydsys :: 0 : { 1 Generate Error : y S i 11"8' x Theta 1} for C :: 1 TO 1b DO Ydsys :: Ydsys t Sigma((JtTheta(tJ: Error :: Ydnoise Ydsys : fOR J :: l TO lb DO BEGIN Psig(J) :: O: {t Psig : P{t} x Signa 1} fOR C DO Psig(JJ :: Psig(JJ t P(J,(J1Sigma1CJ: END: fOR [ :: 1 TO ib DO {I Update Theta 1} Theta([) :: Theta[.() t Psig[()IError: END: {1 System DARftA Step Ahead Estimator Routine 1} Ydptusl :: UinBlll [ :: 2 TOe DO Ydplusl :: Ydplusl t 8[[)1U[l0[) A{()IY[lO-tl: fOR [ 1 DO {1 Reset the values for U{k} and Y{k} 1} BEGIN U[CJ :: U(CtlJ: Y[CJ :: Y(CtlJ: END: Y[e) :: Ydnolse: U(eJ := fd urn: fOR [ :: e DO {t Reset the values for X{k}, A{k}, 1} BEGIN {I 8{1:}, and Signa{i:} 1} Sigma[() :: Yl'Hl: SigmalteJ :: X!CI :: XNEV!CI: A([J :: Theta[CJ: Bltl := Theta!teJ: END: ',yd,' ',Ydnoise,' ',Ydplusl,' ',Error}: t :: t I td: UNTIL {tf 1 epsilon< t}: CLOSHDA TAl: END. 172

PAGE 182

PROGRAn AUPT_EST_leastJquares_With_!alman_filtering: CONST td : 0.05 : epsilon : 0.0001 : Rk : 0.05: {1 Rk is the ll,ll covariance matrix 1} TYPE vector : array[1 eJ of real: vector2 : arrayU .. 1bl of real: matrix : array[l .. e,1 .. of real: matrlx2 : arrayU .. lb,1 .. lb] of real: YAR I, J, ( : integer: t, tf, yd, fd, Uin, Error, Denom, Ydsys, Ydnoise : real: AS, PI:, AI:, API:, tCk : P : matrlx2: A, B, U, Y, BS, X, Xnew, Bk, (1:, Ck, XI:, Xknew : vector: SigP, Psig, Theta, Sigma : vector2: DATA : text: BEGIN All) :: 0.1 At21 := -l[3] :: -0.3 :: 0.3 A[51 :: -0.3 Ubi :: A(?) :: 0.1 :: l {1 PARAftETfRS 1} : Blll :: 2.7-2: : Bl2l :: l.2E-2 : : 8[3) :: l.3E-2 : : B(q) :: -b.bE-3 : : 8[51 :: b.bE-3 : : Blbl ;: -l.3E-2 : : 8171 :: -1.2E-2 : : 8{61 :: -2.7E-2: {1 DARftA ftODEl 1} {1 which Is the transfer 1} {1 function of the system 1} AS[1,l]:: -0.256b6: AS[1,2J:: Q,qbb02: AS(1,3J:: 0.09909: -0.0073: AS[1,5):: -29.66b79: AS[1,b]:: 3.65011: ASl1,7]:: b.30602: AS[1,6J:: 0.6067b: AS[2,l):: o,qbb02: AS[2,2]:: -O.l59b: AS[2,3]:: o,q7069: AS[2,5):: 3.65011: AS[2,6]:: -23.51"677: AS12,1J:: q,36929: AS[2,6]:: 7.3eb37: AS[3,1J:: 0.09909: AS[3,2]:: O,q7069: AS[3,3J:: -O.l6blq: AS(3,q);: O.b06b6: AS[3,SJ:: b.30602: AS[3,b):: q,36929: AS[3,7):: -2q,9b254: AS[3,6]:: l3,q2902: ASI4,1J:: O.ooqa7: AS(q,2):: 0.07255: AS(q,3):: 0.40579: AS[4,4):: 0.51662: ASI4,5):: 0.53916: AS[4,b):: 4.92q25: AS[4,7]:: AS(q,aJ:: -1q,q435b: AS[5,l]:: D.02523: AS[5,2J:: O.OlOb: AS[5,3]:: 0.00111: AS[S,q);: 0.00005: AS(5,5);: -0.25666: AS[5,6):: o,q6602: ASl5,7J:: 0.09909: AS(S,6J:: 0.0073: AS[b,lJ:: O.OlOb: ASlb,2J:: 0.02634: ASlb,3J:: O.OlObq: AS[b,q);: 0.00119: AS[b,5);: o,q6b02: AS(b,b):: -O.l59b: AS[b,7J:: Q,q7069: AS[L,6):: 0.10662: AS[7,1J:: 0.00111: AS(7,2J:: O.OlObq: AS[7,3);: 0.02602: ASf7,q);: 0.01217: AS(7,5J:: 0.09909: AS{7,bJ:: 0.47069: ASl7,7J:: -O.l6bl4: AS[7,6J:: O.b06b6: AS16,1J:= O.OOOOq: ASl6,2):: 0.00079: AS[6,3J:: 0.00611: AS[6,4):: O,Oq10b: AS[6,5J:: 0.00467: AS[6,bl:= 0.07255: ASl6,7J:: 0.40579: AS(6,6J:: O.Slbb2: 8SilJ:: 3.bb0134E-S: BSl2J:: 7.921027-4: BS(3J:: 6.ll32b9E-3: BSl4J:: 2.73713qE-2: BSl5J:: 2.370377-7: BS(1):: b.9b2379E-b: BSI7J:: 1.10ql55E-4: BSl6J:: 7.5492[-q; 173

PAGE 183

fOR ( :: 1 TO 1b DO fOR J :: 1 TO 1b DO P(I:,JJ :: O: fOR [ :: 1 TO lb tO PU:,:J :: 1: fOR [ :: 1 TO ft. DO {1 Initialize Parameters for theta and sigma 1} BEGIN U([J :: 0 : Y([J :: 0 : X([J :: 0 : Si!JIIilUl :: -YI,[J: Sigma([tftl :: Ul,CJ: Theta([) :: l((J: Theta([t5) :: Bl[J: fND: t :: O: Uin :: O: Error :: O: fOR I :: 1 TO a DO fOR J :: l TO 6 DO BEGIN AUI,JJ :: O: PI:(I,JJ :: D: END: fOR I :: 2 TO a DO AI:U,I-11 :: l: fOR I :: 1 TO a DO BEGIN CUll :: O: :: O: XI:: (I) :: O: Pl::[l,rJ :: 0.1: END: :: l: {1 Initialize [alman Observable State natrices, 1} {1, [alman Gain ftatrix Error Covariance 1} {I natr lx PI:, and hlman State Vector XI:. 1} WRIT({'ENTER THE DESIRED Tift Of THE SinULATION: '}; RUDLN{tf}; ASSIGN{DATA, 'UlftNEST. STC I}: REWRITHDAT U: 174

PAGE 184

REPEAT {1 Begin Simulation Routine 1} {1 Simulate the System using State natrix Equations I} If t : 0 THEN fd := -10 ElSE fd :: O: FOR J :: l DO BEGIN XNEW(JJ : = O: FOR ( :: l TO DO XNEV(J) := XNEW(J) I XNEW(J) :: XNEVIJI 1 BS(J)Ifd 1 BS{J)tUin: END: yd :: xm: {1 Adding sensor noise to the system 1} Randcrnize: Ydnoise :: yd O.l{randcrn-O.S}t{random-O.S}t{randomO.S}: { (AlftAN STATE ESTiftATOR ALGORITHft 1} {1 Update the Observable [alman State ftatrices 1} FOR I:: l DO BEGIN :: :: END: {1 STEP l : State Estimates Extrapolation 1} FOR J :: l TO 6 DO BHIN Xl:new[J) :: O: FOR [ := l TO 6 DO Xl:new[J) := Xl:new(J) 1 Xl:new(JJ :: Xl:new(J) 1 Uin I Bk[JJ lfd: END: fOR I := l DO Xk[ll :: Xl:new[IJ: {1 STEP 2 : Error Covariance Estimate Extrapolation 1} fOR I :: 1 TO 6 DO fOR J :: l TO 6 DO BEGIN API:U,JI :: O: fOR ( :: 1 DO APk(I,JI :: API:[I,J) 1 AI:(I,[)IPk[t,JI: END: 175

PAGE 185

FOR I :: 1 DO FOR J :: 1 TO & DO BEGIN Pk(I,JI :: D: fORK :: 1 TO & DO Pk(I,JI :: Pkli,Jl t APk(l,[)IAk([,JJ: END: { STEP 3 : New State Estimate 1} FOR I :: 1 TO & DO Xklll :: Xklii t Kkll)I{Ydnolse Xk[6)}: {1 STEP : New Error Covariance Estimate 1} fOR I :: 1 TO a DO fOR J :: 1 TO a DO KCtli,JJ :: Kk(I]I(k(JJ: fOR I :: l TO a DO KCtli,IJ :: tCtli,IJ t 1: fOR I :: l TO a DO fOR J :: 1 DO BEGIN APk(I,JJ :: D: fOR [ :: l TO 6 DO APk[I,J] :: APk[I,Jl t tCk[I,[]IPk{[,JJ: END: fOR I :: 1 TO 6 DO fOR J :: 1 TO 6 DO Pkli,JI :: APkli,Jl: {1 STEP 5 : New Caiman Gain natrlx 1} fOR I :: 1 TO 6 DO tklii :: Pkli,61/{Pk[8,61 t Rk}: { Return [alman fitter's Output Estimate 1} Ydnolse :: Xk(6J: {1 LEAST SQUARES PARAftETER ESTiftlTOR 1} If t I td epsilon THEN BEGIN {1 Update Theta's Parameters the estimator's 1} {1 Parmeters Using the ftethod of least Squares 1} 176

PAGE 186

fOR J :: l TO 11. DO {1 Generate denominator : l t Sigma' x P x Sigma 1} BEGIN SlgP[J) :: O: {1 SlgP : Sigma' x P 1} fOR [ :: l TO 11. DO SlgP(JJ :: SlgP(J) t Slgma([JtP(C,JI: END: Denom :: 1: fOR [ :: 1 TO 11. DO Denom :: Denom t SigP(tJSigma(tJ: fOR J :: l TO 11. DO {1 Update P{t d} 1} fOR [ :: l TO 11. DO PU ,t:J :: P(J,[J PsiglJIISigPIU/Denom: Ydsys :: 0 : {1 Generate Error : y Sigma' x Theta 1} for t :: l TO l!. DO Tdsys :: Tdsys t Slgma([)ITheta(t): Error :: Ydnolse Ydsys : fOR J := l TO lb 10 BEGIN Psig[J) :: 0; {1 Pslg : P{t} x Sigma 1} fOR t :: l TO lb tO Pslg(J) :: Psig(JJ t P(J,[)ISigma([); END: fOR [ :: l TO lb 10 Theta(() :: Theta[[) t Pslg[[) Error: END: {1 Update Theta 1} fOR [ :: 1 TO 7 DO {1 Reset the values for U{k} and Y{k} 1} BEGIN U([) :: U[[tlJ: Y[tl :: Y([tll: END: Y(6) :: Ydnolse: U[6) :: fd t Uln: fORt :: 1 TO 6 DO BEGIN Sigma([) :: YI,tl: Sigma([t61 :: Ul,tl: Xltl :: XNEVltl: l[[J ::Theta([): { Reset the values for X{k}, l{k}, 1} { B{k}, and Sigma{k} 1} B([J :: Theta([t61: END: ',yd,' ',Ydnolse,' ',Error}: t :: t td: UNTIL {tf t epsilon c t}: CLOS[{DUAl: END. 177

PAGE 187

APPENDIX C EXAMPLE OF A RESTRICTIVE COMPLEXITY CONTROLLER FOR MODAL ESTIMATION

PAGE 188

179 As shown in chapter 2.2, the behavior of the physical system can be modeled by using the modal supe.rposi tion fotni. 1 y(t,r) = 2:'1li(t)\{'k(r) k=O ( c .1 ) At a stationary location (r) along the structure equation (i3) can be expressed as Yk(t) = Mk(t) sin(GUkt + E1k(t)] for mode k. (C.2) Where Mk(t) is a function of time assuming inherent or active control dampening or additional disturbances during simulation. ek(t) is a function of time with random other certain other types of inputs-otherwise it is a constant. Note, (t)k(t) is a constant with independent modal velocity state space feedback. The following statements can be made about equation (C.2) 1) tis known and Yk(t) is measured. 2) Wk should remain fairly constant, but can be inaccurate with modeling errors. 3) Mk(t) and E1k(t) are the unknowns. Depending on the type of control may or may not be a function of time. The prediction problem is given Yk(O), ,yk(t), t, and 6Uk, find Mk(t+d) and(k(t+d) so that y(t+d) can be accurately predicted. An algorithm that can solve this type of problem is an example of a restricted complexity predictor.