**Stability Analysis and Stability Enhancement Based on Virtual Harmonic Resistance for Meshed DC Distributed Power Systems with Constant Power Loads**

#### **Huiyong Hu, Xiaoming Wang, Yonggang Peng \*, Yanghong Xia, Miao Yu and Wei Wei**

College of Electrical Engineering, Zhejiang University, Hangzhou 310027, China; huhuiyong@zju.edu.cn (H.H.); wxm\_15@zju.edu.cn (X.W.); royxiayh@zju.edu.cn (Y.X.); zjuyumiao@zju.edu.cn (M.Y.); wwei@zju.edu.cn (W.W.) **\*** Correspondence: pengyg@zju.edu.cn; Tel.: +86-571-8795-2653

Academic Editor: Birgitte Bak-Jensen Received: 24 November 2016; Accepted: 2 January 2017; Published: 9 January 2017

**Abstract:** This paper addresses the stability issue of the meshed DC distributed power systems (DPS) with constant power loads (CPLs) and proposes a stability enhancement method based on virtual harmonic resistance. In previous researches, the network dynamics of the meshed DC DPS are often ignored, which affects the derivation of the equivalent system impendence. In addition, few of them have considered the meshed DC DPS including multiple sources with voltage-controlled converters and CPLs. To tackle the aforementioned challenge, this paper mainly makes the following efforts. The component connection method (CCM) is employed and expanded to derive the stability criterion of the meshed DC DPS with CPLs. This stability criterion can be simplified to relate only with the network node admittance matrix, the output impendences of the sources, and the input admittances of the CPLs. A virtual harmonic resistance through the second-order generalized integrator (SOGI) is added in the source with the voltage-controlled converter to lower the peak of the source output impendence, which can enhance the stability of the meshed DC DPS. The effectiveness of the proposed stability criterion and stability enhancement method are verified by nonlinear dynamic simulations.

**Keywords:** meshed DC distributed power systems; stability criterion; virtual harmonic resistance; stability enhancement

#### **1. Introduction**

With the increasing penetration of renewable energy generation into the modern electric girds, meshed DC distributed power systems (DPS) have been widely adopted [1–3]. The meshed DC DPS physically is composed of smaller power modules/subsystems. Usually, these smaller power modules/subsystems are designed only based on the stability requirement in its stand-alone operation, and thus it can operate well in the stand-alone application. However, the meshed DC DPS may become unstable due to the interaction among the modules/subsystems and the network [4]. Furthermore, the negative resistance characteristic of the constant power load (CPL) is an important unstable factor.

Many stability methods in previous researches are given to analyze the stability of the single-source-single-load system or the parallel-source-parallel-load system [4–6]. However, the network dynamics of the meshed DC DPS are often ignored, which affects the derivation of the equivalent system impendence. In addition, few of them have considered the meshed DC DPS including multiple sources with voltage-controlled converters and CPLs. A general approach for the stability analysis of the meshed DC DPS is to build its whole state-space model, and to identify the eigenvalues of the state matrix [7]. However, this method requires the detailed models of loads and network dynamics, and the formulation of the system matrices may be very complex. To overcome this

problem, the component connection method (CCM) is introduced for the stability analysis [8]. In the CCM, the dynamics of network and power modules/subsystems are modeled separately by a set of two vector-matrix equations, which is beneficial to reduce the computation burden of formulating the system transfer matrices [9]. For the ac power-electronics-based power system, the CCM is employed to derivate the closed-loop response transfer matrix of the overall power which can predict the instability of the power system [10], but it does not separate the source impedance and load admittance with the network node admittance matrix. Thus, it is difficult to obtain the general stability criterion intuitionally. In this paper, the CCM is also employed and expanded to derive the stability criterion for the meshed DC DPS with the CPLs. The stability criterion can be simplified to relate only with the network node admittance matrix, the output impendence of the source, and the input admittance of the CPLs which can be separately obtained by theoretical calculation or impedance analyzer measurement.

Many stability enhancement methods have been given in previous researches. An adaptive active capacitor converter (AACC) is introduced to effectively stabilize the cascaded system [11], but this method needs additional hardware devices. Another approach is to add a virtual adaptive parallel impedance [12] or adaptive series-virtual-impedance [13] in the input of the load converter to improve the stability in the cascaded DC/DC converter system. However, in the meshed DC DPS, it is not very effective due to the existence of an electromagnetic interference (EMI) filter for the CPL with a switching-mode power converter. The input admittance of the CPLs is mainly determined by the EMI filter, and the virtual impedance is difficult to be added by the load converter control part. Another way is to adjust the source output impedance. Modifying the source output impedance to be zero [14] and adding virtual impedances in the source output filters [15] are both feasible strategies by the source converter control part, but both of them will influence the low and high frequency characteristics of the source output impedance. To overcome the problem, this paper proposes an effective and achievable method, which is to add a virtual harmonic resistance through the second-order generalized integrator (SOGI) in the source converter control part to lower the peak of the source output impendence.

The main contributions of this paper are follows. Firstly, the stability criterion for the meshed DC DPS is derived, which has been simplified to relate only with the network node admittance matrix, the output impendence of the source, and the input admittance of the CPLs. It is easy to implement, as the impedance and the admittance can both be separately obtained by the theoretical calculation or impedance analyzer measurement. Secondly, an effective stability enhancement method by adding a virtual harmonic resistance is proposed, with which only the small range middle frequency characteristic of the source output impedance is modified, while the low and high frequency characteristics of the source output impedance are still unmodified. The results of the nonlinear dynamic simulations verify the effectiveness of the proposed stability criterion and stability enhancement method.

The remainder of this paper is organized as follows. The stability criterion for the meshed DC DPS is given in Section 2. The modeling of the source with voltage-controlled converters and the CPLs are derived, and the proposed stability enhancement method and the stability analysis for the meshed DC DPS are given in Section 3. Then, the nonlinear dynamic simulations are conducted in Section 4. At last, conclusions are drawn in Section 5.

#### **2. Stability Criterion for the Meshed DC DPS**

The general form of DC DPS is given in [4], but it ignores the network dynamics. In [4], any converter in a DC DPS can be classified as either a bus voltage controlled converter (BVCC) or a bus current controlled converter (BCCC). In this study, a four node meshed DC DPS as Figure 1 is considered as an example, which is composed by two sources with voltage-controlled converters and two CPLs. The sources with voltage-controlled converters are controlled with constant voltage, hence they can be treated as BVCCs. The CPLs are controlled to absorb constant power from the meshed DC DPS, hence they can be treated as BCCCs.

**Figure 1.** The block diagram of the component connection method (CCM) applied for the studied meshed DC distributed power systems (DPS).

Although all of the sources and CPLs are independently designed well, there will be complex interactions through the power cables in the meshed DC DPS. This consequently necessitates the use of CCM to analyze the stability of the meshed DC DPS.

The block diagram of the CCM applied for the studied meshed DC DPS is shown in Figure 1, where the CCM decomposes the overall DPS into four subsystems and the network. The two BVCCs are modeled by the Thevenin equivalent circuits [16], and the two BCCCs are modeled by the Norton equivalent circuits [17].

In Figure 1, *Zv*<sup>1</sup> and *Zv*<sup>2</sup> are the closed-loop output impedances of the sources, and *Yc*<sup>1</sup> and *Yc*<sup>2</sup> are the closed-loop input admittances of the CPLs. % *v*<sup>1</sup> *v*<sup>2</sup> *v*<sup>3</sup> *v*<sup>4</sup> &*T* is the node voltage column vector, and % *i*<sup>1</sup> *i*<sup>2</sup> *i*<sup>3</sup> *i*<sup>4</sup> &*T* is the node injection current column vector. % *icable*<sup>1</sup> *icable*<sup>2</sup> *icable*<sup>3</sup> *icable*<sup>4</sup> &*T* is the cable current column vector. *RL*<sup>3</sup> and *RL*<sup>4</sup> are the resistive loads connected to node 3 and node 4, respectively. % *Lcable*<sup>1</sup> *Lcable*<sup>2</sup> *Lcable*<sup>3</sup> *Lcable*<sup>4</sup> &*T* is the inductance column vector of cables, % *Rcable*<sup>1</sup> *Rcable*<sup>2</sup> *Rcable*<sup>3</sup> *Rcable*<sup>4</sup> &*T* is the parasitic resistance column vector of cables, and % *Ccable*<sup>1</sup> *Ccable*<sup>2</sup> *Ccable*<sup>3</sup> *Ccable*<sup>4</sup> &*T* is the capacitance column vector of cables. % *Gv*<sup>1</sup> *Gv*<sup>2</sup> *Gc*<sup>1</sup> *Gc*<sup>2</sup> &*T* is the voltage reference-to-output transfer function column vector of the source and CPLs. % *vv*<sup>1</sup>*ref vv*<sup>2</sup>*ref vc*<sup>1</sup>*ref vc*<sup>2</sup>*ref* &*<sup>T</sup>* is the reference voltage column vector of the source and CPLs, respectively.

Generally, the source is designed to be a voltage source that is stable when unloaded, and the CPL is designed to be stable when supplied by an ideal voltage source. That is, the unterminated behaviors of inverters % *Gv*<sup>1</sup> *Gv*<sup>2</sup> *Gc*<sup>1</sup> *Gc*<sup>2</sup> &*T* are stable, and % *Zv*<sup>1</sup> *Zv*<sup>2</sup> *Yc*<sup>1</sup> *Yc*<sup>2</sup> &*T* are stable, which means that there are no right-half plane poles.

To facilitate the formulation of the nodal admittance matrix, the Thevenin circuits of the BVCCs are converted to Norton circuits [10]. Then, the studied meshed DC DPS is represented as Figure 2.

**Figure 2.** The equivalent Norton circuit of the studied meshed DC DPS.

Applying the node voltage equation to the external dashed box part in Figure 2, it can be derived as

$$
\begin{bmatrix}
\begin{bmatrix}
Z\_{v1}^{-1}G\_{v1}\upsilon\_{v1ref} \\
Z\_{v2}^{-1}G\_{v2}\upsilon\_{v2ref} \\
G\_{\ell1}\upsilon\_{\ell1ref} \\
G\_{\ell2}\upsilon\_{\ell2ref}
\end{bmatrix} = 
\begin{pmatrix}
\begin{bmatrix}
Z\_{v1}^{-1} \\
Z\_{v2}^{-1} \\
\upsilon\_{\ell1} \\
\upsilon\_{\ell2}
\end{bmatrix} + \chi\_{\rm sys} \\
\begin{bmatrix}
\upsilon\_{\ell1} & \upsilon\_{\ell2}
\end{bmatrix}
\end{pmatrix} + \chi\_{\rm sys}
\begin{bmatrix}
\upsilon\_{1} \\
\upsilon\_{2} \\
\upsilon\_{3}
\end{bmatrix}
\end{bmatrix}
\tag{1}
$$

where *Ysys* is the network nodal admittance matrix as

$$\mathbf{Y}\_{sp} = \begin{bmatrix} \mathbb{C}\_{\text{idid}:0} + \mathbb{C}\_{\text{idid}:0} + \mathbf{Y}\_{\text{idid}} + \mathbf{Y}\_{\text{idid}:2} & \mathbf{0} & -\mathbf{Y}\_{\text{idid}:2} & -\mathbf{Y}\_{\text{idid}:1} \\ \mathbf{0} & \mathbb{C}\_{\text{idid}:0} + \mathbb{C}\_{\text{idid}:0} + \mathbf{Y}\_{\text{idid}} + \mathbf{Y}\_{\text{idid}} & -\mathbf{Y}\_{\text{idid}} & -\mathbf{Y}\_{\text{idid}:1} & \mathbf{0} \\ -\mathbf{Y}\_{\text{idid}:2} & -\mathbf{Y}\_{\text{idid}:0} & \mathbf{Y}\_{\text{id}} + \mathbb{C}\_{\text{idid}:0} + \mathbb{C}\_{\text{idid}:0} + \mathbf{Y}\_{\text{idid}:2} + \mathbf{Y}\_{\text{idid}:1} & \mathbf{0} & \mathbf{0} \\ -\mathbf{Y}\_{\text{idid}:1} & -\mathbf{Y}\_{\text{idid}:1} & \mathbf{0} & \mathbf{Y}\_{\text{idid}:1} + \mathbb{C}\_{\text{idid}:0} + \mathbf{Y}\_{\text{idid}:1} & \mathbf{0} & \mathbf{Y}\_{\text{idid}:1} \\ \end{bmatrix} \tag{2}$$

where *YL*<sup>3</sup> and *YL*<sup>4</sup> are the admittances for loads connected to node 3 and node 4, respectively. *Ycable*1~*Ycable*<sup>4</sup> are the admittances for cables and *Ycablei* = 1/(*Rcablei* + *Lcableis*), *i* = 1∼4 .

Then, the node voltage can be derived from (1) as follows:

$$
\begin{aligned}
\begin{bmatrix} v\_1 \\ v\_2 \\ v\_3 \\ v\_4 \end{bmatrix} &= \left( \begin{bmatrix} Z\_{v1}^{-1} \\ & Z\_{v2}^{-1} \\ & & Y\_{c1} \\ & & & Y\_{c2} \end{bmatrix} + Y\_{sys} \right)^{-1} \begin{bmatrix} Z\_{v1}^{-1} \\ & Z\_{v2}^{-1} \\ & & 1 \\ & & & 1 \end{bmatrix} \begin{bmatrix} G\_{v1}v\_{v1ref} \\ G\_{v2}v\_{v2ref} \\ G\_{v1}v\_{v1ref} \\ G\_{v2}v\_{v2ref} \end{bmatrix} \tag{3} \\\\ v &= \left( \begin{bmatrix} Z\_v^{-1} \\ & Y\_c \\ & & Y\_c \end{bmatrix} + Y\_{sys} \right)^{-1} \begin{bmatrix} Z\_v^{-1} \\ & I\_c \\ & & I\_c \end{bmatrix} \begin{bmatrix} G\_v \\ & G\_v \\ & G\_{vc} \end{bmatrix} \begin{bmatrix} v\_{vref} \\ v\_{vref} \\ v\_{cref} \\ v\_{cref} \end{bmatrix} \end{aligned} \tag{4}
$$

where *Zv* is the set of the closed-loop output impedances of the sources, and *Yc* is the set of the closed-loop input admittances of the CPLs. *v* is the set of node voltages, and *i* is the set of the node injection currents. *Gv* is the the set of voltage reference-to-output transfer functions of the source, *Gc* is the the set of voltage reference to output transfer functions of the CPLs, *vvre f* is the set of voltage references of the source, and *vcre f* is the set of voltage references of the CPLs. *Ic* is a unit diagonal

matrix whose dimension is equal to *Yc*. *Iv* is a unit diagonal matrix whose dimension is equal to *Zv*. And they satisfy the following relations:

$$\mathbf{Z}\_{v} = \begin{bmatrix} \mathbf{Z}\_{v1} \\ & \mathbf{Z}\_{v2} \end{bmatrix}, \mathbf{Y}\_{c} = \begin{bmatrix} \mathbf{Y}\_{c1} \\ & \mathbf{Y}\_{c2} \end{bmatrix}, \mathbf{G}\_{v} = \begin{bmatrix} \mathbf{G}\_{v1} \\ & \mathbf{G}\_{v2} \end{bmatrix}, \mathbf{G}\_{c} = \begin{bmatrix} \mathbf{G}\_{c1} \\ & \mathbf{G}\_{c2} \end{bmatrix}, \mathbf{v}\_{\text{var}f} = \begin{bmatrix} \mathbf{v}\_{c1\text{rf}f} \\ \mathbf{v}\_{c2\text{rf}f} \end{bmatrix}, \mathbf{v}\_{\text{var}f} = \begin{bmatrix} \mathbf{v}\_{c1\text{rf}f} \\ \mathbf{v}\_{c2\text{rf}f} \end{bmatrix} \tag{5}$$

$$\mathbf{v}\_{\text{var}} = \begin{bmatrix} \mathbf{v}\_{c1\text{rf}f} \\ \mathbf{v}\_{c2\text{rf}f} \end{bmatrix}$$

$$v = \begin{bmatrix} v\_1 & v\_2 & v\_3 & v\_4 \end{bmatrix}^T,\\ i = \begin{bmatrix} i\_1 & i\_2 & i\_3 & i\_4 \end{bmatrix}^T \tag{6}$$

Applying the node voltage equation to the internal dashed box part in Figure 2, *i* can be derived as:

$$\mathbf{M} = \mathbf{Y}\_{\text{sys}} \mathbf{v} = \mathbf{Y}\_{\text{sys}} \left( \begin{bmatrix} I\_{\text{v}} \\ & Y\_{\text{c}} \end{bmatrix} + \begin{bmatrix} Z\_{\text{v}} \\ & I\_{\text{c}} \end{bmatrix} \mathbf{Y}\_{\text{sys}} \right)^{-1} \begin{bmatrix} \mathbf{G}\_{\text{v}} \\ & \mathbf{G}\_{\text{c}} \end{bmatrix} \begin{bmatrix} v\_{\text{v}ref} \\ v\_{\text{c}ref} \end{bmatrix} \tag{7}$$

In the above equation, *Ysys* is the network nodal admittance matrix for a real physical system, which means *Ysys* is stable and has no right-half plane poles. *Gv* and *Gc* are both stable as in the above analysis. Therefore, the stability criterion of the meshed DC DPS *Tm* is

$$T\_m = \left( \left[ \begin{array}{cc} I\_v \\ & Y\_c \end{array} \right] + \left[ \begin{array}{cc} Z\_v \\ & I\_c \end{array} \right] Y\_{\text{sys}} \right)^{-1} \tag{8}$$

If the stability criterion *Tm* has right-half plane poles, the meshed DC DPS will be unstable. Through the analysis of the stability criterion *Tm*, it can be found that the stability of the meshed DC DPS is only related to the network node admittance matrix *Ysys*, the output impendence of the sources *Zv*, and the input admittance of the CPLs *Yc*, which can be separately obtained by theoretical calculation or impedance analyzer measurement.

#### **3. Modeling, Stability Enhancement Method and Stability Analysis**

#### *3.1. Modeling of the CPL*

The modeling of the CPL connected to node 3 is similar to the CPL connected to node 4. Without loss of generality, taking the CPL connected to node 4 as an example. The CPL usually exploits a switching-mode power converter as the interface with the node, and the compliance with EMI standards usually requires the insertion of an EMI filter between the switching-mode power converter and the node [18–20]. The close-loop circuit of the CPL with an EMI filter is shown as Figure 3.

**Figure 3.** The close-loop circuit of the constant power load (CPL) with an electromagnetic interference (EMI) filter.

Applying switch period average, the small-signal model of the CPL with EMI filter can be achieved as Figure 4.

**Figure 4.** The small-signal model of the CPL with an EMI filter.

The variables showing in capital letters represent the steady-state values of duty cycle, inductor current, capacitor voltage, and output current, and the variables with " ˆ " indicate the small-signal disturbance.

It can be derived from Figure 4 that the small-signal disturbance of the output current ˆ*i*<sup>4</sup> can be expressed as:

$$
\hat{\mathfrak{u}}\_4 = -\mathcal{Y}\_{\mathfrak{c2}}\mathfrak{d}\_4 + \mathcal{G}\_{\mathfrak{c2}}\mathfrak{d}\_{\mathfrak{c2}ref} \tag{9}
$$

$$Y\_{c2} = \frac{a\_4s^4 + a\_3s^3 + a\_2s^2 + a\_1s^1 + a\_0}{\Delta\_{c2}} \tag{10}$$

$$G\_{c2} = -\frac{c\_3s^3 + c\_2s^2 + c\_1s^1 + c\_0}{\Delta\_{c2}}\tag{11}$$

$$
\Delta\_{c2} = b\_5 s^5 + b\_4 s^4 + b\_3 s^3 + b\_2 s^2 + b\_1 s^1 + b\_0 \tag{12}
$$

where *Gc*<sup>2</sup>*<sup>v</sup>* is the voltage controller, *Yc*<sup>2</sup> is the closed-loop input admittance of the CPL, and *Gc*<sup>2</sup> is the voltage reference to output transfer function of the CPL. *a*<sup>0</sup> ∼*a*<sup>4</sup> , *b*<sup>0</sup> ∼*b*5, and *c*<sup>0</sup> ∼*c*<sup>3</sup> are given in Appendix A.

#### *3.2. Modeling of the Source with Voltage-Controlled Converters*

The modeling of the source with voltage-controlled converters connected to node 1 is similar to the source connected to node 2. Without loss of generality, taking the source connected to node 1 as an example, the close-loop circuit of the source with voltage-controlled converters is shown as Figure 5.

**Figure 5.** The close-loop circuit of the source with voltage-controlled converters (where *Gv*<sup>1</sup>*<sup>v</sup>* and *Gv*<sup>1</sup>*<sup>i</sup>* are the voltage and current controllers).

Figure 6 shows the small-signal block diagram of the multi-loop control system of the source with voltage-controlled converters.

**Figure 6.** The small-signal block diagram of the multi-loop control system of the source.

Thus, the dynamic behavior of the close-loop system can be given by

$$
\psi\_1 = -Z\_{v1}\hat{\imath}\_1 + G\_{v1}\vartheta\_{v1ref} \tag{13}
$$

where

$$Z\_{v1} = \frac{G\_{v1i}K\_{v1pwvu}V\_1 + L\_1s + r\_1}{L\_1\mathbb{C}\_1s^2 + (\mathbb{C}\_1\mathbb{G}\_{v1i}K\_{v1pwv}V\_1 + r\_1\mathbb{C}\_1)s + G\_{v1i}G\_{v1v}K\_{v1pwv}V\_1 + 1} \tag{14}$$

$$\mathbf{G}\_{v1} = \frac{\mathbf{G}\_{v1i}\mathbf{G}\_{v1v}\mathbf{K}\_{v1pwm}V\_1}{L\_1\mathbf{C}\_1\mathbf{s}^2 + (\mathbf{C}\_1\mathbf{G}\_{v1i}\mathbf{K}\_{v1pwm}V\_1 + r\_1\mathbf{C}\_1)\mathbf{s} + \mathbf{G}\_{v1i}\mathbf{G}\_{v1v}\mathbf{K}\_{v1pwm}V\_1 + 1} \tag{15}$$

where *Zv*<sup>1</sup> is the closed-loop output impedance of the source with voltage-controlled converters, and *Gv*<sup>1</sup> is the voltage reference-to-output transfer function of the source with voltage-controlled converters.

#### *3.3. Stability Enhancement Method*

The source output impedance has a peak at the resonant frequency. Based on the Nyquist criterion, instability between the source and the network may occur if the network input impedance becomes lower than the source output impedance at this frequency. Therefore, to ensure the system stability and to minimize the potential for inadvertent interactions with the source, it is important to lower the peak of the source output impedance. For this purpose, this paper proposes a method by adding a virtual harmonic resistance through the SOGI to be in parallel with the capacitor.

The transfer function *Gf*(*s*) of SOGI can be expressed as

$$G\_f(\mathbf{s}) = \frac{\mathbf{s}^2 + \omega^2}{\mathbf{s}^2 + k\omega\mathbf{s} + \omega^2} \tag{16}$$

where *ω* is the resonant frequency, and *k* is the frequency coefficient. The frequency-domain characteristic of *Gf*(*s*) is shown as Figure 7. Around the resonant frequency, the amplitude of *Gf*(*s*) is very small, and the amplitude of *Gf*(*s*) in other frequency ranges is 0 dB.

**Figure 7.** The frequency-domain characteristic of *Gf*(*s*).

Figure 8 shows the small-signal control block of the source 1. If a virtual harmonic resistance *Rvh* is required to be added in parallel with the capacitor, it is intuitively to introduce *Rvh* to the block as a dot-dashed part in Figure 8, which is the basic idea of this control method. However, this method cannot be realized by control directly. Therefore, *Rvh* is moved to the output of *Gv*<sup>1</sup>*v*(*s*) as the dashed part in Figure 8.

**Figure 8.** Adding a virtual harmonic resistance through the second-order generalized integrator (SOGI).

*GR*(*s*) is expressed as

$$G\_R(s) = \frac{L\_1s + r\_1 + G\_{v1i}K\_{v1puvu}V\_1}{G\_{v1i}K\_{v1puvu}V\_1R\_{vlu}}.\tag{17}$$

With the parallel virtual harmonic resistance *Rvh* = 1 Ω, the frequency-domain characteristic of the source output impedance is shown as Figure 9.

**Figure 9.** The source output impedance with/without parallel virtual harmonic resistance *Rvh*.

As can be seen, the peak of the source output impedance has been lower, and only the small range middle frequency characteristic of the source output impedance is modified, while the low and high frequency characteristics of the source output impedance are still unmodified.

#### *3.4. Stability Analysis*

All the main electrical parameters of the studied DC meshed DPS as shown in Figure 1 are listed in Appendix B.

The closed-loop output impedances of the sources *Zv* can be represented by using impedance calculated in Equations (5) and (14), the closed-loop input admittances of the CPLs *Yc* can be represented by using admittance calculated in Equations (5) and (10), and the network nodal admittance matrix *Ysys* can be represented by using admittance calculated in Equation (2). Then, the the stability of the meshed DC DPS can be analyzed with the stability criterion *Tm* expressed in Equation (8). If the stability criterion *Tm* has right-half plane poles, the meshed DC DPS will be unstable.

At first, it is necessry to compare the proposed stability criterion *Tm* in this paper and the previous stability criterion *Tm*<sup>2</sup> given in [4] as follows:

$$T\_{m2} = \frac{1}{1 + (Z\_{v1} // Z\_{v2})(Y\_{c1} + Y\_{c2} + Y\_{L3} + Y\_{L4})}.\tag{18}$$

The dominant poles of *Tm* ans *Tm*<sup>2</sup> are displayed as Figure 10, in which the cables are long *Rcablei* = 0.05 Ω , *Lcablei* = 0.5 mH, *Ccablei* = 100 μF, and *i* = 1∼4. As shown in Figure 10, all the dominant poles of *Tm*<sup>2</sup> are both in the left-half plane, which indicates that the meshed DC DPS is stable. However, there are two right-half plane poles 0.15 ± 94.5*i* in the dominant poles of *Tm*, which indicates the meshed DC DPS is unstable and have an oscillation with approximately period 0.066 s. The results of the proposed stability criterion *Tm* and the previous stability criterion *Tm*<sup>2</sup> are not consistent, and the results of proposed stability criterion *Tm* will be confirmed to be right with the nonlinear dynamic simulations in Section 4. Moreover, many poles related to the network dynamics are also lost in *Tm*<sup>2</sup> compared to *Tm*.

**Figure 10.** The dominant poles of *Tm* ans *Tm*<sup>2</sup> when the cables are long.

With the electrical parameters listed in Table B1, when the power of resistive load 4 *PL*<sup>4</sup> is 50 kW, which means *RL*<sup>4</sup> = 5 Ω, the dominant poles of *Tm* are displayed as Figure 11. When the power of resistive load 4 *PL*<sup>4</sup> reduces from 50 kW to 10 kW, which means *RL*<sup>4</sup> = 25 Ω, the dominant poles of *Tm* are displayed as Figure 12. Furthermore, when the power of resistive load 4 *PL*<sup>4</sup> is 10 kW and the absorbed power of CPL2 *Pc*<sup>2</sup> reduces from 100 kW to 50 kW, which means *RL*<sup>4</sup> = 25 Ω and *Rc*<sup>2</sup> = 1.25 Ω, and the dominant poles of *Tm* are displayed as Figure 13.

**Figure 11.** The dominant poles of *Tm* when *PL*<sup>4</sup> = 50 kW, *Pc*<sup>2</sup> = 100 kW.

**Figure 12.** The dominant poles of *Tm* when *PL*<sup>4</sup> = 10 kW, *Pc*<sup>2</sup> = 100 kW.

**Figure 13.** The dominant poles of *Tm* when *PL*<sup>4</sup> = 10 kW, *Pc*<sup>2</sup> = 50 kW.

As can be seen from Figure 11, all of the poles of *Tm* are in the left half plane, which indicates that the meshed DC DPS is stable. But there are two right half plane poles 0.48 ± 97*i* in the Figure 12, which indicates the meshed DC DPS is unstable and have an oscillation with approximately period 0.0647 s. Furthermore, with the power reducation of CPL2, the meshed DC DPS becomes stable as shown in Figure 13.

Applying the proposed stability enhancement method in the sources connected with node 1 and node 2 of the meshed DC DPS when *PL*<sup>4</sup> is 10 kW and *Pc*<sup>2</sup> is 100 kW, the dominant poles of *Tm* are displayed as Figure 14. There are no right half plane poles, which means that the meshed DC DPS is stable. Compared to the poles in Figure 12, the stability of the meshed DC DPS is enhanced by the proposed stability enhancement method.

**Figure 14.** The dominant poles of *Tm* with the proposed stability enhancement method.

#### **4. Nonlinear Dynamic Simulations**

To validate the effectiveness of the proposed stability criterion and stability enhancement method, the meshed DC DPS in Figure 1 is built in the nonlinear time-domain simulation by using Matlab/Simulink. The electrical and controller parameters are all given in Appendix B.

To compare the proposed stability criterion *Tm* and the previous stability criterion *Tm*2, the cables are chosen to be long cables with parameters *Rcablei* = 0.05 Ω , *Lcablei* = 0.5 mH, *Ccablei* = 100 μF, *i* = 1∼4, the nonlinear dynamic simulation results of the node 1 voltage *v*1(*t*) and the cable 1 current *icable*1(*t*) are shown in Figure 15. As can be seen in Figure 15, before 2.5 s, the meshed DC DPS gradually becomes to a stable state with *PL*<sup>4</sup> = 50 kW, and once the power of resistive load 4 *PL*<sup>4</sup> reduces to 15 kW at 2.5 s, the meshed DC DPS becomes to an oscillation state. The oscillation amplitude gradually grows, and the oscillation period is about 0.066 s, which satisfy the theoretical analysis in Figure 10. The results verify the effectiveness of the proposed stability criterion *Tm* compared to the previous stability criterion *Tm*<sup>2</sup> which does not take the network dynamics into account.

**Figure 15.** Comparison between *Tm* and *Tm*2.

To verify the effectiveness of the proposed stability criterion *Tm* with power changes of the resistive load and CPL, and the proposed stability enhancement method, the nonlinear dynamic simulation results are shown as the following Figures 16–18.

As can be seen in Figure 16, before 1.5 s, the meshed DC DPS gradually becomes to stable state with *PL*<sup>4</sup> = 50 k*W* and *Pc*<sup>2</sup> = 100 k*W*, which confirms the theoretical analysis in Figure 11. Then, the power of resistive load 4 *PL*<sup>4</sup> reduces to 10 kW at 1.5 s, the meshed DC DPS becomes to the oscillation state with about period 0.065 s, which confirms the theoretical analysis in Figure 12. The power of CPL2 reduces from 100 kW to 50 kW at 2 s, and the meshed DC DPS becomes to a stable state again, which confirms the theoretical analysis in Figure 13.

**Figure 16.** The node 1 voltage *v*1(*t*) with power changes of the resistive load and CPL.

**Figure 17.** The node 1 voltage *v*<sup>1</sup> without/with the proposed stability enhancement method.

**Figure 18.** The currents of cable 1 *icable*<sup>1</sup> and cable 4 *icable*4.

As can be seen in Figure 17, when the power of CPL2 increases from 50 kW to 100 kW at 3.5 s, the meshed DC DPS becomes to the oscillation state again between 3.5–4.5 s. Then, the vitural harmonic resistances are added in the controllers of source 1 and source 2 at 4.5 s. With the vitural harmonic resistance, the meshed DC DPS becomes to the stable state again, and the node 1 voltage *v*<sup>1</sup> converges again to the steady value 500 V. The simulation results between 4.5–5.5 s confirm the theoretical analysis results in Figure 14. The power of resistive load 4 *PL*<sup>4</sup> increases to 50 kW at 5.5 s, and the meshed DC DPS is maintained at a steady state.

The currents of cable 1 *icable*<sup>1</sup> and cable 4 *icable*<sup>4</sup> are displayed as Figure 18. The steady state and unstable state are similar to the node 1 voltage *v*<sup>1</sup> as shown in Figures 16 and 17. It is worth noting that the steady state current values in 0–1.5 s and 5.5–6 s are the same, which means that with or without the proposed stability enhancement method, the current distributions in the meshed DC DPS are not affected, which verifies that the low frequency characteristic of the source output impedance has not been modified.

#### **5. Conclusions**

This paper discusses the modeling and stability analysis of the meshed DC DPS including multiple sources with voltage-controlled converters and CPLs, and proposes a stability enhancement method. The stability criterion is found to relate only with the network node admittance matrix, the output impendence of the source, and the input admittance of the CPLs, and the three parts can be obtained separately. Virtual harmonic resistance through the second-order generalized integrator (SOGI) for the source proposed in this paper is a decoupling design which is only based on the characteristic of the source output impedance. Nonlinear dynamic time-domain simulation results verify the effectiveness of the proposed stability criterion and stability enhancement method.

Furthermore, the current control sources or loads can be treated as BCCCs, and the input admittances of current control sources or loads are related to the stability analysis. Multiple sources with droop control can be treated as BVCCs, and the output impendence of the droop control sources are related to the stability analysis. So the proposed stability criterion can also deal with the current control sources or loads and multiple sources with droop control. Related analyses are in the future researches.

**Acknowledgments:** The research leading to these results has received funding from the Nation "863" project of China (2014AA052001), the National Natural Science Foundation of China (51377142), the Zhejiang Provincial Natural Science Foundation of China (LY16E070002), and the Key Research and Development Project of Zhejiang Province (2016C01G2010726).

**Author Contributions:** All authors contributed equally in the presented research, and all authors approved the final manuscript.

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

#### **Appendix A**

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎩

⎧ ⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩ *<sup>a</sup>*<sup>0</sup> = *<sup>D</sup>*<sup>2</sup> − *Gc*<sup>2</sup>*<sup>v</sup> Ic*2*Kc*2*pwmRc*2*<sup>D</sup> <sup>a</sup>*<sup>1</sup> = −*Rc*2*Rd f* <sup>2</sup>*Cd f* <sup>2</sup>*Gc*<sup>2</sup>*<sup>v</sup> Ic*2*Kc*2*pwmD* + (*Cd f* <sup>2</sup> + *Cf* <sup>2</sup>)(*Gc*2*vKc*2*pwmRc*2*Vf* <sup>2</sup> + *Rc*<sup>2</sup> + *rc*2)+(*Rc*2*Cc*<sup>2</sup> + *Rd f* <sup>2</sup>*Cd f* <sup>2</sup>)*D*<sup>2</sup> *a*<sup>2</sup> = *Rc*2*Rd f* <sup>2</sup>*Cd f* <sup>2</sup>(*Cf* <sup>2</sup>*Gc*2*vKc*2*pwmVf* <sup>2</sup> + *Cc*2*D*2) + *Cd f* <sup>2</sup>*Cf* <sup>2</sup>*Rd f* <sup>2</sup>(*Rc*<sup>2</sup> + *rc*2)+(*Lc*<sup>2</sup> + *Cc*2*Rc*2*rc*2)(*Cd f* <sup>2</sup> + *Cf* <sup>2</sup>) *a*<sup>3</sup> = (*rc*2*Rc*2*Cc*<sup>2</sup> + *Lc*2)*Rd f* <sup>2</sup>*Cd f* <sup>2</sup>*Cf* <sup>2</sup> + *Cc*2*Lc*2*Rc*2(*Cf* <sup>2</sup> + *Cd f* <sup>2</sup> *a*<sup>4</sup> = *Rc*2*Cc*2*Rd f* <sup>2</sup>*Cd f* <sup>2</sup>*Lc*2*Cf* <sup>2</sup> (A1) *b*<sup>0</sup> = *Gc*2*vKc*2*pwmRc*2*Vf* <sup>2</sup> + *Rc*<sup>2</sup> + *rc*<sup>2</sup> *<sup>b</sup>*<sup>1</sup> = (*Rc*2*Rd f* <sup>2</sup>*Cd f* <sup>2</sup>*Vf* <sup>2</sup> − *Rc*2*D Ic*2*Lf* <sup>2</sup>)*Gc*2*vKc*<sup>2</sup>*pwm* + *Rc*2*Cc*2*rc*<sup>2</sup> + *Rd f* <sup>2</sup>*Cd f* <sup>2</sup>(*Rc*<sup>2</sup> + *rc*2) + *<sup>D</sup>*<sup>2</sup>*Lf* <sup>2</sup> + *Lc*<sup>2</sup> *<sup>b</sup>*<sup>2</sup> = −*Rc*2*Rd f* <sup>2</sup>*Cd f* <sup>2</sup>*D Ic*2*Gc*2*vKc*2*pwmLf* <sup>2</sup> + (*Cd f* <sup>2</sup> + *Cf* <sup>2</sup>)(*Rc*2*Vf* <sup>2</sup>*Gc*2*vKc*2*pwmLf* <sup>2</sup> + *Rc*2*Lf* <sup>2</sup> + *rc*2*Lf* <sup>2</sup>)+(*Rc*2*Cc*<sup>2</sup> + *Rd f* <sup>2</sup>*Cd f* <sup>2</sup>)(*Lf* <sup>2</sup>*D*<sup>2</sup> + *Lc*2) + *Rc*2*Cc*2*Rd f* <sup>2</sup>*Cd f* <sup>2</sup>*rc*<sup>2</sup> *b*<sup>3</sup> = *Rc*2*Rd f* <sup>2</sup>*Cd f* <sup>2</sup>*Cf* <sup>2</sup>*Gc*2*vKc*2*pwmLf* <sup>2</sup>*Vf* <sup>2</sup> + (*D*<sup>2</sup>*Lf* <sup>2</sup> + *Lc*2)*Rc*2*Cc*2*Rd f* <sup>2</sup>*Cd f* <sup>2</sup> + *Cd f* <sup>2</sup>*Cf* <sup>2</sup>*Lf* <sup>2</sup>*Rd f* <sup>2</sup>(*Rc*<sup>2</sup> + *rc*2)+(*Cd f* <sup>2</sup> + *Cf* <sup>2</sup>)(*Cc*2*Lf* <sup>2</sup>*Rc*2*rc*<sup>2</sup> + *Lc*2*Lf* <sup>2</sup>) *b*<sup>4</sup> = *Cd f* <sup>2</sup>*Cf* <sup>2</sup>*Lf* <sup>2</sup>*Rd f* <sup>2</sup>(*Rc*2*rc*2*Cc*<sup>2</sup> + *Lc*<sup>2</sup> + *Cd f* <sup>2</sup> + *Cf* <sup>2</sup>)*Cc*2*Lc*2*Lf* <sup>2</sup>*Rc*<sup>2</sup> *b*<sup>5</sup> = *Cc*2*Cd f* <sup>2</sup>*Cf* <sup>2</sup>*Lc*2*Lf* <sup>2</sup>*Rc*2*Rd f* <sup>2</sup> (A2) ⎧ ⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎩ *c*<sup>0</sup> = *Gc*2*vKc*2*pwmVf* <sup>2</sup>*D* + *Gc*<sup>2</sup>*<sup>v</sup> Ic*2*Kc*<sup>2</sup>*pwm*(*Rc*<sup>2</sup> + *rc*2) *c*<sup>1</sup> = (*Cc*2*Rc*<sup>2</sup> + *Cd f* <sup>2</sup>*Rd f* <sup>2</sup>)*Gc*2*vKc*2*pwmVf* <sup>2</sup>*D* + *Lc*<sup>2</sup> + *Rc*2*Rd f* <sup>2</sup>*Cd f* <sup>2</sup> + *rc*2*Rd f* <sup>2</sup>*Cd f* <sup>2</sup> + *rc*2*Rc*2*Cc*2)*Gc*2*vKc*<sup>2</sup>*pwm Ic*<sup>2</sup> *c*<sup>2</sup> = *Rc*2*Cc*2*Rd f* <sup>2</sup>*Cd f* <sup>2</sup>*Gc*2*vKc*2*pwmVf* <sup>2</sup>*D* + *Rc*2*Cc*2*Lc*<sup>2</sup> + *Rd f* <sup>2</sup>*Cd f* <sup>2</sup>*Lc*<sup>2</sup> + *Rc*2*Cc*2*Rd f* <sup>2</sup>*Cd f* <sup>2</sup>*rc*2)*Gc*<sup>2</sup>*<sup>v</sup> Ic*2*Kc*<sup>2</sup>*pwm* (A3)

```
c3 = Rc2Cc2Rd f 2Cd f 2Lc2Gc2v Ic2Kc2pwm
```
#### **Appendix B**

This Appendix lists all the main electrical parameters of the studied DC meshed DPS as shown in Figure 1.


**Table B1.** Main electrical parameters of the studied DC meshed distributed power systems (DPS).

#### **References**


© 2017 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/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Energies* Editorial Office E-mail: energies@mdpi.com www.mdpi.com/journal/energies

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18