**Short Overview of Early Developments of the Hardy Cross Type Methods for Computation of Flow Distribution in Pipe Networks**

**Dejan Brki´c 1,2,3,\*,**† **and Pavel Praks 1,2,\***


Received: 9 March 2019; Accepted: 16 April 2019; Published: 16 May 2019

**Abstract:** Hardy Cross originally proposed a method for analysis of flow in networks of conduits or conductors in 1936. His method was the first really useful engineering method in the field of pipe network calculation. Only electrical analogs of hydraulic networks were used before the Hardy Cross method. A problem with flow resistance versus electrical resistance makes these electrical analog methods obsolete. The method by Hardy Cross is taught extensively at faculties, and it remains an important tool for the analysis of looped pipe systems. Engineers today mostly use a modified Hardy Cross method that considers the whole looped network of pipes simultaneously (use of these methods without computers is practically impossible). A method from a Russian practice published during the 1930s, which is similar to the Hardy Cross method, is described, too. Some notes from the work of Hardy Cross are also presented. Finally, an improved version of the Hardy Cross method, which significantly reduces the number of iterations, is presented and discussed. We also tested multi-point iterative methods, which can be used as a substitution for the Newton–Raphson approach used by Hardy Cross, but in this case this approach did not reduce the number of iterations. Although many new models have been developed since the time of Hardy Cross, the main purpose of this paper is to illustrate the very beginning of modeling of gas and water pipe networks and ventilation systems. As a novelty, a new multi-point iterative solver is introduced and compared with the standard Newton–Raphson iterative method.

**Keywords:** Hardy Cross method; pipe networks; piping systems; hydraulic networks; gas distribution

#### **1. Introduction**

Hardy Cross solved the problem of distribution of flow in networks of pipes in his article "Analysis of Flow in Networks of Conduits or Conductors" [1] published on 13 November 1936.

Networks of pipes are nonlinear systems since the relation between flow and pressure is not linear. On the contrary, the relation between current and voltage in electrical networks with regular resistors is governed by the linear Ohm's law. Electrical circuits with diodes as well as hydraulic networks are nonlinear systems where resistance depends on current and voltage, i.e., on flow and pressure, respectively [2]. Nonlinear electrical circuits are electrical circuits containing nonlinear components. Nonlinear components can be resistive, capacitive, and inductive.

The distribution of flow in a network of pipes depends on the known inputs and consumptions at all nodes, on the given geometry of pipes, and on network topology. A stable state of flow in a network must satisfy Kirchhoff's laws, which are statements of the conservation of mass and energy. Although in theory an indefinite number of flow distributions that satisfy that the conservation of mass is possible, only one distribution from this set also satisfies the conservation of energy for all closed paths formed by pipes in the network. This state is unique for the given network and in- and outflows [3].

Since the relation between flow and pressure is not linear, Hardy Cross used a relation between an increment of flow and an increment of pressure, as this relation is linear for the given quantity of flow. If, however, the increments are fairly large, this linear relation is somewhat in error, such as for gas compressible flow. However, if the pressure drop in pipes is minor, such as in a municipality network for natural gas distribution, the Hardy Cross method can be used without significant errors [4–6]. Moreover, the Hardy Cross method can also be used for water pipe networks (district heating [7] and cooling networks [8]) and ventilation systems [9,10] (a related formulation is presented in Appendix).

The Hardy Cross method is an iterative method, i.e., a method using successive corrections [4]. Lobaˇcev and Andrijašev in the 1930s, writing in Russian, offered similar methods [11,12]. Probably because of the language barrier and the political situation in Soviet Russia, Hardy Cross was not aware of Lobaˇcev and Andrijašev's contributions.

Today, engineers use the most improved version of the Hardy Cross method (the Δ*Q* method [13]; for Δ*p*, see [14]), which analyzes the whole looped network of pipes simultaneously [15].

As a novel approach presented for the first time here, we tested multi-point iterative methods [16,17] that can be used as a substitution for the Newton–Raphson approach used by Hardy Cross. This approach, however, did not in this case reduce the number of required iterations to reach the final balanced solution.

One example of the pipe network for distribution of gas is analyzed using the original Hardy Cross method [1] in Section 3.1, its related equivalent from Russian literature [11,12] in Section 3.2, the improved version of the Hardy Cross method [15,17,18] in Section 3.3, and finally the approach which uses multi-point iterative methods instead of the commonly used Newton–Raphson method in Section 3.4.

#### **2. Network Piping System; Flow Distribution Calculation**

#### *2.1. Topology of the Network*

The first step in solving a pipe network problem is to make a network map showing pipe diameters, lengths and connections between pipes (nodes). Sources of natural gas supply and consumption rates have to be assigned to nodes. For convenience in locating pipes, code numbers are assigned to each pipe and closed loop of pipes (represented by roman numbers for loops in Figure 1). Pipes on the network periphery are common to one loop and those in the network interior are common to two loops. Figure 1 is an example of a pipe network for distribution of natural gas for consumption in households.

**Figure 1.** The network of pipes for natural gas distribution for domestic consumption.

The next step is to write the initial gas flow distribution through pipes in the network. This distribution must be done according to Kirchhoff's first law. The choice of initial flows is not critical, and the criterion should satisfy Kirchhoff's first law for every node in the network [3]. The total gas flow arriving at a node equals the total gas flow that leaves that node. The same

conservation law is also valid for the whole network in total (except for gas input and output nodes that cannot be changed during calculations; see consumption nodes in Figure 1). The sum of pseudo-pressure drops along any closed path must be approximately zero for the network to be in balance according to Kirchhoff's second law. In this paper, the flow distribution, which satisfies both of Kirchhoff's laws, will be calculated using the Hardy Cross iterative method.

#### *2.2. A Hydraulic Model*

The Renouard formula; Equation (1) best fits a natural gas distribution system built with polyvinyl chloride (PVC) pipes [19,20]. Computed pressure drops are always less than the actual drop since the maximal consumption occurs only during extremely severe winter days [21,22].

$$f = \Delta \overline{p}^2 = p\_1^2 - p\_2^2 = 4810 \cdot \frac{\rho\_{r^\cdot L^\cdot Q}^{1.82}}{D^{4.82}} \tag{1}$$

where *f* is a function of pressure, <sup>r</sup> is relative gas density (dimensionless), ρ*<sup>r</sup>* = 0.64, *L* is the pipe length (m), *D* is the pipe diameter (m), *Q* is flow (m3/s), and *p* is pressure (Pa).

As shown in Appendix A, other formulas are used in the case of waterworks systems [23,24] and ventilation networks [7].

Regarding the Renouard formula (Equation (1)), one has to be careful since the pressure drop function, *f*, does not relate pressure drop, but actually the difference of the quadratic pressure at the input and the output of the pipe. This means that \$ <sup>Δ</sup>*p*<sup>2</sup> <sup>=</sup> *p*2 <sup>1</sup> <sup>−</sup> *<sup>p</sup>*<sup>2</sup> <sup>2</sup> is not actually pressure drop despite using the same unit of measurement, i.e., the same unit is used as for pressure (Pa). The parameter \$ <sup>Δ</sup>*p*<sup>2</sup> can be noted as a pseudo-pressure drop. In fact, the gas is actually compressed, and hence that volume of the gas is decreased, and then such a compressed volume of the gas is conveying with a constant density through the gas distribution pipeline. The operating pressure for a typical distribution gas network is 4 <sup>×</sup> <sup>10</sup><sup>5</sup> Pa abs, i.e., 3 <sup>×</sup> 105 Pa gauge, and accordingly the volume of the gas decreases four times compared to the volume of the gas in normal (or standard) conditions. Pressure in the Renouard formula is for normal (standard) conditions.

The first derivative *f'* of the Renouard relation (Equation (2)), where the flow is treated as a variable, is used in the Hardy Cross method.

$$f' = \frac{\partial f(Q)}{\partial Q} = 1.82 \cdot 4810 \cdot \frac{\rho\_r \cdot L \cdot Q^{0.82}}{D^{4.82}} \tag{2}$$

First assumed gas flow in each pipe is listed in the third column of Table 1. The plus or minus sign preceding flow indicates the direction of flow through the pipe for the particular loop [18,25]. A plus sign denotes counterclockwise flow in the pipe within the loop, while the minus sign clockwise. The loop direction can be chosen to be clockwise or counterclockwise (in Figure 1, all loops are counterclockwise).

#### **3. The Hardy Cross Method; Di**ff**erent Versions**

The Hardy Cross method is presented here: the original approach in Section 3.1, a version of the Hardy Cross method from Russian practice in Section 3.2, the modified Hardy Cross method in Section 3.3, and finally the method that uses multi-point iterative procedures instead the Newton–Raphson method and which can be implemented in all the aforementioned methods.

#### *3.1. The Hardy Cross Method; Original Approach*

The pressure drop function for each pipe is listed in Table 1 (for initial flow pattern, in the fourth column). The sign in front of the pressure drop function shown in the fourth column is the same as for flow from the observed iteration. The fifth column of Table 1 includes the first derivatives of the pressure drop function, where the flow is treated as a variable. The column of the function of pressure drops is computed algebraically, while the column of the first derivatives is estimated numerically for each loop. Flow correction Δ*Q* has to be computed for each loop *x* (Equation (3)).

$$
\Delta Q\_{\rm x} = \left(\frac{\sum \pm f}{\left|f'\right|}\right)\_{\rm x} = \left(\frac{\sum \pm 4810 \cdot \frac{\rho\_r \cdot L \cdot Q^{1.82}}{D^{4.82}}}{\sum \left|1.82 \cdot 4810 \cdot \frac{\rho\_r \cdot L \cdot Q^{0.82}}{D^{4.82}}\right|}\right)\_{\rm x} \tag{3}
$$

For the network from Figure 1, flow corrections for the first iteration in each loop can be calculated using Equation (4).

$$\begin{pmatrix} \left( -f\_1' \right) + \left| f\_7' \right| + \left| f\_8' \right| + \left| f\_9' \right| + \left| -f\_{10}' \right| + \left| -f\_{12}' \right| \\ \left( -f\_2' \right) + \left| -f\_{11}' \right| + \left| f\_{12}' \right| \end{pmatrix} \Delta Q\_{II} = -f\_1 + f\_7 + f\_8 + f\_9 - f\_{10} - f\_{12} \\ \Delta Q\_{II} = -f\_1 + f\_7 + f\_{12} \\ \Delta Q\_{III} = \begin{pmatrix} -f\_2' \end{pmatrix} + \begin{pmatrix} -f\_{11}' \end{pmatrix} + \begin{pmatrix} f\_{12}' \end{pmatrix} \\ \Delta Q\_{II} = -f\_2 - f\_{11} + f\_{12} \\ \left( -f\_3' \right) + \left| f\_{13}' \right| + \left| -f\_{14}' \right| + \left| f\_{15}' \right| + \left| f\_{16}' \right| + f\_{11} - f\_{14} \\ \left( \left| f\_5' \right| + \left| -f\_9' \right| + \left| f\_{13}' \right| + \left| f\_{14}' \right| \right) \Delta Q\_{IV} = f\_5 - f\_9 + f\_{13} + f\_{14} \\ \left( \left| f\_6' \right| + \left| -f\_8' \right| + \left| f\_{13}' \right| \right) \Delta Q\_V = f\_6 - f\_8 + f\_{13} \end{pmatrix} \tag{4}$$

In the second iteration, the calculated correction Δ*Q* has to be added algebraically to the assumed gas flow (the first initial flow pattern). Further, the calculated correction Δ*Q* has to be subtracted algebraically from the gas flow computed in the previous iteration. This means that the algebraic operation for the first correction is the opposite of its sign, i.e., add when the sign is minus, and vice versa. A pipe common to two loops receives two corrections simultaneously. The first correction is from the particular loop under consideration, while the second one is from the adjacent loop, which the observed pipe also belongs to.


**Table 1.** Procedure for the solution of the flow problem for the network from Figure 1 using the modified Hardy Cross method (first two iterations)—First iteration.


**Table 1.** *Cont.*

<sup>a</sup> Pipe lengths, diameters and initial flow distribution are shown in Table 2 and Figure 1; <sup>b</sup> *f* calculated using the Renouard Equation (1); <sup>c</sup> *<sup>f</sup>* calculated using the first derivative of the Renouard, Equation (2), where flow is variable; <sup>d</sup> calculated using the matrix Equation (10) and entering <sup>Δ</sup>*<sup>Q</sup>* with the opposite sign (using the original Hardy Cross method for Iteration 1: Δ*QI*= +0.1126; Δ*QII*= −0.0088; Δ*QIII* = −0.0208; Δ*QIV* = −0.0892; Δ*QV*= −0.0902; using the Lobaˇcev method for Iteration 1: <sup>Δ</sup>*QI* <sup>=</sup> <sup>−</sup>0.1041; <sup>Δ</sup>*QII* <sup>=</sup> <sup>−</sup>0.0644; <sup>Δ</sup>*QIII*<sup>=</sup> <sup>−</sup>0.0780; <sup>Δ</sup>*QIV*= + 0.1069; <sup>Δ</sup>*QV*<sup>=</sup> <sup>−</sup>0.1824); <sup>e</sup> <sup>Δ</sup>*Q*<sup>2</sup> is <sup>Δ</sup>*Q*<sup>1</sup> from the adjacent loop; <sup>f</sup> the final calculated flow in the first iteration is used for the calculation in the second iteration, etc.; <sup>g</sup> if *Q* and *Q*<sup>1</sup> have a different sign, this means that the flow direction is opposite to that in the previous iteration, etc. (this occurs with the flow in pipe 13 between Iteration 3 and 4).

The upper sign after the second correction in Table 1 is plus if the flow direction in the mutual pipe coincides with the assumed orientation of the adjacent loop, and minus if it does not (Figure 2). The lower sign is the sign in front of correction Δ*Q* calculated for the adjacent loop (Figure 2).

**Figure 2.** Rules for the upper and lower sign (correction from the adjacent loop; second correction).

Details of the signs of corrections were reported by Brki´c [18] and Corfield et al. [25].

The algebraic operation for the second correction should be the opposite of its lower sign when its upper sign is the same as the sign in front of flow *Q*, and, as indicated by its lower sign, when its upper sign is opposite to the sign in front of flow *Q*.

The calculation procedure is repeated until the net algebraic sum of pressure functions around each loop is as close to zero as the desired degree of precision demands. This also means that the calculated corrections of flow and the change in calculated flow between two successive iterations is approximately zero. The pipe network is then in approximate balance and the calculation after the Hardy Cross can be terminated.

In the original Hardy Cross method, the corrections for the first iteration are:

$$
\Delta Q\_{I} = \frac{1575448179.8}{13987715480.9} = +0.1126$$

$$
\Delta Q\_{II} = \frac{-8424412.4}{957889226.7} = -0.0088$$

$$
\Delta Q\_{III} = \frac{-170493836.7}{8186058014.8} = -0.0208$$

$$
\Delta Q\_{IV} = \frac{-749453158.7}{8402810812.8} = -0.0892$$

and

$$
\Delta Q\_V = \frac{-325325177.5}{3605869136.3} = -0.09021
$$

#### *3.2. A Version of the Hardy Cross Method from Russian Practice*

As mentioned in the Introduction, two Russian authors, Lobaˇcev [11] and Andrijašev [12], proposed a similar method to Hardy Cross [1]. These two methods are also from the 1930s. It is not clear if Hardy Cross had been aware of the contribution of these two authors from Soviet Russia and vice versa, but most probably the answer to this question is no, for both sides. The main difference between the Hardy Cross and Andrijašev methods is that in the method of Andrijašev contours can be defined to include few loops. This strategy only complicates the situation, while the number of required iterations remains unchanged.

Further on, the Andrijašev method can be seen from the example in the paper of Brki´c [3]. Here, the method of Lobaˇcev is shown in more detail. In the Hardy Cross method, the influence of adjacent loops is neglected. The Lobaˇcev method takes into consideration this influence (Equation (5)):

+ -% % %−*f* 1 % % % + % % % *f* 7 % % % + % % % *f* 8 % % % + % % % *f* 9 % % % + % % %−*f* 10 % % % + % % %−*f* 12 % % % ·Δ*QI* + % % % *f* 12 % % %·Δ*QII* <sup>+</sup> % % % *f* 10 % % %·Δ*QIII* <sup>+</sup> % % % *f* 9 % % %·Δ*QIV* <sup>+</sup> % % % *f* 8 % % %·Δ*QV* <sup>=</sup> <sup>−</sup>*f*<sup>1</sup> <sup>+</sup> *<sup>f</sup>*<sup>7</sup> <sup>+</sup> *<sup>f</sup>*<sup>8</sup> <sup>+</sup> *<sup>f</sup>*<sup>9</sup> <sup>−</sup> *<sup>f</sup>*<sup>10</sup> <sup>−</sup> *<sup>f</sup>*<sup>12</sup> % % % *f* 12 % % %·Δ*QI* <sup>−</sup> -% % %−*f* 2 % % % + % % %−*f* 11 % % % + % % % *f* 12 % % % ·Δ*QII* − % % % *f* 11 % % %·Δ*QIII* <sup>=</sup> <sup>−</sup>*f*<sup>2</sup> <sup>−</sup> *<sup>f</sup>*<sup>11</sup> <sup>+</sup> *<sup>f</sup>*<sup>12</sup> + % % % *f* 10 % % %·Δ*QI* <sup>−</sup> % % % *f* 11 % % %·Δ*QII* <sup>−</sup> -% % %−*f* 3 % % % + % % % *f* 4 % % % + % % % *f* 10 % % % + % % % *f* 11 % % % + % % %−*f* 14 % % % ·Δ*QIII* − % % % *f* 14 % % %·Δ*QIV* <sup>=</sup> <sup>−</sup>*f*<sup>3</sup> <sup>+</sup> *<sup>f</sup>*<sup>4</sup> <sup>+</sup> *<sup>f</sup>*<sup>10</sup> <sup>+</sup> *<sup>f</sup>*<sup>11</sup> <sup>−</sup> *<sup>f</sup>*<sup>14</sup> + % % % *f* 9 % % %·Δ*QI* <sup>−</sup> % % % *f* 14 % % %·Δ*QIII* <sup>−</sup> -% % % *f* 5 % % % + % % %−*f* 9 % % % + % % % *f* 13 % % % + % % % *f* 14 % % % ·Δ*QIV* − % % % *f* 13 % % %·Δ*QV* <sup>=</sup> *<sup>f</sup>*<sup>5</sup> <sup>−</sup> *<sup>f</sup>*<sup>9</sup> <sup>+</sup> *<sup>f</sup>*<sup>13</sup> <sup>+</sup> *<sup>f</sup>*<sup>14</sup> + % % % *f* 8 % % %·Δ*QI* <sup>−</sup> % % % *f* 13 % % %·Δ*QIV* <sup>−</sup> -% % % *f* 6 % % % + % % %−*f* 8 % % % + % % % *f* 13 % % % ·Δ*QV* = *f*<sup>6</sup> − *f*<sup>8</sup> + *f*<sup>13</sup> ⎫ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭ (5)

In Equation (5), signs in front of terms from the left side of the equal signs have to be determined (this is much more complex than in the Hardy Cross method). Thus, in the Lobaˇcev method, if ( *<sup>f</sup>*)*<sup>x</sup>* <sup>&</sup>gt; 0, then the sign in front of -% % % *f* % % % *<sup>x</sup>* has to be positive, and vice versa (for the first iteration, this can be seen in Table 1; *fI* = +1575448179.8 > 0, *fII* = −8424412.4 < 0, *fIII* = −170493836.7 < 0, *fIV* = −749453158.7 < 0, *fV* = −325325177.5 < 0). The sign for the other terms (these terms are sufficient in the Hardy Cross method) are determined using further rules and the scheme in Figure 3.

**Figure 3.** Rules for terms from Lobaˇcev equations, which do not exist in the Hardy Cross method.

As shown in Figure 3 if ( *f*)*<sup>x</sup>* > 0 and if the assumed flow coincides with the loop direction, then the sign of flow in the adjacent pipe is negative and, if the flow does not coincide with the loop direction, then the sign of flow in the adjacent pipe is positive. Conversely, if ( *f*)*<sup>x</sup>* < 0 and if the assumed flow coincides with the loop direction, then the sign of flow in the adjacent pipe is positive and, if the flow does not coincide with the loop direction, then the sign of flow in the adjacent pipe is negative. This procedure determines the signs in the front of the flow corrections (Δ*Q*), which are shown in Figure 3 with black letters (and also in Figure 4 for our pipe network example).

**Figure 4.** Rules for terms from Lobaˇcev equations, which do not exist in the Hardy Cross method applied for the network from Figure 1.

If ( *f*)*<sup>x</sup>* from the adjacent loop is positive, while the loop direction and assumed flow do not coincide, the flow correction from the adjacent loop changes its sign, and, conversely, if ( *f*)*<sup>x</sup>* from the adjacent loop is positive, while the loop direction and assumed flow coincide, the flow correction from the adjacent loop does not change its sign. If ( *f*)*<sup>x</sup>* from the adjacent loop is negative, while the loop direction and assumed flow do not coincide, the flow correction from the adjacent loop does not change its sign, and, conversely, if ( *f*)*<sup>x</sup>* from the adjacent loop is negative, while the loop direction and assumed flow do not coincide, the flow correction from the adjacent loop changes its sign. These four parameters are connected in Figure 3 with the same colored lines. Flow corrections (Δ*Q*) shown in Figure 4 with different colors are used with the related signs in Equation (5). They are chosen in a similar way as explained in the example in Figure 3.

Thus, instead of the simple equations of the original Hardy Cross method, the system of equations has to be solved in the Lobaˇcev method (Equation (6)).


Underlined terms in Equation (6) do not exist in the Hardy Cross method.

In the Lobaˇcev method, corrections for the first iterations are <sup>Δ</sup>*Qx* <sup>=</sup> <sup>Δ</sup>(Δ*Qx*) <sup>Δ</sup> , where Δ for the first iteration is (Equation (7)):


while ΔQx for the first iteration is (Equation (8)).


The correction for the first loop in the first iteration is (Equation (9)).

$$
\Delta Q\_l = \frac{\Delta(\Delta Q\_l)}{\Delta} = \frac{-4.14 \cdot 10^{+47}}{+3.97 \cdot 10^{+48}} = -0.1041\tag{9}
$$

Other corrections in the first iteration are Δ*QII* = −0.0644, Δ*QIII* = −0.0780, Δ*QIV* = +0.1069 and Δ*QV* = −0.1824.

The Lobaˇcev method is more complex compared to the original Hardy Cross method. However, the number of required iterations is not reduced using the Lobaˇcev procedure compared with the original Hardy Cross procedure.

#### *3.3. The Modified Hardy Cross Method*

The Hardy Cross method can be noted in matrix form. The gas distribution network in Figure 1 has five independent loops (Equation (10)).

$$
\begin{bmatrix}
\Sigma[f\_I'] & 0 & 0 & 0 & 0 \\
0 & \Sigma[f\_{II}'] & 0 & 0 & 0 \\
0 & 0 & \Sigma[f\_{III}'] & 0 & 0 \\
0 & 0 & 0 & \Sigma[f\_{IV}'] & 0 \\
0 & 0 & 0 & 0 & \Sigma[f\_V']
\end{bmatrix} \ge \begin{bmatrix}
\Delta Q\_I \\
\Delta Q\_{II} \\
\Delta Q\_{III} \\
\Delta Q\_{IV} \\
\Delta Q\_{V}
\end{bmatrix} = \begin{bmatrix}
\Sigma \pm f\_I \\
\Sigma \pm f\_{II} \\
\Sigma \pm f\_{III} \\
\Sigma \pm f\_{IV} \\
\Sigma \pm f\_V
\end{bmatrix} \tag{10}
$$

Equation (4) provides for each particular loop in the network the same corrections as Equation (10) using matrix calculation. Epp and Fowler [15] improved the original Hardy Cross method [1] by replacing some of the zeroes in the non-diagonal terms of Equation (10). For example, if pipe 8 is mutual for loop I and V, the first derivative of the pressure drop function for the observed pipe, where flow treated as a variable, will be put with a negative sign in the first column and the fifth row, and also in the fifth column and the first row (Equation (11)).

$$
\begin{bmatrix}
\Sigma \left| f\_I' \right| & -f\_{12}' & -f\_{10}' & -f\_9' & -f\_8' \\
\end{bmatrix} \ge \begin{bmatrix}
\Delta Q\_I \\
\Delta Q\_{II} \\
\Delta Q\_{II} \\
\Delta Q\_{IV} \\
\Delta Q\_{IV} \\
\Delta Q\_V
\end{bmatrix} = \begin{bmatrix}
\Sigma \pm f\_I \\
\Sigma \pm f\_{II} \\
\Sigma \pm f\_{II} \\
\Sigma \pm f\_{IV} \\
\Sigma \pm f\_V \\
\Sigma \pm f\_V
\end{bmatrix} \tag{11}
$$

In the modified Hardy Cross method, corrections for the first iterations are shown in Equation (12), and the solutions are listed in Table 1.


This procedure significantly reduces the number of iterations required for the solution of the problem (Figure 5).

**Figure 5.** The number of required iterations for the solution using the original vs. the improved Hardy Cross method.

The first two iterations for the network in Figure 1 are shown in Table 1. Pipe diameters and lengths, as well as the first, assumed, and the final calculated flow distributions for the network in balance are shown in Table 2.

The gas velocity in the network is small (can be up to 10–15 m/s). The network can be the subject of diameter optimization (as in [4]), which can also be done by using the Hardy Cross method (diameter correction Δ*D* should be calculated for known and locked flow, where the first derivative of the Renouard function has to be calculated for diameter as a variable). The network should stay unchanged, even if planned gas consumption on nodes 5, 6, 8 and 10 increases, as pipes 4 and 13 will be useful thanks to increased gas flow.

Similar examples, but for water flow, can be seen in [26]. Optimization of pipe diameters in a water distributive pipe network using the same approach can be seen in [6].


**Table 2.** Pipe diameters and lengths, flows, and velocities of gas within pipes.

<sup>a</sup> Network from Figure 1 (flows are for normal pressure conditions; real pressure in the network is 4 × 105 Pa abs, i.e., 3 × <sup>10</sup><sup>5</sup> Pa); <sup>b</sup> chosen to satisfy Kirchhoff's first law for all nodes (dash arrows in Figure 1); <sup>c</sup> calculated to satisfy Kirchhoff's first law for all nodes and Kirchhoff's second law for all closed path formed by pipes (full errors in Figure 1); <sup>d</sup> the minus sign means that the direction of flow is opposite to the initial pattern for assumed flows.

#### *3.4. The Multi-Point Iterative Hardy Cross Method*

The here described multipoint method can substitute the Newton–Raphson iterative procedure used in all the above described methods. Recently, we successfully used this multipoint method for acceleration of the iterative solution of the Colebrook equation for flow friction modeling [16,17]. On the contrary, for the gas network example in Figure 1, the multipoint method requires the same number of iterations as the original Newton–Raphson procedure.

For the test, we used the three-point method from Džuni´c et al. [27]. Flow corrections Δ*QI* from Equations (10) and (11) from the first loop I should be calculated using the three-point procedure (Equation (13)):

$$\begin{cases} \Delta Q\_{I} = -\frac{(\sum \pm f)\_{i}}{\left(\sum [f']\_{i}\right)\_{i}}, \text{ where } i = 1;\\ \Delta Q\_{I} = -\frac{(\sum \pm f)\_{i}}{(\sum \pm f)\_{i} - 2 \cdot (\sum \pm f)\_{i+1}} \cdot \frac{(\sum \pm f)\_{i+1}}{\left(\sum [f']\_{i}\right)\_{i}}, \text{ where } i = 2;\\ \Delta Q\_{I} = -\frac{(\sum \pm f\_{i})\_{i+2}}{\left(\sum [f']\_{i}\right)\_{i} \left[1 - 2 \cdot \frac{(\sum \pm f)\_{i+1}}{\left(\sum [f']\_{i}\right)\_{i}} - \left(\frac{(\sum \pm f)\_{i+1}}{\left(\sum [f]\_{i}\right)\_{i}}\right)^{2}\right] \left[1 - \frac{(\sum \pm f)\_{i+2}}{\left(\sum [f]\_{i}\right)\_{i+1}}\right]}, \text{ where } i > 2; \end{cases} \tag{13}$$

Formulas of flow corrections Δ*QI* depend on the counter*i.* The algorithm starts from *i* = 1, in which the multipoint method is the same as the original Newton–Raphson procedure (Equation (13a)):

$$
\Delta Q\_I = -\frac{(\sum \pm f\_I)\_i}{\left(\sum |f'\_I|\right)\_i} \tag{13a}
$$

In the second iteration, *i* = 2, flow corrections Δ*QI* have a little bit more complicated form (Equation (13b)):

$$\Delta Q\_l = -\frac{(\sum \pm f\_l)\_i}{(\sum \pm f\_l)\_i - 2 \cdot (\sum \pm f\_l)\_{i+1}} \cdot \frac{(\sum \pm f\_l)\_{i+1}}{\left(\sum |f'\_l|\right)\_i} \tag{13b}$$

The symbol( ±*fI*)*<sup>i</sup>* represents stored values from the first iteration, whereas ( ±*fI*)*i*+<sup>1</sup> represents values from the second iteration.

For the third iteration, *i* = 3, flow corrections Δ*QI* have the most complicated form (Equation (13c)):

$$\Delta Q\_{I} = -\frac{(\sum \pm f\_{\mathbf{x}})\_{i+2}}{\left(\sum |f'\_{I}|\right)\_{i} \cdot \left[1 - 2 \cdot \frac{(\sum \pm f)\_{i+1}}{(\sum \pm f)\_{i}} - \left(\frac{(\sum \pm f)\_{i+1}}{(\sum \pm f)\_{i}}\right)^{2}\right] \cdot \left[1 - \frac{(\sum \pm f)\_{i+2}}{(\sum \pm f)\_{i+1}} \cdot \left[1 - 2 \cdot \frac{(\sum \pm f)\_{i+2}}{(\sum \pm f)\_{i}}\right] \right]}\tag{13c}$$

This iterative process can continue, as the formula from the third iteration is also used for iterations *i* = 4, 5, 6, 7, etc. This procedure should be done for all loops in the network separately (in our case for I, II, III, IV and V). However, to simplify calculations, derivative-free methods can be used [28,29].

#### **4. Conclusions**

Hardy Cross simplified mathematical modeling of complex problems in structural and hydraulic engineering long before the computer age. Moment distributions in indeterminate concrete structures described with differential equations were too complex for the time before computers. Hardy Cross later applied these finding from structural analysis to balancing of flow in pipe networks. He revolutionized how the profession addressed complicated problems. Today, in engineering practice, the modified Hardy Cross method proposed by Epp and Fowler [15] is used rather than the original version of the Hardy Cross method [1]. Methods proposed by Hamam and Brameller [30], as well as by Wood and Charles [31] and Wood and Rayes [32], are used in common practice [33], too. Moreover, the node-oriented method proposed by Shamir and Howard [34] is also based on the Hardy Cross method.

Professional engineers use a different kind of looped pipeline in professional software [35], but. even today, engineers invoke the name of Hardy Cross with awe. When petroleum and natural gas or civil engineers have to figure out what is happening in looped piping systems [36], they inevitably turn to what is generally known as the Hardy Cross method. The original Hardy Cross method is still extensively used for teaching and learning purpose [6]. Here, we introduced into the Hardy Cross method the multi-point iterative approach instead of the Newton–Raphson iterative approach, but it does not affect the number of required iterations to reach the final solution in our case.

The view of Hardy Cross was that engineers lived in the real world with real problems and that it was their job to come up with answers to questions in design tasks, even if initial approximations were involved. After Hardy Cross, the essential idea which he wished to present involves no mathematical relations except the simplest arithmetic.

For example, ruptures of pipes with leakage can be detected using the Hardy Cross method because every single-point disturbance affects the general distribution of flow and pressure [37,38].

This paper has the purpose of illustrating the very beginning of modeling of gas or water pipe networks. As noted by Todini and Rossman [39], many new models have been developed since the time of Hardy Cross.

Some details about the life and work of Hardy Cross are given in Appendix B.

**Author Contributions:** The paper is a product of the joint efforts of the authors, who worked together on models of natural gas distribution networks. P.P. has a scientific background in applied mathematics and programming while D.B. has a background in control and applied computing in mechanical and petroleum engineering. D.B. performed calculations with advice from P.P. who has extensive experience with the implementation of probabilistic gas network modeling tools.

**Funding:** Resources to cover the Article Processing Charge were provided by the European Commission. This article is registered as Pubsy JRC116624 in the internal system for publications of the European Commission, Joint Research Center.

**Acknowledgments:** We acknowledge support from the European Commission, Joint Research Centre (JRC), Directorate C: Energy, Transport and Climate, Unit C3: Energy Security, Distribution and Markets, and we especially thank Marcelo Masera and Ricardo Bolado-Lavin for their scientific supervision and approvals. We also thank John Cawley from IT4Innovations, who as a native speaker kindly checked the correctness of English expressions throughout the paper. This work was partially supported by the Ministry of Education, Science and Technological Development of the Republic of Serbia through the project iii44006 and by the Ministry of Education, Youth and Sports of the Czech Republic through the National Programme of Sustainability (NPS II) project "IT4Innovations excellence in science-LQ1602". Dejan Brki´c acknowledges support from the mentioned Czech's project LQ1602 for funding of his stay in Ostrava at IT4Innovations.

**Conflicts of Interest:** The authors declare no conflict of interest. Neither the European Commission, the VŠB— Technical University of Ostrava, or the Alfatec nor any person acting on their behalf is responsible for any use which might be made of this publication.

#### **Nomenclature**

The following symbols are used in this paper:


#### **Appendix A Hydraulic Models for Water Pipe Networks and for Ventilation Systems**

To relate pressure *p* [40] with flow *Q*, instead of Equation (1), which is used for gas distribution networks in municipalities, for water distribution the Darcy–Weisbach correlation and Colebrook equation are recommended (Equation (A1)) [23,41–47], and for ventilation systems, the Atkinson equation (Equation (A2)) [9]:

$$\begin{split} \frac{1}{\sqrt{\lambda}} &= -2 \cdot \log\_{10} \left( \frac{2.51}{\text{Re}} \cdot \frac{1}{\sqrt{\lambda}} + \frac{\varepsilon}{3.71 \cdot D} \right) \\ \Delta p &= \frac{8 \cdot p \cdot \lambda \cdot \mathcal{Q}^2}{\pi^2 \cdot D^8} \end{split} \tag{A1}$$

$$
\Delta p = \frac{\rho}{2 \cdot \mathbb{C}\_d^2 \cdot \mathbb{A}^2} \cdot \mathbb{Q}^2 \tag{A2}
$$

#### **Appendix B The Life and Work of Hardy Cross**

Hardy Cross (1885–1959) was one of America's most brilliant engineers [48–55]. He received a BSc degree in arts in 1902 and BSc degree in science in 1903, both from Hampden-Sydney College, where he taught English and Mathematics. Hardy Cross was also awarded a BSc degree in 1908 from Massachusetts Institute of Technology and an MCE degree from Harvard University in 1911, both in civil engineering. He taught civil engineering at Brown University from 1911 to 1918. He left teaching twice to become involved in the practice of structural and hydraulic engineering, from 1908 to 1909, and from 1918 to 1921. The most creative years of Hardy Cross were spent at the University of Illinois in Champaign-Urbana where he was a professor of structural engineering from 1921 to 1937. His famous article "Analysis of flow in networks of conduits or conductors" was published in 1936 in Urbana Champaign University Illinois Bulletin; Engineering Experiment Station number 286 [1]. His name is also famous in the field of structural engineering [53–55]. He developed the moment distribution method for statically indeterminate structures in 1932 [56]. This method has been superseded by more powerful procedures, however the moment distribution method made possible the efficient and safe design of many reinforced concrete buildings for the duration of an entire generation. Furthermore, the solution of the here discussed pipe network problems was a byproduct of his explorations in structural analysis. Later, Hardy Cross was Chair of the Department of Civil Engineering at Yale, from 1937 to the early 1950s.

Related to the moment distribution method for statically indeterminate structures developed by Hardy Cross in 1932 [56], it need to be noted that in 1922 and 1923, Konstantin A. Cališev, emigrant from ˇ Soviet Russia, writing in Serbian, offered a similar method of solving the slope deflection equations by successive approximations [57–60]. The method of Hardy Cross is an improved version of the Cališev's method, but with important circumstance that Hardy Cross most probably was not aware of ˇ Cališev's contributions. As noted in [ ˇ 49]: "It was Hardy Cross's genius that he recognized he could bypass adjusting rotations to get to the moment balance at each and every node.", which was the part that Cališev did not developed. ˇ

#### **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* **Multi-Switching Combination Synchronization of Three Fractional-Order Delayed Systems**

**Bo Li 1, Yun Wang 1,\* and Xiaobing Zhou <sup>2</sup>**


Received: 19 July 2019; Accepted: 13 August 2019; Published: 15 October 2019

**Abstract:** Multi-switching combination synchronization of three fractional-order delayed systems is investigated. This is a generalization of previous multi-switching combination synchronization of fractional-order systems by introducing time-delays. Based on the stability theory of linear fractional-order systems with multiple time-delays, we propose appropriate controllers to obtain multi-switching combination synchronization of three non-identical fractional-order delayed systems. In addition, the results of our numerical simulations show that they are in accordance with the theoretical analysis.

**Keywords:** multi-switching combination synchronization; time-delay; fractional-order; stability

#### **1. Introduction**

Fractional calculus has attracted researchers from various fields due to fractional dimensions widely existing in nature and engineering fields [1–3]. Compared to the integer-order dynamical systems, the fractional-order counterparts can exhibit more complex dynamical behaviors. Some of the researches on integer-order dynamical systems can be generalized to fractional-order dynamical systems. Fractional-order dynamical systems have been widely investigated, such as synchronization [4], identification [5], stabilization [6] and approximate entropy analysis [7,8]. Time-delay is a frequently encountered phenomenon in real applications, such as physical, communication, economical, pneumatic and biological systems [9]. Introducing time-delay into a system can enrich its dynamic characteristics and describe a real-life phenomenon more precisely. Thus, the fractional-order delayed differential equation (FDDE) is becoming a hot topic for scientists and engineers, and it has many theoretical and practical applications [10]. Nowadays, the chaotic behavior and synchronization of FDDE attract intensive research interests. In [11], Bhalekar et al. introduced the fractional-order delayed Liu system. The fractional-order delayed financial system was presented in [12], and hybrid projective synchronization between the aforementioned two systems was achieved in [13]. The fractional-order delayed Chen system was proposed in [14], while its adaptive synchronization was investigated in [15]. The fractional-order delayed porous media was proposed in [16]. In [17], a fractional-order delayed Newton–Leipnik system was taken as an example to present intermittent synchronizing delayed fractional nonlinear system.

Due to its wide applications in secure communication, synchronization of fractional-order delayed chaotic systems are extensively investigated [18–21]. However, all the above-mentioned and other synchronization schemes are in traditional drive–response ways, which only have unique drive and response systems. Recently, Luo et al. [22] extended the traditional drive–response synchronization to combination synchronization, which has two drive systems and one response system. Compared to the drive–response synchronization, combination synchronization has stronger anti-decode and anti-attack abilities in image encryption and secure communication, in which origin messages are split into two parts and each part can be embedded into two separate drive systems. There are many works on combination synchronization [23–26]. To further strengthen the security in secure communication, Vincent et al. [27] proposed multi-switching combination synchronization scheme, in which the two drive systems are synchronized with the response system in different states. Based on the nonlinear control technique, Zheng [28] studied multi-switching combination synchronization of three non-identical chaotic systems. Khan [29] investigated adaptive multi-switching combination synchronization among three non-identical chaotic systems. In [30], Ahmad et al. proposed globally exponential multi-switching combination synchronization scheme, and applied it to secure communications. Multi-switching combination synchronization was applied in encrypted audio communication in [31]. The previous work on multi-switching combination synchronization schemes are based on integral-order chaotic systems. To improve the security in these synchronization schemes, based on fractional-order chaotic systems, Bhat et al. [32] extended the work in [27] to study multi-switching combination synchronization among three non-identical fractional-order chaotic systems. Khan et al. [33] investigated multi-switching combination-combination synchronization among a class of four non-identical fractional-order chaotic systems. In multi-switching combination-combination synchronization scheme, the state variables of two drive systems synchronize with different state variables of two response systems simultaneously, which makes the security of this scheme higher than that in [32].

Since time-delay is a frequently encountered phenomenon in real applications, and the time-delay can be used as an additional parameter in synchronization to increase security in secure communication, we consider the work in [32] to investigate multi-switching combination synchronization scheme for non-identical fractional-order delayed systems by introducing time-delays in fractional-order systems.

The rest of this paper is organized as follows. In Section 2, the concept of fractional calculus and the stability theory of linear fractional-order systems with multiple time-delays are briefly introduced. The multi-switching combination synchronization scheme of three non-identical fractional-order delayed systems is analyzed in Section 3. In Section 4, numerical simulations performed using MATLAB are presented. Finally, conclusions are drawn in Section 5.

#### **2. Preliminaries**

Fractional calculus is a generalization of integration and differentiation to non-integer order fundamental operator *aD<sup>r</sup> <sup>t</sup>* , which is defined as

$$\_aD\_t^r = \begin{cases} \frac{d^r}{dt^r} & : r > 0, \\\ 1 & : r = 0, \\\ \int\_a^t (d\tau)^{-r} : r < 0. \end{cases} \tag{1}$$

There are several different definitions for the fractional-order differential operator [34]. Because the Caputo definition is easy to understand and is frequently used in the literature, we apply this definition in this paper, which is

$$\_aD\_t^r f(t) = \frac{1}{\Gamma(n-r)} \int\_a^t \frac{f(\tau)}{t-\tau}^{r-n+1} d\tau \tag{2}$$

where 1 < *r* < *n*.

**Lemma 1** ([35])**.** *Suppose <sup>f</sup>*(*t*) <sup>∈</sup> *<sup>C</sup><sup>α</sup> <sup>a</sup>* ([*a*, *b*]), *D<sup>α</sup> <sup>a</sup> <sup>f</sup>*(*t*) <sup>∈</sup> *<sup>C</sup><sup>β</sup> <sup>a</sup>* ([*a*, *b*]), *α* > 0, *β* > 0, *m* − 1 < *β* < *m*, *n* − 1 < *α* < *n, then C*

$$\prescript{C}{a}{D}\_{t}^{\beta} \left( \prescript{C}{a}{D}\_{t}^{a} f(t) \right) = \prescript{C}{a}{D}\_{t}^{a} \left( \prescript{C}{a}{D}\_{t}^{\beta} f(t) \right) = \prescript{C}{a}{D}\_{t}^{a + \beta} f(t) \tag{3}$$

*Appl. Sci.* **2019**, *9*, 4348

**Lemma 2** ([35])**.** *Consider*

$$\, \_a^C D\_t^a \mathbf{x}^n = \frac{\Gamma(n+1)\mathbf{x}^{n-a}}{\Gamma(n+1-a)} \, \_a^C D\_t^a \mathbf{x}$$

**Lemma 3** ([35])**.** *Let φ*(*t*) ∈ *R be a continuous and derivable function. Then, for any time instant t* ≥ *t*<sup>0</sup>

$$\frac{1}{2}D^q\phi^2(t) = \phi(t)D^q\phi(t),\ \forall q \in (0,1)\tag{5}$$

Given the following n-dimensional linear fractional-order system with multiple time-delays [36]:

$$\begin{cases} D^{q\_1}y\_1(t) = a\_{11}y\_1(t - \tau\_{11}) + a\_{12}y\_2(t - \tau\_{12}) \\ \qquad + \cdots + a\_{1n}y\_n(t - \tau\_{1n}), \\ D^{q\_2}y\_2(t) = a\_{21}y\_1(t - \tau\_{21}) + a\_{22}y\_2(t - \tau\_{22}) \\ \qquad + \cdots + a\_{2n}y\_n(t - \tau\_{2n}), \\ \vdots \\ D^{q\_n}y\_n(t) = a\_{n1}y\_1(t - \tau\_{n1}) + a\_{n2}y\_2(t - \tau\_{n2}) \\ \qquad + \cdots + a\_{nn}y\_n(t - \tau\_{nn}), \end{cases} (6)$$

where *qi* ∈ (0, 1) is the fractional-derivative order, *yi*(*t*) is the state, and *τij* > 0 is the time-delay, the initial value *yi*(*t*) = *φi*(*t*) is given by −*maxτij* = −*τmax* ≤ *t* ≤ 0, *A* = [*aij*] ∈ *Rn*×*<sup>n</sup>* is the coefficient matrix.

Performing Laplace transform on the system in Equation (6) yields

$$
\Delta(\mathbf{s}) \cdot \mathbf{Y}(\mathbf{s}) = b(\mathbf{s}),
\tag{7}
$$

where *Y*(*s*)=(*Y*1(*s*), *Y*2(*s*), ..., *Yn*(*s*))*<sup>T</sup>* is the Laplace transform of *y*(*t*) = (*y*1(*t*), *y*2(*t*), ..., *yn*(*t*))*T*, *b*(*s*)=(*b*1(*s*), *b*2(*s*), ..., *bn*(*s*))*<sup>T</sup>* is the remaining non-linear part, the characteristic matrix of the system in Equation (6) is

$$
\Delta(s) = \begin{pmatrix}
 s^{q\_1} - a\_{11}e^{-s\tau\_{11}} & -a\_{12}e^{-s\tau\_{12}} & \cdots & -a\_{1n}e^{-s\tau\_{1n}} \\
 -a\_{21}e^{-s\tau\_{21}} & s^{q\_2} - a\_{22}e^{-s\tau\_{22}} & \cdots & -a\_{2n}e^{-s\tau\_{2n}} \\
 \vdots & \vdots & \ddots & \vdots \\
 -a\_{n1}e^{-s\tau\_{n1}} & -a\_{n2}e^{-s\tau\_{n2}} & \cdots & s^{q\_n} - a\_{nn}e^{-s\tau\_{nn}}
\end{pmatrix}.\tag{8}
$$

**Theorem 1** ([36])**.** *If all the roots of the characteristic equation* det(Δ(*s*)) = 0 *have negative real parts, then the zero solution of the system in Equation (6) is Lyapunov globally asymptotically stable.*

**Corollary 1** ([36])**.** *If q*<sup>1</sup> = *q*<sup>2</sup> = ··· = *qn* = *β* ∈ (0, 1)*, all the eigenvalues λ of the coefficient matrix A satisfy* | arg(*λ*)| > *βπ*/2*, and the characteristic equation* det(Δ(*s*)) = 0 *has no purely imaginary roots for any τij* > 0, *i*, *j* = 1, 2, ... , *n*, *then the zero solution of the system in Equation (6) is Lyapunov globally asymptotically stable.*

#### **3. Multi-Switching Combination Synchronization Scheme**

Multi-switching combination synchronization among three non-identical fractional-order delayed systems is investigated in this section.

The two drive systems are

$$\begin{aligned} D^{\mathfrak{a}}\mathbf{x}(t) &= \mathbf{x}(t) + \mathbf{x}(t-\tau) + A(\mathbf{x}(t), \mathbf{x}(t-\tau)), \\ \mathbf{x}(t) &= \mathbf{x}(0), t \in [-\tau, 0], \end{aligned} \tag{9}$$

and

$$\begin{aligned} D^{\mu}y(t) &= y(t) + y(t-\tau) + B(y(t), y(t-\tau)), \\ y(t) &= y(0), t \in [-\tau, 0]. \end{aligned} \tag{10}$$

The response system is

$$\begin{aligned} D^{\mathfrak{a}}z(t) &= z(t) + z(t-\tau) + \mathbb{C}(z(t), z(t-\tau)) + \mathbb{U}, \\ z(t) &= z(0), t \in [-\tau, 0], \end{aligned} \tag{11}$$

in which, *α* ∈ (0, 1) is the fractional order, *τ* > 0 is the time-delay, *U* = (*U*1, ..., *Un*) is the controller vector, *<sup>x</sup>* = (*x*1, ..., *xn*)*<sup>T</sup>* ∈ *<sup>R</sup>n*, *<sup>y</sup>* = (*y*1, ..., *yn*)*<sup>T</sup>* ∈ *<sup>R</sup><sup>n</sup>* and *<sup>z</sup>* = (*z*1, ..., *zn*)*<sup>T</sup>* ∈ *<sup>R</sup><sup>n</sup>* are state vectors, and *<sup>A</sup>* : *<sup>R</sup>*2*<sup>n</sup>* → *<sup>R</sup>n*, *<sup>B</sup>* : *<sup>R</sup>*2*<sup>n</sup>* → *<sup>R</sup><sup>n</sup>* and *<sup>C</sup>* : *<sup>R</sup>*2*<sup>n</sup>* → *<sup>R</sup><sup>n</sup>* are continuous vector functions.

Define the error state as *eklm* = *fkzk* − *glxl* − *hmym*(*k*, *l*, *m* = 1, ..., *n*). Then, we have the error state vector

$$\mathbf{e}(t) = \mathbf{F}\mathbf{z} - \mathbf{G}\mathbf{x} - \mathbf{H}\mathbf{y},\tag{12}$$

where *<sup>e</sup>*(*t*) is the vector form of *eklm*, *<sup>F</sup>* = *diag*{ *<sup>f</sup>*1, *<sup>f</sup>*2, ... , *fn*} ∈ *<sup>R</sup>n*×*n*, *<sup>G</sup>* = *diag*{*g*1, *<sup>g</sup>*2, ... , *gn*} ∈ *<sup>R</sup>n*×*<sup>n</sup>* and *<sup>H</sup>* = *diag*{*h*1, *<sup>h</sup>*2, ... , *hn*} ∈ *<sup>R</sup>n*×*<sup>n</sup>* are real scaling matrix. Accordingly, *eklm*(*<sup>t</sup>* − *<sup>τ</sup>*) = *fkzk*(*<sup>t</sup>* − *<sup>τ</sup>*) − *glxl*(*t* − *τ*) − *hmym*(*t* − *τ*).

**Definition 1** ([27])**.** *The systems in Equations (9) and (10) and the system in Equation (11) are defined to be multi-switching combination synchronization if F*, *G*, *H are non-zeros, and k* = *l* = *m*, *k* = *l* = *m*, *k* = *l* = *m*, *k* = *m* = *l, such that:*

$$\lim\_{t \to +\infty} \parallel \varepsilon(t) \parallel = \lim\_{t \to +\infty} \parallel Fz - Gx - Hy \parallel = 0 \tag{13}$$

*where* . *represents the matrix norm.*

**Remark 1.** *If k* = *l* = *m, the systems in Equations (9) and (10) and the system in Equation (11) are defined to be combination synchronization [22].*

**Remark 2.** *If the scaling matrix F* = 0*, G* = 0 *or H* = 0*, the multi-switching combination synchronization mentioned above is simplified to multi-switching hybrid projective synchronization.*

From the systems in Equations (9)–(11), we have the error system as follows

$$D^a e(t) = FD^a z(t) - GD^a x(t) - HD^a y(t) \tag{14}$$

To achieve multi-switching combination synchronization among the above systems, a non-linear controller is constructed:

$$dI = \mathcal{R}c(t) + GA(\mathbf{x}(t), \mathbf{x}(t-\tau)) + HB(y(t), y(t-\tau)) - FC(z(t), z(t-\tau)),\tag{15}$$

where *<sup>K</sup>*˜ <sup>=</sup> *<sup>K</sup>* <sup>−</sup> *<sup>I</sup>*, *<sup>I</sup>* is an n-dimensional unit matrix, *<sup>K</sup>* <sup>=</sup> *diag*{*k*1, *<sup>k</sup>*2, ... , *kn*} is a feedback gain matrix. Substituting the systems in Equations (9)–(11) and (15) into the system in Equation (14), we have

$$D^{\mathfrak{a}}e(t) = (\mathbb{k} + l)e(t) + e(t - \tau) = \mathcal{K}e(t) + e(t - \tau). \tag{16}$$

Thus, the multi-switching combination synchronization between the systems in Equations (9) and (10) and the system in Equation (11) is changed into the analysis of the asymptotical stability of the system in Equation (16).

In light of Corollary 1, a sufficient condition to achieve multi-switching combination synchronization between the systems in Equations (9) and (10) and the system in Equation (11) is obtained as follows.

**Proposition 1.** *Multi-switching combination synchronization between the systems in Equations (9) and (10) and the system in Equation (11) can be achieved if there exists a matrix K* = *diag*{*k*1, *k*2, ... , *kn*} *in the system in Equation (16) such that ki* < −1/ sin(*απ*/2), (*i* = 1, 2, . . . , *n*)*.*

**Proof.** For the system in Equation (16), *A* = *K* + *I* is the coefficient matrix. Since *ki* < −1/ sin(*απ*/2), *α* ∈ (0, 1), the eigenvalues of *A* are *λ<sup>i</sup>* = *ki* + 1 < 0, (*i* = 1, 2, ... , *n*). Then, |arg(*λ*)| > *π*/2 > *απ*/2 holds.

Performing Laplace transform on the system in Equation (16) yields

$$
\Delta(s) \cdot E(s) = s^{a-1} e(0) + e(0) e^{-s\tau} \int\_{-\tau}^{0} e^{-s\tau} \, d\mathbf{x}\_{\prime} \tag{17}
$$

where *<sup>E</sup>*(*s*) is the Laplace transform of *<sup>e</sup>*(*t*), *<sup>e</sup>*(0) = *Fz*(0) <sup>−</sup> *Gx*(0) <sup>−</sup> *Hy*(0), <sup>Δ</sup>(*s*) = *<sup>s</sup><sup>a</sup> <sup>I</sup>* <sup>−</sup> *<sup>K</sup>* <sup>−</sup> *<sup>e</sup>*−*s<sup>τ</sup> <sup>I</sup>* is the characteristic matrix. Then,

$$\det(\Delta(\mathbf{s})) = \left| \mathbf{s}^a I - \mathbf{K} - \mathbf{e}^{-s\tau} I \right| = (\mathbf{s}^a - k\_1 - \mathbf{e}^{-s\tau})(\mathbf{s}^a - k\_2 - \mathbf{e}^{-s\tau}) \dots (\mathbf{s}^a - k\_{\text{il}} - \mathbf{e}^{-s\tau}) = 0. \tag{18}$$

Assume

$$(s^a - k\_i - \varepsilon^{-s\tau}) = 0, i = 1, 2, \dots, n. \tag{19}$$

has a root *s* = *wi* = |*w*|(cos(*π*/2) + *i* sin(±*π*/2)). Thus,

$$|w|^a \left(\cos(a\pi/2) + i\sin(\pm a\pi/2)\right) - k\_i - \cos(\omega \tau) + i\sin(\omega \tau) = 0. \tag{20}$$

From the above equation, we can get

$$\begin{aligned} |w|^a \cos(a\pi/2) - k\_i &= \cos(\omega \tau), \\ |w|^a \sin(\pm a\pi/2) &= -\sin(\omega \tau). \end{aligned} \tag{21}$$

Hence,

$$\left|w\right|^{2a} - 2k\_i \cos(a\pi/2) \left|w\right|^a + k\_i^2 - 1 = 0. \tag{22}$$

Since *ki* < −1/ sin(*απ*/2), *α* ∈ (0, 1), the discriminant of the roots satisfies

$$\begin{split} \Delta &= \left(-2k\_{\bar{i}}\cos(a\pi/2)\right)^{2} - 4(k\_{\bar{i}}^{2} - 1) \\ &= 4(1 - k\_{\bar{i}}^{2}\sin^{2}(a\pi/2)) \\ &< 0. \end{split} \tag{23}$$

Then, Equation (22) has no real solutions, and Equation (18) has no purely imaginary roots.

In light of Corollary 1, the zero solution of the system in Equation (16) is globally asymptotically stable, i.e., multi-switching combination synchronization is obtained between the systems in Equations (9) and (10) and the system in Equation (11).

#### **4. Numerical Examples**

Numerical simulations were carried out to illustrate the above proposed multi-switching combination synchronization scheme. We used the same systems as in [32] with time-delays, which are fractional-order delayed Lorenz, Chen, and Rössler systems, and the numerical simulations were carried out in MATLAB.

The fractional-order delayed Lorenz system [37] was considered as the first drive system

$$\begin{cases} D^{\alpha} \mathbf{x}\_{1} = a\_{1} (\mathbf{x}\_{2} - \mathbf{x}\_{1}), \\ D^{\alpha} \mathbf{x}\_{2} = c\_{1} \mathbf{x}\_{1} - \mathbf{x}\_{2} - \mathbf{x}\_{1} \mathbf{x}\_{3}, \\ D^{\alpha} \mathbf{x}\_{3} = \mathbf{x}\_{1} \mathbf{x}\_{2} - b\_{1} \mathbf{x}\_{3} (t - \tau). \end{cases} \tag{24}$$

The system in Equation (24) exhibits a chaotic attractor, as illustrated in Figure 1. The system in Equation (24) can be rewritten as

$$\begin{aligned} D^{\mathfrak{a}}\mathbf{x}(t) &= \mathbf{x}(t) + \mathbf{x}(t-\tau) + A(\mathbf{x}(t), \mathbf{x}(t-\tau)), \\ \mathbf{x}(t) &= \mathbf{x}(0), t \in [-\tau, 0], \end{aligned} \tag{25}$$

where

$$A(\mathbf{x}(t), \mathbf{x}(t-\tau)) = \begin{pmatrix} a\_1 \mathbf{x}\_2 - (a\_1 + 1)\mathbf{x}\_1 - \mathbf{x}\_1(t-\tau) \\ c\_1 \mathbf{x}\_1 - 2\mathbf{x}\_2 - \mathbf{x}\_1 \mathbf{x}\_3 - \mathbf{x}\_2(t-\tau) \\ \mathbf{x}\_1 \mathbf{x}\_2 - (b\_1 + 1)\mathbf{x}\_3(t-\tau) - \mathbf{x}\_3 \end{pmatrix}. \tag{26}$$

**Figure 1.** Chaotic attractor of Lorenz system with *α* = 0.95, *τ* = 0.4: (**a**) *x*<sup>3</sup> − *x*<sup>1</sup> plane; (**b**) *x*<sup>3</sup> − *x*<sup>2</sup> plane; (**c**) *x*<sup>2</sup> − *x*<sup>1</sup> plane; and (**d**) *x*<sup>3</sup> − *x*<sup>1</sup> − *x*<sup>2</sup> space.

The fractional-order delayed Chen system [14] is the second drive system

$$\begin{cases} D^\mu y\_1 = a\_2 (y\_2 - y\_1)\_\prime \\ D^\mu y\_2 = (c\_2 - a\_2) y\_1 - y\_1 y\_3 + c\_2 y\_2 \\ D^\mu y\_3 = y\_1 y\_2 - b\_2 y\_3 (t - \tau) .\end{cases} \tag{27}$$

*Appl. Sci.* **2019**, *9*, 4348

The system in Equation (27) displays a chaotic attractor, as shown in Figure 2. We rewrite the system in Equation (27) as

$$\begin{aligned} D^a y(t) &= y(t) + y(t - \tau) + B(y(t), y(t - \tau)), \\ y(t) &= y(0), t \in [-\tau, 0], \end{aligned} \tag{28}$$

where

$$B(y(t), y(t-\tau)) = \begin{pmatrix} a\_2(y\_2 - y\_1) - y\_1 - y\_1(t-\tau) \\ (c\_2 - a\_2)y\_1 - y\_1y\_3 + (c\_2 - 1)y\_2 - y\_2(t-\tau) \\ y\_1y\_2 - (b\_2 + 1)y\_3(t-\tau) - y\_3 \end{pmatrix} . \tag{29}$$

**Figure 2.** Chaotic attractor of Chen system with *α* = 0.95, *τ* = 0.4: (**a**) *y*<sup>3</sup> − *y*<sup>1</sup> plane; (**b**) *y*<sup>3</sup> − *y*<sup>2</sup> plane; (**c**) *y*<sup>2</sup> − *y*<sup>1</sup> plane; and (**d**) *y*<sup>3</sup> − *y*<sup>1</sup> − *y*<sup>2</sup> space.

The fractional-order delayed Rössler system is the response system, given by

$$\begin{cases} D^{\mathfrak{a}}z\_1 = -(z\_2 + z\_3) + 0.2z\_1(t - \tau) + \mathcal{U}\_1, \\ D^{\mathfrak{a}}z\_2 = z\_1 + a\_3 z\_2 + \mathcal{U}\_2, \\ D^{\mathfrak{a}}z\_3 = z\_3(z\_1 - m\_3) + b\_3 + \mathcal{U}\_3. \end{cases} \tag{30}$$

where *U*1, *U*<sup>2</sup> and *U*<sup>3</sup> are determined afterwards. Without the controllers, the system in Equation (30) exhibits a chaotic attractor, as illustrated in Figure 3. The system in Equation (30) is rewritten as

$$\begin{aligned} D^a z(t) &= z(t) + z(t - \tau) + \mathbb{C}(z(t), z(t - \tau)) + \mathcal{U}, \\ z(t) &= z(0), t \in [-\tau, 0], \end{aligned} \tag{31}$$

where

$$\mathcal{C}(z(t), z(t-\tau)) = \begin{pmatrix} -(z\_2 + z\_3) - 0.8z\_1(t-\tau) - z\_1 \\ z\_1 + (a\_3 - 1)z\_2 - z\_2(t-\tau) \\ z\_3(z\_1 - m\_3) + b\_3 - z\_3 - z\_3(t-\tau) \end{pmatrix}.\tag{32}$$

For the systems in Equations (24), (27) and (30), there are eight possible switch combination synchronization cases.

When *k* = *l* = *m*, we have *e*123, *e*231, *e*<sup>312</sup> and *e*132, *e*213, *e*321. When *k* = *m* = *l*, we have *e*121, *e*232, *e*<sup>313</sup> and *e*131, *e*212, *e*323. When *k* = *l* = *m*, we have *e*122, *e*233, *e*<sup>311</sup> and *e*133, *e*211, *e*322. When *k* = *l* = *m*, we have *e*112, *e*223, *e*<sup>331</sup> and *e*113, *e*221, *e*332.

**Figure 3.** Chaotic attractor of Rössler system with *α* = 0.95, *τ* = 0.4: (**a**) *z*<sup>1</sup> − *z*<sup>3</sup> plane; (**b**) *z*<sup>2</sup> − *z*<sup>3</sup> plane; (**c**) *z*<sup>1</sup> − *z*<sup>2</sup> plane; and (**d**) *z*<sup>2</sup> − *z*<sup>1</sup> − *z*<sup>3</sup> space.

We randomly pick two cases

$$\begin{cases} \mathbf{e}\_{123} = f\_1 \mathbf{z}\_1 - g\_2 \mathbf{x}\_2 - h\_3 y\_3, \\ \mathbf{e}\_{231} = f\_2 \mathbf{z}\_2 - g\_3 \mathbf{x}\_3 - h\_1 y\_1, \quad \mathbf{case 1} \\ \mathbf{e}\_{312} = f\_3 \mathbf{z}\_3 - g\_1 \mathbf{x}\_1 - h\_2 y\_2. \end{cases} \tag{33}$$

and

$$\begin{cases} \mathbf{e}\_{12} = f\_1 \mathbf{z}\_1 - g\_1 \mathbf{x}\_1 - h\_2 \mathbf{y}\_2, \\ \mathbf{e}\_{223} = f\_2 \mathbf{z}\_2 - g\_2 \mathbf{x}\_2 - h\_3 \mathbf{y}\_3, \quad \mathbf{case 2} \\ \mathbf{e}\_{331} = f\_3 \mathbf{z}\_3 - g\_3 \mathbf{x}\_3 - h\_1 \mathbf{y}\_1. \end{cases} \tag{34}$$

In the following, we analyze these two cases in detail.

#### **Case 1**

From the systems in Equations (24), (27) and (30), we have the error dynamical system

$$\begin{cases} D^{\mathfrak{a}} \mathfrak{e}\_{123} = f\_1 D^{\mathfrak{a}} z\_1 - \mathfrak{g}\_2 D^{\mathfrak{a}} x\_2 - h\_3 D^{\mathfrak{a}} y\_3 \\ D^{\mathfrak{a}} \mathfrak{e}\_{231} = f\_2 D^{\mathfrak{a}} z\_2 - \mathfrak{g}\_3 D^{\mathfrak{a}} x\_3 - h\_1 D^{\mathfrak{a}} y\_1 \\ D^{\mathfrak{a}} \mathfrak{e}\_{312} = f\_3 D^{\mathfrak{a}} z\_3 - \mathfrak{g}\_1 D^{\mathfrak{a}} x\_1 - h\_2 D^{\mathfrak{a}} y\_2 \end{cases} \tag{35}$$

such that

$$\begin{cases} \lim\_{t \to +\infty} \parallel f\_1 z\_1 - g\_2 x\_2 - h\_3 y\_3 \parallel = 0, \\\lim\_{t \to +\infty} \parallel f\_2 z\_2 - g\_3 x\_3 - h\_1 y\_1 \parallel = 0, \\\lim\_{t \to +\infty} \parallel f\_3 z\_3 - g\_1 x\_1 - h\_2 y\_2 \parallel = 0. \end{cases} \tag{36}$$

Substituting the systems in Equations (24), (27) and (30) into the system in Equation (35) yields

$$\begin{cases} D^{a}e\_{123} = f\_1(-(z\_2+z\_3) + 0.2z\_1(t-\tau) + \mathcal{U}\_1) - g\_2(c\_1\mathbf{x}\_1 - \mathbf{x}\_2 - \mathbf{x}\_1\mathbf{x}\_3) \\ \qquad \quad - h\_3(y\_1y\_2 - b\_2y\_3(t-\tau)), \\ D^{a}e\_{231} = f\_2(z\_1 + a\_3z\_2 + \mathcal{U}\_2) - g\_3(\mathbf{x}\_1\mathbf{x}\_2 - b\_1\mathbf{x}\_3(t-\tau)) \\ \qquad \quad - h\_1(a\_2(y\_2 - y\_1)), \\ D^{a}e\_{312} = f\_3(z\_3(z\_1 - m\_3) + b\_3 + \mathcal{U}\_3) - g\_1(a\_1(\mathbf{x}\_2 - \mathbf{x}\_1)) \\ \qquad \quad - h\_2((c\_2 - a\_2)y\_1 - y\_1y\_3 + c\_2y\_2). \end{cases} \tag{37}$$

Here, we obtain the following results.

**Theorem 2.** *Multi-switching combination synchronization between the systems in Equations (24) and (27) and the system in Equation (30) can be achieved with the following controllers*

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *<sup>U</sup>*<sup>1</sup> <sup>=</sup> <sup>1</sup> *f*1 {(*k*<sup>1</sup> − 1)(*f*1*z*<sup>1</sup> − *g*2*x*<sup>2</sup> − *h*3*y*3) + *g*2(*c*1*x*<sup>1</sup> − 2*x*<sup>2</sup> − *x*1*x*<sup>3</sup> − *x*2(*t* − *τ*)) + *h*3(*y*1*y*<sup>2</sup> − (*b*<sup>2</sup> + 1)*y*3(*t* − *τ*) − *y*3) − *f*1(−(*z*<sup>2</sup> + *z*3) − 0.8*z*1(*t* − *τ*) − *z*1)}, *<sup>U</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup> *f*2 {(*k*<sup>2</sup> − 1)(*f*2*z*<sup>2</sup> − *g*3*x*<sup>3</sup> − *h*1*y*1) + *g*3(*x*1*x*<sup>2</sup> − (*b*<sup>1</sup> + 1)*x*3(*t* − *τ*) − *x*3) + *h*1(*a*2(*y*<sup>2</sup> − *y*1) − *y*<sup>1</sup> − *y*1(*t* − *τ*)) − *f*2(*z*<sup>1</sup> + (*a*<sup>3</sup> − 1)*z*<sup>2</sup> − *z*2(*t* − *τ*))}, *<sup>U</sup>*<sup>3</sup> <sup>=</sup> <sup>1</sup> *f*3 {(*k*<sup>3</sup> − 1)(*f*3*z*<sup>3</sup> − *g*1*x*<sup>1</sup> − *h*2*y*2) + *g*1(*a*1*x*<sup>2</sup> − (*a*<sup>1</sup> + 1)*x*<sup>1</sup> − *x*1(*t* − *τ*)) + *h*2((*c*<sup>2</sup> − *a*2)*y*<sup>1</sup> − *y*1*y*<sup>3</sup> + (*c*<sup>2</sup> − 1)*y*<sup>2</sup> − *y*2(*t* − *τ*)) − *f*3(*z*3(*z*<sup>1</sup> − *m*3) + *b*<sup>3</sup> − *z*<sup>3</sup> − *z*3(*t* − *τ*))}. (38)

Supposing *F* = 0 and *G* = 0 or *H* = 0, we have the following results.

**Corollary 2.** *(i) Supposing that fi* = 0*, gi* = 0 *and hi* = 0 (*i* = 1, 2, 3)*, multi-switching hybrid projective synchronization between the systems in Equations (27) and (30) can be achieved with the following controllers*

$$\begin{cases} \mathcal{U}\_{1} = \frac{1}{f\_{1}} \{ (k\_{1} - 1)(f\_{1}z\_{1} - hy\_{3}) + h\_{3}(y\_{1}y\_{2} - (b\_{2} + 1)y\_{3}(t - \tau) - y\_{3}) \\ \qquad \quad - f\_{1}(-(z\_{2} + z\_{3}) - 0.8z\_{1}(t - \tau) - z\_{1}) \}, \\ \mathcal{U}\_{2} = \frac{1}{f\_{2}} \{ (k\_{2} - 1)(f\_{2}z\_{2} - h\_{1}y\_{1}) + h\_{1}(a\_{2}(y\_{2} - y\_{1}) - y\_{1} - y\_{1}(t - \tau)) \\ \qquad \quad - f\_{2}(z\_{1} + (a\_{3} - 1)z\_{2} - z\_{2}(t - \tau)) \}, \\ \mathcal{U}\_{3} = \frac{1}{f\_{3}} \{ (k\_{3} - 1)(f\_{3}z\_{3} - h\_{2}y\_{2}) + h\_{2}((c\_{2} - a\_{2})y\_{1} - y\_{1}y\_{3} + (c\_{2} - 1)y\_{2} - y\_{2}(t - \tau)) \\ \qquad \quad - f\_{3}(z\_{3}(z\_{1} - m\_{3}) + b\_{3} - z\_{3} - z\_{3}(t - \tau)) \}. \end{cases} \tag{39}$$

*(ii) Similarly, supposing that fi* = 0*, gi* = 0 *and hi* = 0 (*i* = 1, 2, 3)*, multi-switching hybrid projective synchronization between the systems in Equations (24) and (30) can be achieved with the following controllers*

$$\begin{cases} l\mathcal{U}\_{1} = \frac{1}{f\_{1}} \{ (k\_{1} - 1)(f\_{1}z\_{1} - g\_{2}x\_{2}) + \mathcal{g}\_{2}(c\_{1}x\_{1} - 2x\_{2} - x\_{1}x\_{3} - x\_{2}(t - \tau)) \\ \qquad \quad - f\_{1}(-(z\_{2} + z\_{3}) - 0.8z\_{1}(t - \tau) - z\_{1}) \}, \\ l\mathcal{U}\_{2} = \frac{1}{f\_{2}} \{ (k\_{2} - 1)(f\_{2}z\_{2} - g\_{3}x\_{3}) + \mathcal{g}\_{3}(x\_{1}x\_{2} - (b\_{1} + 1)x\_{3}(t - \tau) - x\_{3}) \\ \qquad \quad - f\_{2}(z\_{1} + (a\_{3} - 1)z\_{2} - z\_{2}(t - \tau)) \}, \\ l\mathcal{U}\_{3} = \frac{1}{f\_{3}} \{ (k\_{3} - 1)(f\_{3}z\_{3} - g\_{1}x\_{1}) + g\_{1}(a\_{1}x\_{2} - (a\_{1} + 1)x\_{1} - x\_{1}(t - \tau)) \\ \qquad \quad - f\_{3}(z\_{3}(z\_{1} - m\_{3}) + b\_{3} - z\_{3} - z\_{3}(t - \tau)) \}. \end{cases} \tag{40}$$

**Corollary 3.** *Supposing that fi* = 0*, gi* = 0 *and hi* = 0 (*i* = 1, 2, 3)*, the system in Equation (30) can be stabilized to its equilibrium O*(0, 0, 0) *with the following controllers*

$$\begin{cases} l\mathcal{U}\_1 = \frac{1}{f\_1} \{ (k\_1 - 1)(f\_1 z\_1) - f\_1(-(z\_2 + z\_3) - 0.8z\_1(t - \tau) - z\_1) \}, \\ \mathcal{U}\_2 = \frac{1}{f\_2} \{ (k\_2 - 1)(f\_2 z\_2) - f\_2(z\_1 + (a\_3 - 1)z\_2 - z\_2(t - \tau)) \}, \\ \mathcal{U}\_3 = \frac{1}{f\_3} \{ (k\_3 - 1)(f\_3 z\_3) - f\_3(z\_3(z\_1 - m\_3) + b\_3 - z\_3 - z\_3(t - \tau)) \}. \end{cases} \tag{41}$$

#### **Case 2**

From the systems in Equations (24), (27) and (30), we have

$$\begin{cases} D^{\mathfrak{a}}c\_{112} = f\_1 D^{\mathfrak{a}} z\_1 - \mathfrak{g}\_1 D^{\mathfrak{a}} x\_1 - h\_2 D^{\mathfrak{a}} y\_2, \\ D^{\mathfrak{a}} c\_{223} = f\_2 D^{\mathfrak{a}} z\_2 - \mathfrak{g}\_2 D^{\mathfrak{a}} x\_2 - h\_3 D^{\mathfrak{a}} y\_3, \\ D^{\mathfrak{a}} c\_{331} = f\_3 D^{\mathfrak{a}} z\_3 - \mathfrak{g}\_3 D^{\mathfrak{a}} x\_3 - h\_1 D^{\mathfrak{a}} y\_1. \end{cases} \tag{42}$$

such that

$$\begin{cases} \lim\_{t \to +\infty} \parallel f\_1 z\_1 - g\_1 \mathbf{x}\_1 - h\_2 y\_2 \parallel = 0, \\\lim\_{t \to +\infty} \parallel f\_2 z\_2 - g\_2 \mathbf{x}\_2 - h\_3 y\_3 \parallel = 0, \\\lim\_{t \to +\infty} \parallel f\_3 z\_3 - g\_3 \mathbf{x}\_3 - h\_1 y\_1 \parallel = 0. \end{cases} \tag{43}$$

Substituting the systems in Equations (24), (27), and (30) into the system in Equation (42) yields:

$$\begin{cases} D^{\mathfrak{a}}e\_{112} = f\_1(-(z\_2+z\_3) + 0.2z\_1(t-\tau) + \mathcal{U}\_1) - g\_1(a\_1(\mathbf{x}\_2-\mathbf{x}\_1)) \\ \qquad \quad - h\_2((c\_2-a\_2)y\_1 - y\_1y\_3 + c\_2y\_2), \\ D^{\mathfrak{a}}e\_{223} = f\_2(z\_1+a\_3z\_2 + \mathcal{U}\_2) - g\_2(c\_1\mathbf{x}\_1 - \mathbf{x}\_2 - \mathbf{x}\_1\mathbf{x}\_3) \\ \qquad \quad - h\_3(y\_1y\_2 - b\_2y\_3(t-\tau)), \\ D^{\mathfrak{a}}e\_{33} = f\_3(z\_3(z\_1-m\_3) + b\_3 + \mathcal{U}\_3) - g\_3(\mathbf{x}\_1\mathbf{x}\_2 - b\_1\mathbf{x}\_3(t-\tau)) \\ \qquad \quad - h\_1(a\_2(y\_2-y\_1)). \end{cases} \tag{44}$$

Here, we have the following similar results.

**Theorem 3.** *Multi-switching combination synchronization between the systems in Equations (24) and (27) and the system in Equation (30) can be achieved with the following controllers*

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *<sup>U</sup>*<sup>1</sup> <sup>=</sup> <sup>1</sup> *f*1 {(*k*<sup>1</sup> − 1)(*f*1*z*<sup>1</sup> − *g*1*x*<sup>1</sup> − *h*2*y*2) + *g*1(*a*1*x*<sup>2</sup> − (*a*<sup>1</sup> + 1)*x*<sup>1</sup> − *x*1(*t* − *τ*)) + *h*2((*c*<sup>2</sup> − *a*2)*y*<sup>1</sup> − *y*1*y*<sup>3</sup> + (*c*<sup>2</sup> − 1)*y*<sup>2</sup> − *y*2(*t* − *τ*)) − *f*1(−(*z*<sup>2</sup> + *z*3) − 0.8*z*1(*t* − *τ*) − *z*1)}, *<sup>U</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup> *f*2 {(*k*<sup>2</sup> − 1)(*f*2*z*<sup>2</sup> − *g*2*x*<sup>2</sup> − *h*3*y*3) + *g*2(*c*1*x*<sup>1</sup> − 2*x*<sup>2</sup> − *x*1*x*<sup>3</sup> − *x*2(*t* − *τ*)) + *h*3(*y*1*y*<sup>2</sup> − (*b*<sup>2</sup> + 1)*y*3(*t* − *τ*) − *y*3) − *f*2(*z*<sup>1</sup> + (*a*<sup>3</sup> − 1)*z*<sup>2</sup> − *z*2(*t* − *τ*))}, *<sup>U</sup>*<sup>3</sup> <sup>=</sup> <sup>1</sup> *f*3 {(*k*<sup>3</sup> − 1)(*f*3*z*<sup>3</sup> − *g*3*x*<sup>3</sup> − *h*1*y*1) + *g*3(*x*1*x*<sup>2</sup> − (*b*<sup>1</sup> + 1)*x*3(*t* − *τ*) − *x*3) + *h*1(*a*2(*y*<sup>2</sup> − *y*1) − *y*<sup>1</sup> − *y*1(*t* − *τ*)) − *f*3(*z*3(*z*<sup>1</sup> − *m*3) + *b*<sup>3</sup> − *z*<sup>3</sup> − *z*3(*t* − *τ*))}. (45)

**Corollary 4.** *(i) Supposing that fi* = 0*, gi* = 0 *and hi* = 0 (*i* = 1, 2, 3)*, multi-switching hybrid projective synchronization between the systems in Equations (27) and (30) can be achieved with the following controllers*

$$\begin{cases} \mathcal{U}\_{1} = \frac{1}{f\_{1}} \{ (k\_{1} - 1)(f\_{1}z\_{1} - h\_{2}y\_{2}) + h\_{2}((c\_{2} - a\_{2})y\_{1} - y\_{1}y\_{3} + (c\_{2} - 1)y\_{2} - y\_{2}(t - \tau)) \} \\ \qquad \quad - f\_{1}(-(z\_{2} + z\_{3}) - 0.8z\_{1}(t - \tau) - z\_{1}) \}, \\ \mathcal{U}\_{2} = \frac{1}{f\_{2}} \{ (k\_{2} - 1)(f\_{2}z\_{2} - h\_{3}y\_{3}) + h\_{3}(y\_{1}y\_{2} - (b\_{2} + 1)y\_{3}(t - \tau) - y\_{3}) \} \\ \qquad \quad - f\_{2}(z\_{1} + (a\_{3} - 1)z\_{2} - z\_{2}(t - \tau)) \}, \\ \mathcal{U}\_{3} = \frac{1}{f\_{3}} \{ (k\_{3} - 1)(f\_{3}z\_{3} - h\_{1}y\_{1}) + h\_{1}(a\_{2}(y\_{2} - y\_{1}) - y\_{1} - y\_{1}(t - \tau)) \} \\ \qquad \quad - f\_{3}(z\_{3}(z\_{1} - m\_{3}) + b\_{3} - z\_{3} - z\_{3}(t - \tau)) \}. \end{cases} \tag{46}$$

*(ii) Supposing that fi* = 0*, gi* = 0 *and hi* = 0 (*i* = 1, 2, 3)*, hybrid projective synchronization between the systems in Equations (24) and (30) can be achieved with the following controllers*

$$\begin{cases} l\mathcal{U}\_{1} = \frac{1}{f\_{1}} \{ (k\_{1} - 1)(f\_{1}z\_{1} - g\_{1}\mathbf{x}\_{1}) + \mathbf{g}\_{1}(a\_{1}\mathbf{x}\_{2} - (a\_{1} + 1)\mathbf{x}\_{1} - \mathbf{x}\_{1}(t - \tau)) \\ \qquad \quad - f\_{1}(-(z\_{2} + z\_{3}) - 0.8z\_{1}(t - \tau) - z\_{1}) \}, \\ l\mathcal{U}\_{2} = \frac{1}{f\_{2}} \{ (k\_{2} - 1)(f\_{2}z\_{2} - g\_{2}\mathbf{x}\_{2}) + g\_{2}(c\_{1}\mathbf{x}\_{1} - 2\mathbf{x}\_{2} - \mathbf{x}\_{1}\mathbf{x}\_{3} - \mathbf{x}\_{2}(t - \tau)) \\ \qquad \quad - f\_{2}(z\_{1} + (a\_{3} - 1)z\_{2} - z\_{2}(t - \tau)) \}, \\ l\mathcal{U}\_{3} = \frac{1}{f\_{3}} \{ (k\_{3} - 1)(f\_{3}z\_{3} - g\_{3}\mathbf{x}\_{3}) + g\_{3}(\mathbf{x}\_{1}\mathbf{x}\_{2} - (b\_{1} + 1)\mathbf{x}\_{3}(t - \tau) - \mathbf{x}\_{3}) \\ \qquad \quad - f\_{3}(z\_{3}(z\_{1} - m\_{3}) + b\_{3} - z\_{3} - z\_{3}(t - \tau)) \}. \end{cases} \tag{47}$$

**Corollary 5.** *Supposing that fi* = 0*, gi* = 0 *and hi* = 0 (*i* = 1, 2, 3)*, the system in Equation (30) can be stabilized to its equilibrium O*(0, 0, 0) *with the following controllers*

$$\begin{cases} l\mathcal{U}\_1 = \frac{1}{f\_1} \{ (k\_1 - 1)(f\_1 z\_1) - f\_1(-(z\_2 + z\_3) - 0.8z\_1(t - \tau) - z\_1) \}, \\ \mathcal{U}\_2 = \frac{1}{f\_2} \{ (k\_2 - 1)(f\_2 z\_2) - f\_2(z\_1 + (a\_3 - 1)z\_2 - z\_2(t - \tau)) \}, \\ \mathcal{U}\_3 = \frac{1}{f\_3} \{ (k\_3 - 1)(f\_3 z\_3) - f\_3(z\_3(z\_1 - m\_3) + b\_3 - z\_3 - z\_3(t - \tau)) \}. \end{cases} \tag{48}$$

The system parameters are given as *a*<sup>1</sup> = 10, *b*<sup>1</sup> = <sup>8</sup> <sup>3</sup> , *c*<sup>1</sup> = 28, *a*<sup>2</sup> = 35, *b*<sup>2</sup> = 3, *c*<sup>2</sup> = 28, *a*<sup>3</sup> = 0.4, *b*<sup>3</sup> = 0.2, *m*<sup>3</sup> = 10, thus the systems in Equations (24), (27) and (30) exhibit chaotic behaviors, respectively. We assume *f*<sup>1</sup> = *f*<sup>2</sup> = *f*<sup>3</sup> = 1, *g*<sup>1</sup> = *g*<sup>2</sup> = *g*<sup>3</sup> = 1 and *h*<sup>1</sup> = *h*<sup>2</sup> = *h*<sup>3</sup> = 1, and the initial values are (*x*1(0), *x*2(0), *x*3(0)) = (−20, 2, 3),(*y*1(0), *y*2(0), *y*3(0)) = (7, 4.04, 20) and (*z*1(0), *z*2(0), *z*3(0)) = (1, 2, 40), respectively. Multi-switching combination synchronization between the systems in Equations (24), (27) and (30) can be realized with *K* = *diag*{−10, −10, −10}. Figure 4 illustrates synchronization errors *e*123, *e*231, *e*312. Figure 5 shows synchronization states *x*<sup>2</sup> + *y*<sup>3</sup> vs. *z*1, *x*<sup>3</sup> + *y*<sup>1</sup> vs. *z*<sup>2</sup> and *x*<sup>1</sup> + *y*<sup>2</sup> vs. *z*<sup>3</sup> of the drive systems in Equations (24) and (27) and the response system in Equation (30). Figure 6 displays synchronization errors *e*112, *e*223, and *e*331. Figure 7 illustrates synchronization states *x*<sup>1</sup> + *y*<sup>2</sup> vs. *z*1, *x*<sup>2</sup> + *y*<sup>3</sup> vs. *z*<sup>2</sup> and *x*<sup>3</sup> + *y*<sup>1</sup> vs. *z*<sup>3</sup> of the drive systems in Equations (24) and (27) and the response system in Equation (30). In Figures 4–7, we can see that the multi-switching combination synchronization errors converge to zero, i.e., the multi-switching combination synchronizations for Cases 1 and 2 are achieved, respectively.

The feedback gain matrix *K* is an important factor to affect the convergence of the error systems. With the increase of the absolute value of *ki*, the convergence time will be shortened. Thus, we carried out one more simulation with *K* = *diag*{−40, −40, −40}. Figures 8 and 9 illustrate the synchronization errors for Cases 1 and 2, respectively. By comparing Figures 4 and 8, as well as Figures 6 and 9, it is easy to see that convergence times are shortened obviously.

**Figure 4.** Synchronization errors among the systems in Equations (24), (27) and (30) with *k*<sup>1</sup> = *k*<sup>2</sup> = *k*<sup>3</sup> = −10: (**a**) *e*123; (**b**) *e*231; and (**c**) *e*312.

**Figure 5.** Responses for states between the systems in Equations (24) and (27) and the system in Equation (30): (**a**) *x*<sup>2</sup> + *y*<sup>3</sup> vs. *z*1; (**b**) *x*<sup>3</sup> + *y*<sup>1</sup> vs. *z*2; and (**c**) *x*<sup>1</sup> + *y*<sup>2</sup> vs. *z*3.

**Figure 6.** Synchronization errors among the systems in Equations (24), (27) and (30) with *k*<sup>1</sup> = *k*<sup>2</sup> = *k*<sup>3</sup> = −10: (**a**) *e*112; (**b**) *e*223; and (**c**) *e*331.

**Figure 7.** Responses for states between the systems in Equations (24) and (27) and the system in Equation (30): (**a**) *x*<sup>1</sup> + *y*<sup>2</sup> vs. *z*1; (**b**) *x*<sup>2</sup> + *y*<sup>3</sup> vs. *z*2; and (**c**) *x*<sup>3</sup> + *y*<sup>1</sup> vs. *z*3.

**Figure 8.** Synchronization errors among the systems in Equations (24), (27) and (30) with *k*<sup>1</sup> = *k*<sup>2</sup> = *k*<sup>3</sup> = −40: (**a**) *e*123; (**b**) *e*231; and (**c**) *e*312.

**Figure 9.** Synchronization errors among the systems in Equations (24), (27) and (30) with *k*<sup>1</sup> = *k*<sup>2</sup> = *k*<sup>3</sup> = −40: (**a**) *e*112; (**b**) *e*223; and (**c**) *e*331.

#### **5. Conclusions**

We extended previous work [32] to investigate multi-switching combination synchronization among three non-identical fractional-order delayed systems by introducing time-delays. Based on the stability theory for linear fractional-order systems with multiple time-delays, we designed appropriate controllers to obtain multi-switching combination synchronization among three non-identical fractional-order delayed systems. The simulations are in accordance with the theoretical analysis.

On the one hand, when applying multi-switching combination synchronization of fractional-order delayed chaotic systems in secure communications, fractional-order and time-delay can enrich systems' dynamics. On the other hand, the origin information can be separated into two parts and embedded different parts in separate drive systems via combination synchronization scheme. Besides, because the switched states are unpredictable, this synchronization scheme can increase the security of the transmitted information in secure communication. Thus, the communication security will be enhanced, which makes multi-switching combination synchronization of fractional-order delayed chaotic systems able to find better applications in security communication.

**Author Contributions:** Methodology, X.Z.; Writing—Original Draft Preparation, B.L.; Writing—Review and Editing, Y.W.; and Funding Acquisition, X.Z.

**Funding:** This work was supported by the Natural Science Foundations of China under Grant No. 61463050, the NSF of Yunnan Province under Grant No. 2015FB113.

**Acknowledgments:** The authors sincerely thank the anonymous reviewers for their constructive comments and insightful suggestions.

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

#### **References**


c 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*
