TH E U N I V E RS I T Y OF MICHIGAN COLLEGE OF ENGINEERING Computer, Information and Control Engineering Technical Report DESIGN OF DECOUPLED MULTIVARIABLE CONTROL SYSTEMS John R. Pivnichny E. G. Gilbert Project Director ORA Project 026450 supported by: UNITED STATES AIR FORCE AIR FORCE OFFICE OF SCIENTIFIC RESEARCH GRANT NO. AFOSR-69-1767 ARLINGTON, VIRGINIA administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR March 1971

Ackn owle dgement s The author wishes to express his sincere gratitude and appreciation to Professor Elmer G. Gilbert for his constant guidance and helpful suggestions over the period of time in which I was working on the research reported here and during the preparation of the dissertation. The other committee members are also to be thanked for their suggestions and willingness to carefully review the work during the busiest time of the year. Finally, I should like to thank Mrs. Maureen Konchar for an excellent job of typing the final draft. ii

TABLE OF CONTENTS Page Acknowledgements ii List of Tables vi List of Figures vii List of Appendices viii Chapter 1 Introduction 1 1.1 Introduction to Problem 1 1.2 Outline of Dissertation 8 Chapter 2 Previous Work 10 2.1 Morgan's Control Law 10 2.2 Necessary and Sufficient Condition for Decoupling 12 2.3 Weakly Coupled Systems 13 2.4 Decoupling Theorem 16 2.5 Fixed Poles 19 2.6 Summary of Results Available 21 Chapter 3 Theoretical Results 22 3.1 Control Effort and Error Performance Procedure 22 3.1.1 Problem Approach 23 3.1.2 Input Selection 26 3.1.3 Cost Functional 28 3.1.4 Form of Gradient and Hessian 30 3.1.5 Algorithm 33 3.1.6 Independent Inputs 40 iii

3.2 Disturbance Inputs 44 3.2.1 Complete Cancellation Results 3.2.2 Inclusion of Disturbance Inputs in Minimization Procedure 58 3.3 Sensitivity 61 3.3.1 Complete Cancellation of Parameter Variation Effects 61 3.3.2 Calculation of Static Sensitivity and its Gradient 65 3.3.3 Integral Control for Static Sensitivity Improvement 73 Chapter 4 Design Procedure 75 4.1 Design Data 75 4.2 Reduction to State Variable Feedback Problem 78 4.3 Application of Optimization Procedure 79 4.4 Investigation of Disturbances 82 4.5 Check on Parameter Variations 84 4.6 Limitations of Design Procedure 87 Chapter 5 Computational Questions 88 5.1 Input Data - Decoupling Program Results 89 5.2 Numerical Methods for Solution of Lyapunov Equation 92 5.3 Run Times 96 5.4 Core Storage Required and Numerical Accuracy 99 5.5 An Efficient Program 102 5.6 Maximum Problem Size 104 Chapter 6 Example Problems 106 6.1 Motor Generator Example 106 6.2 Aircraft Example 114 6.3 Distillation Column Example 119 6.4 Helicopter Example 130 iv

Chapter 7 Conclusion 137:1X' _reC e 11 0 Appendices 146 Table of Nomenclature 159

LIST OF TABLES Table Page 6.1 Effect of Weighting Parameters on Error and Control Variances, Motor Generator Example 112 6.2 Effect of Filter Bandwidth on Error and Control Variances, Motor Generator Example 112 6.3 Designs Using Second Order Filters, Motor Generator Example 113 6.4 Distillation Column A and B Matrices 121 6.5 Decoupling F and G Matrices for Distillation Column 12 2 6.6 Perturbed A and B Matrices for Distillation Column 126 6.7 Perturbed Distillation Column Variances 127 6.8 Helicopter A and B Matrices 131 6.9 Decoupling F and G Matrices for Helicopter 134 6.10 Comparison of Pole Locations of Closed-Loop Designs 136 vi

LIST OF FIGURES Figure Page 3.1 Steps of Algorithm 35 4.1 Compensator 80 6.1 Motor Generator 107 6.2 Distillation Column 120 vii

LIST OF APPENDICES Appendix Page A - Difficulties of Quadratic Approach 146 B - Modified W Matrix 148 C - First Order Uniqueness Result 150 D - Details of Result V 152 E - Application to Sample Data Systems 154 F - Canonically Decoupled Form 156 viii

CHAPTER 1 INTRODUCTI ON The first section of this introductory chapter describes in general terms the main problem to be considered. Following that is a description of the arrangement of material in the main chapters. 1.1 Introduction to Problem This dissertation is concerned with the design of decoupled multivariable feedback control systems. Multivariable control systems involve multivariable plants, that is, systems having several inputs and several outputs. In general there is coupling within the plant so that one plant input, or control, affects more than just one output. The reason for applying feedback to the plant is to improve the quality of the system response, and in the case of decoupling to eliminate coupling between the closed loop inputs and plant outputs. A large number of practical systems are multivariable in nature, including aircraft and space vehicles, chemical processes, biological systems and electric power generating systems. Consequently, there is a need for techniques to design control systems for these multivariable systems. The classical, frequency domain approach is directed more toward analysis of single input-single output systems and is not well suited to multivariable systems with coupling. Therefore, the modern state variable approach has been used.

The systems mentioned are in general non-linear multivariable systems and a few results are available for such systems, principally by Iwai [Ill, and Nazar and Rekasius N1l]. The greater part of the work to date has assumed that these systems can be described sufficiently well in a small neighborhood of an "operating point" by a set of linear constant coefficient differential (or difference) equations. The same assumption is made throughout this dissertation including the remainder of this section. The modern approach to decoupling was first introduced in 1964 by Morgan [M51] who formulated the problem of decoupling by state feedback, i. e., a control strategy consisting of linear feedback of the system state variables combined with a linear combination of the closed loop system inputs. He was able to define a class of systems which could be decoupled into independent first order subsystems with this control strategy. These results were extended somewhat by Rekasius [R1]. Later, Falb and Wolovich [F1, F21 developed a necessary and sufficient condition that a system can be decoupled and a procedure for choosing a control law which decouples. However, their procedure did not yield all possible closed loop designs. More recently Gilbert [G2] developed complete specifications for the class of control laws which decouple, the class of decoupled closed loop systems, and the correspondence between elements of these two classes. Some

additional results involving pole assignability and system observability by Slivinsky and Schultz ES2] and Mufti [M71 respectively have appeared since 1968. However, these are rather obvious extensions of Gilbert's work. Two limitations of decoupling with Morgan's control law which have been mentioned by Gilbert and others, will now be discussed. The first limitation involves those multivariable systems which cannot be decoupled by Morgan's control law, but which can be decoupled by control laws that include dynamic elements. Such systems are said to have "weak inherent coupling." The second limitation follows from the decoupling requirements itself. As it turns out, some systems can be decoupled only by fixing certain closed loop poles. These poles may have undesirable locations in the complex plane, causing instability or other objectionable effects. Silverman rS1 has given a design procedure for weak inherently coupled systems which is based on his work on inverse systems. It results in a control law which includes dynamic elements and which allows decoupling. Whether or not his dynamic compensator has minimal order is an unanswered question at this time. Wonham and Morse [W3, M6J on the other hand, have results which produce a minimal order dynamic compensator. At the present time their geometric approach has not been translated into a constructive design procedure. The most important result is the following. If a multivariable system has weak inherent coupling, then a

dynamic compensator can be constructed so that the system plus compensator can be decoupled by state variable feedback alone. In addition, Wonham and Morse show that the fixed pole locations can be shifted by additional dynamic compensation and give a way of calculating the minimal order of compensator needed to shift any prescribed set of fixed poles. The compensator poles can also be freely assigned. Howze and Pearson [Hi] also have a procedure for calculating a relatively low order compensator for free assignment of poles. However their results are incomplete; they do not apply to all systems. It should be emphasized that all of the preceding remarks apply to the case where the number of plant inputs is equal to the number of plant outputs. When the number of plant inputs exceeds the number of outputs, much more difficult questions are involved and are still not totally resolved [W3, M6, S1]. In this dissertation it is assumed that the number of plant inputs and outputs are equal. All of the results mentioned above involve a considerable amount of computation so that for systems of moderate or high order (greater than three) hand calculations are are impractical. Thus a good part of the work in this area is concerned with computer programs which make design of such large systems feasible. An example is the recent paper of Gilbert and Pivnichny [G3].

The end result of all this work is that all cases which arise can be reduced to the basic problem of decoupling a plant by state variable feedback. It is not clear at the present time what the best way is to do this reduction, and other researchers are currently working on the question. For our purposes it is sufficient to indicate that reduction is possible, because our concern is with the basic decoupling design problem itself. More specific statements and the currently available details are given in the next chapter. There are many factors which must be examined in order to achieve an acceptable practical design. The following are the most important: 1. The error between outputs and inputs should be small. 2. Control magnitudes must not be too large. 3. Effects of disturbances inputs on the closed loop system behavior must be kept small. 4. Parameter variation effects must also be small, particularly changes in static gains should be very small. 5. Static and dynamic cross-coupling due to parameter variations should be very small. At the present time there is no technique which addresses all of these factors perfectly. The aim of this

dissertation is to present a reasonable way of approaching the design of high order multivariable decoupled systems. A number of complementing procedures are presented which together can be used to achieve satisfactory closed loop performance. The most important of these involves selecting a class of random inputs and solving for those decoupling control law parameters which minimize a cost functional involving output error variances and control variances. This technique has the advantage of simplicity along with the freedom to adjust the weight coefficients in the cost functional so as to achieve desired compromises between the tracking errors and control efforts. Another advantage is that is is possible to include disturbance inputs and measurement errors and minimize their effect on output error. The resulting parameter minimization problem is non-linear in nature and at this time there are no general existence or uniqueness results for its solution. Because of the complexity of the minimization problem, a computational algorithm has been developed which searches numerically for a solution. Other techniques which might be considered for addressing the control amplitude aspect of the design problem are discussed in a later chapter. Results concerning the parameter variation and disturbance input questions are also presented. The main value of any synthesis technique, of course, is its usefulness in the design of practical control systems. For this reason a number of examples have

been worked out to demonstrate the practicality of the proposed techniques. Several of the examples are of reasonable complexity.

1.2 Outline of Dissertation The contents of this dissertation are arranged in the following way. First a precise definition of the decoupling problem is given in Chapter 2. Here the main results which have appeared and are necessary to the understanding of what follows are presented. Some results for weakly coupled systems and alteration of the fixed poles are also included. All these results and the results in the subsequent chapters can be applied to sample data systems. This is discussed briefly in Appendix E. In Chapter 3 the objectives and criteria for a design synthesis are examined. Theoretical results involving the selection of a decoupling feedback law which minimizes a quadratically weighted sum of output error and control effort appear here. The basic steps of an algorithm for minimization are described. Some invariance results which define the class of disturbance inputs and parameter variations which affect only the unobservable part of a decoupled system, and therefore have no effect on the system output are presented. The necessary restrictions on feedback parameters which make this invariance possible are also given. In general the disturbance inputs and plant parameter variations cannot be cancelled completely however. It is shown in this chapter how the effects of such disturbance inputs can be included in the parameter minimization procedure mentioned above. Results concerning the effect of parameter variations which cannot be cancelled completely, on the static

gain of the closed-loop system also appear in Chapter 3. A steady state sensitivity matrix is defined, justification given for its use, and a method of calculating changes in static sensitivity as feedback parameters are varied is given. The question of how to put all these results together into a useful design procedure is answered in Chapter 4. Questions about selection of design requirements, noise characteristics, etc. are investigated. Weak areas in the overall design procedure are pointed out. Chapter 5 contains details concerning the computational aspects of design synthesis. It is shown how the previous work of Gilbert and Pivnichny [G3] forms a basis for the computation required by the design technique of Chapter 4. Important computer programming questions are investigated. Numerical results for several example problems appear in Chapter 6. Here the results of investigation concerning problem dimensionality weights on control and error, input noise characteristics, and computational speed are presented.

CHAPTER 2 PREVIOUS WORK The results which we shall consider apply to those multivariable plants whose operation can be described sufficiently well in a small region about an "operating point" by a set of linear constant-coefficient differential equations of the form: i(t) = Ax(t) + Bu(t)'2.1) y(t) = Cx(t), where t is time; the state x is an n-vector; the control u and output y are both m-vectors. The dot denotes a time derivative, and A, B, and C are constant matrices of appropriate size. The principal design constraint is that the resulting closed-loop system be decoupled. Consequently in this chapter we will review the pertinent decoupling results which are known today. The main objective is to provide a foundation for what follows in the subsequent chapters. Readers who are interested in the details about the historical development of the results should consult the following references [MS, R1, Fl, F2, G2, M6, S2, W31. 2.1 Morgan's Control Law In 1964 Morgan [MS] first proposed that the plant (2.1) be decoupled by means of a state variable feedback control law u(t) = Fx(t) + Gv(t) (2.2) 10

11 where the closed-loop system input v(t) is an m-vector and F and G are constant matrices of sizes m x n and m x m respectively. He stated the following: Definition 2.1 The closed-loop system consisting of plant (2.1) with control law (2.2) x(t) = (A + BF) x(t) + BGv(t) (2.3) y(t) = Cx(t) is dec'oupled if the resulting closed-loop matrix transfer function H(s, F, G) = C(sI - A - BF)-1BG (2.4) is diagonal and non-singular. Whether or not there exists a matrix pair {F, G} for a given plant (2.1) so that control law (2.2) can achieve decoupling is determined from the following result.

12 2.2 Necessary and' Sufficient Condition for Decoupling Let Ci be the i-th row of C, and define the row vector D. D = C AdiB i i and O. if C.B 0 d. =-^ 1 (2.5) j, if C~B - 0 where j is the largest integer from (1,..., n-l) such that C.AkB = 0 for k = 0,..., j-1. Now form the m x m matrix D D D 1 (2.6) then a necessary and sufficient condition for decoupling is given by: Theorem 2.1 A given system S = {A, B, C}can be decoupled by state variable feedback (2.2) if and only if D is non-singular. Proof of Theorem 2.1 is given in [F1] and [G2].

13 2.3 Weakly Coupled Systems The question arises of what to do when the conditions of Theorem 2.1 are not satisfied. In that case we may ask whether some other control law will achieve decoupling. It is obvious that if the open loop matrix transfer function H(s) = C(sI - A) 1B (2.7) is singular, for all s, then no control law can effect decoupling. On the other hand, Wonham and Morse EW3J have shown that those plants for which H(s) is non-singular but have det D = 0 can always be decoupled by the addition of suitable dynamic compensation. Such systems are said to have weak inherent coupling. Consider, for example, a compensator of the form x = A x (t) + B ui (t) (2.8) where x is an n-vector, u (t) is an m-vector, A and B are constant n x n and n x m matrices respectively. If this compensator is interconnected with the plant (2.1) by u (t) - K1 ut) )) K2 t) (2.9) where K1 is m x m, K2 is m x n, then the overall system can be viewed as a new system S = {A, B, C} of order n = n + n with state x -, input u = u and output y = y with A A f BK21 [~l =C A- 2 B - C =[CC 0J

14 We remark that although the new control u is different from the original u, this fact does not cause any great difficulty in design. A detailed discussion appears in Section 4.3. Wonham and Morse's result can now be stated as: Theorem 2.2 Given a plant (2.1) with weak inherent coupling, there exists a minimal order compensator of form (2.8) for which the combined system S can be decoupled by state variable feedback. It is clear then that the presence of weak inherent coupling in a plant need not deter us from satisfying the decoupling requirement of design, provided of course we are willing to add the necessary dynamics. Reduction of Wonham's geometrical results to a computational algorithm for the synthesis of dynamic compensators has not appeared as of yet. An alternative approach, however, is to utilize Silverman's [S1] algorithm for calculating the inverse of a dynamic system (which exists if H(s) is n. s.) [S1]. The way in which this algorithm is used is as follows. A necessary and sufficient condition [S1 for H(s) to be non-singular is the existence of a positive integer a < n, a non-singular differential operator N of degree a, a non-singular constant matrix D and a constant matrix C such that Ny(t) = C x(t) + D u(t) (2.10) a a

15 Silverman's algorithm is a procedure for constructing N, Ca, Da. If these have been determinedI, we let s replace the differential operator d/dt and form N(s) from N. Now let B. = max k such that sk is an element of the j-th column of N(s) and define a non-singular matrix differential operator A corresponding to the matrix function A(s). A(s) = diag (s-Bj) (2.11) Then defining a control law -1 -1 u(t) = D NAu(t) - D 1C x(t) (2.12) a a a we see from (2.10) that Ny(t) = NAu(t) (2.13) as for zero initial condition y(t) A= u(t) (2.14) which is a decoupled closed-loop system, so the conditions of Theorem 2.1 must be satisfied for this system (det D O0). The dynamic compensator NA which does not contain any differentiators can be synthesized i. e., find A, B, K1, K2 such that u(t) = NAu(t), by standard procedures described in EK1] and elsewhere. A disadvantage of the tchnique is that it is unknown at this time whether it yield-!' a minimal order compensator.

16 2.4 Decoupling Theorem Once it has been determined that det D ~ 0 or the original plant (2.1) has been suitably compensated so that the combined system has det D 0, then we may proceed to examine the main decoupling results which are due to Gilbert [G2] and are presented here for use in later parts of the thesis. First of all, an m x n matrix A* is defined Adl + 1 A*= [ 1 (2.15) Cm dm + 1 then we state: Theorem 2.3 If D is non-singular, the following data are uniquely determined from S = {A, B, C}: integers Pi > 0, i = 1,..., m; integers ri > 0, i 1,..., m + 2; polynomials ai (s) sri- ailsri-l iri i = 1,..., m + 2 (if ri = 0, ai. (s) = 1); m x m matrices Gi, i = 1,..., m; 1 m x n matrices Jik' i =1,..., m, k 1,..., pi; m x n matrices Kik, i 1,..., m, k 1, rm + 2 (the Kik are not defined if rm + 2 =0) The class of {F, GI which decouple is given by G - m x. G. i 1 1 and (2.16)

17 p. m rm+2 F -D1 A* + E (ik Tik) Jik + Pik Kik i=1 k=l i=1 k=l where nik ik' k 1,..., ri, ik = k r+l, ik 1ik 1 k 1 Pi and A., 0 Pik are arbitrary real numbers ( X. 0, i = 1,.., m). Let hi (s, F, G) denote the elements of the decoupled (diagonal) closed-loop transfer function matrix H (s, F, G). Then.i a. (s) hi (s, F, G) = 1 (2.17) 1i (s, i) Pi Pi (s, a.) s S i a ip. a. = (.,. il) where.i and the elements of a i may be chosen arbitrarily. The characteristic equation of the closedloop system is: m q (s,F)=det(sI-A-BF)= a+m+2 (s) T i(sa) (2.18) m;2 i-1=l The actual details of how these equations arise and a proof of the results in the decoupling theorem can be found in Gilbert's work [G21]. We can see from (2.17) all possible decoupled subsystem characteristics. Notice that all of the subsystem

18 poles can be freely assigned, but the zeros are fixed. A particular characteristic can be specified by selection of the real numbers., i 1,..., m and.ik i 1,..., m, k = 1,''' Pi' Once these numbers are chosen, then the corresponding control law which produces that closed-loop characteristic can be calculated by equation (2.16). If the open loop system is controllable, which we shall assume in Chapter 3 for simplicity, then r = 0 and the matrix m+ 2 pair {F, G}is uniquely determined. Observe from equation (2.18) that unless rm + = 0 m+1 and rm + 2 = 0 the roots of a + 1 (s) = 0 and the roots of am + 2 (s) = 0 appear as fixed roots of the closed-loop characteristic equation. In that case we may well ask what should be done if any of these roots lie in the right half plane. Clearly it would be impossible to design a stable decoupled closed-loop system with state variable feedback alone and another approach such as the one to be described in Section 2.5 is needed. Determination of all the integers, polynomial coefficients, polynomial roots, and matrix elements in the synthesis theorem above would be impossibly difficult without the aid of a computer program such as the one developed by Gilbert and Pivnichny EG3]. This program calculates all this data from the system matrices A, B, and C and as such forms a basis for the synthesis procedure which is described in Chapter 3.

19 2.5 Fixed Poles In the preceding section, the possibility of undesirable fixed poles appearing as a result of the decoupling requirement was noted. As mentioned in EG3] this phenomena occurs in a large percentage of the sample problems considered and is a real concern in any practical case. Fortunately there are results available which allow us to circumvent the difficulty of undesirable fixed poles. For it is shown in [W3] that the location of the fixed roots can always be shifted by dynamic compensation of the type considered in Section 2.3. Wonham and Morse EW3] have developed the most general results. Their solution is stated in the framework of first separating the fixed poles into those with "good" locations and those with "bad" locations, then choosing a minimal order compensator which simultaneously (1) allows decoupling and (2) results in a decoupled system for which the only fixed poles are those with the "good" locations. Unfortunately, their geometrical results lack computer implementation at this time. Alternative, though less general results are given in Howze and Pearson [H11. Their generality can be seen from the main theorems which are stated here. Define v. as the greatest integer j such that the vectors Ai Bi, j = 0,..., vi - 1 are linearly independent. Theorem 2.4 Assume (A, B) is controllable, S = {A, B, C} can be decoupled by {F, G} and

20 i Pi + p for i = 1,..., m. It is then i m + 1 possible to compensate S by means of a dynamic compensator of order k = (m-l) p such that n + k poles can be placed arbitrarily with state feedback. Theorem 2.5 Let the assumptions of Theorem 2.4 hold except v. < p. + p for one or more i. 1 i m + 1 Then if the fixed roots are all real and distinct, it is possible to compensate S by means of a compensator of order k j v. - n such that n + k poles can be i=1 arbitrarily placed.

21 2.6 Summary of Results Available To summarize the material presented in this chapter we may say that results are available which specify a decoupling feedback law for any plant of the type (2.1) for which det H(s) 3 0. This specification involves using state variable feedback primarily with the addition of only enough dynamic compensation to allow decoupling and assignability of undesirable "fixed" poles. Although all the details involved in selecting suitable dynamic compensation have not been taken care of completely at the present time, and only a sketch is given of available results, it is clear that compensators do exist which allow all cases to be considered as state variable feedback problems. We have consequently arrived at the starting point for the original work in this thesis. It will be assumed in the following that the necessary compensation has been determined, and the task at hand is to choose a decoupling control law from the class (2.16) (by suitable choice of aij i l,...,m; j = 1,...,Pi) which produces a "best" overall design. Specification of what constitutes a best design and techniques for computing it will now be presented.

CHAPTER 3 THEORETI CAL RESULTS This chapter is divided into three main sections. The first deals with the random input technique for design of decoupled systems. A precise statement of the problem is given and the steps of a basic algorithm for its solution are described. Once the basic algorithm is understood, various changes in the individual steps of the algorithm and in the algorithm itself can be made to improve the efficiency of computation. These are discussed in this section and also in Chapter 5. Theoretical results concerning disturbance input effects on decoupled systems are presented in the second main section. Finally methods for understanding, reducing, and in some cases eliminating static gain variations and cross coupling due to plant parameter variations appear in Section 3. 3.1 Control Effort and Error Performance Procedure One approach for the design of multivariable systems is the well known quadratic tracking problem [K1, All and its extension to model following techniques [Y1, T11. Although it is possible to use these to design multivariable systems with reasonable error and control amplitudes, how is a decoupling requirement included? Unlike the technique which follows, there does not appear to be any reasonable way to use such an approach to achieve decoupled multivariable designs except in an approximate way. Appendix A contains a discussion 22

23 of the difficulties one encounters. 3.1.1 Problem Approach Given the plant (2.1), we assume that it can be decoupled, that is, any necessary dynamic compensation has been included and the fixed constants in equations (2.16) for the F and G matrices of a control law have been determined. Let us consider now the synthesis problem of selecting values for the free scalars Xi, aik i = 1,..., m; k = 1,..., Pi. Although these scalars appear in the expression of the decoupled sub-system transfer functions hi.(s) in a simple way, their effect on the plant inputs is not clear. We begin by selecting a random input vector v(t). This vector is generated by passing zero-mean Gaussian white noise through a filter =- Av~ + B w(t) (3.1) v(t) = C v where ~ is a v-vector representing the state of the filter, Av, B and C are constant matrices of sizes v x v, v x m, V V v and m x v respectively, and w(t) is a zero-mean Gaussian white noise m-vector with co-variance matrix E[w(t)w'(T)] = r6(t - T) (3.2) where r is a diagonal matrix and 6(t -T ) is the Dirac delta function. A vector or matrix transpose is indicated by a prime. We assume that the filter is asymptotically stable (eigenvalues of Av have real parts < 0), and the filter

24 consists of m independent single input-single-output subsections which are represented in the combined form of (3.1) only for convenience of notation. Now combine the plant and filter states into the augmented state vector z of size n + v and the augmented system equations can be expressed z = Az + 9w(t) A (3.3) y(t) = Cz where A A + BF BGC | A 1 A] v B, C = IC 0. We define a covariance matrix for the augmented system vector P(t) = E[z(t)z'(t)] then it is well known EA21 that the covariance matrix of a linear constant coefficient dynamical system obeys the Lyapunov differential equation A A A A P(t) = AP(t) + P(t) A' + BrB' (3.4) A with initial condition P(to) = P. For a stable matrix A, this matrix differential equation has a unique symmetric steady state solution matrix P = lim P(t) given by the t+o solution of: A._._A A A AP + PA' = -BrB' (3.5)

25 Define another 2 m vector E(t) consisting of the error e(t) = v(t) - y(t) and control u(t) vectors. v(t) - y(t) -C C (t) = = y(t) Wz(t), W v v (3.6) u(t ) F GCv For convenience let the expected value notation with a bar above indicate the steady state expected value if it exists, and supress the time argument e. g. lim E [i' E] = tom E[E'(t)i(t)]. Now define a cost functional J consisting of a quadratically weighted sum of the steady state error and control variances J- EE['Qo ] Q [ (3.7) where the m x m constant matrices Q and R are positive definite and symmetric. Then it is an easy derivation to show that the cost functional can be expressed in terms of the trace (sum of diagonal elements) of a matrix product J = tr[W'QOWP] (3.8) Observe that the cost is a function of the free scalars Xi, Cik, i = 1,..., m, k = 1..., Pi, which appear in the expression for W and in the Lyapunov equation (3.5) which has P for its solution. The synthesis problem can now be stated precisely as given a plant (2.1) with control law matrices (2.16), given an input filter (3.1), a diagonal covariance matrix r and weighing matrices Q and R, select the values of the free

26 scalars Xi, oik' to minimize the cost functional J of equation (3.8). Because of the highly complex way in which these scalars affect the cost functional, no analytical solution could be found and an iterative computational procedure for the minimization problem was developed. Before going into this procedure, however, we shall discuss in more detail, selection of the random input v(t) and selection of the weighing matrices Q and R. 3.1.2 Input Selection We assume that the filter (3.1) consists of a combination of m completely independent single-input, single output filters. The order, pole location and zero location of the filters may be freely chosen although we do assume that the filters are stable. If the i-th filter element has ordera. and is described by the equations \;i - Avi + bviwi(t) (3.9) vi(t)= cvii whereti is aiu vector, AV is b is a -). column 1nii vetr v. 1 vector, and Cvi is a -. row vector, then one simple way to combine the individual filters is to let

27 A 0 b 1 - ~ -I. - -- I - - i' 0 A 0'b A- — hv- B 2 V' V e II 0 0 Av bv tdm ilf m Cv kl Icy2 v v2 V I V m it is also assumed that the spectral density of each element v(t), i =,..., m of the class of inputs v(t) which the closed-loop decoupled system will be asked to follow is known. The output spectral density Sv (w) of each individual filter on the other hand can be calculated from the well kn own formula i i (3.10) where Sw (w) is a constant, (white noise) and Hi(jw) is the Fourier transform of the filter inpulse response hi(t). One rationale for filter selection, then, is to choose a filter so that Svi(w) approximates the expected content of vi(t). In practice low (first or second) order filters seem to be adequate. Some experimental results on the subject of filter selection and its effect on the synthesis problem are presented in Chapter 5.

28 3.1.3 Co'st''Function al In general, no single procedure for selecting the elements of the weighting matrices Q and R can be stated. We assume that some knowledge of the acceptable control magnitudes and acceptable error magnitudes is available. In order to get started, we may select Q and R to be diagonal matrices with element values chosen to reflect this knowledge of acceptable nagnitudes. Furthermore, the relative weight of the elements of Q compared to R should reflect our tradeoff of error magnitude for control amplitude. This latter consideration is particularly important in situations where the plant (2.1) will be asked to follow inputs with significant high frequency energy components compared to the time constants of the plant itself. In that case, we must accept either large controls or large errors, or a tradeoff where each is significant. Usually the synthesis result using our first choice of Q and R will not be entirely satisfactory. After examining the control and error variances, we may decide some are too large and decide to choose another Q-R pair, this time weighing more heavily (increasing the corresponding elements of Q and R) those errors or controls which turned out to be unacceptable. Another synthesis result can be calculated with confidence that these objectional items will be reduced. Because no statement can be made about exactly how much each will be reduced for a given increase in weight, several tries may be necessary. In practice only a few such

29 tries are necessary before reaching an acceptable design. Some practical results which demonstrate this technique and the effect of changes in the weighting matrices are presented in Chapter 5.

30 3. 1.4 Form of Gradient and Hessian A Newton procedure is used to minimize the cost functional J. In this section, we shall examine the form of the gradient vector and Hessian matrix of J, both of which are the essential ingredients of the Newton procedure. The procedure itself will be described in detail in the next section. Assume we have already obtained the matrices D-1A*; Gi i 1,..., m; Jik' i = 1,..., m, k = 1,... Pi the scalars W ik' i = 1,..., m, k 1,...' Pi of a control law (2.16). Select a starting value for the free scalars xi and aik' i 1,..., m, k = 1,..., Pi. Call these starting values lio and aiko. For convenience we now convert the double subscript ik to a single subscript j. The conversion is given by i-1 j =k + and ca =a j ik The control law (2.16) can then be rewritten in terms of the new variables X., i = 1,..., m andj, j - 1,..., P p F = F + ajF. (3.11) G=G + X.G. 0 i=l1

31 where P- P., a- = aj -' -''i. and -=1 1 1 1 ]O 1 10 F -D-1A* + (a. )F 0~~ r j= t ]o ] J (3.12) m G - 2 iG. o ioGi Note that the new variables oj and x. all have a value of ] 1 zero when aj. and. x. io This property is very important in reducing the complexity of the expressions for the gradient and Hessian by removing cross product terms which would otherwise appear. Corresponding expressions are written for A and W Am AL A + a.A + Ai o j-1 j ] i i (3.13) W=W - +X GW + xiW o j1 ] i=1 1 11 where m A + BD 1A O ~ j-1 io~] i 1 io i + O W~0 =- a m ~ J 1A j=l 1 J i=l 10 1 D-1A* 0 BF. 0 0 BGiC A. -j 1 [ Cv j =l ] i =a L0 1 (3.14) ] F O LO GCJ

32 Again for simplicity of notation we re-name the i variables + i = 1,..., m (3.15) Now taking the gradient of J, that is, the derivative with respect to each of the variables a, j = 1,..., p + m from (3.8), we denote the jth element of the gradient vector evaluated at a. - 0, j 1,..., p + m as J (0) tr+ W QOWjP + iQoWoP1 (3.16) Note that the last two terms are symmetric and thus have the same value under the trace operator. The terms P may be aoj obtained by implicit differentiation of the Lyapunov equation (3.5). A a + a A' -A.P - PA.' (3.17) Evaluating everything at.j - 0, j = 1,..., p + m, we find that if A is stable, then P has already been specified as the solution of (3.5), and A = A is also stable. Then (3.17) also has a unique p. d. solution DP which may be asj used above in evaluating the gradient. We differentiate equation (3.8) twice to obtain elements of the Hessian matrix 2 r 2i (0) 2 tr j(WjQoWiP + Wo'QoWi. —-- + W'QoW. I + tr LwiO 2.] 1 (8 I 1 ] I (3.18)

33 The matrices P, P having been previously determined, we obtain a2P from (3.5) by differentiating it twice 2- 2Ao.o.aC + - - A -A. aP - PA - A. aP - PA.' (3.19) 1T 1 i J We state, as before, that equation (3.19) has a unique p. d. symmetric solution if Ao is a stable matrix. 3.1.5 Algorithm The basic steps of an algorithm for the minization of the cost functional J (o) by means of a Newton procedure will be described in this section. The Newton procedure itself consists of repetitive application of the algorithm. k + i k -1 ka lk -Hk VJk (3.20) where akdenotes the independent vector variable at the k' th iteration step, and Hk and VJk are the Hessian and gradient respectively of the cost functional evaluated at a=a. This procedure is known to converge very rapidly in the vicinity of the optimum a*. Exactly how close the starting point must be to a* is difficult to determine for this problem in light of the fact that existence or uniqueness results are rmt available. In any event, our purpose is to locate a minimum (if it exists) and not to study the fine details of minimization procedures. In line with current practice, the existence of additional local minima

34 can usually be detected by using a wide variety of starting points. Only one minima has ever been found for any of the examples tried. For clarity, a flow diagram of the computational process is shown in Figure 3.1. Step 1 The matrices A and W are calculated from the current O O values for,i' i = 1,..., p + m. This step must be repeated after each Newton iteration in order to guarantee that the variables.j and Xi have a value of zero. Otherwise the expression for the gradient and Hessian will not be valid. Step 2 The Lyapunov equation (3.5) is solved for P. Solution methods for this equation is a subject of current research. There are several techniques available at the present time. For reference see Chapter 5 and EBl, M1, M3, D1, K41. The most obvious is to convert the (n +v) x (n + v) symmetric matrix P into a (n + v) (n + v +1) vector p. Because P is symmetric, only the upper triangular part needs to be determined. Let the i'th element pi of p be given by the j, k'th element [Pjk] of P, i. e, k (k - 1) i = " 2 + J The symmetric matrix BrB' is simularly converted into a vector q. Then equation (3.5) can be written as A p = q (3.21) BIG where ABIG is a square matrix of size (n +v ) (n +v +1). A 2 If A is a stable matrix, then the large matrix ABIG must be

35 |Read Data anul Print It FA A = h.A. +,A iterate W W +. w +C -W O 1+ jwj B Solve AP + PA' = -BEB' J =tr W'Q WP S -A.P -'PA. - LSolve APi + PA' S loop J. = tr I'Q WP. + W'QQW.P + W'QWP} 0 1 L 1 0 i = l,i3.,p j Inver Hiess~ian w ll S -Calculate P AiPj PjA! Figure 31 Steps o flve Algor Pithm S J. - 2tr tQ WP + W'QoWiPj +W'Qowi] + trWIQ oWPi Invert Hessian Calculate New Step| Figure 3.1 Steps of Algcrithm

36 non-singular because equation (3.5) and thus the equivalent system (3.18) is known to have a solution. = AB1q (3.22) BIG Techniques have appeared for forming the matrix ABIG from A with a small number of steps. See EB2] and [Cil. Inversion is a straightforward matrix inverse problem. Other solution methods for equation (3.5) seek to avoid inverting such a large matrix by using an iterative procedure instead. Another item to consider for reduction of the computational burden is to partition equation (3.5) A + BF BGC P P (A + BF) 0 11 12 11 1 (A + BF) 0 A P P (BGC )' A' [7:v 21 2 21 221 2 v v 0O 0 o BrB1 l Writing the partitioned equations: (A+BF) P11 + BGCvP21 + P11 (A + BF) + P21 (BGCv )= 0 1A + BF) P21 + BGCvP2 1 11 v (A + BF) P12 + BGC P + P A' - 0 12v 22 21 (3.23) A P + P (A + BF)' + P22(BGCV)' = 0 212 12 22-BB v 22 22 v

37 The third equation is merely the transpose of the second and presents no new information. In addition, the fourth equation is independent of the a. variables and needs to be solved only once for each filter selection, while the first and second equations need to be solved at each iteration of the aj variables. Step 3 Calculation of the cost J. Simple matrix algebra and a trace operation are all that is needed here. aP Step 4 Determination of P. a 1 1 The right hand side of equation (3.17) is computed using the matrix P determined in Step 2. Then the Lyapunov equation is solved for P.. Note however that if the straightforward matrix inversion technique is used to solve for P, then the same matrix inverse is used to determine P. so there is no lengthy calculation here. Step 5 Calculation of the gradient terms Ji i(o). Again, simple matrix algebra and a trace operation are used. The computation is somewhat longer than in Step 3. Take note that Steps 4 and 5 must be repeated for each element of the gradient vector. Also the Pi matrices determined in Step 4 must be stored for use in Step 6 below. Otherwise they will have to be recomputed. aZP Step 6 Determination of P. =........13 1

38 The right hand side of equation (3.19) is computed using the Pi, i = 1,..., p + m, matrices determined in Step 4. Then the Lyapunov equation is solved for P.. as before, 1J making use of the simplification discussed in Step 4. 2 Step 7 Calculation of the Hessian terms Jik ik Matrix algebra and a trace operation, this time on an even lengthier expression is used. Steps 6 and 7 must be repeated for each element of the Hessian matrix. Because of its symmetry, however, only the upper triangular portion needs to be computed in this way. Step 8 Invert Hessian An inversion routine which takes advantage of the symmetry of the Hessian matrix is preferred here. Step 9 Calculation of the step size k + 1 The new starting vector for the k + 1 iteration o is computed by equation (3.20) from the kth starting vector: k + 1 k H-1 (3.24) (a - H VJ) (3.4) k where H and VJ are evaluated at ok. The matrices A and W 0 0 are recomputed and another iteration is taken. There are many possibilities for a stopping criteria and the one selected for this problem is to specify a maximum number of steps and after each step, a norm of the gradient vector specified by norm- =lJ 1

39 is computed and compared to a specified scalar. If the norm of the gradient is less than this specified constant, then it is assumed that the minimum has been reached, and no further steps are taken. If one uses the method of inverting ABIG to find P, then the Newton procedure is preferred over a pure gradient or conjugate direction method. The reason is that it is difficult to determine the cost J by the method of inverting ABIG, but once A1iG is determined, computation of the gradient vector and Hessian matrix elements does not require much additional effort. The additional computer time required to compute the Hessian and invert it is probably well justified by the more rapid convergence of the Newton procedure, particularly in the vicinity of the optimum. However, if an iterative method [K4] is used to find P, then it is more reasonable to use a gradient procedure because the iterative method would have to be repeated once for each element of the Hessian matrix. Many authors suggest that for the first iteration with a Newton procedure, only a partial step be taken as a precaution against a serious overshoot which is theoretically possible. For the practical examples considered in Chapter 6 this precaution was found to be unnecessary, because overshoots seldom occurred. For those few cases in which it did, the overshoot entered the region of an unstable A matrix. A simple correction procedure which worked well

40 was to cut the initial step in half and use this point as a new starting point. In all cases, this procedure was sufficient and convergence was achieved. Additional questions in regard to efficiency of computation and computer run time requirements are discussed in detail in Chapter 5. 3.1.6 In'dependent Inputs Because the closed-loop system is decoupled, it is possible to reduce the problem as stated in Sections 3.1.1 - 3.1.5 into m independent problems of lower order as the following development demonstrates. The cost functional J from (3.8) can be expressed as: J = tr [W'Q WP] = tr [QoWPW'] (3.25) Also denote the upper left hand part of the error and control variance matrix WPW' as (Xi) where Xi is an m x m error covariance matrix. For a decoupled system, the individual subsystems are independent single input-single output systems; consequently, the errors e.(t) - v.(t) - y.(t), i - 1,... m 1 1 1 are independent and thus Xi is a diagonal matrix. So for decoupled systems, off diagonal elements in Q have no effect on the cost J. We express J equivalently as J Eu'Ru + q.e. 2] (3.26) i=l 1 1

where q. is the ith diagonal element of Q. The control u can be expressed as the sum of control responses due to the individual inputs because of linearity. m i u =- u (3.27) i=l where u is the control response due to vi(t) with v.(t) = 0, j ~ i. Thus m u m u'Ru =, (u )' R (uk) i=l k=l m m i k -E E (U )?Ru i=l k=l i k Now u and u are independent and thus uncorrelated. Hence m i E[u'Ru] - (ul)'Ru i=l and 2 m J i 2 A m i J E (ui)TRu + q.e.:. J (3.28) i=l i= i-l where i - i i 2 J = E[(ui)Rui + qie 3] (3.29) Using the canonically decoupled form EG2] described in Appendix F it follows that if a suitable change of coordinates for state space is made (x = T2x) and the resulting A 1 state vector (x) is partitioned into subvectors x,.., x + 1of dimension P1 Pm + 1' the equations describing i i u and y can be written as

42 i ii x (t) =A.xi(t) + b.(X.v.(t) + a x 1 1' 1 1 x (t) A x m+lt) + ACxi(t) + b{(X.v (t) + 0.xi(t)) (3.30) 1 i i m+l u (t) W.x (t) + W+1x (t) + V. v.(t) + 1 m+l i V. 0 x (t) y (t) = ci.x(t) where the matrices Ai, bi, Am+l, A b, ci are given as partitions of the CD system matrices a, P, C in EG2] which are also shown in (3.34). Since there is a one-to-one correspondence between the elements of the row vector 0. and the elements of a = (aip,.., ai 1) it follows that J is a function of X. and a only. As a result, the overall minimization problem has been separated into m smaller (in the number of search variables) problems. This separation into m sub-systems has the consequence that even if the original problem is solved directly without separation, the Hessian matrix is now known to be block diagonal in blocks of sizes (Pi + 1) x (Pi +1). The blocks are Pi x Pi if the unity static gain constraint is imposed. This simplification means a considerable reduction in the number of elements of the Hessian matrix which need to be calculated and a consequent reduction in the computer time required. The saving is discussed in later chapters.

43 In addition to the reduction in the number of search variables, the m individual problems will each have smaller state dimensiors than the overall problem. A saving in overall computer time may be possible because of the reduced dimensions even though m individual problems must be solved. Questions of this nature are discussed in Chapter 5. If disturbance inputs are present, then it is unlikely that their effects will be confined to just one subsystem. In that case, the preceding results do not apply.

44 3.2 D'isturbance Inputs 3.2.1 Complete Cancellation Re'sults Consider the case when the plant (2.1) has disturbance inputsvhich can be described by the modified set of equations - = Ax + Bu(t) + Nr(t) (3.31) y(t) = Cx where r(t) is a q-vector (elements r.(t), i = 1,..., q) of disturbance input functions. The matrix N is a n x q constant matrix. Definition'3'.1 - The output y(t) is said to be invariant with respect to the (scalar) disturbance input r.(t) if y(t)- 0 for all r.(t) and x(o) = 0 It is assumed that Morgan's control law (2.2) is used for closed-loop control. A necessary and sufficient condition for invariance was first developed by Wang EW1]. This condition is stated below in Theorem 3.1 as it is reported by Cruz and Perkins EC3]. Take note that the theorem applies whether F and G decouple or not. Theorem 3.1 A necessary and sufficient condition for the system (3.30) with control law (2.2) to be invariant with respect to r.(t) is that 1

45 C(A + BF)k v 0= for k = 0,..., n-1 where v is the ith column of the matrix N. For a proof of Theorem 3.1 see Wang [W1]. It is easy to show that a consequence of this condition is that except for the trivial case v a 0, the closedloop system - = (A + BF)x + BGv(t) + vr.(t) (3.32) y(t) = Cx must not be completely observable. Returning now to the additional requirement that F and G decouple, we shall see what effect Theorem 3.1 has when applied to the decoupling problem. The problem we shall consider can be stated as: Given A, B, C what v are permissible so that y is invariant with respect to ri(t) and the system (3.31) with control law (2.2) is decoupled? In general, for a given plant S {A, B, C} the space of permissible v given by =L {v: veR, C(A + BF)kv=O, k=O,..., n-l) (3.33) will be dependent on the free elements aik, i = 1,..., m, k = 1,..., Pi in the expression (2.16) for the decoupling feedback F matrix. If 7is known in terms of the aik' then we can examine the columns of N and determine which of the dis

46 turbance inputs can be made invariant (complete cancellation) and which cannot. For those which can it may not be possible to make them all invariant because of conflicting requirements made on the a ik elements. One has to choose which disturbances are more important. For those which cannot and those which would require an unacceptable feedback F matrix, the results of Section 3.2.2 can be applied. Some results can be seen more easily if we consider the problem transformed into CD form as in Appendix F. A A T T (A- BD A*)T1 T2/1 2 A B = T T BD2 1 (3.34) A C - CT -1 T - 1 2 A T2T1 A A A Note: These particular A, B, C matrices should be distinguished from other matrices with the same notation in Section 3.1. The A A A form of A, B, C is given by equation (F2 ) of Appendix F. Furthermore, in this CD form F can be expressed as F f 0 f2 0 (3.35) f O m

47 A Partition the vector v similarly into m + 1 sub-vectors each of dimension P i = 1,..., m + 2 A v -I v V2 (3.36) Vm+l The following results can now be stated. Result I The allowable v is completely free for all feedback matrices which decouple. Result II The conditions for the allowable v:, i = 1,..., m are independent of each other. th and depend only on the respective i subsystem A., b., c., f.. 1.1 1 1 Proof: From the partitioned form of the matrices given in Appendix F and (3.35) we calculate for any power k the expression A A AA C (A + BF) = C (Al+bfl)k 0 0 0 0 G2(A2+b2f2)k 0 0 k 0 0 c (A +b f) 0 m m m m The location of the many partitions which are zero for any k, k = 0,..., n - 1 establishes Results I and II.

48 Now consider the individual subsystem A, bi, Ci fi, vi. From Result II we can examine each subsystem i = 1,..., m separately, and by the Cayley Hamilton theorem, powers of (Ai + bifi) higher than (Pi - 1) can be expressed in terms of lower powers so we can write the condition of Theorem 3.1 for a subsystem as k ci(Ai + bifi) v. 0 for k 0=,..., pi- (3.37) rather than for k 0,..., n - 1. This simplifies the procedures for finding the permissible vi. Result III The space Hi {vi:vi.ER, ci(Ai+bifi) v. = O, k = 0,..., p. - 1} will have dimension d(ti.) o <d(A.)_<r. and in general will depend on the actual values taken by the free scalars a ik k = 1'' Pi Proof: Form a matrix M. consisting of the row vectors c.(Ai + bifi), k 0,..., p. - 1 C. c.(A. + b.f.) M. c.(A. + b.f.+ ) 1 1 1 1 1.~ p. - 1 Ci(Ai + bifi) 1 1 1 1 1

49 The space.i is given by the null space of this matrix. From the form of Ai bi, ci', fi we calculate the upper left hand (di + 1) x (di + 1) elements (d pi - ri -1) which turn out to represent an identity matrix. Clearly rank M. > d. + 1, thus the null space has dimension 1 1 < ri.. 1 Obviously as a consequence of Result III if r. 0= then v. 0 is the only element ofti, that is if there are no numerator dynamics in hi.(s) then that subsystem will be completely observable for all choices of aik, k 1,...' Pi In fact, the only possibility for vie. when vi. 0 occurs when pole-zero cancellations are made in hi.(s) by a suitable choice of aik, k = 1,.'., Pi. Recall that the numerator dynamics are fixed, i. e. cannot be altered by state variable feedback alone. Conditions on the o ik variables in order to cancel several numerator roots are rather complicated, however, some rather general results can be obtained. It will be easier to prove the results if we make an additional transformation to companion form indicated by the v superscript. Transform v v v v v A., bi' ci' fi' vi A., bi' ci) fi) v.. This transformation which corresponds to a change in coordinates in state space can always be made, see [G1], and for a decoupled subsystem, the matrices will have the following form [G2].

50 v I v O A. L bi 0O LOi0 0 1 V v. = K. v 1 1 i (3.38) c. = [-a. 1 0...0] V 1 [Pi Pi-1' a rri tr.' ar. - where ai. is an ri row vector ai. = [a a' 1 1 1 ri r.-l' 1 1 and sri - rsri-l ri h.(s) = 1 1 sPi s Pia _j -. -.P and Ki is the nonsingular transformation matrix. The conditions for complete cancellation of disturbances can be stated for the the companion form sub-systems v v v k v c. (A. + b.f.) v. 0= for k = 0,..., p. (3.39) 1 1 1 1 1 1 by applying the transformation to (3.37). Silverman ES1] has recently obtained powerful results for a more general disturbance cancellation problem which can be applied here. We state the following result

51 which is based upon the idea in his Theorem 5.2 but as stated here applies only to single input-single output companion form systems (3.38) and is easy to prove by direct expansion without reference to his work. It involves a matrix L which he introduces. Theorem 3.2 There exists a matrix f. such that condition (3.39) is satisfied if and only if Li) 0. The (p. - r.) x Pi matrix L for this problem is given by: Ci L = i 1 (3.40) v.v Pi -ri eiAi Proof: Necessity is established by noting from V V V V the form of ci, Ai, bi, fi just given that V v v k V k i (Ai + b. f) = ciA. for k = O,..., Pi - ri - 1. These are the rows of L. For sufficiency, choose 1 0= 0f rk < i r< Pi then Cv(A. + b.fl ) 0 for k > P r. - Czi + bifi)= for k > Pi i This result establishes clearly th at there'does exist a matrix f. such that d(?. ) = r. and that this space 1 1 1

52.i can be calculated as the null space of L. The column V vector vi must lie in this null space' in order for it to satisfy (3.39). If there are several disturbance inputs, V then the vectors vi corresponding to each disturbance input must lie in this null space, otherwise, they cannot be comv pletely cancelled by any choice of f.. Direct calculation shows that elements in the null space of LB are given by vi., 1 vi, 2 i). -=. free vi, r. vi, r. + 1 = ij v j=- ( 3.41) r vi, r. + 2 1-= jj =1 "i' Pi -a. mi'j Pi-j j-1 Another result originally due to Silverman [S1] can be applied to calculate conditions on the elements of V a matrix f. which satisfies (3.39). 1

53 V V Theorem 3.3 Let v be such that L v. = 0. Then (3.'39') is satisfied if and only if V V V fi(Ai.) v. - 0 for j = 0,..., p. (3.42) Proof of this theorem is lengthy and is omitted. The theorem is established by direct application of Silverman's Theorem (5.3) for the companion form shown. It is clear that simultaneous satisfaction of LSV 0 and (3.42) is equivalent to satisfaction of (3.39). For a general companion form system, the two equivalent expressions are quite difficult to solve in terms of the parameters a.. On the other hand, if specific numerical values are available for a particular companion form system, then Silverman's results are much easier to use in calculating conditions on the elements of f. than (3.39). General solutions have been found for two specific cases. They were developed before Silverman's results became V available directly from (3.39), and application of Lv = 0 and (3.42) do not simplify either case. The results are presented because they do provide general solutions which would otherwise be unavailable and because they apply to the great majority of examples, that is, those for which p.~3, i = 1, 1., m as well as some others. The first case for which a solution can be found easily is if ri = 1, i. e. Pi = di + 2. 1 ~~~~~1

54 Result IV If r. = 1 then 1C-* 1 0 0 0 _i 1 v V V 1 Ci(Ai+bifi) 0 -a1 1 0 0 0 O -al 1 0 ~ ~~~- ~~(3.43) v v v v Pi-1 Ci(A.+b.f.) 1Pi 1 1 -11 1 1 Ifti. 0, then a1 c2 free (3.44) 1 I ~Pi-2 P 1i a cl -i a i2 -...ap. 1C'Pi 1 1 1 2 1''' - and for any Vi Vi ii Pi-l

55 Proof: Obviously, the vector vi must lie in the null space of the Pi x Pi matrix. The partitioned portion is triangular with all l's along the diagonal and thus has rank di + 1 which agrees with Result III. For a ~ 0, we require that the entire matrix be singular. Evaluating its determinant along the last row and requiring it to be zero gives 0 C- + CY (- a) -al(c )2 + Pi pi-1 1 Pi1 +(CT1-l) (-al)Pi-1 and consequently (3.44). From (3.43) with condition (3.44) it can be seen by inspection that for any V V vi-mi condition (3.45) must apply to vi. Note that the expression for the permissible vi is independent of the free values of the ai variables. Furthermore, note that Pi - 1 of the ai variables can be arbitrarily specified since it is only necessary to cancel one numerator root. The other case for which a solution can easily be found is if ri = Pi-l and all of the numerator zeros are cancelled. The resultant h.(s) will then be first order. In that case Result V If ri - 1 then v (3.46) i 1(3i-1 1 V V V 1 1 ci(Ai+bifi) |a|Pi Pi (ai 2-a2 ) (a -cL ) v v p._]. ci(Ai+bifi) Pi 11

56 and a necessary and sufficient condition for (3.46) to have rank = 1 is that (a2 - a2) = -a (al - l ) (a3 - a3) = -a2(a al1 ) al' free (3.47) a-r ( - ( pi-l Pr 1 1 CYPi ria ( 1 al Furthermore for this relationship we have for v v any l aii v i2. = i2 free (3.48) iri ri j= j i ri +1 j Proof: The details of this result which are somewhat involved can be found in Appendix D. The equation (3.47) gives the greatest possible dimension to7i. when ri = Pi-1, that is, complete cancellation of the numerator. Of course, greatest dimension in71i results in the strictest condition on selection of the a ik variables so that only one element of oi can be chosen freely. Notice

57 V that this expression for the permissible vi is also independent of the particular value of a1. The transfer function for condition (3.47) on the aik' s will be hi(s) = 1. - It appears that partial cancellation of the numerator when ri = Pi-1 or for either complete or partial cancellation when 1 < ri < Pi - 1 is a particularly difficult problem because there are many ways in which the cancellation can occur. Another technique which may be considered is the possibility of feedforward control of measurable disturbance inputs by means of a control law q u(t)= Fx + Gv(t) + K Klri(t) (3.49) i=1 where K. is a n x 1 constant matrix and the summation is made over those indices i which are measurable. One simple result by McLane and Davidson [M2] is available here, but in general, the overall problem of using both feedback and feedforward control to cancel disturbance inputs in decoupled systems is an open area for research with no general results available at this time.

58 3.2.2 Inclusion of Disturbance Inputs in Minimization Probcedure Assume we are given a plant with disturbance inputs (3.31). If the effects of the disturbance inputs cannot be eliminated from the output y(t) by the results of Section 3.2.1, or if we choose not to eliminate them completely, then an alternative procedure is to include their effect in the minimization procedure of Section 3.1. We assume the spectral content and magnitude of the disturbance vector r(t) can be expressed as the output of a filter driven by a white noise source. This filter is described by the dynamical system - An% + Bnn(t) (3.50) r(t) = Cn where ~ is an n state vector; A n, B, C are constant matrices of size n x nn, nn x q, and q x n respectively; n(t) is a white noise source q-vector. E [n(t)n'(T)) = r 6(t - T) where r is a positive definite symmetric n x n correlation n n n matrix. It is also assumed that the noise input n(t) to filter (3.50) and noise input w(t) to filter (3.1) are independent and therefore uncorrelated. E[n(t)w' (T)] ] 0

59 All of the comments made about selection of the filter characteristics in Section 3.1.2 apply here to the question of selection of the elements of An, Bn Cn We assume furthermore that we wish to decouple the plant (2.1) with a control law (2.16). Now form the augmented state vector z, and augmented noise input vector ~. z [ ~(t)= [w(t)] +(tR L1 and combining equations (3.31), (3.1), and (3.50) z = A z + B Z(t) (3.51) y(t) = C z where A + BF BGC NCj O A - = 0 ~ A, B = (3.52) 0 0 An' Bn C =[ C 0 0 As before, define i(t) as the error and control vector v(t) - y(t) i(t) = - | = W (t) (3.53) u(t) where -C C 0 W v (3.54) n F GC O ~ ar n

60 We observe that the inclusion of a disturbance input to the minimization problem of Section 3.1 has not made any significant difference in the problem formulation. Equations (3.51) and (3.53) are analogious to (3.3) and (3.6) respectively and the same procedure can be used to solve the problem posed in this section. As before, it may take several tries with different weighting matrices Q and R in the cost functional (3.7) to achieve an overall acceptable design, however, also as before we have the advantage this technique offers in knowing how to modify a previous choice of Q and R in order to reduce an unacceptable control or error magnitude. As we note in Section 3.1.6, the disturbance inputs in general couple to all aik variables and consequently the independent inputs simplification of that section does not apply when disturbance inputs are treated as above.

61 3.3 Sensitivity 3. 3.1 Comple te Cancellation of Parameter Variation Effects Consider the closed-loop system (2.3). If the A matrix changes by 6A from its nominal value A 0 A = A + 6A (3.55) 0 then we would like to know under what conditions the response y(t) for the perturbed system is identical to the response yo(t) of the unperturbed system. Some results for this problem are reported by Cruz and Perkins [C3] however, the development which follows is more general and includes their results. Definition 3.2. The closed-loop system (2.3) is said to be invariant with respect to 6A if for all inputs v(t), t > 0 and all initial conditions x(O), y(t) — y (t), for all t > 0. o The output of the mnperturbed system given by (Ao+BF)t t (A +BF)(t-T) yO(t) = Ce x(o) + JCe BGv(T)dT (3.56) and the perturbed system (A A+BF)tt (A + A+BF)(t-T) (3.57) y(t)=Ce x(O)+ Ce BGv(T)dT 0 Because x(o) and v(t) are arbitrary, we have the necessary and sufficient conditions:

(A +6A+BF)t (A +BF)t C[e o - e o = 0, for all t20 (3.58) and (Ao+6A+BF) (t-T) (Ao+BF) (t-T) C[e - e ]BG = 0, (3.59) for all t<O, O<T< t Clearly (3.58) implies (3.59). Using the series expansion for exponential yields tk k C -. [(A +6A+BF)k - (A +BF)] 0 for all tk0 k=O 0 0 (3.60) there fore k k C[(A +6A+BF) - (A +BF)] 0 for all integers k>0 By mathematical induction and the Cayley-Hamilton theorem, we establish the necessary and sufficient condition that (2.3) is invariant with respect to 6A, C(Ao+BF) kA = 0 for k = 0,..., n-1 (3.61) By following the above steps it is also possible to show that if B changes by 6B from its nominal Bo, but A and C remain fixed, then a necessary and sufficient condition that (2.3) is invariant with respect to 6B is that: C(A+B F)k6BF 0 for k 0,..., n-1 (3.62) Unfortunately, there are no simple conditions if both A and B change simultaneously. Observe that conditions (3.61) and (3.62) on parameter variations are essentially identical to the condition of

63 Theorem 3.1 for invariance to disturbance inputs. Consequently, all the subsequent results of Section 3.2.1 apply directly to this problem. We merely require that the columns of 6A or 6BF simultaneously satisfy conditions on v so that (2.3) is invariant with respect to 6A. Conceivably, the unperturbed system may be invariant to a disturbance input whereas the perturbed system is not. The following result applies to this possibility. Definition 3.3. The closed-loop system (2.3) is said to be simultaneously invariant with respect to a disturbance input ri (t) and a parameter variation 6A if for all inputs v(t), t> 0; ri(t), t> 0 and initial conditions x(o), y(t) = yo(t) for all t > 0 Theorem 3.4. The system (2.3) is simultaneously invariant with respect to r.(t) and 6A if and only if C(A + BF)kv = 0 (3.63) o and k C(A + BF) 6A = 0, for k = 0,..., n-l (3.64) 0 Proof: Necessity is established by following the previous development where y(t) = Ce (A F)t jtCe(A+BF)(t- T)BGv(T)d 0o t (A +BF) (t-T)

64 and noting that C(Ao+BF) v = 0, k = 0,..., n-1 is necessary, from Section 3.2.1 for the definition. Sufficiency is proved by observing that the necessity proofs are also sufficient for the conditions of the definition and for arbitrary v(t), ri(t), x(0).

65 3.3.2 Calculation of Static Sensitivity' and its Gradient In this section, we consider the effect of parameter variations on the static gain of a closed-loop decoupled system. The approach involves definition of a static sensitivity matrix, giving its interpretation, and the computing the effect of changes in the feedback parameters upon it. Some preliminary results will have to be collected first. Extended Leverrier Algorithm Given a sum* G = A + aJ of two square n x n real matrices A and J where a is a real scalar, then the well -1 known expression for (sI-A) can be extended to this sum G yielding the following equations -l - — 1 (sI-G) (sI-A-aJ) q l(s, o)[R(s,o)] n-l n-2 n-3 R(s,a) = Isnl+Rl()s +R2()s +...+Rn _l(a) (3.65) q(s,a) = sn-ql()sn-l-q2 ()s -2_...-qn(a) - 2 The coefficients q.(a) i = 1,..., n and matrices Rj(a) j = 1,..., n-l can be calculated simultaneously by the use of the Leverrier algorithm, which is here extended to the case of a sum of two matrices. *also called a pencil of matrices, see Gantmacher [G1l.

66 Al(a) = A + aJ ql(a) = trA1(a) A2(a) = (A + aJ)R1 (a) q2(a) = 1/2 trA2(a) An (a) = (A + aJ)Rn_2(a) qn-l(a) 1 trA (a) n-i n-2 n-ii n-i n-i An(a) (A + aJ)Rn 1(a) qn(a) = 1/n trA (a) n n R1(a) = A1(a) - ql(a)I (3.66) R2(a) = A2(o) - q2(a)I R (a) = A (a) - q (a)I n-i n-1 n-1 R (a) = A (a) qn(o)I = 0 From the extended Leverrier formulas, one concludes that the elements of Ai(a), and Ri(a) as well as qi(a) i = 1,..., n are polynomials of degree n (or less) in the scalar a. Because the algorithm is recursive, the calculation of the formulas is easily performed on a digital computer with a small amount of storage required. Computation and storage of the constant matrices which are the coefficients of the powers of a in the expansions for Ri(a) and Ai(a) as well as the scalar coefficients in the expansion for qi(a) is straightforward. For use in later results we shall see that only the coefficients of the first and

zero powers of a will be needed. It should be clear from the algorithm expressions that only the coefficients of these two powers need to be stored at each iteration, rather th than the entire n order polynomial coefficients. Lemma 3.1 The matrix G = A + aJ is singular iff qn(G) = 0. Definition 33.'4 (A, J)cR1 is defined as 1,- {a:aoc R, det (A + aJ) = 0} Lemma 3.2 For any two square real matrices A and J of dimension n, the set.2 has either 0, 1,..., n elements or else J-R1. Proof: By previous lemma det (A + aJ) = Oqn (a) = 0. From the Leverrier development qn(~) is a polynomial in a of order n or less and thus will have at most n real roots. However if q (a) ) 0, then A + aJ is singular for all a. Lemma 3.3 If G = A + aJ is non singular, i. e. a/J, then (A + qJ) q (a )R (o) (3.67) Consider the system of closed-loop differential equations = (A + aJ)x + BGv(t) (3.6 y(t) = Cx, x(o) = xo

68 where A, J, B, C, and G are real constant matrices of size n x n, n x n, n x m, p x n, and m x m respectively, a is a scalar parameter, x an n-vector, v an m-vector, and y a p-vector denoting the state, input, and output respectively. Assume A + aJ is a stable matrix. Also assume v(t) is a piecewise continuous function for t > t and v(t) approaches A a constant vector v as t+ +-. Then it is known (see EK3 p.338] and [C21p. 74] that the solution (A + aJ)(t-to) + (A+aJ)(t-T) x(t) e x + e BGv(T)d xto (3.69) A will approach a constant vector x called the steady state A A state vector. The corresponding y Cx is called the steady state output vector. Lemma 3.4 Let A + aJ be stable, then the steady state output vector y is given by A -1 A y = C(A + aJ) BGv (3.70) Under the assumption of this Lemma we have Definition 3.5 The p x m matrix T(A, B, C, a) = C(A + aJ) 1BG (3.71) is called the steady state' gain matrix. Returning to the system of differential equations (3.68) consider the case when the matrices A, B, C change by 6A,SB, 6C respectively from their nominal values of Ao, Bo, Co. That is

69 A = Ao + 6A B = B + 6B (3.72) C = C + 6C Assume that both A + aJ and A + aJ are stable matrices, 0 then the steady state gain matrix T(A, B, C, a) = (C + 6C) (Ao + 6A + aJ) (B +6B)G (3.73) will also be changed from its nominal value To(Ao, B, C a) C (Ao +aJ) B G (3.74) 0 0O 0 0 Definition 3.6 Again, let m = p and also To(Ao, Bo, Co, a) be diagonal and non-singular. The sensitivity matrix M is defined as M = To (T - To) (3. 75) The diagonal elements of M may be interpreted as the percentage change in the diagonal elements of the steady state gain matrix with respect to the corresponding nominal diagonal elements of To. The i, jth(i j) element of M 0.th measures the steady state interaction between the jth input and ith output. -1 Lemma 3.5 Let m = p, if H (s,a) = C(sI-A-aJ) BG is a diagonal and non-singular matrix, then for all a/J(A, J), T(A,B,C, a) is also diagonal and non-singular.

70 Consider the open loop system = Ax+Bu(t) (3.76) y(t) = Cx A is n, n, B is n x m, C is p X n, Let p = m. Assume (3.76) is controllable and can be decoupled so that the control law (2.16) applies with the Kik=0O Temporarily fix the scalars aik ~iko i, k o i, k (3.77) and let a-oko = a (3.78) zioko - so F = F1 + Jioko' F1 -D-1A +ZikJik ik # ioko The the closed-loop transfer function -1 -1 H(s, a) = C(sI-A-BF) BG = C(sI-Ao-cJ) BG (3.79) where A = A + BF, J = BJi (3.80) o ioko Theorem 3.5 Given any decoupled closed-loop system with transfer function (3.79), let the elements of A change by 6A from their nominal Ao A - A + 8A, (3.81) 0

71 so that the corresponding elements of A also change by 6A from their nominal A A = A0 + 6A (3.82) Then dM is given by du dM'dL(A,) R1 (A, a) dc = a ) C ________BG (3.83) qn(A, a) qn(Ac)dRnl dq (A,ca) + L(Aox)C -(A,a)- R (A c) 0....d...... nd B _____ _____ ____ _-___ _"_'" )..... BG qn (A,cx) where L(A,o) = To is a diagonal and non-singular matrix with at most one element which is a function -l of a and that one is of the form r c, w a non zero scalar. Proof: If H (s, ca) is diagonal and non-singular then F has the form of (2.16) and H(s,a) = diag (hi(s, a) is obtained from (2.17) where ~ij and Ai are fixed i - 1,..., m, j = 1,..., ri and independent of aik* From Lemma 3.6, To will be diagonal and non-singular and can be calculated from T = H(o, ac). By inspection of (2.17) it is clear that ir. if io, ko = i, pi hi(s) C=.. (3. 84) s = 0 i a 1pi I1 if io, ko X i, Pi i iPi

72 which demonstrates that at most one element of T (and also T1 because T is diagonal) can be a function of a. All 0 o other elements of To are fixed (constant). The scalar 7, is non-zero because otherwise H(o,a) would be singular. Now by definition -1 M T (T-T) o o (3.85) -1 L(Ao, a)C(A + aJ) BG-I By Lemma 3.3 M = L (Ao,a)C' R (A, a)BG-I (3.86)' (A,a) n-l Taking the derivative gives the desired result (3.83). Application of Theorem 3.5 in the overall design of decoupled systems is discussed in Chapter 4. The evaluation of (3.83) is easier than might be expected. R 1('' ) and qn(', ) as previously mentioned are polynomials in a of order n or less, so the indicated differentiation can easily be calculated. Furthermore if we evaluate the derivative at a = 0, then only the coefficients of the first and zero powers of a in R 1('') and qn(',') are required, resulting in a reduction in computational effort as pPeviously mentioned.

73 3.3.3 Integral1 Cont'r'ol'for Static Sensit'ivity'Improvement The idea of integral control is invariably considered when one has stringent requirements on allowable variations in static gain. A general statement can be made regarding the case of integral control in multivariable systems, whether or not they are decoupled. The'author feels that this result must be known by other control theorists, however, it does not seem to have appeared previously in the form presented here. Frequency domain techniques are used to establish the result. Let the (m x m) matrix transfer function Q(s) describe a multivariable system with frequency domain input V(s) and output Y(s). Then Y(s) = Q(s)V(s) (3.87) In practice Q(s) might be a decoupled feedback system obtained by using the control law (2.16), i. e., Q(s) = H(s, F, G). We require Q(O) to be defined and nonsingular. Add integral control by letting V(s) = A1 + A2)E(s) (3. 88) E(s) = W(s) - Y(s) where A and A are m x m real constant matrices representing 1 2 the integration and proportional constants respectively of the controller. The controller operates on the error E(s) between the closed-loop input W(s) and the output Y(s). It is assumed

74 that Ai is nonsingular. An analysis of the combined system leads to: Y(s) Q(s)l 1 + A (V1() - Y(s) (3.89) [I+Q@) (1 A1 + A2)] Y(s) Q(s) (1 A1+A )V1(s) The assumption is made that the matrices A1 and A2 are selected so that the combined system is stable. Then if the input has a constant steady state A value v, we use the Final Value Theorem EK3] to calculate the A steady state output y and find: A A y = V (3.90) indicating that the static gain is independent of parameter variations in Q, Al, and A2. If Q(s) were diagonal it would be natural to specify A1 and A as diagonal. This would sim2 plify the choice of A1 and A2 since then the design problem separates into m problems each with only two parameters to consider.

CHAPTER 4 DESIGN PROCEDURE In this chapter we show how the theoretical results of Chapter 3 can be consolidated into an over-all design procedure for fairly high order systems. The nature of the approach involves trial and error design methods and by no means is a truly automated procedure. However it is systematic in the sense that after a particular trial is arrived at, its overall performance is evaluated as explicitly as possible within current theoretical limitations to determine whether the design requirements are met. If they are not, then systematic guidelines are given which, if followed, lead to an improved design. 4.1 Des'i gn' Data Assume that for a particular application, we have the following information: 1. A linear model of the plant 2. The decoupling requirement 3. The expected power spectrum Sv. () of the inputs v.(t), i - 1,..., m under actual operating conditions. It is further assumed that the actual inputs are uncorrelated, E[v.(t)vj(T)]- 0, i 0 j. I J 4. Control magnitude constraints.io' i = 1,..., m, where u.(t) _<.., i = 1,...,m, for 1 10 all time t. 75

76 5. Desired error performance - For each i = 1,..., m the variance'of the steady state error Ev - yi)2] < cio where the set of o io, i- 1,.., m are positive scalars. 6. The power spectrum S (w) of disturbance inputs r.(t), i = 1,..., k. The disturbance 1 inputs are assumed to be independent of and therefore uncorrelated with the plant inputs vi(t), i = 1,..., m. However, the disturbance inputs may be correlated with each other in which case the correlation matrix or equivalent cross spectral density matrix is assumed to be known. The maximum allowable steady state variance of each output E(yi)2, i =,..., m caused by these disturbances with zero inputs v.(t) 0, i = 1,..., m is specified. 7. Known parameter variations given by the triples (6A., 6Bi, 6C.), i 1,..., k which are explained in detail later in this chapter. 8. Allowable static and dynamic cross-coupling due to parameter variations. Restrictions on static cross coupling are given as limits on the magnitude of the off-diagonal elements of the static sensitivity matrix M defined in Section 3.3.2. Dynamic cross coupling is specified by limits on the allowable magnitude of steady state output variances

77 EL2] when the inputs v.(t), i = 1,..., m are excited one at a time by their normal input functions. 9. Allowable static and dynamic gain variations due to parameter variations. Static gain variations are specified as limits on the magnitude of the diagonal elements of the static sensitivity matrix M. Dynamic gain variations are restricted by requiring that the error performance specified in item 5 above be satisfied for all parameter variations specified in item 7.

78 4.2'Reduction'to State: Variab'le Feedback Problem We determine if decoupling is possible, and if not, suitable dynamic compensation is determined, and an augmented system S = {A, B, C}is formed. The computer program described in Gilbert and Pivnichny EG3] is used to determine all of the fixed data given by Decoupling Theorem 2.3. In particular, the location in the complex plane of any "fixed poles" are calculated and examined. Because these poles appear in the "unobservable part" of the closedloop decoupled system, they have no effect on the output, and the information above gives no guide to evaluating whether their locations are acceptable.' Obviously from stability considerations, those which lie on or to the right of the imaginary axis are unacceptable. In addition, although no specific rules are given in Chapter 3, any poles located close to the imaginary axis in the left half plane should be examined carefully because parameter variations given by item 7 above ray cause their location to shift onto or across the imaginary axis. Furthermore, if a complex pair of roots is located near the imaginary axis, the corresponding "part" of the closed-loop system will have a highly oscillatory response which is usually undesirable. Once the indesirable fixed poles are'identified, then additional compensation is determined so that only the acceptable fixed poles remain fixed, i. e., all other poles of the augmented system can be altered by state variable feedback while maintaining the decoupling requirement item 2 above.

79 4. 3 Application of'Opt'imization'Procedure In order to apply the results in Section 3.1, an input filter is selected which when driven by zero mean Gaussian white noise produces a signal similar in magnitude and spectral content to the information given in item 3 above. Initial weighting matrices Q and R are chosen from the error performance and control magnitude constraint information. If the plant has been augmented by a pre-compensator, as required in the previous section, then if the results of Section 3.1 are applied directly, the matrix R will weight the inputs u to the compensator rather than the desired control inputs u to the plant. See Figure 4.1. A modification in the definition of the matrix W defined by (3.6) is used in this case so that the proper variables are weighted by the matrix R and appear directly in the cost functional J. The details of the modification appear in Appendix B. Then we define for the augmented plant S ={A,B,C}. W -C C - ~ v. (4.1) KF + K K G C 1 2 1 v where {F, G} are a decoupling pair for the plant S and K2 EO K2] is an m x n matrix. The cost functional J is minimized by a computer program using the algorithm described in Section 3.1.5. When a minimum is reached, the error variance matrix Eree'] and control variance matrix E[uu' ] can be examined by having the

compensator pldnt V, ~,,X cK I,~J I I 4 Figure 4.1 Compensated Plant

81 computer print the matrix WPW' evaluated at the optimum values of the variables ai, i = 1,..., p + m. The error variance matrix will be given by the upper left hand m x m partition and for a decoupled system will be diagonal. The control variance matrix on the''other hand will be'given by the lower right hand m x m partition and in general will not be diagonal, but since it is a variance matrix, it must be positive definite and symmetric. By examining the diagonal elements of these variance matrices we'can compare their magnitude with the allowable error and control magnitude specification above and decide whether another minimization should be run with stronger weights on certain error or control elements as discussed in Section 3.1.3. Because of the rature of the problem, particularly the lack of existence and uniqueness results, it is not possible to know if other weights will allow all error and control constraints to be met simultaneously. That is, there is no way to tell whether a particular set of constraints represents an impossible design problem or not. Although this lack of knowledge represents a weakness of the design approach, in practice the difficulty is most likely to occur only for systems which allow only small control amplitudes yet require exceptional error performance.

82 4.4 Investigation of Disturbances Disturbance inputs given in item 6 above are examined using the results of Section 3.2.1 to determine whether they will disturb the output. If so, we determine whether their effects can be eliminated by a proper choice of the ai variables, determine what conditions on the variables will cause this elimination and then decide whether to alter the variables from the optimal calculated values, so that they satisfy the complete cancellation conditions. In some cases, a small shift in the ai variables may be sufficient. However even a small shift in the a. might make large changes elsewhere, particularly in the control variances, so these must be checked with the shifted a. values. Another possibility is to do a constrained minimization design with the a. variables constrained so as to eliminate the disturbances. In that case, the basic algorithm of Section 3.1.5 which utilizes a Newton procedure cannot be used, and another procedure would have to be devised. An evaluation of the effect of disturbance inputs which are not eliminated can be made in the following way. Assume an optimal set of a. variables has been determined 1 without considering the effect of disturbances. Select a disturbance input variance matrix rn and disturbance filter given by An, Bn, Cn, whose output r(t) has approximately the same amplitude and spectral content as the disturbance inputs in item 6. Following the development in Section 3.2.2 evaluate the covariance matrix WnPWn with zero inputs, i.e.,

83 set r 0, at the optimum.i values. This matrix represents the output and control variances, similar to WPW', due to the disturbance inputs alone. In general, that portion of WnPW'n representing the outputs will not be diagonal as before. A scalar measure of these variances' can be defined, as before Jn =E [ Q] [on [ O R[n] (L4..2) where Qn and R are positive definite symmetric m m matrices, and J is calculated by the equivalent trace operation. n J trEW'Q WP Q (4.23) n on Such a scalar measure can be used to compare various designs for their response to disturbance inputs. In many cases, we shall let Q= Q, R = R, although this is not necessn n ary for comparison purposes. If it is decided that the cost J is too large, then the disturbance input effects can be included in the minimization procedure as described in Section 3.2.2. As mentioned there we assume that the disturbances and ordinary inputs are independent so that J and J are additive when both inputs and n disturbances are present.

84 4. 5'Che'ck' on Parame'ter' Variations The effects of known plant parameter variations are examined and compared to the information in items 8 and 9 above to determine whether the closed-loop design has acceptable sensitivity performance, that is, we calculate what effect a variation in plant parameters will have if it causes the A, B, and C matrices to change by 6A, 6B, 6C respectively from their nominal values Ao, Bo, C. o A A + 6A 0 B = B + 6B (4.4) o0 C C + 6C 0o A suggested method of determining the triples (6A, 6Bi, 6Ci ) for parameter variations is the following. The linear model (2.1) is usually obtained from a set of nonlinear differential equations which are linearized about an "operating point" the static solution, which in turn depends upon the parameters of the process. Suppose there are k of these parameters. Then we perturb each of them one at a time, i = 1,..., k by a fixed amount indicative of the anticipated variation, compute a new static solution, linearize about it, and determine which elements of A, B, and C have changed from the nominal A, B, C. For i 1,..., k the triple (6Ai, 6Bi, 6C ) represents this change. The effects of these variations are calculated by the following procedure which is repeated for i = 1,..., k.

85 The steady state sensitivity matrix M defined in Section 3.3.2 is calculated and its elements are examined. If any of the elements are too large compared to the data in items 8 and 9, then the gradient of this matrix with respect to each of the (free) (i variables is calculated from Theorem 3.5. It may turn out that a small shift in one or several of these variables will reduce Ithe magnitude of the unacceptable element of M without causing an unacceptable increase in the error or control magnitude variances from the nominal Ao, Bo, Co, which will have to be recalculated. What effect any change in the ai variables has on previously obtained disturbance input calculations is determined by repeating them for the new values of the m. variable s. If this procedure does not result in an acceptable design, then the addition of integral control should be considered. The results of Section 3.3.3 apply. The effect of a parameter variation (4.4) on the dynamic performance of the closed-loop system is evaluated by calculating the cost J(3.7) and the variance matrix WPW' for the perturbed A, B, C matrices. The change in the cost and the variances can be evaluated using the information in items 8 and 9. In general, the perturbed system will not be decoupled and the magnitude of the off diagonal elements of that part of WPW' which represents the error variances, gives a measure of the severity of cross coupling introduced by the parameter variation under consideration.

86 If either the static or. dynamic effects due to parameter variations are unacceptable, then the possibility of complete'cancellation of parameter variations should be checked using the results of Section 3.3.1. Because of the rather stringent conditions for this cancellation to occur, however, it is rot likely that they will be' met. Nevertheless it is worth trying; the calculation is not very difficult. Other than complete cancellation, there is at the present time no systematic procedure for altering the closedloop system while maintaining decoupling if the dynamic interaction caused by parameter variations in unacceptably large. In theory, a scalar measure of dynamic interaction can be defined in terms of these off diagonal elements and this scalar measure can be added as an additional term in the cost functional J (3.7). A new minimization problem can then be formulated which includes this additional term in the cost. It is not clear, however, what the functional form of this scalar measure should be. Even if one is defined, the problem of determining expressions for the gradient and Hessian of J is a formidable task.

87 4*. 6 Limitations of Desi'gn Pr:oc:edure The lack of any clearly defined procedure to handle dynamic interaction due to parameter variations in decoupled systems is one of the weakest areas in the overall design procedure. Another weak spot is the lack of general existence and uniqueness results for the minimization procedure. Indeed it can be shown' that for the first order example'of Appendix C, with certain choices of parameters, no minimum of J exists such that the system is stable. Consequently, a certain amount of good engineering judgement must be exercised when using the procedures just outlined. If one permits some interaction for parameter changes, then an argument can be made that it should be permitted in the choice of F and G for nominal system. In answer, one can argue that there is at present no other way to guarantee that the interaction will be small without resorting to large control amplitudes as are required by the quadratic tracking problem, see Appendix A. Use of decoupling matrices for F and G also very greatly reduces the number of free elements in F and G which have to be selected. This is an important computational advantage.

CHAPTER 5 COMPUTATIONAL QUESTIONS In this chapter, we examine computational questions which arise when the techniques of Chapters 3 and 4 are implemented on a digital computer. The main considerations are core storage requirements, run time, and round off errors. The minimization procedure of Section 3.1 was programmed in Fortran IV language. Extensive use was made of the subroutine call feature in order to keep the programs flexible. One subroutine, for example,' called INNER computes the trace of the product of two matrices, another called SOLVE, solves the Lyapunov equation A'P + PA = -Q. For matrix manipulation, extensive use is made of the well known Fortran IV SSP subroutine written by IBM Corp. [Z1 ]. Most of the computing was done on the University of Michigan computer system which consists of an IBM System 360 Model 67 time sharing system with dual processors. For comparison purposes, this system has an effective internal CPU cycle time of 0.2 microseconds. When references are made later to run times, they were determined for this computer configuration. Some of the examples in Chapter 6 were run on an IBM System 360 Model 65 computer. As we shall see, because of the relatively large number of computations involved, it is necessary to use a large computer such as one of these two just mentioned in order to keep run times down to reasonable levels for examples of order greater than 5. 88

89 5.1 Input Data - Decoupling Program Ressults Input Data for the minimization program consists of the following: System Data - n; m; pi, i = 1,..., m + 1; A, B, C, D A,, i =,..., m, k- 1,..., Pi; G., i - 1,..., m Input Filter Data - v; A; B; C V V V Weighting Data -Q; R; r Starting Data - aik, i 1,..., m; k = 1,..., pi Stopping Data - maximum number of steps; norm of gradient For each example, we assume the integers n, and m and the matrices A, B, C are given. The additional system data is determined by the use of the decoupling program described in [G31 and E 1. Filter and weighting data are selected as previously described. For each example, these determine the final closed-loop design, consequently a capability for repeating minimization runs for different values of Filter and Weighting data is included. Because of the absence of uniqueness results, capability for repeated runs with various starting data is also included. Stopping data provides some control over the total run time. Some of the largest problems may require up to

90 30 seconds CPU time per step. For these problems, a good practice is to allow only two or three steps for the initial run, then stop and observe the results to determine if everything is going all right. If so, then iteration is continued starting where the previous iteration left off. There'are two conditions for which intervention should be made to change the starting data. 1. A step causes the closed-loop system to be unstable for the next iteration. In this' case the procedure will be invalid theoretically for the next step because equation (3.4) may not have a constant solution for the covariance matrix P(t) and even if it does, J is not meaningful if A is unstable. It sometimes occurs that if the procedure is continued anyway, the next and every succeeding step results in a stable closed-loop system. This is pure chance however, and cannot be counted upon to occur. A reasonable procedure is to check for stability at each iteration. When a step causes instability at the next iteration, then the step size is reduced by 1/2 and another stability check is made.' Repeated 1/2 step size reductions are taken until a stable closed-loop system results. Stability determination in general involves finnding the roots of the denominators of each subsystem and checking to see that all have negative real parts. Although this can be programmed,

91 it turns out that most decoupled systems have first and second order subsystems which' are known to be stable if all the coefficients a,. are negative. This condition is easily checked manually, particularly since instability is most likely to occur on the'first iteration, and manually teducing the step size just once has always been sufficient for all examples tried. 2. The cost J increases from its value at the preceding iteration. This increase is known to be a possibility for the Newton method if the initial point is not sufficiently close to a local minimum. If it occurs, then again the usual procedure is to take steps of 1/2, 1/4, 1/8,... size until a reduction in J is attained. For the examples described later, an increase very rarely occurred, and when it did it was always on the first step. Reducing the initial step to 1/2 size by manual intervention was sufficient in every case, so the convenience of a programmed step size reduction was not found to be necessary.

92 5.2 Numerical Methods for Solution of Lyapuhov Equation A- -A Solution of the Lyapunov equation A'P + PA -Q requires a longer computer time than any other step in Figure 3.1. For this reason, we shall go into some detail on just how this solution is determined numerically. Formation of ABI.G - The direct method used by- the author. Although Fortran computer listings have appeared by Chen and Shiek ECl], and Bingulac EB21 which claim to be efficient algorithms for the formation of ABIG from A, the direct method to be described here is superior to either. The program coding is shorter, fewer interations are required, and in the author's opinion, the direct method will execute faster on practically any computer. Consider the j, kth element of the Lyapunov matrix equation. A A n A _ n A,[A l.[ PI + [P]Il A] = [-Q] k (5.1) 1=1 j, lk 11 1, 1 k,l In converting from double to single subscript, this th scalar equation will be represented by the s- element of the vector equation. ABIGP =-q where we choose is k (k - 1)/2 + j, that is P [P]j.k qs- [ Q]j,k

93 A A The coefficients [A]j..l [A 1= *..., n,l k 21,l will appear in the s-th row of ABIG. The proper column is A A selected so that the elements [A].j and [Alk 1 will be coefficients of [P k and [Pj1 respectively when the product th of the s- row of ABIG with the vector p is taken. Now [P1] k and [P]I 1 are the elements i and i2 of the vector p where i j(j - 1)/2 + 1 if 1 < j 1(1- 1)/2 + j if j < 1 (5.2) i2 = 1(1 - )/2 + k if k < 1 k(k - 1)/2 + 1 if 1 < k Matrix ABIG is initially set to zero. Then for each 1 = 1..., n;k = 1,..., n; j= 1, A A ~..,k [A]k 1 and [A]. are added to [ABIG], i and (ABIG),i2 This procedure constitutes the direct algorithm used by the author. The number of iteration steps required is n2(n + 1)/2 each with two additions and the subscript arithmetic (5.2). Also there are (n(n +- )/2) iterations required to zero ABIG initially, but these are executed very fast. On 2 the other hand, Chen's algorithm requires (n(n + i)/2) iterations all with subscript arithmetic. Bingulac's approach requires (n(n + 1)/2)2 iterations to zero ABIG' n3 iterations with one addition and some subscript arithmetic, and n2(n + 1)/2 iterations with one multiplication and some subscript arithmetic.

94 Inversion,of ABiG BIG Inversion of ABIG is handled as a straightforward matrix inversion problem. Because of the large size of ABIG(91 x 91 for one example), a very fast inversion routine must be used if the minimization procedure is to be numerically feasible. A routine INV written by Miss Sarah Lightfoot of the University of Michigan Computer Science Staff in IBM System/360 Assembly Language was selected. This routine is considerably faster than the SSP subroutine MINV. With this or any matrix inversion routine using the Gaussian elimination method, the processing time t required increases as the cube of the dimension N of the matrix [F3]. N3 t =a —- (5.3) 3 because the inner loop of the routine must be executed approximately N /3' times. For the routine INV executing on -6 the computer configuration described earlier a = 15.8 x 10 seconds. For example, less than four seconds of CPU time is required to invert a 90 x 90 matrix. Because of the good results using routine INV, no attempt was made to determine whether there are any inversion techniques available which can take advantage of the sparce nature of ABIG. Perhaps some saving in computer time or storage requirements is possible.

95 The iterative techniques for solution of the Lyapunov equation of Davidson and Man [D1], and Konar EK3] suffer from a lack of any information on how the required number of iterations increases with dimension n, or what a reasonable stopping criteria might be. Both'authors indicate that the convergence'rate'of their algorithms is A strongly dependent upon the location of eigenvalues of A. Konar also proposes iterative improvement EF3] which involves additional iterations. Further research into these techniques may disclose a distinct computer time advantage over inversion of ABIG for large systems. The reduction of core storage requirements, although less important, is well defined [D1].

96 5.3'Run Times The actual CPU time required for the execution of each step in Figure 3.1 was determined by the use of a timer subroutine. Two fairly large examples show the effect of dimensionality on execution times. Example 6.3 has dimension n = 8, m = 3, v - 3, thus N = 66 whereas example 6.4 has n = 9, m = 4, v = 4, N = 91. A summary of the times required is shown in Table 5.1. Table 5.1 Execution Time For Steps of Minimization Program Distillation Step Column Helicopter 1. Read data and print it.790 sec. 1.010 sec. 2. Compute A and W.170.216 3. Solve Lyapunov equation 4.5 8.0 4. Compute cost J.241.415 5. Compute one element of gradient vector.717.870 6. Compute one element of Hessian matrix.975 1.429 7. Invert Hessian.05.05 Total for one iteration 17.0 sec. 27.0 sec. The time reported for calculation of the gradient and Hessian is merely the time required for calculation of a single term of either. Consequently, the total time required to calculate the gradient is p times the figure shown, where p is the number of variables. For examples 6.3 and 6.4, p 4, and 5 respectively.

97 Because of its symmetry, Hessian calculation requires p(p + 1)/2 times the figure shown for calculation of an individual element. However if we can make use of the separation results in Section 3.1.-6, then only elements of the'non-zero blocks of' H need to be calculated. Since each' of the blocks is symmetric, there are m p.(p + 1)/2 (5.4) i=1 i i elements to calculate. For instance, with this simplification only five elements of H rather than ten need to be calculated at each iteration for example 6.3. Inversion of H, although not a large contributor to execution time, is somewhat faster if the individual blocks are inverted, rather than the entire'matrix. A reduction in the time required to compute the cost J may be possible if the following approach is used. From (3.8) which is repeated here for reference, J = tr[W'Q WP], the matrix product W'Q W is computed after each new step 0 because W is a function of the search variables oi, i = 1,., p. Then the matrix inner product or trace of W'QoWxP is taken. The alternative to forming W and then computing the matrix product W'QoW is to compute the product from the equivalent expression

98 P P W'QoW =zi > o.cx T* T(5.5) w1 j-1 i j' Where the constant T.. matrices would be determined just once,, 13 for each example. This same approach can also be used when computing the gradient and Hessian elements. Whether or not it requires less CPU time depends upon the integers m, n, v, and p. The savings if any did not seem to justify the added coding effort required to initially form T., i 1,..,p j = i..., p or the additional memory required to store them.

5.4 Core Storage Re quirre'd an'd Nume rical' Accur'acy A large quantity of core storage is required by the program written. Table 5.2'lists all of the arrays which are stored along with the' maximum dimensions specified in this version of the program. Table'5.2 Matrix Storage Requirements Matrix Size Storage A 10 x 10 400 bytes B 10 x 5 200 C 5 5x 10 200 F 55 x 10 2200 G 5 x 5 100 Av 10 x 10 400 Bv 10 x 5 200 C 5 x 10 200 v Q 5 x 5 100 R 5 x 5 100 r 5 20 avector 10 40 Xvector 5 20 A A 20 x 20 1600 W 10 x 20 800 P 20 x 20 1600 Pi5 i = 1,...,p 10(20 x 20) each 16000 Wi, i - 1,...,p 10(10 x 20) each 8000 Aj, i =1,...,p 10('20 x 20) each 16000 Total.1'8','180' bytes

100 In addition to all of.these, A must be stored. BIG If the inversion is done in double precision and ABIG is for instance 90 x 90, then 64',800 bytes of storage will be required, for this matrix alone. The inverse is developed in the same storage locations. With different coding, small reductions in storage may be possible, some at the expense of longer computational time, such as not storing the Pi. Wi, and A. matrices and generating them when necessary, in general, these figures accurately represent the large quantity of storage required by the' minimization procedure design method for moderately large control problems. All of the computations except inversion of ABIG are done with single precision accuracy. For the IBM System/ 360 computers, six hexadecimal or about 7 decimal significant figures are used. With this precision, no roundoff problems were encountered on any of the examples run. The searching nature of the minimization procedure has the ability to avoid excessive accummulation of small errors. For example rather than using data from previous steps, a new Hessian matrix is computed at each step. The inversion of ABIG is done in double precision simply because the inversion subroutine INV was written this way and probably would not run any faster in single precision anyway. For comparison, example 6.3 which has a 91 x 91 ABIG matrix was run using the single precision version of subroutine

101 MINV to determine in a practical sense what differences would show up. In addition to the greatly increased time required by this subroutine, some elements of the inverse matrix differed from that produced by INV, generally in the last or last two significant digits for the larger (in magnitude) elements. Elements with large negative exponents representing what would be zero elements generally had more negative exponents for the double precision inverse. The overall effect of these differences on the actual calculation of cost J, gradient vector, or Hessian matrix was nevertheless very slight. Both approaches reached the same minimum location to within five significant figures in all search variables, in exactly the same number of steps.

102 5.5 An Efficient Program The main concern in writing the' computer program has been to get some examples computed and demonstrate feasibility of the technique. A well designed program requires more effort and it may be asked what features a highly efficient computer program would have if we assume minimal computer time is our goal rather than minimal storage. In order of decreasing importance, the following five changes should result in a considerable reduction of overall time. 1. Separation into individual sub-problems (Section 3.1.6). The savings depends on how the problem separates. A tenth order problem with m 2, p - 1, P2 = 2, P3 7 separates into an 8th and a 9th order problem and there probably is no savings at all. On the other hand a tenth order problem with m = 3, P1 = P2 = p3 = 2 p4 = 4 th separates into three 6 order problems and there should be considerable savings. 2. Improved method for solution of Lyapunov equation. Because this is the most time consuming step, a reduction here would be significant. For problems of order 10, even if this time were reduced to zero, the maximum savings would be about 50%. The most likely prospects for improvement are the results of Man EM31 and Konar EK3], and the separation procedure of equation (3.23). 3. Re-organization of calculations. The technique of (5.5) can be used to reduce the time required by steps 2, 4, 5, 6 of Table 5.1.

103 4. Programming in assembly language.' Although Fortran is a convenient language'to use it generally does not produce the most efficient machine code for array calculations. Re-programming some'of the. computation steps or the entire program in assembly language can cause a significant reduction in overall computer time. 5. Improvements to the search algorithm itself. There are many possibilities here, such' as fixing the Hessian after a reasonable number of iteration steps. Since only about ten iterations have been required for the examples of Chapter 6 improvements of greater than 50% probably are not possible.

104 5.6 Max'imum Prob'Iem''Size The question of just how large a problem can be solved using the minimization procedure is best answered as follows. Assume the problem is separated into m subproblems then the question is what is the largest subproblem which can be'solved? Three basic'limitations on problem size are computer run time, storage size, and roundoff errors. Computer time limitations are primarily a matter of cost. Using the most time consuming step in the iteration procedure, inversion of ABIG, as a guide in determining run time per iterath tion we find that the cost increases as the 6 power of the number of states n in the'subproblem. For n = 20 about two minutes of time per iteration is required. Considering that ten or more iterations may be necessary and that we might wish to repeat the procedure a few times with differing weights, several hours of computer time can easily be consumed. At $7.00 per minute current rates, the cost may be prohibitive. For n = 10, only 1/64 of this figure, or just a few seconds is needed. If the Lyapunov equation is solved by some other faster technique, then larger problems may be solved. However no figures' are available on computer time for any other method, or indeed what is more important, how computer time increases with dimension n.

105 At the present time, unless one is willing to spend a lot of money, computer time limitations appear to impose a problem size limit somewhere between n=10 and n=20. In modern computers, core storage is large. Nevertheless it should be noted that a n=20 problem will require on the order of half a million bytes. While this is easily within the reach of today's large computers and time sharing medium size computers, not everyone has such a system available. The roundoff error limitation is difficult to determine because it depends in an extremely complex way on the problem data itself and how it is scaled. As a comparison even the effect of scaling on the solution of a linear system of equations is poorly understood at this time, see [F3 p. 38]. Current practice is to scale the rows and columns of matrices by powers of 10 so that the elements with largest absolute value in each row and each column are about equal. Because this problem is so complex, experience with example problems probably gives the best answer. The examples of Chapter 6 with n=8 and n=9 were run successfully which indicates that problems with subsystems of this order are possible. To date these are the largest problems for which decoupling feedback matrices have been calculated. Indeed it may turn out that computing the solution of the decoupling problem may be more difficult because of roundoff errors than design of the closed-loop system itself.

106 CHAPTER 6 EXAMPLE PROBLEMS A variety of examples were studied and they are discussed in this chapter. A second order example representing a motor-generator combination was used for a considerable amount of investigation concerning effects of filter and weighting parameter selection. Also, to demonstrate versatility of the optimization procedure, it was applied to increasingly higher order systems including: 1. A fourth order aircraft landing problem with two inputs 2. An eighth order distillation column with three inputs 3. A ninth order helicopter with four inputs 6.1 Motor Generator Example Cruz, Perkins, and Gonzales [P1] develop state equations for a dc to ac rotary converter. The converter, shown in Figure 6.1, consists of a dc motor driving an ac generator. Their equations are linearized about a nominal operating point yielding the linear system described by the A, B, and C matrices shown below where Ul - dc motor field current U2 - ac generator field current

107 + I I + V V'V V U. = V- Vo Y, = (- (o 2 igure 61 Figure 6.1 Mot or Generator Example

108 y = x - shaft angular velocity 1 1 Y2 x 2 - output voltage -4 - 4 0 A - 2 B I | [ C 2 -2 _44 v -4 2 0~[ 11 The control system is required to be decoupled. It is known that the inputs will be unit step functions. Furthermore, the static operating values are subject to variation causing a change in the parameters of the linearized model. By use of the Gilbert and Pivnichny computer program [G31, it is determined that this system can be decoupled into two first order subsystems with single input - single output transfer functions. 1 hi(s) = h2(s) s -.2 2 s- 2 The corresponding decoupling pair {F, G) of control law (2.2) is chosen from the class. F = -1 -.25a -.5 G.25X 0 -1 -.5 1 1 +.S a.5X.51 L1 2 For simplicity, we assume a unity static gain constraint as described in Section 3.1,thus 1 -a-' 2 -a2'

109 A number of minimizations were computed with various filter and weighting parameters in order to get a "feel" for their effect on the problem. This example was used for these investigations because the order is low and run times are short; less than 0.4 seconds of CPU time is required per iteration. As we expect, the number of iterations required to reach optimum is strongly dependent on the starting values selected for the search variables. Nevertheless for this example the procedure converged (norm of gradient -5 < 10 ) in less than ten iterations for starting values which lie within one or two orders of magnitude of the optimum. It is not difficult to select starting values within these limits. If the open loop poles are known, then it may be assumed for most designs that the closedloop poles will not lie many orders of magnitude away. For if they did, then it is likely that the design would have unreasonably large control or error magnitudes and in either case would be an unacceptable design. While there may be exceptions, for instance in poorly scaled problems, (see Section 5.6), "this rule of thumb" of selecting starting values so that the closed-loop ~2 poles are within 10 magnitude of the open loop poles has worked well for all examples tried.

110 If starting values happen to be selected close to the optimum, such as using the optimal values of aik from a previous entry in Table 6.1 as the starting values for the next entry, then convergence is very rapid requiring only two or three iterations. Table 6.1 shows the effect of increasing weights on the controls while maintaining the weight on error variances fixed. Identical first order filters are used for these designs with break frequency wo of one radian per second. That is Av v v 0 _W L0 2w L The factor of two in the elements of B is used v to normalize the filter gain so that its output variances are unity when driven by zero mean Gaussian white noise with unit variance, r = diag (1,1), E(vi) = 1.0, i 1, 2. Notice the tradeoff which can be made between control variances and error variances by suitable adjustment of the weighting parameters. Because the closed-loop system is being asked to follow input signals with a bandwidth similar to the bandwidth of the open loop system (open loop poles are located at s = -2, -6), both error variances and control variances cannot be made small simultaneously. For Table 6.2, the weighting matrices are held constant and w is varied thus changing the bandwidth of the signals which the system is asked to follow. Here we see

111 the error variances decrease as the filter bandwidth is decreased without any significant increase in control variance. Some minimizations were run with second order filters to see what effect their more rapid falloff of high frequency gain would have. The results are shown in Table 6.3. Identical second order filters with no zeros, and poles as indicated were used for each input element. As with the first order filter case, filter gains were normalized so that E )~2 2 2 E(v E(v2) 1.0. It appears that the error variances are sensitive to high frequency components in the input. Use of the underdamped filters requires considerably more control effort for a similar error variance when compared to the critically damped filters, probably because its output has a greater proportion of high frequency components. In general, the second order filter designs have much less error variance than the first order design with comparable bandwidths. This is not surprising since first order filters put much more energy into very high frequencies.

112 Table 6.1 Effect of Weighting Parameters on Error and Control Variances Motor Generator Example, First Order Filters Break Frequency of 1.0 Radian/Second Control Control Total R Matrix Errors Effort Effort Control Variables with v2=0 with vl- Effort E(1-y)2 E 2 E() 1 E(v - (u u E(u ) a [ 1 2 2 2 2 el 0.441 1.01.458 1.12 -4.158. 0..409 1.28 1.37 1.88 -5.0.2 0.475.97.456 1.07 -3.435 0.2.399 1.20 1.36 1.82 -4.916.680.76.430.87 -1.166 L 0.2.512.835 1.12 1.40 -2.808 Note: Q = Table 6.2 Effect of Filter Bandwidth on Error and Control Variances Motor Generator Example, First Order Filters Break ~Break 2 2 - 2 2 Frequency E(vl-Yl)2 E(v2-Y2) E(u2) a1 a2 Rad/Sec 1.0.1939.1667 1.2239 3.5190 -4.1575 -5.0.5.0984.0842 1.2598 2.9565 -4.5843 -5.4371.2.0395.0338 1.2602 2.4366 -4.8636 -5.7194.1.0198.0170 1.2564 2.2278 -4.9606 -5.8170 Note: Q = R =[ for all cases above Filter gains are selected so that E(v1)2 E(v2) 1.0

113 Table 6. 3 Designs Using Second Order Filters Motor Generator Example 2nd Order Filter Poles at -1, -1 R Matrix E(v -Y ) E(v2 y )2 1 E 2 a 11 2 20 1 00071.00039.65314 1.2329 -25.545 -34.557 0. 1.2 0.00446.00196.6457 1.2060 -9.595 -14.969 0 2 3 0 [.08743.02583.54168 1.0037 -1.391 -3.400 0.5 2nd Order Filter Poles at -.707 ~ j.707 [: 2i.0056.0021 1.2997 2.4495 -13.002 -20.770

114 6.2 Aircraft Example Another application of this technique arises in the design of an aircraft system for decoupled control of pitch motion and airspeed during approach to landing. This problem is considered by Smith [S3] for a typical KC 135 making use of the independent elevator and engine thrust variables. The linearized model is a fourth order system described by the matrices 0 1.0 0 0.00162 -1.2358 -1.1179.07281 A. 00550.95810 -.64900 -.24630 - 12736 0.07280 -. 03323 0 1 -1.00800 6.72000 1 0 0 0 B = -.03419 0 C 0.79210 O 0 0 1 with Y1 X - pitch angle in radians 2 -angle of attack, radians 3- angle of attack, radians y x - perturbation velocity, fraction of nominal 2 4 uI - elevator angle, radians u2 - engine thrust, 106 lb. The aircraft is assumed to be on a 2.5 degree glide path, at an airspeed of 253 ft./sec., with its gear down and 40 degree flaps. From the decoupling program [G3], it is determined that the system can be decoupled by state variable feedback into the subsystems with transfer functions:

115 h (s) = Al 2 s - a11S - r12 h2(S) =S - s - a21 there is one fixed root at s = -.59. The corresponding class of decoupling pairs - {F, G} are given by: F 1.07031 -.99206a12 -1.22599 —.99206a11 -1.72175.16079 0.09191.35191 + 8.41645a.04195 + 1.2624a21 2 0O 1.26247A2 Because of this structure, any disturbance inputs which enter only the unobservable "part" of the closed-loop system will have no effect on the output. From the results of Section 3.2.1 we find for the CD system A 0o v = [0, V31 a scalar is free and transforming back to the original system v = T2 v where T2 1 is given by the decoupling program. = O v 31

116 In addition, from Section 3.3.1 we can state that any parameter variation of the form: 0 0 0 0 SA- 0 0 0 0, Sa. i 1,..., 4 free 0 0 0 0 6a1 6a2 6a3 6a4 will also have no effect on the output. One must however, check the poles of the closed-loop system for such a variation, to determine that the "fixed" root at s = -.590 has not moved into the right half complex plane, indicating an unstable system. Smith does not give information concerning the type of inputs the aircraft will be expected to follow. We expect pilot inputs to contain significant high frequency components and arbitrarily selected first order filters with 0.2 radian per second break frequencies. Although this may seem very low, recall that first order filters put a significant amount of energy at frequencies well above the break frequency. The open loop system has poles located at s = -.944 i j.993, -.0144 ~ j.147. Filter gains are selected so that E(vi)2 E(v2)2 = 1.0. These input amplitudes are chosen merely for convenience in comparing error and control amplitudes to input amplitudes rather than to represent the actual signal amplitudes anticipated. As we have shown in Section 3.1.6, the problem can be separated, with the effect of each input considered individually. From linearity, the amplitudes scale.

117 With inputs characterized by these filters, and weighting matrices given by: 1 0 1 0 Q = 0 1 R.1 a unity static gain constraint 1 = -a12' 2 = -21 is imposed. The minimization procedure converged in eight iterations of 1.9 seconds CPU time each, yielding: (v1 - 1)2.2475 E(v2 - y2 = 4876 E(ul) = 4.1045 2 E(u2 ) =.0569 11 - = -4.24913 12 = -3. 34053 a21 = -.21015 It appears that large elevator deviations will be required to follow such inputs with this design. Separating the control variances by input, we find both inputs contribute strongly to the u1 variance. due to vl due to v2 E(ul)2 2.0573 2.0472 E(u2 )2.0199.0370

118 If these are unacceptably large, then the design of course can be repeated with a stronger weight on the u1 variance. Some typical amplitudes as used by Smith are ~10 elevator angles ~=.0175 radian and 2 000 lb. thrust =.002 for u2. If we assume Ev12 =.0175 then the following error and controls results. vl _ yl2 =.5o E(v1 - 2 50 2 2 E(ul) 1.4 E(u2) = 2460 lb. which seem to be acceptable values.

119 6.3 Distillation Column Example A design for the binary distillation column originally described in Gilbert and Pivnichny EG3] was computed using the minimization program. This eighth order system with three inputs in shown in Figure 6.2. Equations describing the column operation are linearized about an operating point and the resulting system of linear differential equations are represented by the matrices shown in Table 6.4 where: Y1 = xl - concentration of liquid on top tray x 2 - concentration of liquid on tray 3 x 3 - concentration of liquid on feed tray X4 - concentration of liquid on tray 1 Y2 = x5 - concentration of liquid in reboiler X6 Y overflow rate of tray 3 x 7 - overflow rate of feed tray y3 = X4 - overflow rate of tray 1 u1 - reflux ratio u2 - feed rate u3 - reboiler evaporation rate As reported [G3] the system can be decoupled into two first order and one second order subsystems with the control law matrices shown in Table 6.5. hl(s) = 1 S - 0 11

120 CONDENSER V L | T / I, u, XT.IL LXB COLUMN XB B REBOILER Figure 6.2 Binary Distillation Column

121 Table 6.4 Distillation Column A and B Matrices -.09.12 0 01 0 0 0 0.06 -.18.12 0 0 -.134 0 0 0.06 -.28.12 0.134 -.096 0 0 0.16 -.28.12 0.096 -.084 0 0 0.16 -.2 0 0.084 0 0 0 0 0 -.2 0 0 0 0 0 0 0.2 -.2 0 0 0 0 0 0 0.2 -.2,077 0 -.067.154 0 -.033 0.1 -.018 0 0 -.025 0 0 -,099.2 0 0 0.2 0 0 0 0

122 Table 6.5 Decoupling F and G Matrices for Distillation Column 1.17 + 13a11 -1.56 0 1.41 -1.76 - 8.8a21 F = | 0 0 0 0 0 0 0 0 1.62 -2.02 - 10.1a21 0 0.738 -1' 2 + 5a 32 -1 + 25a31 + 5032 0 0.848 13X1 -8.8X2 0 A: = 0 0 25X3 0 -0.1X2 0 2

123 h (s) = 2X3 2 h (s) 3 3 2 s -a31 S -32 Four "fixed" poles, all stable, are located at s = -.08, -.14, -.5~ j.0336. For this pEoblem, first order filters with fairly low break frequencies (.02 rad/sec) were selected because it is known the usual inputs do not have a considerable high frequency content. Filter parameters also were constrained to produce unity output variances when driven by unity variance Gaussian white noise. Unity gain was assumed for each subsystem, A1 -11' 2 = -X21, X3 = -a32 and weighting parameters were selected for overall acceptable control and error variances. With Q = diag (1.0 1.0 1.0), R = diag (0.1 0.1 0.1) the optimization procedure yields the following results which are listed along with the arbitrary choice previously reported in [G3]. Minimization Gilbert S Pivnichny Variables all = -.24388 -.08 a21= -.17286 -.08 31 = -.88988 -.08 a32 = -.13632 -.02 32

124 Mitnimization Ggi'lbert g Pivnichny Control variances E(u2) = 2.58078 2.09422 E(uN) = 2.27182 2.12783 E(u ) = 2.05858 1.96914 Error variances E(v1 Y1).97579.2 E(v2 - y2)2 =.10371.2 E(v3 - y3)2 =.13759.31818 Both designs appear to have reasonable control and error amplitudes. The design produced by the minimization procedure has less than half the error variance with slightly larger control magnitudes required. The sources of control effort requirements vl,v2, v3 were introduced individually and the results found to be. due to vl due to v2 due to v3 E(ul).77191 1.29396.51490 E(u2)2.44937 1.16841.65402 E(u3)2.026 82 1.30868.72311 Neone- of the inputs requires an unusually large control. A source of variation in performance for distillation columns is non-uniformity of the feed composition. Such changes occur slowly compared to the column dynamics and thus fall into the category of parameter variations as discussed

125 in Section 3.3. An increase in the composition of the ore volatile component in the feed of 10% above its nominal value of 0.1 is considered to determine how sensitive the results above are to such changes. A new operating point based upon the higher feed composition was calculated as the steady state solution of the system of nonlinear differential equations which describe column operation. Linearization about this new operating point yielded the slightly perturbed A and B matrices shown in Table 6.6. Note that a change in the operating point of a nonlinear system may be viewed as a parameter variation in the corresponding linear system. Using the perturbed A and B matrices along with the F and G matrices determined as optimal for the unperturbed system and the previously described filters and weighting constants, error and control variances were computed as described in Chapter 4. The computer time required to compute these variances is less than required for one iteration cycle and involves doing steps 1 - 3 of Table 5.1 and then computing the matrix product WPW'.

126 Table 6.6 Perturbed A and B Matrices for Distillation Column -.09.12 0 0 0 0 0 0.06 -.18.12 0 0 -.148 0 0 0.06 -.28.12 0.148 -.106 0 0 0.16 -.28.12 0.106 -.092 0 0 0.16 -.2 0 0.092 0 0 0 0 0 -.2 0 0 0,0 0 0 0.2 -.2 0 0 0 0 0 0 0.2 -.2.085 0 -.0741.17 0 -.0373 0.11 -.02 0 0 -.028 0 0 -.108.2 0 0 0.2 0 0 0 0

127 The magnitude of the cross coupling which is evident from the error variance matrices due to each input separately is not severe as shown below in Table 6.7. Table 6.7 Perturbed Distillation Column Variances Due to v B[(v-y)'(v-y)] =.07130 -.00042 0 -.00042.00007 0. 0 0 0 Due to v2.00012 -.00171 0 -.00171.10125 0 0 0 0 Due to va.00003.00005.00059.00005.00008.00115.00059.00115.1376 0 Controls due to vl due to v2 due to v3 1Ul 2 3 E(u1)2.69483 1.22277.52963 E(u2)2.39262 1.10074.64930 E(u3)2.02584 1.23274.73257

128 In addition, the small change in control variances shows that control magnitude demands have'not altered significantly so that overall, we expect this design is relatively insensitive to the 10% increase in feed composition. The'static gain variation due to this perturbation and the sensitivity matrix M defined in Section 3.3.2 were computed. Recall that the unperturbed system has unity static gain (To = I); we find 1.0119 -.0101.0015 T = -, 084 1.0480.0024 0 0 1.0018 M= f.00119 -.0101.0015 -.0084.00480.002 4 0 0.0018 The elements of M appear to be acceptably small in magnitude so that no additional effort should be required. The closedloop poles for the perturbed system were calculated and all are located in the left half complex plane. The individual subsystems do not have any numerator dynamics and consequently the only possibility for complete disturbance input cancellation occurs if disturbances enter the fixed unobservable part of the system which for this example has dimension 4,. Consequently A = column (0, 0, 0, 0, v41, v42' v143, v44) where v14 - V44 are free scalars. Transforming back from the CD to the original system, we have

129 v = column (0, v1 v4, v44, v43, 0, v42, 0, 0). Following the results of Section 3.3.'1 we may state from this disturbance'input result that parameter variations in the second, third, fourth, and sixth rows of A will not affect the output, provided the system remains stable, and the elements ofB and C remain fixed. Unfortunately the perturbation above causes changes in the elements of B so the theory does not apply to this perturbation.

130 6.4 Helicopter Example Previous efforts by Murphy and Narenda EM8] and Wolovich and Shirley [W2] to design control systems for the Sikorsky SH3D Sea King helicopter provide useful comparisons to our results for this example. The vehicle is required to spend a large amount of the total time hovering at an altitude of 40 feet, maintaining a fixed position with respect to ground. The hover dynamics for small deviations can be represented as a ninth order system x = Ax + Bu with four inputs, where A and B are given in Table 6.8 and x1- longitudinal velocity, normalized x2 - lateral velocity, normalized x - vertical velocity, normalized x4 - roll angle, radians x - roll rate X6 - pitch angle, radians x - pitch rate x 8 - yaw angle, radians x - yaw rate u1 - longitudinal cyclic pitch, radians U2 - lateral cyclic pitch, radians u3 - tail rotor collective, radians U4 - main rotor collective, radians

131 Table 6.8 Helicopter A and B Matrices -.016 -.0047 0 0.0047 -.033 -.0007.05 -.00001.00004 -.3242 0 0 0 0 0 A = 2.61 -7.25 -.163 0 0 0 0 0 1.97.736 1.0 0 0 0 0 0.016 5.59 -.193 0 -.001 -.05.0025 0 0 -.0025 0 -.0024 0.00 09 0 0.00018 0 0 1 0 0 0 0 -1.96 0 -1.94 0.01 0 0 1 0 0.548 0 -.542 0 0 0 0 0 1 -.0043 0 -.0083 0 -.303.05 0 0.005 0.05.022.01 0 0 0 -4.24 0 0 0 0 B = 0 21.81.3475 -2.13 0 0 0 0 -6.15 0 0.69 0 0 0 0 0.174 -7.48 5.12

132 Decoupling is cited as a desirable handling quality by both authors, however, there is some disagreement as to which variables should be considered as outputs. Murphy and Narenda use a quadratic regulator type of approach and choose to weight xl, x2, x3, x8. Wolovich and Shirley use a frequency domain approach and consider x3, x4, x6, x8 as the most important output variables. It was decided to use Murphy and Narenda's outputs Y1 = xl, y2 = x2 3 = x3, Y4 = X8 because their approach has more similarity to the optimization procedure. In fact their choice of weighting coefficients provides the starting point for our selection of Q and R matrices. In their case, inputs are chosen to be solutions of a homogeneous second order differential equation of the form d2Vl dv 2 + al - + av = O dt2 dt or v 0 z 6;, I i-a a V o 111 1 By varying a, a and 9 (o) a wide range of commands is geneo 1 1 rated with the condition of v (0) = O. Although no specific values are given, the shape of their pilot command graphs seems to indicate that they have chosen ao and al to produce overdamped responses with predominant time constants from two to six seconds.

133 We have consequently chosen our input filters to be first order systems with time constants of 50 seconds which should guarantee that the high frequency components of the input are sufficiently attenuated. The similarity of the two approaches ends here because Murphy and Narenda have no way to guarantee a decoupling requirement as we do, and in fact their resulting feedback matrix 2.59 1.19 -0.12 0.10 0 F Narenda - -1.16 2.66 0.0L4 0.20 0.07 -0.05 -0.58 -0.99 -0.01 0 -0. 02 0.14 -2.15 -0.01 0 -0.28 -0.20 0.01 0 0.11 -0.01 0.02 0.01 0. 01 0 -0. 22 -0. 28 0 0 0.09 0.15 cannot possibly yield a decoupled control system. Comparison with the class of decoupling pairs {F, G} given in Table 6.9 verifies this fact. As it turns out, the system considered by Murphy and Narenda can be decoupled into the following subsystems: h (s) = 1 1 1 h2(s) = X2 s - a21 h3(s) = X3 s - <31 h (s) = 4 41 442

134 Table 6.9 Decoupling F and G Matrices for Helicopter.32 + 20 a11.094.0765 +.236a 0 31 -.094.328 + 19.8o.4044 + 1.17a -.99 ~~F -=~~ ~21 31 -.00006.755 +.46a21 -.5398 - 1.59a31 -.023 -.00003.00008 -.7646 - 2.3631 0.002 1.0 -.05 0 0.05 0.048.058a42 -.00017 +.058a.0006 0.0003 -.132a42 -.0405 -.132a41.00001 0.0004 0 0 20RX 0.236,3 0 0 19.8X2 1.17X3.058X G = 0.46X2 -1.59X3 -.132X 0 0 -2.36X 3 0 There are four fixed poles, all stable, located at s = -.437 ~j 4.70, -.117 ~j 2.44. Although these are lightly damped, they are comparable to the open loop poles located at s 1.31 ~ j.651, -.048 ~ j 4.14, -.083 ~ j.316, -.324, -.305,0.0. Using the weighting matrices from Murphy and Narenda Q = diag (10, 10, 10, 0.5), R = I, and filters as described with unity output variances and assuming the unity gain constraint for the system, x1 = -~11' ~2 - -21' X3 -31'

135 4 = -a42 a minimization was started at all = -.1, 21 = -.1, 31 = -.5, 41 = -1.0, 042 -1.0. The procedure converged to the following results in seven iterations. variables al = -.12593 11 Fa = -.18279 21 a31 = -. 7996 3 o41 = -2.62629 a42 = -3.74346 variances E(vl - y) =.137 E(u) 2.774 E(v2 - Y2) = 099 E(u22.487 E(v3 - y3) = 024 E(u) =.852 E(v4 - y4) =.021 E(u4) =.657 Closed-loop system poles are located as shown below in Table 6.10 where they are compared to the closed-loop poles of Murphy and Narenda's and Wolovich and Shirley's designs. Iteration time as reported in Chapter 5 is about 27 seconds without the independent input simplification. Computation of the conditions for complete cancellation of disturbance inputs yields v = column (0, 0, 0, 0, O0, v51, v52, v53' v54) and transforming back to the original system, v = column (0, 0, 0, v53, v54' v51' 5 52' 0, 0).

136 Table 6.10 Comparison of Pole Locations of Closed-Loop Designs Author's Design Murphy & Narenda Wolovich & Shirley -.126 -.460 i j.734 -2.121 ~ j 2.121 -.183 -.822 + j 1.168 -2.121 ~ j 2.121 -.800 -1.002 + j.252 -1.0 -1.313 i j 1.42 -1.522 -10.0 -.437 7 j 4.70 -1.418 0.0 -.117 i j 2.44 -2.455 0.0 0.0 Reports of pilot evaluation by Wolovich and Shirley indicate that "decoupling is a desirable design objective to achieve if one does not sacrifice stability in the process". Although Murphy and Narenda's design is asymptotically stable, it lacks the desirable decoupling whereas Wolovich and Shirley's design is the other way around. The design presented here contains both these features without compromising on either. As such, it represents a very desirable alternative to be considered.

CHAPTER 7 CONCLUSION In summary we have presented a design technique for multivariable feedback control systems which involves decoupling, minimization of error and control amplitudes, and reduction of the effects of disturbance inputs and plant parameter variations. Actually a number of results which apply to these various aspects have been presented and although one way of making use of all of them in design is discussed in Chapter 4, there is open area for improvements in determining ways to best put all of the results together for a unified synthesis method. The error and control variance minimization procedure of Section 3.1 represents the most reasonable approach to the control amplitude factor in decoupled systems where its nearest competitor, the quadratic tracking theory has certain difficulties. We conclude from the discussion of computer calculations and example problems that this minimization technique can be used for the design of fairly large systems without using unreasonably large amounts of computer time. While the basic algorithm and computer program clearly demonstrate feasibility of the procedure, there is ample room for additional study and implementation of improvements such as those identified in Chapters 3 and 5. These are particularly important if one is interested in 137

138 designing systems with very high order (greater than 10) subsystems. A more efficient program is clearly possible and should be a desirable objective now that the method has been shown to work for the examples of Chapter 6. There is room for additional work on better ways to specify compensators for plants which have weak inherent coupling or are unstable under the constraint of decoupling. The approach taken here assumes they have been determined if needed. An important implication of the results of Chapter 3 is that they help to develop intuition in how to handle multivariable systems and in how to tell what is happening within them. The examples of Chapter 6 show how this intuition is used in the application of these theoretical results. Two techniques are presented for handling disturbance inputs. The first involves determining whether the inputs will affect the output. When possible the control law is structured so they will not. These results depend on a part of the system being unobservable, which the experience of Chapter 6 shows is a common occurance in decoupled systems. Otherwise the effects of disturbance inputs can be reduced by including them in the error term of the cost functional J. The approach to parameter variations differs from that developed by Cruz and Perkins [C4], in that it is assumed we have information about how the elements of A, B, and C vary. We attempt to minimize only the effects of these specific variations. Particular emphasis is placed upon static gain

139 variations. Changes in system dynamics due to parameter variations remains a weak and difficult area in terms of the results which are available.

References Al M. Athans and P. L. Falb, Optimal Control, McGraw Hill, New York, 1966. A2 M. Athans and E. Tse,' "A Direct Derivation of the Optimal Linear Filter Using the Maximum Principle", IEEE Trans. on Auto. Control, vol. AC12, no. 6, December 1967, pp. 690-698. B1 S. Barnett, "Comments on "A Theorem on the Lyapunov Matrix Equation'", IEEE Trans. on Auto. Control, vol. AC15, no. 2, April 1970, pp. 279-280. B2 S. P. Bingular, "An Alternative Approach to Expanding PA + A'P = -Q." IEEE Trans. on Auto. Control, vol. AC15, no. 1, February 1970, pp. 135-139. C1 C. F. Chen and L. S. Shieh, "A Note on Expanding PA + A P = -Q." IEEE Trans. on Auto. Control, vol. AC13, no. 1, February 1968, pp. 122-123. C2 E. A. Coddington and N. Levinson, Theory of Ordinary Differential' Equations, McGraw-Hill, New York, N. Y. 1955. C3 J. B. Cruz and W. R. Perkins,"Conditions for Signal and Parameter Invariance in Dynamical Systems," IEEE Trans. on Auto. Control, vol. AC11, no. 3, July 1966, pp. 614-615. 140

141 C4 J. B. Cruz, Jr. and W. R. Perkins, "A New Approach to the Sensitivity Problem in Multivariable Feedback System Design, "IEEE Trans. on Auto. Control, vol. AC9, no. 3, July 1964, pp. 216-223. D1 E. J. Davidson and F. T. Man, "The Numerical Solution of A'Q + QA = -C." IEEE Trans. on Auto. Control, vol. AC13, no. 4, Aug. 1968, pp. 448-449. D2 P. M. DeRusso, R. J. Roy, and C. M. Close, State Vari'ab le's for En gi nee'rs, Wiley, New York, 1966. Fl P. L. Falb and W. A. Wolovich, " Decoupling in the Design of Multivariable Control Systems," IEEE Trans. on Auto. Control, vol. AC-12, no. 6, Dec. 1967, pp. 651-659. F2 P. L. Falb andW. A. Wolovich, "On the Decoupling of Multivariable Systems, "Proc. 1967 JACC, Philadelphia, Pa., pp. 791-796. F3 G. Forsythe and C. B. Moler, Computer Solution of Linear AlgZebraic Systems, Prentice Hall, Englewood Cliffs, N. J., 1967. G1 F. R. Gantmacher, The Theory of Matrices, Chelsea, New York, 1969. G2 E. G. Gilbert, "The Decoupling of Multivariable Systems by State Variable Feedback," SIAM Journal on Control, vol. 7, no. 1, 1969.

142 G3 E. G. Gilbert and J. R. Pivnichny, "A Computer Program for the Synthesis of Decoupled Multivariable Feedback Systems," IEEE Trans. on Auto. Control, vol. AC14, no. 6, Dec. 1969, pp. 652-659. Hl J. W. Howze and J. B. Pearson, "Arbitrary Pole Placement in Decoupled Linear Systems," to appear Il Z. Iwai, "Studies on the Noninteracting Control of Multivariable Systems," Ph.D. Thesis 1970. K1 R. E. Kalman, "Contributions to the Theory of Optimal Control," Bol. Soc. Mat. Mexicana, Second Ser., vol. 5, 1960, pp. 102-119. K2 R. E. Kalman and B. L. Ho, "Effective Construction of Linear State Variable Models from Input/output Functions," Proc. Third Allerton Conf. 1965, pp. 449-459. K3 W. Kaplan, Operational Methods For Linear Systems, Reading, Mass., Addison Wesley, 1952. K4 A. F. Konar, "Improved Iterative Algorithm for Solving Lyapunov Equation," to appear L1 C. K. Liu and N. J. Bergman, "On Necessary Conditions for Decoupling Multivariable Control Systems," IEEE Trans. on Auto. Control, vol. AC15, no. 1, Feb. 1970, pp. 131-133. M1 E. C. Ma, "A Finite Series Solution of the Matrix Equation AX-XB = C", SIAM J. Appl. Math., vol. 14, no. 3, May 1966, pp. 490-495.

143 M12 P. J. Maclane and E. J. Davidson, "Disturbance Localization and Decoupling in Stationary Linear Multivariable Systems", IEEE Trans. on Auto. Control, vol. AC15, no. 1, February 1970, pp. 133-134. M3 F. T. Man, "A Theorem on the Lyapunov Matrix Equation", IEEE Trans. on Auto. Control, vol. AC14, no. 3, June 1969, p. 306. M4 B. S. Morgan, Jr., "Sensitivity Analysis and Synthesis of Multivariable Systems," IEEE Trans. on Auto. Control, vol. ACll, no. 3, July 1966, pp. 506-512. M5 B. S. Morgan, Jr., "The Synthesis of Linear Multivariable Systems by State Variable Feedback", IEEE Trans. on Auto. Control, vol. AC9, October 1964, pp. 405-411. M6 A. S. Morse and W. M. Wonham, "Decoupling and Pole Assignment by Dynamic Compensation", SIAM Journal Control, vol. 8, no. 1, 1970, pp. 1-18. M7 J. H. Mufti, "On the Observability of Decoupled Systems", IEEE Trans. on Auto. Control, vol. AC14, no. 1, February 1969, pp. 75-77. M8 R. D. Murphy and K. S. Narenda, "Design of Helicopter Stabilization Systems Using Optimal Control Theory", J. Aircraft, vol. 6, no. 2, March-April 1969, pp. 129-136.

144 N1 S. Nazar and Z. V. Rekasius, "On the Decoupling of Nonlinear Multivariable Systems and Sufficient Conditions for a Special Class'", to appear in IEEE Trans. on Auto. Control. N2 G. C. Newton, Jr., L. A. Gould, and J. F. Kaiser, Analytical~ Design' of Linear? Feedb'ack Controls, Wiley, New York, 1957. P1 W. R. Perkins, J. B. Cruz, and R. L. Gonzales, "Design of Minimum Sensitivity Systems", IEEE Trans. on Auto. Control, vol. AC13, no. 2., April 1968, pp. 159-167. P2 W. A. Porter, "Minimizing System Sensitivity Through Feedback", Journal of Franklin Institute, vol. 286, no. 3, Sept. 1968, pp. 225-240. R1 Z. V. Rekasius, "Decoupling of Multivariable Systems by Means of State Feedback", Proc. Third Allerton Conf., Monticello, Ill., 1965, pp. 439-448. R2 Z. V. Rekasius, "Optimal Linear Regulators with Incomplete State Feedback," IEEE Trans. on Auto. Control, vol. AC12, no. 3, June 1967, pp. 296-299. S1 L. M. Silverman and H. J. Payne, "Input-Output Structure of Linear Systems with Application to the Decoupling Problem", to appear in 1971. S2 C. R. Slivinsky and D. G. Schultz, "State Variable Feedback Decoupling and Design of Multivariable Systems", to appear.

145 S 3 F. B. Smith, Jr., "Application of Mark IV-SOC to Multivariable Control Problems, Part III", Air Force Flight Dynamics Laboratory, AFFDL-TR-68-10, Part IV, September 1968. Tl J. S. Tyler, Jr. and F. B. Tuteur, "The Use of a Quadratic Performance Index to Design Multivariable Control Systems", IEEE Trans. on Auto. Control, vol. ACll, no. 1, January 1966, pp. 84-92. W1 P. K. C. Wang, "Invariance, Uncontrollability, and Unobservability in Dynamical Systems, "IEEE Trans. on Auto. Control, vol. AC10, no. 3, July 1965, pp. 366-367. W2 W. A. Wolovich and R. S. Shirley, "A Frequency Domain Approach to Handling Qualities Design', Proc. 1970 JACC, Atlanta, Ga., pp. 297-305. W3 W. M. Wonham and A. S. Morse, "Decoupling and Pole Assignment in Linear Multivariable Systems, to appear in SIAM Journal on Control, 1971. Y1 E. E. Yore, "Optimal Decoupling Control Applied to the Lateral Direction Axes of a Sweep Wing Aircraft, "Air Force Flight Dynamics Lab., Wright-Patterson, Ohio, Tech. Rept. TR68-10, February 1968. Z1 IBM Corp., "System/360 Scientific Subroutine Package (360A-cm-03X) Version III Programmer's Manual," Form #H20-0205-3, 1968.

Appendix A Difficulties of Quadratic Tracking Approach Suppose one tries to design a closed-loop decoupled system with a control law (2.2) using the quadratic tracking approach of [All, [K1]. Because there is no guarantee that such an approach produces decoupled systems, we may attempt to modify it in two ways. First we impose a constraint on the elements of the F and G matrices such that they must satisfy the decoupling conditions given by (2.16). If we solve such a constrained tracking problem it turns out that unless certain restrictions are placed upon the weighting terms in the tracking problem cost functional, the application of a constraint to the tracking problem causes its solution, F and G, to be functions of the system's initial condition [R2]. The author has found a very complicated set of conditions which because' of their complexity are almost impossible to use. Iwai[Il1 has also developed somewhat similar conditions which suffer the same drawback. The other possibility is to ignore any constraint and note that cross coupling effects contribute to the error term in the tracking problem cost functional. If we weight these errors heavily enough, then it follows that cross. coupling effects can be made arbitrarily small. In general it turns out that the contribution of cross coupling to the error 146

147 terms is of about the same magnitude as tracking errors themselves. A heavy weight on these error terms invariably leads to large control amplitudes. Furthermore, solution of this unconstrained problem requires determination of all the elements of F and G involving many more variables than just the oik and Xi variables of (2.16).

Appendix B Modified W Matrix for Augmented Plants Given the plant (2.1) - = Ax + Bu(t) y(t) = Cx and compensator (2.8) x = Ax + Bu(t) with interconnection (2.9) u(t) = K u(t) + K2x we form the augmented system with stare vector x = [x x = Ax + B u(t) y = Cx where A = A BK2 B = BK C =EC 0]. We apply state variable feedback to the augmented plant S ={A, B, C} u(t) = Fx + Gv(t) and use an input filter given by (3,1) = AvP + Bvw(t) v(t) = Cv~ 148

149 Now combining the augmented plant states and the filter states, z =x we have z = Az + Bw(t) y =z where A+;BF BGC = O A B The error (v - y) as before is given by: v - y= Cv - Cx But the control u is now u = K1u + K2x = K1Fx + K1GC + K x 2 1 v 2 Consequently the error and control vector 5 is expressed as i = = Wz where (4.1) -C C W = - v K1F+K2 K1GCV K2 = [o K2]

Appendix C First Order Case Uniqueness Result Given a plant (2.1) with n = 1 and A = [-a], B = [1], F = [a + a], G = [-a] where a > 0 and a is free. For a stable closed-loop system a<O. The closed-loop transfer function is h(s) = -a S-a We select a stable first order filter with Av =[-], B = [1], Cv = [1], c>O; and weighting matrices Q = [1], R = rJ], r= [1], r > 0. Then the cost functional J given by equation (3.7) for this first order case is J = -ar' 2+ra2o -a 2ac(a-Oa) and its gradient J =' J-4..2.r. 2+4'3ra-2' 2ra2 +2'2 -y 4 (cy-a )2 - To find the mimina (or maxima) of J, we set its gradient equal to zero. The sum of the (two) roots of the numerator'43r a1 + a2 = — 4ar = +a Since a is positive (stable filter) both roots of the numerator cannot be negative, thus there can be only one minimum or maximum of J for which the closed-loop system will be stable. 1

151 To see that it has to be a minimum, we note that lim J(a) = -a lim J(a) = +o Because there are at most two stationary points, the one corresponding to a negative a must be a minimum. We also observe that no stable minimum may exist. For example let r = 1, a l, then the stationary points are located at a= 0, +a.

Appendix D Details of Result V If the matrix Ci V V V ~i(Ai+bifi) V V V V Ci(Ai+bifi )2 v V V V P.-l Ci(Ai+bifi) is expanded for the case i = ri+l, then we find the resulting matrix r -1 -C a 1 Pi 1- 1 +c3-Ca3 +a 2-~2 must have rank = 1 to cancel all the numerator roots. We require that all the rows be multiples of each other. Let ( a2x _) c (a' -cx ) 3(32 _ 2 1 1 1 (o3-a3) = -~2( 1-15 152

153 (apic- ar.) ar0 1 aPi - - (a (-al) Pi ri ri 1 1 then (row 2) = (o-a l)'(row 1) and (row 3) = (-) (row 2), etc. v v If rank = 1, then civi = 0 is a n + s condition for v vi yielding Vi,1 v. I vi,2 Ufree l,ri r. 1jj

Appendix E Application to Sample Data Systems Consider a linear constant coefficient sample data system which is described by: x(k + 1) = A x(k) + Bsu(k) y C x(k) where x(k) is an n-vector denoting the state at sampling instant k, u(k) is an m-vector control at index k, and y(k) is an m-vector output. As, Bs, Cs are constant matrices of appropriate size. A closed-loop sample data control law is applied u(k) = Fx(k)+Gv(k) where v(k) is an m-vector closed-loop input sequence and F and G are constant m x n and m x m matrices respectively. Then by association with (2.1) and (2.2) with obvious minor changes, all the following results of Chapter 3 apply. In this case the transfer functions will be z-transforms (function of z rather than s) and for stability the roots of det [zI-A-BF] must lie inside the unit circle. This is necessary so that an analog of the Lyapunov equation will have a solution. In place of a Gaussian white noise source (3.2) we can define a Gaussian sequence with covariance matrix E[w(j)w' (k)] = r:j k 154

155 for all integers j,k 2 0. Any reader who is interested in application to sample data systems should consult [Al1 and [D2 ].

Appendix F Canonically De oupled Form Gilbert [G23 has shown that any system S = {A, B, C} for which det D O can be transformed to a canonically decoupled (CD) form through a change in coordinates. We have assumed that S is controllable. Let x = T2x, then A = T2 (A-BD A*)T2A B = T BD'1 A -1 C = CT2 A where the CD state vector x can be partitioned into m + 1 sub-vectors of size pi,i = 1,..., m + 1 and the CD matrices A A A' A, B, C have the partitioned form as shown below. The f 0 0 A 0 f 0 P - f2 (F1) 0 0 f 0 m where the p -vectors f. i = 1,..., m have a 1-1 corres1 1 pondence with the ai vectors. This correspondence is shown in the transformation to companion form (3.38). 156

157 A 0...0 0 1 0 A 0...0 0 Ai is Pi x Pi A A.. ~.. ~ Ai is Pm+l x Pi 0 0 0...A 0 AC A0 AC.AC A 1 2 3 m m+1l b 0.0.. 0 b.0.0 bi is Pi x 1 B =|.i (F2) 1b is Pm+ x 1 b b. b rc o. o' c 0..b 1 O c2..0 0 C. is 1 x pi A2 O 0 C 0 m where pi > 0, i = 1,..., m, Pp r 0. M+1 m+l For i = 1,.., m the matrices Ai, bi, and ci have the partitioned form: Ai is (p1i -1) x (pi -r. ~ T j1 |D Ti is ri x (Pi-ri) i is r. x r.

158 ~~0 O8Bi 0 ai2 it 0 hi i='i=10. b = * a.*, c. = [100..0]. 1 ~1 ir where r. > 0, i = 1, *.., m. 1

TABLE OF NOMENCLATURE A,B,C plant matrices (2.1) x,u,y plant state,input,output vectors (2.1) F, G control law matrices (2.2) v closed loop input (2.2) H(s) transfer function matrix (2. 4) hi(s) diagonal element of H(s) (2.17) s Laplace transform variable (A,J) set defined in section 3.3.2 Da constant matrix (2.10) A, A(s) differential operator,corresponding transfer function (2.11) A1 A2 constant matrices (3.88) N Operator (2.12) N disturbance input matrix (3.31) q(s,F) closed-loop characteristic equation (2.18) ai(s),'i(s) polynomials in s (2.17) aik,Xi free constants (2.17) AV,BV,Cv filter matrices (3,1),w filter state and input vectors (3.1) E[ ] expected value (3.2) E[ ] steady state expected value (3.7) det determinant of matrix tr trace of matrix r diagonal covariance matrix (3.2) rn general covariance matrix z augmented state vector (3.3) 159

160 P(t),P covariance matrix of z,its steady state (3.5) combined error and control vector (3.6) W,Wn constant matrices (3.6) J scalar cost (3.7) D,A,Jik decoupling matrices (2.16) M sensitivity matrix Def. 3.4 i null space (3.33) v denotes companion form (3.38) A denotes CD form CD canonically decoupled L Theorem 3.2 An9,BnC, disturbance input filter matrices (3.50) q,n(t) filter state and input vectors (3.50) 6A,6B,6C perturbations to A,B,C (3.72) A A y,v steady state vectors (3.70) To,T steady state gain matrices (3.73) L(Ao,a) T1 (3.83) V(s),Y(s),Q(s) frequency domain input,output,plant(3.87) Qo,Q,R weighting matrices (3.7) SVi(W) spectral density of vi (3.10) Hk Hessian matrix (3.20) VJk gradient of J at k'th iteration (3.20) p,q vector form of P,BrB (3.21) A see (3.21) BIG T1,T2 similarity transformation matrices (3.34) 0i row vector, part of F (3.30) r(t) disturbance input vector (3.31) v column of N (3.32)

UNIVERSITY OF MICHIGAN 3 9015 03695 1989