Prediction of limit cycle oscillations in piecewise linear systems with multiple piecewise nonlinearities

Mathematical modelings of many electric and mechanical systems involve piecewise linear system. Piecewise linear system can possess a periodic solution called a limit-cycle oscillation (LCO), which can seriously undermine the system performance. Therefore, how to analyze LCO and its parameters in piecewise linear systems is one of the primary concerns for the control and system engineers. This work presents a novel framework to predict and analyze LCO of piecewise linear systems, focused on systems with multiple piecewise nonlinearities. On top of the well-known piecewise linear analysis, the Floquet theory is applied to identify LCO parameters and determine the stability of the LCO. The introduction of Floquet theory to piecewise linear systems is allowed through transforming piecewise nonlinearities to corresponding equivalent analytic functions. In addition, the establishment of switching equation provides another necessary condition to predict LCO parameters. An example of a realistic ﬂight control system is taken to demonstrate the effectiveness and efﬁciency of authors’ framework.

PNs because the bandwidth of each linear transfer function near each PN can be different.
By piecewise linear systems approach, we can predict the exact frequency of LCO along with the stability of LCO of the piecewise system in which the system dynamics are divided into a finite number of linear affine systems by the PNs [4]. Assuming the unimodal sinusoidal signal as an input to a PN, this approach calculates half the LCO period based on the switching time between surfaces defined by the PN. The stability of the LCO is also verified by the surface Lyapunov function, and in the example of a linear system with a hysteresis relay, the effectiveness of this framework is demonstrated [4]. However, we virtually cannot apply this approach to the system with multiple PNs, because the switching order between surfaces defined by each PN is generally not clear. Therefore, the clarification of the switching order needs to be found in order to successfully apply the piecewise linear approach to the LCO prediction of general piecewise linear systems.
Analysis of LCO for more general nonlinear differential equations also led to methods that identify and describe the LCO. Harmonic balance method (HBM) and energy balance method (EBM) presents the LCO amplitude-period relationship in a simple and accurate framework [5,6,13,14]. But HBM along with EBM are applicable only to conservative systems, which is generally not the case for piecewise linear systems. Furthermore, the solution for each coefficients of harmonic terms are not necessarily feasible and HBM lacks the ability to analyze the stability of LCO. We also have studies on the center manifold theorem and its branches for the analysis of the LCO of general nonlinear differential equation of up to three dimensional state space, which yield excellent analytic results only if calculation of Lyapunov exponent is feasible [2,15,16].
The main contribution of this paper is the construction of a framework for predicting and analyzing limit cycle oscillations in piecewise linear systems with multiple PNs, which includes the followings.
• Equivalent analytic functions (EAF) are designed corresponding to common PNs. The point-wise convergence of them to the original PN function is proven. • The perturbed linear time-periodic (LTP) system evaluated at a periodic solution of piecewise linear systems is established so that Floquet theory may be applied. • The primary characteristic multiplier (PCM) of the perturbed LTP system is calculated leading to the first necessary condition to predict an LCO, or the  equation. • A framework to analyze the switching order of piecewise linear systems is organized. This involves the classification of switching configurations, switching time identification, and the establishment of switching equation, or the  equation. • The absolute stability of every Lur'e type system corresponding to each switching configuration is proven based on Popov Criteria and KYP lemma [17]. Showing that at least one of the systems is not absolutely stable leads to the construction of the  equation.
• From the  and the  equations, LCO parameters of a practical piecewise linear system with multiple PNs (YF-12) are predicted along with the stability of the LCO.
Based on our approach, we are able to add more knowledge on how switchings between specific switching plains take place in the presence of LCO, and therefore understand how to adjust gains in order to either suppress or promote LCO phenomena. Furthermore, compared to other LCO stability analysis which barely checks the stability of an LCO, this research provides information on detailed dynamic modes that mostly contributes to the stability or instability of LCO. The following section will introduce the background necessary to describe it. They include Floquet theory and absolute stability of Lur'e type system. In the following section we present the main result of this work. Firstly, how to represent PNs in a piecewise linear system and how it can be incorporated into the Floquet analysis will be discussed. Then a framework is introduced to analyze the LCO of piecewise linear system based on the Floquet analysis and piecewise linear analysis. In the presence of multiple PNs, a more sophisticated approach will be demonstrated, especially in the analysis of a non-trivial switching order. Following section demonstrates this framework with an example of a realistic flight control system of the YF-12, which is represented by a piecewise linear system with multiple PNs. Finally, we present conclusive remarks along with suggestions for future studies.

Floquet theory
G. Floquet developed a classical result in the analysis of an LTP system in 1883 and this result still remains one of the most important theorems in the area of LTP systems [18]. The essence of Floquet theory is that the fundamental matrix of an LTP system (P(t )) is decomposed as a product of a timeperiodic matrix and a matrix exponential, of which the eigenvalues reveal the characteristics of the periodic solution.
Theorem 1. [19] In an LTP system whose state space representation isẋ(t ) = A(t )x(t ) with the period of the coefficient matrix A(t ) T , the fundamental matrix P(t ) of this system is decomposed as where H(t ) is a T-periodic matrix. Furthermore, the following holds for the state transition matrix of this system (t, t 0 ).
Proof. (for the proof of (1), refer to the proof of Theorem 3.2.1 in [19]) □ The proof of (2) comes from the basic properties of the state transition matrix. The state transition matrix (T, 0) in (2) is called as the monodromy matrix of the LTP systemẋ(t ) = A(t )x(t ).

Definition 1.
Characteristic multiplier [20] The characteristic multipliers of an LTP system are the eigenvalues of the monodromy matrix of the LTP system. The next theorem tells us about the eigenvalues of the monodromy matrix in relation to the properties of the periodic solution of nonlinear systems. Consider an autonomous nonlinear systemẋ where f : ℝ n ↦ ℝ n , f ∈ C 1 and the autonomous system above has a periodic solution of the least period T . Now we linearize the system above around its periodic solution (t ).̇x where is a Jacobian matrix of the system of (3) evaluated at the periodic solution (t ). Since (t ) is a periodic solution with period T, D(t ) is a T-periodic matrix thus the system of (4) is an LTP system. Then the following holds for the LTP system of (4). Theorem 2. [20] One of the characteristic multipliers of the LTP system of (4) is unity. In addition, if all of the rest of the characteristic multipliers reside within the unit circle in the complex plane, the periodic solution (t ) of the system of (3) is a stable limit cycle. Otherwise, (t ) is an unstable limit cycle.
Proof. (Refer to the proof of Lemma 10.2 and Theorem 11.1 of [20]) □

Absolute stability of Lur'e type system
The type of system that includes piecewise linear systems of our interest is called a Lur'e type system and is defined as below [21]. Figure 1. m is the number of system inputs. The nonlinear input function i is in sector [0, k i ], that is, it is within the sector bounded by y i axis and the line of slope k i > 0, as shown in Figure 2. Therefore, a PN is a special type of such s. In short, The notion of absolute stability is introduced to describe the stability of Lur'e type systems.

Definition 2.
Absolute stability [22,23] The system of (5) is said to be absolutely stable within the sector [0, k i ](i = 1, … , m) if the system is globally asymptotically stable about its equilibrium x e = 0.
We can check the absolute stability of the system (5) with Popov's criterion as below: Popov criterion [24] System (5) in which every nonlinearity is in sector [0, asymptotically stable (necessary condition) and if there exists a positive real number such that In particular, the condition of (7) is equivalent to the statement "the transfer function (1 + i )G (i ) + 1 k is positive real (PR)". In case there are multiple (m) PNs, (7) becomes where A, B, C, D, P, L, and W are of appropriate dimensions.
In other words, if we can find a positive definite matrix P that satisfies (9), the absolute stability of the system 5 is guaranteed by the KYP lemma and the Popov criterion.

MAIN RESULTS
The basic framework to analyze the LCO parameters-LCO period (frequency) and LCO amplitude is as follows. Firstly, we express the system both in state space representation and a cascaded series of linear transfer functions and nonlinearities. The later form, with the transformation of piecewise nonlinearities into corresponding equivalent analytic functions (EAF), yields the LTP by which to extract the  equation based on Floquet theory. Secondly, select the approximate range of LCO parameters with numerical integration and DF analysis. Both of them have limitations in refining the results into accurate LCO parameters, but it is important to narrow down the search space for the true LCO parameters. This significantly speeds up the process because our technique is based on Floquet theory and piecewise linear analysis with numerical calculations. Then we establish the necessary conditions on the LCO parameters.
Since there are two unknowns-LCO period and LCO amplitude, we also need two algebraic equations, one of which we obtain from Floquet analysis and the other from a switching equation derived from the piecewise linear analysis. We call the equation built on Floquet analysis the  equation, and another one from the switching equation the  equation, respectively. However, the two algebraic equations are not necessarily solved in closed form, so we need to iterate those equations in a certain range of the unknowns in order to find the answer numerically. That is why we needed the range of the two parameters with the preliminary work. The detailed procedure is described in Table 1.

Nonlinear EAFs for PNs
Since PNs usually consist of multiple sections, they are potentially discontinuous. Therefore, linear systems with PNs are, in -Incorporate b) and d).
h) Obtain the approximate range of T and A from (f) and (g).

STEP3
i) Assume an input signal to a specific PN. fact, also a Lur'e type system as in the equation below [25].
The sector property of a PN plays a critical role in preventing Floquet theory from being used in the analysis of LCO for piecewise linear systems, because PNs are nondifferentiable. Fortunately, we can build the EAF that enables us to make corresponding PN differentiable. The construction of EAF starts from an observation that g(u) in the following expression reproduces a signal identical to a simple relay. (11) where r (u) is a simple relay. Hence we define the EAF for a simple relay as follows.
Note that as grows larger, g(u) becomes closer to r (u). Thus, we claim that g(u) converges to r (u) point-wisely, but not uniformly.
Theorem 3. Convergence of g(u) to r (u) in (11) g → r pointwise on u ∈ (−∞, ∞), while g ↛ r uniformly on u ∈ (−∞, ∞) where | ⋅ |  is Lebesgue measure and E is a set whose (Lebesgue) measure is zero [26]. Now define a set F as then there exists a natural number n such that F ∉ E. In this case, Thus it is a contradiction to (12). □ Practically, with finitely large enough g(u) becomes close enough to r (u). Next, we check the differentiability of this EAF for a simple relay. By the appearance of a simple relay, the most dubious point is where u = 0. Checking the differentiability at u = 0 leads us to dg du , and with respect to t (time) as well, only if u = u(t ) is differentiable in time. Similarly, we can also express other common PNs with corresponding EAF. Observing that the input derivatives of most PNs look like a combination of reversed relays as shown in Figure 3, we can construct EAFs corresponding to PNs that are common in control systems.

Initial search on LCO parameters
For preliminary search of LCO parameters-LCO frequency and LCO amplitude, we utilize describing function and numerical integration.

Describing function analysis
Consider a single-feedback linear system with multiple nonlinearities shown in Figure 4. In this system, G 1 (s) and G 2 (s) are linear transfer functions and N 1 (A, ) and N 2 (A, ) are DFs of corresponding PNs. We now utilize the singularity of 'the system state flow' matrix [3].
by the singularity condition, we obtain However, this argument is based on the assumption that PNs are transformed into a completely linear component. Therefore this could lead to the violation of either of the two assumptions on which the validity DF analysis holds [3]. Nevertheless, DF analysis can still provide an estimation of the LCO parameters with an accuracy reasonable enough for a preliminary search.

Numerical integration
Numerical integration is the most straightforward and convenient way to track the propagation of the state trajectory of even quite a complicated differential equation of motion. The same argument is applied to the analysis of LCO because propagation of the state for a sufficient amount of time may converge to a periodic solution. However, numerical integration has various sources of error such as, machine error, mathematical truncation error, noise in function evaluation, under or overflow errors in calculation, and so on [27]. Furthermore, in the course of integrating in forward direction-that is, in a direction in which time increases as integrating steps go by, depending on the property of Jacobian along the integration trajectory, errors stemming from those sources can positively propagate only to deteriorate the reliability of the numerical integration [28,29]. Notwithstanding, we may still utilize the numerical integration as a useful tool to approximately locate the LCO parameters, together with DF analysis.

Floquet analysis: The  equation
Since the exact analytic expression of piecewise nonlinearities are available with corresponding EAFs, we can model a piecewise linear system as a differentiable nonlinear system, such as (3). Since the existence of LCO suggests that the piecewise linear system has a periodic solution, Floquet theory is applicable. By incorporating the LCO parameters in question into the expression of the characteristic multipliers of the system, we extract a necessary condition through which to establish the  equation. As described in Section 2.1, Theorem 2 implies that a periodic solution to the system guarantees one of the characteristic multipliers is unity (one). This characteristic multiplier is very important for the construction of the  equation.

Definition 3. Primary characteristic multiplier (PCM)
The primary characteristic multiplier is a characteristic multiplier of an LTP system of which the real part is the greatest among all of the LTP's characteristic multipliers.
The calculation of PCM is important especially for the system with a stable LCO because, according to Theorem 2, PCM should necessarily be equal to unity for a stable LCO. The next step is to linearize the system around the hypothetical periodic solution (t ) to obtain corresponding LTP system.
where x ∈ ℝ n is small perturbation around the periodic solu- is periodic as well. Assuming that the period of (t ) is T , (15) is an LTP system. Now we need to obtain the expression for the characteristic multipliers of the LTP system of (15). Before that, however, the expression for (t ) should be determined because (t ) appears explicitly in D(t ). Considering that the input to the EAFs usually are the output of the linear part of (10), we assure once again the input signal to each PN is a unimodal sinusoid signal. Therefore, the periodic solution (t ) constitutes this input signal (t ) ∈ ℝ m , where m is the number of PNs in the system. where where A p s and p s are the input amplitude and the input phase angle for the pth PN, respectively, and without loss of generality, we may put 1 = 0. The rest of A p s and p s are determined through the reasoning based on the transfer functions between each p (t )s. In other words, if we define G pq (i ) (p ≠ q, p, q = 1, … , m) as the transfer function between the input to the pth PN and the input to the qth PN, then It follows that A q = A q ( ) and q = q ( ). In some configurations, however, this G pq (i ) involves describing function. For example, if linear part and PN appears in series, is a function of both A p and . The introduction of describing function into G pq is possible only on the assumption that the open loop bandwidth of transfer function L q (i ) is small compared to current frequency in concern-. Otherwise, we need to construct a composite describing function that binds N p (A p , ), L q (i ), and N q (A q , ) [3]. Figure 5 illustrates how G pq is obtained in each configuration.
To obtain the first necessary condition-the  equation, we firstly need to obtain the expression for the monodromy matrix (T, 0). Since the closed form solution of the fundamental matrix does not generally exist for linear time varying systems including LTP system, we need to turn to a numerical approach [30]. With sufficient number of small intervals that divide the overall time interval, we can approximate the state transition matrix of the system [31].
where N is the number of intervals, Δ = T N , and t i = iΔ. Now applying (19) we can calculate the characteristic multipliers by obtaining the eigenvalue of (T, 0). From the arguments in (5) and (16), the equation of motion of piecewise linear system iṡ Therefore, in the presence of an LCO with its expression for the periodic solution (t ), the Jacobian matrix D(t ) in (15) becomes where ′ (H (t )) ∈ ℝ m×n is a spatial derivative of (Hx(t )) with respect to Hx(t ) evaluated at x(t ) = (t ). Since (t ) contains T and A as discussed in (17) and (18), D(t ) also contains T and A. It also follows that the monodromy matrix (T, 0) is also a function of T and A from (19). The requirement that the PCM should be equal to a unity and the fact that (T, 0) is a function of T and A lead to By definition PCM has the greatest real part of all the characteristic multipliers, therefore making PCM unity guarantees stable LCO, if one exists. This leads to the  equation where F (T, A) stands for PCM as a function of the expected LCO period (T = 2 ) and LCO amplitude (A). Finally, we check the stability of the LCO by inspecting the locations of all the characteristic multipliers in the complex plain. Recent studies have shown the stability of LTP with matrix polynomial approach which works well, but the procedure is much simpler and intuitive based on Floquet theory [32]. All we have to do is to see if they are all within a unit circle in a complex plane. If they are we may say that the corresponding LCO is a stable one, and otherwise the LCO is unstable.

Switching equation: The  Equation
Besides the  equation, we need one more necessary condition to predict the LCO parameters-the  equation, which is build upon a switching equation.

Switching feasiblity
Depending on A p s with p = 1, 2, … , m, switching in a certain region or direction may or may not happen in the presence of LCO. The possible switchings in one PN, in combination with the possible switchings in other PNs form a finite number of switching configurations, and therefore each configuration has its own state space representation. If we succeed in proving that the system of a certain switching configuration is absolutely stable, we may say that there is no LCO in this system. This is because, according to the definition of the absolute stability of Lur'e type system in Definition 2, the absolute stability implies the global asymptotic stability towards the zero equilibrium, while the existence of an LCO implies the existence of other invariant set than the stable equilibrium. To this end, we need to i. Find the zero equilibrium of the system representing each switching configuration and check its stability ii. Try to show the transfer function M(i ) in (8) of the system is positive real (PR), or find a solution to (9) (positive definite matrix P) We continue eliminating one configuration by another until we have switching configurations that we cannot prove their absolute stability. Proving a certain system is not absolutely stable does not guarantee the existence of an LCO, but this at least tells us that we cannot prove the non-existence of an LCO, either.

Switching order
Given some switching configurations that are not absolutely stable, we develop how to find the switching orders depending on the pair of expected LCO parameters-T LCO and A LCO . We start from defining switching times for each PN. Firstly, we define the linear or affine section of each PN-ranging from p1 to pn p for p-th PN. Then, considering that y p s are all in a unimodal sinusoidal shape, we may index each switching time in each PN. For example, if the pth PN is one shown in Figure 6, the sequence of switching times are t p3,p4 , t p4,p5 , t p5,p4 , t p4,p3 , t p3,p2 , t p2,p1 , t p1,p2 , t p2,p3 . In other words, t pi,p j implies a switching time at which the linear affine region i is switched to the linear affine region j of pth PN. Assuming a unimodal sinusoidal signal as potential LCO, it is straightforward to express each switching time in terms of LCO amplitude (A) and LCO frequency(period) (T ).
After identifying the switching times in all of the PN of the system, we arrange those switching times in an ascending order. One remark is that, the order needs to be determined only after we take the phase differences between each y p into account, because there usually are nonzero phase difference between the input to each PNs, as discussed in (18).

Switching equation
After identifying the switching order we sort all the switching times up to n s over the range of T and assign indices from 1 FIGURE 7 Definitions of switching time, corresponding linear affine system index, and state variables(x( i )). Note that x(0) = x(T ) and system 0 =system n s to n s , that is, the phase-driven index, in ascending order-1 , 2 , … , n s , as in Figure 7. Also, we define each linear affine system corresponding to each switching time interval based on the combination of the output of each PN. The ith linear affine system is defined as one that governs the system dynamics during Now we are ready to derive the switching equation. State transition from system i to system i+1 is formulated by the standard solution of linear time invariant systems.
in other words, Given this, the following claim will lead to the switching equation.
Claim 1. Recurrence equation: , n > 2, Proof. This claim is proven by a mathematical induction. The cases of n = 1 and n = 2 are straightforward. Then try to propagate from n to n + 1 in the equation above, then we obtain .
Since A 0 = A n s and x 0 = x T , the switching equation is finally derived.
The  equation is obtained from the reasoning that without loss of generality, one element(usually the first one) of the output vector y = Hx in (20) is put to zero. This also corresponds to the assumption we make in (18) that the with x 0 replaced by 0 . Since each switching time interval Δt p in (26) is a function of LCO parameters as shown in (23), 0 is also expressed in the LCO parameters T and A. Therefore, (27) is a scalar equation with T and A, namely the  equation.

YF-12 flight control system
YF-12 is a former model of the reconnaissance plane SR-71 back in 1970s and the flight test program was run by NASA [10]. The pilot of the program reported a pilot induced oscillation (PIO) depending on the fuel storage condition allegedly resulting from the aeroelasticity of airframe. This PIO could essentially be the LCO that we are trying to analyze. The system diagram is in Figure 8 and the detailed specification of the system is in (29) [10].
where P (s) is a pilot model, A(s) an elevator actuator model, G p (s) and G h (s) transfer functions of the airframe, and H (s) a lead-lag compensator of pitch rate feedback. Based on the transfer functions, the state space model is established as follows [10].
where (Hx(t )) = where y i = H i x, i = 1, 2 and PNs g 1 and g 2 are depicted in Figure 9. Their corresponding EAFs and input-derivatives are as follows.
, g 2 (y 2 ) = 1 2 [ where d 5 = 5, d 25 = 25, and d 2.5 = 2.5. As shown in Figure 8, inputs of both the PNs are all coming after passing through a sufficient number of linear systems. Therefore we assume those inputs to be unimodal sinusoidal signals.
where A 1 and A 2 are the amplitudes of y 1 and y 2 , respectively, and 21 is the phase difference between y 2 and y 1 . The relationship between the parameters of y 1 and y 2 is clarified by the transfer function reasoning as below.

DF analysis
This DF analysis is based on the observation that every signal in the loop, that is, y 1 , u 1 , y 2 , and u 2 are all nonzero in the presence of an LCO (Figure 8). where N 1 and N 2 are the DFs of g 1 and g 2 , respectively. Thus, Since z in the equation above is a non-trivial solution in the presence of an LCO, det (N) = 0. This leads to (37) From Equations (37) and (35), we see that there are three equations (amplitude and phase condition from (37) and amplitude condition from (35)

Numerical integration
With an appropriate initial condition, the YF-12 system shows a stable circulating behavior as shown in Figure 10 and we observe that T LCO ≈ 1.318 (sec) and A LCO ≈ 45.426 (deg).

Floquet analysis
Assuming a periodic solution ( (t )) to this system, we linearize the equation of motion in (30). where Then we obtain the  equation based on the argument that PCM of the monodromy matrix is equal to unity, because the magnitudes of other characteristic multipliers should be within a unit circle (sphere) to guarantee the stability of the LCO. Since the input derivatives h 1 and h 2 contain the sinusoidal expression with LCO parameters as in the (33), (T, 0) is a function of A and T . Thus (33), (34), and (35) lead to the  equation The next step is to calculate the  equation through all the possible range of T LCO and A LCO . Based on the approximation of LCO parameters in the previous section we have narrowed down their feasible range. If we fix the range as from 1.2 to 1.5 for T LCO and 44 through to 52 for A LCO , then we obtain the  equation as depicted in Figure 11(a). Now we check the stability of the potential LCOs on the PCM curve in Figure 11(a).
To this end, we display all the characteristic multipliers corresponding to each point in the PCM cureve. As depicted in Figure 11(b), since all the characteristic multipliers are within a unit cirvle, we conclude that the LCO that we are going to identify is a stable one.

Switching feasiblity
All the possible switching configurations for YF-12 flight control system are shown in Table 2, each of which has its own state space representation.
Taking an example of configuration two in Table 2, we havė x = F 2 x + G 2 (y), y = Hx, where F is from (30). If we succeed in proving that the system (41) is absolutely stable, we may say that there is no LCO in this system.
Claim 2. The Lur'e type system corresponding to the switching configuration two in Table 2 is absolutely stable.
Proof. To find the equilibrium, we solve the following algebraic equation Then the equilibrium set is We now check the stability of this equilibrium. The Jacobian matrix around this equilibrium is Since y 2 (x e ) = 0, J(x e ) = F and its eigenvalues are −4, −10, 0, −0.785 ± 15.665i, −0.750 ± 1.854i, −34, the equilibrium is stable with only one mode marginally stable.
Since we have nonzero equilibrium x e , we need to check the absolute stability of the following system with zero equilibrium, because proving the absolute stability based on PR or KYP lemma implies the absolute stability towards the zero equilibrium.
or, equivalentlẏ Define 2 (H 2 x) ≡ 2 (H 2 x − H 2 x e ) + 2 (H 2 x e ) then 2 (H 2 x) = 2 (H 2 x) because H 2 x e = 0. Therefore, the equivalent zero equilibrium system is identical to the original system of configuration two. Now we need to check if this configuration is absolutely stable. We try to find the matrix P in (9), which is equivalent to the following LMI problem [33].
Therefore, by KYP lemma, the system with following transfer function M(i ) is PR.
Hence, by the Popov criterion, the system of the switching configuration two is absolutely stable. □ Likewise, we can check the switching feasibility for each configuration in Table 2. Each of them is proven to be absolutely stable except the last configuration, which cannot be proven to be PR. One of the necessary conditions to check the PR of a system is to check if the Nyquist plots of the transfer function M(i ) of (8) of the system are entirely in the right half complex plane for some 1 > and 2 > 0 [17]. In other words, if we prove that for every 1  where M(i, j )(s) is the ith row and j th column entry of M(s) and N i j (s), D i j (s) are corresponding nominator and denominator, respectively. Then we want to confirm if the real part of at least one entry of M(i, j )(i ) always has negative real part for some . It is simple to show that M(2, 1)(i ) always has negative real part when ∈ [1.4608, 3.2816]. Therefore, the Lur'e type system corresponding to configuration six is not PR, and not absolutely stable.

Switching order
We start from defining switching times for each PN. Define the linear or affine section of each PN-ranging from (1) to (5) for FIGURE 12 Linear or affine sections for g 1 (left) and g 2 (right)

FIGURE 13
Switching time for g 1 (left) and g 2 (right) g 1 and from (a) to (c) for g 2 , as shown in Figure 12. In addition, the output signal versus the input sinusoidal to g 1 and g 2 should look like Figure 13, where t i j stands for the switching time that transits from linear affine section i to j . i and j could be 1,2,3,4, or 5 for g 1 and a,b, or c for g 2 , respectively. The switching order in a single PN is straightforward-for g 1 , it is (3) ). Likewise, all of the twelve switching times-t 34 , t 45 , t 54 , t 43 , t 32 , t 21 , t 12 , t 23 , t bc , t cb , t ba , t ab are expressed in terms of T , A 1 , and 21 = 21 (T ), and are arranged in ascending order.

FIGURE 14
The frequency response of y 2 over y 1 over the expected LCO period range

Switching equation
Now we assign the twelve switching times the phase driven indices from 1 to 12 in ascending order-1 , 2 , … , 12 . Also, we define a few more terms based on the phase-driven index. For all i = 1, 2, … , 12, (Δt 0 = 1 , Δt 12 = T − 12 ), system i ≡ the linear affine system corresponding to Δt i .
The switching times segregate the whole system into a finite number of linear affine systems based on the combination of indexes. For example, if g 1 belongs to the section '1' and g 2 in section 'a' in Figure 12, respectively, the system index is '1a' as specified below.
Likewise, all linear affine systems are defined in Table 3. Once a switching order is established for a specific T and A pair extracted from the feasible ranges determined through Section 4.2, we can construct the switching equation, in such a way that each linear affine system governs the system dynamics during corresponding switching time interval Δt i .  Now we derive the switching equation using (26).
By (34), x 7 (0)) = 0 at phase zero. Since (44) is expressed in A and T , x 7 (0) is also a function of A and T . Therefore, this initial phase condition plays for another necessary condition, the  equation, along with (40).
The  equation is depicted in the T -A plane in Figure 15. A blurred line on the right side of the solid line in the right part of Figure 15 is due to the singularity of the surface, that is, the evaluation of x 7 (0) oscillates so radically that we cannot guarantee the continuity of the surface around this part. Still, the solid line is convincing and can be entrusted as the  equation.

Determination of the LCO parameters
Now that we have two necessary conditions ((40) and (45)) with two unknowns (T LCO and A LCO ), we may determine them. Figure 16 depicts the two necessary conditions in the T -A plane, and thus the intersection of the two solid curves should make the solution-T LCO =1.356 (sec) and A LCO =51.57 (deg). In Table 4 the results are summarized as to the LCO parameters for the YF-12 flight control system with different approaches. The characteristic multipliers at T LCO = 1.356 and A LCO = 51.57 1.798e−7, and 9.873e−18, respectively, hence we confirm this LCO is a stable one. Inspection on the characteristic multipliers indicates that the most dominant one other than the PCM is 0.4806 and its corresponding dynamic mode should involve corresponding eiginvector.

CONCLUSION
In this research a new framework is established to analyze the limit-cycle oscillation (LCO) of piecewise linear systems with multiple piecewise nonlinearities(PN). While PN in the piecewise linear system prevented the application of analytic methods such as Floquet theory to the analysis of LCO, the EAFs of corresponding PNs are developed to practically incorporate the PNs into the Floquet analysis. In addition, the switching function based on the exact switching order of linear affine systems that organize the piecewise linear system plays another necessary condition along with one obtained from the Floquet analysis. Finding the exact switching order is possible through proving each switching configuration is absolutely stable but only one, which possess a stable LCO. The stability of LCO is also analyzed with Floquet theory, by inspecting the location of characteristic multipliers of the perturbed linearized system evaluated at the expected periodic solution. Through this procedure, we gained independence on the open loop bandwidths of each subloop that each contains piecewise nonlinearity, thus this independence negated the inaccuracy shown in existing methods such as DF. The example of the YF-12 flight control system demonstrates the effectiveness of this framework. Future studies include, but are not limited to the followings. Firstly, we need to explore the multi-modal sinusoidal assumption as a periodic solution. Especially in presence of multiple feedback loops such as the YF-12 system, we have a high chance of having another nonnegligible harmonic components due to the different loop bandwidths. Another challenge will be finding the sufficient condition for the existence of LCO in piecewise linear systems. Furthermore, finding the region of attraction to LCO will be a substantial contribution to the effort to prevent adverse LCOs in the system. Parametric studies will also add perfection to the analysis of LCO properties, including robustness analysis with the system parameters such as pilot time constant.