.93;

mm

.
EL

{ME .

J

' L: 33

"3 V9

’5

mm

5.

Ant? '4'”

S

m

u

I."

W0

 

 

 

 

 

 

a»;
’J

u I
~. .firzéa

 

 

 

 

. .L, .«
. . . .
L. 1x.:.1‘;3;.:=

r '51th 2'24?

137.11
.Pv"
v

.v'tf‘: rc-

.

..a,~v
l._
- .

.
=. - .4.

”3:36

«‘iu -.

IT

 

.35
fluid: .3;

a
It

Date

MIIGANCH

\ \llllll ll Lilli Wall ill2 \lllllll

312903

 

 

 

 

 

 

l

 

This is to certify that the

thesis entitled

DESIGN OF AN ACTIVE HYDRAULIC MOUNT
WITH PROGRAMMABLE MECHANICAL IMPEDANCE

presented by

Jerome J. Palazzolo

has been accepted towards fulfillment

of the requirements for

Master of Mechanical
Science degree in Engineering

 

 

 

Major rofessor

July 18, 1994

 

0-7639

MS U is an Affirmative Action/Equal Opportunity Institution

 

 

LIBRARY

M'Chigan State
University

 

 

 

PLACE I! RETURN BOXtomnovombdnckwlm youncord.
To AVOID FINES Mum on or baton dd- duo.

DATE DUE ‘ DATE DUE DATE DUE

       

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU I'M mum WM Oppommlty Institwon
W1

DESIGN OF AN ACTIVE HYDRAULIC MOUNT
WITH PROGRAMMABLE MECHANICAL IMPEDANCE

By

Jerome J. Palazzolo

A THESIS

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

MASTER OF SCIENCE

Department of Mechanical Engineering

1994

ABSTRACT

DESIGN OF AN ACTIVE HYDRAULIC MOUNT
WITH PROGRAMMABLE MECHANICAL IMPEDANCE

By

Jerome J. Palazzolo

The problem of prescribing the mechanical impedance of a one degree of freedom
system is considered. The system consists of a mass that is driven by an external force
and a hydraulic actuator. The goal is to make the system behave as if it were a different
sized mass driven solely by the external force. A two degree of freedom control system
that consists of a feedback controller and a prefilter is proposed as a solution. A
mathematical model of the plant was developed by representing plant nonlinearities with
conic sector bounds. Then, a feedback controller was designed using H2 Optimal Control
Theory. The prefilter, which specifies the size of the emulated mass and compensates for
the amplitude reduction and phase lag of the closed loop system, completed the design.
The control system performance was analyzed using frequency domain criteria from
Robust Control Theory. Time domain simulations were then used to investigate how

well the frequency domain analysis predicted the time domain behavior.

Copyright

Copyright by

JEROME JOSEPH PALAZZOLO

1994

Dedication

To my parents, Jerome and Linda

iv

Acknowledgments

First, I would like to sincerely thank my faculty advisor, Dr. Philip M.
FitzSimons. His advice and friendship since the senior year of my undergraduate study
has allowed me to expand my interest in the fields of Dynamic Systems and Control
Theory. His enthusiasm for our project has taught me that research, although
challenging, is one of the most rewarding tasks undertaken by scientists and engineers.

I would also like to extend my graditude to my other committee members, Dr.
Clark J. Radcliffe and Dr. Ronald C. Rosenberg. Their support of my intellectual growth
began when I was an undergraduate and their suggestions for improvement will always
be appreciated.

The funding for this project was provided by the Rotorcraft Aeroelasticity Group
at the NASA Langley Research Center. The advice given to me by the group has been
appreciated. In particular, Paul Mirick and William Yeager, Jr. have been valuable
sources of information during this project.

Finally, I would like to give a special thanks to my family and friends. I will
always be grateful for their support of my educational pursuits. Although the support
they have given me has not always been technical, I feel that it has made me a better

person and therefore a better engineer.

TABLE OF CONTENTS

LIST OF TABLES .......................................................................................................... vii
LIST OF FIGURES ......................................................................................................... viii
CHAPTER 1. INTRODUCTION .................................................................................. 1
CHAPTER 2. MODELING THE HYDRAULIC SYSTEM ......................................... 7
2.1. Servovalve Dynamics .................................................................................. 7
2.2. Hydraulic Actuator Dynamics ..................................................................... 9
2.3. Mathematical Model of the System ............................................................ 11
CHAPTER 3. DESIGN OF A ROBUST CONTROL SYSTEM ................................... 16
3.1. Feedback Controller Design ........................................................................ 16
3.2. Closed Loop Transfer Characteristics and the Prefilter Design .................. 19
3.3. The MA Structure ....................................................................................... 20
3.4. Robust Stability and Performance ............................................................... 24
3.5. Implementing the Tests for Robust Stability and Performance .................. 26
3.6. MatLab Numerical Results .......................................................................... 29
CHAPTER 4. TIME DOMAIN SIMULATIONS ......................................................... 36
CHAPTER 5. CONCLUSIONS AND RECOMMENDATIONS ................................. 47
APPENDIX A. MATLAB CODE ................................................................................. 49
APPENDD( B. FORTRAN CODE ................................................................................ 58
APPENDIX C. TIME DOMAIN RESULTS ................................................................. 81
BIBLIOGRAPHY ........................................................................................................... 85

vi

LIST OF TABLES

Table 4.1 Outline of Time Domain Simulations ............................................................ 4O

vii

LIST OF FIGURES

Figure 1.1. Schematic Diagram of ARES I ..................................................................... 2

Figure 1.2. Schematic Diagram of ARES 1.5 ................................................................. 3

Figure 1.3. Schematic Diagram of ARES II .................................................................... 4

Figure 1.4. Schematic of the Simple One Degree of Freedom Problem ......................... 6

Figure 1.5. Two Degree of Freedom Controller Topology ............................................. 6

Figure 2.1. Schematic Diagram of a Servovalve ............................................................. 8

Figure 2.2. Schematic Diagram of the Hydraulic Actuator ............................................. 10
Figure 2.3. Conic Sector Bounds of a Nonlinear Function of One Variable .................. 12
Figure 3.1. The Internal Model Control Structure .......................................................... 17
Figure 3.2. The Classical Control Structure .................................................................... 17
Figure 3.3. The MA Structure ......................................................................................... 21
Figure 3.4. Incorporating Uncertainty into a State Space Model .................................... 22
Figure 3.5. Robust Stability Test for Method 1, wc=100 Hz .......................................... 31
Figure 3.6. Robust Stability Test for Method 1, we =200 Hz .......................................... 31
Figure 3.7. Robust Stability Test for Method 2, mc=100 Hz .......................................... 33
Figure 3.8. Robust Performance Test for Method 2, mc=100 Hz, M=50 lb ................... 33
Figure 3.9. Robust Performance Test for Method 2, wc=100 Hz, M=150 lb ................. 34
Figure 3.10. Robust Stability Test for Method 2, wc=200 Hz ........................................ 34
Figure 3.11. Robust Performance Test for Method 2, wc=200 Hz, M=50 lb ................. 35
Figure 3.12. Robust Performance Test for Method 2, (05200 Hz, M=150 lb ............... 35

viii

Figure 4.1 Enforcing a Saturation Limit in the Derivative of State 3 ............................. 37

Figure 4.2. Time Trace with we = 200 Hz, M = 501b, f = 1 Hz ..................................... 41
Figure 4.3. Time Trace with we = 200 Hz, M = 50 lb, f = 10 Hz ................................... 41
Figure 4.4. Time Trace with we = 200 Hz, M = 50 lb, f = 25 Hz ................................... 43
Figure 4.5. Time Trace with we = 200 Hz, M = 501b, f = 50 Hz ................................... 43
Figure 4.6. Time Trace with we = 200 Hz, M = 150 lb, f = 1 Hz ................................... 44
Figure 4.7. Time Trace with we = 200 Hz, M = 150 lb, f = 10 Hz ................................. 44
Figure 4.8. Time Trace with we = 200 Hz, M = 150 lb, f = 25 Hz ................................. 45
Figure 4.9. Time Trace with we = 200 Hz, M = 150 lb, f = 50 Hz ................................. 45
Figure C]. Time Trace with we = 100 Hz, M = 501b, f = 1 Hz .................................... 81
Figure C.2. Time Trace with we = 100 Hz, M = 501b, f = 10 Hz .................................. 81
Figure C.3. Time Trace with we = 100 Hz, M = 501b, f = 25 Hz .................................. 82
Figure C.4. Time Trace with we = 100 Hz, M = 501b, f = 50 Hz .................................. 82
Figure C.5. Time Trace with we = 100 Hz, M = 150 lb, f = 1 Hz .................................. 83
Figure C.6. Time Trace with we = 100 Hz, M = 150 lb, f = 10 Hz ................................ 83
Figure C.7. Time Trace with we = 100 Hz, M = 150 lb, f = 25 Hz ................................ 84
Figure C.8. Time Trace with we = 100 Hz, M = 150 lb, f = 50 Hz ................................ 84

ix

CHAPTER 1. INTRODUCTION

In order to develop safe and reliable helicopter rotor systems, a significant amount
of time and resources must be invested in the analysis and testing of scale rotor
assemblies. At the NASA Langley Research Center (LaRC), the Rotorcraft
Aeroelasticity Group (RAG) conducts much of its research in the Transonic Dynamics
Tunnel (TDT) facility. This facility is ideally suited for scale rotor tests because of its
unique ability to use a heavy gas as its working medium. By scaling the model
appropriately, the results can be readily extrapolated to full scale. The primary scaled
rotor test stand used by the RAG has been the Aeroelastic Rotor Experimental System I
(ARES I) displayed in Figure 1.1. The rotor system is mounted on a fixed balance that
has one pitching degree of freedom at the longeron. The longeron is mounted on the
model stand which, in turn, is bolted to the floor of the TDT. Although the rotor system
has been isolated from the helicopter fuselage in this configuration, ARES I is used to
investigate rotor dynamics, operating loads, and the performance of the rotor. Another
test stand, ARES 1.5, is available to investigate specific rotor/body coupling. It allows
for pitch and roll of the test balance relative to the longeron and uses elastomeric springs
and rotary viscous dampers to control the stiffness and the damping of the pitch and roll
axes. A schematic of this test stand is shown in Figure 1.2.

A third generation test stand, ARES II, is currently being developed to extend the
experimental capabilities of the RAG and is displayed in Figure 1.3. It is unique in that
the rotor system is mounted on a six degree of freedom balance. The longeron, which
serves as a fixed table, has a degree of freedom at the pitch pivot and sets the reference

angle of the coordinate system of the test balance. This pitching degree of freedom

1

2
allows the RAG to test the rotor system over various forward flight regimes. By sensing

the forces and moments generated by the rotor at the test balance, it is possible to make
the platform emulate a scale fuselage by modeling the dynamics of the actual helicopter.
This will enable the RAG to investigate the full rotor/body coupling dynamics of scaled
rotorcraft. In order to emulate a scale helicopter fuselage, the positioning system of the
test balance must have a high power rating so that the platform is able to operate over a
wide frequency bandwidth. Due to this requirement, six hydraulic actuators connect the
test balance to the longeron. The length of these actuators determine the position and
orientation of the test balance relative to ground. This six degree of freedom parallel

linkage is commonly known as 3 Stewart Platform [Nguyen 1991].

  
   
  

Balance

 

 

Pitch Pivot
‘ Longeron
Model Stand
Q.)

Figure 1.1. Schematic Diagram of ARES I

 

 

 

 

 

Roll Damper .... Pitch Damper

Roll Spring

 

 

Figure 1.2. Schematic Diagram of ARES 1.5

Hydraulic
A Actuators

 

Hydraulic
Actuators

 

 

Pitch Pivot

 

Model Stand

Figure 1.3. Schematic Diagram of ARES II

5
The success of ARES II hinges on the design of an actuation system with

programmable mechanical impedance characteristics. To demonstrate this concept using
ARES II hardware, the problem of prescribing the mechanical impedance of a simple one
degree of freedom system will be considered. In Figure 1.4, the mass (m) is driven by the
external force (F,,,) and the force due to a hydraulic actuator (Fm). The goal is to
prescribe the signal driving the hydraulic actuator such that the system behaves as if it
were another mass (M) driven solely by the external force. The proposed solution to this
problem (Figure 1.5) is to use a two degree of freedom control design strategy. First, a
high gain feedback controller, C, will be designed using H2 Optimal Control Theory to
precisely position the actuator over a wide frequency range. A prefilter, PF, that specifies
the size of the emulated mass and compensates for the amplitude reduction and phase lag
introduced by the nominal closed loop system will then be designed to improve the
system's performance. After the control system has been developed, Robust Control
Theory will be used to ensure the closed loop system performs robustly over the
frequency range of interest. Time domain simulations will then be used to investigate
how well the frequency domain analysis predicted the time domain behavior. Before the
robust controller and prefilter can be determined, a high fidelity model of the plant, P,

must be developed.

 

 

 

 

 

 

 

 

1:“act '7 m Fext
M I FCXI
WW

Figure 1.4. Schematic of the Simple One Degree of Freedom Problem

F ext

 

 

 

 

Figure 1.5. Two Degree of Freedom Controller Topology

CHAPTER 2. MODELING THE HYDRAULIC SYSTEM

Determination of the nominal plant model will be carried out in three phases. In
the first phase, the servovalve dynamics will be modeled. Next, a model of the hydraulic
actuator will be developed. In particular, the method of conic sector bounding will be
used to define a set of linear systems that capture all of the nonlinear behavior of the
actual system model. Once these mathematical models are complete, the results will be
combined to define the plant model of the overall system in state space form. From these

results, an equivalent transfer function representation will be defined.

2.1. Servovalve Dynamics

The servovalve being modeled, shown in Figure 2.1, has flapper force feedback.
An input current to the torque motor causes the flapper to change its position. This
results in a pressure differential across the spool forcing it to move from the null position.
The null position is defined to be the center position of the spool. Since the spool valve
under consideration is critically lapped, any movement of the spool from the null position
will result in fluid flow through the valve. The flapper also creates a feedback torque that
opposes the input torque of the motor. This feedback causes the flapper to return to the
null position whereas the required torque balance is achieved by spool displacement.

A relationship between the input current and the spool position of the servovalve

will now be derived. The resultant torque on the flapper can be given by [Watton 1989]

M — kwxe = z, (2.1)

where i is the input current

k1 is the torque motor current gain

7

   
    

Torque Motor

Flapper

Spool

Figure 2.1. Schematic Diagram of a Servovalve

kw is the feedback wire stiffness
x, is the spool position
Tfis the resultant flapper torque
The equation of motion for the first stage of the flapper can be written in the form
x}. + 2493, + a): = (w: /k,)r, (2.2)
where kfis the net stiffness of the flapper
xf is the displacement of the flapper
w" is the natural frequency of the first stage of the flapper
C is the damping ratio of the first stage of the flapper
The relationship between the position of the flapper and the velocity of the spool is
defined by
kzxf = Aejre (2.3)
where A, is the spool end area

k2 is the hydraulic amplifier flow gain

9
Equations (2.1), (2.2), and (2.3) can now be combined to determine the relationship

between the spool position and the input current. The resulting equation is
365 + (12555 + alfce + aoxe = boi (2.4)
where do: (k2wfkw / Aekf)
at]: w:
(12: 2§we

b0=(k,w3k, /A,k,)

2.2. Hydraulic Actuator Dynamics
The next phase of the modeling process to be considered is the dynamic behavior
of the hydraulic actuator. This system is shown in Figure 2.2. The input current to the
servovalve causes the spool position to change which, in turn, creates a pressure
differential across the hydraulic actuator. This pressure differential, accompanied by
fluid flow, results in a change in the position of the hydraulic actuator. A simple force
balance on the actuator results in the following equation:
mite = —ae5ce + PlA, - Per (2.5)
where A1 is the piston side area of the actuator
A2 is the rod side area of the actuator
m is the load mass
P1 is the pressure on the piston side of the actuator
P2 is the pressure on the rod side of the actuator

xa is the position of the actuator

ad is a damping term due to the viscosity of the hydraulic fluid

10

 

 

 

 

 

 
   
 

 

 

 

Torque Motor /
Servovalve
Pl, Ql 1D2 , Q2 Load Mass
Hydraulic _____ m
Actuator

 

 

 

W

Figure 2.2. Schematic Diagram of the Hydraulic Actuator

The amount of fluid flow on sides 1 and 2 of the actuator depends on the spool position

of the servovalve. This relationship is given by [Watton 1989]

Ql = kexe(u(xe)./Pe — Pl + u(-x5)\/F,) (2.6)
Q, = k,x_,(u(x,),/172 + u(-xe)1/Pe — P2) (2.7)

{0: x < O} '
u(x) = (2.8)

l: x 2 0

where kx relates the spool displacement to volume flow rate

P3 is the supply pressure (assumed constant)

Q1 is the volume flow rate on side 1 of the actuator

Q2 is the volume flow rate on side 2 of the actuator
The switching function, u(x), is needed to model the role reversal that occurs when the
servovalve goes from extending the load mass to retracting it. Namely, one of the supply
lines in Figure 2.1 is cut-off by the spool valve to allow the actuator to move in one

direction.

11
The volume flow rates also depend upon the velocity of the actuator and the

compressibility of the hydraulic fluid. These relations are given by [Watton 1989]

 

QI = Alia +[Vh + A'xa )PI (29)
l3
, (V, +A,(1e -xe)].
Q2 = Azx. - [3 P2 (2.10)

where 1,, is the stroke length of the actuator

V), is the volume of the fluid in the hose from the servovalve to the actuator
[3 IS the bulk modulus of the hydraulic flu1d (4175—17)

The bulk modulus is a measure of the compressibility of a fluid. This fluid property
arises when it is determined that the fluid cannot be considered incompressible. Since the
control system developed must operate over a wide frequency bandwidth, the bulk
modulus of the hydraulic fluid must be considered to maintain accurate set-point tracking

and to ensure that no ill effects arise due to the hydraulic “spring” in the system.

2.3. Mathematical Model of the System
The equations of motion for this system are nonlinear and time invariant. In total,
seven states are needed to define the dynamic behavior of the system. They are defined

as

T . . .. T
x=[x1 x2 x3 x4 x5 x6 x7] =[xee xe, P1 P2 xs x3 x5] (2.11)
The nonlinear set of state space equations is

1- x2 -

—(ae,/m)x2 +(A1/m)x3 - (AZ/m)x4
l111(351’112’153’3‘5)

x: 'I’z(xl,x2,x4,x5) (2.12)

x6

x7

 

 

-—a2x7 - alx6 - “0355 + boi

12

 

 

 

 

Figure 2.3. Conic Sector Bounds of a Nonlinear Function of One Variable

where '1" = (fly—[4J2 + keex5 (u(x5 )1/Pe - x3 + u(-x5 NZ» (2.13)
h l 1

 

fl
'2”, = ( V,, + A2(l _ x] )](A,x2 — keex5 (u(xs )Jx—4 + u(—x5 )ije — x4 )) (2.14)

a

Because of the complicated nature of the functions ‘1’; and '15, it was decided that a
simple linearization may not capture all of the dynamic behavior of the system. In order
to ensure that all of the system dynamics are included in the mathematical description of
the system, the method of conic sector bounding was used [FitzSimons 1994]. For a
simple nonlinear function of one variable, conic sector bounding is equivalent to finding
the upper and lower lines whose swept area contain the complete nonlinear behavior of

the function. These bounding lines are displayed in Figure 2.3. The goal of the conic

sector bound analysis is to represent the functions 'PI and ‘I’z in a parametric linear form,

13
i.e., the nonlinear function is replaced by a linear function that has a range of parameter

values. The linear form of the two nonlinear functions can be defined by

'P, = auxl + aux2 + c1331:3 + a35x5 (2.15)

'1’2 = aelxl + (142x2 + aeex4 + “45155 (2.16)

where 2.; 5 ae. S (7,}. were obtained from the conic sector bound method

The nominal parameter values are defined by

ae. = —2— (2.17)

whereas the amount that the upper and lower bounds differ from the nominal parameter

values is given by
Aa.. = ——"——"—’- (2.18)

After replacing ‘1’, and '2”, with 52’, and ‘I’z in Equation (2.12), the state space form of

the nominal plant can be written as
x = Apx + Bpi (2.19)
y = Cpx + Dpi (2.20)

where

 

O l 0 O O O O
O 022 a23 (124 O O O
031 “32 033 0 “35 0 0
AI) = 041 042 O (144 (145 O O
O 0 O O 0 l 0
O O O O 0 O l
O O O 0 a75 (176 077‘

 

B=[0 0 0 0 0 Obnlr

14

and a22 = —ae / m
a23 = A1 / m
“24 " TA2 / m
“75 = T“o
“76 = T“l
“77 = T“2
“71 = “o

The state space equations for the nominal plant can also be expressed by an

equivalent transfer function. This is accomplished by applying the familiar equation

T(s) = Cp(sI — Ap)"Be, + DP (2.21)
The software package Mathematica was used to determine the symbolic representation of

the transfer function. The results are as follows:

nls+n0
S7 +d6S6 +£15.95 +d4S4 +d3S3 +d2S2 +dls+d0

 

T(s) = (2.22)

where nl = b71(a23a35 + (124045)
"0 = Tb71(“23“35“44 + “24“45“33)
d6 = 51 T “77
d5 = 62 T “7761 T “76
“4 = 53 T “7752 T “7651 T “75
“3 = “4 T “7753 T “7662 T “7561
“2 = T“77“4 T “7653 T “7552
“1 = T“7654 T “7553

“0 = T617564

15
and 6, = —(a22 + a33 + a“)
“2 = T(“23“32 + “24“42 T “22“33 T “22“44 T “33“44)
53 = “23“32“44 + “24“42“33 T “22“33“44 T “23“31 T “24“41
54 = “23“31“44 + “24“4l“33
Now that a high fidelity model of the hydraulic system has been developed, the design of

the two degree of freedom control system may begin.

CHAPTER 3. DESIGN OF A ROBUST CONTROL SYSTEM

Since the plant model has been determined, the design of the two degree of
freedom controller outlined in Chapter 1 can begin. First, H2 Optimal Control Theory
will be used to define a closed loop controller (C) that satisfies criteria for internal
stability and nominal performance. Once the controller has been designed, a prefilter
(PF) will be defined that allows the hydraulic mount to have programmable mechanical
impedance. In order to compensate for the amplitude reduction and phase lag in the
output of the nominal system, knowledge of the closed loop system will also be used in
the prefilter design. Once the two degree of freedom control system has been determined,

tests for robust stability and robust performance will be conducted.

3.1. Feedback Controller Design

One of the building blocks of Robust Control Theory is the Internal Model
Control (IMC) structure. This structure is shown in Figure 3.1. The error signal for the
IMC controller depends on how the actual plant (P) varies from the nominal plant (P).
Although the controller will be defined with the IMC structure, it can easily be converted
to an equivalent classical controller, C, as shown in Figure 3.2 with the control signal

disturbance, u’, set equal to zero. The relationship between C and Q is given by

Q
=__,,_ 3.1
1_PQ ( )

When designing a robust controller with the IMC structure, it is important to first

guarantee internal stability and nominal performance. Internal stability is guaranteed

16

 

 

 

 

 

 

 

 

%
I o
I
I

 

"Ut

 

 

 

 

 

Figure 3.1. The Internal Model Control Structure

 

 

 

 

 

 

 

 

 

Figure 3.2. The Classical Control Structure

when bounded signals injected at any point of the control system generate bounded
responses at any other point. As shown in Figure 3.2, the control system has two
independent inputs (r, u’) and two independent outputs (y, u). For the classical feedback
structure, internal stability is achieved if and only if all elements of the transfer matrix
relating the independent inputs to the independent outputs have all their poles in the open

left-half plane (LHP) [Morari 1989]. The following transfer matrix can be found with

 

block diagram algebra:
PC P
y r
= 1+ PC 1+ PC (3 2)
C - PC u’ '

 

1+PC 1+PC

18
Before checking if the condition for internal stability has been met, the classical

controller must be defined. The first step in the design process is to determine the
controller that optimizes the performance of the nominal plant. The performance criteria
is to keep the H2 norm of the error signal, e, small with respect to the input signal, r.
[Morari 1989]. The difficulty in specifying nominal performance criteria is quantifying
what size error signal can be considered small and what type of inputs make up the set of
external inputs. The criteria for nominal performance of a single-input single-output
system is guaranteed by the design of an H2 optimal controller. The condition that the
IMC controller must satisfy in order to be considered H2 optimal can be given by [Morari

1989]

ngnllellz =ngnll(1—f»€1)wll, (3.3)

where ||e||2 is the integral-square-error (ISE) jjfez(t)dt

p is the transfer function that describes the nominal plant
a is the transfer function that describes the H2 optimal controller
w is the weighting applied to the input signal
Since it has been assumed that the conic sector bounds will be chosen such that the
transfer function of the nominal plant is stable and minimum phase, the nominal
performance criteria is satisfied when
a=i4 (in
Due to the nature of the controller defined in Equation (3.4), the H2 optimal controller is
independent of the weighted set of specified inputs to the system. This is shown by
substituting Equation (3.4) into Equation (3.3) since the minimal ISE will be zero with
any input signal weighting chosen.
The first step towards achieving internal stability is to make the H2 optimal
controller physically realizable by augmenting the controller, 6;, with an IMC filter, f.

Since the transfer function defined in Equation (2.22) has six more poles than zeros, its

l9
inverse will have six more zeros than poles. In order to make the controller physically

realizable, a sixth-order filter is used to make the transfer function of q proper, i.e., the

number of poles equals the number of zeros. The IMC controller is therefore given by

q = 51f (3.5)
1
where f = -———.
(1s + 1)
and A is the filter tuning parameter

The filter tuning parameter is used to trade-off system performance and robustness by
limiting the bandwidth of the controller. By Equation (3.1), the equivalent classical

controller can be given by

.. f
“111—31 (3'6)

Now that the classical closed loop controller has been defined, the poles of the individual
transfer functions defined in Equation (3.2) can be found. If all of the real parts of these
poles are less than zero, the conditions for internal stability and nominal performance will

have been met.

3.2. Closed Loop Transfer Characteristics and the Prefilter Design

Now that a classical controller has been defined for the hydraulic system, a
prefilter can be designed that will not only allow the user to specify the mechanical
impedance of the system, but will also compensate for the amplitude reduction and phase
lag of the nominal closed loop system. For the simple one degree of freedom problem
displayed in Figure 1.4, a mass, m, is acted upon by an external force and a hydraulic
actuator. The goal is for the system to behave as if it were another sized mass, M, being
forced solely by the external measurable force. The equation of motion that describes the

dynamic behavior of this system is given by

(3.7)

20
Therefore, the transfer function relating the external force to the position of the emulated

mass, M, is defined as

X(s) = 1
Fe,,(s) Ms2

 

(3.8)

From Equation (3.2), the transfer characteristic from the input, r, to the output, y,

is given by

PC

If it is assumed that the actual plant, P, can be approximated by the nominal plant, P,

over the frequency range of interest, Equation (3.9) simplifies to

1

CL = f = ——
(S) (is +1)6

(3.10)

Since the denominator of Equation (3.8) is of second order, the nominal closed loop
amplitude reduction and phase lag can be compensated by a second order approximation
to f '1, namely

f“ z 15.12..~2 +6its+l (3.11)
Although the effect of Equation (3.11) is more pronounced at lower frequencies, it will
still provide some compensation for amplitude reduction and phase lag at the higher
operating frequencies of the control system. In order to obtain the desired effects of
programmable mechanical impedance and improved closed loop performance, the

prefilter was defined as

1512.92 + 615 +1
Ms2

 

PF(S) = (3.12)

3.3. The MA Structure
With Robust Control Theory, the MA structure displayed in Figure 3.3 is used to

test the system for robust stability and robust performance. The power of Robust Control

Theory lies in its ability to check whether or not the closed loop system meets

21

>7
1

 

 

[Q

 

 

 

 

 

:17

 

 

 

 

 

Figure 3.3. The MA Structure

requirements for stability and performance when uncertainty exists in the model of the
plant. By satisfying the requirements for internal stability and nominal performance, the
only thing that has been proven is that the controller can perform well if the actual plant
is equivalent to the nominal plant. However, the MA structure can be used to test the
controller's performance for an entire family of plants that is described by the model
uncertainty.

The model uncertainty is given by the range of conic sector bound parameters,

i.e., Q.)- S “a S aej, that can be used to describe the behavior of the hydraulic system.

Recall, the model uncertainty was divided into nominal and differential values in

Equations (2.17) and (2.18)

alj + gij
2
21' —a
_. 'J —U
2

In order to represent the uncertainty that exists in the parameters of the A-matrix of the

state space representation of the nominal plant, the MA structure shown in Figure 3.4 was

investigated.

22

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

C(s)

 

 

 

 

 

 

 

 

Figure 3.4. Incorporating Uncertainty into a State Space Model

Before defining the transfer characteristics of the MA structure, several aspects of
Figure 3.4 must be described. In order to test the system for robust stability, the
uncertainty that exists in the matrix Ap is structured into a diagonal matrix, AA. The
amount of uncertainty to be considered will be varied by specifying the percentage, (1) , of
the differential parameter value that will be considered in the analysis. Also, the
weighting matrices B11 and Cu will be used to normalize the uncertainty matrix so that
each diagonal entry in AA varies between [-1, +1]. In particular, these relationships and

the symbols given on Figure 3.4 can be given by

F -

 

 

A3, 0 o o 0 o o o
0 A32 0 0 0 0 0 0
0 0 A33 0 o 0 0 o
0 0 0 A35 0 0 0 0

AA: 0 0 o 0 A41 0 0 0 (3-13)
0 o o o 0 A42 0 0
0 0 0 o 0 0 A44 0
_ 0 0 o o o o 0 A45

 

 

' O 0 0 O
0 0 0 0
“31d “32.1 “33d “35d
B): = O O O O
0 0 0 O
O O 0 O
_ O 0 O O
P“31d 0 0
0 “32d 0
O O a33d
O 0 0
C“ = “41d 0 0
0 “42d 0
O O O
_ O O 0
where AU e[—l, +1]
“7,21 = q) . Aaej
and e = r’ T y
r = Far
r’ = xdcs

x = Apx+Bpi+Buul
y=xa

”O = Cu"

23

O O O O
O 0 O 0
O O O O

 

“41d “42d “44d “454
0 0 0 0
0 0 o 0
0 0 0 0 j

0 0 0 0 ‘

0 0 0 0

0 0 0 0

0 6235,, o 0

0 0 0 0

0 o 0 0

a44d O O O

0 (145,, 0 0 j

 

(3.14)

(3.15)

In order to define the state space representation that relates the inputs of M (111 and r) to

its outputs (110 and e) the state space representations of the plant, controller, and prefilter

must be formally defined. The dynamics of the plant are described by

x, = Apartl + Bei + Ben,

“o
x

a

H

C
C

p

P

l...

(3.16)

24
In a similar manner, the controller and the prefilter are given by the following state space

representations:
x =A x +B e
.2 c 2 c (3.17)
1=ch2+Dce
x =A x +B r
3 r 3 r (3.18)

r' = CfX3 + Dfr
After manipulating Equations (3.16), (3.17), and (3.18) so that the states only depend on

the inputs it! and r, the state space representation of M can be given by

x, Ap-Bpnccp BpCe BpDch x, BlLl Benchf

x, = —Bccp Ac 13ccf x2 + o 13,1)f [”1]
e
x 0 o A x o B
3 r 3 r (3.19)
uo_c, 0 0 “+0 0 111
e"-cp0c,"2 onfr
X3

When considering robust stability and performance, it is convenient to think of the

individual transfer functions that make up the M structure, i.e.,
[”0] = MPH] ___ [M11 M12 :l[u1:l (3.20)
e r M21 M22 r

3.4. Robust Stability and Performance

In order to describe the family of plants in the least conservative manner, the
development of the representation of the MA structure has been conducted such that the
uncertainty can be represented in a structured manner. The conditions for robust stability
for structured uncertainty are provided in Chapter 11.2 of [Morari 1989]. Before
discussing the requirements for robust stability, several terms need to be defined.

The singular values of a complex (nx m) matrix A are the non-negative square

roots of the eigenvalues of AHA , i.e.,

0,.(A) = ,1,(A”A) (3.21)

25

where i = 1, 2, min(n, In)

and H denotes the complex conjugate transpose of A

In the development of the conditions for robust stability, the maximum singular value, 6 ,
is used to specify the norm of a matrix induced by the vector 2-norm. An alternate

definition of the maximum singular value is

 

A
3(A)=max| Xlz
x¢0 IXI2

(3.22)
Therefore, the geometric interpretation of the maximum singular value is that it is the
least upper bound on the magnification of a vector, x, by a matrix (or operator), A.

In order to define the structured singular value, the set of structured uncertainties

of interest must be defined. For the hydraulic system, we have

X1) = {All = diag(A3,,A32,A33,A35,A41,A42,A44,A45) I 6(Aij) S v} (3.23)

The structured singular value, u(M), is defined such that p-1(M) is equal to the smallest

5(Au) needed to make det(I — MA“) = 0. If no Au that destabilizes the system exists,

then u(M) E 0. A formal definition is given by [Morari 1989]
pT](M) = min{v I det(I — MA“) = O, for some AILI e Xv} (3.24)
I)

The system is considered to be robustly stable for all perturbations Au 6 X ”:1 if and only
if [Morari 1989]
u(M11(iw)) < 1 We (3.25)

According to the definition of the structured singular value, any uncertainty that

destabilizes the system must satisfy v >1: 6(Aem) > 1. Therefore, the uncertainty that
exists in the system, Au 6 Xv=1, cannot destabilize the control system.

The conditions for robust performance are provided in Chapter 11.3 of [Morari
1989]. With the MA structure, robust performance is guaranteed when the closed loop
system is augmented with a performance perturbation block and the modified system

satisfies the conditions for robust stability. As shown in Figure 3.3, the block

26

perturbation for the uncertainty, A“, is augmented with a perturbation, Ap’ for the
performance specifications

A = diag(Au, AP) (3.26)
The system meets the conditions for robust performance for all perturbations A 6 X2, ”:1

if and only if [Morari 1989]
u(M(iw)) <1 Vw (3.27)

6(Aeee) S v}

Satisfying the condition for robust performance is equivalent to having a robustly stable

where X2. e = {A = diag(Au, AP)

 

system whose transfer function from the input, r, to the output, e, is bounded by 1 at all

frequencies. That is,

“(Mll) < 1
u(M) <1 c» and (3.28)

llT L. <1

ml
However, if the structured singular value of the MA structure is greater than one over a

certain frequency range, alternative performance limits are given by

u(M) < fie u(M,,)<1 for Au 6 X1243"
and <=> and (3.29)
[3. >1 lITmIL. < 13.

Equation (3.29) states that the only way to ensure that the system is robustly stable and

the transfer function is bounded by ,8" is to reduce the plant model uncertainty by [3“.

3.5. Implementing the Tests for Robust Stability and Performance

With the conditions for robust stability and performance defined, the most
efficient way to determine the structured singular value was to use the u-Analysis and
Synthesis Toolbox developed for the software package MatLab. This toolbox is a
collection of functions that were developed for the analysis and synthesis of robust

control systems. In particular, the MatLab function mu calculates the structured singular

27
value of the MA structure at distinct frequency points. However, several aspects of the
description of the system had to be addressed before using the toolbox.

The MatLab scripts that were written to test the system for internal stability,
robust stability, and robust performance are included in Appendix A. The script,
hssvpfm, is the main program that sets up the system and calls the various functions to
determine whether the system meets the specifications previously outlined. After the
constants of the actuator and the servovalve are defined from data provided by the
manufacturers, the results of the conic sector bound analysis are read from data files to
specify all of the other system parameters. The program then defines the state space
representation of the nominal plant, Equation (3.16). The classical controller is then
defined by using Equation (3.6). The function, intstab.m, was then written to test the
nominal system for internal stability. If the transfer matrix defined by Equation (3.2)
contains any poles that are not in the open LHP, the program terminates before the tests
for robust stability and performance can begin.

One of the problems encountered while trying to describe the system with its state
space representation and the equivalent transfer function was numerical ill-conditioning.
This problem arose because of the wide range of parameter values that are used to
describe the system. Specifically, the smallest value is given by the torque motor current
gain, k1=O.O25 in-lbs/mA, while the largest value is given by the bulk modulus of the
hydraulic fluid, B=250,000 psi. This results in a transfer function whose coefficients vary
by seventeen orders of magnitude. Once the state space representation of the system is
defined, an ordered balanced realization of the system was completed by using the
MatLab function, obalreal. This function defines a similarity transformation matrix, T,
that gives an equivalent state space representation of the system. This transformation can

be given by

[A, B, c, D] = [TAT-1, TB, CT“, 1)] (3.30)

28

However, the states of this system are ordered, resulting in a system whose numerical
conditioning is superior to its original representation. After balanced realizations are
conducted on the state space representations of the plant, the controller, and the prefilter;
the M-structure is defined and its frequency response is determined at specified
frequencies.

In order to calculate the structured singular value, 1.1, efficiently, the u-Analysis
and Synthesis Toolbox determines its upper and lower bounds. These bounds are easier
to calculate then the actual structured singular value and, as it will be shown, usually
converge to a tight description of the structured singular value. The lower bound is
provided by the spectral radius of M, p(M), while the upper bound is provided by
3(M). The optimal bounds are determined by using transformations that do not affect
p(M), but reduce the gap between p(M) and 3(M). First, the spectral radius is defined
by [Morari 1989]

p(M) = maxlle-(M)| (3.31)
l
The transformations of interest are defined as [Morari 1989]

Q={QeA=Q”Q=I.}

(3.32)
D = {D e A : diag(d,~) where d,- e 91"}

With these transformations defined, the bounds on the structured singular value are given

The MatLab function mu contains an algorithm that determines the upper and lower
bounds defined by Equation (3.33). After defining the MA structure in MatLab, the test
for robust stability was completed by finding u(M“) over a specified frequency range

with an uncertainty defined by A“. Similarly, robust performance was determined by

p(M) with the uncertainty A.

29
3.6. MatLab Numerical Results

Before conducting the tests for robust stability and performance, the uncertain
elements in the plant description must be specified by a conic sector bounding method.
Once their nominal values are determined, the nominal plant and classical controller can
be defined and the structured singular value can be determined by the software package
MatLab. The first method involves a set of conic sector bounds that only considered the
variation of the nonlinear function with two of its dependent variables. In particular, this

parametric linear form is given by

A

‘1’, = aux2 + a35x5 (3.34)
‘I’2 = aux2 + “45x5 (3.35)
The slope values that bound the nonlinear functions were approximated by defining the
proposed operating ranges for the actuator position and the pressures of the system and
considering their effect on Equations (2.13) and (2.14). With the nominal plant defined,
the only parameter that needed to be defined was the detuning parameter of the filter, 2..

Since the IMC filter is low-pass, the detuning parameter can be defined by specifying the

cut-off frequency of the filter, we. In particular,
II = — (3.36)

The closed loop system was tested for robust stability and performance for cut-off
frequencies of 100 Hz and 200 Hz. The performance of the control system can be judged
over a range of emulated masses because the prefilter has been designed to achieve the
goal of programmable mechanical impedance. Since the actual mass was equal to 100
pounds, emulated masses of 50 pounds and 150 pounds were considered in this study.

In the MatLab p—Analysis and Synthesis Toolbox , the uncertainty can be specified
as real, complex, or mixed. Due to the increased difficulty of the pure real structured
singular value problem, MatLab warns that some convergence problems may occur.

Indeed, this was the case for the hydraulic system. However, when the uncertainty is

30

considered to be complex, the results of the analysis are conservative because the true
uncertainty is strictly real. For this reason, only twenty-five percent of area swept out by
the conic sector bound was included in this analysis. Although a portion of the nonlinear
function may be left out with this subset, it was hoped that adding the complex
component to the uncertainty would help compensate for the decrease in real-valued
parameter variation. These assumptions will be investigated when the frequency domain
results are compared to the time domain simulations.

The first case to be considered will be the control system with a cut-off frequency
of 100 Hz. The results for the robust stability of the system are shown in Figure 3.5.
Although the structured singular value for the MA structure is less than 1 for most of the
frequency range of interest, the maximum structured singular value is approximately 7.7.
Since the requirement for robust stability is that the maximum structured singular value is
less than one for all frequencies, the effect of increasing the filter cut-off frequency to 200
Hz was investigated (Figure 3.6). Although the structured singular values are
approximately equal over most of the frequency range, the peak has now risen to almost
9.5. For this reason, it was decided to investigate the method of conic sector bounding
that includes variations in all four dependent variables, i.e., Equations (2.15) and (2.16).

In order to determine a set of conic sector bounds with respect to four coordinate
directions, a simple investigation of the nonlinear equations was not considered feasible.
With this under consideration, we developed a constrained optimization to determine the
optimal set of conic sector bounds for a function of four variables. This technique is
related to previous results [FitzSimons 1994]. The constraints in this optimization
ensured that the conic sector bounds captured the entire nonlinear function over a grid of

the dependent variables. The cost function to be minimized for this optimization,

2 o
(p: 27,}.(Aae) , was a measure of the area swept out between the maxrmum and

minimum parameter values in each of the four coordinate directions. The weights, ye,

were used to specify the degree of importance of each of the coordinate directions.

31

SW sewn: Values

 

Frequency [rad/sec)

Figure 3.5. Robust Stability Test for Method 1, we=100 Hz

 

      

 

 

 

 

 

 

 

   

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

A ...., ., . .4 .. “I “”1

70‘ 10" 10° 10' 10’ 10‘

 

 

 

 

 

 

10'
Frequency [radlsoc]

Figure 3.6. Robust Stability Test for Method 1, we=200 Hz

 

32
After the optimal set of conic sector bounds were determined, the closed loop

system was again tested for robust stability. Figure 3.7 displays the structured singular
values for a filter cut-off frequency of 100 Hz. With the optimal set of conic sector
bounds, the maximum structured singular value was less than 0.5. Therefore, the
conditions for robust stability had been satisfied. Then the system was tested for robust
performance for emulated masses of 50 lb and 150 lb. These results are displayed in
Figures 3.8 and 3.9. Although the system is robustly stable, the condition for robust
performance is violated for frequencies less than approximately 0.6 rad/sec. This means
that the error resulting from a unit force input in the lower frequency range is not norm-
bounded by 1. However, Equation (3.29) could be used to establish alternative bounds on
the transfer function Tm. It is also important to remember that the structured singular
value provides a least upper bound on the worst case performance of the control system
and may be conservative.

By increasing the cut-off frequency to 200 Hz, the maximum structured singular
value increased slightly from the last test case (Figure 3.10). However, the advantage of
increasing the filter cut-off frequency for this system becomes apparent when the
conditions for robust performance are considered. The robust performance results for an
emulated mass of 50 lb are shown in Figure 3.11, whereas the results for 150 lb are given
in Figure 3.12. At the low frequencies, the structured singular value decreased from 3 to
approximately 2.2. However, the control systems with cut-off frequencies of 100 Hz and
200 Hz seem to perform approximately the same at the higher frequencies. That is, the
maximum gain from the input force to the error signal is approximately equal to 0.5 for

both systems.

33

Structured ShalerVeluee

  

10‘
Frequency [redieec]

Figure 3.7. Robust Stability Test for Method 2, we=100 Hz

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Frequency [red/sec]

10‘ 1o" 10"

Figure 3.8. Robust Performance Test for Method 2. we=100 Hz, M=50 lb

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

22“.]
I

10 1
Frequency [redeec]

 

 

Figure 3.9. Robust Performance Test for Method 2. we=100 Hz, M=150 lb

summsmuvm

 

Frequency [redIeec]

Figure 3.10. Robust Stability Test for Method 2. we=200 Hz

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

3 __, E _, E 533:2, §§::::::,' :35: : ::::‘
10‘ 10" 10° 10’ 10’ 10

10'
Frequency [redisec]

Figure 3.11. Robust Performance Test for Method 2, we=200 Hz, =50 lb

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

433.1 221' ' 2 322222' 222 23'.‘ 3
" 10° 10' 10' 10 10
Frequency [redieecl

Figure 3.12. Robust Performance Test for Method 2. we=200 Hz. M=150 lb

CHAPTER 4. TIME DOMAIN SIMULATIONS

Chapter 3 dealt with the stability and performance of the closed loop system in the
frequency domain. However, the performance of most systems is easier to quantify if the
system is investigated in the time domain. For this reason, it seemed beneficial to
simulate the prefiltered response of the system. The dynamics of the plant, the controller,
and the prefilter are given by a set of nonlinear ordinary differential equations. When
simulating the dynamics of the hydraulic system, the actual nonlinear equations defined
in Equation (2. 12) were used. The dynamics of the controller and the prefilter are defined
in Equations (3.17) and (3.18). Since the conditions for robust stability and robust
performance were proven for the conic sector bounds of the nonlinear functions, the
controller should also perform robustly with the nonlinear system in the time domain.

In order to solve this set of ordinary differential equations, the Livermore Solver
for Ordinary Differential Equations (LSODE) was used [Hindmarsh 1983]. This software
is written in FORTRAN and is linked with a file that defines both the derivatives of the
states and the Jacobian matrix of the partial derivatives of the states. Appendix B
contains the FORTRAN code, hpfsim.f, which describes the dynamic behavior of the
prefiltered closed loop system. The results of the integration are written to a file in a data
format that is read by MatLab to generate the actual time histories of interest.

Although the dynamics of the hydraulic system are given by Equation (2.12),
there are certain limits that exist in the operating range of the states. For instance, the
load pressures should lie between the supply pressure and the return pressure. This

results in the servovalve having saturation limits. That is, the location of the spool in the

36

37

'1’,(x,,x2,33,x4)

'Pi(x1,x2,_x3,x4)

 

 

 

 

_ _ 1’
£3,101 £3 x3 x3,tol PS

Figure 4.1 Enforcing a Saturation Limit in the Derivative of State 3

servovalve (Figure 2.1) determines the pressure differential which causes the fluid to flow
from the supply lines to the return lines in the hydraulic system. Once the spool moves
all of the way to the right or left, the resulting load pressures should be bounded by the
supply pressure and the return pressure under normal operating conditions. In order to
enforce these operating limits in terms of the pressures, the derivatives of the pressures
will be saturated when the pressures approach the supply pressure or the return pressure.
This is shown qualitatively in Figure 4.1.

Recall, the derivatives of the pressures were defined in Equations (2.13) and

(2.14). Specifically

'11,(x, ,x2,x3,x5) = (__fl__](—A,x2 + kxx5(u(x5)1/Ps - x3 + u(—x5 )Jg»

VII '1' Alxl

 

W2(X1,X2,X4,x5) = [ Vh + A51 — xl)](Alx2 — kxx5(U(x5)\/E ‘1' “("XS )‘IPS —' X4 ))

38
Since the equations for both of the derivatives are similar, the saturation limits will only

be formally developed for the derivative of state 3, '1’]. In order to define the limiting
pressures shown in Figure 4.1, the lower bound pressure tolerance, 53,1011 is defined and
then the other pressures are given by

£3 = 253.101

x, = Pe — 253,0, (4.1)

f3.101 = P5 T £3,101

The state is saturated by setting its derivative equal to zero. However, numerical

problems may arise in LSODE if the transition from the actual derivative value to zero is
too sudden, i.e., discontinuous. Therefore, a tolerance band was set up such that the
derivative ramps up from a value of zero to the actual function value at the lower bound
pressure. Similarly, the derivative ramps down from the function value at the upper
bound pressure to a value of zero. For this case, the saturation limits are also dependent
upon the nominal value of the derivative. When the pressure lies within the lower
tolerance band, its derivative should only be saturated if the derivative is negative.
Similarly, the derivative in the upper tolerance band should only be saturated if its value
is positive. That is, the saturation only occurs if the operating pressure is approaching a

pressure limit. Therefore, the corresponding saturated derivative is given by

 

 

f I- -‘
0: x3 < £3,101
‘1’ <0' x —x
1 - 3 -3 .
1+ ‘I’,(x,,x2,x3,x5). 53.1015 x3 553
_ 53,101 J
\PIJ'GI =< T _ T i (42)
x --x
1— 3 3‘Px,x,3r‘,x:x-Sx Sf
11;! >0: [ x 1( 1 2 3 5) 3 3 3.101
—3.tol
e _ O J‘3 > x3 rol

 

 

 

 

Since the saturation limits are enforced in the derivative of the state, they must also be
taken into account when considering the partial derivatives that define the Jacobian

matrix of the system. In particular,

39

 

 

 

 

 

 

 

 

f r- _\
0° x3 < £3,101
W, <0: 1+X3-53 ”1(xlex29539x5). < <
- £3.10! - x3 — £3
N! £310! ax
__;‘_‘1 =1 _-, ' ‘ g t (4 3)
1 ° $3,“)! axi
, _ 0: x3 > fBJoIJ
where i=1, 2, 5
and
- ~1
' 0' x3 < £3,101
‘1’ <0: 1
‘ ‘P(x,x,x,x): x Sx Ex
1 1 2 _3 5 —3.rol 3 _3
”I sat l. 13-"),
—'—— =1 r i > (4.4)
8x3 -1 _ _ _
111 > 0. ‘I’,(x,,x2,x3,x5): x3 S x3 5 x3301
I ' £3,101
1 e 0 x3 > x3 to”,

 

 

 

 

With these saturation limits in place, the FORTRAN code was ready to be linked with
LSODE so that the time response of the prefiltered closed loop system could be
simulated.

Although simulations were conducted for both cut-off frequencies of 100 Hz and
200 Hz, the same observations can be made from either set. Therefore, the results of the
control system whose filter’s cut-off frequency was 200 Hz will be discussed in detail
while the results for 100 Hz can be found in Appendix C. Although the control system
that was defined by the two-dimensional conic sector bound method failed to meet the
robust stability criteria in the frequency domain, it was still of interest to see how the
system performed in the time domain. Therefore, all figures will display three time
traces. The time history of the reference position will be shown by a dotted line, the
actuator position from the optimal four-dimensional conic sector bound method by a solid

line, and the actuator position from the two-dimensional conic sector bound method by a

dashed line.

40
Before simulating the control system, the operating limits of the actuator had to be

considered. It was assumed that the actuator would initially lie at its mid-stroke position
and the pressures would be set such that the actuator was in equilibrium. The three
operating limits identified arose from kinematic, force, and power limitations of the
actuator. Although the actuator has a stroke length of 2.1 inches, a kinematic limit of 1.8
inches was established in case the control system was too aggressive. In order to
determine the force and power limitations of the actuator, a representative from a
servovalve manufacturer was contacted. Assuming that the forcing was sinusoidal, the
proposed force limit of 400 lb and power limit of 1 Hp were found to be reasonable. A
set of simulations was then developed and is given in Table 4.1. The table provides the
forcing frequency and amplitude for each emulated mass and the resultant displacement
from the initial actuator position.

Figure 4.2 displays the time response of the system for an emulated mass of 50 lb
and a forcing frequency of 1 Hz. Although both controllers are tracking the reference
position, the controller designed from the optimal conic sector bounds is experiencing
“ringing.” The ringing or high frequency component of the actuator motion is a

phenomenon that has been noticed with the actual hydraulic system when the controller

Table 4.1 Outline of Time Domain Simulations

 

 

 

 

 

 

 

 

 

Emulated Mass of 50 lb Emulated Mass of 1501b
Frequency, Hz Force, lb Deflection, in Force, lb Deflection, in
1 4.6 0.9 13.8 0.9
10 328 0.64 400 0.261
25 400 0.125 400 0.042
50 400 0.031 400 0.010

 

 

 

 

 

WWW]

 

 

 

 

 

to I I I I I I I I I

e
e e e - u . . e o
. . . . . - . .
n . e . e u u , . e 1
‘Gb ------------- I .................................................... ' ................. .q
C e . e e e ' 1 e . e
I l l I l . I . I O
I . . , e u - e . .
h e e o e a e

d
.
.
.
.
O
O
C
.
C
.
I
.
C
C
.
.
.
.
.
I
.
.
O
.
.
.
O
.
.
U
U
C
l
.
.
.
.
.
C
.
.
.
.
.
U
.
.
D
.
.
O
.
I
O
.
.
.
.
.
.
.
.
.
O
.
.
.
.
.
O
O
I
.
.
I
.
I
U
.
C
D
l

I C I I I I I I '
O I I D D I Q I I
0.3 ...................... . .......... , ....... , ..... , ....................... , ........ ..
I O O ‘ I O O 9 .
O I O I O I I O C
C i I I e D e O O
I O I ‘

 

 

I e . . . . . . . .
0. - L 1 1 1 1 1 1 1 4 1

0.5 0.55 0.5 0.55 0.7 0.75 0.8 0.85 0.9 0.95 1
Time [e]

 

Figure 4.3. Time Trace with we = 200 Hz, M = 501b, f = 10 Hz

42
gains are set too aggressively. When the forcing frequency is increased to 10 Hz, the

ringing phenomenon no longer occurs with the control system defined by the optimal
conic sector bound method (Figure 4.3). As in the previous simulation, the controller
defined by the two-dimensional conic sector bound method performs better than the
optimal method. Figure 4.4 displays the results when the forcing frequency is equal to
25 Hz. Once again, the control system defined from the two-dimensional method is
tracking the reference signal better than the optimal method. The actuator position from
the optimal method is also displaying a high frequency component that may indicate the
system is undergoing a gentler ringing than the one shown in Figure 4.2. When the time
histories of a forcing frequency of 50 Hz are considered (Figure 4.5), it becomes apparent
that the controller developed from the two-dimensional method performs satisfactorily at
all frequencies, whereas the optimal method results in a controller that may be too
aggressive.

Since the control system was able to emulate a mass of 50 lb, the next step was to
investigate how well the controller performed when emulating a mass of 150 lb. Figure
4.6 displays the time traces when the forcing frequency equals 1 Hz. The results are very
similar to the smaller emulated mass and the control system from the optimal method is
experiencing the ringing phenomenon. Figure 4.7 displays the time histories when the
forcing frequency is increased to 10 Hz. Unlike the emulated mass of 50 lb, the control
system from the optimal method displays a high frequency ringing. The simulation
results for forcing frequencies of 25 Hz and 50 Hz are given in Figures 4.8 and 4.9. The
ringing phenomenon from the optimal method is becoming more pronounced and results
in a control system whose performance is not acceptable. Since the deflections from the
nominal actuator position decrease with increasing forcing frequency, the amplitude due

to the ringing becomes more significant.

43

 

 

a ”51111 q q 5
I O
O
i I
H e
I O
05.5" a
.5L -
... 555555 .
. 555 .
.1555
'5.
D '-
. V"
v ee
5.5555 ....
555
55
55
5 .
.
C
.5 .
5 . - ......................... u eeeeeeeeeeee O eeeeeeeeeeeeeeeeeeeeeeee 1%
5. . . o
n 11111.1
9'55 1
. . .555.
I 'n"
.'
.
. v.”
5 ...
. 55 0
. 555....V5
. 555
|
- “
. ‘ed:..ei..1
5
.‘.‘.
55
5
T ......... u ......... .555

 

I
'.
r
t r
I‘
l
0.85

.1
1
1
I.
0.9
Time [s]

5
5
5
. .55
5555

5
5
n 5.555 .
. . 555 .
. 5555 .o. u
55

 

 

1.2

Figure 4.4. Time Trace with we = 200 Hz, M = 501b, f = 25 Hz

 

 

 

 

 

 

A . .
0 e n
. 5’ 1.. . .
t.’ n .
e ‘ e
I eeeeeeeeeee W eeeeeeeeeeee u eeeeeeee + eje’Ikehe-eeeOeee-.eu eeeeeeeeeee
- .e1u' .
. . .J
u .1.)
‘2-
. ‘ eee
553-4....
5‘ eeeeee e .
555 ..... .
II N .
I u .
5". . .
‘0' . .
55..” 5. .
a- .
n .1 1 .
o i .
. .al\
. \..U.
‘ Wee.
‘Je‘eue-ee
5 4.4 ..... .
5 5 .
.6
. .0“
, I o
I‘ n ..
. . 'l
. Pee .
. . 5.55.5 .
m m 1.3“.
e e \e-e
. . ”5.5.. .
p L- k25efl.de.b..
2 A m m m
er. I
1 U 5

Time [e]

Figure 4.5. Time Trace with we = 200 Hz, M = 501b, f = 50 Hz

 

 

 

 

 

 

 

I
O I I O C O I I I
O O I I l I I I O
‘ .. ........ 1. ...... 1. ......... -. .- .......... ' ........ - .......... - ...... .1 ......... r ...... .1
I ‘ l O I I I U 0

Actual! Position [In]

. e e .
e e e e e I e
e e e n 'e u e e e
. . T. e e - e n e
09.1., ......... - ....... . .......... - ....... - .......... - ....... - .......... - ....... .; ........ .
. . e . e u e e -
e e g e u e e e e
e u . e n e e e
~ I . . . . -

 

 

Figure 4.7. Time Trace with we = 200 Hz, M = 150 lb, f= 10 Hz

45

 

 

 

 

 

 

 

Time [s]

 

 

 

 

 

 

 

 

 

 

f=25Hz

=200Hz,M=1501b,

Figure 4.8. Time Trace with we

 

 

 

 

 

 

 

 

 

 

0.95
Tlrne [e]

 

 

1.1

 

0.8

' Figure 4.9. Time Trace with we = 200 Hz, M = 150 lb, f = 50 Hz

46
Although the frequency domain results implied that the controller designed with

the optimal conic sector bound method would be superior to the two-dimensional method,
the time domain results show that the opposite is true. It is suspected that the peak
displayed in Figures 3.5 and 3.6 may be indicative of a resonance that occurs in the actual
system near 150 Hz. It is important to remember that the method of conic sector
bounding allows you to capture all of a function’s behavior with a set of linear equations,
but the nominal parameter values of one set of conic sector bounds may be superior to
another. In this case, it appears that the two-dimensional conic sector bound method
resulted in a transfer function that more closely predicts the nonlinear behavior of the
system. Since the control system defined by the optimal conic sector bound method does
not display a sudden peak in the structured singular value analysis, the controller may be
acting more aggressively at this frequency than the two-dimensional case. This

aggressiveness was shown in the time domain by the ringing phenomenon.

CHAPTER 5. CONCLUSIONS AND RECOMMENDATIONS

By representing the actual nonlinear system with an uncertain set of linear
systems, a plant description that was readily amenable to the application of established
robust control design methodologies was determined. This model was used to design a
controller that satisfied criteria for nominal performance and internal stability. By using
Robust Control Theory, the frequency domain analysis predicted that the controller
designed with the two-dimensional conic sector bound method would not perform as well
as the controller designed from the optimal four-dimensional conic sector bound method.
However, the time domain results found the opposite to be true. To understand this,
recall that Robust Control Theory provides guaranteed bounds on the worst case
performance of the control system. In this case, the time domain results showed that the
frequency domain predictions for robust stability and robust performance were
conservative. However, by using both methods to analyze the hydraulic system, a control
system design that allows the mechanical impedance of a one degree of freedom system
to be prescribed has been determined.

With the first phase of the project near completion, several areas can be suggested
for future work. First, the control system that has been designed will be implemented on
a test stand at the NASA Langley Research Center to compare the performance of the
actual system with the time domain simulations. Next, the results of the one degree
freedom system will be expanded to the six degree of freedom hydraulic system on
ARES 11. Therefore, the Rotorcraft Aeroelasticity Group will be able to emulate different

scale helicopter fuselage characteristics with a single hydraulic mount. Another area of

47

48

investigation may be to determine if it is possible to find a less conservative set of
frequency domain results. This would allow for better correlation between time domain

and frequency domain results.

APPENDIX A. MATLAB CODE

49

% Name: Jerome J. Palazzolo
% File: hssvpf.m
% Date: July 1, 1994

% This file determines the structured singular value of the

% M-Delta structure using the mu-synthesis toolbox in MatLab.
% Optimal conic sector bounds were provided by Dr. FitzSimons
clear

% Define the state space model of the nominal plant: [aap,bbp,ccp,ddp]

% Define nominal system parameters
mas = (1/386)* 100;

1h = 12;
ps = 3000;
plim = 300;

psuppc = ps-plim;
pslowc = plim;

disp(' ');
disp(['Pressure Range: (',num2str(pslowc),',',num2$tr(psuppc),')‘]);

% Define the intermediate constants

aal = 0.771;
a2 = 0.466;
aas = 0.0532;
bet = 250000;
k1 = 0.025;
k2 = 175;

kf = 93;

kw = 13.4;
kx = 85;

la = 2.1;

wn = 2*pi*730;
xi = 0.4;

% Define the equilibrium parameter values
c1 = la/2;

c3 = ps/(l+aa1/aa2);

c4 = c3*aa1/aa2;

% Define the average volume of fluid in the hydraulic line
% that connects the actuator to the servovalve

dh = 0.234;

vh = lh*pi*dh*dh/4;

v1= vh + aa1*c1;

v2 = vh + aa2*(la-c1);

% Calculate the constants that describe the dynamics of
% the servovalve.

a0 = k2*wn*wn*kw/(aas*kf);

a1 = wn*wn;

a2 = 2*xi*wn;

b0 = k2*wn*wn*kl/(aas*kf);

50

% Define the operating range of the hydraulic actuator

xmin = 0.05; %input('Enter the minimum actuator position [in]: ');
xmax = 2.05; %input('Enter the maximum actuator position [in]: ');
disp(' ');

disp(['Actuator Range: (',num23tr(xmin),',',num25tr(xmax),')']);

% Enter the conic sector bound results from the following
% files: csb2#3.m, csb2#4.m, where # is the case (a,b,c,...)
casenum = '3C';

csb3c3

csb3c4

% Define the coefficients of the state space matrices in terms
% of the system constants. This defines the nominal system.
a22=-0.40/mas; % damping based upon viscosity of hydraulic fluid
a23= aal/mas;

a24=~aa2lmas;

a31=(a3lmax+a31min)/2;

a32=(a32max+a32min)/2;

a33=(a33max+a33min)/2;

a35=(a35max+a35min)/2;

a4l=(a41max+a4lmin)/2;

a42=(a42max+a42min)/2;

a44=(a44max+a44min)/2;

a45=(a45max+a45min)/2;

a75=-a0;

a76=-a1;

a77=-a2;

b7 1=b0;

% Define the state space matrices of the nominal plant
aap = [ O, l, 0, 0, 0, O, 0;
0, a22, a23, a24, 0, 0, 0;
a31, a32, 333, 0, a35, 0, 0;
a41, a42, 0, a44, a45, 0, 0;
0, 0, 0, 0, 0, 1, O;
0, 0, 0, 0, 0, 0, l;
O, 0, 0, 0, a75, a76, a77 ];
bbp=[0;0;0;0;0;0;b71];
ccp=[1,0,0,0,0,0,0];
ddp = 0;
% Define the state space of the controller: [aac,bbc,ccc,ddc]

% Since the IMC controller is the inverse of the nominal
% plant augmented with filter dynamics, the first step will
% be to describe the transfer function of the nominal plant.
nlp = b71*(a23*a35+a24*a45);

n0p = -b7l*(a23*a35*a44+a24*a45*a33);

nump = [n1p, nOp ];

51

nomzer = roots(nump);
disp(' ');

% Check to see if the nominal system's zeros are all <= 0.
disp(['Nominal plant zero: ',num231r(nomzer)]);
if ( (nomzer > 0) l (nomzer == NaN) ),
disp(' ');disp('Sorry, nominal plant zeros can only be <= 0')
return
end

dconl= -(a22+a33+a44);

dcon2= -(a23*a32+a24*a42-a22*a33-a22*a44-a33*a44);
dcon3= a23*a32*a44+a24*a42*a33-a22*a33*a44-a23*a31-a24*a41 ;
dcon4= a23*a31*a44+a24*a41*a33;

d6p = dconl-a77;

d5p = dcon2-a77*dcon1-a76;

d4p = dcon3-a77*dcon2-a76*dcon1-a75;

d3p = dcon4-a77*dcon3-a76*dcon2-a75*dcon1;

d2p = -a77*dcon4-a76*dcon3-a75*dcon2;

dlp = -a76*dcon4-a75*dcon3;

dOp = -a75*dcon4;

denp = [ 1, d6p, dSp, d4p, d3p, d2p, dlp, dOp ];

% Define the transfer function of the controller in terms

% of the nominal plant and the filter dynamics.

wc = input('Enter the filter cut-off frequency [Hz]: ');

lam = l/(2*pi*wc);

numc = denp;

fltr = [lam"6,6*lam"5,15*lam"4,20*lam"3,15*lam"2,6*lam,0];
denc = conv(nump,fltr);

% Since the nominal plant and the controller have been
% defined, we can now check the internal stability of
% the system.
instchk = intstab(nump,denp,numc,denc);
if instchk = 0,
disp(' '); disp('Sorry, internal stability conditions violated');
return
end
if instchk == 1,
disp(' '); disp('Intemal stability conditions satisfied');
end
disp(' ');disp('Controller minreal');
[betc,alpc]=minreal(numc,denc);

% Use the transfer function of the controller to define the

% state space representation in controller canonical form.

% First, scale the matrices with tml to improve the conditioning
% for the balanced realization

[aac l,bbc1 ,ccc 1 ,ddc] = thss(betc,alpc);

tml = diag([le4, 1e7, 1e11, 1e14, 1e17, 1e20]);

aac = tml *aac1*inv(tm1);

bbc = tml*bbc1;

ccc = ccc 1 *inv(tml);

52

% Conduct a balanced realization to improve the numerical
% conditioning of the system.

atmpc = aac;

[ablac,bbalc,cbalc,ggc,ttc] = obalreal(atmpc,bbc,ccc);

aac = ablac;

bbc = bbalc;

ccc = cbalc;

% Define the prefilter to the system. The prefilter to the

% system not only allows you to prescribe the mechanical

% impedance of the system, but it also reduces the nominal

% amplitude reduction and phase lag of the closed loop system.
numpf = [15*lam"2, 6*lam, 1];

denpf= conv([l 0.1],[1 0.l]);

masd = (1/386)*input('Enter the emulated mass [lb]: ');

denpf( 1, l )=masd;

disp(' ');disp('Prefilter minreal');

[numpf,denpf] = minreal(numpf,denpf);

% Once again, conduct a balanced realization to improve
% the numerical conditioning of the system.

[aaf, bbf, ccf, ddf] = thss(numpf, denpf);
[aaf,bbf,ccf,ggf,ttf]= obalreal(aaf,bbf,ccf);

% Define the M-Delta uncertainty scaling matrices

perunc3 = (1/100)*input(‘Enter the uncertainty in x3dot [%]: ');
perunc4 = (1/100)*input('Enter the uncertainty in x4dot [%]: ');
a31d = abs(a31max-a31min)/2;

a32d = abs(a32max-a32min)/2;

a33d = abs(a33max-a33min)/2;

a35d = abs(a35max-a35min)/2;

a4ld = abs(a41max-a41min)/2;

a42d = abs(a42max-a42min)/2;

a44d = abs(a44max-a44min)/2;

a45d = abs(a45max-a45min)/2;

a31del = sqrt(perunc3*a31d);

a32del = sqrt(perunc3*a32d);

a33del = sqrt(perunc3*a33d);

a35del = sqrt(perunc3*a35d);

a41del = sqrt(perunc4*a4ld);

a42del = sqrt(perunc4*a42d);

a44del = sqrt(perunc4*a44d);

a45del = sqrt(perunc4*a45d);

bmu = [ zeros(2,8);
a3 ldel, a32del, a33del, a35del, O, 0, 0, O;
0, 0, 0, 0, a41del, a42del, a44del, a45del;
zeros(3,8) ];

cmu = [ a31del, zeros(1,6);
0, a32del, zeros(l,5);
0, 0, a33del, zeros(l,4);
zeros(l,4), a35del, 0, 0;

53

a41del, zeros(l,6);

0, a42del, zeros(l,5);
zeros(1,3), a44del, zeros(1,3);
zeros(l,4), a45del, 0, 0 ];

% Matrices defined to this point

% Nominal plant: [aap,bbp,ccp,ddp]

% Controller: [aac,bbc,ccc,ddc]

% Interconnections: (bmu, cmu)

% We are now ready to define the M-Delta structure.

% First consider the open loop system whose inputs are MUin and
% the controller output, and whose outputs are MUout and the

% actuator position.

tm2 = diag([1,1,1e-3,1e-3,1e3,l,le-3]);

aam = tm2*aap*inv(tm2);

bbm = tm2*[bmu, bbp];

ccm = [cmu; ccp]*inv(tm2);

ddm = zeros(5,5);

% Conduct a minimal realization of the open loop system.
disp(' ');disp('Open Loop (OL) minreal');
[aam,bbm,ccm,ddm] = minreal(aam,bbm,ccm,ddm);

% Conduct a balanced realization of the open loop system

% to improve the numerical conditioning.

shft = 0.01;

disp(' ');disp(['Amount of A-shifting in OL(s): ',num2str(shft)]);
atmp = aam - shft*eye(max(size(aam)));

[abla,bbal,cbal,gg,tt] = obalreal(atmp,bbm,ccm);

abal = abla + shft*eye(max(size(aam)));

% Extract the new [aap,bbp,ccp,ddp] and (bmu,cmu).
aap = abal;

bmu = bbal( 1 :max(size(aap)), 1 :8);

bbp = bba1( 1 :max(size(aap)),9);

cmu = cbal( 1 :8,l:max(size(aap)));

ccp = cbal(9,1:max(size(aap)));

% Define M-structure.

szp = max(size(aap));

szc = max(size(aac));

szf = max(size(aafl);

amtotl = [ aap-bbp*ddc*ccp, bbp*ccc, bbp*ddc*ccf;...
-bbc*ccp, aac, bbc*ccf;
zeros(szf,szp), zeros(szf,szc), aaf];

bmtotl = [bmu, bbp*ddc*ddf;
zeros(szc,8), bbc*ddf;...
zeros(szf,8), bbt];

cmtotl = (emu, zeros(8,(szc+szf)); -ccp, zeros(1,szc), ccf];

dmtot = zeros(9,9); dmtot(9,9)=ddf;

% Conduct minimal and balanced realizations of the M—structure.
disp(' '); disp('M-structure minreal');

54

[amtot,bmtot,cmtot,dmtot]=minreal(amtot l ,bmtot 1 ,cmtot l ,dmtot);
conamt=cond(amtot);

shft2 = 0;

disp(' ');disp(['Amount of A-shifting in M(s): ',num23tr(shft2)]);
atmp2 = amtot - shft2*eye(max(size(amtot)));
[amtot,bmtot,cmtot,ggt,ttt]=obalreal(atmp2,bmtot,cmtot);

amtot = amtot + shft2*eye(max(size(amtot)));

diSP(' ');

disp('cond(amtota)=');

disp(conamt);

disp('cond(amtot)=');

disp(cond(amtot));

disp('eig(amtot)=‘);

disp(eig(amtot));

disp('Stability alone: complex uncertainty - (8x8)');

disp('Performance: complex uncertainty, complex performance - (9x9)');
test = input('Robust Stability (0) / Robust Performance (1): ');

% Use the SYSTEM data structure and convert it to
% the VARYING format at discrete frequency points
sysmat = pck(amtot,bmtot,cmtot,dmtot);

disp(' ');minfo(sysmat)

wrad = logspace(-2, 4, 500);

sysfrsp = frsp(sysmat, wrad);

% Define the Delta structure diag(Delta_U, Delta_P)
% Real uncertainty: [~ones(8,l), ones(8,l)]
% Complex uncertainty: [ones(8,l), ones(8,l)]
% Performance uncertainty is always complex.
if (test == 0),

stadel = [ones(8,l),ones(8,l)];

else

stadel = [ones(8,l),ones(8,l)];

end

perdel=[1,1];

sysdel = [stadel; perdel];

if (test == 0),

% Consider Robust Stability

m1 lfrsp = sel(sysfrsp, [1:8], [1:8]);

[rsbnds,rsdvec,rssens,rspvec,rsgvec] = mu(m11frsp, stadel);

vplot('liv,m', rsbnds);

title(['Case ',casenum,': Wc=',...
num25tr(wc),'[Hz]; Uncertainty - x3=',...
num25tr(100*perunc3),'%, x4=',...
num25tr(100*perunc4),'%; Mass — Actual=', num2str(386*mas),...
'[lb], Emulated=', num23tr(386*masd), '[lb]']);

xlabel('Frequency [rad/sec]‘);

ylabel('Robust Stability SSV');

grid;

maxssv=max(rsbnds( 1 :max(size(wrad)), 1 :2))

55

else

% Consider Robust Performance

[rpbnds,rpdvec,rpsens,rppvec,rpgvec] = mu(sysfrsp, sysdel);

vplot('liv,m', rpbnds);

title([‘Case ',casenum,': Wc=',...
num25tr(wc),'[Hz]; Uncertainty - x3=',...
num28tr(100*perunc3),'%, x4=',...
num2str(100*perunc4),'%; Mass - Actual=', num2str(386*mas),...
'[lb], Emulated=', num2str(386*masd), '[lb]']);

xlabel('Frequency [rad/secl');

ylabel('Robust Performance SSV');

grid;

maxssv=max(rpbnds( 1 :max(size(wrad)), 1 :2))

end;

56

% Name: Jerome J. Palazzolo
% File: csb3c3.m
% Date: July 1, 1994

% Conic Sector Bound results from optimization

disp(' ')
disp('Case 3C Conic Sector Bounds, x3dot')

%a31max= 2.7460e6; needed to make a31max=-a31min

a31max= 7.7540e6;

a31min= —7.7540e6;

a32max= (1/40)*(-5.8 l62e6);
a32min= (1/40)*(-l.5610e7);
a33max= (1/1000)*(0);

a33min= (l/lOOO)*(-6.7498e6);
a35max= (1/0.0175)*(2.1093e7);
a35min= ( l/0.0175)*(7.3144e6);

% Name: Jerome J. Palazzolo
% File: csb3c4.m
% Date: July 1, 1994

% Conic Sector Bound results from optimization

disp(' ')
disp('Case 3C Conic Sector Bounds, x4dot')

%a41max= 3.3936e6; needed to make a41max=-a41min

a41max= 8.1215e6;

a41min= -8.1215e6;

a42max= (1/40)*(9.4970e6);
a42min= (1/40)*(3.7817e6);
a44max= (1/1000)*(5.8900e2);
a44min= (1/1000)*(-6.3890e6);
a45max= (1/0.0175)*(-1.1977e7);
a45min= (1/0.0175)*(-2.2187e7);

57

% Name: Jerome J. Palazzolo
% File: intstab.m
% Date: July 1, 1994

% This script tests the nominal plant and controller for
% internal stability.

function intstl = intstab(nump,denp,numc,denc);

% Pre-load the numerator of the plant to match dimensions
nump = [zeros( 1,6),nump];

% Define the transfer characteristics of the transfer
% matrix that defines internal stability

% [1/(1+pc)]* [1pc p 1 =[[a b]

% I 0 -pc 11 [c -a l

numa = conv(nump,numc);

numb = conv(nump,denc);

numc = conv(numc,denp);

den = numa + conv(denp,denc);

% Conduct minimal realizations to cancel out common poles
% and zeros of the transfer functions.

[num 1 ,den 1 ]=minreal(numa,den);
[num2,den2]=minreal(numb,den);
[num3,den3]=minreal(numc,den);

rts=[roots(denl); roots(den2); roots(den3)];

% The system is internally stable if all of the poles
% are in the open left-half plane (LHP, eig's < 0)
chk = max(real(rts));

if (chk < 0)
% Stable
intstl = 1;

else
% Unstable
intstl = 0;

end

return

APPENDIX B. FORTRAN CODE

58

c Name: Jerome J. Palazzolo
c File: hpfsim.f

c Date: July 1, 1994

c

c This fortran program will be linked to LSODE in order to simulate the H2
c optimal controller for a hydraulic actuator. In this case, the derivatives

c of the pressures saturate to guarantee { (ptol) <= Pi <= (Ps-ptol) }. These

c saturation limits are enforced with the functions: ydot3, ydot4, pd3i, pd4i.
c The input signal will be prefiltered in order to simulate different masses.

c...MAIN PROGRAM of Hydraulic Servo-actuator Simulation
c ........ Size of rwork (mf=21): 22+9*neq+neq"2
c ........ Size of iwork (mf=21): 20+neq

c.Variable type declaration
implicit none
external fex, jex, ssform
real*8 atol, rtol, rwork(500), t, tout, y(l6), dt, tf
real*8 sscon(14), wc, masl, masZ
integer i, iout, iwork(50), itol, itask, iopt, lrw
integer liw, mf, neq, istate, i1
character* 10 fname

c.Common variable type declaration
real*8 a22,a23,a24,a31,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aa] ,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b7 l ,
+ bet,cc1 1,cc12,cc13,cc14,cc15,cc16,cc17,curi,dc1 1,dh,frqr,
+ kl ,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

common a22,a23,a24,a31,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aa 1 ,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b7 l ,
+ bet,cc1 1,cc12,cc13,ccl4,cc15,ccl6,cc17,curi,dc1 1,dh,frqr,
+ k1,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

c.Open data file to record the integration results
write (*,*) 'Input the file name for the Matlab data file: '
read (*,2) fname

2 format(A10)
open (unit=1, file=fname, status='new')

c.Define the number of equations, i.e., the number of states

c ................. y( l) = actuator position
c ................. y(2) = actuator velocity
c ................. y(3) = pressure 1

c ................. y(4) = pressure 2

c ................. y(5) = spool position

c ................. y(6) = spool velocity

c ................. y(7) = spool acceleration

c ......... y(8) to y( 14) = controller states
c ....... y( 15) to y(16) = prefilter states
neq = 16

59

c.Define the user inputted parameters
masl = 100.d0
mas = masl/386.d0
write(*,*) 'Enter emulated mass [lbm]: '
read (*,*) masZ
masd = mas2/386.d0
1h = 12.d0
ps = 3000.d0
ptol = 50.d0
write (*,*) 'Enter forcing signal amplitude [lbf]: '
read (*,*) ampr
write (*,*) 'Enter forcing signal frequency [Hz]: '
read (*,*) frqr
write (*,*) 'Enter the filter cut-off frequency [Hz]: '
read (*,*) wc

c.Define the intermediate constants
aal = 0.771d0
aa2 = 0.466d0
aas = 0.0532d0
bet = 250000.d0
k1 = 0.025d0
k2 = 175.d0
kf = 93.d0
kw = 13.4d0
kx = 85.d0
1a = 2.1d0
pi = datan(l.d0)*4.d0
pslowc = 0.d0
psuppc = P5 ,
wn = 2.d0*p1*730.d0
xi = 0.4d0

c.Define the average line volume.
db = 0.234d0
vh = lh*pi*dh*dh/4.d0

c.Calculate the linear coefficients of the state space matrices for the system.

c.a3* and a4* are the results of the conic sector bound optimization.
a22=—0.40d0/mas
a23= aal/mas
a24=-aa2/mas
a31= (7.7540d6-7.7540d6)/2.d0
a32= (1 .d0/40.d0)*(-5.8 l62d6-1.5610d7)/2.d0
a33= (1.d0/1000.d0)*(0.d0-6.7498d6)/2.d0
a35= (1.d0/0.0175d0)*(2.1093d7+7.3144d6)/2.d0
a41= (8.1215d6-8.1215d6)/2.d0
a42= (l.d0/40.d0)*(9.4970d6+3.78 l7d6)/2.d0
a44= ( l .dO/1000.d0)*(5.8900d2-6.3890d6)/2.d0
a45= (l.d0/0.0175d0)*(-1.1977d7-2.2187d7)/2.d0
a75=-k2*wn*wn*kw/(aas*kf)
a76=-wn*wn
a77=-2.d0*xi*wn

60
b7l= k2*wn*wn*kl/(aas*kf)

c.Define the state space representation constants (non-zero/non-unity) for
c.the nominal controller, C, with the given lambda.

lam = 1.d0/(2.d0*pi*wc)

call ssform(sscon)

ac72=sscon( 1)

ac73=sscon(2)

ac74=sscon(3)

ac75=sscon(4)

ac76=sscon(5)

ac77=sscon(6)

ccl 1=sscon(7)

cc 12=sscon(8)

cc 13=sscon(9)

cc 14=sscon( 10)

cc 15=sscon( 1 1)

cc 16=sscon( 12)

cc17=sscon(13)

dc 1 1=sscon( 14)

c.Zero out the y-array
do 3 i1 = 1,neq
y(il)=0.d0
3 continue

c.Define the initial values for the system
y(1)=1a/2.d0
y(3)=ps/(1.d0+aa1/aa2)
y<4>=(aa1/aa2)*y(3)
y( l 6) = -ampr/(masd*2.d0*pi*frqr)

c.Define the integrator time controls
t=0.d0
dt=0.001d0
tf=1.d0
write(*,*) ' '
write(*,*) 'Initial simulation settings have been set.’
write(*,*) ' '

c.Echo initial settings to the data file
write (1,4) 'Initial actuator position: ', y(l)
write (1,4) 'Initial actuator velocity: ', y(2)
write (1,4) 'Initial pressure on side 1: ', y(3)
write ( 1,4) 'Initial pressure on side 2: ', y(4)
write (1,4) 'Initial spool position: ', y(5)
write (1,4) 'Initial spool velocity: ', y(6)
write (1,4) 'Initial spool acceleration: ', y(7)
write (1,4) 'Forcing signal amplitude: ', ampr
write (1,4) 'Forcing signal frequency: ', frqr
write (1,4) 'Filter cut-off frequency: ', we
write (1,4) 'Initial time: ', t
write (1,4) 'Integrator time step: ', dt
write (1,4) 'Final time: ', tf

61
4 format('% ', A30, e125)

c.Write the beginning array specification for the data file.
write (l,*) ' '
write (1,*) 'amp =' , ampr, ';'
write (l,*) 'bet =' , bet, ';'.
write (1, *) 'frq— -' ,,';frqr
write (1, *) 'lh= -' ,lh , ';'
write (1, *)' mas— -' ,masl, ';'
write (l,*) 'masd=' , mas2, ';'
write (l,*) 'wc =' , wc , ';'
write (l,*) 'ps =' , ps , ';'
write (1,*) 'ptol=' , ptol
write (l,*) 'timhis = [ '

c.Set the lsode tolerance parameters for the state variables
itol = 1
rtol = 1.d-9
atol = 1.d-9

c.Set lsode parameters
itask = l
istate = l
iopt = 1
do 5 i=5,10
iwork(i)=0
rwork(i)=0.d0
5 confinue
iwork(6)=1000
rwork(6)=0.001d0

lrw = 500
liw = 50
mf = 21
c.Call LSODE to solve the Dif EQ set
iout- — 0
write(l 20) 0 d0.y(1).y(2).)’(3) y(4).Y(5) Y(6) 3(0)
+ 0. d0, 1 0.5d0,1.005d

10 continue
iout = iout + 1
tout = dfloat(iout)*dt
call lsode(fex,neq,y,t,tout,itol,rtol,atol,itask,istate,
+ iopt,rwork,lrw,iwork,liw,jex,mf)
if (istate .lt. 0) then
write( 1.25) t,)'(1).Y(2).Y(3).Y(4),Y(5).Y(6).y(7),curi,
+ ref,refc
write(l,*) ' '
write(l,*) 'istate = ',istate
goto 15
endif
if (t .lt. tf) then
write(1.20) t.y(1).y(2).y(3>.y(4).y(5).y<6).y(7).curi.
+ ref,refc

62

else
write(l.25) t.y(1).Y(2).Y(3),Y(4),y(5).Y(6).y(7).curi.
+ ref,refc
goto 15
endif
goto 10

15 continue
20 format(el 1.3, 10(',' ,e12.5), ';')
25 format(el 1.3, 10(',' ,e12.5), '];')

if (istate .lt. 0) then
write(6,90) istate
endif
9O format(l/l22h error halt.. istate =,i3)
close(l)
end

 

szmm=-___-—======——_.=============: == =============

c...Subroutine fex - finds the derivatives of the state vector
subroutine fex(neq, t, y, ydot)

c.Variable declaration
real*8 t, y( 16), ydot(l6)
external ydot3, ydot4
real*8 ydot3, ydot4
integer neq

c.Declare parameters as double precision
real*8 v1,v2,conl,force

c.Common variable type declaration
real*8 a22,a23,a24,a3 l,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aa] ,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b7 1 ,
+ bet,ccl 1,cc l 2,cc l 3,cc l4,cc 15,cc16,cc17,curi,dc1 1 ,dh,frqr,
+ kl ,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

common a22,a23,a24,a3 1,a32,a33,a35,a41 ,a42,a44,a45,a75,a76,
+ a77,aa1,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,arnpr,b71,
+ bet,cc1 1,cc12,cc13,cc14,cc 15,cc 16,cc17,curi,dcl 1,dh,frqr,
+ k1,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

c.Determine the line volumes
v1 = vh+aal*y(l)
v2 = vh+aa2*(la-y(1))

c.Define the current reference signal
force = ampr*dsin(2.d0*pi*frqr*t)
ref=la/2.d0 + y(15)
refc=la/2.d0 + y(15) + 6.d0*lam*y( 16) +

63
+ 15.d0*lam*lam*force/masd

c.Define the control action (current, i) including saturation.
curi = ccll*y(8)+ccl2*y(9)+cc13*y(10)+ccl4*y(11)+cc15*y(12)
curi = curi+cc 16*y(l3)+cc l7*y( 14)+dcl 1*(refc-y( 1))
if (curi .lt. -10.d0) then
curi=-10.d0
else
if (curi .gt. 10.d0) then
curi=10.d0
else
curi=curi
endif
endif

c.Calculate the derivatives of the state vector
ydot(l) = y(2)
ydot(2) = a22*y(2)+a23*y(3)+a24*y(4)+force/mas

c. Derivatives of y3 and y4 are saturated. Call functions ydot3 and ydot4.
ydot(3) = ydot3(>’(2). y(3). y(5), v1)
ydot(4) = ydot4(y(2), y(4), y(5), v2)

ydot(S) = Y(6)
ydot(6) = y(7)
ydot(7) = a75*y(5)+a76*y(6)+a77*y(7)+b71*curi

ydot(S) = y(9)

ydot(9) = y(10)

ydot(10)= y(l l)

ydot(11)= y(l2)

ydot( 12) = y(13)

ydot(l3) = y(14)

conl = ac72*y(9)+ac73*y(10)+ac74*y(1 1)+ac75*y(12)

ydot( 14) = con 1+ac76*y( 13)+ac77*y( 14)+(refc-y( l ))

ydot(15) = y(16)
ydot( 16) = force/masd

c.Return to the call statement from LSODE
return
end

C:==——=:$=
c...Subroutine jex - defines the Jacobian matrix

subroutine jex(neq, t, y, ml, mu, pd, nrpd)
c.Variable declaration

real*8 t, y(l6), pd(16,16)

integer nrpd, ml, mu

c.Declare intermediate constants and elements as double precision

64

external pd31, pd32, pd33, pd35, pd41, pd42, pd44, pd45
real*8 pd3 1, pd32, pd33, pd35, pd41, pd42, pd44, pd45
real*8 v1,v2

integer i1,i2

c.Common variable type declaration
real*8 a22,a23,a24,a31,a32,aB3,a35,a41,a42,a44,a45,a75,a76,
+ a77,aa l ,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b7 l ,
+ bet,ccl 1,cc12,ccl3,cc14,cc15,cc16,ccl7,curi,dcl 1,dh,frqr,
+ kl ,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

common a22,a23,a24,a31,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aa],aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b7l,
+ bet,cc1 l,cc12,cc13,ccl4,cc15,cc16,ccl7,curi,dcl 1,dh,frqr,
+ kl ,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

c.Define the line volumes
v1= vh+aa1*y(l)
v2 = vh+aa2*(la-y(1))

c.Zero out the J acobian matrix.
do 500i1=1, neq
do 490 i2 = l, neq
pd(i1,12)=0.d0
490 continue
500 continue

c.Define the elements of the Jacobian (partial derivatives) that do not depend

c.upon the sign of y(5).
pd(l,2) = 1.d0
pd(2,2) = a22
pd(2,3) = a23
pd(2,4) = a24
pd(5,6) = 1.d0
pd(6,7) = 1.d0
pd(7,1) = -b7l*dcll
pd(7,5) = a75
pd(7,6) = a76
pd(7,7) = a77
pd(7,8) = b71*ccl 1
pd(7,9) = b71*ccl2
pd(7,10) = b71*cc13
pd(7,11) = b71*cc14
pd(7,12) = b71*cc15
pd(7,13) = b71*ccl6
pd(7,14) = b71*ccl7
pd(7,15) = b71*dc11
pd(7,16) = 6.d0*b71*dc1 1*lam
pd(8,9) = 1.d0
pd(9,10) = 1.d0
pd(10,11)= 1.d0
pd(11,12)=1.d0

65

pd(12,13)= 1.d0
pd(13,14)= 1.d0
pd(14,1) = -l.d0
pd(14,9) = ac72
pd(14,10)= ac73
pd(14,11)= ac74
pd(14,12)= ac75
pd(14,13)= ac76
pd(14,14)= ac77
pd(14,15)= 1.d0
pd(14,16)= 6.d0*lam
pd(15,16)= 1.d0

c.Define the elements of the Jacobian (partial derivatives) that depend
c.upon the saturated derivatives, namely pd3], pd32, pd33, pd35, pd41,
c.pd42, pd44, and pd45.

Pd(3.l)=Pd31(Y(2). y(3). y(5). vl)

Pd(3.2)=Pd32(y(2). y(3). y(5). v1)

Pd(3.3)=Pd33(Y(2). y(3). y(5). V 1)

pd(3.5)=pd35(y(2). y(3). y(5). v1)

Pd(4.1)=pd41(y(2). y(4). Y(5). v2)

pd(4,2)=pd42(y(2), y(4). Y(5). v2)

Pd(4.4)=pd44()’(2). y(4). y(5). v2)

pd(4.5)=pd45(y(2). y(4), y(5). v2)

c.Retum to the call statement from LSODE
return
end

 

C======m——__=::==__—-===__—-===__-—===__——========== : == === =

c...Subroutine ssform(sscon) - defines the nominal controller, c.
subroutine ssform(sscon)

c.Variable declaration
real*8 sscon(l4)

c.Declare intermediate constants and elements as double precision
real*8 n0p,nlp,dcon1,dcon2,dcon3,dcon4,d0p,dlp,d2p,d3p,d4p
real*8 d5p,d6p,fac,dc l ,dc2,dc3,dc4,d05,dc6

c.Common variable type declaration
real*8 a22,a23,a24,a31,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aa1,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,cc1 1,cc 12,cc 13,cc l4,cc 15,cc16,cc l7,curi,dcl 1,dh,frqr,
+ k1,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

common a22,a23,aZ4,a31,a32,a33,a35,a4],a42,a44,a45,a75,a76,
a77,aa],aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
bet,ccl 1,cc12,ccl3,cc14,cc15,ccl6,ccl7,curi,dc1 1,dh,frqr,
kl,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
ptol,ref,refc,vh,wn,xi

++++

66

c.Define the transfer function of the nominal plant
nlp = b71*(a23*a35+a24*a45)
nOp = -b71*(a23*a35*a44+a24*a45*a33)
dcon1= -(a22+a33+a44)
dcon2= -(a23*a32-I-a24*a42-a22*a33-a22*a44-a33*a44)
dcon3= a23*a32*a44+a24*a42*a33-a22*a33*a44-a23*a31-a24*a41
dcon4= a23*a3 1*a44+a24*a41*a33
d6p = dconl-a77
d5p = dcon2-a77*dcon 1 -a76
d4p = dcon3-a77*dcon2—a76*dcon1-a75
d3p = dcon4-a77*dcon3-a76*dcon2-a75*dcon1
d2p = -a77*dcon4-a76*dcon3-a75*dcon2
dlp = ~a76*dcon4-a75*dcon3
dOp = -a75*dcon4

c.Define the controller constants
fac = l.d0/(n1p*lam**6)
dc6 = fac*(6.d0*nlp*lam**5 + n0p*lam**6)
dc5 = fac*(15.d0*n1p*lam**4 + 6.d0*n0p*lam**5)
dc4 = fac*(20.d0*nlp*lam**3 + 15.d0*n0p*lam**4)
dc3 = fac*(15.d0*nlp*lam**2 + 20.d0*n0p*lam**3)
dc2 = fac*(6.d0*nlp*lam + 15.d0*n0p*lam**2)
del = fac*6.d0*n0p*lam

c.Define the coefficients for the state space form of C.
sscon( 1) - -dc1

sscon(2) = -dc2
sscon(3) = -dc3
sscon(4) = -dc4
sscon(5) = -dc5

sscon(6) — -dc6

sscon(7) = fac*dOp
sscon(8) = fac*(dlp-dcl)
sscon(9) = fac*(d2p-dc2)
sscon(IO) = fac*(d3p-dc3)
sscon(l 1) = fac*(d4p-dc4)
sscon(12) = fac*(d5p-dc5)
sscon(13) = fac*(d6p-dc6)
sscon(14) = fac

c.Retum to the main program

return
end

C======================================= __
c...Function ydot3 - determines the saturated derivative at state 3.
function ydot3(x2, x3, x5, voll)

c.Variable declaration
real*8 ydot3, x2, x3, x5, voll

c.Declare intermediate constants

67

real*8 dptol, yd3nom, yd310w, yd3upp
real*8 p11, p12, pul, pu2

c.Common variable type declaration
real*8 a22,a23,a24,a31,a32,a33,a35,a4],a42,a44,a45,a75,a76,
+ a77,aa1,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b7l,
+ bet,cc1 1,cc12,ccl3,cc14,cc15,ccl6,cc l7,curi,dcl 1,dh,frqr,
+ k 1 ,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

common a22,a23,a24,a3 1,a32,a33,a35,a4],a42,a44,a45,a75,a76,
+ a77,aal ,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,cc1 1,cc12,ccl3,ccl4,cc15,ccl6,cc17,curi,dcl 1,dh,frqr,
+ kl,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,

+ ptol,ref,refc,vh,wn,xi

c.Define intermediate constants
dptol = ptol/2.d0
p11 = dptol
p12 = ptol
pul = ps-ptol
pu2 = ps-dptol

c.Evaluate the saturated derivative.

if (x5 .gt. 0.d0) then

yd3nom = bet*(-aal *x2+kx*x5*dsqrt(dabs(ps-x3)))/voll
yd310w = bet*(-aa1*x2+kx*x5*dsqrt(ps-ptol))/vol1

yd3upp = bet*(-aa1*x2+kx*x5*dsqrt(ptol))/voll

else
yd3nom = bet*(-aa1*x2+kx*x5*dsqrt(dabs(x3)))/vol1
yd310w = bet*(-aa1 *x2+kx*x5*dsqrt(ptol))/vol1
yd3upp = bet*(-aal *x2+kx*x5*dsqrt(ps-ptol))/vol1

endif

ydot3 = yd3nom

if ( (x3 .lt. p12) .and. (yd3nom .lt. 0.d0) ) then
if (x3 .lt. p11) then

ydot3 = 0.d0
else
ydot3 = (l.d0+(x3-pl2)/dptol)*yd3low
endif
endif

if( (x3 .gt. pul) .and. (yd3nom .gt. 0.d0) ) then
if (x3 .gt. pu2) then
ydot3 = 0.d0
else
ydot3 = ( 1.d0-(x3-pul)/dptol)*yd3upp
endif
endif

c.Retum to the function call
return

68

end

 

C==================$======m========::

 

c...Function ydot4 - determines the saturated derivative at state 4.
function ydot4(x2, x4, x5, v012)

c.Variable declaration
real*8 ydot4, x2, x4, x5, v012

c.Declare intermediate constants
real*8 dptol, yd4nom, yd4low, yd4upp
real*8 p11, p12, pul, pu2

c.Common variable type declaration
real*8 a22,a23,a24,a3 1,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aa1,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,cc1 1,cc12,ccl3,cc14,cc15,cc16,cc17,curi,dc1 1,dh,frqr,
+ k1,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

common a22,a23,a24,a3 1,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aal,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,cc1 1,cc12,cc13 ,cc l4,cc 15,cc16,cc l7,curi,dc 1 1,dh,frqr,
+ kl ,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

c.Define intermediate constants
dptol = ptol/2.d0
p11 = dptol
p12 = ptol
pul = ps-ptol
pu2 = ps—dptol

c.Evaluate the saturated derivative.

if (x5 .gt. 0.d0) then
yd4nom = bet*(aa2*x2-kx*x5*dsqrt(dabs(x4)))/v012
yd4low = bet*(aa2*x2-kx*x5*dsqrt(ptol))/v012
yd4upp = bet*(aa2*x2-kx*x5*dsqrt(ps-ptol))/v012

else
yd4nom = bet*(aa2*x2-kx*x5*dsqrt(dabs(ps-x4)))/v012
yd4low = bet*(aa2*x2-kx*x5*dsqrt(ps-ptol))/vol2
yd4upp = bet*(aa2*x2-kx*x5*dsqrt(ptol))/v012

endif

ydot4 = yd4nom

if ( (x4 .lt. p12) .and. (yd4nom .lt. 0.d0) ) then
if (x4 .lt. pll) then
ydot4 = 0.d0
else
ydot4 = (1.d0+(x4-p12)/dptol)*yd4low
endif

69
endif

if ( (x4 .gt. pul) .and. (yd4nom .gt. 0.d0) ) then
if (x4 .gt. pu2) then
ydot4 = 0.d0
else
ydot4 = ( 1.d0-(x4-pu1)/dptol)*yd4upp
endif
endif

c.Retum to the function call
return
end

 

C ‘1'.

 

c...Function pd3] - determines the saturated partial derivative pd3].
function pd31(x2, x3, x5, voll)

c.Variable declaration
real*8 pd3], x2, x3, x5, voll

c.Declare intermediate constants
real*8 pd31nom, pd3llow, pd31upp, conl
real*8 p11, p12, pul, pu2, dptol, yd3nom

c.Common variable type declaration
real*8 a22,a23,a24,a3 1,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aa1,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,ccl 1,cc l2,cc13,ccl4,cc15,cc16,cc17,curi,dc1 1,dh,frqr,
+ k1,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

common a22,a23,a24,a3 l ,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aal,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,ccl 1,cc12,cc13,ccl4,cc15,cc16,ccl7,curi,dc1 1,dh,frqr,
+ k 1 ,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

c.Define intermediate constants
dptol = ptol/2.d0
p11 = dptol
p12 = ptol
pul = ps-ptol
pu2 = ps-dptol

c.Evaluate the nominal first derivative
if (x5 .gt. 0.d0) then
yd3nom = bet*(-aal *x2+kx*x5*dsqrt(dabs(ps-x3)))/vol1
else
yd3nom = bet*(-aa1*x2+kx*x5*dsqrt(dabs(x3)))/vol1
endif

70

c.Evaluate the saturated partial derivatives
conl = bet*aa1/(vol 1 *vol 1)
if (x5 . gt. 0.d0) then
pd3 1nom=conl*(aa1*x2-kx*x5*dsqrt(dabs(ps-x3)))
pd3 llow=con1*(aa1*x2-kx*x5 *dsqrt(ps-ptol))
pd31upp=con1*(aa1*x2-kx*x5*dsqrt(ptol))
else
if (x5 .lt. 0.d0) then
pd3 lnom=conl *(aal *x2-kx*x5*dsqrt(dabs(x3)))
pd3llow=con1*(aa1*x2-kx*x5*dsqrt(ptol))
pd3 lupp=con l *(aal *x2-kx*x5 *dsqrt(ps-ptol))
else
pd3 lnom=con1*aal *x2
pd3 llow=pd3 lnom
pd31upp=pd3 lnom
endif
endif

pd31 = pd31nom

if ( (x3 .lt. p12) .and. (yd3nom .lt. 0.d0) ) then
if (x3 .lt. p11) then

pd31 = 0.d0
else
pd31 = (1.d0+(x3-pl2)/dptol)*pd3 llow
endif
endif

if ( (x3 .gt. pul) .and. (yd3nom .gt. O.d0) ) then
if (x3 .gt. pu2) then
pd31 = 0.d0
else
pd31 = (1.d0-(x3-pu 1)/dptol)*pd3 lupp
endif
endif

c.Retum to the function call
return
end

 

C:==-_—_-==-——_—_=== I--"-‘ _ _--—---- _-

 

 

nun-i — ---——w—‘ —-

c...Function pd32 - determines the saturated partial derivative pd32.
function pd32(x2, x3, x5, voll)

c.Variable declaration
real*8 pd32, x2, x3, x5, voll

c.Declare intermediate constants
real*8 pd32nom, pd3210w, pd32upp
real*8 p11, p12, pul, pu2, dptol, yd3nom

c.Common variable type declaration

71

real*8 a22,a23,a24,a31,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aa] ,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b7 1,
+ bet,cc1 l,cc12,ccl3,ccl4,cc15,cc16,ccl7,curi,dcl 1,dh,frqr,
+ k1,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

common a22,a23,a24,a31,a32,a33,a35,a4l,a42,a44,a45,a75,a76,
+ a77,aa1,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,cc1 1,cc12,cc13,ccl4,cc15,cc16,cc l7,curi,dc 1 1,dh,frqr,
+ k1,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

c.Define intermediate constants
dptol = ptol/2.d0
pll = dptol
p12 = ptol
pul = ps-ptol
pu2 = ps-dptol

c.Evaluate the nominal first derivative
if (x5 .gt. 0.d0) then
yd3nom = bet*(-aa1*x2+kx*x5*dsqrt(dabs(ps-x3)))/vol1
else
yd3nom = bet*(-aal *x2+kx*x5*dsqrt(dabs(x3)))/vol1
endif

c.Evaluate the saturated partial derivatives
pd32nom= -bet*aa1/voll
pd32low= pd32nom
pd32upp: pd32nom

pd32 = pd32nom

if ( (x3 .lt. p12) .and. (yd3nom .lt. 0.d0) ) then
if (x3 .lt. pll) then

pd32 = O.d0
else
pd32 = (1.d0+(x3-p12)/dptol)*pd3210w
endif
endif

if ( (x3 .gt. pul) .and. (yd3nom .gt. 0.d0) ) then
if (x3 .gt. pu2) then
pd32 = 0.d0
else
pd32 = ( l .d0-(x3-pu1)/dptol)*pd32upp
endif
endif

c.Retum to the function call
return
end

 

 

C:================: : =====:

72

c...Function pd33 - determines the saturated partial derivative pd33.
function pd33(x2, x3, x5, voll)

c.Variable declaration
real*8 pd33, x2, x3, x5, voll

c.Declare intermediate constants
real*8 pd33nom, yd310w, yd3upp
real*8 p11, p12, pul, pu2, dptol, yd3nom

c.Common variable type declaration
real*8 a22,a23,a24,a31,a32,a33,a35,a4],a42,a44,a45,a75,a76,
+ a77,aal,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,ccl 1,cc 12,cc13,cc14,cc15,cc16,cc17,curi,dcl 1,dh,frqr,
+ kl,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

common a22,a23,a24,a31,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aal,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,ccl 1,cc12,cc13,cc14,cc15,ccl6,ccl7,curi,dc1 1,dh,frqr,
+ k1,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

c.Define intermediate constants
dptol = ptol/2.d0
pll = dptol
p12 = ptol
pul = ps-ptol
pu2 = ps-dptol

c.Evaluate the first derivatives

if (x5 .gt. 0.d0) then

yd3nom = bet*(-aa1*x2+kx*x5*dsqrt(dabs(ps-x3)))/vol1
yd310w = bet*(-aa1*x2+kx*x5*dsqrt(ps-ptol))/vol1

yd3upp = bet*(-aa1*x2+kx*x5*dsqrt(ptol))/vol1

else
yd3nom = bet*(-aa1*x2+kx*x5 *dsqrt(dabs(x3)))/voll
yd310w = bet*(-aal *x2+kx*x5*dsqrt(ptol))/vol1
yd3upp = bet*(-aa1*x2+kx*x5*dsqrt(ps-ptol))/vol1

endif

c.Evaluate the saturated partial derivatives
if (x5 .gt. 0.d0) then
pd33nom=—0.5d0*bet*kx*x5/(dsqrt(dabs(ps-x3))*vol 1)
else
if (x5 .lt. 0.d0) then
pd33nom=0.5d0*bet*kx*x5/(dsqrt(dabs(x3))*vol1)
else
pd33nom=0.d0
endif
endif

73
pd33 = pd33nom

if ( (x3 .lt. p12) .and. (yd3nom .lt. 0.d0) ) then
if (x3 .lt. p11) then

pd33 = 0.d0
else
pd33 = (1.d0/dptol)*yd3low
endif
endif

if ( (x3 .gt. pul) .and. (yd3nom .gt. 0.d0) ) then
if (x3 .gt. pu2) then
pd33 = 0.dO
else
pd33 = (-1.d0/dptol)*yd3upp
endif
endif

c.Return to the function call
return
end

C:-"====-______—-—====== :======_—== ==

 

c...Function pd35 - determines the saturated partial derivative pd35.
function pd35(x2, x3, x5, voll)

c.Variable declaration
real*8 pd35, x2, x3, x5, voll

c.Declare intermediate constants
real*8 pd35nom, pd3510w, pd35upp, conl
real*8 p11, p12, pul, pu2, dptol, yd3nom

c.Common variable type declaration
real*8 a22,a23,a24,a31,a32,a33,a35,a4],a42,a44,a45,a75,a76,
+ a77,aal,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,ccl 1,cc12,cc13,cc14,cc15,cc16,cc17,curi,dc1 1,dh,frqr,
+ k1 ,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

common a22,a23,a24,a31,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aa1,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b7 l ,
+ bet,cc1 1,cc12,ccl3,cc l4,cc 15,cc 16,cc17,curi,dc1 1,dh,frqr,
+ k l ,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

c.Define intermediate constants
dptol = ptol/2.d0
p11 = dptol
p12 = ptol
pul = ps—ptol
pu2 = ps-dptol

74

c.Evaluate the nominal first derivative
if (x5 .gt. 0.d0) then
yd3nom = bet*(-aa1*x2+kx*x5*dsqrt(dabs(ps-x3)))/voll
else
yd3nom = bet*(-aal *x2+kx*x5*dsqrt(dabs(x3)))/vol1
endif

c.Evaluate the saturated partial derivatives
conl=(dsqrt(dabs(ps-x3))+dsqrt(dabs(x3)))/2.d0
if (x5 .gt. 0.d0) then
pd35nom=bet*kx*dsqrt(dabs(ps-x3))/vol1
pd3510w=bet*kx*dsqrt(ps-ptol)/vol 1
pd35upp=bet*kx*dsqrt(ptol)/vol1
else
if (x5 .lt. 0.d0) then
pd35nom=bet*kx*dsqrt(dabs(x3))/vol1
pd3510w=bet*kx*dsqrt(ptol)/voll
pd35upp=bet*kx*dsqrt(ps-ptol)/vol 1
else
pd35nom=bet*kx*con l/voll
pd3510w=pd35nom
pd35upp=pd35nom
endif
endif

pd35 = pd35nom

if ( (x3 .lt. p12) .and. (yd3nom .lt. 0.d0) ) then
if (x3 .lt. pll) then

pd35 = 0.d0
else
pd35 = (1.d0+(x3-p12)/dptol)*pd3510w
endif
endif

if ( (x3 .gt. pul) .and. (yd3nom .gt. 0.d0) ) then
if (x3 .gt. pu2) then
pd35 = 0.d0
else
pd35 = (1.d0-(x3-pu 1 )/dptol)*pd35upp
endif
endif

c.Return to the function call
return
end

 

C:===========m=========m==m=:======: =2
c...Function pd4l - determines the saturated partial derivative pd41.

function pd4l(x2, x4, x5, v012)

75

c.Variable declaration
real*8 pd4l, x2, x4, x5, v012

c.Declare intermediate constants
real*8 pd41nom, pd4llow, pd4lupp, conl
real*8 p11, p12, pul, pu2, dptol, yd4nom

c.Common variable type declaration
real*8 a22,a23,a24,a31,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aa],aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,ccl 1,cc12,cc13,ccl4,cc15,ccl6,ccl7,curi,dc1 1,dh,frqr,
+ k1,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

common a22,a23,a24,a31,a32,a33,a35,a4l,a42,a44,a45,a75,a76,
+ a77,aa],aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,cc1 1,c012,cc13,cc14,cc15,cc16,cc17,curi,dcl 1,dh,frqr,
+ k l ,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

c.Define intermediate constants
dptol = ptol/2.d0
p11 = dptol
p12 = ptol
pul = ps-ptol
pu2 = ps-dptol

c.Define the nominal first derivative
if (x5 .gt. 0.d0) then
yd4nom = bet*(aa2*x2-kx*x5*dsqrt(dabs(x4)))/vol2
else
yd4nom = bet*(aa2*x2-kx*x5*dsqrt(dabs(ps-x4)))/vol2
endif

c.Evaluate the saturated partial derivatives
conl = bet*aa2/(vol2*v012)
if (x5 . gt. 0.d0) then
pd4 lnom=con1 *(aa2*x2-kx*x5*dsqrt(dabs(x4)))
pd4llow=con1 *(aa2*x2-kx*x5 *dsqrt(ptol))
pd4 1 upp=con 1 *(a32*x2-kx*x5*dsqrt(ps-ptol))
else
if (x5 .lt. 0.d0) then
pd4 1nom=con1 *(aa2*x2-kx*x5 *dsqrt(dabs(ps-x4)))
pd4llow=con1 *(aa2*x2-kx*x5 *dsqrt(ps-ptol))
pd4 1upp=conl *(aa2*x2-kx*x5 *dsqrt(ptol))
else
pd41nom=con l *aa2*x2
pd4 1 low=pd4 l nom
pd4 1upp=pd4 l nom
endif
endif

pd4l = pd41nom

76

if ( (x4 .lt. p12) .and. (yd4nom .lt. 0.d0) ) then
if (x4 .lt. p11) then

pd41 = 0.d0
else
pd4l = ( l .dO+(x4-pl2)/dptol)*pd4 1 low
endif
endif

if ( (x4 .gt. pul) .and. (yd4nom .gt. 0.d0) ) then
if (x4 .gt. pu2) then
pd4] = 0.d0
else
pd4] = ( 1.dO-(x4-pu1)/dptol)*pd41upp
endif
endif

c.Retum to the function call
return
end

 

C==m==——-__..===::==========___—-==-—-.__.==___——-=====

 

c...Function pd42 - determines the saturated partial derivative pd42.
function pd42(x2, x4, x5, v012)

c.Variable declaration
real*8 pd42, x2, x4, x5, v012

c.Declare intermediate constants
real*8 pd42nom, pd4210w, pd42upp
real*8 p11, p12, pul, pu2, dptol, yd4nom

c.Common variable type declaration
real*8 a22,a23,a24,a31,a32,a33,a35,a4],a42,a44,a45,a75,a76,
+ a77,aa],aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,cc1 1,cc12,cc 13,cc 14,cc15,cc16,cc17,curi,dc1 1,dh,frqr,
+ k1,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

common a22,aZ3,a24,a31,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aa1,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,cc1 1,cc12,cc 13,cc l4,cc15,cc16,cc17,curi,dc 1 1,dh,frqr,
+ k1,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

c.Define intermediate constants
dptol = ptol/2.d0
pll = dptol
p12 = ptol
pul = ps-ptol
pu2 = ps-dptol

c.Define the nominal first derivative

77

if (x5 .gt. 0.d0) then

yd4nom = bet*(aa2*x2-kx*x5*dsqrt(dabs(x4)))/v012
else

yd4nom = bet*(aa2*x2-kx*x5*dsqrt(dabs(ps-x4)))/v012
endif

c.Evaluate the saturated partial derivatives
pd42nom=bet*aa2/v012
pd42low=pd42nom
pd42upp=pd42nom

pd42 = pd42nom

if ( (x4 .lt. p12) .and. (yd4nom .lt. 0.d0) ) then
if (x4 .lt. p11) then

pd42 = 0.d0
else
pd42 = (1.dO+(x4-p12)/dptol)*pd4210w
endif
endif

if ( (x4 .gt. pul) .and. (yd4nom .gt. 0.d0) ) then
if (x4 .gt. pu2) then
pd42 = 0.d0
else
pd42 = (1.d0-(x4-pu l )/dptol)*pd42upp
endif
endif

c.Retum to the function call
return
end

 

C========$==m============m==m===========

c...Function pd44 - determines the saturated partial derivative pd44.
function pd44(x2, x4, x5, v012)

c.Variable declaration
real*8 pd44, x2, x4, x5, v012

c.Declare intermediate constants
real*8 pd44nom, yd4low, yd4upp, yd4nom
real*8 p11, p12, pul, pu2, dptol

c.Common variable type declaration
real*8 a22,a23,a24,a31,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aa1,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b7 1,
+ bet,cc1 1,cc12,cc13,cc14,cc15,ccl6,ccl7,curi,dcl 1,dh,frqr,
+ kl ,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

common a22,a23,a24,a31,a32,a33,a35,a4],a42,a44,a45,a75,a76,

78

a77,aa1,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b7 l,
bet,ccl 1,cc12,cc13,cc14,cc15,cc16,cc17,curi,dcl 1,dh,frqr,
k1 ,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
ptol,ref,refc,vh,wn,xi

++++

c.Define intermediate constants
dptol = ptol/2.d0
p11 = dptol
p12 = ptol
pul = ps-ptol
pu2 = ps-dptol

c.Define the first derivatives

if (x5 .gt. O.d0) then
yd4nom = bet*(aa2*x2-kx*x5*dsqrt(dabs(x4)))/v012
yd4low = bet*(aa2*x2-kx*x5*dsqrt(ptol))/v012
yd4upp = bet*(aa2*x2-kx*x5*dsqrt(ps-ptol))/v012

else
yd4nom = bet*(aa2*x2-kx*x5*dsqrt(dabs(ps-x4)))/v012
yd4low = bet*(aa2*x2-kx*x5*dsqrt(ps-ptol))/v012
yd4upp = bet*(aa2*x2-kx*x5*dsqrt(ptol))/v012

endif

c.Evaluate the saturated partial derivatives
if (x5 .gt. 0.d0) then
pd44nom=-0.5d0*bet*kx*x5/(dsqrt(dabs(x4))*vol2)
else
if (x5 .lt. 0.d0) then
pd44nom=0.5d0*bet*kx*xS/(dsqrt(dabs(ps-x4))*v012)
else
pd44nom=0.d0
endif
endif

pd44 = pd44nom

if( (x4 .lt. p12) .and. (yd4nom .lt. 0.d0) ) then
if (x4 .lt. pll) then

pd44 = 0.d0
else
pd44 = (l.d0/dptol)*yd4low
endif
endif

if ( (x4 .gt. pul) .and. (yd4nom .gt. 0.d0) ) then
if (x4 .gt. pu2) then
pd44 = 0.d0
else
pd44 = (-1.d0/dptol)*yd4upp
endif
endif

c.Retum to the function call
return

 

79

end

c===__.-—======$3========mm============:===

 

c...Function pd45 - determines the saturated partial derivative pd45.
function pd45(x2, x4, x5, v012)

c.Variable declaration
real*8 pd45, x2, x4, x5, vol2

c.Declare intermediate constants
real*8 pd45nom, pd4510w, pd45upp, conl
real*8 pll, p12, pul, pu2, dptol, yd4nom

c.Common variable type declaration
real*8 a22,a23,a24,a31,a32,a33,a35,a4],a42,a44,a45,a75,a76,
+ a77,aa1,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b71,
+ bet,ccl 1,cc12,cc13,ccl4,cc 15,cc16,cc17,curi,dcl 1,dh,frqr,
+ k1,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

common a22,a23,a24,a31,a32,a33,a35,a41,a42,a44,a45,a75,a76,
+ a77,aa l ,aa2,aas,ac72,ac73,ac74,ac75,ac76,ac77,ampr,b7 1,
+ bet,ccl 1 ,cc12,cc13,cc l4,cc15,cc16,cc17,curi,dcl 1,dh,frqr,
+ k 1 ,k2,kf,kw,kx,la,lam,lh,mas,masd,pi,ps,psuppc,pslowc,
+ ptol,ref,refc,vh,wn,xi

c.Define intermediate constants
dptol = ptol/2.d0
p11 = dptol
p12 = ptol
pul = ps-ptol
pu2 = ps-dptol

c.Define the nominal first derivative
if (x5 .gt. 0.d0) then
yd4nom = bet*(aa2*x2-kx*x5*dsqrt(dabs(x4)))/v012
else
yd4nom = bet*(aa2*x2-kx*x5*dsqrt(dabs(ps—x4)))/v012
endif

c.Evaluate the saturated partial derivatives
con1=(dsqrt(dabs(ps-x4))+dsqrt(dabs(x4)))/2.d0
if (x5 .gt. 0.d0) then
pd45nom=-bet*kx*dsqrt(dabs(x4))/v012
pd45]ow=-bet*kx*dsqrt(ptol)/vol2
pd45upp=-bet*kx*dsqrt(ps-ptol)/vol2
else
if (x5 .lt. 0.d0) then
pd45nom=-bet*kx*dsqrt(dabs(ps-x4))/v012
pd45low=-bet*kx*dsqrt(ps-ptol)lv012
pd45upp=-bet*kx*dsqrt(ptol)/v012
else

80

pd45nom=-bet*kx*con 1/v012
pd4510w=pd45nom
pd45upp=pd45nom
endif
endif

pd45 = pd45nom

if ( (x4 .lt. p12) .and. (yd4nom .lt. 0.d0) ) then
if (x4 .lt. p11) then
pd45 = O.d0
else
pd45 = (l.dO+(x4-p12)/dptol)*pd4510w
endif
endif

if ( (x4 .gt. pul) .and. (yd4nom .gt. 0.d0) ) then
if (x4 .gt. pu2) then
pd45 = 0.d0
else
pd45 = (l .dO-(x4-pu1)/dptol)*pd45upp
endif
endif

c.Retum to the function call

return
end

 

 

APPENDIX C. TIME DOMAIN RESULTS

81

 

 

 

'
‘55
n 5. .
55
n 555 -
.5.l ...........................

“
. 555-.
5.5I
.5.5 ..................

 

 

 

' 0.95

0.8

0.85

0.8

0.75

Time [s]

0.7

Trace with we = 100 Hz. M: 501b, f=10 Hz

0.65

0.8
1me

T

.1 0.55
C.2.

 

 

 

Figure C.l. Time Trace with we = 100 Hz, M = 501b, f = 1 Hz
'-
l
I
l
l
l
l
I
i
I
I

18
16"
1.4---------
.2
1

 

Figure

82

 

 

 

 

M .\ .
‘ .-
. . . IIII.‘ I...“ .....
. n -l-‘ .......... .
. . “I|| ........... .
II .....

. tl‘ ......
. 8‘. .....
\.“ ...‘. a

r ...... . .. . .............................................. c ........ l e

' ....
. t‘-
‘1’- '
. ’l-
p ’
.‘t’.
. . v.
. . . ‘ ..
. . "4..‘.o.a..
. . “'.. ..... .
. . . “.‘ ...... .
. . I‘ ............. .
.......
. ItII .......
. ‘ ...
| .....
“ .. .
II N
. "
' .
".t“ . 9
.l....... .h ................ leen‘J. ....... . ........... .. ........ I e
I." o
.l’ .
. ,.’
. ..’V.
. II .
. .....
. . --..e-‘e -----
u . “.| 4.... .......
‘II‘..|..-‘. .........
. t‘ ......
‘ e.
. “ ...-- .
‘ a.
I O H
l
. " n
. I .
,';-.8805 .
v ‘ I
. . I'll,
’
. . IO,

1 ........................................... . .................... . A“
. . .‘ e
. ‘u ...

. .“.‘..u‘...... o
u '8‘- .........
. ‘.|I|'. ..........
. ......... n
. .....
‘ ...
. ‘ ...... a
‘ ..
‘ o e
ll H
8'...
'15-!“ .
'
.""
'e‘ .
. 'l-’ .
I ,’I
. . I
a u
l I ‘ O.
. u “ o
. . ‘IIIU.I.I ....... .0
.-
. p . illnin. ..... a a.
5 2 5 1 w 1 % 30
ac. 4| . el . 0

Time [5]

wc=100Hz,M=501b,f=25Hz

Figure C.3. Time Trace wi

 

 

d d N u .d
\ u .
. I e .
I O C. I
’ u u
U II. n
. 4.:
.o..”
. a... I
..k.
‘ .
v“
\.
\ .
a u
I’.
v" e
'.
......"
..... I

...................

. I
..&
‘ e
.II
I
....J
\J
. \ .
.4.”
u... I
C.“
\ ..

I.
0.95
Time [5]

 

8.9

 

1.25

Figure C.4. Time Trace with a)‘ = 100 Hz, M = 50 lb, f = 50 Hz

83

 

 

 

 

 

0.5 0.6 0.7 0.8 0.9

Time [s]

0.4

0.3

0.2

0.1

=1Hz

Figure C.5. Time Trace with we = 100 Hz, M = 1501b, f

 

 

 

 

 

u a .. -eris . W q
. .8 . .
.‘8.8 . .
\.
. ‘-
88"-
8 ........................... I'lst.r. ............................... L
"n"
'8'
'
8“
8
‘
u 8888 P
I. ............................. DI .............................
I |
‘-
. 88.818
n‘
\
. '8
8
- "'e"
8 ..... . ...................... ‘8! ................................. l
8"
c 8"
8
. 8'
c '-
. 8“
. 888
. 88
. 888
.l ..... . ....................... I I I ........................... I.
. .8
. 8.8.!
. .8.‘
c‘
.' "
. '8'."
. '0
.I ......................... 8.8 ................................ L
. '8..."
.l"
. ’.'o
u 88 ~
‘
888
. .888
8 eeeee u ....................... 8.. 8 ............................. 8
n 8'
. \‘.8‘
\.\
'
. 0'
‘l."8
'
.l .......................... 8 8“ .......................... . ..... l
"'
:8" .
2" .
’.(e
8 o
8 88
8888 .
r ............................... .8. 8 ............ u ............. u ..... l
l - u l I
. 8.888 .
. 8 . .
I8 . .
‘ .
. .
' . e
a.
8' .
8.
.I . 8' ... . . I.
.......................... . e8 .............................-......
. , ‘f‘o'.’ . .
e . . ‘5" .
. . 8'. .
n . o l a?
I I I A.
I I I ‘\ I.
. . I?“ .
. . 8
888
— b b I-“ b n
1 1 1 0 0

0.95

0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9
Time [5]

0.5

Figure C.6. Time Trace with we = 100 Hz, M = 1501b, f = 10 Hz

 

84

 

 

q #1 1 3.1.1.4...- - e1
. 88.8.... . u
“‘8 ..... .
I ....
,'nln. e
'8'rr .
8.8.8
‘88; .
. 88".
. u'
. 1..
U H 88.8...8...:
H 883-885-8- eeeee .
. 888. ...... .
8‘. ......

.
‘ nil
. 88 ......

_'\
I \
\
._\
l
0.95

ole
. .8
. 8 .
u 8 ee-
. 8 ....
. 113424.... .
. U 88.]. ..4-....” o
. . 88. .....
8IJ....
- ‘ .e o
8 ..
..\\ ..
. l
I, ...... s
. 11.01:... 91.
. . t.ll.'n'l
........... ...u... .... ...-.... "t.n..........-.......-..e.u .8 O
11. 0
. '8', 8 ...........
. -14. n
I I I I
. . .
. . 8
‘
n o ‘ ..-
‘ 'IOe
n ” 1.1.1.1242...
88. ..... .
. 88:... .
8 ..
. 8 ...I
.qii.... .
I .... .
,'k .....
,8'.t‘
- It
’I.
8“
8'
.8,
I
............................................. ............. .l.
e u‘ I
U 888.” ...... o
8 ....
‘n ..... .
888.‘.-8.....o
. 8 ......
88....
. 88 ..... .
8 ..
. 8 ...e .
8 ..
.\ .. .
Or I
' ...
' ...
11..1

o
e
on
e
en
en .
.-

 

 

 

1.12

0.95

Time [5]

we=100Hz,M=1501b,f=50Hz

 

 

1.1
-_1
'..\
Figure C.7. Time Trace with we = 100 Hz, M = 1501b, f = 25 Hz

11"
__1.08"'””"
11
1.09""

Figure C.8. Time Trace wi

1.08....m.

BIBLIOGRAPHY

BIBLIOGRAPHY

FitzSimons, P. M., paper to be presented at the 1994 ASME Winter Annual Meeting.

Hindmarsh, A. C., "ODEPACK, A Systematized Collection of ODE Solvers." Scientific
Computing, 1983, pp. 55-64.

Morari, M. and E. Zafiriou, Robust Process Control, Prentice-Hall, New Jersey (1989).

Nguyen, C. C, S. Antrazi, and Z. Zhou, "Analysis and Design of a Six-Degree-of-
Freedom Stewart Platform-Based Robotic Wrist." (submitted to NASA Goddard
Space Flight Center, February 1991).

Watton, J ., Fluid Power Systems - Modeling, Simulation, and Microcomputer Control,
Prentice-Hall, New Jersey (1989).

85

"lullnu“