# Robust Stability Clearance of Flight Control Law Based on Global Sensitivity Analysis

Authors: Liuli Ou; Lei Liu; Shuai Dong; Yongji Wang
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

To validate the robust stability of the flight control system of hypersonic flight vehicle, which suffers from a large number of parametrical uncertainties, a new clearance framework based on structural singular value (μ) theory and global uncertainty sensitivity analysis (SA) is proposed. In this framework, SA serves as the preprocess of uncertain model to be analysed to help engineers to determine which uncertainties affect the stability of the closed loop system more slightly. By ignoring these unimportant uncertainties, the calculation of μ can be simplified. Instead of analysing the effect of uncertainties on μ which involves solving optimal problems repeatedly, a simpler stability analysis function which represents the effect of uncertainties on closed loop poles is proposed. Based on this stability analysis function, Sobol’s method, the most widely used global SA method, is extended and applied to the new clearance framework due to its suitability for system with strong nonlinearity and input factors varying in large interval, as well as input factors subjecting to random distributions. In this method, the sensitive indices can be estimated via Monte Carlo simulation conveniently. An example is given to illustrate the efficiency of the proposed method.


---

## Body

## 1. Introduction

In the past decades, clearance of flight control system has been paid great attention by the air force of many countries. Rational clearance improves not only the reliability and safety of flight control system, but also feedback valuable information to designers. Clearance methods split neatly into two types based on the different clearance principle. One is the so-called analytical-model-based method (AMBM) [1] and another is simulation-data-based method (SDBM) [2, 3]. AMBM analyses stable margin, robust stability and some other criteria based on advanced control theories which include (1) structural singular value (μ) analysis [4, 5]; (2) polynomial-based clearance [6]; (3) v-gap analysis; (4) bifurcation and continuation method [7]; and (5) optimization-based clearance [8]. AMBM has been improved and extended to more criteria and more types of air vehicles [9–14]. SDBM demands high fidelity simulation platform and combines experiment design techniques with decision sciences. Generally speaking, AMBM is suitable for early stages of the whole controller designing process, whereas SDBM is suitable for late stages. These two types of methods are complementary.For hypersonic flight vehicle (HFV) whose maximum flight speed may reach to about 20 Mach, the principal criterion to validate is robust stability of its attitude control system because of the existence of uncertain coefficients. The uncertainties are mainly caused by the following three factors. First, the aerodynamic coefficients of HFV, obtained via wind tunnel test and computational fluid mechanics tools, often suffer from much more serious deviations than those of supersonic and subsonic vehicles. Second, it is still unknown how the external environment affects the HFV due to lack of flight experience. Finally, erosion and corrosion may damage the aerodynamic configuration of the vehicle and hence change the aerodynamic coefficients. Many advanced methods [15–17] have the potential to improve the ability and performance of the attitude control system, which may also lead to more complicated controllers and require more rational clearance.Among these model-based clearance methods mentioned above,μ analysis is the favourable method for robust stability criterion due to its nonconservatism. Calculation of μ with pure real uncertainties [18] had always been a challenge until the optimization-based approach [19] was proposed. Hybrid optimization [8], which combines sequential quadratic programming (SQP) with intelligent optimization, can produce an exact estimation of μ. However, the computational burden is still heavy due to iteratively sweeping frequencies. In particular, when the number of uncertainties is big, it will take a long time to get μ values at a series of frequency points.In order to reduce the computational burden ofμ, engineers usually ignore some unimportant uncertain factors according to their experience. How to find those minor uncertain factors of the HFV when we have little experience? The answer is sensitivity analysis (SA). SA methods can be broken down into local SA [20, 21] and global SA [22–25]. Global SA has been applied widely in industrial production, ecological engineering, hydrology systems, and so on [26] due to its two main advantages. First, it is suitable for input factors varying in broad range and system model with strong nonlinearity. Second, it can consider input factors obeying certain random distributions.In this work, we proposed a new clearance framework for robust stability clearance of the control system of HFV based onμ theory and global uncertainty SA (GUSA), which is shown in Figure 1. The GUSA serves as the preprocess of the μ analysis to help engineers to determine those minor uncertainties. Because measuring the direct effect of uncertainties on μ involves computing μ repeatedly, a novel indirect GA index based on our newly defined stability analysis function is adopted. This function represents the effect of uncertainties on dominant closed loop poles. With this new GA index, we can extend the widely used Sobol’s method [23] and apply it to analyse the effect of uncertainties on stability of the flight control system. It is worth noting that the framework is suitable for not only μ analysis but also other robust analysis approaches. Furthermore, GUSA can provide valuable information about the effect of aerodynamic coefficients on stability for subsequent design and validation.Figure 1
The new clearance process.The paper is organized as follows. In Section2, the uncertain model of the attitude control system of HFV is established and μ framework is presented. The Sobol’s method is extended and applied to GUSA in Section 3. Section 4 gives an example to illustrate the efficiency of the proposed method. In Section 5, some conclusions are stated.


## 2. Problem Description

The motion equations of the hypersonic flight vehicle are given as follows [27]:
(1)α˙=CAqSmsin⁡α-CNqSmcos⁡⁡αmvcos⁡⁡β-ωxcos⁡⁡αsin⁡βcos⁡⁡β+ωysin⁡αtan⁡β+ωz-gyvcos⁡⁡ββ˙=CAqSmcos⁡⁡αsin⁡β+CNqSmsin⁡αsin⁡β+CZqSmcos⁡⁡βmv+gz+ωxsin⁡α+ωycos⁡⁡αvγ˙=ωxω˙x=CmxqSmlkJx+ClpqSmlk2ωx2vJx+(Jy-Jz)ωyωzJxω˙y=CmyqSmlkJy+CnrqSmlk2ωy2vJy+(Jz-Jx)ωzωxJyω˙z=CmzqSmlkJz+CmqqSmlk2ωz2vJz+(Jx-Jy)ωxωyJz,
whereα,β,γ
= attack, sideslip, and roll angle;ωz,ωy,ωx = pitch, yaw, and roll rate;CA,CN,CZ = aerodynamic coefficients;Cmx,Cmy,Cmz = aerodynamic moment coefficients;Clp,Cnr,Cmq = damping moment coefficients;Jx,Jy,Jz = inertia moment in x, y, z axes;Sm,lk = wing area and span;q = dynamic pressure;m = mass;v = velocity; andδφ,δψ,δγ = elevator, rudder, and aileron deflection.These aerodynamic coefficients are nonlinear functions of the deflection anglesδφ, δψ, and δγ. Consider
(2)Ci=fi(δφ,δψ,δγ),i=A,N,Z,mx,my,mz,lp,nr,mq
whose nominal values are given in the form of look-up tables. Ci, q, Jx, Jy, and Jz are uncertain coefficients varying in some certain intervals around their nominal values.Linearizing (1) and (2) along the nominal flight trajectory according to the methods introduced in [4, 18], the uncertain model P can be written as
(3)x˙p=Apxp+Bpu,
where xp=[αωzβωyγωx]T is the state vector, u=[δφδψδγ]T is the control input, and Ap∈R6×6 and Bp∈R6×3 are uncertain matrices in the affine form as follows:
(4)Ap=Ap0+∑i=1kσiApiBp=Bp0+∑i=1kσiBpi,
where k is the number of normalized uncertainties after combination and linearization and σi∈[-1,1] are uncertain parameters that represent the uncertainties of these coefficients. The uncertain model (4) is not unique via different modelling methods. Because of different designing and manufacturing environment for different vehicles, uncertainties may obey certain random distribution. Consider
(5)σi~R,
for instance, the normal distribution σi~N(μ,υ2) or the uniform distribution σi~U(a,b).Suppose that a state feedback controllerC, which can stabilize the nominal system of (3), is
(6)x˙c=Acxc+Bc1xp+Bc2ru=Ccxc+Dc1xp+Dc2r,
where xc∈Rd×1 is the state vector of controller, r=[αcβcγc]T is the command attitude angles to be tracked, Ac∈Rd×d, Bc1∈Rd×6, Bc2∈Rd×3, Cc∈R3×d, Dc1∈R3×6, and Dc2∈R3×3.Then, the closed loop system, as shown in Figure2, can be written as
(7)[x˙px˙c]=(A0+ΔA)[xpxc]+(B0+ΔB)r,
where
(8)A0=[Ap0+Bp0Dc1Bp0CcBc1Ac]B0=[Bp0Dc2Bc2]ΔA=∑i=1kσiAi=∑i=1k(σi[Api+BpiDc1BpiCc00])ΔB=∑i=1kσiBi=∑i=1k(σi[BpiDc20]).Figure 2
Closed loop system of the attitude control of HFV.It is clear that the stability of the system (7), suffering from affine interval uncertainties σi, is determined by the system matrix A0+ΔA. In μ framework, the linear fraction representation (LFR) of (7), denoted as F(M,Δ), should be modelled first according to the approach presented in [4, 18]. Consider
(9)[x˙px˙c]=[Ap0+Bp0Dc1Bp0CcBc1Ac][xpxc]+[U1U2⋯Uk]wz=[V1TV2T⋯VkT]T[xpxc]w=Δz=diag⁡(σiIri)z}M,
where M represents known part of the system with input w and output z, Δ=diag⁡{σiIri} represents uncertainties of diagonal structure, rank⁡(·) is the rank of a matrix, and Ui∈R(n+d)×ri and Vi∈Rri×(n+d) are obtained by full rank decomposition. Consider
(10)UiVi=[Api+BpiDc1BpiCc00]ri=rank⁡([Api+BpiDc1BpiCc00]).At a specific frequency pointω, structural singular value of F(M,Δ) (9) is defined as
(11)μΔ(M(jω))=1min⁡⁡{∥Δ∥∞∣det⁡(I-M(jω)Δ)=0}>0,
where ∥·∥∞ is the infinite norm of a matrix. Define μ*=max⁡ω⁡μΔ(M(jω)). Since ∥Δ∥∞≤1, the sufficient condition of robust stability is max⁡⁡{∥Δ∥∞}<1/μ*; that is, μ*<1. Furthermore, the system stays stable even if the uncertain intervals extend to σi∈[-1/μ*,1/μ*]. For simplicity, we denote μΔ(M(jω)) as μ hereinafter.According to definition (11), if the problem
(12)min⁡⁡γ-s.t.0(1)∥Δ∥∞≤γ-s.t.0(2)det⁡(I-M(jω)Δ)=0
has an optimal solution γ-*, then the estimation of μ is 1/γ-*. By combining global intelligent algorithm (such as particle swarm optimization (PSO) and differential evolution (DE) [28]) with local SQP, this problem can be iteratively solved [8]. In order to get the estimation of μ*=max⁡ω⁡μ, we need to solve the optimal problem (12) at a series of frequency points.When the number of uncertainties is large, the dimension of the searching space is also large. In order to cover the whole space uniformly, the number of particles of PSO (or analogy of other intelligent algorithms) should exponentially increase with respect to the number of uncertainties. Then, the computational burden will increase exponentially. Therefore, our aim is to find thoseσis which affect the stability dramatically and ignore remaining unimportant σis when calculating μ*. If we measure the effect of uncertainties on μ* directly, we need to calculate μ* repeatedly. So we introduce an indirect stability analysis function.As we know, a system is stable if and only if all of its closed loop poles lie in the left half of the complex plan. Therefore, we define the nonlinear stability analysis function as follows:(13)f(σi)=max⁡j{Re[λj(A0+ΔA)]},
where λj(·) is the jth eigenvalue of a matrix and Re[·] is the real component of a complex. The effect of σi on the stability of the closed loop system can be represented by (13) and SA can be invoked. It is clear that the stability analysis function (13) is also suitable for other robust stability clearance method.All local SA methods need to calculate the differential or derivative; that is, they are based on the analytical model. However, the nonlinear function (13) does not satisfy this requirement. So we will use the global Sobol’s method to deal with this problem in the upcoming sections.


## 3. Global Uncertainty Sensitivity Analysis

### 3.1. Sobol’s Method

Consider the nonlinear model in the form(14)f:Rk⟶Rx⟶y=f(x),
where x=[x1,x2,…,xk] are k independent input factors and y is the scalar output. If x is defined over Ω={x∣0≤xi≤1,i=1,2,…,k}, then the function can be decomposed as the following form [23]:
(15)f=f0+∑ifi+∑i∑j>ifij+∑i1∑i2>i1∑i3>i2fi1i2i3+⋯+f12⋯k,
where fi=fi(xi), fij=fij(xi,xj), and so on. The uniqueness of (15) is guaranteed by
(16)∫01fi1,i2,…,is(xi1,xi2,…,xis)dxiw=0,
where 1≤i1<i2<⋯<is≤k and iw∈{i1,i2,…,is}.Denotexi1i2⋯is=[xi1,xi2,…,xis], 1≤i1<i2<⋯<is≤k, and x~i1i2⋯is as the complement of xi1i2⋯is; then we have
(17)x~i=[x1,…,xi-1,xi+1,…,xk]x~ij=[x1,…,xi-1,xi+1,…,xj-1,xj+1,…,xk]
and so on. Define the corresponding domain as Ωi1i2⋯is and Ω~i1i2⋯is; then we have
(18)Ωi1i2⋯is={xi1i2⋯is∣0≤xl≤1,l=i1,i2,…,is}Ω~i1i2⋯is={x~i1i2⋯is∣0≤xl≤1,l≠i1,i2,…,is}.These functions in (15) are obtained by
(19)f0=∫Ωf(x)dxfi=∫Ω~if(x)dx~i-f0fij=∫Ω~ijf(x)dx~ij-fi-fj-f0
and similarly for higher orders, where ∫Ωi1⋯irf(x)dxi1⋯ir=∫01∫01⋯∫01f(x)dxi1dxi2⋯dxir.The variance ofy is
(20)V=V(y)=∫Ωf2(x)dx-f02.
And the partial variance is defined as
(21)Vi1i2,…,is=∫Ωi1i2⋯isfi1i2⋯is2dxi1i2⋯is.
Then the decomposition (15) leads to the following decomposition of the variance of y(22)V=∑i=1kVi+∑1≤i<j≤kVij+⋯+V12⋯k.The normalized sensitivity index is defined as(23)Si=ViVSij=VijVSi1i2,…,is=Vi1i2,…,isV,
where Si are the so-called first-order sensitivity indices and Sij and Si1i2,…,is are the second-order and the sth-order indices, respectively. Si1i2,…,is indicates sensitivity of y to the interaction among xi1,xi2,…,xis. So the total sensitivity indices STi for the input xi is defined as the sum of all indices relating to it. For instance, if k=3, we have ST1=S1+S12+S13+S123. It is easy to find that
(24)STi=1-V~iV
and 0≤STi≤1, where V~i is the (k-1)th indices of the remaining k-1 factors except the ith one.Because these global sensitivity indices are based on variance, enough samples are needed. The best and most widely used method is Monte Carlo (MC) simulations.SupposingX∈RN×k is the matrix of N groups of k input factors
(25)X=(xij)N×k=[X1,X2,…,Xk]=[x1T,x2T,…,xNT]T,
where Xi∈RN×1 and xi∈R1×k. Y∈RN×1 is the corresponding output
(26)Y=f(X)=[f(x1),f(x2),…,f(xN)]T=[y1,y2,…,yN]T,
and P and Q are two samples of X generated by MC. Then we define that
(27)PQ(i)=[P1⋯Pi-1QiPi+1⋯Pk]
is the combined matrix with column i from Q and all other k-1 columns from P, where Pi and Qi are samples of Xi. Similarly, we define
(28)QP(i)=[Q1⋯Qi-1PiQi+1⋯Qk].The estimation below can be obtained [24, 25] as follows:
(29)f^0=12N∑j=1N[f(pj)+f(qj)]V^=12N∑j=1N[f2(pj)+f2(qj)]-f^02V^i=1N∑j=1Nf(qj)f((PQ(i))j)-f^02V^~i=1N∑j=1Nf(qj)f((QP(i))j)-f^02S^i=V^iV^,S^Ti=1-V^~iV^,
where pj, qj, (PQ(i))j, and (QP(i))j are, respectively, the jth row of P, Q, PQ(i), and QP(i).


### 3.2. Uncertainties Sensitivity Analysis Based on Sobol’s Method

It is worth noting that the input factors of the nonlinear model considered in Sobol’s method are limited to0≤xi≤1. However, the normalized uncertainties of the flight control system is σi∈[-1,1]. One manner is to renormalize the uncertainties to the interval [0,1]. Here we will introduce an easier method.Considering the nonlinear model (14), if its domain of definition is
(30)Π={x∣-0.5≤xi≤0.5,i=1,2,…,k},
we can get the below theorem.Theorem 1.
The nonlinear model (14) defined over Π can be decomposed uniquely to (15) with its items (34) if
(31)∫-1/21/2fi1,i2,…,is(xi1,xi2,…,xis)dxiw=0
hold for all 1≤i1<i2<⋯<is≤k and iw∈{i1,i2,…,is}.Proof.
From (31), we can obtain that the members in (15) are orthogonal. That is, for any
(32)1≤i1<i2<⋯<is≤k1≤j1<j2<⋯<jr≤k,
if (i1,i2,…,is)≠(j1,j2,…,jr), there exists
(33)∫Πfi1i2⋯isfj1j2⋯jrdx=0.
Then, we have
(34)∫Πf(x)dx=f0∫Π~if(x)dx~i=fi(xi)+f0∫Π~ijf(x)dx~ij=fij(xi,xj)+fi(xi)+fj(xj)+f0
and so on for high order terms. Obviously, (34) is similar to (19).Remark. In fact, Sobol’s method can be extend to any nonlinear function defined over
(35)Γ={x∣x_≤xi≤x-,i=1,2,…,k},
where x--x_=1.The uncertain matrixΔA of the system matrix in (7) can be rewritten as
(36)ΔA=[∑iσi2(2Api+2BpiDc1)∑iσi22BpiCc00]=∑i2σ-i[Api+BpiDc1BpiCc00],
where σ-i=σi/2∈[-0.5,0.5] are the normalized uncertainties. It is clear that the uncertain space of σ-i is Π defined in (30).Substituting (7) and (36) into the nonlinear function (13), we can get
(37)f(σ-)=max⁡j{[λj(A0+∑iσ-i[2Api+2BpiDc12BpiCc00])]Re[λj(∑iσ-i[2Api+2BpiDc12BpiCc00]A00000000000000000000+∑iσ-i[2Api+2BpiDc12BpiCc00])]},
where σ-=(σ-1,σ-2,…,σ-k), σ-∈Π. The Sobol’s method can be applied to analyse the effect of uncertainty terms σ-i on the real part of the closed loop poles. The sample matrices of σ- can be generated randomly according to the distribution of σ.Generally, the whole clearance process can be summarized into the following.Step  1:
establish the uncertain model of the attitude control system and normalize the uncertainties to[-0.5,0.5].Step  2:
construct stability analysis function as (37) according to the specific uncertain model.Step  3:
generate the sample matricesP and Q of N groups of uncertainties σ-.Step  4:
calculate output samples of the stability analysis functions (37) according to MC simulations with input samples P and Q.Step  5:
calculate the SA index according to (29) and ignore the minor uncertainties with relatively small SA index.Step  6:
calculate the singular structure value according to (12).


## 3.1. Sobol’s Method

Consider the nonlinear model in the form(14)f:Rk⟶Rx⟶y=f(x),
where x=[x1,x2,…,xk] are k independent input factors and y is the scalar output. If x is defined over Ω={x∣0≤xi≤1,i=1,2,…,k}, then the function can be decomposed as the following form [23]:
(15)f=f0+∑ifi+∑i∑j>ifij+∑i1∑i2>i1∑i3>i2fi1i2i3+⋯+f12⋯k,
where fi=fi(xi), fij=fij(xi,xj), and so on. The uniqueness of (15) is guaranteed by
(16)∫01fi1,i2,…,is(xi1,xi2,…,xis)dxiw=0,
where 1≤i1<i2<⋯<is≤k and iw∈{i1,i2,…,is}.Denotexi1i2⋯is=[xi1,xi2,…,xis], 1≤i1<i2<⋯<is≤k, and x~i1i2⋯is as the complement of xi1i2⋯is; then we have
(17)x~i=[x1,…,xi-1,xi+1,…,xk]x~ij=[x1,…,xi-1,xi+1,…,xj-1,xj+1,…,xk]
and so on. Define the corresponding domain as Ωi1i2⋯is and Ω~i1i2⋯is; then we have
(18)Ωi1i2⋯is={xi1i2⋯is∣0≤xl≤1,l=i1,i2,…,is}Ω~i1i2⋯is={x~i1i2⋯is∣0≤xl≤1,l≠i1,i2,…,is}.These functions in (15) are obtained by
(19)f0=∫Ωf(x)dxfi=∫Ω~if(x)dx~i-f0fij=∫Ω~ijf(x)dx~ij-fi-fj-f0
and similarly for higher orders, where ∫Ωi1⋯irf(x)dxi1⋯ir=∫01∫01⋯∫01f(x)dxi1dxi2⋯dxir.The variance ofy is
(20)V=V(y)=∫Ωf2(x)dx-f02.
And the partial variance is defined as
(21)Vi1i2,…,is=∫Ωi1i2⋯isfi1i2⋯is2dxi1i2⋯is.
Then the decomposition (15) leads to the following decomposition of the variance of y(22)V=∑i=1kVi+∑1≤i<j≤kVij+⋯+V12⋯k.The normalized sensitivity index is defined as(23)Si=ViVSij=VijVSi1i2,…,is=Vi1i2,…,isV,
where Si are the so-called first-order sensitivity indices and Sij and Si1i2,…,is are the second-order and the sth-order indices, respectively. Si1i2,…,is indicates sensitivity of y to the interaction among xi1,xi2,…,xis. So the total sensitivity indices STi for the input xi is defined as the sum of all indices relating to it. For instance, if k=3, we have ST1=S1+S12+S13+S123. It is easy to find that
(24)STi=1-V~iV
and 0≤STi≤1, where V~i is the (k-1)th indices of the remaining k-1 factors except the ith one.Because these global sensitivity indices are based on variance, enough samples are needed. The best and most widely used method is Monte Carlo (MC) simulations.SupposingX∈RN×k is the matrix of N groups of k input factors
(25)X=(xij)N×k=[X1,X2,…,Xk]=[x1T,x2T,…,xNT]T,
where Xi∈RN×1 and xi∈R1×k. Y∈RN×1 is the corresponding output
(26)Y=f(X)=[f(x1),f(x2),…,f(xN)]T=[y1,y2,…,yN]T,
and P and Q are two samples of X generated by MC. Then we define that
(27)PQ(i)=[P1⋯Pi-1QiPi+1⋯Pk]
is the combined matrix with column i from Q and all other k-1 columns from P, where Pi and Qi are samples of Xi. Similarly, we define
(28)QP(i)=[Q1⋯Qi-1PiQi+1⋯Qk].The estimation below can be obtained [24, 25] as follows:
(29)f^0=12N∑j=1N[f(pj)+f(qj)]V^=12N∑j=1N[f2(pj)+f2(qj)]-f^02V^i=1N∑j=1Nf(qj)f((PQ(i))j)-f^02V^~i=1N∑j=1Nf(qj)f((QP(i))j)-f^02S^i=V^iV^,S^Ti=1-V^~iV^,
where pj, qj, (PQ(i))j, and (QP(i))j are, respectively, the jth row of P, Q, PQ(i), and QP(i).


## 3.2. Uncertainties Sensitivity Analysis Based on Sobol’s Method

It is worth noting that the input factors of the nonlinear model considered in Sobol’s method are limited to0≤xi≤1. However, the normalized uncertainties of the flight control system is σi∈[-1,1]. One manner is to renormalize the uncertainties to the interval [0,1]. Here we will introduce an easier method.Considering the nonlinear model (14), if its domain of definition is
(30)Π={x∣-0.5≤xi≤0.5,i=1,2,…,k},
we can get the below theorem.Theorem 1.
The nonlinear model (14) defined over Π can be decomposed uniquely to (15) with its items (34) if
(31)∫-1/21/2fi1,i2,…,is(xi1,xi2,…,xis)dxiw=0
hold for all 1≤i1<i2<⋯<is≤k and iw∈{i1,i2,…,is}.Proof.
From (31), we can obtain that the members in (15) are orthogonal. That is, for any
(32)1≤i1<i2<⋯<is≤k1≤j1<j2<⋯<jr≤k,
if (i1,i2,…,is)≠(j1,j2,…,jr), there exists
(33)∫Πfi1i2⋯isfj1j2⋯jrdx=0.
Then, we have
(34)∫Πf(x)dx=f0∫Π~if(x)dx~i=fi(xi)+f0∫Π~ijf(x)dx~ij=fij(xi,xj)+fi(xi)+fj(xj)+f0
and so on for high order terms. Obviously, (34) is similar to (19).Remark. In fact, Sobol’s method can be extend to any nonlinear function defined over
(35)Γ={x∣x_≤xi≤x-,i=1,2,…,k},
where x--x_=1.The uncertain matrixΔA of the system matrix in (7) can be rewritten as
(36)ΔA=[∑iσi2(2Api+2BpiDc1)∑iσi22BpiCc00]=∑i2σ-i[Api+BpiDc1BpiCc00],
where σ-i=σi/2∈[-0.5,0.5] are the normalized uncertainties. It is clear that the uncertain space of σ-i is Π defined in (30).Substituting (7) and (36) into the nonlinear function (13), we can get
(37)f(σ-)=max⁡j{[λj(A0+∑iσ-i[2Api+2BpiDc12BpiCc00])]Re[λj(∑iσ-i[2Api+2BpiDc12BpiCc00]A00000000000000000000+∑iσ-i[2Api+2BpiDc12BpiCc00])]},
where σ-=(σ-1,σ-2,…,σ-k), σ-∈Π. The Sobol’s method can be applied to analyse the effect of uncertainty terms σ-i on the real part of the closed loop poles. The sample matrices of σ- can be generated randomly according to the distribution of σ.Generally, the whole clearance process can be summarized into the following.Step  1:
establish the uncertain model of the attitude control system and normalize the uncertainties to[-0.5,0.5].Step  2:
construct stability analysis function as (37) according to the specific uncertain model.Step  3:
generate the sample matricesP and Q of N groups of uncertainties σ-.Step  4:
calculate output samples of the stability analysis functions (37) according to MC simulations with input samples P and Q.Step  5:
calculate the SA index according to (29) and ignore the minor uncertainties with relatively small SA index.Step  6:
calculate the singular structure value according to (12).


## 4. Examples

In this section, a numerical example is given to reveal the effectiveness of the proposed method. Consider the nominal flight trajectory of a HFV depicted as Figure3, from which 30 operation points for robust stability clearance have been selected. At first, we will illustrate the application of the new clearance process with one single point. Then, we will analyse all these operation points together.Figure 3
Nominal trajectory of the HFV.At the 9th point (Ma=11.8670, h=41928), the nominal model of the attitude control equation in the form of (3) is
(38)Ap0=[-0.013810000-7.64330000000-0.00820.972500.2300-5.338500000000100-33.8537000]Bp0=[-0.000500-2.5109000-0.0001-0.00010-0.2711-1.24410000-1.1006-52.7017].
Because the pitch motion has no coupling with yaw and roll motion, the controller can be designed separately. For convenience, a low order controller that can stabilize the nominal system has been designed as follows:
(39)u=[L100L2]r-[K100K2]xp,
where
(40)L1=-99.9845L2=[-3.36640.6189-0.2088-9.9808]K1=[-98.1456-7.9508]K2=[-1.6586-3.19710.61890.0781-1.74282.4028-9.9808-0.6311].Then the system matrix of the closed loop system is(41)A=[Aclo100Aclo2]+ΔA,
where
(42)Aclo1=[-0.05870.9964-254.0734-19.9634]Aclo2=[-0.00850.9726-0.00120.2329-7.95622.1227-12.2491-0.76400001-127.5277123.1129-525.3240-33.1741]ΔA=[000000Δa21Δa22000000000000Δa43Δa44Δa45Δa4600000000Δa63Δa64Δa65Δa66]Δa21=-1.5287σ1+49.286σ4Δa22=3.9927σ4Δa43=-2.6692σ2+0.0899σ5+0.4336σ6Δa44=0.1733σ5-0.5979σ6Δa45=-0.0336σ5+2.4834σ6Δa46=-0.0042σ5+0.1570σ6Δa63=-6.7707σ3+1.4604σ7+64.2940σ8Δa64=2.8150σ7-88.6422σ8Δa65=-0.5449σ7+368.2036σ8Δa66=-0.0688σ7+23.2820σ8σi~U(-1,1),i=1,2,…,8.Calculate the global sensitivity indices of the function (37) with respect to σ-i=σi/2~U(-0.5,0.5) according to the steps presented in Section 3. In order to illustrate the consistency and repeatability of the proposed method, 10 MC simulations have been conducted. The obtained global sensitivity indices STi are depicted in Figures 4 and 5. It is evident that (1) the results of 10 simulations are nearly consistent; (2)  ST1, ST3, and ST4 are extremely close to 0 and obviously smaller than other indices; and (3) when the sample number N increases from 20000 to 50000, the results converge accordingly. Varying σi in the interval [-1,1] singly and fixing σj=0, j≠i, the distribution of the dominant closed loop poles is shown as Figure 6. This can directly validate the SA results.Figure 4
The global uncertainty sensitivity indices (N=20000).Figure 5
The global uncertainty sensitivity indices (N=50000).Figure 6
Dominant poles distributions for systems with one uncertainty.SinceSTi≈0, i=1,3,4, we can ignore σi, i=1,3,4 in the analysis of robust criterion. μ is calculated via hybrid optimization method at 136 frequency points ω=[0.01,0.02,…,0.09,0.1,0.2,…,3,4,5,…,100] in the MATLAB environment. Contrastive results are shown in Figures 7 and 8. When only considering the remaining 5 major uncertainties (σi,i=2,5,6,7,8), the μ curve is very close to the situation of considering all 8 uncertainties. The computational time, however, decreases from 3.3387 s to 2.2903 s when the number of uncertainties decreases from 8 to 5. When taking into account all the operation points together, the proposed clearance frame can save a large amount of time. Because μ*≈1.1>1, the controller cannot stabilize the system in the whole uncertain space. However, we can get the valuable information that the approximate allowable uncertain space is |σi|<1/1.1≈0.9 and feed it back to designers.Figure 7
Structural singular value of the uncertain systems.Figure 8
Structural singular value in low frequency interval of the uncertain systems.The global sensitivity indicesSTi of all operation points are plotted in Figure 9. Note that STi≈0, i=1,3,4, holds for all operation points. That is, because the flight trajectory in Figure 3 is relatively smooth and the uncertain model varies slowly. Therefore, we just need to analyse a few operation points when the trajectory is smooth. A comparison of computational time of μ is also depicted in Figure 10. With the proposed clearance framework, we can save more than 30 s totally while the GUSA just takes 3.6645 s.Figure 9
GSA of all operation points.Figure 10
Computational time ofμ for all operation points.In summary, GUSA provides a quantitative index for engineers to measure the effect of uncertainties on dominant closed loop poles. It is worth noting that the GUSA index is relative index. It means that we can order the effect of all uncertainties based on the GUSA index, but we cannot find a threshold value to ignore unimportant uncertainties. There is a tradeoff between precision and computational time when ignoring unimportant uncertainties. For instance, if we need to reduce more computational time, we can ignore one more uncertaintyσ6.


## 5. Conclusion

In this paper, a newμ framework based on GUSA is proposed. GUSA serves as the preprocessing step of μ analysis and obviously release the computational burden. Based on our proposed stability analysis function, Sobol’s method is extended and applied to determine the major uncertainties because of its maturity and flexibility. An example is given to illustrate the efficiency of the proposed method. In the potential applications, this work provides an easy-to-follow procedure for engineers to digest and expand.



---

*Source: 153602-2014-10-13.xml*