\(\(H\^∞\ Control\ Design\ on\ Frequency\ Domain\)\(
\)
      \)

Béla Paláncz
Budapest Technical University, Hungary
palancz@epito.bme.hu
http://geo2.fmt.bme.hu/palancz
Budapest, 2006

Abstract

In this case study we demonstrate, how to employ the Mathematica Application, Control System Professional Suite, (CSPS) together with the Professor Helton's package, OPTDesign to develope a robust controller for a nonlinear plant. First,  CSPS is used to linearize the nonlinear system in state space then to transform it into the frequency domain. In the second step, the OPTDesign package is applied to design optimal (H^∞) controller on frequency domain considering different constrains formulated by Helton and Merino (1998). In addition a new technique is suggested to consider the disturbance model as a constrains inequality. In the last step, we use CSPS again to transform the Laplace representation of  the controller dynamics back to state space and to find a minimal state space representation of it. This state space model of the controller was connected to the nonlinear plant model, and closed loop simulation of the controlled plant dynamics has been achieved.

Off[Syntax::stresc];

Off[General::spell];

Introduction

Let us consider a water storage vessel with a time dependent inlet flow, Qin[t]. The water level, h[t] can be controlled by the outlet flow, Qout[t] regulated by a control valve, c[t] as actuator.

Fig.1 Level control problem

The dynamical model of the controlled system can be developed as follows:

Model equations

The mass balance:

deq = ((Qin[t] - Qout[t])/A) == h'[t]

The outlet flow rate can be expressed by the valve "constant" and the water level:

Qout[t_] = c[t] (Sqrt[h])[t]

Combining them, we get the mathematical model of the system, which is a nonlinear, non-autonomous ordinary differential equation for the state-variable h[t].

deq

Nonlinear model

The general form of the nonlinear model is

\(x\& .  = f \((x, u)\)\)
y = g (x, u)

In our case

f = deq[[1]]

and

x = {h[t]}; y = {h[t]}; u = {Qin[t], c[t]}; g = y;

To design the control for the system , it should be linearized. Now, we load the CSPS application. However, in order to remove any packages easily, first we save the original content of the $ContextPath,

BasicContextPath=$ContextPath;

<< ControlSystems`

Linearization

The linearization should be carried out at a steady state, namely at (h0, Qin0, c0):

ControlObjectSS = 
      Linearize[f, g, {{h[t], h0}}, {{Qin[t], Qin0}, {c[t], c0}}]

From the state representation, one can get the equation form of the linearized system:

ControlObjectSS // EquationForm

Here we should keep in mind, that the variables x and u represent the deviation from the steady state. Our model in matrix form is:

TraditionalForm[%]

or its representation on frequency domain:

ControlObjectTF = 
      TransferFunction[s, ControlObjectSS] // Simplify

The numerical values are

h0 = 9. ; A = 10. ; c0 = 3. ; Qin0 = 9. ;

The transfer function between the level h[t] and the inlet flow rate Qin[t], as disturbance:

pd[s_] = ((2 (Sqrt[h0]))/(c0 + 2 A (Sqrt[h0]) s))

The transfer function between the level h[t] and the control valve position, c[t] as actuator:

p[s_] = -((2 h0)/(c0 + 2 A (Sqrt[h0]) s))

Now, we remove ControlSystems application and load the OPTDesign package, which is available on http://www.math.ucsd.edu/~anopt/

Remove[t, s, s$];

$ContextPath = BasicContextPath;

<< OPTDesign`

Concept of the control design

Consider the closed-loop system depicted in Fig.2

Fig.2 Concept of the control design

Let us consider C(s) and P(s) the transfer functions for the compensator and plant, respectively, then the closed-loop transfer function is,

\(T \((s)\)\  = \ \(P \((s)\)\ C \((s)\)\)\/\(1 + P \((s)\)\ C \((s)\)\)\)

Roughly speaking, to design a system one must choose all of the system's parts that have not been given in advance. This is done according to some criteria or requirements. In this case,  the most obvious way is to specify C(s). However, now we use another way, namely we specify T(s) instead of C(s). If T(s) is specified, the one can calculate C(s) as

C (s) = ((T (s))/((1 - T (s)) P (s)))

Under these circumstanses we refer to T(s) as a variable, as a designable transfer function.
The designable transfer function cannot take an arbitrary form. The set of admissible functions to which the designable T(s) must belong is a subset of  ℛ(ℋ^∞).  These rational functions are stable, proper and real (the coefficients are real). A rational function is stable if every pole of the function has negative real part and the relative degree is not positve. The relative degree of a rational function F(s) = ((N (s))/(D (s))) is,

d(F(s)):= degree(D(s)) - degree(N(s))

A function F(s) is called proper if d(F(s)) ≥ 0 and strictly proper if  d(F(s)) > 0. If  F(s) proper and stable, then it is bounded,

|F(i ω)| ≤ M for all ω∈R

\(\(\(The\ smallest\ such\ M\ denoted\ \((\(|\)\(F\)\(|\))\)\_∞\)\(,\)\)\)

\(\((\(|\)\(F\)\(|\))\)\_∞ := sup | F \((i\ ω)\) | \ 
        for\ \ all\ \ ω ∈ R\)

In addition a designable transfer function should ensure internal stability, namely there should be no cancellation of left-half-plane(unstable) poles and zeros in the product P(s) C(s).
The usual convention is to denote the set of all rational functions that are stable and real by  ℛ(ℋ^∞). Note that functions in  ℛ(ℋ^∞) can be described as proper rational functions with real coefficents and with no poles in the closed rigth half-plane RHP.

Performance Requirements

A frequency domain performance requirement is an inequality that the designable transfer function must satisfy, and an interval in ω for which the inequality is required to hold.

Disk inequalities

All of  the frequency domain performance requirements treated here can be written in the form of a disk inequality

|K(i ω) - T(i ω)| ≤ R(i ω)                for  (Subscript[ω, a ])≤ ω ≤ (Subscript[ω, b])

which T must satisfy. K and R are fixed functions that embody the desired specifications of the system. K is called the center of the disk and R is called the radius. Disk inequalities are easy to plot as regions in 3-D space and correspond to the inside of  a tubelike domain. If the frequency ω is fixed, then a disk inequality is represented in the complex plane by a region consisiting of a solid disk,

(Subscript[S, ω]) = {z : 
          K (i ω) - z  ≤ 
            R (i ω)}

Gain-phase margin

The gain-phase margin is a composite measure of stability representing the classical  gain margin and phase margin together. Define the gain -phase margin as

m = inf |1 + P C (iω)|   for ω

Graphically m is just the distance of the Nyquist plot of PC to the point -1 in complex plaine. One may get a constraint in form of disk inequality, see Fig. 3.,

\(\(\(|\)\(T \((i\ ω)\) - 
        1\)\(|\)\(\(≤\)\(\ \)\(1\/m\ \ \ \ \ for\ all\ ω\)\)\)\)

Fig. 3. The gain margin g, phase margin ϕ, and gain - phase margin m

If a given stable system has m near 0, then it is close to the unstable case, which is undesirable. If (Subscript[α, g]) denotes the largest value of 1/m considered to be acceptable, then

\(\(\(|\)\(1\  - 
        T \((i\ ω)\)\)\(|\)\(\(≤\)\(\ \ \ \)\(α\_g\ \ for\ \ \ all\ \ \ \ ω\)\)\)\)

Disturbance rejection

Systems are effected by disturbances, and one of the purposes of feedback control is to minimize the effect of them. The equivalent form of Fig.2. can be considered as

Fig. 4. The control scheme with the disturbance transfer function

where Qin(t) is the disturbance - generally \(\(d\^*\)\)(t). Let us suppose that the disturbance has the  transfer function (Subscript[p, d]), then

\(d \((s)\) = \(p\_d\) \((s)\)\ \(d\^*\) \((s)\)\)

The transfer function from disturbance d(s) to plant is

(Subscript[T, d  ->  h])  =  (1/(1 + P C)) = 1 - T

Now, we want to restrict the size of

\(T\_\(\(d\^*\) -> h\) = \ \(h\/\(d\^*\) = \(\(h\/d\) 
            d\/\(d\^*\) = \((1 - T)\)\ p\_d\)\)\)

namely

\(\(\(|\)\(T\_\(\(d\^*\) -> 
              y\)\)\(|\)\) = \(\(|\)\(\((1 - 
              T)\)\ p\_d\)\(|\)\(\(≤\)\(\ \)\(c\)\)\)\)

therefore the corresponding disk inequality constraint for T is

\(\(\(|\)\(1 - 
        T\)\(|\)\(\(≤\)\(c\/\(\(|\)\(p\_d\)\(|\)\)\)\)\)\)

where c is a contant value.

Bandwidth

Bandwidth of a system is commonly defined as the frequency (Subscript[ω, b]) at which  |T(i ω) | falls below a given constant times the low - frequency value of the input. This constant usually taken to be (1/(Sqrt[2])). A particular bandwith is required to ensure that the system at high frequency is not upset by noise, plant uncertainty, etc. The disk inequality form of bandwidth constraint is

\(\(\(|\)\(T \((i\ ω)\)\)\(|\)\(\(<\)\(\ \)\(α\_b\ \ \ \ \ \ for\)\)\(\ \ \)\(|\)\(ω\)\(|\)\(\ \)\(\(>\)\(\ \)\(ω\_b\)\)\)\)

Closed-loop roll-off

The open loop transfer function L = P C rolls off at high frequency, because the plant and the compensator are strictly proper

|P(i ω) C(i ω)| -> 0               if  ω -> ∞

Considering the closed loop transfer function

T  =  ((P C)/(1  -  P C))

for large frequency

\(\(\(|\)\(\ \)\(T\ \((i\ ω)\)\)\(|\)\(\(≤\)\(\ \ \)\(\(\(|\)\(P\ \((i\ ω)\)\ C \((i\ ω)\)\)\(|\)\)\/\(\(\(1\)\(\ \)\(-\)\) | \ 
              P \((i\ ω)\)\ C \((i\ ω)\) | \)\)\(\ \)\(≈\)\)\(\ \ \ \ \ \)\(|\)\(P\ \((i\ ω)\)\ C \((i\ ω)\)\)\(|\)\(\ \ \ \ \ \)\(for\ large\  ∣ ω ∣ \)\)\)

To obtain an inequality, C(i ω) must be eliminated. Therefore let

\(C \((i\ ω)\) ≤ \ α\_r\/\((\(|\)\(ω\)\(|\))\)\^n\)

\(If\ the\ designer\ specify\ n\ and\ α\_r\ \ as\ well\ as\ the\ the\ roll - off\ frequency\ ω\_r, \ 
    the\ constraint\ for\ \  | 
      T \((i\ ω)\) | \ is, \)

\(\(\(|\)\(\ \)\(T\ \((i\ ω)\)\)\(|\)\(\(≤\)\(α\_r\ \(\(|\)\(P \((i\ ω)\)\)\(|\)\)\/\((\(|\)\(ω\)\(|\))\)\^n\ \ \ \ \ \ \ \ \ \ \ \ \ \ for\)\(\ \)\(∣\)\(ω\)\(∣\)\(\ \)\(>\)\(\ \)\(ω\_r\)\)\)\)

Requirements envelope

The regions defined by the constraints can be put together providing a performance requirements envelope, see Fig. 5.,

Fig. 5. Requirements envelope

In our case, the parameters for the performance requirements are chosen as it follows,

c = 0.9; αr = 185. ; 
    n = 2; ωd = 0.275; ωb = 2. ; ωr = 4. ; αg = 
      2.5; αb = 0.9;

Then the disk inequalities can be written for the center,

k0[ω_] = 
        Which[0 ≤ Abs[ω] < ωb, 
          1.0, ωb ≤ Abs[ω], 0];

and for the radius

r0[ω_] = 
        Which[ 0 ≤  Abs[ω] < ωd, 
          c/Abs[pd[I ω]], 
                           ωd ≤ 
            Abs[ω] < ωb, αg, 
                           ωb ≤ 
            Abs[ω] < ωr, αb, 
                           ωr ≤ 
            Abs[ω],  αr Abs[
                p[I ω]]/(ω^n)];

Now we remove the Controlsystems     application and load the OPTDesign package,

The detailed description of these packages can be found in Helton and Merino (1998).
The requirements envelope can be drawn,

p1 = 
        EnvelopePlot[Radius -> r0, Center -> k0, 
          FrequencyBand -> {0. , ωr + 1}];

Fig. 6 The original and smoothed requirements envelopes

For further computations, the smoothed form of the envelope will be considered, see Fig. 6. A 3D plot of the requirement envelope can be produced with

EnvelopePlot3D[Radius -> r0, Center -> k0];

Fig. 7 3-D plot of the  requirements envelope

The next figure shows the discrete grid points, where the value of the T(i ω) function will be computed.

EnvelopeLogPlot[Radius -> r0, Center -> k0, 
        FrequencyBand -> {0. , ωr + 1}, 
        Discrete -> True];

Fig. 8  Plot of the discrete profile of the requirements envelope

Optimization

Concepts of the Optimization

Let us consider the closed loop transfer function of our system,

\(T \((s)\) = \(P \((s)\)\ C \((s)\)\)\/\(1\  - \ P \((s)\)\ C \((s)\)\)\)

or

P (s) C (s)  = ((T (s))/(1  -  T (s)))

where P(s) and C(s) are proper, real rational functions of  s, the system internally stable if  T(s) has no RHP poles, and there is no RHP pole-zero cancellation in the product PC.
The design problem is the following:

Given
1. a plant,  P(s)
2. a set of performance requirements, 𝒫
3. a set ℐ = {T(s) ∈ ℛ(ℋ^∞)} where T(s) is the closed loop transfer function of an internally stable system with plant P(s).

Determine T ∈ ℐ  if there exist that satisfy all the performance requirements in 𝒫. Such a T is called feasible.

Optimal design means that find a feasible T that yields the best performance.  If many satisfactory T exist, tighten or loosen the specification, and repeat the design.

Generating a performance function

As we have seen the performance requirements can be expressed in from of disk inequalities, which can be written in the following way, too:

\(\(1\/\((R \((i\ ω)\))\)\^2\) \((\(|\)\(K \((i\ ω)\)\  - \ T \((i\ ω)\)\)\(|\))\)\^2 ≤ \ 
      1\ \ \ for\ \ \ \ \(∀ \ ω\)\)

Now, for a given T calculate the largest value of the left-hand side of this, see Fig. 9.,

\(γ \((T)\) = \ \(sup\+ω\) \(1\/\((R \((i\ ω)\))\)\^2\) \((\(|\)\(K \((i\ ω)\)\  - \ T \((i\ ω)\)\)\(|\))\)\^2\)

Fig. 9. The largest values of different T functions

If a particular T is avaiable, then one can check if the original performance requirements are satisfied. In case of γ(T) ≤ 1, T is satisfactory.

For designable T and given center and radius functions, K and R respectively, the expression

\(\(1\/\((R \((i\ ω)\))\)\^2\) \((\(|\)\(K \((i\ ω)\)\  - \ T \((i\ ω)\)\)\(|\))\)\^2\)

is called the circular performance function, or simply performance function. The number γ(T) is the worst-case performance of T, but we will call it simply, the performance function of T.

Finding T with the best performance

\(A\ feasible\ function\ \(T\^*\)\ with\ property\ that\)

\(γ \((\(T\^*\))\) ≤ γ \((T)\)\ \ \(∀ \ 
          T\ feasible\)\)

is called optimal. The performance of the optimal designable transfer function is

\(\(γ\^*\)\  = \ γ \((T)\)\)

Optimization

The process of finding the optimal \(TraditionalForm\`\(T\^*\)\)and the corresponding optimal performance \(TraditionalForm\`\(γ\^*\)\)is called optimization. More precisely, we say that \(TraditionalForm\`\(T\^*\)\)is a solution to

\(\(γ\^*\)\  = \ 
      inf\+T\ \(sup\+ω\) \(1\/\((R \((i\ ω)\))\)\^2\) \((\(|\)\(K \((i\ ω)\)\  - \ T \((i\ ω)\)\)\(|\))\)\^2\)

There exist A, B ∈ ℛ(ℋ^∞), that if  T ∈ℐ, then

T = \(A\  + \ 
          B\ \(T\^+\)\ \ \ \ \ \ \ \ for\ \ \ some\ \ \(T\^+\)\) ∈ ℛ(ℋ^∞)

Design with OPTDesign

OPTDesign[p, Radius -> r0, Center -> k0, Nsmth -> 100, 
      Ngrid -> 512, Cdecay -> 2]

The output parameter \(\(γ\^*\)\) (gammastar) is less than 1. One can  put together the numerical values of  the function T with the grid on which it is defined,

Tpairs = OPTDParametrize[T];

\(Then\ the\ computed\ numerical\ values\ of\ \ the\ pairs, \ {T\_i, ω\_i}\ \ \ can\ be\ displayed\ together\ with\ the\ requirement\ envelopes, \ see\ Fig \( \(.10\)\(.\)\), \)

EnvelopePlot3D[Radius -> r0 , Center -> k0 , ClosedLoop -> T ];

Fig. 10. The numerical values of the optimal T inside the constraining envelope

Producing a Rational Compensator

The approximation of these points with a rational function

(Subscript[T, rat]) (s)  =  
      A (s)  +  B (s) (Subscript[T, 1  rat]) (s)

with stable function  (Subscript[T, 1  
            rat]) (s).

The proper rational functions A and B depend on the plant and incorporate automaticaly the internal stability requirements on T.

Trat[s_] = RationalModel[s, DegreeOfDenominator -> 5]

The compensator can be computed from the closed - loop transfer function,

C  =  (1/P) (T/(1  -  T))

therefore the transfer function of the compensator is,

comp1[s_] = 1/p[s] Trat[s]/(1 - Trat[s]) // Simplify

PlotZP[comp1[s], s];

Fig.11 The poles (x) and zeros (o) of the transfer function of the compensator

Now, no pole-zero pair cancellation is necessary. The poles and zeros of the closed loop transfer function

PlotZP[Trat[s], s];

Fig.12 The poles (x) and zeros (o) of the closed - loop transfer function

Analysis of the  result

In addition of Fig.10. further plots can be produced in order to illustrate of the result. We can plot the nonsmoothed solution in 2D.

p2 = 
        Plot[1/r0[ω] Abs[
              k0[ω] - Trat[I ω]], {ω, 0, 
            10}, PlotStyle -> RGBColor[1, 0, 0], 
          PlotRange -> {0, 1}];

Fig.13 Plot of the original (nonsmoothed) performance of Trat

The maximum value is

1/r0[ω] Abs[
          k0[ω] - Trat[I ω]] /. ω -> 
        0

One may put the figures  Fig.6 and Fig. 13 together

Show[{p1, p2}];

Fig.14. The nonsmoothed solution in the requirements envelope in 2D

The performance also can be analyzed in the Bode plots,

Tdisc = Discretize[Trat[s]];

BodeMagnitude[Tdisc, 
        PlotLabel -> "BodeMagnitudePlot"];

Fig.15 The Bode plot (magnitude) of the closed-loop transfer function

BodePhase[Tdisc, 
        PlotLabel -> "Phase Bode Plot"];

Fig.16 The Bode plot (phase) of the closed-loop transfer function

In addition one may display the Nyquist or the Nichols plots, too,

Nyquist[Tdisc, PlotLabel -> "Nyquist Plot"];

Fig.17 The Nyquist plot  of the closed-loop transfer function

Nichols[Tdisc, PlotLabel -> "Nichols Plot"];

Fig.18 The Nichols diagram of the closed-loop transfer function

There are certain diagnostic parameters, which indicate, how good our result is. We mentioned the γ already.

gammastar

Another parameter is for flatness diagnostics,which measures the "nonflatness".If it is small,then the diagnosis OK.

flat

The parameter grAlign should be theoretically zero,

grAlign

The numerical error of the computation can be characterized by the parameter ned, which should a small positive number,

ned

Quality of the disturbance rejection

We have seen that the transfer function from disturbance d(s) to plant is (1 - T(s))  (Subscript[p, d])(s). In our case

Tdy[s_] = (1 - Trat[s]) pd[s] // FullSimplify

The magnitude of this transfer function can indicate the effectivness of the disturbace rejection. If  it is less than 1, then the disturbace will be damped,

Plot[Abs[Tdy[I ω]], {ω, 0, 6}, 
        PlotRange -> All];

Fig .  19 The magnitude of the \(T\_\(\(d\^*\)\  -> \ h\)\) transfer function

The maximum is

Abs[Tdy[I ω]] /. ω -> 0

Now, we save the transfer function of the compensator is

DumpSave["foo.mx", comp1];

Simulating the nonlinear closed-loop dynamics

In order to  simulate the dynamics  nonlinear closed-loop we load again the ControlSystems application. But first, we should remove the OPTDesign application by restarting  Mathematica and continue the compuation from this point!

Off[General::spell];

<< ControlSystems`

Let us load the transfer function of the compensator

<< foo.mx

The transfer  function of the compensator,

K = TransferFunction[s, -comp1[s]]

SK = StateSpace[K]; EquationForm[SK]

One may try to reduce the states looking for  minimal realization

SKM = MinimalRealization[SK]; EquationForm[SKM]

It means, that the number of the states can not be reduced, which is in harmony with Fig. 11. Consequently we consider the original state space model, which is in traditional form

TK = TraditionalForm[SK]

Now, considering the Fig.2, the compensator dynamics can be represented in state space as

x'(t) = AA x(t) + BB (h(t) - h0)

c(t ) - c0 = CC x(t) +  DD (h(t) - h0)

These matrices can be deduced from TK,

AA = TK[[1, 1]]; MatrixForm[AA]

BB = TK[[1, 2]]; MatrixForm[BB]

CC = TK[[1, 3]]; MatrixForm[CC]

and DD is zero, therefore

Length[TK[[1]]]

We shall create the system of the differential equations describing the state of the compensator dynamics,

n = Length[Transpose[AA]]

X = Table[(Subscript[x, i])[t], {i, 1, n}];

deqx = 
      Map[# == 0 &, 
        Flatten[Table[
            D[X[[i]], t] - AA[[i]] . X - 
              BB[[i]] (h[t] - h0), {i, 1, n}]]]

The nonlinear plat model

deq=((-c[t] (Sqrt[h])[t] + Qin[t])/A)==h'[t];

so the output of the compensator model should be substituted into the nonlinear plant model

deqh = 
      Map[# == deq[[2]] &, 
        deq[[1]] /. c[t] -> c0 + CC . X]

The disturbance function, in our case, Qin(t) should be defined. Let us consider a unit step function  at t = 10

Qin[t_] = Qin0 + UnitStep[t - 10];

Th numerical values are for the plant model and for the initial (steady state) values,

h0 = 9. ; A = 10. ; c0 = 3. ; Qin0 = 9. ;

The initial conditions for the controller state variables are

incx = Map[(# /. t -> 0) == 0 &, X]

Then the system of differential equations, deqx and deqh can be solved simultanuesly,

sol = 
      NDSolve[Join[deqx, deqh, {h[0] == h0}, incx], 
        Join[X, {h}], {t, 0, 150},  MaxSteps -> 50000]

We plot now only the level as function of time

Plot[Evaluate[(h[t] /. sol)] - h0, {t, 0, 150}, 
        PlotStyle -> {RGBColor[0, 0, 1]}, Frame -> True, 
        PlotRange -> All];

Fig. 20 The level change as state variable of the plant after step change in the input flow rate (disturbance)

Here is the actuator dynamics, c(t)

Plot[Evaluate[c0 + CC . X] /. sol, {t, 0, 150}, 
        PlotStyle -> {RGBColor[0, 0, 1]}, Frame -> True, 
        PlotRange -> All];

Fig. 21 The  change of the control valve position as actuator during the control operation

Remark

Employing higher order approximation of the numerically evaluated closed loop transfer function, which provides smaller approximation error (see in Producing Rational Compensator), we get smaller deviation from the steady state, see Figures 22-23.

Figs.22-23 The response of the level in case of different order systems representing controller dynamics

Reference

Helton J . W .  and Merino O(1998).  Classical Control Using (H^∞) Methods,
 Theory,Optimization,and Design .  SIAM, Philadelphia, LCCCN : 98 - 86612