**Transient-Flow Modeling of Vertical Fractured Wells with Multiple Hydraulic Fractures in Stress-Sensitive Gas Reservoirs**

#### **Ping Guo 1,\*, Zhen Sun 1, Chao Peng 2, Hongfei Chen <sup>3</sup> and Junjie Ren 1,4**


Received: 18 February 2019; Accepted: 28 March 2019; Published: 31 March 2019

**Abstract:** Massive hydraulic fracturing of vertical wells has been extensively employed in the development of low-permeability gas reservoirs. The existence of multiple hydraulic fractures along a vertical well makes the pressure profile around the vertical well complex. This paper studies the pressure dependence of permeability to develop a seepage model of vertical fractured wells with multiple hydraulic fractures. Both transformed pseudo-pressure and perturbation techniques have been employed to linearize the proposed model. The superposition principle and a hybrid analytical-numerical method were used to obtain the bottom-hole pseudo-pressure solution. Type curves for pseudo-pressure are presented and identified. The effects of the relevant parameters (such as dimensionless permeability modulus, fracture conductivity coefficient, hydraulic-fracture length, angle between the two adjacent hydraulic fractures, the difference of the hydraulic-fracture lengths, and hydraulic-fracture number) on the type curve and the error caused by neglecting the stress sensitivity are discussed in detail. The proposed work can enrich the understanding of the influence of the stress sensitivity on the performance of a vertical fractured well with multiple hydraulic fractures and can be used to more accurately interpret and forecast the transient pressure.

**Keywords:** gas flow; stress-sensitive porous media; multiple hydraulic fractures; vertical fractured well

#### **1. Introduction**

Gas flow in porous media has recently attracted much attention and stimulated great interest in the development of oil and gas reservoirs [1–15]. In the past decades, many gas reservoirs with stress-sensitive permeability have been discovered and developed around the world. During the development process of stress-sensitive reservoirs, the continuously decreasing formation pressure leads to the decrease of the formation permeability. Therefore, it is critical to investigate the effect of the stress sensitivity of the permeability on the production. Recently, lots of models have been established to study the performance of various wells in stress-sensitive reservoirs [16–18].

On the other hand, in order to obtain economic benefit, hydraulic fracturing has been widely applied in the development of low-permeability gas reservoirs. A great variety of seepage models of vertical fractured wells have been proposed and the characteristic of pressure response has been studied in detail. Gringarten et al. [19] established a seepage model of a vertical fractured well with two symmetrical infinite-conductivity hydraulic fractures, which can be used to identify the linear flow. Later, by considering the effect of the fluid flow within the hydraulic fractures, a seepage

model of a vertical fractured well with two symmetrical finite-conductivity hydraulic fractures was established [20,21], which can identify the bilinear flow. Some researchers pointed out the existence of two asymmetrical hydraulic fractures in practice and established some seepage models of vertical fractured wells with two asymmetrical hydraulic fractures [22–27]. In fact, massive hydraulic fracturing usually creates fractured network or several main hydraulic fractures near the wellbore instead of two symmetrical/asymmetrical hydraulic fractures [28–31]. Recently, seepage models of vertical fractured wells with multiple hydraulic fractures have received more attention. Shaoul et al. [32] proposed a numerical model of a vertical fractured well with multiple hydraulic fractures and investigated the pressure behavior of a vertical fractured well. Restrepo and Tiab [33] and Wang et al. [34] established semi-analytical models of a vertical fractured well with multiple infinite-conductivity hydraulic fractures. Ren and Guo [35] and Luo and Tang [36] proposed semi-analytical models of a vertical fractured well with multiple finite-conductivity hydraulic fractures. Compared with the pressure behavior of vertical fractured wells with two hydraulic fractures, the pressure response of vertical fractured wells with multiple hydraulic fractures is very complex. Although the effect of relevant parameters on the pressure response of vertical fractured wells with multiple hydraulic fractures has been studied, and the characteristics of the type of curves is not clearly understood. More important, little work has focused on vertical fractured wells with multiple hydraulic fractures in stress-sensitive gas reservoirs. In particular, the effect of relevant parameters on the error caused by neglecting the stress sensitivity has not been recognized clearly.

In this paper, we propose a seepage model of vertical fractured wells with multiple finiteconductivity hydraulic fractures. Type curves are identified and analyzed. The effects of relevant parameters on type curves and the error caused by neglecting the stress sensitivity are discussed in detail. This work can enrich the understanding of the influence of stress sensitivity on the performance of vertical fractured wells with multiple hydraulic fractures.

#### **2. Physical Model**

Complex geometry patterns of hydraulic fractures for vertical fractured wells can be described by the multiple-fracture model which is shown in Figure 1. Some assumptions of the present model are given as follows:

(1) The gas reservoir is a homogeneous and laterally infinite formation with a constant thickness. The top and bottom boundaries of the gas reservoir are assumed to be impermeable.

(2) The finite-conductivity hydraulic fractures emanate from the wellbore in the horizontal plane (as shown in Figure 1b) and completely penetrate the formation in the vertical plane. The length of each hydraulic fracture may be different.

(3) The vertical fractured well produces only from the hydraulic fractures and the production rate of the vertical fractured well keeps constant.

(4) The permeabilities of both the reservoir and hydraulic fractures change with pressure.

(5) The gas reservoir has a constant temperature and uniform initial pressure.

**Figure 1.** Vertical fractured well with multiple hydraulic fractures and its model. (**a**) Complex geometry pattern, (**b**) Multiple fracture model.

#### **3. Mathematical Model**

A semi-analytical model of vertical fractured well with multiple finite-conductivity hydraulic fractures in stress-sensitive gas reservoirs is established and solved by incorporating various methods and techniques. The nomenclatures of the present model are listed in Appendix A. The flow chart of solution methodology is shown in Figure 2.

**Figure 2.** Flow chart of the solution methodology.

#### *3.1. Flow Model for Stress-Sensitive Gas Reservoirs*

The reservoir permeability changes with reservoir pressure and can be described as follows [37]:

$$k = k\_i e^{-\gamma(\psi - \psi)},\tag{1}$$

where *ψ* is the pseudo-pressure, which is defined as

$$
\Psi = \int\_0^p \frac{2p}{\mu Z} \mathrm{d}p.\tag{2}
$$

The governing equation of gas flow in stress-sensitive porous media is given as follows [37]:

$$
\frac{\partial^2 \psi}{\partial r^2} + \frac{1}{r} \frac{\partial \psi}{\partial r} + \gamma \left(\frac{\partial \psi}{\partial r}\right)^2 = \epsilon^{\gamma(\psi\_i - \psi)} \frac{\phi \mu c\_t}{k\_i} \frac{\partial \psi}{\partial t}.\tag{3}
$$

Both gas viscosity *μ* and total compressibility *c*<sup>t</sup> in Equation (3) are functions of pressure. In order to simplify the proposed model and obtain an analytical solution, many research papers [26,37,38] treat the values of *μ* and *c*<sup>t</sup> in Equation (3) as constants, which could make the simulation results be accurate enough for engineering requirements. Following the suggestion by Wang [37], the values of *μ* and *c*<sup>t</sup> are evaluated at the initial condition, i.e., *μ* = *μ*(*p*i) and *c*<sup>t</sup> = *c*t(*p*i) in this paper.

With the assumption of the uniform initial pressure, the initial condition is expressed as

$$
\psi(r, t=0) = \psi\_{\text{i.}} \tag{4}
$$

The gas reservoir is assumed to be laterally infinite, and thus, the outer boundary condition is

$$
\psi(r \to \infty, t) = \psi\_{\text{i.}} \tag{5}
$$

The inner boundary condition for a point source can be written as

$$\lim\_{r \to 0} e^{-\gamma(\psi\_i - \psi)} r \frac{\partial \psi}{\partial r} = \frac{p\_{\rm sc} q(r = 0, t) T}{\pi k\_i h T\_{\rm sc}}.\tag{6}$$

Introducing the dimensionless variables listed in Table 1, Equations (3)–(6) are rewritten as

$$\frac{\partial^2 \psi\_{\rm D}}{\partial r\_{\rm D}^2} + \frac{1}{r\_{\rm D}} \frac{\partial \psi\_{\rm D}}{\partial r\_{\rm D}} - \gamma\_{\rm D} \left(\frac{\partial \psi\_{\rm D}}{\partial r\_{\rm D}}\right)^2 = \varepsilon^{\gamma\_{\rm D} \psi\_{\rm D}} \frac{\partial \psi\_{\rm D}}{\partial t\_{\rm D}},\tag{7}$$

$$
\psi\_{\rm D}(r\_{\rm D}, t\_{\rm D} = 0) = 0,\tag{8}
$$

$$
\psi\_{\rm D}(r\_{\rm D}, t\_{\rm D} = 0) = 0,\tag{9}
$$

$$\lim\_{r\_{\rm D}\to 0} e^{-\gamma\_{\rm D}\psi\_{\rm D}} r\_{\rm D} \frac{\partial \psi\_{\rm D}}{\partial r\_{\rm D}} = -q\_{\rm D}(r\_{\rm D} = 0, t\_{\rm D}). \tag{10}$$

**Table 1.** Definitions of dimensionless variables.


The point source model including Equations (7)–(10) shows strong non-linearity which can be alleviated by introducing the following expression [16]:

$$\psi\_{\rm D}(r\_{\rm D}, t\_{\rm D}) = -\frac{1}{\gamma\_{\rm D}} \ln[1 - \gamma\_{\rm D} \mathbf{j}\_{\rm D}(r\_{\rm D}, t\_{\rm D})].\tag{11}$$

The perturbation technique, Laplace transform, and superposition principle [39] were employed to obtain the pressure distribution in the gas reservoir with a vertical fractured well with multiple hydraulic fractures producing at a constant rate (see Appendix B):

$$\tilde{\zeta}\_{\rm D0}(r\_{\rm D},\theta,s) = \sum\_{j=1}^{n} \left[ \int\_{0}^{L\_{\rm fD}} \overline{q}\_{fD}(\mathbf{a},\theta\_{\rm fj},s) \mathbb{K}\_{\rm 0}(\sqrt{s}\sqrt{r\_{\rm D}^{2} + a^{2} - 2r\_{\rm D}\mathbf{a}\cos\left(\theta - \theta\_{\rm fj}\right)}) d\mathbf{a} \right].\tag{12}$$

#### *3.2. Flow Model for Stress-Sensitive Hydraulic Fractures*

The dimensionless model for gas flow within stress-sensitive hydraulic fractures is given as (see Appendix C)

$$\frac{\partial^2 \zeta\_{\rm fD}(\mathbf{x}\_{\rm iD}, t\_{\rm D})}{\partial \mathbf{x}\_{\rm iD}^2} - \frac{2\pi}{\mathbb{C}\_{\rm fD}} q\_{\rm fD}(\mathbf{x}\_{\rm iD}, t\_{\rm D}) = 0, \ (0 < \mathbf{x}\_{\rm iD} < L\_{\rm fD}), \tag{13}$$

$$\left.\frac{\partial\zeta\_{\rm fD}(x\_{\rm iD},t\_{\rm D})}{\partial x\_{\rm iD}}\right|\_{x\_{\rm iD}=0} = -\frac{2\pi}{\mathcal{C}\_{\rm fD}}Q\_{\rm iD}(t\_{\rm D}).\tag{14}$$

$$Q\_{\rm iD}(t\_{\rm D}) = \int\_{0}^{L\_{\rm iD}} q\_{fD}(\mathbf{x}\_{\rm iD}, t\_{\rm D}) d\mathbf{x}\_{\rm iD} \tag{15}$$

$$\sum\_{j=1}^{n} Q\_{j\mathbf{D}}(t\_{\mathbf{D}}) = 1.\tag{16}$$

According to the relationship between Cartesian coordinate system (*xi*, *yi*) and polar coordinate system (*r*, *θ*) (see Figure 3), Equations (13)–(16) can be rewritten in the polar coordinate system as follows:

$$\frac{\partial^2 \zeta\_{\rm fD}(r\_{\rm D}, \theta\_{\rm fi}, t\_{\rm D})}{\partial r\_{\rm D}^2} - \frac{2\pi}{\mathcal{C}\_{\rm fD}} q\_{\rm fD}(r\_{\rm D}, \theta\_{\rm fi}, t\_{\rm D}) = 0, \ (0 < r\_{\rm D} < L\_{\rm fD}), \tag{17}$$

$$\left.\frac{\partial\xi\_{\rm fD}(r\_{\rm D},\theta\_{\rm fi},t\_{\rm D})}{\partial r\_{\rm D}}\right|\_{r\_{\rm D}=0} = -\frac{2\pi}{\mathcal{C}\_{\rm fD}}Q\_{\rm iD}(t\_{\rm D}),\tag{18}$$

$$Q\_{\rm iD}(t\_{\rm D}) = \int\_{0}^{L\_{\rm iD}} q\_{f\rm D}(r\_{\rm D}, \theta\_{\rm fi}, t\_{\rm D}) \mathrm{d}r\_{\rm D} \tag{19}$$

$$\sum\_{j=1}^{n} Q\_{j\mathbf{D}}(t\_{\mathbf{D}}) = 1.\tag{20}$$

**Figure 3.** Schematic of different coordinate systems.

Integrating Equation (17) from 0 to *x* with respect to *r*<sup>D</sup> yields

$$\frac{\partial \mathfrak{J}\_{\rm fD}(\mathbf{x}', \theta\_{\rm fi}, t\_{\rm D})}{\partial r\_{\rm D}} - \frac{\partial \mathfrak{J}\_{\rm fD}(0, \theta\_{\rm fi}, t\_{\rm D})}{\partial r\_{\rm D}} = \frac{2\pi}{\mathcal{C}\_{\rm fD}} \int\_{0}^{\mathbf{x}'} q\_{fD}(\mathbf{x}'', \theta\_{\rm fi}, t\_{\rm D}) d\mathbf{x}''. \tag{21}$$

Substituting Equations (18) and (19) into Equation (21) yields

$$\frac{\partial \mathfrak{J}\_{\rm fD}^{\rm r}(\mathbf{x}', \theta\_{\rm fi}, t\_{\rm D})}{\partial r\_{\rm D}} + \frac{2\pi t}{\mathsf{C}\_{\rm fD}} \int\_{0}^{L\_{\rm fD}} q\_{f\rm D}(r\_{\rm D}, \theta\_{\rm fi}, t\_{\rm D}) \mathrm{d}r\_{\rm D} = \frac{2\pi}{\mathsf{C}\_{\rm fD}} \int\_{0}^{\mathbf{x}'} q\_{f\rm D}(\mathbf{x}'', \theta\_{\rm fi}, t\_{\rm D}) \mathrm{d}\mathbf{x}''. \tag{22}$$

Integrating Equation (22) from 0 to *r*<sup>D</sup> with respect to *r*<sup>D</sup> yields

$$\mathbf{f}\_{\rm SWID}^{\rm x}(\mathbf{t\_D}) - \mathbf{f}\_{\rm fD}^{\rm r}(r\_{\rm D}, \theta\_{\rm fi}, t\_{\rm D}) = \frac{2\pi}{\mathbf{C}\_{\rm fD}} \left[ r\_{\rm D} \int\_{0}^{\mathbf{L}\_{\rm fD}} q\_{f\rm D}(r\_{\rm D}, \theta\_{\rm fi}, t\_{\rm D}) \mathrm{d}r\_{\rm D} - \int\_{0}^{r\_{\rm D}} \int\_{0}^{\mathbf{x'}} q\_{f\rm D}(\mathbf{x''}, \theta\_{\rm fi}, t\_{\rm D}) \mathrm{d}x'' \mathrm{d}x' \right], \tag{23}$$

where

$$
\underline{\chi}\_{\rm wHD}(t\_{\rm D}) = \underline{\chi}\_{\rm dD}(0, \theta\_{\rm fi}, t\_{\rm D}).\tag{24}
$$

#### *3.3. Coupled Discrete Model*

The flow model for gas reservoirs can associate with the flow model for hydraulic fractures by the following relationship:

$$
\Psi\_{\rm fD}(r\mathrm{D}\_{\prime}\theta\_{\mathrm{fi}},t\mathrm{D}) = \Psi\mathrm{D}(r\mathrm{D}\_{\prime}\theta\_{\mathrm{fi}},t\mathrm{D}), \ (0 < r\mathrm{D} < L\_{\mathrm{fiD}}).\tag{25}
$$

Substituting Equations (11) and (A37) into Equation (25), one can derive

$$
\delta\_{\rm fD}^{\mathbb{Z}}(\mathbf{r\_D}, \theta\_{\rm fi}, t\_{\rm D}) = \delta\_{\rm D}^{\mathbb{Z}}(\mathbf{r\_D}, \theta\_{\rm fi}, t\_{\rm D}). \tag{26}
$$

The zero-order perturbation solution of the *ξ*<sup>D</sup> is accurate enough in practice, and thus we can obtain

$$
\delta\_\mathrm{D}(r\_\mathrm{D}, \theta\_\mathrm{fi}, t\_\mathrm{D}) = \delta\_\mathrm{D0}(r\_\mathrm{D}, \theta\_\mathrm{fi}, t\_\mathrm{D}). \tag{27}
$$

Combining Equations (26) and (27) yields that

$$
\check{\varsigma}\_{\rm fD}(r\_{\rm D}, \theta\_{\rm fi}, t\_{\rm D}) = \check{\varsigma}\_{\rm D0}(r\_{\rm D}, \theta\_{\rm fi}, t\_{\rm D}).\tag{28}
$$

Substituting Equation (28) into Equation (23) leads to

$$\mathbf{f}\_{\rm wHD}(t\_{\rm D}) - \mathbf{f}\_{\rm D0}(r\_{\rm D}, \theta\_{\rm fi}, t\_{\rm D}) = \frac{2\pi}{\mathbf{C}\_{\rm dD}} \left[ r\_{\rm D} \int\_{0}^{l\_{\rm dD}} q\_{fD}(r\_{\rm D}, \theta\_{\rm fi}, t\_{\rm D}) \mathrm{d}r\_{\rm D} - \int\_{0}^{r\_{\rm D}} \int\_{0}^{\mathbf{x'}} q\_{fD}(\mathbf{x''}, \theta\_{\rm fi}, t\_{\rm D}) \mathrm{d}\mathbf{x''} \mathrm{d}\mathbf{x'} \right]. \tag{29}$$

Taking the Laplace transform of Equation (29), one can obtain

$$\mathcal{J}\_{\text{wHD}}(\mathbf{s}) - \mathcal{J}\_{\text{DD}}(r\_{\text{D}}, \theta\_{\text{fi}}, \mathbf{s}) = \frac{2\pi}{\mathcal{C}\_{\text{HD}}} \Big[ r\_{\text{D}} \int\_{0}^{L\_{\text{HD}}} \overline{q}\_{fD}(r\_{\text{D}}, \theta\_{\text{fi}}, \mathbf{s}) \mathrm{d}r\_{\text{D}}} - \int\_{0}^{r\_{\text{D}}} \int\_{0}^{\mathbf{x'}} \overline{q}\_{fD}(\mathbf{x''}, \theta\_{\text{fi}}, \mathbf{s}) \mathrm{d}x'' \mathbf{d}x' \Big]. \tag{30}$$

Substituting Equation (12) into Equation (30), one can derive

$$\begin{split} \overline{\mathbb{J}}\_{\text{wHD}}(s) &- \sum\_{j=1}^{n} \left[ \int\_{0}^{\text{L}\_{\text{f}} \text{p}} \overline{q}\_{fD}(a, \theta\_{\text{f}}, s) \mathbb{K}\_{\text{0}}(\sqrt{s}\sqrt{r\_{\text{D}}^{2} + a^{2} - 2r\_{\text{D}}a\cos\left(\theta\_{\text{f}i} - \theta\_{\text{f}j}\right)} \right] da \right] \\ &= \frac{2\pi}{\text{L}\text{D}} \left[ r\_{\text{D}} \int\_{0}^{\text{L}\_{\text{f}} \text{p}} \overline{q}\_{fD}(r\_{\text{D}}, \theta\_{\text{f}i}, s) \text{d}r\_{\text{D}} - \int\_{0}^{\text{rD}} \int\_{0}^{\text{x'}} \overline{q}\_{fD}(\mathbf{x''}, \theta\_{\text{f}i}, s) \text{d}\mathbf{x''} \mathbf{d}\mathbf{x'}} \right]. \end{split} \tag{31}$$

Employing the Laplace transform of Equations (19) and (20), one can obtain

$$\sum\_{i=1}^{n} \int\_{0}^{L\_{\rm fD}} \overline{q}\_{fD}(r\_{\rm D}, \theta\_{\rm fi}, s) \mathrm{d}r\_{\rm D} = \frac{1}{s}. \tag{32}$$

The coupled model consists of Equations (31) and (32), which can be solved by the numerical discrete method [40]. Each hydraulic fracture is discretized into some segments, and the flow rate in each segment is considered to stay the same at a certain time. The discrete schematic of a hydraulic fracture of a vertical fractured well is shown in Figure 4.

**Figure 4.** Discrete schematic of the *i*th hydraulic fracture of a vertical fractured well.

And then, the discrete forms of Equations (31) and (32) can be obtained as follows

*i*=1

*j*=1

$$
\overline{\mathbf{T}}\_{\text{wHD}}(s) - \sum\_{i=1}^{n} \sum\_{j=1}^{N\_i} \left[ \overline{q}\_{\text{f},j\text{D}}(s) \int\_{\overline{r}\_{i,\text{D}}}^{r\_{i,j}} \mathbf{K}\_0(\sqrt{s} \sqrt{r\_{\text{m}k,\text{rD}}^2 + a^2 - 2r\_{\text{m}k,\text{rD}}a \cos(\theta\_{\text{f}\text{k}} - \theta\_{\text{f}\text{i}})} \right] da \right] \tag{33}
$$

$$
\frac{2\pi}{\varepsilon\_{\text{dD}}} \left\{ \sum\_{i=1}^{n-1} \overline{q}\_{\text{f}\text{k},\text{jD}}(s) \cdot \left[ (v - i) \cdot \Delta L\_{\text{f\text{E}D}}^2 \right] + \frac{\Delta L\_{\text{f\text{E}D}}^2}{8} \cdot \overline{q}\_{\text{f}\text{k},\text{rD}}(s) - r\_{\text{m}k,\text{rD}} \sum\_{i=1}^{n} \overline{q}\_{\text{f}\text{k},\text{jD}}(s) \Delta L\_{\text{f\text{E}D}} \right\} = 0,
$$

$$
\sum\_{i=1}^{n} \sum\_{j=1}^{N\_i} \left[ \overline{q}\_{\text{f}\text{i},j\text{D}}(s) \Delta L\_{\text{f\text{E}D}} \right] = \frac{1}{s},\tag{34}
$$

where 1 ≤ *k* ≤ *n*, 1 ≤ *v* ≤ *Nk* and Δ*L*f*i*<sup>D</sup> = *L*f*i*D/*Ni*.

Equations (33) and (34) compose *<sup>n</sup>* ∑ *i*=1 *Ni* + 1 linear equations with *<sup>n</sup>* ∑ *i*=1 *Ni* + 1 unknowns. The unknowns are *ξ*wHD(*s*) and *q*f*k*,*v*D(*s*) (1 ≤ *k* ≤ *n*, 1 ≤ *v* ≤ *Nk*), which can be obtained by solving the linear equations, and then the wellbore storage and the skin near the wellbore can be taken into account by the following formula [41]

$$\mathfrak{J}\_{\rm wD}(s) = \frac{s\mathfrak{J}\_{\rm wHD}(s) + S\_{\rm f}}{s + s^2 \mathbb{C}\_{\rm D} \Big[s\mathfrak{J}\_{\rm wHD}(s) + S\_{\rm f}\big]}. \tag{35}$$

With the aid of the numerical inversion method [42], *ξ*wD(*s*) in Laplace space is transformed into *ξ*wD(*t*D) in real space

$$\zeta\_{\rm wD}^{\rm T}(t\_{\rm D}) = \frac{\ln 2}{t\_{\rm D}} \sum\_{i=1}^{N} V\_i\_{i\rm wD}^{\rm T}(s), \ (N = 8), \tag{36}$$

where

$$s = \frac{i \ln 2}{t\_{\rm D}},\tag{37}$$

$$V\_i = (-1)^{\frac{N}{2} + i} \sum\_{k=\frac{i+1}{2}}^{\min\left(i, \frac{N}{2}\right)} \frac{k^{\frac{N}{2} + 1} (2k)!}{\left(\frac{N}{2} - k\right)! k! (k - 1)! (i - k)! (2k - i)!}. \tag{38}$$

Finally, the dimensionless bottom-hole pseudo-pressure of a vertical fractured well with constant production rate is obtained by *ψ*wD(*t*D) = −ln[1 − *γ*<sup>D</sup> · *ξ*wD(*t*D)]/*γ*D.

#### **4. Model Validation and Application**

Because there is no analytical inversion solution of *ξ*wD(*s*), it is difficult to discuss the accuracy of the Laplace transform by comparing the analytical inversion solution with the numerical inversion solution [43]. In order to validate the proposed model and show how this model was used in practice, the drawdown test data of a vertical fractured well with four fracture wings were collected from the published literature [44]. The drawdown test data in the literature [44] were generated by the reservoir simulator which is based on the implicit finite difference method. The vertical fractured well produces at a constant production rate of 0.1639 m3/s. Basic parameters of the vertical fractured well are listed in Table 2. The proposed model is employed to simulate the bottom-hole pressure under the constant-rate-production condition. The simulated bottom-hole pseudo-pressure and bottom-hole pressure are compared with the results published in the literature [44]. It is shown from Figure 5 that the simulated bottom-hole pseudo-pressure and bottom-hole pressure by the proposed model are in good agreement with the bottom-hole pressure data published in the literature [44], which validates the proposed model.

Furthermore, the proposed model can be used to predict the bottom-hole pressure under the constant-rate-production condition. Based on the proposed model, it was easy to obtain the bottom-hole pressure of the vertical fractured well at any given time. For example, when the vertical fractured well produces at a constant production rate of 0.1639 m3/s, the bottom-hole pressures are 41.132 MPa, 37.676 MPa, 31.484 MPa at *t* = 100, 1000, 10000 hours, respectively.


**Table 2.** Parameters of a vertical fractured well with four fracture wings in the published literature [44].

**Figure 5.** Comparisons of bottom-hole pseudo-pressure and bottom-hole pressure between our simulation results and the test data [44]. (**a**) bottom-hole pseudo -pressure, (**b**) bottom-hole pressure

#### **5. Results and Analysis**

#### *5.1. Type Curves and Flow Regimes*

The bottom-hole pseudo-pressure of a vertical fractured well with multiple hydraulic fractures in stress-sensitive gas reservoirs was calculated by the proposed model. Type curves for the bottom-hole pseudo-pressure were plotted and the characteristics of the bottom-hole pseudo-pressure behavior was analyzed. The established model in this work was suitable for multiple hydraulic fractures with arbitrary fracture number, arbitrary fracture length, and arbitrary fracture orientation. In order to investigate the effect of the hydraulic-fracture distribution on the pseudo-pressure response, without loss of generality, we mainly focused on the vertical fractured well with four fracture wings which is shown in Figure 6. The first and second fracture wings were assumed to have the equal length (*X*f), and the third and fourth fracture wings were set to be the equal length (*Y*f).

Figure 7 shows the type curves for the bottom-hole pseudo-pressure of a vertical fractured well with multiple hydraulic fractures in stress-sensitive gas reservoirs. The type curves for the bottom-hole pseudo-pressure consist of the pseudo-pressure curve (*x*−axis: *t*D/*C*D, *y*−axis: *ψ*wD) and the pseudo-pressure derivative curve (*x*−axis: *t*D/*C*D, *y*−axis: *ψ* wD · *t*D/*C*D). It is observed that there are seven possible flow regions in the type curves as follows:

**Figure 6.** Schematic of a vertical fractured well with four fracture wings.

**Figure 7.** Type curves for the bottom-hole pseudo-pressure of a vertical fractured well with multiple hydraulic fractures in stress-sensitive gas reservoirs with different values of the *γ*<sup>D</sup> (*n* = 4, *L*f1 = *L*f2 = *<sup>X</sup>*<sup>f</sup> <sup>=</sup> 20 m, *<sup>L</sup>*f3 <sup>=</sup> *<sup>L</sup>*f4 <sup>=</sup> *<sup>Y</sup>*<sup>f</sup> <sup>=</sup> 80 m, *<sup>θ</sup>* <sup>=</sup> <sup>30</sup>◦, *<sup>C</sup>*fD <sup>=</sup> <sup>5</sup> <sup>×</sup> 104, *<sup>S</sup>*<sup>f</sup> <sup>=</sup> <sup>10</sup><sup>−</sup>3, *<sup>C</sup>*<sup>D</sup> <sup>=</sup> 0.1).

(1) Wellbore storage period (WSP): The produced gas comes from the storage gas in the wellbore, and the gas in the reservoir has not flowed into the wellbore. Both the pseudo-pressure curve and the pseudo-pressure derivative curve exhibit as straight lines with the slope being one.

(2) Transitional flow after wellbore storage period (TFAWSP): the gas in the reservoir begins to flow into the wellbore. This flow period is controlled by both the wellbore storage and skin near the wellbore. The pseudo-pressure derivative curve appears as a "hump".

(3) Bilinear flow period (BFP): Linear flows take place within the hydraulic fractures and perpendicular to the hydraulic-facture surfaces in the reservoir simultaneously. The pseudo-pressure derivative curve exhibits as a straight line with the slope being 1/4.

(4) Transitional flow after bilinear flow period (TFABFP): The duration of this flow period is dependent on the difference of the hydraulic-fracture lengths. The shorter the hydraulic-fracture length is, the earlier the BFP ends. Therefore, this flow period will occur when the BFP of the shorter hydraulic fracture ends and the longer hydraulic fracture still lies in the BFP. The pseudo-pressure derivative curve appears as a straight line with the slope in the range of 1/4–1/2.

(5) Linear flow period (LFP): When the BFP of each hydraulic fracture ends, the LFP will appear. In this period, the linear flow perpendicular to the hydraulic-facture surfaces in the reservoir dominates in the area near the wellbore. The pseudo-pressure derivative curve appears as a straight line with the slope being 1/2.

(6) Transitional flow after linear flow period (TFALFP): When the pressure wave continues to travel in the reservoir, the interference between the pressure waves from different hydraulic fractures will take place, which makes the slope of the pseudo-pressure derivative curve be greater than 1/2. It should be noted that the stress sensitivity begins to obviously affect the type curves in this period. As the magnitude of the *γ*<sup>D</sup> increases, the pseudo-pressure and its derivative curves are shifted up.

(7) Pseudo-radial flow period (PRFP): Compared with the previous flow period (i.e., the TFALFP), much slower growth of the pseudo-pressure drop was observed. The pseudo-pressure derivative without the effect of the stress sensitivity (i.e., *γ*<sup>D</sup> = 0) keeps a value of 0.5, while the magnitude of the pseudo-pressure derivative with the effect of the stress sensitivity (i.e., *γ*<sup>D</sup> > 0) increases with the time. Furthermore, the pseudo-pressure and its derivative increase with increasing the value of the *γ*<sup>D</sup> at a fixed time. Therefore, if the stress sensitivity of the gas reservoir is stronger, the larger pressure drop is needed to remain the constant rate of produced well.

In order to quantify the effect of the stress sensitivity on the pseudo-pressure, the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity is introduced as follows

$$\delta = \frac{|\psi\_{\rm wDl} - \psi\_{\rm wDnl}|}{\Psi\_{\rm wDl}} \times 100\% \,\tag{39}$$

where *ψ*wDnl and *ψ*wDl are the bottom-hole pseudo-pressures with and without the effect of the stress sensitivity, respectively.

Figure 8 shows the effect of the dimensionless permeability modulus (*γ*D) on the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity. With the increase of the time, the relative difference *δ* first increases very slowly and finally increases rapidly. Considering the flow regimes, it was found that the impact of the stress sensitivity on the pseudo-pressure was negligibly small during and before the LFP, while after the LFP, this impact became more and more obvious with the increase of the time. Furthermore, the *δ* became larger with increasing the magnitude of the *γ*<sup>D</sup> at a fixed time.

**Figure 8.** Effect of the *γ*<sup>D</sup> on the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity (*n* = 4, *L*f1 = *L*f2 = *X*<sup>f</sup> = 20 m, *L*f3 = *L*f4 = *Y*<sup>f</sup> = 80 m, *θ* = 30◦, *<sup>C</sup>*fD <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>4, *<sup>S</sup>*<sup>f</sup> <sup>=</sup> <sup>10</sup><sup>−</sup>3, *<sup>C</sup>*<sup>D</sup> <sup>=</sup> 0.1).

#### *5.2. Sensitivity Analysis*

Besides the dimensionless permeability modulus (*γ*D), other parameters may have effects on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity. Therefore, it is critical to conduct the sensitivity analysis of relevant parameters.

Figures 9 and 10 show the effects of the fracture conductivity coefficient (*C*fD) on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity, respectively. It is obvious that the impact of the *C*fD on the type curve takes place in the early periods (i.e., from the TFAWSP to the LFP), where the values of the pseudo-pressure and its derivative increase with the decrease of the magnitude of the *C*fD. Furthermore, it was found from Figure 9 that as the value of the *C*fD increases, the duration of the BFP becomes shorter and the start time of the LFP becomes earlier. As shown in Figure 10, the *δ* was little affected by the *C*fD in all periods, and this was because the main impact of the *C*fD on the pseudo-pressure occurred in the early periods, but the effect of the stress sensitivity on the pseudo-pressure was negligibly small in the early periods.

**Figure 9.** Type curves for the bottom-hole pseudo-pressure of a vertical fractured well with multiple hydraulic fractures in stress-sensitive gas reservoirs with different values of the *C*fD (*n* = 4, *L*f1 = *L*f2 = *X*<sup>f</sup> = 50 m, *L*f3 = *L*f4 = *Y*<sup>f</sup> = 50 m, *θ* = 90◦, *S*<sup>f</sup> = 10<sup>−</sup>3, *C*<sup>D</sup> = 0.1, *γ*<sup>D</sup> = 0.1).

**Figure 10.** Effect of the *C*fD on the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity (*n* = 4, *L*f1 = *L*f2 = *X*<sup>f</sup> = 50 m, *L*f3 = *L*f4 = *Y*<sup>f</sup> = 50 m, *θ* = 90◦, *S*<sup>f</sup> = 10<sup>−</sup>3, *C*<sup>D</sup> = 0.1, *γ*<sup>D</sup> = 0.1).

Figures 11 and 12 show the effects of the length of the hydraulic fractures (*L*f) on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity, respectively. It is seen from Figure 11 that the *L*<sup>f</sup> has an important effect on the magnitudes of the pseudo-pressure and its derivative in the intermediate and late flow periods (i.e., from the BFP to the PRFP). Decreasing the *L*<sup>f</sup> reduced the duration of the BFP and resulted in the increase of the pseudo-pressure and its derivative during and after the BFP. Furthermore, it was found from Figure 12 that the *δ* increased with decreasing the magnitude of the *L*<sup>f</sup> in the late flow period, indicating that the pseudo-pressure obtained by the conventional model without the effect of the stress sensitivity will result in a much bigger error when the length of the hydraulic fractures becomes shorter.

**Figure 11.** Type curves for the bottom-hole pseudo-pressure of a vertical fractured well with multiple hydraulic fractures in stress-sensitive gas reservoirs with different values of the *L*<sup>f</sup> (*n* = 4, *L*<sup>f</sup> = *L*f1 = *<sup>L</sup>*f2 <sup>=</sup> *<sup>L</sup>*f3 <sup>=</sup> *<sup>L</sup>*f4 <sup>=</sup> *<sup>X</sup>*<sup>f</sup> <sup>=</sup> *<sup>Y</sup>*f, *<sup>θ</sup>* <sup>=</sup> <sup>90</sup>◦, *<sup>C</sup>*fD <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>4, *<sup>S</sup>*<sup>f</sup> <sup>=</sup> <sup>10</sup><sup>−</sup>3, *<sup>C</sup>*<sup>D</sup> <sup>=</sup> 0.1, *<sup>γ</sup>*<sup>D</sup> <sup>=</sup> 0.1).

**Figure 12.** Effect of the *L*<sup>f</sup> on the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity (*<sup>n</sup>* <sup>=</sup> 4, *<sup>L</sup>*<sup>f</sup> <sup>=</sup> *<sup>L</sup>*f1 <sup>=</sup> *<sup>L</sup>*f2 <sup>=</sup> *<sup>L</sup>*f3 <sup>=</sup> *<sup>L</sup>*f4 <sup>=</sup> *<sup>X</sup>*<sup>f</sup> <sup>=</sup> *<sup>Y</sup>*f, *<sup>θ</sup>* <sup>=</sup> <sup>90</sup>◦, *<sup>C</sup>*fD <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>4, *S*<sup>f</sup> = 10<sup>−</sup>3, *C*<sup>D</sup> = 0.1, *γ*<sup>D</sup> = 0.1).

Figures 13 and 14 show the effects of the angle between the two adjacent hydraulic fractures (*θ*) on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity, respectively. As shown in Figure 13, with the decrease of the *θ*, the start time of the TFALFP becomes earlier, and the magnitudes of the pseudo-pressure and its derivative in the TFALFP increase. This is because decreasing the *θ* results in the earlier and stronger interference between the hydraulic fractures, which makes the pressure drop increase to remain the constant rate of the produced well. Whereas the impact of the *θ* on the type curve is relatively unobvious in the PRFP. Furthermore, it is seen from Figure 14 that the *δ* increases with decreasing the magnitude of the *θ* during and after the TFALFP, indicating that the effect of the stress sensitivity on the pseudo-pressure becomes greater when the *θ* decreases.

**Figure 13.** Type curves for the bottom-hole pseudo-pressure of a vertical fractured well with multiple hydraulic fractures in stress-sensitive gas reservoirs with different values of the *θ* (*n* = 4, *L*f1 = *L*f2 = *<sup>X</sup>*<sup>f</sup> <sup>=</sup> 50 m, *<sup>L</sup>*f3 <sup>=</sup> *<sup>L</sup>*f4 <sup>=</sup> *<sup>Y</sup>*<sup>f</sup> <sup>=</sup> 50 m, *<sup>C</sup>*fD <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>4, *<sup>S</sup>*<sup>f</sup> <sup>=</sup> <sup>10</sup><sup>−</sup>3, *<sup>C</sup>*<sup>D</sup> <sup>=</sup> 0.1, *<sup>γ</sup>*<sup>D</sup> <sup>=</sup> 0.1).

**Figure 14.** Effect of the *θ* on the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity (*<sup>n</sup>* <sup>=</sup> 4, *<sup>L</sup>*f1 <sup>=</sup> *<sup>L</sup>*f2 <sup>=</sup> *<sup>X</sup>*<sup>f</sup> <sup>=</sup> <sup>50</sup> m, *<sup>L</sup>*f3 <sup>=</sup> *<sup>L</sup>*f4 <sup>=</sup> *<sup>Y</sup>*<sup>f</sup> <sup>=</sup> <sup>50</sup> m, *<sup>C</sup>*fD <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>4, *S*<sup>f</sup> = 10<sup>−</sup>3, *C*<sup>D</sup> = 0.1, *γ*<sup>D</sup> = 0.1).

Introducing the coefficient of *β* = |*X*<sup>f</sup> − *Y*f|/(*X*<sup>f</sup> + *Y*f) to quantify the difference of the hydraulic-fracture lengths, the effects of *β* on the type curve, and the relative differences between the pseudo-pressures with and without the effect of the stress sensitivity are shown in Figures 15 and 16, respectively. It is obvious that the difference of the hydraulic-fracture lengths becomes larger as the magnitude of the *β* increases. The range of the *β* is from zero to one. As shown in Figure 15, with the decrease of the *β*, the end of the BFP takes place later, the duration of the TFABFP becomes shorter, and the start time of the LFP occurs earlier. The reason for the existence of the TFABFP is that the difference of the hydraulic-fracture lengths leads to the inconsistency of the end times for the BFP for each hydraulic fracture. If *β* is equal to zero, the TFABFP may not even appear. Furthermore, the difference of the hydraulic-fracture lengths makes the pseudo-pressure and its derivative increase

in the TFABFP. It is interesting to find in Figure 15 that during the TFALFP, the pseudo-pressure and its derivative increase as the *β* decreases. This may be because that for the distribution of the hydraulic fractures shown in Figure 6, the interference between hydraulic fractures enhances with the decrease of the *β* in the TFALFP. As shown in Figure 16, the *β* has an effect on the error caused by neglecting the stress sensitivity during and after the TFALFP. Decreasing the *β* results in the increase of the error.

**Figure 15.** Type curves for the bottom-hole pseudo-pressure of a vertical fractured well with multiple hydraulic fractures in stress-sensitive gas reservoirs with different values of the *β* (*n* = 4, *L*f1 = *L*f2 = *X*f, *<sup>L</sup>*f3 <sup>=</sup> *<sup>L</sup>*f4 <sup>=</sup> *<sup>Y</sup>*f, *<sup>X</sup>*<sup>f</sup> <sup>+</sup> *<sup>Y</sup>*<sup>f</sup> <sup>=</sup> 100 m, *<sup>θ</sup>* <sup>=</sup> <sup>90</sup>◦, *<sup>C</sup>*fD <sup>=</sup> <sup>5</sup> <sup>×</sup> 104, *<sup>S</sup>*<sup>f</sup> <sup>=</sup> <sup>10</sup><sup>−</sup>3, *<sup>C</sup>*<sup>D</sup> <sup>=</sup> 0.1, *<sup>γ</sup>*<sup>D</sup> <sup>=</sup> 0.1).

**Figure 16.** Effect of the *β* on the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity (*n* = 4, *L*f1 = *L*f2 = *X*f, *L*f3 = *L*f4 = *Y*f, *X*<sup>f</sup> + *Y*<sup>f</sup> = 100 m, *θ* = 90◦, *<sup>C</sup>*fD <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>4, *<sup>S</sup>*<sup>f</sup> <sup>=</sup> <sup>10</sup><sup>−</sup>3, *<sup>C</sup>*<sup>D</sup> <sup>=</sup> 0.1, *<sup>γ</sup>*<sup>D</sup> <sup>=</sup> 0.1).

Figures 17 and 18 show the effects of the hydraulic-fracture number (*n*) on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity, respectively. The distribution of the hydraulic fractures was assumed to be uniform (i.e., the angle between the two arbitrary adjacent hydraulic fractures was equal). As shown in Figure 17, it is obvious that the *n* affects the pseudo-pressure behavior in all periods except for the WSP. Increasing the *n* leads to the decrease of the pseudo-pressure and its derivative. Due to the increase of the *n*, the angle between the two adjacent hydraulic fractures was reduced, and thus the interference between hydraulic fractures occurred earlier and the effect of this interference on the pseudo-pressure was

more significant. Figure 18 indicates that a smaller value of the *n* leads to a bigger error caused by neglecting the stress sensitivity during and after the TFALFP.

**Figure 17.** Type curves for the bottom-hole pseudo-pressure of a vertical fractured well with multiple hydraulic fractures in stress-sensitive gas reservoirs with different values of the *n* (*L*f1 = *L*f2 = *X*<sup>f</sup> = 80 m, *<sup>L</sup>*f3 <sup>=</sup> *<sup>L</sup>*f4 <sup>=</sup> *<sup>Y</sup>*<sup>f</sup> <sup>=</sup> 80 m, *<sup>θ</sup>* <sup>=</sup> <sup>360</sup>◦/*n*, *<sup>C</sup>*fD <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>4, *<sup>S</sup>*<sup>f</sup> <sup>=</sup> <sup>10</sup><sup>−</sup>3, *<sup>C</sup>*<sup>D</sup> <sup>=</sup> 0.1, *<sup>γ</sup>*<sup>D</sup> <sup>=</sup> 0.1).

**Figure 18.** Effect of *n* on the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity (*L*f1 = *L*f2 = *X*<sup>f</sup> = 80 m, *L*f3 = *L*f4 = *Y*<sup>f</sup> = 80 m, *θ* = 360◦/*n*, *<sup>C</sup>*fD <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>4, *<sup>S</sup>*<sup>f</sup> <sup>=</sup> <sup>10</sup><sup>−</sup>3, *<sup>C</sup>*<sup>D</sup> <sup>=</sup> 0.1, *<sup>γ</sup>*<sup>D</sup> <sup>=</sup> 0.1).

#### **6. Conclusions**

We have established a non-linear model of a vertical fractured well with multiple hydraulic fractures in gas reservoirs by incorporating the pressure-dependent permeability. Both transformed pseudo-pressure and perturbation techniques were employed to linearize the proposed model. Superposition principle and numerical discrete methods were used to obtain the semi-analytical solution. Type curves for pseudo-pressure were plotted and discussed in detail. The influence of relevant parameters on the type curve and the relative difference between the pseudo-pressures with and without the effect of the stress sensitivity were analyzed. Some main conclusions are listed as follows:

(1) The type curve of a vertical fractured well with multiple hydraulic fractures was identified by seven possible flow regions. Compared with the conventional model of a vertical fractured well with two symmetrical hydraulic-fracture wings, two transitional flow regimes (i.e., the TFABFP and TFALFP) were observed, which were caused by the difference of the hydraulic-fracture lengths and the interference between hydraulic fractures, respectively.

(2) As the time increased, the pressure drops increased, which made the permeability decrease. The impact of the stress sensitivity on the pseudo-pressure increased with the increase of the time. Furthermore, the influence of the stress sensitivity on the pseudo-pressure was negligibly small in the early and intermediate flow periods but becomes very significant in the late flow period.

(3) Some relevant parameters, such as dimensionless permeability modulus, fracture conductivity coefficient, hydraulic-fracture length, angle between the two adjacent hydraulic fractures, the difference of the hydraulic-fracture lengths, and hydraulic-fracture number, not only affected the type curve, but also have an influence of the error caused by neglecting the stress sensitivity.

**Author Contributions:** Data curation, Z.S.; Formal analysis, P.G.; Funding acquisition, J.R.; Investigation, Z.S.; Methodology, P.G.; Project administration, P.G.; Visualization, H.C.; Writing—original draft, Z.S.; Writing—review and editing, C.P.

**Funding:** We acknowledge the support provided by the Open Fund of the State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation (PLN201627) and the Scientific Research Starting Project of SWPU (No. 2017QHZ031).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Symbol Description**

*Nomenclature*



#### **Appendix B. Pressure Distribution in a Gas Reservoir with a Vertical Fractured Well with Multiple Hydraulic Fractures Producing at a Constant Rate**

With the aid of Equation (11), Equations (7)–(10) can be rewritten as

$$\frac{\partial^2 \zeta\_\mathrm{D}}{\partial r\_\mathrm{D}^2} + \frac{1}{r\_\mathrm{D}} \frac{\partial \zeta\_\mathrm{D}}{\partial r\_\mathrm{D}} = \frac{1}{1 - \gamma\_\mathrm{D} \zeta\_\mathrm{D}} \frac{\partial \zeta\_\mathrm{D}}{\partial t\_\mathrm{D}},\tag{A1}$$

$$
\xi\_{\mathcal{D}}(r\_{\mathcal{D}}, t\_{\mathcal{D}} = 0) = 0,\tag{A2}
$$

$$\mathbf{q}\_{\rm D}(\mathbf{r}\_{\rm D}, \mathbf{t}\_{\rm D} = 0) = 0,\tag{A3}$$

$$\lim\_{r\_{\rm D} \to 0} r\_{\rm D} \frac{\partial \xi\_{\rm D}}{\partial r\_{\rm D}} = -q\_{\rm D}(r\_{\rm D} = 0, t\_{\rm D}). \tag{A4}$$

The perturbation technique of *γ*<sup>D</sup> was employed to obtain that

$$\mathbf{\upvarphi}\_{\rm D} = \mathbf{\upvarphi}\_{\rm D0} + \gamma \mathbf{\upvarphi}\_{\rm D1}^{\mathbf{x}} + \gamma \mathbf{\upvarphi}\_{\rm D}^{2} \mathbf{\upvarphi}\_{\rm D2} + \cdots,\tag{A5}$$

$$-\frac{1}{\gamma\_{\rm D}}\ln[1-\gamma\_{\rm D}\mathbb{1}\_{\rm D}] = \mathbb{s}\_{\rm D} + \frac{1}{2}\gamma\_{\rm D}\mathbb{s}\_{\rm D}^{2} + \cdots \ , \tag{A6}$$

$$\frac{1}{1 - \gamma\_{\rm D}\xi\_{\rm D}} = 1 + \gamma\_{\rm D}\xi\_{\rm D} + \gamma\_{\rm D}^2\xi\_{\rm D}^2 + \cdots \ . \tag{A7}$$

Considering a small value of *γ*<sup>D</sup> in general, the zero-order perturbation solution was accurate enough for engineering requirements. So, Equations (A1)–(A4) can be rewritten as

$$\frac{\partial^2 \zeta\_{\rm D0}}{\partial r\_{\rm D}^2} + \frac{1}{r\_{\rm D}} \frac{\partial \zeta\_{\rm D0}}{\partial r\_{\rm D}} = \frac{\partial \zeta\_{\rm D0}}{\partial t\_{\rm D}},\tag{A8}$$

$$\, \, \, \, \, \, \, \_{\rm D0}^{x} (r\_{\rm D}, t\_{\rm D} = 0) = 0,\tag{A9}$$

$$\,^\delta \mathfrak{g}\_{\mathcal{D}0}(r\_{\mathcal{D}} \to \infty, t\_{\mathcal{D}}) = 0,\tag{A10}$$

$$\lim\_{r\_{\rm D} \to 0} r\_{\rm D} \frac{\partial \zeta\_{\rm D0}}{\partial r\_{\rm D}} = -q\_{\rm D}(r\_{\rm D} = 0, t\_{\rm D}). \tag{A11}$$

Laplace transform was employed to solve the above model and it is defined by

$$\overline{f}(s) = \int\_0^\infty f(t\_\mathcal{D}) e^{-st\_\mathcal{D}} \mathbf{d}t\_{\mathcal{D}\_\prime} \tag{A12}$$

where *f* is an arbitrary variable in real space.

Taking the Laplace transform of Equations (A8)–(A11), one can obtain

$$\frac{\partial^2 \zeta\_{\rm D0}^{\rm T}}{\partial r\_{\rm D}^2} + \frac{1}{r\_{\rm D}} \frac{\partial \zeta\_{\rm D0}^{\rm T}}{\partial r\_{\rm D}} = s\_{\rm D0}^{\rm T} \tag{A13}$$

$$\,^\sharp \zeta\_{\rm D0}(r\_{\rm D} \to \infty, s) = 0,\tag{A14}$$

$$\lim\_{r\_{\rm D} \to 0} r\_{\rm D} \frac{\partial \xi\_{\rm D0}}{\partial r\_{\rm D}} = -\overline{q}\_{\rm D}(r\_{\rm D} = 0, s). \tag{A15}$$

The solution of Equations (A13)–(A15) in Laplace space can be derived as

$$
\overline{\zeta}\_{\rm D0}(r\_{\rm D}, s) = \overline{q}\_{\rm D}(r\_{\rm D} = 0, s) \mathbb{K}\_0(r\_{\rm D}\sqrt{s}).\tag{A16}
$$

According to Equation (A16), the superposition principle [39] was employed to obtain the pressure distribution in a gas reservoir with a vertical fractured well with multiple hydraulic fractures producing at a constant rate

$$\overline{\zeta}\_{\rm D0}(r\_{\rm D}, \theta, s) = \sum\_{j=1}^{n} \left[ \int\_{0}^{L\_{\rm fD}} \overline{q}\_{f\rm D}(\mathfrak{a}, \theta\_{\rm f}, s) \mathbb{K}\_{\mathbb{O}}(\sqrt{s}\sqrt{r\_{\rm D}^{2} + a^{2} - 2r\_{\rm D}a\cos\left(\theta - \theta\_{\rm f}\right)}) da \right]. \tag{A17}$$

#### **Appendix C. Derivation of Dimensionless Model for Gas Flow within Stress-Sensitive Hydraulic Fractures**

Considering the effect of the pressure-dependent permeability of hydraulic fractures, the governing equation for gas flow in the *i*th hydraulic fracture is

$$
\frac{\partial^2 \psi\_{\rm f}}{\partial \mathbf{x}\_i^2} + \frac{\partial^2 \psi\_{\rm f}}{\partial y\_i^2} + \gamma \left[ \left(\frac{\partial \psi\_{\rm f}}{\partial \mathbf{x}\_i}\right)^2 + \left(\frac{\partial \psi\_{\rm f}}{\partial y\_i}\right)^2 \right] = \epsilon^{\gamma(\varphi\_{\rm f} - \varphi\_{\rm f})} \frac{\not p\_{\rm f} \mu \epsilon\_{\rm ft}}{k\_{\rm ft}} \frac{\not \partial \psi\_{\rm f}}{\partial t} \,, \tag{A18}
$$

where 0 <sup>&</sup>lt; *xi* <sup>&</sup>lt; *<sup>L</sup>*f*i*, <sup>−</sup>*w*<sup>f</sup> <sup>2</sup> <sup>&</sup>lt; *yi* <sup>&</sup>lt; *<sup>w</sup>*<sup>f</sup> 2 .

The volume within the hydraulic fracture was very small, and thus the compressibility of the hydraulic fracture was neglected. Then Equation (A18) was reduced to

$$
\frac{
\partial^2 \psi\_\mathbf{f}
}{
\partial x\_i^2
} + \frac{
\partial^2 \psi\_\mathbf{f}
}{
\partial y\_i^2
} + \gamma \left[ \left(\frac{
\partial \psi\_\mathbf{f}
}{
\partial x\_i
}\right)^2 + \left(\frac{
\partial \psi\_\mathbf{f}
}{
\partial y\_i
}\right)^2 \right] = 0.\tag{A19}
$$

The flow-rate density of the *i*th hydraulic fracture is given as

$$q\_{\rm f}(\mathbf{x}\_{i},t) = \frac{hT\_{\rm sc}k\_{\rm i}e^{-\gamma(\psi\_{i}-\psi)}}{2p\_{\rm sc}T} \left[ \left. \frac{\partial\psi(\mathbf{x}\_{i},y\_{i},t)}{\partial y\_{i}} \right|\_{y\_{i}=\frac{w\_{\rm i}}{2}} - \left. \frac{\partial\psi(\mathbf{x}\_{i},y\_{i},t)}{\partial y\_{i}} \right|\_{y\_{i}=-\frac{w\_{\rm i}}{2}} \right],\tag{A20}$$

where 0 < *xi* < *L*f*i*.

The interface boundary condition between the wellbore and the *i*th hydraulic fracture is

$$\left| h \int\_{-\frac{w\_i}{2}}^{\frac{w\_i}{2}} \frac{T\_{sc} k\_{fi} e^{-\gamma(\psi\_i - \psi\_t)}}{2p\_{sc}T} \frac{\partial \Psi\_t(\mathbf{x}\_i, y\_i, t)}{\partial \mathbf{x}\_i} d\mathbf{y}\_i \right|\_{\mathbf{x}\_i = 0} = Q\_i(t), \tag{A21}$$

where *Qi*(*t*) is the production rate from the *i*th hydraulic fracture, which is expressed as

$$Q\_i(t) = \int\_0^{L\_{\text{fi}}} q\_\mathbf{f}(\mathbf{x}\_i, t) \, \mathbf{d}x\_i. \tag{A22}$$

*Appl. Sci.* **2019**, *9*, 1359

The total production rate of the vertical fractured well is the sum of the production rate of each hydraulic fracture, and thus we can obtain

$$\sum\_{i=1}^{n} Q\_i(t) = Q\_{\text{sc}}.\tag{A23}$$

The gas velocity should be continuous at the interface between the hydraulic fracture and reservoir, and thus the boundary conditions of the fracture surface are given as

$$\frac{T\_{\rm sc}k\_{\rm f\bar{t}}e^{-\gamma(\psi\_{i}-\psi\_{\rm t})}}{2p\_{\rm sc}T}\frac{\partial\psi\_{\rm l}(\mathbf{x}\_{i},y\_{i},t)}{\partial y\_{i}}\bigg|\_{y\_{i}=\frac{w\_{\rm l}}{2}} = \left.\frac{T\_{\rm sc}k\_{\rm i}e^{-\gamma(\psi\_{i}-\psi)}}{2p\_{\rm sc}T}\frac{\partial\psi(\mathbf{x}\_{i},y\_{i},t)}{\partial y\_{i}}\right|\_{y\_{i}=\frac{w\_{\rm l}}{2}}\tag{A24}$$

$$\frac{T\_{\rm sc}k\_{\rm fi}e^{-\gamma(\psi\_{i}-\psi\_{i})}}{2p\_{\rm sc}T}\frac{\partial\psi\_{\rm f}(\mathbf{x}\_{i},y\_{i},t)}{\partial y\_{i}}\bigg|\_{y\_{i}=-\frac{w\_{\rm f}}{2}} = \frac{T\_{\rm sc}k\_{\rm f}e^{-\gamma(\psi\_{i}-\psi)}}{2p\_{\rm sc}T}\frac{\partial\psi(\mathbf{x}\_{i},y\_{i},t)}{\partial y\_{i}}\bigg|\_{y\_{i}=-\frac{w\_{\rm f}}{2}}.\tag{A25}$$

Introducing the following expressions [16]

$$
\psi\_{\mathbf{f}} = \psi\_{\mathbf{i}} + \frac{1}{\gamma} \ln[1 - \gamma \cdot \mathbb{g}\_{\mathbf{f}}],\tag{A26}
$$

$$
\psi = \psi\_{\rm i} + \frac{1}{\gamma} \ln[1 - \gamma \cdot \xi],
\tag{A27}
$$

Equations (A19)–(A21), (A24), and (A25) can be rewritten as

$$\frac{\partial^2 \xi\_l(\mathbf{x}\_i, y\_i, t)}{\partial \mathbf{x}\_i^2} + \frac{\partial^2 \xi\_l(\mathbf{x}\_i, y\_i, t)}{\partial y\_i^2} = 0,\\ \left(0 < \mathbf{x}\_i < L\_{\text{fi}} - \frac{w\_{\text{f}}}{2} < y\_i < \frac{w\_{\text{f}}}{2}\right),\tag{A28}$$

$$q\_{\rm f}(\mathbf{x}\_{i},t) = -\frac{hT\_{\rm sc}k\_{\rm i}}{2p\_{\rm sc}T} \left[ \frac{\partial \mathfrak{J}(\mathbf{x}\_{i},y\_{i},t)}{\partial y\_{i}} \Big|\_{y\_{i}=\frac{w\_{\rm i}}{2}} - \frac{\partial \mathfrak{J}(\mathbf{x}\_{i},y\_{i},t)}{\partial y\_{i}} \Big|\_{y\_{i}=-\frac{w\_{\rm i}}{2}} \right], (0 < \mathbf{x}\_{i} < L\_{\rm fi}), \tag{A29}$$

$$\left.h\int\_{-\frac{w\_i}{2}}^{\frac{w\_i}{2}} \frac{T\_{sc}k\_{fi}}{2p\_{sc}T} \frac{\partial \xi\_t(\mathbf{x}\_i, y\_i, t)}{\partial \mathbf{x}\_i} \mathrm{d}y\_i\right|\_{\mathbf{x}\_i=0} = -Q\_i(t),\tag{A30}$$

$$\frac{T\_{\rm sc}k\_{\rm fi}}{2p\_{\rm sc}T} \frac{\partial \mathfrak{J}\_{\rm f}^{\rm r}(\mathbf{x}\_{i\prime}, y\_{i\prime}, t)}{\partial y\_{i}} \bigg|\_{y\_{i} = \frac{w\_{\rm i}}{2}} = \frac{T\_{\rm sc}k\_{\rm i}}{2p\_{\rm sc}T} \frac{\partial \mathfrak{J}\_{\rm r}^{\rm r}(\mathbf{x}\_{i\prime}, y\_{i\prime}, t)}{\partial y\_{i}} \bigg|\_{y\_{i} = \frac{w\_{\rm i}}{2}}\tag{A31}$$

$$\frac{T\_{\rm sc}k\_{\rm fi}}{2p\_{\rm sc}T} \frac{\partial \mathfrak{J}\_{\rm f}^{\rm x}(\mathbf{x}\_{i\prime}y\_{i\prime},t)}{\partial y\_{i}}\bigg|\_{y\_{i}=-\frac{w\_{\rm i}}{2}} = \frac{T\_{\rm sc}k\_{\rm i}}{2p\_{\rm sc}T} \frac{\partial \mathfrak{J}\_{\rm f}^{\rm x}(\mathbf{x}\_{i\prime}y\_{i\prime},t)}{\partial y\_{i}}\bigg|\_{y\_{i}=-\frac{w\_{\rm i}}{2}}.\tag{A32}$$

The width of the hydraulic fracture was very small, and thus the pressure changes within the *i*th hydraulic fracture in the *yi* direction can be neglected. Then the average pressure within the *i*th hydraulic fracture in the *yi* direction is introduced as follows:

$$\zeta\_{\rm f}^{\rm r}(\mathbf{x}\_{i\prime},t) = \frac{1}{w\_{\rm f}} \int\_{-\frac{w\_{\rm f}}{2}}^{\frac{w\_{\rm f}}{2}} \zeta\_{\rm f}(\mathbf{x}\_{i\prime},y\_{i\prime},t) \mathbf{d}y\_{i\prime} \ (0 < \mathbf{x}\_{i} < L\_{\rm fi}). \tag{A33}$$

With the aid of Equations (A31)–(A33), Equation (A28) becomes

$$\frac{\partial^2 \zeta\_l(\mathbf{x}\_i, t)}{\partial \mathbf{x}\_i^2} + \frac{k\_i}{w\_l k\_{\text{fi}}} \left[ \left. \frac{\partial \zeta(\mathbf{x}\_i, y\_i, t)}{\partial y\_i} \right|\_{y\_i = \frac{w\_i}{2}} - \left. \frac{\partial \zeta(\mathbf{x}\_i, y\_i, t)}{\partial y\_i} \right|\_{y\_i = -\frac{w\_i}{2}} \right] = 0, \ (0 < \mathbf{x}\_i < L\_{\text{fi}}). \tag{A34}$$

Substituting Equation (A29) into Equation (A34) yields

$$\frac{\partial^2 q\_\mathbf{f}(\mathbf{x}\_i, \mathbf{t})}{\partial \mathbf{x}\_i^2} - \frac{2p\_{\text{sc}}T}{w\_l k\_\text{fl} hT\_{\text{sc}}} q\_\mathbf{f}(\mathbf{x}\_i, \mathbf{t}) = \mathbf{0}, \ (0 < \mathbf{x}\_i < L\_{\text{fi}}). \tag{A35}$$

According to Equation (A33), Equation (A30) is expressed as

$$\left. \text{div}\_{\mathbf{f}} \frac{T\_{\text{sc}} k\_{\text{fi}}}{2p\_{\text{sc}}T} \frac{\partial f\_{\text{f}}(\mathbf{x}\_{i}, t)}{\partial \mathbf{x}\_{i}} \right|\_{\mathbf{x}\_{i} = \mathbf{0}} = -Q\_{i}(t). \tag{A36}$$

Introducing the dimensionless variables, Equations (A26), (A35), (A36) are rewritten as follows:

$$
\psi\_{\rm fD} = -\frac{1}{\gamma\_{\rm D}} \ln[1 - \gamma\_{\rm D} \cdot \xi\_{\rm fD}],
\tag{A37}
$$

$$\frac{\partial^2 \mathbf{f}\_{\mathrm{fD}}(\mathbf{x}\_{\mathrm{iD}}, t\_{\mathrm{D}})}{\partial \mathbf{x}\_{\mathrm{iD}}^2} - \frac{2\pi}{\mathbf{C}\_{\mathrm{fD}}} q\_{\mathrm{fD}}(\mathbf{x}\_{\mathrm{iD}}, t\_{\mathrm{D}}) = 0, \ (0 < \mathbf{x}\_{\mathrm{iD}} < L\_{\mathrm{fD}}), \tag{A38}$$

$$
\left.\frac{\partial \xi\_{\rm fD}(x\_{\rm iD}, t\_{\rm D})}{\partial x\_{\rm iD}}\right|\_{x\_{\rm iD}=0} = -\frac{2\pi}{\mathcal{C}\_{\rm fD}} Q\_{\rm iD}(t\_{\rm D}).\tag{A39}
$$

The dimensionless expressions of Equations (A22) and (A23) are given as

$$Q\_{\rm iD}(t\_{\rm D}) = \int\_{0}^{L\_{\rm iD}} q\_{f\rm D}(x\_{\rm iD}, t\_{\rm D}) \mathrm{d}x\_{\rm iD} \tag{A40}$$

$$\sum\_{j=1}^{n} Q\_{j\mathbf{D}}(t\mathbf{p}) = 1.\tag{A41}$$

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Parametrical Non-Complex Tests to Evaluate Partial Decentralized Linear-Output Feedback Control Stabilization Conditions from Their Centralized Stabilization Counterparts**

#### **Manuel De la Sen 1,\* and Asier Ibeas <sup>2</sup>**


Received: 18 March 2019; Accepted: 22 April 2019; Published: 26 April 2019

**Abstract:** This paper formulates sufficiency-type linear-output feedback decentralized closed-loop stabilization conditions if the continuous-time linear dynamic system can be stabilized under linear output-feedback centralized stabilization. The provided tests are simple to evaluate, while they are based on the quantification of the sufficiently smallness of the parametrical error norms between the control, output, interconnection and open-loop system dynamics matrices and the corresponding control gains in the decentralized case related to the centralized counterpart. The tolerance amounts of the various parametrical matrix errors are described by the maximum allowed tolerance upper-bound of a small positive real parameter that upper-bounds the various parametrical error norms. Such a tolerance is quantified by considering the first or second powers of such a small parameter. The results are seen to be directly extendable to quantify the allowed parametrical errors that guarantee the closed-loop linear output-feedback stabilization of a current system related to its nominal counterpart. Furthermore, several simulated examples are also discussed.

**Keywords:** Output-feedback; centralized control; decentralized control; closed-loop stabilization

#### **1. Introduction**

Control systems are very important in real world applications and, therefore, they have been investigated exhaustively concerning their properties of stability, stabilization, controllability control strategies etc. See, for instance, [1–4] and references therein. Some extra constraints inherent to some systems, like solution positivity in the case of biological systems or human migrations or the needed behavior robustness against parametrical changes of disturbance actions add additional complexity to the related investigations and need the use of additional mathematical or engineering tools for the research development, [5–7]. A large variety of modeling and design tools have to be invoked and developed in the analysis depending on the concrete systems under study and their potential applications as, for instance, the presence of internal and external delays, discretization, dynamics modeling based on fractional calculus, the existence of complex systems with interconnected subsystems, [8–13], hybrid coupled continuous/digital tandems, nonlinear systems and optimization and estimation techniques [14–19] as well as robotic and fuzzy-logic based systems, [20,21]. In particular, decentralized control is a useful tool for controlling dynamic systems by cutting some links between the dynamics coupling a set of subsystems integrated in the whole system at hand. It is claimed to keep the main properties related to the use of centralized control such as stability, controllability, observability, etc. In summary, a centralized controller keeps all the information on the system and coupling links as

available to the control designer while decentralized control ignores some of such information or even cuts on occasions some of coupling signals between the various subsystems integrated in the whole system at hand. It can be pointed out that the stability studies are often performed trough Lyapunov theory which requires to find a Lyapunov function (see [20,21] and some references therein). It turns out that, if the neglected couplings are strong and are not taken into account by the controller, the stabilization and other properties such as the controllability can become lost. The use of decentralized stabilization and control tools is of interest when the whole system has physically separated subsystems that require the implementation of local control actuators but the control has to be global for the whole system. An ad-hoc example provided in [2–4,19] where decentralized control is of a great design interest is the case of several coupled cascade hydroelectric power plants allocated in the same river but separated far away from each other. It has to be pointed out that the term "decentralized control" versus "centralized control" refers to the eventual cut of links of the shared information between tandems of integrated subsystems, or coupling signals between them, to be controlled rather than to the physical disposal of the controller. In other words, if the whole controller keeps and uses all the information on measurable outputs and control components design to implement the control law, such a control is considered to have a centralized nature even if its various sub-control stations are not jointly allocated. It is a common designer´s basic idea in mind for complex control designs to try to minimize the modeling designs and computational loads without significantly losing the system´s performance and its essential properties. For instance, in [8], the dynamic characteristic of a discrete-time system is given as an extended state space description in which state variables and output tracking error are integrated while they are regulated independently. The proposed robust model predictive control is much simpler than the traditional versions since the information of the upper and lower bounds of the time-varying delay are used for design purposes. On the other hand, in [9], a control law might be synthesized for a hydropower plant with six generation units working in an alternation scheme. To assess the behavior of the controlled system, a model of such a nonlinear plant is controlled by a fractional proportional/integral/derivative control device through a linearization of its set points, the fractional part being relevant in the approach on the controller derivative actions. In addition, a set of applied complex control problems are studied, for instance, in [10–16] with the aim of giving different ad hoc simplification tools to deal with the appropriate control methodologies. In particular, a decentralized control approach is proposed in [16].

In this paper, the decentralized control design versus its decentralized control counterpart, under eventual output linear feedback, are studied from the point of view of the amount of information that can be lost or omitted in terms of the total or partial knowledge of the coupled dynamics between subsystems necessary in the decentralized case to keep the closed-loop stability. The study is made by using the information on the worst-case deviation, in terms of norms, between the respective matrices of open-loop dynamics and the respective controller gains under which the closed-loop stability is kept. This paper is organized as follows. The problem statement is given in Section 2 while the main stabilization results of the paper are provided in Section 3. The proofs of some of the results of Section 3, which are technically involved, but conceptually simple, are distributed in various technical auxiliary that are given in Appendices A and B. It is claimed to give a non-complex method to test the feasibility of the implementation of decentralized control and conditions for its design, which be a fast and simpler stability test compared to Lyapunov stability theory [20,21], for instance, under a partial removal of information or physical cuts of links of coupling dynamics between the various subsystems or state, control and output components. Section 4 discusses several examples and, finally, the concluding remarks end the paper.

*Notation*

$$\overline{n} = \{1, 2, \cdots, n\},$$

$$\mathcal{R}\_{+} = \{z \in \mathcal{R} \,:\, z > 0\} \,; \, \mathcal{R}\_{0+} = \{z \in \mathcal{R} \,:\, z \ge 0\}.$$

*sp*(*A*) and det(*A*) are the spectrum and determinant of *<sup>A</sup>* <sup>∈</sup> *<sup>R</sup>n*×*n*, respectively. For *<sup>A</sup>* <sup>∈</sup> *<sup>R</sup>q*×-, being in general rectangular, *A* denotes any unspecified norm of *A*, *A*<sup>2</sup> denotes the -<sup>2</sup> or spectral norm of a matrix *A*, ρ(*A*) denotes its spectral radius, and .∞ denotes the *H*∞-norm of a stable real rational transfer matrix or function, *Iq* denotes the *<sup>q</sup>*th identity matrix, and *<sup>i</sup>* <sup>=</sup> <sup>√</sup> −1 is the complex imaginary unit.

Let *<sup>A</sup>* <sup>∈</sup> *<sup>R</sup>n*×*<sup>n</sup>* be symmetric. Then, <sup>λ</sup>max(*A*) and <sup>λ</sup>min(*A*) are, respectively, the maximum and minimum eigenvalues of *A*.

*Mn*×*<sup>n</sup> <sup>E</sup>* is the set of Metzler matrices (any off-diagonal entry is non-negative) of *n*th order.

*Zn*×*<sup>n</sup>* is the set of *Z*-matrices (any off-diagonal entry is non-positive) of *n*th order.

*Mn*×*<sup>n</sup>* is the set of *M*-matrices (*Z*-matrices which are stable or critically stable) of *n*th order.

Assume that *A* = -*Aij* , *B* = -*Bij* <sup>∈</sup> *<sup>R</sup>n*×*n*. Then, the notations *<sup>A</sup>B*, *<sup>A</sup> <sup>B</sup>* and *<sup>A</sup> <sup>B</sup>*, are, respectively, equivalent to *B*≺*A*, *B* ≺ *A* and *B* ≺≺ *A*, meaning that *Aij* ≥ *Bij*, *Aij* ≥ *Bij* (and *B* - *A*) and *Aij* > *Bij*; ∀*i*, *j* ∈ *n*, respectively. In particular, *A*0, *A* 0 and *A* 0 are reworded as *A* is non-negative, positive and strictly positive, respectively, and *A*≺0, *A* ≺ 0 and *A* ≺≺ 0 are reworded as *A* is non-positive, negative and strictly negative, respectively.

#### **2. Problem Statement**

Consider the following linear and time-invariant system under linear output-feedback centralized control: .

$$\dot{\mathbf{x}}\_{\mathcal{L}}(t) = A\_{\mathcal{L}} \mathbf{x}\_{\mathcal{L}}(t) + B\_{\mathcal{C}} \boldsymbol{u}\_{\mathcal{C}}(t);\ \mathbf{x}\_{\mathcal{C}}(0) = \mathbf{x}\_{\mathcal{C}0} \tag{1}$$

$$\mathbf{y}\_{\mathfrak{c}}(t) = \mathbf{C}\_{\mathfrak{c}} \mathbf{x}\_{\mathfrak{c}}(t) + D\_{\mathfrak{c}} u\_{\mathfrak{c}}(t) \tag{2}$$

$$u\_{\mathcal{C}}(t) = \mathcal{K}\_{\mathcal{C}} y\_{\mathcal{C}}(t) \tag{3}$$

where *xc*(*t*) <sup>∈</sup> *<sup>R</sup><sup>n</sup>* is the state vector; *uc*(*t*) <sup>∈</sup> *<sup>R</sup><sup>m</sup>* is the centralized control vector; *yc*(*t*) <sup>∈</sup> *<sup>R</sup><sup>p</sup>* is the output; *Ac*, *Bc*, *Cc* and *Dc* are the system, control, output and input–output interconnection matrices, respectively, of orders being compatible with the dimensions of the above vectors; and *Kc* <sup>∈</sup> *<sup>R</sup>m*×*<sup>p</sup>* is the control matrix. If the system runs in a decentralized control context, we have:

$$\dot{\mathbf{x}}\_d(t) = A\_d \mathbf{x}\_d(t) + B\_d \boldsymbol{u}\_d(t) \; ; \; \mathbf{x}\_d(0) = \mathbf{x}\_{d0} \tag{4}$$

$$\mathbf{y}\_d(t) = \mathbb{C}\_d \mathbf{x}\_d(t) + D\_d \mathbf{u}\_d(t) \tag{5}$$

$$u\_d(t) = \mathcal{K}\_d y\_d(t) \tag{6}$$

where *xd*(*t*) <sup>∈</sup> *<sup>R</sup><sup>n</sup>* is the state vector; *ud*(*t*) <sup>∈</sup> *<sup>R</sup><sup>m</sup>* is the centralized control vector; *yd*(*t*) <sup>∈</sup> *<sup>R</sup><sup>p</sup>* is the output; *Ad*, *Bd*, *Cd* and *Dd* are the system, control, output and input–output interconnection matrices, respectively, of orders being compatible with the dimensions of the above vectors; and *Kd* <sup>∈</sup> *<sup>R</sup>m*×*<sup>p</sup>* is the control matrix.

Basically, the differences between centralized and decentralized controls are as follows:


centralized ones and, roughly speaking, to a have a "more diagonal" or "sparser" structure than their centralized counterparts *Bc*, *Cc* and *Dc*. If the parameterization of the system (or dynamics) matrix is available to the designer, then some off-diagonal block matrices of *Ac* could be zeroed or simply re-disposed in a more sparse structure to construct *Ad*.

(3) The only strictly necessary condition for the system to be subject to partially (or, respectively, fully) decentralized control is that some (or, respectively, all) of the off-diagonal entries of *Kd* are forced to be zero even if the system, control, output and interconnection matrices remain identical in Equations (4) and (5) with respect to Equations (1) and (2).

**Assumption 1.** *The system in Equations (1) and (2) is linear output-feedback stabilizable via some centralized control law (Equation (3)).*

Note that Assumption 1 does not hold if the open-loop system in Equations (1) and (2) has unstable or critically stable fixed modes that cannot be removed via linear feedback.

**Proposition 1.** *If Assumption 1 holds, then there exists a centralized stabilizing controller gain Kc* <sup>∈</sup> *<sup>R</sup>m*×*<sup>p</sup> such that the matrices* (*Im* <sup>−</sup> *KcDc*) *and* - *Ip* − *DcKc are non-singular, thus the closed-loop centralized control system is solvable and given by:*

$$\dot{\mathbf{x}}\_{\varepsilon}(t) = \left(A\_{\varepsilon} + B\_{\varepsilon}(I\_{\mathfrak{m}} - K\_{\mathfrak{c}}D\_{\mathfrak{c}})^{-1}K\_{\mathfrak{c}}\mathbb{C}\_{\varepsilon}\right)\mathbf{x}\_{\varepsilon}(t);\ \mathbf{x}\_{\varepsilon}(0) = \mathbf{x}\_{\varepsilon 0} \tag{7}$$

$$y\_c(t) = \left(I\_p - D\_c K\_c\right)^{-1} \mathbb{C}\_c x\_c(t) \tag{8}$$

*and asymptotically stable for any given xc*<sup>0</sup> <sup>∈</sup> *<sup>R</sup><sup>n</sup> under the generated control law:*

$$
\mu\_{\mathfrak{c}}(t) = \left(I\_{\mathfrak{m}} - K\_{\mathfrak{c}}D\_{\mathfrak{c}}\right)^{-1} K\_{\mathfrak{c}} \mathbb{C}\_{\mathfrak{c}} \mathbf{x}\_{\mathfrak{c}}(t) \tag{9}
$$

*that is, the polynomial p*(*s*) = det- *sIn* <sup>−</sup> *Ac* <sup>−</sup> *Bc*(*Im* <sup>−</sup> *KcDc*)−1*KcCc is Hurwitz.*

**Proof.** The replacement of Equation (3) into Equations (1) and (2) yields Equation (7)–(9). Since the (1) and (2) is linear output stabilizable, a stabilizing controller gain *Kc* <sup>∈</sup> *<sup>R</sup>m*×*<sup>p</sup>* has to exist such that (7)–(9) are solvable and the closed-loop dynamics is stable. -

**Assumption 2.** *The system in Equations (4) and (5) is linear output-feedback stabilizable via some decentralized control law (Equation (6)).*

In the same way as Proposition 1, we get the following result:

**Proposition 2.** *If Assumption 2 holds, then there exists a decentralized stabilizing controller gain Kd* <sup>∈</sup> *<sup>R</sup>m*×*<sup>p</sup> such that the matrices* (*Im* <sup>−</sup> *KdDd*) *and* - *Ip* − *DdKd are non-singular, thus the closed-loop decentralized control system is solvable and given by:*

$$\dot{\mathbf{x}}\_d(\mathbf{t}) = \left(\mathbf{A}\_d + B\_d (I\_m - \mathbf{K}\_d \mathbf{D}\_d)^{-1} \mathbf{K}\_d \mathbf{C}\_d\right) \mathbf{x}\_d(\mathbf{t});\ \mathbf{x}\_d(\mathbf{0}) = \mathbf{x}\_{d0} \tag{10}$$

$$y\_d(t) = \left(I\_p - D\_d K\_d\right)^{-1} \mathbb{C}\_d \mathbf{x}\_d(t) \tag{11}$$

*and asymptotically stable for any given xd*<sup>0</sup> <sup>∈</sup> *<sup>R</sup><sup>n</sup> under the control law:*

$$
\mu\_d(t) = (I\_m - K\_d D\_d)^{-1} K\_d \mathbb{C}\_d \mathbf{x}\_d(t) \tag{12}
$$

*that is, the polynomial p*(*s*) = det- *sIn* <sup>−</sup> *Ad* <sup>−</sup> *Bd*(*Im* <sup>−</sup> *KdDd*)−<sup>1</sup>*KdCd is Hurwitz.* **Proposition 3.** *Assume that Ad* = *Ac, Bd* = *Bc, Cd* = *Cc and Dd* = *Dc, and that the system in Equations (1) and (2) is not linear output-feedback stabilizable via some centralized control law (Equation (3)). Then, it is not stabilizable under any linear output-feedback decentralized control (Equation (6)) either.*

**Proof.** Obviously, if there is no completely free-design matrix *Kc* that stabilizes Equations (1) and (2), then there is no *Kd* with at least a forced zero off-diagonal entry that stabilizes it since *Kd* has extra design constraints related to *Kc*. -

It can be pointed out that decentralized control has also been proved to be useful in applications. For instance, an integral-based event-triggered asymptotic stabilization of a continuous-time linear system is studied in [17] by considering actuator saturation and observer-based output feedback are considered. In the proposed scheme, the sensors and actuators are implemented in a decentralized manner and a type of Zeno-free decentralized integral-based event condition is designed to guarantee the asymptotic stability of the closed-loop systems. On the other hand, two decentralized fuzzy logic-based control schemes with a high-penetration non-storage wind-diesel system are studied in [18] for small power system with high-penetration wind farms. In addition, several examples concerning decentralized control are described in [4] to illustrate the theoretical design analysis. A typical described case is that of tandems of electrical power system with a tandem disposal on the same river which are not physically nearly allocated. The next section discusses some simple sufficiency-type conditions which ensure that, provided that the system is stabilizable under linear output-feedback centralized control, it is also stabilizable under decentralized control in two cases: (a) the system matrix remains identical but the other parameterization matrices can eventually vary; and (b) the system matrix can vary as well in the decentralized case with respect to the centralized one. A result elated to the maintenance of the stability of a matrix under an additive matrix perturbation is summarized through a set of sufficiency-type conditions simple to test in Theorem A1. Theorem A2 proves sufficiency-type for the stability of the matrix function *<sup>A</sup>*(*t*) <sup>=</sup> *<sup>A</sup>*<sup>0</sup> <sup>+</sup> *<sup>A</sup>*(*t*) with *<sup>A</sup>*<sup>0</sup> stable and *<sup>A</sup>*(*t*) being time-varying. Appendix B includes calculations and auxiliary results to quantify the tolerance to cut some dynamics links between subsystems, state components or control centers or components while keeping the closed-loop stability of the whole coupled system. The results of Appendices A and B are used in the proofs of the main results in the next section.

#### **3. Main Results**

The first set of technical results which follow are concerned with centralized and decentralized control stabilizability.

**Assertion 1.** *A necessary and su*ffi*cient condition for the system to be linear state-stabilizable via some centralized control law is that rank*[*sIn* − *Ac Bc*] = *n for all Res* ≥ 0*.*

**Proof.** Assume that *rank*[*sIn* − *Ac Bc*] = *n*<sup>1</sup> < *n* for some *Res* ≥ 0. Then, there is some Laplace transform *x*ˆ*T*(*s*), *u*ˆ*T*(*s*) *T* = *Lap xT*(*t*), *uT*(*t*) *<sup>T</sup>* - 0 such that [*sIn* − *Ac Bc*] *x*ˆ*T*(*s*), *u*ˆ*T*(*s*) *<sup>T</sup>* = 0 for some *Res* ≥ 0 and [*sIn* − *Ac Bc*] *In Kc <sup>x</sup>*ˆ(*s*) = 0 for any *Kc* <sup>∈</sup> *<sup>R</sup>n*×*<sup>n</sup>* and some *<sup>x</sup>*ˆ(*s*) - 0 with *Res* ≥ 0 since *rank* [*sIn* − *Ac Bc*] *In Kc* <sup>≤</sup> min *rank*[*sIn* <sup>−</sup> *Ac Bc*] , *rank In Kc* ≤ min(*n*<sup>1</sup> , *n*) ≤ *n*<sup>1</sup> < *n* for some *Res* ≥ 0. Therefore, the closed-loop system has some unstable or critically stable solution for any given (centralized) control gain. This proves the necessary part. Sufficiency follows since, if *rank*[*sIn* <sup>−</sup> *Ac Bc*] = *<sup>n</sup>*, then *<sup>x</sup>*ˆ(*s*) <sup>≡</sup> 0 for all *Res* <sup>≥</sup> 0 and some *Kc* <sup>∈</sup> *<sup>R</sup>n*×*<sup>n</sup>* which can be found so that min *rank*[*sIn* <sup>−</sup> *Ac Bc*] , *rank In Kc* = min(*n*, *n*) = *n*. -

**Assertion 1.** is a particular adapted ad-hoc test for stabilizability of the celebrated Popov–Belevitch –Hautus rank controllability test [6]. Note that, if there exist unstable or critically stable fixed modes (i.e., those present in the open-loop system that cannot be removed via feedback control), then neither centralized nor decentralized stabilizing control laws can be synthesized. Note that the stabilizability rank test of Assumption 1 can only be evaluated for the critically and unstable eigenvalues of *Ac* instead for all the open right-hand complex plane. In all the remaining points of such a plane, the test always gives a full rank of the tested matrix. The parallel controllability test should always be applied in the same matrix to any eigenvalues of *Ac*.

**Assertion 2.** *A necessary condition for the system to be linear state-stabilizable via some partially or totally decentralized control is that it be stabilizable via centralized control (i.e., Assertion 1 holds).*

**Proof.** It is obvious from Assertion 1 since any gain *Kd* used for centralized or decentralized is sparser than a centralized gain counterpart so that the proof follows from Assertion 1. -

**Assertion 3.** *A necessary condition for the system to be linear output-stabilizable via some partially or totally decentralized control is that it be stabilizable via centralized control (i.e., that Assertion 1 holds).*

**Proof.** It is obvious from Assertions 1 and 2 that, when replacing *Kc* <sup>→</sup> (*Im* <sup>−</sup> *KcDc*)−1*KcCc* and *Kd* <sup>→</sup> (*Im* <sup>−</sup> *KdDd*)−<sup>1</sup>*KdCd* (see Equations (9) and (12)), the second replacement happens under sparser parameterizations. -

Now, consider the closed-loop system matrices from Equations (7) and (10) for the case *A* = *Ac* = *Ad*.

$$A\_{\mathfrak{C}} = A + B\_{\mathfrak{C}} (I\_{\mathfrak{M}} - K\_{\mathfrak{C}} D\_{\mathfrak{C}})^{-1} K\_{\mathfrak{C}} \mathbb{C}\_{\mathfrak{C}}; \\ A\_{\mathrm{dc}} = A + B\_{\mathfrak{d}} (I\_{\mathfrak{m}} - K\_{\mathfrak{d}} D\_{\mathfrak{d}})^{-1} K\_{\mathfrak{d}} \mathbb{C}\_{\mathfrak{d}}$$

with its parametrical error being:

$$\overline{A}\_{d\varepsilon} = A\_{c\varepsilon} - A\_{d\varepsilon} = B\_{\varepsilon}(I\_m - K\_{\varepsilon}D\_{\varepsilon})^{-1}K\_{\varepsilon}C\_{\varepsilon} - B\_{d}(I\_m - K\_{d}D\_{d})^{-1}K\_{d}C\_{d}$$

A first main technical result follows:

**Theorem 1.** *If Assumption 1 holds, assume also that Kc* <sup>∈</sup> *<sup>R</sup>m*×*<sup>p</sup> is a centralized linear output-feedback stabilizing controller gain such that the resulting closed-loop system matrix Acc* <sup>∈</sup> *<sup>R</sup>n*×*<sup>n</sup> has a stability abscissa* (−ρ*cc*) < 0*. Then, the following properties hold:*

*(i) Adc* <sup>∈</sup> *<sup>R</sup>n*×*<sup>n</sup> is a closed-loop stability matrix under a linear output-feedback stabilizing controller gain Kd* <sup>∈</sup> *<sup>R</sup>m*×*<sup>p</sup> if any of the subsequent su*ffi*ciency-type conditions holds:*

*(1) The H*∞*-norm of* (*sIn* <sup>−</sup> *Acc*)−1*Adc satisfies* (*sIn* <sup>−</sup> *Acc*)−1*Adc*∞ <sup>&</sup>lt; <sup>1</sup>*, (2) Adc*<sup>2</sup> <sup>&</sup>lt; 1/ sup <sup>ω</sup>∈*R*0<sup>+</sup> (*i*ω*In* <sup>−</sup> *Acc*)−12*. Other alternative su*ffi*ciency-type conditions to Conditions 1 and 2 for the stability of Adc are:* 

$$\begin{array}{l}(3)\ \rho\big(A\_{\rm c}^{-1}A\_{\rm dc}\big)<1,\\(4)\ \|A\_{\rm c}^{-1}\overline{A}\_{\rm dc}\|\_{2}<1,\\(5)\ \|\overline{A}\_{\rm dc}\|\_{2}<1/\|A\_{\rm c}^{-1}\|\_{2},\ \text{that is,}\ \lambda\_{\rm max}\left(\overline{A}\_{\rm dc}^{T}\overline{A}\_{\rm dc}\right)<\lambda\_{\rm min}\left(A\_{\rm c}^{T}A\_{\rm c}\right).\\\ \text{in the following particular cases:}\\(a)\ \begin{array}{l}(a)\ A\_{\rm c}<0 \text{ and }\overline{A}\_{\rm dc}>A\_{\rm c};\text{and}\\(b)\ A\_{\rm c}=\left(A\_{\rm c}\overline{c}\_{ij}\right)\in\mathcal{M}\_{\rm E}^{\rm p\times n}\text{ and }\overline{A}\_{\rm dc}=\left(\overline{A}\_{\rm c}\overline{c}\_{ij}\right)\text{fullths }\overline{A}\_{\rm dcij}\le A\_{\rm c};j,\forall i,j\in\overline{n}.\end{array}\end{array}$$

*(ii) Assume that Property (i) holds and that the number of inputs and outputs are identical, i.e., p* = *m, and decompose both the controller gain matrices as sums of their diagonal and o*ff*-diagonal parts leading to Kc* = *Kcd* + *Kcod and Kd* = *Kdd* + *Kdod, thus K* <sup>=</sup> *Kc* <sup>−</sup> *Kd* <sup>=</sup> (*Kcd* <sup>−</sup> *Kdd*) <sup>+</sup> (*Kcod* <sup>−</sup> *Kdod*)*. Then, the system is stabilizable under partially decentralized control linear output-feedback control in the sense that Equations (4) and (5), is asymptotically stable under a control law (Equation (6)), if Kd* <sup>∈</sup> *<sup>R</sup>m*×*<sup>p</sup> is such that, if there*

*is at least one non-diagonal zero entry in at least one of its rows in the o*ff*-diagonal controller error matrix K od* <sup>=</sup> *Kcod* <sup>−</sup> *Kdod. If <sup>K</sup> od* <sup>=</sup> <sup>0</sup>*, then the system is stabilizable under decentralized control.*

**Proof.** Property (i) is a direct translation of the results of Theorem A1 in Appendix A to the closed-loop system matrices. Property (ii) holds if Property (i) holds with an off-diagonal controller error matrix between the centralized a decentralized controller gain that has at least one non-diagonal zero at some row (or its identically zero) so that a feedback from some crossed output to some of the inputs is not provided to the control law for closed-loop stabilization the stabilization. -

The following result follows for the time-varying case from Theorem 1 and Theorem A2:

**Corollary 1.** *Assume that Adc*(*t*) *and then <sup>A</sup>dc*(*t*) *are everywhere piecewise-continuous time-varying. Then, Theorem 1 still holds if Condition 1 is replaced with kcc* <sup>ρ</sup>*cc* sup 0≤τ<*t Adc*(τ) <sup>&</sup>lt; <sup>1</sup>*;* <sup>∀</sup>*<sup>t</sup>* <sup>∈</sup> *<sup>R</sup>*0<sup>+</sup> *with kcc* <sup>≥</sup> <sup>1</sup> *and* <sup>ρ</sup>*cc* <sup>&</sup>gt; <sup>0</sup> *being real constants such that e*−*Acct* <sup>≤</sup> *kcce*−ρ*cct ;* ∀*t* ∈ *R*0+*.*

**Remark 1.** *Theorem 1 (ii) has been stated for the case m* = *p. Note that the case m* > *p (i.e., there are more inputs than outputs) is irrelevant for the stabilization from the strict algebraic point of view since the* (*m* − *p*) *extra inputs would be redundant. In the case that m* ≤ *p, Theorem 1 (ii) might be directly generalized to a subsystem´s decomposition philosophy if a number <sup>q</sup>* <sup>≤</sup> *<sup>m</sup> of subsystems of inputs and outputs* - *uT* <sup>1</sup> , *uT* <sup>2</sup> , ··· , *uT q T and* - *yT* <sup>1</sup> , *yT* <sup>2</sup> , ··· , *<sup>y</sup><sup>T</sup> q <sup>T</sup> with ui* <sup>∈</sup> *<sup>R</sup>mi , yi* <sup>∈</sup> *<sup>R</sup>pi ;* <sup>∀</sup>*<sup>i</sup>* <sup>∈</sup> *n with p* <sup>=</sup> *<sup>q</sup> <sup>i</sup>*=<sup>1</sup> *pi and m* <sup>=</sup> *<sup>q</sup> <sup>i</sup>*=<sup>1</sup> *mi.*

**Remark 2.** *Theorem 1 can be easily generalized to cases when some dynamics transmission links between state, input or output components (or subsystems, in general) can be suppressed by manipulation. In more general cases, it is possible to extend Theorem 1 to combinations of the subsequent situations with the matrix decompositions having the same sense (in the various modified contexts) as that of Theorem 1 (ii):*

• **Case 1**. Suppression of some transmission links between the coupled open-loop dynamics by examining the decompositions: *Ac* <sup>=</sup> *Acd* <sup>+</sup> *Acod*, *Ad* <sup>=</sup> *Add* <sup>+</sup> *Adod*, and *<sup>A</sup>* <sup>=</sup> *Ac* <sup>−</sup> *Ad* <sup>=</sup> (*Acd* − *Add*) + (*Acod* − *Adod*).

(*a*) If there is at least one non-diagonal zero entry in at least one of its rows in the off-diagonal controller error matrix *<sup>A</sup>od* <sup>=</sup> *Acod* <sup>−</sup> *Adod* which is not a corresponding zero in *Acod*; and (*b*) if there is at least one non-diagonal zero entry in at least one of its rows in the off-diagonal controller error matrix *K od* <sup>=</sup> *Kcod* <sup>−</sup> *Kdod*, then the closed-loop system is stabilizable under a partial decentralized control even if some links of the dynamics between crossed components are cut if Theorem 1 (ii) holds. If only Condition a is addressed, then the system is stabilizable by centralized control when cutting certain transmission links between coupled dynamics in the open-loop system. This idea can be extended to total decentralized control for a purely diagonal open-loop system´s dynamics under full zeroing of the off-diagonal corresponding error dynamics. It can be also generalized to the "ad hoc" decompositions between subsystems. Other cases with similar interpretations in the new contexts are:


**Problem 1.** *Find a stabilizing decentralized family of control gains by assuming that Ad* = *Ac, such that <sup>A</sup>* <sup>=</sup> *Ac* <sup>−</sup> *Ad* <sup>=</sup> <sup>0</sup> *and Assumption 1 holds with Kc being a stabilizing centralized controller gain.*

The following more general result for the eventual case *<sup>A</sup>* - 0 (that is *Ac* - *Ad* eventually), follows from Theorem 1, Theorem A1 and Theorem A2 and Lemmas B1, B2 and B3:

**Theorem 2.** *Define the following error matrices between the centralized and decentralized system parameterizations:*

$$
\widetilde{A} = A\_{\mathfrak{c}} - A\_{\mathfrak{d}'} ; \widetilde{B} = B\_{\mathfrak{c}} - B\_{\mathfrak{d}'} ; \widetilde{\mathbb{C}} = \mathbb{C}\_{\mathfrak{c}} - \mathbb{C}\_{\mathfrak{d}'} ; \widetilde{D} = D\_{\mathfrak{c}} - D\_{\mathfrak{d}'} ; \widetilde{K} = K\_{\mathfrak{c}} - K\_{\mathfrak{d}} \tag{13}
$$

*such that <sup>A</sup>* ≤ σ*A*ε*, B* ≤ σ*B*ε*, <sup>C</sup>* ≤ σ*c*ε*, <sup>D</sup>* ≤ σ*D*<sup>ε</sup> *and <sup>K</sup>* ≤ σ*K*<sup>ε</sup> *for some* <sup>ε</sup> <sup>∈</sup> *<sup>R</sup>*0<sup>+</sup> *and given* σ*A,* σ*B,* σ*C,* σ*<sup>D</sup>* <sup>∈</sup> *<sup>R</sup>*+*. Assume that:*

*(1) Assumption 1 holds;*


$$\overline{\varepsilon} = \frac{\sqrt{\left(\overline{\sigma}\_{\mathrm{D}} \|\mathcal{K}\_{\mathrm{c}}\| + \overline{\sigma}\_{\mathrm{K}} \|\mathcal{D}\_{\mathrm{c}}\|\right)^{2} + 4\overline{\sigma}\_{\mathrm{D}}\overline{\sigma}\_{\mathrm{K}} / \|(I\_{\mathrm{m}} - \mathcal{K}\_{\mathrm{c}}\mathcal{D}\_{\mathrm{c}})^{-1}\|} - (\overline{\sigma}\_{\mathrm{D}} \|\mathcal{K}\_{\mathrm{c}}\| + \overline{\sigma}\_{\mathrm{K}} \|\mathcal{D}\_{\mathrm{c}}\|)} \,\,\,\,\widetilde{\varepsilon}\_{1} = \frac{\left(\overline{\sigma}\_{\mathrm{D}} \|\mathcal{K}\_{\mathrm{c}}\| + \overline{\sigma}\_{\mathrm{K}} \|\mathcal{D}\_{\mathrm{c}}\|\right)}{2\|\mathcal{K}\_{\mathrm{c}}\|\overline{\sigma}\_{\mathrm{D}} + \overline{\varepsilon}\overline{\sigma}\_{\mathrm{D}} \overline{\sigma}\_{\mathrm{K}} + \|\mathcal{D}\_{\mathrm{c}}\| \overline{\sigma}\_{\mathrm{K}}\|} \,\,\,\,\widetilde{\varepsilon}\_{2} = 1 / \left(\overline{\overline{a}\_{\mathrm{dc}}} \sup\_{\omega \in \mathcal{R}\_{\mathrm{0}+}} \|(i\omega I\_{\mathrm{n}} - A\_{\mathrm{c}c})^{-1}\|\_{2}\right).$$

*where*

$$\widetilde{a}\_{\mathrm{dc}} = (1 - \|K\_{\mathrm{c}}D\_{\mathrm{c}}\|)^{-1} \times \left[\widetilde{\sigma}\_{B} \|K\_{\mathrm{c}}\mathbb{C}\_{\mathrm{c}}\| + \left(\|B\_{\mathrm{c}}\| + \|K\_{\mathrm{c}}\mathbb{C}\_{\mathrm{c}}\|\right) \left(\|K\_{\mathrm{c}}\| \widetilde{\sigma}\_{\mathrm{C}} + \|\operatorname{C}\_{\mathrm{c}}\| \widetilde{\sigma}\_{\mathrm{K}} + \operatorname{C}\left[\|K\_{\mathrm{c}}\| \widetilde{\sigma}\_{\mathrm{D}} + \|D\_{\mathrm{c}}\| \widetilde{\sigma}\_{\mathrm{K}}\right]\right)\right]$$

*where the non-negative real constant C is given in Equation (A17). Then, the following properties hold: (i) If* σ*<sup>A</sup>* <sup>=</sup> <sup>0</sup> *(that is, Ac* <sup>=</sup> *Ad), then Adc is stable and* <sup>ε</sup> <sup>∈</sup> [0 , <sup>ε</sup><sup>∗</sup> ]*. (ii) If <sup>A</sup>*<sup>2</sup> <sup>≤</sup> σ*A*ε*, then Adc is stable and* ε ∈ [0 , ε∗] *where* ε∗ = min (1 , ε, ε<sup>1</sup> , ε <sup>2</sup>) *and* ε <sup>2</sup> = 1/ ⎛ ⎜⎜⎜⎜⎝ (*adc* <sup>+</sup>σ*A*) sup ω∈*R*0<sup>+</sup> (*i*ω*In* <sup>−</sup> *Acc*)−1<sup>2</sup> ⎞ ⎟⎟⎟⎟⎠*. (iii) If <sup>A</sup>*(*t*) *is piecewise continuous and bounded, then Property (ii) holds by replacing <sup>A</sup>*<sup>2</sup> <sup>≤</sup> σ*A*<sup>ε</sup> *by* sup 0≤*t*<∞ *<sup>A</sup>*(*t*)<sup>2</sup> <sup>≤</sup> σ*A*ε*.*

**Remark 3.** *Some quantified results are given in Lemmas B.2 and B.3 to modify* ε2*(and hence* ε <sup>2</sup>*) in Theorem 2 by considering the second power of* ε *in the calculations of the disturbed parameterization guaranteeing the closed-loop stability in the decentralized case.*

**Remark 4.** *If the corresponding parametrical error matrices of Equation (13) have some zero o*ff*-diagonal entries (or o*ff*-diagonal block matrices in the more general case that the system is described by coupled subsystems), then we have at least a partial closed-loop stabilization under decentralized control or, eventually, cut coupled dynamic links to the light of the various Cases 1–5 described after Remark 2 such that closed-loop stability is preserved.*

**Remark 5.** *Theorem 2 also applies to the case of state-feedback control by replacing the output matrices Cc* , *Cd* <sup>→</sup> *In and fixing* σ*<sup>C</sup>* <sup>=</sup> <sup>0</sup>*.*

**Remark 6.** *Theorem 2 also applies directly to the cases where Equations (1)–(3) are a given nominal asymptotically stable closed-loop system configuration and Equations (4)–(6) are a perturbed one whose closed-loop asymptotic stability maintenance related to its nominal counterpart is a suited objective and which is not necessarily of partial of compete decentralized type.*

#### **4. Simulation Examples**

This Section contains some numerical simulation examples to illustrate the theoretical results introduced in Section 3.

**Example 1.** *Consider the interconnected linear system with less inputs than outputs given by, [19]:*

$$\begin{aligned} \dot{\mathbf{x}}\_1(t) &= A\_{11}\mathbf{x}\_1(t) + A\_{12}\mathbf{x}\_2(t) + B\_1\boldsymbol{u}\_1(t) \\ \dot{\mathbf{x}}\_2(t) &= A\_{21}\mathbf{x}\_1(t) + A\_{22}\mathbf{x}\_2(t) + B\_2\boldsymbol{u}\_2(t) \\ y\_1(t) &= C\_1\mathbf{x}\_1(t) \\ y\_2(t) &= C\_2\mathbf{x}\_2(t) \end{aligned}$$

*with x*1(*t*) *<sup>T</sup>* = *x*11(*t*) *x*12(*t*) *, x*2(*t*) *<sup>T</sup>* = *x*21(*t*) *x*22(*t*) *and matrices defined by:*

$$A\_{11} = \begin{bmatrix} 1 & 0 \\ 0 & 2 \end{bmatrix}, A\_{12} = A\_{21} = \begin{bmatrix} -1 & 0 \\ 0 & -1 \end{bmatrix}, A\_{22} = \begin{bmatrix} -1 & 0 \\ 0 & -2 \end{bmatrix}, B\_1 = \begin{bmatrix} 1 \\ 2 \end{bmatrix}, B\_2 = \begin{bmatrix} 2 \\ 1 \end{bmatrix}$$

$$\mathbf{C}\_1 = \mathbf{C}\_2 = I\_2$$

This system can be cast into the form of Equations (1) and (2) by composing the matrices:

$$A\_{\mathfrak{c}} = \begin{bmatrix} 1 & 0 & -1 & 0 \\ 0 & 2 & 0 & -1 \\ -1 & 0 & -1 & 0 \\ 0 & -1 & 0 & -2 \end{bmatrix}, B\_{\mathfrak{c}} = \begin{bmatrix} 1 & 0 \\ 2 & 0 \\ 0 & 2 \\ 0 & 1 \end{bmatrix}, C\_{\mathfrak{c}} = I\_{4\mathfrak{c}}, D\_{\mathfrak{c}} = 0$$

Note that matrix *Ac* is unstable with eigenvalues given by {2.36, .1.41,−1.41,−2.23}. A static feedback output controller of the form of Equation (3) can be designed for this system, which leads to the following gain:

$$K\_c = \begin{bmatrix} -9.3 & 8.2 & -0.015 & 0.02 \\ 0.01 & -0.01 & 0.375 & 0.25 \end{bmatrix}$$

that places the closed-loop poles at {−0.84,−1.1,−2.8,−3.3} and thus stabilizes the closed-loop system. The static feedback gain *Kc* corresponds to a centralized controller as it can be readily observed. The question that arises now is whether a decentralized controller defined by:

$$K\_{\rm d} = \begin{bmatrix} -9.3 & 8.2 & 0 & 0\\ 0 & 0 & 0.375 & 0.25 \end{bmatrix}$$

is enough to stabilize the system or not. Note that *Kd* is a block-diagonal matrix with zero off-block-diagonal entries. Theorem 1 enables us to guarantee the asymptotic stability of the above system when the decentralized controller *Kd* is used. Therefore, *Ad* = *Ac*, *Bd* = *Bc* and *Cd* = *Cc* while the feedback gain *Kd* is restricted to the proposed particular structure. In this way, consider now *Acc* <sup>=</sup> *Ac* <sup>+</sup> *BcKcCc*, *Adc* <sup>=</sup> *Ac* <sup>+</sup> *BcKdCc* and *<sup>A</sup>dc* <sup>=</sup> *Bc*(*Kc* <sup>−</sup> *Kd*)*Cc*. Condition 2 of Theorem 1 *(i)* yields:

$$0.055 = \|\tilde{A}\_{\text{dc}}\|\_{2} < \frac{1}{\sup\_{\omega \in \mathbb{R}\_{0+}} \|\left(\|\omega I\_{4} - A\_{\text{cc}}\right)^{-1}\|\_{2}} = 0.07$$

Consequently, we can conclude from Theorem 1 that the closed-loop system controlled by the decentralized static output gain *Kd* is asymptotically stable. Thus, we have been able to easily analyze the stability of the decentralized case from the stability property of the centralized one. Figure 1 shows the trajectories of the closed-loop system when the gain *Kd* is deployed with initial conditions given by *x*1(*t*) *<sup>T</sup>* = <sup>−</sup><sup>3</sup> <sup>−</sup><sup>4</sup> , *x*2(*t*) *<sup>T</sup>* = 5 6 . It can be observed in Figure 1 that all the states converge to zero as predicted by Theorem 1.

**Figure 1.** States evolution when the decentralized controller *Kd* is employed.

**Example 2.** *Consider the linear system with the same number of inputs and outputs composed of two identical pendulums THAT are coupled by a spring and subject to two distinct inputs, as displayed in Figure 2, [19].*

**Figure 2.** Two inverted pendulums coupled by a spring.

The mathematical model of such interconnected system is given by:

$$\begin{aligned} \dot{x}\_1(t) &= A\_{11}x\_1(t) + A\_{12}x\_2(t) + B\_1u\_1(t) \\ \dot{x}\_2(t) &= A\_{21}x\_1(t) + A\_{22}x\_2(t) + B\_2u\_2(t) \\ y\_1(t) &= C\_1x\_1(t) \\ y\_2(t) &= C\_2x\_2(t) \end{aligned}$$

with *x*1(*t*) *<sup>T</sup>* = θ1 . θ1 , *x*2(*t*) *<sup>T</sup>* = θ2 . θ2 and matrices defined by:

$$A\_{11} = A\_{22} = \begin{bmatrix} 0 & 1 \\ \frac{\mathcal{S}}{\mathcal{I}} - \frac{kx^2}{m^2} & -\mu \end{bmatrix},\\ A\_{12} = A\_{21} = \begin{bmatrix} 0 & 0 \\ \frac{kx^2}{m^2} & 0 \end{bmatrix},\\ B\_1 = B\_2 = \begin{bmatrix} 0 \\ \frac{1}{m^2} \end{bmatrix}$$

$$\mathbf{C}\_1 = \mathbf{C}\_2 = \begin{bmatrix} 1 & 0 \end{bmatrix}$$

where *g* represents the gravity, μ accounts for the friction, *m* = *m*<sup>1</sup> = *m*<sup>2</sup> are the masses of both pendulums, *k* is the spring constant and the meanings of the geometrical parameters are shown in Figure 2. This linear model corresponds to the linearization of the pendulum nonlinear equations around the up-right position equilibrium point. The following values were used in simulation, [19]:

$$\frac{\mathcal{S}}{l} = 1, \frac{1}{ml^2} = 1, \,\mu = 1, \,\frac{k}{m} = 2, \,\frac{a}{l} = 0.5$$

This system can be cast into the form of Equations (1) and (2) as:

$$A\_{\mathfrak{c}} = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0.5 & -1 & 0.5 & 0 \\ 0 & 0 & 0 & 1 \\ 0.5 & 0 & 0.5 & -1 \end{bmatrix}, B\_{\mathfrak{c}} = \begin{bmatrix} 0 & 0 \\ 1 & 0 \\ 0 & 0 \\ 0 & 1 \end{bmatrix}, \mathbf{C}\_{\mathfrak{c}} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix}, D\_{\mathfrak{c}} = 0$$

A static output feedback controller can be designed for this system to achieve its asymptotic stability. In this way, the feedback gain

$$K\_{\odot} = \begin{bmatrix} 2.92 & 0.65 \\ 0.64 & 2.80 \end{bmatrix}$$

places the closed-loop poles at {−0.5 ± 1.5*i*,−0.5 ± 1.4*i*} with negative real parts. Now, we implement a decentralized controller with feedback gain given by:

$$\mathbf{K}\_d = \begin{bmatrix} 2.92 & 0\\ 0 & 2.80 \end{bmatrix}$$

Theorem 2 is now used to analyze the stability of the closed-loop system when this controller is employed. This case is of practical importance and corresponds to the situation when the local controller has only available for control purposes the information regarding the local output, and not the output of the complete system. Thus, the centralized and decentralized systems are the same and only the static feedback gain changes. Theorem 2 conditions are applied with σ*<sup>A</sup>* <sup>=</sup> σ*<sup>B</sup>* <sup>=</sup> σ*<sup>C</sup>* <sup>=</sup> σ*<sup>D</sup>* <sup>=</sup> 0, σ*<sup>K</sup>* <sup>=</sup> *Kc* <sup>−</sup> *Kd*<sup>2</sup> <sup>=</sup> 0.65 while the stability condition for this special case (see Appendix B) reads:

$$0.65 = 0.65 \times 1 \times 1 = \|K\_c - K\_d\|\_2 \|B\_c\|\_2 \|C\_c\|\_2 < 1$$

Accordingly, the closed-loop system attained with the decentralized controller is asymptotically stable and all the outputs will converge to zero asymptotically. Figure 3 displays the evolution of both angles from initial conditions *x*1(*t*) *<sup>T</sup>* = 0.5 <sup>−</sup>0.5 , *x*2(*t*) *<sup>T</sup>* = 0.15 0.5 , where it can be observed that both pendulums are stabilized in the up-right position with the decentralized controller.

**Figure 3.** Evolution of the angles of both pendulums.

**Example 3.** *Consider the linear interconnected system given by:*

$$\begin{aligned} \dot{\mathbf{x}}\_{\mathfrak{c}}(t) &= A\_{\mathfrak{c}} \mathbf{x}\_{\mathfrak{c}}(t) + B\_{\mathfrak{c}} \boldsymbol{\mu}\_{\mathfrak{c}}(t) \\ y\_{\mathfrak{c}}(t) &= \mathsf{C}\_{\mathfrak{c}} \mathbf{x}\_{\mathfrak{c}}(t) + D\_{\mathfrak{c}} \boldsymbol{\mu}\_{\mathfrak{c}}(t) \end{aligned}$$

*with matrices defined by:*

$$A\_{\mathfrak{c}} = \begin{bmatrix} -1 & 0.4 & 0.3 \\ 0.2 & -2 & 0.1 \\ -0.1 & 0.2 & -3 \end{bmatrix}, B\_{\mathfrak{c}} = \begin{bmatrix} 1 & 0.1 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 2 \end{bmatrix}, \mathbf{C}\_{\mathfrak{c}} = \begin{bmatrix} 1 & 0.1 & 0 \\ 0 & 1 & 0.1 \\ 0 & 0 & 1 \end{bmatrix}, D\_{\mathfrak{c}} = 0.1l\_{\mathfrak{c}}$$

This system is controlled by the static output feedback gain given by:

$$K\_c = \begin{bmatrix} 0.19 & 0.05 & 0.04 \\ 0.05 & -0.01 & 0.03 \\ -0.02 & 0.02 & -0.46 \end{bmatrix}$$

which places the closed-loop poles at {−1.17,−2,−2.1}. The decentralized system is now given by:

$$A\_d = A\_{c'} B\_d = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 2 \end{bmatrix} \\ \text{'C}\_d = I\_{3\prime} D\_d = 0$$

The decentralized system corresponds to the case when some transmission links have been suppressed from the original open-loop coupled dynamics, as considered in Remark 2. The following decentralized gain iz employed to stabilize the decentralized system in Equations (4)–(6) parameterized by the above matrices:

$$K\_d = \begin{bmatrix} 0.19 & 0 & 0 \\ 0 & -0.01 & 0 \\ 0 & 0 & -0.46 \end{bmatrix}$$

Theorem 2 is used to analyze the stability of the decentralized closed-loop system. To this end, we calculate:

$$\|A\_{\mathfrak{c}} - A\_{d}\|\_{2} = \overline{\sigma}\_{A} = 0$$

$$0.1 = \|B\_{\mathfrak{c}} - B\_{d}\|\_{2} \le \overline{\sigma}\_{B}\varepsilon = 2.7 \times 0.07 = 0.19$$

$$0.1 = \|\mathbf{C}\_{\mathfrak{c}} - \mathbf{C}\_{d}\|\_{2} \le \overline{\sigma}\_{\mathfrak{C}}\varepsilon = 2.7 \times 0.07 = 0.19$$

$$0.1 = \|D\_{\varepsilon} - D\_{d}\|\_{2} \le \overline{\sigma}\_{D}\varepsilon = 2.7 \times 0.07 = 0.19$$

$$0.07 = \|K\_{\varepsilon} - K\_{d}\|\_{2} \le \overline{\sigma}\_{K}\varepsilon = 2.7 \times 0.07 = 0.19$$

With these values, we can compute ε = 0.31, ε<sup>1</sup> = 0.24, ε<sup>2</sup> = 0.071 so that 0.07 = ε < ε<sup>∗</sup> = *min*(1, ε, ε1, ε2) = 0.071. Since *KcDc*<sup>2</sup> = 0.05 < 1, we are in conditions of applying Theorem 2 *(i)* and we can conclude that the decentralized closed-loop system is asymptotically stable. In this way, the presented results allow establishing the stability of the decentralized system by a simple method based on the stability and design of the centralized system. Figure 4 shows the state variables evolution from the initial state *x*(*t*) *<sup>T</sup>* = <sup>5</sup> <sup>−</sup>5 1 . As shown in Figure 4, the state variables converge to zero asymptotically, as concluded from Theorem 2.

**Figure 4.** Evolution of the state space variables when the decentralized controller *Kd* is used.

#### **5. Concluding Remarks**

This paper is devoted to formulating sufficiency-type linear-output feedback decentralized closed-loop stabilization conditions if the continuous-time linear dynamic system can be stabilized under linear output-feedback centralized stabilization. The developed stability tests are conceptually simple to evaluate and they rely on the quantification in terms of worst-case norms of interconnection and open-loop system dynamics matrices and the corresponding control gains in the decentralized case compared to the centralized counterpart. The tolerances of the various parametrical matrix errors have been quantified by considering the first or second powers of a small parameter. Such a parameter is a design factor to characterize in the worst-case for the allowed tolerances to the perturbed parameterization norms. Simulated examples are discussed to illustrate the obtained results. The decentralized control design versus its decentralized control counterpart, under eventual output linear feedback, has been studied from the point of view of the amount of information that can be lost or omitted in terms of the total or partial knowledge of the coupled dynamics between subsystems necessary in the decentralized case to keep the closed-loop stability. A foreseen related future work relies on the application of the method to some applied control problems such as consensus protocols under decentralized control and continuous-discrete hybrid controller designs.

**Author Contributions:** Author M.D.l.S. contributed to the manuscript in the aspects of conceptualization; methodology; formal analysis; investigation; writing original draft preparation; writing review and editing; supervision and fund acquisition. Author A.I. contributed to the manuscript in the various aspects of preparation of the numerical experiments; software implementation; visualization and original draft preparation.

**Funding:** This research was co-funded by the Spanish Government and the European Fund of Regional Development FEDER, grant numbers DPI2015-64766-R and DPI2016-77271-R (MINECO/FEDER, UE) and the University of the Basque Country UPV /EHU, grant PGC 17/33.

**Acknowledgments:** The authors are grateful to the Spanish Government for grants DPI2015-64766-R and DPI2016-77271-R (MINECO/FEDER, UE), and to UPV /EHU for grant PGC 17/33. They also thank the referees for their useful comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A** *Auxiliary Stability Results on Perturbed Matrices under Constant and Time-Varying Perturbations of Stability Matrices*

**Theorem A1.** *Assume that <sup>A</sup>*<sup>0</sup> <sup>∈</sup> *<sup>R</sup>n*×*<sup>n</sup> is a stability matrix with stability abscissa* (−ρ*c*) <sup>&</sup>lt; 0*. Then, <sup>A</sup>* <sup>=</sup> *<sup>A</sup>*<sup>0</sup> <sup>+</sup> *A is a stability matrix if any of the subsequent su* ffi*ciency-type conditions holds:*


**Proof.** Note that

$$\begin{split} \det(sI\_{\mathbb{H}} - A) &= \det(sI\_{\mathbb{H}} - A\_0 - \overline{A}) = \det\Big( (sI\_{\mathbb{H}} - A\_0) \left( I\_{\mathbb{H}} - (sI\_{\mathbb{H}} - A\_0)^{-1} \overline{A} \right) \Big) \\ &= \det(sI\_{\mathbb{H}} - A\_0) \det\Big( I\_{\mathbb{H}} - (sI\_{\mathbb{H}} - A\_0)^{-1} \overline{A} \Big); \forall s \in \mathbb{C} \end{split} \tag{A1}$$

and then det(*sIn* <sup>−</sup> *<sup>A</sup>*) <sup>=</sup> det(*sIn* <sup>−</sup> *<sup>A</sup>*0) det- *In* <sup>−</sup> (*sIn* <sup>−</sup> *<sup>A</sup>*0)−1*A* - 0 Then, for all *s sp*(*A*), and also for all *<sup>s</sup>* <sup>∈</sup> *sp*(*A*0) if the *<sup>H</sup>*∞-norm of (*sIn* <sup>−</sup> *<sup>A</sup>*0)−1*A*, which exists since *<sup>A</sup>*<sup>0</sup> is a stability matrix, satisfies (*sIn* <sup>−</sup> *<sup>A</sup>*0)−1*A*∞ <sup>&</sup>lt; 1, which is guaranteed if *A*<sup>2</sup> <sup>&</sup>lt; 1/ sup ω∈*R*0<sup>+</sup> (*i*ω*In* <sup>−</sup> *<sup>A</sup>*0)−12. Then, *<sup>A</sup>* is a stability matrix if Conditions 1 or 2 holds. On the other hand, if *A*<sup>0</sup> and *A* are negative (implying that *A* ≺ −*A*0), or if they are both Metzler-stable (implying for all off-diagonal entries that *<sup>A</sup>ij* ≥ −*A*0*ij* ; <sup>∀</sup>*i*, *<sup>j</sup>*( *i*) ∈ *n*), then their dominant abscissa (perhaps multiple) eigenvalue is real and negative since *A*<sup>0</sup> being a stability matrix is claimed to guarantee that *A* is stable. Since *A*<sup>0</sup> is a stability matrix, it is non-singular with eigenvalues with negative real parts. Then, by the continuity of the eigenvalues with respect to the matrix entries, *A* = *A*<sup>0</sup> - *In* + *A*−<sup>1</sup> <sup>0</sup> *<sup>A</sup>* is a stability matrix if ρ - *A*−<sup>1</sup> <sup>0</sup> *<sup>A</sup>* ≤ *A*−<sup>1</sup> <sup>0</sup> *<sup>A</sup>* <sup>2</sup> ≤ *A*−<sup>1</sup> <sup>0</sup> 2 *A*<sup>2</sup> <sup>&</sup>lt; <sup>1</sup> leading to the sufficiency of Conditions 3–5 for the stability of *A* if is *A*<sup>0</sup> stable. The last, sufficient condition comes directly by upper-bounding Condition 4 by norm product and it is equivalent to *A*<sup>2</sup> <sup>=</sup> <sup>λ</sup>1/2 max- *<sup>A</sup>TA* <sup>&</sup>lt; 1/*A*−<sup>1</sup> <sup>0</sup> <sup>2</sup> <sup>=</sup> 1/λ1/2 max- *A*−<sup>1</sup> <sup>0</sup> *<sup>A</sup>*−*<sup>T</sup>* 0 = λ1/2 min- *AT* <sup>0</sup> *A*<sup>0</sup> . -

**Theorem A2.** *Assume that <sup>A</sup>*<sup>0</sup> <sup>∈</sup> *<sup>R</sup>n*×*<sup>n</sup> is a stability matrix and that <sup>A</sup>*(*t*) <sup>=</sup> *<sup>A</sup>*<sup>0</sup> <sup>+</sup>*A*(*t*)*, where <sup>A</sup>* : *<sup>R</sup>*0<sup>+</sup> <sup>→</sup> *<sup>R</sup>n*×*<sup>n</sup> is piecewise-continuous and bounded. Then, <sup>A</sup>* : *<sup>R</sup>*0<sup>+</sup> <sup>→</sup> *<sup>R</sup>n*×*<sup>n</sup> is stable if <sup>k</sup>*<sup>0</sup> <sup>ρ</sup><sup>0</sup> sup 0≤τ≤*t A*(τ) <sup>&</sup>lt; <sup>1</sup>*;* <sup>∀</sup>*<sup>t</sup>* <sup>∈</sup> *<sup>R</sup>*0+*, where* (−ρ0) < 0 *is the stability abscissa of A*<sup>0</sup> *and k*<sup>0</sup> = *k*0(*A*0) ≥ 1 *is a real constant satisfying eA*0*<sup>t</sup>* ≤ *<sup>k</sup>*0*e*−ρ0*<sup>t</sup> ;*∀*t* ∈ *R*0+*.*

**Proof.** Consider the linear time-varying system:

$$\dot{\mathbf{x}}(t) = \left(A\_0 + \tilde{A}(t)\right)\mathbf{x}(t), \; \mathbf{x}(0) = \mathbf{x}\_0; \; \forall t \in \mathbb{R}\_{0+} \tag{A2}$$

where *<sup>A</sup>* : *<sup>R</sup>*0<sup>+</sup> <sup>→</sup> *<sup>R</sup>n*×*<sup>n</sup>* is piecewise-continuous and bounded [1]. Such a system is globally asymptotically stable if and only if *<sup>A</sup>* : *<sup>R</sup>*0<sup>+</sup> <sup>→</sup> *<sup>R</sup>n*×*<sup>n</sup>* is a stability matrix. The state-trajectory solution of Equation (A2) satisfies:

$$\begin{split} \|\mathbf{x}(t)\| &\leq \|e^{A\_0 t}\| \|\mathbf{x}\_0\| + \int\_0^t \|e^{A\_0(t-\tau)}\| \, \|\widetilde{A}(\tau)\| \|\mathbf{x}(\tau)\| \, d\tau \leq k\_0 e^{-\rho\_0 t} \|\mathbf{x}\_0\| + \frac{k\_0}{\rho\_0} \sup\_{0 \leq \tau \leq t} \|\mathbf{x}(\tau)\| \\ &= K\_0 e^{-\rho\_0 t} \|\mathbf{x}\_0\| + \frac{k\_0}{\rho\_0} \sup\_{0 \leq \tau \leq t} \left(\|\widetilde{A}(\tau)\|\right) \|\mathbf{x}(t)\| \, \forall t \in \mathcal{R}\_{0+} \end{split} \tag{A3}$$

Let *t* = *t*(*t*) be defined for each *t* ∈ *R*0<sup>+</sup> as *t* = *z* = max 0≤τ≤*t* τ : *x*(*z*)≥*x*(*t*) . Then,

*x*(*t*) ≤ sup 0≤τ≤*t x*(τ) = sup 0≤τ≤*t x*(τ) = *x*(*t*) ≤ *k*<sup>0</sup> −ρ0*t x*0 + *k*0 ρ0 sup 0≤τ≤*t A*(τ) *x*(*t*); <sup>∀</sup>*<sup>t</sup>* <sup>∈</sup> *<sup>R</sup>*0<sup>+</sup> (A4)

Since 1 > *<sup>k</sup>*<sup>0</sup> <sup>ρ</sup><sup>0</sup> sup 0≤τ≤*t A*(τ); <sup>∀</sup>*<sup>t</sup>* <sup>∈</sup> *<sup>R</sup>*0+, one gets:

$$\begin{split} \|\mathbf{x}(t)\| \le & \|\mathbf{x}(t\tau)\| = \sup\_{0 \le \tau \le t\tau} \|\mathbf{x}(\tau)\| \le \left(1 - \frac{k\_0}{\rho\_0} \sup\_{0 \le \tau < \infty} \|\widetilde{A}(\tau)\|\right)^{-1} k\_0 e^{-\rho\_0 t} \|\mathbf{x}\_0\| \\ \le & M = k\_0 \Big(1 - \frac{k\_0}{\rho\_0} \sup\_{0 \le \tau < \infty} \|\widetilde{A}(\tau)\|\Big)^{-1} \|\mathbf{x}\_0\| \; \forall t \in \mathbb{R}\_{0+} \end{split} \tag{A5}$$

Therefore, *x*(*t*) is bounded for any *t* ∈ *R*0<sup>+</sup> if *x*<sup>0</sup> is finite. Now, assume the following cases.

*Case a*: For any *Ts* ∈ *R*+, the sequence ⎧ ⎪⎪⎨ ⎪⎪⎩ sup *nTs*≤*t*<(*n*+*mn*)*Ts x*(*t*) ⎫ ⎪⎪⎬ ⎪⎪⎭ ∞ *n*=*k* is strictly decreasing for some finite positive integer *k* = *k*(*Ts*) and some positive sequence {*mn*} of bounded integer numbers which satisfies *mn*+<sup>1</sup> > *mn* − 1 for *n* ≥ 0. As a result, ⎧ ⎪⎪⎨ ⎪⎪⎩ sup *nTs*≤*t*<(*n*+*mn*)*Ts x*(*t*) ⎫ ⎪⎪⎬ ⎪⎪⎭ ∞ *n*=*k* → 0 as *n* → ∞ for any given *Ts* > 0. Then, one gets from Equation (A3) that *x*(*t*) → 0, as *t* → ∞ since *x*(*nTs* + *t*)≤ *k*0 ⎛ ⎜⎜⎜⎜⎝ 1 + <sup>1</sup> <sup>ρ</sup><sup>0</sup> sup *nTs*≤τ<*t x*(τ) ⎞ ⎟⎟⎟⎟⎠ ; ∀*t* ∈ (*nTs* , (*n* + 1)*Ts*). The result is proved for this case. ∞

*Case b*: For some *Ts* ∈ *R*+, a sequence ⎧ ⎪⎪⎨ ⎪⎪⎩ sup *nTs*≤*t*<(*n*+*mn*)*Ts x*(*t*) ⎫ ⎪⎪⎬ ⎪⎪⎭ *n*=0 can be built, with {*mn*} ∞ *<sup>n</sup>*=<sup>0</sup> → ∞ as *m* → ∞ satisfying sup *x*(*t*) = sup *x*(*t*) ≥ sup *x*(*t*) ; *n* ≥ 0

(*n*+*mn*)*Ts*≤*t*<(*n*+*mn*)*Ts nTs*≤*t*<(*n*+*mn*+*mn*+*mn* )*Ts nTs*≤*t*<(*n*+*mn*)*Ts* (note that the above inequality cannot be strict as *n* → ∞ since it has already been proven that *x*(*t*) < +∞; ∀*t* ∈ *R*0+). However, then one gets from Equation (A3) for some *tn*+<sup>1</sup> ∈ [*n* + *mn* , *n* + *mn* + *mn*+*mn* ) since 1 > *<sup>k</sup>*<sup>0</sup> <sup>ρ</sup><sup>0</sup> sup *A*(τ); <sup>∀</sup>*<sup>t</sup>* <sup>∈</sup> *<sup>R</sup>*0+:

0≤τ≤*t*

$$\begin{array}{lcl}\displaystyle\sup\_{t\in T\_{s}\leq t\leq (n+m\_{n})T\_{s}}\|\mathbf{x}(t)\|\leq&\sup\_{nT\_{s}\leq t\leq (n+m\_{n}+m\_{n+m\_{n}})T\_{s}}\|\mathbf{x}(t)\|\leq&K\_{0}e^{-p\rho\left(t\_{n+1}-nT\_{s}\right)}\|\mathbf{x}(nT\_{s})\|\|\\\\\displaystyle+\frac{K\_{0}}{\rho\_{0}}\sup\_{nT\_{s}\leq t\leq (n+m\_{n}+m\_{n+m\_{n}})T\_{s}}\left(\|\widetilde{A}(\tau)\|\right)\sup\_{nT\_{s}\leq t\leq (n+m\_{n}+m\_{n+m\_{n}})T}\|\mathbf{x}(t)\|\\\\ \leq&K\_{0}e^{-p\rho\left(t\_{n+1}-nT\_{s}\right)}\|\mathbf{x}(nT\_{s})\|\|+&\sup\_{nT\_{s}\leq t\leq (n+m\_{n}+m\_{n+m\_{n}})T}\|\mathbf{x}(t)\|\end{array}\tag{A6}$$

If *n* → ∞, *mn* → ∞, then (*tn*+<sup>1</sup> − *mn*+*mn* ) → ∞, thus the following contradiction arises:

$$0 = \lim\_{n \to \infty} \sup \left( \sup\_{nT\_s \le t < (n + m\_n + m\_{n+m\_n})T\_s} \|\mathbf{x}(t)\| - \sup\_{nT\_s \le t < (n + m\_n + m\_{n+m\_n})T\_s} \|\mathbf{x}(t)\| \right) < 0.$$

Thus, Case b is not possible and the whole result follows from Case a. -

Note that the stability abscissa of *A*0, that is, (−ρ0) < 0 is not smaller than the dominant eigenvalue abscissa.

#### **Appendix B** *Calculations for solving Problem 1*

Assume that the matrix (*Im* − *KcDc*) is non-singular with *KcDc*<sup>2</sup> < 1 and *A* = *Ac* = *Ad*. Then, one gets from Equations (10) and (13) that:

$$\begin{split} & \quad A\_{d} + B\_{d} (I\_{m} - K\_{d} D\_{d})^{-1} K\_{d} \mathbb{C}\_{d} \\ &= A + \left( B\_{\text{c}} - \widehat{B} \right) \left[ I\_{m} - K\_{\text{c}} D\_{\text{c}} - \widetilde{\Delta}\_{0} \right]^{-1} \left( K\_{\text{c}} D\_{\text{c}} - \widetilde{\Delta}\_{1} \right) \\ &= A + \left( B\_{\text{c}} - \widehat{B} \right) \left[ \left( I\_{m} - K\_{\text{c}} D\_{\text{c}} \right) \left( I\_{m} - \left( I\_{m} - K\_{\text{c}} D\_{\text{c}} \right)^{-1} \widetilde{\Delta}\_{0} \right) \right]^{-1} \left( K\_{\text{c}} C\_{\text{c}} - \widetilde{\Delta}\_{1} \right) \\ &= A + \left( B\_{\text{c}} - \widehat{B} \right) \left( I\_{m} - \left( I\_{m} - K\_{\text{c}} D\_{\text{c}} \right)^{-1} \widetilde{\Delta}\_{0} \right)^{-1} \left( I\_{m} - K\_{\text{c}} D\_{\text{c}} \right)^{-1} \left( K\_{\text{c}} C\_{\text{c}} - \widetilde{\Delta}\_{1} \right) \\ &= A + \left( B\_{\text{c}} - \widehat{B} \right) \left( I\_{m} + \widetilde{\Delta}\_{2} \right) \left( I\_{m} - K\_{\text{c}} D\_{\text{c}} \right)^{-1} \left( K\_{\text{c}} C\_{\text{c}} - \widetilde{\Delta}\_{1} \right) \end{split} \tag{A7}$$

provided that *<sup>D</sup>* and *<sup>K</sup>* are such that (*Im* <sup>−</sup> *KdDd*)−<sup>1</sup> <sup>=</sup> - *Im* <sup>−</sup> *KcDc* <sup>−</sup> Δ<sup>0</sup> <sup>−</sup><sup>1</sup> exists (note that this always holds if *<sup>D</sup>* <sup>=</sup> <sup>0</sup>*p*×*<sup>m</sup>* and *<sup>K</sup>* <sup>=</sup> <sup>0</sup>*m*×*<sup>p</sup>* from Assumption 1), where:

$$
\widetilde{\Delta}\_0 = K\_d D\_d - K\_\mathfrak{c} D\_\mathfrak{c} = \widetilde{\mathcal{K}} (\widetilde{D} - D\_\mathfrak{c}) - K\_\mathfrak{c} \widetilde{D} = \left( \widetilde{\mathcal{K}} - K\_\mathfrak{c} \right) \widetilde{D} - \widetilde{\mathcal{K}} D\_\mathfrak{c} \tag{A8}
$$

$$
\widetilde{\Lambda}\_1 = \widetilde{K} \{ \mathbf{C}\_{\varepsilon} - \widetilde{\mathbf{C}} \} + K\_{\varepsilon} \widetilde{\mathbf{C}} = \left( K\_{\varepsilon} - \widetilde{K} \right) \widetilde{\mathbf{C}} + \widetilde{K} \mathbf{C}\_{\varepsilon} \tag{A9}
$$

$$
\widetilde{\Delta}\_2 = \left[ I\_{\rm ll} - \left( I\_{\rm ll} - K\_{\rm c} D\_{\rm c} \right)^{-1} \widetilde{\Delta}\_0 \right]^{-1} - I\_{\rm m} \tag{A10}
$$

and note that

$$\|\Lambda\_0\| \le \varepsilon \left[ \|K\_{\varepsilon}\| \overline{\sigma}\_D + \varepsilon \overline{\sigma}\_D \overline{\sigma}\_K + \|D\_{\varepsilon}\| \overline{\sigma}\_K \right] \tag{A11}$$

$$\|\|\Delta\_1\|\| \le \varepsilon \left[ \|\|K\_{\mathbf{c}}\|\|\overline{\sigma}\_{\mathbf{C}} + \varepsilon \overline{\sigma}\_{\mathbf{C}} \overline{\sigma}\_K + \|\|\mathbf{C}\_{\mathbf{c}}\|\|\overline{\sigma}\_K\right] \tag{A12}$$

and one gets from Banach´s Perturbation Lemma [7] that

$$\|\|\tilde{\Delta}\_2\|\| \le 1 + \frac{1}{1 - \|(I\_m - K\_c D\_c)^{-1}\widetilde{\Delta}\_0\|} \le 1 + \frac{1}{1 - \|(I\_m - K\_c D\_c)^{-1}\| \|\|\widetilde{\Delta}\_0\|\|}\tag{A13}$$

provided that Δ0 <sup>&</sup>lt; δ<sup>0</sup> <sup>=</sup> 1/(*Im* <sup>−</sup> *KcDc*)−1. Equivalently, if

$$q(\varepsilon) = \widetilde{\sigma}\_D \widetilde{\sigma}\_k \varepsilon^2 + \left(\widetilde{\sigma}\_D \|\mathcal{K}\_\mathbf{c}\| + \widetilde{\sigma}\_K \|D\_\mathbf{c}\|\right) \varepsilon - 1/\|(I\_m - \mathcal{K}\_\mathbf{c} D\_\mathbf{c})^{-1}\|\tag{A14}$$

since *q*(ε) is a convex parabola with zeros ε<sup>1</sup> < 0 and ε = ε<sup>2</sup> > 0, Equation (A14) holds, guaranteeing that Δ0 <sup>&</sup>lt; δ0, if <sup>ε</sup> <sup>∈</sup> [0 , <sup>ε</sup>), where:

$$\overline{\varepsilon} = \frac{\sqrt{\left(\widetilde{\sigma}\_{D} \|\mathbf{K}\_{\rm c}\| + \widetilde{\sigma}\_{\rm K} \|\mathbf{D}\_{\rm c}\|\right)^{2} + 4\widetilde{\sigma}\_{\rm D}\widetilde{\sigma}\_{\rm K} / \|(\mathbf{I}\_{\rm m} - \mathbf{K}\_{\rm c}\mathbf{D}\_{\rm c})^{-1}\|} - \left(\widetilde{\sigma}\_{\rm D} \|\mathbf{K}\_{\rm c}\| + \widetilde{\sigma}\_{\rm K} \|\mathbf{D}\_{\rm c}\|\right)}\tag{A15}$$

#### Before continuing with the calculations, we give the following auxiliary result:

**Lemma B1.** *If* (*Im* <sup>−</sup> *KcDc*) *is non-singular with KcDc*<sup>2</sup> <sup>&</sup>lt; <sup>1</sup> *and* Δ0 <sup>&</sup>lt; δ<sup>0</sup> <sup>=</sup> 1/(*Im* <sup>−</sup> *KcDc*)−1*, equivalently if* <sup>ε</sup> <sup>∈</sup> [0 , <sup>ε</sup>)*, with* <sup>ε</sup> *defined in Equation (A15), then* Δ2 ≤ *<sup>C</sup>*Δ0 <sup>&</sup>lt; <sup>1</sup> *with a norm-dependent real constant C* <sup>≥</sup> <sup>1</sup> <sup>2</sup>δ<sup>0</sup> *if* <sup>ε</sup> <sup>∈</sup> [0 , <sup>ε</sup>1) *with*

$$\overline{\varepsilon}\_{1} = \frac{\|(I\_{m} - K\_{\varepsilon}D\_{\varepsilon})^{-1}\|}{2[\|K\_{\varepsilon}\|\overline{\sigma}\_{D} + \varepsilon \overline{\sigma}\_{D} \overline{\sigma}\_{K} + \|D\_{\varepsilon}\| \overline{\sigma}\_{K}]} \tag{A16}$$

*Appl. Sci.* **2019**, *9*, 1739

**Proof.** One gets from Equation (A10) and Banach´s Perturbation Lemma [7] that, if Δ2 ≤ *<sup>C</sup>*Δ0 for some *C* ∈ *R*+, then:

$$\frac{1}{\|\mathbf{1} - \mathbf{C}\| \|\widetilde{\boldsymbol{\Delta}}\_{\mathcal{O}}\|} \ge \|\{\widetilde{\boldsymbol{\Delta}}\_{2} + I\_{m}\}^{-1}\| = \|\boldsymbol{I}\_{m} - \left(\boldsymbol{I}\_{m} - \boldsymbol{K}\_{\boldsymbol{c}}\boldsymbol{D}\_{\boldsymbol{c}}\right)^{-1} \widetilde{\boldsymbol{\Delta}}\_{\mathcal{O}}\| \ge 1 - \|(\boldsymbol{I}\_{m} - \boldsymbol{K}\_{\boldsymbol{c}}\boldsymbol{D}\_{\boldsymbol{c}})^{-1}\| \|\boldsymbol{\widetilde{\Delta}}\_{\mathcal{O}}\| $$

provided that *<sup>C</sup>* <sup>&</sup>lt; 1/Δ0. One gets that the above inequality holds if 1/Δ0 <sup>&</sup>gt; *<sup>C</sup>* <sup>≥</sup> <sup>1</sup> <sup>2</sup>δ<sup>0</sup> ≥ (*Im*−*KcDc*)−1 <sup>1</sup>+Δ0 (*Im*−*KcDc*)−1 and, one gets from Equation (A11) that

$$\|\overline{\Delta}\_2\| \le \varepsilon \, \mathbb{C} \left[ \|\mathbb{K}\_{\varepsilon}\| \widetilde{\sigma}\_D + \varepsilon \overline{\sigma}\_D \widetilde{\sigma}\_K + \|D\_{\varepsilon}\| \, \widetilde{\sigma}\_K \right] < 1 \tag{A17}$$

if ε < ε1. -

Now, rewrite the system matrices of closed-loop dynamics of Equations (7) and (10), equivalently Equation (A7), with *A* = *Ac* = *Ad* and its incremental value as follows:

$$A\_{\mathfrak{c}\mathfrak{c}} = A + B\_{\mathfrak{c}} (I\_m - K\_{\mathfrak{c}} D\_{\mathfrak{c}})^{-1} K\_{\mathfrak{c}} \mathbb{C}\_{\mathfrak{c}} \tag{A18}$$

$$A\_{dc} = A + B\_d \left(I\_m - K\_d D\_d\right)^{-1} K\_d \mathbb{C}\_d \tag{A19}$$

*<sup>A</sup>dc* <sup>=</sup> *Acc* <sup>−</sup> *Adc* <sup>=</sup> *Bc*(*Im* <sup>−</sup> *KcDc*)−1*KcCc* <sup>−</sup> *Bd*(*Im* <sup>−</sup> *KdDd*)−<sup>1</sup>*KdCd* = *Bc*(*Im* <sup>−</sup> *KcDc*)−1*KcCc* <sup>−</sup> - *Bc* <sup>−</sup> *<sup>B</sup>* -*Im* <sup>+</sup> Δ<sup>2</sup> (*Im* − *KcDc*) −1 - *KcCc* <sup>−</sup> Δ<sup>1</sup> = *Bc*(*Im* <sup>−</sup> *KcDc*)−1*KcCc* <sup>−</sup> *Bc* - *Im* <sup>+</sup> Δ<sup>2</sup> (*Im* − *KcDc*) −1 - *KcCc* <sup>−</sup> Δ<sup>1</sup> +*B* - *Im* <sup>+</sup> Δ<sup>2</sup> (*Im* − *KcDc*) −1 - *KcCc* <sup>−</sup> Δ<sup>1</sup> <sup>=</sup> *Bc*(*Im* <sup>−</sup> *KcDc*)−1Δ<sup>1</sup> −*Bc* Δ2(*Im* <sup>−</sup> *KcDc*)−1*KcCc* <sup>+</sup> *Bc* Δ2(*Im* <sup>−</sup> *KcDc*)−1Δ<sup>1</sup> <sup>+</sup>*B*(*Im* <sup>−</sup> *KcDc*) −1 *KcCc* <sup>−</sup> *B*(*Im* <sup>−</sup> *KcDc*) −1 Δ1 <sup>+</sup>*B*Δ2(*Im* <sup>−</sup> *KcDc*) −1 *KcCc* <sup>−</sup> *B*Δ2(*Im* <sup>−</sup> *KcDc*) −1 Δ1 (A20)

Now, the following technical result follows directly from Equations (A20), (A11), (A12) and (A17), the norm upper-bounding values of the control, output interconnection and controller matrix errors and Lemma B1:

**Lemma B2.** *The following properties hold for any* ε ∈ [0 , εˆ) *with* εˆ = min(ε , ε1) *calculated from Equations (A15) and (A16):*

*(i)*

*Adc*≤*Bc*(<sup>1</sup> − *KcDc*)−<sup>1</sup> Δ1 - <sup>1</sup> <sup>+</sup> Δ2 <sup>+</sup> *KcCc*Δ2 +ε (1 − *KcDc*) −1 σ*B* - <sup>1</sup> <sup>+</sup> Δ2 -*KcCc* <sup>+</sup> Δ1 <sup>≤</sup> <sup>ε</sup>*Bc*(<sup>1</sup> − *KcDc*)−<sup>1</sup> <sup>×</sup>[(*Kc*σ*<sup>C</sup>* <sup>+</sup> <sup>ε</sup>σ*C*σ*<sup>K</sup>* <sup>+</sup> *Cc*σ*K*)(<sup>1</sup> <sup>+</sup> *KcCc* <sup>+</sup> <sup>ε</sup>*<sup>C</sup>* [*Kc*σ*<sup>D</sup>* <sup>+</sup> <sup>ε</sup>σ*D*σ*<sup>K</sup>* <sup>+</sup> *Dc*σ*K*])*<sup>C</sup>* [*Kc*σ*<sup>D</sup>* <sup>+</sup> *Dc*σ*<sup>K</sup>* <sup>+</sup> <sup>ε</sup>σ*D*σ*K*]] +ε (1 − *KcDc*) −1 σ*B*(<sup>1</sup> <sup>+</sup> <sup>ε</sup>*<sup>C</sup>* [*Kc*σ*<sup>D</sup>* <sup>+</sup> <sup>ε</sup>σ*D*σ*<sup>K</sup>* <sup>+</sup> *Dc*σ*K*])(*KcCc* <sup>+</sup> <sup>ε</sup> [*Kc*σ*<sup>C</sup>* <sup>+</sup> <sup>ε</sup>σ*C*σ*<sup>K</sup>* <sup>+</sup> *Cc*σ*K*] ) (A21)

*(ii) If, furthermore,* ε ≤ <sup>1</sup>*, then* ε*<sup>r</sup>* ≤ ε *for any real r* ≥ <sup>1</sup> *so that one gets from Equation (A21) by taking the upper-bound* ε<sup>2</sup> *for* ε<sup>3</sup> *that*

$$\begin{split} \|\overline{\boldsymbol{A}}\_{\mathrm{id}}\| &\leq \varepsilon \|\boldsymbol{B}\_{\mathrm{c}}\| (1 - \|\boldsymbol{K}\_{\mathrm{c}}\boldsymbol{D}\_{\mathrm{c}}\|)^{-1} \Big( \|\overline{\boldsymbol{A}}\_{10}\| + \|\boldsymbol{K}\_{\mathrm{c}}\boldsymbol{C}\_{\mathrm{c}}\| \|\overline{\boldsymbol{A}}\_{20}\| \Big) + \varepsilon^{2} \|\overline{\boldsymbol{A}}\_{20}\| \|\boldsymbol{B}\_{\mathrm{c}}\| (1 - \|\boldsymbol{K}\_{\mathrm{c}}\boldsymbol{D}\_{\mathrm{c}}\|)^{-1} \\ &+ \varepsilon \left( 1 - \|\boldsymbol{K}\_{\mathrm{c}}\boldsymbol{D}\_{\mathrm{c}}\| \right)^{-1} \overline{\boldsymbol{\sigma}\_{\mathrm{B}}} \|\boldsymbol{K}\_{\mathrm{c}}\boldsymbol{C}\_{\mathrm{c}}\| + \varepsilon^{2} \left( 1 - \|\boldsymbol{K}\_{\mathrm{c}}\boldsymbol{D}\_{\mathrm{c}}\| \right)^{-1} \overline{\boldsymbol{\sigma}\_{\mathrm{B}}} \Big( \|\widetilde{\boldsymbol{A}}\_{10}\| + \|\boldsymbol{K}\_{\mathrm{c}}\boldsymbol{C}\_{\mathrm{c}}\| \|\widetilde{\boldsymbol{A}}\_{20}\| \Big) \\ &+ \varepsilon^{3} \left( 1 - \|\boldsymbol{K}\_{\mathrm{c}}\boldsymbol{D}\_{\mathrm{c}}\| \right)^{-1} \overline{\boldsymbol{\sigma}\_{\mathrm{B}}} \|\widetilde{\boldsymbol{A}}\_{20}\| \|\widetilde{\boldsymbol{A}}\_{10}\| \leq \mu\varepsilon + \upsilon\varepsilon^{2} \leq (\mu + \upsilon)\varepsilon \end{split} \tag{A22}$$

*where* Δ*<sup>i</sup>* <sup>=</sup> <sup>ε</sup>Δ*i*<sup>0</sup> *for i* <sup>=</sup> 1, 2 *and it has been used that* <sup>ε</sup><sup>3</sup> <sup>≤</sup> <sup>ε</sup>2*, with*

$$\begin{split} u &= \|\mathbb{B}\_{\mathbf{c}}\| (1 - \|\mathbb{K}\_{\mathbf{c}}D\_{\mathbf{c}}\|)^{-1} \Big( \|\widetilde{\Lambda}\_{10}\| + \|\mathbb{K}\_{\mathbf{c}}\mathbb{C}\_{\mathbf{c}}\| \|\widetilde{\Lambda}\_{20}\| \Big) + \left(1 - \|\mathbb{K}\_{\mathbf{c}}D\_{\mathbf{c}}\|\right)^{-1} \widetilde{\sigma}\_{\mathbf{B}} \|\mathbb{K}\_{\mathbf{c}}\mathbb{C}\_{\mathbf{c}}\| \\ v &= \left(1 - \|\mathbb{K}\_{\mathbf{c}}D\_{\mathbf{c}}\|\right)^{-1} \widetilde{\sigma}\_{\mathbf{B}} \|\widetilde{\Lambda}\_{20}\| \|\widetilde{\Lambda}\_{10}\| + \|\widetilde{\Lambda}\_{20}\| \|\mathbb{B}\_{\mathbf{c}}\| (1 - \|\mathbb{K}\_{\mathbf{c}}D\_{\mathbf{c}}\|)^{-1} + (1 - \|\mathbb{K}\_{\mathbf{c}}D\_{\mathbf{c}}\|)^{-1} \widetilde{\sigma}\_{\mathbf{B}} \Big( \|\widetilde{\Lambda}\_{10}\| + \|\mathbb{K}\_{\mathbf{c}}\mathbb{C}\_{\mathbf{c}}\| \|\widetilde{\Lambda}\_{20}\| \Big) \end{split} \tag{A23}$$

*(iii) If the upper-bound* ε *is used for* ε<sup>2</sup> *and* ε3*, one gets that*

*Adc* ≤ <sup>ε</sup>*Bc*(<sup>1</sup> − *KcDc*)−<sup>1</sup> <sup>×</sup>[(*Kc*σ*<sup>C</sup>* <sup>+</sup>σ*C*σ*<sup>K</sup>* <sup>+</sup> *Cc*σ*K*)(<sup>1</sup> <sup>+</sup> *<sup>C</sup>* [*Kc*σ*<sup>D</sup>* <sup>+</sup>σ*D*σ*<sup>K</sup>* <sup>+</sup> *Dc*σ*K*]) <sup>+</sup> *KcCc<sup>C</sup>* [*Kc*σ*<sup>D</sup>* <sup>+</sup>σ*D*σ*<sup>K</sup>* <sup>+</sup> *Dc*σ*K*]] +ε (1 − *KcDc*) −1 σ*B*(<sup>1</sup> <sup>+</sup> *<sup>C</sup>* [*Kc*σ*<sup>D</sup>* <sup>+</sup>σ*D*σ*<sup>K</sup>* <sup>+</sup> *Dc*σ*K*])(*KcCc* <sup>+</sup> [*Kc*σ*<sup>C</sup>* <sup>+</sup>σ*C*σ*<sup>K</sup>* <sup>+</sup> *Cc*σ*K*] ) (A24)

*which also yields that Adc* ≤ <sup>ε</sup>*adc* <sup>+</sup> *<sup>o</sup>*(ε) *in the case that* ε < <sup>1</sup> *after grouping all the additive contributions of terms in* ε*<sup>i</sup> for i* ≥ 2 *into an additive bounded term, which converges to zero as* ε → 0*, where*

$$\begin{split} \widetilde{\boldsymbol{a}}\_{\rm dc} &= \left(1 - \|\boldsymbol{\mathcal{K}}\_{\rm c} \boldsymbol{D}\_{\rm c}\|\right)^{-1} \\ &\times \left[\widetilde{\boldsymbol{\sigma}}\_{\rm B} \|\boldsymbol{\mathcal{K}}\_{\rm c} \boldsymbol{\mathcal{C}}\_{\rm c}\| + \left(\|\boldsymbol{\mathcal{B}}\_{\rm c}\| + \|\boldsymbol{\mathcal{K}}\_{\rm c} \boldsymbol{\mathcal{C}}\_{\rm c}\|\right) \left(\|\boldsymbol{\mathcal{K}}\_{\rm c} \|\widetilde{\boldsymbol{\sigma}}\_{\rm C} + \|\boldsymbol{\mathcal{C}}\_{\rm c} \|\widetilde{\boldsymbol{\sigma}}\_{\rm K} + \boldsymbol{\mathcal{C}}\_{\rm c} \left[\|\boldsymbol{\mathcal{K}}\_{\rm c} \|\widetilde{\boldsymbol{\sigma}}\_{\rm D} + \|\boldsymbol{\mathcal{D}}\_{\rm c} \|\widetilde{\boldsymbol{\sigma}}\_{\rm K}\right]\right)\right] \end{split} \tag{A25}$$


Now, the following technical result follows directly from Lemma B2 and Theorem A1 (i):

**Lemma B3.** *Define* ε<sup>2</sup> = 1/ ⎛ ⎜⎜⎜⎜⎝*adc* sup ω∈*R*0<sup>+</sup> (*i*ω*In* <sup>−</sup> *Acc*)−1<sup>2</sup> ⎞ ⎟⎟⎟⎟⎠ *from Equation (A25) and assume that KcDc*<sup>2</sup> < 1 *and that* ε<sup>2</sup> < min (ε , ε1, 1)*. Then, Adc is stable if Acc is stable and* ε ∈ [0 , ε2]*. By using Equations (A22) and (A23), a better bound of the maximum allowable Adc*<sup>2</sup> *is found as* ε ∈ [0 , ε3] *with* ε<sup>3</sup> < min (ε , ε20, 1) *and* ε<sup>20</sup> = √ *u*2+4*v*ω−*u* <sup>2</sup>*<sup>v</sup> .*

**Proof.** Note from Theorem A1 (i) that the *<sup>H</sup>*∞-norm of (*sIn* <sup>−</sup> *Acc*)−1*Adc* satisfies (*sIn* <sup>−</sup> *Acc*)−1*Adc*∞ <sup>&</sup>lt; 1, which is guaranteed if ε < ε2, then *Adc* is stable since *Acc* is stable. The result follows by taking also into account, in addition, the constraints in Equations (A15) and (A16) of Lemma B.1 by using <sup>ε</sup><sup>3</sup> <sup>≤</sup> <sup>ε</sup><sup>2</sup> <sup>≤</sup> <sup>ε</sup>. If the second power of <sup>ε</sup> is considered and the third one is upper-bounded as <sup>ε</sup><sup>3</sup> <sup>≤</sup> <sup>ε</sup>2, we examine the stability constraint *u*ε + *v*ε<sup>2</sup> < ω = 1/ sup ω∈*R*0<sup>+</sup> (*i*ω*In* <sup>−</sup> *Acc*)−1<sup>2</sup> by building the convex parabola

<sup>θ</sup>(ε) <sup>=</sup> *<sup>v</sup>*ε<sup>2</sup> <sup>+</sup> *<sup>u</sup>*<sup>ε</sup> <sup>−</sup> ω < 0 whose negative and positive zeros are <sup>ε</sup>1,2 <sup>=</sup> <sup>−</sup>*u*<sup>±</sup> √ *u*2+4*v*ω <sup>2</sup>*<sup>v</sup>* . Hence, the second part of the result. -

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Review*
