Category A PRACTICAL. APPROACH TO. ROBUSTNESS. ANALYSIS

The #oo missile autopilot was proved to exhibit good robust stability and performance properties in the presence of uncertainties in the aero­dynamic model. The robustness margin obtained when analyzing the robust stability property inside the left half plane is very good, as well as the robustness margin corresponding to the robust stability property inside a truncated sector. The results obtained when defining the per­formance in the frequency domain (subsection 1.3) appear disappointing. However, as indicated in chapter 9, the poor quality of these results is due to the use of classical p tools: it will come out that the results ob­tained when defining the performance in the frequency domain and when using skewed p tools appear very close to those obtained when defining the performance through a truncated sector.

Concerning the methods, see first the table above for a comparison of the estimates of the robustness margins, obtained with the different computational methods. In the case of the missile problem, whose di­mension is not too large, the best results are obtained by combining the upper bound by Zadeh and Desoer and the lower bound by Dailey. The mixed /j upper bound by Fan et al also gives good results.

15 -25 ■ ft.

frequency in racVs

Figure 6.5. Missile autopilot – template for robust performance – nominal sensitivity

function in dashdot line – template of (Balas and Packard, 1992) in dashed line (with

factor 0.9) – new template in solid line.

Robust performance is analyzed in the frequency domain, using the sensitivity function S (see section 2. of chapter 2 for the definition of S). A fictitious performance block is added to the model perturbation, which now contains 4 non repeated real scalars and a full complex block.

Figure 6.6. Missile autopilot – robust performance (P3) – mixed и lower and upper bounds in solid and dashed lines (the two plots nearly coincide) – complex и upper bourn in (lash-clotted ine.

The template on the sensitivity function S is essentially the same as the one in (Balas and Packard, 1992), except that the low frequency performance is relaxed and that the template is multiplied by factor 0.9 (see Figure 6.5). The low frequency performance is relaxed in order to focus on the performance at medium frequencies. On the other hand, the template is multiplied by 0.9, so as to choose the worst allowable performance as uic = 4 rad/s and tr«^ = 0.25s (the bandwidth of the nominal sensitivity function S iscuc = 5 rad/s, thus leading to a nominal closed loop rise time tr w 0.20s).

Figure 6.6 presents the mixed p lower and upper bounds in solid and dashed lines and the complex p upper bound in dash-dotted line. The mixed p lower and upper bounds nearly coincide at all frequencies. The maximal value of the mixed p upper bound is 0.87 at ui = 11.00 rad/s (the result is nearly non conservative). The corresponding uncertainty in the stability derivatives is thus 5/0.87 « 5.8%.

Note finally that a peak at medium frequency is here again obtained (around 11 rad/s). This strongly suggests that the degradation of the stability and performance properties is essentially due to a decrease of the damping ratios of some closed loop poles at medium frequencies.

Robust stability inside a truncated sector is studied. The minimal degree of stability (resp. the minimal damping ratio) is chosen as 0.3 (resp. 0.4). The nominal degree of stability (resp. the nominal damping ratio) is 0.49 (resp. 0.62).

Figure 6.3. Missile autopilot – robust stability inside a truncated sector (P2) – the real /і upper bound by Zadeli and Desoer is represented in dashed line and the real lower bound by Dailey in solid line.

We proceed in the same way as in the previous subsection. The same 5 methods are used. Here again, all these methods give nearly the same result (around 0.12) at ш = 0. We then focus on the real // upper bound

by Zadeh and Desoer (dashed line – see Figure 6.3) and the real lower bound by Dailey (solid line). Note here again the good accuracy of the p interval at nearly all frequencies, and the discontinuity of the p bounds by Zadeh and Desoer and by Dailey at и = 0: p w 0.12 at и = 0, while p « 0.06 at very low frequencies. The real s. s.v. is thus discontinuous at the zero frequency.

Missile autopilot – robust stability inside a truncated sector (P2) – the

classical complex ц upper bound is represented in solid line, the real и upper bound by Jones in dashed line, the mixed ц upper bound by Fan et al in dash-dotted line and the real fi upper bound by Zadeh and Desoer in dotted line.

The maximal value of the upper bound by Zadeh and Desoer is 0.25 at u) = 13.6 rad/s, while the maximal value of the lower bound by Dailey is 0.24 at з = 1 3.5 rad/s. The accuracy of the estimate of the robustness margin is thus very good (less than 5 %). The corresponding uncer­tainties in the stability derivatives are 5/0.25 = 20%. The Hoo missile autopilot consequently presents good robust performance properties in

the presence of uncertainties in the aerodynamic model.

When analyzing robust stability inside the left half plane, two peaks were obtained at ш = 0 (the real ц was discontinuous at this frequency) and at medium frequencies ( a a 1 4.4 sa4/s( . In the same way, when analyzing robust stability inside a truncated sector, two peaks are here again obtained at u> = 0 and at medium frequencies (со и 13.6 rad/s).

As a final point, Figure 6.4 presents the results obtained with the four fj, upper bounds (classical complex /r upper bound in solid line, real ц upper bound by Jones in dashed line, mixed ц upper bound by Fan et al in dash-dotted line, real ц upper bound by Zadeh and Desoer in dotted line). The same comments can be done as for Figure 6.2.

5 methods of chapter 5 are used, namely the mixed /і upper bound by Fan et al (subsection 2.3), the classical complex д upper bound (subsec­
tion 2.2), the real fj, upper bound by Zadeh and Desoer (subsection 1.1), the real Ц upper bound by Jones (subsection 1.3), and the real /л lower bound by Dailey (subsection 1.2). Note first that all these methods give nearly the same result (around 0.13) at u> = 0.

Figure 6.1. Missile autopilot – robust stability (PI) – the real p upper bound b

Zadeh and Desoer is represented in dashed line and the real p lower bound by Dailey in solid line.

Figure 6.1 presents the real л upper bound by Zadeh and Desoer (dashed line) and the real p lower bound by Dailey (solid line). Note the good accuracy of the p interval at nearly all frequencies. The maximal value of the s. s.v. is obtained as 0.13 at ш = 0. The result is non conser­vative. The corresponding uncertainties in the stability derivatives are 5/0.13 « 40%. This means that the closed loop missile is stable despite simultaneous uncertainties of ± 40 % in the four stability derivatives. Note also the discontinuity of the p bounds by Zadeh and Desoer and by Dailey at u> = 0: д и 0.13 at w = 0, while p « 0.05 at very low
frequencies. This means that the exact value of the real s. s.v. is also discontinuous at the zero frequency, as already mentioned in (Ferreres et al., 1996a).

Missile autopilot – robust stability (PI) – the classical complex /і upper

bound is represented m solid line, the real /і upper bound by Jones in dashed line, the mixed Ц upper bound by Fan et al in dasli-dotted line and the real ц upper bound by Zadeh and Desoer in dotted line.

Figure 6.2 presents the results obtained with the four p upper bounds (classical complex /x upper bound in solid line, real p upper bound by Jones in dashed line, mixed p upper bound by Fan et al in dash-dotted line, real p upper bound by Zadeh and Desoer in dotted line). This Figure illustrates the decreasing conservatism of these bounds. The best results are obtained with Zadeh and Desoer’s and Fan et al’s methods. More precisely, better results are obtained with Zadeh and Desoer’s method at low frequencies, equivalent results are obtained at medium frequencies, and better results are finally obtained with Fan et al’s method at high
frequencies. The classical complex // upper bound, the real // upper bound by Jones and the mixed //upper bound by Fan et al are continuous at = = 0 , unlike the real /і upper bound by Zadeh and Desoer.

The aim of this chapter is to evaluate the applicability of the com­putational methods, which were presented in the previous chapter. See (Ferreres et al., 1996a) for a related work in the context of a single-axis and three-axes missile autopilot. See also chapter 5 (section 4.) for a summary of the computational methods.

The first section evaluates the robust stability and performance prop­erties of the longitudinal missile autopilot, robust performance being defined either as robust pole location inside a truncated sector, or as ro­bust shaping of the sensitivity function. The lateral flight control system is analyzed in the same way in the second section. Since this controller was synthesized with a modal approach, robust performance is defined in this context as robust pole location inside a truncated sector.

1. THE MISSILE AUTOPILOT

The aim of this section is to analyze the local stability and performance properties of the #oo autopilot in the presence of parametric uncertain­ties in the 4 stability derivatives Ma , Ms, Za, Z$. The high frequency bending mode is not taken into account (see chapter 9), so that the model perturbation only contains 4 non repeated real scalars. The weights in the stability derivatives are chosen as 5 %. See chapters 2 (subsection 2.2) and 4 (section 1.) for details concerning the description of the linearized missile model and the building of the interconnection structure.

More sophisticated methods can be used to move the poles of A + BA*RC through the imaginary axis. In the spirit of (Magni and Doll, 1997), a solution is to introduce an additional model perturbation Ar. The aim is then to find the minimal size model perturbation Дд + Ад, which moves one pole of A + B(A*R -t – Ад)С through the imaginary axis. This problem can be easily recast as a simple Linear Programming problem.

First remember that the model perturbation is diagonal:

Ад = diag(6iIkl, S2lk2i – ■ ■ ,firlkr)

Ад = diag(6lIkl,6*2Ik2,…,6*rIkr) (5.26)

Let A о the closed loop pole, which is the closest to the imaginary axis. Using well-known results on eigenvalues derivatives (Kato, 1980), and assuming that the magnitude of the additional parametric uncertainties Si is sufficiently small, the first order variation of Ao is computed as an affine function of the Sfs:

Г

*Ao = ^2

1=1

In order to impose that Ao is moved onto the imaginary axis, one imposes 7£e(Ao + <5Ao) = 0. Remember then that the infinity norm of the model perturbation A*r + Ад is to be minimized, so that the problem reduces to the minimization of scalar 7 under the constraints:

-7 < s* + й <_7, i = l,…,r TZe(Xo + ^2i a(Si) = 0

It may be necessary in practice to modify the above problem, when the eigenvalue Ao is not sufficiently close to the imaginary axis. In this case indeed, the accuracy of the first order development may not be sufficient, to directly move the pole onto the imaginary axis. A solution consists in partitioning the real segment between Ao and the imaginary axis, and to iteratively perform the migration on each sub-segment. At the end of this process, the method of the previous subsection is applied, to ensure an exact pole placement on the imaginary axis.

Remarks:

(/’) In the same way as in (Magni and Doll,

of the model perturbation A*R + Дд could be minimized, instead of the infinity norm. The results are however more conservative than those ob­tained with the LP method above, since the s. s.v. handles the infinity norm.

(//) It would be possible to impose that the imaginary axis is crossed at a given point juj, in order to compute a lower bound of ^AR(M(juj)) at a given frequency u>. In this case, two equality constraints (7?.e(Ao + £Ao) = 0 and 2m(Ao + <SAo) = u>) instead of a single one (7£e(Ao + £Ao) = 0) are to be considered. Nevertheless, the approach becomes numerically sensitive, since it appears difficult to move the closed loop poles onto

some regions of the imaginary axis (especially at very low and very high frequencies).

(/77) It is worth emphasizing that the aim is not to compute a lower bound of n at a given frequency, but to detect the peak values on the H plot, in order to directly compute a lower bound of the maximal real s. s.v. over the frequency range. Assume that a ц lower bound was com­puted for the regularized ц problem at a frequency wo, which is far from the critical frequency ш*. Using especially the LP method above, the imaginary axis can be crossed at a frequency <I>o which is very close to to* (see Figure 5.4). Indeed, the imaginary axis is not constrained to be crossed at a given point jco, and the norm of the model perturba­tion Дд + Дд is to be minimized. As a consequence, it is generally observed that the imaginary axis is crossed at frequencies which corres­

pond to the peak values on the p plot (remember that p is homogeneous to the inverse of the robustness margin, so that the peak values on the p plot correspond to model perturbations which are of minimal size). The method consequently detects the critical frequencies on the p plot: see also chapter 10 (subsection 6.4).

(iv) The following table summarizes the three methods presented above.

4. SUMMARY

A great deal of work was devoted to the problem of computing the s. s.v.. Existing computational algorithms can be sorted following various cri­teria:

■ Nature of the structured model perturbation: two large categories of methods can be considered, corresponding to the special case of a real model perturbation and to the general case of a mixed model perturbation [7].

■ Nature of the result: the algorithm provides the exact value of p, a p lower bound or a p upper bound.

■ Computational requirement: in polynomial time algorithms (resp. ex­ponential time algorithms), the computational amount is a polynomial (resp. exponential) function of the size of the problem. Nevertheless, it must be noted that the computational amount can also largely differ for two methods inside the same category. Remember finally that an algorithm, which computes the exact value of p, is necessarily exponential time, since this problem is NP hard.

As an illustration, the characteristics of the computational methods, which were presented in this chapter, are summarized in the above table. In practice, a trade-off is to be achieved between the accuracy of the res­ult, the computational requirement and the size of the problem. As a simple example, the exact value of jj, can be computed in the case of small size problems, without an excessive amount of computation. How­ever, computing the exact value of /x in the case of medium or large size problems would require an excessive computational amount.

More generally, in the case of small dimension problems, exponen­tial time algorithms can be used. However, it is necessary for large size problems to compute an interval of the s. s.v. j with polynomial-time al­gorithms. Even if the gap between the ц bounds can not be guaranteed a priori, good results can be nevertheless obtained in realistic examples, as illustrated in chapter 6. Finally, concerning more specifically the com­putational methods, the following two points are recalled:

■ The mixed fx upper bound by Fan et al is obtained as the solution of a quasi-convex optimization problem, namely an LMI problem.

■ Conversely, the mixed /x lower bound by Young and Doyle is ob­tained as the solution of a non convex optimization problem: the idea is more precisely to obtain the lower bound as the limit of a fixed point iteration Xk+i = f(xic)- However, the associated power algorithm is not guaranteed to converge, and the final result depends on the initialization of the fixed point iteration. Nevertheless, good results are generally obtained, except in the case of a purely real model perturbation. In this specific context, it is more interesting to use Dailey’s method if the size of the problem is sufficiently small (this method is indeed exponential-time). Otherwise, section 3. proposes a polynomial-time method, which directly computes a lower bound of the maximal s. s.v. over the frequency range.

Nevertheless, it can be remarked that the real part A*R of A* is prob­ably a good initial guess in the search of a model perturbation Ад, which renders the matrix I — М()ш)Дд singular. There is however a technical difficulty, which is to decide when a matrix X is singular. An obvious solution is to compute the magnitude of the determinant det(X), the minimal singular value g_{X), or the condition number and to de­cide that X is singular when one of these quantities is lower than a given value. The value of this threshold is however difficult to determine, and it may depend on the problem data.

If the aim is to compute a lower bound of the maximal s. s.v. over the frequency range, rather than a lower bound of n&R(M(ju>)) at a fixed frequency w, a natural solution is to remember thatMO’w) is the value of the transfer matrix M(s) at s = ju. Let (А, В, C, 0) a state-space repres­entation of M(s), which is here again assumed to be strictly proper just for the ease of notation. It can be expected that one pole of A + BA*RC is close to the point ju of the imaginary axis, if a lower bound of the regularized s. s.v. fj,^(H(ju>)) was computed.

A very simple solution is thus to increase the size of the model per­turbation until one pole of A + aBA*RC crosses the imaginary axis at the point jCj (for a > 1). Since one pole of A + BA*RC is expected to be close to the point jio of the imaginary axis, the value of a, for which a pole of A + aBA*RC crosses the imaginary axis, is possibly close to 1.

In an obvious way, a lower bound of R(M(ju)) was computed in­stead of a lower bound of цаR(M(ju)). But more importantly, a model perturbation aA*R was obtained, which brings one closed loop pole on the imaginary axis. The inverse of the size of otA*R is thus a lower bound of the maximal s. s.v. over the frequency range.

The method is thus the following. A frequency gridding is first chosen, as usually in ц analysis. A regularized fi lower bound is computed at each point of this gridding, and the real part A*R of the augmented model per­turbation, which is provided by the power algorithm, is extracted. The value of a is then increased from the initial value of 1 until one pole of A + aBARC crosses the imaginary axis.

At each point of the frequency gridding, a lower bound of the maximal s. s. v. over the frequency range was thus obtained as the inverse of the size of the model perturbation aAR. The best estimate of the maximal s. s.v. over the frequency range is finally chosen as the highest value of this lower bound over the frequency gridding.

Remarks:

(i) The problem is slightly more complex in practice. There exist indeed two ways of optimizing a. In the above approach, a is increased from its nominal value of 1, until the target pole crosses the imaginary axis (.i. e. the pole of A + BA*rC, which is the closest to the point ju of the imaginary axis, if a lower bound of the regularized s. s.v. /хд (H(jw)) was computed). It is however possible that some other poles of A 4- BA*RC were found to be strictly unstable. In this case, an alternative is to decrease a until all unstable poles cross the imaginary axis. This ap­proach provides a smaller destabilizing perturbation, but the associated frequency w is no longer guaranteed to be close to the initial value oj.

(.ii) The method can be readily extended to the problem of robust stabil­ity inside a region fi of the complex plane (typically a truncated sector).

A well-known solution for improving the convergence properties of the power algorithm by (Young and Doyle, 1990) is to regularize the real p problem, by adding a small amount e of complex uncertainty to each real uncertainty (Packard and Pandey, 1993). Let Ад and M(juj) the original data of the real ju problem, which thus consists in computing (a lower bound of) paR(M(ju>)). A model perturbation Ac is introduced, with the same structure as Ад, except that the real scalars become complex.

Let then the augmented model perturbation A = diag(An, Ac) and let:

The issue is now to compute a lower bound of /ід (H(ju)). The higher the value of e, the larger the amount of complex uncertainties, and the better the convergence properties of the power algorithm, which is applied to this regularized problem. e is classically chosen as 5 % or 10 %.

However, the lower bound obtained for this regularized p problem is not a lower bound for the original real p problem. The power algorithm provides indeed a critical model perturbation A* = diag(AR, A^), which renders the matrix I – H(juj)A* singular. No model perturbation Ад was found, which would render the matrix I — M{ju))A. R singular. It is especially worth emphasizing that the matrix I — M(ju})AR is not a priori singular.

3.1 BACKGROUND

The aim of this section is to compute a reliable upper bound of the robustness margin (i. e. a lower bound of the maximal s. s.v. over the fre­quency range) in the context of a purely real model perturbation, which contains a large number of parametric uncertainties (Ferreres and Bian – nic, 1998b). Note indeed that the primary aim is to compute an interval of the robustness margin, rather than an interval of the s. s.v. n(M(juj)) as a function of ui.

This is not an easy problem, since the real ) lower bound of (Dailey, 1990) is exponential time. On the other hand, the mixed д lower bound of (Young and Doyle, 1990) is polynomial time, but the convergence prop­erties of the power algorithm are very poor in the context of a purely real model perturbation. The result provided by this power algorithm is consequently not reliable in this specific context.

A state-space approach was proposed in (Magni and Doll, 1997): con­sider the interconnection structureM{s) – Д, and let (А, В, C, 0) a state – space representation ofM(s), which is assumed to be strictly proper just for the ease of notation. The idea is to interpret the model perturbation Д as a fictitious feedback gain which moves the poles of the closed loop, whose state matrix is A + ВАС, from the left half plane through the imaginary axis. The norm of Д is to be minimized during the process of migration of the poles towards the imaginary axis.

The algorithm provides a model perturbation A*, which brings one closed loop pole on the imaginary axis at jui*. An upper bound of the robustness margin is thus obtained as the size of A*, and its inverse is a lower bound of the s. s.v. )(M (M*) ) (and thus a lower bound of the maximal s. s.v. over the frequency range).

Good results are generally obtained with this technique. There are however some technical difficulties. A key issue is especially to choose

which poles of the nominal closed loop (i. e. which poles of the state matrix A) are to be moved towards the imaginary axis.

The exact value of the mixed s. s.v. can be obtained as the solution of a non convex optimization problem, whose global maximum coincides with the exact value of p. The idea of the power algorithms by (Packard et al., 1988; Young and Doyle, 1990) is to write the necessary conditions of optimality as f (x) = x. Power iterations x(k + 1) = f ‘(x(k)) are then used to asymptotically solve f (x) = x. Note that this algorithm only provides (at least a priori) a local maximum of the non convex optimiz­ation problem, i. e. a p lower bound.

The optimization problem is first detailed in the following two propos­itions. As a preliminary, the real spectral radius of a complex matrix A is defined as the magnitude of the largest real eigenvalue of A:

Pr(A) = sup(|A| / det(A — XI) = 0 and А є R) (5.18)

Pr(A) is zero if matrix A has no real eigenvalue.

PROPOSITION 2..2

p(M) = гам pr(AM)

Proof: a structured model perturbation Д is searched, for which there exists a non-zero vector x satisfying:

(I – AM)x = 0

Д can be rewritten as:

so that /? is an eigenvalue of matrix A°M. Since the smallest size destabilizing model perturbation is searched, the maximal value of (3 is thus to be obtained, i. e. the maximal value of the real spectral radius of matrix AqM over Д° Є BA.

Proposition 2..3 Let (see equation (5.10)):

Q = {A / SI є [-1,1], 6f6Cj = 1 and AcqHAcq = 1} (5.23)

Then:

ц{М) = шxpR(QM)

It is worth emphasizing the coherence between the definition of Q and the counterexamples of (Ackermann, 1992; Holohan and Safonov, 1993). When considering complex uncertainties, it suffices to search the destabil­izing model perturbation over the unit sphere. However, as proved in (Ackermann, 1992; Holohan and Safonov, 1993), no assumption can be made about the real perturbations <5[, which are thus simply assumed to lie inside their unit ball (i. e. S Є [—1,1]).

The issue is now to solve the optimization problem, noting as a pre­liminary that a local minimum will be necessarily obtained in the general case (since the problem is non convex) and that the computational bur­den should remain minimal. As said above, a classical solution (Packard et al., 1988; Young and Doyle, 1990) is to rewrite the necessary con­ditions of optimality under the form f (x) = x, where f is a vectorial function of vector x. A fixed point method finds a limit x* of the series xk+l = f(xk)’

When the series converges, x* satisfies indeed the necessary condi­tions of optimality and a local minimum has been obtained, i. e. a /і lower bound. Nevertheless, the value of x* depends on the initial condi­tion жо of the series, because of the non convexity of the problem. On

the other hand, the series does not necessarily converge: a limit-cycle may especially appear inside the power algorithm.

In practice, the algorithm by (Young and Doyle, 1990) generally presents good convergence properties. The computational burden is low and the quality of the ц lower bound is usually good. However, it is worth em­phasizing that this power algorithm exhibits poor convergence properties in the case of a purely real model perturbation (see below). Note finally that the power algorithm by (Young and Doyle, 1990) is not detailed here, since a generalization of this power algorithm to the skewed ц problem will be presented in section 3. of chapter 8.