- Email: [email protected]
A CAD Program for Suboptimal Linear Regulators
A CAD PROGRAM FOR SUBOPTIMAL LINEAR REGULATORS P.
J.
Fleming
Inter-University Institute of Engineering Control, School of Electronic Engineering Science, University College of North Wales, Bangor, Gwynedd, Wales
Abstract. A number of suboptimal linear regulator design strategies is discussed and it is shown that they may be incorporated within a general problem formulation. It is therefore possible to adopt a common numerical approach to the solution of these problems and a gradient minimization technique is selected as the solution method following a study of four alternative schemes. This method requires the solution of a pair of Liapunov matrix equations during each iterative step and it is found that these equations are best solved by application of the Bartels-Stewart algorithm. An interactive CAD program has been developed and it functions both as a design tool and as a research aid. As a design tool it offers access to a number of existing design strategie&, provides added flexibility in the choice of controller configuration and enables different strategies to be merged to form an effective overall design procedure. As a research aid the program is easily modified to include new controller cOufigurations or test new performance index specifications; it may also be extended to handle linear discrete-time systems and linear stochastic systems. Keywords. Computer-aided system design; linear systems; multivariable control systems; optimal control; numerical methods. INTRODUCTION
this dependence are suggested and incorporated in the program.
One of the main criticisms of linear optimal regulator design is that the resulting feedback controller requires access to the full set of plant state variables or to a state reconstructor. This has led to the development of alternative suboptimal design strategies for linear time-invariant systems, in which the feedback configuration is prespecified and tailored to suit the particular plant under consideration e.g. output feedback, dynamic compensation. Allied to these developments a number of algorithms has evolved for tackling sensitivity reduction and model-following problems in linear regulator design.
Since the program package is to be used in an interactive mode it is important to devise an efficient solution procedure and a gradient minimization approach is selected following a study of four alternative schemes. This method requires the solution of a pair of Liapunov matrix equations during each iterative step. A numerical comparison of the chief algorithms for solving these equations reveals that the Bartels-Stewart algorithm (Bartels and Stewart 1972) is well suited for this application. The program structure is illustrated by means of a block diagram and the steps involved in setting up the problem are clearly defined. Various features which aid the design task are highlighted. The program may be easily adapted to include new design strategies and planned extensions for handling discrete-time systems and stochastic systems are discussed.
It is possible to adopt a common approach to the solution of these problems and to this end a general problem formulation is described which incorporates the various design strategies. Apart from the obvious software advantages in unifying the methods in this way the developed program package enables different design philosophies to be merged to form an effective overall design procedure, e.g. model-following via output feedback control. The resulting controllers are initial-condition dependent and, in the absence of knowledge of initial conditions, various policies for reducing
SUBOPTIMAL LINEAR REGULATOR DESIGN OPTIONS The basic linear optimal regulator problem may be stated thus:259
260
P. J. Fleming
Problem:Given the linear time-invariant plant description A
x + B u , x (0) p -p p -p -p
Cp -p x
=
Joo {xT Q x -p
+ uT R
P -p
-p
u} dt. p-p
(2)
In this section the performance indices contain the control weighting term
(3)
It is well known that the solution to this problem is a constant gain feedback controller of the form
uT R u This expression is only valid -p p-p for state feedback and output feedback and should Controller Option 11 be selected this term is replaced by the compensator gain weighting terms given in Eq. (4). I Sensitivity reduction. D'Angelo 1967) Problem:-
K
the PI penalize
PI Options
o
u* (t) -p
~n
(1)
find the control, u (t), which minimizes the -p performance index (PI) J
second and third terms compensator gains.
x.
Given the plant description
(2), find the controller, u , which
(1),
s -p
(Hendricks and
-p
minimizes the PI. This formulation has certain drawbacks in connection with the implementation of the optimal controller and also the specific form of the PI, Eq.(3). The following suboptimal design options are therefore presented as an alternative to the basic design and form the basis of the CAD program. (For a full description of the symbols used - see Appendix 1).
J
0
where
X-
-s
A x + A x + B u + Bp dU-p /da, p -s s -p s -p ~s(O) = O.
Controller Options I Output feedback.
T T T } dt {x Q x + x Q x-s + U-p R -p -s s p U -p p -p
f
(Levine and Athans 1970)
Here A and B are functions of a scalar p p time-invariant parameter, a, and ~s denotes the differential trajectory sensitivity vector, dX /da. The PI contains a sensitiv-p ity weighting term and hence minimization of J will tend to reduce plant sensitivity. If Controller Option 11 is selected it is necessary to generate an additional dynamic equation in dX / da.
Problem:-
Given the plant description (1), (2), find the controller, u , where -p
u = KO Lp v , _p
which minimizes the PI, Eq.(3).
-c
Note here that the constant-gain feedback controller is constrained to feed back only the available plant outputs.
Il Virtual model-following. (Hirzinger 1975)
11 Dynamic compensation.
Problem:Given the plant description (1\ (2), find the controller u , which minimizes the PI -p
(Johnson and
Athans 1970) Problem:-
Given the plant description
(1),
(2), find the controller, u , where
u -p
A
X-
-c
J
=
-p
+ B x c Y..p c -c
C Y..p + D ~c' ~c(O) c c
Joo
{Yp-~m)T
Qp(Y..p-Y.. m) +
~~ Rp~p }
dt
o
where
-m
Ax m-m
Y.. m
C x
X
~m (0)
x -m
(5 ) o
0
and
(6 )
m-m
which minimizes the PI J
f
T T Q x + T (AT R A + C R Cc)Y..p Y..p -p p -p c p c c c
{x
0
T T (B R B + DT Rc Dc) x } dt. -c -c c P c c
+ x
(4)
The compensator description resembles that of a Luenberger observer, but differs in the respect that the order of the compensator is prespecified by the designer. The
The mode I, Eqs. (5), (6) may be the result of an earlier design or may simply represent some desirable response characteristics. Minimization of J will tend to ensure that the output of the plant behaves like that of the model. III Implici t model-following. (Tyler 1964) Problem:(I),
Given the plant description
(2), find the controller u , which
-p
A CAD Program for Suboptimal Linear Regulators
261
minimizes the PI
(x-p -x-m ) Qp (x-p -x-m ) +
J
where A is a symmetric matrix of Lagrangian mul tip Hers .
uT R u} dt -p p-p
o
x
where
Clearly the solution is dependent on the initial condition matrix, X ' and hence on o
A x
m -m
-m
and X The model-following term here is less precise than the corresponding term in PI Option 11 since it is expressed by plant and model derivative vectors. This approach demands that the order of the model equal that of the plant and excludes the use of Controller Option 11. The advantages of implicit model-following over virtual model-following are that the order of the overall problem is reduced and that it is possible in this case to set R =0 (see Markland 1970). p
m o
since ?Sc(O) = ?Ss (0)
T = -m x -m x o
=0
since
0
(see Appendix 2) •
If x
l.S not known the following alternative -Po policies may be adopted and these are incorporated l.n the program:T
(i)
Set X = E{x x } , when there is Po -Po-Po statistical information available.
(ii)
Set X = I to minimize ~ve J (see Po -Po Kleinman and Athans 1968).
GENERAL PROBLEM FORMULATION (iii) The following problem formulation was devised in order to e xpress each of the above options in a general form.
P
ll
, where P
ll
is the nxn
submatrix of P; this tends to ml.nl.ml.ze the "worst-case" cost, i.e.
Problem:-
Given the system A?S , ?S(O)
max x -Po
x
r
x
T
(see Fleming 1973) .
-0
Alternatives (i) - (iii) are not applicable for PI Option 11.
find the controller gain matrix, K, which minimizes the PI J
J
x xT -Po-Po
[Q + !Ci (KTRK + TKTRKT) Cl] ?S dt
NUMERICAL SOLUTION
o
Eqs. (10), (11) and (12) suggest the following solution methods:-
OO
r
J
~
T
(8)
Q?S dt.
(A)
o
(B)
Appendix 2 gives details o f how to construct A, Bl , B2 , Cl' C2 ' Q, R, T, K and ?S for the individual opti ons and for composite design strategies whi ch employ more than one option. It is well known that from Eqs. (7) and (8) we have J = trace {PX } , where Xo
o
- -T and PA+A P+Q = O.
x x
T
-0-0
(9)
AAT
+
AA
+ X = 0 o
K(i+l) = _2R- l (B Tp(i)A(i)C T BTp(i)A(i)C T) 1 1+ 2 2 . (C1 A(i)ci + TC1A(i)ciT)-1, where p(i) and A(i) are found by solving Eqs. (10) and (11) for K = K(i-l). (C)
Minimize J, Eq. (9) , with respect to K.
(D)
~inimize
(10)
Using gradient matrix operations (Athans and Schweppe 1965) and the method of Lagrangian multipliers (Mendel 1974) it can be shown that the nece ssary conditions for a solution to the problem are Eq. (10), (11)
Solution of the three nonlinear algebraic matrix equations for K, P and A. It e rate on
J, Eq. (9), with respect to K using gradient information, where
aJ
aK
and P and A are found by solving Eqs. (10) and (11) for a fixed value of K. In method (A) the nonlinear equations may be
P. J. Fleming
262
solved by function mLnLmization techniques where the minimizing variables are K, P and A. However, since method (C) involves fewer minimizing variables, (A) is rejected. Convergence of method (B) may be slow and is not guaranteed for all designs. Following tests on a number of problems method (D) was preferred to (C) because of the superior convergence properties of gradient search algorithms, especially since the gradient can be evaluated with modest extra computational effort. Experience has shown that quasiNewton methods perform effectively for these types of problem.
A quasi-Newton method in general minimizes a function F(~) of P independent variables, = Gl' k2' ... , kJT.
Solution of the Liapunov Matrix Equation Since a pair of these equations must be solved at each iteration of the search routine, a comparative study of three leading methods was undertaken, testing 24 matrices with orders ranging from 3 to 16. The results are summarized in Table 1. TABLE 1 Matrix Order, n
Algorithm
~
point.
From a starting
point supplied by the user it generates a sequence of points which converge to a local minimum. A typical iteration begins at
~(i) with the gradient g(~(i»
3 4 6 8 10 12 16
Computing Times for Solving the Liapunov Matrix Equation Method I
II
III
31
40 79 280 530 1050 1930 4260
14 26 60 120 190 580 940
77
220 625 770 2240 4180
and H(i) (~(i», a positive definite approximation to the Hessian matrix. The matrix equation
Computing times are in milliseconds and the results were computed on a DEC-10 computer.
is solved to give an estimate, s, of the Newton direction and the functi~n is minimized along this line. An improved . f the so 1 · k(i+l). estLmate 0 utLon, _ ,LS f oun d and the Hessian approximation is updated. The process is terminated when 11 g 11 is suitably small. -
All three methods have execution times proportional to n 3 . Method I is a widely used iterative technique developed by Smith (1968). Heinen (1971) proposed Method 11 which is an eigenvector transformation technique. From Table 1 we see that the fastest algorithm is Method III (Bartels and Stewart 1972). In this method A is transformed to upper Hessenberg form and thence to Schur form. In this way the problem is reduced to the solution of a succession of linear equations of order ~4.
In our problem the function to be minimized is J, the gain matrix, K, is stacked columnwise to form k and 3J/ 3K is similarly stacked to form g. This operation affords the opportunity 5f simplifying K by fixing specific elements at zero or non-zero values and redefining ~ so that it contains only the free elements of K. This design aspect will be discussed later. At each stage of the search J and 3J/3K are computed for the current estimate of K and this requires the solution of Eqs. (10) and (11) which are now Liapunov matrix equations. Each value of K, both at the starting point and at subsequent estimates, must stabilize A. The tendency of the minimizing routine will be to generate stabilizing valuet.~f K, however should a destabilizing K, (K=K L ), be encountered the method will fail and one recommended procedure is to restart the search from (K(i)+K(i-l»/2 with a reduced step length. On completion of a successful search a local minimum will have been found. In most designs this will coincide with the global minimum, however, if it is suspected that "better" solutions exist the routine should be re-started from a different starting
Method III offers a further advantage in this application. Since a Schur form of A is obtained for the solution of P, Eq. (10), this transformation stage may be bypassed for the corresponding solution of A, Eq. (11), and replaced by some minor program re-coding. On average this reduces the computing time for the solution of A by 60%. Although Method III is not quite as accurate as either Method I or 11 it has proved to be sufficiently accurate for this application. PROGRAM DETAILS The program is intended f or use ln an interactive conversational mode and was developed on the U.C.N.~. DEC-10 computer using a Tektronix 4010 VDU terminal. The program language is FORTRAN-IO. The organization of the program is illustrated in Fig.l and a brief description of the function of the main subroutines is given below. For full details of operation and the conversational mode used, see Fleming and Jones (1979). SUPER is the controlling routine which directs the flow of data and allows the user to select the appropriate function block.
A CAD Program for Suboptimal Linear Regulators
263
I
I
~ Fig. 1
Program structure
SETSYS interrogates the user on the design strategy to be adopted, handles the input of plant, sensitivity and model matrices and constructs the general system matrices (see Appendix 2). Data may be input from user terminal or file (DATABASE). Throughout the program, adequate provision is made for checking data before proceeding to the next phase. SETX~
interrogates the user about initial condition data, selects the appropriate policy and constructs Xo'
SETQR handles the input or modification of PI weighting matrices and constructs Q and R subject to the appropriate design strategy (see Appendix 2). SETK~ determines the fixed and free elements of K and constructs a matrix operator, S, to be used for packing the minimizing vector ~ from K and unpacking aJ/aK to form g. The initial value of K is input or modified and checked to ensure that A is a stable matrix.
SETMIN allows the user access to certain parameters of the minimization routine, namely the accuracy criterion, the maximum step length per iteration, the frequency of printing (see MONIT) and maximum number of iterations.
monitoring the progress of the minimization procedure and outputs either to user terminal or line printer. COST allows the user to examine, on request, the contribution of individual PI components e.g.
J u R u dt, thus enabling the designer o -p p-p
to make more effective use of the PI weighting matrices. PLOT uses the GHOST package (Calderbank and Prior 1977) to produce graphical output and, upon user request, displays plant, sensitivity, model, compensator and control responses. The solution of Eq. (7) may be stored on file (DATABASE) and specific graphical output may be retained using a "hard copy unit" (HCU). DATABASE offers the user the facility to store data permanently on named disc files for use in a later program phase or execution. Upon entry to the program the user is led through a set sequence of operations, as indicated by the dotted line in Fig. 1, to enter E04DDF for the first time. Upon return to SUPER from E04DDF the user is then at liberty to select successive operations. If SETSYS is re-entered the set sequence is initiated again. DESIGN APPLICATIONS
A~~
EXTENSIONS
E~4DDF
is a Numerical Algorithms Group (NAG) Library gradient minimization routine (see Gill, Murray and Pit field 1972). This routine interrogates FUNCT and GRAD and communicates with MONIT. A "soft fail" option is employed which returns program control to SUPER in the event of failure, thus allowing re-entry to E~4DDF. FUNCT and GRAD compute the PI and g, respectively, for a specific value of k using the packing operator matrix, S. K Is checked to ensure that A is stable by interrupting the solution of Eq. (l0) to derive the eigenvalues from the Hessenberg form of A using the QR algorithm. MONIT is a facility for
Practical applications for most of the individual design options which have been discussed may be found in the originating papers. There have also been examples which demonstrate the value of combining certain design options in order to satisfy a specific objective. For example, Michael and Farrar (1973) have designed an output feedback controller in a model-following context for the regulation of the longitudinal modes of a VTOL flight vehicle to yield a response which closely matches the optimal system response. Fleming (1978) has also shown that PI Options I and III may be successfully used in the second stage of a
P. J. Fleming
264
2-stage design process to desensitize an existing controller. The program facility which a l lows the designer to fix specific elements of the controller gain matrix may be effectively used in various ways . For example the device enab l es the designer to place constraints on the compensator structure, in order to set up , say , conventional P+I controllers , or to fix the elements of the compensator matrix, D , in order to establish certain compens~tor eigenvalues. Also , when measurement noise is present, it is good practice to fix A =0 so that noise is not transmitted directly to the system. Further, when using either Control Option I or 11 it is often found that certain gain elements are redundant or play an insignificant part in the controlling process and so, by application of this device , these gains may be set to zero, thus simplifying the controller structure. Linear Time- Invariant Discrete-time Systems There is a close correspondence between the solution for the discrete-time case and that for the continuous-time case. For the discrete - time case, the general problem becomes that of choosing K to minimize
subject to the constraining system equation
developed program is less straightforward. A number of design approaches have been discussed in the literature, however the formulation proposed by Mendel and Feather (1975)is promising in this case, since, in the absence of disturbances, it reduces to the deterministic formulation expressed 1n Eqs. (6) and (7). Numerical problems encountered by Mendel and Feather are currently being investigated.
CONCLUDING REMARKS A suboptimal linear regulator CAD program has been developed and fully tested. It is not meaningful to quote typical execution times for the numerical solution method since these are dependent on a variety of factors such as the degree of accuracy required and the choice of initial starting point. Nevertheless the number of iterations is roughly proportional to P, the number of minimizing variables and execution times per iteration vary as n 3 , where n is the order of A. Although the program adopts an "algori thmic approach" to control sys tern design its effectiveness may be enhanced by the designer in two ways: the intelligent use of PI weighting matrices in the iterative design process and the judicious choice of parameters which govern the minimization process. Finally, difficulties may arise in finding a stab i lizing K for open-loop unstable plants. Choi and Sirisena (1977) have cast this problem as a constrained optimization problem; this and other methods are being evaluated with a view to implementation in the program.
A~,
where the system matrices are defined as in Appendix 2. The gradient matrix 1S now
~~
2(B~PAAC~+B~PAAC~)
+
RK(CIAC~+TC1AC~T), (14)
where and
ATpA - P + Q = 0 AAi..T - A+ X
o
O.
(15) ( 16)
The discrete-time alternative is therefore easily incorporated in the pr ogram since Eqs. ( 15) and (16) can be transformed into Liapunov matrix equations by means of the bilinear transformation,
and Eq. (13) is equivalent to Eq . (14) when P is repla c ed by PA. Linear Time-Invariant Stochastic Systems Assimilation of stochastic systems into the
REFERENCES Athans, M. and F . C. Schweppe (1965) . Gradient matrices and matrix calculations. Tech. Note 1965-53, M.I . T. Lincoln Labs . , Lexington, Mass. Bartels, R. H. and G. W. Stewart (1972) . Solution of the matrix equation AX+XB=C. Commun. ACM, 15, 820-826. Calderbank, v. J-:-and \{ . A. J. Prior (1977). The GHO~.Ehl-sal output system. Culham Laboratory, UKAE (HMSO) . Choi, S. S. and H. R. Sirisena (1977). Computation of optimal output feedback controls for unstable linear multivariable systems. ~rans . Autom. Control, AC-22, 134-136. Fleming, P. J. (1973). Trajectory sensitivi:.!l--Eed\jction i..~ th~ optimal linear regulator. Ph.D. Thesis, Queen's University, Belfast. Fleming, P . J . (1978) . Desensitizing constant gain feedback linear regulators. IEEL..IEanS~l!..t.9~~0...!1.~rol, AC.:1) , 933936 . Fleming, P . J. and R. K. Jones (1979).
A CAD Program for Suboptimal Linear Regulators
Suboptimal linear regulator CAD program: User manual. V.C.N.W. Gill, P. E., W. Murray and R. A. Pitfield (1972). The implementation of two revised Quasi-Newton algorithms for vnconstrained optimization. Report NAC 11, National Physical Laboratory. Heinen, J. A. (1971). A technique for solving the extended Liapunov matrix equation. Froc. IEEE, ~, 295-296. Hendricks, T. C. and H. D'Angelo (1967). An optimal fixed control structure design with minimal sensitivity for a large elastic booster. Proc. 5th Ann. ~U .,.J;.9.nf..:..~I.L0.lZf.\!~ t and Sys tern Theory, 142-151. Hirzinger, G. (1975). Decoupling multivariable systems by optimal control techniques. Int. J. Control, 1l.. 157-168. Johnson, T. L. and M. Athans (1970). On the design of optimal constrained compensators for linear control systems. IEEE Trans. Autom. Control, AC-15, 658660. Kleinman, D. L. and M. Athans (1968). The design of suboptimal linear time-varying systems. IEEE Trans. Autom. Control, AC-13, 150-159. Levine, W. S. and M. Athans (1970). On the determination of the optimal constant output feedback gains for linear multivariable systems. IEEE Trans. Autom. Control, AC-15, 44-48. Markland,C.~1970). Optimal modelfollowing control system synthesis techniques. Proc. lEE, lll, 1705-1713. Mendel, J.M. (1974). A concise derivation of optimal constant limited state feedback gains, IEEE Trans. Autom. Control, AC-19, 447-448. Mendel, J:M.~J. Feather (1975). On the design of optimal time-invariant compensators for linear stochastic time-invariant systems. IEEE Trans. Autom. Control, AC-20, pp.653-656. Michael, G.J. and F.A. Farrar (1973). The use of model-following methods to generate suboptimal fe edback control, Int.J.Control, 17, 1325-1333. Smith, R.A. (1968)~Matrix equation XA+BX=C. SIAM J. Appl. Math., ~, 198-Z01. Tyler , J.S. (1964). The characteristics of model-following systems as synthesized by optimal control. IEEE Trans. Autom. Control, AC-9, 485-498. APPFNr:IX 1 List of Symbols
~m
model state vector, £xl
~m
model output vector, qxl
A plant system matrix, nxn p
B
plant control matrix, nxm
C
plant output matrix, rxn
p
p
Ks state feedback matrix, mxn Ko output feedback matrix, rnxr Q positive semidefinite weighting matrix, p nxn R
p
x
compensator state vector, sxl
x -s
t~aject0ry
-c
sensitivity vector, nxl
positive definite weighting matrix, mxm
~~}"ompen'"'o"
~~
m"',i"",
Rc positive definite weighting matrix, sxs As sensitivity matrix (= aAp/ aa ), nxn Bs sensiti v ity matrix (= aBp/ aa), nxm Am model system matrix,
~ x£
Cm model output matrix, qx£ . APPENDIX Z Construction of Matrices f or General Problem Formulation The full augmented state vector (ord e r N) is x =[ xTT x xT ( ax l aa )T xTJT . -p -c -s -c -m The orders of the matrices are:- A(NxN), Bl (NxM) , B (NxM) , Cl (RxN) , C (RxN) , K (MxR), 2 Z Q (NxN), R(MxM) and T (RxR). Firstly, either operation (a), (b) or (c) is carried out according to co ntroll er configuration. (a) Full state f e edback. N=n, M=m, R=n; A=A p ' Bl=Bp' BZ=O, Cl=I n , C2 =O, K=K s ' Q=Q p ' R=R T=I E=I C=I. n'
p'
n'
n
(b) Controller Opti on I. N=n, M=m, R=r; A=A , Bl=B , BZ=O, Cl=C , CZ=O, K=K o ' Q=Q p ' P P A P R=R T=I E=I C=C. p' n' n' p (c) Controller Option 11. N=n+s, M=m+s, R=r+s; A = C = 1
x plant state vector, nxl -p u plant control vector, rnxl -p ~p plant output vector, rxl
Z65
Q R
~p ~,J . IJ
op
u·
ETQ E
~,'
Bl
[:p
C
0
z
~J
K =
'
B = 0 Z
[~: ;~
where
The following steps are performed sequentially t o accommodate the adoption of any of the PI
P. J. Fleming
266
~
If PI Option 11 is not to be adopted, go to step 4; otherwise set
options. ~
If PI Option III is not to be adopted,
N
'T '
go to step 2; otherwise set Q ~ A QA, where A = (A +B KC -L). (Note:- PI Option III p P P is not used in conjunction with Controller Option 11.) Go to step 2. Step 2 If PI Option I is not to be adopted, ~ = 0 and go to Step 3; otherwise, set s N ~ 2N, A~
~
Cl
E
+
~A/ao. ~, ".~:,/,"] [Cl
~
step 3.
ol, On,N/2] ,
C 2 C
~
[0
• '2 •
C~
[~ ~J
[~J Q~
ETQsE
and go to
A+
Cl Q+
·
~o o J
~
+
N+£,
[~ ~J @l
[Q
l
'T
Q2
oJ QJ
'Q1 h were
C'TQC'
'
Q
Q =_ETCTQ C and 0 = CTQ C ,and go to
,
3 2 P P m mP m step 4. T ~ Set up X ' where Xo ~(O) ~(O). If o either alternative policy (i), (ii) or (iii) is employed, Xm o. o
J.
Fleming
Inter-University Institute of Engineering Control, School of Electronic Engineering Science, University College of North Wales, Bangor, Gwynedd, Wales
Abstract. A number of suboptimal linear regulator design strategies is discussed and it is shown that they may be incorporated within a general problem formulation. It is therefore possible to adopt a common numerical approach to the solution of these problems and a gradient minimization technique is selected as the solution method following a study of four alternative schemes. This method requires the solution of a pair of Liapunov matrix equations during each iterative step and it is found that these equations are best solved by application of the Bartels-Stewart algorithm. An interactive CAD program has been developed and it functions both as a design tool and as a research aid. As a design tool it offers access to a number of existing design strategie&, provides added flexibility in the choice of controller configuration and enables different strategies to be merged to form an effective overall design procedure. As a research aid the program is easily modified to include new controller cOufigurations or test new performance index specifications; it may also be extended to handle linear discrete-time systems and linear stochastic systems. Keywords. Computer-aided system design; linear systems; multivariable control systems; optimal control; numerical methods. INTRODUCTION
this dependence are suggested and incorporated in the program.
One of the main criticisms of linear optimal regulator design is that the resulting feedback controller requires access to the full set of plant state variables or to a state reconstructor. This has led to the development of alternative suboptimal design strategies for linear time-invariant systems, in which the feedback configuration is prespecified and tailored to suit the particular plant under consideration e.g. output feedback, dynamic compensation. Allied to these developments a number of algorithms has evolved for tackling sensitivity reduction and model-following problems in linear regulator design.
Since the program package is to be used in an interactive mode it is important to devise an efficient solution procedure and a gradient minimization approach is selected following a study of four alternative schemes. This method requires the solution of a pair of Liapunov matrix equations during each iterative step. A numerical comparison of the chief algorithms for solving these equations reveals that the Bartels-Stewart algorithm (Bartels and Stewart 1972) is well suited for this application. The program structure is illustrated by means of a block diagram and the steps involved in setting up the problem are clearly defined. Various features which aid the design task are highlighted. The program may be easily adapted to include new design strategies and planned extensions for handling discrete-time systems and stochastic systems are discussed.
It is possible to adopt a common approach to the solution of these problems and to this end a general problem formulation is described which incorporates the various design strategies. Apart from the obvious software advantages in unifying the methods in this way the developed program package enables different design philosophies to be merged to form an effective overall design procedure, e.g. model-following via output feedback control. The resulting controllers are initial-condition dependent and, in the absence of knowledge of initial conditions, various policies for reducing
SUBOPTIMAL LINEAR REGULATOR DESIGN OPTIONS The basic linear optimal regulator problem may be stated thus:259
260
P. J. Fleming
Problem:Given the linear time-invariant plant description A
x + B u , x (0) p -p p -p -p
Cp -p x
=
Joo {xT Q x -p
+ uT R
P -p
-p
u} dt. p-p
(2)
In this section the performance indices contain the control weighting term
(3)
It is well known that the solution to this problem is a constant gain feedback controller of the form
uT R u This expression is only valid -p p-p for state feedback and output feedback and should Controller Option 11 be selected this term is replaced by the compensator gain weighting terms given in Eq. (4). I Sensitivity reduction. D'Angelo 1967) Problem:-
K
the PI penalize
PI Options
o
u* (t) -p
~n
(1)
find the control, u (t), which minimizes the -p performance index (PI) J
second and third terms compensator gains.
x.
Given the plant description
(2), find the controller, u , which
(1),
s -p
(Hendricks and
-p
minimizes the PI. This formulation has certain drawbacks in connection with the implementation of the optimal controller and also the specific form of the PI, Eq.(3). The following suboptimal design options are therefore presented as an alternative to the basic design and form the basis of the CAD program. (For a full description of the symbols used - see Appendix 1).
J
0
where
X-
-s
A x + A x + B u + Bp dU-p /da, p -s s -p s -p ~s(O) = O.
Controller Options I Output feedback.
T T T } dt {x Q x + x Q x-s + U-p R -p -s s p U -p p -p
f
(Levine and Athans 1970)
Here A and B are functions of a scalar p p time-invariant parameter, a, and ~s denotes the differential trajectory sensitivity vector, dX /da. The PI contains a sensitiv-p ity weighting term and hence minimization of J will tend to reduce plant sensitivity. If Controller Option 11 is selected it is necessary to generate an additional dynamic equation in dX / da.
Problem:-
Given the plant description (1), (2), find the controller, u , where -p
u = KO Lp v , _p
which minimizes the PI, Eq.(3).
-c
Note here that the constant-gain feedback controller is constrained to feed back only the available plant outputs.
Il Virtual model-following. (Hirzinger 1975)
11 Dynamic compensation.
Problem:Given the plant description (1\ (2), find the controller u , which minimizes the PI -p
(Johnson and
Athans 1970) Problem:-
Given the plant description
(1),
(2), find the controller, u , where
u -p
A
X-
-c
J
=
-p
+ B x c Y..p c -c
C Y..p + D ~c' ~c(O) c c
Joo
{Yp-~m)T
Qp(Y..p-Y.. m) +
~~ Rp~p }
dt
o
where
-m
Ax m-m
Y.. m
C x
X
~m (0)
x -m
(5 ) o
0
and
(6 )
m-m
which minimizes the PI J
f
T T Q x + T (AT R A + C R Cc)Y..p Y..p -p p -p c p c c c
{x
0
T T (B R B + DT Rc Dc) x } dt. -c -c c P c c
+ x
(4)
The compensator description resembles that of a Luenberger observer, but differs in the respect that the order of the compensator is prespecified by the designer. The
The mode I, Eqs. (5), (6) may be the result of an earlier design or may simply represent some desirable response characteristics. Minimization of J will tend to ensure that the output of the plant behaves like that of the model. III Implici t model-following. (Tyler 1964) Problem:(I),
Given the plant description
(2), find the controller u , which
-p
A CAD Program for Suboptimal Linear Regulators
261
minimizes the PI
(x-p -x-m ) Qp (x-p -x-m ) +
J
where A is a symmetric matrix of Lagrangian mul tip Hers .
uT R u} dt -p p-p
o
x
where
Clearly the solution is dependent on the initial condition matrix, X ' and hence on o
A x
m -m
-m
and X The model-following term here is less precise than the corresponding term in PI Option 11 since it is expressed by plant and model derivative vectors. This approach demands that the order of the model equal that of the plant and excludes the use of Controller Option 11. The advantages of implicit model-following over virtual model-following are that the order of the overall problem is reduced and that it is possible in this case to set R =0 (see Markland 1970). p
m o
since ?Sc(O) = ?Ss (0)
T = -m x -m x o
=0
since
0
(see Appendix 2) •
If x
l.S not known the following alternative -Po policies may be adopted and these are incorporated l.n the program:T
(i)
Set X = E{x x } , when there is Po -Po-Po statistical information available.
(ii)
Set X = I to minimize ~ve J (see Po -Po Kleinman and Athans 1968).
GENERAL PROBLEM FORMULATION (iii) The following problem formulation was devised in order to e xpress each of the above options in a general form.
P
ll
, where P
ll
is the nxn
submatrix of P; this tends to ml.nl.ml.ze the "worst-case" cost, i.e.
Problem:-
Given the system A?S , ?S(O)
max x -Po
x
r
x
T
(see Fleming 1973) .
-0
Alternatives (i) - (iii) are not applicable for PI Option 11.
find the controller gain matrix, K, which minimizes the PI J
J
x xT -Po-Po
[Q + !Ci (KTRK + TKTRKT) Cl] ?S dt
NUMERICAL SOLUTION
o
Eqs. (10), (11) and (12) suggest the following solution methods:-
OO
r
J
~
T
(8)
Q?S dt.
(A)
o
(B)
Appendix 2 gives details o f how to construct A, Bl , B2 , Cl' C2 ' Q, R, T, K and ?S for the individual opti ons and for composite design strategies whi ch employ more than one option. It is well known that from Eqs. (7) and (8) we have J = trace {PX } , where Xo
o
- -T and PA+A P+Q = O.
x x
T
-0-0
(9)
AAT
+
AA
+ X = 0 o
K(i+l) = _2R- l (B Tp(i)A(i)C T BTp(i)A(i)C T) 1 1+ 2 2 . (C1 A(i)ci + TC1A(i)ciT)-1, where p(i) and A(i) are found by solving Eqs. (10) and (11) for K = K(i-l). (C)
Minimize J, Eq. (9) , with respect to K.
(D)
~inimize
(10)
Using gradient matrix operations (Athans and Schweppe 1965) and the method of Lagrangian multipliers (Mendel 1974) it can be shown that the nece ssary conditions for a solution to the problem are Eq. (10), (11)
Solution of the three nonlinear algebraic matrix equations for K, P and A. It e rate on
J, Eq. (9), with respect to K using gradient information, where
aJ
aK
and P and A are found by solving Eqs. (10) and (11) for a fixed value of K. In method (A) the nonlinear equations may be
P. J. Fleming
262
solved by function mLnLmization techniques where the minimizing variables are K, P and A. However, since method (C) involves fewer minimizing variables, (A) is rejected. Convergence of method (B) may be slow and is not guaranteed for all designs. Following tests on a number of problems method (D) was preferred to (C) because of the superior convergence properties of gradient search algorithms, especially since the gradient can be evaluated with modest extra computational effort. Experience has shown that quasiNewton methods perform effectively for these types of problem.
A quasi-Newton method in general minimizes a function F(~) of P independent variables, = Gl' k2' ... , kJT.
Solution of the Liapunov Matrix Equation Since a pair of these equations must be solved at each iteration of the search routine, a comparative study of three leading methods was undertaken, testing 24 matrices with orders ranging from 3 to 16. The results are summarized in Table 1. TABLE 1 Matrix Order, n
Algorithm
~
point.
From a starting
point supplied by the user it generates a sequence of points which converge to a local minimum. A typical iteration begins at
~(i) with the gradient g(~(i»
3 4 6 8 10 12 16
Computing Times for Solving the Liapunov Matrix Equation Method I
II
III
31
40 79 280 530 1050 1930 4260
14 26 60 120 190 580 940
77
220 625 770 2240 4180
and H(i) (~(i», a positive definite approximation to the Hessian matrix. The matrix equation
Computing times are in milliseconds and the results were computed on a DEC-10 computer.
is solved to give an estimate, s, of the Newton direction and the functi~n is minimized along this line. An improved . f the so 1 · k(i+l). estLmate 0 utLon, _ ,LS f oun d and the Hessian approximation is updated. The process is terminated when 11 g 11 is suitably small. -
All three methods have execution times proportional to n 3 . Method I is a widely used iterative technique developed by Smith (1968). Heinen (1971) proposed Method 11 which is an eigenvector transformation technique. From Table 1 we see that the fastest algorithm is Method III (Bartels and Stewart 1972). In this method A is transformed to upper Hessenberg form and thence to Schur form. In this way the problem is reduced to the solution of a succession of linear equations of order ~4.
In our problem the function to be minimized is J, the gain matrix, K, is stacked columnwise to form k and 3J/ 3K is similarly stacked to form g. This operation affords the opportunity 5f simplifying K by fixing specific elements at zero or non-zero values and redefining ~ so that it contains only the free elements of K. This design aspect will be discussed later. At each stage of the search J and 3J/3K are computed for the current estimate of K and this requires the solution of Eqs. (10) and (11) which are now Liapunov matrix equations. Each value of K, both at the starting point and at subsequent estimates, must stabilize A. The tendency of the minimizing routine will be to generate stabilizing valuet.~f K, however should a destabilizing K, (K=K L ), be encountered the method will fail and one recommended procedure is to restart the search from (K(i)+K(i-l»/2 with a reduced step length. On completion of a successful search a local minimum will have been found. In most designs this will coincide with the global minimum, however, if it is suspected that "better" solutions exist the routine should be re-started from a different starting
Method III offers a further advantage in this application. Since a Schur form of A is obtained for the solution of P, Eq. (10), this transformation stage may be bypassed for the corresponding solution of A, Eq. (11), and replaced by some minor program re-coding. On average this reduces the computing time for the solution of A by 60%. Although Method III is not quite as accurate as either Method I or 11 it has proved to be sufficiently accurate for this application. PROGRAM DETAILS The program is intended f or use ln an interactive conversational mode and was developed on the U.C.N.~. DEC-10 computer using a Tektronix 4010 VDU terminal. The program language is FORTRAN-IO. The organization of the program is illustrated in Fig.l and a brief description of the function of the main subroutines is given below. For full details of operation and the conversational mode used, see Fleming and Jones (1979). SUPER is the controlling routine which directs the flow of data and allows the user to select the appropriate function block.
A CAD Program for Suboptimal Linear Regulators
263
I
I
~ Fig. 1
Program structure
SETSYS interrogates the user on the design strategy to be adopted, handles the input of plant, sensitivity and model matrices and constructs the general system matrices (see Appendix 2). Data may be input from user terminal or file (DATABASE). Throughout the program, adequate provision is made for checking data before proceeding to the next phase. SETX~
interrogates the user about initial condition data, selects the appropriate policy and constructs Xo'
SETQR handles the input or modification of PI weighting matrices and constructs Q and R subject to the appropriate design strategy (see Appendix 2). SETK~ determines the fixed and free elements of K and constructs a matrix operator, S, to be used for packing the minimizing vector ~ from K and unpacking aJ/aK to form g. The initial value of K is input or modified and checked to ensure that A is a stable matrix.
SETMIN allows the user access to certain parameters of the minimization routine, namely the accuracy criterion, the maximum step length per iteration, the frequency of printing (see MONIT) and maximum number of iterations.
monitoring the progress of the minimization procedure and outputs either to user terminal or line printer. COST allows the user to examine, on request, the contribution of individual PI components e.g.
J u R u dt, thus enabling the designer o -p p-p
to make more effective use of the PI weighting matrices. PLOT uses the GHOST package (Calderbank and Prior 1977) to produce graphical output and, upon user request, displays plant, sensitivity, model, compensator and control responses. The solution of Eq. (7) may be stored on file (DATABASE) and specific graphical output may be retained using a "hard copy unit" (HCU). DATABASE offers the user the facility to store data permanently on named disc files for use in a later program phase or execution. Upon entry to the program the user is led through a set sequence of operations, as indicated by the dotted line in Fig. 1, to enter E04DDF for the first time. Upon return to SUPER from E04DDF the user is then at liberty to select successive operations. If SETSYS is re-entered the set sequence is initiated again. DESIGN APPLICATIONS
A~~
EXTENSIONS
E~4DDF
is a Numerical Algorithms Group (NAG) Library gradient minimization routine (see Gill, Murray and Pit field 1972). This routine interrogates FUNCT and GRAD and communicates with MONIT. A "soft fail" option is employed which returns program control to SUPER in the event of failure, thus allowing re-entry to E~4DDF. FUNCT and GRAD compute the PI and g, respectively, for a specific value of k using the packing operator matrix, S. K Is checked to ensure that A is stable by interrupting the solution of Eq. (l0) to derive the eigenvalues from the Hessenberg form of A using the QR algorithm. MONIT is a facility for
Practical applications for most of the individual design options which have been discussed may be found in the originating papers. There have also been examples which demonstrate the value of combining certain design options in order to satisfy a specific objective. For example, Michael and Farrar (1973) have designed an output feedback controller in a model-following context for the regulation of the longitudinal modes of a VTOL flight vehicle to yield a response which closely matches the optimal system response. Fleming (1978) has also shown that PI Options I and III may be successfully used in the second stage of a
P. J. Fleming
264
2-stage design process to desensitize an existing controller. The program facility which a l lows the designer to fix specific elements of the controller gain matrix may be effectively used in various ways . For example the device enab l es the designer to place constraints on the compensator structure, in order to set up , say , conventional P+I controllers , or to fix the elements of the compensator matrix, D , in order to establish certain compens~tor eigenvalues. Also , when measurement noise is present, it is good practice to fix A =0 so that noise is not transmitted directly to the system. Further, when using either Control Option I or 11 it is often found that certain gain elements are redundant or play an insignificant part in the controlling process and so, by application of this device , these gains may be set to zero, thus simplifying the controller structure. Linear Time- Invariant Discrete-time Systems There is a close correspondence between the solution for the discrete-time case and that for the continuous-time case. For the discrete - time case, the general problem becomes that of choosing K to minimize
subject to the constraining system equation
developed program is less straightforward. A number of design approaches have been discussed in the literature, however the formulation proposed by Mendel and Feather (1975)is promising in this case, since, in the absence of disturbances, it reduces to the deterministic formulation expressed 1n Eqs. (6) and (7). Numerical problems encountered by Mendel and Feather are currently being investigated.
CONCLUDING REMARKS A suboptimal linear regulator CAD program has been developed and fully tested. It is not meaningful to quote typical execution times for the numerical solution method since these are dependent on a variety of factors such as the degree of accuracy required and the choice of initial starting point. Nevertheless the number of iterations is roughly proportional to P, the number of minimizing variables and execution times per iteration vary as n 3 , where n is the order of A. Although the program adopts an "algori thmic approach" to control sys tern design its effectiveness may be enhanced by the designer in two ways: the intelligent use of PI weighting matrices in the iterative design process and the judicious choice of parameters which govern the minimization process. Finally, difficulties may arise in finding a stab i lizing K for open-loop unstable plants. Choi and Sirisena (1977) have cast this problem as a constrained optimization problem; this and other methods are being evaluated with a view to implementation in the program.
A~,
where the system matrices are defined as in Appendix 2. The gradient matrix 1S now
~~
2(B~PAAC~+B~PAAC~)
+
RK(CIAC~+TC1AC~T), (14)
where and
ATpA - P + Q = 0 AAi..T - A+ X
o
O.
(15) ( 16)
The discrete-time alternative is therefore easily incorporated in the pr ogram since Eqs. ( 15) and (16) can be transformed into Liapunov matrix equations by means of the bilinear transformation,
and Eq. (13) is equivalent to Eq . (14) when P is repla c ed by PA. Linear Time-Invariant Stochastic Systems Assimilation of stochastic systems into the
REFERENCES Athans, M. and F . C. Schweppe (1965) . Gradient matrices and matrix calculations. Tech. Note 1965-53, M.I . T. Lincoln Labs . , Lexington, Mass. Bartels, R. H. and G. W. Stewart (1972) . Solution of the matrix equation AX+XB=C. Commun. ACM, 15, 820-826. Calderbank, v. J-:-and \{ . A. J. Prior (1977). The GHO~.Ehl-sal output system. Culham Laboratory, UKAE (HMSO) . Choi, S. S. and H. R. Sirisena (1977). Computation of optimal output feedback controls for unstable linear multivariable systems. ~rans . Autom. Control, AC-22, 134-136. Fleming, P. J. (1973). Trajectory sensitivi:.!l--Eed\jction i..~ th~ optimal linear regulator. Ph.D. Thesis, Queen's University, Belfast. Fleming, P . J . (1978) . Desensitizing constant gain feedback linear regulators. IEEL..IEanS~l!..t.9~~0...!1.~rol, AC.:1) , 933936 . Fleming, P . J. and R. K. Jones (1979).
A CAD Program for Suboptimal Linear Regulators
Suboptimal linear regulator CAD program: User manual. V.C.N.W. Gill, P. E., W. Murray and R. A. Pitfield (1972). The implementation of two revised Quasi-Newton algorithms for vnconstrained optimization. Report NAC 11, National Physical Laboratory. Heinen, J. A. (1971). A technique for solving the extended Liapunov matrix equation. Froc. IEEE, ~, 295-296. Hendricks, T. C. and H. D'Angelo (1967). An optimal fixed control structure design with minimal sensitivity for a large elastic booster. Proc. 5th Ann. ~U .,.J;.9.nf..:..~I.L0.lZf.\!~ t and Sys tern Theory, 142-151. Hirzinger, G. (1975). Decoupling multivariable systems by optimal control techniques. Int. J. Control, 1l.. 157-168. Johnson, T. L. and M. Athans (1970). On the design of optimal constrained compensators for linear control systems. IEEE Trans. Autom. Control, AC-15, 658660. Kleinman, D. L. and M. Athans (1968). The design of suboptimal linear time-varying systems. IEEE Trans. Autom. Control, AC-13, 150-159. Levine, W. S. and M. Athans (1970). On the determination of the optimal constant output feedback gains for linear multivariable systems. IEEE Trans. Autom. Control, AC-15, 44-48. Markland,C.~1970). Optimal modelfollowing control system synthesis techniques. Proc. lEE, lll, 1705-1713. Mendel, J.M. (1974). A concise derivation of optimal constant limited state feedback gains, IEEE Trans. Autom. Control, AC-19, 447-448. Mendel, J:M.~J. Feather (1975). On the design of optimal time-invariant compensators for linear stochastic time-invariant systems. IEEE Trans. Autom. Control, AC-20, pp.653-656. Michael, G.J. and F.A. Farrar (1973). The use of model-following methods to generate suboptimal fe edback control, Int.J.Control, 17, 1325-1333. Smith, R.A. (1968)~Matrix equation XA+BX=C. SIAM J. Appl. Math., ~, 198-Z01. Tyler , J.S. (1964). The characteristics of model-following systems as synthesized by optimal control. IEEE Trans. Autom. Control, AC-9, 485-498. APPFNr:IX 1 List of Symbols
~m
model state vector, £xl
~m
model output vector, qxl
A plant system matrix, nxn p
B
plant control matrix, nxm
C
plant output matrix, rxn
p
p
Ks state feedback matrix, mxn Ko output feedback matrix, rnxr Q positive semidefinite weighting matrix, p nxn R
p
x
compensator state vector, sxl
x -s
t~aject0ry
-c
sensitivity vector, nxl
positive definite weighting matrix, mxm
~~}"ompen'"'o"
~~
m"',i"",
Rc positive definite weighting matrix, sxs As sensitivity matrix (= aAp/ aa ), nxn Bs sensiti v ity matrix (= aBp/ aa), nxm Am model system matrix,
~ x£
Cm model output matrix, qx£ . APPENDIX Z Construction of Matrices f or General Problem Formulation The full augmented state vector (ord e r N) is x =[ xTT x xT ( ax l aa )T xTJT . -p -c -s -c -m The orders of the matrices are:- A(NxN), Bl (NxM) , B (NxM) , Cl (RxN) , C (RxN) , K (MxR), 2 Z Q (NxN), R(MxM) and T (RxR). Firstly, either operation (a), (b) or (c) is carried out according to co ntroll er configuration. (a) Full state f e edback. N=n, M=m, R=n; A=A p ' Bl=Bp' BZ=O, Cl=I n , C2 =O, K=K s ' Q=Q p ' R=R T=I E=I C=I. n'
p'
n'
n
(b) Controller Opti on I. N=n, M=m, R=r; A=A , Bl=B , BZ=O, Cl=C , CZ=O, K=K o ' Q=Q p ' P P A P R=R T=I E=I C=C. p' n' n' p (c) Controller Option 11. N=n+s, M=m+s, R=r+s; A = C = 1
x plant state vector, nxl -p u plant control vector, rnxl -p ~p plant output vector, rxl
Z65
Q R
~p ~,J . IJ
op
u·
ETQ E
~,'
Bl
[:p
C
0
z
~J
K =
'
B = 0 Z
[~: ;~
where
The following steps are performed sequentially t o accommodate the adoption of any of the PI
P. J. Fleming
266
~
If PI Option 11 is not to be adopted, go to step 4; otherwise set
options. ~
If PI Option III is not to be adopted,
N
'T '
go to step 2; otherwise set Q ~ A QA, where A = (A +B KC -L). (Note:- PI Option III p P P is not used in conjunction with Controller Option 11.) Go to step 2. Step 2 If PI Option I is not to be adopted, ~ = 0 and go to Step 3; otherwise, set s N ~ 2N, A~
~
Cl
E
+
~A/ao. ~, ".~:,/,"] [Cl
~
step 3.
ol, On,N/2] ,
C 2 C
~
[0
• '2 •
C~
[~ ~J
[~J Q~
ETQsE
and go to
A+
Cl Q+
·
~o o J
~
+
N+£,
[~ ~J @l
[Q
l
'T
Q2
oJ QJ
'Q1 h were
C'TQC'
'
Q
Q =_ETCTQ C and 0 = CTQ C ,and go to
,
3 2 P P m mP m step 4. T ~ Set up X ' where Xo ~(O) ~(O). If o either alternative policy (i), (ii) or (iii) is employed, Xm o. o