Application of the Filippov Method to PV‐fed DC‐DC converters modeled as hybrid‐DAEs

In this article, a method to study the nonlinear dynamics of a PV‐fed boost converter is developed. The model for the solar panel introduces an algebraic constraint into the system and the resulting mathematical model is a set of hybrid differential algebraic equations (DAEs). The behavior of the system is observed to exhibit undesired operation as parameter values vary. Conventional mathematical tools developed to analyze DC‐DC converters are not directly applicable to systems modeled as a set of hybrid DAEs. In order to find the monodromy matrix to assess the eigenvalues of the system, we modify the Filippov method to calculate the saltation matrix. The derived formula is general and versatile such that it can be applied to any PV‐fed DC‐DC converter irrespective of the topology or control algorithm employed. Two case studies are investigated: current mode control and maximum power point tracking. Each case study serves to demonstrate different considerations that must be taken into account when conducting stability analysis of hybrid‐DAEs.


INTRODUCTION
A large-scale redesign process is currently under way in power systems across the globe. Conventional sources of energy based on fossil fuels are being replaced by renewable sources of energy such as photovoltaic (PV) solar panels, wind turbines, and batteries. This redesign is primarily motivated by growing environmental concerns as well as the economic aspects of nations maximizing their available energy resources. 1 PV solar panels have gained an increase in market share due to falling prices, the free availability of sunlight, and their low rate of emissions. 2,3 PV panels have a nonlinear I-V characteristic curve. Hence, as the solar irradiance varies, the output voltage of the PV panel will vary according to the I-V curve. The output voltage of the PV panel is nonlinear with a maximum power point (MPP) voltage in the range of 15 to 40 V. 4 The output voltage of each PV panel must be stepped up to meet the demands of the interconnecting device. Such applications include connecting to the AC grid through a DC-AC inverter, charging a storage battery or electric vehicle, or contributing to the 380 V DC distribution system. 5 The conventional method for stepping up the output voltage of multiple PV panels is to connect several PV panels together in series. As the combined output of the array of panels is variable, some form of power electronics is used to condition and regulate the output voltage. However, shading and other unideal operating conditions can drastically reduce the efficiency of series connected PV devices. 6 The approach under study in this article is to connect a DC-DC converter in series with each PV panel. Each PV module operates independently from one another and can be connected to the load, battery, AC-grid, or DC-link separately. The control algorithm for each DC-DC converter can be optimized as environmental conditions vary in order to operate at the MPP.
The desired behavior of DC-DC converters is a stable periodic motion around a predefined value with a frequency that is equal to the external clock. However, as parameters vary, such systems exhibit a wide range of nonlinear behaviors such as subharmonic oscillations. [7][8][9][10] These oscillations can manifest themselves through a series of bifurcations, which can increase the current/voltage ripple, add extra harmonics, and increase the switching losses. 10 Thus, bifurcations can greatly impact on the efficiency of the system. This bifurcation behavior can lead to high ripples in the state variables or control variables and high stress in the switches. For example, a period-doubling bifurcation can double the current ripple, which can damage the DC-DC converter, the PV source, or the load. 11 Hence, accurate techniques to predict and avoid their occurrence in PV applications is essential in order to maximize the efficiency of the overall system. A variety of different control algorithms for DC-DC converters have been proposed in the literature to avoid such undesirable behavior. [12][13][14][15] The typical method for predicting and avoiding fast-scale instabilities in switching circuits is through numerical analysis by discrete-time modeling 16,17 or Floquet theory. 18 Further details of the later approach can be found in Section 3. However, there is limited research into the prediction of these undesirable behaviors when the input is nonlinear and the system is modeled as a set of hybrid differential algebraic equations (DAEs) as opposed to a set of piecewise-linear ordinary differential equations (ODEs). This work aims to address this issue, that is, to develop a method to determine the stability of PV-fed DC-DC converters that are modeled as a set of hybrid DAEs. Herein lies the novel contribution of this work: • To develop an algorithmic approach that is mathematically accessible, to assess the stability of PV-fed DC-DC converters modeled as a set of hybrid DAEs. The derived formulas are versatile in that they can be applied to any PV-fed DC-DC converter topology directly. This can be found in Section 4.
• The derived algorithm is verified through brute-force computer based simulations in Section 5. The circuit studied in this work is also studied in Reference 11, where experimental results are presented but mathematical analysis is not performed. This research conducts that quantitative analysis by extending the Filippov method to enable its application to systems modeled as a set of hybrid DAEs. By considering the same system model with same parameters, the experimental results of Reference 11 can be related to this work.
• Two case studies are considered in Section 5. The first is chosen in order to demonstrate that the method for the calculation of the transition matrix for DAE modeled systems must be modified from that for ODE modeled systems in order to find the correct monodromy matrix.
• The second case study illustrates when the algebraic constraint is included in the switching condition that the formulas derived in this research must be applied in order to find the correct eigenvalues of the system.
• In Section 2, we compare the results observed and novelty of this body of research to that found in the literature.

LITERATURE SURVEY
The research relating to PV-fed DC-DC converters is sparse and limited. In this section, we provide a brief overview of the available literature and identify how it differs from the novelty and results of this research. There exists a variety of methods to predict bifurcations in DC-DC converters. These include trajectory sensitivity analysis, 19 an auxiliary vector method, 20 among others. However, they have not found widespread acceptance or application in power electronics design. This may be due to the mathematical complexity in their application. In recent years, the Filippov method has gained the most traction for studying the nonlinear dynamics of DC-DC converters most probably due to the mathematical simplicity of its application. However, little research has been conducted into quantitative analysis of PV-fed DC-DC converters or DC-DC converters modeled as a set of hybrid DAEs. Hence, the focus of this research is on the extension of the Filippov method to PV-fed DC-DC converters modeled as a set of hybrid DAEs. This is not the first article to study the nonlinear dynamics of PV-fed DC-DC converters, a new method for studying the nonlinear dynamics of DC-DC converters and predicting unstable modes of operation was proposed in References 21 and 22. In References 4,22,23, this new method is applied to a PV-fed boost converter and quadratic boost converter. However, the system differs slightly from that studied in this research as the output of the PV panel does not feed directly into the DC-DC converter. Instead, there is an input capacitor in order to smooth the input voltage and avoid the injection of high current ripples. This enables an additional state variable to be created and the system can be modeled as a hybrid ODE.
In this work, the PV panel is connected directly to the DC-DC converter and the DC-link capacitor is omitted. The resulting mathematical model is a hybrid DAE. The rationale for this can be found in Section 5. The results of References 4,22,23 are only applicable to hybrid-ODEs. Hence, the results of this work are more general as they can be applied to DAEs and ODEs with no adaptation to the equations to model the system.
The application of the Filippov method to PV-fed DC-DC converters is not a novel concept. In Reference 11, stability analysis of a PV-fed boost converter operating under different control algorithms is performed. The authors identify over what range of values of R, I ref , and I ph give a stable period-1 output through simulations and bifurcation diagrams. The results are experimentally verified. For peak current mode control (CMC) and average CMC, the eigenvalues of the system are found using the conventional Filippov method. As outlined in this research, the conventional method is applicable as the algebraic variable is not present in the switching manifold, that is, n T y = 0 for these control strategies. However, only qualitative results are presented when the maximum power point tracking (MPPT) control algorithm is considered. There is no quantitative analysis performed. The stability of the system is determined through bifurcation diagrams, time-domain simulations, and experimental waveforms. The eigenvalues of the system are not found. This is because the standard Filippov method cannot be employed as the algebraic constraint is present in the switching condition for MPPT control.
The approach taken in this research is to extend the work of References 11 and 18. Both of these works consider quantitative analysis of hybrid ODEs and the derived formulas only have a mathematical foundation in such systems. This work considers hybrid-DAEs and extends the mathematical foundation developed in Reference 18 to such systems. A new formula to calculate the saltation matrix is derived to include the effect of the algebraic term on the calculation of the stability of the system. This enables the Filippov method to be applied to hybrid DAEs, as well as, hybrid ODEs.
By considering the same circuit, parameters, and operating conditions, the experimental results of Reference 11 can be extended to verify the validity of proposed method developed in this work.
The main author in References 4,23 contributed to References 24-26 where they study the bifurcation behavior of a PV-fed boost converter. In a similar manner to Reference 11, the conventional Filippov method is applied in order to find the eigenvalues of the monodromy matrix. However, there is an input capacitor filtering the output of the PV-panel. Hence, the approach taken in References 24-26 is not as general as that developed in this research.
In Reference 27, an MPPT controlled boost converter fed by a PV-panel is studied. In all the articles discussed in this literature survey, the temperature and solar irradiance is fixed when studying the effect of other parameters on the nonlinear dynamics of the system or the solar irradiance is the bifurcation parameter with all other parameters fixed. However, in Reference 27, the authors calculate the maximum power point for each value of solar irradiance and fix v in at the MPP. For this reason, the results of Reference 27 are for a specific converter, at a specific temperature, and operating under specific conditions. The algorithm presented in this article is generalized and can be applied to any DC-DC converter with a PV source with any operating conditions.
In all the research articles discussed here, the system is modeled and analyzed as an ODE. Hence, it is important to consider the work carried out in the field of hybrid DAE analysis. The work of Federico Bizarri and Angelo Brambilla in Reference 28 studies mixed analog/digital circuits that are modeled as hybrid DAEs. In a similar manner that is undertaken in this research, they derive a generalized version of the saltation matrix that is applicable to hybrid DAEs. However, the systems that they study undergo an impact event, that is, G 1 ≠ G 2 unlike this work where G 1 = G 2 . Furthermore, the application of the derived saltation matrix is to analog mixed signal circuits and not PV-fed DC-DC converters or similar.
The proposed algorithmic approach offers a powerful technique that can be applied to any DC-DC converter modeled as a set of hybrid DAEs, a set of piecewise-linear ODEs, or a set of hybrid ODEs. The approach can be applied regardless of the index of the DAE, the number of algebraic constraints present in the DAE, or the topology of the DC-DC converter. The boost converter is studied in this work as experimental results are presented in Reference 11. Hence, the same system model and parameters are used. The approach is valid if isolated DC-DC converters are considered, such as the forward converter or other transformerless high-gain DC-DC converters such as the quadratic boost converter.

OVERVIEW OF FILIPPOV THEORY
In this section, a brief overview of the Filippov method is presented. The differences between studying a hybrid-DAE and piecewise-linear ODE are highlighted along with the relevant notation used throughout this work. The desirable behavior of a DC-DC converter is a stable periodic motion around a predefined value with a frequency that is equal to the external clock. As parameters vary, the system may undergo a bifurcation and the stability of the periodic motion may be lost. This can greatly increase the current ripple and therefore damage the DC-DC converter itself, the PV-source, or the interconnecting load. 11 Operating in the unstable region may lead to a degradation in the efficiency of PV array or reduce the lifespan of the PV module. Hence, the prediction of these stability bounds is imperative for PV arrays.
One popular method to analytically determine the stability bounds of DC-DC converters is the Filippov method. The stability of a general orbit, say x(t), is assessed by placing a small perturbation at t = t 0 and monitoring the evolution of the perturbation Δx(t). The evolution is related to the initial perturbation by the fundamental solution matrix. When the vector field that governs the original orbit is linear time invariant (LTI), the fundamental solution matrix is given by the exponential matrix: where Φ(t, t 0 ) is termed the state transition matrix. In this work, the system is modeled as a nonlinear hybrid DAE. Hence, the fundamental solution matrix is not termed the state transition matrix as it refers to the evolution of both the perturbed state variables and the perturbed algebraic variables. Instead, we refer to it as the transition matrix, Γ(t, t 0 ). The transition matrix contains two submatrices and is of the following form: where N D is the number of differential variables in the system and N A is the number of algebraic variables. We will still refer to Φ(t, t 0 ) as the state transition matrix as it describes the evolution of the state variables with respect to the state variables. We will refer to Ψ(t, t 0 ) as the algebraic transition matrix as it relates the evolution of the state variables to the algebraic variables. Hence, if N A = 0, then Γ(t, t 0 ) = Φ(t, t 0 ). Furthermore, the transition matrix is not given by the exponential matrix as the system is not LTI. Instead, it must be calculated numerically. Section 4.2 provides a mathematical reasoning as to why Γ(t, t 0 ) is constructed in this form and presents a method for numerically calculating Φ(t, t 0 ), Ψ(t, t 0 ), and Γ(t, t 0 ).
If the orbit is periodic, the stability can be quantitatively determined by assessing the eigenvalues of the fundamental solution matrix evaluated at t = t 0 + T, where T is period of the orbit under study. The fundamental solution matrix obtained at t = T is termed the monodromy matrix.
DC-DC converters switch between two or more vector fields. As a result, the stability of (1) cannot be assessed directly as the switching action must be taken into account. The saltation matrix, S, is a form of the state transition matrix that relates the perturbation before and after switching.
where h(x,t) = 0 represents the switching manifold, n is the vector normal to the switching manifold and n T its transpose, and f represents the right-hand side of the differential equations, and the subscript notation 1 and 2 indicate if the function is evaluated before or after switching. However, (3) is the equation to find the saltation matrix for a piecewise-linear ODE.
In Section 4.1, we will present a reformulated version of (3) applicable to hybrid DAEs. Since we are interested in the stability of a periodic oscillation exhibited by converters, we need to calculate the transition matrix over an entire clock cycle. DC-DC converters are switched-mode power supplies that operate by opening and closing switches in order to step up/down the input voltage. In this work, we analyze the DC-DC step-up boost converter fed by a PV panel. This has one switch that opens and closes depending on the value of the control signal. If we assume that the switch is closed (ON) at the beginning of the clock cycle from t = 0 to t = dT (where d is the duty cycle) and then opens (OFF) when some switching condition is met, h(x,t) = 0, and remains open from t = dT to t = T. The monodromy matrix is given by: S 2 relates to the resetting event at the end of the switching cycle. Ramp compensation is typically used for DC-DC converters and the falling edge of this signal leads to S 2 being the identity matrix. For this reason, S 2 is omitted from future discussion in this work. Hence, (4) becomes: When a PV-fed DC-DC converter is considered, the system is modeled as a nonlinear hybrid DAE. As a result, the formula for the calculation of the saltation matrix must be adapted accordingly. This adaptation is presented in the next section.

MONODROMY MATRIX FOR DAES
In this section we will outline a method to calculate the monodromy matrix for a hybrid-DAE. First we will derive the method used to calculate the saltation matrix before briefly describing the method to find the transition matrices.

Saltation matrix for DAEs
The PV-fed DC-DC converter is modeled as a hybrid DAE of the form: where x are the differential state variables and y the algebraic variables, and the dot notation indicates the derivative with respect to time. It is assumed that f and g are differentiable and their partial derivative matrices are referred to as f x , f y , g x , g y , and g t . The initial condition at t = t 0 for the state variables is given by x 0 . The initial condition for the algebraic variables, y 0 , must be chosen in order to satisfy the constraint 0 = g(x 0 ,y 0 ,t 0 ). We can reformulate the DAE as an equivalent ODE by writing: For brevity, the dependence on (x,y,t) is omitted. Letting G = −g −1 y (g xẋ + g t ), the equivalent ODE is given by: Let S ode be the saltation matrix for the reformulated ODE given in (6): The saltation matrix, S ode , relates the differential and algebraic variables before and after switching: However, instabilities manifest themselves through the differential variables. The algebraic variables are a constraint on the system. Hence, we are not interested in the relationship between Δx 2 and Δy 1 . However, due to the algebraic constraints in (6), a perturbation of Δx 1 to x 1 leads to an perturbation in y 1 of −g −1 y g x Δx 1 . Hence, [ Δx 2 In order to derive an expression for the saltation matrix that is directly applicable to a hybrid DAE, let us consider each component of (7) separately. First, consider the bottom line of the fraction: As n T is the derivative of the switching manifold with respect to the state and algebraic variables, we can define n T x as the first 1 × N D entries of n T and n T y as the remaining 1 × N A entries. Hence: we can write: Rewrite (8) as: where S D relates the perturbation of the state variables before and after switching and S A captures the effect on the algebraic variables. Hence, where .
If both matrix elements of S D are postmultiplied by , we get the following: Hence, we can rewrite (9) as: Now consider: where S A is given by: y 1 ,t 1 ) .
However, G 1,2 = −g −1 y g x f 1,2 . If we postmultiply by , we get the following: Hence, we can rewrite (11) as: Hence, we can modify the expression of S ode in (7) and, as the perturbation of the algebraic variables do not have an impact on stability, rewrite (8) to derive an expression for the saltation matrix that is directly applicable to a DAE: where S Φ and S Ψ are given by (10) and (12), respectively. The resulting equations derived in this section, (10) and (12), can be directly applied to the hybrid DAE system model. Herein lies the novelty of this work. The system model does not need to be reformulated as an ODE. The derived equations take the reformulation of the DAE to an ODE about the switching point into account to enable the direct application of the adapted Filippov method to the hybrid DAE. In order to determine the stability of a periodic orbit the saltation matrices and transition matrices must be calculated. In the next section, we present an adapted method from Reference 29 to find the transition matrix of a hybrid DAE.

Calculation of the transition matrix
Consider an ODE with the initial condition x = x 0 and in the following form: where u is the input to the system and A and B are N D × N D and N D × 1 matrices, respectively. If we apply a perturbation at t = t 0 , the distance between the perturbed orbit and the original orbit will evolve in time satisfying the variational equation: 30Φ where P = J for an ODE and J is the Jacobian matrix of the ODE. If the system is LTI, Φ(x 0 , t) = e A(t−t 0 ) . Otherwise, Φ must be calculated numerically. Hence, we solve both (15) and (16) in parallel so as to obtain the solution of the state variables, which can be used to form the Jacobian matrix at each time step as the Jacobian matrix is now time-dependent. In this research, the system is modeled as a nonlinear hybrid DAE. For this reason, (16) cannot be solved by setting P = J. Consider a linearized form of the DAE: If we assume that g y is nonsingular, we can write y = −g −1 y g x x. Hence: Setting P = (f x − f y g −1 y g x ) allows (16) to be solved for a DAE. Hence: In general, (18) can be solved in parallel with (17) using any numerical integration method. The sensitivity of y with respect to x 0 is given by: The formulation of the transition matrix is given by: .
It is important to note that the calculation of (18) is computationally expensive as finding g −1 y at every time step is not trivial for large systems. A computationally efficient method for estimating both Φ and Ψ can be found in Reference 30. However, the focus of this research is to develop an algorithm for a PV-fed DC-DC converter, which typically has very few differential and algebraic variables. Hence, this research is not concerned with optimizing the algorithm for computational efficiency. Instead, we are concerned with reducing the algebraic complexity of the proposed approach and, for this reason, (18) and (19) are used to find Γ in this work.
We will now illustrate the adapted method by finding the eigenvalues of a PV-fed boost converter.

PV-FED BOOST CONVERTER
The output voltage of a PV panel is in the range of 15 to 40 V. However, applications at grid level require a much higher voltage, such as the 380 V DC-link voltage. Hence, the output voltage of the PV panel must be stepped up. For some applications, the required high-voltage gain could be as high as 20. 4 Conventional DC-DC converter topologies, such as the boost converter, cannot achieve this due to parasitic resistances in the switching devices and reactive components. 4 However, the purpose and novelty of this article lies in the development of an algorithm to determine the monodromy matrix for PV solar panels operating with DC-DC converters in series where the system is modeled as a set of hybrid DAEs. Hence, the topology of the converter is not important as all relevant DC-DC converters can be modeled as a set of ODEs. The PV panel supplying the input to the boost converter is what introduces the algebraic term and hence is the reason we model the system as a hybrid DAE. For this reason, for simplicity but without loss of generality, the boost converter is chosen.
In this section, we will introduce the mathematical model used to simulate the boost converter and PV module. We will then consider two control algorithms for the PV-fed boost converter: 1. Average CMC: this control algorithm is chosen because n T y = 0. This case study will demonstrate that while the saltation matrix widely found in literature for ODEs can be applied to DAEs when n T y = 0, the conventional methods for finding the transition matrix cannot. Hence, the system cannot be strictly treated as an ODE. 2. MPPT: this control strategy is chosen as n T y ≠ 0. Hence, the adapted saltation matrix given in (14) must be applied to the system. Otherwise, the calculation of the eigenvalues is wrong and may lead to incorrect design or overdesigned components.
All simulations are obtained using Matlab 2017b using the ode15s solver with the tolerance set to 1 × 10 −10 with a maximum step size, Δt, of T/100, where T is the switching period of the boost converter. All simulations were executed on a 64-bit Windows 10 operating system running on an 4 core 1.90 GHz Intel Core with 16 GB of RAM.
It is important to note that there is no experimental verification presented in this research. However, in Reference 11, a PV-fed boost converter is considered. The results of Reference 11 are verified experimentally. For this reason, we have chosen the same parameter values for the circuit and control system. Hence, by extension, the experimental results of Reference 11 verify the approach undertaken in the proceeding analysis.

Mathematical models for simulation
The circuit diagram for the PV-fed boost converter, using the single diode model, 31 is shown in Figure 1. The diode equation is: where A = q∕( kT e ) and I 0 is the saturation current of the diode, q is the charge of an electron =1.6 × 10 −19 C, k is the Boltzmann's constant =1.38 × 10 −23 J/K, T e is the absolute temperature, and is the diode ideality factor. Applying Kirchhoff's current law: Since v pv = v in + i L R s and substituting in (20), we can derive: PV devices have a nonlinear I-V characteristic. 31 Hence, the resulting equation to find the input voltage to the boost converter, v in , is nonlinear. Equation (21) is an algebraic constraint that causes the overall system to be modeled as a hybrid-DAE.
It is important to note that in Figure 1, and the resulting mathematical model, there is no DC-link capacitor between the PV panel and the boost converter. A real-world system would include this capacitor and the resulting system would be modeled as a set of hybrid ODEs as in References 24-26 and not a set of hybrid DAEs. The reason for not including the DC-link capacitor is 2-fold: 1. Power systems and microgrids are often modeled as a set of DAEs. 32 The inclusion of a switched mode power supply changes the system model to a set of hybrid-DAEs. An example of this is Reference 33 where a PV panel is connected F I G U R E 1 PV-fed boost converter to a battery. The battery is modeled as a DAE and hence the resulting system is a set of hybrid-DAEs. There are a variety of different applications and loads where the resulting mathematical model for the entire system is a set of hybrid-DAEs. Hence, for simplicity without loss of generality, an resistive load is considered and the DC-link capacitor removed in order to model the system as a set of hybrid-DAEs. 2. The work carried out in this article extends that of Reference 11. The experimental work in Reference 11 serves to verify the proposed approach, where the DC-link capacitor is not included. Hence, it is not included in this work.
For the boost converter, the switch is closed at the start of the switching period (ON) for a time dT and then opens and remains open (OFF) for the rest of the switching period, that is, (1 − d)T. The control signal, v con , is compared against a ramp signal, v ramp . When v con = v ramp , the converter switches from ON to OFF. Hence, the switching manifold is given by h(x,y,t) = v con − v ramp = 0, which determines when switching occurs. The state equations for the boost converter are expressed as:[ where Hence, combining (21) and (22), the equations to model the system illustrated in Figure 1 gives a nonlinear hybrid DAE:⎡ Hence, the differential variables, x, are i L and v o and the algebraic variable, y is v in . Unless otherwise stated, the parameters used in this work are as follows (Table 1):

Average current mode control
CMC algorithms for switching power devices can give rise to many issues such as poor noise immunity and peak-to-average current errors which the low current loop gain cannot correct. 34 Average CMC is a popular method employed to overcome the limitations of other current sensing techniques. In average CMC, the inductor current is compared against a reference current to create an error signal. The error signal is then adjusted by a proportional integral (PI) control algorithm to produce a control voltage: The control signal is compared against a ramp voltage: where V L is the lower value of the ramp signal and V U is the upper value. Using this control strategy, at the start of the switching period the switch is open. When the ramp signal is equal to the control signal, the switch is closed and remains closed for the rest of the switching period. Figure 2 shows a bifurcation diagram with R as the bifurcation parameter ranging from 40 to 46 Ω. The system is simulated for 100 switching cycles for each value of R and the output sampled at the end of each switching cycle. It is clear, that for values of R ≤ 43.5 Ω, the system operates with a stable period-1 orbit. When R = 44 Ω, the system has exceeded the critical load value and undergoes a period-doubling bifurcation. The system operates in this unstable region for all values of R ≥ 44 Ω.
Bifurcation diagrams are qualitative tools for assessing the stability of systems. The Filippov method is a quantitative tool that can achieve the same results. As the system is a nonlinear hybrid DAE, we must apply the adapted Filippov method given by (10) and (12) to find the relevant saltation matrices. It is important to note that no adaptation to the DAE needs to be made. The equations derived in (10) and (12) can be applied directly to the DAE model of the system.
] , As has been shown in the literature, the integrator is effectively isolated from the rest of the system dynamics. It introduces an eigenvalue very close to +1 and it does not contribute to the systems stability. It must be ignored when estimating the stability margin. 18 For these reasons, it has been omitted from the calculation of the saltation matrix and transition matrix. It has been included in the simulation as it plays a role in deciding the switching condition.
The saltation matrix is found using the above reduced formulas and the transition matrix is numerically calculated by solving (18) and (19). Since the formulation of the transition matrix and the saltation matrix are: The monodromy matrix will take the form: .
There are N D + N A eigenvalues of this matrix. However, since the last N A columns are zero, there are N A zero eigenvalues. The Ψ ON S Ψ Ψ OFF block of the monodromy matrix does not factor into the calculation of the eigenvalues. Hence, it is sufficient to only consider the following reduced monodromy matrix for stability analysis: Using (23) and (26), the eigenvalues of the monodromy matrix are determined to be: It is clear that the eigenvalues of the system move outside the unit circle across = −1, which indicates a period-doubling bifurcation. This matches Figure 2 and serves to verify the derived equations.
We will compare the results obtained if the system is assumed to be a two-dimensional piecewise-linear ODE as opposed to the hybrid DAE approach proposed in this article.

Conventional method
DC-DC converters are systems typically described by:̇x where j denotes the different topologies of the DC-DC converter, and A ON/OFF are given by the respective f x . The state transition matrix is given by Φ(t, t 0 ) = e A(t−t 0 ) if the system is LTI. It is important to note that v in does not appear in the state transition matrix as, for a LTI system, the input will affect both the perturbed orbit and unperturbed orbit equally as Bv in does not depend on the state variables. As shown in Section 5.2.1Ψ ON/OFF and S Ψ do not have an effect on the eigenvalues of the monodromy matrix. It is sufficient to consider For these reasons, it would be natural to assume that the fluctuations of v in do not need to be taken into account for the state transition matrix and, hence, the monodromy matrix of a PV-fed boost converter would be given by: (28) where d is the duty cycle and S Φ calculated according to (23). However, for R = 43.5 Ω, (28)  . This indicates that the system is unstable. However, according to Figure 2, it is clear that the system operates with the desired stable period-1 orbit.
The reason these assumptions are not valid is because the system is modeled as a set of DAEs. It is a nonlinear system and cannot be described by (27). The algebraic constraint can cause the perturbation to shrink or grow over time. Hence, it can have an effect on the stability of the system and this must be taken into account over the entire switching period. Hence, Φ(t, t 0 ) is not given by e (f x (t−t 0 )) . A far better estimation is given by Φ(t, t 0 ) = e (f x −f y g −1 y g x )(t−t 0 ) .  (18), the eigenvalues correctly predict the bifurcation point. Calculating the state transition matrix according to Φ(t, t 0 ) = e (f x (t−t 0 )) gives incorrect eigenvalues of the system and the bifurcation point is not correctly calculated.

Conclusions: CMC control
The conclusions drawn from this case study are: • Conventional ODE-based algorithms to calculate the state-transition matrix, and hence the monodromy matrix, are not applicable for PV-fed DC-DC converters. As the system is not LTI, the algebraic constraint must be taken into account.
• While it is important to include the algebraic terms when calculating the transition matrix, the perturbation of the algebraic variables does not need to be calculated when considering the stability of the system. Hence, S Ψ does not need to be calculated, only S Φ does. For very large simulations, this will reduce the computational burden as calculating the transition matrix is not trivial and computationally expensive. 30 • If n T y = 0, then the reduced form of the saltation matrix which closely resembles that widely found in literature (such as in Reference 18) can be used.
In the next section, we will find the eigenvalues of the monodromy matrix for a PV-fed boost converter with MPPT.

Maximum power point tracking
MPPT is a technique used to maximize the power extraction of PV panels regardless of the environmental conditions. As shown in Reference 11, the MPPT algorithm is far more efficient compared with conventional control techniques. Unlike these conventional techniques, MPPT aims to maximize the output of the PV panel. In this algorithm, the input to the boost converter, that is, the output of the PV cell, is compared against a suitable reference voltage, v ref . The resulting signal, v con , is then amplified by a proportional controller, K p , to create the control signal v con : The control signal is then compared with a ramp signal. Thus, the switching surface is given by: At the start of the switching period, the switch is closed (ON). When v con = v ramp , the switch opens (OFF) and remains open for the rest of the switching period. Figure 3 illustrates a bifurcation diagram of the PV-fed boost converter with MPPT. It follows a similar pattern to Figure 2. For low values of R, the system operates with a stable period-1 orbit. The critical load value is at R ≈ 37 Ω as the system undergoes a period-doubling bifurcation for values greater than this. A further bifurcation point can be seen at R = 41.95 Ω.
In order to analytically determine the bifurcation point, we find S Φ , Φ ON , and Φ OFF and determine the eigenvalues of (26). ] , The MPPT control algorithm differs from CMC in that the algebraic variable, v in , is compared against a reference voltage as opposed to a state variable. Hence, the algebraic variable forms part of the switching manifold in (29) and n T y ≠ 0. For this reason, the conventional formulation of the Filippov method, shown in (23), cannot be applied. The resulting eigenvalues of the monodromy matrix will be incorrect and this may lead to an inaccurate prediction of the bifurcation points of the system. This case study highlights the relevance and requirement of this work for maximizing the power output of PV solar panels. Figure 4 shows the time-domain response for the inductor current, i L , with R = 37 Ω and R = 37.1 Ω. It is clear, that the system operates with the desired period-1 orbit for R = 37 Ω hence, we determine that the system is stable. An undesirable period-2 orbit can be seen for R = 37.1 Ω hence, the system is deemed to be unstable. The system loses stability at R ≈ 37.1 Ω. We will now find the eigenvalues for both cases.
In order to find the reduced monodromy matrix given in (26), we find Φ ON and Φ OFF by solving the variational equation for R = 37 Ω and R = 37.1 Ω. The saltation matrix is calculated in two manners. The first uses the adapted method to find S Φ given by (10). The second method uses the conventional formula to find S Φ that is widely given in literature such as Reference 18. This method is equivalent to setting n T y = 0 of (10) as this ignores the effect of the algebraic constraint. The eigenvalues for the system with R = 37 Ω and R = 37.1 Ω are presented in Table 3.
When R = 37 Ω, using both the adapted method and the conventional method for finding S Φ results in both sets of eigenvalues being placed inside the unit circle. Both methods indicate that the system is stable and operating with a period-1 orbit as shown in Figure 4A. However, it is important to note that the adapted method shows that one of the eigenvalues is very close to = −1. Hence, the system is on the cusp of instability, which may result in a period-doubling bifurcation if any of the system parameters vary. The conventional method places the eigenvalues further away from the edge of the unit circle. The distance an eigenvalue is from the edge of the unit circle can be used as a metric to determine the relative stability of a system. This concept is similar to gain and phase margin and can be used to indicate how far away from the instability the system is. If a control engineer was to interpret the results using the conventional method, they may believe that the system is robust to noise or parameter variations. However, R shifting from 37 to 37.1 Ω causes the system to undergo a period-doubling bifurcation. These types of bifurcations can greatly increase the current ripple and damage the DC-DC converter, the solar panel, or the interconnecting devices. Hence, the incorrect calculation of may lead to designs that damage devices. For this reason, it is imperative that we employ the correct method for finding the eigenvalues of the monodromy matrix.
When R = 37.1 Ω, the adapted method for finding S Φ correctly shows that the eigenvalues now lie outside the unit circle and that the system is unstable. However, the conventional method no longer gives the correct indication of stability, instead it indicates that the eigenvalues lie inside the unit circle. Figure 4B clearly shows that the system is operating with a period-2 orbit and hence, is unstable.
This case study has served to show the importance of considering the effect of the algebraic terms in the calculation of the saltation matrix and, hence, the monodromy matrix. If n T y ≠ 0, the algebraic terms have a significant effect on the stability of the system. This is due to the fact that they feature in the switching manifold. Furthermore, this section has also served to show the simplicity of the proposed approach. While the derivation of the new formulas to calculate the saltation matrix is mathematically complicated, its application does not greatly differ from the approach of applying the standard Filippov method. Hence, the mathematical burden is not greatly increased and the approach retains the main benefit of the Filippov method, which is the simplicity of its application. However, it is important to note that the transition matrix must be calculated numerically as part of the simulation and cannot be calculated in isolation unlike the state transition matrix for equivalent linear time-invariant ODEs.

CONCLUSION
This study, and the available literature discussed in Section 2, has shown that PV-fed DC-DC converters may undergo bifurcations, leading to unstable and undesired outputs. Traditionally, these systems have been modeled as a set of ODEs where the effect of the algebraic constraint imposed by the PV panel has been neglected during the stability analysis or an additional state variable has been created through the inclusion of additional input capacitors allowing the algebraic constraint to be expressed in the form of an ODE.
In this work, we propose modeling the system as a set of DAEs. However, the typical methods for studying DC-DC converters have been developed around modeling such systems as ODEs. Hence, in order to study PV-fed DC-DC converters modeled as a set of DAEs, we demonstrate how the Filippov method can be adapted accordingly. The main advantage of the proposed approach is that it is a general method as, regardless of the converter topology or the number of algebraic constraints, the derived formulas to perform stability analysis can be applied while maintaining the main benefit of the Filippov method which is the simplicity of its application. This reformulation of the Filippov method simplifies the analysis as no additional reformulation of the system equations needs to be performed.
We consider two case studies with the purpose of demonstrating two separate scenarios. In the first, we consider CMC. In CMC, there is no algebraic term in the switching manifold. Hence, the reformulated formulas to find the saltation matrix reduces to the conventional formula for the Filippov method found in the literature. The calculation of the saltation matrix can be treated as if it were an ODE. However, since the system is a DAE, it is shown that the calculation of the transition matrix differs and must be done so numerically.
In the second case study, we study a PV-fed boost converter with MPPT. It is shown in Section 5.3, since the algebraic term is present in the switching manifold, the equation to find the saltation matrix differs to that when the system is an ODE. It is shown that the application of the derived formula, presented in Section 4.1, retains the main benefit of the Filippov method, which is the simplicity of its application compared with other available methods outlined in Section 2. The approach was verified through brute-force simulations and experimental results can be seen in Reference 11.
The derived formula for the saltation matrix can be applied to any PV-fed DC-DC converter modeled as a hybrid DAE of any index with any number of algebraic constraints nor is it limited to the boost converter. The approach can be applied to any DC-DC converter, such as the isolated forward converter or the transformerless quadratic boost converter among many others, once the system model is a set of hybrid DAEs.

PEER REVIEW INFORMATION
Engineering Reports thanks Abdelmajid Abouloifa, Abdelali El Aroudi, and Meng Huang for their contribution to the peer review of this work.

CONFLICTS OF INTEREST
The authors declare no potential conflict of interest.