MODERN CONTROL SYNTHESIS AND ANALYSIS IN A
MATHEMATICA/CONTROL SYSTEM PROFESSIONAL ENVIRONMENT

David J. Dana-Bashian
The Boeing Company
3370 Miraloma Avenue
Anaheim, CA  92803
Attn.:  Dept. 797, M/S DC 27
david.j.dana-bashian@boeing.com (e-mail)

Copyright, The Boeing Company, All Rights Reserved

Abstract

For a given plant, an important objective of control theory is to design a controller whose corresponding closed-loop system satisfies a given set of specifications.  Another important objective of control theory is to reduce the number of states of a given plant model to be no greater than a prespecified positive integer while satisfying a given set of specifications.  One means to accomplish this task is to transform the specifications into a performance index, minimize the index with respect to some measure, and compute the resulting controller.  One may define either optimization problem using (H^2), (H^∞), or mixed-norm (H^2)/(H^∞) criteria in a finite-dimensional linear-time-invariant setting.  Haddad, Bernstein, and others [1]-[5] analyze a number of aspects of each of these problems.  Please consult the References section.

In 1996 Wolfram Research, Inc., introduced Control System Professional (CSP) to solve a wide range of control problems, including classical control design and analysis and introductory modern control design.  This section describes additions the author has made to CSP, all of which harmonize with current CSP conventions and usage.

Two major additions are full- and reduced-order (H^2), (H^∞), and (1-, 2-) (H^2)/(H^∞) controller design routines and (H^2) and (H^2)/(H^∞) model reduction routines.  Both sets of routines, based on the optimal projection equations of Haddad and Bernstein, are about as close to "no-guesswork" as possible.  Compactness of the routines results from using the so-called augmented plant approach, which transparently allows cross-coupling and obviates the need for separate computation of estimator and regulator gains and assembly of a controller or a reduced-order model.  Auxilliary routines allow the user to analyze a system for which individual states of a system are controllable/uncontrollable, observable/unobservable, and/or to prove that a result satisfies 2-norm and/or ∞-norm specifications.

Another major addition is incorporation of a numerical time-domain simulation routine which is considerably more flexible, efficient, and compact than the corresponding CSP routine.  The routine takes advantage of the form of the input corresponding to the response requested to cut the time to compute the response.  The routine transparently allows the user to output only a uniformly-spaced fraction of the points computed, saving the computer both time and space.  The routine also obviates the need to add states to compute the response to a time-delayed input, thus avoiding output spread over many screens' worth of space.  All told, the user may choose from 60 different types of responses using a single command.

Introduction

Anyone in the control community has a continuing need to obtain controllers which satisfy a given set of specifications and to analyze the resulting closed-loop performance to see how well one meets those specifications.  The development includes code suitable for modern control and model reduction and associated analysis.

The modern control portion represents a compact method for designing controllers satisfying 2-norm and/or ∞-norm constraints.  The compactness results from using the so-called augmented plant approach, commonplace in modern control papers, for “building in” a priori specifications.  This approach, which also transparently allows cross-coupling, obviates the need for separate computation of estimator and regulator gains and assembly of the controller.  Satisfaction of the resulting equations results in corresponding closed-loop systems which are guaranteed robustly stable and nominally performing to the degree of the desired specifications.  The development allows the user to design continuous controllers satisfying full- or reduced-order (H^2), (H^∞), or (1-, 2-) (H^2)/(H^∞) criteria or to design discrete controllers satisfying full- or reduced-order (H^2) criteria.  The model reduction portioin allows the user to reduce continuous models satisfying (H^2) or (H^2)/(H^∞) criteria or discrete models satisfying (H^2) criteria.

Other listed additions in the development represent convenient ways to analyze systems.  Proving a system satisfies 2-norm and/or ∞-norm constraints requires computing these system norms.  The 2-norm computation is straightforward from the definitions of matrix trace and controllability and observability gramians.  The ∞-norm computation uses the quadratically convergent method of Boyd and Balakrishnan.  The damping computation is straightforward.  The Popov-Belevitch-Hautus (PBH) rank tests provide an alternative means to the Kalman decomposition to compute the controllability and observability and corresponding ranks of systems.  In particular, PBH rank tests also tell which individual states are controllable/uncontrollable, observable/unobservable.

The numerical time-domain simulation routine takes advantage of the form of the input corresponding to the response requested to cut the time to compute the response.  Options for the routine include a numerical plot sampling period (if the user wants to output only a uniformly-spaced fraction of the points computed), an initial state vector, and a time delay.  The routine can accommodate time-delayed inputs without adding states to compute the response.  This feature avoids output spread over many screens' worth of space.  Computing even the continuous response to a continuous impulse presents no difficulties to the development, even if that impulse has a time delay, even if the system feedthrough matrix is nonzero.  All told, the user may choose from 60 different types of responses using a single command:  (discrete or continuous) x (unforced or impulse or step or ramp or general) x (output or state) x (the three aforementioned options).

I.  Mathematica Command Synopsis

Let a, b, c, d, d1, d2, e1, e2, e1inf, einf, e2inf stand respectively for matrices A, B, C, D, (Subscript[D, 1]), (Subscript[D, 2]), (Subscript[E, 1]), (Subscript[E, 2]), (Subscript[E, 1 ∞]), (Subscript[E, ∞]), (Subscript[E, 2 ∞]), and let ed and ee stand respectively for the model reduction disturbance and estimation weighting matrices (Subscript[E, d]) and (Subscript[E, e]), respectively.  Let gam stand for the upper bound for the so-called "γ-iteration" of (H^∞) theory.  Let system or systemi stand for any of the CSP constructs for its state-space description, such as StateSpace[a,b,c,d] or, if d is zero, StateSpace[a,b,c].  Then the general synopsis for the new main commands for CSP is as follows.

H2Controller[system, d1, d2, e1, e2]

compute a continuous controller from continuous system meeting (H^2) specifications only

HinfinityController[system, gam, d1, d2, e1inf, einf, e2inf]

compute a continuous controller from continuous system meeting (H^∞) specifications only

Options for HinfinityController:

Iterate:  performs γ-iteration with upper bound gam if True; default False

GammaFractionTolerance:  minimum change in γ-iteration; default 0.01

GammaTolerance:  minimum change in observer γ-iteration; default 0.0001

H2HinfinityController[system, d1, d2, e1, e2, e1inf, einf, e2inf]

compute a continuous controller from continuous system meeting (H^2) and/or (H^∞) specifications

Options for H2HinfinityController (continuous case):

Iterate:  same as above

GammaFractionTolerance:  same as above

GammaTolerance:  same as above

Gam:   γ, as one value or as iteration upper bound; default Infinity

DbetauTol:  minimum change in controller authority; default 1

ControllerOrder:  positive integer upper bound for controller; default n, the number of states of a

InputFileName:  name beginning external file output; default "testfilename"

H2HinfinityController[system, d1, d2, e1, e2]

compute a discrete controller from discrete system meeting (H^2) specifications

Options for H2HinfinityController (discrete case):

ControllerOrder:  same as above

InputFileName:  save as above

H2HinfinityReducedOrderModel[system, ed, ee]

compute a reduced-order model from continuous system meeting (H^2) and optional (H^∞) specifications

or from discrete system meeting (H^2) specifications

Options for H2HinfinityReducedOrderModel:

Iterate:  same as above

GammaTolerance:  same as above (applies only to continuous model reduction)

ModelOrder:  positive integer upper bound less than n for reduced model

InputFileName:  same as above

NumericalResponse[system, u, ToString[type], timevector]

compute the type response to system from input u

type is of the form XYZ, where X is Discrete or Continuous, Y is Unforced or Impulse or Step or Ramp

or General, Z is Output or State

u is the matrix of inputs if Y is General, the input number else

system and type must be both Continuous or both Discrete

if X is Discrete, then timevector is nonexistent if Y is General, {t, tend} else, where t is a symbol, time

ranges between 0 and tend in increments of deltat, where deltat is the sampling rate of system,

tend being an integral multiple of delta

if X is Continuous, then timevector is of the form {t, tend, deltat}, where t is a symbol, time ranges

between 0 and tend in increments of deltat, tend being an integral multiple of deltat; if in addition

Y is General, input u must be consistent with timevector

Options for NumericalResponse:

NumericalPlotSampling Period

outputs a uniformly spaced fraction of points computed

InitialConditions

same usage as in CSP

TimeDelay

same usage as in CSP

The general synopsis for the new auxilliary commands is as follows.

Norm[matrix, type] or Norm[number]

compute the norm of a matrix; the norm of a number is its absolute value

choices for type are 1, 2, Infinity, "Frobenius"; default type 2

Grm[a,b] or Grm[Transpose[a],Transpose[c]]

compute the controllability Gramian of a and b or the observability Gramian of a and c using Laub's method

SVD[matrix]

compute the "non-economy size" singular value decomposition of a real matrix

for example, if Dimensions[matrix]={m,n}, m>n, SVD[matrix] produces {u,s,v} with Dimensions[u]={m,m}, Dimensions[s]={m,n}, Dimensions[v]={n,n}, matrix = u.s.Conjugate[Transpose[v]]; in contrast, {u,s,v} from SingularValues[matrix] has the "economy size" Dimensions[u]={n,m}, Dimensions[s]=n, Dimensions[v]={n,n}

SystemNorm[system, norm]

compute the norm of system, where norm is 2 or Infinity

DampingRatios[system] or Damping Ratios[a]

compute the damping ratios of system, which are those of a

NaturalFrequencies[system] or NaturalFrequencies[a]

compute the natural frequencies of system, which are those of a

Damp[system] or Damp[a]

list containing the list of damping ratios and the list of natural frequencies

PBHRank[system, ToString[type]]

compute the Popov-Belevitch-Hautus properties of a system

choices for type:

Controllable, Observable, ControllableRank, ObservableRank,

UncontrollableStates, UnobservableStates, ControllableStates, ObservableStates

UnityFeedback[system1, system2]

compute the system which results from negatively feeding back the negative outputs of system1 followed by system2 (which is not the same as positively feeding back the positive outputs of system 1 followed by system2), where the gain in the feedback loop is 1 (i.e., unity)

Sensitivity[system1, system2]

compute the sensitivity "Inverse[IdentityMatrix[n]+system1.system2]"; n=number of states in system1.system2

currently, system1 and system2 must both be single-input-single-output (i.e., n=1)

SimultaneousMargin[system1, system2]

compute the minimum distance from the critical point -IdentityMatrix[n] to system1.system2, n=dimension of  system1.system2; better as a measure of nearness to instability than both gain or phase margin

currently, system1 and system2 must both be single-input-single-output (i.e., n=1)

For demonstration information and sample runs, see the appendices which follow.

Appendix A

DEMONSTRATION:  CONTINUOUS CONTROL DESIGN

The coupled rotating disk problem [1] is an eight-state, one-input-one-output, neutrally stable plant with one pair of nonminimum phase zeros.  Except for the two pure integrators, the damping of all the poles is approximately 0.02.  This plant has a fairly significant reputation for being very difficult to construct a controller satisfying its specifications.  As such, this example provides an excellent opportunity to observe how relatively easy designing a reduced-order controller can be using this program.

Step-by-step instructions to construct and compare (H^2), (H^∞), and (H^2)/(H^∞) controllers for the coupled rotating disk problem:

Set current location to the desired subdirectory (or folder).  One should have "write-access" to that location.  Further, one should have access to data file "handb.m" as well.

Enter Mathematica.

Type (as indicated)

$Line = 0;

<< "ControlSystems`Master`"; << "H2Hinfinity`Master`";

$Line = 1;

<< "handb."; system = StateSpace[a, b, c, d]

$Line = 2;

These statements examine the controllable and observable properties of system using the PBHRank tests.  Since PBHRank finds system controllable and observable, the ranks are full, and all states appear.

{PBHRank[system, "Controllable"], 
      PBHRank[system, "Observable"], 
      PBHRank[system, "ControllableRank"], 
      PBHRank[system, "ObservableRank"]}

$Line = 3;

PBHRank[system, "ControllableStates"]

$Line = 4;

PBHRank[system, "ObservableStates"]

$Line = 5;

This statement displays the remaining matrices.  For the purposes of this demonstration, e2inf will be {{0.},{1.*(10^(-4)) }}.

e2inf = {{0. }, {1. *10^(-4)}}; {d1, d2, e1, e2, e1inf, einf, 
      e2inf}

$Line = 0;

The statements below ask for the central (H^∞) and (H^2) controllers.  To compute the (H^∞) controller, the program first computes the minimum γ for the observer equation, then uses this γ as a lower bound to find the minimum γ for the entire augmented plant.  Apparently, this γ-iteration improvement does not appear in any currently publicly available software.  The letter T signifies a passed test, F a failed test.  The first controller computed is an (H^2) controller (for γ=∞).  The program decreases γ if all tests pass, increases γ if not, by a simple bisection.  If there is a failed test, all tests beyond that point for the particular γ have no evaluation (signified I for indeterminate).  The iteration continues until the increment in γ falls below a preset small positive number.  In this case, all the tests pass for all values of γ tried.  All full-order controllers are unstable for γ<0.22.  The nonminimum-phase plant and unstable controller kinfinity give rise to unusual results for closedloophinfinity, later.

kinfinity = 
      HinfinityController[system, Infinity, d1, d2, e1inf, einf, e2inf, 
        Iterate -> True];

\(\(\($Line = 7; \)
\)\)

k2 = H2Controller[system, d1, d2, e1, e2];

$Line = 8;

These statements close the loop using the unity feedback construction.

closedloophinfinity = UnityFeedback[system, kinfinity];

$Line = 9;

closedlooph2 = UnityFeedback[system, k2];

$Line = 10;

These statements allow comparison of the eigenvalues of the closed loops.

Eigenvalues[closedloophinfinity[[1]]]

$Line = 11;

Eigenvalues[closedlooph2[[1]]]

$Line = 12;

These statements allow comparison of the 2-norms and ∞-norms of the closed loops.

SystemNorm[closedloophinfinity, 2]

$Line = 13;

SystemNorm[closedloophinfinity, Infinity]

$Line = 14;

SystemNorm[closedlooph2, 2]

$Line = 15;

SystemNorm[closedlooph2, Infinity]

$Line = 16;

These statements compute the damping ratios ([[1]]) and natural frequencies ([[2]]) of the open and closed loops.

Damp[system]

$Line = 17;

Damp[closedloophinfinity]

$Line = 18;

Damp[closedlooph2]

$Line = 19;

These statements allow comparison of the time-domain responses of the closed loops.  The computational software consists of Mathematica code, together with an Installed time-domain simulation executable.  The code takes advantage of the form of the input corresponding to the response requested to cut the time to compute the response.  The notation

(1)	"ContinuousStepOutput" signifies the request to compute the output response to a continuous unit-step input

(2)	{t,100,1} signifies the computation of the response will be from time t = 0 sec to t = 100 sec in uniform

increments of 1 sec with the appropriate hold (zero-order hold in this case)

(3)	NumericalPlotSamplingPeriod->10 signifies the output back to the Mathematica session will have

uniform sampling every 10 sec

(4)	TimeDelay->5 signifies the unit-step input will not occur until 5 sec after the input initial state vector

(denoted by Initial Conditions, whose meaning is the same as that in CSP)

All told, the user may choose from 60 different types of responses using this single command:  (discrete or continuous) x (unforced or impulse or step or ramp or general) x (output or state) x (the 3 options of NumericalPlotSamplingPeriod, InitialConditions, TimeDelay).

NumericalResponse[closedloophinfinity, 
      "ContinuousStepOutput", {t, 100, 1}, 
      NumericalPlotSamplingPeriod -> 10, TimeDelay -> 5, 
      InitialConditions -> {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 
          15, 16}]

$Line = 20;

NumericalResponse[closedlooph2, 
      "ContinuousStepOutput", {t, 100, 1}, 
      NumericalPlotSamplingPeriod -> 10, TimeDelay -> 5, 
      InitialConditions -> {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 
          15, 16}]

$Line = 21;

This statement asks for the (H^2)/(H^∞) controller.  To compute this controller, the program first computes the minimum γ for the observer ("Q") equation, then uses this γ to improve convergence characteristics of the solution process.  The program computes the full-order (H^2)/(H^∞) controller by iterating on the (H^2) regulator ("P") and (H^∞) regulator ("(OverHat[Q])") equations until the residual errors P_E and (OverHat[Q])_E are both below their respective prespecified tolerances.  A value of 1 for β (controller authority) indicates an attempt to use all the driving term possible using this method.  Since γ is high enough, the resulting controller is stable.

kh2hinfinity = 
      H2HinfinityController[system, d1, d2, e1, e2, e1inf, einf, e2inf, 
        InputFileName -> "handb", Gam -> 100, ControllerOrder -> 6];

$Line = 22;

This statement closes the loop using the unity feedback construction.

closedlooph2hinfinity = UnityFeedback[system, kh2hinfinity];

$Line = 23;

This statement computes the eigenvalues of the closed loop.

Eigenvalues[closedlooph2hinfinity[[1]]]

$Line = 24;

These statements compute the 2-norm and ∞-norm of the closed loop.  Comparison with results already supplied show that the closed loop exhibits much more 2-norm behavior than ∞-norm behavior, consistent with the high γ specification Gam->100.  To impose more ∞-norm behavior, one would reduce maximum γ level specification.

SystemNorm[closedlooph2hinfinity, 2]

$Line = 25;

SystemNorm[closedlooph2hinfinity, Infinity]

$Line = 26;

This statement computes the damping ratios ([[1]]) and natural frequencies ([[2]]) of the closed loop.

Damp[closedlooph2hinfinity]

$Line = 27;

This statement computes a time-domain response of the closed loop.  A comparison of this response to results already supplied may not be fair because the initial state vector must necessarily be different, since the closed loop has but 12 instead of 16 states.

NumericalResponse[closedlooph2hinfinity, 
      "ContinuousStepOutput", {t, 100, 1}, 
      NumericalPlotSamplingPeriod -> 10, TimeDelay -> 5, 
      InitialConditions -> {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}]

$Line = 28;

Appendix B

DEMONSTRATION:  DISCRETE MODEL REDUCTION PROBLEM

The file myles110.m contains state space matrices for a 2-input, 2-output, wing and vortex example, of order 110.  The user supplies matrices ed and ee, where the intensity of the disturbance is v=ed.Transpose[ed] and the model-reduction error-weighting matrix is r=Transpose[ee].ee.

Step-by-step instructions to construct a reduced-order model of order 6 for the above plant (follows as the previous set of instructions, except that the date file is myles110.m):

$ContextPath =. ; Remove["Global`*"]; $Line = 0;

<< "ControlSystems`Master`"; << "H2Hinfinity`Master`";

$Line = 1;

<< "myles110."; system = StateSpace[a, b, c, d, Sampled -> True];

$Line = 2;

ed = {{1, 0}, {0, 1}}; ee = {{1, 0}, {0, 1}};

$Line = 3;

reducedordermodel = 
      H2HinfinityReducedOrderModel[system, ed, ee, ModelOrder -> 6, 
        InputFileName -> "myles110"]

$Line = 4;

One may compute the output responses to, say, a discrete impulse input for, say, ten samples beyond the initial condition.  Note that Out[5] and Out[6], the output responses to the first input, are fairly close to each other, as are Out[7] and Out[8], the output responses to the second input, to each other.  One may find the singular values of reducedordermodel and system each over a wide range of frequencies, say, a log scale from (10^(-6)) rad/sec to (10^6) rad/sec, are very close to each other.  The time-domain responses and the singular values argue that reducedordermodel faithfully represents system.

NumericalResponse[reducedordermodel, 1, 
      "DiscreteImpulseOutput", {t, 10}]

$Line = 5;

NumericalResponse[system, 1, "DiscreteImpulseOutput", {t, 10}]

$Line = 6;

NumericalResponse[reducedordermodel, 2, 
      "DiscreteImpulseOutput", {t, 10}]

$Line = 7;

NumericalResponse[system, 2, "DiscreteImpulseOutput", {t, 10}]

$Line = 8;

References

D.S. Bernstein and W.M. Haddad, LQG Control with an (H^∞) Bound:  A Riccati Equation Approach, IEEE Transactions on Automatic Control 34 (1989) 293-305.

W.M. Haddad and D.S. Bernstein, Generalized Riccati Equations for the Full- and Reduced-Order Mixed-Norm (H^2)/(H^∞) Standard Problem, Systems and Control Letters 14 (1990) 185-197.

W.M. Haddad and D.S. Bernstein, Combined (H^2)/(H^∞) Model Reduction, Proc. 1989 American Control Conference 2660-2665.

D.S. Bernstein, L.D. Davis, and D.C. Hyland, The Optimal Projection Equations for Reduced-Order, Discrete-Time Modeling, Estimation, and Control, J. Guidance 9 (1986) 288-293.

D.J. Dana-Bashian, On Pure (H^∞) Reduced-Order Dynamic Controllers, Systems and Control Letters 37 (3) (05 July 1999) 177-179.

Vita

David Dana-Bashian, senior engineer scientist at The Boeing Company (before 04 August 1997, McDonnell Douglas Corporation), has researched nonclassical controller design approaches for over 10 years.  He specializes in transforming equations in the corresponding literature into powerful numerical software tools.  He holds a B.S., M.S., and B.S. in Chemistry, all from Purdue University, and is a member of the honorary chemical society Phi Lambda Upsilon.  He co-organized (and in the second case also co-chaired) the two invited sessions "Case Studies in Robust Control Design:  General Applications and Flight Control Applications" at the IEEE Conference on Decision and Control, December 1990.  He also authored the paper "Generic Controller Design for Autonomous Vehicles" at the Annual National Symposium of the Association for Unmanned Vehicle Systems in May 1994.