A Vector-Controlled Distributed Generator Model for a Power Flow Based on a Three-Phase Current Injection Method

Abstract

This paper proposes a vector-controlled distributed generator (DG) model for a power flow based on a three-phase current injection method (TCIM). In order to represent the DG models in the power flow, steady-state phase current output equations are formulated. Using these equations, the TCIM power flow formulation is modified to include the DG models. In the proposed power flow, a DG-connected bus is modeled as either a load bus (PQ bus) or a voltage-controlled bus (PV bus), depending on the control mode of the reactive power. However, unlike conventional bus models, the values of the DG-connected bus models are represented by three-phase quantities: three-phase active and reactive power output for a PQ bus, and three-phase active power and positive-sequence voltage for a PV bus. In addition, a method is proposed for representing the reactive power limit of a voltage-control-mode DG by using the q-axis current limit. Utilizing a modified IEEE 13-bus test system, the accuracy of the proposed method is verified by comparison to the power systems computer aided design (PSCAD) model. Furthermore, the effect of the number of DGs on the convergence rate is analyzed, using the IEEE 123-bus test system.

1. Introduction

Recently, integration of distributed generators (DGs) into distribution systems has been increasing, owing to their various environmental and economic advantages. Moreover, various techniques have been proposed for utilizing DGs as system operational resources, such as voltage and reactive power (Volt/VAR) control [1,2,3]. However, the characteristics of DGs (especially inverter-interfaced DGs) are quite different from those of conventional generators. Therefore, to achieve stable and effective operation of a distribution system with DGs, the operating characteristics of the DGs should be taken into account.

Power flow studies are the backbone of power system design, analysis, and operation [4]. Three-phase power flow methods are generally adopted for a distribution system, since the operating conditions are inherently unbalanced. Various DG models for unbalanced three-phase power flow studies have been proposed. In [5,6,7,8,9], a three-phase connected DG that operates in constant reactive power mode was modeled as a negative constant power load. However, this is not appropriate, since the phase outputs of a DG depend on the phase voltages under unbalanced operating conditions [10]. Consequently, power flow results obtained with these models may have errors. DG models considering unbalanced operating conditions to improve the accuracy of the power flow results were presented in [10,11,12,13,14]. Models for synchronous generator-based DGs were proposed in [10,11,12], and voltage-source converter (VSC)-coupled DG models (including output filters) were introduced in [12,13]. In [14], the power flow model for an induction generator-based wind turbine generating system was proposed. In [10], the phase-frame approach was used (i.e., the variables are represented by a-b-c phase quantities), while sequence-frame analysis was adopted in [11,12,13,14] (i.e., the variables are represented by positive-negative-zero sequence quantities). In [15], it was reported that a steady-state voltage-controlled model could be represented by a controllable three-phase balanced voltage source, and a current-controlled VSC model by a controllable three-phase balanced current source.

In [16], phase-frame-based power flow, based on a three-phase current injection method (TCIM), was proposed for robust and fast unbalanced three-phase power flow studies. In this method, the power flow problem was formulated using current mismatch equations written in rectangular coordinates, and solved via the Newton-Raphson method. Voltage control device models for TCIM power flow were presented in [17]. In [18], a step-size optimization factor was proposed to improve the convergence characteristics of TCIM power flow. According to [19], the TCIM power flow method is superior to the forward/backward sweep method, especially in cases involving a large-scale system, heavy loading conditions, and/or a system with numerous voltage control devices.

In this paper, a modified TCIM power flow incorporating vector-controlled DGs is proposed. The remainder of the paper is divided into four sections, organized as follows. In Section 2, a steady-state phase output current model for a DG is formulated. The implementation of TCIM power flow with the DG model formulation is described in Section 3. A method for representing the reactive power limit of a voltage-controlled bus (PV bus) is also proposed. In Section 4, the accuracy of the proposed power flow is verified by comparing its results to those obtained from the power systems computer aided design (PSCAD) model, utilizing a modified IEEE 13-bus test system. Moreover, the effect of the number of DGs on the convergence rate of the power flow is analyzed, using the IEEE 123-bus test system. Finally, Section 5 contains concluding remarks.

2. Power Flow Model for a DG

Because a vector-controlled DG controls three-phase outputs rather than single-phase outputs, a-b-c phase outputs should be modeled to represent the DG in the TCIM power flow. In this section, DGs are first classified into three types according to the topology of the output filter and the point of current control. Steady-state phase current output equations are formulated for each type.

2.1. DG Overview

Vector control theory is widely used to control the outputs of VSC-coupled DGs, since the active and reactive power can be controlled independently. Figure 1a shows a general power control block for a vector-controlled DG [20]. The three-phase active and reactive power reference values, P_ref and Q_ref, are used to calculate the reference values of the d- and q-axis currents, i_d,ref and i_q,ref. As another option, the dc link voltage can be used for active power control, and the three-phase ac voltage can be used for reactive power control. Figure 1b shows three different DG configurations determined by the topology of the output filter [21]. For a DG with an L filter, the current controller controls the grid-side current, I_DG. On the other hand, either the VSC output current (I_VSC) or I_DG can be used as a control variable for a DG equipped with an LC or LCL filter [20,22,23].

Figure 1. Overview of a vector-controlled DG: (a) power control block and (b) output filter topology.

Figure 1. Overview of a vector-controlled DG: (a) power control block and (b) output filter topology.

2.2. Steady-State Phase Current Output Models

In the steady state, the a-b-c phase currents at the current control point are balanced, since the d- and q-axis currents are maintained at specific values by the controller. Therefore, the steady-state output of a DG can be modeled as a balanced three-phase current source (BTCS) [15]. Table 1 lists three types of steady-state DG models, according to the topology of the output filter and the controlled current. Because the steady-state phase current output equations of the Type 1 and 2 DG models can be formulated by modifying those of the Type 3 DG model, the current output model for the Type 3 DG is first formulated, and the modifications required to represent the other types are discussed at the end of the section.

Table 1. DG types and corresponding equivalent circuits.

Property	Type 1 DG	Type 2 DG	Type 3 DG
Output filter	L, LC, LCL	LC	LCL
Controlled current	IDG	IVSC	IVSC
Equivalent circuit

Table 1. DG types and corresponding equivalent circuits.

Property	Type 1 DG	Type 2 DG	Type 3 DG
Output filter	L, LC, LCL	LC	LCL
Controlled current	IDG	IVSC	IVSC
Equivalent circuit

The phase current outputs of a BTCS are given by:

$$I_{BTCS}^a=|I_{BTCS}|e^{j\theta }$$

(1)

$$I_{BTCS}^b=|I_{BTCS}|e^{j(\theta -\frac{2}{3}\pi )}=I_{BTCS}^a{e}^{-j\frac{2}{3}\pi }$$

(2)

$$I_{BTCS}^c=|I_{BTCS}|e^{j(\theta +\frac{2}{3}\pi )}=I_{BTCS}^a{e}^{j\frac{2}{3}\pi }$$

(3)

where |I_BTCS| and θ represent the magnitude and phase angle, respectively, of the phase a current phasor. From the equivalent circuit shown in Table 1, the output current of phase s for a Type 3 DG connected to bus k can be written as:

$$I_{DG,k}^s=\frac{R_d+j{X}_C}{R_{eq,k}+j{X}_{eq,k}}{I}_{BTCS}^s-\frac{V_k^s}{R_{eq,k}+j{X}_{eq,k}}$$

(4)

where V_k^s is the voltage of phase s for bus k; and R_eq,k and X_eq,k are the equivalent resistance and reactance of the output filter, defined by:

$$R_{eq,k}=R_d+R_2$$

(5)

$$X_{eq,k}=X_C+X_{L2}$$

(6)

The three-phase complex power output of the DG is the sum of the phase power outputs:

$$P_{DG,k}^{3\varphi }+j{Q}_{DG,k}^{3\varphi }=V_k^a{(I_{DG,k}^a)}^{*}+V_k^b{(I_{DG,k}^b)}^{*}+V_k^c{(I_{DG,k}^c)}^{*}$$

(7)

where the superscript * denotes the complex conjugate. By substituting Equations (1)–(3) into Equation (4), and then substituting Equations (1)–(3) into Equation (7), the three-phase complex power output can be written as:

$$P_{DG,k}^{3\varphi }+j{Q}_{DG,k}^{3\varphi }=3\frac{R_d-j{X}_C}{R_{eq,k}-j{X}_{eq,k}}{V}_k^{+}{(I_{BTCS}^a)}^{*}-\frac{|V_k^a{|}^2+|V_k^b{|}^2+|V_k^c{|}^2}{R_{eq,k}-j{X}_{eq,k}}$$

(8)

where V_k⁺ is the positive sequence voltage, defined by:

$$V_k^{+}=\frac{V_k^a+V_k^b{e}^{j\frac{2}{3}\pi }+V_k^c{e}^{-j\frac{2}{3}\pi }}{3}$$

(9)

By rearranging Equation (8), the current output of the BTCS for phase a can be calculated from:

$$I_{BTCS}^a=\frac{(P_{DG,k}^{3\varphi }-j{Q}_{DG,k}^{3\varphi })(R_{eq,k}+j{X}_{eq,k})+|V_k^a{|}^2+|V_k^b{|}^2+|V_k^c{|}^2}{3{(V_k^{+})}^{*}(R_d+j{X}_C)}$$

(10)

Finally, the phase current outputs can be formulated as functions of the three-phase active and reactive power outputs, phase voltages, and impedances of the output filter by substituting Equation (10) into Equations (1)–(3), and then substituting Equations (1)–(3) into Equation (4):

$$I_{DG,k}^a=\frac{P_{DG,k}^{3\varphi }-j{Q}_{DG,k}^{3\varphi }}{3{(V_k^{+})}^{*}}+\frac{(R_{eq,k}+j{X}_{eq,k})(|V_k^a{|}^2+|V_k^b{|}^2+|V_k^c{|}^2)}{3{(V_k^{+})}^{*}}-\frac{V_k^a}{R_{eq,k}+j{X}_{eq,k}}$$

(11)

$$I_{DG,k}^b=\frac{P_{DG,k}^{3\varphi }-j{Q}_{DG,k}^{3\varphi }}{3{(V_k^{+})}^{*}}{e}^{-j\frac{2}{3}\pi }+\frac{(R_{eq,k}+j{X}_{eq,k})(|V_k^a{|}^2+|V_k^b{|}^2+|V_k^c{|}^2)}{3{(V_k^{+})}^{*}}{e}^{-j\frac{2}{3}\pi }-\frac{V_k^b}{R_{eq,k}+j{X}_{eq,k}}$$

(12)

$$I_{DG,k}^c=\frac{P_{DG,k}^{3\varphi }-j{Q}_{DG,k}^{3\varphi }}{3{(V_k^{+})}^{*}}{e}^{j\frac{2}{3}\pi }+\frac{(R_{eq,k}+j{X}_{eq,k})(|V_k^a{|}^2+|V_k^b{|}^2+|V_k^c{|}^2)}{3{(V_k^{+})}^{*}}{e}^{j\frac{2}{3}\pi }-\frac{V_k^c}{R_{eq,k}+j{X}_{eq,k}}$$

(13)

For the Type 1 DG, only the first term in each of the Equations (11)–(13) will be used, since there are no shunt or series impedances. In other words, the Type 1 DG injects balanced three-phase currents even when the voltages are unbalanced. The current output equations of the Type 2 DG have the same formulation as Equations (11)–(13), except that R_eq,k is equal to R_d and X_eq,k is equal to X_C.

3. Implementation of the TCIM Power Flow

In a power flow, a DG-connected bus can be modeled as either a load bus (PQ bus) or a PV bus, depending on the reactive power control mode of the DG. However, unlike conventional bus models, the given values of the models are three-phase quantities: three-phase active and reactive power outputs for a PQ bus, and three-phase active power and positive-sequence voltage for a PV bus. In this section, a method is proposed for implementing a TCIM-based power flow incorporating DG output models. A power flow method for the Type 2 and 3 DG models (which have the same phase current output equations) is first proposed. Modifications of the proposed method for the Type 1 DG model are presented at the end of the section.

3.1. Basic Equations

Taking into account the phase current outputs of the DGs, the current mismatch equations in rectangular form for phase s of bus k, presented in [16], are modified as follows:

$$\Delta I_{r,k}^s={\displaystyle \sum _{i=1}^n{\displaystyle \sum _{t=\{a,b,c\}}(G_{ki}^{st}{V}_{r,i}^t-B_{ki}^{st}{V}_{m,i}^t)}}-\frac{P_k^s{V}_{r,k}^s+Q_k^s{V}_{m,k}^s}{{(V_{r,k}^s)}^2+{(V_{m,k}^s)}^2}-I_{DG,r,k}^s$$

(14)

$$\Delta I_{m,k}^s={\displaystyle \sum _{i=1}^n{\displaystyle \sum _{t=\{a,b,c\}}(G_{ki}^{st}{V}_{m,i}^t+B_{ki}^{st}{V}_{r,i}^t)}}-\frac{P_k^s{V}_{m,k}^s-Q_k^s{V}_{r,k}^s}{{(V_{r,k}^s)}^2+{(V_{m,k}^s)}^2}-I_{DG,m,k}^s$$

(15)

where:

$\Delta I_{r,k}^s$, $\Delta I_{m,k}^s$

real and imaginary parts of current mismatch;

$G_{ki}^{st}$, $B_{ki}^{st}$

real and imaginary parts of bus admittance matrix;

$V_{r,k}^s$, $V_{m,k}^s$

real and imaginary parts of phase voltage;

$P_k^s$, $Q_k^s$

active and reactive power injections of load;

$I_{DG,r,k}^s$, $I_{DG,m,k}^s$

real and imaginary parts of output current of the DG.

The last terms of Equations (14) and (15) are newly introduced to represent the current injection of the DG, and are calculated from Equations (11)–(13) as follows:

$$I_{DG,r,k}^a=\frac{{\alpha }_k(P_{DG,k}^{3\varphi }+R_{eq,k}{\gamma }_k)+{\beta }_k(Q_{DG,k}^{3\varphi }+X_{eq,k}{\gamma }_k)}{{\alpha }_k^2+{\beta }_k^2}-\frac{R_{eq,k}{V}_{r,k}^a+X_{eq,k}{V}_{m,k}^a}{R_{eq,k}^2+X_{eq,k}^2}$$

(16)

$$I_{DG,m,k}^a=\frac{-{\alpha }_k(Q_{DG,k}^{3\varphi }+X_{eq,k}{\gamma }_k)+{\beta }_k(P_{DG,k}^{3\varphi }+R_{eq,k}{\gamma }_k)}{{\alpha }_k^2+{\beta }_k^2}-\frac{R_{eq,k}{V}_{m,k}^a-X_{eq,k}{V}_{r,k}^a}{R_{eq,k}^2+X_{eq,k}^2}$$

(17)

$$I_{DG,r,k}^b=-\frac{(-{\alpha }_k+\sqrt{3}{\beta }_k)(P_{DG,k}^{3\varphi }+R_{eq,k}{\gamma }_k)-(\sqrt{3}{\alpha }_k+{\beta }_k)(Q_{DG,k}^{3\varphi }+X_{eq,k}{\gamma }_k)}{2({\alpha }_k^2+{\beta }_k^2)}-\frac{R_{eq,k}{V}_{r,k}^b+X_{eq,k}{V}_{m,k}^b}{R_{eq,k}^2+X_{eq,k}^2}$$

(18)

$$I_{DG,m,k}^b=\frac{(-\sqrt{3}{\alpha }_k-{\beta }_k)(P_{DG,k}^{3\varphi }+R_{eq,k}{\gamma }_k)+({\alpha }_k-\sqrt{3}{\beta }_k)(Q_{DG,k}^{3\varphi }+X_{eq,k}{\gamma }_k)}{2({\alpha }_k^2+{\beta }_k^2)}-\frac{R_{eq,k}{V}_{m,k}^b-X_{eq,k}{V}_{r,k}^b}{R_{eq,k}^2+X_{eq,k}^2}$$

(19)

$$I_{DG,r,k}^c=\frac{(-{\alpha }_k-\sqrt{3}{\beta }_k)(P_{DG,k}^{3\varphi }+R_{eq,k}{\gamma }_k)+(\sqrt{3}{\alpha }_k-{\beta }_k)(Q_{DG,k}^{3\varphi }+X_{eq,k}{\gamma }_k)}{2({\alpha }_k^2+{\beta }_k^2)}-\frac{R_{eq,k}{V}_{r,k}^c+X_{eq,k}{V}_{m,k}^c}{R_{eq,k}^2+X_{eq,k}^2}$$

(20)

$$I_{DG,m,k}^c=\frac{(\sqrt{3}{\alpha }_k-{\beta }_k)(P_{DG,k}^{3\varphi }+R_{eq,k}{\gamma }_k)+({\alpha }_k+\sqrt{3}{\beta }_k)(Q_{DG,k}^{3\varphi }+X_{eq,k}{\gamma }_k)}{2({\alpha }_k^2+{\beta }_k^2)}-\frac{R_{eq,k}{V}_{m,k}^c-X_{eq,k}{V}_{r,k}^c}{R_{eq,k}^2+X_{eq,k}^2}$$

(21)

$${\alpha }_k=V_{r,k}^a-\frac{1}{2}{V}_{r,k}^b-\frac{\sqrt{3}}{2}{V}_{m,k}^b-\frac{1}{2}{V}_{r,k}^c+\frac{\sqrt{3}}{2}{V}_{m,k}^c$$

(22)

$${\beta }_k=V_{m,k}^a-\frac{1}{2}{V}_{m,k}^b+\frac{\sqrt{3}}{2}{V}_{r,k}^b-\frac{1}{2}{V}_{m,k}^c-\frac{\sqrt{3}}{2}{V}_{r,k}^c$$

(23)

$${\gamma }_k=\frac{{(V_{r,k}^a)}^2+{(V_{m,k}^a)}^2+{(V_{r,k}^b)}^2+{(V_{m,k}^b)}^2+{(V_{r,k}^c)}^2+{(V_{m,k}^c)}^2}{R_{eq,k}^2+X_{eq,k}^2}$$

(24)

To represent a PV bus, a three-phase voltage mismatch equation should be defined. Because the DG generally controls the magnitude of the positive-sequence voltage [9], the voltage mismatch equation for bus k is defined as:

$$\Delta V_k^{+}=|V_k^{+}|-V_{ref,k}^{+}$$

(25)

where V⁺_ref,k is the reference value for the magnitude of the positive-sequence voltage; and |V_k⁺| is the magnitude of the positive-sequence voltage calculated in the power flow:

$$|V_k^{+}|=\frac{\sqrt{{\alpha }_k^2+{\beta }_k^2}}{3}$$

(26)

To solve the power flow problem, one must find the unknown voltages that make the mismatch values defined in Equations (14), (15) and (25) equal to zero.

In order to solve the power flow problem via the Newton-Raphson method, the current mismatch vector (ΔI_mr) and voltage vector (V_rm) are defined as:

$$\Delta I_{mr}={\left[\begin{array}{ccccccc}\Delta I_{m,2}^a& \Delta I_{r,2}^a& \Delta I_{m,2}^b& \Delta I_{r,2}^b& \cdots & \Delta I_{m,n}^c& \Delta I_{r,n}^c\end{array}\right]}^t$$

(27)

$$V_{rm}={\left[\begin{array}{ccccccc}{V}_{r,2}^a& V_{m,2}^a& V_{r,2}^b& V_{m,2}^b& \cdots & V_{r,n}^c& V_{m,n}^c\end{array}\right]}^t$$

(28)

Because bus 1 is assumed to be a slack bus, the current mismatch and voltages for bus 1 are not included in Equations (27) and (28). It should be noted that the order of the real and imaginary parts is reversed in the current mismatch vector and voltage vector [16].

As in [24], a three-phase reactive power vector is introduced to represent a PV bus. The three-phase voltage mismatch vector (ΔV⁺) and the three-phase reactive power vector (Q^3φ) are defined as:

$$\Delta V^{+}={\left[\begin{array}{ccc}\Delta V_k^{+}& \cdots & \Delta V_l^{+}\end{array}\right]}^t$$

(29)

$$Q^{3\varphi }={\left[\begin{array}{ccc}{Q}_{DG,k}^{3\varphi }& \cdots & Q_{DG,l}^{3\varphi }\end{array}\right]}^t$$

(30)

In the hth iteration of the Newton-Raphson method, the state variable vectors V_rm and Q^3φ are updated using the mismatch vectors ΔI_mr and ΔV⁺, as follows:

$${\left[\begin{array}{c}{V}_{rm}\\ Q^{3\varphi }\end{array}\right]}^{h+1}={\left[\begin{array}{c}{V}_{rm}\\ Q^{3\varphi }\end{array}\right]}^h-{\left[J\right]}^{-1}{\left[\begin{array}{c}\Delta I_{mr}\\ \Delta V^{+}\end{array}\right]}^h$$

(31)

where J is the Jacobian matrix.

3.2. Structure of the Jacobian Matrix

The Jacobian matrix is the sensitivity matrix between the state variable vectors and the mismatch vectors. The Jacobian matrix can be partitioned as follows:

$$J=\left[\begin{array}{cc}\frac{\partial \Delta I_{mr}}{\partial V_{rm}}=A& \frac{\partial \Delta I_{mr}}{\partial Q^{3\varphi }}=B\\ \frac{\partial \Delta V^{+}}{\partial V_{rm}}=C& \frac{\partial \Delta V^{+}}{\partial Q^{3\varphi }}=O\end{array}\right]$$

(32)

Because there is no explicit equation for the relationship between the three-phase voltage mismatch and the three-phase reactive power, the corresponding submatrix should be calculated via the chain rule, as follows:

$$\frac{\partial \Delta V^{+}}{\partial Q^{3\varphi }}=\left[\frac{\partial \Delta V^{+}}{\partial V_{rm}}\right]{\left[\frac{\partial \Delta I_{mr}}{\partial V_{rm}}\right]}^{-1}\left[\frac{\partial \Delta I_{mr}}{\partial Q^{3\varphi }}\right]$$

(33)

However, if the submatrix obtained from Equation (33) is used, the Jacobian matrix does not have full rank, and thus it is not invertible. In order to avoid this problem, the submatrix is replaced by the zero matrix, as specified in Equation (32).

Figure 2 shows a simple radial system with three DGs and the corresponding Jacobian structure. Since there are four buses (excluding the slack bus), and each bus has six terms pertaining to current mismatches and rectangular voltages, ΔI_mr and V_rm are 24 × 1 column vectors. Since two DGs (DG 1 and DG 3) operate in voltage control mode, Q^3φ and ΔV⁺ are 2 × 1 column vectors. Therefore, the Jacobian matrix is a 26 × 26 matrix, as shown in Figure 2b.

In the Jacobian matrix, the submatrices of the partition A, B, and C are a 6 × 6 matrix, a 6 × 1 column vector, and a 1 × 6 row vector, respectively. Because the three-phase reactive output of the DG affects only the current mismatches of the DG-connected buses, B₃ and B₅ have nonzero values, and the other vectors in B are zero vectors. Similarly, C₃ and C₅ are the only nonzero vectors in C.

Figure 2. Jacobian structure: (a) five bus system and (b) corresponding Jacobian matrix.

Figure 2. Jacobian structure: (a) five bus system and (b) corresponding Jacobian matrix.

3.3. Calculation of the Jacobian Matrix

The elements of the matrix A can be calculated from the partial derivatives of Equations (14) and (15) with respect to the rectangular voltages. For the elements of the off-diagonal block matrices, Equation (12) of the original TCIM [16] can be used without any modification, since these terms are not related to the current outputs of the DG. The diagonal terms, on the other hand, should be modified as follows:

$$A_{kk}=A_{kk,O}-A_{kk,DG}$$

(34)

The first term, A_kk,O, is identical to the matrix of Equation (13) in [16], and the second term, A_kk,DG, represents the partial derivatives of the current outputs of the DG:

$$A_{kk,DG}=\left[\begin{array}{cccccc}\frac{\partial I_{DG,m,k}^a}{\partial V_{r,k}^a}& \frac{\partial I_{DG,m,k}^a}{\partial V_{m,k}^a}& \frac{\partial I_{DG,m,k}^a}{\partial V_{r,k}^b}& \frac{\partial I_{DG,m,k}^a}{\partial V_{m,k}^b}& \frac{\partial I_{DG,m,k}^a}{\partial V_{r,k}^c}& \frac{\partial I_{DG,m,k}^a}{\partial V_{m,k}^c}\\ \frac{\partial I_{DG,r,k}^a}{\partial V_{r,k}^a}& \frac{\partial I_{DG,r,k}^a}{\partial V_{m,k}^a}& \frac{\partial I_{DG,r,k}^a}{\partial V_{r,k}^b}& \frac{\partial I_{DG,r,k}^a}{\partial V_{m,k}^b}& \frac{\partial I_{DG,r,k}^a}{\partial V_{r,k}^c}& \frac{\partial I_{DG,r,k}^a}{\partial V_{m,k}^c}\\ \frac{\partial I_{DG,m,k}^b}{\partial V_{r,k}^a}& \frac{\partial I_{DG,m,k}^b}{\partial V_{m,k}^a}& \frac{\partial I_{DG,m,k}^b}{\partial V_{r,k}^b}& \frac{\partial I_{DG,m,k}^b}{\partial V_{m,k}^b}& \frac{\partial I_{DG,m,k}^b}{\partial V_{r,k}^c}& \frac{\partial I_{DG,m,k}^b}{\partial V_{m,k}^c}\\ \frac{\partial I_{DG,r,k}^b}{\partial V_{r,k}^a}& \frac{\partial I_{DG,r,k}^b}{\partial V_{m,k}^a}& \frac{\partial I_{DG,r,k}^b}{\partial V_{r,k}^b}& \frac{\partial I_{DG,r,k}^b}{\partial V_{m,k}^b}& \frac{\partial I_{DG,r,k}^b}{\partial V_{r,k}^c}& \frac{\partial I_{DG,r,k}^b}{\partial V_{m,k}^c}\\ \frac{\partial I_{DG,m,k}^c}{\partial V_{r,k}^a}& \frac{\partial I_{DG,m,k}^c}{\partial V_{m,k}^a}& \frac{\partial I_{DG,m,k}^c}{\partial V_{r,k}^b}& \frac{\partial I_{DG,m,k}^c}{\partial V_{m,k}^b}& \frac{\partial I_{DG,m,k}^c}{\partial V_{r,k}^c}& \frac{\partial I_{DG,m,k}^c}{\partial V_{m,k}^c}\\ \frac{\partial I_{DG,r,k}^c}{\partial V_{r,k}^a}& \frac{\partial I_{DG,r,k}^c}{\partial V_{m,k}^a}& \frac{\partial I_{DG,r,k}^c}{\partial V_{r,k}^b}& \frac{\partial I_{DG,r,k}^c}{\partial V_{m,k}^b}& \frac{\partial I_{DG,r,k}^c}{\partial V_{r,k}^c}& \frac{\partial I_{DG,r,k}^c}{\partial V_{m,k}^c}\end{array}\right]$$

(35)

The elements of A_kk,DG can be obtained from the partial derivatives of Equations (16)–(21), as follows:

$$\begin{array}{ll}{A}_{kk,DG}& =\frac{1}{2{({\alpha }_k^2+{\beta }_k^2)}^2}\left[\begin{array}{cccccc}{f}_1& f_2& f_3& f_4& f_5& f_6\\ -f_2& f_1& -f_4& f_3& -f_6& f_5\\ f_3& f_4& f_5& f_6& f_1& f_2\\ -f_4& f_3& -f_6& f_5& -f_2& f_1\\ f_5& f_6& f_1& f_2& f_3& f_4\\ -f_6& f_5& -f_2& f_1& -f_4& f_3\end{array}\right]-\frac{1}{R_{eq,k}^2+X_{eq,k}^2}\left[\begin{array}{cccccc}-X_{eq,k}& R_{eq,k}& 0& 0& 0& 0\\ R_{eq,k}& X_{eq,k}& 0& 0& 0& 0\\ 0& 0& -X_{eq,k}& R_{eq,k}& 0& 0\\ 0& 0& R_{eq,k}& X_{eq,k}& 0& 0\\ 0& 0& 0& 0& -X_{eq,k}& R_{eq,k}\\ 0& 0& 0& 0& R_{eq,k}& X_{eq,k}\end{array}\right]\\ & +\left[\begin{array}{cccccc}2h_1{V}_{r,k}^a& 2h_1{V}_{m,k}^a& 2h_1{V}_{r,k}^b& 2h_1{V}_{m,k}^b& 2h_1{V}_{r,k}^c& 2h_1{V}_{m,k}^c\\ 2h_2{V}_{r,k}^a& 2h_2{V}_{m,k}^a& 2h_2{V}_{r,k}^b& 2h_2{V}_{m,k}^b& 2h_2{V}_{r,k}^c& 2h_2{V}_{m,k}^c\\ -(h_1+\sqrt{3}{h}_2)V_{r,k}^a& -(h_1+\sqrt{3}{h}_2)V_{m,k}^a& -(h_1+\sqrt{3}{h}_2)V_{r,k}^b& -(h_1+\sqrt{3}{h}_2)V_{m,k}^b& -(h_1+\sqrt{3}{h}_2)V_{r,k}^c& -(h_1+\sqrt{3}{h}_2)V_{m,k}^c\\ (\sqrt{3}{h}_1-h_2)V_{r,k}^a& (\sqrt{3}{h}_1-h_2)V_{m,k}^a& (\sqrt{3}{h}_1-h_2)V_{r,k}^b& (\sqrt{3}{h}_1-h_2)V_{m,k}^b& (\sqrt{3}{h}_1-h_2)V_{r,k}^c& (\sqrt{3}{h}_1-h_2)V_{m,k}^c\\ -(h_1-\sqrt{3}{h}_2)V_{r,k}^a& -(h_1-\sqrt{3}{h}_2)V_{m,k}^a& -(h_1-\sqrt{3}{h}_2)V_{r,k}^b& -(h_1-\sqrt{3}{h}_2)V_{m,k}^b& -(h_1-\sqrt{3}{h}_2)V_{r,k}^c& -(h_1-\sqrt{3}{h}_2)V_{m,k}^c\\ -(\sqrt{3}{h}_1+h_2)V_{r,k}^a& -(\sqrt{3}{h}_1+h_2)V_{m,k}^a& -(\sqrt{3}{h}_1+h_2)V_{r,k}^b& -(\sqrt{3}{h}_1+h_2)V_{m,k}^b& -(\sqrt{3}{h}_1+h_2)V_{r,k}^c& -(\sqrt{3}{h}_1+h_2)V_{m,k}^c\end{array}\right]\end{array}$$

(36)

where:

$$f_1=-4{\alpha }_k{\beta }_k(P_{DG,k}^{3\varphi }+{\gamma }_k{R}_{eq,k})+2({\alpha }_k^2-{\beta }_k^2)(Q_{DG,k}^{3\varphi }+{\gamma }_k{X}_{eq,k})$$

(37)

$$f_2=2({\alpha }_k^2-{\beta }_k^2)(P_{DG,k}^{3\varphi }+{\gamma }_k{R}_{eq,k})+4{\alpha }_k{\beta }_k(Q_{DG,k}^{3\varphi }+{\gamma }_k{X}_{eq,k})$$

(38)

$$f_3=(\sqrt{3}{\alpha }_k^2+2{\alpha }_k{\beta }_k-\sqrt{3}{\beta }_k^2)(P_{DG,k}^{3\varphi }+{\gamma }_k{R}_{eq,k})-({\alpha }_k^2-2\sqrt{3}{\alpha }_k{\beta }_k-{\beta }_k^2)(Q_{DG,k}^{3\varphi }+{\gamma }_k{X}_{eq,k})$$

(39)

$$f_4=-({\alpha }_k^2-2\sqrt{3}{\alpha }_k{\beta }_k-{\beta }_k^2)(P_{DG,k}^{3\varphi }+{\gamma }_k{R}_{eq,k})-(\sqrt{3}{\alpha }_k^2+2{\alpha }_k{\beta }_k-\sqrt{3}{\beta }_k^2)(Q_{DG,k}^{3\varphi }+{\gamma }_k{X}_{eq,k})$$

(40)

$$f_5=-(\sqrt{3}{\alpha }_k^2-2{\alpha }_k{\beta }_k-\sqrt{3}{\beta }_k^2)(P_{DG,k}^{3\varphi }+{\gamma }_k{R}_{eq,k})-({\alpha }_k^2+2\sqrt{3}{\alpha }_k{\beta }_k-{\beta }_k^2)(Q_{DG,k}^{3\varphi }+{\gamma }_k{X}_{eq,k})$$

(41)

$$f_6=-({\alpha }_k^2+2\sqrt{3}{\alpha }_k{\beta }_k-{\beta }_k^2)(P_{DG,k}^{3\varphi }+{\gamma }_k{R}_{eq,k})+(\sqrt{3}{\alpha }_k^2-2{\alpha }_k{\beta }_k-\sqrt{3}{\beta }_k^2)(Q_{DG,k}^{3\varphi }+{\gamma }_k{X}_{eq,k})$$

(42)

$$h_1=\frac{-{\alpha }_k{X}_{eq,k}+{\beta }_k{R}_{eq,k}}{({\alpha }_k^2+{\beta }_k^2)(R_{eq,k}^2+X_{eq,k}^2)}$$

(43)

$$h_2=\frac{{\alpha }_k{R}_{eq,k}+{\beta }_k{X}_{eq,k}}{({\alpha }_k^2+{\beta }_k^2)(R_{eq,k}^2+X_{eq,k}^2)}$$

(44)

The nonzero elements of the submatrix B corresponding to bus k are calculated from the partial derivatives of Equations (14) and (15) with respect to the three-phase reactive power output, as follows:

$$\begin{array}{ll}{B}_k& ={\left[\begin{array}{cccccc}\frac{\partial \Delta I_{m,k}^a}{\partial Q_{DG,k}^{3\varphi }}& \frac{\partial \Delta I_{r,k}^a}{\partial Q_{DG,k}^{3\varphi }}& \frac{\partial \Delta I_{m,k}^b}{\partial Q_{DG,k}^{3\varphi }}& \frac{\partial \Delta I_{r,k}^b}{\partial Q_{DG,k}^{3\varphi }}& \frac{\partial \Delta I_{m,k}^c}{\partial Q_{DG,k}^{3\varphi }}& \frac{\partial \Delta I_{r,k}^c}{\partial Q_{DG,k}^{3\varphi }}\end{array}\right]}^t\\ & =\frac{1}{2({\alpha }_k^2+{\beta }_k^2)}{\left[\begin{array}{cccccc}2{\alpha }_k& -2{\beta }_k& -{\alpha }_k+\sqrt{3}{\beta }_k& \sqrt{3}{\alpha }_k+{\beta }_k& -{\alpha }_k-\sqrt{3}{\beta }_k& -\sqrt{3}{\alpha }_k+{\beta }_k\end{array}\right]}^t\end{array}$$

(45)

Similarly, the nonzero elements of the submatrix C for bus k are given by:

$$\begin{array}{ll}{C}_k& =\left[\begin{array}{cccccc}\frac{\partial \Delta V_k^{+}}{\partial V_{r,k}^a}& \frac{\partial \Delta V_k^{+}}{\partial V_{m,k}^a}& \frac{\partial \Delta V_k^{+}}{\partial V_{r,k}^b}& \frac{\partial \Delta V_k^{+}}{\partial V_{m,k}^b}& \frac{\partial \Delta V_k^{+}}{\partial V_{r,k}^c}& \frac{\partial \Delta V_k^{+}}{\partial V_{m,k}^c}\end{array}\right]\\ & =\frac{1}{6\sqrt{{\alpha }_k^2+{\beta }_k^2}}\left[\begin{array}{cccccc}2{\alpha }_k& 2{\beta }_k& -{\alpha }_k+\sqrt{3}{\beta }_k& -\sqrt{3}{\alpha }_k-{\beta }_k& -{\alpha }_k-\sqrt{3}{\beta }_k& \sqrt{3}{\alpha }_k-{\beta }_k\end{array}\right]\end{array}$$

(46)

3.4. Representation of the Reactive Power Limit

In general, the reactive power limit of a vector-controlled DG is represented by the q-axis current limit, instead of being expressed in terms of the three-phase reactive power output itself [25,26]. Figure 3 illustrates the current vector of the BTCS and the voltage vector of a DG-connected bus in a rotating reference frame.

Figure 3. Voltage and current vectors in a rotating reference frame.

Figure 3. Voltage and current vectors in a rotating reference frame.

As discussed in Section 2, the phase voltages are normally unbalanced, while the currents are balanced. Thus, the difference of angle between the voltage and current vectors, θ, will change continuously. In power flow studies, the average value of the angle (θ_av), given by Equation (47), can be used. For the voltage angle, the angle of the positive-sequence voltage will be used, since the negative- and zero-sequence terms have zero average.

$${\theta }_{av}=\arg \left(I_{BTCS}^a\right)-\arg \left(V_k^{+}\right)$$

(47)

Since the current output of the BTCS is balanced, the magnitude of the current vector, |i_BTCS|, is equal to the peak value of the phase current. Therefore, the d- and q-axis components of i_BTCS can be expressed in terms of θ_av and the magnitude of the current phasor of phase a, |I^a_BTCS|, as follows:

$$i_{BTCS,d}=\sqrt{2}|I_{BTCS}^a|\cos ({\theta }_{av})$$

(48)

$$i_{BTCS,q}=\sqrt{2}|I_{BTCS}^a|\sin ({\theta }_{av})$$

(49)

When the q-axis current is bounded by a specific limit, i_q,limit, as shown in Figure 3, the magnitude of the limited current vector, |i_BTCS,limited|, and the angular change due to the limitation, φ, are obtained from:

$$|i_{BTCS,limited}|=\sqrt{i_{BTCS,d}^2+i_{q,limit}^2}$$

(50)

$$\phi ={\tan }^{-1}\left(\frac{i_{BTCS,q}}{i_{BTCS,d}}\right)-{\tan }^{-1}\left(\frac{i_{q,limit}}{i_{BTCS,d}}\right)$$

(51)

Because the limited current is still balanced, the angular change of the current phasor is also φ. Therefore, the current phasor of phase a is given by:

$$I_{BTCS,limited}^a=\frac{|i_{BTCS,limited}|e^{j(\arg (I_{BTCS}^a)-\phi )}}{\sqrt{2}}$$

(52)

By substituting Equation (52) into Equation (8) and taking the imaginary part, the three-phase reactive power output of the DG with limited q-axis current can be calculated from:

$$Q_{DG,limited}^{3\varphi }=\mathrm{Im}\left\{3\frac{R_{d,k}-j{X}_{C,k}}{R_{eq,k}-j{X}_{eq,k}}{V}_k^{+}{(I_{BTCS,limited}^a)}^{*}-\frac{|V_k^a{|}^2+|V_k^b{|}^2+|V_k^c{|}^2}{R_{eq,k}-j{X}_{eq,k}}\right\}$$

(53)

If the reactive power of a voltage-control-mode DG is not limited, the DG-connected bus is modeled as a PV bus. However, if the q-axis current of a DG is restricted by its limit, the bus model should be changed to a bounded PQ-bus model, meaning that the reactive power of the PQ bus is not constant, and is calculated via Equation (53) for every iteration. The transition between a PV-bus model and a bounded PQ-bus model is determined by i_BTCS,q calculated via Equation (49) and ΔV_k⁺ calculated via Equation (25), as shown in Figure 4. To avoid the numerical instability caused by frequent changes between PV-bus and PQ-bus models, the bus model is only changed if the transition condition occurs in two consecutive iterations. Note that in the frame of reference shown in Figure 3, a positive q-axis current represents negative reactive power (i.e., reactive power absorption), and vice versa. Accordingly, a lower-bounded PQ-bus model (i_q,limit = i_q,min) indicates that the DG connected to that bus supplies the maximum reactive power. Therefore, if the voltage is higher than the reference value (i.e., ΔV_k⁺ > 0), the bus model should be changed to a PV-bus model. For an upper-bounded PQ-bus model (i_q,limit = i_q,max), the opposite condition is applied.

Figure 4. Transition between a PV-bus model and bounded PQ-bus models.

Figure 4. Transition between a PV-bus model and bounded PQ-bus models.

3.5. Power Flow Procedure

Step 1:

Determine the bus model for the voltage-control-mode DG-connected buses, based on the rules shown in Figure 4;

Step 2:

Calculate the three-phase reactive power outputs of the bounded PQ buses, using Equation (53);

Step 3:

Calculate the current and voltage mismatches, using Equations (14), (15) and (25);

Step 4:

Test the convergence: if the absolute values of all mismatches are within the convergence tolerance, terminate the power flow; otherwise, go to Step 5;

Step 5:

Calculate the Jacobian matrix;

Step 6:

Update the state variable vector using Equation (31), and go to Step 1.

3.6. Modifications for Type 1 DG

For a Type 1 DG, some of the foregoing equations are modified.

(1)

The second terms of Equations (16)–(21) are eliminated;

(2)

γ_k in Equationa (16)−(21) and Equations (37)–(42) is set equal to zero;

(3)

The last two matrices of Equation (36) are eliminated; and

(4)

Equation (53) is modified to

$$Q_{DG,limited}^{3\varphi }=3\times \mathrm{Im}\left\{V_k^{+}{(I_{BTCS,limited}^a)}^{*}\right\}$$

(54)

4. Case Studies and Results

To verify its accuracy and convergence performance, the proposed method was implemented using MATLAB. For all mismatches, the convergence tolerance was set at 1.0 × 10⁻⁶ pu on a system base of 100 kVA. The Type 3 DG, which is the most complex type, was used in all simulations.

4.1. Verification of Accuracy

The accuracy of the proposed power flow method was verified by comparing its results to those obtained from the PSCAD model. In addition, power flow studies using conventional simplified DG model [5,6,7,8,9] were conducted to show the improvement in accuracy by the proposed method. In the simplified DG model, the three-phase power outputs of the DG are equally divided into each phase, regardless of the voltage of the DG-connected bus.

Figure 5 shows the modified IEEE 13-bus test system, with a 1-MW DG connected to bus 680 through a transformer [27]. The parameters of the DG are summarized in Table 2. The DG can supply three-phase reactive power in the range −500 kVAr to 600 kVAr at the rated voltage. In the PSCAD model, the DG was modeled with a three-wire two-level VSC using IGBT switches. For the switch control scheme, the sinusoidal pulse width modulation (SPWM) method was adopted. The simulation was performed with time step of 1 μs.

Figure 5. Modified IEEE 13-bus test system.

Figure 5. Modified IEEE 13-bus test system.

Table 2. Parameters of the DG.

Property	Value	Unit
Maximum q-axis current, iq,max	650	A
Minimum q-axis current, iq,min	−650	A
System frequency, fsys	60	Hz
Switching frequency, fsw	10	kHz
Converter side inductance, L1	0.25	mH
Grid side inductance, L2	0.11	mH
Capacitance, C	279	μF
Damping resistance, Rd	0.17	Ω

Table 2. Parameters of the DG.

Property	Value	Unit
Maximum q-axis current, iq,max	650	A
Minimum q-axis current, iq,min	−650	A
System frequency, fsys	60	Hz
Switching frequency, fsw	10	kHz
Converter side inductance, L1	0.25	mH
Grid side inductance, L2	0.11	mH
Capacitance, C	279	μF
Damping resistance, Rd	0.17	Ω

4.1.1. Reactive Power Control-Mode DG

In this simulation, the DG controlled the active and reactive power outputs in accordance with the reference values. We tested five different three-phase active power outputs (200 kW, 400 kW, 600 kW, 800 kW and 1000 kW), with the DG operating at a unity power factor (i.e., the reactive power output was set equal to zero).

Figure 6 shows the active and reactive power outputs of each phase, measured in the PSCAD model. In all cases, the active power outputs of the phases differed because the system line and load conditions were unbalanced. Furthermore, although the three-phase reactive power (sum of the three phases) was zero, one phase supplied reactive power while the others absorbed it. In addition, the mismatches between phases increased as the DG active power output increased. When the active power output was 1000 kW, the largest active and reactive power mismatches were 22.16 kW and 18.81 kVAr, respectively.

Figure 6. Phase power outputs measured in the PSCAD model: (a) active power and (b) reactive power.

Figure 6. Phase power outputs measured in the PSCAD model: (a) active power and (b) reactive power.

The simulation using the proposed power flow produced almost the same results as the PSCAD model. The maximum phase output power difference between the two models was less than 0.13 kVA. The maximum differences in the phase voltage magnitude and angle were less than 2.0 × 10⁻⁵ pu (0.002%) and 0.023 degrees, respectively, as shown in Figure 7.

Figure 7. Maximum difference between the power flow results with the proposed DG model and the PSCAD simulation results: (a) voltage magnitude and (b) voltage angle.

Figure 7. Maximum difference between the power flow results with the proposed DG model and the PSCAD simulation results: (a) voltage magnitude and (b) voltage angle.

In all cases, the power flow converged in five iterations, which is the same as the number of iterations required for the original system. In other words, adding a DG did not increase the number of iterations in this scenario.

Figure 8 shows the comparison between the PSCAD model and the power flow using simplified DG model. The maximum phase voltage magnitude and angle differences were 2.7 × 10⁻³ pu (0.27%) and 0.112 degrees, respectively, which are significantly larger than the error of the proposed DG model. Moreover, phase a and c voltages calculated from the power flow were higher than those obtained from the PSCAD model, while the phase b voltages were reversed, as shown in Figure 9. From the above results, it is concluded that the proposed method can more accurately represent the output power characteristics of DG under unbalanced operating conditions.

Figure 8. Maximum difference between the power flow results with the simplified DG model and the PSCAD simulation results: (a) voltage magnitude and (b) voltage angle.

Figure 8. Maximum difference between the power flow results with the simplified DG model and the PSCAD simulation results: (a) voltage magnitude and (b) voltage angle.

Figure 9. Voltage profile of the nodes on the path from the secondary side of the voltage regulator to the DG-connected node when the active power output of the DG was 1000 kW.

Figure 9. Voltage profile of the nodes on the path from the secondary side of the voltage regulator to the DG-connected node when the active power output of the DG was 1000 kW.

4.1.2. Voltage Control-Mode DG

In this simulation, the DG controlled the positive-sequence voltage of bus 681 instead of the reactive power output. The three-phase active power output of the DG was set at 400 kW, and three cases with different reference voltages (0.98 pu, 1.01 pu, and 1.04 pu) were tested. Figure 10a compares the reference voltage, positive-sequence voltage, and phase voltages measured in the PSCAD model. The magnitudes of the phase voltages differed appreciably from each other, due to the unbalanced operating conditions. In Case 2 (reference voltage = 1.01 pu), the DG was able to control the voltage of bus 681 at the reference value. However, in the other two cases, the voltage could not be controlled, due to the q-axis current limit. In Case 1 (reference voltage = 0.98 pu), the three-phase reactive power output of the DG was limited to −497 kVAr, and the positive-sequence voltage was 0.991 pu, since the q-axis current was restricted by the upper limit. In Case 3 (reference voltage = 1.04 pu), the q-axis current was restricted by the lower limit, and the positive-sequence voltage and reactive power output of the DG were 1.036 pu and 626 kVAr, respectively.

The results of proposed power flow method matched well with the PSCAD simulation results. The maximum three-phase reactive power difference between the two methods was 0.25 kVAr, which is 0.05% of the rated reactive power limit. The differences in voltage magnitude were also very small. The maximum differences between the two methods were less than 3.0 × 10⁻⁵ pu (or 0.003%), as shown in Figure 10b. Six iterations were required for the first case, and five iterations for each of the other cases.

Figure 10. Results for the voltage control-mode DG: (a) reference voltage and voltages measured in the PSCAD model and (b) maximum voltage magnitude difference.

Figure 10. Results for the voltage control-mode DG: (a) reference voltage and voltages measured in the PSCAD model and (b) maximum voltage magnitude difference.

4.2. Application to a Large Distribution System

To investigate the effect of the number of DGs on the convergence rate of the proposed power flow method, we used the IEEE 123-bus test system [27]. The power flow problem was solved repeatedly, adding DGs to the three-phase buses one at a time, and observing the effect on the required number of iterations. It was assumed that each DG was connected to a bus through a 150 kVA 4.16/0.69 kV transformer, and generated 60 kW of active power. The main parameters of the DG are summarized in Table 3. In order to analyze the effects of operating mode and current limit on the convergence rate, we simulated the following three cases.

Case 1:

Each DG operates in reactive power control mode. The reactive power reference value is set at 0 kVar. The DG-connected buses are modeled as PQ buses;

Case 2:

Each DG operates in voltage control mode. The reference voltage is 1.03 pu. A q-axis current limit is not imposed (i.e., the current is not bounded). Therefore, the DG-connected buses are always modeled as PV buses;

Case 3:

The operating mode and reference voltage are the same as in Case 2, but a q-axis current limit is included. The DG-connected buses can be modeled as either PV buses or bounded PQ buses, depending on the q-axis current and voltage mismatch.

Table 3. Parameters of the DG.

Property	Value	Unit
AC voltage, Vac	690	V
Maximum q-axis current, iq,max	100	A
Minimum q-axis current, iq,min	−100	A
Grid side inductance, L2	0.3	mH
Capacitance, C	27.8	μF
Damping resistance, Rd	0.9	Ω

Table 3. Parameters of the DG.

Property	Value	Unit
AC voltage, Vac	690	V
Maximum q-axis current, iq,max	100	A
Minimum q-axis current, iq,min	−100	A
Grid side inductance, L2	0.3	mH
Capacitance, C	27.8	μF
Damping resistance, Rd	0.9	Ω

Figure 11 shows the number of iterations required for convergence in each of the three cases, according to the number of DGs. In Cases 1 and 2, in which there was no current limit, the power flow converged in at most six iterations, regardless of the number of DGs. In other words, even when as many as 62 DGs were added, only one additional iteration was required for convergence. In Case 3, the average iteration number was 7.8 and the maximum number of iterations was 11, when the number of installed DGs was 51. In this case, the q-axis current was bounded in 30 of the 51 DGs, as shown in Figure 12. Thus, the convergence rate of the proposed power flow method was not substantially degraded by the addition of DG models, even when q-axis current limits were imposed on many of the DGs.

Figure 11. Number of iterations required for convergence.

Figure 11. Number of iterations required for convergence.

Figure 12. Numbers of PV buses and bounded PQ buses for Case 3.

Figure 12. Numbers of PV buses and bounded PQ buses for Case 3.

5. Conclusions

We proposed a TCIM-based power flow model for vector-controlled DGs. Three kinds of equivalent circuits for DGs were considered, classified according to output filter topology and current control target. Based on these equivalent circuits, steady-state phase current output models were formulated in terms of the three-phase active and reactive power outputs of the DG, phase voltages of the DG-connected bus, and impedances of the output filter. Finally, a method was proposed for implementing a TCIM power flow including the DG models. In the proposed power flow, a DG-connected bus was modeled as either a PQ bus or a PV bus. The active and reactive power and the voltage of a DG-connected bus were represented by three-phase quantities. The q-axis current limit was used to represent the reactive power limit of a voltage-control-mode DG. If the q-axis current was not bounded, a DG-connected bus was modeled as a PV bus. Otherwise, the bus was modeled as an upper-bounded or lower-bounded PQ bus. Rules for changing the bus model were also presented. In the case studies, the accuracy of the proposed method was verified by comparing the results with those of the PSCAD model. When the proposed power flow method was used, the calculated voltage magnitudes for the test system differed from those of the PSCAD model by less than 3.0 × 10⁻⁵ pu (0.003%). In addition, the effect of the number of DGs on the convergence rate of the power flow was analyzed by increasing the number of DGs in the IEEE 123-bus test system. It was proved that the convergence rate of the proposed power flow method was not substantially degraded by the addition of DG models, even when the q-axis current was limited in many of the DGs.