Polynomial
Solution of the SISO Mixed Sensitivity H-infinity Problem |
Signals, Systems and Control
Department
Faculty of Mathematical Sciences
University of Twente
The Netherlands
To download the application in .pdf format, click here. |
To obtain satisfactory high-frequency roll-off nonproper weighting functions may be needed. These cannot be directly handled in the conventional state space solution of the H-infinity problem. Nonproper weighting functions present no problems in the frequency domain solution of the H-infinity problem. The frequency domain solution may be implemented in terms of polynomial matrix manipulations.
We present the mixed sensitivity problem and its solution for single-input-single-output plants. Step by step the Toolbox function
mixeds is developed that implements the algorithm. A simple but relevant example is used for illustration. IntroductionFig. 1. Single-degree-of-freedom feedback loop
The sensitivity matrix S and input sensitivity matrix U of the feedback system are defined as
The mixed sensitivity problem is the problem of minimizing the infinity norm
ofwith suitably chosen weighting matrices
and , and V a suitably chosen shaping matrix. In the SISO case the infinity norm is given byThe problem may be reduced to a standard H-infinity problem by considering the block diagram of Fig. 2, which includes the weighting and shaping filters.
The diagram of Fig. 2 defines a standard problem whose generalized plant has the transfer matrix
Fig. 2. Mixed sensitivity configuration
The closed-loop transfer matrix from w to
is precisely the function H whose infinity-norm we wish to minimize.
There are many ways to solve this problem, but only the frequency domain solution (Kwakernaak, 1996) allows the generalized plant to have a nonproper transfer matrix G. To enhance robustness at high frequencies it usually is necessary to make the weighting filter
nonproper, which in turn makes G nonproper.For simplicity we develop the solution for the SISO case only. Suppose that
where the numerators and denominators are (scalar) polynomials. Note that P and V have the same denominators d this makes partial pole placement possible (Kwakernaak, 1993).
Numerical Examplewhere we let c = r = 1. Note that
is nonproper and, hence, the generalized plant G is nonproper. The various polynomials are defined by the following command lines:The partitioning is needed for the solution of the
problem. The following code lines define the various polynomial matrices:The * denotes conjugation, that is,
. In the next subsection it is seen how a spectral factorizationof
may be obtained. The spectral factor is a square rational matrix such that both and have all their poles in the open left-half plane. The matrix J is a signature matrix of the formwhere the two unit matrices do not necessarily have the same dimensions.
Given this spectral factorization, all compensators whose norm is less than the number
are of the form , whereU is a rational stable matrix such that
. In particular one may chose . Spectral FactorizationIt requires the following steps.
Then the desired spectral factor is
The first spectral factorization may actually be done analytically, because we have
Inspection shows that if each of the polynomials
, and is strictly Hurwitz then the desired spectral factorization may be rendered in slightly modified form asThe computation of the left-to-right conversion
may easily be coded.
The result is
-0.71 0.71 + s
0.71 + 1.7s + s^2 0.71 + 0.71s
-0.71 -0.71 + s
Lam =
0.71 + s 0.71
-0.71 0.71 + s
The second spectral factorization now takes the form
It is not difficult to write the necessary code lines
% Spectral factorization
Jgamma = eye(3); Jgamma(3,3)= -gamma^2;
DelDel = Del'*Jgamma*Del;
[Gam,J] = spf(DelDel)
Gam =
-1.3e+002 - 50s - s^2 -2.2e+002 - 0.046s
1.3e+002 + 48s + 0.17s^2 2.2e+002 - 3.8s
J =
1 0
0 -1
where for simplicity we choose
. It is advantageous to implement this computation as a right-to-left conversionso that the compensator transfer function is
Again this may be coded straightforwardly:
The output is
2.4e+002 + 1.6e+002s + 50s^2 + s^3
y =
63 + 1.8e+002s - 0.66s^2
Computation
of the Closed Loop Poles We use it to test closed-loop stability and to compute the closed-loop poles.
The closed-loop characteristic polynomial is
% and closed-loop poles
phi = d*x+n*y;
clpoles = roots(phi)
clpoles =
-46.7980
-0.9727 + 0.6271i
-0.9727 - 0.6271i
-0.7071 + 0.7071i
-0.7071 - 0.7071i
We first run the script with
gamma = 4. The closed-loop poles all have negative real parts, and, hence, the closed-loop system is stable.Next, we run the macro with
gamma = 3.5. The script returns4.9995
-0.9590 + 0.5600i
-0.9590 - 0.5600i
-0.7071 + 0.7071i
-0.7071 - 0.7071i
The closed-loop system is unstable. Hence,
gamma has been chosen too small.Note that four of the five closed-loop poles do not change much with
gamma. The fifth pole is very sensitive to changes in gamma.We test this dependence by running the script several times for different values of
gamma without showing the output. Table 1 shows the results. Apparently as gamma decreases the fifth pole crosses over from the left- to the right-half complex plane, but does so through infinity.For the final run we take
gamma = 3.9515, close to the optimal value but such that the closed-loop system is stable. The corresponding numerator polynomial y and the denominator polynomial x of the compensator arex =
1.5e+004 + 1e+004s + 3e+003s^2 + s^3
Table 1. Dependence of the fifth pole on
gammagamma |
the fifth pole |
4 |
46.798 |
3.5 |
4.9995 |
3.75 |
11.369 |
3.875 |
30.297 |
3.9375 |
173.86 |
3.9688 |
127.80 |
3.95 |
315.38 |
3.951 |
3153.8 |
3.9515 |
2987.6 |
In fact, we may cancel the leading coefficients and simplify y and x. This amounts to cancelling the compensator pole-zero pair and eliminating the corresponding closed-loop pole that pass through infinity as
gamma passes through the optimal value.Recalculation of the closed-loop poles after this cancellation confirms that the large closed-loop pole has disappeared. These calculations are performed by typing a few simple command lines
y =
1.3 + 3.8s
x =
5.1 + 3.4s + s^2
phi = d*x+n*y;
clpoles = roots(phi)
clpoles =
-0.9703 + 0.6201i
-0.9703 - 0.6201i
-0.7075 + 0.7072i
-0.7075 - 0.7072i
We type in the command lines
Fig. 3. Magnitude plots of the sensitivity functions
we rearrange the factorization as
L is diagonal in the form
where
and are diagonal nonnegative-definite but not unit matrices. If the factorization is close to noncanonical then L is close to nonsingular. The large numbers disappear.The alternative factorization is obtained by an option in the
spf command. The computation of the compensator is not affected, so we only need to change the code line that contains the spf command toRerunning the script with this modification for
gamma = 3.9515 shows that the large numbers have disappeared. We also see that instead of the large closed-loop poles and corresponding large pole and zero of the compensator we now have a closed-loop pole, compensator pole and zero at 1:y =
1.3 + 5.1s + 3.8s^2
x =
5.1 + 8.4s + 4.4s^2 + s^3
rootsx = roots(x), rootsy = roots(y), clpoles
rootsx =
-1.6787 + 1.5027i
-1.6787 - 1.5027i
-1.0000
rootsy =
-1.0000
-0.3476
clpoles =
-0.7071 + 0.7071i
-0.7071 - 0.7071i
-1.0002
-0.9715 + 0.6197i
-0.9715 - 0.6197i
Cancellation of the common root of y and x leads to the same optimal compensator that was previously obtained:
ans =
5.1 + 3.4s + s^2
ans =
1.3 + 3.8s
This search algorithm has been implemented in the macro
mixeds. Calling the routine in the form gmin = 3.5; gmax = 4; accuracy = 1e-4;executes the search while showing the intermediate
results. The search stops if gmax-gmin is less than the input parameter accuracy. This is the output:
----- -----------
4 stable
3.5 unstable
3.75 unstable
3.875 unstable
3.9375 unstable
3.96875 stable
3.95313 stable
3.94531 unstable
3.94922 unstable
3.95117 stable
3.9502 unstable
3.95068 unstable
3.95093 stable
3.95081 stable
3.95074 stable
Cancel root at -0.999999
y =
1.3 + 3.8s
x =
5.1 + 3.4s + s^2
gopt =
3.9507
Rather than relying on one of the polynomial division routines we write a few dedicated
code lines. Suppose that the numerator y has a root z. Then by polynomial
division we have for the denominator x
where the remainder r is a constant. By substituting s = z we see that actually r = x(z). Expanding the polynomials x and q as
it follows that the quotient q and the remainder r may recursively be computed as
Note that this way of computing r = x(z) is nothing else than Horner's algorithm. If the remainder r is small then we cancel the factor sz. By the same algorithm the factor sz may be cancelled from the numerator y.
These are the necessary code lines. A tolerance
tolcncl is used to test if the remainder is small.% Divide xo(s) by s-z
q = zeros(1,degxo);
q(degxo) = xo(degxo+1);
for j = degxo-1:-1:1
q(j) = xo(j+1)+z*q(j+1);
end
% If the remainder is small then cancel
% the factor s-z both in xo and in yo
if abs(xo(1)+z*q(1)) < tolcncl*norm(xo,1)
xo = q; degxo = degxo-1;
p = zeros(1,degyo);
p(degyo) = yo(degyo+1);
for j = degyo-1:-1:1
p(j) = yo(j+1)+z*p(j+1);
end
yo = p; degyo = degyo-1;
if show
disp(sprintf('\nCancel root at %g\n',z));
end
end
end
includes a four-dimensional optional tolerance parameter
which allows to fine tune the macro. For a description of the various tolerances consult the manual page for
mixeds.