Research on random fuzzy power flow calculation of AC/DC hybrid distribution network based on unified iterative method

Correspondence Ping Dong, College of Electric Power, South China University of Technology, Guangzhou, 510640, China Email: epdping@scut.edu.cn Abstract In recent years, with the promotion of policies and technologies, a large number of distributed generations are connected to the urban distribution network in various forms. At the same time, the flexible DC transmission technology has become mature and gradually applied to the distribution network. In this condition, the AC/DC hybrid distribution networks with flexible regulation ability will be the trend of urban distribution network. This paper accurately considers the loss of voltage source converter, and proposes a unified iterative power flow method for AC/DC hybrid distribution network, which takes the voltage source converter control modes into account. Because of the uncertainties of the distributed generations and loads, a random fuzzy power flow model is proposed, in which cumulant method is applied in the random stage. In the fuzzy stage, the fuzzy simulation technology is applied. Finally, an improved PG&E 69-Bus system is taken as an example to verify the effectiveness of the unified iterative method and the random fuzzy power flow calculation, and put forward some suggestions for the operation of distribution network with distributed generations.


INTRODUCTION
The types and the amount of the power loads increased recently with the higher penetration of distributed generations (DGs) such as photovoltaic and wind power generation in the distribution network. In addition, with the development of power electronic technologies, the performance of voltage source converter (VSC) has been greatly improved, and it has been gradually applied in the distribution network. Therefore, AC/DC hybrid distribution network based on VSC becomes the trend of the urban distribution network in the future. Compared with the traditional AC distribution network, the AC/DC hybrid distribution network has the functions of controlling power flow and optimizing network operation. The access of VSC changed the calculation method of deterministic power flow in the distribution network.
In the power flow solution of AC/DC hybrid network, according to the coupling relationships between AC and DC, AC/DC power flow can be divided into two categories: alternating iterative method (AIM) and unified iterative method (UIM).
The AIM means that the AC/DC hybrid system is divided into AC and DC subsystems, which are decoupled and equivalent at the AC and DC sides according to the different control modes of VSC. Information is exchanged between AC and DC subsystems through converter equations, and the system state including voltages and angles is obtained when all subsystems converge.
In [1], the AIM is adopted, but the VSC loss is ignored. Taking the VSC loss into account, Sun et al. [2] and Zheng et al. [3] use the AIM to calculate the power flow, however, the former thinks that the loss accounts for a certain proportion of transmission power, and the latter equates the loss with resistance to simplify calculation. In fact, the loss will affect the accuracy of power flow calculation results so the loss of VSC should be modelled accurately. The AIM has a complex calculation process and poor convergence. In [4], the AIM is improved, that is, the converter is incorporated into the AC side, which simplifies the calculation process, but does not improve the convergence performance of VSC under the condition of constant AC side voltage.
The UIM takes the AC/DC hybrid system as a unified whole, which writes the balance equations according to the types of AC and DC buses and VSC control modes, and solves it with a Jacobian matrix. It has the characteristic of strong convergence. In some literature, the UIM is used to calculate the AC/DC power flow. In [5], the UIM is involved without the consideration of the VSC control modes. Even using the UIM, the redundant variables make the Jacobian matrix complex which burdens the calculation [6]. And the VSC losses are not accurately considered in [7][8][9].
When it comes to renewable energy sources, the deterministic power flow cannot describe the state of AC/DC distribution network due to the stochastic characteristic. In order to discuss the influence of the probabilistic uncertainties, in 1974 borkowska put forward the concept of probabilistic power flow (PPF) [10]. Since then, the theory of PPF has gradually enriched.
In general, there are several types of PPF methods: analytical methods, approximate methods and Monte Carlo simulation (MCS) [11,12]. The cumulant method (CM) [13][14][15][16] represents the analytical method. And the point estimate method (PEM) [17][18][19][20][21] and the unscented transformation (UT) represent the typical approximate methods. The MCS deals uncertainties with numerous deterministic power flow to obtain accurate probabilistic states, which takes a long time in computation. Generally, it is recognized as the standard reference for other methods. Compared to MCS, CM and PEM have less calculation time at the cost of accuracy loss, which is usually involved in many literature.
There have been many studies on the PPF in the AC distribution network, but research works on PPF in AC/DC network are rare. In [22], MCS is used to calculate the PPF in AC/DC hybrid network. The UT is used in [23]. As an improved sampling method, Latin hypercube sampling is developed in [24]. In [25], PEM is employed. On account of the VSC, there are few research works adopting CM to calculate the PPF in the hybrid network.
PPF with deterministic distribution parameters is unable to fully describe the uncertainties of renewable energy sources. Literature focusing on multiple uncertainties of renewable energy sources is rare. In [26], wind speed model is established considering randomness and fuzziness. In [27], by introducing fuzzy theory, random fuzzy power flow (RFPF) is put forward in the AC distribution network, which is feasible to analyze two uncertainties.
This paper proposes a unified iterative RFPF calculation method for AC/DC hybrid distribution network, which takes accurately the VSC loss into account. First, the unified iterative model is established considering different VSC control modes, then the double uncertainty of DGs and loads in the distribution network is modelled, and the RFPF calculation based on CM and fuzzy simulation technology is developed. In addition, since multi wind farms are connected to one power grid, the Cholesky decomposition is adopted to address the correlation of the variable of different wind farms. Finally, the PG&E 69-Bus system is taken as an example to verify the effectiveness of the UIM and the RFPF calculation.

Model of VSC
The AC/DC hybrid distribution network is mainly composed of AC network, DC network and VSC. VSC is not only a power interaction device for DC network and AC network, but also a key device for controlling system power flow. The equivalent circuit is shown in Figure 1.
Regardless of the filter, the AC network is connected to the DC network through a power transformer, a phase reactor and a VSC in the condition of steady state. The transmission power of the AC side can be expressed as where P s and Q s are the active power and reactive power injected into the VSC by the AC network; U s and δ s are the voltage and phase angle of the AC network; P c and Q c are the active power and reactive power at the AC side of the VSC; U c and δ c are the voltage and phase angle at the AC side of the VSC; R and X are equivalent resistance and reactance of transformer and phase reactor; U d and I d are DC voltage and current, and δ = δ s − δ c , α = arctan(X/R), Y = 1∕ √ R 2 + X 2 . The power and voltage on both sides of the VSC meet the following constraints: where μ is the utilization rate of DC voltage whose value is usually √ 3∕2; M is the modulation index of the VSC, 0 ≤ M ≤ 1. P loss,VSC is the VSC loss.
The VSC loss can be expressed as where I c is the current at the AC side of the VSC; A is the fixed loss of VSC, B is the coefficient of the loss in direct proportion to I c , C is the coefficient of the loss in direct proportion to I c 2 . A, B and C are given by: where S N is the rated power of VSC and U Bdc is the rated voltage of DC side. The subscripts r and i represent the rectifier and inverter states respectively.

Unified iterative power flow equations
VSC can control the active and reactive physical quantities through M and δ. The active physical quantities include the active power P s at the AC side, the voltage U d at the DC side, the frequency f and the current I d at the DC side. The reactive physical quantities include the reactive power Q s at the AC side and the voltage U s at the AC side. Generally, VSC has four control modes, such as U dc −Q s constant mode, P s −Q s constant mode, U dc −U s constant mode and P s −U s constant mode. The control modes will affect the balance equations of AC/DC network, which should be analyzed in detail. The balance equations of the buses not directly connected with VSC in the network are the same as those in the common AC network and the common DC network, which will not be discussed here.

U dc −Q s constant mode
Suppose the VSC number is m (1 ≤ m ≤ N VSC ), N VSC is the number of VSC in the network. The subscript sm represents the AC bus of the AC network, the subscript cm represents the AC bus of the VSC, and the subscript dcm represents the DC bus of VSC. The equivalent injected reactive power of the AC bus is known, whose value is -Q VSCm , and Q VSCm is the reactive power controlled by VSC. However, the equivalent active power is unknown. According to (1), the balance equations of the AC bus of VSC can be expressed as follows: .
(6) And the balance equations of VSC are given by In (6) and (7), j(1≤j≤N AC ) represents the AC bus number, i(1≤i≤N DC ) represents the DC bus number, P loss,VSCm is the VSC loss, N AC and N DC are the number of buses in the AC and DC networks, respectively. Since the DC bus of VSC is a DC slack bus, there is no need to write the balance equations of the bus.
It can be seen that for the various control modes of VSC, as long as the balance equations of AC and DC buses connected with them are properly modified, and the corresponding balance equations are added according to the control mode, the numbers of balance equations and variables in the network are equal, and the network can be solved. The balance equations of the remaining three control modes are obtained in a similar way.

P s −Q s constant mode
AC bus: . (8) DC bus: Balance equations of VSC:

U dc -U s constant mode
AC bus: Balance equations of VSC:

P s -U s constant mode
AC bus: DC bus: Balance equations of VSC: .
(15) The number of variables of AC/DC buses is N, and there are 2N VSC variables of VSC in total. Newton Raphson method is used to solve the network, the Jacobian matrix is J, and the order is N + 2N VSC , then the modified equation is as follows: After the equations are solved, according to (3), we have It means the AC bus voltage and DC bus voltage of VSC should meet this constraint otherwise the power flow should be recalculated.
In the above-mentioned AC/DC power flow calculation model based on VSC control mode, the modulation index M and DC current I dc are used as intermediate variables, so that the number of solving variables decreased, which is suitable for large-scale distribution network power flow calculation.

Random fuzzy models of wind power, solar power and load
In the past, wind speed and solar radiation are recognized as probabilistic variables, so the probabilistic theory is involved in power flow. But actually, other uncertainties should also be considered. Through simulation, in [25], the parameters of Weibull function of wind speed are found fuzzy, which means fuzzy and random characteristic coexist in renewable energy sources.
A fuzzy variable ξ is the function from the possibility space (Θ, P(Θ), Pos) for the real line R. By introducing fuzzy theory, the random fuzzy variable can be defined as a function from the possibility space (Θ, P(Θ), Pos) for the set of random variables. That means, for any θ ∈ P(Θ), ξ(θ) is a random variable. So, in the following, by applying fuzzy numbers to the parameters of the probabilistic distribution function, we can get the random fuzzy models.

Wind power
This paper uses three-parameter Weibull distribution to simulate the distribution of wind speed. For Weibull distribution, the parameters are represented by fuzzy numbers, then the wind speed becomes a random fuzzy variable, and the expression is as follows: where v is the wind speed; v 0 is the position parameter, and ξ k and ξ c are the shape and scale parameter, respectively. The wind power output function at speed v is expressed as where P w is the output of the wind turbine; P wN is the rated power of the wind turbine; v is the wind speed; and v ci , v r , v co are the cut-in, rated, and cut-out wind speeds respectively. The values of k 1 and k 2 are expressed as According to statistics, wind speed is mostly between v ci and v r, so the function relationship between wind power output and wind speed is approximately linear. Then the random fuzzy output of the wind turbine is expressed as The wind turbine usually absorbs the reactive power, which can be calculated by the active power and power factor.

Solar power
For photovoltaic output, the model of light intensity should be studied first. It is generally considered that the light intensity conforms to beta distribution, and the expression is where r and r max represent the actual and maximum values of light intensity, respectively, ξ a and ξ b are shape parameters of light intensity distribution function and Γ represents gamma function. Given the total area of the solar modules A and efficiency η, the total active power P PE is evaluated by the following function: Similarly, the random fuzzy model of a solar power output is expressed as where P max is the maximum power of the solar power.

Load
The load power can usually be represented by the normal distribution, and their random fuzzy models can be formulated by where P Load and Q Load are the active and reactive load powers, respectively; ξ μP and ξ μQ are the means of the active power and reactive power, respectively; ξ σP and ξ σQ are standard deviations of the active power and reactive power, respectively.

Random fuzzy power flow method
Based on the linearization of power flow equations, CM, Gram-Charlier series expansion and fuzzy simulation technology, the RFPF is conducted in this paper.

Cumulant Method and linearized power flow equation
The cumulant γ k of a random variable can be obtained by where β k is the k-order central moment of the random variable. The matrix form of the linearized power flow equation is where ΔW is the increment of bus injection power, ΔX is the increment of bus voltage and phase angle and J 0 is the Jacobian matrix at the operating point. It can be seen from the additivity of cumulant that if there are several random variables in the injected power of a bus, ΔW (k) of the injected power increment of the bus is equal to the sum of the cumulants of each random variable, that is For the bus connected with VSC, its power cumulant is different from that of the ordinary bus. Take the AC bus in the AC network connected to a VSC with U dc -Q s constant mode as an example, the equivalent injection active power of this bus is Suppose there are N buses in the DC part of the VSC, and the injection power of each bus are P iDC,1 , P iDC,2 …P iDC,N , we have P is = P iDC,1 + P iDC,2 + ⋯ + P iDC,N − P loss,VSC − P loss,DC (30) where P loss,VSC , P loss,DC are the loss of VSC and line loss of DC network, respectively. Considering their small values, the cumulants of injected power of the AC bus of VSC can be approximately expressed as The cumulants of the buses connected to the other control modes VSC and the cumulants related to the VSC balance equations can be obtained by referring to the above process.
According to the linear properties of cumulant, the cumulants of output variables are obtained by

Expansion of Gram-Charlier series
According to the theory of Gram-Charlier series, if the μ and σ of random variable x are known, the formula z = (x -μ)/σ can be used to standardize it. The probability density function f (z) can be expressed as where H i (z) is Hermite polynomials of order i. This paper takes the first seven orders to calculate f (z). After obtaining f (z), the PDF of x can be obtained according to x = zσ + μ as follows:

Random fuzzy power flow calculation process
The complete calculation process is shown in Figure 2.

CASE STUDIES
The method proposed in this paper was applied to the improved PG&E 69-Bus system. The improved PG&E 69-Bus system is shown in Figure 3. The DC network operates in a unipolar circuit that has two DC lines. The resistance of the modified DC lines is the same as the original AC lines, and the other parameters of the system are shown in Table 1. This paper uses WTs and PVs are used; the basic parameters of wind speed density distribution function are c = 10.7, k = 4, v 0 = 3 m/s. The cut-in, rated and cut-out speeds of the WTs are 4, 15 and 25 m/s, respectively, and the power factor is 0.9, absorbing reactive power. The basic parameters of the solar irradiance density distribution function are a = 1.5, b = 1. The basic values of load standard deviations are 0.1 times of the mean values. The network structure is shown in Figure 3.

Model validation of different fuzzy numbers
There are many forms of fuzzy numbers in fuzzy theory, such as triangular fuzzy numbers and trapezoidal fuzzy numbers. In order to illustrate the impact of the fuzzy numbers, we use different fuzzy numbers in the models, as shown in Table 2. The interval (k 1 , k 2 , k 3 ) for triangular fuzzy number means the three values of all the triangular fuzzy numbers are k 1 , k 2 and k 3 times of their base values, respectively. And so are the trapezoidal fuzzy numbers.
There are two PVs in buses 35 and 52, and one WT generator in bus 15. The rated power of all the DGs is 0.6 MW, and the total penetration is 48.4%. The mean values and standard deviations of voltage and phase angles of network nodes under different conditions are extracted, as shown in Figure 4.
From Figure 4 we can see that, given the same range of fluctuations, the form of the fuzzy numbers has no effect on the result. For the sake of simplicity, the triangular fuzzy numbers are adopted in our models.

Comparison between pure AC network and AC/DC network
In order to verify the advantages of AC/DC hybrid distribution network, this example compares the voltages and network loss of AC/DC distribution network and AC distribution network. Disconnect branches 35-39 and 27-54 in Figure 3. All VSCs in  Figure 5 (the numbers in brackets are the ratio of network loss to load). It can be seen from Figure 5 that due to the functions of VSC to control power flow and voltage, not only the voltages of the modified system are improved, but also the total loss of the network is reduced, which proves the advantages of the operation of the AC/DC hybrid power grid.

Comparison of two iterative methods
In order to verify the effectiveness of the UIM proposed in this paper, without adding DGs, the UIM and the AIM are used to calculate the network in Figure 3. The control modes and parameters of each VSC are shown in Table 3, and the results are shown in Figure 6. It can be seen from Figure 6 that the results of the two iterative methods are almost the same, the maximum deviation of voltages is 5 × 10 -5 p.u., and the maximum deviation of phase angles is 1.08 × 10 -3• , which proves the effectiveness of the UIM proposed in this paper. In terms of the number of iterations, AIM needs more total number of iterations to achieve convergence of all the subsystems. Compared to the UIM with the same calculation accuracy, AIM needs more iteration. In addition, the number of iterations of the UIM will increase approximately linearly with the accuracy.

The comparison between three PPF methods
When the fluctuation intervals of all fuzzy parameters in RFPF calculation are set to (1,1,1), it is actually simplified into PPF based on CM and Gram-Charlier series expansion. As a comparison, CM, MCS and three-point estimate method (3-PEM) are adopted to calculate PPF. The results are shown in Table 4, where μ V and μ δ represent the mean values of voltages and phase angles of all buses, σ V and σ δ represent the standard deviations of voltages and phase angles, ε mean and ε max represent the mean and maximum relative errors, respectively.
From Table 4, we can see the CM have little difference from the other two PPF methods. Because of the short calculation time of CM, we think it is an effective method to apply in AC/DC network.

4.5
The influence of the penetration and fluctuation of DGs on the network 1) There are three DGs in buses 15, 35 and 52. And there are different DGs in different cases as shown in Table 5.
All the intervals for triangular fuzzy numbers are set to (0.8,1,1.2). The mean values and standard deviations of voltages and phase angles under different cases are shown in Figures 7 and 8. And the voltage PDF of bus 15 is shown in Figure 9. The comparison between RFPF and PPF is shown in Table 6. In Table 6, we calculate RFPF in the condition of small fluctuation to study the compatibility between RFPF and PPF under different DG penetration (case A-D) and different DG types (cases E-G). And the mean relative errors of all cases in Table 6 seems to be relatively small, which states the RFPF and PPF are perfectly compatible. The RFPF method is able to handle the fuzziness and randomness.  Tables 9 and 10.
For bus 15, when V ≥ 1.07 or V ≤ 0.97, it is considered that the bus voltage exceeds the limit, as shown in Figure 12.
When we fix the DG penetration, we need to calculate only once PPF for Tables 8 and 10 because the cases in Tables 8  and 10        there are some differences between Tables 9 and 10. The mean values in Table 9 are bigger than Table 10 for the corresponding cases, which indicates the fluctuations of DG power have a greater influence on the global states of the network than that of loads. From the point of view of maximum values in Tables 9  and 10, we see that the maximum relative errors in Table 10 are mostly bigger than those in Table 9, which means a high fluctuation of loads will have a greater impact on some buses of the network. Generally, the DGs should be concerned more than loads. From Figure 12, we know that when the penetration is high, the big fluctuation range of WT results in a high probability of node voltage exceeding limit.

The influence of different access schemes of DG on the network
This section explores the impact of different DG access schemes on the network. Four schemes are set as shown in Table 11. In the schemes, the fluctuation intervals of WTs, PVs and loads are set to (0.8,1,1.2), and the total penetration of DGs is fixed as 64.6%. The loads near the DG access bus are shown in Table 12. The PDF of voltage of buses 15, 35 and 52 under different schemes are shown in Figure 14. And the comparison between RFPF and PPF is shown in Table 13.
It can be seen from Figure 13 and Table 13 that given the same penetration, access schemes of multiple DGs are

Analysis for correlated wind sources
The above-mentioned section does not cover cases where variables are correlated. Actually, the correlation is common in the distribution network. There are many ways to describe correlation, one of which is the correlation coefficient matrix. Given a correlation coefficient matrix C W of several variables and marginal distribution function F, we can address the correlation in RFPF.

1)
Step 1: Generate sample Es of independent standard normal distribution variable E. Then transform C W to C Q using where ρ q,ij and ρ w,ij are members of C W and C Q , respectively. And D(ρ w,ij ) is the conversion factor obtained from (36). 2) Step 2: Get sample Qs for the normally distributed variable Q whose correlation coefficient matrix is C Q .
where G Q is the lower triangular matrix obtained by Cholesky decomposition of C Q .

3)
Step 3: Through the principle of equal probability conversion, the sample Ws of the input variable W with the correlation coefficient matrix C W can be generated where Φ is the cumulative distribution function of the standard normal distribution variables. where B = G −1 , and G is the lower triangular matrix obtained by Cholesky decomposition of C W .

5)
Step 5: Use the cumulant of Y to calculate the cumulant of W, and the cumulant of output variables will be obtained. So, the characteristic of correlation is handled.
In this section, there are two WTs in buses 15 and 17 whose rated powers are 1 MW. The fluctuation range of WT parameters is set to (0.9,1,1.1), and loads are set with no fluctuation. The correlation coefficient matrix of wind speed in two sites is given: ] .
Take bus 15 and branch 1-2 as the example to illustrate the influence of the correlation coefficient. Some results are given in Figures 14 and 15.
From Figures 14 and 15 show that the correlation coefficient of wind speeds ρ has no influence on mean voltage and mean active power of branch, while the standard deviation of voltage and branch active power have the linear relationship with ρ. That means a strong correlation may lead to wild fluctuation of voltage and branch power. So, when there are several wind sites with strong correlation, it is vital to concentrate more on the influence of the wind sites.