SUM-OF-SQUARES BASED STABILITY ANALYSIS FOR SKID-TO-TURN MISSILES WITH THREE-LOOP AUTOPILOT Hyuck-Hoon Kwon*, Min-Won Seo*, Dae-Sung Jang*, and Han-Lim Choi *Department of Aerospace Engineering, KAIST, Daejeon 305-701, Republic of Korea E-mail : {hhkwon, mwseo, dsjang, hanlimc}@lics.kaist.ac.kr Abstract In recent years, sum-of-squares optimization has been employed to estimate region of attraction for nonlinear systems. We present region of attraction for missiles with three-loop autopilot in boost phase and give variation of region of attraction according to gainscheduling. Comparing the region of attraction between them can be used for qualitative insight to check the stability and effect of gainscheduling method. 1 Introduction Conventional autopilot design for tactical guided missiles employs acceleration and rate feedback such as three-loop structure to track the guidance command [1]. The gains of the linear autopilot are obtained at several trim points indexed by scheduling variable for gainscheduling method. However, linear design method about equilibrium point ensure stability and performance in only small area and may even lose stability between trim points due to the nonlinear characteristics of missiles. In the design of linear control, control designers don’t know the exact or approximate region of attraction, and they should obtain as many trim points as possible to guarantee the stability for entire flight regime. Recently, significant development has been put on the nonlinear stability analysis based on sumof-squares optimization [2-4]. They use sum-ofsquares relaxation to check non-negativity and efficient solve them using semi-definite programming algorithm [5]. If nonlinear system can be expressed by polynomial functions, then sum-of-square optimization is used to construct Lyapunov function or to give approximate region of attraction based on nonlinear dynamics [6-8]. However, it can be applied to only polynomial systems. Then nonlinear systems with non-polynomial terms should be transformed or approximated to polynomial form. In this paper, regions of attraction for missile in boost phase are calculated using sum-of-squares optimization tool. Standard three-loop autopilot is applied to the missile and is scheduled by Mach number. At first, the region of attraction is obtained at each trim point based on the result of linear control design. And then, gains by gain-scheduling method are obtained between trim points and used for estimating gainscheduled region of attraction. It can gives insight for stability of nonlinear system and show the effect of gain-scheduling based on nonlinear dynamics, not linear one. Generally, validity for gain-scheduling can be computed using stability margins at some check points between trim points. Since it is also based on linearized systems at check points, it cannot include the exact nonlinear nature. Comparison between the region of attraction at check points gives valuable results to check the stability and effect of gain-scheduling. This paper is organized as follows : Section 2 describes the Lyapunov stability, region of attraction and V-s iteration method to apply SOS optimization. Section 3 introduces pitch dynamics of missiles, especially short-period mode. Section 4 include the result of linear control design such as calculation of trim points, design of three-loop autopilot. Section 5 describes simulation results for estimating and comparing region of attraction. Finally, Section 6 presents conclusion. 1 Hyuck-Hoon Kwon, Min-Won Seo property in practice. The existence of a Lyapunov function is a sufficient condition for stability of the zero equilibrium as shown in the following theorem. 2 Lyapunov Stability and Region of Attraction Consider the autonomous nonlinear system ẋ = f(x) (1) where x(t) ∈ ℝ is state vector and f ∶ ℝ → ℝ is continuous such that f(0) = 0 , i.e., the origin is an equilibrium point of (1), and f is locally Lipschitz. Let (; x ) denote the solution to (1) at time t with the initial state x = x(0) . Without loss of generality, any equilibrium point can be shifted to the origin 0 via a change of variables and we may always assume that the equilibrium point of interest occurs at the origin. Definition 1 (Lyapunov stability). The equilibrium point 0 of (1) is · Stable if, for any ϵ > 0 , there exists δ = δ(ϵ) > 0 such that ‖x(0)‖ < δ ⇒ ; x(0) < ϵ, ∀t > 0 · · unstable, if it is not stable. asymptotically stable, if it is stable, and δ can be chosen such that ‖x(0)‖ < δ ⇒ lim (; x(0)) = 0 → · (2) (3) globally asymptotically stable if it is stable, and for all x(0) ∈ ℝ , lim (; x(0)) = 0 → (4) The equilibrium point 0 is stable if all solutions with the initial conditions in a neighborhood of the origin remain near the origin for all time. Asymptotically stable means that all solutions starting at nearby points not only remain close enough but also finally converge to the equilibrium point. And globally asymptotically stable means that asymptotically stable for all initial state x(0) ∈ ℝ . Since arbitrarily small perturbations of the initial state about the equilibrium point lead to arbitrarily small perturbations in the corresponding solution trajectories of (1), the stability is an important Theorem 2 ([10, Theorem 4.1]) Let ⊂ be a domain containing the equilibrium point = 0 of the system (1). Let : → ℝ be a continuously differentiable function such that (0) = 0, () > 0 \{0} (5) ̇ ≔ ∇ ∗ ≤ 0 (6) then the origin is stable. Moreover, if ̇ ≔ ∗ < 0 \{0} (7) then the origin is asymptotically stable. It is commonly known as a Lyapunov function satisfying conditions (5) and (6) in Theorem 2. And globally asymptotic stability of system (1) can be verified by using Lyapunov functions stated as follows. Theorem 3 ([10, Theorem 4.2]) Let the origin be an equilibrium point for (1). If there exists a continuously differentiable function : ℝ → ℝ such that (0) = 0, () > 0 ∀ ≠ 0, (8) ‖‖ → ∞ ⇒ () → ∞, (9) ̇ () < 0 ∀ ≠ 0 (10) then the origin is globally asymptotically stable. Remark that () satisfying condition (9) is radially unbounded. Although globally asymptotic stability is very desirable, it is difficult to achieve in many applications. Quite often, determining a given system has an asymptotically stable is not sufficient. It is important to determining how far from the origin the trajectory can be and converge to the 2 SUM-OF-SQUARES BASED STABILITY ANALYSIS FOR SKID-TO-TURN MISSILES WITH THREE-LOOP AUTOPILOT origin as approaches infinity. This gives rise to the following definition. Definition 4 (region of attraction). The region of attraction (ROA) Ω of the equilibrium point 0 of (1) is defined as Ω = {x ∈ ℝ | lim (; x) = 0} → The ROA is the set of all points x such that any trajectory starting at initial state x(0) will converges to the equilibrium point. In the literature, the terms “attraction basin” and “domain of attraction” are also used. Definition 5 (positively invariant set). A set M ⊂ ℝ is called a positively invariant set of the system (1), if x(0) ∈ implies (; x(0)) ∈ for all ≥ 0 . Namely, if a solution belongs to a positively invariant set M at some time instant, then it belongs to for all future time. It is difficult to analytically find the exact ROA for nonlinear systems if not impossible. In general, Lyapunov functions can be used to compute estimates of the ROA. Lyapunov based approaches to ROA estimation using a characterization of invariant subsets of the ROA rely on statements of the following lemma. For c > 0 and a function : ℝ → ℝ, define the csublevel set Ω, of as Ω, ≔ {x ∈ ℝ |V(x) < c}. Lemma 1 ([13]) Let ∈ ℝ be positive. If there exists a function : ℝ → ℝ such that , is bounded, and (0) = 0, () > 0 for all ∈ ℝ (11) In order to find the estimate invariant subset of the ROA, we describe P with a semi-algebraic set P ≔ { ∈ |() ≤ } where ∈ , a fixed positive definite convex polynomial, and maximizing subject to P ∈ Ω, satisfying the constraints (11)-(13). We pose the following optimization to search for Lyapunov function ∗ () ≔ max ,∈ s.t. (11)-(13), P ∈ Ω, Here denotes over which the maximum is defined the set of candidate Lyapunov functions. It is shown that a characterization of the invariant subsets of the ROA with regard to the sublevel sets Lyapunov functions. Since the optimization problem in (14) is an infinitedimensional problem, is restrained all polynomials of some fixed degree. Using simple generalizations of the S-procedure [11], sufficient conditions for set containment constraints are obtained. Usually, we take l (x) of the form l (x) = ∑ ϵ , where ϵ are positive small real number. Then, the constraint V − l ∈ SOS (15) And V(0) = 0 are sufficient conditions for (11) and (12). Additionally, if s ∈ SOS, then −[(β − p)s + (V − c)] ∈ SOS (16) means the set containment P ⊆ Ω, , and if s , s ∈ SOS, then (12) −[(c − V)s + ∗ + l ] ∈ SOS , \{0} ⊂ { ∈ ℝ | ∗ < 0}, (14) (17) (13) then for all x(0) ∈ , , the solution of (1) exists, satisfies (; x(0)) ∈ , for all ≥ 0, and → (; x(0)) = 0 , i.e,. , is a positively invariant region contained in the equilibrium point’s domain of attraction. is a sufficient condition for (13). Using these sufficient conditions, a lower bound on β∗ (ν) can be defined as 3 Hyuck-Hoon Kwon, Min-Won Seo ∗ () ≔ max ,∈, ∈ (18) s.t. (15)-(17), (0) = 0, s ∈ , > 0 Here, the set ν is specified finite-dimensional subspaces of polynomials. Even though β∗ depends on these subspaces, it will not always be definitively notated. Remark that β∗ () ≤ ∗ () because conditions (15)-(17) are only sufficient conditions. The optimization problem in (18) is bilinear due to existing the product terms βs in (16) and Vs and ∗ in (17). However, there are so many approaches to solve BMI problem like v-s iteration [2], BMI [8], and using simulation data [14]. In this paper v-s iteration approach is used to find Lyapunov function and estimate ROA. 3 Analysis of a controller for missile pitch dynamics 3.1 Missile longitudinal dynamics A nonlinear short period mode in the pitch dynamic model of a SRAAM (Short-Range Airto-Air Missile) is given by − ̇ = { + 2 ̇ = − The autopilot an automatic control system is to guarantee stability and follow command given by guidance law. In most skid-to-turn configuration, tracking acceleration normal to the missile longitudinal axis is desired. In this paper, assume that general three-loop autopilot is applied to control SRAAM in the boost phase as shown in Figure 1. Three-loop autopilot is comprised of a rate loop, a synthetic stability loop, and an accelerometer feedback loop. The three autopilot gains K , w and K must be chosen to fulfill some designer-chosen criteria and the gain K is obtained from other gains therefore the achieved acceleration will correspond the commanded acceleration. Figure 1. Structure of Three-loop Autopilot The linear control input based on Three-loop structure is shown below. = ∆ + K ∆ (20) where = ∫{( − ) + } . The closed-loop system is combined nonlinear pitchaxis model of (19) and control input of (20). (19) , − − + + } where α is the angle of attack, q is the pitch rate, δq is the control surface deflection, V is the missile velocity, Q, S, D, m are the dynamic pressure, reference area, length, missile mass and I is moment of inertia about the body frame. The aerodynamic force and moment coefficients (C , C , C , C , C ) are given in terms of an aerodynamic table. 3.2 State-space description and Three-loop autopilot design () − () ∆̇ = { () + ∆ 2 ∆̇ = ∆ − − ( , − ) () + () + () )} ∆̇ = () + () + − The model is valid for an equilibrium point characterized by the altitude h = 2km and 4 (21) SUM-OF-SQUARES BASED STABILITY ANALYSIS FOR SKID-TO-TURN MISSILES WITH THREE-LOOP AUTOPILOT number of Mach M . The signal ∆ = − , ∆ = − and ∆ = − represent perturbation from the equilibrium values , and . The terms C (), C (), C (), C () are polynomial functions of and were obtained with polynomial fitting on data of the aerodynamics table as shown in Figure 2. CM -0.3019 -0.1537 -0.0906 -0.0571 -0.0567 -0.0607 Angle of attack and deflection angle at trim points according to Mach number are shown in Figure 3 and Figure 4. As shown in the figures, nonlinear characteristics of missile dynamics is evident from Mach 0.85 to Mach 1.20, which is general transonic region. Therefore, the stability analysis in this area is more important than any other area and it will be analyzed through estimation of region of attraction in the following chapter. 0 0 2-Polyfit -10 1.5668 1.1093 0.8362 0.6468 0.5352 0.4332 1.60 2.00 2.40 2.80 3.20 3.60 -20 CM 0 -30 -40 -50 -60 -70 7 0 10 20 30 40 50 AOA(deg) 60 70 80 90 6 0 5 a TRIM -20 3 dq -40 CM 4 2 -60 1 -80 0 0.5 -100 CM 0 10 20 30 40 50 AOA(rad) 60 70 1.5 2 2.5 3 3.5 4 Mach Figure 3. Angle of Attack at Trim Points 3-Polyfit -120 1 dq 80 90 Figure 2. Polynomial Fitting on Aero. Data 0 -0.5 4 Linear Controller Design For the design of linear controller, trim points should be calculated. Trim points for the levelflight are shown in Table 1. (The gravitational effect is neglected) Table 1. Calculation of Trim Points Mach 0.85 0.95 1.05 1.20 AOA (deg) 5.8254 4.6989 3.7948 2.8231 Deflection (deg) -2.0721 -2.0523 -2.1989 -1.5251 dTRIM -1 -1.5 -2 -2.5 0.5 1 1.5 2 2.5 3 3.5 4 Mach Figure 4. Deflection Angle at Trim Points The structure of general three-loop autopilot is already shown in Figure 1. Then, the gains of 5 Hyuck-Hoon Kwon, Min-Won Seo three-loop autopilot are calculated at trim points to meet the stability margin requirement. In this paper, we obtain at least 6dB gain margin and 45deg phase margin as follows : range for one side become larger, but admissible range for the other side decreases. At each trim point, the Lyapunov functions are calculated as follows : Table 2. Gains of Three-loop Autopilot V . = 4.6847x + 0.3791 + 0.0460 Ω . ≔ {x ∈ ℝ |V. < c . = 1.0042} V . = 5.2914x + 0.4702 + 0.0513 Ω . ≔ {x ∈ ℝ |V . < c . = 1.0033} V . = 4.9698x + 0.3599 + 0.0538 Ω . ≔ {x ∈ ℝ |V . < c . = 1.0042} V . = 4.7734x + 0.2314 + 0.0516 Ω . ≔ {x ∈ ℝ |V . < c . = 1.0014} Mach KDC 1.0735 1.0725 1.0645 1.0599 1.1122 1.1318 1.1486 1.1666 1.1609 1.1507 0.85 0.95 1.05 1.20 1.60 2.00 2.40 2.80 3.20 3.60 KA 0.0481 0.0437 0.0444 0.0274 0.0167 0.0114 0.0084 0.0064 0.0058 0.0055 wI 8.1835 8.6233 8.4221 8.4447 10.0004 10.5640 10.8591 11.0255 10.8589 10.5670 KR 0.2605 0.2024 0.1426 0.1236 0.0827 0.0661 0.0541 0.0445 0.0384 0.0335 5 Simulation Results In this paper, we concentrate on the analysis of trim points in the transonic region due to limited space. 5.1 Estimation of Region of Attraction Simulation results from Mach 0.85 to Mach 1.20 are shown in the figure 5. It appears that each trim point has different region of attraction. M=0.85 M=0.95 M=1.05 M=1.20 5 4 5.2 ROA Effect of Gain-Scheduling Three-loop autopilot gain set is calculated at each trim point as shown in Table2. In the region between trim points, gain-scheduling is generally used for linear control methodology. Variation of region of attraction by gainscheduling reveals the difference between nonlinear analysis and linear analysis. Figure 6 indicate the results at Mach 0.95 for three cases. Blue line is the result for the proper gain set calculated at Mach 0.95. Red line means the result for linear interpolation gain set. Finally, green line is the result for gain set calculated at Mach 1.05. As shown in the figure, blue region is larger than red and green region and red region is larger than green one. It means that interpolation method is better than using improper gain set in view of ROA. 3 2 Proper gain set at Mach 0.95 Gain-scheduled gain set at Mach 0.95 Improper gain set at Mach 1.05 q(rad/sec) 5 1 4 0 3 -1 2 q(rad/sec) -2 -3 -4 1 0 -1 -5 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 -2 AOA(rad) -3 Figure 5. Region of Attraction at Trim Points -4 -5 Region of attraction at Mach 0.95 is smaller than that at Mach 0.85 for angle of attack and pitch rate. For Mach 1.05 and 1.20, admissible -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 AOA(rad) Figure 6. Region of Attraction at Mach 0.95 6 SUM-OF-SQUARES BASED STABILITY ANALYSIS FOR SKID-TO-TURN MISSILES WITH THREE-LOOP AUTOPILOT Then, the results of linear analysis for three cases are as follows : Table 3. Linear Results of Mach 0.95 Gain Set Gain Margin Phase Margin 0.95 1.05 Scheduled 12.3 dB 15.3 dB 12.4 dB 46.0 deg 44.9 deg 46.5 deg Gain margin by proper gain set is the smallest among them and phase margin by improper gain set is the smallest one. Gain margin is expected to be similar to the result of ROA. However, the simulation results display that linear analysis don’t express the nature of nonlinear systems properly. The figure 7 is the simulation results for Mach 1.05. It is similar to those in figure 6, but region of attraction by gain-scheduling method is a little bit large in some axis. It means that stability margin in each state of nonlinear system is different Proper gain set at Mach 1.05 Gain-scheduled gain set at Mach 1.05 Improper gain set at Mach 1.20 5 4 3 q(rad/sec) 2 1 0 -1 -2 -3 -4 -5 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 The results by linear analysis in Table 4 are hard to compare the nonlinear one in figure 7. It represents that linear dynamics at the trim points is sometimes highly nonlinear or results from linear analysis don’t capture the nonlinear nature in terms of stability. In conclusion, control designers should consider nonlinear analysis in addition to traditional linear one to guarantee the stability of nonlinear systems. 6 Conclusion In this paper, region of attraction for missile with three-loop autopilot is calculated based on sum-of-squares optimization. For simplicity, the short period mode in pitch dynamics of missiles is considered. Aerodynamic data in look-up table is transformed into 2nd or 3rd polynomial functions for use of sum-of-squares method. The region of attraction at transonic trim points are calculated and compared with the results of gain-scheduled gain set and improper gain set. Simulation results generally show that the region of attraction with better gain set is larger than any other region of attraction in most direction. However, Nonlinear region of attraction seems to be hard to understand based on the results of linear control during highly nonlinear regime. As a result, although stability margin on linear dynamics reflect stability features on nonlinear dynamics in a measure, it doesn’t contain nonlinear property completely. Therefore, nonlinear analysis need to be carried out with linear analysis for a better understanding of nonlinear systems. 0.5 AOA(rad) Figure 7. Region of Attraction at Mach 1.05 The results of linear analysis for Mach 1.05 are as follows : Acknowledgments This work was supported by “Core Technology Research for the Next Generation of Precesion Guided Munitions” Project funded by LIG Nex1. Table 4. Linear Results of Mach 1.05 Gain Set Gain Margin Phase Margin 1.05 1.20 Scheduled 11.4 dB 11.9 dB 9.53 dB 46.1 deg 48.6 deg 44.2 deg References [1] ZARCHAN P. Tactical and Strategic Missile Guidance. Fifth edition, AIAA, 2007. [2] Z. Jarvis-Wloszek, Lyapunov Based Analysis and Controller Synthesis for Polynomial Systems using 7 Hyuck-Hoon Kwon, Min-Won Seo Sum-of-Squares Optimization, Ph.D. dissertation, University of California, Berkeley, 2003. [3] A. Chakraborty, P. Seiler, and G. Balas, “Nonlinear Region of Attraction Analysis for Flight Control Verification and Validation,” Control Engineering Practice, Vol.19, Issue 4, April 2011, Pages 335-345. [4] A. Ataei-Esfahani, Algorithms for Sum-of-SquaresBased Stability Analysis and Control Design of Uncertain Nonlinear Systems, Ph.D. dissertation, The Pennsylvania State University, 2012. [5] P. Parrilo, Structured Semidefinite Programs and Semialgebraic Geometry Methods in Robustness and Optimization, Ph.D. Thesis, California Inst. Of Technology, Pasadena, CA, 2000 [6] Antonis P. and Stephen P., “A Tutorial on Sum of Squares Techniques for Systems Analysis,” American Control Conference, June 8-10, 2005. [7] Antonis P. and Stephen P., “On the Construction of Lyapunov Functions using the Sum of Squares Decompositionm,” Proceedings of the 41st IEEE Conference on Decision and Control, Dec., 2002. [8] Stephen P., Antonis P., and Fen W., “Nonlinear Control Synthesis by Sum of Squares Optimization: A Lyapunov-based Approach,” 5th Asian Control Conference, 2004. [9] Sturm J. Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones. Optimization Methods and Software, 11, 625-653, 1999. [10] H. K. Khalil. Nonlinear Systems. Prentice Hall, 3rd edition, 2002. [11] P. Parrilo. Semidefinite programming relaxations for semialgebraic problems. Mathematical Programming Ser. B, 96:2:293-321, 2003. [12] Tan W, Nonlinear Control Analysis and Synthesis Using Sum-of-Squares Programming, Ph.D. Thesis, University of California, Berkeley, 2006. [13] Vidyasagar. M. Nonlinear Systems analysis. Prentice Hall, 2nd edition, 1993. [14] U. Topcu, A. Packard, and P. Seiler, “Local stability analysis using simulations and sum-of-squares programming,” Automatica, vol44, no.10, pp.26692675, 2008. Copyright Statement The authors confirm that they, and/or their company or organization, hold copyright on all of the original material included in this paper. The authors also confirm that they have obtained permission, from the copyright holder of any third party material included in this paper, to publish it as part of their paper. The authors confirm that they give permission, or have obtained permission from the copyright holder of this paper, for the publication and distribution of this paper as part of the ICAS 2014 proceedings or as individual off-prints from the proceedings. 8