- Email: [email protected]
NONLINEAR RE-ENTRY FLIGHT CONTROL SYSTEM ROBUSTNESS EVALUATION USING INTERVAL ANALYSIS
NONLINEAR RE-ENTRY FLIGHT CONTROL SYSTEM ROBUSTNESS EVALUATION USING INTERVAL ANALYSIS Sinar Juliana, ∗,1 Q.P. Chu, ∗,2 J.A. Mulder ∗,3 ∗
Faculty of Aerospace Engineering, Delft University of Technology, PO Box 5058, 2600 GB Delft, The Netherlands
Abstract: A stability robustness analysis is performed on a re-entry vehicle flight control system. The mathematical model of this system is nonlinear and may contain uncertain parameters, for which a set of models has to be examined instead of a singleton. The robustness evaluation is performed using interval analysis, which allows one to certify uncertainty regions affecting the closed-loop system. The celebrated Lyapunov function is taken as the stability criterion for the closed-loop system with uncertainties, since the existence of this function provides sufficient stability condition. The result shows the feasibility of the method for c nonlinear system analysis. Copyright 2007 IFAC Keywords: closed-loop system, interval, Lyapunov, robustness, stability
1. INTRODUCTION A re-entry flight vehicle can be modeled mathematically as a highly parameter dependent nonlinear system, for the purpose of controller design and analysis. Because of the nature of the model, nonlinear control techniques are suitable for this type of system to achieve desired stability and performance. In practice, the model may contain uncertainties which may be caused by modeling errors, disturbances from the environment, unknown parameter variation, and measurement errors. As a consequence, the controller may not be able to perform adequately and eventually may cause damage to the system. Uncertainties in re-entry flight models may arise from various sources. For instance, the knowledge about re-entry aerodynamics is still not mature; hence no accurate models are available. Nonethe1 2 3
Researcher Associate Professor Professor, Head of Control and Simulation Division
less, one has to make sure that the designed controllers are robust. To this objective, this paper presents a method to evaluate the stability robustness of re-entry flight control systems in the presence of parametric uncertainties. The Lyapunov stability is chosen as the criterion for the evaluation, for the existence of a Lyapunov function provides sufficient condition for stability in the system. Moreover, the criterion is applicable to both linear and nonlinear input-output systems. A closed-loop system with parametric uncertainties can be considered as a set of model whose numerical values may vary within some unknown but bounded intervals. Therefore, robustness analysis must be applied to this set, instead of to a single model. In this paper, an interval analysis algorithm based on set inversion is proposed to perform this evaluation. Using this approach, all elements in the set can be analysed simultaneously. In (Verde and Corraro, 2002), a similar algorithm was proposed for clearance of flight control law; however, the scope was limited to
linear systems only. In this paper, the algorithm is applied to analyse nonlinear control systems. For completeness, this algorithm is described in Section 3, based on an extensive discussion in (Jaulin et al., 2001; Jaulin and Walter, 1993). To illustrate the method, a nonlinear re-entry flight control system is taken as a case study (Juliana et al., 2004). The controller is designed using the nonlinear dynamic inversion (NDI) technique, which is sensitive to possible model errors (Slotine and Li, 1991). The error can propagate through the time, which eventually can cause system instability. The system under consideration is assumed to contain uncertain parameters in its aerodynamic models. A modification of the NDI technique allows the description of the parametric uncertainties to be included in the model (Juliana et al., 2005). The Lyapunov function evaluation in interval analysis enables one to perform the analysis on a continuous system, omitting the drawback of the classical grid point approach. The result of the robustness analysis, described in Section 4, shows the feasibility of the method to be applied to general nonlinear systems with uncertain parameters. The method is guaranteed to give effective result, for no critical region can be missed in the analysis. The result is validated by performing worst-case simulations of the closed-loop system.
2. LYAPUNOV STABILITY For a general nonlinear autonomous system model, one can construct a positive definite function, which has continuous partial derivatives in the neighbourhood of the equilibria, and whose time derivative along any state trajectory of the system is negative semi-definite. Consider a general nonlinear system model, x˙ = f (x),
(1)
where x is the state vector of the system, and f is a nonlinear continuous function. The system is autonomous, as the function f does not depend explicitly on time. For this system, construct a positive definite function V (x) : Rn → R, V (x) > 0, (2)
The Lyapunov theorem for local stability is defined by Equation (2) and (3) (La Salle and Lefschetz, 1961). When the time derivative of a positive-definite Lyapunov function is negative semi-definite, the system is said to be locally stable. When this derivative is negative definite, the system is said to be locally asymptotically stable. For global stability, one has to consider the whole state-space instead of the neighbourhood of the equilibria. Moreover, the Lyapunov function V (x) has to be radially unbounded. The existence of a Lyapunov function for the system in Equation (1) provides sufficient condition of stability in the system. However, there is generally no effective method to construct a Lyapunov function. Therefore, trial-and-error methods are usually employed to find it. Fortunately, for mechanical systems, the energy dissipation principle can be used to form a suitable choice of the Lyapunov function, as will be shown by a case study.
3. INTERVAL ANALYSIS For the computation with intervals of numbers, the rules are defined differently from the operations for single numbers. The basic rules for arithmetic operations, set operations and inclusion tests can be found in various textbooks and papers on interval analysis such as (Hansen and Walster, 1992; Jaulin et al., 2001). These operations are the basics to build the algorithm Set Inverter via Interval Analysis (SIVIA), which is used to perform the stability robustness evaluation in this paper. This section is dedicated to describe the algorithm.
3.1 Subpavings and Set Inversion The notation of subpavings is defined in the following (Jaulin et al., 2001). Take a set of interest X in Rn . This set is covered with sets of non-overlapping boxes of Rn , which are called the subpavings. The set X is also included between its inner and outer approximations, i.e. X ⊂ X ⊂ X,
(4)
for which the knowledge of the inner and outer approximations provides some valuable information about X.
which has continuous partial derivatives in the neighbourhood of the equilibria. If its time derivative along any state trajectory of the system in Equation (1) is negative semi-definite, i.e.
Set inversion is the basic building block in the algorithm presented in the next section. Let f be a possibly nonlinear function from Rn to Rm , and let Y be a subset of Rm . Set inversion is the computation of the reciprocal image
V˙ (x) ≤ 0,
X = {x ∈ Rn |f (x) ∈ Y} = f −1 (Y),
(3)
then V (x) is called the Lyapunov function for the system in Equation (1).
m
(5)
of a regular subpaving Y ⊂ R by a function f : Rn → Rm . The algorithm serves as a method
to compute two subpavings X and X of Rn , such that Equation (4) is satisfied.
3.2 Set Inversion via Interval Analysis The regular subpavings, such as in Equation (4) for X, for any function f with a convergent inclusion function [f ], can be obtained with the set inversion algorithm. This algorithm is illustrated in Figure 1.
of the parameter space of p, for which the derivative of the Lyapunov function is negative semidefinite. In other words, the stability domain can be defined as Sp ≡ p ∈ Rnp V˙ (p) ≤ 0 = V˙ −1 (]−∞, 0]) . (6) The characterization of Sp can, therefore, be cast into the framework of set inversion and performed by the SIVIA algorithm.
4. CONTROL SYSTEM APPLICATION The algorithm is applied to analyse the stability robustness of a nonlinear re-entry flight control system. As a case study, a re-entry capsule equipped with NDI controller is taken from (Juliana et al., 2004). The block diagram of the closed-loop system is shown in Figure 2. While the inner loop feedback control is designed using NDI, in the outer feedback loop, the Proportional, Integral, and Derivative (PID) technique is applied. The vehicle is spun around its body axis to provide gyroscopic rigidity. The controller is designed to keep the spinning rate p constant, and to keep angle of attack α and sideslip angle β zero. Fig. 1. Illustration of set inversion via interval analysis As shown in Figure 1, four cases may be encountered (denoted by Roman numerals): (1) Case I: [f ] ([x]) has a non-empty intersection with Y, but is not entirely in Y. The set [x] then may contain part of the solution set; it is said to be undetermined. If the magnitude of [x] is larger than a predefined small number , the algorithm bisects the interval. The test is further performed recursively on the newly generated intervals. (2) Case II: [f ] ([x]) has an empty intersection with Y. This means that [x] does not belong to the solution set. (3) Case III: [f ] ([x]) is entirely in Y. Therefore, [x] belongs to the solution, and stored in X and X. (4) Case IV: the interval [x] is undetermined, while its width is smaller than . As the interval is small enough, it can be stored in the outer approximation X of X. The stability robustness evaluation can be cast into a set inversion problem. As shown by Equation (3), the derivative of the positive-definite Lyapunov function has to be negative semi definite for the system to be locally stable. When a system contains parametric uncertainties, denoted as p, the Lyapunov function may also contain these uncertainties. The robustness evaluation amounts to find the stability domain Sp , which is a subset
Fig. 2. Feedback control system with NDI and PID
4.1 Closed-loop system with uncertainties In the following, the model of the nonlinear reentry flight control system in the presence of parametric uncertainties is described. The model of the vehicle rotational motion can be written as a set of first-order nonlinear differential equations. The first-order derivatives of the outputs y = T p α β are written in Equation (7) through (9), where q and r are pitch rate and roll rate, Ixx is the moment of inertia in x-axis, MCx is the control moment about x-axis, and MAx is the aerodynamic moment about x-axis. The expressions for the second order derivative of α and β are given in Equation (10) and (11). Equation (7), (10) and (11) are used to perform the dynamic inversion. In the model construction, the magnitude of the aerodynamic moments is assumed to be not fully known, with the uncertainties assumed to be bounded. In a linear fractional representation form, the true aerodynamic moments can be expressed as functions of predicted aerodynamic moments and the uncertainties, as shown
MCx + MAx , Ixx α˙ = −p cos α tan β + q − r sin α tan β, β˙ = p sin α − r cos α, p˙ =
(7) (8) (9)
α ¨ = α˙ (p sin α tan β − r cos α tan β) − β˙ p cos α sec2 β + r sin α sec2 β − cos α tan β (MCx + MAx ) /Ixx (MCy + MAy ) + (Izz − Ixx ) · pr sin α tan β (MCz + MAz ) + (Iyy − Ixx ) · pq sin α tan β + − , (10) Iyy Izz cos α (MCz + MAz ) + (Iyy − Ixx ) · pq cos α sin α (MCx + MAx ) − . (11) β¨ = α˙ (p cos α + r sin α) + Ixx Izz ˜ A, (12) MA = I3×3 + ∆MA · M where I3×3 is a 3-by-3 identity matrix, and ⎡ ⎤ δMAx 0 0 T ˜A = M ˜ Ax M ˜ Ay M ˜ Az T . ∆MA = ⎣ 0 δMAy 0 ⎦ ; MA = MAx MAy MAz ; M 0 0 δMAz Fig. 3. Mathematical model of the re-entry vehicle with aerodynamic uncertainties in Equation (12). The terms δMAx , δMAy , and δMAz represent the lump of all unknowns in the model of aerodynamic moments. Their bounds show the relative values of the uncertainties with respect to the predicted aerodynamic moments. For instance, ˜ Ax 1 + δM MAx = M . (13) Ax
bounds are the first guess of the uncertainties, which can be shrinked or enlarged based on the result of the analysis. Table 1. Uncertain parameters
In Equation (13), the true moment MAx is equal ˜ Ax only if δM = 0. to the predicted moment M Ax Otherwise, the magnitude of MAx will be larger ˜ Ax , depending on δM . or smaller than M Ax As a result of imperfect feedback linearization, the closed-loop system model still contains the uncertain nonlinear terms given by Equation (12). This model can be expressed as ⎡ ⎤ ˜ Ax /Ixx δMAx M ⎡ ⎤ ⎡ ⎤ p˙ ν1 ⎢ ⎥ ˜ ⎥ ⎣α ¨ ⎦ = ⎣ν2 ⎦+ A1 · ⎢ ⎢δMAy MAy /Iyy ⎥ , ⎣ ⎦ ν3 β¨ ˜ /I δ M MAz
with A1
Az
(14)
zz
⎤ ⎡ cos β 0 0 1 ⎣ − cos α sin β cos β − sin α sin β ⎦, = cos β sin α cos β 0 − cos α cos β
and ν1 = (KP 1 + KI1 /s) (pref − p) , ν2 = (KP 2 + KD2 s) (αref − α) , ν3 = (KP 3 + KD3 s) (βref − β) , where the terms KP , KI , and KD are the proportional, integral, and derivative gains, respectively. The stability robustness evaluation is performed to the closed-loop system model in Equation (14) in the presence of the aerodynamic uncertainties, whose values are given in Table 1. The 10%
parameter
bound
δMAx δMAy
[−0.05, 0.05] [−0.05, 0.05]
δMAz
[−0.05, 0.05]
4.2 Stability robustness analysis The form of the closed-loop system may be simplified by defining a new state vector, T ∂q q˙ = = p α˙ β˙ , (15) ∂t and positive definite matrices KI and KII , whose non-zero elements are the PID controller gains, ⎡ ⎤ KP 1 0 0 KI = ⎣ 0 KD2 0 ⎦ (16) 0 0 KD3 ⎡
⎤ KI1 0 0 KII = ⎣ 0 KP 2 0 ⎦ . 0 0 KP 3
and
(17)
Substituting Equation (15) through (17) into (14), and removing the static term from the model, the closed-loop system can be rewritten as ˜ A, q¨ = −KI q˙ − KII q + A2 In−1 ∆MA M with
(18)
⎤ 1 0 0 A2 = ⎣ − cos α tan β 1 − sin α tan β ⎦ ; and sin α 0 − cos α ⎡
⎡
In−1
⎤ 1/Ixx 0 0 = ⎣ 0 1/Iyy 0 ⎦ . 0 0 1/Izz
Table 2. Parameter variation parameter
A Lyapunov function can be constructed as the sum of the artificial kinetic and potential energy functions of the system. The artificial kinetic energy for the system in Equation (18) can be written as a quadratic function of the body rate, 1 T q˙ q, ˙ (19) 2 and the artificial potential energy of the system can be written as 1 T q KII q . Ep = (20) 2 The Lyapunov function is, therefore, 1 T q˙ q˙ + q T KII q , (21) V = 2 which is positive semi-definite. Its first-order derivative can be written as Ek =
V˙ = q˙T q¨ + q˙T KII q,
(22)
which by substituting Equation (18) becomes ˜ A. V˙ = −q˙T KI q˙ + q˙T AIn−1 ∆MA M
(23)
The values of V˙ are functions of the flight conditions, the states, and the uncertain parameters.
interval
α
[−1, 0) ∪ (0, 1] deg
α˙ β β˙
[−5.6, 0) ∪ (0, 5.6] deg/s [−1, 0) ∪ (0, 1] deg
p V ρ
[−5.6, 0) ∪ (0, 5.6] deg/s [−1, 1] deg/s [1041, 1175] m/s [2.4, 2.7] · 10−2 kg /m3
δMAx is regarded as a fix interval during the analysis. This reduces the computational workload. The stable domain SP for the closed-loop system is a 10-dimensional space, whose dimensions are given by the parameters P in Equation (25), in which Equation (24) is satisfied. Based on the analysis result, for all parameter variations in Table 2, stability robustness is confirmed for uncertain parameters δMAy × δMAz shown in Figure 4. The shaded area in the plot represents the cleared, i.e. stable, parameter combinations. The figure shows that the system is always stable when the uncertain parameters δMAy and δMAz are both positive. In other words, when the model of the aerodynamic moments is smaller than the true aerodynamic moments, the control system is stable.
The system is robustly stable in the presence of the uncertainties ∆MA , if V˙ in Equation (23) is negative semi-definite. Equivalently, the stability domain of the system can be defined as SP ≡ P ∈ RnP V˙ (P ) ≤ 0 = V˙ −1 (]−∞, 0]) , (24) where V˙ (P ) is given by Equation (23), with P = δMAx δMAy δMAz V ρ α α˙ β β˙ p , (25) where V represents velocity, and ρ represents air density. The characterization of SP can be performed by the SIVIA algorithm. The algorithm is implemented using the Interval Toolbox (Zemke, 2002) in the MATLABTM environment. In the stability robustness evaluation, the states of the rotational motion are allowed to vary around the equilibrium point, which is the origin. In this manner, the intervals of the aerodynamic uncertainties, as well as the variations in the states and flight conditions are included. Table 2 gives the variations in the states and flight conditions as part of the vehicle re-entry flight trajectory, which is taken in this case study. The range of Mach number is between M = 3.5 and M = 3.9. Note: Preliminary analysis shows that the stability robustness of the system is insensitive to the uncertainty in the roll aerodynamic moment, δMAx . Together with the parameters in Table 2,
Fig. 4. Analysis result in δMAy × δMAz space The same steps can be performed to analyse other flight regions in the re-entry trajectory. By covering the whole re-entry flight regime with the analysis, the result will show the overall stability robustness of the control system in the whole trajectory.
4.3 Nonlinear Simulation To verify the results, re-entry flight simulations of the vehicle with nonlinear control system are performed. During the simulations, the true aerodynamic moments may have different values than the ones used in the controller. The simulations are performed in SIMULINKTM environment, using General Simulator for Atmospheric Re-entry
In this section, the worst case simulation of the reentry flight is performed by taking the unstable δMAy and δMAz combinations in the model of the controller. For this simulation, the initial conditions are given in Table 3. Table 3. Simulation initial conditions parameter
initial condition
angle of attack, α sideslip angle, β
2 deg 2 deg
aerodynamic roll angle, σ roll rate, p
0 deg 20 deg s−1
pitch rate, q
0 deg s−1
yaw rate, r velocity, V flight path angle, γ heading, χ altitude, h latitude, δ longitude, τ
0 deg s−1 6650 m s−1 -0.2 deg 120 deg 26 · 105 m 77 deg 118 deg
Figure 5 shows the simulation of the system with uncertainty in the worst-case scenario, namely δMAx δMAy δMAz = 0.05 0.05 −0.02 . (26) These worst-case parameters are located in the unstable region of the parameter space (see Figure 4). The system is unstable, as shown by the errors of the outputs p and β which diverge. The errors in the aerodynamic moments propagate throughout the simulation time and eventually caused system instability. This simulation verifies the clearance result shown in the previous section. 5. CONCLUSION The stability robustness evaluation of a re-entry flight control system model with uncertainties in aerodynamic moments, performed using interval analysis, has successfully located the worst-case parameter combinations. These are the ones which give worst stability robustness to the system. Using this approach, one can analyse the nonlinear parameter-dependent model as a continuous system. Nonlinear simulations have been performed
ref
p − p , deg/s
ref
α − α , deg
ref
(1) the nominal system, where the model matches the true system, (2) the system with stable δMAy and δMAz assigned to the model, and (3) the system with unstable δMAy and δMAz assigned to the model.
0.5
β − β , deg
Dynamics (GESARED), a six-degree-of-freedom re-entry simulation toolbox developed in Delft University of Technology (Costa et al., 2002). Three types of simulations can be performed to validate the result of the analysis, i.e. with different values of the aerodynamic moments assigned to the controller. The simulations can be classified as follows:
0 −0.5 400 0.5
410
420
430
440
450
410
420
430
440
450
410
420
430 time, s
440
450
0 −0.5 400 40 20 0 −20 −40 400
Fig. 5. Response errors in the simulation to verify the results of the analysis, which shows that the algorithm has performed the evaluation with reliable results.
REFERENCES Costa, R. R., J. A. Silva, S. F. Wu, Q. P. Chu and J. A. Mulder (2002). Atmospheric Reentry Modeling and Simulation. Journal of Spacecraft and Rockets 39(4), 636–639. Hansen, E. and G.W. Walster (1992). Global Optimization using Interval Analysis. Marcel Dekker, Inc.. USA. Jaulin, L. and E. Walter (1993). Set Inversion via Interval Analysis for Non-linear Bounded Error Estimation. Automatica 29(4), 1053– 1064. Jaulin, L., M. Kieffer, O. Didrit and E. Walter (2001). Applied Interval Analysis. SpringerVerlag. London, GB. Juliana, S., Q. P. Chu, J. A. Mulder and T. J. van Baten (2004). Flight Control of Atmospheric Re-entry Vehicle with Non-linear Dynamic Inversion. In: AIAA Paper 2004-5330. American Institute of Aeronautics and Astronautics. USA. Juliana, S., Q.P. Chu, J.A. Mulder and T.J. van Baten (2005). The Analytical Derivation of Non-linear Dynamic Inversion Control for Parametric Uncertain Systems. In: AIAA Paper 2005-5849. American Institute of Aeronautics and Astronautics. USA. La Salle, J. and S. Lefschetz (1961). Stability by Liapunov’s Direct Method. Academic Press. USA. Slotine, J-J. E. and W. P. Li (1991). Feedback Linearization. In: Applied Non-linear Control. Prentice Hall. USA. pp. 207–271. Verde, L. and F. Corraro (2002). A Polynomialbased Clearance Method. In: Advanced Techniques for Clearance of Flight Control Laws. Springer-Verlag. Germany. pp. 77 – 88. Zemke, J. (2002). b4m, A Free Interval Arithmetic Toolbox for MATLAB based on BIAS, version 1.02.004. Technische Universitaet Hamburg Harburg. Hamburg, Germany.
Faculty of Aerospace Engineering, Delft University of Technology, PO Box 5058, 2600 GB Delft, The Netherlands
Abstract: A stability robustness analysis is performed on a re-entry vehicle flight control system. The mathematical model of this system is nonlinear and may contain uncertain parameters, for which a set of models has to be examined instead of a singleton. The robustness evaluation is performed using interval analysis, which allows one to certify uncertainty regions affecting the closed-loop system. The celebrated Lyapunov function is taken as the stability criterion for the closed-loop system with uncertainties, since the existence of this function provides sufficient stability condition. The result shows the feasibility of the method for c nonlinear system analysis. Copyright 2007 IFAC Keywords: closed-loop system, interval, Lyapunov, robustness, stability
1. INTRODUCTION A re-entry flight vehicle can be modeled mathematically as a highly parameter dependent nonlinear system, for the purpose of controller design and analysis. Because of the nature of the model, nonlinear control techniques are suitable for this type of system to achieve desired stability and performance. In practice, the model may contain uncertainties which may be caused by modeling errors, disturbances from the environment, unknown parameter variation, and measurement errors. As a consequence, the controller may not be able to perform adequately and eventually may cause damage to the system. Uncertainties in re-entry flight models may arise from various sources. For instance, the knowledge about re-entry aerodynamics is still not mature; hence no accurate models are available. Nonethe1 2 3
Researcher Associate Professor Professor, Head of Control and Simulation Division
less, one has to make sure that the designed controllers are robust. To this objective, this paper presents a method to evaluate the stability robustness of re-entry flight control systems in the presence of parametric uncertainties. The Lyapunov stability is chosen as the criterion for the evaluation, for the existence of a Lyapunov function provides sufficient condition for stability in the system. Moreover, the criterion is applicable to both linear and nonlinear input-output systems. A closed-loop system with parametric uncertainties can be considered as a set of model whose numerical values may vary within some unknown but bounded intervals. Therefore, robustness analysis must be applied to this set, instead of to a single model. In this paper, an interval analysis algorithm based on set inversion is proposed to perform this evaluation. Using this approach, all elements in the set can be analysed simultaneously. In (Verde and Corraro, 2002), a similar algorithm was proposed for clearance of flight control law; however, the scope was limited to
linear systems only. In this paper, the algorithm is applied to analyse nonlinear control systems. For completeness, this algorithm is described in Section 3, based on an extensive discussion in (Jaulin et al., 2001; Jaulin and Walter, 1993). To illustrate the method, a nonlinear re-entry flight control system is taken as a case study (Juliana et al., 2004). The controller is designed using the nonlinear dynamic inversion (NDI) technique, which is sensitive to possible model errors (Slotine and Li, 1991). The error can propagate through the time, which eventually can cause system instability. The system under consideration is assumed to contain uncertain parameters in its aerodynamic models. A modification of the NDI technique allows the description of the parametric uncertainties to be included in the model (Juliana et al., 2005). The Lyapunov function evaluation in interval analysis enables one to perform the analysis on a continuous system, omitting the drawback of the classical grid point approach. The result of the robustness analysis, described in Section 4, shows the feasibility of the method to be applied to general nonlinear systems with uncertain parameters. The method is guaranteed to give effective result, for no critical region can be missed in the analysis. The result is validated by performing worst-case simulations of the closed-loop system.
2. LYAPUNOV STABILITY For a general nonlinear autonomous system model, one can construct a positive definite function, which has continuous partial derivatives in the neighbourhood of the equilibria, and whose time derivative along any state trajectory of the system is negative semi-definite. Consider a general nonlinear system model, x˙ = f (x),
(1)
where x is the state vector of the system, and f is a nonlinear continuous function. The system is autonomous, as the function f does not depend explicitly on time. For this system, construct a positive definite function V (x) : Rn → R, V (x) > 0, (2)
The Lyapunov theorem for local stability is defined by Equation (2) and (3) (La Salle and Lefschetz, 1961). When the time derivative of a positive-definite Lyapunov function is negative semi-definite, the system is said to be locally stable. When this derivative is negative definite, the system is said to be locally asymptotically stable. For global stability, one has to consider the whole state-space instead of the neighbourhood of the equilibria. Moreover, the Lyapunov function V (x) has to be radially unbounded. The existence of a Lyapunov function for the system in Equation (1) provides sufficient condition of stability in the system. However, there is generally no effective method to construct a Lyapunov function. Therefore, trial-and-error methods are usually employed to find it. Fortunately, for mechanical systems, the energy dissipation principle can be used to form a suitable choice of the Lyapunov function, as will be shown by a case study.
3. INTERVAL ANALYSIS For the computation with intervals of numbers, the rules are defined differently from the operations for single numbers. The basic rules for arithmetic operations, set operations and inclusion tests can be found in various textbooks and papers on interval analysis such as (Hansen and Walster, 1992; Jaulin et al., 2001). These operations are the basics to build the algorithm Set Inverter via Interval Analysis (SIVIA), which is used to perform the stability robustness evaluation in this paper. This section is dedicated to describe the algorithm.
3.1 Subpavings and Set Inversion The notation of subpavings is defined in the following (Jaulin et al., 2001). Take a set of interest X in Rn . This set is covered with sets of non-overlapping boxes of Rn , which are called the subpavings. The set X is also included between its inner and outer approximations, i.e. X ⊂ X ⊂ X,
(4)
for which the knowledge of the inner and outer approximations provides some valuable information about X.
which has continuous partial derivatives in the neighbourhood of the equilibria. If its time derivative along any state trajectory of the system in Equation (1) is negative semi-definite, i.e.
Set inversion is the basic building block in the algorithm presented in the next section. Let f be a possibly nonlinear function from Rn to Rm , and let Y be a subset of Rm . Set inversion is the computation of the reciprocal image
V˙ (x) ≤ 0,
X = {x ∈ Rn |f (x) ∈ Y} = f −1 (Y),
(3)
then V (x) is called the Lyapunov function for the system in Equation (1).
m
(5)
of a regular subpaving Y ⊂ R by a function f : Rn → Rm . The algorithm serves as a method
to compute two subpavings X and X of Rn , such that Equation (4) is satisfied.
3.2 Set Inversion via Interval Analysis The regular subpavings, such as in Equation (4) for X, for any function f with a convergent inclusion function [f ], can be obtained with the set inversion algorithm. This algorithm is illustrated in Figure 1.
of the parameter space of p, for which the derivative of the Lyapunov function is negative semidefinite. In other words, the stability domain can be defined as Sp ≡ p ∈ Rnp V˙ (p) ≤ 0 = V˙ −1 (]−∞, 0]) . (6) The characterization of Sp can, therefore, be cast into the framework of set inversion and performed by the SIVIA algorithm.
4. CONTROL SYSTEM APPLICATION The algorithm is applied to analyse the stability robustness of a nonlinear re-entry flight control system. As a case study, a re-entry capsule equipped with NDI controller is taken from (Juliana et al., 2004). The block diagram of the closed-loop system is shown in Figure 2. While the inner loop feedback control is designed using NDI, in the outer feedback loop, the Proportional, Integral, and Derivative (PID) technique is applied. The vehicle is spun around its body axis to provide gyroscopic rigidity. The controller is designed to keep the spinning rate p constant, and to keep angle of attack α and sideslip angle β zero. Fig. 1. Illustration of set inversion via interval analysis As shown in Figure 1, four cases may be encountered (denoted by Roman numerals): (1) Case I: [f ] ([x]) has a non-empty intersection with Y, but is not entirely in Y. The set [x] then may contain part of the solution set; it is said to be undetermined. If the magnitude of [x] is larger than a predefined small number , the algorithm bisects the interval. The test is further performed recursively on the newly generated intervals. (2) Case II: [f ] ([x]) has an empty intersection with Y. This means that [x] does not belong to the solution set. (3) Case III: [f ] ([x]) is entirely in Y. Therefore, [x] belongs to the solution, and stored in X and X. (4) Case IV: the interval [x] is undetermined, while its width is smaller than . As the interval is small enough, it can be stored in the outer approximation X of X. The stability robustness evaluation can be cast into a set inversion problem. As shown by Equation (3), the derivative of the positive-definite Lyapunov function has to be negative semi definite for the system to be locally stable. When a system contains parametric uncertainties, denoted as p, the Lyapunov function may also contain these uncertainties. The robustness evaluation amounts to find the stability domain Sp , which is a subset
Fig. 2. Feedback control system with NDI and PID
4.1 Closed-loop system with uncertainties In the following, the model of the nonlinear reentry flight control system in the presence of parametric uncertainties is described. The model of the vehicle rotational motion can be written as a set of first-order nonlinear differential equations. The first-order derivatives of the outputs y = T p α β are written in Equation (7) through (9), where q and r are pitch rate and roll rate, Ixx is the moment of inertia in x-axis, MCx is the control moment about x-axis, and MAx is the aerodynamic moment about x-axis. The expressions for the second order derivative of α and β are given in Equation (10) and (11). Equation (7), (10) and (11) are used to perform the dynamic inversion. In the model construction, the magnitude of the aerodynamic moments is assumed to be not fully known, with the uncertainties assumed to be bounded. In a linear fractional representation form, the true aerodynamic moments can be expressed as functions of predicted aerodynamic moments and the uncertainties, as shown
MCx + MAx , Ixx α˙ = −p cos α tan β + q − r sin α tan β, β˙ = p sin α − r cos α, p˙ =
(7) (8) (9)
α ¨ = α˙ (p sin α tan β − r cos α tan β) − β˙ p cos α sec2 β + r sin α sec2 β − cos α tan β (MCx + MAx ) /Ixx (MCy + MAy ) + (Izz − Ixx ) · pr sin α tan β (MCz + MAz ) + (Iyy − Ixx ) · pq sin α tan β + − , (10) Iyy Izz cos α (MCz + MAz ) + (Iyy − Ixx ) · pq cos α sin α (MCx + MAx ) − . (11) β¨ = α˙ (p cos α + r sin α) + Ixx Izz ˜ A, (12) MA = I3×3 + ∆MA · M where I3×3 is a 3-by-3 identity matrix, and ⎡ ⎤ δMAx 0 0 T ˜A = M ˜ Ax M ˜ Ay M ˜ Az T . ∆MA = ⎣ 0 δMAy 0 ⎦ ; MA = MAx MAy MAz ; M 0 0 δMAz Fig. 3. Mathematical model of the re-entry vehicle with aerodynamic uncertainties in Equation (12). The terms δMAx , δMAy , and δMAz represent the lump of all unknowns in the model of aerodynamic moments. Their bounds show the relative values of the uncertainties with respect to the predicted aerodynamic moments. For instance, ˜ Ax 1 + δM MAx = M . (13) Ax
bounds are the first guess of the uncertainties, which can be shrinked or enlarged based on the result of the analysis. Table 1. Uncertain parameters
In Equation (13), the true moment MAx is equal ˜ Ax only if δM = 0. to the predicted moment M Ax Otherwise, the magnitude of MAx will be larger ˜ Ax , depending on δM . or smaller than M Ax As a result of imperfect feedback linearization, the closed-loop system model still contains the uncertain nonlinear terms given by Equation (12). This model can be expressed as ⎡ ⎤ ˜ Ax /Ixx δMAx M ⎡ ⎤ ⎡ ⎤ p˙ ν1 ⎢ ⎥ ˜ ⎥ ⎣α ¨ ⎦ = ⎣ν2 ⎦+ A1 · ⎢ ⎢δMAy MAy /Iyy ⎥ , ⎣ ⎦ ν3 β¨ ˜ /I δ M MAz
with A1
Az
(14)
zz
⎤ ⎡ cos β 0 0 1 ⎣ − cos α sin β cos β − sin α sin β ⎦, = cos β sin α cos β 0 − cos α cos β
and ν1 = (KP 1 + KI1 /s) (pref − p) , ν2 = (KP 2 + KD2 s) (αref − α) , ν3 = (KP 3 + KD3 s) (βref − β) , where the terms KP , KI , and KD are the proportional, integral, and derivative gains, respectively. The stability robustness evaluation is performed to the closed-loop system model in Equation (14) in the presence of the aerodynamic uncertainties, whose values are given in Table 1. The 10%
parameter
bound
δMAx δMAy
[−0.05, 0.05] [−0.05, 0.05]
δMAz
[−0.05, 0.05]
4.2 Stability robustness analysis The form of the closed-loop system may be simplified by defining a new state vector, T ∂q q˙ = = p α˙ β˙ , (15) ∂t and positive definite matrices KI and KII , whose non-zero elements are the PID controller gains, ⎡ ⎤ KP 1 0 0 KI = ⎣ 0 KD2 0 ⎦ (16) 0 0 KD3 ⎡
⎤ KI1 0 0 KII = ⎣ 0 KP 2 0 ⎦ . 0 0 KP 3
and
(17)
Substituting Equation (15) through (17) into (14), and removing the static term from the model, the closed-loop system can be rewritten as ˜ A, q¨ = −KI q˙ − KII q + A2 In−1 ∆MA M with
(18)
⎤ 1 0 0 A2 = ⎣ − cos α tan β 1 − sin α tan β ⎦ ; and sin α 0 − cos α ⎡
⎡
In−1
⎤ 1/Ixx 0 0 = ⎣ 0 1/Iyy 0 ⎦ . 0 0 1/Izz
Table 2. Parameter variation parameter
A Lyapunov function can be constructed as the sum of the artificial kinetic and potential energy functions of the system. The artificial kinetic energy for the system in Equation (18) can be written as a quadratic function of the body rate, 1 T q˙ q, ˙ (19) 2 and the artificial potential energy of the system can be written as 1 T q KII q . Ep = (20) 2 The Lyapunov function is, therefore, 1 T q˙ q˙ + q T KII q , (21) V = 2 which is positive semi-definite. Its first-order derivative can be written as Ek =
V˙ = q˙T q¨ + q˙T KII q,
(22)
which by substituting Equation (18) becomes ˜ A. V˙ = −q˙T KI q˙ + q˙T AIn−1 ∆MA M
(23)
The values of V˙ are functions of the flight conditions, the states, and the uncertain parameters.
interval
α
[−1, 0) ∪ (0, 1] deg
α˙ β β˙
[−5.6, 0) ∪ (0, 5.6] deg/s [−1, 0) ∪ (0, 1] deg
p V ρ
[−5.6, 0) ∪ (0, 5.6] deg/s [−1, 1] deg/s [1041, 1175] m/s [2.4, 2.7] · 10−2 kg /m3
δMAx is regarded as a fix interval during the analysis. This reduces the computational workload. The stable domain SP for the closed-loop system is a 10-dimensional space, whose dimensions are given by the parameters P in Equation (25), in which Equation (24) is satisfied. Based on the analysis result, for all parameter variations in Table 2, stability robustness is confirmed for uncertain parameters δMAy × δMAz shown in Figure 4. The shaded area in the plot represents the cleared, i.e. stable, parameter combinations. The figure shows that the system is always stable when the uncertain parameters δMAy and δMAz are both positive. In other words, when the model of the aerodynamic moments is smaller than the true aerodynamic moments, the control system is stable.
The system is robustly stable in the presence of the uncertainties ∆MA , if V˙ in Equation (23) is negative semi-definite. Equivalently, the stability domain of the system can be defined as SP ≡ P ∈ RnP V˙ (P ) ≤ 0 = V˙ −1 (]−∞, 0]) , (24) where V˙ (P ) is given by Equation (23), with P = δMAx δMAy δMAz V ρ α α˙ β β˙ p , (25) where V represents velocity, and ρ represents air density. The characterization of SP can be performed by the SIVIA algorithm. The algorithm is implemented using the Interval Toolbox (Zemke, 2002) in the MATLABTM environment. In the stability robustness evaluation, the states of the rotational motion are allowed to vary around the equilibrium point, which is the origin. In this manner, the intervals of the aerodynamic uncertainties, as well as the variations in the states and flight conditions are included. Table 2 gives the variations in the states and flight conditions as part of the vehicle re-entry flight trajectory, which is taken in this case study. The range of Mach number is between M = 3.5 and M = 3.9. Note: Preliminary analysis shows that the stability robustness of the system is insensitive to the uncertainty in the roll aerodynamic moment, δMAx . Together with the parameters in Table 2,
Fig. 4. Analysis result in δMAy × δMAz space The same steps can be performed to analyse other flight regions in the re-entry trajectory. By covering the whole re-entry flight regime with the analysis, the result will show the overall stability robustness of the control system in the whole trajectory.
4.3 Nonlinear Simulation To verify the results, re-entry flight simulations of the vehicle with nonlinear control system are performed. During the simulations, the true aerodynamic moments may have different values than the ones used in the controller. The simulations are performed in SIMULINKTM environment, using General Simulator for Atmospheric Re-entry
In this section, the worst case simulation of the reentry flight is performed by taking the unstable δMAy and δMAz combinations in the model of the controller. For this simulation, the initial conditions are given in Table 3. Table 3. Simulation initial conditions parameter
initial condition
angle of attack, α sideslip angle, β
2 deg 2 deg
aerodynamic roll angle, σ roll rate, p
0 deg 20 deg s−1
pitch rate, q
0 deg s−1
yaw rate, r velocity, V flight path angle, γ heading, χ altitude, h latitude, δ longitude, τ
0 deg s−1 6650 m s−1 -0.2 deg 120 deg 26 · 105 m 77 deg 118 deg
Figure 5 shows the simulation of the system with uncertainty in the worst-case scenario, namely δMAx δMAy δMAz = 0.05 0.05 −0.02 . (26) These worst-case parameters are located in the unstable region of the parameter space (see Figure 4). The system is unstable, as shown by the errors of the outputs p and β which diverge. The errors in the aerodynamic moments propagate throughout the simulation time and eventually caused system instability. This simulation verifies the clearance result shown in the previous section. 5. CONCLUSION The stability robustness evaluation of a re-entry flight control system model with uncertainties in aerodynamic moments, performed using interval analysis, has successfully located the worst-case parameter combinations. These are the ones which give worst stability robustness to the system. Using this approach, one can analyse the nonlinear parameter-dependent model as a continuous system. Nonlinear simulations have been performed
ref
p − p , deg/s
ref
α − α , deg
ref
(1) the nominal system, where the model matches the true system, (2) the system with stable δMAy and δMAz assigned to the model, and (3) the system with unstable δMAy and δMAz assigned to the model.
0.5
β − β , deg
Dynamics (GESARED), a six-degree-of-freedom re-entry simulation toolbox developed in Delft University of Technology (Costa et al., 2002). Three types of simulations can be performed to validate the result of the analysis, i.e. with different values of the aerodynamic moments assigned to the controller. The simulations can be classified as follows:
0 −0.5 400 0.5
410
420
430
440
450
410
420
430
440
450
410
420
430 time, s
440
450
0 −0.5 400 40 20 0 −20 −40 400
Fig. 5. Response errors in the simulation to verify the results of the analysis, which shows that the algorithm has performed the evaluation with reliable results.
REFERENCES Costa, R. R., J. A. Silva, S. F. Wu, Q. P. Chu and J. A. Mulder (2002). Atmospheric Reentry Modeling and Simulation. Journal of Spacecraft and Rockets 39(4), 636–639. Hansen, E. and G.W. Walster (1992). Global Optimization using Interval Analysis. Marcel Dekker, Inc.. USA. Jaulin, L. and E. Walter (1993). Set Inversion via Interval Analysis for Non-linear Bounded Error Estimation. Automatica 29(4), 1053– 1064. Jaulin, L., M. Kieffer, O. Didrit and E. Walter (2001). Applied Interval Analysis. SpringerVerlag. London, GB. Juliana, S., Q. P. Chu, J. A. Mulder and T. J. van Baten (2004). Flight Control of Atmospheric Re-entry Vehicle with Non-linear Dynamic Inversion. In: AIAA Paper 2004-5330. American Institute of Aeronautics and Astronautics. USA. Juliana, S., Q.P. Chu, J.A. Mulder and T.J. van Baten (2005). The Analytical Derivation of Non-linear Dynamic Inversion Control for Parametric Uncertain Systems. In: AIAA Paper 2005-5849. American Institute of Aeronautics and Astronautics. USA. La Salle, J. and S. Lefschetz (1961). Stability by Liapunov’s Direct Method. Academic Press. USA. Slotine, J-J. E. and W. P. Li (1991). Feedback Linearization. In: Applied Non-linear Control. Prentice Hall. USA. pp. 207–271. Verde, L. and F. Corraro (2002). A Polynomialbased Clearance Method. In: Advanced Techniques for Clearance of Flight Control Laws. Springer-Verlag. Germany. pp. 77 – 88. Zemke, J. (2002). b4m, A Free Interval Arithmetic Toolbox for MATLAB based on BIAS, version 1.02.004. Technische Universitaet Hamburg Harburg. Hamburg, Germany.