NASA/TM-2002-212129

Multirate Flutter Suppression System Design for

the Benchmark Active Controls Technology

Wing

Part II: Methodology Application Software Toolbox

Gregory S. Mason and Martin C. Berg

University of Washington, Seattle, Washington

Vivek Mukhopadhyay

Langley Research Center, Hampton, Virginia

December 2002

https://ntrs.nasa.gov/search.jsp?R=20030011272 2018-05-15T15:38:46+00:00Z
The NASA STI Program Office... in Profile

Since its founding, NASA has been dedicated to the

advancement of aeronautics and space science. The

NASA Scientific and Technical Information (STI)

Program Office plays a key part in helping NASA

maintain this important role.

The NASA STI Program Office is operated by

Langley Research Center, the lead center for NASA's

scientific and technical information. The NASA STI

Program Office provides access to the NASA STI

Database, the largest collection of aeronautical and

space science STI in the world. The Program Office is

also NASA's institutional mechanism for

disseminating the results of its research and

development activities. These results are published by

NASA in the NASA STI Report Series, which

includes the following report types:

TECHNICAL PUBLICATION. Reports of

completed research or a major significant phase

of research that present the results of NASA

programs and include extensive data or

theoretical analysis. Includes compilations of

significant scientific and technical data and

information deemed to be of continuing

reference value. NASA counterpart of peerreviewed formal professional papers, but having

less stringent limitations on manuscript length

and extent of graphic presentations.

TECHNICAL MEMORANDUM. Scientific

and technical findings that are preliminary or of

specialized interest, e.g., quick release reports,

working papers, and bibliographies that contain

minimal annotation. Does not contain extensive •

analysis.

CONTRACTOR REPORT. Scientific and

technical findings by NASA-sponsored

contractors and grantees.

CONFERENCE PUBLICATION. Collected

papers from scientific and technical

conferences, symposia, seminars, or other

meetings sponsored or co-sponsored by NASA.

SPECIAL PUBLICATION. Scientific,

technical, or historical information from NASA

programs, projects, and missions, often

concerned with subjects having substantial

public interest.

TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific and

technical material pertinent to NASA's mission.

Specialized services that complement the STI

Program Office's diverse offerings include creating

custom thesauri, building customized databases,

organizing and publishing research results ... even

providing videos.

For more information about the NASA STI Program

Office, see the following:

• Access the NASA STI Program Home Page at

http://www.stLnasa.gov

E-mail your question via the Internet to

help@sti.nasa.gov

Fax your question to the NASA STI Help Desk

at (301) 621-0134

Phone the NASA STI Help Desk at

(301) 621-0390

Write to:

NASA STI Help Desk

NASA Center for AeroSpace Information

7121 Standard Drive

Hanover, MD 21076-1320
NASA/TM-2002-212129

Multirate Flutter Suppression System Design for

the Benchmark Active Controls Technology

Wing

Part II." Methodology Application Software Toolbox

Gregory S. Mason and Martin C. Berg

University of Washington, Seattle, Washington

Vivek Mukhopadhyay

Langley Research Center, Hampton, Virginia

National Aeronautics and

Space Administration

Langley Research Center

Hampton, Virginia 23681-2199

December 2002
Available from:

NASA Center for AeroSpace Information (CASI)

7121 Standard Drive

Hanover, MD 21076-1320

(301) 621-0390

National Technical Information Service (NTIS)

5285 Port Royal Road

Springfield, VA 22161-2171

(703) 605-6000
iii

CONTENTS

p_ggg

Section One: Introduction ........................................................................................................................................... 1

1.0 Overview ....................................................................................................................................................... 1

1.1 Software Requirements ................................................................................................................................. 1

Section Two: Background ............................................................................................................................................ 3

2.0 Overview ....................................................................................................................................................... 3

2.1 A Multirate Feedback System ...................................................................................................................... 3

2.1.1 The Sampling Schedule .............................................................................................................. 4

2.1.2 Digital Processor Gains .............................................................................................................. 4

2.2 The Equivalent Time-Invariant System ('ETIS) ........................................................................................... 4

2.3 Optimizing the Digital Processor Gains ....................................................................................................... 6

Section Three: M.File Documentation ....................................................................................................................... 7

3.0 Standard Variable Definitions ...................................................................................................................... 7

3.1 Modeling M-Files ......................................................................................................................................... 9

3.1.1 mr2etis ........................................................................................................................................ 9

3.1.2 mr2ptv ................................................................................................................................ •........ 9

3.1.3 ptv2etis ..................................................................................................................................... 10

3.1.4 sr2etis ........................................................................................................................................ 10

3.1.5 stack .......................................................................................................................................... 10

3.1.6 unstack ...................................................................................................................................... 11

3.2 Modeling M-Files ....................................................................................................................................... 12

3.2.1 mrsim ........................................................................................................................................ 12

3.2.2 mrfeedback ............................................................................................................................... 13

3.3 Synthesis M-Files ....................................................................................................................................... 14

3.3.1 Input Variables ......................................................................................................................... 14

3.3.2 mropt_imt ................................................................................................................................. 17

3.3.3 mropt_optim ....................................................................................................... 18

3.3.4 mropt extract ............................................................................................................................ 19

3.3.5 calc_LQGcost ........................................................................................................................... 19

Section Four: Synthesis Example .............................................................................................................................. 21

4.0 Problem Description ................................................................................................................................... 21

4.1 Compensator Design ................................................................................................................................... 22

4.2 Design Analysis .......................................................................................................................................... 26

References ................................................................................................................................................................... 27

Appendix A: Global Workspace Variables .............................................................................................................. 29

Appendix B: Support M-Files ................................................................................................................................... 31

Appendix C: M-File mkm .......................................................................................................................................... 37

Appendix D: Script mropt_mkm .............................................................................................................................. 39

SECTION ONE

INTRODUCTION

1.0 OVERVIEW

The design and analysis of even a simple multirate compensator can be a complex task. In this document we

describe some Matlab M-Files which aid in multirate design. They include M-Files for modeling multirate systems,

for computing optimum values of a multirate compensator' s gains, and for generating time domain simulations.

The remainder of the document is divided into three sections. In Section Two we review the basics of multirate

modeling and our optimization algorithm. We also present the key concepts and notation which are utilized in

Section Three. Section Three describes each M-File, detailing its inputs and outputs. The M-Files in this section

are divided into three categories: modeling, simulation, and synthesis. Finally, in Section Four we conclude with a

multirate design example.

1.1 SOFTWARE REQUIREMENTS

In addition to the standard Maflab toolbox routines, this software uses M-Files from both the Control_Toolbox

and the Optimization_Toolbox. These two toolboxes must be present for the multirate software to operate.
2

Blank
SECTION TWO

BACKGROUND

2.0 OVERVIEW

In the following paragraphs we present the notation used in the remaining sections and review some key

concepts which will be helpful to those using this software. A detailed explanation of the theory behind our M-Files

can be found in References 1-5.

2.1 A MULTIRATE FEEDBACK SYSTEM

A multirate sampled-data system consists of a continuous plant in feedback with a multirate compensator. A

block diagram of such a system is shown in Figure 2.1. The vector Ys in this figure represents the continuous plant

sensor output. Each element of Ys can be sampled at an independent rate. The vector _- represents the sampled

value ofy s available to the digital processor. (In our double index notation, the discrete signal p(m,n) results from

sampling the continuous signal p(t) at the times t = (raN+ n)T; where the integer N is the period of repetition; T is

the sampling time; m =0, 1.... ; and n=0, 1..... N-1.) The digital processor obtains the current value off and

combines it with the current processor state vector, _-, using the state space structure shown in Figure 2.1. Each

element of the processor state vector, _-, can be updated at an independent rate. The continuous output from the

compensator, represented by the vector u, is formed by holding the output from the digital processor,m, with a zeroorder-hold. Each element of the vector _ can be held at an independent rate to form u. The vectors v and w

represent the discrete sensor and continuous process noise respectively.

Conceptually, one can divide the multirate compensator into two parts, the "sampling schedule" and the digital

processor gains. This is the approach used in the synthesis and analysis software. The "sampling schedule" is a

description of when each compensator input is sampled and when each output and processor state is updated, while

the digital processor gains determine the dynamics of the digital processor. In the following paragraphs we discuss

each in detail.

Continuous Plant

Yc(O

Ys(O

_m,n+ 1) = Az(n)zT-mn)+Bz(n)y(_n)

-u(m,n) = CzXn)zTrn, n) + Dz(n)y(_,n)

[ Zero Order Hold(s)

[ Digital Processor

- Multirate Compensator

Figure 2. l Closed-loop Multirate System

pier(s)

v(m,n)
2.1.1 The Sampling Schedule

In general, the sampling and updating of Ys, z-, and _ in Figure 2.1 can occur at any time. We, however, require

that these sample and update activities occur only at integer multiples of some fixed time, called the shortest time

period (STP). The actual value of the STP is arbitrary, but is often a function of the hardware and software used to

implement the control law. We also require that the sampling and updating activities of the sensors, states and

outputs repeat themselves after some fixed period of time. The period of repetition of the sampling schedule is

called the basic time period (BTP). Finally, we define

BTP

the integer N - sT P and the value T=STP (2.1)

In our double index notation, the first index (m) in p(m,n) indicates the integer number of BTP's which have elapsed

when the sample/update occurred and the second index (n) indicates the integer number of STP's which have

elapsed within the current BTP when the sample/update occurred.

We can represent the sampling schedule for the multirate compensator graphically, as shown in Figure 2.2.

The figure shows a time line for each sampler, processor state, and zero-order-hold. The time line is divided into

one STP increments. On the left side of the time line is a description of the signal or state, being sampled or

updated. On the right side is a description of the particular activity represented by the time line, e.g. state update,

sampler, or zero-order-hold. Circles on each time line indicate when a sample or update activity associated with that

particular signal or state takes place. Usually the sampling schedule is shown for only one BTP since the sampling

schedule repeats itself every BTP.

In most applications, the sampling/updating activities for a given sensor, output or state will be periodic within

the BTP, as is shown in Figure 2.2. However, the sampling/updating activities do not have to be periodic within the

BTP. The only requirement is that the sampling/updating activities have some period of repetition (the BTP) and

that they occur at integer multiples of the STP. Once the STP and BTP have been selected, the designer can

arbitrarily specify sampling/updating activities at any multiple of the STP with in one BTP. An example of a

multirate sampling schedule in which the sampling/updating activities are not periodic within the BTP is shown in

Figure 2.3. A sampling policy like this might be used to multiplex multiple inputs through a single analog to digital

converter.

2.1.2 Digital Processor Gains

The processor gains are the values of the matrices A z, Bz, Cz, and D z in Figure 2.1. Like the sampling schedule,

they can be periodically time-varying with a period of repetition of one BTP. Generally, these matrices are free

design parameters which can be adjusted by the designer to improve the performance of the multirate compensator.

The synthesis software discussed later in this document can be used to calculate optimum values for these gains.

22 THE EQUIVALENT TIME-INVARIANT SYSTEM (ETIS)

A multirate compensator with the structLtre discussed in Section 2.1 can be modeled as a periodically timevarying single-rate compensator by appending appropriate hold states to the digital processor model. This new

compensator has the form

x(m,n+l) = A(n)x(m,n)+ B(n)u(m,n) (2.2a)

ys(m,n) = C(n)x( m,n) + D (n)u(m,n) (2.2b)

for n = 0,1,...N-I, and m = 0,1 ....

This compensator has a sampling period of one STP and a period of repetition of one BTP (or NT). y s(m,n)

represents the values of y(t) sampled every STP; u(t) is formed by holding u(m,n) with a zero-order-hold.
Zl _1 I VI i ,9 I ', ,9 I I ,9 I I ? State Update

Ysensor ? I ,, ,91 ,, ,9,, _, ,£_, ,, ,9,, _, ,£ Sampter

Ucontrol _ F ', : I ? ', ', _, ', ? '_ I I I ? ZOH

0 3T 6T 9T 12T 15T

STP _._ I_- BTF _._l

? = sample/update activity

N=15

Figure 2.2. Example multirate sampling schedule with periodic sampling/updating activity

Ysensor#1 I ,_ ? ,_ ', ? ? ?

Ysensor #2 £, ,' ,' I ? I I a

0 T 2T 3T 4T 5T 6T 71"

STP 14--

"_ BTI_

Sampler

_) = sample/update activity

Sampler N=4

Figure 2.3. Example multirate sampling schedule with aperiodic sampling activity

The periodically time-varying single-rate compensator can be further transformed into a single-rate Equivalent

Time-lnvariant System (ETIS) with the form shown below.

x(m+l,0) = AEX(m,O) + BEuE(m,0) (2.3a)

yE(m,0) = CEX(m,0) + DEuE(m,0 ) (2.3b)

where

Lys(m ,N- 1)J Lu(m _N- l)J

We use the subscript E to denote vectors and matrices strictly associated with the ETIS. Notice that the ETIS

input/output vectors are composite vectors containing the input/output values of the multirate (or periodically-time

varying) system at N sampling times. Consequently, an ETIS is always MIMO even if the original system is SISO.

If the multirate system has p inputs, q outputs and a sampling period of one STP then the ETIS is a single-rate linear

time-invariant system with Np inputs, Nq outputs and a sampling period of one BTP.

The ETIS is fundamental to the analysis of multirate systems. It allows one to evaluate the performance and

stability of complex systems comprised of multi.rate, periodically time-varying and/or single-rate components using

only modified linear time-invariant single-rate techniques. For example, to evaluate the stability of the system in

Figure 2.1, we would first transform the multirate compensator into its ETIS with a given value for N. Then we

would discretize the plant at the STP of the compensator using a zero-order-hold and transform the resulting singlerate system into an ETIS using the BTP of the compensator. Next, the plant and compensator ETIS's could be

combined in feedback just as if they were traditional single-rate systems. Finally, we could determine the stability

of the original multirate sampled-data system from the eigenvalues of its closed-loop ETIS.
2.3 OlWI'IMIZING TIlE DIGITAL PROCESSOR GAINS

The synthesis software discussed later calculates optimum values of the digital processor gains, A z, B z, C z, and

D z, by minimizing a quadratic cost function which reflects the performance of the closed-loop system in Figure 2.1.

The cost function has the form

r

2"= lira EffyT(t) u (t)]IM T MIVYc(t)ll t---_*_ t" Q2JL u($) ]J (2.5)

where Y is the cost associated with the closed-loop system shown in Figure 2.1. The vector Yc is the continuous

criterion output and u is the continuous conlrol input Q1, Q2 and M are the cost function weighting matrices which

are selected by the designer.

This cost function has the same form as that minimized by a continuous time LQR design. Thus the cost

associated with the optimized multirate compensator and that of an LQR design can be compared directly. The

designer can also use this fact to help select appropriate values for Q1, Q2 and M.

To improve the robustness of the compensator, the optimization algorithm can optimize the digital processor

gains for several different plant perturbations simultaneously. The resulting compensator will stabilize the plant at

each perturbation and provide overall optimum performance. This is accomplished by minimizing the new cost

function of Eqn. (2.6) which is the sum of the costs associated with each plant perturbation.

J= Ji = YmaE yT(t) uiT(t)

i=I = t--_** /_ J[. u/(t) .JJ

(2.6)

Here 7/is the cost associated with the ith plant perturbation and there are Np plant perturbations.

Optimum values of Az, Bz, Cz, and D z, occur when

9J 9J OJ OJ

--,_z =0, _-z =0, _-_z =0, and _-_z=O (2.7)

Our algorithm uses a gradient type numerical search to determine values of the digital processor gains such that the

conditions in Eqn. (2.7) are satisfied. Because the synthesis software uses an iterative process to determine optimum

values for the digital processor gains, the user must provide the software with an initial guess for Az, Bz, C z, and D z.

The compensator corresponding to these values must stabilize every plant perturbation considered in (2.6).

Obtaining a suitable stabilizing initial guess can sometimes be a difficult problem. We refer the interested reader to

Reference 1.
SECTION THREE

M-FILE DOCUMENTATION

3.0 STANDARD VARIABLE DEFINITIONS

Many of the M-Files discussed in this document require similar input variables or provide similar outputs. To

simplify the documentation of these M-Files a set of standard variables are used throughout this document. They are

defined in Table 3.1 with Matlab variables and functions bolded.

Table 3.1 Standard Variable Definitions

Variable Description

pit Plant matrices in the form

plant is

nplt

cmp

or

pit = [F, G; H, J] where the state-space representation of the

Jc(t) = Fx(t) + Gu(t) (3. la)

y(t) = nx(t)+ Ju(t) (3.1 b)

x(m,n+l) = Fx(m,n) + Gu(m,n)

y(m,n) = Hx(mm) + Ju(m,n)

depending on the value of stp defined later.

Number of states in (3. I), or equivalently the number of rows in F

(3.1c)

(3.1 d)

Multirate compensator gain matrices. Given the state-space representation of the digital

processor

2(m, n + 1) = Az(n)5(m, n) + Bz(n)y(m, n) (3.2a)

_(m,n) = Cz(n)_(m,n)+ Dz(n)y(m,n ) (3.2b)

where m=0,1 .... and n----O,l,...,N-1, so that the gain matrices, Az, Bz, Cz and D z, are

periodically time-varying, crop has the form

-[[-Az(O) "z<l l,Az(N-I)

cmp-LLCz(O ) Dz(O)j, LCz(1) Dz(1)j...,LCz(N_I) Dz(N-I)]] (3.3)

IrA z, Bz, Cz, and Dz are constant, then

=F/=(o)8z(0) 1

crop LCz (0) Dz(O)] (3.4)

and it is assumed that Az(O)=Az(1) ..... A z(N-1), Bz(O)=Bz(1 )..... B z(N-1 ), etc.

The software automatically deduces from the size of crop, the value of ncmp and the

number of compensator inputs whether the digital processor gains are periodically timevarying or not.

Continued onfollowing page...
Table 3.1 Standard Variable Definitions (continued from previous page)

Variable Description

nemp

stp

stppbtp

SU

sy

SZ

Number of digital processor states in Eqn. (3.2), or equivalently the number of rows in Az.

If stp > 0 then stp is STP, the shortest sample/update period of the compensator. See

Section 2.1.1, and pit describes a continuous plant.

If stp = 0 then the plant matrices pit describe a discrete system whose sampling rate is

one STP of the compensator, and pit describes a discrete plant.

Number of STP's per BTP. stppbtp = N. See Section 2.1.1.

Sampling schedule for the compensator output, su has two forms

Case I: su is a matrix, su must have as many columns as there are compensator outputs,

and must have N rows (N = BTP/STP). Each element of su has a value of 1 or 0. A 1

in the iah column andj _ row of su indicates that the i _ compensator output is updated

at the 3_ STP within the current BTP. A 0 indicates no update and the previous value

is held. For example, given a two output system with the following sampling schedule

Output#1 _ ZOH

? = update activity

Output #2 _ ZOH N=3

STP_ _--- BTP

the corresponding value of su is

(3.5)

Case II: su is a row vector. The i fi column of su specifies the sampling period of the /th

compensator output hold in number of STP's. For example, su=[1 4] specifies that the

first output is updated every STP and the second output is updated every fourth STP.

This form assumes that the updating is synchronized on the first STP of the BTP.

Sampling schedule for the compensator input, sy has the same form as su except that it

specifies when the plant output (the compensator input) is sampled.

Sampling schedule for the compensator states, sz has the same form as su except that it

specifies when the compensator's digital processor states are updated.
3.1 MODELING M-FILES

3.1.1

Format:

Descrip_on:

mr2etis

[ae,be,ce,de]= mr2etis(emp, nemp, su, sy, sz, stppbtp)

Converts a multirate system into an ETIS with the form:

I y(m,0) ]

v(m, n + 1) !

x(m + 1,0) = aex(m, O) + be " .

[y(m, N- 1)J

I u(m,O) I [ y(m,O) ]

uCm, !n+l) =eex(m'0)+ael . ly(m,n+l)| " I

[u(m, N - 1) [y(m, N- 1)J

Inputs:

Outputs:

cmp, ncmp, su, sy, sz, stppbtp The compensator description. See Table 3.1.

ae, be, ee, de The ETIS matrices in Eqn. (3.6)

(3.6a)

(3.6b)

3.1.2

Format:

Description:

Inputs:

Outputs:

mr2ptv

[P, nu, ny, nz]=mr2ptv (cmp; ncmp, su, sy, sz, stppbtp)

Converts a multirate system into a periodically time-varying system of the form

x(m,n +1) = A(n)x(m,n) + B(n)y(rn,n)

u(m,n) = C(n)x(m,n) + D(n)y(m,n)

cmp, ncmp, su, sy, sz, stppbtp The compensator description. See Table 3.1.

P The periodically time-varying system matrices in Eqn. (3.7), where

p=IID(O) C<Oqp<'>+'<'>]"

LLB(O) AC0)J'LB(1) AC1)J' ""LB(N-1) ACN-1)JJ

nu Number of inputs u in Eqn. (3.7)

ny Number of outputs y in Eqn. (3.7)

nz Number of states x in Eqn. (3.7)

(3.7a)

(3.7a)

(3.8)
10

3.1.3

Format:

Description:

Inputs :

Ou_u_:

ptv2etis

[ae, be, ce, de]=ptv2etis(P, nu, ny, nz, st-ppbtp)

Converts a periodically time-varying system into its ETIS of the form of Eqn. (3.6)

P, nu, ny, nz The periodically time-varying system matrices. See Section 3.1.2

stppbtp See Table 3.1

ae, be, ee, de The ETIS matrices. See Eqn. (3.6)

3.1.4

Format:

Description:

Inputs ".

Output:

sr2etis

[ae, be, ee, de]=sr2etis(a, b, c, d, stppbtp)

Converts the single-rate system in Eqn. (3.9) into an ETIS with an N of stppbtp, and the form of

Eqn. (3.6)

x(m,n +1) = a(n)x(m,n) + b(n)y(m,n) (3.9a)

u(m,n) = c(n)x(m,n) + d(n)y(m,n) (3.9b)

a, b, e, d The single-rate systems matrices. See Eqn. (3.9)

stppbtp The desired N of the new ETIS. See Table 3.1

ae, be, ce, de The ETIS matrices. See Eqn. (3.6)

3.1.5

Format:

Description:

Inputs:

Outputs:

stack

[y]= stack (w, stppbtp)

Converts a traditional discrete-time vector w into an ETIS vector y

w A matrix whose columns contain the values of the vector y(m,n) for Nk successive samples times,

where k is an integer. The matrix w has the form

w= [y(0,0), y(0,1) ..... y(0,N-1) y(1,0), y(1,1) ..... y(1,N-1) ..... y(k-l,N-1)] (3.10)

y is written using the double index notation of Section 2.1.

stppbtp See Table 3.1.

y An ETIS vector of the values ofyE (the ETIS of y) with the form

Y=[YE(O,O),yE(1,O) ..... YE(k_l,O)]=ll y(0,1)i y(1,1)! ... y(k-l,1)i (3.11)

[[_y(0,N - 1) [y(1,N-1) [.y(k-l,N-1)J
11

3.1.6

Format:

Description:

Inputs:

Outputs:

unstack

[w]=unstack(y, stppbtp)

Converts an ETIS vector y into a traditional discrete-time vector w.

y An ETIS vector of the form of Eqn. (3.11)

stppbtp See Table 3.1

w A vector of the values ofy(m,n) as in Eqn. (3.10)

NOTE: Additional M-Files used in support of those described in this section are described in Appendix B.
12

3.2 SIMULATIONM-FILES

3.2.1

Format:

Description:

mrsim

[y, t] =mrsim (plt, nplt, cmp, ncmp, su, sy, sz, stppbtp, simpstp, inputl, outputl, ur, stp, XO)

Generates a time response simulation of a continuous or discrete plant in feedback with a multirate

compensator. The routine assumes positive feedback as shown below with the continuous inputs

and outputs of interest sampled or updated every stp/simpstp seconds. This configuration allows

the user to simulate the inter-sample response of the muir_irate system.

., J- yp

( l Plant Tf

Compensator Multirate _

"If = $tp/simpstp

Inputs:

Outputs:

pit, nplt A description of the plant, see Table 3.1.

crop, ncmp, su, sy, sz, stppbtp,A description of the multirate compensator, see Table 3.1.

simpstp This routine can calculate the response of the closed-loop system between the

sample/update activities when the plant is continuous (stp > 0). simpstp is the number of

simulation steps per STP of the compensator, simpstp must be an integer greater than zero and

simpstp = 1 if sip = 0.

inputl A vector specifying which plant inputs are connected to the compensator outputs. For

example inputl = [1 4] indicates that the compensator outputs are connected to the first and

fourth plant inputs.

outputl A vector specifying which plant outputs are connected to the compensator inputs. For

example outputl = [2 4] indicates that the compensator inputs are connected to the second and

fourth plant outputs.

u¢ The discrete input vector containing the values of the input at every (stp/slmpstp) time for

k-stppbtp-simpstp sample times, ur has the form

t=r = [Up(0,0), Up(0,1) ..... up(k-l,stppbtp.simpstp- 1)] (3.12)

Note that in Eqn. (3.12) the sample/update period is one stp/simpstp. This sample period should

not be confused with the sample period of the multirate compensator which is simply one stp.

The second index in the double index notation in Eqn. (3.12) takes on values between zero and

simpstp- stppbtp.

stp see Table 3.1.

X0 A vector containing the initial conditions of the plant state vector. If X0 is omitted, the initial

conditions are assumed to be zero.

y A vector containing the response of the plant every (sttgsimpst p) time. y has the form

y = [yp(0,0), yp(0,1) ..... yp(k-l,stppbtp.simpstp- 1)] (3.13)
13

ComlTzenls :

3.2.2

Format:

Description:

Inputs;

Outputs :

Note that the second index in the double index notation corresponds to the second index of

Eqn. (3.12)

t a vector containing the times corresponding to the columns of y. When stp= 0 then t indicates the

number of sampling periods.

This M-File uses the ETIS of the system to compute the time domain response. Hence, if k in

Eqn. (3.12) is not a multiple of simpstp.N where (N=BTP/STP) the system response is only

simulated to the nearest multiple of simpstp. N

mrfeedback

[ae, be, ce, de] =mrfeedback(plt, nplt, cmp, ncmp, su, sy, sz, stppbtp, inputl, outputl, stp)

Creates a closed-loop ETIS system from a discrete or continuous plant, depending on stp, and a

multirate compensator using positive feedback.

pit, nplt The plant description. See Table 3.1.

cmp, ncmp, su, sy, sz, stppbtp The multirate compensator description. See Table 3.1.

input1 A vector specifying which plant inputs are connected to the compensator outputs. For

example input1 = [ 1 4] indicates that the compensator outputs are connected to the first and

fourth plant inputs.

output1 A vector specifying which plant outputs are connected to the compensator inputs. For

example output1 = [2 4] indicates that the compensator inputs are connected to the second and

fourth plant outputs.

stp See Table 3.1. If stp is omitted mrfeedback assumes stp=0.

ae, be, ce, de The system matrices of the ETIS closed-loop system with the form of Eqn. (3.6). The

inputs and outputs of this system are the inputs and outputs of the discrete plant
14

3.3 SYNTtIESlSM-FILES

Thesynthesisalgorithmisimplementedasaseriesof scriptEAes. They allow the designer to optimize the

performance of the closed-loop system shown in Figure 3.1 by calculating values of Az, Bz, C z, and Dz, such that the

cost function of Eqn. (2.6) is minimized. There are three distinct steps to the optimization process: 1) entering the

data which describes the problem; 2) discretizing and preprocessing the data; and 3) optimizing the compensator

gain values. Each step of the optimization process is described in detail in the following paragraphs.

Although the synthesis algorithm was designed to solve the multirate sampled-data problem, it is fairly general

and can be applied to a variety of related problems as well. It can compute optimum compensator gains for a

sampled-data system consisting of a continuous plant in feedback with a discrete compensator as shown in

Figure 2.1, or for a discrete system consisting of a discrete plant in feedback with a discrete compensator. In either

case, the compensators may be single-rate or multirate and may have either time-invariant or periodically timevarying digital processor gains.

3.3.1 Input Variables

The synthesis algorithm looks for specific Maflab "workspace" variables to def'me the optimization problem.

(A "workspace" variable is defined from the >>prompt as opposed to being passed as an argument to a function.)

These variables define things such as the plant dynamics, the cost function and the compensator structure and are

defined in Table 3.2. The user must assign values to these variables before beginning the optimization. This can

either be done manually or automatically from a script. An example of a script which defines these variables is

provided in Appendix D.

: u(m,n)

Zero Order Hold(s)

W

Plant

Continuous or Discrete

_(rn;n+ I)= Az(n)z_mn) + Bz(n)y(_n)

-u(m,n) = Cz(n)zTmn)+Dz(n)y(_n)

Digital Processor

Figure 3.1 Closed-loop Multirate System

Yc

Ys

_rn, n)

'4_ampl_r(s)

_(mn)
15

Table 3.2 Matlab Workspace Variables Defining Optimization Problem

Variable Description

r_

pit

The number of plant perturbations. Np=l indicates only the nominal pant is defined. See

Section 2.3

State-space matrices of the plant(s). Given the i_ continuous plant perturbation

xi(t ) = Fixi(t)+ Gi_ ui(t) l (3.14a)

(wi(t)J

Ysi(t)J (3.14b)

or the i_ discrete plant perturbation, depending on the value of stp defined later,

xi(m, n + 1) = Fixi(m, n) + Giui(m, n) + wi(m, n) (3.14c)

Ysi(m,n)J

where u is the control input; wi is zero mean white noise input.; Yci is the criterion output;

and Ysi are the outputs which are sampled by the compensator. Then

G2, _P

pit=[[H11 jGIII, IHF_ j21...F FNp,LHN p JNpl ) (3.15)

For convenience we assumed in Eqns. (3.14) that the input vector is partitioned into

the control input ui and the disturbance input wi. The software, however, can accommodate

a more general form for the input vector. Using the variables ueol and neol, defined later,

the user can specify which inputs correspond to control and noise inputs respectively. Any

row of the input vector can be either control, disturbance or both. (We interpret an input

vector specified as both control and disturbance as two separate inputs, one a control input

and the other a disturbance input whose distribution vectors have the same values.)

Similarly, the outputs Yci and Ysi are specified by crow and srow respectively, and

any row of the output vector can be criterion output, sensor output or both.

NOTE: For a discrete plant, the process noise distribution matrix is unity.

A vector of the PSD values of the white noise disturbance wi. If

WiS ( _) = E{w i(t - z)wT(t)} (3.16a)

or if Wi(_(k) + S(j)) = E{wi(m, n)wT(m + k, n + j)} (3.16b)

depending on the value of stp, then

wxx = [W 1 W2 ... WNp]

continued on following page...
16

Table 3.2 Maflab Workspace Variables DeEming Optimization Problem (continued from previous page)

Variable Description

vIg

where

then

The PSD values of the discrete sensor noise covariance. If the output, sampled every STP,

is

yst(k ) = Hixi(k ) + vi(k ) (3.17)

VirS(j ) = E{vi(k)vT(k + j)}

w= [vl v2 ... VNp]

neol A vector specifying which columns of Gi correspond to w i. If neol = [1 4], then the plant

disturbance vector w is comprised of the first and fourth plant inputs.

ueol A vector specifying which columns of G i correspond to ui. If ueoi = [2 4], then the plant

control input vector u is comprised of the second and fourth plant inputs. (Again, we

interpret an input vector specified as both control and disturbance as two separate inputs,

one a control input and the other a disturbance input whose distribution vectors have the

same values.)

srow A vector specifying which rows of Q correspond to Ysi. If stow = [1 3], then the first and

third continuous plant outputs are sampled and connected to the compensator.

crow A vector specifying which rows of Q correspond to Yci. If srow = [1 3], then the first and

third continuous plant outputs are used to calculate the cost given in Eqn. (3.18).

qla, q2a, maThe cost function weighting matrices for either a continuous or discrete cost function,

depending on the value of st1) defined later. The cost function has the form

{[ T ]IQIi MilFYci(t)ll (3.18a) Jcontinuous = lira E yT(t) ui (t) MT Q2iJ[ ui(t) JJ

f --q,_ *.= 1

or

1 N stp_p- 1

LM?

where i corresponds to a particular plant condition.

qla = [Qll Q12... OlNp]

q_ = [Q21022... L_Np]

ma = [N1 N2 ... NNV]

qa and ra must be symmetric and positive semi-definite.

If ma is left undefined its value is assumed to be zero.

(3.18b)

continued on following page...
17

Table 3.2 Matlab Workspace Variables Defining Optimization Problem (continued from previous page)

Variable Description

emp

cmpfree

SU, by, SZ

S_pb_

stp

The compensator gain matrices as defined in Table 3.1. The digital processor gains can be

either periodically time-varying or time-invariant.

IMPORTANT: The user must provide an initial guess for crop which stabilizes all plant

perturbations.

This matrix specifies which of the gains in crop are fixed, empfree has the same

dimensions as crop, but its elements are either 0 or 1. A 0 indicates that the corresponding

element in crop is fixed and should not be optimized. A 1 indicates that the corresponding

element in crop is free and can be optimized.

The sampling schedule description. See Table 3.1. By selecting appropriate values for su,

sy, sz and Stppbtp the designer can specify whether the compensator is multirate or singlerate.

The same as stppbtp in Table 3.1. Note that the first letter of Stppbtp is upper case,

indicating this is a global workspace variable. See Section 3.3.2.

STP as defined in Table 3.1.

Ifstp > 0 the algorithm assumes that the plant (pit), the processor noise (wxx) and the

cost weighting matrices (qla, q2a, and ma) are all defined for a continuous

plant in feedback with the discrete compensator.

If stp = 0 the algorithm assumes that the plant (pit), the processor noise (wxx) and the

cost weighting matrices (qla, q2a, and ma) are all defined for a discrete plant

in feedback with the discrete compensator, and that the sampling period of the

discrete plant is one STP of the mullirate compensator

3.3.2

Fo rrruzt :

Description:

Inputs •

Outputs :

Comment.

mropt_init

mropt_init

This script operates on the workspace variables defined in Table 3.2. It discretizes (if necessary) the

plank process noise and cost function; and defines a set of global workspace variables used by the

optimization script (mropt__opfim). This script must be rerun after making any changes in the

workspace variables defined in Table 3.2, with exception to changes in the values of crop or

empfree.

mropt__init is a script which uses the workspace variables defined in Table 3.2.

mropt__init save it's output to the global workspace variables defined in Appendix A.

The global workspace variables defined by mropt_init are used by both the optimization script

mropt_optim and by the M-File calc_LQGeost. All global variables generated by mropt_init

begin with a capital letter. This helps to differentiate them from other workspace variables. A

brief description of these global variables is given in Appendix A.
18

3.3.3

Format:

Description:

Inputs:

Ou_u_:

Comments:

mropt_optim

mropt_optim

mropt_optim computes optimum values of the digital processor gains such that the cost function

defined in Eqn. (3.18) is minimized.

mropt_optim uses the global workspace variables generated by mropt_init in addition to the three

workspace variables defined in Table 3.3.

The script returns a value of the cost function, the gradient of the cost function with respect to the

free compensator parameters, and the optimum values of the digital processor gains. These

values are automatically displayed on the computer screen and simultaneously stored in the

workspace variables defined in Table 3.4.

1) As already noted, if the user modifies the problem definition by changing the values of the

variables defined in Table 3.2, the script mropt_init must be rerun before running mropt_optim.

The exceptions are changes to the variables defined in Table 3.3.

2) mropt_optim can be computationally intensive. The user can abort the optimization with a

Control-C key sequence. Also see mropt_extraet, Section 3.3.4

Table 3.3 Input Workspace Variables for mropt_oplim

Variable

cmp

cmpfa-ee

dseale

Description

This variable contains the compensator gains and is the same variable as described in

Table 3.2. These are the starting values for the optimization and must stabilize all plant

perturbations.

This matrix specifies which of the gains in crop are fixed and is the same variable as

described in Table 3.2.

During optimization the compensator gains are scaled using a variable named dseale to help

improve the numerical accuracy of the search. They are scaled as follows

Xscaled = inv (dseale).Xunscaled

where Xunscaled is an unscaled vector containing the free compensator gains and Xscaled

are those values scaled. Typically dseale is a diagonal matrix chosen so that all the

elements in Xscaled have the same magnitude. The user can either manually set dseale or

have mropt optim set it automatically. If dseale exists and has the correct dimensions

then mropt_optim uses that value of dseale, otherwise a new value for dseale, based on the

current values in crop, is automatically calculated. Typing "elaar d.mc,,lo" before

running mropt_optim forces the routine to scale the problem based on current values in

emp.
19

Table3.4OutputWorkspace Variablesformropt optim

Variable

cmpf

jcost

djdp

Description

This contains the optimized values of emp after mropt_optim has completed the

optimization.

This contains the value of the cost function after mropt_optim has completed the

optimization.

This contains a scaled value of the gradient of the cost function with respect to the

compensator values after mropt optim has completed the optimization.

3.3.4

Format:

Description:

Inputs:

Outputs:

Comments:

mropt_extract

mropt_extract

After aborting the script mropt_optim, usually with a Control-C key sequence, the script

mropt_extraet can extract the value of the "optimized" compensator at the last iteration.

None. Can only be rerun after aborting ml'opt_oplim.

The optimized values of the digital processor gains at the last iteration before mropt_optim was

aborted are displayed on the computer and stored in the workspace variable cmpf.

Sometimes it is necessary to abort mropt_optim before is has completed the optimization. The user

might wish to rescale the free parameters by changing dseale, or might decide to free up

additional gains in cmp by changing cmpfree. The user can restart the optimization from the last

iteration by setting emp=onpf and rerunning mropt_opfim.

3.3.5

Format:

Description:

Inputs:

Outputs:

Comments:

cale_LQGcost

[jcost ]= calc_LQGcost (j)

Calculates the discrete LQG cost assuming a sampling period of one stp, This routine uses the

global workspace variables computed by mropt_init. Use this routine to find the minimum value

of the cost function attainable by mropt_oplim

j = 0 or is undefined, then the routine calculate the sum of the LQG costs for all Np plant conditions

corresponding to the cost Jdiscrete given in Eqn. (3.18).

0<j<Np then the routine calculates the LQG cost for only the j#_ plant condition

jcost The value of the LQG cost function

cale_LQGcost uses the Controls_Toolbox routines dlqr and dlqe. They require that ra in Table 3.2

be posilive-defmite.

NOTE: Additional M-Files used in support of those described in this section are described in Appendix B.
20

Blank
21

SECTION FOUR

SYNTrIESIS EXAMPLE

4.0 PROBLEM DESCRIPTION

In this section we present an example design problem which utilizes the M-Files discussed in the previous

section. Our objective is to design a multirate compensator for a lightly damped mass-spring-mass system.

The mass-spring-mass (MKM) system is shown in Figure 4.1. It consists of two masses connected by a spring

and damper. The control inputs are the force inputs u 1 and u2; the sensor outputs are the displacements Xl and x2;

and the disturbance input is w. The sensor outputs are corrupted by discrete sensor noise (not shown on figure) with

covariance v. Nominal values for the plant parameters, along with those for one known plant perturbation, are given

in Table 4.1. These parameter values result in a system with two poles at the origin associated with Xl, and two

lightly damped poles associated with x 2. The open-loop poles axe listed in Table 4.1. A M-File which defines the

system matrices for the MKM system is given in Appendix C.

The goal of our design is to increase the damping coefficient on the lightly damped poles to _=0.707 without

shifting their natural frequency, and to move the two poles at the origin so that the Xl response has a frequency of

0.825 rad/sec with a damping coefficient _--0.707.

,,_---. xl -------._

/

/

/

/

u1 = M

,t--- x2---_

k

W

/////////////////

Figure 4.1 Mass-Spring-Mass System

m

k

b

E{ wT(t-z)w(t)}

z{ vY(k+j)v(k)}

Open-Loop Poles

Table 4.1 Mass-Spring-Mass Parameter Values

Nominal Plant Perturbed Plant

M 1.0 kg 1.0 kg

0.1 kg 0.2 kg

0.01 N/m 0.01 N/m

0.01 N.s/m 0.01 N.s/m

0.18(,r)N 2 0.15('r)N 2

0.0016(j)m 2 0.0018(j')m 2

0, 0, -0.055+ 3.316i 0, 0,-0.030+_2.449i
22

4.1 COMPENSATOR DESIGN

Our compensator design followed the preceding steps. All computer inputs and outputs in the following

paragraphs are in this font. Bolded text should be input by the user, plain text is returned by the

computer, instructional comments are in italic.

Step 1: Select the weighting matrices for the quadratic cost function given in FAIn. (2.6) which characterize

the desired closed-loop performance

The cost function given in Eqn. (2.6) is the sum of the costs associated with each plant perturbation. In our case

we have two plants - the nominal and the perturbed plant. We determined the weighting matrices for our cost

function by designing two continuous time LQ regulators which placed the closed-loop poles for the nominal and

perturbed plants in the desired locations. The cost functions are

Jnominal = E{5.5x? + 2.2*2 2 + 10000u? + I0u22}/10.6 (4.1a)

Jperturbed = E{6.5x? + 4.8i 2 + 10000u12 + 10u2}/8.0 (4.1b)

where the costs have been scaled so that the LQR cost for each is unity. The cost function characterizing the

performance of our multirate compensator is the sum of the these two cost functions. By minimizing (,/nominal +

Jperturbed) the optimization routine will attempt to find a single compensator whose performance approaches that of

the LQR designs for both the nominal and the perturbed plant. Of course we do not expect such a compensator to

exist. Instead the resulting compensator will represent a compromise between good performance at the nominal

plant condition and good performance at the perturbed condition.

Step 2: Select an appropriate compensator structure and sampling schedule based on the desired closedloop performance

The closed-loop LQR system has a fast mode associated with Xl and Ul and a slow mode associated with x2

and u2. We selected a multirate compensator which capitalized on this structure. It consists of two coupled firstorder loops, one from x 1 to u 1 and another x 2 to u2 . A block diagram of this structure is shown in Figure 4.2. The

x 1 to Ul loop is sampled and updated every 0.8 seconds while the x2 to u2 loop is sampled and updated four time as

often, every 0.2 seconds. These sampling rates are approximately 10 time the desired closed-loop bandwith for each

loop. In addition, the compensator accounts for a one-half sampling period computational delay. It's sampling

schedule is shown in Figure 4.3.

Ul ! Mass-Spring-Mass U2 4th Order Plant

1st Order Compensator

State update = 0.8 sec

1st Order Compensator

State update = 0.2 sec

Figure 4.2 Block Diagram of Closed-Loop Multirate System

,'7 x1 Command

x2 Command
23

Xl ,91 X I I I I I_ Sampler

?l?r ?',,_? Samp.er

Ul I_; ', I ', ', ', ', ZOH

,,2 j?..?,,?,.?., zo.

Z-slow_l I t I I I I t State Update

_.,,?I ?'.?'.?I? StateUpdate

STP--O., .eo_ _-- L-- BTP--O.8.ec

= sample/update activity

Figure 4.3 Sampling Schedule for Multirate Compensator

Step 3: Design a compensator with the sampling schedule and structure selected in Step 2 that stabilizes all

plant perturbations. This will be the starting point for the optimization.

We designed a compensator using successive loop closures and root locus for the nominal plant. The

compensator has the same structure and sampling schedule as the multirate design discussed in Step 2 except there is

no input coupling and we did not account for computational delay. This compensator does however stabilize both

the nominal and perturbed plants even with the computational delay and so can be used as the starting point for the

optimization.

Refer to Reference 1 for an introduction to successive loop closures.

Step 4: Load the problem defmition into Matlab's workspace variables def'med in Table 3.2.

We defined the workspace variables using the Matlab script provided in Appendix D. To load the workspace

variables for this example, type

>>au:opt_mkm

Step 5: Initialize the workspace variables and generate the necessary global variable using mropt_init

Type

>>aucoI_t_4 nit see Section 3.3.2

At this point you can compute the minimum value of the cost function as follows (see Section 3.3.5)

>>calc_r.QQcoat: see Section 3.3.5

>> ans =

2. O108e+O0 this is the minimum cost our design could achieve
24

Step 6: Calculate optimum values of the digital processor gains using mropt_opfim

A partial output of the optimization follows

>>mropt_optim see Section 3.3.3

Automatically

Calculate gradient only?

selecting dscale a diagonal dscale is automatically selected if the variable

dscale does not exist, or is not compatible with the

current compensator

(y or n)n if we replied y the algorithm would compute only

the gradient of the cost with respect to the free

digital processor gains and the value of the cost

for the current compensator gains in crop

f-COUNT=

FUNCTION =

the total number of function evaluations

the value of the cost function at the current iteration. Refer to Matlab FMINU

documentation

f-COUNT

2

-> Plant 1

-> Plant 1

-> Plant 1

-> Plant 1

-> Plant 1

-> Plant

-> Plant

-> Plant

-9.9655e-01

-1.6399e-02

1.0030e+00

-I.0676e-02

-2.2977e-02

-9.4363e-01

-1.3703e-03

1.0529e+00

1.0030e+00

1.0036e+00

1.0529e+00

1.0085e+00

FUNCTION STEP-SIZE GRAD/SD

176 0.01 -9.48e+07

unstable while calculating J

unstable while calculating J

unstable while calculating J

unstable while calculating J

unstable while calculating J

1 unstable while calculating J

1 unstable while calculating J

1 unstable while calculating J

LINE-SEARCH

The algorithm encountered a destabilizing

solution. It automatically adjusts the step

size and continues the optimization

these are the values of the scaled digital processor gains

at the current iteration
25

313 4.62 14.2 2.45e-13

Optimization Terminated Successfully

Gradient less than options(2)

NO OF ITERATIONS=3!4

-8.0984e-01

-1.6700e+0i

1.8697e+00

2.2849e+01

7.4219e-01

-7.3479e+00

-8.6882e+00

7.9308e+00

1.0706e+00

3.8424e-03

1.0367e+00

1.4037e+00

int_st

**** Final Results  

Gain Gradient

-1.6197e-01 -4.6527e-06

-1.6700e+01 6.6547e-04

1.8697e+00 -2.8675e-04

2.2849e+01 4.5050e-04

7.4219e-0! !.6576e-04

-5.5109e+00 -8.9955e-04

-8.6882e+00 -4.8353e-05

7.9308e+00 -8.5171e-04

8.5648e-02 -i.0838e-04

1.9212e-03 6.9640e-04

6.2202e-01 3.4133e-03

1.4037e-01 3.2564e-04

these are the unscaled digital processor gain values

and their scaled gradient values

Final cost = 4.62482 _eopHmum value of_ecostfunc6on

Optimized compensator gains

cmpf =

1.9212e-03

0

8.5648e-02

0

0 1.8697e+00 -8.6882e+00

1.4037e-01 2.2849e+01 7.9308e+00

0 -1.6197e-0! 7.4219e-01

6.2202e-01 -1.6700e+01 -5.5109e+00

Elapsed time: 1516.05 sec
26

4.2 DESIGN ANALYSIS

We looked at two criteria when evaluating our multirate design: I) the final cost; and 2) the step response to a

command input.

4.2.1 Final Cost

The final cost for the multirate compensator is calculated automatically by the optimization routine

mropt_optim. The LQG cost is computed by the function calc_LQGcost in step 5. For our system

Jmultirate = 4.6 and JLQR = 2.0

4.2.2 Step Response

We obtained the response to a unit command step to x 1 as follows

>>[a,b, c, d] -_; get the nominalplant

>>b- [b zerom (4, 1) ] ; c= [c; [1 0 0 0] ] ; create a reference input

>>d-[d zeroB(4,1);O 0 -1];

>>ur- [zeros { 2,100) ; ones (1,100) ] ; generate a step input function for the input

and compute the time domain response, see Section 3.2.]

>>[y, t] =mrllim( [a b;c dJ,4,c_l)f,2,11u, sy, sz,SgppbtD, l, [I 2], [5 3],ur, lltp);

assuming cmpf, su, sy, sz, Stppbtp, and stp were previously defined and

that cmpf contains the optimized values of the digital processor gains

>>plot(t,y([1 3],:))

We similarly obtained the step response for the LQR design. The results are shown m Figure 4.4.

As expected the performance of the multirate design does not match that of the LQ regulator. The LQ regulator

was optimized for only the nominal plant. The multirate designs on the other hand represent a compromise,

stabilizing both plants and providing "optimum" performance for both. In addition, the LQ design was a continuous

full-state design whereas the multirate design was discrete second order and used only two state measurements.

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

; °,-" "°"-... °..

,-" ..0

2

_ i

D 0 "-.... ..... . ....... ......... .

i/ / --

I \...: .......... Multirate

.1 I J i i t 1 _ i

0 1 2 3 4 5 6 7 8 9

Time (Seconds)

Figure 4.4 Responseof MKM systemto command step responseto xI
.

.

.

.

.

27

References

Mason, G. S., "Multirate Compensator Synthesis and Analysis," Ph.D. Dissertation, University of Washington,

Seattle, Washington, 1992.

Mason, G. S., and Berg, M. C., "Reduced Order Mul_ate Controller Synthesis," AIAA Journal of Guidance,

Control and Dynamics, Vol. 15, No. 3, May-June. 1992, pp. 700-706.

Mason, G. S., and Berg, M. C., "Multirate Flutter Suppression System Design for a Model Wing," AIAA

Journal of Guidance, Control and Dynamics, Vol. 17, No. 6, Nov. -Dec. 1994, pp. 1267-1274.

Mason, G. S., Berg, M. C., and Yang, G. S., "A New Multirate Sample-Data Control Law Structure and Synthesis

Algorithm," AIAA Journal of Guidance, Control and Dynamics, Vol. 15, No. 5, Sep.-Oct. 1992, pp. 1183-1191.

Mason, G. S., Berg, M. C., and Mukhopadhyay, V., "'Multirate Flutter Suppression System Design for the

Benchmark Active Controls Technology Wing: Part I : Theory and Design Procedure, NASA TM 2002-

212128, December 2002.
28

Blank
29

APPENDIX A

GLOBAL WORKSPACE VARIABLES

The following global workspace variables are defined by mropt_init. Note that the Etrst letter of every global

variable is capitalized.

For each plant perturbation, the discrete closed-loop system can be written in the form:

xi(m,n+l) = Fix(m,n) + Giu(m,n) + Ww i (A.1)

yi(m,n) = Hix(m,n) + Vwi (A.2)

ui (re,n) = IS 1(n)P(n)S2(n) + S3(n)] Ys(re,n) (A. 3)

where

and w is a discrete zero mean white noise input with

E{wi(m, n)wT (m + j, n + k)} = {_(k) + 6(j)}/_"

(A.4)

(A.5)

The subscript i indicates the ith plant condition and the overbarred matrices are the discretized plant matrices. The

matrices S1, $2, and $3 are periodically time-varying switching matrices which model the multirate sampling

scheme, and P is a matrix of the compensator's digital processor gains.

The discrete cost function for the system in Eqns. (A. 1-A.3) can be written as

1 NpN-1 _frxi(m,n)qTFQli Mi]Fx/(m,n)]

J= lira NZZ /5_/ui(m.n)|L L "J LMiT Q2iJLui(rn'n)J (A.6) m--** * i=l n=O

See Reference 1 for more information.

The global Matlab workspace variables are then defined as follows.

Plant Description:

Np = number of plant conditions

MF = [F1, F2 ..... FNp ]

MG= [ G 1, G2 ..... GNp ]

MW= [ W 1, W 2..... WNp]

M]:I= [HI, H2 ..... HNp]

MV= IV1, V2 ..... VNp ]

MR= [R1, R 2 ..... RNp ]

NcF = number of columns in Fi

NeG= number of columns in Gi

NeW = number of columns in Wi

Nell= number of columns in Hi

NoV = number of columns in I_

NoR = number of columns in Ri
30

Compensator Description:

Nzhcmp = number of states in the periodically time-varying model of the multirate compensator

Stppbtp = the global version of stppbtp

Slk = [SI(0), SI(0) ..... SI(N-1)] Nsl = the number of columns in S1

S2k = [$2(0), $2(0) ..... S2(N-I)] Ns2 = the number of columns in $2

S3k = [$3(0), $3(0) ..... S3(N-1)] Ns3 = the number of columns in $3

Ptv: abs(Ptv) = number of columns in the compensator gain matrix [Az, Bz; Cz, Dz]. Ptv is positive if the

digital processor gains are time-invariant and negative if they are periodically time-varying.

Cost Function Description:

MQI= [Qll, Q12 .... Q1Np]

MQ2= [Q21, Q22 .... a2Np]

MM= [M1, M2 .... MNp ]

NeQI= number of columns in Qli

NeQ2= number of columns in Q2i

NcM= number of columns in M i

Miscellaneous Defmitions:

Ppn, Ppm, and Ppinit : used internally

Ppfree = index to the "free" digital processor gains

Fd = closed-loop state transitions matrix of the last computed ETIS system

Xe_global = value of the compensator gains at the last iteration

Unstable: If Unstable =0 the compensator gains at the current iteration stabilize the plants

Otherwise they destabilize the plant
31

APPENDIX B

SUPPORT M-FILES

B.1 MODELING SUPPORT M-FILES

B.I.1 mrspfit

Format: [a, b, c, d] = mrspfit(P, n)

Description: Separates the composite gain matrix P into its constituents, where

P:[::]

lnputs: P see F_xln.(B.1)

n the number of columns in a

Outputs: a, b, c,d the gain maWices in Eqn. (B. 1)

B.I_ mrmakesk

Format: [slk, nslk, s2k, ns2k, s3k, ns3k]=mrmakesk(su, s'y, sz, stppbtp)

Description:

Inputs:

Outputs:

B.1.3

Format:

Description:

Inputs:

Outputs:

03.1)

The periodically time-varying representation of a multirate compensator in Eqn. (2.2) can be written

as

Cz(n ) u(m,n) Fy<,,,.,,>I=I,:,<,,>c<,,>Tu<,,..,,>1 r.... ,<,,>],,<,,>+ s,<,,>}[=<,,,. ,,>] Lx(m,n+l)J LB(.)A(,,)J[_(,,,,,,)J--1"'t"'L_(,,)

where A z, Bz, Cz and D z are the distal processor gains and S1, $2 and $3 are switching matrices

which model the compensator's sampling schedule, mrmakesk creates the periodically timevarying switching matrices given the sampling schedule description.

su, sy, sz, stppbtp see Table 3.1

slk = [Sl(0) SI(1) ... Sl(stppbtp-1)]

s2k = [ff2(O) ff2(1) ... S2(stppbtp-1)]

s3k = [$3(0) 53(1) ... S3(stppbtp-1)]

nslk, ns2k and ns3k the number of columns in SI(-), $2(.) and $3(-) respectively

mrgetsk

[sl, s2, s3]=mrgetsk(k, slk, nslk, s2k, ns2k, s3k, ns3k)

Extracts the individual switching matrices for a given STP from slk, s2k and s3k

k an integer specifying for which STP the switching matrices are to be exlracted

slk nslk, s2k, ns2k, s3k, ns3k see Section B.1.2.

sl, s2, s3 switching matrices corresponding to S10t), S2(k) and S3(k) respectively. See

Section B. 1.2.
32

B.2 SYNTHESIS SUPPORT M-FILES

B.2.1

Format:

Description:

Inputs:

Comments:

calc_djdp

djdp= calc djdp(pguess, dscale)

Used by the routine mropt_optim to calculate the derivative of the cost J with respect to the free

compensator parameters. See Section 2.3.

pguess a vector of the values of the free digital processor gains at the current iteration

dscale see Table 3.3.

This routine uses the global workspace variables generated by mropt_init

B.2.2

Format:

Description:

Inputs:

Comments:

c_cj

djdp= calc_.j (pguess, dscale)

Used by the routine mropt_optim to calculate the cost J at the current iteration

pguess a vector of the values of the free digital processor gains at the current iteration

dscale see Table 3.3.

This routine uses the global workspace variables generated by mroptinit

B.2.3

Format:

Description:

Inputs:

Output:

Comments:

cale_L

[laee]= calc_L (pk, f, g, w, h, ql, q2, m)

Called by calc_djdp to calculate the steady-state values of a Lagrange multiplier. (The multiplier is

used to adjoin the stability equality constraints to the cost function.)

pk gains for the periodically time-varying representation of the multirate compensator

pk = [P(O), P(1) ..... P(N- 1)] and P(i) = FD(i) C(i)]

LB(i) A(i)J 03.3)

D(i), C(i), B(i), and A(i) are the state space matrices for the periodically time varying

representation of the compensator. See Section 2.2 and Section B. 1.2

t; g, w, h plant description matrices for the current plant perturbation corresponding to Fi, Gi, W and

Hi in Eqn. (A.4).

ql, q2, m cost weighting matrices for the current plant perturbation corresponding to Q1 i, Q2i and

Mi in Eqn. (A.6).

laec periodically time varying lagrange multipliers where

lacc = [A(1) A(2) ... A(N-1) A(O)]

and A(i) is the lagrange multiplier for the i th STP

This routine uses the global workspace variables generated by mropt_init
33

B.2.4 ealc_P

Format: [pk]= calc_P(p)

Description: Called by calc_.j and ealc_djdp to compute the periodically time-var3'ing representation of the

multirate compensator

Inputs: p the digital processor gains in the form

p = [p(O),p(1),...,p(N-1)] where p(i) =I Dz(i) Cz(i)l

LBz(i) az(i) j

Outputs: pk matrix of the compensator gains for the periodically time-varying representation of the multirate

compensator

pk = [e(1), e(2), P(3) .... P(N)] where P(i)= Sl(_)p(i)S2(i)+ S3(i)

See Eqn. (A.2) or Section B. 1.2.

Comments: This routine uses the global workspace variables generated by mropt_init

B.2.5 caic_X

Format: [xaec]= calc X(pk, f, g, w, h, v, r)

Description: Called by calc_djdp and calc_j to calculate the steady-state covariance values of the closed-loop

system

Inputs: pk matrix of the compensator gains for the periodically time-varying representation of the multirate

compensator, see Section B.2.4.

l; g, w, h plant description matrices for the current plant perturbation corresponding to Fi, Gi, W, H i

and V in Eqn. (A.4) respectively

r process and sampling noise covariances for the current plant perturbation corresponding to Ri in

Eqn. (A.5).

Comments: This routine uses the global workspace variables generated by mropt_init

B.2.6 disc_cost

[qld, q2d, md]=disc_cost (pit, nplt, stp, ql, q2, m)

Computes the discrete cost function weighting matrices such that the cost associated with a

continuous plant in feedback with a single-rate compensator is identical to the cost associated

with the discretized plant in feedback with the single-rate compensator. The continuous cost

function is given in Eqn. (2.6); the discrete cost function is given by

1 _S_I {rx(m,n)lTrQld MdlF_(_,")II

Jd#crete m--r- n=O l [ufm'n)3 LM_ Q2dJLufm'n)JJ

= I.Lm _. _... £ 03.4)

disc_cost assumes the control input is updated using a zero-order-hold.

pit, nplt, stp see Table 3.1

ql, q2, m the continuous cost function weighting matrices in Eqn. (2.6)

qld, q2d, md the discrete weighting matrices in Eqn. (13.4)

Format:

Description:

[rtputs :

Outputs:
34

B.2.7

Format:

Description:

Inputs:

Outputs :

disc_noise

[w]=disc_noise(wxx, pit, nplt, stp)

Computes the covariance of a discrete-time white noise process disrarbance such that the expected

value of the states of a continuous plant, and the expected value of the states of the discretized

plant are equivalent at the sampling instances. Thus, given the continuous plant

x(t) = Ax(t)+ Gww(t), where E(w(t)wT(t + r)) = ¢3('r)wxx

and the plant discretized with a zero-order-hold at a sampling period T

xd(m, n + 1) = Adxd(m, n) + wd(m, n), where E(wd(m, n)wT(m + k, n + j)) = (6(k) + c$(j))w

then E{xd(m, n)x_(m, n)} = E{x((mN + n)T)xT((mN + n)T)}

wxx PSD of the continuous processor noise

pit, nplt, stp see Table 3.1. The inputs to pit are assumed to be only the process noise

w PSD of the discrete process noise

B.2.8 dlyap2

Format: x= dlyap2 (a, c, method )

Description: Solves the discrete lyapunov equation

x= a-x-aT + c

Inputs:

Outputs:

a, c system matrices in Eqn. (B.5)

if method = 'Bartels' then the algorithm solves Eqn. (I3.5) using Bartels method

otherwise Eqn. 03.5) is solved using partial sums

x the solution to the lyapunov equation

03.5)

B.2.9

Format:

Description:

Inputs:

Outputs:

Comments:

get_cost

[ql, q2, m]= get_cost(i)

Extracts the cost function weighting matrices Q1, Q2 and M for the i th plant perturbation from the

global workspace variables MQ1, MQ2, and MM. See Appendix A.

i an index to the desired plant perturbation

ql, q2, m the cost weighting matrices corresponding to Q1, Q2 and M respectively

This routine uses the global workspace variables generated by mropt_init
35

B.2.10

Format:

Description:

Inputs:

Outputs :

Comments:

get_pit

[f, g, w, h, v, r]=get_plt(i)

E,..a,.,.. ..... the i., Fi ' t-z_, W, Hi , _z Ri _:,.r Lhe ;th ,.!o,,, p._,,.r, ot;,.,, ¢,-,.m th,, glnh.1

workspace variables MF, MG, MW, MH, MV and MR. See Appendix A.

i an index to the desired plant perturbation

f, g, w, h, v, r plant description matrices for the ith plant perturbation corresponding to Fi, Gi, W, Hi,

V, and Ri respectively

This routine uses the global workspace variables generated by mropt_init

B.2.11

Format:

Description:

Inputs:

Outputs:

get_._ppfree

ppfree = get__ppfree (cmpfree, ncmp)

Called by mropt__optim to determine which compensator gains will be optimized

empfree, nemp see Table 3.1.

ppfree vector whose elements point to the flee elements in crop

B.2.12

Format:

Description:

Inputs:

Outputs:

Comments:

get__sk

[sl, s2, s3]=get_sk(k)

Extracts the switching matrices Sl(k), S2(k) and S3(k) for the k th STP from the global workspace

variables Slk, S2k, and S3k. See Appendix A.

k indicates the k th STP of the current BTP. NOTE: k = 1, 2..... stppbtp

sl, s2, s3 switching matrices for the kth STP corresponding to Sl(k), S2(k) and S3(k) respectively

This routine uses the global workspace variables generated by mropt_init

B.Z13

Format:

Description:

Inputs:

Outputs:

make_noise

[r]= make_noise (rw, Iv)

Creates a compound noise covariance matrix from rw and rv

rw, rv the process and sensor noise covariance matrices respectively

r compound noise matrix of the form
36

B.2.14 makeplt

Format:

Description:

Inputs:

Output:

B.2.15

Format:

Description:

If, g, w, h, v, nzcmp ]= make_pit (pltc, nplt, ncmp, stp)

.... t,._ the compound plant _•_,a'_° .... s Fi, Gi, W and V for the current plant perturbation. See

Appendix A.

pltc has the same form as pit in Table 3.1 except it describes only a single plant perturbation

nplt, ncmp, stp See Table 3.1.

f, g, w, h, v compound plant matrices Fi, Gi, W and V for the current plant perturbation respectively

nzcmp the number states in the periodically time-varying representation of the multirate

compensator

mropt_check

mropt_check

mropt_check performs elementary error checking on the data in Table 3.2 and set the variable

err=-i if it finds an error in the data. mropt_check is called by mroptinit.

None

None

B.2.16

Format:

Description:

Inputs:

Outputs:

mropt_global

mropt_global

Called by mropt__init to define the global workspace variables

None

None

B.2.17

Format:

Description:

Inputs:

Outputs :

Comments:

mropt_fminu

Ix, OPTIONS] = mropt_fminu(FUN, x, OPTIONS, GRADFUN, P1, P2, P3, P4, 1'5....

P6, PT, PS, P9, P10)

A modified version of FMINU fi-om the Optimization_Toolbox. This version automatically reduces

the search step size when a destabilizing solution is encountered.

See FMINU

See FMIN-U

Requires the Optimization_Toolbox
37

APPENDIX C

M-FILE MKM

The following M-File creates the state space matdccs for the Mass-Spring-Mass system of Section 4.0

function [a,b,c,d]=mkm(ml,m2,k,b)

% [a,b, c, d] =mkm (ml, m2, k,b)

if nargin-=4

ml=l;

m2=.l;

k=l;

b=O.Ol;

end;

%nominal plant

t=(i/ml)+(i/m2);

a=[O I 0 0

0 0 k/ml b/ml

0 0 0 1

0 0 -k*t -b*t];

b-_-- [0 0

llml 0

0 0

-llml llm2];

c=eye (4) ;

d=zeros (4,2) ;
38

Blank
39

APPENDIX D

SCRIPT MROPT_MKM

The following script defines the workspace variables in Table 3.2 for the example problem of Section 4.0

% Script: mropt_mkm.m

% Example script for creating input data for mropt

%

% Creates the following data required by mropt_init

% plt,nplt,wxx,rv,qla,q2a,ma

% ucol,ncol,srow, crow

% sz,su,sy, cmp,stp,Stppbtp,cmpfree

% G. Mason, U.W. 1992

Continuous plant description

Np=2;

nplt=4;

plt= [] ;

% The number of plant conditions

% The number of states in the plants

% Clear the variable which holds the plant data

% Plant #I

[a,b,c,d]=mkm(l,0.1,1,0.01) ;

plt = [a b; c d] ;

% Plant #2

[a,b, c, d] =mkm(l, 0.2, I, 0.01) ;

plt = [plt [a b;c d] ];

% Pointer to rows and columns in the plant

ucol=[l 2]; % control input

ncol=2; % process noise input

srow=[l 3]; % sensor output from the plant

crow=-[l 2 3 4]; % criterion output from the plant

% Continuous process noise PSD

% wxx is a compound matrix like plt, it contains wxx for plant #i and #2

wxx=[.l .i];

% Discrete sensor noise PSD

% with same compound form as plt

rv=eye(2)*0.001;

rv= [rv rv] ;

% ..... Continuous Cost weighting matrices

% Output criterion weights

% Use le-8 instead of 0 so that qal & qa2 are positive definite

% The synthesis algorithm would accept semi-definite qla & q2a

% but calc_LQGcost will not

qlal=diag([5.5 le-8 le-8 2.2]);

qla2=diag([6.5 le-8 le-8 4.8]);

qla=[qal/10.6 qa2/8.0];

% Control input weights

q2a=diag([le4 lel]);
4O

q2a=[ra/10.6 ra/8.0] ;

%Cross weighting (ma)

%if there is no cross weighting it can be left undefined

clear ma

% ..... Compensator description .....

% Digital processor gains from successive loop closures design

cmp= [ 0.5 0 1 0

0 0.I 0 1

0.08 0 -0.2 0

0 0.6 0 -0.75] ;

% Sampling schedule

su=[0 0;I i;0 0;0 I;0 0;0 I;0 0;0 i];

% Compensator output updating w/ delay

sy=[8 2]; % Sensor input sampling, multiplexed

sz=[8 2]; % Compensator state update

stp = .i; % The shortest sampling period

Stppbtp = 8; % The number of stp's in one BTP

% Free compensator gains

cmpfree=[l 0 1 1

0 1 1 1

1 0 1 1

0 1 1 i];

% clear temporary variables

clear a b c d
,m

Form Approved

REPORT DOCUMENTATION PAGE OMB No. 0704-0188

The publ c reporting burden Ior this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existingdata sources,

gathering and maintaining the data needed, and completing and reviewing the collection ol information. Send comments regarding this burden estimate or any other aspect of this

collect on o n ormation nclud ng suggestions lor reducing this burden, to Department of Delense Washington Headquarters Services, Directorate for Information Operations and

Reports (0704-0188) 1215 Jellerson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision ollaw, no person

shall be subiect to any penalty for failingto comply with a collection Of inlormation if it does not display a currently valid OM8 control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY) 2. REPORT "_fPE

12-2002 Technical Memorandum

4. TITLE AND SUBTITLE

Multirate Flutter Suppression System Design for the Benchmark Active

Controls Technology Wing

Part 11." Methodology Application Software Toolbox

6. AUTHOR(S)

Mason, Gregory S.; Berg, Martin C.; and Mukhopadhyay, V,

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

3. DATES COVERED (From - To)

5a. CONTRACT NUMBER

5b. GRANT NUMBER

iSc. PROGRAI_I ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

706-17-51-03

8. PERFORMING ORGANIZATION

REPORT NUMBER

L-18249

: 101 SPONSOR/MONITOR'S ACRONYM(S)

NASA

11. SPONSOR/MONITOR'S REI:;ORT

NUMBER(S)

NASA/TM-2002-212129

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 63

Availability: NASA CASI (301) 621-0390 Distribution: Nonstandard

13. SUPPLEMENTARY NOTES

An electronic version can be found at http://techreports.larc.nasa.gov/ltrs/or http://techreports.larc.nasa.gov/cgi-bin/NTRS

14. ABSTRACT

To study the effectiveness of various control system design methodologies, the NASA Langley Research Center initiated the

Benchmark Active Controls Project. In this project, the various methodologies were applied to design a flutter suppression

system for the Benchmark Active Controls Technology (BACT) Wing. This report describes the user's manual and software

toolbox developed at the University of Washington to design a multirate flutter suppression control law for the BACT wing.

15. SUBJECT TERMS

Cybernetics; Digital Control Systems; Multirate sampling; Flutter suppression; Control law design and optimization; Metlab

toolbox

16. SECURITY CLASSIFICATION OF:

a. REPORT b. ABSTRACT c. THIS PAGE

U U U

17. LIMITATION OF

ABSTRACT

UU

18. NUMBER

OF

PAGES

46

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa._ov)

19b. TELEPHONE NUMBER (Includo area code)

(301) 621-0390

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Sld. Z39,18