OPTIMIZATION-BASED CONTROLLER TUNiNG USING THE Q-PARAMETERIZATIONByAlan F. LynchB.A.Sc. University of Toronto, Toronto, Ont. 1991A THESIS SUBMITTED IN PARTIA. FULFILLMENT OFTHE REQUIREMENTS FOK PHE DEGREE OFMASTER OF APPLIED SCIENCEinPHE FACULTY F C DU E S iDlESDEPARTMEN OF ELI’:i IICAL I GINEEBINCWe accept this t s as nformingto the requil standardTHE UNIVERSITY OF BRITISH COLUMBIADecember 1993© Alan F. Lynch, December 31, 1993In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the Universityof British Columbia, I agree that the Library shall make it freely available for reference and study. Ifurther agree that permission for extensive copying of this thesis for scholarly purposes may be granted bythe head of my department or by his or her representatives. It is understood that copying or publicationof this thesis for financial gain shall not be allowed without my written permission.Department of Electrical EngineeringThe University of British Columbia2075 Wesbrook PlaceVancouver, CanadaV6T 1W5Date: 3 / /AbstractWe consider the problem of tuning a closed-loop linear controller for a continuous-time linear, lumped,exponentially stable plant using measurements. A two-step tuning method is proposed based on an“internal model” controller structure. The motivation for developing the method is to reduce the effectof modeling error on the design and to reduce the amount of computation required compared to a similardesign performed off-line.The first step is to solve a controller design problem formulated as a sequence of convex optimizationproblems in which the cost and constraint functionals and their descent directions are computed directlyfrom plant measurements. The designer is comfortable in specifying desired performance by adjustingthe weights of a weighted-max cost function composed of a wide range of time and frequency domainperformance functionals.To implement the “internal model” controller, a plant model is obtained in the second step. Theidentification objective makes the plant model depend on the desired closed-loop performance. To ensurea robust design, the model satisfies a worst-case error bound which depends on information about theplant, the experiment duration, the noise level, and the order of the model.11Table of ContentsAbstract iiList of Figures viiList of Tables xNotation xiAcknowledgements xiii1 Introduction 11.1 Feedback system 11.2 Control design step 21.2.1 Background on linear controller design methods 21.2.2 Using plant measurements to perform an optimization-based design 81.3 Control-oriented identification step 91.4 Statement of problem 101.5 Contributions 111.6 Thesis outline 112 Control design step of the tuning method 132.1 Closed-loop transfer matrix 132.2 Parameterization of the controller 132.2.1 Internal model controller structure Stable plant 142.2.2 Implementation of the design parameter Q for optimization-based controller design 152.3 Controller design problem 172.3.1 Performance functionals 172.3.2 Performance specifications 172.3.3 Form of the optimization problem 181112.3.4 Tuning procedure.192.4 Performing a design on P 202.5 Optimization algorithms 232.6 Internal model controller structure Unstable plant 243 Computing the performance functionals in the control design step 273.1 The step response functional 273.1.1 Using measurements to evaluate the step response functional 293.1.2 Using measurements to evaluate a subgradient of the step response functional . . 303.1.3 Computational cost 303.1.4 A step response functional example 313.2 Norms of particular response functionals 313.2.1 Using measurements to evaluate the norm of a particular response functional andits subgradient 353.2.2 Computational cost 353.2.3 An example of using the “L norm” functional 353.3 The “peak gain norm” functional 373.3.1 Using measurements to evaluate the “peak gain norm” functional 383.3.2 Using measurements to evaluate a subgradient of the “peak gain norm” functional 393.3.3 Computational cost 403.3.4 An example of using the “peak gain norm” functional 403.4 The “H2 norm” functional 413.4.1 Using measurements to evaluate the “H2 norm” functional 423.4.2 Using measurements to evaluate the subgradient of the “H2 norm” functional . . 433.4.3 Computational cost 433.4.4 An example of using an “H2 norm” functional 453.5 The “H norm” functional 453.5.1 Using measurements to evaluate the “H, norm” functional 473.5.2 Using measurements to evaluate a subgradient of the “H,,, norm” functional . . 473.5.3 Computational cost 483.5.4 An example of using the “Hc,, norm” functional 48iv4 Identification step of the tuning method4.1 Identification method4.2 Deriving the bound for the unweighted model error4.2.1 Transformation to the z-plane via the bilinear transformation4.2.2 Assumed knowledge about the plant4.2.3 Approximation error of a linear spline interpolant4.2.4 Error bound for truncating a Taylor series4.2.5 An expression for the model error bound50505152525354545 Simulation examples5.1 Example 1: Stable plant5.1.1 An analytic solution5.1.2 Open-loop tuning and another interpretation of the design problem5.1.3 Results of the identification step5.2 Example 2 Unstable plant, the ACC’ benchmark problem5.2.1 Statement of the design problemPerformance specificationsResults of the control design stepA Some rules for computing subgradientsA.1 Basic definitionsA.2 Rules for computing subgradientsB A two-step design for strongly stabiizable plants 91‘ACC stands for the American Control Conference.5.2.25.2.35656575964717172766 Conclusions/Future work6.1 Conclusions6.1.1 Control design step6.1.2 Identification step6.2 Future work8585858687898989VC Cost of evaluating the performance functionals 930.1 The step response functional 930.1.1 Evaluating the functional using a plant model 930.1.2 Evaluating the functional using measurements 970.2 Norms of a particular response functional 970.2.1 Evaluating the functional using a plant model 970.2.2 Evaluating the functional using measurements 100C.3 The “peak gain norm” functional 1000.3.1 Evaluating the functional using a plant model 1000.3.2 Evaluating the functional using measurements 1020.4 The “H2 norm” functional 1030.4.1 Evaluating the functional using the plant model 1030.4.2 Evaluating the functional using measurements 1050.5 The “H, norm” functional 1090.5.1 Evaluating the functional using the plant model 1090.5.2 Evaluating the functional using measurements 111D Cost of computing a time domain response in MATLAB 114Fl Loop Transfer Recovery (LTR) design method 115E.1 Linear Quadratic Gaussian (LQG) design method 115E.2 LTR design procedure 118F Obtaining the frequency response data for the identification step 119Bibliography 121viList of Figures1.1 Closed-loop system 21.2 Open-loop system 81.3 Tradeoff between uncertainty and achievable performance 92.4 Controller structure for a stable plant 142.5 An equivalent open-loop system 142.6 Parameterization of the design parameter Q 152.7 Iterations of the ellipsoid algorithm 252.8 Controller structure for a strongly stabilizable plant 263.9 Open-loop system 273.10 The step response functional 283.11 Computational cost : The step response functional 313.12 Computational cost The step response functional 323.13 A step response for the step response functional example 333.14 Computational cost The norm of a particular response functional 363.15 A step response for the norm of a particular response functional example 373.16 Actuator effort for the norm of a particular response functional example 383.17 Computational cost : The “peak gain norm” functional 403.18 A step response for the “peak gain norm” functional example 413.19 Actuator effort for the “peak gain norm” functional example 423.20 Computational cost The “H2 norm” functional 443.21 Computational cost The “H2 norm” functional 453.22 A response due to white-noise excitation for the “H2 norm” functional example 463.23 Computational cost The “H® norm” functional 483.24 A Bode magnitude plot for the “H norm” functional example 495.25 Example 1: An augmented plant used in the LTR design method 58vii5.26 Example 1: A weighting transfer function used in the LTR design method 595.27 Example 1: Performing the analytic design method 605.28 Example 1: Another interpretation of the analytic design problem 615.29 Example 1 : Optimal cost function values versus the number of controller parameters. . 635.30 Example 1: Number of functional evaluations to obtain the solution 645.31 Example 1: A tradeoff curve 655.32 Example 1: Assumed information about the plant 665.33 Example 1 : Worst-case model error bounds 675.34 Example 1: A step response of the optimal controller 695.35 Example 1: Bode plots of the open-loop gain 705.36 Example 2: A benchmark problem 715.37 Example 2: Describing uncertainty in the plant 735.38 Example 2 : Measuring an impulse response 755.39 Example 2: A gain margin specification 765.40 Example 2 : Bode magnitude plots of the optimal solution 795.41 Example 2 : Bode magnitude plots of the optimal solution 805.42 Example 2 : Bode plots the optimal open-loop gain 815.43 Example 2 : Time domain responses of the optimal solution 825.44 Example 2: A tradeoff curve 835.45 Example 2: Convergence of the optimization algorithm 84A.46 An example of a convex nondifferentiable functional 90C.47 Cost of evaluating the step response functional 95C.48 Cost of evaluating the norm of a particular response functional 98C.49 Cost of evaluating the “peak gain norm” functional 101C.50 Cost of evaluating the “H2 norm” functional 104C.51 Relative cost of evaluating the “H2 norm” functional 109C.52 Cost of evaluating the “H, norm” functional 112C.53 Relative cost of evaluating the “H norm” functional 113E.54 LQG controller structure 116viiiE.55 An augmented plant used in the LTR design method 117F.56 Identification experiment setup 119ixList of Tables5.1 Example 1: An analytic solution.615.2 Example 1: An optimal controller 685.3 Example 2 Optimal values of the objectives 775.4 Example 2: Optimal controller parameters 78xNotationNotation Meaningarg mm A minimizer of the argument.arg inf An infimizer of the argument.C The complex numbers.Co+ {s e CIRs> O}, the open right-half plane.B,, {s e Clisi < p}.B B1 = {s e Cilsi < 1}, the open unit disk.FIX The expected value of a random variable X.Hba The closed-loop transfer function from signal a to signal b.The imaginary part of the complex number s.inff The infimum of a function..1L1 One-sided inverse Laplace transform operator.max The maximum of a function or set.mm The minimum of a function or set.N The nonnegative integers.R, R The Real numbers and nonnegative Real numbers, respectively.RN The vector space of n-component real vectors.The space of functions {f: C0+—+ Cif is analytic in C0÷ and sup5€C+ If(s)I <oo}.The space of functions {f B —* Cif is analytic in B and supEB If(z)L <oo}.RH,,, The space of proper real-rational transfer functions.RH2 The space of strictly proper real-rational transfer functions.Js The real part of a complex number s.sgn(.) The signum function.sup The supremum of a function.The complex conjugate of the complex number s.AH The complex conjugate transpose of the matrix A with complex entries.xiVf(x) The gradient of the function f with respect to x at a point x.(V, llll) The vector space V endowed with a normlull2 The L2 norm of a function jj’° uQ)2dt.lIHll2 The H2 norm of a transfer function. If H C RH2 then1/2llHll2= (acr IH(iw)l2dw)I Iull The L norm of a sequence supkEN lu(k)l or function supR+ luQ)lll-”ll sup3€c÷ IH(s)I or sup2D IH(z)l, depending on whether H is a belongs of Hor H,D.supjj IH(z)I.xiiAcknowledgementsI would like thank my supervisor Dr. S.E. Salcudean for his patience and support. Many of the ideaspresented in the thesis are a result of our discussions.I would like to thank NSERC and ASI for providing me with the financial support to obtain myM.A.Sc.As well, I thank my friends and the students I have met at UBC while pursuing my degree. Inparticular, I am indebted to Peter Racksel, Brian Perkins, John Ru, and Joseph Yan for their help.Finally I thank Christine for her encouragement and for being able to put up with me.X1HChapter 1IntroductionA standard approach to designing a feedback controller applies a design method to a plant model derivedfrom laws of physics or experimental data. A particularly powerful design method solves a nonlinearprogramming problem where the cost and constraint functionals are determined by performance specifications (see for example [PMS84, GD83]). However, this design method is typically performed on a loworder model of the plant to reduce the complexity of the design. The effects of ignoring the higher orderdynamics of the plant on the performance of the system are uncertain. Hence, this design method hastwo drawbacks: computationally intensive optimization problems for high order plant models and lowperformance due to design on a simplified model. The tuning1 method proposed in this thesis exploitsthe flexibility of an optimization-based design method while addressing its shortcomings by performinga design using plant to a model) and then identifying the plant accountingfor desired closed-loop performance.In this chapter we first present the feedback system considered throughout the thesis followed bybackground on linear controller design methods stressing the benefits of a parameter optimization-basedmethod. Next, motivation for designing using measurements from the plant as opposed to a model isdiscussed and an outline of an optimization-based tuning method consisting of a design and identificationstage is presented. Finally, contributions of the thesis are given and its organization outlined.1.1 Feedback systemThe unity-feedback system considere shown in Figure 1.1, where P is a plant, K a controller, r areference signal, d a plant input disturbanc” y. a system output, iii a measurement noise, e an errorsignal, u a controller output, and w a plant state disturbance.We restrict our attention to plants and controllers which are linear, time-invariant (LTI), lumped,single-input, single-output (SISO), and operate in continuous-time. The time domain input-output‘Tuning means designing the controller using measurements from the system and not the plant’s transfer function as isoften the case in off-line design.1Chapter 1. Introduction 2+ + mFigure 1.1: Closed-loop system.properties of such systems are described by the convolution relationy = g * it,where the system’s impulse response is denoted by g and its (scalar valued) input and output are denotedby it and y respectively. Since the systems are lumped and LTI they have real-rational transfer functions.Throughout the thesis the upper case letter denotes both the system and its transfer function. The “real”plant with unknown transfer function F(s), is assumed strictly proper and exponentially stable. Thelatter restriction means there exists a a < 0 such thatwhere p denotes the impulse response of the plant. Finally, the notation Hba is used to denote thesystem’s closed-loop transfer function from signal a t ;signal b.1.2 Control design step1.2.1 Background on linear controller design methodsThe control design step of the tuning method is based on parameter optimization-based design methods.It is only recently, with development in optimization techniques and algebraic control theory, that theChapter 1. Introduction 3flexibility of such design methods can be fully exploited. This flexibility is understood by comparingthe methods with synthetic and analytic design procedures. The division of methods into the classes ofsynthetic, analytic, and parameter optimization-based is due to [BBN9O].Synthetic methodsSynthetic methods, discussed in [NGK57, 11or63], are the oldest design methods in which the controlleris constructed in sections which are added until design specifications are met. Although easy to use andunderstand, many of the design steps of these methods involve trial and error or applying rules of thumb.For example, consider the synthetic loop shaping technique discussed in [Hor63] which constructs theopen-loop gain L = KP to achieve performance specifications and solves for K using K = L/P. If thedesign problem is to obtain a stabilizing controller which asymptotically tracks reference signals withenergy uniformly concentrated over [0, 1] rad/s, and if the plant is minimum phase and has a relativedegree of 1, a suitable choice for L is the low-pass filterL(s) = j—, (1.1)which attenuates frequencies outside the operating bandwidth of 1 rad/s [DFT9O]. Tracking performanceis improved by multiplying (1.1) by a constant gain chosen small enough to ensure stability. Thedesign must be repeated with some other choice of transfer functions to construct L if unsatisfactoryperformance is obtained. Hence, in cases where performance specifications are difficult to achieve, therough guidelines available to obtain L are inadequate and the method fails. Although synthetic methodsprovide a simple way for designing a controller with understandable structure, each section havingan identifiable purpose, their application is awkward and relies on the experience of the designer andundemanding design specifications.Analytic methodsAnalytic methods, the second type of design method, produce controllers that minimize some specificcost function. For example, the LQG method described in [BH75] minimizes the weighted-sum of the“steady-state” mean-squared value of the state and output in the presence of white Gaussian processand measurement noise. The designer chooses spectral factors of the noise and weights on the outputand state so that the controller produced meets the performance specifications. Although the LQGmethod guarantees a controller which has desirable stability margin and minimizes the cost function (aChapter 1. Introduction 4measure of the system’s sensitivity to noise), it is rare that design specifications are directly related tothe physical interpretation of the cost function (the mean-squared value of the weighted state and outputdue to white Gaussian noise inputs). Consequently, desirable performance is obtained by applying theLQG method in an iterative manner with the weighting matrices and noise spectral factors acting astuning parameters. Loop Transfer Recovery (LTR) is an example of one iterative LQG design methodwhich gives the designer a systematic procedure for loop shaping by solving a sequence of LQG problems[DS81, Mac89],The problem of weight selection is also common to fl-optimal control theory, introduced by Zames[Zam8l], which considers the problem of minimizing the H, norm of the closed-loop transfer matrix.The recent popularity of the H,-optimal method is due to its treatment of robustness specificationswhich limit the maximum variation in the performance of the closed-loop system due to a set of plantperturbations.To illustrate the difficulty of weight selection for H-optimal control methods, consider a specialcase of the standard problem which minimizesmax{IW1(jw)S(jw)12+1W2(jw)T(jw)1}, (1.2)wcRwhere l41 and l42 are stable weighting trans’ f’.v:ions chosen by the designer and T = Hyr, S = Herare the complementary sensitivity and sensitivity functions respectively. This problem is known asthe mixed sensitivity problem and was first solved for the multi-input, multi-output (MIMO) case in{Kwa85, VJ84]. A controller is said to provide robust performance if it has a cost function (1.2) lessthan one. In this case the nominal performance specification II14”i5II < 1 holds for all plants in themultiplicative perturbation setf P(iw)-PUw) <1W2(jw)I,VwP(jw) —where P is a stable nominal plant rndel. ‘T.e.staudaid; approach is to chose IWi(iw)I (1W2iw)I) largeat low (high) frequencies to force IS(jw;I (ITuwI) t:. ‘;esniall there. The difficulty lies in incorporatinga wide range of engineering specifications such as corv’aints on. the magnitudes of closed-loop frequencyresponses and on time domain responses into the weights Wi and W2. In fact, there is no general way ofdetermining the weights directly from the performance specifications and, as with the LQG procedure,the weighting transfer functions end up being tuned until desirable performance is obtained. For example,[DFT9O] considers the problem of weight selection for time domain specifications on a flexible beam.Chapter 1. Introduction 5The authors show how the framework of the mixed sensitivity problem cannot deal with constraintson the step response overshoot and peak value of the plant input precisely and an iterative procedurefor tuning the transfer function design weights is required. Although tuning procedures are sometimessuccessful in practice, they may not be adequate if the optimization problem is sensitive to the choice ofweights [PS89]. Moreover, the design weights should not be tuned on-line since it is sometimes difficultto predict their effects on performance.Hence, analytic method are more methodical than synthetic methods because they produce controllerswhich guarantee certain desirable feedback properties such as stability. On the other hand, they aredifficult to apply when performance specifications are not related in an obvious way to the choice ofdesign weights. Parameter optimization-based design methods and their advantages over analytic andsynthetic methods are discussed next.Parameter optimization methodsMuch work has been done on parameter optimization-based design, see for example [ZA73, Edm79,DF81, 0D83, PMS84, 0D85, NT86, PS89, BB91]. The method involves selecting a controller structureand performance index which is minimized over the controller parameters. Since the designer is freeto select the cost function, more complex and realistic design specifications can be dealt with. Forexample, the cost function can readily combine functionals measuring system performance in both thetime and frequency domain. In addition to the benefits of a general cost function, recent developmentsin numerical optimization techniques and the emergence of the Q-parameterization theory have madeparameter optimization-based methods particularly attractive.Recent development in optimization algorithms and design methodsOnly recently have optimization algorithms been developed which accommodate the nondifferentiablefunctionals specifying performance of a control system. For example,= sup y(z,t)I (1.3)measures the peak value of the output response y and is not differentiable with respect to the controllerparameter x due to the supremum operation. Also, in many cases performance is specified by semiinfinite inequalities such as= sup IHer(x,iW)I — K(w) 0, (1.4)eRChapter 1. Introduction 6which is equivalent to an infinite system of inequalities (indexed by w) in the finite-dimensional controllerparameter. The optimization problem involving functionals such as (1.3) and (1.4) is called semi-infiniteand because the functionals are not differentiable, a nondifferentiable optimization algorithm must beused {PMS84, Sal86, Po187].Although flexibility in the choice of cost and Donstraint functionals allows an optimization-baseddesign to deal directly with a wide range of specifications, the inherent difficulty in stating a realisticcontrol problem requires a interactive procedure for its solution {NT86]. For instance, in many applications objectives compete with each other and the designer must decide on a satisfactory tradeoff.This compromise is obtained by performing a sequence of optimization runs with cost functions emphasizing different objectives. Additionally, specifications cannot always be quantified precisely: becausethere is insufficient information about the system or because the specifications are flexible. An interactive procedure allows the designer to explore different tradeoffs between the objectives before decidingwhat constitutes an optimal design. For example, consider the problem of meeting specifications on therise-time of the output step response and the peak value of the actuator output u. If the designer isflexible with the rise-time specification and wants to know how its slight violation can reduce the peakvalue of u, then the solution is obtained by exploring the tradeoff curve between the objectives untilthe best compromise is obtained, It is important to note that although computer-designer interactionis a key element to successfully applying an optimization-based design method, it is only recently thatmethodologies have been developed, e.g. the CAD system DELIGHT described in [NT86j.Formulation of a convex optimization problemNumerical algorithms applied to nonconvex2 optimization problems can yield only local minima, anddetermining whether these minima are global is generally computationally intractable when there aremore than a few design parameters. Only when the problem is convex can we guarantee that a globalsolution to the problem can be found. We believe this confidence in the solution to be an importantattribute of the tuning method presented in this thesis. Although the idea of formulating a controllerdesign problem as a convex optimization problem is old, only recently with the development of theso-called Q-parameterization theory has the importance of convexity been stressed [Sal86, PS89, BB91j.The original form of the Q-parameterization theory is due independently to [YBJ76, Kuc79] andparameterizes the set of internally stable closed-loop transfer matrices using a factorization in terms of2An optimization problem is convex if the cost and constraint functions are convex.Chapter 1. Introduction 7polynomial matrices. The present form of the theory used in this thesis is due to [DLMS8O, Vid85] andparameterizes the set of internally stable closed-loop transfer matrices using a factorization in termsof stable transfer matrices. This parameterization has been widely used in controller design methods:[GD83, 0D85] describes an optimization-based design method, [MV87] presents a globally convergentadaptive control scheme, [NJM88] describes a controller design method detecting failures in the plant, and[MZ89] gives a robust control design method for the process control industry. The Q-parameterizationtheory has also been used in identification to parameterize all plants stabilized by a known controller:[11L87] uses it to parameterize the input-output dynamics of the plant, and [Han89] uses it to simplifythe analysis of closed-loop identification,Most commonly used controller structures (e.g. proportional-integral-derivative (PID)) do not necessarily lead to performance functionals that are convex in the controller parameter. Consider a controllerthat is convex in its parameters (a PID controller for example), the closed-loop transfer matrix is usuallya nonconvex function of the controller parameter due to its linear fractional dependence on K. Hence,functions important for design such as the magnitude of the frequency response of a closed-loop transferfunction at a specific frequency are not necessarily convex in the controller parameter [BBB+ 88]. On theother hand, by introducing the Q-parameterization ‘e transform the closed-loop transfer matrix into an‘ affine function of the free parameter Q which ranges over the set of stable rational transfer functions. For- example, Hy,. = PQ is a linear function of Q. Consequently, most important performance functionalsassociated with control design problems, including Hc,o, H2, and L1 norm functionals, are convex in Q.This observation is made explicitly in [PS 89].As mentioned above, both the ideas of performing the design on Q as opposed to K and formulatinga control design problem as convex programming problem are old. The idea that H5,. is a linear functionof Q is due to [NGK57], who suggest a design is simplified in terms of Q with the controller obtainedusing K Convex programming was first used to minimize the difference between a referenceinput and the system output of a discrete-time system in [Feg64].To summarize, new optimization techniques allow a wide range of important performance functionalsto minimized in an interactive design environment which enables the designer to solve realistic controlproblems. Furthermore, the development of the Q-parameterization theory has re-emphasized the importance of solving a convex problem and making the designer confident that a (approximately Pareto3)3Pareto optimal solutions of an miiLtiobjective optimization problem are the desirable solutions that lie on the boundarybetween unachievable and achievable performance.Chapter 1. Introduction 8optimal design can be found.1.2.2 Using plant measurements to perform an optimization-based designWhile parameter optimization is emerging as a flexible design method it has two drawbacks: low performance due to modeling error in the plant and computationally intensive optimization problems.The traditional approach to designing a controller is based on a plant model. The design is performedrobustly due to error in the nominal model and resulting performance is typically conservative. Thisthesis suggests that the control design step he performed directly from measurements on the plant toreduce the effect of modeling errors on performance.Further justification for using measurements as opj.osed to a model to directly perform the designis given in {GD85] in the form of an example which is repeated here. Consider the Routh-Hurwitzmethod and the Nyquist stability criterion for determining the stability of the closed-loop. Whereasthe Nyquist stability criterion can be checked directly using the measured frequency response of theplant, the Routh-Hurwitz method requires the characteristic polynomial of the closed-loop system. It isusually difficult to obtain the characteristic polynomial and the effect of the approximations made in itscalculation are unpredictable.In addition, measurements not only offer a more direct way of performing a design, they also reducethe cost of performing optimization for many important cost functions. This is particularly importantwhen the plant model is high order or when many controller parameters are used.Hence, this thesis exploits the benefits of numerical optimization to solve a design problem on-line inwhich performance functionals of the “Youla parameter” Q and their descent directions are computedfrom measurements of an open-loop system shown in Figure 1.2. The motivation for this approachis to reduce the amount of computation required to perform the design (compared to an identicaloptimization-based design performed on a plant model) and reduce the effect of modeling errors on thesolution._H__Figure 1.2: The open-loop system used to compute the performance functionals from measurements.Chapter 1. Introduction 9The optimal solution to this on-line design problem is denoted Q°,,t. Due to the controller structureused a transfer function expression of the plant is required to obtain the expression for K; the secondstep of the tuning method is therefore to perform identification.1.3 Control-oriented identification stepThe relationship between uncertainty in the plant and achievable performance of the control systemis shown in Figure 1.3 [HJN91b]. As the level of uncertainty increases, the achievable performancedecreases. This implies that to guarantee the level of performance designed for in the control designstep, uncertainty in the plant must be reduced to a certain level by performing identification. Evidently,this procedure is only possible if the measure of plant model error is compatible with the control designstep.Limit ofperformancesystem IsunstableFigure 1.3: Tradeoff between uncertainty in the plant and achievable performance. After performing thecontrol design step, the uncertainty in the plant must be reduced by identification to a level determinedby the desired closed-loop performance.Many of the standard approaches to system identification are discussed in {Lju87j. These methodsusually model uncertainty in the form of stochastic additive noise. As well, most methods ignore uncertainty due to unmodeled dynamics of the plant [SD92]. However, for the control design step to guaranteerobust performance the identification method presented here must not only account for additive noisebut also provide H norm bound on the unmodeled dynamics.a4?iLargest acceptable level ofuncertainty required to achieve a desiredlevel of performanceUncertainty system isunstableChapter 1. Introduction 10The upper bound on the model error derived for the identification method presented in this thesisis motivated by the series of papers [HJN92, HJNOO, HJN91a, 11JN89, IIJN91c] which discuss “control-oriented identification in H,”. (This series of papers also refers to other similar relevant work on“identification in H”.) Identification methods are “control-oriented” if they provide bounds on theremaining uncertainty which are compatible with robust control design methods. “Identification inmeans the measure of the model error is the H norm, making the method compatible withcontroller design methods discussed in [Zam8l, Fra87j and the control design step discussed previously.In addition, the above papers do not assume any structure for the additive noise except that it is boundedin magnitude.Typically, the behaviour of a plant must be known accurately only over some finite interval offrequencies. Hence, it is logical to weight the model error strongly over the bandwidth of the closed-loopsystem. This is accomplished by mapping a set of frequency response data of the plant into a model Pwhich minimizes the error—By using the weight the model error is small at frequencies where the Nyquist plot of PK nearsthe —1 point since Q = 1+PK This is important because the performance of the system depends highlyon the accuracy of the model near the critical —1 point. The idea of weighting the model error with Q(which is obtained in a different manner than that presented here) appears in [BYM92] which describesan iterative control-identification design method. Note that as with the design for Q, the optimizationproblem is convex in the design parameter P (provided P is stable).1.4 Statement of problemLoosely speaking, the tuning method is performed in two steps:1. The control design step : Ideally we assume= argmin{(P,Q) I Q eexists, where the performance functional and its descent direction are computed from measurements of the open-loop system PQ.2. The identification step: Ideally we assumePop: = argmin{Q0(F— )II I P € RH00}Chapter 1. Introduction 11exist, where frequency response data of the plant is used to obtain the solution.The closed-loop controller is implemented as K(F0,Q°) = Q0(1 —1.5 ContributionsThis thesis presents a new two-step parameter optimization-based tuning method, roughly outlined inSection 1.4, with the following properties:• Since the optimization problem in the design parameter Q is convex, assuming the estimatedperformance functionals are exact, an arbitrarily good approximation to the optimal (with respectto designer specified criteria) solution is obtained.• Computation of performance functionals and their descent directions from the open-loop systemis computationally efficient because measurements are used.• The control design step for Q is performed directly from measurements of the plant and avoids theeffect of modeling errors on the design.• It solves realistic control problems: since a wide rauge of performance functionals can be incorporated into the cost and constraint functiorials inJ because an interactive design procedure allowsthe designer to understand how the objectives ampete.• The identification step minimizes a frequency domain cost function which makes the plant modeldependent on the desired closed-loop performance.• The stability of the system is guaranteed because the control design step is performed on the stableopen-loop system in Figure 1.2 and because the identification step provides a model which satisfiesa worst-case error bound.iS Thesis outlineChapter 1.2.1 describes the control design step of the tuning method by first describing how the controlleris parameterized and specifying what is meant by the “controller design problem”. An interactiveoptimization-based procedure for solving the controller design problem is presented. Finally, someremarks about the optimization algorithms used are given, and the extension of the method to stronglystabilizable plants is discussed.Chapter 1. Introduction 12As opposed to their vague description in Chapter 2, Chapter 3 precisely defines the performance functionals and their subgradients used in the tuning method and how they are computed from measurementsof the open-loop system in Figure 1.2. Furthermore, it compares the cost of computing functionals usingmeasurements with the cost of computing them from a plant model.Chapter 4 describes the identification step for the plant model required for closed-loop operation andpresents an expression for the bound on the unmodeled part of the plant.Chapter 5 presents two simulation examples using MATLAB4which illustrate the properties of thetuning method.Chapter 6 presents the conclusions and discusses possible future work that could be performed.Appendix A presents some elementary results from convex analysis. Appendix B gives details on atwo-step design procedure to extend the “internal model” controller structure to strongly stabilizableunstable plants. Appendix C presents the details of how the computational costs of the performancefunctionals in Chapter 3 are obtained. Appendix D derives the cost of computing a time domain responseusing MATLAB. Appendix E explains the Loop Transfer Recovery design method used in Chapter 3, andAppendix F explains how the frequency response data assumed available for the identification methodis obtained.4MATLAB stands for matrix laboratory, interactive software for scientific computing.Chapter 2Control design step of the tuning methodThis chapter describes the control design step of the two-step tuning method outlined in Section 1.4.The closed-loop transfer matrix of the feedback system described first, followed by a description ofthe “internal model” controller structure, the controller design problem, and the optimization-basedprocedure for its solution. Finally, there are comments on the benefits derived from using the plant toperform the design, the optimization algorithms used, and how the controller structure is applied tounstable plants.2.1 Closed-loop transfer matrixThe system’s closed-loop transfer matrix contains all ci :sd-1oop transfer functions that are of interestto the designer and is defined as[H(P, K) = (1 + PK’ rip —K K , (2.5)-P -1 1y dj =H(P,K) rne rThe vector-valued signals [y, u, ejT and [d, rn, nT are the system’s exogenous inputs and regulated outputsrespectively. The tuning method allows the designer to constrain any of the elements of (2.5).2.2 Parameterization of the controllerThis section presents the “internal model” structure used for the controller, and how the design parametercontained in the structure is itself parameterized.13Chapter 2. Control design step of the tuning method 14_____- 1 JFigure 2.5: Equivalent open-loop system from r to y when P = F, m = 0, d = 0, and w = 0 in Figure 2.4.2.2.1 Internal model controller structure: Stable plantThe structure of the controller, which is referred to as Internal Model Control (IMC) in process control[MZ89J, is shown in Figure 2.4. The expression for AK(l, Q) = Q(1 — PQF’, Q E RH, (2.6)where P is a stable model of the plant and Q is a free design parameter. The feature of the IMCstructure is that when P is exact, there is no measuiiient noise (m = 0), no plant disturbances (d = 0and w = 0), the feedback signal fi is equal to f2 and the system is effectively open-loop. Under theseconditions the equivalent system from r to y is shown in Figure 2.5 where it is obvious that any stabletransfer function Hg,. can be obtained for some choice of Q.As well, substituting (2.6) into (2.5) with P = F, we can express the closed-loop transfer matrix as+Figure 2.4: The IMC controller K(P, Q) = Q(1 — PQy’ for a stable plant.Chapter 2. Control design step of the tuning method 15I TL IFigure 2.6: The design parameter Q in the IMC controller is parameterized linearly in the controllerparameters xk. This corresponds to the Nth Ritz approximation of Q.the affine function of Q“(1- PQ) —PQ FQH(P,K(P,Q)-Q Q . (2.7)F(1-PQ) PQ-1 1-PQIt can be shown, that a stable Q is a necessary and sufficient condition for the system to be (internally)stable [MZ89].2.2.2 Implementation of the design parameter Q for optimization-based controller designSince the design parameter Q belongs to an infinit’-dnnensional space of functions, and it is evidentthat we can not numerically optimize over all RHc,o functions. We must approximate RH by somefinite-dimensional linear subspace spanned by functiops. Q’ of the formQN(x, s) := XkBk(S), (2.8)for some Bk E RH and controller parameters xk C B., k = 1,. . ., N. The parameterization of Q isshown in Figure 2.6 and is known as the Nth Ritz approximation [BB91].+Chapter 2. Control design step of the tuning method 16Examples of basis functions areB(s) = ()k>oklN (2.9)B(s) = (2.10)B3— ,_(s+IpI)((s_p)(s_p*))k_l2k_i(S)— V P ((s+p)(s+p*))kBk(s) -(s- IpI)((sp)(sp*)), ?p> 0,k = 1,. . ., N/2. (2.11)The functions (2.9) are used in [PS89] and (2.10) are the Laplace transforms of the Laguerre functions.A more general form of (2.11) are used by Kautz in [Kau54] and allow QN to contain complex poles.It is important to note that the sets {QN(x, s)Ix e RN, N N} with basis functions (2.9) aredense in (RH, II.II). This implies that any element of RH can be approximated arbitrarily well,at all frequencies, by sums of the form (2.8) with N large enough. Intuitively this result result can beunderstood by considering the bilinear transformationf(z)= p (•4-) (2.12)which maps the closed left half-plane into tii. clos ‘iit disk. The composition of (2.12) with QN,QN(xr(z)) =yields a polynomial in z of order N. Weierstrass’ Approximation Theorem says that a continuousfunction on a compact set (the unit circle) can be approximated uniformly by a polynomial function[Str84]. Since (2.12) preserves the H norm QN can approximate any element of RH by choosingN large enough. When Q is known to be strictly proper then (2.10) and (2.11) are more appropriatechoices of basis functions than (2.9). It can be shown that the sets {QN(z, s)Ix e RN, N e N} withbasis functions (2.10) are dense in (RH2,11.112) using the reasoning given above. Also, [Wah9l] containsan explanation of this result for the dis. rete-time equivalent of (2.11).It is important to note that the parameterization (2.8) typically requires a large number of controllerparameters. This is because QN has a fixed pole at s = —p and only its zeros can determine itsshape. Little work is known to have been performed on methods for the choice of basis functionsand it is a possible area of future work. When faced with a high order QN the designer has twooptions: direct implementation or perform model order reduction. Due to advanced filter technology,it is possible to directly implement high ordered QN. For example, controllers are implemented usingChapter 2. Control design step of the tuning method 17many-tap finite impulse response (FIR) filters in {BBB88] (where N = 30) and {0B90] (where N = 90).Alternatively, model order reduction can be performed on QN. If the basis functions (2.9) or (2.10) areused, QN composed with a bilinear transformation results in an FIR model for which simple model orderreduction methods exist [HJN9O, GKL89]. Hence, due to its simple structure it is possible to implementhigh order QN; the reader should therefore not be alarmed by the number of controller parameters usedin the designs presented later in this thesis.2.3 Controller design problemHaving specified the controller parameterization, performance functionals and specifications are introduced to describe the controller design problem. A detailed description of the performance functionalsis left until Chapter 3. As well, an interactive optimization-based procedure for the controller designproblem’s solution is presented.2.3.1 Performance functionalsPerformance functionals, convex in the controller parameter x, quantify some aspect of closed-loop;d- performance with lower values corresponding jr bc or tighter designs. The performance functional/ 1‘ \1/2RMS(X)= ( lim—J u(x,t)2dt) , (2.13)T-+c,oT jwhere u is the controller output signal due to a step reference input (with m = d = w = 0), measuresthe root-mean-squared of u and is constrained to prevent potential damage to the actuator causedby excessive heating. It is important to note that the performance functionals considered are convexbecause of the controller parameterization chosen. This is because the closed-loop transfer matrix (2.5)is a convex function of Q which is itself parameterized linearly in the controller parameters. Therefore(2.5) is a convex function of x and consequently functionals measuring some closed-loop quantity suchas the RMS value of u in (2.13) are also convex in a,.2.3.2 Performance specificationsOne method, due to [ZA73], of formulating a controller design problem is as a set of functional inequalitiesthat must be satisfied. The designer adjusts the constraints until desirable performance is achieved. AChapter 2. Control design step of the tuning method 18more flexible approach is to divide the performance specifications into two classesbkQc) ak,ke {1,...,n}, (2.14)I4(x) bk,kE{l,...,m},ak,bkER, (2.15)where b and are performance functionals [NT86]. The inequalities (2.14) are soft constraints or specifications about which the designer is flexible; the inequalities (2.15) are hard constraints or specificationsthat must be met. An obvious example of a hard constraint is stability of the system, whereas an example of a soft constraint involves the gain margin. Small violations of the specification that the positivegain margin be greater than 5dB typically do not affect the quality of the design and can be tolerated.As opposed to the formulation in [ZA73], the approach of dividing the performance functionals into twocategories addresses the designer’s flexibility with many specifications.Hence, the controller design problem is to find the set of parameters x which satisfy the inequalities(2.14) and (2.15). (Assuming the engineering specification can be formulated as these inequalities.)For every design there are two possibilities: either the inequalities (2.15) cannot be satisfied and noacceptable design exists, or (2.15) are satisfied and (2.14) are not satisfied and the designer must decidet on a best solution. The design procedure discusstH .Section 2.3.4 allows the designer to solve thiscontroller design problem.As mentioned in Section 1.2.1, in many cases it is difficult to formulate performance specifications asfunctional inequalities (2.14) and (2.15). For instance, specifications are often given vaguely in words suchas “the step response must contain little oscillation”. Additionally, in many cases there is insufficientknowledge about the system to specify exact levels of performance a priori. The interactive designprocedure described in Section 2.3.4 gives the designer a way of obtaining a set of inequalities that bestdescribe the optimal solution.2.3.3 Form of the optimization problemThe general form of the optimization problem sol ve&x determine whether the specifications (2.14) and(2.15) are achievable ismm max{’fr)I &‘(x) <O}. (2.16)The functional ‘ is defined as= max{(x), k.e {1, . . .,Chapter 2. Control design step of the tuning method 19the functional “ is defined as= max{(x), k C {1, ..and the normalized objectives and constraints are given by— #k—Gk ,— Vk—Gkk Bk—Gk’ k — Bk—GkBy chosing the weights Gk and Bk in (2.17) as “good” and “bad” values respectively for the performancefunctionals, the designer assigns a relative importance among the objectives c5k Although other choicesfor the normalization exist, the reason for using the one in (2.17) is because it tends to improve theobjectives throughout the optimization the way the designer wants them to trade off [NT86j. This canbe understood by quantifying the designer’s degree of satisfaction which is large (small) when the designeris satisfied (dissatisfied) with the current values of the objectives. The normalization (2.17) provides twolevels of designer satisfaction at which the value of the normalized objectives are equal (0 (1) when theobjective reaches its “good” (“bad”) value). Note that because of the way the weights Gk and Bk werechosen, it is likely that the level of designer satisfaction is linear in the normalized objective. That is, ifallthe normalized objectives take on the :tiue 1/3, the level of designer satisfaction is typically 33% awayir from its “good” value. A linear relationship between Ui;: normalized objectives and the designer’s degreeof satisfaction means that objectives decrease the way the designer wants them to trade off throughoutthe optimization. On the other hand, for the case when the relationship is nonlinear some interactiveadjustment of the weights may be necessary when the design specifications (2.14) are not satisfied. Forexample, if the current weights leads to an optimal objective value that is too large then its good weightis decreased (or equivalently its bad weight decreased)It should be noted that the designer is comfortable with adjusting the good and bad weights since theyhave a predictable effect on the resulting design. This contrasts with the nonintuitive weight selectionwhich can be required for the analytic dsign methods discussed in Section 1.2.1.2.3.4 Tuning procedureThis section presents the procedure for solving the controller design problem specified in Section 2.3 bysolving a sequence of constrained minimax problems of the form (2.16).ProcedureChapter 2. Control design step of the tuning method 201. Select N, the number of controller parameters, and the basis functions used in the function2. If the design has hard constraints, solve= argmin{’(x), x Cwhere Gk = bk, Bk = 2bk, k = 1, . . ., in. The procedure stops if the hard constraints are notachievable. That is b’(x0t) > 0. If a’vconstraints are satisfied and there are no soft constraintsthen the procedure stops. Otherwise let a =3. Starting at x0, solve= argmin{’(x)b’(x) 0, r e RN}, (2.18)where initially Gk = 0k, Bk = 2ak, Ic = 1,. . . , n. If ‘(x0) 0 then all soft constraints aremet and the procedure stops. Otherwise select the objective with the most undesirable value.(Typically it is the functional with the largest normalized value which competes with at least oneobjective.)4. Reassess what “good” and “bad” values of the selectee objective are and adjust its good and badweights. Repeat Step 3 with x0 = x00 until the b-t tradeoff between objectives is achieved.Note that if the performance of the system is undesirable it is possible that N in Step 1 was chosentoo small, and the procedure is repeated with larger N. Also, in Step 3 it is not always necessary tofind a truly optimal solution. If the designer observes a feasible point with a negative cost function orsatisfactory performance then optimization can be terminated. As well, if in Step 3 performance has notimproved significantly after a few iterations the designer can stop the optimization run and select theworst objective. In this case = Tnop in Step 4, where is current value of the design vector.2.4 Performing a design on PAs discussed in Section 1.2.2, one drawback of the traditional approach to controller design is that plantmodeling errors lead to conservative performance. To motivate the tuning method which computesperformance functionals from plant measurements, this section shows how the controller parameterizationin Section 2.2.1 and a optimization-based design with an exact plant model P gives higher or equalperformance to a design with an incorrect plant model.Chapter 2. Control design step of the tuning method 21The performance functional was introduced in Section 2.3.1 as a function of the set of controllerparameters x. In this section it is a function of the plant model P used to perform the design, F, andQ. Hence the notation (P, F, Q) is the evaluation of the performance of the system involving F, f(contained in the controller), and Q.LetP1 : inf çi(P,P,Q),QERH= infA(F,P,Q),QERfLand definePk mm {(ppQN(x))IIIQNII — bN O},xEwhere {bN} E R+ and bN oo (the sequence diverges). It is a fact that sequences of optimizationproblems Pj with increasing number of controller parameters N can approximate arbitrarily well theinfinite-dimensional problem P1 [PS89]. Equivalently, if x is a solution that achieves the minimum inproblem P, then {QN(x)} is a sequence of minimizers for F1, i.e.urn (P, p ON(XN )) = 71.N—The following theorem shows that the performance of any system containing P can be achieved withknowledge of an exact plant model P = P.Theorem 1 For any P,P E RH and Q E A, where A = {Q E RHIH(P,K(P,Q)) is internallystable } there exists an R RH such that(P,, Q) = (P, F, R).Proof: Let P,P e RH, and let R=1? is stable since Q E A. Substituting R into (2.7)for Q we getPh PQ) -PQ PQH(P,K(P,R))=(1—Q(P-P))’-QP-Q Q , (2.19)-P(1-PQ) FQ-1 1—PQhence (P,P,Q) = cb(P,F,R). .Chapter 2. Control design step of the tuning method 22Letting Q0i = arginf{q5(F,F,Q)Q C RH0}and Qopt2 = arginf{(F,P,Q)IQ C RH} (assuming these infima are attained), we have shown thatqS(F, F, Qopt2) = b(F, F, It) 6(F, F, Q0i) = 7iEquivalently, when designing for F, using an exact plant model always gives better or equal performanceto a design performed on an incorrect model.To illustrate this idea, consider a simple design problem (with analytic solution) performed on theexact plant modelF1: inf (F,F,Q),whereqS(F,F,Q) IIWi1terIIoo,1—sF(s)= (1+s)2’P(s) = F(s),Wi(s) = (1+s)tDefining Q0i = argirif{qS(F,F,Q)IQ C RH } then qS(F, P,Qopn) = 1/2. Now consider the designproblem on the incorrect modelF2: inf sF,F,Q),QE RH0.where( 1.5—s— (1.5+s)(1+s)and defining Qopt2 = arginf{q5(F,F,Q)Q C RH,}. The closed loop performance of K(F,Q02)isinferior to that of K(F, Q0i) which agrees with the conclusion obtained above. That is,(F, P, Qopt2) = l4i 1 — FQOP:2 = 2/3> 1/2.1-Q02(F-F)To summarize, since performance functionals of the’ ‘ontrol design step are computed ftom measurements of F and assuming the values of the estimated I unctionals are exact, the tuning method has aclear advantage in performance over the design performed off-line on an incorrect model.Chapter 2. Control design step of the tuning method 232.5 Optimization algorithmsTwo optimization algorithms were investigated to solve the optimization problem (2.16): the ellipsoidalgorithm, and a modified sequential quadratic programming (SQP) method. This section discusses theadvantages and disadvantages of both algorithms.Ellipsoid algorithmRecently the ellipsoid algorithm, which is thoroughly discussed in [GLS88 and the references therein, isreceiving much attention in control system lzarature (see for example [BB91, Z093, BCFB93]). Roughlyspeaking, the ellipsoid algorithm generates a sequence of ellipsoids decreasing in volume; the minimizeris contained in the intersection of this sequence. At each iteration the algorithm requires the value ofthe objective and its subgradient. (See Appendix A for the definition of a subgradient.)The ellipsoid algorithm has the following merits:• It is simple to implement; the number of lines of code required are a fraction of that required forother algorithms [EK85].• The functionals minimized need only be convex so that the algorithm is guaranteed to converge.• When the functionals are nonconvex exam eles have shown that in many cases the algorithm findsthe global solution [EKS5].However, the algorithm has disadvantages:• It requires a large number of functional and subgradient evaluations to obtain a small error in thesolution [EK85].• It is not a descent method; although at each iteration the distance to the solution decreases,functional values may increase. This is because the negative of the subgradient of a function isa descent direction for the distance to the minimizer and not a descent direction for the functionitself. In addition, particularly on the first few iterations when the initial ellipsoid is large, thevalue of objective can vary wildly.Sequential quadratic programming methodThe modified sequential quadratic programming (SQP) method implemented in [GraOO] has the followingbenefit:Chapter 2. Control design step of the tuning method 24• It requires less functional and gradient evaluations than other algorithms to obtain small errors inthe solution [EK85].The algorithm has the disadvantage:• The cost and constraint functional must be of the formmax{cbk(x),1c C {1,...,n}},where c5k must have continuous first and second derivatives for convergence to be guaranteed[Gra9O]. This implies that the convergence of a problem minimizing the L1 norm (which is non-differentiable) of a response of the system canuo be guaranteed using the SQP method.An exampleThe behaviour of the two algorithms is illustrated with the following simple convex problemmm max{fri — 1)2, (z2 — 1)2g(x) < 0},xeR2gfr) = x2 + ni — 1. (2.20)The solution to (2.20) is x0 = [0.5, 0.5]. The ellipsoid method, with an initial sphere of radius 106centered at [0,0], (using deep-cuts in both the constraint and objective iterations) finds the solution in23 iterations. Figure 2.7 shows several iterations .$gorithm with the values of the cost functionlabeled. Note that for the second iteration the value ol the objective functional is large (1.09 x io). Onthe other hand, the SQP method converges in only 5 iterations with a maximum objective value over alliterations of 8. The absolute error in the cost and constraints functionals is 10—6 and the absolute errorin the solution is 10i.Simulations have shown that for more useful control system design examples the SQP method stillconverges in less iterations than the ellipsoid algorithm. For this reason it is used in the simulation examples presented. However, either algorithm can be usEd and when discussing the performance functionalsin Chapter 3 the ellipsoid method allows a wider range of functionals to be considered.2.6 Internal model controller structure: Unstable plantAs mentioned in Section 2.2.1, when the plant model is exact (P = F), then a necessary and sufficientcondition for an IMC controller to provide internal stability is for F and Q to be stable. Hence the IMCChapter 2. Control design step of the tuning method 25C%IxFigure 2.7: Three iterations of the ellipsoid algorithm for the example problem (2.20). Note the largevalue of the cost function at the second iteration.structure shown in Figure 2.4 can not be used for unstable P. However, the IMC controller can stillbe used in a two-step design procedure if a nominal stabilizing stable controller exists (i.e. the plantis strongly stabilizable). (A condition for checki. bether a plant is strongly stabilizable is given in[DFT9O].)The first step of the design procedure is to obtain a nominal stable controller K0 that stabilizes P.The second step is to modify the dynamics of the stabilized system using the IMC controller structureas shown in the Figure 2.8. The expression for the controller isK(P, Q) = K0 + Q(1 - H(P)Q)1,Q e RII,whereH(P) = P(i + KP)1. (2.21)Although not obvious, because the nominal controller is stable the two-step design procedure is over allstabilizing controllers when P = P (see [Mac89] or Appendix B for an explanation).Finally, it is important to note that when P = P the closed-loop transfer matrix in terms of P, K0-500 0 500 1000x(1)Chapter 2. Control design step of the tuning method 26Figure 2.8: The IMC controller K(P, Q) = K0 + Q(1 — H(P)Q)1 for a strongly stabilizable plant.and Q isH(P,K(P,Q)) = (1+PK0)2P((1 + PK0)- PQ) -PK0(1 + PK0) - PQ PK0(1 + PK0)+ PQ-PK0(1 + P1<0) - PQ -(Ko(1 + P1<0)+ Q) Ko(1 + P1<0)+ Q—P((1 + PK0) - PQ) “t1 + PK0 — PQ) 1+ PK0 - PQ+an affine function of Q as for the stab1 plant.Chapter 3Computing the performance functionals in the control design stepThis chapter describes the performance functionals and their subgradients’ in detail. It gives the stepsfor computing the functionals and their subgradients for all entries of the closed-loop transfer matrix(2.5) from measurements of the open-loop system shown in Figure 3.9. As well, the amount of computation required to evaluate the functionals using a model is compared with that using measurements.Throughout the chapter it is assumed that the responses measured are exact, and to reduce the notationQ will be used instead of QN.U,HQix) H U3Figure 3.9: Open-loop system used to evaluate performance functionals and their subgradients usingmeasurements.3.1 The step response functionalConsider the step response functional= max max{s(x,t) a(t),,BQ) — s(x,t),0}, (3,22)where s is the unit step response of the transfer function H which is an entry of (2.5) {BB91]. Thefunctions a and 3 are piecewise continuous upper and lower bound functions chosen by the designer.Figure 3.10 gives an example of a step response and typical bounding functions a and fi such thatt(x) <0.‘In this chapter it is assumed that the ellipsoid algorithm is used. For a definition of a subgradient see Appendix A.27Chapter 3. Computing the performance functionals in the control design step 28I I I I I I I P Ict(t)‘/ /zY/O.8 0E it0.2p3(t)0%/%0 1 2 3 4 5 8 7 8 9 10Time (sees)Figure 3.10: Typical bounding functions a and fi for the step response functional 4st, and a step responses satisfying (x) 0.The functional is used to meet specifications involving overshoot, rise-time, “steady-state” error,and settling-time of the step responses of the entries of (2.5). As well, it can be easily generalized toinputs other than steps, and in Chapter 5 it is used to bound the amplitude of the system’s impulseresponse.A subgradient of can be determined using the rules in Appendix A. Letting to be a time at whichthe maximum occurs, if st(x) = s(to, x) — a(to) > 0 or the upper bound is violated then(st(x) = V(s(x, to) — a(to))= vr(L’ (:‘) (to))L1 (fjjfJ) (to)= . , (3.23)L1 (“) (to)where Hk is the partial derivative of H with respect to xk. If sj(x) = 9(t0) — sfr,to) > 0, or the lowerChapter 3. Computing the performance functionals in the control design step 29bound is violated then= V(i3(to) — s(x,to))[ —L (4ü) (to)=. (3.24)—L’ (!J21±I) (to)Finally, if &tfr) = 0 then= 0. (3.25)It is important to note the possibility there be several maximizers of (3.22). That is, there may beseveral times at which the maximum of the functional is attained. In this case any of the times can beused to compute a subgradient as described above.3.1.1 Using measurements to evaluate the step response functionalThe following steps are used to evaluate (3.22) from measurements of the open-loop system. Step 1 isperformed before optimization begins, and Steps 2—5 are performed at each iteration of the optimizationalgorithm.1. A step input is applied at ‘i (see Figure 3.9) with Q = 1 and the step response of P is measuredat u3 for 0 <t <tmar.2. A step input is applied at u1 and the step responses of Q and PQ are measured at it2 and it3 for0 tmaa,.3. The response measured at it3 in Step 2 is applied at it1 with Q = 1, and the step response of P2Qis measured at it3 for 0 <t tmar.4. The step responses of P(l TQ), —P(1 — PQ), —FQ,—Q, PQ — 1, and 1 — PQ are determinedfor 0 t trfl by adding £he step respon htained in Steps 1—3.5. The functional is computed using (3.22) where s is one of the step responses (depending on whichclosed-loop transfer function is being considered) obtained in Steps 2 and 4. A time, to, at whichthe maximum occurs is determined.Chapter 3. Computing the performance functionals in the control design step 303.1.2 Using measurements to evaluate a subgradient of the step response functionalAssuming a non-zero subgradient exists, the following steps are used to evaluate a subgradient of (3.22).Steps 1—3 are performed before optimization begins, and Step 4 is performed at each iteration of theoptimization algorithm.1. Let k = 1. The kth controller parameter is set to 1 and the remaining controller parameters areset to zero, a step input is applied 4: (see Figure 3.9), and the step responses of Bk and PBkare measured at u2 and u3 respectively for U I < trn.2. The response measured at it3 in Step 1 is applied to it1 with Q = 1, and the step response ofP2Bkis measured at it3 for 0 t < t,flax.3. lv = lv + 1. Repeat Steps 1 and 2 until lv = N + 1.4. If the upper bound is violated, (3.23) is computed by selecting one of the vectors of step responses(depending on which closed-loop transfer function is being considered) computed in Steps 1 and2 at to. If the lower bound is violated, (3.24) is computed by selecting the negative of one of thevectors of step responses (depending on which closed-loop transfer function is being considered)computed in Steps 1 and 2 at to.3.1.3 Computational costUsing a model the cost, measured in floating point operations (flops), associated with computing (3.22)for ft = H5,. and N = 20 before optimization begins versus the number of states in a plant model isshown in Figure 3.11. In comparison, using measurements requires no computation.At each iteration of the optimization algorithm, the relative cost of computation versus the numberof controller parameters is shown in Figure 3.12. Since the relative cost is greater than 1 for all valuesof parameters, the cost of computing (3.22; is reducd i.; using measurements. The reason for thiscomputational advantage is that when using nrasi; :qinents the only computation required is to subtractthe measured step response from the bounding functions a and jY. It will be shown that particularlyfor performance functionals involving time domain response, using measurements substantially reducesthe amount of computation required. The details of how Figures 3.11 and 3.12 were generated are inAppendix C.Chapter 3. Computing the performance functionals in the control design step 310C0a)a.0C0a0)C0LIFigure 3.11: The cost in flops of computing the step responses of PBk, ic = 1,.. . , 20 versus the numberof states in the plant model. These responses are required to compute the functional before optimizationbegins when using a plant model.3.1.4 A step response functional exampleThis example illustrates how rise-time, overshoot, and “steady-state” tracking error specifications on theoutput step response can be satisfied using 4. The plant isF(s) = , (3.26)(52+065+1)(s+1)and the specifications are a rise-time of less than 4 seconds, no overshoot, and a “steady-state” errorof .1%. The uncompensated and compensated plant’s step response, and the bounding functions areshown in Figure 3.13 using N = 20 and the basis functions (2.9) with p = 1 in Q. Clearly, the desiredperformance is achieved.3.2 Norms of particular response functionalsDifferent norms of signals can be combined with different inputs giving rise to a wide range of usefulperformance functionals. The following is a list of functionals and their subgradients that are norms ofresponses to a particular input, and can be used in the tuning method. The method used to evaluateChapter 3. Computing the performance functionals in the control design step 32U)SCUa)Figure 3.12: The relative cost of computing the step response functional at each iteration of the optimization algorithm using a model relative to using measurements versus the number of controller parametersN. The relative cost is greater than one for all N.the functional and its subgradient from measurements of the open-loop system in Figure 3.9 is similarto that used for the step response functional and is only outlined briefly.Throughout this section a particular input is assumed applied to the system and the resulting systemoutput is denoted by y. As well, it is assumed that y is such that the functional exists.The “L norm” functionalThe “L norm” functional is used to measure peak values of a signal and is defined as#iGr)= sup Iy(x,t)I. (3.27)tcR+Although (3.27) is not differentiable, a subgradiellt can be specified using the rules in Appendix A as(p1x) = sgn(y(x,t0))Vy(,x) (3.28)2 4 6 8 10 12 14Number of controller parameters16 18 20where t0 is a time at which the maximum in (3.27) is attained.Chapter 3. Computing the performance functionals in the control design stepa)-DaE33Figure 3.13: Output step response of compensated plant (solid), uncompensated plant (dash-dot) andbounding functions (dash) for Example 3.1.4. The step response meets the design specifications on therise-time, overshoot and “steady-state” error.Time (secs)Chapter 3. Computing the performance function als in the control design step 34The “L2 norm” functionalThe L2 norm of a signal can be interpreted as the square root of the total energy of the signal. The “L2norm” functional is defined asp2(X)=Y(, t)2dt) (3.29)which is differentiable, and its subgradient is equal to its gradientp2(X)=fY(x,t)vY(x,t)dt. (3.30)p2(X) jThe “Li norm” functionalThe “L1 norm” functional is defined asp3(X)= j y(,t)Idt. (3.31)Although (3.31) is not differentiable, a subgradient can be determined as3(x) = jsgn(y(x,t))Vxy(x,t)dt (3.32)[BB91].Note that for persisting output signals the values of (3.29) and (3.31) are infinite, and we define theaverage-absolute value and the RMS value functionals for signals that do not decay to zero.Average-absolute value functionalThe average-absolute value functional is4(X) = urn j y(x, t)Idt. (3.33)A subgradient of (3.33) is=1 j sgn(y(x, t))Vy(x, t)dt (3.34)[BB91].RMS value functionalThe RMS value of a signal is/ 1 T\i/2RMS(X) = lim f Y(Xt)2dt) , (3.35)Too oChapter 3. Computing the performance function als in the control design step 35the eventual average size of the signal and can be interpreted as a measure of average power of a signal.The functional (3.35) is differentiable and its subgradient is equal to its gradient1 1 fT(RMsfr) = lim—y(x,t)Vy(x,t)dt. (3.36)IRMS(X) T—*c’o T.‘o3.2.1 Using measurements to evaluate the norm of a particular response functional andits subgradientThe steps used to evaluate (3.27), (3.29), (3.31), (3.33), or (3.35) from measurements of the open-loopsystem are almost identical to those for the step response functional (described Section 3.1.1) and arenot repeated. The differences being that a particular input is applied instead of a step, and one of (3.27),(3.29), (3.31), (3.33) or (3.35) is used to evaluate the functional.Similarly, the procedure for evaluating a subgradient of the norm of a particular response functionalis almost identical to that of a step response functional and is not repeated. The differences being thatthe particular input is applied instead of a step, and one of (3.28), (3.30), (3.32), (3.34), or (3.36) is usedto evaluate the subgradient in Step 4.3.2.2 Computational costThe functional (3.35) is used to illustrate the computational cost of evaluating the norm of a particularresponse functional. Before optimization begins, the cost is approximately equal to that required forthe evaluation of the step response functional. Hence, Figure 3.11 shows the cost versus the number ofstates in the model. Also, as with the step response functional, it is important to note no computationis required when using measurements.At each iteration of the optimization algorithm the cost of evaluation using a plant model relativeto using measurements versus the number of controller parameters is shown in Figure 3.14. As with thestep response functional, the relative cost is greater than one for all values of design parameters becausethe particular response is obtained without any computation ftom measurements. The details of howFigure 3.14 was generated are in Appendix C.3.2.3 An example of using the “L norm” functionalTo illustrate how the norm of a particular response functional is used, consider Example 3.1.4 in whicha constraint was imposed on the output step response. Since there was no specification on the size ofChapter 3. Computing the performance functionals in the control design stepU,0VCua)36Figure 3.14: The relative cost of computing the functional (3.35) at each iteration of the optimizationalgorithm using a plant model relative to using measurements versus the number of controller parametersN. The relative cost is greater than one for all N.the actuator effort a large peak value for u resulted. This example shows how the “L, norm” functionalconstrains the peak value of the actuator output. The design specifications are1(x) = < 0, 2(X) = ifr) 1, (3.37)where the bound functions associated with the basis functions for Q, and N are the same as inExample 3.1.4.It is assumed that the designer is flexible with the specifications (i.e. the inequalities (3.37) are softconstraints) and, in an attempt to reduce 2 to 1, the designer selects the weights G1 = .001, B1 = .002,= 1, and B2 = 2. Optimization yields minimum raw functional values of 9.98 x 1O and .998for 4i and 2 respectively. At this point the designer explores how #2 can be decreased further whileaccepting more violation of the constraint on The weights G2 and B2 are set to .5 and 1 respectivelywhich result in minimum functional values of 1.99 x 1O and .994 for and #2 respectively. Clearly,relatively large increases in 4 are required to achieve relatively small reductions in c2• In a last attemptto reduce #2 B1 is increased to .01; the resulting values of i and #2 are .00982 and .979 respectively. Atthis point the procedure stops since further increase in is considered unacceptable. The best outputNumber of controller parametersChapter 3. Computing the performance function als in the control design step 37a,-oza.EFigure 3.15: The output step response for the example in Section 3.2.3 (solid) and bounding functions(dash-dot). Note the step response’s violation of the lower bound function at “steady-state” (a 1%tracking error).step response is shown in Figure 3.15 and should be compared with that in Figure 3.13. Note the 1%violation of the lower bounding function at “steady-state”. The step response of the actuator output isshown in Figure 3.16 along that in Example 3.1.4. Note the peak value of u is reduced dramatically bya factor of about 15 at t = 0.3.3 The “peak gain norm” functionalConsider the performance functional ç6pkgn= j I(x,r)dr, (3.38)where Ii is the impulse response of an element of (2.5) [BB91J. It can be shown that 4pkgn is equal toIIftwIIoosupIIwII0 IIwIIthe peak gain of the transfer function H, which can be interpreted as the worst-case peak value of thesystem output for all input signals with peak values less than or equal to one. Note that (3.38) is the “L1norm” of a closed-loop impulse response and therefore a special case of the functional p2 (see (3.31)).10Time (secs)Chapter 3. Computing the performance functionals in the control design step 3816141210’•8’6i2 ,:4:34 5Time (secs)Figure 3.16: Controller output step response in example in Section 3.2.3 (solid) and for Example 3.1.4(dash). Note the large decrease in the actuator output peak value due to the constraint onA subgradient of (3.38) is(pkgn(X) = jVIi(sr)Idrf’° sgn(ui(c, r))hi(r)dr(3.39)j° sgn(I(x, T))hN(T)dTwhere hk is the partial derivative of Ii with respect to x.3.3.1 Using measurements to evaluate the “peak gain norm” functionalThe following steps are used to evaluate (3.38) from measurements of the open-loop system. Step 1 isperformed before optimization begins, and Steps 2—5 are performed at each iteration of the optimizationalgorithm.1. An input of “unit-area” (i.e. f u(t)dt = 1) with large magnitude and short duration is appliedat u1 (see Figure 3.9) with Q = 1, and the impulse response of P is measured at u3.Chapter 3. Computing the performance functionals in the control design step 392. An input of “unit-area” with large magnitude and short duration is applied at u1, and the impulseresponses of Q and PQ are measured at u2 and u3.3. The response measured at u3 in Step 2 is applied to u1 with Q = 1, and the impulse response ofP2Q is measured at u3.4. The impulse responses of P(1 — PQ), —P(l — PQ), —PQ,—Q, PQ — 1, and 1— PQ are computedby adding the impulse responses obtained in Steps 1—3.5. The functional is computed by performing the integration in (3.38) where h is one of the impulseresponses (depending on which closed-loop transfer function is being considered) obtained in Steps 2and 4.3.3.2 Using measurements to evaluate a subgradient of the “peak gain norm” functionalThe following steps are used to evaluate the subgradient (3.39) of the “peak gain norm” functional.Steps 1—4 are performed before optimization begins, and Step 5 is performed at each iteration of theoptimization algorithm.1. Let k = 1. The kth controller parameter is set to 1 and the remaining controller parameters areset to zero. An input of “unit-area” with large magnitude and short duration is applied at Ui(Figure 3.9) and the step responses of Bk and PBk are measured at U2 and t13 respectively.2. The response measured at u3 in Step 1 is applied at u1 with Q = 1, and the impulse response ofP2Bk is measured at u3.3. The impulse responses of —Bk , —PBk, and —P2Bk are computed directly from the responses inSteps 1 and 2.4. k = Jr + 1. Repeat Steps 1—3 until Jr = N + 1.5. The integration in (3.39) is performed, where h is the impulse responses used in the computationof the functional (obtained in Section 3.3.1) and Vh is one of the vectors measured in in Steps 1—3(depending on which closed-loop transfer function is being considered).Chapter 3. Computing the performance functionals in the control design step 403.3.3 Computational costBefore optimization begins, the cost of computing the impulse responses is approximately equal to thecost required for the step response functional. Hence, Figure 3.11 shows the cost versus the number ofstates in the model.The cost of evaluating (3.38) using measurements relative to using a plant model at each iterationof the optimization algorithm versus the number of controller parameters N is shown in Figure 3.17.As with the previous two functionals, computational savings are possible for more than one controllerparameter.0)0C.,a,>a,Figure 3.17: The relative cost of computing (3.38) at each iteration of the optimization algorithmusing a plant model relative to using measurements versus the number of parameters in Q. Note thatcomputational savings are possible for more than one controller parameter.3.3.4 An example of using the “peak gain norm” functionalThis example illustrates the usefulness of the “peak gain norm” functional in reducing the oscillation ina response. The problem involves the design specifications (3.37) of Example 3.2.3 and2 4 6 8 10 12 14 16 18 20Number of controller parameters= pkgn(2J) .3, (3.40)Chapter 3. Computing the performance functionals in the control design step 41a constraint on the peak gain of Hur. The basis functions for Q and N are the same as in Example 3.2.3.After a few design iterations the design weights are chosen as 01 = .001, B1 = .01, 02 = .5, B2 = 1,03 = .3, and B3 = 1. The resulting optimal output step response is shown in Figure 3.18 which showsincreased oscillation due to the new constraint on.The controller output is shown in Figure 3.19(solid line) which has considerably less oscillation than the the output for Example 3.2.3 (dash-dot line).At the optimal solution only c5i and 4cj are competing; hence 2(x0t) = 1.01 could be reduced furtherif desired.ci,VaE.<Figure 3.18: The output step response in Example 3.3.4 (solid line). Note the increased oscillation dueto the constraint on peak gain of Hur. The bounding functions associated with 4’ are the dash-dotlines.3.4 The “H2 norm” functionalConsider the “H2 norm” performance functional= (FEw1 + w2(1=(Efr,jw)Wi(jw)+W2(fw)F2dw), (3.41)where W1 and W2 are stable weighting transfer functions selected by the designer [BB9 1). The functionalmeasures the H2 norm of the transfer function HW1 + W2. A physical interpretation of iUw1w2Time (secs)Chapter 3. Computing the performance functionals in the control design step 4210Figure 3.19: The controller output step response in Example 3.3.4 (solid) and the controller output inExample 3.2.3 (dash-dot). Note the decrease in oscillation in the response due to the constraint on thepeak gain of Hur.is the RMS value of the output of H for an input modeled as stationary stochastic process with powerspectral density S(w) = Wi(jw)Wi(jw)*.If we assume ,bH2(x) 0, then (3.41) is differentiable with respect to x, and its gradient isH3(X)= 1 J Vft(x,jw)Wi(jw)+W2(j )Id2bH(x)21r= 1 J (((xjw)Wi(jw) +2(jw))*V(Hfr,jw)Wi(jw) +W2(jw)))dw—oof°°, R((fI(x,jw)Wi(jw) +W(jw))*Hi(jw)Wi(jw))dw= 1 : (3.42)1H2(X)rf° R((ft(x,jw)W1(jw +2(jw))*H(jw)W1j ))dL,Jwhere “i” denotes the complex conjugate and Hk is the partial derivative of H with respect to xk.3.4.1 Using measurements to evaluate the “112 norm” functionalThe following steps are used to evaluate (3.41) from measurements of the open-loop system. Steps 1and 2 are performed before optimization begins, and Steps 3 and 4 are performed at each iteration ofTime (secs)Chapter 3. Computing the performance function ads in the control design step 43the optimization algorithm.1. An input with enough frequency content is applied to input u1 (see Figure 3.9) with Q = 1; theresponse at u3 is measured, and the frequency response of P is computed using the “tn-spectrum”average technique discussed in Appendix C.2. Using the expression for Q and the weighting function W1, the frequency responses of Wi F,W1PBk, W1Bk, 117P, and W1P2Bk, k = 1, . . . , N are obtained.3. The frequency responses of the weighted closed-loop transfer functions W1P( 1 — PQ) + W2,-W1F(1 - PQ) + W2, W(1 - PQ) + 1472, W1(PQ - 1) + P172, W1PQ + W2, -WPQ + P172,W1Q + 142, and — W1Q + 1472 are formed by multiplying and adding the frequency responsesobtained in Step 2, the controller parameter x, and the weighting function 1472. For example,Wi(jw)P(jw)Q(x,jw) + W2(jw) = Wi(jw) (xkP(iw)Bk(iw)) + W2(jw).4. The integration in (3.41) is performed where ftfr, jw)Wi (jw) + W2(jw) is one of the frequencyresponses computed in Step 3 (depending on which closed-loop transfer function is being considered).3.4.2 Using measurements to evaluate the subgradient of the “H2 norm” functionalThe following steps are performed at each iteration of the optimization algorithm to evaluate (3.42) frommeasurements of the open-loop system.1. At each iteration of the optimization algorithm one of the frequency responses of the weightedclosed-loop transfer functions (depending on which closed-loop transfer function is being considered) computed in Section 3.4.1, Step 3 is multiplied by its derivative, one of the vectors of frequencyresponses computed in Section 3.4.1, Step 2.2. The subgradient is computed by performing the integration in (3.42) using one of the vectors offrequency responses computed in Step 1.3.4.3 Computational costFigure 3.20 shows the cost required before optimization begins using a plant model relative to usingmeasurements versus the number of states in the plant model. Computational savings are possible whenChapter 3. Computing the performance functionals in the control design step 44the number of frequency points used is large and the number of states in the plant model is large. Forexample, savings are possible when the model has more than 12 states and the response is computed at370 frequency points.At each iteration of the optimization algorithm, the relative cost versus the number of controllerparameters N is shown in Figure 3.21. Because of the simple method available for computing the functional with a plant model there is no computational advantage to using measurements. It is important tonote that note that (3.41) could be computed by integrating the impulse response of the system squared.Although this method would require less computation, in practice measuring the impulse response response of a system is more difficult than obtaining its frequency response. Hence, at the cost of morecomputation (3.41) is computed in the frequency domain. The details of how Figures 3.20 and 3.21 weregenerated are in Appendix C.I I I I I I I I I///1.0 //1.61.47al20 —0.80.6£ 6 8 10 12 14 16 18 20Number of states in plant modelFigure 3.20: The cost of computing the frequency response of PBk, k = 1, . . . ,20 at 370 (dash) and 92(solid) points using a plant model relative to using measurements before optimization begins. Note thatcomputational savings by using measurements are possible for 370 frequency points and when there aremore than about 12 states in the model.Chapter 3. Computing the performance functionals in the control design step 450.20.15Co0C)a)>. 01 Va)0.05 V V2 4 6 8 10 12 14 16 18 20Number of controller parametersFigure 3.21: The relative cost of computing (3.41) at each iteration of the optimization algorithm usinga plant model relative to using measurements versus the number of parameters in Q. There is nocomputational advantage to using measurements.3.4.4 An example of using an “H2 norm” functionalConsider the plantF(s) = , (3.43)(s2--.06s+9)(s+1)and the problem of minimizing the RMS value of the output due to a white-noise input applied at theplant input disturbance d. With N = 60 and the basis functions for Q as in Example 3.3.4 a minimumvalue of IIHydII2 = 1.29 is obtained. Figure 3.22 shows a sample of the “steady-state” output responsedue to a white noise excitation applied at d. The RMS value of the response is labeled with a dashedline. For comparison the theoretical minimum is 1.27 [DFT9Oj.3.5 The “Hc,, norm” functionalConsider the “H norm” performance functionalc5Hinf(x) = sup IW1(jw)if(x,jw) + W2(jw) (3.44)wcRChapter 3. Computing the performance fun ctionals in the control design step 46a)Va.2Figure 3.22: A sample of the “steady-state” output response to a white-noise excitation applied to theplant input disturbance d (solid curve) for the optimal solution in Example 3.4.4. The RMS value of thesignal is 1.29 and is labeled with a dotted line.[BB91]. By the Maximum Modulus Theorem, q5Hinf is the H,, norm of HW1+ W2 ([Fra87]) and allowsus to place hard bounds on the magnitudes of the frequency responses of the system. The H,,, normcan be interpreted as the worst-case RMS value of the output for all input signals with RMS values lessthan or equal to one.Although q5Hinf is not differentiable because of the supremum operation, it is convex and an expression for a subgradient can be determined using the rules in Appendix A.A subgradient of qSH1.r is(Hinf(X) = VrWi(jwo)H(x,jwo) + W2(jwo)I= Vr((Wi(jwo)&fr,jwo) + (jwo))(Wi(jwo)H(x,jwo) += 1V((Wi(jwo)H(x,jwo) +2(jwo))(Wi(jwo)Hfr,jwo) + W2(jwo))*)2çbjjf (a’)IR((Wi(jwo)k(x, iwo) + (jwo))*Wi(jwo)iti(jwo))= 1, (3.45)#Hinf (a’)Jt((Wi(jwo)ftQr, iwo) +2(jwo))*Wl(jwo)ftN(jwo))Time (secs)Chapter 3. Computing the performance functionals in the control design step 47where w0 is any frequency such thatc5Hjnf(x) = i’V1(j0)ft( w + W2(jwo)I. To obtain (3.45) it wasassumed that IWi(iwo)ftCiwo) + W2(jwo)I 0. Also, ilk is the partial derivative of H with respect to3.5.1 Using measurements to evaluate the “lI functionalThe following steps are used to evaluate (3.44) from measurements of the open-loop system. Steps 1and 2 are performed before optimization begins, and Steps 3 and 4 are performed at each iteration ofthe optimization algorithm. Note that the method is similar to the one in Section 3.4.1.1. An input with enough frequency content is applied to input u1 (see Figure 3.9) with Q = 1; theresponse at 123 is measured, and the frequency response of P is computed.2. Using the expression for Q and the weighting function Wi, the frequency responses of W1P,W1PBa, WiBi, and W1P2Bi, ic = 1, . . ., N are obtained.3. The frequency responses of the weighted closed-loop transfer functions Wi P( 1 — PQ) + W2,-W1P(1 - PQ) + W2, Wi(1 - PQ) + W2, W1(PQ - 1) + W2, W1PQ + W2, -W1PQ + W2,Wi Q + W2 and — Wi Q + W2 are formed using multiplication and addition from the frequencyresponses obtained in Step 2, the controller parameter x and the weighting function W2.4. The maximum magnitude of one of the weighted closed-loop frequency responses (depending onwhich closed-loop transfer function is being considered) obtained in Step 3 is computed and thelowest frequency at which the maximum occurs is w0.3.5.2 Using measurements to evaluate a subgradient of the “H norm” functionalThe following steps are used to evaluate (3.45) from measurements of the the open-loop system. Notethat the method is similar to the one in Section 3.4.2.1. At each iteration of the optimization algorithm, the conjugate of one of the weighted closed-loopfrequency responses (depending on which closed-loop transfer function is being considered) at wo(computed in Section 3.5.1, Step 3) is multiplied by its derivatives, one of the vectors of frequencyresponses computed in Section 3.5.1, Step 2.2. The subgradient is computed using (3.45), where (Wi(jwo)H(x,jwo) +w2(jw0))*W1H(jw iskth component of the vector of frequency responses computed in Step 1.Chapter 3. Computing the performance functionals in the control design step3.5.3 Computational cost48The relative cost required before optimization begins is shown in Figure 3.23. Computational savingsare possible when there is more than one state in the plant model.Figure 3.23: The relative cost of using a plant model as opposed to measurements to compute FEk,k = 1,. . . , 20 at 370 frequency points versus the number of states in the model. Computational savingsare possible for more than one state in the plant model.Because the method for evaluating the “II, norm” functional using the model and the measurementsare the same, the relative cost of computation during optimization is one. The details of how Figure 3.23was generated is in Appendix C.3.5.4 An example of using the “H norm” functionalThis example considers the problem of minimizing the H0 norm of Hd for the plant (3.43). The basisfunctions for Q, and N are the same as in Example 3.4.4. A minimum value of IIHyczIIoo = 1.00 isobtained. The resulting Bode magnitude plot of Hd is shown in Figure 3.24 with the minimum value ofthe functional labeled with a dashed line. Figure 3.24 illustrates clearly how 6Hinf can be used to placehard bounds on the magnitude of the frequency responses of the system’s closed-loop transfer functions.Chapter 3. Computing the performance functionals in the control design step 49102ici 1 0 I 2 3 4io2 10 10 10 10 10 10Frequency (rad/s)Figure 3.24: The Bode magnitude plot of H2d for Example 3.5.4 The maximum magnitude of thefrequency response is one and is labeled with a dashed line.Chapter 4Identification step of the tuning method4.1 Identification methodAs discussed previously, the IMC controller structure requires a model of the plant. As well, in order toguarantee desired closed-loop performance a quantitative bound on the uncertainty in the unmodeleddynamics of the plant is required. This chapter presents a parameter optimization-based method forcontrol-oriented identification in H,. As mentioned in Chapter 1, “control-oriented” implies the modeland the measure of its “error” are compatible with robust control design methods. “Identification inmeans uncertainty is measured in the frequency domain by a disk in the complex plane containingP(jw) and centered at P(jw) for all frequencies. This additive perturbation model of uncertaintyimplies the identification method is compatible with R-optimal controller design methods described in[Zam8l, Fra87]. More importantly, it is compatible with the design procedure described in Section 2.3.4which accommodates the H, type objective as a special case.As an example of how the tuning procedure guarantees robust stability, assume Q00 is designed tosatisfyW2HyrIIoo <1,where W2 is a bound on the relative uncertainty in the plant,P(jw)—P(jw)P w) <IW2Ow)I,Vw.If the identified model satisfiesIF—PII < ‘/IIQ0t ,then the system is ensured to be stable. Similarly, it can be shown that a specification for robustperformance can be guaranteed if the plant model satisfies a worst-case error bound.Therefore, the identification objective minimized isIIQ0(P—50Chapter 4. Identification step of the tuning method 51where Q° is the optimal design parameter obtained in the control design step described in Chapter 2. Weighting the model error by Q0,t makes the identification dependent on desired closed-loopperformance. This can be understood by considering the expression for H5,. in terms of F, P and Qo,,t,H5,. = (1— —As Q0(P — F) tends to zero then H51. tends to FQ0, the desired closed loop transfer function.Identification of the plant can be formulated as a convex optimization problem= argmin{max{cit(jwk)(F(jwk)— P”(v,jwk)) I v C R’,k C {1, . .where jiF isP”(vs)=Evk (Z2),A>O, (4A6)and where the frequencies wk, Ic = {1,. . . , n} and order of the model F are to be determined by thedesired worst-case model error using Theorem 2 which is presented in this chapter. The choice of theparameter A affects how many terms are required to obtain a certain accuracy in the solution. Simulationshave shown that a wise choice for A is the inverse of the plant’s dominant time constant.4.2 Deriving the bound for the unweighted model errorThe derivation of the model error bound is motivated by work done by ilelmicki et al. in the series ofpapers [11JN92, HJN9O, HJN91a, HJN89, HJN91c]. These papers present identification algorithms whichmap a set of noisy frequency response data of the plant into a model for which worst-case error boundsare provided. In particular this thesis uses results of the algorithms in [11JN89, llJN9la, 11JN92], whichare described next, to obtain an error bound for the optimization-based identification method presentedhere.The algorithm in [11JN89] involves computing a cubic spline to noisy frequency response data. Thecoefficients of the Fourier series of the spline are then used as the parameters of an FIR model. Anexpression for an Hc, model error bound in the presence of unknown bounded noise is derived. Since anFIR model is also used in the identification method presented, results in [11JN89] are useful in quantifyingthe model error.In [HJN91a] the model parameters are obtained in two steps. The first step is to compute thecoefficients of a truncated Fourier series of a linear spline interpolating a set of noisy plant frequencyChapter 4. Identification step of the tuning method 52response data. This step is identical to that in {HJN89] except that a linear as opposed to a cubic splineis used. The second step is to map the truncated Fourier series into a transfer function model by solvinga Nehari problem (see [FraS7] for a discussion of the Nehari problem). The results in [HJN91a] fordiscrete-time plants are extended to the continuous-time case in [11JN92] by composing the continuous-time plant with a bilinear transformation. The method presented in this chapter follows the work in[HJN92] by deriving the error bound in discrete-time.4.2.1 Transformation to the z-plane via the bilinear transformationMost of the results for identification in II have been obtained for discrete-time plants hence it isconvenient to compose P with a bilinear transformationPbl(z) = P1Ø(z)),whereb(z) = A—,A>o, (4.47)maps the closed unit disk onto the closed right half plane. This transformation is made because it iseasier to obtain a bound for the model error by working with Pb! instead of P.As well, it is assumed that a stable discrete-time system has a transfer function which is analytic onthe unit disk. That is, the z-transform of an impulse response h is taken asH(z) =to avoid using the cumbersome notation “r1”.4.2.2 Assumed knowledge about the plantIn order for the model error to be bounded in Chapter 1 it was assumed that the plant is exponentiallystable. For such plants there exist [RJNO2I:• p> 1 such that Pb! is analytic in D,, = {s C CJJsJ <p}. Equivalently, Pb! has no poles in the diskD.• M > 0 such that IPb!II,,p M, where IPb!IIcp = sup€D IPb!(z)I.• a continuous function (w) > 0 such that IP(iw)I <(w), Yw and lim. (w) = 0.Chapter 4. Identification step of the tuning method 53It is assumed that p, M, and 4 are a priori information about the plant which has the followingphysical interpretation: the parameter p is a measure of relative stability of P1q, is a bound on thehigh frequency “roll-off” rate of the frequency response of P, and M is the steady-state gain of Ps,, dueto an exponentially decaying complex sinusoidal input ri(k) = p_kei.4.2.3 Approximation error of a linear spline interpolantThe samples of the plant’s frequency response areP(jwk),k = 0,...,n— 1,PD(Ok)= 0,k=n,...,1/2, , (4.48)(P(.jwl_k))*, k = (1/2) + 1,.. .,l — 1.where wk = A tan(irk/l). Note that choosing the particular spacing of the frequency response data leadsto evenly spaced data on the unit circle after transforming the jw-axis using (4.47).A linear spline to the data (4.48) is denoted g1PD and is defined asg’PD(O)=PD(Ok)+(6—6k) [PD(9k)_PD(Ok+1)] o <e <9k+1,k=O,...,l. (4.49)The error between a function FbI and a linear spline g1PD is given byIJF -g’PDI <min{IIPiII0,()2(IIFj1[ + IIIIc)}.Using Cauchy’s integral formula we can express IIP’Il i terms of M and p, and a bound between gPDand FbI can be obtained and is presented as Lemma 1.Lemma 1 Let gP denote the linear spline interpolant to the exact extended frequency response dataas in (4.48), thenIIFbI — g’FD <max 4Mir , sup (Atan(6/2)) . (4.50)I l(p — 1) ee[2/l,] JProof: See [HJN92, Lemma 5.1] and [HJIN91a, Fact 5.3]. •RemarksNote that in (4.50) that if 1 is chosen such that I > 2n(1 + n), a > 0, the second entry in the maxtends to zero as n tends to cc since 27rn/1 —* ir and limw_c,o .(w) = 0. Also, the first entry in the maxtends to 0 since 1 tends to cc.Chapter 4. Identification step of the tuning method 54Also, for small a (1 >> 2n) the frequency response data points are closely spaced and the term 1M)in (4.50) is small. For large a (1 2n) we specify a large bandwidth over which the data is collected,and the term SUp9E[2flh1,,] (Atan(9/2)) in (4.50) is small. For a given number of data points n andthe optimal value of a is determined by minimizing (4.50) numerically.4.2.4 Error bound for truncating a Taylor seriesWhen the plant model (4.46) is transformed into the z-plane, we obtain the FIR modelP(v,z) = pF(v(z)) = EVkzk.The error in truncating the Taylor series of FbI will be required to obtain the bound on the model errorand is provided in Lemma 2.Lemma 2 Le TPj be he Nth order Taylor series truncation of FbI thenIIT1—Pb1(Ic NMp (p—i)Proof: See [RJN89, Theorem 3.1]. .4.2.5 An expression for the model error boundTo be of practical use, an identification method should provide a model which has desirable noiserejection properties. Noise arises from many sources including sensor noise corrupting measurements orunavoidable roundoff errors when implementing an algorithm numerically. To account for its effect, thefrequency response data is assumed corrupted by a term k which is unknown but bounded in magnitudeby M. The data in (4.48) becomes= PD(Ok)+e(Ok),kE 0,...,l, (4.51)whereC(Ok)t 0,k=n,...,l/2, . (4.52)An identification method for obtaining the data (4.51) and M6 is discussed in [HJN92] and Appendix F.Chapter 4. Identification step of the tuning method 55Theorem 2 Let P’ be a model of the form (4.46) where the Yk are determined by solvingP argmin{max{PD(jwk) — PF(v,jwk) I v E R’,k E {1, . . .,n}}}, (4.53)where l,n, ), and F are free parameters, thenll — P1100 M6 + FM+ max{ 4M , sup (Atan(9/2))}. (4.54)p (p — 1) l(p — 8E[2rn/1,r]Proof:rnax{I(PD(wk) - PF(wk)I k e {i,. .,i}} = IIgFD - FF1100=IIgP— P1100 + g’eco + lip —M 14M7rF +Me+maxjp (p—i) l(p—i)sup (Atan(O/2)). (4.55)9E[2irn/i,r] JLemmas 1 and 2 were used in (4.55) and the fact that IIg’eiI,x, <E.•RemarksNote that as n tends to oo and a > 0 is chosen so that 1 > 2n(1 + n), the error (4.54) tends toM6 + pF(p_l) Furthermore, as the number of model parameters F is increased, the error tends to M.Hence the bound on the error in the frequency response data Me places a lower bound on the uncertaintyin the identified model.Also, the weighted model errorlIQ0(P — 1t)II00, (4.56)where is the solution to (4.53) is bounded by,, ñF ,,.U opt 00 ‘Jopt 00using the sub-multiplicative property of the U00 norm. Hence Theorem 2 can be used to bound theerror in a model obtained by minimizing the weighted objective (4.56).Chapter 5Simulation examplesThis chapter presents two design examples illustrating the tuning method. An analytic method is usedto solve the first problem so that the performance of the tuning method can be compared with an exactsolution. The second example considers a benchmark problem involving a wider range of performancefunctionals and an unstable plant. Since the identification of unstable plants is beyond the scope ofthis thesis, only the control design step of the tuning method is presented. Throughout the chapter thefunctionals estimated from the open-loop system are assumed exact. As well, for display all figures arerounded to 3 digits.5.1 Example 1 Stable plantThe plant considered in this section is9F(s) (1+s)(s2+.8s+9)which is exponentially stable and has a lightly damped complex pole pair at s — .40 ± j2.97. Thedesign specifications are:• an output step response that is well damped and asymptotically tracks a reference step input (withd = 0, w = 0, and in = 0).• a (positive) gain margin greater than 10dB.One of the features of the tuning method is that because the control design step solves a convexoptimization problem it obtains “good” solutions to the design problem. This is shown in this exampleby comparing the method’s performance with a solution obtained using an analytic metbod. Since tbeanalytic method chosen can be interpreted as minimizing a cost composed of “H2 norm” functionals(see (3.41)), it is straightforward to compare its solution to one obtained using the tuning method.56Chapter 5. Simulation examples 575.1.1 An analytic solutionThe Loop Ttansfer Recovery (LTR) design method, discussed in Appendix E and [DS81, Mac89], is usedto obtain the analytic solution to the design problem. This procedure gives the designer a systematicway of adjusting the design weights of an Linear Quadratic Gaussian (LQG) problem (which is alsopresented in Appendix E) to shape the system’s frequency responses while guaranteeing a design withgood stability margin. Since the LTR method solves a series of LQG problems, the controller has thesame structure as that in the LQG problem. The LQG controller structure is shown in Figure E.54,where I<f is the Kalman filter gain matrix, K, is the optimal state feedback matrix, and A, B, and Care state-space matrices of the augmented plant with input Up and output Yp shown in Figure 5.25. InFigure 5.25, the measurement noise m and v are white-noise and uncorrelated in time, and the stabletransfer transfer function Wd is the spectral factor of the process noise d. It can be shown that the state-space equations (see (E.103)) of the augmented plant are of form required by the LQG method. Thereasons for considering this augmented plant: it gives the designer an understandable way of adjustingK1 by choosing Wd so that— C(sI — A)1K, (5.57)is a desirable open-loop gain for the system, and it gives us another interpretation the problem (discussedin Section 5.1.2) allowing us to compare the analytic solution with that of the tuning method.Designing Kj: the Kalman filter gain matrixAs mentioned previously, the choice of the spectral factor Wd allows us to shape (5.57) so that it is adesirable open-loop gain. At frequencies where the spectral factor of the process noise Wd is large relativeto that of the measurement noise, K1 is increased. Hence, to cancel the plant’s resonance Wd containsa “notch” at 3 rad/s, and to meet the tracking specification it is chosen large at low frequency. Afterexperimenting with the values of the poles and zeros, a Wd with the above specifications is obtained andits Bode plot is shown in Figure 5.26. The transfer function of Wd is1.23s2 + 1.23s + 10.6Wd(s)= s + 2.00s2 + 10.6s + .0211(5.58)The Kalman filter gain K1, is computed by solving (E.100) with W = I and V = I, and bysubstituting into (E.99). The Bode magnitude plot of the resulting —C(sI — A)’K1 is the solid line inFigure 5.27. It is important to note that the gain margin of (5.57) was not considered in the selection ofChapter 5. Simulation examples 58Wd because it is guaranteed to be infinite. The next step of the LTR procedure is to compute K, theoptimal state feedback matrix, so that the frequency response of —C(sI — A)1K is recovered as theopen-loop gain for a sufficiently wide range of frequencies.Designing K,, : the optimal state feedback matrixThe standard approach to computing K,, is to set Q = I and iteratively decrease p in R = p1, solve (E.98)and obtain K,, by substitution in (E.97) until KP has converged to —C(sI — A)1Kf over a sufficientlywide range of frequencies to satisfy the performance specifications. To illustrate this procedure, theBode magnitude plot of (5.57) (solid) and the open-loop gain for decreasing values of p is shown inFigure 5.27. In this example, the final value of p is obtained when the closed-loop step response showsno oscillation, asymptotic tracking is achieved, and the positive gain margin is greater than 10dB. Notethat p should not be decreased more than is necessary since, as shown in Figure 5.27, it unnecessarilyincreases the open-loop gain at high frequencies. Although a value of p could be determined the standardway described above, to facilitate the comparison of the solutions of the tuning method with the analyticmethod we defer its selection till Section 5.1.2.Figure 5.25: The augmented plant used in the LTR design method. The inputs m and v are white-noiseand uncorrelated in time.Chapter 5. Simulation examples 59a)0DC0)‘UFigure 5.26: The Bode magnitude plot of Wd, the spectral factor of the process noise d. The lowfrequency gain is large to meet the tracking specification, and a notch at 3 rad/s removes oscillationfrom the output step response.5.1.2 Open-loop tuning and another interpretation of the design problemWhen solving a multiobjective optimization problem, where the cost function is composed of severalobjectives, the designer is interested in solutions that lie on the boundary between achievable andunachievable performance. This set of points are known as the Pareto optimal points of the problem. Itis well-known that any Pareto point is the solution to a weighted minimax problem {BS8O]. For example,Figure 5.28 shows a Pareto optimal point for the problem involving the objectives = RMS(y)2 and= RMS(u)2,where u and y are the responses due to noise at m and i’. The point is obtained byminimizingA1RMS(y)2+A2RMS(u), (5.59)with A1 = 3.30 and A2 = 1. The resulting optimal objective values are = 14.7 and 2 = 147. Notethat the dashed line in Figure 5.28 has a slope —A2/A1 and is ratio of the change in 4’2 over the changein q on the tradeoff curve about the Pareto point. A controller minimizing the cost (5.59) is also aFrequency (rad/s)Chapter 5. Simulation examples 6011102101a) 0D 10CC)c -110102111 1 0 1 02 1 01 100 101 102Frequency (radis)Figure 5.27: The procedure for synthesizing the optimal state feedback controller K. The Bode magnitude plot of the open-loop gain KP converges to that of —C(sI — A)—1K as p is decreased.solution to the problem minimizingmax{)3RMS(y2,.)4RMS(u)2}, (5.60)with weights )3 = 10 and )i4 = 1.Hence, if we select good and bad weights G1, G2, B1, and B2 for the problem that minimizesf RMS(y) — G1 RMS(u)2— G2 5 61max’ G1—B ‘ G2—B . )to determine the optimal p for the LTR procedure described in Section 5.1.1, we can use the same weightsto perform the control design step of the tuning method which solvesmm max{q(x),qS(x)}, (5.62)ERwwhereçSi(x) = 421(z) + 422(2), 2(X) = 1)H23(Z) + 4a4(x),= IIWdHydII2, H2(X) = II1tyinII2, H23(X) = II14’dUudII2, ‘H24(Z) = 1111um112,‘ ‘ \‘\ ‘\ \\‘‘ \‘ \\\ \Chapter 5. Simulation examples 61and the performance of both the designs can be sensibly compared. The specific choice of good andbad weights and the optimal values of the objectives are given in Table 5.1, and the expression for thecontroller isK (— 9.78s5 + 38.8s4 + 238s + 498s2 + 1220s + 1084Opt\S)— s6 + 11.0s5 + 77.2s4 + 331s + 793s2 + 1464s+ 11.2(5.63)This choice of weights corresponds to a value of p .007c’JC’)c’J140Figure 5.28: The Tradeoff curve between the objectives i and qf2. The dashed area is the region ofachievable performance and a Pareto optimal point, RMS(y)2 = 14.7 and RMS(u)2 = 147, is labeled.This point is the minimizer of the objective (5.59) for the weights . = 3.3 , = 1 and (5.60) for theweights.= 10 , .\4 = 1. The dashed line has a slope equal to —1/A2= —3.30 and equals the ratio ofthe changes in the objectives about the Pareto point on the tradeoff curve.Table 5.1: Optimal values of the cost function for the analytic solution and design weights used (theobjective is (5.61) ).I RMS(y)2 I RMS(u)2 G1 I I B1 B2 I max{RM)_1 RMS(u)2—G3 IG—B ‘ G2—B j I[ 1.02 251 1 250 2 300 .0168 Ii’ RMS(y)2Chapter 5. Simulation examples 62Results of the control design stepThe basis functions considered are (2.10) and (2.11) since Hgm= Q must be strictly proper for its H2norm to be finite. The range of controller parameters is [2, 201.Convergence to the exact solutionTo investigate how well the exact solution can be approximated, the number of controller parameters isincreased for several choices of p in the basis functions (2.10) and (2.11). The value of the minimum costfunction versus the number of controller parameters N is shown in Figure 5.29 where the dashed linelabels the minimum objective value of .0168 for the exact analytic solution. Except for the choice p = 100with the basis functions (2.10) and p = 1 + jlO with the basis functions (2.11), a good approximationto the exact solution is obtained with few controller parameters.Since we know the exact solution to Q, some explanation can be given for the effect of the choiceof p on the rate of convergence observed in Figure 5.29. The poles of the exact solution lie in the setS {s E C — .4 < Js < —3, —5 < s < 5}; we expect faster convergence for choices of —p in or closeto S. Although in general this behaviour is observed it is evident that “errors” in the imaginary part ofp lead to slower convergence than “errors” in its real part.To get some idea of the computational complexity of the design, the number of functional evaluationsrequired for an error less than .01% in the cost function is shown in Figure 5.30 for the same choice of basisfunctions as in Figure 5.29. For all choices of p and basis functions, the number of functional evaluationsincreases with the number of controller parameters, and the average number of functional evaluationsis about 45 for N = 20. This amount of computation is clearly much higher than that required by theanalytic method which requires the solution of a few Lyapunov equations. Hence, the purpose of thisexample is to illustrate that good approximations to the exact solution can be obtained by tuning and notto convince the reader of its computational advantage over an off-line design. It is important to rememberthat especially when specifications are simple analytic methods can be computationally efficient designmethods. On the other hand, if the design problem involved more demanding specifications then thetuning method is more flexible. Also, the example is unrealistic because the most convenient wayto ensure asymptotic tracking with the tuning method is to use the step response functional (3.22).However, choosing this functional prevents us from comparing the performance of the tuning method.Chapter 5. Simulation examples 630a).000a)NCu202Figure 5.29: The optimal value of the cost function for problem (5.62) versus the number of controllerparameters N. The choice of basis in QN affects the rate of convergence to the exact value of the problemwhich is labeled with a dashed line. The solid lines are for basis functions (2.10), the dash-dot lines arefor basis functions (2.11), and the values of p are labeled.Tradeoff curve between l’i and ‘2To illustrate how the choice of good and bad weights determines the the relative importance betweenthe objectives, Figure 5.31 shows the tradeoff curve obtained by varying the weights with basis functions(2.9), N = 20, and p = 10. The values of three choices of weights used to obtain particular points onthe tradeoff curve are labeled. By increasing G2 and decreasing G1 (with B1 = 2G1 and B2 = 2G)the designer trades off an increase in 4 with a decrease in This tradeoff has the identical effectof increasing p in Section 5.1.1. In other words, it decreases the tracking performance of system whileincreasing robustness to high frequency plant perturbations by decreasing the value of the open-loopgain at all frequencies.101Controller parameters, NChapter 5. Simulation examples 64U,C0CuzCu>a,C00CLIFigure 5.30: The number of functional evaluations of the optimization algorithm required to obtainsolutions of (5.62) to within .01% of their optimal values versus the number of controller parameters N.The solid lines are for basis functions (2.10), the dash-dot lines are for the basis functions (2.11), andthe values of p are labeled.5.1.3 Results of the identification stepWe assume given the correct plant information p = 1.15, M = 0, (w)= (1+(/7))’ and M = .7. Aswell, arbitrarily A = 3 in (4.47). Figure 5.32(a) shows the Bode magnitude plots of , the upper boundon the “roll-off” rate of the plant, and P. Figure 5.32(b) shows a physical interpretation of M and p asa bound on the impulse response of Fbi. The optimal design parameter, denoted by Q, uses the basisfunctions (2.10), N = 20, and p = 1.Before the plant model is obtained the worst-case error bounds are computed using Theorem 2 todetermine an appropriate number of model parameters, and bandwidth of the plant’s frequency responsedata. First we consider the dependence of worst-case error on the number of model parameters F and thenumber of data points n. An optimization problem is solved, as discussed in Section 4.2.3, to determinean optimal bandwidth and spacing of the frequency response data for a range of F and a.10 12Controller parameters, NChapter 5. Simulation examples 65CMDU)eq-a-80Figure 5.31: The tradeoff curve between the objectives and By adjusting the good and badweights, the designer trades off a reduction in the RMS value of the output with an increase in the RMSvalue of the actuator signal due to noise at in and d.Figure 5.33 shows the value of the error bound times IIQII°o 1.29 versus the number of modelparameters F for a = 200, 500, 1000, and 2000. The system is guaranteed to be stable when the weightedmodel error is less than one. Hence, stability is guaranteed with F = 14 model parameters and a = 200data points.Note that for small F the error bound is almost independent of a and is dominated by the “truncationerror” term in (4.54). On the other hand, when 20 < F < 50 the model error is approximatelyequal to the “spline error” term max{!f), supsE[2,1l #(A tan(6/2))} and is less sensitive to increasesin model parameters F.The actual error of an identified model denoted is obtained with “noise-free” frequency responsedata; its weighted model error is the solid line in Figure 5.33. Although in comparison to the actual modelerrors the bounds are conservative, they provide the designer confidence the controller will stabilize the0 10 20 30 40 50 60 70RMS(y)2plant.Chapter 5. Simulation examples 66‘5Va2,m20K0(s) =1 — P,(v, s)Q(x, s)100Frequency (radis)(a)100No. of Samples(b)0C0)‘UFigure 5.32: (a) shows the upper bound on the “roll-off” rate of the plant (dash) and the Bodemagnitude plot of F, (b) shows the physical interpretation of the information M = .7 and p = 1.15. Theimpulse response of Pb! is the solid line.Closed-loop performanceThe identified model and are used to implement the controller(5.64)The parameters of and are in Table 5.2 and the step response is shown in Figure 5.34 withthe open-loop step response of the plant for comparison. The Bode plots of open-loop gain of both theanalytic and tuned solutions with gain and phase margin labeled are shown in Figure 5.35. Although themagnitude plots of both design are indistinguishable (with acceptable gain margins of 15dB), the phaseplots differ by 180 degrees at DC because of an error in the DC magnitude of F. The consequence isthat the tuned system is conditionally stable.It is evident that the controller obtained from the two-step tuning procedure meets the design specifications and yields a good approximation to the optimal analytic solution of the design problem. Indeed,the value of the performance objective for the controller (5.64) is 23% away from its “good” value. Hence,this example underlines one of the main advantages of the tuning method. Because global optimizationis used to determine both P and Q we are guaranteed to find “good” solutions to the design problem.This guarantee cannot be made for methods based on local optimization such as PID tuning (see [SDL91]Error boundsn=200—-•----.. n=10001 I I I I I0 5 10 15 20 25 30 35 40 45 50Plant model parameters FFigure 5.33: The worst-case model error bound times IQ IIo versus the number of model parametersF (dash). The system is stable for values of error less than one (dash). For all values of n the system isguaranteed to be stable for more than 14 model parameters. Actual error times IIQII of identifiedplant model (solid).67I I I I I I IChapter 5. Simulation examples1011000i5-l1000102Actual m errorSystem is stablefor example) where the designer can never be confident that the global minimum has been found.Chapter 5. Simulation examples 68Table 5.2: Parameters of the controller K0 in (5.64).2; V1.10 0.1120.707 -0.3860.699 0.4400.850 -0.07540.940 -0.1480.936 -0.03270.848 0.1630.695 0.00620.508 -0.1080,315 -0.00420.142 0.09240.0043 0.0099-0.0900 -0.0602-0.142 0.0028-0.157 0.0584-0.146 0.0102-0.120 -0.0297-0.0877 0.0097-0.0541 0.041-0.0185 -0.00171Chapter 5. Simulation examples 69— I Istep responseI’ Ig ‘I II VI0.8 ,AjTuned step response0.6E<I0.40.2CI0 5 10 15 20 25 30Time (secs)Figure 5.34: The output step response of the tuned system which meets the design specifications (solid)and the open-loop step response of the plant (dash).Chapter 5. Simulation examples15dB100L,9flfl I I I I Iio 1o3 102 101 100 101 102Frequency (radls)(b)102Figure 5.35: (a) and (b) show the Bode magnitude and phase plots respectively of the open-loop gainfor the exact (solid) and tuned (dash) solutions. Note that whereas the magnitude plots are indistinguishable, the phase plots differ by 180 degrees at DC. Both controllers have gain margins of about 15dB satisfying the design specification.70(a)0VCCUa0)ci,VcicciCU010Frequency (radls)Chapter 5. Simulation examples 715.2 Example 2: Unstable plant, the ACC’ benchmark problemThe following is an example of the control design step of the tuning method for an unstable plant. Thedesign problem is the ACC benchmark robust control design problem described in [BW9O].5.2.1 Statement of the design problemThe plant considered is shown in Figure 5.36 which models an uncertain system with a non-collocatedsensor and actuator in which the controller exerts a force on Mass 1 and the output to be controlled isthe position of Mass 2.The state-space form of the plant isii:12:13:14= [o 1 0 ojx,1 0 x1o i x2o o X3o o :34where the values of the parameters k, m1, and in2 are 1.25 kg/s, 1 kg, and 1 kg respectively.The design problem is to obtain a controller:• such that the closed-loop system is stable for plants with.5 <k 2, (5.65)Ux2=ywFigure 5.36: The benchmark problem.0000—k/mi k/mik/in2 —k/in20 0o o+ ttp+ W1/mi 00 1/in2‘ACC stands for the American Control Conference.Chapter 5. Simulation examples 72and m1 = 1 and in2 = 1,• for w = unit impulse, the system output y has a settling-time of about 15 seconds for the actualsystem in1 = rn2 = 1, and k = 1.25,• which gives reasonable system performance/robustness,• which minimizes the controller effort,• which provides robustness to the measurement noise at in characterized realistically by the designer.As with the first example, some of the specifications are given exactly whereas others are vague andleft to the interpretation of the designer. The particular interpretation chosen by the designer determinesthe level of realism of the problem.Nominal controllerThe two-step design procedure described in Section 2.6 is used since the plant is unstable and stronglystabilizable having a single unstable zero at oo [DFT9O]. Throughout the thesis it is assumed that thetransfer function of P is unknown hence we must assume the nominal controller K0 is given. The nominalcontroller is obtained using the pole placement method on the actual plant, and the set of eigenvaluesof the A matrix of H=(see (2.21)) is taken as {—1 ± j, —2 ± j}. The resulting controller is— —156s + 105s2 + —288s —80Ko(s)—4 + 12s + 63.5s2 + 186s + 302’which happens to be stable.5.2.2 Performance specificationsThis section describes the performance functionals used in the design.Robustness to variation of plant parametersSince the stability of the system is of utmost importance to the designer, it is dealt with first. Thevariation of the of the plant with the spring constant k can represented in perturbation feedback form[BB91]. That is, by a nominal plant perturbed by an internal feedback transfer function A. The blockdiagram of the perturbed plant ppett is shown in Figure 5.37 and the set of plants{pPert(A)IIIAII 1}, (5.66)Chapter 5. Simulation examples 73describes the variation of the spring constant. By expressing the uncertainty in perturbation feedbackform it is possible to obtain a sufficient condition for robust stability in terms of a bound on the Rnorm of the closed-loop transfer function from u1 to Yl (see Figure 5.37) using the Small Gain Conditiondescribed in [BB91].1up IyI IFigure 5.37: The variation in the spring constant .5 < Ic < 2 can be represented as varying the feedbackperturbation A inside the perturbed plant pPt•The closed loop transfer function from it1 to Pi isW1 + W2Hyr,whereWi(s)= s2+25and—.8s2W2(s)=52+25Hence, the performance specification for robust stability iskHinf(X) = .7511W1 + W2Hyr(x)IIco 1.Since stability is a hard constraint we define+I II I=Chapter 5. Simulation examples 74following the convention in Section 2.3.4 of using the letter b for objectives associated with hard constraints.Impulse disturbance rejectionThe definition of the settling-time of an impulse response is taken asinf{rI I(OI .01 ,Vt r).The easiest way to ensure a settling-time of 15 seconds for the impulse response from w to y is to usethe step response functional (3.22) with an impulse as opposed to step input. Hence, we defineifr)= tcWmax{h(x, 1) — a(t), /3(t) — h(x, t), 0}, (5.67)where h is the impulse response of H5, and the upper and lower bounds a and /3 are shown in Figure 5.43(a). These bounds also constrain the magnitude of the response for I < 15 seconds.Note that the procedures given in Chapter 3 for evaluating performance functionals and their subgradients are for stable plants. To illustrate how measurements are used to compute the impulse responseof II for an unstable plant, consider the system in Figure 5.38. First a force with large amplitude andshort duration and “unit-area” (see Section 3.3.1) is applied to Mass 2 (at input w) and its responsemeasured at u3. Next, the measured response is applied to u1 and the response at u3 (the positionof Mass 2) is measured. Finally, the measured outputs are subtracted to obtain the impulse responseof Hence, the procedures for obtaining closed-loop responses from an open-loop system given inChapter 3 differ only slightly to that used for the unstable plant.Tracking performancePerformance of the system is in part specified by constraining the output step response using the stepresponse functional (3.22) which is denoted as 4’2. The bounding functions a and /3 (which are not thesame as those in (5.67)) are shown in Figure 5.43(b). If these bounds are not violated the step responsehas a settling-time of 20 seconds, a “steady-state” tracking error of .1%, an overshoot of less than 1.5m, and an undershoot less than 1 m.ControUer effortIt is important that a realistic design should limit the size of the actuator signal. The approach takenhere is measure the size of the actuator effort with a weighted-sum of RMS gains of the closed-loopChapter 5. Simulation examples 75transfer functions with output u. That is43(X) = IIHurIIoo + IIHudIIc,o + (5.68)Note that IIHumIIoo is not included in (5.68) since it equals II1tyrIIoo.An alternative approach to measuring the size of the actuator signal using (5.68), is to use the timedomain functionals in Sections 3.2 and 3.3 which measure norms of particular responses.Rejection of measurement noiseThe RMS value of the system output due to white measurement noise is measured by ç5, an “H2 norm”functional discussed in Section 3.4. That is4)4(X) = IIHymII2.Specification on Gain marginThe specification of a negative gain margin equal to 20 log LdB and a positive gain margin equal to20 log UdB, where 0 < L 1 < U, can be satisfied by ensuring robust stability for the perturbed plantsetUIU2w+Figure 5.38: The system used to measure the impulse response of Hd. The response at u3 is measuredto an impulse signal applied at w. This response is then applied to u1 and the output at u measured.{GIG=cP,L c < U} (5.69)Chapter 5. Simulation examples 76The set (5.69) can be represented in perturbation feedback form shown in Figure 5.39 [BB91]. That is,the perturbed plant set (5.69) is a subset of{pPert(A)IIIII, < max{1— L, U — 1}}.Using Figure 5.39, the feedback seen by the perturbation is Hence, applying the Small GainCondition, the inequality.11 145(x)= I[yrIIco <min1— L’ U—iensures a positive gain margin of 20 log UdB and a negative gain margin of 20 log LdB.p pee ()+ ,_______up IyCDI+ IFigure 5.39: The perturbation feedback form used for the gain margin specification.5.2.3 Results of the control design stepThe basis functions (2.9) with p = 1 and N = 20 are chosen. After several design iterations in whichthe good and bad weights of associated with the objectives were varied, the performance specificationsare determined by the inequalities‘1 .01, 42 .01, < 50, < 1.5,qf 1. (5.70)The choice of good and bad weights used to obtain the optimal solution are in Table 5.3 along withthe raw and normalized optimal values of the objectives. Note that although the specifications (5.70)Chapter 5. Simulation examples 77are not obtained, performance is deemed satisfactory. This illustrates the importance of specifyingperformance as soft constraints. Due to the vagueness in the statement of the problem, small violationsof the inequalities (5.70) are permissible since they do not affect the overall quality of the design. On theother hand, because an unstable design is useless, hard constraints are necessary to ensure the robuststability specification is met.Table 5.4 gives the parameters of Q. The Bode magnitude plots of the closed-loop system are inshown Figures 5.40 and 5.41, the open-loop gain is shown in Figure 5.42, and the step responses of u andy for a unit step applied at the reference input are shown in Figures 5.43(a-c) along with the impulseresponse of All figures show the effect of varying the spring constant k; the dashed curves are fork = .5 and the dotted curves for k = 2. Although the system is stable for all plants in the perturbedset (5.66), a small stability margin is observed for k = .5. That is, the peak magnitude of the frequencyresponse of Her is about 20 (see Figure 5.40(b)).Note that in Figure 5.43(a-b) the impulse response meets the settling-time specification and that theoutput step response lies within the desired bounds. Additionally, the actuator output step response hasa peak value of 70 N, and the (positive) gain and phase margins are 3.3dB and 25.6 degrees respectively.As well, note that in Table 5.3 ç is the only functional which is not competing at the optimal solution.Hence, if desired the designer could further reduce the effect of measurement noise by decreasing G4 orB4. This example shows clearly how the general form of the cost function allows the designer to dealdirectly with a wide range of performance specifications. This is compared with other tuning methodswhich minimize some crude measure of performance.Table 5.3: Optimal values of the objectives and the choice of good and bad weights.Gk Bk Raw value Normalizedvaluei 4/3 8/3 1.33 0.00ci .01 .02 .0227 1.272 .01 .02 .0227 1.2750 70 75.3 1.271.5 2 1.32 -.3251 2 2.27 1.27One of the features of the control design step is that it gives the designer an easy way to investigatethe relationships between the objectives by adjusting the normalizing weights in the weighted-max costfunction. For example, suppose the designer is unsatisfied with the 3.3dB positive gain margin obtainedChapter 5. Simulation examples 78Table 5.4: The optimal controller parameters.Xl_lo Xfl_IJ-.775 6.94.164 4.561.06 2.493.46 1.006.13 .1218.79 -.25910.7 -.30311.4 -.22210.9 -.1039.23 -0.0362for the optimal design and wants to know how it can be increased by increasing the controller effort. Byadjusting the weights G3, B3, G5, and B5 and constraining the values of the objectives , , and 44to their “good” values given in Table 5.3, the designer explores the tradeoff curve shown in Figure 5.44.This figure shows that a positive gain margin can be guaranteed greater than 4.23dB if the designeraccepts a value of 325 for (which is roughly equivalent to an controller output step response with apeak value of 300 N).As mentioned previously, a central problem with performing a design is specifying the performancespecifications quantitatively. The tradeoff procedure described above illustrates clearly how the tuningmethod simplifies this problem. By adjusting good and bad weights which have obvious physical interpretation the designer is able to reassess his or her notion of optimal performance throughout the designprocess. For the above example, the designer has no idea of what a realistic specification on the gainmargin is to ensure a small enough actuator effort. Additionally, even when exact design specificationsare known, in many cases because they are unachievable and some best compromise must be obtained.Such a compromise is obtained by performing tradeoff analyses such as the one described above.Finally, we reemphasize that tuning is performed on-line optimizing actual closed-loop performance.This is opposed to an off-line design which is inevitably corrupted by modeling errors. Figure 5.45illustrates typical convergence observed by the designer during the procedure. The step response of Hg,.is “pulled” between upper and lower bounding functions associated with the step response functional‘p2.Chapter 5. Simulation examples 79100 101Frequency (racils)(b)(a)1010 .1 • Ii::10.20 I I I I102 io1 102 io3102Frequency (radls)1Figure 5.40: (a) shows the Bode magnitude plot of the complimentary sensitivity function. (b) showsthe Bode magnitude plot of the sensitivity function. The effect of the variation of the spring constantis shown by the dashed curves for k = .5 and the dotted curves for k = 2. The small stability margin ofthe design for k = .5 is indicated by the large peak value of dashed curve in (b).Chapter 5. Simulation examples 80(a)1010 I . .... Ia)I::: II I102 10_i 10° 101 102 3Frequency (rad/s)(b)1020 Ia)1020102 10i 100 101 102 to3Frequency (rad/s)(c)10 I I I I.jj2 i0Frequency (rad/s)Figure 5.41: (a) shows the Bode magnitude plot of (b) shows the Bode magnitude plot of Hdy.(c) shows the Bode magnitude plot of Hur. The effect of the variation of the spring constant is shownby the dashed curves for k = .5 and the dotted curves for k = 2.Chapter 5. Simulation examples 81VCCu0a)ci0)0)Cu0(a)Figure 5.42: (a) and (b) show the Bode magnitude and phase plots respectively of the open-loop gainof the system. For k = 1.25 the systems’ positive gain margin and phase margin are 3.30dB (w = .565rad/s) and 25.9 degrees (w = .292 rad/s). The effect of the variation of the spring constant is shown bythe dashed curves for k = .5 and the dotted curves for k = 2.Frequency (radfs)(b)102 10.1 10° 101 102 ioFrequency (radis)10 15 20 25 30 35 40Time (secs)(b)Figure 5.43: (a) shows the output step response and bounding functions associated with the step responsefunctional 42, (b) shows the impulse response of and bounding functions associated with qfi, (c)shows the actuator effort for a unit step applied at r. Note that the specifications of asymptotic trackingand and settling-time are satisfied. The effect of the variation of the spring constant is shown by thedashed curves for k = .5 and the dotted curves for Ic = 2.Chapter 5. Simulation examples 82(a)2aU)zci20LITime (secs)(c)5Time (secs)Chapter 5. Simulation examples 833.232.82.62.42.221.81.614 I I I I0 50 100 150 200 250 300 35003Figure 5.44: The tradeoff curve between controller effort measured by 4 and a measure of the systems’gain margin measured by.The objectives , qf2, and qo are constrained to be less than or equalto their “good” values in Table 5.3. Increases in positive gain margin can be traded off with increasedactuator effort.G3=15, G5=4.5, B3=25, B5=4.6G =55, G =2.5, B =65, B =2.63 5 3 5U,G3=75, G5=1.5, B3=85, B5=1.6Chapter 5. Simulation examplesa0a0a0aCl)Figure 5.45: (a)-(d) show the typical convergence observed for the tuning method. An output stepresponse converges to within bounding functions (associated with q2) for the first five iterations of theoptimization algorithm.84(a)I I I I I I I I Io--•-5 I I I I I I I Io 5 10 15 20 25 30 35 40 45 50Time (secs)(b)20J7”—”’•1 I I I I—2 I I I I I I0 5 10 15 20 25 30 35 40 45 50Time (secs)(c)•I I I I II I I I I I I-20 5 10 15 20 25 30 35 40 45 50Time (secs)(d)II I-2 I I I I I I I I0 30 35 40 45 505 10 15 20 25Time (secs)Chapter 6Conclusions/Future work6.1 ConclusionsIn conclusion, this thesis has presented a new optimization-based on-line tuning method for linear compensators with “internal model” structure, The method comprises of a control design step for thecontroller parameter Q and an identification step to obtain the internal model of the plant. This chaptersummarizes the key points of the method and discusses possible future work.6.1.1 Control design stepA standard approach to designing a feedback compensator is to first obtain a model of the plant andthen to apply one of the well-established design techniques discussed in Section 1.2.1. The difficultywith this approach is that the achievable performance depends on the quality of the model used. Forexample, in many cases models ignore higher order dynamics because they are impossible to model orbecause they make the design process too complex. In an effort to overcome this limitation, the tuningmethod proposed uses plant measurements to perform the control design step.Another fundamental difficulty with most design methods is that they require that real world engineering specifications be transcribed into some standard form. For example, recall that the LQG methodrequired that the designer select weighting matrices in the cost function and spectral factors of fictitiousnoise inputs. In most cases these design parameters cannot be related to the given specifications andtheir choice is difficult. This difficulty is addressed with parameter optimization design methods whichallow a wide range of performance specifications. Furthermore, with recent advances in optimizationtechniques an even wider selection of measures of performance is possible. The tuning method presentedexploits this advantage of parameter optimization-based design methods and allows the designer to specify desired performance with inequalities involving H, H2, L1 norm type functionals, and a range oftime domain objectives.Of course, even empowered with a general cost function, it is still difficult sometimes to specify85Chapter 6. Conclusions/Future work 86desired performance precisely (it may not even be known). This is because for many realistic designproblems desired performance cannot be quantified and the designer must experiment with a range ofspecifications. The tuning method addresses this problem by proposing an iterative design procedurewherein the designer adjusts the weights associated with a weighted-max cost function. The weightshave obvious physical interpretation and the designer is comfortable in adjusting their values. Thisis contrasted with the iterative procedure based on other design methods (such as LQG), where the“tuning” weights have little physical interpretation.The tuning method uses an “internal model” controller structure to formulate the design problem asa convex optimization problem. This implies that, assuming the performance functionals measured areexact, an arbitrarily “good” approximation to the optimal solution is obtained by solving a sequenceof convex optimization problems with increasing number of controller parameters. This is opposed totuning with other controller structures (such as that of the commonly used PID controller) where thedesigner is never certain whether performance can be improved. Although the price paid for convexityis a higher order compensator, due to its simple structure this is not believed to be a major drawback.A problem with using parameter optimization-based design methods is that they sometimes requirethe solution of computationally costly optimization problems. The tuning method addresses this problemby using plant measurements to evaluate the performance functionals. A comparison of the cost ofcomputing the performance functionals in the tuning method shows that it is almost always advantageousto use measurements as opposed to a plant model. Hence, even if a good plant model is available thereis still reason to apply the method on-line to reduce the amount of computation.Finally, for the case where the plant is unstable it is possible to apply the control design step provideda stable stabilizing controller can be found.6.1.2 Identification stepThe second step of the tuning method identifies a plant model required by the internal model controllerstructure.Recently much attention is being paid to identifying a plant for the purpose of controller design.Most standard identification methods emphasize obtaining the nominal model but pay less attentionto the description of its uncertainty. On the other hand, most modern control design methods rely ona description of the plant’s uncertainty. For example, the H-optimal design method design methodassumes the plant lies in some H,,, frequency domain bound. This thesis presents a control-orientedChapter 6. Conclusions/Future work 87numerical optimization-based identification method that depends on desired closed-loop performanceand provides a description of the uncertainty in plant which is compatible with robust control methodsand the control design step of the tuning method.The identified model is obtained by minimizing the 11cc norm of a model error weighted by theoptimal design parameter obtained in the control design step. This weighting ensures that the desiredclosed-loop performance is reflected in the plant model. Since an FIR model structure is used, theoptimization problem solved is convex and as in the control design step the designer is guaranteed aglobal solution to the design problem.In addition to providing a nominal model, the identification method provides a worst-case bound onthe model error. This bound is computed from assumed information about the plant, the experimentduration, the noise level, and the model order. Such a bound allows performance of the first design stepto be guaranteed.6.2 Future workThis section presents possible future work that could be performed on the tuning method.• When performing the control design step, the data required was assumed free from noise. Thisis not a realistic assumption in practice. An analysis of the effect of noise (realistically modeled)should be performed to determine the sensitivity properties of the tuning method.• The effectiveness of the control design step relies on ability of the designer to interact with thecomputer. Hence it is important that a user-friendly interface to the existing software be developed.This software should allow the designer to adjust the good and bad weights during optimizationand should provide a graphical summary of the performance functionals. As well, work should beperformed on a compiler which allows the designer to input performance functionals in convenientway. This idea is similar to the QDES software proposed in [BBB 88] for off-line controller design.• An adaptive control scheme based on the parameter optimization-based control and identificationshould be developed. The general form of the cost function would allow more realistic measures ofperformance to be used (compared to pole placement or LQR design methods used in traditionaladaptive control techniques).Chapter 6. Conclusions/Future work 88• One shortcoming of the identification step is the method used to obtain the frequency responsedata and the bound on its error (Appendix F). A popular way of computing the frequency responseis to use the “tn-spectrum” average method described in Appendix C.4.2. It would therefore bedesirable to obtain a hard bound (i.e. an M) on the error of the frequency response data identifiedusing this method.• When using the Ritz approximation for identification it is possible to use information about theplant and noise level to derive a worst-case expression for the error for the model (This was notpossible when designing Q since in most cases the designer has little information about the optimalsolution.) Although an FIR plant model structure was used, it is well known that the choice ofbasis functions affects the number of parameters required to obtain a given error in the nominalmodel (see [Wah9lj and the references therein). Hence in an effort to reduce the number of modelparameters, expressions for the worst-case model error using other popular basis functions wouldbe useful.• The tuning strategy should be applied to an actual system to fully demonstrate its usefulness.Appendix ASome rules for computing subgradientsThe ellipsoid algorithm is a simple optimization algorithm which requires only the value of the objectiveand a subgradient at each iteration. A subgradient of a convex function exists even where it is nondifferentiable. This appendix defines the subgradient of a convex function and specifies rules determiningits expression. An elementary treatment of convex analysis is given in [BB91] and a more thoroughreference is {C1a83].A.1 Basic definitionsA functional 4 : RN —* R is convex if for all x1, x2 in RN and A e [0, 1] then#Ptxi + (1 — A)x2) <A#(xi) + (1 — A)c6(x2).A subgradient of a convex functional is defined as any g satisfying(x) (xo) + gTfr — xo),Vx.An example of convex functional 4 is shown in Figure A.46. Although 4 is not differentiable at x0,subgradients can be computed. The dotted tangent lines correspond to two subgradients gi and g.A.2 Rules for computing subgradientsThe following is a set of rules for computing a subgradient of certain general functionals. These rulesare applied in Chapter 3 to the performance functionals specifying the control system performance.• If the functional is differentiable at x, then its subgradient is its derivative at z.• If w 0 and q5 is convex then a subgradient of wq5 at x is wg where g is a subgradient of at x.• If #(x) = çbkfr) and j, i = 1, . . . N are convex then =1Ok, where g is the subgradientof k at x, is a subgradient of 4’ at x. Note that the peak gain norm functional (the L1 norm89Appendix A. Some rules for computing subgradients 90slope g2slope g1Figure A.46: An example of a convex functional which is not differentiable at x0, The dotted tangentlines correspond to subgradients g and g.of the system’s impulse response) (3.38), can be interpreted as in integral of a family of convexfunctionals and the summation can be replaced by an integral in the above rule to obtain a wiseguess for the subgradient (3.39).• Take q’ = sup{qj(x) i I}, where I is any index set and çb is convex. A subgradient of 4j0(x),where i0 e I such that 4j0(x) = q(x), is a subgradient of q at x.• Take = max{(x) Ii e I}, where I is any finite index set and qj is convex. A subgradient of4j(z), where i0 E I such that qj(x) = q(x), is a subgradient of at x.4(x)—— VV——VxoAppendix BA two-step design for strongly stabilizable plantsThe IMC controller structure can be extended to strongly stabilizable unstable plants using a two-stepdesign procedure. The first step is to design any stable stabilizing controller K0, and the second stepis to modify the dynamics of the stabilized controller using the IMC structure. The expression for thecontroller isK(F,Q) = K0 + Q(1 — H(P)Q),Q€ RH0,, (B.71)whereH(F) = F(1 +K0F)1 (B.72)(assuming P = P). This appendix shows that as Q varies over RH, in (B.71) we parametrize the setof all stabilizing controllers. The explanation is due to [Mac89j.LetF = NM1, (B.73)with N, M coprime and K0 = NK0Mj be a stable controller that stabilizes F with Nç0,MK0 coprime.ThenX = NV’ (B.74)Y = MK0V1, (B.75)where V = NNK0 + MMK0,satisfyMY+NX=1.Note that because K0 was assumed stable, Y is invertible in RH0,.To show that the two-step design is over all stabilizing controllers we show (B.71) is of the formK(P,Q) = (X + MQ)(Y - NQ)1, (B.76)where X and Y are defined in (B.74) and (B.75) respectively; N and M are defined in (B.73).K(F, Q) = (K0 + MQY1)(1 -91Appendix B. A two-step design for strongly stabilizable plants 92= (K0 + MQY’)(1 + NQY1(1 — NQY1)= K0 +K0NQY’(l - NQY1)+ MQY’(1 - NQY1)’= K0 + (K0NQ + MQ)(1 - NQY)Y= K0 + (XY1N+ M)Q(1 -= K0 + (XN + MY)Y’Q(1 -YQN)’Y= K0 +Y1Q(1 -1QN)’Y= Ko+R(1—YRN). (B.77)Where R is defined as R Y2Q in (B.77). Note that because K0 is stable so is R. Since (B.72) canbe rewritten asH(P) = NM1(1+NM’XY’)= N(M + NXY’)’= NY’(MY+NX)= NY1,and (B.71) can be expressed asK = K0 + Q(1 — NY’Q)1. (B.’78)Comparing (B.77) and (B.78) and since there is a one-to-one relationship between R and Q, (B.71) isequivalent to (B.76) and the two-step design procedure is over all stabilizing controllers.Appendix CCost of evaluating the performance functionalsOne of the advantages of using the tuning method is that it avoids the calculation of time and frequencyresponses from the plant model:1. These can be costly when the order of the system is high, or when many time or frequency pointsof the response are required,2. They can be inaccurate due to the approximations made when obtaining the plant model.In Chapter 3 a comparison of the cost of evaluating performance functionals using a model and usingmeasurements from the plant was made. This Appendix presents the details of how the comparison wasmade.For simplicity the set of plant models considered isF(s)= (i)k’ k = 1,...,20, (C.79)and the closed-loop transfer function is 11y = FQ. For each functional two costs are considered: thecost before optimization begins, and the cost at each iteration of the optimization algorithm.C.1 The step response functionalThis section considers the step response functional=max max{s(x, 1) — a(t), /3(1) — s(x, 1), 0). (C.80)C.1.1 Evaluating the functional using a plant modelThe easiest way to evaluate (C.80) is to compute the step response s(x, 1), 0 < I < tmar• To avoidcomputing s at each iteration of the optimization algorithm we determine its expression in term of thecontroller parameter x. Letting p, bk and qN denote the impulse responses of Bk, QN, and F respectively93Appendix C. Cost of evaluating the performance functionals 94thens(x,t)= jp(r)*qN(x,r)dr= J j p(r cr)qN (x, cr)ddr= xkjjp(T—)bk(o)dodr=where the kth component of c (an N by 1 vector) is defined asfOGck(t) = I I p(r — 0)qk(0)du. (C,81)Jo JoThe value of c at N5 evenly spaced points on the interval [0, tmaj is evaluated before optimizationbegins, generating an N by N3 matrix. This is equivalent to evaluating the step response of PBk Ntimes at N8 points, which costs roughly N((Na + i)(6 + j + 1/3) + N(2N + 3N0 + 1)) floatingpoint operations (flops). The dimension of A, the state-space matrix of PBk, is N0, and j = max{0, 1 +floor(log2(AIOG))}. For example, with N0 = 40, N = 20, and hAil00 = 20 the cost is approximately 100Mflops. (See Appendix D for how the expression for the cost of a time domain response is obtailled.)The amount of computation obtained by simulation for evaluating ck, k = 1, . . . ,20 at N3 = 1000time points versus the number of states in the plant model (C.79) is shown in Figure C.47. Note that therough estimate obtained above for the cost of evaluating c agrees with the value obtained by simulation(128 Mflops at 20 model states in Figure C.47).Appendix C. Cost of evaluating the performance fun ctionals 95Cl)C0Cua)0C00CCu0L18 10 12 14Number of states in plant modelFigure C.47: The number of flops versus the number of states in the plant model to evaluate ck,k = 1,. . . ,20 in (C.81) at N3 = 1000 time points.7x 10Appendix C. Cost of evaluating the performance functionals 96The MATLAB program used to generate Figure C,47 is:clearflops_step_accum=[j;% Plant transfer function denominator.denom=l;% 1000 time pts.t=(0:999)*.Ol;% State space form of B_i.N=20;p=l; % p is arbitrary.Aq=-ones (N) *2*p;Aq=tril (Aq, 0) +p*eye(N);Bq=ones(N,l);Cq(l:N)=_2*p*ones(l,N);Dq=1;for k=1:20,% Plant P_k=l/(l+s)’k.denom=conv(denom, [1 11);[Ap,Bp,Cp,Dpj=tf2ss(l,denom);% PB_i.[Apq,Bpq,Cpq,Dpqj=series(Ap,Bp,Cp,Dp,Aq,Bq,Cq,Dq);flops(0)% Compute step response of PB_i and count the number of operations.s=step (Apq, Bpq, Cpq, Dpq, l,t);% Multiply by 20 to account for i=1,...,20.flops_step_accuin [flops_step_accum flops*20];endplot(l:20, flops_step_accum)xlabel(’Number of states in plant model’)ylabel ( ‘Floating point operations’)During optimization the evaluation of (C.80) reciuires (2N — 1)NL3 flops for evaluating cTx and2N5 flops to add bounding functions a and 9.Appendix C. Cost of evaluating the performance functionals 97C.1.2 Evaluating the functional using measurementsWhen using measurements the step response is measured directly from the open-loop system, and nocomputation is required before optimization begins. The only computation is to add the boundingfunctions which requires 2N3 flops. Hence, the relative cost (during optimization) of using a modelversus using measurements isRelative Cost = (2N + 1)N3 = 2N + 1> 1.2N8The relative cost of evaluating (C.80) during optimization increases linearly with N, and computationalsavings by using measurements are possible for all N.C.2 Norms of a particular response functionalConsider the RMS value functiollalRMS(X) = (lirn j y(x, t)2d1), (C.82)where y is the output response to a square wave input at r (with m = d = w = 0).C.2.1 Evaluating the functional using a plant modelAs with the step response functional (C.80), the easiest way to evaluate (C.82) is to first compute theresponse y and then use the expression for the functional directly. To determine the cost of evaluationbefore optimization begins, the expression for the response in terms of the controller parameter isobtained. Denoting the input u, the output response y in terms of x isy(x,t)=(p(t) * qN(xt)) * u(t)= j (jp( — r)qN(z, r)dr) u(t —Xk j (fp( — r)bk(T)dr) u(t —= dT(t)x, (C.83)where the kth component of d is defined asdk(t) = j (j p( — T)bk(T)dT) u(t —Appendix C. Cost of evaluating the performance functionals 98The value of d is evaluated at N3 evenly spaced points on the interval [0, tmavj , where tma is chosento be the time at which the system is in approximately in “steady-state”. The cost of evaluating d isroughly the same as for c in Section C.1.1 since both require computation of time domain responses.Figure C.48 shows the cost of computing d versus the number of states in the plant model.Cl)C0Cua)0C0Q.CCu0U-Figure C.48: The number of flopsof d in (C.83) at 1000 time points.versus the number of states in the plant model for the computation2 4 6 8 10 12 14Number of states in plant model16 18 20Appendix C. Cost of evaluating the performance functionals 99The MATLAB program used to generate Figure C.48 is:clearflops_part_accum= [1;% Plant transfer function denominator.denom=l;% 1000 time pts.t=(0:999) *01;% State space form of B_i.N=20;p=l; % p is arbitraryAq=-ones (N) *2*p;Aq=tril(Aq,0)÷p*eye(N);Bq=ones (N, 1);Cq(1:N)=_2*p*ones(l,N)Dq=1;% Input is a square wave.fs=l;u=square(t*2*pi*fs);for k=1:20,% Plant P_k=1/(1+s)”k.denom=conv(denom, [1 1]);[Ap,Bp,Cp,Dp]=tf2ss(1,denom);% PB_i.[Apq,Bpq,Cpq,Dpq]=series(Ap,Bp,Cp,Dp,Aq,Bq,Cq,Dq);flops (0)% Compute response of PB_i and count the number of operations.s=lsim(Apq,Bpq,Cpq,Dpq,u,t);% Multiply by 20 to account for i=1, . . . ,20.flops_part_accum = [flops_part_accum flops*20];endplot (1:20, flops_part_accum)xlabel(’Number of states in plant model’)ylabel ( ‘Floating point operations’)At each iteration of the optimization algorithm, the evaluation of (C.82) requires (2N — 1)N flopsfor evaluating dTx, and 2N3 + 1 flops for the integration.Appendix C. Cost of evaluating the performance function als 100C.2.2 Evaluating the functional using measurementsSince the cost of evaluating the functional from measurements is the cost of performing the integration,the relative cost using the plant model compared to using measurements (during optimization) is(2N + 1)N3 + 1Relative Cost = > 1.2N3 + 1The relative cost during optimization increases linearly with N and is greater than one for all N.C.3 The “peak gain norm” functionalConsider the “peak gain norm” functionalpkgn(X)= j (x, r)Idr, (C.84)where Ii is the impulse response of Hyr.C.3.1 Evaluating the functional using a plant modelAs with the step response functional and the norm of a particular response functional, to avoid computingthe impulse response at each iteration its expression in terms of x is obtained. Letting p and q’ denotethe impulse responses of P and Q’ respectively, the impulse response of PQ in term of xp(t)*qN(x,t) = Jp(r)q”1(x,t_r)dr= p(r)bk(t—r)drk=1 0= gT(t)x (C.85)where the kth component of g is defined asgk(t)=— r)dr. (C.86)The integral in (C.86) is approximated by summing the magnitude of the impulse response at N3evenly spaced points on the interval [0, Imax], where tmax is the time at which the impulse response isapproximately zero. Using N3 = 1000, the cost of computing gy, Ic = 1, . . . , 20 is shown in Figure C.49.Appendix C. Cost of evaluating the performance functionals 101C’,0(50.000.0>CC’0U-8 10 12 14 16 18 20Number of states in plant modelFigure C.49: The number of flops versus the number of states in the plant model for the computationof g, k = 1, . . . ,20 in (C.86) at 1000 time points.7x 102 4 6Appendix C. Cost of evaluating the performance function als 102The MATLAB program that produced Figure C.49 is:clearflops_impulse_accum= [1;% Plant transfer function denominator.denomDl;% 1000 time pts.t(0:999) *01;% State space form of B_i.N=20;p=l; % p is arbitraryAq=-ones (N) *2*p;Aq=tril (Aq,0)+p*eye(N);Bqtones (N, 1);Cq(l:N)=_2*p*ones(l,N);Dqzl;for k=l:20,% Plant P_k=l/(l+sV’k.denom=conv(denom, dt);[Ap,Bp,Cp,DpJ=tf2ss(numer,denorn);PB_i.[Apq,Bpq,Cpq,DpqJ=series(Ap,Bp,Cp,Dp,Aq,Bq,Cq,Dq);flops(0)% Compute impulse response of PB_k and count the number of operations.s=impulse(Apq,Bpq,Cpq,Dpq,l,t);flops_impulse_accum = [flops_impulse_accum flops*201;endplot (1:20, flops_impulse_accum)xlabel(’Nuinber of states in plant model’)ylabel(’Floating point operations’)During optimization, the cost of evaluating (C.84) is 2NN3 — N3 — 1 flops. This is because themultiplication in (C.85) requires 2N8(N — 1) flops, and the integration requires — 1 flops.C.3.2 Evaluating the functional using measurementsDuring optimization, no computation for the impulse response is required since it is obtained by measurement. Hence, the cost of evaluating the functional is N8 —1 (the cost of performing the integration).Appendix C. Cost of evaluating the performance functionals 103The relative cost is therefore2N3N—N—1Relative Cost = > 1,2N3 — 1which increases linearly with N and is greater than 1 for N 2.C.4 The “H2 norm” functionalConsider the “H2 norm” functional (3.41) with Wi(s) = 1, W2(s) 0, and ft = PQ:= ( J IP(iw)Q(x jw)I2dw). (C.87)C.4.1 Evaluating the functional using the plant modelAn efficient method exists for evaluating the “H2 norm” functional using a model, First PQN is expressedin terms of the controller parameterpQN(x)= ZkPBk= (xkGk) (sI—A)1B,where A, B and Ck are the state space matrices of PBk (D = 0 since P is strictly proper). Hence1pQN12= /xTEx, (C.88)whereC1E = Wcrb [ Gj1’ ... C j,CNand Wctrb is the controllability Gramian, the solution to the Lyapunov equationAWctrb + Wctr&AT + BBT = 0 (C.89)[BB91]. The cost of computing E in (C.88) versus the number of plant model states is shown inFigure C.50.Appendix C. Cost of evaluating the performance functionals 1048 10 12Number of states in plant modelFigure C.50: The number of flops versus the number of states in the plant model for the computationof E in (C.88).2Appendix C. Cost of evaluating the performance functionals 105The MATLAB program used to generate Figure C.50 is:clearflops_h2_accuni= [1;% Plant transfer function denominator.denom=l;% State space form of B_i.N=20;p=l; % p is arbitraryAq=-ones (N-i) *2*p;Aq=tril (Aq, 0) +p*eye (Ni);Bq=ones (N-i, 1);Cq(l:N_l)=(_2*p)*ones(i,N_l)Dq=1;for k=l:20,% Plant P_k=l/(l+s)’ .denomzconv(denom, [1 1]);[Ap,Bp,Cp,Dp]=tf2ss(l,denom);% PB_i.[Apq,Bpq,Cpq,Dpq]=series(Ap,Bp,Cp,Dp,Aq,Bq,Cq,Dq);flops (0)% Compute the controllability grarnmian.W = gram(Apq,Bpq);[N_apq,N_apq]=size(Apq);E = sqrt((ones(N_apq,l)*Cpq)*W*(ones(N_apq, l)*Cpq)I);flops_h2_accum = [flops_h2_accurn flops];endplot (1:20, flops_h2_accum)xlabel(’Nurnber of states in plant model’)ylabel (‘Floating point operations’)At each iteration of the optimization algorithm, the cost of computing (C.87) is 2N flops for performing two multiplications of a vector and a matrix, and taking the square root in (C.88).C.4.2 Evaluating the functional using measurementsThe method to evaluate (C.87) using measurements requires the frequency response of the plant whichis obtained by applying a sequence of square waves to the plant input u and measuring the output y. Ifthe number of points measured at a sampling frequency f3 is M, and the frequency of the square wave isPf3/M, where P is the number of periods of square wave applied; the estimate of the transfer functionAppendix C. Cost of evaluating the performance functionals 106of the plant is obtained by computing the “tn-spectrum” averageP(k)P(wk)= IJ.12irkf3Wk=k = 1, . . ., M/2 — 1. (C.90)Where F is the cross power spectral density, and P, P are auto power spectral densities of theinput and output respectively. That is,P(k) = U(k)12P(k) IY(k)12P(k) = Y(k)*U(k),where U and Y are the Discrete Fourier Transform of the sampled input IZ and output weighted by aHanning window respectively. That is,M-1U(k) = ii(n)w(n) exp_i(2M)Y(k)k = 0,1,...,M—1,1/ /2irnw(n)= 2lcosMlThe frequency response of Bk, obtained using the state-space method (described in Section C.5.1),is multiplied by the estimated frequency response of the plant. The cost of evaluating the frequencyresponse of FBk using measurements relative to using a plant model versus the number of states in theplant model is shown in Figure C.51 using 370 frequency points (dashed) and 92 frequency points (solid).It is evident that computational savings are possible with 370 frequency points and when there are morethan about 12 states in the model. The MATLAB programs used to generate the cost of computing thefrequency response of FBk using measurements are:Appendix C. Cost of evaluating the performance function als 107clearH=[J;f=[];flops_acc=0;% The number of periods of square wave.P=10;% Sampling frequency.fs=.256;% Number of samples.M=1024*4;% call the function ‘compute_trispectrum_average five times.for k=l:5,[Z, f_o, flps] =compute_trispectrum_average (fs, P,M);if k > 1,% Throw away overlapping frequency points.ind=find(f_o > f(length(f)));H=[H; Z(ind));f=[f f_o(ind)];elseH=Z;f=f o;endflops_acc=flops_acc÷flps;% Increase sampling frequency.fsfs*lO;end% Interpolate to obtain logarithmcally spaced points.f_inter=logspace(-4, 3,370);H_inter_mag=interp2(f,abs(H) ,f_inter);H_inter_ph=interp2(f,unwrap(angle(H) ) ,f_inter);H=H_inter_mag. *exp( ((H_inter_ph) *sqrt (-1)));flops_acc=flops;% State speace form for B_i.N=20;p=l; % p is arbitraryAq=-ones(N) *2*p;Aq=tril(Aq,0)÷p*eye(N);Bq=ones (N, 1);Cq(l:N)=_2*p*ones(l,N);Dq=l;% Compute frequency reponse of PB_i and count the number of operations.flops(0)% Compute frequency response of PB_k and count the number of operations.% Balance A[t,a] = balance(Aq);b = t \ Bq;c = Cq * t;% Reduce A to Hesenburg form[pa] = hess(a);% Apply similarity transformations from Hessenberg reduction to B and C:b = p * b;c = c *g = 1tifr(a,b,sqrt(_l)*f_inter);% Compute frequency response of PBk.for k=l:N,PB(k,l:length(f_inter)) = (c *endflops_acc=flops+flops_acc;Appendix C. Cost of evaluating the performance functionals 108function [Z, f_o, flps]=compute_trispectrum_average(fs, P,M)fsq=P*fs/M;% P periods of square wave.% Increase sampling frequency by inc.inc=4;t=0: (1/ (inc*fs)) : (((M/fs)-(l/ (fs))));% input square wave.u=. 5*square(2*pi*fsq*t);[a,b,c,d]=tf2ss(l,[l 1]);y=lsim(a,b,c,d,u,t);% decrease sampling rateu=decimate(u, mc);y=decimate(y, mc);flops(0);% Compute transfer function using trispectrum% average.p=spectrum(u,y,length(uH;[mT,nT]=size(p);Z=p(2: (mT) ,4)f_base=(l: ( (M/(2) )—l) )*(fs/M)% .8 times bw due to filtering in ‘decimate’.% take only the frequency components at which square% wave has energy.f_o=f_base(P:2*P:length(f_base) *.8)Z=Z (P: (2*P) : (length(Z) * .8));flps= flops;At each iteration of the optimization algorithm the computation of the frequency response of pQNrequires N8(2N — 1) flops, and the integral requires 2N3 + 1 flops. The relative cost is therefore2NRelative Cost = (2N + 1)Np:s + 1Since N3 is typically larger than the number of controller parameters N, the relative cost is less than1. Hence, evaluating (C.87) using a plant model is more efficient than using measurements duringoptimization.Appendix C. Cost of evaluating the performance functionals 109/40 /1.0 //,/1.0 /11.401.2cLlE0.80.6 --:1O 12 14 16 18 20Number of states in plant modelFigure C.51: The cost of estimating the frequency response of PBk at 370 (dashed) and 92 (solid) pointsfrom measurements relative to using a plant model versus the number of states in the plant model.C.5 The “H norm” functionalConsider the “H norm” functional (3.44) with Wi(s) = 1, W2(s) = 0, and H = PQ4)Hinf(X) = sup IP(jw)QN(z,jw)I. (C.91)C.5.1 Evaluating the functional using the plant modelThe method used to evaluate (C.91) is to discretize an interval [0, Wma], where is a sufficientlylarge frequency, by N3 points, and to find the maximum magnitude of the frequency response of PQ”on the interval.Expressing pQN as a function of x,p(s)QN(s)= XkPBk=1sB, (C.92)where A, B and Ck are the state space matrices of PBk. It is obvious that by computing the frequencyAppendix C. Cost of evaluating the performance function als 110response of FBk before the optimization begins, the frequency response of pQN can be obtained bymultiplication at each iteration of the optimization algorithm. The frequency response of PB,, is computed from its state-space form. In MATLAB this is done by first computing A, the Flessenberg formof A, and applying the similarity transformation T associated with the llessenberg transformation to Band Ck. Then the frequency response is computed usingPBk(s) = — A)1B, (C.93)where A = THAT, B = THB, and Ck = CkT (the superscript H denotes the complex conjugatetranspose). Including all matrix multiplications the cost of evaluating (C.93) at N3 frequency pointsrequires about 6N + 2(2N — Na) + (Na3 + N + 2Na — 1)N5 flops, where Na is the dimension ofA. The cost of obtaining the flessenberg form of A is roughly N3, and the cost of the matrix divisionin (C.93) is about N3 + N2 [Mac89j, For example, with Na = 40, N8 = 370, the total cost is about25 Mflops. The cost computing the frequency response of PBk using simulation versus the number ofstates in the plant model is in Figure C.52. Note that the cost obtained in simulation at 20 model statesis 26 Mflops which agrees with the estimate obtained above.Appendix C. Cost of evaluating the performance function als 111The MATLAB program used to generate Figure C.52 is:clearflops_hin_accum= [1;% Plant transfer function denominator.denon=l;% 370 frequency pts.f=logspace(—4,3,370)% State space form of B_i.N=20;p=l; % p is arbitraryAq=-ones (N) *2*p;Aq=tril (Ag, 0) +p*eye (N);Bq=ones(N,l);Cq(l:N)=_2*p*ones(l,N);Dq=l;for k=l:20,% Plant P_k=l/(l+s)”k.denon=conv(denom, [1 1]);[Ap,Bp,Cp, Dpi =tf2ss(l,denon);% PB_i.[Apq,Bpq,Cpq,Dpqi=series(Ap,Bp,Cp,Dp,Aq,Bq,Cq,Dq);flops(0)% Compute frequency response of PB_k and count% the number of operations.% Balance Aft,ai = balance(Apq);b = t \ Bpq;c = Cpq * t;% Reduce A to I-lesenburg form[p,ai = hess(a);% Apply similarity transformations from Hessenberg% reduction to B and C:b=p1 * b;c = c * p;p = ltifr(a,b,aqrt(_l)*f);% Conpute frequency response of PBk.for k=l:N,PB(k,l:length(f)) = (c * g);endflops_hin_accum = [flops_hin_accum flops);endplot (1:20, flops_hin_accum)xlabel(’Nunber of states in plant model’)ylabel (‘Floating point operations’)At each iteration of the optimization algorithm, (2N— 1)N flops are required to perform the matrixmultiplication in (C.92).C,5.2 Evaluating the functional using measurementsThe method to evaluate (C.91) using measurements is similar to the one used in Section C.4.2 forthe “H2 norm” functional. The cost of estimating the frequency response of PBk at 370 points fromAppendix C. Cost of evaluating the performance functionals 112C’)C0C’)El)a0C0al0)CC’)01.U-Figure C.52: The number of flops versus the number of states in the plant model to evaluate the frequencyresponse of PBk, k = 1, . . . ,20, at 370 freqnency points.measurements relative to using the plant model versus the number of states in the plant model is shownin Figure C.53. The same programs are used to compute the cost of obtaining the frequency response ofPBk from measurements in Section C.4.2; the relative cost is the cost obtained in Section C.5.1 dividedby cost using measurements.Since at each iteration the multiplication in (C.92) must be performed regardless of whether measurements or the plant model is being used, the relative cost of evaluating the functional during optimizationis one.2 4 6 8 10 12 14 16 18 20Number of states in plant modelAppendix C. Cost of evaluating the performance functionals 113Figure C.53: The relative cost of computing the frequency response of PBk, Ic = 1, . . . ,20 at 370frequency points using the plant model and using measurements versus the number of states in the plantmodel.2 4 6 8 10 12 14 16 18 20Number of states in plant modelAppendix DCost of computing a time domain response in MATLABThis appendix outlines the cost of computing a time domain response using MATLAB.First the continuous-time system is discretized using a zero-order hold. This involves computing theexponential of a matrix. The state-space matrices of the discretized system are Ad and Bd, given byAd = expfhBd= J exprAdrB,0where h is the sampling interval. MATLAB uses a scaling and squaring with diagonal Pade approximation method to compute the exponential of a matrix. This method is described in {GL89] and requiresroughly (6 + j + )N flops, wherej = max{O, 1 + floor(log2(AII))},N = b + a, a is the dimension of A, and nb the number of columns in B.The state-space equations of the discretized system arex(k + 1) = Adx(k) + Bdu(k)y(k) = Cx(k) + Du(k). (D.94)Hence, the state history can be obtained in (2N + N) length(u) flops and the response obtained fromthe state in (2N + 1)length(u) flops. Hence, including the cost of calculating the exponential of a matrixthe total cost of computing the response is (6 + j + 1/3)(na + i) + (2n + 3fla + 1) length(u) flops. Forexample, if the plant has 10 states, IIAII,0 = 1, and the number of time points is 100 then the cost ofcomputing a time response is roughly 33000 flops.114Appendix ELoop Transfer Recovery (LTR) design methodE.1 Linear Quadratic Gaussian (LQG) design methodThe plant is assumed described byth = Ax+B’Up+WprocYp = CX+Vsens, (E.95)where and V8, are independent zero-mean Gaussian stochastic processes with spectral densitymatrices W > 0 and V > 0 respectively.The LQG method finds the controller which minimizes the costlirnE{x(t)TQx(t) + u(t)TRu(t)}, (E.96)where Q and R are weighting matrices satisfying Q > 0 and R 0. The standard assumptions arethat the systems with state-space realizations (A, B,Q112C, 0) and (A, W1/2C, 0) are stabilizable anddetectable.The LQG-optimal controller is given byK(s) = —K,(sI—A+BK+KjC)1j,where the optimal state feedback matrix is= R_iBTP, (E.97)and P is the solution to the Lyapunov equationATP + PA - PBRl BTP + Q =0. (E.98)The Kalman filter gain matrix isK1 = P1 CT V_i , (E.99)115Appendix E. Loop Transfer Recovery (LTR) design method 116where P1 is the solution to the Lyapunov equationAP1 + P1AT - i1 CT V’CP1+ W =0.The block diagram of the LQG controller structure is shown in Figure E.54.(E.100)The difficulty with the LQG method is that only rules of thumb exist to select the weights W, V, Q, R(see [BH75] for example). Design methods such Loop Transfer Recovery (LTR) have been developedwhich systematically adjust the design parameters required for the LQG method [BBB+88].In this thesis it is convenient to consider the augmented plant shown in Figure E.55 which allows usto apply the LTR design method in a more flexible way. First we show how this plant can be describedby state-space equations of the form (E.95).The process noise d is filtered white-noise with power spectral density Wd having a realization[Ad, Bd, Cd, Dd]. The state space equations for the process noise are= Ad(+Bdvd =where v is Gaussian white noise. The state-space equations for the measurement noise in are= Amn+Bmp(E.101)Figure E.54: The LQG controller structure.in = Cml)+9, (E.102)Appendix P2. Loop Transfer Recovery (LTR) design method 117where S and p are Gaussian white-noise, and the transfer function Wm in Figure E.55 has the realization[Am, Bm, Cm, 0]. As well, it is assumed p, S and v are independent. Combining (E.101), (E.102) andthe state-space equations of the plant yieldsth ABCd 0 x BDd 0 Bii( = 0 Ad 0 ( + Bd 0 + 0 UpIA7) 0 0 Am 7] 0 Bm 0Yp = [c 0 Cm] (E.103)Since (E.103) is of the form (E.95) the LQG method can be applied. There are two reasons for augmentingthe plant: it makes the selection of the LQG design weights more intuitive (discussed in the next section)and it gives us a another interpretation for the LTR cost function allowing a comparison between thesolution obtained by the analytic design and that obtained using the tuning method (see the discussionin Chapter 5).Figure E.55: The augmented plant used in the LTR design method.Appendix F. Loop Transfer Recovery (LTR) design method 118L2 LTR design procedureThis section discusses how the LTR design procedure is used in the first example in Chapter 6 (involvinga SISO plant). The procedure is:1. Choose the spectral factor transfer matrices Wd, Wm so that —G(sI — A)1K is a desirableopen-loop gain by increasing Wd relative to Wm at frequencies where K1 is desired large. EitherWd or Wm can be chosen to shape —C(sI — A)1. In the example problem in Chapter 5, iii istaken as white noise (Am, Bm, and Cm are zero) and the spectral factor of d is used to shape—C(sI — A)1K.The positive gain margin of —C(sI — A)1 is guaranteed to be infinite and itsphase margin greater than 60 degrees.2. Obtain the optimal state feedback matrix of an LQG-optimal controller K by setting Q = I andR = p1 in (E.98). The parameter p is decreased until KP (K uses K just determined and K1from Step 1) converges to —C(sI — A)’K1 over a sufficiently wide range of frequencies.Appendix FObtaining the frequency response data for the identification stepThis section describes the experiments performed on the system shown in Figure F.56 to obtain thedata (4.51). The additive noise term v is included to represent the effect of sensor noise, nonzero initialconditions, plant output disturbances or finite-length identification experiments. The noise term isassumed to be unknown but bounded in magnitude for all time by ij.Figure F.56: Identification experiment setup.To obtain a bound on the error in the estimate of the frequency response it is assumed the plant’simpulse response belongs to the set{ f IfQ)Ietdt M’ (F.104)The input applied to the system isuQ) = ae’, 0 t (n — 1)T.The nth sample of the output of the systemy(nT)= (p * u)(nT) + v(nT),where T is the sampling interval. An estimate of the plant’s frequency response is— y((n — 1)T)—UV+y+119Appendix F. Obtaining the frequency response data for the identification step 120and a bound in the error in P(jw) islip — iioo= :+ f p(t)e_iwdt—_________< sup I p()e_iwtdt — çn_1)T p(T)ceiw((n_l)T_T)dT + 7)— WER+ J0 ae3w(11)Too (n—1)T= p(t)edt— f p(r)e_Tdr — aeiw(n_l)T< J p(t)e_tdt + (F.105)(n—1)TMe11)7’+ (F.106)To go from (F.105) to (F.106), (F.104) andJ lp(t)le1Tdt <(n—1)TJ p(t)idt(n—i )Twere used. The measurement of n samples of the output is repeated for each frequency as describedin (4.51). Hence, a bound for k, the error in the frequency response data (4.51), is given by (F.106)a function of the time domain bound on the additive noise, the experiment duration, the amplitude ofthe input, and information about the plant. Note that in practice it is not always possible to apply asinusoidal input to the plant.Bibliography[BB91] S. Boyd and C.H. Barratt. Linear Controller Design Limits of Performance. Prentice-Hall,1991.[BBB’88] S. Boyd, V. Balakrishna, C. Barratt, N. Khraishi, X. Li, D. Meyer, and S. Norman. A newCAD method and associated architectures for linear controllers. Proc. IEEE, 33(3):268—283,March 1988.{BBN9O] S. Boyd, C. Barratt, and S. Norman. Linear controller design: Limits of performance viaconvex optimization. Proc. IEEE, 78(3):529—574, March 1990.[BCFB93] S. Boyd, L. El Chaoui, E. Feron, and V. Balakrishnan. Linear matrix inequalities in systemsand control theory. draft, June 1993.[BH75] A.E. Bryson and Y.C. Ho. Applied Optimal Control. Hemisphere Publishing, 1975.[BS8OJ R.K. Brayton and R. Spence. Sensitivity and Optimization, Elsevier, 1980.[BW9OJ D.S. Bernstein and B. Wie. A benchmark problem for robust controller synthesis. In Proc.American Control Conf, pages 961—962, 1990.[BYM92I D.S. Bayard, Y. Yam, and E. Mettler. A criterion for joint optimization of identification androbust control. IEEE Trans. Aut. Control, 37(7):986—991, July 1992.[C1a83] F.H. Clarke. Optimization and Nonsmooth Analysis. Canadian Math. Soc. Series. John Wiley& Sons, 1983.[DF81] E. J. Davison and I.J. Ferguson. The design of controllers for the multivariable robust servomechanism problem using parameter optimization methods. IEEE Trans. Aut. Control,26(2):93—110, February 1981.[DFT9O] J. Doyle, B. Francis, and A. Tannenbaum. Feedback Control Theory. Macmillan, 1990.[DLMS8O] C.A. Desoer, R.W. Liu, J. Murray, and R. Saeks. Feedback system design: The fractionalrepresentation approach to analysis and synthesis. IEEE Trans. Aut. Control, 25(3):399—412,June 1980.[DS81] J. Doyle and G. Stein. Multivariable feedback design. Concepts for a classical/modern synthesis. IEEE Trans. Aut. Control, 26(1):4—16, January 1981.[Edm79] J.M. Edmunds. Control system design and analysis using closed-loop nyquist and bode arrays.ml. J. Control, 30(5):773—802, November 1979.[EK85] J .G. Ecker and M. Kupferschmid. A computational comparison of the ellipsoid algorithmwith several nonlinear programming algorithms. SIAM J. Control and Opt., 23(5):657—673,September 1985.[Feg64] K.A. Fegley. Designing sample-data control systems by linear programming. IEEE Trans. onAppi. and Ind., 83:198—200, May 1964.121Bibliography 122[Fra87] B. Francis. A Course in H Control Theory. Springer-Verlag, 1987.[GD83] C.L. Gustafson and C.A. Desoer. Controller design for linear multivariable feedback systemswith stable plant, using optimization with inequality constraints. mt. J. Control, 37(5):881—907, May 1983.[GD85] C.L. Gustafson and C.A. Desoer. A CAD methodology for linear multivariable feedback systems based on algebraic theory. ml. J. Control, 41(3):653—675, March 1985.[GKL89] G. Gu, P.P. Khargonekar, and E.B. Lee. Approximations of infinite dimensional systems.IEEE Trans. Ant. Control, 34(6):610—618, June 1989.[GL89] G. Golub and C. Van Loan. Matrix Computations. John Hopkins Univ. Press, second editionedition, 1989.[GLS88] M. Grötschel, L. Lovsz, and A. Schrijver. Geometric Algorithms and Combinatorial Optimization. Volume of Algorithms and Combinatorics. Springer-Verlag, 1988.[Gra9O] A. Grace. Optimization Toolbox for use with MATLAB. The Math Works, Inc., 1990.[11an89] F,R. Hansen. A Fractional Representation Approach to Closed-loop System Identification andExperiment Design. PhD thesis, Stanford University, California, 1989.{HJN89] A.J. Helmicki, C.A. Jacobson, and C.N. Nett. H identification of stable LSI systems: Ascheme with direct application to controller design. In Proc. American Control Conf, pages1428—1434, 1989.[HJN9O] A.J. Helmicki, C.A. Jacobson, and C.N. Nett. Identification in H: A robustly convergent,nonlinear algorithm. In Proc. American Control Conf, pages 386—391, 1990.[HJN91a] A.J. Helmicki, C.A. Jacobson, and C.N. Nett. Control oriented system identification: A worst-case/deterministic approach in H. IEEE Trans. Aut. Control, 36(10):1163—1176, October1991.[HJN91b] A.J. Helmicki, C.A. Jacobson, and C.N. Nett, Distributed Parameter Control Systems: NewTrends and Applications, volume 10, pages 301—332. Marcel Dekker, 1991.[HJN91c] A.J. Helmicki, C.A. Jacobson, and C.N. Nett. Fundamentals of control-oriented system identification and their application in H. In Proc. American Control Conf, pages 89—99, 1991.[HJN92] A.J. Helmicki, C.A. Jacobson, and C.N. Nett. Worst-case/deterministic identification in H:The continuous-time case. IEEE Trans. Ant. Control, 37(5):1163—1176, May 1992.[HL87] Q. Haung and R.W. Liu. A necessary and sufficient condition for stability of a perturbedsystem. IEEE Trans. Aut. Control, 32(4):337—340, April 1987.[Hor63] I.M. Horowitz. Synthesis of Feedback Systems. Academic Press, 1963.[Kau54] W.H. Kautz. Transient synthesis in the time domain. IRE Transactions on Circuit Theory,1(3):29—39, July 1954.{Kuc79] V. Kucera. The Polynomial Equation Approach. John Wiley & Sons, 1979.[Kwa85] H. Kwakernaak. Minimax frequency domain performance and robustness optimization of linearsystems. IEEE Trans. Ant. Control, 30(10):994—1004, October 1985.Bibliography 123[Lju87] L. Ljung. System Identification: Theory for the User. Prentice-Hall, 1987.[Mac89] J.M. Maciejowski. Multivariable Feedback Design. Addison-Wesley, 1989.[MV87] C,C.H. Ma and M. Vidyasagar. Direct globally convergent adaptive regulation of stable multivariable plants. IEEE Trans. Aut. Control, 32(6):543—547, June 1987.[MZ89] M. Moran and E. Zafiriou. Robust Process Control. Prentice-Hall, 1989.[NGK57] G.C. Newton, L.A. Gould, and J.F. Kaiser. Analytical Design of Feedback Controls. JohnWiley & Sons, second edition, 1957.[NJM88] C.N. Nett, C.A. Jacobson, and A.T. Miller. An integrated approach to controls and diagnostics:The 4-parameter controller. In Proc. Joint Auto. Cont. Cont. Conf., pages 824—835, 1988.[NT86] W.T. Nye and A.L. Tits. An application-oriented, optimization-based methodology for interactive design of engineering systems. mt. J. Control, 46(6):1693—1721, June 1986.[0B90] C. Oakley and C. Barratt. End-point controller design for an experimental two-link flexiblemanipulator using convex optimization. In Proc. American Control Conf, pages 1752—1759,1990.[PMS84] E. Polak, D.Q. Mayne, and D.M. Stimler. Control system design via semi-infinite optimization:A review. Proc. IEEE, 72(12):1777—1794, December 1984.[Po187] E. Polak. On the mathematical foundations of nondifferentiable optimization in engineeringdesign. SIAM Review, 29(1):21—89, March 1987.[PS89] E. Polak and S.E. Salcudean. On the design of linear multivariable feedback systems via constrained nondifferentiable optimization in H spaces. IEEE Trans. Aut. Control, 34(3):268—276, March 1989.[5a186] S.E. Salcudean. Algorithms for Optimal Feedback Control. PhD thesis, University of California,Berkeley, 1986.[5D92] R.S. Smith and J.C. Doyle. Model validation: A connection between robust control and identification. IEEE Trans. Aut. Control, 37(7):942—952, July 1992.[SDL91j C.M. Seaman, A.A. Desrochers, and G.F. List. Identification and Multi-Objective ControllerThning of a Plastic Injection Molding Machine. In Proc. American Control Conf., pages 262—267, 1991.[Str84] K. R. Stromberg. An Introduction to Classical Real Analysis. Wadsworth Press, 1984.[Vid85] M. Vidyasagar. Control System Synthesis: A Factorization Approach. MIT Press, 1985.[VJ84] M. Verma and E. Jonckheere. L-compensation with mixed sensitivity as a broadband matching problem. Systems and Control letters, 4:125—129, 1984.[Wah9l] B. Wahlberg. Identification of resonant systems using Kautz filters. In Proc. IEEE Conf onDecision and Control, pages 2005—2010, 1991.[YBJ76] D.C. Youla, J.J. Bongiorno, and H.A. Jabr. Modern Wiener-Hopf design of optimal controllers— part I: The single input-output case. IEEE Trans. Aut. Control, 21(1):3—13, February 1976.Bibliography 124[ZA73J V. Zakian and U. Al-Naib. Design of dynamical and control systems by method of inequalities.Proc. lEE, 120:1421—1427, 1973.[Zam8l] G. Zames. Feedback and optimal sensitivity: Model reference transformations, weighted seminorms, and approximate inverses. In Proc. 17th Allerton Conf, pages 744—752, 1981.[Z093] G. Zames and J.G. Owen. Duality theory for MIMO robust disturbance rejection. IEEE Trans.Aut. Control, 38(5):743—752, May 1993.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Optimization-based controller tuning using the Q-parameterization
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Optimization-based controller tuning using the Q-parameterization Lynch, Alan F. 1993
pdf
Page Metadata
Item Metadata
Title | Optimization-based controller tuning using the Q-parameterization |
Creator |
Lynch, Alan F. |
Date Issued | 1993 |
Description | We consider the problem of tuning a closed-loop linear controller for a continuous-time linear, lumped, exponentially stable plant using measurements. A two-step tuning method is proposed based on an “internal model” controller structure. The motivation for developing the method is to reduce the effect of modeling error on the design and to reduce the amount of computation required compared to a similar design performed off-line. The first step is to solve a controller design problem formulated as a sequence of convex optimization problems in which the cost and constraint functionals and their descent directions are computed directly from plant measurements. The designer is comfortable in specifying desired performance by adjusting the weights of a weighted-max cost function composed of a wide range of time and frequency domain performance functionals. To implement the “internal model” controller, a plant model is obtained in the second step. The identification objective makes the plant model depend on the desired closed-loop performance. To ensure a robust design, the model satisfies a worst-case error bound which depends on information about the plant, the experiment duration, the noise level, and the order of the model. |
Extent | 2337429 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-02-24 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0064876 |
URI | http://hdl.handle.net/2429/5019 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1994-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1994-0093.pdf [ 2.23MB ]
- Metadata
- JSON: 831-1.0064876.json
- JSON-LD: 831-1.0064876-ld.json
- RDF/XML (Pretty): 831-1.0064876-rdf.xml
- RDF/JSON: 831-1.0064876-rdf.json
- Turtle: 831-1.0064876-turtle.txt
- N-Triples: 831-1.0064876-rdf-ntriples.txt
- Original Record: 831-1.0064876-source.json
- Full Text
- 831-1.0064876-fulltext.txt
- Citation
- 831-1.0064876.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0064876/manifest