A New Control Strategy for Photovoltaic System Connected to the Grid via Three-Time-Scale Singular Perturbation Technique with Performance Analysis

Abstract

This chapter addresses the problem of controlling single-phase grid-connected photovoltaic system through a full bridge inverter with L-filter. The control objectives are threefold: (i) forcing the voltage in the output of photovoltaic panel to track a reference. This reference has been obtained from the maximum power point tracking strategy; (ii) guaranteeing a tight regulation of the DC-link voltage; and (iii) ensuring a satisfactory power factor correction (PFC) at the grid such as the currents injected must be sinusoidal with the same frequency and the same phase as the grid voltage. The considered control problem entails several difficulties including: (i) the high dimension and strong nonlinearity of the system; (ii) the changes in atmospheric conditions. The problem is dealt with by designing a synthesized nonlinear multi-loop controller using singular perturbation technique, in which a three-time-scale dynamics is artificially induced in the closed-loop system. A formal analysis based on the three-time-scale singular perturbation technique and the averaging theory is developed to proved that all control objectives are asymptotically achieved up to small harmonic errors (ripples). The performance of the proposed approach and its strong robustness with respect to climate changes are evaluated based on the various simulations results carried out under Matlab/Simulink software.

1. Introduction

Due to dramatic increase in energy consumption and thrust to reduce carbon and greenhouse gas emissions from the traditional electric power generation systems, photovoltaic (PV) power generators have gained a great popularity in recent years. Indeed, photovoltaic systems produce electric power without harming the environment, transforming a free inexhaustible source of energy, solar radiation, into electricity. Furthermore, the major advantage of the photovoltaic systems is to meet the basic power requirement of non-electrified remote areas, where grid power has not yet reached. Also, there are other advantages such as the declining cost and prices of solar modules. On the other hand, the importance of PV systems in the solar industry makes these systems more efficient and reliable, especially for utility power in distributed generation (DG) at medium and low voltages power systems [1]. All these considerations assure a promising role for PV generation systems in the near future.

On the other hand, many technical problems, such as untimely failures, could be found on electronic systems related in particular to the transfer and conversion of this energy to the network. Today, most conversion systems often suffer from low yields in real production sites. To meet the requirements of the new international standards on expected performance on associated conversion systems, it is important to make a research effort to solve the many control problems associated with the static power converter and bring this area to a degree of sufficient maturity to make them industrial products in their own right. One of the difficulties caused by the use of a photovoltaic conversion chain is focused on the problem of non-perfect control of the chain between the photovoltaic generator itself and the continuous or alternative type of load.

The efficiency of a PV plant is affected mainly by three factors: the efficiency of the PV panel, the efficiency of the static power converter and its control, and the efficiency of the maximum power point tracking (MPPT) algorithm.

PV grid-connected systems represent the most important field applications of solar energy [2, 3, 4]. In general, the power converter interface, from PV module (the DC source) to the load or to the grid, consists of two-stage converters: The first-stage DC/DC converter is usually used to boost the PV voltage and to implement the maximum power point technique. While the second stage is used to convert this power into high-quality AC voltage, with power factor correction (PFC) respecting to the power supply grid (i.e. sinusoidal and in phase with the AC supply voltage).

Maximum power point tracking (MPPT) is required to match the PV array power to the environmental changes achieving to extract the maximum power output from a solar cell [5]. To this end, different MPPT techniques have been proposed such as incremental conductance [6], perturbation and observation (P&O) [7, 8], the hill-climbing, and some other special methods, such as neural networks, fuzzy logic technique [9]. Among all available techniques, a simple and effective MPPT of incremental conductance algorithm is applied to attain the maximum power of PV array in different solar irradiance and temperature condition parameters.

In order to provide a stable controller of DC/DC and DC/AC converters, many linear control methods have been proposed using many methods, such as a fuzzy-proportional integral controlled [10], a simple PR controller [11] where the performances have been illustrated by experimental result. However, in both proposed controllers, the problem of maximizing PV power transfer is not accounted for in the controller design. In contrast to linear control, nonlinear approaches can optimize the dynamic performance of system, such as sliding mode [12], fuzzy-sliding mode [13], feedback linearization [14], singular perturbation technique [15], and many others works [16, 17, 18]. In light of the previous descriptions, no theoretical analysis is made to formally prove that the closed-loop control performances are actually achieved.

In this chapter, a multi-loop nonlinear controller is designed and developed via singular perturbation technique (Chapter 11 in Refs. [19, 20]), as was shown in Refs. [21, 22], where three-time-scale dynamics is artificially induced in the closed-loop system. The control objectives are threefold: (i) achieving the MPPT for the PV array; (ii) ensuring a tight regulation of the DC-link voltage; and (iii) ensuring a grid connection with unity power factor (PF). These objectives must be met despite changes of the climatic variables (temperature and radiation). A theoretical stability analysis, for the closed-loop system, is provided using the three-time-scale singular perturbation technique [23, 24] and averaging technique (Chapter 10 in Refs. [19, 25]. The three-time-scale analysis allows to construct a suitable composite Lyapunov function candidate for the closed-loop photovoltaic system, and the stability properties of the resulting subsystems are analyzed providing mathematical expressions for the upper bounds of the singularly perturbed parameters.

Compared to previous works, the contribution of the new nonlinear controller enjoys several interesting features including the following:

• Several control objectives are simultaneously taken into account such as: MPPT, DC regulation, and PFC, whereas only some of these objectives have been tackled in previous works [10, 11].

• A theoretical analysis will prove, using three-time-scale singular perturbation and averaging technique, that the desired multiple objectives are achieved. Such a formal analysis was missing in the previous works [12, 13, 14, 15].

• The nonlinearity of the controlled system was preserved in the controller design in order to keep all the properties of the studied system, whereas it is partly or totally ignored in previous controllers [14].

• By including of three-time-scale dynamics in the full-order closed-loop system can ensure to achieve desired properties, such as robust zero steady-state error of the reference input realization, desired output performance specifications (overshoot, settling time), and insensitivity of the output transient behavior with respect to parameter variations and external disturbances.

The content of this chapter is outlined as follows: in Section 2, the grid-connected PV system is described and modeled. Section 3 is devoted to the cascade nonlinear controller design and its performances are formally analyzed in Section 4. The global performance of the closed-loop photovoltaic system will be illustrated by numerical simulation using MATLAB/SIMULINK tool in section 5. A conclusion and a reference list end the chapter.

2. System description

This section describes the modeling of photovoltaic system connected to the grid. The power circuit topology used in the proposed single phase grid connected to the photovoltaic array is shown in Figure 1. It consists of the following components: (i) a photovoltaic array which consists of an arrangement of $N{}_s$-series and $N{}_p$-parallel strings; (ii) an input capacitor $C_{\mathit{pv}}$ and a DC-DC boost converter used to increase the voltage level and achieve MPPT for photovoltaic array; (iii) a DC link capacitor $C_{\mathit{dc}}$; and (iv) a single-phase full-bridge inverter including four power semiconductors with $L$ filter that is used to provide the energy to the grid and ensure power factor correction.

Figure 1.

Single phase grid-connected PV system.

Typical (Ip-Vp) characteristics of solar cells arranged in $N{}_p$-parallel and $N{}_s$-series can be found in many places (see, e.g. [26]). The PV array module considered in this paper is of type KC200GT. In this chapter, a simple and effective MPPT of incremental conductance algorithm is applied to attain the maximum power of PV array in different solar irradiance and temperature condition parameters.

By analyzing the circuit and applying the well-known Kirchhoff laws, the system of Figure 1 can be described by the following set of differential equations:

$$\frac{{\mathit{di}}_i}{\mathit{dt}}=\frac{v_{\mathit{pv}}}{L_i}-\left(1-u_1\right)\frac{v_{\mathit{dc}}}{L_i}-\frac{r_i}{L_i}{i}_i$$E1

$$\frac{{\mathit{dv}}_{\mathit{pv}}}{\mathit{dt}}=\frac{1}{C_{\mathit{pv}}}\left(i_{\mathit{pv}}-i_i\right)$$E2

$$\frac{{\mathit{di}}_g}{\mathit{dt}}=-\frac{v_g}{L_o}-\frac{r_g}{L_g}{i}_g+u_2\frac{v_{\mathit{dc}}}{L_o}$$E3

$$\frac{{\mathit{dv}}_{\mathit{dc}}}{\mathit{dt}}=-\frac{1}{C_{\mathit{dc}}}{u}_2{i}_g+\left(1-u_1\right)\frac{i_i}{C_{\mathit{dc}}}$$E4

where $v_{\mathit{pv}}$ and $i_{\mathit{pv}}$ are, respectively, the photovoltaic generator voltage and current. $v_{\mathit{dc}}$ and $i_{\mathit{dc}}$ are, respectively, DC link voltage and current. $i_i$ designates the input current chopper, $C_{\mathit{dc}}$ is DC link capacitor, $r_i$ and $r_g$ are, respectively, the equivalent series resistances (ESR) of input inductance $L{}_i$ and the filter inductance $L_g$. $v_g$ and $i_g$ are, respectively, the voltage and current of the grid. Here, the grid voltage is defined by $v_g=E_{g}sin\left({\omega }_{g}t\right)$, where $E_g$ and ${\omega }_g$ denote the constant amplitude and the constant angular frequency. The switching functions ${\mu }_1$ and ${\mu }_2$ are defined by:

$${\mu }_1=\left\{\begin{array}{l}1\text{if}{s}_0\text{is}\mathit{ON}\\ 0\text{if}{s}_0\text{is}\mathit{OFF}\end{array}\right.$$E5

$${\mu }_2=\left\{\begin{array}{l}1\text{if}\left(s_1,{\mathrm{s}}_4\right)\text{is}\mathit{ON}\text{and}\left(s_2,{\mathrm{s}}_3\right)\text{is}\mathit{OFF}\\ -1\text{if}\left(s_1,{\mathrm{s}}_4\right)\text{is}\mathit{OFF}\text{and}\left(s_2,{\mathrm{s}}_3\right)\text{is}\mathit{ON}\end{array}\right.$$E6

The instantaneous model (1)–(4) cannot be used directly for the development of continuous control laws since it involves, as input variables, the binary signal ${\mu }_1$ and ${\mu }_2$. To overcome this inconvenience, the average model is used [27]. Therefore, the state variables $i_i$, $v_{\mathit{pv}}$, $i_g$, and $v_{\mathit{dc}}$ are replaced by their average values $x_1$, $x_2$, $x_3$, and $x_4$ over a cutting period. The control inputs $u_1$ and $u_2$ denote the average values of ${\mu }_1$ and ${\mu }_2$, respectively.

$$\frac{{\mathit{dx}}_1}{\mathit{dt}}=\frac{1}{L_i}{x}_2-\frac{1}{L_i}{x}_4-\frac{r_i}{L_i}{x}_1+\frac{1}{L_i}{x}_4{u}_1$$E7

$$\frac{{\mathit{dx}}_2}{\mathit{dt}}=\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}-\frac{1}{C_{\mathit{pv}}}{x}_1$$E8

$$\frac{{\mathit{dx}}_3}{\mathit{dt}}=-\frac{r_g}{L_g}{x}_3-\frac{v_g}{L_g}+\frac{1}{L_g}{x}_4{u}_2$$E9

$$\frac{{\mathit{dx}}_4}{\mathit{dt}}=-\frac{1}{C_{\mathit{dc}}}{x}_3{\mu }_2-\frac{1}{C_{\mathit{dc}}}{x}_1{u}_1+\frac{1}{C_{\mathit{dc}}}{x}_1$$E10

3. Controller design

3.1 Input inductor current regulation and PFC objectives

3.1.1 Control law design

The first control objective is to enforce the photovoltaic voltage $x_2$ to track, as closely as possible, the optimal point $V_M^{\prime }$ (called Regulator 2). However, it is well-known that the boost converter has a non-minimum phase feature. Such an issue is generally dealt with by resorting a cascaded loop design strategy that starts with the input current loop (Regulator 1), as it is shown in Figure 2. More specifically, the controller makes the input inductor current $x_1$ to track a reference signal $x_1^{\ast }$, the latter is determined from (Regulator 2).

Figure 2.

Schematic diagram of the proposed controller for single phase grid-connected PV system.

In parallel with the input current controller, the network current controller (Regulator 3) will be designed for power factor correction requirement that amounts to forcing the network current $x_3$ to match the reference signal of the form $x_3^{\ast }=\beta sin\left({\omega }_{g}t\right)$. It means that the grid current $x_3$ should be sinusoidal and in phase with the AC grid voltage $v_g$, with $\beta$ is a signal. In fact, the latter is allowed to (and actually will) be time-varying but it must converge to a positive constant value. For both objectives (i.e. the input current $x_1$ regulation and PFC objectives), let us consider Eqs. (7), (9) in the following form:

$${\dot{X}}_{1,3}=M_{1}U+M_2$$E11

where $X_{1,3}=\left(\begin{array}{c}{x}_1\\ x_3\end{array}\right)$, $M_1=\left(\begin{array}{c}\frac{x_4}{L_i}\\ \frac{x_4}{L_g}\end{array}\right)$, $M_2=\left(\begin{array}{c}-\frac{r_i}{L_i}{x}_1+\frac{x_2}{L_i}-\frac{x_4}{L_i}\\ -\frac{r_g}{L_g}{x}_3-\frac{v_g}{L_g}\end{array}\right)$, and $U=\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)$.

Remark 1. Under condition that the DC-link voltage remains all the time positive, it can be showed that $0<M_{1min}\le M_1\le M_{1\mathit{Max}}<\infty$ and $\left|M_2\right|\le M_{2\mathit{Max}}<\infty$ are satisfied.

Let us introduce the following current tracking errors:

$$e_i=x_i^{\ast }-x_i\left(i=1,3\right)$$E12

Then, the reference model can be constructed in the following form

$$\frac{{\mathit{dx}}_i}{\mathit{dt}}=\frac{e_i}{T_i}\stackrel{\mathrm{\Delta }}{=}{D}_i^{\prime }\left(x_i^{\ast },x_i\right)\left(i=1,3\right)$$E13

where $T_1$ and $T_3$ denote the time constants and they are selected based on the desired settling time, respectively, for the currents $x_1$ and $x_3$. Based on (13), the realization errors of the desired behaviors of ${\dot{x}}_1$ and ${\dot{x}}_3$, namely ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_3$, are given by

$${\mathrm{\Omega }}_i=D_i^{\prime }\left(x_i^{\ast },x_i\right)-{\dot{x}}_i\left(i=1,3\right)$$E14

Therefore, the control problem $\underset{t\to \infty }{lim}{e}_i\left(t\right)=0\left(i=1,3\right)$ corresponds to the insensitivity condition defined by

$${\mathrm{\Omega }}_i=0\left(i=1,3\right)$$E15

Doing so, the behaviors of ${\dot{x}}_1$ and ${\dot{x}}_3$ with prescribed dynamics of expression (13) will be fulfilled. Replacing in Eq. (14)${\dot{x}}_1$ and ${\dot{x}}_3$ by their expressions (7), (9), and since the requirement (15), one gets

$$D^{\prime }-M_2-M_{1}U=0$$E16

Hence, there is an isolated root called the inverse dynamic solutions for $U$ given by

$$U^{\mathit{id}}=\left(\begin{array}{c}{u}_1^{\mathit{id}}\\ u_2^{\mathit{id}}\end{array}\right)=\left(\begin{array}{c}\frac{L_i}{x_4}\left(\frac{e_1}{T_1}-\left(-\frac{r_i}{L_i}{x}_1+\frac{1}{L_i}{x}_2-\frac{1}{L_i}{x}_4\right)\right)\\ \frac{L_g}{x_4}\left(\frac{e_3}{T_3}+\frac{x_4+v_g+r_g{x}_3}{L_g}\right)\end{array}\right)$$E17

The control variable, namely $U$, has emerged in Eq. (16). At this point, an appropriate control law with the first derivative in feedback has to be found, so that the ($e_1$, $e_3$)-systems are made asymptotically stable. As the objective is to drive the error to zero, it is natural to choose the Lyapunov functions candidate

$$\left(\begin{array}{c}{V}_3\\ V_4\end{array}\right)=\frac{1}{2}\left(\begin{array}{c}{\left({\mathrm{\Omega }}_1\left(u_1\right)\right)}^2\\ {\left({\mathrm{\Omega }}_3\left(u_2\right)\right)}^2\end{array}\right)$$E18

From expressions (14) and (16), it can be easily checked that the following time-derivatives

$$\frac{d}{\mathit{dt}}\left(\begin{array}{c}{V}_3\\ V_4\end{array}\right)=\left(\begin{array}{c}\frac{\partial V_3}{\partial u_1}\frac{{\mathit{du}}_1}{\mathit{dt}}\\ \frac{\partial V_4}{\partial u_2}\frac{{\mathit{du}}_2}{\mathit{dt}}\end{array}\right)=\left(\begin{array}{c}\frac{{\mathrm{\Omega }}_1\partial {\mathrm{\Omega }}_1}{\partial u_1}\frac{{\mathit{du}}_1}{\mathit{dt}}\\ \frac{{\mathrm{\Omega }}_3\partial {\mathrm{\Omega }}_3}{\partial u_2}\frac{{\mathit{du}}_2}{\mathit{dt}}\end{array}\right)=-\left(\begin{array}{c}\frac{{\mathrm{\Omega }}_1{x}_4}{L_i}\frac{{\mathit{du}}_1}{\mathit{dt}}\\ \frac{{\mathrm{\Omega }}_3{x}_4}{L_g}\frac{{\mathit{du}}_2}{\mathit{dt}}\end{array}\right)$$E19

are made negative-definite using the following control laws:

$$\frac{d}{\mathit{dt}}\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}\frac{k_1}{{\epsilon }_1\epsilon {}_2}{\mathrm{\Omega }}_1\\ \frac{k_3}{{\epsilon }_1\epsilon {}_2}{\mathrm{\Omega }}_3\end{array}\right)$$E20

At this point, $k_1$ and $k_3$ are any design parameters, ${\epsilon }_1$ and $\epsilon {}_2$ are small positive parameters. In view of $M_1>0$, it is easily seem from expressions (11) and (14) that

$$\frac{d}{\mathit{dt}}\left(\begin{array}{c}{V}_3\\ V_4\end{array}\right)=-\left(\begin{array}{c}\frac{k_1}{{\epsilon }_1\epsilon {}_2}\frac{x_4{\mathrm{\Omega }}_{\mathit{1}}^{\mathit{2}}}{L_i}\\ \frac{k_3}{{\epsilon }_1\epsilon {}_2}\frac{x_4{\mathrm{\Omega }}_3^{\mathit{2}}}{L_g}\end{array}\right)<0\left({\mathrm{\Omega }}_{\mathrm{i}}\ne 0\right)$$E21

for a positive values of $k_i>0\left(i=1,3\right)$. Therefore, from expressions (14) and (20), the discussed nonlinear control laws are formulated as follows

$${\epsilon }_1{\epsilon }_2\frac{d}{\mathit{dt}}\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}{k}_1\left(\frac{e_1}{T_1}-\frac{{\mathit{dx}}_1}{\mathit{dt}}\right)\\ k_3\left(\frac{e_3}{T_3}-\frac{{\mathit{dx}}_3}{\mathit{dt}}\right)\end{array}\right)$$E22

3.1.2 Singular perturbation system of the inner current loops

Consider the closed-loop system of inner loops composed of the Eqs. (7)–(10) and the control laws (22), which can be rewritten in the following form

$${\epsilon }_1{\epsilon }_2\frac{\mathit{dZ}}{\mathit{dt}}=h\left(X,Z,t\right)$$E23

$$\frac{\mathit{dX}}{\mathit{dt}}=f^{\prime }\left(X,Z,t\right)$$E24

with $Z={\left(z_1{z}_2\right)}^T={\left(\begin{array}{cc}{u}_1& u_2\end{array}\right)}^T$, $X={\left(\begin{array}{cccc}{x}_1& x_2& x_3& x_4\end{array}\right)}^T$, $h\left(X,Z,t\right)=\left(\begin{array}{c}{k}_1\left(\frac{e_1}{T_1}-\frac{x_4}{L_i}{z}_1+\frac{r_i}{L_i}{x}_1-\frac{x_2}{L_i}+\frac{x_4}{L_i}\right)\\ k_3\left(\frac{e_3}{T_3}-\frac{x_4}{L_g}{z}_2+\frac{v_g}{L_g}+\frac{r_g}{L_g}{x}_3\right)\end{array}\right)$, $f^{\prime }\left(X,Z,t\right)=\left(\begin{array}{c}{z}_1\frac{x_4}{L_i}-\frac{r_i}{L_i}{x}_1+\frac{x_2}{L_i}-\frac{x_4}{L_i}\\ -\frac{x_1}{C_{\mathit{pv}}}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}\\ -\frac{v_g}{L_g}-\frac{r_g}{L_g}{x}_3+\frac{x_4}{L_g}{z}_2\\ \frac{x_1}{C_{\mathit{dc}}}-\frac{x_1}{C_{\mathit{dc}}}{z}_1-\frac{x_3}{C_{\mathit{dc}}}{z}_2\end{array}\right)$.

Now, we go to the fast time ${\tau }_1=t/{\epsilon }_1{\epsilon }_2$. Then for ${\epsilon }_1=0$, the ultra-fast dynamic subsystem (UFDS) is given by

$$\frac{\mathit{dZ}}{d{\tau }_1}=h\left(X,Z,t\right)$$E25

$$\frac{\mathit{dX}}{d{\tau }_1}=0$$E26

After the rapid decay of transients in expression (25), the steady state (more precisely, quasi-steady state) tends toward an equilibrium $Z^e=\tilde{h}\left(X\right)$. The manifold defined by $Z^e$ is called the slow manifold, which is given by

$$Z^e=\left(\begin{array}{c}{z}_1^e\\ z_2^e\end{array}\right)=\left(\begin{array}{l}\frac{L_i}{x_4}\left(\frac{e_1}{T_1}+\frac{x_4-x_2+r_i{x}_1}{L_i}\right)\\ \frac{L_g}{x_4}\left(\frac{e_3}{T_3}+\frac{v_g+r_g{x}_3}{L_g}\right)\end{array}\right)$$E27

Remark 2.

1. During the fast transient in expression (25), the variables $X$ are treated as the frozen parameters.

2. The equilibrium point given by expression (27) involves a division by the DC link voltage $x_4$, from a practical point of view this division is not a problem because the DC link voltage remains all the time positive for the power converter to work correctly.

By substituting of this equilibrium $Z^e$ into Eq. (24), the slow dynamic subsystem (SDS) of inner loops takes place on the slow manifold, according to the equation

$$\frac{\mathit{dX}}{\mathit{dt}}=\left(\begin{array}{c}\frac{e_1}{T_1}\\ -\frac{x_1}{C_{\mathit{pv}}}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}\\ \frac{e_3}{T_3}\\ -\frac{L_o{x}_3}{C_{\mathit{dc}}{x}_4}\left(\frac{e_3}{T_3}+\frac{v_g+r_g{x}_3}{L_g}\right)-\frac{L_i{x}_1}{C_{\mathit{dc}}{x}_4}\left(\frac{e_1}{T_1}+\frac{r_i{x}_1-x_2}{L_i}\right)\end{array}\right)$$E28

Proposition 1. Consider the closed-loop system composed of Eqs. (23) and (24). For ${\epsilon }_1\to 0$, the system takes the singular perturbation form where the UFDS is defined by equation (25), while the SDS of inner loop is defined by equation (28). Under the considerations given by Remark 2, one has the following properties

1. If the design parameters $k_i$ ($i=1,3$) are positives, the UFDS (25) will be exponentially stable, and $Z$ converge exponentially fast to $Z^e$.

2. The behaviors of $x_i$ ($i=1,3$) are prescribed by the stable reference equations of the form ${\mathit{dx}}_i/\mathit{dt}=\left(x_i^{\ast }-x_i\right)/T_i$. Then, the requirements $\underset{t\to \infty }{lim}{e}_i=\underset{t\to \infty }{lim}\left(x_i^{\ast }-x_i\right)=0$ are maintained.

3.2 MPPT and DC bus voltage regulation objective

The second step consists in completing the inner control loops by outer control loops for PV voltage (Regulator 2) and DC-link voltage (Regulator 4). The aim is now to enforce the photovoltaic voltage $x_2$ and the DC-link voltage $x_4$ to track, respectively, the optimal point $x_2^{\ast }=V_M^{\prime }$ and a given reference voltage $x_4^{\ast }=V_{\mathit{dc}}^{\prime }$, such that tuning laws for the ratio $\beta$ and $x_1^{\ast }$ must be designed. According to the three-time-scale design methodology that is employed in this work, the general formulation of the three-time-scale singular perturbed systems requires the system to possess three different time scales. To this end, the voltage loops will be slow compared to the transients of the current loops. Therefore, the design parameters for voltage loops, in particular (${\epsilon }_{\mathit{2}}$, $T_2$, and $T{}_4$) must satisfy: $0<{\epsilon }_1{\epsilon }_{\mathit{2}}\ll {\epsilon }_{\mathit{2}}\ll 1$, and $\left(T_1,T_3\right)<\left(T_{\mathit{2}},T_4\right)$. In addition, the steady states for the current $x_1$ and $x_3$ yield, respectively, $x_1^{\ast }$ and $x_3^{\ast }$. Therefore, the SDS of inner loop given by equation (28) will be reduced to

$$\frac{{\mathit{dX}}_r}{\mathit{dt}}=\frac{d}{\mathit{dt}}\left(\begin{array}{c}{x}_2\\ x_4\end{array}\right)=\left(\begin{array}{c}\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}-\frac{x_1^{\ast }}{C_{\mathit{pv}}}\\ \frac{x_2{x}_1^{\ast }-v_g{x}_3^{\ast }-r_i{\left(x_1^{\ast }\right)}^2-r_g{\left(x_3^{\ast }\right)}^2}{C_{\mathit{dc}}{x}_4}\end{array}\right)$$E29

Now, the first step is to establish the relation between the ratio $\beta$ (which acts as the control input of the outer loop) and the DC-link voltage $x_4$ (representing the output of the outer loop).

3.2.1 Relation between $\left(\beta \right)$ and $\left(x_4\right)$, and the control law

The relation between the ratio $\beta$ and the DC-link voltage $x_4$ is the subject of the following proposition.

Proposition 2. We consider the second equation of expression (29) and the power factor correction requirement defined by $x_3^{\ast }=\beta sin\left({\omega }_{g}t\right)$.

1. The relation between $x_4$ and $\beta$ is described in the following results

   $$\frac{{\mathit{dx}}_4}{\mathit{dt}}=\frac{x_2{x}_1^{\ast }-r_i{\left(x_1^{\ast }\right)}^2}{C_{\mathit{dc}}{x}_4}-\frac{v_g^2\beta }{C_{\mathit{dc}}{x}_4{E}_g}-\frac{r_g{v}_g^2{\beta }^2}{C_{\mathit{dc}}{x}_4{E}_g^2}$$E30

2. Therefore, the squared-voltage $x_{r,4}=x_4^2$ varies, in response to the tuning ratio β, according to the following first-order time-varying nonlinear equation:

$$\frac{{\mathit{dx}}_{r,4}}{\mathit{dt}}=2\frac{x_2{x}_1^{\ast }-r_i{\left(x_1^{\ast }\right)}^2}{C_{\mathit{dc}}}-\frac{2v_g^2\beta }{C_{\mathit{dc}}{E}_g}-\frac{2r_g{v}_g^2{\beta }^2}{C_{\mathit{dc}}{E}_g^2}$$E31

The second step is to establish control laws for the outer loops, in which $Y={\left(\begin{array}{cc}{y}_1& y_2\end{array}\right)}^T={\left(\begin{array}{cc}{x}_1^{\ast }& \beta \end{array}\right)}^T$ represent the new control inputs, while $X_r={\left(\begin{array}{cc}{x}_{r,2}& x_{r,4}\end{array}\right)}^T={\left(\begin{array}{cc}{x}_2& x_4^2\end{array}\right)}^T$ represent the new output variables. To this end, introduce the following tracking errors

$$e_j=x_j^{\ast }-x_{r,j}\left(j=2,4\right)$$E32

Then, let the desired behavior of $X_r$ be assigned by

$$\frac{{\mathit{dx}}_{r,j}}{\mathit{dt}}=\frac{e_j}{T_j}\stackrel{\mathrm{\Delta }}{=}{D}_j\left(x_j^{\ast },x_j\right)\left(j=2,4\right)$$E33

The error of the desired dynamic realization it follows

$${\mathrm{\Omega }}_j=D_j\left(x_{r,j}^{\ast },x_{r,j}\right)-{\dot{x}}_{r,j}\left(j=2,4\right)$$E34

Then, the insensitivity condition is given by

$${\mathrm{\Omega }}_j=0$$E35

Similar to the previous subsection and bearing in mind the fact that $\beta$ and their first derivative must be available, and in order to meet the requirement (35), we should apply the control law given by the following structure

$${\epsilon }_2\frac{d}{\mathit{dt}}\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)=\left(\begin{array}{c}{k}_2\left(e_2/T_2-{\mathit{dx}}_{2,r}/\mathit{dt}\right)\\ k_4\left(e_4/T_4-{\mathit{dx}}_{r,4}/\mathit{dt}\right)\end{array}\right)$$E36

3.2.2 Singular perturbation system of the outer voltage loops

Combining expressions (31) and (36), and the first equation of expression (29), one obtains

$${\epsilon }_2\frac{\mathit{dY}}{\mathit{dt}}=g\left(X_r,Y,t\right)=\left(\begin{array}{c}{k}_2\left(\frac{e_2}{T_2}+\frac{y_1-i_{\mathit{pv}}}{C_{\mathit{pv}}}\right)\\ k_4\left(\frac{e_4}{T_4}-f_2^{″}\right)\end{array}\right)$$E37

$$\frac{{\mathit{dX}}_r}{\mathit{dt}}=f^{″}\left(X_r,Y,t\right)=\left(\begin{array}{c}\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}-\frac{1}{C_{\mathit{pv}}}{y}_1\\ f_2^{″}\end{array}\right)$$E38

with $f_2^{″}=2\left(\frac{x_{r,2}{y}_1-r_i{y}_1^2}{C_{\mathit{dc}}}-\frac{E_g{v}_g^2{y}_2+r_g{v}_g^2{y}_2^2}{C_{\mathit{dc}}{E}_g^2}\right)$. The fast dynamic subsystem (FDS) is obtained by transforming the slow time-scale $t$ to the fast time-scale ${\tau }_2=t/{\epsilon }_2$, then, by setting ${\epsilon }_2=0$

$$\frac{\mathit{dY}}{d{\tau }_2}=g\left(X_r,Y,t\right)$$E39

$$\frac{{\mathit{dX}}_r}{d{\tau }_2}=0$$E40

Notice that expression (39) has an isolated equilibrium at $Y^e=\left(\begin{array}{cc}{y}_1^e& y_2^e\end{array}\right)$, which will be determined (in the mean) in Appendix. As the FDS (39) is nonlinear, the stability properties of its equilibrium can be checked through the analysis of the Jacobian matrix of the linearized version defined as follows

$$A_F=\left[\begin{array}{cc}\frac{k_2}{C_{\mathit{pv}}}& 0\\ k_4{A}_F^{\prime }& k_4{A}_F^{″}\end{array}\right]$$E41

with $A_F^{\prime }=\frac{4r_i{y}_1^e}{C_{\mathit{dc}}}-\frac{2x_{r,2}}{C_{\mathit{dc}}}$, $A_F^{″}=\frac{2v_g^2}{C_{\mathit{dc}}{E}_g}+\frac{4r_g{v}_g^2{y}_2^e}{C_{\mathit{dc}}{E}_g^2}$. Taking into account that $y_2^e$ is positive, we conclude that all eigenvalues of $A_F$ satisfy $Re\left({\lambda }_{1,2}\right)<0$ for the negative values of $k_2$ and $k_4$. Therefore, $A_F$ is Hurwitz matrix. By substituting of the equilibrium $Y^e$ into expression (38), the reduced SDS of outer loops takes place on the slow manifold, according to the equation

$$\frac{{\mathit{dX}}_r}{\mathit{dt}}=\left(\begin{array}{c}{e}_2/T_2\\ e_4/T_4\end{array}\right)$$E42

Proposition 3. Consider the system closed loop composed of expressions (37) and (38). For ${\epsilon }_2\to 0$, this system takes the singular perturbation form, where the FDS is given by expression (39) and the reduced SDS of outer loop is given by expression (42). One has the following properties

1. If the design parameters $k_2$ and $k_4$ are negative, the FDS (39), will be exponentially stable and $Y$ converge exponentially fast to $Y^e$.

2. The behaviors of $x_{r,j}$ ($j=2,4$) are prescribed by the stable reference equations of the form ${\mathit{dx}}_{r,j}/\mathit{dt}=\left(x_j^{\ast }-x_{r,j}\right)/T_j$. Then, the requirements $\underset{t\to \infty }{lim}{e}_j=\underset{t\to \infty }{lim}\left(x_j^{\ast }-x_{r,j}\right)=0$ ($j=2,4$) are maintained.

4. Control system analysis

The objective of the global stability of closed-loop system can be analyzed in the following theorem. It is shown that the control objectives are achieved (in the mean) with an accuracy that depends on the network frequency ${\omega }_g$ and the small parameters ${\epsilon }_i^{\otimes }\left(i=1,2\right)$.

Theorem. Consider the overall control system composed of the Pv panel, boost DC-DC converter and DC-AC inverter, described by the model (7)–(10), in closed loop with the multi-cascade multi-loop composed of:

• The inner regulators (23), where (${\epsilon }_1{\epsilon }_2$, $k_1$, $k_3$, $T_1$, $T_3$) are the design parameters;

• The outer regulators (37), where (${\epsilon }_2$, $k_2$, $k_4$, $T_2$, $T_4$) are the design parameters.

Then, one has the following property

1. The augmented state vector $Z={\left(\begin{array}{cc}{u}_1& u_2\end{array}\right)}^T={\left(\begin{array}{cc}{z}_1& z_2\end{array}\right)}^T$, $Y={\left(\begin{array}{cc}{x}_1^{\ast }& \beta \end{array}\right)}^T={\left(\begin{array}{cc}{y}_1& y_2\end{array}\right)}^T$, and $X^{\prime }={\left(\begin{array}{cccc}{x}_1^{\prime }& x_2^{\prime }& x_3^{\prime }& x_4^{\prime }\end{array}\right)}^T={\left(\begin{array}{cccc}{x}_1& x_2& x_3& x_4^2\end{array}\right)}^T$ undergoes the following state equations

   $${\epsilon }_1{\epsilon }_2\dot{Z}\left(t\right)=h\left(X^{\prime },Y,Z,t\right)$$E43
   $${\epsilon }_2\dot{Y}\left(t\right)=g\left(X^{\prime },Y,Z,t\right)$$E44
   $${\dot{X}}^{\prime }\left(t\right)=f\left(X^{\prime },Y,Z,t\right)$$E45

   with

   $h=\left(\begin{array}{c}{k}_1\left(\frac{e_1}{T_1}-\frac{z_1\sqrt{x_4^{\prime }}}{L_i}+\frac{r_i{x}_1^{\prime }-x_2^{\prime }+\sqrt{x_4^{\prime }}}{L_i}\right)\\ k_3\left(\frac{e_3}{T_3}-\frac{\sqrt{x_4^{\prime }}}{L_g}{z}_2+\frac{v_g+r_g{x}_3^{\prime }}{L_g}\right)\end{array}\right)$, $f=\left(\begin{array}{c}-\frac{r_i}{L_i}{x}_1^{\prime }+\frac{x_2^{\prime }}{L_i}-\frac{\sqrt{x_4^{\prime }}}{L_i}+z_1\frac{\sqrt{x_4^{\prime }}}{L_i}\\ -\frac{y_1}{C_{\mathit{pv}}}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}\\ \frac{\sqrt{x_4^{\prime }}}{L_o}{z}_2-\frac{v_g}{L_g}-\frac{r_o}{L_g}{x}_3^{\prime }\\ 2\frac{x_2^{\prime }{y}_1-r_i{y}_1^2}{C_{\mathit{dc}}}-\frac{2v_g^2{y}_2}{C_{\mathit{dc}}{E}_g}-\frac{2r_g{v}_g^2{y}_2^2}{C_{\mathit{dc}}{E}_g^2}\end{array}\right)$, $g=\left(\begin{array}{c}{k}_2\left(\frac{e_2}{T_2}+\frac{y_1}{C_{\mathit{pv}}}-\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}\right)\\ k_4\left(\frac{e_4}{T_4}+\frac{2v_g^2{y}_2}{C_{\mathit{dc}}{E}_g}+\frac{2r_g{v}_g^2{y}_2^2}{C_{\mathit{dc}}{E}_g^2}-2\frac{x_2^{\prime }{y}_1-r_i{y}_1^2}{C_{\mathit{dc}}}\right)\end{array}\right)$.

2. Let the control design parameters be selected, such that the following inequalities hold $T_1<T_2$, $T_3<T_4$, $k_1>0$, $k_2<0$, $k_3>0$, $k_4<0$, and $0<{\epsilon }_1{\epsilon }_{\mathit{2}}\ll {\epsilon }_{\mathit{2}}\ll 1$.

Then, there exist positive constants ${\rho }^{\otimes }$, ${\epsilon }_n^{\otimes }\left(n=1,2\right)$, such that for and $0<\rho <{\rho }^{\otimes }0<{\epsilon }_n<{\epsilon }_n^{\otimes }\ll 1\left(n=1,2\right)$ the system (43)–(45) has an asymptotically stable $\left(2\pi /{\omega }_g\right)$ periodic solution $X\left(t,{\epsilon }_n,\rho \right)$, $Y\left(t,{\epsilon }_n,\rho \right)$, $Z\left(t,{\epsilon }_n,\rho \right)$ that continuously depends on ${\epsilon }_n\left(n=1,2\right)$ and $\rho =1/{\omega }_g$. Furthermore, when (${\epsilon }_n=0$ and $\rho =0$), one has

$$\underset{{\epsilon }_n\to 0,\rho \to 0}{lim}X\left(t,{\epsilon }_n,\rho \right)=X_0^{\ast },\underset{{\epsilon }_n\to 0,\rho \to 0}{lim}Y\left(t,{\epsilon }_n,\rho \right)=Y_0^{\ast },\underset{{\epsilon }_n\to 0,\rho \to 0}{lim}Z\left(t,{\epsilon }_n,\rho \right)=Z_0^{\ast },$$E46

where $X_0^{\ast }={\left(\begin{array}{cccc}{x}_{1,0}^{\ast }& x_{2,0}^{\ast }& x_{3,0}^{\ast }& x_{4,0}^{\ast }\end{array}\right)}^T$, $Y_0^{\ast }={\left(\begin{array}{cc}{y}_{1,0}^{\ast }& y_{2,0}^{\ast }\end{array}\right)}^T$, and $Z_0^{\ast }={\left(\begin{array}{cc}{z}_{1,0}^{\ast }& z_{2,0}^{\ast }\end{array}\right)}^T$ with $x_{2,0}^{\ast }=V_M^{\prime }$, $x_{4,0}^{\ast }=V_{\mathit{dc}}^{\prime }$, $z_{1,0}^{\ast }=1-\frac{x_{2,0}^{\ast }}{\sqrt{x_{4,0}^{\ast }}}+\frac{r_i{x}_{1,0}^{\ast }}{\sqrt{x_{4,0}^{\ast }}}$, $z_{2,0}^{\ast }=\frac{r_o{x}_{3,0}^{\ast }}{\sqrt{x_{4,0}^{\ast }}}$, $y_{1,0}^{\ast }=i_{\mathit{pv},0}$, $y_{2,0}^{\ast }=\frac{-E_g+2\sqrt{{\left(E_g/2\right)}^2+2r_o\left(x_{2,0}^{\ast }{i}_{\mathit{pv},0}-r_i{\left(i_{\mathit{pv},0}\right)}^2\right)}}{2r_o}$.

See Appendix for the proof.

5. Simulation and discussion of results

The experimental setup is described by Figure 2 and the nonlinear controller, developed in Section 3, including the control laws (23) and (37), will now be evaluated by simulation in MATLAB/SIMULINK platform using the electromechanical characteristics of Table 1.

Parameters	Symbol	Values
Network	$E/f$	$220\sqrt{2}\mathrm{V}/50\mathrm{Hz}$
Boost	$C_{\mathit{pv}}$; $L_i$; $r_i$	$0.3\mathrm{mF}$; $6\mathrm{mH}$; $20\mathrm{m\Omega }$
L-filter	$L_o$; $r_o$	$5\mathrm{mH}$; $50\mathrm{m\Omega }$
PWM switching frequency	$f_{\mathit{PWM}}$	$14\mathrm{kHz}$
DC capacitance	$C_{\mathit{dc}}$	$7\mathrm{mF}$

Table 1.

PV system and single-phase grid characteristics.

The numerical values used for the design parameters are chosen as follows: ${\epsilon }_1=5.81\times {10}^{-2}$, ${\epsilon }_2=7.21\times {10}^{-5}$, $T_1=2.1\times {10}^{‐4}$, $T_2=9.3\times {10}^{‐3}$, $T_3=9.51\times {10}^{‐5}$, $T_4=1.43\times {10}^{‐2}$, $k_1=7.1\times {10}^{‐7}$, $k_2=-2.56\times {10}^{-6}$, $k_3=2.3\times {10}^{‐6}$, $K_4=8.1\times {10}^{‐8}$. These values have proved to be suitable based on several trials respecting the singular perturbation technique. In this simulation, we consider the KC200GT type of PV array module with $N_s=54$ and $N_p=6$.

The performances of the proposed controller are illustrated by Figures 3–5. Figure 3 shows that the DC-link voltage $x_4$ is well regulated and quickly settles down after each change in the signal reference (stepping from $600V$ to $750V$ at $t=0.4s$). The wave frame of the output current $x_3$ is showed in Figure 4. The current is sinusoidal and in phase with the network voltage complying with the PFC requirement. This is further demonstrated by Figure 5, which shows that the ratio $\beta$ takes a constant value after transient periods following the changes in reference signals.

Figure 3.

DC-link voltage.

Figure 4.

Power factor correction checking.

Figure 5.

Tuning parameter β.

5.1 Radiation variation effect

Figure 6a–c illustrates the resulting closed-loop control performances in presence of radiation changes. Specifically, the radiation takes a low, medium, and high value (equal to 800, 1000 and step to 600 W/m² at times 0, 0.4, and 0.7 s, respectively), meanwhile the temperature is kept constant, equal to $298.17K$ (i.e. $25C^{\circ }$). Figure 6a shows that the captured PV voltage varies between 182, 212, and 138 V. These values correspond to the maximum points. Figure 6b shows that the DC-link voltage regulation is recovered after a short transient period following each change of the irradiation. Figure 6c shows that the current amplitude changes whenever the radiation varies. It is seen that the output current $x_3$ and the grid voltage $v_g$ are actually sinusoidal and in phase. Hence, the converter connection to the supply network is done with a unitary power factor.

Figure 6.

Radiation variation effect: (a) Photovoltaic voltage. (b) DC-link voltage. (c) PFC checking.

5.2 Temperature variation effect

The perfect MPPT is illustrated by Figure 7a. Here, the temperature steps from 298.15 to 318.15 K, then to 308.15 K while the radiation λ is kept constant equal to $1000\mathrm{W}/{\mathrm{m}}^2$. Figure 7b shows that the DC-link voltage $x_4$ is tightly regulated: it quickly settles down after each change in the temperature. Figure 7c illustrates the current amplitude changes whenever the temperature varies. The current remains (almost) sinusoidal and in phase with the network voltage complying with the PFC requirement.

Figure 7.

Temperature variation effect: (a) Photovoltaic voltage. (b) DC-link voltage. (c) PFC checking.

6. Conclusion

In this chapter, an advanced controller is developed for PV grid-connected system. The latter is described by fourth-order nonlinear averaged model. The multi-loops nonlinear controller has been designed and developed using three-time singular perturbation technique and averaging theory.

Using the theoretical analysis (via three-time-scale singular perturbation technique and averaging theory) and simulation, it is proved that the controller does meet the performances for which it was designed, namely: (i) Maximum power point tracking of PV array; (ii) tight regulation of the DC bus voltage; (iii) perfect power factor in the grid; and (iv) global asymptotic stability of the all system.

Several simulation results have been made that illustrate the high performances of the proposed controller in ideal operating conditions (in the presence of meteorological constant) and its robustness against radiation and temperature change.

Appendix (Proof)

Part 1: Eqs. (43)–(45) are immediately obtained from expressions (37), (38), (23), and the first and the third equation of expression (24).

Part 2: The stability of the time-varying system (43)–(45) will now be performed in two steps using the averaging theory (e.g. Chapter 10 in Refs. [19, 25]) and the singular perturbation theory (e.g., Chapter 11 in Refs. [19, 23, 24]). The next step consists in using the averaging theory. To this end, let us introduce the time-scale change $\overline{t}={\omega }_{g}t$. Using this time, one gets ${\epsilon }_1{\epsilon }_2\dot{\overline{Z}}\left(\overline{t}\right)=\rho \overline{h}\left({\overline{X}}^{\prime },\overline{Z},\overline{t}\right)$, ${\epsilon }_2\dot{\overline{Y}}\left(\overline{t}\right)=\rho \overline{g}\left({\overline{X}}^{\prime },\overline{Y},\overline{t}\right)$, and ${\dot{\overline{X}}}^{\prime }\left(\overline{t}\right)=\rho \overline{f}\left({\overline{X}}^{\prime },\overline{Y},\overline{Z},\overline{t}\right)$ with $\rho =1/{\omega }_g$. In view of expressions (43)–(45), it is seem that the functions $\overline{h}\left({\overline{X}}^{\prime },\overline{Z},\overline{t}\right)$, $\overline{g}\left({\overline{X}}^{\prime },\overline{Y},\overline{t}\right)$, and $\overline{f}\left({\overline{X}}^{\prime },\overline{Y},\overline{Z},\overline{t}\right)$ as functions of $\overline{t}$ are periodic with period–2π, let us introduce the averaged functions:

$${\epsilon }_1{\epsilon }_2{\dot{Z}}_0\left(\overline{t}\right)=\rho \underset{t\to \infty }{lim}\frac{1}{2\pi }\int \overline{h}\left({\overline{X}}^{\prime },\overline{Z},\overline{t}\right)d\overline{t}={\mathit{\rho h}}_{\mathit{0}}\left(X_0,Z_{\mathit{0}}\right)$$E47

$${\epsilon }_2{\dot{Y}}_0\left(\overline{t}\right)=\rho \underset{\rho \to 0}{lim}\frac{1}{2\pi }\int \overline{g}\left({\overline{X}}^{\prime },\overline{Y},\overline{t}\right)d\overline{t}={\mathit{\rho g}}_{\mathit{0}}\left(X_0,Y_{\mathit{0}}\right)$$E48

$${\dot{X}}_0\left(\overline{t}\right)=\rho \underset{\rho \to 0}{lim}\frac{1}{2\pi }\int \overline{f}\left({\overline{X}}^{\prime },\overline{Y},\overline{Z},\overline{t}\right)d\overline{t}={\mathit{\rho f}}_{\mathit{0}}\left(X_0,Y_{\mathit{0}},Z_{\mathit{0}}\right)$$E49

Since the systems here studied present equilibrium different from zero and in order to satisfy this requirement, a change of variables is introduced such that defines the new system in terms of its error dynamics. Therefore, the error dynamics are defined by introducing: ${\tilde{Z}}_0={\left(\begin{array}{ll}{\tilde{z}}_{1,0}& {\tilde{z}}_{2,0}\end{array}\right)}^T={\left(\begin{array}{ll}\left(z_{1,0}-z_{1,0}^{\ast }\right)& \left(z_{2,0}-z_{2,0}^{\ast }\right)\end{array}\right)}^T$, ${\tilde{Y}}_0={\left(\begin{array}{ll}{\tilde{y}}_{1,0}& {\tilde{y}}_{2,0}\end{array}\right)}^T={\left(\begin{array}{ll}\left(y_{1,0}-y_{1,0}^{\ast }\right)& \left(y_{2,0}-y_{2,0}^{\ast }\right)\end{array}\right)}^T$, and ${\tilde{X}}_0=\left(\begin{array}{cccc}{\tilde{x}}_{1,0}& {\tilde{x}}_{2,0}& {\tilde{x}}_{3,0}& {\tilde{x}}_{4,0}\end{array}\right)=\left(\begin{array}{cccc}\left(x_{1,0}-x_{1,0}^{\ast }\right)& \left(x_{2,0}-x_{2,0}^{\ast }\right)& \left(x_{3,0}-x_{3,0}^{\ast }\right)& \left(x_{4,0}-x_{4,0}^{\ast }\right)\end{array}\right)$, where the constant $x_{1,0}^{\ast }$, $x_{2,0}^{\ast }$, $x_{3,0}^{\ast }$, $x_{4,0}^{\ast }$, $y_{1,0}^{\ast }$, $y_{2,0}^{\ast }$, $z_{1,0}^{\ast }$, and $z_{2,0}^{\ast }$ represent the desired average values of the state variables. Then, since expressions (43)–(45) and according to expressions (47)–(49), the system can be rewritten into its error-dynamics formulation thus defining the closed-loop error dynamics as:

$${\epsilon }_1{\epsilon }_2{\dot{\tilde{Z}}}_0=\rho {\tilde{h}}_0\left({\tilde{X}}_0,{\tilde{Z}}_0\right)=\rho \left(\begin{array}{c}{k}_1\left(-\frac{{\tilde{x}}_{1,0}}{T_1}-{\tilde{f}}_{1,0}\left({\tilde{X}}_0,{\tilde{Z}}_0\right)\right)\\ k_3\left(-\frac{{\tilde{x}}_{3,0}}{T_3}-{\tilde{f}}_{3,0}\left({\tilde{X}}_0,{\tilde{Z}}_0\right)\right)\end{array}\right)$$E50

$${\epsilon }_2{\dot{\tilde{Y}}}_0=\rho {\tilde{g}}_0\left({\tilde{X}}_0,{\tilde{Y}}_0\right)=\rho \left(\begin{array}{c}{k}_2\left(-\frac{{\tilde{x}}_{2,0}}{T_2}-{\tilde{f}}_{2,0}\left({\tilde{X}}_0,{\tilde{Y}}_0\right)\right)\\ k_4\left(-\frac{{\tilde{x}}_{4,0}}{T_4}-{\tilde{f}}_{4,0}\left({\tilde{X}}_0,{\tilde{Y}}_0\right)\right)\end{array}\right)$$E51

$${\dot{\tilde{X}}}_0=\rho {\tilde{f}}_0\left({\tilde{X}}_0,{\tilde{Y}}_0,{\tilde{Z}}_0\right)=\rho \left(\begin{array}{c}{\tilde{f}}_{1,0}\left({\tilde{X}}_0,{\tilde{Z}}_0\right)\\ {\tilde{f}}_{2,0}\left({\tilde{X}}_0,{\tilde{Y}}_0\right)\\ {\tilde{f}}_{3,0}\left({\tilde{X}}_0,{\tilde{Z}}_0\right)\\ {\tilde{f}}_{4,0}\left({\tilde{X}}_0,{\tilde{Y}}_0\right)\end{array}\right)$$E52

where

$${\tilde{f}}_{1,0}\left({\tilde{X}}_0,{\tilde{Z}}_0\right)=\frac{\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)-\sqrt{\left({\tilde{x}}_{4,0}+x_{4,0}^{\ast }\right)}}{L_i}-\frac{r_i\left({\tilde{x}}_{1,0}+x_{1,0}^{\ast }\right)}{L_i}+\frac{\left({\tilde{z}}_{1,0}+z_{1,0}^{\ast }\right)\sqrt{\left({\tilde{x}}_{4,0}+x_{4,0}^{\ast }\right)}}{L_i}$$E53

$${\tilde{f}}_{2,0}\left({\tilde{X}}_0,{\tilde{Y}}_0\right)=\frac{i_{\mathit{pv},0}}{C_{\mathit{pv}}}-\frac{\left({\tilde{y}}_{1,0}+y_{1,0}^{\ast }\right)}{C_{\mathit{pv}}}$$E54

$${\tilde{f}}_{3,0}\left({\tilde{X}}_0,{\tilde{Z}}_0\right)=-\frac{r_g\left({\tilde{x}}_{3,0}+x_{3,0}^{\ast }\right)}{L_g}+\frac{\sqrt{\left({\tilde{x}}_{4,0}+x_{4,0}^{\ast }\right)}}{L_g}\left({\tilde{z}}_{2,0}+z_{2,0}^{\ast }\right)$$E55

$${\tilde{f}}_{4,0}\left({\tilde{X}}_0,{\tilde{Y}}_0\right)=-\frac{E_g\left({\tilde{y}}_{2,0}+y_{2,0}^{\ast }\right)+r_g{\left({\tilde{y}}_{2,0}+y_{2,0}^{\ast }\right)}^2}{C_{\mathit{dc}}}+2\frac{\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)\left({\tilde{y}}_{1,0}+y_{1,0}^{\ast }\right)-r_i{\left(\left({\tilde{y}}_{1,0}+y_{1,0}^{\ast }\right)\right)}^2}{C_{\mathit{dc}}}$$E56

In order to get stability results regarding the system of interest (43)–(45), it is sufficient (thanks to averaging theory) to analyze the averaged system (50)–(52). Now, the asymptotic stability of the resulting three-time-scale photovoltaic single phase grid system (50)–(52) is discussed, which is based on the sequential (double) time-scale analysis similar to the one presented in Refs. [23, 24], it is an extension of the two-time-scale analysis presented [19].

The use of theory of the three-time-scale singular perturbations for the stability analysis is based on the idea that, for $0<{\epsilon }_{\mathit{1}}{\epsilon }_{\mathit{2}}\ll {\epsilon }_{\mathit{2}}\ll 1$, the trajectories in ${\tilde{X}}_0$, ${\tilde{Y}}_0$, and ${\tilde{Z}}_0$ of the system (50)–(52) can be approximated by three models: the slow dynamic subsystem (SDS) of full system, the fast dynamic subsystem (FDS), and the ultra-fast dynamic subsystem (UFDS). We can thus find Lyapunov functions for each one of the singularly perturbed subsystems.

For the UFDS, it is necessary to ensure that the dynamic of expression (50) does not to shift from the quasi-steady-state equilibrium ${\tilde{Z}}_0={\tilde{h}}_0^{\prime }\left({\tilde{X}}_0\right)$. Then, the associated Lyapunov function candidate is obtained by introducing a change of variables ${\widehat{Z}}_0={\tilde{Z}}_0-{\tilde{h}}_0^{\prime }\left({\tilde{X}}_0\right)$, so that its equilibrium is centered at zero. By letting ${\tau }_1=\overline{t}/{\epsilon }_1{\epsilon }_2$, the UFSD in function of ${\widehat{Z}}_0$ is defined as follows

$$\frac{d{\tilde{Z}}_0}{d{\tau }_1}=\rho {\widehat{h}}_0\left({\tilde{X}}_0,{\widehat{Z}}_0\right)=\rho \left(\begin{array}{c}-\frac{k_1\sqrt{\left({\tilde{x}}_{4,0}+x_{4,0}^{\ast }\right)}}{L_i}{\widehat{z}}_{1,0}\\ -\frac{k_3\sqrt{\left({\tilde{x}}_{4,0}+x_{4,0}^{\ast }\right)}}{L_g}{\widehat{z}}_{2,0}\end{array}\right)$$E57

in which ${\tilde{X}}_0$ is treated as a fixed parameter. Thus, the associated Lyapunov function can be defined by

$$V_U=\frac{1}{2\rho }\left(P_{U1}{\widehat{z}}_{1,0}^2+P_{U2}{\widehat{z}}_{2,0}^2\right)$$E58

In view of $k_i>0\left(i=1,3\right)$ and $\left({\tilde{x}}_{4,0}+x_{4,0}^{\ast }\right)>0$, it is clear that the solutions $P_{U1}$ and $P_{U2}$ can be chosen as

$$P_{U1}=\frac{L_i}{2k_1\sqrt{\left({\tilde{x}}_{4,0}+x_{4,0}^{\ast }\right)}}{q}_{U1}$$E59

$$P_{U2}=\frac{L_g}{2k_3\sqrt{\left({\tilde{x}}_{4,0}+x_{4,0}^{\ast }\right)}}{q}_{U2}$$E60

with $q_{U1}$ and $q_{U2}$ are positive. Similar to the ultra-fast subsystem, it is necessary to ensure that the dynamic of expression (51) does not to shift from the equilibrium ${\tilde{Y}}_0={\tilde{g}}_0^{\prime }\left({\tilde{X}}_0\right)$. Then, by introducing a change of variables ${\widehat{Y}}_0={\tilde{Y}}_0-{\tilde{g}}_0^{\prime }\left({\tilde{X}}_0\right)$, so that its equilibrium is centered at zero and by letting ${\tau }_2=\overline{t}/{\epsilon }_2$, the ultra-fast subsystem can be rewritten in the following form

$$\frac{d{\tilde{Y}}_0}{d{\tau }_2}=\rho {\widehat{g}}_0\left({\tilde{X}}_0,{\widehat{Y}}_0\right)=\rho \left(\begin{array}{c}\frac{k_2}{C_{\mathit{pv}}}{\widehat{y}}_{1,0}\\ {\psi }^{″}\left({\tilde{X}}_0,{\widehat{Y}}_0\right)\end{array}\right)$$E61

with ${\psi }^{″}=2k_4\left[\frac{r_i{\widehat{y}}_{1,0}^2}{C_{\mathit{dc}}}+\frac{r_g{\widehat{y}}_{2,0}^2}{2C_{\mathit{dc}}}+\frac{{\psi }^{\prime }{\widehat{y}}_{2,0}}{C_{\mathit{dc}}}+\left(\frac{2r_i{C}_{\mathit{pv}}}{C_{\mathit{dc}}}\left(\frac{{\tilde{x}}_{2,0}}{T_2}+\frac{i_{\mathit{pv},0}}{C_{\mathit{pv}}}\right)-\frac{\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)}{C_{\mathit{dc}}}\right){\widehat{y}}_{1,0}\right]$.

Since expression (61), it is seem that the fast subsystem is nonlinear. According to the proposition 4 (part i), it is shown that this subsystem can be made asymptotically stable by letting $k_j>0$ (for $j=2,4$). Moreover, we can find a Lyapunov function, which takes the following quadratic form

$$V_F=\frac{1}{2\rho }\left(p_{F1}{\widehat{y}}_{1,0}^2+p_{F2}{\widehat{y}}_{2,0}^2+2p_{F3}{\widehat{y}}_{1,0}{\widehat{y}}_{2,0}\right)$$E62

with the solutions to the associated Lyapunov given as

$$p_{F1}=-\frac{C_{\mathit{pv}}{q}_{F1}}{2k_2}-\frac{2{\left(r_i{C}_{\mathit{pv}}\left(\frac{\tilde{x}}{2,0}+\frac{i_{\mathit{pv},0}}{C_{\mathit{pv}}}\right)-\frac{\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)}{2}\right)}^2{q}_{F2}}{k_2{\psi }^{\prime }\left(\frac{{\psi }^{\prime }}{C_{\mathit{pv}}}+\frac{C_{\mathit{dc}}{k}_2}{2k_4{C}_{\mathit{pv}}^2}\right)}$$E63

$$p_{F2}=-\frac{C_{\mathit{dc}}}{4k_4{\psi }^{\prime }}{q}_{F2}$$E64

$$p_{F3}=\frac{C_{\mathit{dc}}\left(r_i{C}_{\mathit{pv}}\left(\frac{\tilde{x}}{2,0}+\frac{i_{\mathit{pv},0}}{C_{\mathit{pv}}}\right)-\frac{\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)}{2}\right)}{k_4{\psi }^{\prime }\left(2{\psi }^{\prime }+\frac{C_{\mathit{dc}}{k}_2}{k_4{C}_{\mathit{pv}}}\right)}{q}_{F2}$$E65

in which ${\tilde{X}}_0$ is treated as a fixed parameter. Finally, the new slow dynamic subsystem is obtained by substituting of the ultra-fast subsystem equilibrium ${\tilde{Z}}_0={\tilde{h}}_0^{\prime }\left({\tilde{X}}_0\right)$ and the fast subsystem equilibrium ${\tilde{Y}}_0={\tilde{g}}_0^{\prime }\left({\tilde{X}}_0\right)$ into expression (52)

$${\dot{\tilde{X}}}_0={\left(\begin{array}{cccc}{\dot{\tilde{x}}}_{1,0}& {\dot{\tilde{x}}}_{2,0}& {\dot{\tilde{x}}}_{3,0}& {\dot{\tilde{x}}}_{4,0}\end{array}\right)}^T=-\rho {\left(\begin{array}{cccc}{\tilde{x}}_{1,0}/T_1& {\tilde{x}}_{2,0}/T_2& {\tilde{x}}_{3,0}/T_3& {\tilde{x}}_{4,0}/T_4\end{array}\right)}^T$$E66

It is easy to define the associated Lyapunov function for SDS as follows

$$V_S\left({\tilde{X}}_0\right)=\frac{1}{2\rho }{\tilde{X}}_0^T{P}_s{\tilde{X}}_0=\frac{1}{2\rho }\left(p_{S1}{x}_{1,0}^2+p_{S2}{x}_{2,0}^2+p_{S3}{x}_{3,0}^2+p_{S4}{x}_{4,0}^2\right)$$E67

where $p_{S1}=T_1{q}_{S1}/2$, $p_{S2}=T_2{q}_{S2}/2$, $p_{S3}=T_3{q}_{S3}/2$, $p_{S4}=T_4{q}_{S4}/2$, and $q_{S1}$, $q_{S2}$, $q_{S3}$, $q_{S4}$ are the positive constants.

Based on these Lyapunov function candidates, the double application of the standard two-time-scale stability analysis is divided in two stages: in the first stage, the stability analysis focusses on proving the stability properties of the degenerated ${\sum }_{\mathit{SF}}$-subsystem (slow-fast subsystem). The results obtained will be used in the second stage in order to prove the stability properties for the full ${\sum }_{\mathit{SFU}}$-system (slow-fast-ultra-fast subsystem).

$\sum \mathit{SF}$Stability analysis

In the first stage, the standard method for two-time-scale systems is applied in which the previously derived Lyapunov functions for the slow and fast subsystems, that is, $V_S$ and $V_F$, respectively, must satisfy certain inequalities.

• Isolated equilibrium of the origin for the${\sum }_{\mathit{SF}}$-subsystem

The origin $\left({\tilde{X}}_0=0,{\tilde{Y}}_0=0\right)$ is an isolated equilibrium of the ${\sum }_{\mathit{SF}}$-subsystem

$$0={\widehat{g}}_0\left(0,0\right)$$E68

$$0={\tilde{f}}_0\left(0,0,{\tilde{h}}_0\left({\tilde{X}}_0\right)\right)$$E69

moreover, ${\tilde{Y}}_0={\tilde{g}}_0\left({\tilde{X}}_0\right)$ is the root of $0={\widehat{g}}_0\left({\tilde{X}}_0,{\tilde{Y}}_0\right)$ which vanishes at ${\tilde{X}}_0=0$, and $‖{\tilde{g}}_0\left({\tilde{X}}_0\right)‖<{\vartheta }_1\left(‖{\tilde{X}}_0‖\right)$ where ${\vartheta }_1\left(\right)$ is a $\kappa$ function.

• Reduced system condition for the${\sum }_{\mathit{SF}}$-subsystem

Using the Lyapunov function (67) and substituting $\rho {\tilde{f}}_0\left({\tilde{X}}_0,{\tilde{Y}}_0,{\tilde{h}}_0\left({\tilde{X}}_0\right)\right)$ yields

$${\left(\frac{\partial V_S\left({\tilde{X}}_0\right)}{\partial {\tilde{X}}_0}\right)}^T\rho {\tilde{f}}_0\left({\tilde{X}}_0,{\tilde{Y}}_0,{\tilde{h}}_0\left({\tilde{X}}_0\right)\right)\le -{\alpha }_1{\varphi }_1^2\left({\tilde{X}}_0\right)$$E70

where ${\alpha }_1\le 1$ and ${\varphi }_1\left({\tilde{X}}_0\right)=\sqrt{\frac{q_{S1}}{2}{\tilde{x}}_{1,0}^2+\frac{q_{S2}}{2}{\tilde{x}}_{2,0}^2+\frac{q_{S3}}{2}{\tilde{x}}_{3,0}^2+\frac{q_{S4}}{2}{\tilde{x}}_{4,0}^2}$.

• Boundary-layer system condition for the${\sum }_{\mathit{SF}}$-subsystem

Since expression (62)–(65), it is seem that

$${\left(\frac{\partial V_{\mathit{UF}}\left({\tilde{X}}_0,{\widetilde{Ƶ}}_0\right)}{\partial {\widetilde{Ƶ}}_0}\right)}^T\rho \widehat{H}\left({\tilde{X}}_0,{\widetilde{Ƶ}}_0\right)\le -Q_{F2}{\widehat{y}}_{2,0}^2-Q_{F1}{\widehat{y}}_{1,0}^2$$E71

where $Q_{F1}$ and $Q_{F2}$ are defined as follows

$$Q_{F1}=\frac{q_{F1}}{2}-\frac{r_i{q}_{F2}}{2{\psi }^{\prime }}\left(\left|\frac{4\left(r_i{C}_{\mathit{pv}}\left(\frac{\tilde{x}}{2,0}+\frac{i_{\mathit{pv},0}}{C_{\mathit{pv}}}\right)-\frac{\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)}{2}\right)}{\left(2{\psi }^{\prime }+\frac{C_{\mathit{dc}}{k}_2}{k_4{C}_{\mathit{pv}}}\right)}{\widehat{y}}_{1,0}\right|+\left|{\widehat{y}}_{2,0}\right|\right)$$E72

$$Q_{F2}=\frac{q_{F2}}{2}-\frac{r_g{q}_{F2}}{4{\psi }^{\prime }}\left(\left|\frac{4\left(r_i{C}_{\mathit{pv}}\left(\frac{\tilde{x}}{2,0}+\frac{i_{\mathit{pv},0}}{C_{\mathit{pv}}}\right)-\frac{\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)}{2}\right)}{\left(2{\psi }^{\prime }+\frac{C_{\mathit{dc}}{k}_2}{k_4{C}_{\mathit{pv}}}\right)}{\widehat{y}}_{1,0}\right|+\left|{\widehat{y}}_{2,0}\right|\right)$$E73

For physical point of view and domain of working principle, it is supposed that all physical state variables are bounded in domain of interest, where ${\tilde{X}}_{0min}\le {\tilde{X}}_0\le {\tilde{X}}_{0\mathrm{Max}}$, ${\tilde{Y}}_{0min}\le {\tilde{Y}}_0\le {\tilde{Y}}_{0\mathrm{Max}}$, and ${\tilde{Z}}_{0min}\le {\tilde{Z}}_0\le {\tilde{Z}}_{0\mathrm{Max}}$. Therefore, $Q_{F1}$ and $Q_{F2}$ can take positives minimum possible values $Q_{F1,min}>0$ and $Q_{F2,min}>0$. This can be done by ensuring the appropriate selection of $q_{F1}$, $q_{F2}$, $k_2$, and $k_4$ as follows

$$\left\{\begin{array}{l}\frac{k_2}{k_4}>\frac{C_{\mathit{pv}}}{C_{\mathit{dc}}}\left(\frac{4\left|r_i{C}_{\mathit{pv}}\left(\frac{\tilde{x}}{2,0\mathit{Max}}+\frac{i_{\mathit{pv},0}}{C_{\mathit{pv}}}\right)-\frac{\left({\tilde{x}}_{2,0\mathit{Max}}+x_{2,0}^{\ast }\right)}{2}\right|}{\left(\frac{2{\psi }_{\min }^{\prime }}{r_g}-{\widehat{y}}_{2,0\mathrm{Max}}\right)}{\widehat{y}}_{1,0\mathrm{Max}}-2{\psi }_{\min }^{\prime }\right)\\ \frac{q_{F1}}{q_{F2}}>\frac{r_i}{{\psi }_{\min }^{\prime }}\left(\frac{4\left|\left(r_i{C}_{\mathit{pv}}\left(\frac{\tilde{x}}{2,0\mathit{Max}}+\frac{i_{\mathit{pv},0}}{C_{\mathit{pv}}}\right)-\frac{\left({\tilde{x}}_{2,0\mathit{Max}}+x_{2,0}^{\ast }\right)}{2}\right)\right|}{\left(2{\psi }_{\min }^{\prime }+\frac{C_{\mathit{dc}}{k}_2}{k_4{C}_{\mathit{pv}}}\right)}{\widehat{y}}_{1,0\mathrm{Max}}+{\widehat{y}}_{2,0\mathrm{Max}}\right)\end{array}\right.$$E74

Therefore, the boundary-layer system condition for the ${\sum }_{\mathit{SF}}$-subsystem is defined as follows

$${\left(\frac{\partial V_F}{\partial {\tilde{Y}}_0}\right)}^T\rho {\widehat{g}}_0\left({\tilde{X}}_0,{\tilde{Y}}_0\right)\le -{\left(\sqrt{Q_{F2min}{\widehat{y}}_{2,0}^2+Q_{F\mathit{1}min}{\widehat{y}}_{1,0}^2}\right)}^2\le -{\alpha }_2{\varphi }_2^2\left({\widehat{Y}}_0\right)$$E75

where ${\alpha }_2\le 1$ and ${\varphi }_2\left({\widehat{Y}}_0\right)=\sqrt{Q_{F2min}{\widehat{y}}_{2,0}^2+Q_{F1min}{\widehat{y}}_{1,0}^2}$.

Now, we consider the composite Lyapunov function candidate of the ${\sum }_{\mathit{SF}}$ -subsystem given as follow

$$V_{\mathrm{SF}}=\left(1-{\eta }_1\right)V_{\mathrm{S}}\left({\tilde{X}}_0\right)+{\eta }_1{V}_{\mathrm{F}}\left({\widehat{Y}}_0\right)$$E76

with $0<{\eta }_1<1$. The derivative of $V_{\mathrm{SF}}$ presents new terms which represents the effect of the interconnection between the slow and fast dynamics. These interconnections are assumed to satisfy the following conditions

$${\left(\frac{\partial V_s\left({\tilde{X}}_0\right)}{\partial {\tilde{X}}_0}\right)}^T\rho \left({\tilde{f}}_0\left({\tilde{X}}_0,{\tilde{Y}}_0,{\tilde{h}}_0\left({\tilde{X}}_0\right)\right)-{\tilde{f}}_0\left({\tilde{X}}_0,{\tilde{g}}_0\left({\tilde{X}}_0\right),{\tilde{h}}_0\left({\tilde{X}}_0\right)\right)\right)\le {\theta }_1{\varphi }_2\left({\widehat{Y}}_0\right){\varphi }_1\left({\tilde{X}}_0\right)$$E77

$${\left(\frac{\partial V_F\left({\widehat{Y}}_0\right)}{\partial {\tilde{X}}_0}\right)}^T\rho {\tilde{f}}_0\left({\tilde{X}}_0,{\tilde{g}}_0\left({\tilde{X}}_0\right),{\tilde{h}}_0\left({\tilde{X}}_0\right)\right)\le {\theta }_2{\varphi }_2\left({\widehat{Y}}_0\right){\varphi }_1\left({\tilde{X}}_0\right)$$E78

$${\left(\frac{\partial V_F\left({\widehat{Y}}_0\right)}{\partial {\tilde{X}}_0}\right)}^T\left[\rho {\tilde{f}}_0\left({\tilde{X}}_0,{\tilde{Y}}_0,{\tilde{h}}_0\left({\tilde{X}}_0\right)\right)-\rho {\tilde{f}}_0\left({\tilde{X}}_0,{\tilde{g}}_0\left({\tilde{X}}_0\right),{\tilde{h}}_0\left({\tilde{X}}_0\right)\right)\right]\le {\gamma }_1{\varphi }_2^2\left({\widehat{Y}}_0\right)$$E79

where the constants ${\theta }_1$, ${\theta }_2$, and ${\gamma }_1$ are non-negative.

Therefore, from the singular perturbation theory (e.g. Theorem 11.3 in Ref. [19]), it follows that the derivative of V_SF is negative-definite for

$${\epsilon }_2<{\epsilon }_2^{\otimes }\mathit{\text{for}}\left({\eta }_1<{\eta }_1^{\otimes }\right)$$E80

$$\text{with}{\epsilon }_2^{\otimes }=\frac{{\alpha }_1{\alpha }_2}{{\alpha }_1{\gamma }_1+{\theta }_1{\theta }_2}\text{and}{\eta }_1^{\otimes }=\frac{{\theta }_1}{{\theta }_1+{\theta }_2}.$$E81

$\sum \mathit{SFU}$ Stability Analysis:

The stability of full system is analyzed now. In this step, the results obtained in above section will be used in order to prove the asymptotic stability properties of the full ${\sum }_{\mathit{SFU}}$-system, which, for convenience, is first rewritten as

$${\dot{\tilde{\chi }}}_0={\tilde{F}}_0\left({\tilde{\chi }}_0,{\tilde{Z}}_0\right)$$E82

$${\epsilon }_1{\epsilon }_2{\dot{\tilde{Z}}}_0={\widehat{h}}_0\left({\tilde{\chi }}_0,{\tilde{Z}}_0\right)$$E83

where ${\tilde{\chi }}_0={\left(\begin{array}{cccccc}{\tilde{x}}_{1,0}& {\tilde{x}}_{2,0}& {\tilde{x}}_{3,0}& {\tilde{x}}_{4,0}& {\tilde{y}}_{1,0}& {\tilde{y}}_{2,0}\end{array}\right)}^T$, similarly to the ${\sum }_{\mathit{SF}}$ general asymptotic stability analysis presented in the above section, the ${\sum }_{\mathit{SFU}}$-system is treated like a two-time-scale singularly perturbed system, where the ${\sum }_{\mathit{SF}}$-subsystem is treated as the new slow augmented reduced order. The stability of the full system implicates that the previously derived Lyapunov functions for the new slow and new fast subsystems, that is, $V_{\mathrm{SF}}\left({\tilde{\chi }}_0\right)$ and $V_U\left({\widehat{Z}}_0\right)$, respectively, must satisfy certain inequalities.

• Isolated equilibrium of the origin for the${\sum }_{\mathit{SFU}}$-subsystem

The origin $\left({\tilde{\chi }}_0=0,{\tilde{Z}}_0=0\right)$ is an isolated equilibrium of the ${\sum }_{\mathit{SFU}}$-subsystem

$$0={\tilde{F}}_0\left(0,0\right)$$E84

$$0={\widehat{h}}_0\left(0,0\right)$$E85

Moreover, ${\tilde{Z}}_0={\tilde{h}}_0\left({\tilde{\chi }}_0\right)$ is the unique root of $0={\widehat{h}}_0\left({\tilde{\chi }}_0,{\tilde{Z}}_0\right)$, which vanishes at ${\tilde{\chi }}_0=0$, and $‖{\tilde{h}}_0\left({\tilde{\chi }}_0\right)‖<{\vartheta }_{\mathit{2}}\left(‖{\tilde{\chi }}_0‖\right)$ where ${\vartheta }_2\left(\right)$ is a $\kappa$ function.

• Reduced system condition for the${\sum }_{\mathit{SFU}}$-subsystem

Using Lyapunov function candidate $V_{\mathrm{SF}}\left({\tilde{\chi }}_0\right)$ given by expression (76), it is easily shown that

$${\left(\frac{\partial V_{\mathrm{SF}}\left({\tilde{\chi }}_0\right)}{\partial {\tilde{\chi }}_0}\right)}^T={\left(\begin{array}{cccccc}\frac{\partial V_{\mathrm{SF}}\left({\tilde{\chi }}_0\right)}{\partial {\tilde{x}}_{1,0}}& \frac{\partial V_{\mathrm{SF}}\left({\tilde{\chi }}_0\right)}{\partial {\tilde{x}}_{2,0}}& \frac{\partial V_{\mathrm{SF}}\left({\tilde{\chi }}_0\right)}{\partial {\tilde{x}}_{3,0}}& \frac{\partial V_{\mathrm{SF}}\left({\tilde{\chi }}_0\right)}{\partial {\tilde{x}}_{4,0}}& \frac{\partial V_{\mathrm{SF}}\left({\tilde{\chi }}_0\right)}{\partial {\tilde{y}}_{1,0}}& \frac{\partial V_{\mathrm{SF}}\left({\tilde{\chi }}_0\right)}{\partial {\tilde{y}}_{2,0}}\end{array}\right)}^T$$E86

Recalling that ${\widehat{Y}}_0={\tilde{Y}}_0-{\tilde{g}}_0\left({\tilde{X}}_0\right)$, and since (61)–(65), one gets

$$\begin{array}{l}\frac{\partial p_{F1}}{\partial {\tilde{x}}_{2,0}}=\left[\frac{r_g{C}_{\mathit{pv}}{\left(2r_i{C}_{\mathit{pv}}\left(\frac{\tilde{x}}{T_2}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}\right)-\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)\right)}^2\left(\frac{2C_{\mathit{pv}}}{T_2}\left(1-\frac{r_i{C}_{\mathit{pv}}}{T_2}\right){\tilde{x}}_2+C_{\mathit{pv}}\left(\frac{x_{2,0}^{\ast }}{T_2}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}-\frac{2r_i{i}_{\mathit{pv}}}{T_2}\right)\right)}{k_2{{\psi }^{\prime }}^3{\left({\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)}^2}\right.\\ \left.\times \left(2{\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)-\frac{2C_{\mathit{pv}}\left(2r_i{C}_{\mathit{pv}}\left(\frac{{\tilde{x}}_{2,0}}{T_2}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}\right)-\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)\right)\left(\frac{2r_i{C}_{\mathit{pv}}}{T_2}-1\right)}{k_2{\psi }^{\prime }\left({\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)}\right]q_{F2}=A_1{q}_{F2}\end{array}$$

$$\frac{\partial p_{F1}}{\partial {\tilde{x}}_{4,0}}=\frac{r_g{C}_{\mathit{dc}}{C}_{\mathit{pv}}{\left(2r_i{C}_{\mathit{pv}}\left(\frac{{\tilde{x}}_{2,0}}{T_2}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}\right)-\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)\right)}^2\left(2{\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)}{2k_2{T}_4{{\psi }^{\prime }}^3{\left({\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)}^2}{q}_{F2}=A_2{q}_{F2}$$

$$\frac{\partial p_{F2}}{\partial {\tilde{x}}_{2,0}}=\frac{C_{\mathit{dc}}{C}_{\mathit{pv}}{r}_g\left(\frac{2}{T_2}\left(1-\frac{r_i{C}_{\mathit{pv}}}{T_2}\right){\tilde{x}}_{2,0}+\left(\frac{x_{2,0}^{\ast }}{T_2}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}-\frac{2r_i{i}_{\mathit{pv}}}{T_2}\right)\right)}{4k_4{{\psi }^{\prime }}^3}{q}_{F2}=A_3{q}_{F2},$$

$$\frac{\partial p_{F2}}{\partial {\tilde{x}}_{4,0}}=\frac{r_g{C}_{\mathit{dc}}^2}{4T_4{k}_4{{\psi }^{\prime }}^3}{q}_{F2}=A_4{q}_{F2},$$

$$\frac{\partial p_{F3}}{\partial {\tilde{x}}_{4,0}}=-\frac{C_{\mathit{dc}}^2{r}_g\left(2r_i{C}_{\mathit{pv}}\left(\frac{{\tilde{x}}_{2,0}}{T_2}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}\right)-\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)\right)\left(2{\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)}{4k_4{T}_4{{\psi }^{\prime }}^3{\left({\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)}^2}{q}_{F2}=A_6{q}_{F2},$$

$$\begin{array}{l}\frac{\partial p_{F3}}{\partial {\tilde{x}}_{2,0}}=\left[-\frac{C_{\mathit{dc}}{r}_g\left(2r_i{C}_{\mathit{pv}}\left(\frac{{\tilde{x}}_{2,0}}{T_2}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}\right)-\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)\right)\left(2\frac{C_{\mathit{pv}}}{T_2}\left(1-\frac{r_i{C}_{\mathit{pv}}}{T_2}\right){\tilde{x}}_{2,0}+C_{\mathit{pv}}\left(\frac{x_{2,0}^{\ast }}{T_2}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}-\frac{2r_i{i}_{\mathit{pv}}}{T_2}\right)\right)}{4k_4{{\psi }^{\prime }}^3\left({\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)}\right.\\ \left.\times \left(2{\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)+\frac{\left(\frac{2r_i{C}_{\mathit{pv}}}{T_2}-1\right)C_{\mathit{dc}}}{4k_4{\psi }^{\prime }\left({\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)}\right]q_{F2}=A_5{q}_{F2}\end{array}$$

$$p_{F1}\frac{\partial {\widehat{y}}_{1,0}}{\partial {\tilde{x}}_{2,0}}=\frac{C_{\mathit{pv}}^2}{2k_2{T}_2}{q}_{F1}+\frac{C_{\mathit{pv}}^2{\left(2r_i{C}_{\mathit{pv}}\left(\frac{{\tilde{x}}_{2,0}}{T_2}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}\right)-\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)\right)}^2}{2k_2{T}_2{\psi }^{\prime }\left({\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)}{q}_{F2}=A_{13}{q}_{F1}+A_7{q}_{F2}$$

$$p_{F2}\frac{\partial {\widehat{y}}_{2,0}}{\partial {\tilde{x}}_{2,0}}=\frac{C_{\mathit{dc}}\left(2\frac{C_{\mathit{pv}}}{T_2}\left(1-\frac{r_i{C}_{\mathit{pv}}}{T_2}\right){\tilde{x}}_{2,0}+C_{\mathit{pv}}\left(\frac{x_{2,0}^{\ast }}{T_2}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}-\frac{2r_i{i}_{\mathit{pv}}}{T_2}\right)\right)}{4k_4{{\psi }^{\prime }}^2}{q}_{F2}=A_8{q}_{F2}$$

$p_{F3}\frac{\partial {\widehat{y}}_{1,0}}{\partial {\tilde{x}}_{2,0}}=-\frac{C_{\mathit{pv}}{C}_{\mathit{dc}}\left(2r_i{C}_{\mathit{pv}}\left(\frac{{\tilde{x}}_{2,0}}{T_2}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}\right)-\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)\right)}{4T_2{k}_4{\psi }^{\prime }\left({\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)}{q}_{F2}=A_9{q}_{F2}$, $p_{F3}\frac{\partial {\widehat{y}}_{1,0}}{\partial {\tilde{x}}_{4,0}}=0$

$$p_{F3}\frac{\partial {\widehat{y}}_{2,0}}{\partial {\tilde{x}}_{2,0}}=-\frac{C_{\mathit{dc}}\left(2r_i{C}_{\mathit{pv}}\left(\frac{{\tilde{x}}_{2,0}}{T_2}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}\right)-\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)\right)\left(2\frac{C_{\mathit{pv}}}{T_2}\left(1-\frac{r_i{C}_{\mathit{pv}}}{T_2}\right){\tilde{x}}_{2,0}+C_{\mathit{pv}}\left(\frac{x_{2,0}^{\ast }}{T_2}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}-\frac{2r_i{i}_{\mathit{pv}}}{T_2}\right)\right)}{4k_4{{\psi }^{\prime }}^2\left({\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)}{q}_{F2}=A_{10}{q}_{F2}$$

$p_{F1}\frac{\partial {\widehat{y}}_{1,0}}{\partial {\tilde{x}}_{4,0}}=0$, $p_{F3}\frac{\partial {\widehat{y}}_{2,0}}{\partial {\tilde{x}}_{4,0}}=-\frac{C_{\mathit{dc}}^2\left(2r_i{C}_{\mathit{pv}}\left(\frac{{\tilde{x}}_{2,0}}{T_2}+\frac{i_{\mathit{pv}}}{C_{\mathit{pv}}}\right)-\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)\right)}{4T_4{k}_4{{\psi }^{\prime }}^2\left({\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)}{q}_{F2}=A_{12}{q}_{F2}$, $p_{F2}\frac{\partial {\widehat{y}}_{2,0}}{\partial {\tilde{x}}_{4,0}}=\frac{C_{\mathit{dc}}^2}{4k_4{T}_4{{\psi }^{\prime }}^2}{q}_{F2}=A_{11}{q}_{F2}$.

Substituting the expressions of $\rho {\tilde{F}}_0\left({\tilde{\chi }}_0,{\tilde{h}}_0\left({\tilde{\chi }}_0\right)\right)$ and by using absolute value version of Young’s inequality $\left|\mathit{ab}\right|\le \left|\left(a^2+b^{\mathit{2}}\right)/2\right|$, it shows that

$$\begin{array}{l}{\left(\frac{\partial V_1\left({\tilde{\chi }}_0\right)}{\partial {\tilde{\chi }}_0}\right)}^T\rho F_0\left({\tilde{\chi }}_0,{\tilde{h}}_0\left({\tilde{\chi }}_0\right)\right)\le -{\mathrm{{\rm M}}}_1{\tilde{x}}_{1,0}^2-{\mathrm{{\rm M}}}_2{\tilde{x}}_{2,0}^2-{\mathrm{{\rm M}}}_3{\tilde{x}}_{3,0}^2\\ -{\mathrm{{\rm M}}}_4{\tilde{x}}_{4,0}^2-{\mathrm{{\rm M}}}_5{\widehat{y}}_{1,0}^2-{\mathrm{{\rm M}}}_6{\widehat{y}}_{2,0}^2\end{array}$$E87

where ${\mathrm{{\rm M}}}_1=\left(1-{\eta }_1\right)\frac{q_{S1}}{2}$, ${\mathrm{{\rm M}}}_2=\left(1-{\eta }_1\right)\frac{q_{S2}}{2}-\frac{\left(1-{\eta }_1\right)}{2}{\left(\frac{T_2}{2}{q}_{S2}\right)}^2-\frac{{\eta }_1}{2}{\left(-\frac{1}{T_2}\right)}^2$, ${\mathrm{{\rm M}}}_3=-\left(1-{\eta }_1\right)\frac{q_{S3}}{2}$, ${\mathrm{{\rm M}}}_4=\left(1-{\eta }_1\right)\frac{q_{S4}}{2}-\frac{\left(1-{\eta }_1\right)}{2}{\left(\frac{T_4}{2}{q}_{S4}\right)}^2-\frac{{\eta }_1}{2}{\left(-\frac{1}{T_4}\right)}^2$, $M_5=\frac{q_{F1}}{2}-M_5^{\prime }-\frac{\left(1-{\eta }_1\right)}{2C_{\mathit{pv}}^2}$, and $M_6=\frac{q_{F2}}{2}-M_6^{\prime }-\frac{{\eta }_1}{2}{\left(-\frac{2{\psi }^{\prime }}{C_{\mathit{dc}}}-\frac{r_g}{C_{\mathit{dc}}}{\widehat{y}}_{2,0}\right)}_{\mathit{Max}}^2$

$$\begin{array}{l}{M}_5^{\prime }={\eta }_1{\left|\frac{A_{13}{q}_{F1}+\left(\frac{A_1}{2}{\widehat{y}}_{1,0}+A_5{\widehat{y}}_{2,0}+\left(A_7+A_{10}\right)\right)q_{F2}}{C_{\mathit{pv}}}\right|}_{\mathit{Max}}+{\eta }_1{q}_{F2}{\left|\frac{\frac{A_3}{2}{\widehat{y}}_{2,0}+A_8+A_9}{2C_{\mathit{pv}}}\right|}_{\mathit{Max}}\\ +{\eta }_1\frac{q_{F2}}{2}{\left(\frac{2\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)}{C_{\mathit{dc}}}+\frac{2r_i}{C_{\mathit{dc}}}{\widehat{y}}_{1,0}\right)}_{\mathit{Max}}^2+\frac{{\eta }_1}{2}\left(\left(\frac{A_1}{2}{\widehat{y}}_{1,0}+A_5{\widehat{y}}_{2,0}+\left(A_7+A_{10}\right)\right)q_{F2}+\right.\\ {\left.A_{13}{q}_{F1}\right)}_{\mathit{Max}}^2+{\eta }_1{q}_{F2}{\left|\left(\frac{A_2}{2}{\widehat{y}}_{1,0}+A_6{\widehat{y}}_{2,0}+A_{12}\right)\left(\frac{2\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)}{C_{\mathit{dc}}}+\frac{2r_i}{C_{\mathit{dc}}}{\widehat{y}}_{1,0}\right)\right|}_{\mathit{Max}}\end{array}$$

$$\begin{array}{l}+{\eta }_1\frac{q_{F2}}{2}{\left|\left(\frac{A_2}{2}{\widehat{y}}_{1,0}+A_6{\widehat{y}}_{2,0}+A_{12}\right)\left(-\frac{2{\psi }^{\prime }}{C_{\mathit{dc}}}-\frac{r_g}{C_{\mathit{dc}}}{\widehat{y}}_{2,0}\right)\right|}_{\mathit{Max}}+{\eta }_1{q}_{F2}\\ \times {\left|\left(\frac{A_4}{2}{\widehat{y}}_{2,0}+A_{11}\right)\left(\frac{2\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)}{C_{\mathit{dc}}}+\frac{2r_i}{C_{\mathit{dc}}}{\widehat{y}}_{1,0}\right)\right|}_{\mathit{Max}}+{\eta }_1{\left(\left(\frac{A_2}{2}{\widehat{y}}_{1,0}+A_6{\widehat{y}}_{2,0}+A_{12}\right)q_{F2}\right)}_{\mathit{Max}}^2\\ +{\eta }_1{q}_{F2}{\left|-\frac{r_i}{2{\psi }^{\prime }}{\widehat{y}}_{2,0}+\frac{\left(2r_i{C}_{\mathit{pv}}\left(\frac{{\tilde{x}}_{2,0}}{T_2}+\frac{i_{\mathit{pv},0}}{C_{\mathit{pv}}}\right)-\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)\right)r_i}{2{\psi }^{\prime }\left({\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)}{\widehat{y}}_{1,0}\right|}_{\mathit{Max}}\end{array}$$

$$\begin{array}{l}{M}_6^{\prime }={\eta }_1\frac{q_{F2}}{2}{\left|\frac{\frac{A_3}{2}{\widehat{y}}_{2,0}+A_8+A_9}{2C_{\mathit{pv}}}\right|}_{\mathit{Max}}+\frac{{\eta }_1}{2}{\left(\left(\frac{A_3}{2}{\widehat{y}}_{2,0}+A_8+A_9\right)q_{F2}\right)}_{\mathit{Max}}^2\\ +\frac{{\eta }_1}{2}{\left(\left(+\frac{A_4}{2}{\widehat{y}}_{2,0}+A_{11}\right)q_{F2}\right)}_{\mathit{Max}}^2+{\eta }_1{q}_{F2}{\left|\left(\frac{A_4}{2}{\widehat{y}}_{2,0}+A_{11}\right)\left(-\frac{2{\psi }^{\prime }}{C_{\mathit{dc}}}-\frac{r_g}{C_{\mathit{dc}}}{\widehat{y}}_{2,0}\right)\right|}_{\mathit{Max}}\\ +{\eta }_1\frac{q_{F2}}{2}{\left|\left(\frac{A_2}{2}{\widehat{y}}_{1,0}+A_6{\widehat{y}}_{2,0}+A_{12}\right)\left(-\frac{2{\psi }^{\prime }}{C_{\mathit{dc}}}-\frac{r_g}{C_{\mathit{dc}}}{\widehat{y}}_{2,0}\right)\right|}_{\mathit{Max}}\\ +{\eta }_1\frac{q_{F2}}{2}{\left|\left(+\frac{A_4}{2}{\widehat{y}}_{2,0}+A_{11}\right)\left(\frac{2\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)}{C_{\mathit{dc}}}+\frac{2r_i}{C_{\mathit{dc}}}{\widehat{y}}_{1,0}\right)\right|}_{\mathit{Max}}\\ +{\eta }_1{q}_{F2}{\left|-\frac{r_g}{4{\psi }^{\prime }}{\widehat{y}}_{2,0}+\frac{r_g\left(2r_i{C}_{\mathit{pv}}\left(\frac{{\tilde{x}}_{2,0}}{T_2}+\frac{i_{\mathit{pv},0}}{C_{\mathit{pv}}}\right)-\left({\tilde{x}}_{2,0}+x_{2,0}^{\ast }\right)\right)}{4{\psi }^{\prime }\left({\psi }^{\prime }+\left(C_{\mathit{dc}}{k}_2/2k_4{C}_{\mathit{pv}}\right)\right)}{\widehat{y}}_{1,0}\right|}_{\mathit{Max}}\end{array}$$

Now, it needs to ensure that $-{\mathrm{{\rm M}}}_1{\tilde{x}}_{1,0}^2-{\mathrm{{\rm M}}}_2{\tilde{x}}_{2,0}^2-{\mathrm{{\rm M}}}_3{\tilde{x}}_{3,0}^2-{\mathrm{{\rm M}}}_4{\tilde{x}}_{4,0}^2-{\mathrm{{\rm M}}}_5{\widehat{y}}_{1,0}^2-{\mathrm{{\rm M}}}_6{\widehat{y}}_{2,0}^2<0$. This can be done by ensuring the appropriate selection of $q_{S1}$, $q_{S2}$, $q_{S3}$, $q_{S4}$, $q_{\mathrm{F}1}$, and $q_{\mathrm{F}2}$, such that $q_{S1}$ and $q_{S3}$ have sufficient large positive values, $q_{\mathrm{F}1}$ and $q_{\mathrm{F}2}$ have sufficient small positive values, and $q_{S2}$ and $q_{S4}$ are limited as follows $\frac{1-\sqrt{1-\frac{{\eta }_1}{\left(1-{\eta }_1\right)}}}{\frac{T_2^2}{2}}<q_{S2}<\frac{1+\sqrt{1-\frac{{\eta }_1}{\left(1-{\eta }_1\right)}}}{\frac{T_2^2}{2}}$, $\frac{1-\sqrt{1-\frac{{\eta }_1}{\left(1-{\eta }_1\right)}}}{\frac{T_4^2}{2}}<q_{S4}<\frac{1+\sqrt{1-\frac{{\eta }_1}{\left(1-{\eta }_1\right)}}}{\frac{T_4^2}{2}}$ (for ${\eta }_1<\frac{1}{3}$). Therefore, the reduced system condition for the ${\sum }_{\mathit{SFU}}$ system is defined as follows

$${\left(\frac{\partial V_1\left({\tilde{\chi }}_0\right)}{\partial {\tilde{\chi }}_0}\right)}^T\rho F_0\left({\tilde{\chi }}_0,{\tilde{h}}_0\left({\tilde{\chi }}_0\right)\right)\le -{\alpha }_3{\varphi }_3^2\left({\tilde{\chi }}_0\right)$$E88

where ${\alpha }_3\le 1$ and ${\varphi }_3\left({\tilde{\chi }}_0\right)=\sqrt{{\mathrm{{\rm M}}}_1{\tilde{x}}_{1,0}^2+{\mathrm{{\rm M}}}_2{\tilde{x}}_{2,0}^2+{\mathrm{{\rm M}}}_3{\tilde{x}}_{3,0}^2+{\mathrm{{\rm M}}}_4{\tilde{x}}_{4,0}^2+{\mathrm{{\rm M}}}_5{\widehat{y}}_{1,0}^2+{\mathrm{{\rm M}}}_6{\widehat{y}}_{2,0}^2}$.

• Boundary-layer system condition for the ${\sum }_{\mathit{SFU}}$-system

Using the Lyapunov function given by expression (58)–(60), and substituting $\rho {\widehat{h}}_0\left({\tilde{X}}_0,{\widehat{Z}}_0\right)$ of expression (57), one obtains

$$\left(\frac{V_u\left({\tilde{\chi }}_0,{\tilde{Z}}_0\right)}{\partial {\tilde{Z}}_0}\right)\rho {\widehat{h}}_0\left({\tilde{\chi }}_0,{\tilde{Z}}_0\right)=\frac{q_{U1}}{2}{\widehat{z}}_{1,0}^2+\frac{q_{U2}}{2}{\widehat{z}}_{2,0}^2\le -{\alpha }_4{\varphi }_4^2\left({\widehat{Z}}_0\right)$$E89

where ${\alpha }_4\le 1$ and ${\varphi }_4\left({\widehat{Z}}_0\right)=\sqrt{\frac{q_{U1}}{2}{\widehat{z}}_{1,0}^2+\frac{q_{U2}}{2}{\widehat{z}}_{2,0}^2}$.

Now, we consider the composite Lyapunov function candidate of the full system given as follow

$$V_{\mathrm{SFU}}\left({\tilde{\chi }}_0,{\widehat{z}}_0\right)=\left(1-{\eta }_2\right)V_{\mathrm{SF}}\left({\tilde{\chi }}_0\right)+{\eta }_2{V}_{\mathrm{U}}\left({\widehat{z}}_0\right)$$E90

with $0<{\eta }_2<1$. The derivative of $V_{\mathrm{SFU}}\left({\tilde{\chi }}_0,{\widehat{z}}_0\right)$ presents new terms, which represents the effect of the interconnection between the slow and fast dynamics. These interconnections are assumed to satisfy the following conditions:

$${\left(\frac{\partial V_{\mathrm{SF}}\left({\tilde{\chi }}_0\right)}{\partial {\tilde{\chi }}_0}\right)}^T\rho \left({\tilde{F}}_0\left({\tilde{\chi }}_0,{\tilde{Z}}_0\right)-{\tilde{F}}_0\left({\tilde{\chi }}_0,{\tilde{h}}_0\left({\tilde{\chi }}_0\right)\right)\right)\le {\theta }_3{\varphi }_4\left({\widehat{Z}}_0\right){\varphi }_3\left({\tilde{\chi }}_0\right)$$E91

$${\left(\frac{\partial V_u\left({\widehat{Z}}_0\right)}{\partial {\tilde{\chi }}_0}\right)}^T\rho {\tilde{F}}_0\left({\tilde{\chi }}_0,{\tilde{h}}_0\left({\tilde{\chi }}_0\right)\right)\le {\theta }_4{\varphi }_4\left({\widehat{Z}}_0\right){\varphi }_3\left({\tilde{\chi }}_0\right)$$E92

$${\left(\frac{\partial V_u\left({\widehat{Z}}_0\right)}{\partial {\tilde{\chi }}_0}\right)}^T\rho \left[{\tilde{F}}_0\left({\tilde{\chi }}_0,{\tilde{Z}}_0\right)-{\tilde{F}}_0\left({\tilde{\chi }}_0,{\tilde{h}}_0\left({\tilde{\chi }}_0\right)\right)\right]\le {\gamma }_2{\varphi }_4^2\left({\widehat{Z}}_0\right)$$E93

the constants ${\theta }_3$, ${\theta }_4$, and ${\gamma }_2$ are non-negative. Therefore, from the singular perturbation theory (e.g. Theorem 11.3 in Refs. [19, 23, 24]), it follows that the derivative of V_SFU is negative-definite when

$${\epsilon }_2<{\epsilon }_2^{\otimes }\mathit{\text{for}}\left({\eta }_1<{\eta }_1^{\otimes }\right)$$E94

with ${\epsilon }_2^{\otimes }=\frac{{\alpha }_1{\alpha }_2}{{\alpha }_1{\gamma }_1+{\theta }_1{\theta }_2}$ and ${\eta }_1^{\otimes }=\frac{{\theta }_1}{{\theta }_1+{\theta }_2}$.

Therefore, it can be inferred that the equilibrium $X_0^{\ast }={\left(\begin{array}{cccc}{x}_{1,0}^{\ast }& x_{2,0}^{\ast }& x_{3,0}^{\ast }& x_{4,0}^{\ast }\end{array}\right)}^T$, $Y_0^{\ast }={\left(\begin{array}{cc}{y}_{1,0}^{\ast }& y_{2,0}^{\ast }\end{array}\right)}^T$, and $Z_0^{\ast }={\left(\begin{array}{cc}{z}_{1,0}^{\ast }& z_{2,0}^{\ast }\end{array}\right)}^T$ of full system, is asymptotically stable for all ${\epsilon }_i<{\epsilon }_i^{\otimes }\left(i=1,2\right)$. Therefore, Part 2 immediately follows from the averaging theory (e.g. Theorem 10.4 in Ref. [19]). Theorem is established.