175
pages

Voir plus
Voir moins

Mastering Complexity:

an Overview

erospace systems are certainly part of those of which the complexity

Claude Barrouil Ahas been dramatically and continuously increasing for several decades.

Scientific Director Not only are the vehicles increasingly more sophisticated but, since they are

of the Information Processing

interacting increasingly with other so-called “intelligent” systems, human or and System Branch

machines, the overall system is even more complex.

Here, the word “complexity” is used in the common sense: it just means

“too big to fit in the head of a human being”. We undertake projects that

are always aimed beyond the previous ones and we now face challenges

that cannot be addressed without the aid of information technologies, which

extend human capabilities. We are at a point where the matter is not “simply”

to invent new systems, but to invent codes that will help to invent new

complex systems.

This fourth Aerospace Lab issue is dedicated to various techniques used

to address some of the most difficult issues in complex aerospace system

design. It is not focused on a sub-class of aspects; it is a gallery of articles

that will encompass the main issues to be addressed when considering

advanced aerospace systems. There are two groups of articles: those related

to embedded system concepts and those related to concept design aid.

Aerospace System

Embedded System System Design

Automation Operator Engineering Evaluation

Control Decision

Articles Articles Articles Articles Articles

[1] [5] [8] [10] [13]

[2] [6] [9] [11] [14]

[3] [7] [12]

[4]

Issue 4 - May 2012 - Mastering Complexity: an Overview

AL04-00 1Embedded system information processing components Motion control is obtained by the cooperation between three basic

functions: sensing, state estimation and control signal computation.

Flexible aircraft control Likewise, robot “intelligence” is based on the cooperation between

perception, situation assessment and decision. Research in artificial

Modern aircraft (A/C) flight qualities result from the permanent inte- intelligence has been underway for more than 30 years and we are

raction between aeroelastic phenomena and close-to-actuator control still far from being able to make robots as smart as C-3PO or R2D2.

laws, even in manually piloted mode. New materials allow A/C to be However, knowledge has been acquired and efficient methods exist

lighter and larger; in return they offer multiple structure flexible modes for decision making in realistic contexts, in particular on-line

decisionwhich slowly slip as fuel weight decreases and as the flight point making under uncertainty and partial observability, as presented in [5].

changes. It is essential for the control laws to be robust to these varia- For a drone scenario, a scene understanding layer must provide the

tions, i.e., that they provide good flight qualities in all configurations. decision layer with situation assessment information. Which objects

The Theory of Control offers methods for computing robust laws but are in the scene? Where? What is the relative location of the drone with

a naïve approach fails with realistic flexible aircraft, since their model respect to the local scene? As explained in [6] scene understanding

dimension may be as large as several hundred. Recent advances in is a very difficult task. As it is unrealistic to describe a priori all of the

control law synthesis for high-dimensioned systems are presented in objects that may be encountered, learning techniques must be

consi[1]. Robust control laws are optimized given an A/C nominal model dered. How what is expected to be encountered should be encoded,

and an associated model of uncertainty. Of course, the smaller the depending on the sensors used, is also a very hard problem.

uncertainty, the more performing the control is. Therefore, the quality

of the nominal model is essential and this is the flight test purpose, to Fusing heterogeneous information

collect data to identify a precise A/C dynamic model, including

structure flexible modes. This task is complicated by the model dimension The "system of systems" military concept refers to a set of

autonoand by the need to process flight test data very quickly, since the mous systems that must coordinate in order to perform a given action,

flight test campaigns are expensive. [2] gives an overview of the most none of these being able to perform it alone. Situation assessment is

efficient techniques. However, the model resulting from flight test data one of such actions and requires information provided by several types

processing is generally of unnecessarily large dimension for robust of sources to be merged, some of them being human. In the future,

law synthesis purposes and model reduction must be considered. UAVs will be part of systems of systems. They must thus participate

The Linear Fractional Representation (LFR) presented in [3] is espe- in situation assessment based on the signal that they collect with their

cially well suited for reducing the model order. own sensor, or with the sensors of the other UAVs, as well as on high

level information provided by humans. [7] addresses the question of

Close-to-environment flight control fusing human reports for intelligence purposes. It shows a

methodology that has proven to be compliant with NATO recommendations.

For most existing automatic aerial platforms, the navigation task relies

only on proprioceptive sensors (inertial, pressure, pitot sensors) pos- Evolution of the pilot role

sibly coupled with GPS limiting their cases of use to obstacle-free

areas. Flying low among obstacles, possibly without GPS, can be Few systems are really “autonomous”. Except for “fire and forget”

addressed using optical sensors delivering measurements relative to weapons, there is(are) always a (or several) human operator(s)

the surrounding environment. As shown in [4], managing the platform somewhere in the loop. Drones and satellites are monitored from

from this kind of sensors can take on several forms, from designing control stations. Increasingly more pilot assistance will be introduced

piloting or guidance laws compatible with non-metric low level visual in the cockpits of aircraft and helicopters: the role of the pilot is

chanmeasurements, to inferring 3D information or GPS-like measure- ging and it can be foreseen that at some future time transport aircraft

ments by computer vision and scene understanding. will basically not differ from drones, except that they will have an

operator on-board (maybe). Then, two difficult problems arise.

Closed-loop on-line decision making

First, since the automated systems are certified to perform safe

maBefore the numeric age, only analogue electronics or electromecha- nagement and since the operator is licensed to be competent, how

nical devices could implement control law correctors fed with conti- should a conflict between human and machine be managed? This

nuous signals and Boolean information. Then, computers allowed problem analysis and modeling is presented in [8].

more general symbolic data to be dealt with also: on-board

automa(1)ted reasoning became possible, “artificial intelligence” would follow Second, since the automated system makes decisions without

referand robots could be imagined. Note, however, that torpedoes were ring every time to the operator, for the sake of reducing the workload

already simple, but genuine, operational robots. Drones and missiles for instance, the operator may lose the sensation that he/she is still in

are the target applications in the aerospace field, but it is clear that control of the system. He/she could feel disconnected, loose his/her

classic civil A/C and helicopters will benefit from intelligence capabili- situational awareness and the result may be catastrophic. “Sense of

ties, by even more efficient new pilot assistance. Automated decision agency” (= feeling of being an agent) is a new formalism, presented

is a key issue for space systems. in [9], suited for modeling and analyzing this problem.

(1) Of course, this sort of “intelligence” should not be compared to human, or animal, intelligence.

Issue 4 - May 2012 - Mastering Complexity: an Overview

AL04-00 2Designing aerospace systems Multidisciplinary Design Optimization

So far, we have considered issues about aerospace system real- Vehicle architecture design benefits from the progress of

multidiscitime information processing. We now address issues about how plinary design optimization techniques (MDO). It allows fast

dimenaerospace systems, software and avionics should be designed sioning of entirely new concepts or, simply, evaluation of the

potention the one hand, and about how vehicle architecture should be alities of a new subsystem concept by tuning the other subsystems

designed on the other hand. accordingly. Although several generic frameworks exist, achieving a

realistic MDO environment requires the merging of several

compeFormal methods for software verification tences and practical experience, as explained in [12].

The software volume on an A300 (first flight: 1972) was about 2 Large distributed simulation techniques

million lines of code. For an A380 (first flight: 2005) it exceeds

100 million lines of code. It is clear that correctness cannot be Formal property assessment is not always possible. In those cases,

checked without computer assistance. This is not specific to the use of simulation is necessary. However, Monte-Carlo techniques

aerospace; other sorts of critical systems (automobiles, trains are not realistic when a combinatorial explosion of cases is to be

and nuclear power plants) encounter the same concerns. The explored and/or the system is time-dependent. Article [13] presents

research effort, summarized in [10], is aimed at developing ge- various techniques for coping with these problems, in particular when

neric tools based on formal methods, to describe the system at the architecture of the system can be formally described. When the

several level of abstraction, including the safety requirements, system is too large, or when it includes heterogeneous sub-systems,

and to perform automated verification. typically when dealing with systems of systems, simulation remains

the only means to assess the emerging behavior resulting from the

Avionics challenges interaction of several entities. Such simulation should be evolutionary,

i.e., easy to modify when components are modified or, even, when

Avionic system implementation also dramatically changed components are added. Distributed simulation is a key technique,

from the 80’s, when sub-systems had their own dedicated [14] presents operational methods based on the HLA standard.

computer and interacted through dedicated links. Modern

design is based on Integrated Modular Avionics (IMA), which 14 articles are certainly much too few to properly deal with the entire

allows computing resource sharing with no-interference insu- problem addressed above, but a number of references will allow the

rance. As shown in [11], IMA provides flexibility, which should reader to go further. However, important problems have not been

adlead to reconfiguration capabilities. [11] also addresses the dressed. For instance: automated system certification, fault detection,

problem of using a new multi-core generation of processors isolation and recovery and system of systems formal modeling. This

in future avionics. is a good reason to think of a future Aerospace Lab issue

References

[1] C. Roos, C. Döll, J.M. Biannic - Flight Control Laws : Recent Advances in the Evaluation of their Robustness Properties. Aerospace Lab Issue 4,

May 2012.

[2] A. Bucharles, C. Cumer, G. hardier, B. Jacquier, A. Janot, T. Le Moing, C. Seren, C. Toussaint, P. Vacher - An Overview of Relevant Issues for Aircraft

Model Identification. Aerospace Lab Issue 4, May 2012.

[3] C. Poussot-Vassal, F. Demourant - Dynamical Medium(Large)-Scale Model Reduction and Interpolation with Application on Aircraft Systems.

Aerospace Lab Issue 4, May 2012.

[4] M. Sanfourche, J. Delaune, G. Le Besnerais, H. de Plinval, J. Israel, P. Cornic, A. Treil, Y. Watanabe, A. Plyer - Perception for UAV: Vision-Based

Navigation and Environment Modeling. Aerospace Lab Issue 4, May 2012.

[5] G. Verfaillie, C. Pralet, V. Vidal, F. Teichteil, G. Infantes, C. Lesire - Synthesis of Plans or Policies for Controlling Dynamic Systems.

Aerospace Lab Issue 4, May 2012.

[6] S. Herbin, F. Champagnat, J. Israel, F. Janez, B. Le Saux, V. Leung, A. Michel - Scene Understanding from Aerospace Sensors: what can be Expected

? Aerospace Lab Issue 4, May 2012.

[7] J. Besombes, L. Cholvy, V. Dragos - A Semantic-Based Model to Assess Information for Intelligence. Aerospace Lab Issue 4, May 2012.

[8] C. Tessier, F. Dehais - Authority Management and Conflict Solving in Human-Machine Systems. Aerospace Lab Issue 4, May 2012.

[9] B. Berberian, P. Le Blaye, N. Maille, J.-C. Sarrazin - Sense of Control in Supervision Tasks of Automated Systems. Aerospace Lab Issue 4, May 2012.

[10] V. Wiels, R. Delmas, D. Doose, P.-L. Garoche, J. Cazin, G. Durrieu - Formal verification of Critical Aerospace Software. Aerospace Lab Issue 4,

May 2012.

[11] P. Bieber, F. Boniol, M. Boyer, E. Noulard, C. Pagetti - New Challenges for Future Avionic Architectures. Aerospace Lab Issue 4, May 2012.

[12] S. Defoort, M. Balesdent, P. Klotz, P. Schmollgruber, J. Morio, J. Hermetz, C. Blondeau, G. Carrier, N.Bérend - Multidisciplinary Aerospace System

Design: Principes, Issues and Onera Experience. Aerospace Lab Issue 4, May 2012.

[13] J. Morio, H. Piet-lahanier, F. Poirion, J. Marzat, C. Seren, S. Bertrand, Q. Brucy, R. Kervarc, S. merit, P. Ducourtieux, R. Pastel - An Overview of

Probabilistic Performance Analysis Methods for Large Scale and Time-Dependent Systems. Aerospace Lab Issue 4, May 2012.

[14] P. Carle, R. Kervarc, R. Cuisinier, N. Huynh, J. Bedouët, T. Rivière, E. Noulard - Simulation of Systems of Systems. Aerospace Lab Issue 4, May 2012.

Issue 4 - May 2012 - Mastering Complexity: an Overview

AL04-00 3Mastering Complexity

C. Roos, C. Döll, Flight Control Laws: J.-M. Biannic

(Onera) Recent Advances in the Evaluation

E-mail: clement.roos@onera.fr

of their Robustness Properties

his paper reviews a set of robustness analysis tools developed by the authors Tduring the last decade to evaluate the robustness properties of high-dimensional

closed-loop plants subject to numerous time-invariant uncertainties. These tools are

used to compute both upper and lower bounds on the robust stability margin, the

worst-case H performance level, as well as the traditional gain, phase, modulus and

∞

time-delay margins. The key idea is to solve the problem on just a coarse frequency

grid and to perform a fast validation on the whole frequency range, which results in

guaranteed but conservative bounds on the aforementioned quantities. Some heuristics

are then applied to determine a set of worst-case parametric configurations leading to

over-optimistic bounds. A branch and bound scheme is finally implemented, so as to

tighten these bounds with the desired accuracy, while still guaranteeing a reasonable

computational complexity. The proposed algorithms are successfully assessed on a

challenging real-world application: a flight control law validation problem.

Introduction The paper is organized as follows. The problem is first stated in the

next section, where all of the key ingredients for µ-analysis are

recalDespite recent progress in computer-aided control design techniques, led. The robustness margins computation section is then devoted

the development of flight control laws remains a challenging task. to the characterization of (skew-)µ upper bounds (with a particular

Even the most sophisticated approaches are still based on simplified emphasis on computational aspects), thanks to which guaranteed

models and fail to take all of the requirements into account during robustness margins are obtained. Next, in the worst case analysis

the first design phase. As a result, the validation process remains section, computational approaches are proposed to determine some

a major and possibly time-consuming issue. There is consequently enhanced (skew-)µ lower bounds, which are used not only to

quanan obvious need for robustness analysis tools which can be used tify the conservatism, but also to identify some worst-case

configurato perform a reliable and fast preliminary evaluation of the designed tions useful for design purposes. Extensions of the above results are

controllers before the certification phase. described in tools extensions section to compute worst-case gain,

phase, modulus and delay margins. A branch and bound scheme is

Among robustness analysis techniques (Lyapunov-based, SOS- also proposed in this section in order to reduce the conservatism.

based [2], IQC-based [14], …), µ-analysis is now considered as The algorithms are finally illustrated on a challenging aircraft control

a classical and close to mature approach, which has proved to be application in the application to flight control laws validation section.

useful in many applications (see e.g. [5] and included references). Please note that this is a review paper, which presents the

contribuNevertheless, some technical difficulties remain in specific fields tions of the systems control group of Onera to the field of µ-analysis.

such as robust stability and performance analysis of high-order plants Some sections are thus covered quite briefly, but numerous

refeinvolving largely repeated uncertainties. Indeed, despite recent achie- rences are provided throughout the paper.

vements [19], [20], enhancements are still required to further reduce

the conservatism while limiting the required computational time. In

this context, the objective of this paper is to describe a set of algo- Problem statement and motivations

rithms and tools thanks to which robustness analysis becomes more

reliable and less time-consuming. The ultimate goal is to reduce the Let us consider the standard interconnections of figure 1. M(s) is

number of iterations in a control design process and to avoid expen- a stable real-valued linear time-invariant (LTI) plant representing the

sive Monte Carlo simulation campaigns. nominal closed-loop system. is a block-diagonal time-invariant

Issue 4 - May 2012 - Flight Control Laws: Recent Advances in the Evaluation of their Properties

AL04-01 1operator, which contains all model uncertainties. For the sake of sim- µωMj ( ) provides a guaranteed but conservative value of the ( )∆ i

plicity (see nevertheless remark 1), only parametric uncertainties are robustness margin when Ω is restricted to the single frequency ω ,

i

considered in this paper, which means that is a matrix of the form: while a lower bound µωMj ( ) , usually associated with a worst-( )i∆ ∆case parametric configuration , leads to an over-optimistic value.

∆ = diag δδII, , (1)( )1 n Nn1 N

Computing these bounds over the whole frequency range Ω is a

chalwhere the real scalars are said to be repeated if n > 1. Let lenging problem with an infinite number of both frequency-domain

i iN

nn= . The set of real nn matrices with the same structure constraints and optimization variables. It is usually solved on a finite ∑ kk=1

ωas in (1) is denoted by , and B = ∆∈ ∆ :1σ ∆ < , where frequency grid ( ) and an estimate of the robustness margin { ( ) } i∆ iM∈ 1,[ ]

σ (•) denotes the largest singular value. is then obtained as:

11

≤≤k (4)

max

max µωMj ( ) max µωMj ( ){ ( )} { ( )}∆ i i∆iM∈[1, ] iM∈[1, ]

However, a crucial problem appears in this procedure: the grid must

contain the most critical frequency point for which the maximal value

of µωMj ( ) is reached. If not, the upper bound on k can be ( )∆ i max

very poor, notably in the case of flexible systems, whose plot often

exhibits very high and narrow peaks. Even worse, the lower bound

M(s) M(s)

can be over-evaluated, i.e. be larger than the real value of k .

Unformaxe y

tunately, the aforementioned critical frequency is usually unknown.

The same difficulty arises when worst-case H performance is

consi∞

Figure 1 - Standard interconnections for robust stability (left) and robust dered. Problem 2 is indeed equivalent to a specific skew- problem

performance (right) analysis [4]. Similarly to problem 1, it is thus commonly solved on a frequency

grid using polynomial-time algorithms, leading to both lower [16] and

Two main issues are addressed in this paper: robust stability and (possibly under-estimated) upper [8], [9] bounds on .

max

worst-case H performance.

∞

To overcome the above difficulty, [22] suggests to consider

freProblem 1 (robust stability) Compute the maximum value k such quency as an additional parametric uncertainty, but this strategy

max

that the interconnection of figure 1 (left) is stable ∀∆ ∈k B , usually leads to a computational burden when applied to high-di-max ∆

as well as a destabilizing perturbation ∆∈ ∆ such that σ ∆=k . mensional systems. In this context, some alternative methods are ( ) max

proposed in this paper to compute both tight and reliable bounds on

Problem 2 (worst-case H performance) Assuming that the intercon- either k or :

∞ max max

nection of figure 1 (right) is stable ∀∆ ∈B , compute the highest • Either a µ or a skew-µ upper bound is first computed at some ∆

value of the H norm of the transfer matrix F Ms ( ),∆ from e frequency, for which nothing has been assessed yet. This bound is ( )umax ∞

to y when takes all possible values in B , as well as the corres- slightly increased, and a frequency interval on which it remains valid ∆

ponding value ∆ of ∆. is computed. Such a strategy is repeated until the whole frequency

range has been investigated, leading to either a lower bound on k

max

The most efficient technique to answer these two problems is cer- or an upper bound on . The latter is guaranteed over the whole

max

tainly µ-analysis [3]. This is especially true when high-dimensional frequency range, and not only on a frequency grid as is the case of

systems are considered. The underlying theory [17], [5] is not detai- most existing methods (see the Robustness margins computation

led in this paper but a few useful definitions are recalled below. section);

• Some heuristics are then proposed in the Worst case analysis

Definition 1 Let ω be a given frequency. If no matrix ∆∈ ∆ makes section to determine a worst-case parametric configuration ∆ , such

i

I −∆Mjω singular, then the structured singular value (s.s.v.) that the interconnection between M(s) and ∆ has an eigenvalue on ( )i

µωMj ( ) is defined as µωMj ( ) = 0 . Otherwise: the imaginary axis. Either an upper bound on k or a lower bound on ( ) ( )∆ i ∆ i max

is thus obtained. Unlike in most existing methods, frequency is

max

an optimization parameter, which is used to detect critical frequencies

µωMj = (2)( ( ))∆ i and usually leads to more accurate bounds;

−1 • A branch and bound algorithm is finally described in the Accura-min k ∈ℜ : ∃∆ ∈kB ,det I −Mjω ∆ = 0( ( ) ){ }+ ∆ i cy improvements section. It can be used to compute bounds with the

The robustness margin k is then obtained as the inverse of the desired accuracy and at a reasonable computational cost.

max

maximal value µωMj ( ) over the frequency range of interest Ω ( )∆ i

(usually equal to ℜ ): Note also that extensions of the aforementioned methods are pro-+

posed in the Unstructured margins section to solve the worst-case

−1 unstructured margins problem recalled below.k = max µωMj { ( ( ))} (3)max ∆ iω ∈Ωi

Problem 3 (worst-case unstructured margins) With reference to

fiThe exact computation of µωMj ( ) is known to be NP hard in gure 1 (left) and assuming that k > 1, compute the smallest values ( )∆ i max

the general case, but both upper [25] and lower [24], [21] bounds of the gain, phase, modulus and time-delay margins when ∆ takes all

can be determined using polynomial-time algorithms. An upper bound possible values in B .

∆

Issue 4 - May 2012 - Flight Control Laws: Recent Advances in the Evaluation of their Properties

AL04-01 2

∀ω ∈I ωRobustness margins computation Then condition (6) holds true ( ) where:i

I (ω )= [,ω+ δ ω+ δ ] (8)ii ωωii− i+

Computation of a guaranteed stability margin

Proof Let

The classical way [25] to compute an upper bound on µωMj ( ) Hj ()ω =( )∆

for a given frequency point ω = ω is to introduce two scaling ma- −1i 11D ω Mj ω + ω D ω( ) ( ( )) ( )−−i i itrices D(ω ) and G(ω ), which belong to specific sets D and G reflec- 44F ω − jG ωω Fi i ( ) ( ) ( )i i iβting the block-diagonal structure and the real/complex nature of ∆.i

I ωProposition 1 Let be a positive scalar. If there are some scaling The bounds defining ( ) are obtained by searching for both positive ii

*matrices D()ω ∈D and G()ω ∈G such that: and negative ω of smallest magnitude such that I −Hj ()ωωHj () i i

*

−1 becomes singular, i.e. det I −=Hj ()ωωHj () 0 . A state-space ( ) 11 D(ω )Mj ( ωω)D( )−−i ii

44 σωF ( ) −≤jG(ω ) F (ω ) 1 (5) representation (A ,B ,C ,D ) of H(s) is given by equation (7). A i i i H H H H β *

i state-space representation (A ,B ,C ,D ) of I −Hj ()ωωHj () is X X X X

then given by:

2

where F ωω=IG+ and:( ) ( )ii AB0 − HH= , =XX ** *−−CC A C Dnn× * H H H • D= DC∈ ,D =D > 0:∀∆∈ ∆,D∆ =∆D { } *C = DC B , D =I−D D

X H X H H

nn× * *• G= GC∈ ,:G =G ∀∆∈ ∆,G∆ =∆ G { }

Some standard manipulations finally conclude the proof:

µ Mjω ≤ βthen ( ( )) .∆ i

*

det I −HjωωHj = 0( ) ( )( )

β ← 1+ εβLet us then slightly increase this upper bound, i.e. set ( ) , i i −1 −1⇔+det I C ( jωI −A ) BD = 0to enforce a strict inequality in condition (5). The key idea is now ( X X XX )

I ωω ∋to compute the largest frequency interval ( ) for which the ii −1 −1⇔ det I +−jωI A BD C = 0( )( )X XX Xincreased upper bound and the scaling matrices remain valid, i.e.

such that ∀ω ∈I ω :( ) −1i ⇔ det jωI −−A BD C = 0( ( ))X XX X

−1 ⇔ det (ωIj+=H) 0 11 D(ω )Mj ( ωω)D( )−−ii

44 σωF ( ) −≤jG(ω ) F (ω ) 1 (6)i i i β

i In this context, the following algorithm is proposed to compute a gua-

ranteed robustness margin for a high-dimensional uncertain LTI plant

As is shown in proposition 2, the determination of I ω boils down [20], [10], [6]. It mainly consists of a repeated treatment on a list of ( )i

to a standard eigenvalues computation [20]. intervals.

Proposition 2 Let (A ,B ,C ,D ) be a state-space representation of Algorithm 1 (computation of a lower bound on k )

M M M M max

M(s). Build the Hamiltonian-like matrix: 1 - Initialization:

• Define an initial value for the µ upper bound, either by a µ

maxA 0 B HH **H= + X DC B lower bound computation (see Section Computation of a destabili- * ** HH H −CC −−A C D H H H zing perturbation) or by setting = 0.

max

• Define the frequency range Ω = [ω , ω ] on which the µ

min max−1

*where X= I−DD and: upper bound is to be computed. Let I= {I }= ωω , be the ( ) {[ ]}HH 1 min max

initial list of frequency intervals to be investigated.

AB HH 2 - While I ≠∅ , repeat:= (7) C D •ChooseanintervalI ∈ I and a frequency ω ∈I .i i i HH

•Compute , D(ω ), G(ω ) such that (5) holds.

i i iII00

• Set β ←+max (1 εβ) ,β and apply Proposition 2 to −1 ( )ii max 11 A − jωI BDM i M−− compute I(ω ). 44 iFF −1DC DD D − jβG MM i ββ ←• Set and update the intervals in I by eliminating the max i

β β ii frequencies contained in I(ω ).

i

3 - A lower bound on k is given by k =1/ β .LB maxmax

Define δ and δ as follows:ω ωi− i+

The proposed algorithm is not based on a frequency grid to be defined

δ= max λλ∈ℜ : det Ij+ H= 0{ ( ) }ω −i− a priori, with the risk of missing a critical frequency. On the contrary,

= −ω if jH has no positive real eigenvalue it relies on a list of frequency intervals which is updated automatically i

during the iterations. By this approach, the robustness margin is gua-δ= min λλ∈ℜ : det Ij+ H= 0{ ( ) }ω +i+

ranteed over the whole frequency range and no tricky initialization is

= ∞ if jH has no negative real eigenvalue required.

Issue 4 - May 2012 - Flight Control Laws: Recent Advances in the Evaluation of their Properties

AL04-01 3

Remark 1 Algorithm 1 can be directly applied to the general case Remark 2 Algorithm 1 can be further extended to general skew-µ

where ∆(s) is a block-diagonal LTI operator not only composed of problems, where ∆ is structured and composed of mixed

real/comc

real scalars (corresponding to parametric uncertainties), but also of plex uncertainties.

complex scalars and unstructured transfer matrices (representing

neglected dynamics). Worst case analysis

Computation of a guaranteed H performance level Computation of a destabilizing perturbation

∞

The notion of robust performance is of practical importance. Indeed, it The objective is now to compute an upper bound on k to evaluate the

max

is often desirable to quantify the performance degradations, which are conservatism of the lower bound determined in the Computation of a

induced by model uncertainties and appear before instability. The fol- guaranteed stability margin section. This is equivalent to computing a

lowing proposition is a direct consequence of the main loop theorem µ lower bound on the whole frequency range. Constructive

polynomial[17] and is used to reformulate the worst-case H performance pro- time heuristics exist [24], [21], which provide some worst-case values

∞

blem as a specific skew-µ problem. With reference to figure 2, a ficti- of ∆. They usually give fast and accurate results when ∆ contains

tious complex block ∆ is added between the output y and the input e. some complex uncertainties, but they suffer from two drawbacks. First,

c

Its size is to be maximized under the constraint that the interconnec- convergence problems are often encountered in the purely real case,

tion is stable ∀∆ ∈ B . and the resulting lower bound is then equal to 0. Second, the frequency ∆

is fixed. The problem thus has to be solved on a frequency grid, with

the risk of missing the most critical parametric configurations even if a

fine grid is used. In this context, the key idea of the method described

below and initially proposed in [6] is to directly obtain a tight µ lower

bound over the whole frequency range rather than at a fixed frequency.

In this perspective, the real µ problem considered in this paper is first M(s)

regularized by adding a small amount of complex uncertainty to

e y

each real uncertainty [18]: a perturbation ∆ is defined with the same

c

structure as ∆, except that the real scalars become complex. The

method of [24] or [21] is then applied at a given frequency ω , usually

i

with good convergence properties, to the following problem:

c Mjω εωMj( ) ( )i i Mjω =( ) ai

Figure 2 - Augmented system created by addition of a fictitious complex blockεωMj ( ) Mj ( ω ) (9)i i

∆= diag∆∆,( )acProposition 3 The following statements are equivalent:

F Ms ( ),∆• ( )is stable ∀∆ ∈ B and The resulting µ lower bound is not a lower bound for the original u ∆

0 00γγ= max F Ms ( ),∆≤( ) , real µ problem: a perturbation ∆= diag∆ ,∆ has been obtained, ( )max u ac∞

∆∈B 0∆ IM− jω ∆which renders the matrix ( ) singular, but it cannot be a ia

pp× 0I −∆Mjω • the size σ ∆ of the smallest perturbation ∆∈C such claimed that ( ) is itself singular. Nevertheless, if is ( ) c ic

that det I −Mj ( ω)diag(∆∆, ) = 0 for some ∆∈ B and some small enough, an eigenvalue of the interconnection of figure 1 (left) ( )c ∆ 0

ω ∈ℜ is larger than 1/, is usually located near the point jω of the imaginary axis.+ i

0 00µ diag II γωM j diag II γ ≤1,∀ω ∈ℜ ∆= diag δδII, ,• ( ) , Starting from this good initial guess , ( ( ) ( )) ( )∆ np,, np + 1 n Nna 1 N

pp×where ∆ = diag ∆,C . the last step is to move through the imaginary axis to obtain a ( )a 0

destabilizing perturbation for the real µ problem. More precisely, in

Proof [17] the spirit of [13], a solution is to introduce an additional perturbation

∆ = diag δδII, ,, which acts as a fictitious feedback gain. ( )1 n Nn1 N

0 Similarly to the previous section, µ is replaced with its up- The problem is then to find the smallest perturbation ∆ +∆ , which ∆a

per bound µ . For a given ω, the smallest value of such that brings the eigenvalue on the imaginary axis. As shown in [13], [6], ∆ i 0a

µ diag II γωM j diag II γ <1 can then be com- it can be recast as a simple linear programming (LP) problem:( )∆ ( ( np,,) i ( np ))a

0 puted: -νδ≤ +≤δ ν

i i

minν s.t. • either directly using an LMI characterization [8][9], (10)

δ i ℜ λ + αδ = 0( )0 ∑ ii • or iteratively using a dichotomy or a fixed-point algorithm together

with the formulation (5). where:

∂∆

The latter is usually preferred when high-dimensional systems are • α= uB++tD Cv Dw( ) ( )i ∂δianalyzed, since computational time is much lower. Algorithm 1

(especially step 2b) can thus be slightly modified to compute an upper • (A,B,C,D) is a state-space representation of M(s)

bound on the robust H performance level . • u and v are the left and right eigenvectors associated to the

UB ∞ max

−1

00eigenvalue of the closed-loop matrix A= A+B I−∆ D ∆ C ,( )0 0

−1 −1

00 00 • t=uB I−∆ D∆ and w= I−∆ D ∆ Cv( ) ( )

Issue 4 - May 2012 - Flight Control Laws: Recent Advances in the Evaluation of their Properties

AL04-01 4

corresponding to gain, phase, modulus or time-delay varia-In the above formulation, the equality constraint means that the real

tions are introduced either at the input or at the output of the part of must be equal to 0 once the additional perturbation ∆ is

0

considered open-loop system (s), depending on whether applied. An upper bound k on k is obtained as the minimum value

UB max

input or output margins are to be computed (see figure 3). of v such that (10) holds.

Note that either SISO or MIMO margins can be evaluated,

M

δδ= diag(1,,1, ,1, ,1)containing either one ( or several) Note that (10) is a linearized version of the problem to be solved. In M Mi ,

δ = diag(δδ, , ) scalar uncertainties. The expression practice, it may then be necessary to modify it if is not sufficiently ( )M M ,1 Mq ,0

of ,i is given hereafter for each margin:close to the imaginary axis. In this case indeed, the accuracy of the M

ˆˆδ =1,+ δδ∈ℜfirst order development may not be sufficient, so as to directly move • gain margin: Mi , i i

ˆˆthe eigenvalue onto the imaginary axis. A solution consists in partitio- • modulus margin: δ =1,+∈δδ C

Mi , i i

jϕining the real segment which separates from the imaginary axis, and • phase margin: δϕ=e ,∈ℜ

0 Mi , i

−τ sito iteratively perform the migration on each sub-segment. More pre- • time-delay margin: δτ=e ,∈ℜ

Mi , i +

cisely, problem (10) is solved N times, the second constraint being

N −k For gain and modulus margins, the expression of is polynomial k k k M,ireplaced at iteration k with ℜ+λ αδ =ℜ(λ ) , where ( )00∑ ii and can be written directly in linear fractional form. For phase and Nk k−1 k−11 k− 1λ= λ + αδ jϕ −τ sif k>1 and λλ= otherwise. i i00 ∑ ii e etime-delay margins, the non-rational elements and must be

transformed first in order to write the interconnection of figure 3 as a

Note that the aim here is not to compute a µ lower bound at a given linear fractional representation similar to the one of figure 1 (left). For

jϕiefrequency, but to directly obtain the highest possible lower bound this purpose, the phase variation is replaced using the bilinear

ˆˆ ˆover the whole frequency range as well as the associated frequency transformation with (1− jjδδ) / (1+ ), where δ ∈ℜ .ii i

*. Indeed, assume that a bound has been computed for the

regulajϕierized problem at a frequency , which is far from *. Using the LP This new element, similarly to , has a unitary modulus and its

i

ˆmethod above, the imaginary axis can however be crossed very close phase variation covers the whole phase range [-π, π] when δ ∈ℜ . i

−τ sieto j*, since no constraint is imposed on the imaginary part of . For the time-delay margin, the substitution of is more delicate

0

Such a behavior is generally observed in practice. It is thus sufficient because of the Laplace variable s. Nevertheless, it can be replaced by

ˆ ˆto apply the algorithm on a coarse frequency grid to obtain a tight a static rational complex function f δ , δ ∈ℜ , where the varia-( ) ii

ˆupper bound on k . tion range of δ depends on . This dependence can then be treated imax

using some results from [23]. The whole process is detailed in [11].

Remark 3 The size of ∆ in equation (9) is twice the size of the

inia

∆tial matrix ∆. But despite this, computing a µ lower bound usually

remains much faster than computing a µ upper bound.

Worst-case H performance analysis

∞

(s)K(s)

M

In the spirit of the µ lower bound algorithm in the previous section,

a two-step procedure is implemented at each point of a rough

i Input margin

frequency grid:

Output margin

• The unit ball B is investigated by iteratively:∆

σωF Mj ,∆ - computing the gradient of ( ( ) ) ,( )u i

- performing a line search to maximize this quantity (which boils M

down to computing the eigenvalues of a Hamiltonian-like matrix),

until the problem is roughly solved at . Figure 3 - Introduction of a fictitious uncertainty

i M

• Using the value of ∆ computed at step 1 as an

initialization, a quadratic programming problem, which locally maximizes The interconnection of figure 3 can now be written as in figure 1 (left),

ˆσωF Mj ,∆( ( ) ) with respect to both ∆ and ω, is repeatedly where ∆ contains both ∆ and the δ . Problem 3 is finally reformula-( )u i

solved until convergence. ted as follows: assuming that k > 1, compute the maximum value

max

ˆof the δ such that the interconnection between M(s) and ∆ is stable i

∀∆ ∈ BA lower bound of is finally obtained, as well as the associated . This is equivalent to a skew-µ problem. Both upper and Σ ∆LB max ∑

value ∆ of ∆ [19]. lower bounds on the various margins can thus be obtained using the

algorithms described in the Robustness margins computation and

Worst case analysis sections.

Tools extensions

Accuracy improvements

Unstructured margins

Conservatism is defined in this paper as the relative gap between

Structured robustness analysis considers that the uncertain- the lower and the upper bounds on a given quantity x, which can be

ties' effects on the system's behavior are perfectly identified, any of the stability margins or the performance levels considered in

which can be unrealistic. It is thus often desirable to com- the previous sections:

pute worst-case unstructured margins (gain, phase, modulus, x −x

UB LBη =time-delay) to evaluate the effective robustness of a system. (11)

x

LBIn this perspective, some additional fictitious uncertainties

M

Issue 4 - May 2012 - Flight Control Laws: Recent Advances in the Evaluation of their Properties

AL04-01 5 sometimes reaches unacceptable values, notably in the presence of cannot be guaranteed over the whole parametric domain, since

highly repeated real parametric uncertainties. A well-known technique the upper bound is larger than 1. Much better results are

obtaito ensure that it remains below a specified threshold is to use a ned with the branch and bound algorithm, since robust stability

to1

branch and bound algorithm [1], [15]. The idea is to partition the real can be guaranteed as soon as is less than 15%. Moreover,

to1

parametric domain into more and more subsets until the relative gap the additional computational cost induced by the use of branch

between the highest lower bound and the highest upper bound com- and bound is quite low thanks to the efficient strategy introduced

puted on all of the subsets becomes less than . This algorithm is before used to progressively validate the frequency domain (see

to1

known to converge for uncertain systems with only real uncertainties Section Accuracy improvements). Note that all µ upper and lower

[15], i.e. conservatism can be reduced to an arbitrarily small value. bounds are computed using the algorithms described in Sections

However, it usually exhibits an exponential growth of computational Computation of a guaranteed stability margin and Computation

complexity as a function of the number of real uncertainties. Spe- of a destabilizing perturbation respectively. Thus, the only tuning

cifying a threshold is thus used to handle the trade-off between parameter in this stability analysis is the threshold .

to1 to1

the accuracy of the bounds and the computational time.

1.5

µNevertheless, in order to alleviate the computational burden, a stra- UB

µtegy based on the progressive validation of the frequency range is LB1.25 200

proposed here. Assume for example that a subset D of the

parameN

tric domain and a frequency domain Ω are considered at step N of

N 1

the branch and bound procedure. Algorithm 1 is applied to compute 150

a frequency domain Ω Ω such that det I −Mj ( ω )∆≠ 0 ( )iv,N N

0.75∀ω ∈Ωholds ∀∆ ∈ D and . During the next step, the analysis N vN,

100performed on each subset of D then only considers the frequencies

N

in Ω which have not been validated at step N, i.e. which are not 0.5

N

contained in Ω . Consequently, after a few steps, the analysis is only No Branch and Bound

v,N

50restricted to very narrow frequency intervals corresponding to critical 0.25

frequencies. This results in a drastic reduction of the computational

CPU Timeload induced by a classical branch and bound procedure.

0 0

0 5 10 15 20 25 30 35 40 45

Conservatism (%)

to1 Application to flight control laws validation

Figure 4 - µ bounds and CPU time versus conservatism

The algorithms described in the previous sections are now evaluated

on a realistic application. All calculations are performed on a 3GHz Worst-case H performance analysis

∞

PC with 3GB RAM.

Worst-case H performance is evaluated for the transfer function

∞

Description of the model between the vertical wind velocity and the vertical load factor. In order

to identify the secondary peaks of the frequency response, the

anaA high fidelity model composed of 22 states is considered here. It lysis is performed on three contiguous frequency intervals. A skew-µ

describes both the rigid and the flexible closed-loop longitudinal dy- problem is thus solved on each one of these intervals. Figure 5 shows

namics of a civilian passenger aircraft. It is parameterized by 4 real the bounds on obtained with = 20% as well as the frequency

max to1

parameters characterizing the aircraft's mass configuration: center responses of the uncertain system computed on a fine parametric

and outer tanks filling levels and , embarked payload and grid. Results are very accurate.

CT OT PL

position of the center of gravity . The model is written in linear

CG

fractional form as shown in figure 1 using the LFR Toolbox for Matlab Bode Diagram

60

[12]. As the effects of the parameters on the system behavior are

modeled very accurately, the size of ∆ is quite large: 50

∆ = diag(δ I ,δ II,δδ , I ) 40CT 48 OT 28 PL 15 CG 4

30∆ is normalized, which means that the whole parametric domain is

covered when ∆ takes all possible values in B .∆

20

Robust stability analysis 10

Robust stability is first analyzed in order to check whether sta- 0

Upper Boundbility can be guaranteed over the whole parametric domain. For

Lower Bound

-10this purpose, several µ upper and lower bounds are computed,

0 5 10 15 20 25 30 35 40 and the results are illustrated in figure 4. The bounds are first

computed without branch and bound. A relative gap of about 40% Frequency (rad/sec)

is obtained and the computational time is very reasonable

considering the large size of the model. Nevertheless, robust stability Figure 5 - Bounds on and frequency responses on a fine parametric grid

max

Issue 4 - May 2012 - Flight Control Laws: Recent Advances in the Evaluation of their Properties

AL04-01 6

µ

Magnitude (dB)

CPU Time (s)Worst-case unstructured margins Conclusion and prospects

Worst-case SISO unstructured margins are finally computed. With re- Several µ-analysis based tools developed by the systems control

ference to figure 3, the open-loop system (s) is composed of actua- group of Onera are reviewed in this paper. They are used to

comtors, open-loop aircraft and sensors models in a feedback loop with pute both upper and lower bounds on the robust stability margin,

a dynamic controller K(s). The input of (s) is the elevator deflec- the worst-case H performance level, as well as the worst-case

∞

tion, while the outputs are the pitch rate and the vertical load factor. gain, phase, modulus and time-delay margins. Unlike most existing

Figure 6 shows the bounds on the gain, modulus and phase margins methods, these bounds are guaranteed over the whole frequency

obtained with = 25%, and the Nyquist responses of the uncertain range, and not only on a finite frequency grid. Moreover, an efficient

to1

system on a fine parametric grid. Once again, results are quite satis- branch and bound scheme can be used to obtain bounds with the

factory and conservatism is efficiently mastered. desired accuracy, while still guaranteeing a reasonable computa-

tional complexity. These algorithms which will form the basis of a

next release of the Skew Mu Toolbox for Matlab [7] should enable to Nyquist Diagram

considerably improve the flight control systems validation process in 1 2 dB 0 dB the near future

0.8 4 dB

0.6 6 dB

0.4 -10 dB10 dB

0.2

20 dB

0

-0.2

M optimistic-0.4 g

M guaranteedg

-0.6 M optimisticm

M guaranteedm-0.8

M optimistic

-1 M guaranteed

-1 -0.8 -0.6 -0.4 -0.2 0 0.2

Real Axis

Figure 6 - Upper (optimistic) and lower (guaranteed) bounds on the worst

case gain (M ), modulus (M ) and phase (M ) margins and Nyquist

resg m

ponses on a fine parametric grid

Issue 4 - May 2012 - Flight Control Laws: Recent Advances in the Evaluation of their Properties

AL04-01 7

Imaginary Axis