- Email: [email protected]
Robustness Assessment of a Load Factor Flight Control Law Using Dilation Integrals*
20th 20th IFAC IFAC Symposium Symposium on on Automatic Automatic Control Control in in Aerospace Aerospace 20th IFAC Symposium on Automatic Control in August 21-25, 2016. Sherbrooke, Quebec, Canada 20th IFAC Symposium on Automatic Control in Aerospace Aerospace 20th IFAC Symposium on Automatic Control in Aerospace August 21-25, 2016. Sherbrooke, Quebec, Canada August 21-25, 2016. Sherbrooke, Quebec, Canada Available online at www.sciencedirect.com August Canada August 21-25, 21-25, 2016. 2016. Sherbrooke, Sherbrooke, Quebec, Quebec, Canada
ScienceDirect IFAC-PapersOnLine 49-17 (2016) 094–099
Robustness Assessment Assessment of of a a Load Load Factor Factor Robustness Robustness Assessment of a Load Factor Robustness Assessment of a Load Factor Flight Control Law Using Dilation Flight Control Law Using Dilation Flight Control Law Using Dilation Flight Control Law Using Dilation Integrals Integrals Integrals Integrals K. Bessadi Bessadi ∗∗∗ D. D. Saussi´ Saussi´ e ∗∗∗ L. L. Saydy Saydy ∗∗∗ K. e ∗ ∗ ∗ K. Bessadi Saussi´ e ∗ D. ∗ L. K. Bessadi D. Saussi´ e L. Saydy K. Bessadi D. Saussi´ e L. Saydy Saydy ∗ ∗ ∗ Polytechnique Montr´ eeal, QC, H3T 1J4, Canada QC, H3T 1J4, Canada ∗ Polytechnique Montr´ ∗ Polytechnique Montr´ eeeal, al, QC, H3T 1J4, Canada ∗ Polytechnique Montr´ al, QC, H3T 1J4, Canada (e-mail: {karim.bessadi,david.saussie,lahcen.saydy}@polymtl.ca) Polytechnique Montr´ al, QC, H3T 1J4, Canada (e-mail: {karim.bessadi,david.saussie,lahcen.saydy}@polymtl.ca) (e-mail: {karim.bessadi,david.saussie,lahcen.saydy}@polymtl.ca) (e-mail: {karim.bessadi,david.saussie,lahcen.saydy}@polymtl.ca) (e-mail: {karim.bessadi,david.saussie,lahcen.saydy}@polymtl.ca)
Abstract: This article presents a performance assessment approach of aa load factor flight Abstract: This article presents aa robust robust performance assessment approach of factor flight Abstract: This article presents robust performance assessment approach of aaa load load factor flight Abstract: This article presents a robust performance assessment approach of load factor flight control law with the goal of ascertaining whether the closed-loop system’s poles all lie within Abstract: This article presents a robust performance assessment approach of load factor flightaa control law with the goal of ascertaining whether the closed-loop system’s poles all lie within control law with the goal of ascertaining whether the closed-loop system’s poles all lie within control law with the goal of ascertaining whether the closed-loop system’s poles all lie within aaa given subregion of the open left half-plane despite center of gravity and mass variations. The controlsubregion law with of thethe goal of ascertaining whether the closed-loop system’s polesvariations. all lie within given open left half-plane despite center of gravity and mass The given subregion of the open left half-plane despite center of gravity and mass variations. The given subregion of the open left half-plane despite center of gravity and mass variations. The method combines a guardian map approach with dilation integrals to tackle the underlying given subregion of the open left half-plane despite center of gravity and mass variations. The method combines a guardian map approach with dilation integrals to tackle the underlying method combines guardian map approach with dilation to tackle underlying method polynomial combines a a positivity guardian test mapby approach withlatter dilation integrals to tackle the the underlying difficult relaxing the to aa integrals so-called ‘practical’ positivity one. method combines a guardian map approach with dilation integrals to tackle the underlying difficult polynomial positivity test by relaxing the latter to so-called ‘practical’ positivity one. difficult polynomial positivity test by relaxing the latter to a so-called ‘practical’ positivity one. difficult polynomial positivity test by relaxing the latter to a so-called ‘practical’ positivity one. The study is carried carried positivity out for for a a longitudinal longitudinal uncertain model of a F-16 F-16 aircraft aircraft with positivity varying center center difficult polynomial test by relaxing the latter to of a so-called ‘practical’ one. The study is out uncertain model a with varying The study is carried out for a longitudinal uncertain model of a F-16 aircraft with varying center The study is carried out for a longitudinal uncertain model of a F-16 aircraft with varying center of gravity position and mass with the purpose of testing such an approach. The study position is carriedand outmass for a with longitudinal uncertain modelsuch of a an F-16 aircraft with varying center of gravity the purpose of testing approach. of of gravity gravity position position and and mass mass with with the the purpose purpose of of testing testing such such an an approach. approach. of gravity position and mass with the purpose of testing such an approach. © 2016, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: dilation dilation integrals, integrals, generalized generalized stability, stability, guardian guardian maps, maps, robustness robustness analysis, analysis, flight flight Keywords: Keywords: dilation integrals, generalized stability, guardian maps, robustness analysis, flight Keywords: dilation generalized control, polynomial positivity. Keywords: dilation integrals, integrals, generalized stability, stability, guardian guardian maps, maps, robustness robustness analysis, analysis, flight flight control, polynomial positivity. control, polynomial positivity. control, control, polynomial polynomial positivity. positivity. 1. INTRODUCTION INTRODUCTION under consideration consideration is is deemed deemed positive positive from from aa practical practical 1. under 1. INTRODUCTION under consideration 1. under consideration is deemed positive from practical viewpoint. 1. INTRODUCTION INTRODUCTION under consideration is is deemed deemed positive positive from from aaa practical practical viewpoint. viewpoint. viewpoint. viewpoint. The present work deals broadly with the generalized In Section Section 2, 2, we we briefly briefly introduce introduce guardian guardian maps maps to to study study The present work deals broadly with the generalized In The present work deals broadly with the generalized The work broadly with generalized stability of parametrized of uncertain systems. In Section 2, briefly guardian maps The present present work deals dealsfamilies broadly with the theLTI generalized In Section 2, we briefly introduce guardian maps to study the stability of families of parametrized mastability of parametrized families of uncertain LTI systems. In Section 2, we we briefly introduce introduce guardianstate-space maps to to study study the stability of families of parametrized state-space mastability of parametrized families of uncertain LTI systems. stability of families of LTI ensures More specifically, given a controller C the stability of families of parametrized state-space mastability of parametrized parametrized families of uncertain uncertain LTI systems. systems. 0 which the stability of families of parametrized state-space matrices with respect to a given region Ω ⊂ C. Examples which ensures More specifically, given a controller C the stability of families of parametrized state-space ma0 trices with respect to a given region Ω ⊂ C. Examples 1 which ensures More specifically, given a controller C 0 which ensures More specifically, given controller C 1a the generalized stability of the system for a nominal trices with respect to a given region Ω ⊂ C. Examples 0 which ensures More specifically, given a controller C trices with respect to a given region Ω ⊂ C. Examples 0 for a nominal of guardian maps are given for simple regions of interest, the generalized stability of the system 1 trices with respect to a given region Ω ⊂ C. Examples 1 of guardian maps are given for simple regions of interest, k the generalized stability system for a 1 of the generalized stability of the system for a nominal k value of an uncertain parameter qq 000 ∈ we of maps for regions of the generalized stability of the thevector system forR a,, nominal nominal of guardian maps are given for simple regions of interest, as guardian well as as aa necessary necessary and sufficient sufficient condition provided in value of an uncertain parameter vector ∈ R we wish wish k of guardian maps are are given given for simple simple regions provided of interest, interest, 0 k as well and condition in value of an uncertain parameter vector q ∈ R wish 0 ensures k ,, we value of an uncertain parameter vector q ∈ R we wish to ascertain whether the same controller robust as well as a necessary and sufficient condition provided in value of an uncertain parameter vector q ensures ∈ R , werobust wish the as well as a necessary and sufficient condition provided in form of a positivity test involving a single inequality. to ascertain whether the same controller as well asofa necessary and sufficient condition provided in the form a positivity test involving a single inequality. to ascertain whether the same controller ensures robust to ascertain whether the same controller ensures robust stability with respect to Ω when the uncertain parameter the form of a positivity test involving a single inequality. to ascertain whether the same controller ensures robust the form of a positivity test involving a single inequality. In Section 3, the dilation integral method is presented. stability with respect to Ω when the uncertain parameter the form of a positivity test involving a single inequality. In Section 3, the dilation integral method is presented. k stability with respect to Ω when the uncertain parameter stability with within respecta to Ω when when theQ uncertain parameter aa broad vector varies given subset ⊂ In Section 3, dilation integral method presented. stability with respect to Ω the uncertain parameter In Section 3, the dilation integral method is presented. An upper bound Vol(Q is obtained dilation For vector varies within a given subset Q ⊂R Rkkkk ... For In Section 3, the the on dilation integral method is isvia presented. bad ) An upper bound on Vol(Q is obtained via dilation For aaa broad broad vector within a Q bad ) For broad vector varies within a given subset Q ⊂ R class of ofvaries regions and uncertain uncertain LTI systems, systems, the guardian An upper bound on Vol(Q ) is obtained via dilation Forguardian broad integral vector varies within a given given subset subset Q⊂ ⊂R R ..the bad An upper bound on Vol(Q ) is obtained via dilation computations and a conditioning parameter inclass regions and LTI bad An upper bound on Vol(Q ) is obtained via dilation bad integral computations and a conditioning parameter inclass of regions and uncertain LTI systems, the guardian class of regions and uncertain LTI systems, the guardian map approach (Saydy et al., 1990) yields a necessary integral computations and a conditioning parameter inclass of regions (Saydy and uncertain LTI systems, the guardian integral integral computations and a conditioning parameter introduced. In Section4, we combine the two approaches to map approach et al., 1990) yields a necessary computations and a conditioning parameter introduced. In Section4, we combine the two approaches to map approach (Saydy et al., 1990) yields a necessary map approach (Saydy et al., 1990) yields a necessary and sufficiant condition for robustness in the form of troduced. In Section4, we combine the two approaches to map approach (Saydy et al., 1990) yields a necessary troduced. In Section4, we combine the two approaches to analyse the robust stability of a longitudinal flight control and sufficiant condition for robustness in the form of troduced. In Section4, we combine the two approaches to analyse the robust stability of a longitudinal flight control and sufficiant condition for robustness in the form of and sufficiant condition for robustness in the form of a positivity test of a single multivariable polynomial, analyse the robust stability of a longitudinal flight control and sufficiant condition for robustness in the form of analyse the robust stability of a longitudinal flight control system. More specifically, we analyse the robust stability a positivity test of a single multivariable polynomial, analyse the robust stability of a longitudinal flight control system. More specifically, we analyse the robust stability a positivity test of a single multivariable polynomial, a positivity test of single namely, decide whether or not: multivariable More specifically, we the robust stability a positivity test of a a or single multivariable polynomial, polynomial, system. system. More specifically, we analyse the robust stability of aa nominal factor flight control law to see whether namely, decide whether not: system. Moreload specifically, we analyse analyse the robust stability of nominal load factor flight control law to see whether namely, decide whether or not: namely, decide whether or not: of a nominal load factor flight control law to see whether namely, decide whether or not: of a nominal load factor flight control law to see whether the closed-loop poles preserve a certain negativity margin of a nominal load factor flight control law to see whether ν(q) > 0 ∀q ∈ Q the closed-loop poles preserve a certain negativity margin ν(q) > 0 ∀q ∈ Q the closed-loop poles preserve a certain negativity margin ν(q) > 0 ∀q ∈ Q the closed-loop poles preserve a certain negativity margin and a certain damping, despite the variations of center of ν(q) > >0 0 ∀q ∀q ∈ ∈Q Q the closed-loop poles preserve a certain negativity margin ν(q) and a certain damping, despite the variations of center of where ν(q) is a polynomial constructed from the guardian and a certain damping, despite the variations of center of where ν(q) is a polynomial constructed from the guardian and a certain damping, despite the variations of center of gravity and mass. Various assessments using the dilation and a certain damping, despite the variations of center of where ν(q) is a polynomial constructed from the guardian gravity and mass. Various assessments using the dilation where ν(q) is a polynomial constructed from the guardian map associated with the problem under study (cf. Section where ν(q) is a polynomial constructed from the guardian gravity and mass. Various assessments using the dilation map associated with the problem under study (cf. Section gravity and mass. Various assessments using the dilation integrals are presented. gravity and mass. Various assessments using the dilation map associated with the problem under study (cf. Section integrals are presented. map associated with with the the problem problem under under study study (cf. (cf. Section Section integrals are presented. 2). map associated 2). integrals are are presented. presented. integrals 2). 2). 2). Owing to the NP-hardness of tackling such a proposition Owing to the NP-hardness of tackling such a proposition Owing to NP-hardness of such proposition Owing to the the NP-hardness of tackling tackling such a a the proposition 2. GUARDIAN in general, Barmish et al. (2009) introduced concept Owing to the NP-hardness of tackling such a proposition 2. GUARDIAN MAPS MAPS in general, Barmish et al. (2009) introduced the concept 2. in general, Barmish et al. (2009) introduced the 2. GUARDIAN MAPS in general, Barmish et al. (2009) introduced the concept of general, ‘pratical’Barmish positivity to alleviate alleviate the complexity complexity issue 2. GUARDIAN GUARDIAN MAPS MAPS in et al. (2009) introduced the concept concept of ‘pratical’ positivity to the issue of ‘pratical’ positivity to alleviate the complexity issue of ‘pratical’ positivity to alleviate the complexity issue associated with precisely deciding true positivity of a of ‘pratical’with positivity to alleviate the complexity issue We are here interested in stability sets of the form: associated precisely deciding true positivity of a are here interested in stability sets of the form: associated with precisely deciding true positivity of aaa We associated with precisely In deciding trueet positivity ofthe We are here interested in stability sets of the form: multivariable polynomial. (Barmish al., 2009), associated with precisely deciding true positivity of We are here interested in stability sets of the form: n×n We are here interested in stability sets of the multivariable polynomial. In (Barmish et al., 2009), the n×n : σ(A) ⊂ Ω} form: S(Ω) = {A ∈ R (1) multivariable polynomial. In (Barmish et al., 2009), the multivariable polynomial. In (Barmish et al., 2009), the S(Ω) = {A ∈ R : σ(A) ⊂ Ω} (1) n×n authors introduced the method of dilation integrals, multivariable polynomial. In (Barmish et al., 2009), the n×n authors introduced the method of dilation integrals, the S(Ω) = {A ∈ R ⊂ Ω} (1) n×n :: σ(A) S(Ω) = {A ∈ R σ(A) ⊂ Ω} (1) authors introduced the method of dilation integrals, the S(Ω) = {A ∈ R : σ(A) ⊂ Ω} (1) authors introduced the method of dilation integrals, the purpose introduced of which which is is to to decide whether a given givenintegrals, polynomial, Ω is an open subset of the complex plane, and σ(A) authors thedecide method of dilation the where purpose of whether a polynomial, where Ω is an open subset of the complex plane, and σ(A) purpose of which is to decide whether a given polynomial, where Ω is an open subset of the complex plane, and σ(A) purpose of which is to decide whether a given polynomial, where Ω is an open subset of the complex plane, and σ(A) such as ν(q), is practically positive for (almost) all q ∈ denotes the spectrum of A. Such sets S(Ω) will be referred purpose of which is to decidepositive whetherfor a given polynomial, where Ωthe is an open subset the sets complex plane, σ(A) such as ν(q), is practically (almost) all qq ∈ denotes spectrum of A.ofSuch S(Ω) will beand referred such as is positive for all ∈ the of Such sets S(Ω) will such as ν(q), ν(q), is practically practically positive for (almost) (almost) all ∈ denotes denotes the spectrum spectrum of A. A.sets, Suchand setsthus S(Ω)represent will be be referred referred Q. More specifically, the method yields a means to to as generalized stability the set such as ν(q), is practically positive for (almost) all qqtest ∈ denotes the spectrum of A. Such sets S(Ω) will be referred Q. More specifically, the method yields a means to test to as generalized stability sets, and thus represent the set Q. More specifically, the method yields a means to test to as generalized stability sets, and thus represent the set Q. More specifically, the method yields a means to test to as generalized stability sets, and thus represent the set whether the set of parameters for which the polynomial of all matrices which are stable relative to Ω, i.e., which Q. More specifically, the method yields a means to test to as generalized stability sets, and thus represent the set whether the set of parameters for which the polynomial of all matrices which are stable relative to Ω, i.e., which whether the set of parameters for which the polynomial of all matrices which are stable relative to Ω, i.e., which whether the set of parameters for which the polynomial of all matrices which are stable relative to Ω, i.e., which is not positive (so-called Q ) has a volume below a have all their eigenvalues in Ω. Guardian maps are defined whether the set (so-called of parameters which the polynomial bad for of all matrices which are stable relative to Ω, i.e., which is not positive Q ) has a volume below a have all their eigenvalues in Ω. Guardian maps are defined is (so-called Q ))) has a below have all their eigenvalues in Ω. Guardian maps are defined bad is not not positive positive (so-called Qbad has case a volume volume below aaa next. have all certain tolerable threshold , in which the polynomial bad is not positive (so-called Q has a volume below have all their their eigenvalues eigenvalues in in Ω. Ω. Guardian Guardian maps maps are are defined defined bad certain tolerable threshold , in which case the polynomial next. certain tolerable threshold , in which case the polynomial next. certain tolerable tolerable threshold threshold , , in in which which case case the the polynomial polynomial next. next. n×n certain n×n into C, Definition 1. (Saydy et al., 1990) Let ν map R Definition 1. (Saydy et al., 1990) Let νν¯ map R C, n×n n×n into Definition 1. (Saydy et al., 1990) Let map R into C, Definition 1. et al., 1990) Let νν¯ map we say that guards S if for all A ∈ S(Ω): Definition 1.νν (Saydy (Saydy et al., 1990) Let map R Rn×n into into C, C, work was supported by NSERC under Grant RGPIN-2014 This we say that guards S if for all A ∈ S(Ω): ¯ This work was supported by NSERC under Grant RGPIN-2014 ¯ we say that ν guards S if for all A ∈ S(Ω): ¯ we say that ν guards S if for all A ∈ S(Ω): This work was supported by NSERC under Grant RGPIN-2014 we say that ν guards S if for all A ∈ S(Ω): 03942 and RGPIN-2012-122106. This work was supported by NSERC under Grant RGPIN-2014This work was supported by NSERC under Grant RGPIN-201403942 and RGPIN-2012-122106. ν(A) = = 00 ⇔ ⇔A A∈ ∈ ∂S(Ω) ∂S(Ω) (2) 1 ν(A) (2) 03942 and RGPIN-2012-122106. In of 03942 and RGPIN-2012-122106. 1 ν(A) = ⇔ A ∈ ∂S(Ω) (2) 03942 andsense RGPIN-2012-122106. In the the sense of placing placing closed-loop closed-loop poles poles within within a a given given region region ν(A) = 000 ⇔ A ∈ ∂S(Ω) (2) 1 ν(A) = ⇔ A ∈ ∂S(Ω) (2) 1 ¯ In the sense of placing closed-loop poles within a given region 1 Here S S¯ denotes denotes closure closure of of the the set set S S and and ∂S ∂S its its boundary. boundary. Ω In ⊂ C. In C. the sense sense of of placing placing closed-loop closed-loop poles poles within within a a given given region region the Here Ω ⊂ ¯ ¯ Here S denotes closure of the set S and ∂S its boundary. Ω ⊂ C. ¯ Here S denotes closure of the set S and ∂S its boundary. ⊂ C. Ω Here S denotes closure of the set S and ∂S its boundary. Ω ⊂ C. Copyright 2016 IFAC 94 2405-8963 Copyright©© ©2016, 2016IFAC IFAC (International Federation of Automatic Control) 94 Hosting by Elsevier Ltd. All rights reserved. Copyright © 2016 IFAC 94 Copyright 2016 IFAC 94 Copyright ©under 2016responsibility IFAC 94 Control. Peer review © of International Federation of Automatic 10.1016/j.ifacol.2016.09.017
K. Bessadi et al. / IFAC-PapersOnLine 49-17 (2016) 094–099
A(q1 , q2 ) =
• Hurwitz stability :
Let Ω = C◦− (Fig. 1(a)). Then S(Ω) is guarded by:
(3) νH : A → det(A) det(A I) where denotes the bialternate product (cf. Appendix A).
A corresponding guardian maps is given by νΩ = νσ νξ where (Saydy et al., 1990): (8) νσ (A) = det(A + σI) det((A + σI) I) 2 2 νξ (A) = det(A) det(A I + (1 − 2ξ )A A) (9) We obtain : 5 νΩ (A) = − (q1 + 2)(9q1 − 6q2 + 16) 4 × (q1 + 6q2 )(25q12 + 62q1 + 12q2 + 36) (10)
• Negativity margin:
Let Ω be the σ-shifted open left half-plane (Fig. 1.b). Then S(Ω) is guarded by : νσ : A → det(A + σI) det((A + σI) I) (4) • Damping ratio:
Let Ω be the conic sector with inner angle 2θ (Fig. 1.c). Then S(Ω) is guarded by: νξ : A → det(A) det(A2 I + (1 − 2ξ 2 )A A) with ξ = cos θ.
For q 0 = [−2.75, −2.5] , the matrix A(q 0 ) has eigenvalues −3.875 ± j1.654 and is therefore nominally Ω−stable. Furthermore, νΩ (A(q 0 )) < 0. Therefore, according to Theorem 1, the whole family is Ω−stable within Q1 if, and only if, νΩ (A(q)) changes no sign in Q1 , that is νΩ (q) = (q1 + 2)(9q1 − 6q2 + 16) × (q1 + 6q2 )(25q12 + 62q1 + 12q2 + 36) > 0 (11) for all q ∈ Q1 .
(5)
• Schur stability :
Let Ω be the circle of radius ω (Fig. 1.d). Then S(Ω) is guarded by: νω : A → det(A + ωI) det(A − ωI) det(A A − ω 2 I I) (6)
In Figs 2.a and 2.b, setting νΩ (r) = 0 divides the parameter space into many components that are either stable or unstable. One needs to test only one vector in each component to see what prevails (Saussi´e et al., 2006). The green component, being nominally Ω-stable, corresponds in fact to the largest set containing the nominal value for which Ω-stability is preserved. In particular, we have robust Ω-stability within Q1 but not with Q2 = [−3, −2]× [−4, −2] as the guardian map vanishes within the latter parameter set.
A systematic way of constructing guardian maps for various Ω regions can be found in (Saydy et al., 1990). Im
Im
Im
ω 0 Re
(a)
−
σ
Re
(b)
θ
(c)
Re
0 1 (7) 6q2 + q1 5q1 + 6 We wish to study the stability of this family relative to the σ−shifted left half-plane with σ = 2, and to the conic √ sector with θ = π4 ) (ξ = 22 ), i.e., relative to Ω being the intersection of those in Figs 1.b and 1.c.
2.1 Some examples of guardian maps
Im
95
Re
(d)
2
2
1
1
0
q2
Fig. 1. Useful stability sets Let {A(q) : q ∈ Q} be a continuous family of real n × n matrices depending on an uncertain parameter q known to belong to pathwise connected subset Q ⊂ Rk . Let us also assume that this family is nominally stable relative to a given Ω, that is A(q 0 ) ∈ S(Ω) for some q 0 ∈ Q. While the stability of this family amounts in general to testing a number of inequalities for each A(q) (depending on n and Ω), the result below states that only one amongst these needs to be tested for robust stability. Theorem 1. (Saydy et al., 1990) Let S(Ω) be guarded by the map νΩ . The family {A(q) : q ∈ Q} is stable relative to Ω if and only if:
0
q2
IFAC ACA 2016 August 21-25, 2016. Quebec, Canada
-1
-1
-2
-2
-3
-3
-4
-4 -3
-2.5
-2
-1.5
-1
-0.5
0
-3
-2.5
-2
q1
(a) Stable set Q1
-1.5
-1
-0.5
0
q1
(b) Unstable set Q2
Fig. 2. Stable and unstable components 3. DILATION INTEGRAL METHOD Here, we give a brief overview of the dilation integral method of (Barmish et al., 2009) Let f (q) be a multivariable polynomial and let Qbad = {q ∈ Q : f (q) ≤ 0}, the subset of Q for which f (q) is not positive. In order for f (q) > 0, it is thus equivalent Qbad = ∅. However, as this is computationnally very hard in general, the dilation integral method proposes to judge the positivity from a practical point of view by attempting to show that the percentage of positivity violation as bad ) typified by the volume ratio Vol(Q Vol(Q) is below a certain value deemed tolerable.
(1) it is nominally stable, i.e. A(q0 ) ∈ S(Ω) for some q0 ∈ Q ; and (2) ν(A) := νΩ (A(q))νΩ (A(q0 )) > 0, for all q ∈ Q. 2.2 Example Let A(q1 , q2 ) be a family of matrices depending on two parameters −3 ≤ q1 ≤ −2.5 and −3 ≤ q2 ≤ −2 (Q1 = [−3, −2.5] × [−3, −2]) : 95
IFAC ACA 2016 96 August 21-25, 2016. Quebec, Canada
K. Bessadi et al. / IFAC-PapersOnLine 49-17 (2016) 094–099
with increasing k, the symbolic computation can be intractable in a decent time while quadrature formulas can efficiently alleviate the burden.
The dilation integral associated with the pair (f (q), Q) is given by: 1 k (1 − αf (q)) dq (12) k (α) = Vol(Q) Q with α ≥ 0 and k being a positive even integer.
3.1 The conditioner In order to stop the iterations and declare the pair (f, Q) to be non positive from a practical point of view, a new parameter is introduced (Shcherbakov and Barmish, 2003). The conditioner θ is defined as the maximum percentage variation of f (q) about the midpoint of its range. That is, with
For every α positive, the inequality below gives an upper bound for the violation ratio: Vol(Qbad ) ≤ k (α) (13) Vol(Q) It turns out that this upper bound achieves its minimum so that: Vol(Qbad ) . ≤ k = min k (α) (14) α≥0 Vol(Q) and we have the following properties: Theorem 2. (Barmish et al., 2009) For the given pair (f, Q) with associated dilation integrals 1 . k k (α) = (1 − αf (q)) dq (15) Vol(Q) Q defined for positive even integers k, their minima, . k = min k (α)
fmin = minf (q), and fmax = maxf (q) q∈Q
α≥0
To estimate the conditioner’s value, one relies on the fact that . 1/k θ k = k ≤ θ (17)
are attained. Moreover, the following conditions hold: (1) If (f, Q) is a positive pair, then k → 0. (2) If k → 0, then (f, Q) is a non-negative pair with Vol(Qbad ) = 0. Example 1. Let us consider the polynomial (Barmish et al., 2009)
So to stop the computations when one is unable to reach values of k below the acceptable threshold , we use the sequence θk as a revealing indicator of an imminent risk of losing positivity when its value is close to 1. Theorem 3. (Barmish et al., 2009) If (f, Q) is a positive pair, it follows that lim θk = θ (18)
f (x) = 1 + x21 x22 (x21 + x22 − 3) which as it turns out is positive for any Q = [−r, r]2 with 0 < r < 1.
k→∞
For r = 0.75 and 6 (α), we obtain the following: 0.75 0.75 1 6 6 (α) = (1 − αf (x)) dx1 dx2 . (2 · 0.75)2 −0.75 −0.75 ≈ 0.7054α6 − 4.40175α5 + 11.5030α4 − 16.13601α3 + 12.8394α2 − 5.5096α + 1.0.
If (f, Q) is not a positive pair, then lim θk = 1 k→∞
(19)
Moreover, in both cases, the sequence of estimates θk is non-decreasing. Example 2. Let us consider the multivariate polynomial f (q) = q12 (q12 + 1) q22 ((q1 + 0.2)2 + (q2 − 0.3)2 − 0.22 ) with (q1 , q2 ) ∈ [−0.5, 0.5] × [−0.5, 0.5].
The minimum of this function is 6 = 0.00011355 which in turn ensures that positivity is violated on a set whose volume is at worst equal to 0.00011355. If greater accuracy is required, similar integrations may be carried out to try to achieve smaller values (Table 2).
For this example, it is easy to compute the violation ratio: Vol(Qbad ) = 0.22 π ≈ 0.1256 Vol(Q) The results obtained with the dilation integral are given in the table below: Table 2. Utility of θk
Table 1. Convergence de k k 2 4 6 8 10 12 14 16
q∈Q
(fmax − fmin ) (16) θ= | fmax + fmin | In case of nominally positive polynomials f (q), (i.e., f (q0 ) > 0 for some q0 ∈ Q), it is clear that if 0 ≤ θ < 1 then positivity is ensured over the set Q. But if θ → 1, then it would mean that there are parameter configurations where the risk of violating the positivity becomes imminent (Barmish et al., 2009).
k 0.0149 0.0011010 0.0001135 1.36 × 10−5 1.77 × 10−6 2.44 × 10−7 3.45 × 10−8 5.30 × 10−9
k 2 4 6 8 10 12 14 16 .. . 26
Remark 1. To compute symbolic expressions k (α) of the dilatation integrals and the corresponding minimum k quadrature formulas such as Smolyak’s may be used (Smolyak, 1963; Dabbene and Shcherbakov, 2007). Indeed, 96
k 0.85 0.80 0.77 0.75 0.73 0.71 0.70 0.69 .. . 0.65
θk 0.92 0.94 0.95 0.96 0.97 0.97 0.97 0.97 .. . 0.98
IFAC ACA 2016 August 21-25, 2016. Quebec, Canada
K. Bessadi et al. / IFAC-PapersOnLine 49-17 (2016) 094–099
Note that the sequence k is decreasing. However since the conditioner is getting closer to 1, we can stop the calculations and declare that positivity is not practically ensured.
Finally, the closed-loop state-space model can be written as: x˙ = Acl x + Bcl nzc nz = Ccl x and the matrices are: where x = α q δe xi δec δ˙ec A B 0 0
4. FLIGHT CONTROL ROBUSTNESS ANALYSIS
2×1
This section illustrates the robustness analysis of a load factor flight control law extracted from (Lhachemi et al., 2015). For a fixed flight condition, the uncertainties are the mass variation ∆m (in kg) and the center of gravity location ∆x (in m) around nominal values.
Acl =
(23) Bcl = [0 0 0 1 0 0] Ccl = [Cnz Dnz 0 0 0] (24) where Cnz and Dnz denote the row vectors corresponding to nz output. The dependencies in ∆x and ∆m have been dropped to alleviate notation. Nevertheless the jacobian matrix Acl is still a quadratic polynomial matrix in ∆x and ∆m. 4.2 Robustness analysis Our first objective is to study the stability of Acl (∆x, ∆m) relative to the σ−shifted left half-plane on the set Q1 = [−0.15, 0.15] × [−1000, 1000]. Vanishing the corresponding guardian map for σ = {0, 1.5, 1.75, 2} leads to figure 4. σ=0
where Aij are constant coefficient matrices. Numerical values can be found in Appendix B. The elevator dynamics is modeled as a first order transfer function: δe 1 = Fact (s) = δ ec τs + 1 where δec denotes the commanded elevator deflection with time constant τ = 0.05.
Kp
−
+ −
Fro
Fact
δe
1000
1000
500
500
0
0
-500
-500
-1000
-1000
-0.1
0
0.1
-1500 -0.2
0.2
-0.1
∆x (m) σ = 1.75
nz A/C q
Kq
1500
1000
1000
500
500
0
-500
-1000
-1000 0
∆x (m)
0.2
0.1
0.2
0
-500
-0.1
0.1
σ=2
1500
-1500 -0.2
0
∆x (m)
∆m (kg)
+
1500
-1500 -0.2
∆m (kg)
Ki /s
δ ec
σ = 1.5
1500
∆m (kg)
∆m (kg)
0≤i,j≤2 i+j≤2
−
(22)
We consider the longitudinal model of a F-16 aircraft linearized at Mach 0.9 and altitude 5000 m. The classical short period approximation is given by: α˙ α q˙ A(∆x, ∆m) B(∆x, ∆m) q (20) n z C(∆x, ∆m) D(∆x, ∆m) δe q where the state variables are the angle of attack α and the pitch rate q, the input variable is the elevator deflection δe and the output measurements are q and the load factor nz . The state-space matrices (A, B, C, D) are approximated by quadratic polynomial matrices in ∆x and ∆m. For instance, A(∆x, ∆m) = Aij ∆xi ∆mj (21)
+
2×2
01×2 −1/τ 0 1/τ 0 −Cnz −Dnz 0 01×2 0 0 0 1 01×2 2 2 D 2 2 − Kp Kq ωro C −Kp ωro nz Ki ωro −ωro −2ξωro
4.1 Model
nz c
97
0.1
0.2
-1500 -0.2
-0.1
0
∆x (m)
Fig. 4. Stable et unstable sets for different values of σ
Fig. 3. Load factor control law
The nominal case (∆x, ∆m) = (0, 0) is stable for the four values of σ. For σ = {0, 1.5}, the guardian map does not vanish on Q1 , then Acl (∆x, ∆m) is stable for all (∆x, ∆m) ∈ Q1 . Whereas, for σ = {1.75, 2}, the guardian map does vanish on Q1 and one can then deduce that the family is unstable.
The control law depicted in Fig. 3 is comprised of: • a PI controller with integrator Ki /s on the tracking error nzc − nz and proportionnal feedback Kp on the load factor nz • a proportionnal feedback Kq on the pitch rate q • a roll-off filter Fro (s) to limit the controller bandwidth: 2 ωro Fro (s) = 2 2 s + 2ξro ωro s + ωro
We wish to retrieve these results by studying the positivity of the guardian map on Q1 with the dilation integral method. Table 5 gathers the results for σ = {0, 1.5}. In both cases, εk seems to be decreasing to 0. If we can tolerate a violation percentage less than 0.01%, we can conclude to the practical robustness of the family. Note that the conclusion is rapidly reached for σ = 0 (only 3 iterations are necessary).
The controller synthesis in (Lhachemi et al., 2015) yields the following values for this specific flight condition: Kp = 1.491, Ki = 4, Kq = −9.94, ωro = 54, ξro = 0.75 97
IFAC ACA 2016 98 August 21-25, 2016. Quebec, Canada
K. Bessadi et al. / IFAC-PapersOnLine 49-17 (2016) 094–099
Table 3. Dilation integrals with σ = {0, 1.5} k 2 4 6 8 10 12
σ=0 εk 0.02 9.00 × 10−4 5.38 × 10−5 3.70 × 10−6 2.79 × 10−7 2.22 × 10−8
σ = 1.5 εk 0.0911 0.0167 3.9 × 10−3 1.06 × 10−3 3.15 × 10−4 9.93 × 10−5
θk 0.141 0.173 0.194 0.209 0.221 0.230
than 0.01% for k = 10, leading to a practical robustness assessment. Indeed, Figure 5 confirms that Q3 is well within the stable component. The upper-left and bottomright corners of Q1 are unstable, which is translated into an increasing εk and above all, a conditioner θk close to 1; the dilation integral method ensures the instability of Q1 . With regards to Q2 (for which the upper-left corner is unstable), εk does not seem to decrease under 4% but the conditioner increasing to 1 leads us to a practical instability conclusion.
θk 0.302 0.360 0.397 0.425 0.446 0.464
For the unstable cases, Table 4 shows that for σ = 1.75, εk hardly decreases while the conditioner θk seems to come close to 1. From a practical point of view, the family is unstable relative to σ = 1.75; this is indeed confirmed in Figure 4 where the upper-left corner of Q1 is unstable. The same conclusion is reached for σ = 2 as the conditioner tends rapidly to 1.
Table 5. Dilation integrals with ξ = 0.62 k 2 4 6 8 10 12 14 16 18 20
Table 4. Dilation integrals with σ = {1.75, 2} k 2 4 6 8 10 12 14 16 18 20 22
σ = 1.75 εk θk 0.2196 0.469 0.1010 0.564 0.0600 0.626 0.0412 0.671 0.0312 0.707 0.0252 0.736 0.0213 0.760 0.0188 0.780 0.0170 0.797 0.0158 0.813 0.0177 0.833
σ=2 εk θk 0.817 0.904 0.861 0.963 0.872 0.977 0.877 0.984 0.880 0.987 0.882 0.989
θk 0.5 0.7154 0.8098 0.8579 0.8868 0.9060 0.9196 0.9298 0.9377 0.9440
εk 0.1112 0.0549 0.0422 0.0401 0.0407 0.0413 0.0417 0.0420 0.0422 0.0425
Q2
θk 0.3335 0.484 0.5900 0.6688 0.7261 0.7668 0.7970 0.8203 0.8389 0.8540
Q3 εk 0.0385 0.0044 0.0007 0.00014 3.00 × 10−5 6.00 × 10−6 1.60 × 10−6 3.98 × 10−7 1.52 × 10−7 6.26 × 10−8
θk 0.1963 0.2587 0.3003 0.3302 0.3530 0.3710 0.3857 0.3981 0.4181 0.54935
An approach for analysing the generalized stability of uncertain LTI systems has been presented. This approach is based upon guardian maps and the dilation integral method to evaluate the ‘practical’ positivity of multivariate polynomials. It has been applied to the robustness assessment of a flight control law against mass and center of gravity variations. While the simulation results appear to be of some help for the simple two-parameter case we chose to tackle, more studies need to be done involving aircraft models with a higher number of parameters (e.g., flight altitude and airspeed, aerodynamic coefficients, etc.). Unlike the two-parameter case for which we were able to easily confront the dilation approach results with the true ones obtained from two-dimensional vanishing guardian maps, such comparisons will not be easily available for say four parameters or more. The notion of practical positivity and hence practical robustness, while appealing from a complexity standpoint, does bring about a trade-off to be made between various instability risk levels that one can tolerate. A related issue not being addressed by dilation integrals and practical positivity is the lack of information of the location of Qbad relative to nominal. Nevertheless, such an approach could be used to discard rapidly subsets of the parameter space where stability, or lack thereof, is ensured.
Once again, the nominal case (∆x, ∆m) = (0, 0) is stable, i.e., all the eigenvalues of Acl have a damping ratio greater than 0.62. From Figure 5, one can then conclude that Q1 and Q2 are unstable sets, while Q3 is a stable set. ξ = 0.62
1500 Stability boundary Q1 boundary Q2 boundary Q3 boundary
500
∆m (kg)
Q1
5. CONCLUSION
Our second objective is to study the stability of Acl relative to a damping constraint of ξ = 0.62 with respect to three different sets, Q1 = [−0.15, 0.15] × [−1000, 1000], Q2 = [−0.1, 0.1] × [−750, 750] and Q3 = [−0.05, 0.05] × [−500, 500]. In the following, this peculiar choice of damping ratio ensures interesting conclusions on the dilation integral results.
1000
εk 0.2498 0.2619 0.2821 0.2936 0.3010 0.3059 0.3096 0.3124 0.3145 0.3163
0
-500
REFERENCES -1000
-1500 -0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
Barmish, B.R., Shcherbakov, P.S., Ross, S.R., and Dabbene, F. (2009). On positivity of polynomials: The dilation integral method. Automatic Control, IEEE Transactions on, 54(5), 965–978. Dabbene, F. and Shcherbakov, P. (2007). Fast computation of dilation integrals for robustness analysis. In American Control Conference, 2007. ACC’07, 2660– 2665. IEEE.
0.2
∆x (m)
Fig. 5. Robustness analysis against damping ξ = 0.62 The dilation integral method gives a straightforward conclusion for the set Q3 as the volume of violation εk is less 98
IFAC ACA 2016 August 21-25, 2016. Quebec, Canada
K. Bessadi et al. / IFAC-PapersOnLine 49-17 (2016) 094–099
Lhachemi, H., Saussi´e, D., and Zhu, G. (2015). A structured H∞ -based optimization approach for integrated plant and self-scheduled flight control system design. Aerospace Science and Technology, 45, 30–38. Saussi´e, D., Saydy, L., and Akhrif, O. (2006). Longitudinal flight control design with handling quality requirements. Aeronautical Journal, 110(1111), 627–637. Saydy, L., Tits, A.L., and Abed, E.H. (1990). Guardian maps and the generalized stability of parametrized families of matrices and polynomials. Mathematics of Control, Signals and Systems, 3(4), 345–371. Shcherbakov, P.S. and Barmish, B. (2003). On the conditioning of robustness problems. In Decision and Control, 2003. Proceedings. 42nd IEEE Conference on, volume 2, 1932–1937. IEEE. Smolyak, S.A. (1963). Quadrature and interpolation formulas for tensor products of certain classes of functions. In Dokl. Akad. Nauk SSSR, volume 4, 240–243. Stephanos, C. (1900). Sur une extension du calcul des substitutions lin´eaires. Journal de Mathematiques Pures et Appliquees, 73–128.
F11
F20
F02
Appendix A. BIALTERNATE PRODUCT Let A and B be n × n matrices, and let V n be the n(n−1) 2 tuple consisting of pairs of integers (p, q), p = 2, 3 . . . n, q = 1, . . . , p − 1, listed lexicographically. That is, V n = [(2, 1), (3, 1), (3, 2), (4, 1), (4, 2), (4, 3), . . . , (n, n − 1)] Denote by Vin the ith entry of V n and 1 b b a a f ((p, q), (r, s)) = det pr ps + det pr ps bqr bqs aqr aqs 2
where the dependence of f on A and B is kept implicit for simplicity. The bialternate product of A and B, denoted A B, is a n(n−1) × n(n−1) matrix whose ijth entry is 2 2 given by (A B)ij = f Vin , Vjn Additional results and properties are available in (Saydy et al., 1990; Stephanos, 1900). Appendix B. NUMERICAL VALUES Let Fij be the concatenation of the coefficient matrices, i.e., Aij Bij Fij = (B.1) Cij Dij For the flight condition considered in this paper, the numerical values are then:
F00
F10
F01
−1.0170 × 100 1.1740 × 100 = −2.4045 × 101 0 2.3938 × 10−3 −2.8172 × 101 = 1.3688 × 101 0 1.2008 × 10−4 −3.7252 × 10−6 = 2.7762 × 10−3 0
9.4144 × 10−1 −9.5194 × 10−3 −9.4450 × 10−1 −1.2721 × 100 −9.1252 × 10−1 3.7316 × 10−1 1 0 −4.7356 × 10−3 −5.6373 × 10−3 0 −1 −1.6400 × 10 −2.7633 × 10 6.5520 × 10−1 −1.4142 × 10−3 0 0 −6 6.0565 × 10 1.0318 × 10−6 −3.6929 × 10−6 −9.7736 × 10−7 1.4166 × 10−4 2.4296 × 10−5 0 0
99
−4.4109 × 10−7 2.8030 × 10−4 = −1.4088 × 10−4 0 −1.2494 × 10−1 7.3117 × 10−1 = −3.2277 × 100 0 −1.3059 × 10−8 1.4762 × 10−11 = −3.0173 × 10−7 0
99
2.7216 × 10−8 −8.3410 × 10−9 −9.5814 × 10−6 −2.0562 × 10−6 5.0961 × 10−6 7.6601 × 10−7 0 0 −7.1967 × 10−3 −1.2285 × 10−3 −1 −1 −1.3233 × 10 −1.3019 × 10 −1.0458 × 10−1 3.2321 × 10−2 0 0 −10 −6.4278 × 10 −1.1114 × 10−10 −1.2111 × 10−10 2.3641 × 10−11 −1.4795 × 10−8 −2.5790 × 10−9 0 0
ScienceDirect IFAC-PapersOnLine 49-17 (2016) 094–099
Robustness Assessment Assessment of of a a Load Load Factor Factor Robustness Robustness Assessment of a Load Factor Robustness Assessment of a Load Factor Flight Control Law Using Dilation Flight Control Law Using Dilation Flight Control Law Using Dilation Flight Control Law Using Dilation Integrals Integrals Integrals Integrals K. Bessadi Bessadi ∗∗∗ D. D. Saussi´ Saussi´ e ∗∗∗ L. L. Saydy Saydy ∗∗∗ K. e ∗ ∗ ∗ K. Bessadi Saussi´ e ∗ D. ∗ L. K. Bessadi D. Saussi´ e L. Saydy K. Bessadi D. Saussi´ e L. Saydy Saydy ∗ ∗ ∗ Polytechnique Montr´ eeal, QC, H3T 1J4, Canada QC, H3T 1J4, Canada ∗ Polytechnique Montr´ ∗ Polytechnique Montr´ eeeal, al, QC, H3T 1J4, Canada ∗ Polytechnique Montr´ al, QC, H3T 1J4, Canada (e-mail: {karim.bessadi,david.saussie,lahcen.saydy}@polymtl.ca) Polytechnique Montr´ al, QC, H3T 1J4, Canada (e-mail: {karim.bessadi,david.saussie,lahcen.saydy}@polymtl.ca) (e-mail: {karim.bessadi,david.saussie,lahcen.saydy}@polymtl.ca) (e-mail: {karim.bessadi,david.saussie,lahcen.saydy}@polymtl.ca) (e-mail: {karim.bessadi,david.saussie,lahcen.saydy}@polymtl.ca)
Abstract: This article presents a performance assessment approach of aa load factor flight Abstract: This article presents aa robust robust performance assessment approach of factor flight Abstract: This article presents robust performance assessment approach of aaa load load factor flight Abstract: This article presents a robust performance assessment approach of load factor flight control law with the goal of ascertaining whether the closed-loop system’s poles all lie within Abstract: This article presents a robust performance assessment approach of load factor flightaa control law with the goal of ascertaining whether the closed-loop system’s poles all lie within control law with the goal of ascertaining whether the closed-loop system’s poles all lie within control law with the goal of ascertaining whether the closed-loop system’s poles all lie within aaa given subregion of the open left half-plane despite center of gravity and mass variations. The controlsubregion law with of thethe goal of ascertaining whether the closed-loop system’s polesvariations. all lie within given open left half-plane despite center of gravity and mass The given subregion of the open left half-plane despite center of gravity and mass variations. The given subregion of the open left half-plane despite center of gravity and mass variations. The method combines a guardian map approach with dilation integrals to tackle the underlying given subregion of the open left half-plane despite center of gravity and mass variations. The method combines a guardian map approach with dilation integrals to tackle the underlying method combines guardian map approach with dilation to tackle underlying method polynomial combines a a positivity guardian test mapby approach withlatter dilation integrals to tackle the the underlying difficult relaxing the to aa integrals so-called ‘practical’ positivity one. method combines a guardian map approach with dilation integrals to tackle the underlying difficult polynomial positivity test by relaxing the latter to so-called ‘practical’ positivity one. difficult polynomial positivity test by relaxing the latter to a so-called ‘practical’ positivity one. difficult polynomial positivity test by relaxing the latter to a so-called ‘practical’ positivity one. The study is carried carried positivity out for for a a longitudinal longitudinal uncertain model of a F-16 F-16 aircraft aircraft with positivity varying center center difficult polynomial test by relaxing the latter to of a so-called ‘practical’ one. The study is out uncertain model a with varying The study is carried out for a longitudinal uncertain model of a F-16 aircraft with varying center The study is carried out for a longitudinal uncertain model of a F-16 aircraft with varying center of gravity position and mass with the purpose of testing such an approach. The study position is carriedand outmass for a with longitudinal uncertain modelsuch of a an F-16 aircraft with varying center of gravity the purpose of testing approach. of of gravity gravity position position and and mass mass with with the the purpose purpose of of testing testing such such an an approach. approach. of gravity position and mass with the purpose of testing such an approach. © 2016, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: dilation dilation integrals, integrals, generalized generalized stability, stability, guardian guardian maps, maps, robustness robustness analysis, analysis, flight flight Keywords: Keywords: dilation integrals, generalized stability, guardian maps, robustness analysis, flight Keywords: dilation generalized control, polynomial positivity. Keywords: dilation integrals, integrals, generalized stability, stability, guardian guardian maps, maps, robustness robustness analysis, analysis, flight flight control, polynomial positivity. control, polynomial positivity. control, control, polynomial polynomial positivity. positivity. 1. INTRODUCTION INTRODUCTION under consideration consideration is is deemed deemed positive positive from from aa practical practical 1. under 1. INTRODUCTION under consideration 1. under consideration is deemed positive from practical viewpoint. 1. INTRODUCTION INTRODUCTION under consideration is is deemed deemed positive positive from from aaa practical practical viewpoint. viewpoint. viewpoint. viewpoint. The present work deals broadly with the generalized In Section Section 2, 2, we we briefly briefly introduce introduce guardian guardian maps maps to to study study The present work deals broadly with the generalized In The present work deals broadly with the generalized The work broadly with generalized stability of parametrized of uncertain systems. In Section 2, briefly guardian maps The present present work deals dealsfamilies broadly with the theLTI generalized In Section 2, we briefly introduce guardian maps to study the stability of families of parametrized mastability of parametrized families of uncertain LTI systems. In Section 2, we we briefly introduce introduce guardianstate-space maps to to study study the stability of families of parametrized state-space mastability of parametrized families of uncertain LTI systems. stability of families of LTI ensures More specifically, given a controller C the stability of families of parametrized state-space mastability of parametrized parametrized families of uncertain uncertain LTI systems. systems. 0 which the stability of families of parametrized state-space matrices with respect to a given region Ω ⊂ C. Examples which ensures More specifically, given a controller C the stability of families of parametrized state-space ma0 trices with respect to a given region Ω ⊂ C. Examples 1 which ensures More specifically, given a controller C 0 which ensures More specifically, given controller C 1a the generalized stability of the system for a nominal trices with respect to a given region Ω ⊂ C. Examples 0 which ensures More specifically, given a controller C trices with respect to a given region Ω ⊂ C. Examples 0 for a nominal of guardian maps are given for simple regions of interest, the generalized stability of the system 1 trices with respect to a given region Ω ⊂ C. Examples 1 of guardian maps are given for simple regions of interest, k the generalized stability system for a 1 of the generalized stability of the system for a nominal k value of an uncertain parameter qq 000 ∈ we of maps for regions of the generalized stability of the thevector system forR a,, nominal nominal of guardian maps are given for simple regions of interest, as guardian well as as aa necessary necessary and sufficient sufficient condition provided in value of an uncertain parameter vector ∈ R we wish wish k of guardian maps are are given given for simple simple regions provided of interest, interest, 0 k as well and condition in value of an uncertain parameter vector q ∈ R wish 0 ensures k ,, we value of an uncertain parameter vector q ∈ R we wish to ascertain whether the same controller robust as well as a necessary and sufficient condition provided in value of an uncertain parameter vector q ensures ∈ R , werobust wish the as well as a necessary and sufficient condition provided in form of a positivity test involving a single inequality. to ascertain whether the same controller as well asofa necessary and sufficient condition provided in the form a positivity test involving a single inequality. to ascertain whether the same controller ensures robust to ascertain whether the same controller ensures robust stability with respect to Ω when the uncertain parameter the form of a positivity test involving a single inequality. to ascertain whether the same controller ensures robust the form of a positivity test involving a single inequality. In Section 3, the dilation integral method is presented. stability with respect to Ω when the uncertain parameter the form of a positivity test involving a single inequality. In Section 3, the dilation integral method is presented. k stability with respect to Ω when the uncertain parameter stability with within respecta to Ω when when theQ uncertain parameter aa broad vector varies given subset ⊂ In Section 3, dilation integral method presented. stability with respect to Ω the uncertain parameter In Section 3, the dilation integral method is presented. An upper bound Vol(Q is obtained dilation For vector varies within a given subset Q ⊂R Rkkkk ... For In Section 3, the the on dilation integral method is isvia presented. bad ) An upper bound on Vol(Q is obtained via dilation For aaa broad broad vector within a Q bad ) For broad vector varies within a given subset Q ⊂ R class of ofvaries regions and uncertain uncertain LTI systems, systems, the guardian An upper bound on Vol(Q ) is obtained via dilation Forguardian broad integral vector varies within a given given subset subset Q⊂ ⊂R R ..the bad An upper bound on Vol(Q ) is obtained via dilation computations and a conditioning parameter inclass regions and LTI bad An upper bound on Vol(Q ) is obtained via dilation bad integral computations and a conditioning parameter inclass of regions and uncertain LTI systems, the guardian class of regions and uncertain LTI systems, the guardian map approach (Saydy et al., 1990) yields a necessary integral computations and a conditioning parameter inclass of regions (Saydy and uncertain LTI systems, the guardian integral integral computations and a conditioning parameter introduced. In Section4, we combine the two approaches to map approach et al., 1990) yields a necessary computations and a conditioning parameter introduced. In Section4, we combine the two approaches to map approach (Saydy et al., 1990) yields a necessary map approach (Saydy et al., 1990) yields a necessary and sufficiant condition for robustness in the form of troduced. In Section4, we combine the two approaches to map approach (Saydy et al., 1990) yields a necessary troduced. In Section4, we combine the two approaches to analyse the robust stability of a longitudinal flight control and sufficiant condition for robustness in the form of troduced. In Section4, we combine the two approaches to analyse the robust stability of a longitudinal flight control and sufficiant condition for robustness in the form of and sufficiant condition for robustness in the form of a positivity test of a single multivariable polynomial, analyse the robust stability of a longitudinal flight control and sufficiant condition for robustness in the form of analyse the robust stability of a longitudinal flight control system. More specifically, we analyse the robust stability a positivity test of a single multivariable polynomial, analyse the robust stability of a longitudinal flight control system. More specifically, we analyse the robust stability a positivity test of a single multivariable polynomial, a positivity test of single namely, decide whether or not: multivariable More specifically, we the robust stability a positivity test of a a or single multivariable polynomial, polynomial, system. system. More specifically, we analyse the robust stability of aa nominal factor flight control law to see whether namely, decide whether not: system. Moreload specifically, we analyse analyse the robust stability of nominal load factor flight control law to see whether namely, decide whether or not: namely, decide whether or not: of a nominal load factor flight control law to see whether namely, decide whether or not: of a nominal load factor flight control law to see whether the closed-loop poles preserve a certain negativity margin of a nominal load factor flight control law to see whether ν(q) > 0 ∀q ∈ Q the closed-loop poles preserve a certain negativity margin ν(q) > 0 ∀q ∈ Q the closed-loop poles preserve a certain negativity margin ν(q) > 0 ∀q ∈ Q the closed-loop poles preserve a certain negativity margin and a certain damping, despite the variations of center of ν(q) > >0 0 ∀q ∀q ∈ ∈Q Q the closed-loop poles preserve a certain negativity margin ν(q) and a certain damping, despite the variations of center of where ν(q) is a polynomial constructed from the guardian and a certain damping, despite the variations of center of where ν(q) is a polynomial constructed from the guardian and a certain damping, despite the variations of center of gravity and mass. Various assessments using the dilation and a certain damping, despite the variations of center of where ν(q) is a polynomial constructed from the guardian gravity and mass. Various assessments using the dilation where ν(q) is a polynomial constructed from the guardian map associated with the problem under study (cf. Section where ν(q) is a polynomial constructed from the guardian gravity and mass. Various assessments using the dilation map associated with the problem under study (cf. Section gravity and mass. Various assessments using the dilation integrals are presented. gravity and mass. Various assessments using the dilation map associated with the problem under study (cf. Section integrals are presented. map associated with with the the problem problem under under study study (cf. (cf. Section Section integrals are presented. 2). map associated 2). integrals are are presented. presented. integrals 2). 2). 2). Owing to the NP-hardness of tackling such a proposition Owing to the NP-hardness of tackling such a proposition Owing to NP-hardness of such proposition Owing to the the NP-hardness of tackling tackling such a a the proposition 2. GUARDIAN in general, Barmish et al. (2009) introduced concept Owing to the NP-hardness of tackling such a proposition 2. GUARDIAN MAPS MAPS in general, Barmish et al. (2009) introduced the concept 2. in general, Barmish et al. (2009) introduced the 2. GUARDIAN MAPS in general, Barmish et al. (2009) introduced the concept of general, ‘pratical’Barmish positivity to alleviate alleviate the complexity complexity issue 2. GUARDIAN GUARDIAN MAPS MAPS in et al. (2009) introduced the concept concept of ‘pratical’ positivity to the issue of ‘pratical’ positivity to alleviate the complexity issue of ‘pratical’ positivity to alleviate the complexity issue associated with precisely deciding true positivity of a of ‘pratical’with positivity to alleviate the complexity issue We are here interested in stability sets of the form: associated precisely deciding true positivity of a are here interested in stability sets of the form: associated with precisely deciding true positivity of aaa We associated with precisely In deciding trueet positivity ofthe We are here interested in stability sets of the form: multivariable polynomial. (Barmish al., 2009), associated with precisely deciding true positivity of We are here interested in stability sets of the form: n×n We are here interested in stability sets of the multivariable polynomial. In (Barmish et al., 2009), the n×n : σ(A) ⊂ Ω} form: S(Ω) = {A ∈ R (1) multivariable polynomial. In (Barmish et al., 2009), the multivariable polynomial. In (Barmish et al., 2009), the S(Ω) = {A ∈ R : σ(A) ⊂ Ω} (1) n×n authors introduced the method of dilation integrals, multivariable polynomial. In (Barmish et al., 2009), the n×n authors introduced the method of dilation integrals, the S(Ω) = {A ∈ R ⊂ Ω} (1) n×n :: σ(A) S(Ω) = {A ∈ R σ(A) ⊂ Ω} (1) authors introduced the method of dilation integrals, the S(Ω) = {A ∈ R : σ(A) ⊂ Ω} (1) authors introduced the method of dilation integrals, the purpose introduced of which which is is to to decide whether a given givenintegrals, polynomial, Ω is an open subset of the complex plane, and σ(A) authors thedecide method of dilation the where purpose of whether a polynomial, where Ω is an open subset of the complex plane, and σ(A) purpose of which is to decide whether a given polynomial, where Ω is an open subset of the complex plane, and σ(A) purpose of which is to decide whether a given polynomial, where Ω is an open subset of the complex plane, and σ(A) such as ν(q), is practically positive for (almost) all q ∈ denotes the spectrum of A. Such sets S(Ω) will be referred purpose of which is to decidepositive whetherfor a given polynomial, where Ωthe is an open subset the sets complex plane, σ(A) such as ν(q), is practically (almost) all qq ∈ denotes spectrum of A.ofSuch S(Ω) will beand referred such as is positive for all ∈ the of Such sets S(Ω) will such as ν(q), ν(q), is practically practically positive for (almost) (almost) all ∈ denotes denotes the spectrum spectrum of A. A.sets, Suchand setsthus S(Ω)represent will be be referred referred Q. More specifically, the method yields a means to to as generalized stability the set such as ν(q), is practically positive for (almost) all qqtest ∈ denotes the spectrum of A. Such sets S(Ω) will be referred Q. More specifically, the method yields a means to test to as generalized stability sets, and thus represent the set Q. More specifically, the method yields a means to test to as generalized stability sets, and thus represent the set Q. More specifically, the method yields a means to test to as generalized stability sets, and thus represent the set whether the set of parameters for which the polynomial of all matrices which are stable relative to Ω, i.e., which Q. More specifically, the method yields a means to test to as generalized stability sets, and thus represent the set whether the set of parameters for which the polynomial of all matrices which are stable relative to Ω, i.e., which whether the set of parameters for which the polynomial of all matrices which are stable relative to Ω, i.e., which whether the set of parameters for which the polynomial of all matrices which are stable relative to Ω, i.e., which is not positive (so-called Q ) has a volume below a have all their eigenvalues in Ω. Guardian maps are defined whether the set (so-called of parameters which the polynomial bad for of all matrices which are stable relative to Ω, i.e., which is not positive Q ) has a volume below a have all their eigenvalues in Ω. Guardian maps are defined is (so-called Q ))) has a below have all their eigenvalues in Ω. Guardian maps are defined bad is not not positive positive (so-called Qbad has case a volume volume below aaa next. have all certain tolerable threshold , in which the polynomial bad is not positive (so-called Q has a volume below have all their their eigenvalues eigenvalues in in Ω. Ω. Guardian Guardian maps maps are are defined defined bad certain tolerable threshold , in which case the polynomial next. certain tolerable threshold , in which case the polynomial next. certain tolerable tolerable threshold threshold , , in in which which case case the the polynomial polynomial next. next. n×n certain n×n into C, Definition 1. (Saydy et al., 1990) Let ν map R Definition 1. (Saydy et al., 1990) Let νν¯ map R C, n×n n×n into Definition 1. (Saydy et al., 1990) Let map R into C, Definition 1. et al., 1990) Let νν¯ map we say that guards S if for all A ∈ S(Ω): Definition 1.νν (Saydy (Saydy et al., 1990) Let map R Rn×n into into C, C, work was supported by NSERC under Grant RGPIN-2014 This we say that guards S if for all A ∈ S(Ω): ¯ This work was supported by NSERC under Grant RGPIN-2014 ¯ we say that ν guards S if for all A ∈ S(Ω): ¯ we say that ν guards S if for all A ∈ S(Ω): This work was supported by NSERC under Grant RGPIN-2014 we say that ν guards S if for all A ∈ S(Ω): 03942 and RGPIN-2012-122106. This work was supported by NSERC under Grant RGPIN-2014This work was supported by NSERC under Grant RGPIN-201403942 and RGPIN-2012-122106. ν(A) = = 00 ⇔ ⇔A A∈ ∈ ∂S(Ω) ∂S(Ω) (2) 1 ν(A) (2) 03942 and RGPIN-2012-122106. In of 03942 and RGPIN-2012-122106. 1 ν(A) = ⇔ A ∈ ∂S(Ω) (2) 03942 andsense RGPIN-2012-122106. In the the sense of placing placing closed-loop closed-loop poles poles within within a a given given region region ν(A) = 000 ⇔ A ∈ ∂S(Ω) (2) 1 ν(A) = ⇔ A ∈ ∂S(Ω) (2) 1 ¯ In the sense of placing closed-loop poles within a given region 1 Here S S¯ denotes denotes closure closure of of the the set set S S and and ∂S ∂S its its boundary. boundary. Ω In ⊂ C. In C. the sense sense of of placing placing closed-loop closed-loop poles poles within within a a given given region region the Here Ω ⊂ ¯ ¯ Here S denotes closure of the set S and ∂S its boundary. Ω ⊂ C. ¯ Here S denotes closure of the set S and ∂S its boundary. ⊂ C. Ω Here S denotes closure of the set S and ∂S its boundary. Ω ⊂ C. Copyright 2016 IFAC 94 2405-8963 Copyright©© ©2016, 2016IFAC IFAC (International Federation of Automatic Control) 94 Hosting by Elsevier Ltd. All rights reserved. Copyright © 2016 IFAC 94 Copyright 2016 IFAC 94 Copyright ©under 2016responsibility IFAC 94 Control. Peer review © of International Federation of Automatic 10.1016/j.ifacol.2016.09.017
K. Bessadi et al. / IFAC-PapersOnLine 49-17 (2016) 094–099
A(q1 , q2 ) =
• Hurwitz stability :
Let Ω = C◦− (Fig. 1(a)). Then S(Ω) is guarded by:
(3) νH : A → det(A) det(A I) where denotes the bialternate product (cf. Appendix A).
A corresponding guardian maps is given by νΩ = νσ νξ where (Saydy et al., 1990): (8) νσ (A) = det(A + σI) det((A + σI) I) 2 2 νξ (A) = det(A) det(A I + (1 − 2ξ )A A) (9) We obtain : 5 νΩ (A) = − (q1 + 2)(9q1 − 6q2 + 16) 4 × (q1 + 6q2 )(25q12 + 62q1 + 12q2 + 36) (10)
• Negativity margin:
Let Ω be the σ-shifted open left half-plane (Fig. 1.b). Then S(Ω) is guarded by : νσ : A → det(A + σI) det((A + σI) I) (4) • Damping ratio:
Let Ω be the conic sector with inner angle 2θ (Fig. 1.c). Then S(Ω) is guarded by: νξ : A → det(A) det(A2 I + (1 − 2ξ 2 )A A) with ξ = cos θ.
For q 0 = [−2.75, −2.5] , the matrix A(q 0 ) has eigenvalues −3.875 ± j1.654 and is therefore nominally Ω−stable. Furthermore, νΩ (A(q 0 )) < 0. Therefore, according to Theorem 1, the whole family is Ω−stable within Q1 if, and only if, νΩ (A(q)) changes no sign in Q1 , that is νΩ (q) = (q1 + 2)(9q1 − 6q2 + 16) × (q1 + 6q2 )(25q12 + 62q1 + 12q2 + 36) > 0 (11) for all q ∈ Q1 .
(5)
• Schur stability :
Let Ω be the circle of radius ω (Fig. 1.d). Then S(Ω) is guarded by: νω : A → det(A + ωI) det(A − ωI) det(A A − ω 2 I I) (6)
In Figs 2.a and 2.b, setting νΩ (r) = 0 divides the parameter space into many components that are either stable or unstable. One needs to test only one vector in each component to see what prevails (Saussi´e et al., 2006). The green component, being nominally Ω-stable, corresponds in fact to the largest set containing the nominal value for which Ω-stability is preserved. In particular, we have robust Ω-stability within Q1 but not with Q2 = [−3, −2]× [−4, −2] as the guardian map vanishes within the latter parameter set.
A systematic way of constructing guardian maps for various Ω regions can be found in (Saydy et al., 1990). Im
Im
Im
ω 0 Re
(a)
−
σ
Re
(b)
θ
(c)
Re
0 1 (7) 6q2 + q1 5q1 + 6 We wish to study the stability of this family relative to the σ−shifted left half-plane with σ = 2, and to the conic √ sector with θ = π4 ) (ξ = 22 ), i.e., relative to Ω being the intersection of those in Figs 1.b and 1.c.
2.1 Some examples of guardian maps
Im
95
Re
(d)
2
2
1
1
0
q2
Fig. 1. Useful stability sets Let {A(q) : q ∈ Q} be a continuous family of real n × n matrices depending on an uncertain parameter q known to belong to pathwise connected subset Q ⊂ Rk . Let us also assume that this family is nominally stable relative to a given Ω, that is A(q 0 ) ∈ S(Ω) for some q 0 ∈ Q. While the stability of this family amounts in general to testing a number of inequalities for each A(q) (depending on n and Ω), the result below states that only one amongst these needs to be tested for robust stability. Theorem 1. (Saydy et al., 1990) Let S(Ω) be guarded by the map νΩ . The family {A(q) : q ∈ Q} is stable relative to Ω if and only if:
0
q2
IFAC ACA 2016 August 21-25, 2016. Quebec, Canada
-1
-1
-2
-2
-3
-3
-4
-4 -3
-2.5
-2
-1.5
-1
-0.5
0
-3
-2.5
-2
q1
(a) Stable set Q1
-1.5
-1
-0.5
0
q1
(b) Unstable set Q2
Fig. 2. Stable and unstable components 3. DILATION INTEGRAL METHOD Here, we give a brief overview of the dilation integral method of (Barmish et al., 2009) Let f (q) be a multivariable polynomial and let Qbad = {q ∈ Q : f (q) ≤ 0}, the subset of Q for which f (q) is not positive. In order for f (q) > 0, it is thus equivalent Qbad = ∅. However, as this is computationnally very hard in general, the dilation integral method proposes to judge the positivity from a practical point of view by attempting to show that the percentage of positivity violation as bad ) typified by the volume ratio Vol(Q Vol(Q) is below a certain value deemed tolerable.
(1) it is nominally stable, i.e. A(q0 ) ∈ S(Ω) for some q0 ∈ Q ; and (2) ν(A) := νΩ (A(q))νΩ (A(q0 )) > 0, for all q ∈ Q. 2.2 Example Let A(q1 , q2 ) be a family of matrices depending on two parameters −3 ≤ q1 ≤ −2.5 and −3 ≤ q2 ≤ −2 (Q1 = [−3, −2.5] × [−3, −2]) : 95
IFAC ACA 2016 96 August 21-25, 2016. Quebec, Canada
K. Bessadi et al. / IFAC-PapersOnLine 49-17 (2016) 094–099
with increasing k, the symbolic computation can be intractable in a decent time while quadrature formulas can efficiently alleviate the burden.
The dilation integral associated with the pair (f (q), Q) is given by: 1 k (1 − αf (q)) dq (12) k (α) = Vol(Q) Q with α ≥ 0 and k being a positive even integer.
3.1 The conditioner In order to stop the iterations and declare the pair (f, Q) to be non positive from a practical point of view, a new parameter is introduced (Shcherbakov and Barmish, 2003). The conditioner θ is defined as the maximum percentage variation of f (q) about the midpoint of its range. That is, with
For every α positive, the inequality below gives an upper bound for the violation ratio: Vol(Qbad ) ≤ k (α) (13) Vol(Q) It turns out that this upper bound achieves its minimum so that: Vol(Qbad ) . ≤ k = min k (α) (14) α≥0 Vol(Q) and we have the following properties: Theorem 2. (Barmish et al., 2009) For the given pair (f, Q) with associated dilation integrals 1 . k k (α) = (1 − αf (q)) dq (15) Vol(Q) Q defined for positive even integers k, their minima, . k = min k (α)
fmin = minf (q), and fmax = maxf (q) q∈Q
α≥0
To estimate the conditioner’s value, one relies on the fact that . 1/k θ k = k ≤ θ (17)
are attained. Moreover, the following conditions hold: (1) If (f, Q) is a positive pair, then k → 0. (2) If k → 0, then (f, Q) is a non-negative pair with Vol(Qbad ) = 0. Example 1. Let us consider the polynomial (Barmish et al., 2009)
So to stop the computations when one is unable to reach values of k below the acceptable threshold , we use the sequence θk as a revealing indicator of an imminent risk of losing positivity when its value is close to 1. Theorem 3. (Barmish et al., 2009) If (f, Q) is a positive pair, it follows that lim θk = θ (18)
f (x) = 1 + x21 x22 (x21 + x22 − 3) which as it turns out is positive for any Q = [−r, r]2 with 0 < r < 1.
k→∞
For r = 0.75 and 6 (α), we obtain the following: 0.75 0.75 1 6 6 (α) = (1 − αf (x)) dx1 dx2 . (2 · 0.75)2 −0.75 −0.75 ≈ 0.7054α6 − 4.40175α5 + 11.5030α4 − 16.13601α3 + 12.8394α2 − 5.5096α + 1.0.
If (f, Q) is not a positive pair, then lim θk = 1 k→∞
(19)
Moreover, in both cases, the sequence of estimates θk is non-decreasing. Example 2. Let us consider the multivariate polynomial f (q) = q12 (q12 + 1) q22 ((q1 + 0.2)2 + (q2 − 0.3)2 − 0.22 ) with (q1 , q2 ) ∈ [−0.5, 0.5] × [−0.5, 0.5].
The minimum of this function is 6 = 0.00011355 which in turn ensures that positivity is violated on a set whose volume is at worst equal to 0.00011355. If greater accuracy is required, similar integrations may be carried out to try to achieve smaller values (Table 2).
For this example, it is easy to compute the violation ratio: Vol(Qbad ) = 0.22 π ≈ 0.1256 Vol(Q) The results obtained with the dilation integral are given in the table below: Table 2. Utility of θk
Table 1. Convergence de k k 2 4 6 8 10 12 14 16
q∈Q
(fmax − fmin ) (16) θ= | fmax + fmin | In case of nominally positive polynomials f (q), (i.e., f (q0 ) > 0 for some q0 ∈ Q), it is clear that if 0 ≤ θ < 1 then positivity is ensured over the set Q. But if θ → 1, then it would mean that there are parameter configurations where the risk of violating the positivity becomes imminent (Barmish et al., 2009).
k 0.0149 0.0011010 0.0001135 1.36 × 10−5 1.77 × 10−6 2.44 × 10−7 3.45 × 10−8 5.30 × 10−9
k 2 4 6 8 10 12 14 16 .. . 26
Remark 1. To compute symbolic expressions k (α) of the dilatation integrals and the corresponding minimum k quadrature formulas such as Smolyak’s may be used (Smolyak, 1963; Dabbene and Shcherbakov, 2007). Indeed, 96
k 0.85 0.80 0.77 0.75 0.73 0.71 0.70 0.69 .. . 0.65
θk 0.92 0.94 0.95 0.96 0.97 0.97 0.97 0.97 .. . 0.98
IFAC ACA 2016 August 21-25, 2016. Quebec, Canada
K. Bessadi et al. / IFAC-PapersOnLine 49-17 (2016) 094–099
Note that the sequence k is decreasing. However since the conditioner is getting closer to 1, we can stop the calculations and declare that positivity is not practically ensured.
Finally, the closed-loop state-space model can be written as: x˙ = Acl x + Bcl nzc nz = Ccl x and the matrices are: where x = α q δe xi δec δ˙ec A B 0 0
4. FLIGHT CONTROL ROBUSTNESS ANALYSIS
2×1
This section illustrates the robustness analysis of a load factor flight control law extracted from (Lhachemi et al., 2015). For a fixed flight condition, the uncertainties are the mass variation ∆m (in kg) and the center of gravity location ∆x (in m) around nominal values.
Acl =
(23) Bcl = [0 0 0 1 0 0] Ccl = [Cnz Dnz 0 0 0] (24) where Cnz and Dnz denote the row vectors corresponding to nz output. The dependencies in ∆x and ∆m have been dropped to alleviate notation. Nevertheless the jacobian matrix Acl is still a quadratic polynomial matrix in ∆x and ∆m. 4.2 Robustness analysis Our first objective is to study the stability of Acl (∆x, ∆m) relative to the σ−shifted left half-plane on the set Q1 = [−0.15, 0.15] × [−1000, 1000]. Vanishing the corresponding guardian map for σ = {0, 1.5, 1.75, 2} leads to figure 4. σ=0
where Aij are constant coefficient matrices. Numerical values can be found in Appendix B. The elevator dynamics is modeled as a first order transfer function: δe 1 = Fact (s) = δ ec τs + 1 where δec denotes the commanded elevator deflection with time constant τ = 0.05.
Kp
−
+ −
Fro
Fact
δe
1000
1000
500
500
0
0
-500
-500
-1000
-1000
-0.1
0
0.1
-1500 -0.2
0.2
-0.1
∆x (m) σ = 1.75
nz A/C q
Kq
1500
1000
1000
500
500
0
-500
-1000
-1000 0
∆x (m)
0.2
0.1
0.2
0
-500
-0.1
0.1
σ=2
1500
-1500 -0.2
0
∆x (m)
∆m (kg)
+
1500
-1500 -0.2
∆m (kg)
Ki /s
δ ec
σ = 1.5
1500
∆m (kg)
∆m (kg)
0≤i,j≤2 i+j≤2
−
(22)
We consider the longitudinal model of a F-16 aircraft linearized at Mach 0.9 and altitude 5000 m. The classical short period approximation is given by: α˙ α q˙ A(∆x, ∆m) B(∆x, ∆m) q (20) n z C(∆x, ∆m) D(∆x, ∆m) δe q where the state variables are the angle of attack α and the pitch rate q, the input variable is the elevator deflection δe and the output measurements are q and the load factor nz . The state-space matrices (A, B, C, D) are approximated by quadratic polynomial matrices in ∆x and ∆m. For instance, A(∆x, ∆m) = Aij ∆xi ∆mj (21)
+
2×2
01×2 −1/τ 0 1/τ 0 −Cnz −Dnz 0 01×2 0 0 0 1 01×2 2 2 D 2 2 − Kp Kq ωro C −Kp ωro nz Ki ωro −ωro −2ξωro
4.1 Model
nz c
97
0.1
0.2
-1500 -0.2
-0.1
0
∆x (m)
Fig. 4. Stable et unstable sets for different values of σ
Fig. 3. Load factor control law
The nominal case (∆x, ∆m) = (0, 0) is stable for the four values of σ. For σ = {0, 1.5}, the guardian map does not vanish on Q1 , then Acl (∆x, ∆m) is stable for all (∆x, ∆m) ∈ Q1 . Whereas, for σ = {1.75, 2}, the guardian map does vanish on Q1 and one can then deduce that the family is unstable.
The control law depicted in Fig. 3 is comprised of: • a PI controller with integrator Ki /s on the tracking error nzc − nz and proportionnal feedback Kp on the load factor nz • a proportionnal feedback Kq on the pitch rate q • a roll-off filter Fro (s) to limit the controller bandwidth: 2 ωro Fro (s) = 2 2 s + 2ξro ωro s + ωro
We wish to retrieve these results by studying the positivity of the guardian map on Q1 with the dilation integral method. Table 5 gathers the results for σ = {0, 1.5}. In both cases, εk seems to be decreasing to 0. If we can tolerate a violation percentage less than 0.01%, we can conclude to the practical robustness of the family. Note that the conclusion is rapidly reached for σ = 0 (only 3 iterations are necessary).
The controller synthesis in (Lhachemi et al., 2015) yields the following values for this specific flight condition: Kp = 1.491, Ki = 4, Kq = −9.94, ωro = 54, ξro = 0.75 97
IFAC ACA 2016 98 August 21-25, 2016. Quebec, Canada
K. Bessadi et al. / IFAC-PapersOnLine 49-17 (2016) 094–099
Table 3. Dilation integrals with σ = {0, 1.5} k 2 4 6 8 10 12
σ=0 εk 0.02 9.00 × 10−4 5.38 × 10−5 3.70 × 10−6 2.79 × 10−7 2.22 × 10−8
σ = 1.5 εk 0.0911 0.0167 3.9 × 10−3 1.06 × 10−3 3.15 × 10−4 9.93 × 10−5
θk 0.141 0.173 0.194 0.209 0.221 0.230
than 0.01% for k = 10, leading to a practical robustness assessment. Indeed, Figure 5 confirms that Q3 is well within the stable component. The upper-left and bottomright corners of Q1 are unstable, which is translated into an increasing εk and above all, a conditioner θk close to 1; the dilation integral method ensures the instability of Q1 . With regards to Q2 (for which the upper-left corner is unstable), εk does not seem to decrease under 4% but the conditioner increasing to 1 leads us to a practical instability conclusion.
θk 0.302 0.360 0.397 0.425 0.446 0.464
For the unstable cases, Table 4 shows that for σ = 1.75, εk hardly decreases while the conditioner θk seems to come close to 1. From a practical point of view, the family is unstable relative to σ = 1.75; this is indeed confirmed in Figure 4 where the upper-left corner of Q1 is unstable. The same conclusion is reached for σ = 2 as the conditioner tends rapidly to 1.
Table 5. Dilation integrals with ξ = 0.62 k 2 4 6 8 10 12 14 16 18 20
Table 4. Dilation integrals with σ = {1.75, 2} k 2 4 6 8 10 12 14 16 18 20 22
σ = 1.75 εk θk 0.2196 0.469 0.1010 0.564 0.0600 0.626 0.0412 0.671 0.0312 0.707 0.0252 0.736 0.0213 0.760 0.0188 0.780 0.0170 0.797 0.0158 0.813 0.0177 0.833
σ=2 εk θk 0.817 0.904 0.861 0.963 0.872 0.977 0.877 0.984 0.880 0.987 0.882 0.989
θk 0.5 0.7154 0.8098 0.8579 0.8868 0.9060 0.9196 0.9298 0.9377 0.9440
εk 0.1112 0.0549 0.0422 0.0401 0.0407 0.0413 0.0417 0.0420 0.0422 0.0425
Q2
θk 0.3335 0.484 0.5900 0.6688 0.7261 0.7668 0.7970 0.8203 0.8389 0.8540
Q3 εk 0.0385 0.0044 0.0007 0.00014 3.00 × 10−5 6.00 × 10−6 1.60 × 10−6 3.98 × 10−7 1.52 × 10−7 6.26 × 10−8
θk 0.1963 0.2587 0.3003 0.3302 0.3530 0.3710 0.3857 0.3981 0.4181 0.54935
An approach for analysing the generalized stability of uncertain LTI systems has been presented. This approach is based upon guardian maps and the dilation integral method to evaluate the ‘practical’ positivity of multivariate polynomials. It has been applied to the robustness assessment of a flight control law against mass and center of gravity variations. While the simulation results appear to be of some help for the simple two-parameter case we chose to tackle, more studies need to be done involving aircraft models with a higher number of parameters (e.g., flight altitude and airspeed, aerodynamic coefficients, etc.). Unlike the two-parameter case for which we were able to easily confront the dilation approach results with the true ones obtained from two-dimensional vanishing guardian maps, such comparisons will not be easily available for say four parameters or more. The notion of practical positivity and hence practical robustness, while appealing from a complexity standpoint, does bring about a trade-off to be made between various instability risk levels that one can tolerate. A related issue not being addressed by dilation integrals and practical positivity is the lack of information of the location of Qbad relative to nominal. Nevertheless, such an approach could be used to discard rapidly subsets of the parameter space where stability, or lack thereof, is ensured.
Once again, the nominal case (∆x, ∆m) = (0, 0) is stable, i.e., all the eigenvalues of Acl have a damping ratio greater than 0.62. From Figure 5, one can then conclude that Q1 and Q2 are unstable sets, while Q3 is a stable set. ξ = 0.62
1500 Stability boundary Q1 boundary Q2 boundary Q3 boundary
500
∆m (kg)
Q1
5. CONCLUSION
Our second objective is to study the stability of Acl relative to a damping constraint of ξ = 0.62 with respect to three different sets, Q1 = [−0.15, 0.15] × [−1000, 1000], Q2 = [−0.1, 0.1] × [−750, 750] and Q3 = [−0.05, 0.05] × [−500, 500]. In the following, this peculiar choice of damping ratio ensures interesting conclusions on the dilation integral results.
1000
εk 0.2498 0.2619 0.2821 0.2936 0.3010 0.3059 0.3096 0.3124 0.3145 0.3163
0
-500
REFERENCES -1000
-1500 -0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
Barmish, B.R., Shcherbakov, P.S., Ross, S.R., and Dabbene, F. (2009). On positivity of polynomials: The dilation integral method. Automatic Control, IEEE Transactions on, 54(5), 965–978. Dabbene, F. and Shcherbakov, P. (2007). Fast computation of dilation integrals for robustness analysis. In American Control Conference, 2007. ACC’07, 2660– 2665. IEEE.
0.2
∆x (m)
Fig. 5. Robustness analysis against damping ξ = 0.62 The dilation integral method gives a straightforward conclusion for the set Q3 as the volume of violation εk is less 98
IFAC ACA 2016 August 21-25, 2016. Quebec, Canada
K. Bessadi et al. / IFAC-PapersOnLine 49-17 (2016) 094–099
Lhachemi, H., Saussi´e, D., and Zhu, G. (2015). A structured H∞ -based optimization approach for integrated plant and self-scheduled flight control system design. Aerospace Science and Technology, 45, 30–38. Saussi´e, D., Saydy, L., and Akhrif, O. (2006). Longitudinal flight control design with handling quality requirements. Aeronautical Journal, 110(1111), 627–637. Saydy, L., Tits, A.L., and Abed, E.H. (1990). Guardian maps and the generalized stability of parametrized families of matrices and polynomials. Mathematics of Control, Signals and Systems, 3(4), 345–371. Shcherbakov, P.S. and Barmish, B. (2003). On the conditioning of robustness problems. In Decision and Control, 2003. Proceedings. 42nd IEEE Conference on, volume 2, 1932–1937. IEEE. Smolyak, S.A. (1963). Quadrature and interpolation formulas for tensor products of certain classes of functions. In Dokl. Akad. Nauk SSSR, volume 4, 240–243. Stephanos, C. (1900). Sur une extension du calcul des substitutions lin´eaires. Journal de Mathematiques Pures et Appliquees, 73–128.
F11
F20
F02
Appendix A. BIALTERNATE PRODUCT Let A and B be n × n matrices, and let V n be the n(n−1) 2 tuple consisting of pairs of integers (p, q), p = 2, 3 . . . n, q = 1, . . . , p − 1, listed lexicographically. That is, V n = [(2, 1), (3, 1), (3, 2), (4, 1), (4, 2), (4, 3), . . . , (n, n − 1)] Denote by Vin the ith entry of V n and 1 b b a a f ((p, q), (r, s)) = det pr ps + det pr ps bqr bqs aqr aqs 2
where the dependence of f on A and B is kept implicit for simplicity. The bialternate product of A and B, denoted A B, is a n(n−1) × n(n−1) matrix whose ijth entry is 2 2 given by (A B)ij = f Vin , Vjn Additional results and properties are available in (Saydy et al., 1990; Stephanos, 1900). Appendix B. NUMERICAL VALUES Let Fij be the concatenation of the coefficient matrices, i.e., Aij Bij Fij = (B.1) Cij Dij For the flight condition considered in this paper, the numerical values are then:
F00
F10
F01
−1.0170 × 100 1.1740 × 100 = −2.4045 × 101 0 2.3938 × 10−3 −2.8172 × 101 = 1.3688 × 101 0 1.2008 × 10−4 −3.7252 × 10−6 = 2.7762 × 10−3 0
9.4144 × 10−1 −9.5194 × 10−3 −9.4450 × 10−1 −1.2721 × 100 −9.1252 × 10−1 3.7316 × 10−1 1 0 −4.7356 × 10−3 −5.6373 × 10−3 0 −1 −1.6400 × 10 −2.7633 × 10 6.5520 × 10−1 −1.4142 × 10−3 0 0 −6 6.0565 × 10 1.0318 × 10−6 −3.6929 × 10−6 −9.7736 × 10−7 1.4166 × 10−4 2.4296 × 10−5 0 0
99
−4.4109 × 10−7 2.8030 × 10−4 = −1.4088 × 10−4 0 −1.2494 × 10−1 7.3117 × 10−1 = −3.2277 × 100 0 −1.3059 × 10−8 1.4762 × 10−11 = −3.0173 × 10−7 0
99
2.7216 × 10−8 −8.3410 × 10−9 −9.5814 × 10−6 −2.0562 × 10−6 5.0961 × 10−6 7.6601 × 10−7 0 0 −7.1967 × 10−3 −1.2285 × 10−3 −1 −1 −1.3233 × 10 −1.3019 × 10 −1.0458 × 10−1 3.2321 × 10−2 0 0 −10 −6.4278 × 10 −1.1114 × 10−10 −1.2111 × 10−10 2.3641 × 10−11 −1.4795 × 10−8 −2.5790 × 10−9 0 0