Three-Phase Unbalanced Optimal Power Flow Using Holomorphic Embedding Load Flow Method

Abstract

Distribution networks are typically unbalanced due to loads being unevenly distributed over the three phases and untransposed lines. Additionally, unbalance is further increased with high penetration of single-phased distributed generators. Load and optimal power flows, when applied to distribution networks, use models developed for transmission grids with limited modification. The performance of optimal power flow depends on external factors such as ambient temperature and irradiation, since they have strong influence on loads and distributed energy resources such as photo voltaic systems. To help mitigate the issues mentioned above, the authors present a novel class of optimal power flow algorithm which is applied to low-voltage distribution networks. It involves the use of a novel three-phase unbalanced holomorphic embedding load flow method in conjunction with a non-convex optimization method to obtain the optimal set-points based on a suitable objective function. This novel three-phase load flow method is benchmarked against the well-known power factory Newton-Raphson algorithm for various test networks. Mann-Whitney U test is performed for the voltage magnitude data generated by both methods and null hypothesis is accepted. A use case involving a real network in Austria and a method to generate optimal schedules for various controllable buses is provided.

1. Introduction

In recent years, with the integration of distributed generators, electric storage, electrical vehicles, and demand response units, the role of distribution systems is changing. Distributed energy units (DERs) are posing problems mainly in the low-voltage networks with their intermittency and uncontrollability. New innovative solutions are required to maintain grid security. Management of low-voltage distribution networks are challenging since they contain large array of devices which need to be controlled, and monitoring systems are limited. The above DERs along with loads should be run in a sustainable fashion since it is one of the biggest challenges. Various methods to control the DERs are presented in [1].

A so-called advanced distribution management system (ADMS) has come into existence, evolving from the transmission network’s supervisory control and data acquisition systems (SCADA). This is possible with the increase in smart meters and monitoring devices in the network which provides data acquisition abilities [2]. ADMS provides functionalities such as load flow analysis, optimal power flow, monitoring and control capabilities similar to SCADA systems [3]. This must in theory, host advance functionalities such as adaptive protection leading to self-healing, real-time monitoring, dynamic network reconfiguration and control [4,5,6,7,8]. This will provide intelligence to the grid with topology processor, state estimation, load and generation modeling [3]. The grid needs to be operated optimally and the power flows should be optimal to reduce losses, increase security, and maximize economic benefit. Energy balance should be maintained for secure operation of the network to maintain the frequency and voltage within its limits.

Optimal Power Flow (OPF) is one of the most fundamental functionalities of ADMS. In the literature, various OPF algorithms can be found. The authors in [9] describe an OPF algorithm to control active, reactive power, and transformer taps. The objective is to minimize system costs and losses. This method is based on Newton-Raphson load flow. Feasible power flow is solved, and the optimum is close to the load flow solution. Therefore, Jacobian information is used to calculate the optimum in a linear fashion. In [10], a non-linear programming technique is used to provide solution to OPF problem the objectives being economic dispatch and generation cost minimization. Same as before, load flow is performed to determine a feasible solution. Fletcher-Powell method is used to minimize the objective function. A general economic dispatch problem is implemented in [11]. This approach is similar to that in [9,10]. In [12], an OPF method for power system planning is provided. It used generalized reduced gradient technique to find the optimum. Hessian Matrix-based OPF method is illustrated in [13]. It combines non-linear programming, Newton-based methods and uses Hessian matrix load flow to minimize the quadratic objective. In [14], an OPF algorithm using Newton’s method with Hessian in place of Jacobian matrix and Lagrangian multipliers is provided. It provides good convergence when compared to its predecessors. An OPF problem which includes steady-state security is presented in [15]. It is an extension of [9] which includes exact constraints on outage contingencies. In [16], a solution to the optimal dispatch problem using Jacobian matrix is implemented. It provides rapid convergence which can be used in online control. An OPF algorithm based on reduced gradient method is proposed in [17]. It is used to minimize generator loading and optimize voltage levels. Load flow equations are represented as equality constraints. The authors in [18] have described an OPF method using reduced Hessian matrix with systematic constraint handling. It provides accurate solution, good convergence, and description about acceleration factor is provided. In [19], modified recursive quadratic programming (MRQP)-based OPF is implemented. MRQP is based on [11]. An algorithm to solve large OPF problem is presented in [20]. It decomposes the original large problem into set of subproblems which are constrained linearly using augmented Lagrangian.

In 1991, a landmark paper [21] classified various OPF techniques into two categories. Class A describes a series of algorithms which uses ordinary load flow to get an intermediate solution and this solution is under normal circumstances is close to the optimal load flow solutions. Using Jacobian matrix and various other sensitivity relations, optimization is performed iteratively. At each iteration, new load flow is performed. The optimal solution of this class strongly depends on the accuracy of load flow solution. With a load flow solution, set of voltages, and phase angles, Jacobian matrix and set of incremental power flow equations are available or can be extended. If a load flow solution exists, it already satisfies all the constraints. The optimization problem is solved separately by incorporating the sensitivity relation from before to arrive at an optimal one. In [9] an implementation of Class A algorithm is provided. Another example of such implementation is provided in parts one and two in [22] and [23] respectively using linear programming.

Class B refers to the class of algorithms which depend on exact optimal conditions and therefore use load flow equations as equality constraints. The optimal solution is dependent on detailed formulation of the OPF problem with the entire search space. This does not need a load flow solution. However, these kinds of problems are non-convex in nature. Therefore, convex relaxation or non-convex solvers are needed to compute the optimum which poses its own difficulties. It deals with the optimality conditions from Lagrangian function and comprises of derivatives of constraints and objective functions. Since the Hessian matrix is sparse and remains constant, it further increases the simplicity of this method and ease at which the optimum is achieved. Constraint handling is one of the biggest challenges of this class of algorithms. Using a heuristic method, constraints are handled as penalty terms which requires refactoring at every step and therefore, leads to degradation of the performance.

The above two classes of algorithms has various advantages and disadvantages. Performance of Class A directly depends in the performance of load flow techniques such as Newton-Raphson, Gauss-Seidel and widely used Fast Decoupled Method. It is shown in [24] that the above-mentioned methods have convergence and robustness issues. This may result in inaccurate load flow solutions. If the load flow does not result in a so-called high voltage or operable solution, Class A algorithms fails. In Class B algorithms, getting a global operable solution is challenging since it needs convex relaxation or heuristic techniques and the operable solution is difficult to achieve by respecting all the constraints.

The authors in this paper present a third class of algorithms, a Class C. This class combines Class A and Class B. It uses a reliable load flow described in Section 2 method wrapped around a heuristic to determine the optimal solution. The load flow provides accurate operable voltage and phase angle solution at every step and the heuristic uses this as equality constraints as described in Class B. Class C algorithms present various advantages. Operable voltage and phase angle solution is obtained at each iteration with the help of THELM. THELM always finds a solution, if it exists, irrespective of initial conditions whereas, Newton-Raphson load flow method leads to a non-convergent solution at very low or high loading conditions [24]. Since THELM is used in Class C method, the results are high voltage and operable. Global OPF solution can be obtained with a non-convex solver.

The following contributions and structure of the paper is as follows,

• Load Flow Solution to three-phase unbalanced distribution network using Holomorphic Embedding Load Flow Method (HELM) is described in Section 2.
• Benchmarking of HELM against established Newton-Raphson load flow solver from DIgSILENT PowerFactory [25] which is a well-known power system simulation and analysis software. This is discussed in Section 2.
• OPF using Distributed Genetic Algorithm, a Class C algorithm is described in detail in Section 3
• Simulation of OPF is performed to generate active and reactive power schedules at controllable nodes (see Section 5). This algorithm is applied to a real network in Austria.

2. Three-Phase Unbalanced Load Flow Method

A solution to the load flow problem is mostly obtained using numerical iterative methods such as Gauss-Seidel with its slow convergence and improved Newton-Raphson method, which provides better convergence [26,27]. Newton-Raphson method is computationally expensive since it must calculate Jacobean matrix at each iteration step in-spite of using sparse matrix techniques [28]. Various decoupled methods have been implemented which exploits the weak link between active power and voltage, in which Jacobian matrix needs to be calculated only once. One such method is Fast Decoupled Load Flow method which is widely used in the community [29]. The above-mentioned iterative techniques face similar problems with no guaranteed convergence since it depends on the initial conditions. This is due to the fact that load flow equations are non-convex in nature with multiple solutions. It is difficult to control the way these iterative methods converge to an operable solution [24]. In the literature, multiple implementations to improve the convergence of such traditional algorithm have been illustrated with limited success [30,31,32,33,34,35,36].

To use load flow methods in near or real-time applications, the physical models should fully deterministic and solved with reliability. HELM is one such candidate which can full fill these requirements [24].

Three-Phase Holomorphic Embedding Load Flow Method

Power flow equations, for example, the load bus equation described in Equation (1) is inherently non-analytical. Holomorphic principles can be applied to such equations by means of embedding a complex variable $\alpha$ such that the resulting problem is analytic in nature.

$${\displaystyle \sum _{k\in \Omega }}{Y}_{ik}{V}_k=\frac{S_i^{\ast }}{V_i^{\ast }},i\in {\Omega }_{PQ}$$

(1)

Voltage of the slack bus is assumed to be $V_0$ = 1.0 pu. and Bus 00 (see Appendix A) is always set to be slack bus.

Holomorphic embedding can be done in various methods. Equation (2) represents the simplest form. Bus voltages are the functions of the demand scalable complex variable $\alpha$.

$${\displaystyle \sum _{k\in \Omega }}{Y}_{ik}{V}_k\left(\alpha \right)=\frac{\alpha S_i^{\ast }}{V_i^{\ast }\left({\alpha }^{\ast }\right)},i\in {\Omega }_{PQ}$$

(2)

The research work in [24] suggests that the operable voltage solution can be obtained by analytic continuum of Equation (2) at $\alpha =1$ using the unique solution which exists when $\alpha =0$

$${\displaystyle \sum _{k\in \Omega }}{Y}_{ik}{V}_k\left(\alpha \right)=\frac{\alpha S_i^{\ast }}{{\overline{V}}_i\left({\alpha }^{\ast }\right)},i\in {\Omega }_{PQ}$$

(3)

$${\displaystyle \sum _{k\in \Omega }}{Y}_{ik}^{\ast }{\overline{V}}_k\left(\alpha \right)=\frac{\alpha S_i}{V_i\left(\alpha \right)},i\in {\Omega }_{PQ}$$

(4)

Equations (3) and (4), represent a set of polynomial equations and by using the Grobner bases, $V_i$ and ${\overline{V}}_i$ are holomorphic except for finite singularities.

$${\overline{V}}_i\left(\alpha \right)={\left(V_i\left({\alpha }^{\ast }\right)\right)}^{\ast },i\in \Omega$$

(5)

According to [24], if Equation (5) holds good, then Equations (3) and (4) can be reduced to Equation (2). Equation (5) is referred to as reflecting condition.

Since voltages of from Equation (2) for $\alpha =0$ as discussed above, it can be extended to power series described in Equation (6) and (7) at $\alpha =0$.

$$V_i\left(\alpha \right)={\displaystyle \sum _{n=0}^{\infty }}{V}_i\left[n\right]{\alpha }^n,i\in \Omega$$

(6)

$$\frac{1}{V_i\left(\alpha \right)}=W_i\left(\alpha \right)={\displaystyle \sum _{n=0}^{\infty }}{W}_i\left[n\right]{\alpha }^n,i\in \Omega$$

(7)

Equation (9) is obtained by substituting 7 into 2 and power series coefficients can be calculated to a desired degree.

$${\displaystyle \sum _{k\in \Omega }}{Y}_{ik}^{\ast }{\displaystyle \sum _{n=0}^{inf}}{V}_k\left[n\right]\left({\alpha }^n\right)=\alpha S_i^{\ast }{W}_i^{\ast }\left[n\right]{\alpha }^n$$

(8)

The following steps are involved to calculate voltages.

• For $\alpha =0$, solve Equation (9) to obtain a linear equation where the left-hand side of the equation represents the slack bus at which $V_0\left[\alpha \right]=1$.

  $${\displaystyle \sum _{k\in \Omega }}{Y}_{ik}{V}_k\left[0\right]=0,i\in {\Omega }_{PQ}$$

  (9)
• The reduced Y bus matrix is assumed to be non-singular. Equation (10) can be obtained based on the non-singularity assumption.

  $$W_i\left[0\right]=\frac{1}{V_i\left[0\right]}$$

  (10)
• Remaining power series coefficients can be obtained to the desired $n^{th}$ degree by equating the coefficients from Equation (11)

  $${\displaystyle \sum _{k\in \Omega }}{Y}_{ik}{V}_k\left[0\right]=S_i^{\ast }{W}_i^{\ast }[n-1],i\in {\Omega }_{PQ}n\ge 1$$

  (11)
  $W_i[n-1]$ are calculated using the lower order coefficients described in Equation (12).

  $$W_i[n-1]=-\frac{{\sum }_{m=0}^{n-2}{V}_i[n-m-1]W_i\left[m\right]}{V_i\left[0\right]}$$

  (12)
• Pade approximations which are particular kind of rational approximations are used for analytical continuum to determine the voltages at $\alpha =1$.

Based on the fundamentals of HELM discussed above, various research work dealing with enhancing or improving the method is available. One of the major deficiencies of the HELM described in [24] is that the PV/Generator bus is not defined. A PV bus model was presented in [37]. Ref. [38] presents an improved PV bus model and the major contribution of this paper is to provide alternative models capable of solving general networks. The authors have provided four methods with various parameters for PV bus. In the literature, three-phase formulation of HELM is lacking. In this paper, method four developed in [38] is extended to a novel three-phase unbalanced formulation which can be seen below. Equation (13) represents a general form of three-phase unbalanced HELM. Network models including various device models such as loads, generators, transformers are derived from the models developed in [39]. The seed solution, non-singularity of matrix A in Equation (14) and the reflective conditions of holomorphic functions are taken, as is, from [38]. Three-phase unbalanced form for a multi-bus system for PQ and PV bus types is presented below.

$$\left[\begin{array}{cccccc}{A}_1^a& A_1^b& A_1^c& A_2^a& A_2^b& A_2^c\\ A_{P{Q}_3}^a& A_{P{Q}_3}^b& A_{P{Q}_3}^c& A_{P{Q}_4}^a& A_{P{Q}_4}^b& A_{P{Q}_4}^c\\ A_{P{V}_3}^a& A_{P{V}_3}^b& A_{P{V}_3}^c& A_{P{V}_4}^a& A_{P{V}_4}^b& A_{P{V}_4}^c\end{array}\right]\left[\begin{array}{c}Re\left\{V^a\left[n\right]\right\}\\ Re\left\{V^b\left[n\right]\right\}\\ Re\left\{V^c\left[n\right]\right\}\\ Im\left\{V^a\left[n\right]\right\}\\ Im\left\{V^b\left[n\right]\right\}\\ Im\left\{V^c\left[n\right]\right\}\end{array}\right]=\left[\begin{array}{c}{r}_{1,n-1}\\ r_{P{Q}_{2,n-1}}\\ r_{P{V}_{2,n-1}}\end{array}\right]$$

(13)

Is of the form,

$$Ax=b$$

(14)

where the matrix A can be further clarified as,

$$\begin{array}{cc}\hfill A_{1_{ij}}^P=G_{ij}^P+\delta i,jRe\left\{y_i^P\right\},& i,j\in \Omega ,P\in a,b,c\hfill \\ \hfill A_{2_{ij}}^P=B_{ij}^P-\delta i,jIm\left\{y_i^P\right\},& i,j\in \Omega ,P\in a,b,c\hfill \\ \hfill A_{P{Q}_{3_{ij}}}^P=B_{ij}^P-\delta i,jIm\left\{y_i^P\right\},& i,j\in {\Omega }_{PQ},P\in a,b,c\hfill \\ \hfill A_{P{V}_{3_{ij}}}^P=2\delta i,j,& i,j\in {\Omega }_{PV},P\in a,b,c\hfill \\ \hfill A_{P{Q}_{4_{ij}}}^P=G_{ij}^P+\delta i,jIm\left\{y_i^P\right\},& i,j\in {\Omega }_{PQ},P\in a,b,c\hfill \\ \hfill A_{P{V}_{4_{ij}}}^P=0,& i,j\in {\Omega }_{PV},P\in a,b,c\hfill \end{array}$$

(15)

where G and B are the conductance and susceptance, respectively. $\delta i,j=1$ if $i=j$, else 0.

The right-hand side matrix elements are defined as follows,

$$r_{1,n-1,i}={\delta }_{n,1}(P_i-Re\left\{y_i\right\})-Re\left\{{\displaystyle \sum _{m=1}^{n-1}}{V}_i^{\ast }\left[m\right]\sum _{k\in \Omega }\sum _{p\in P}{Y}_{ik}{V}_k[n-m]\right\},i\in \Omega$$

(16a)

$$r_{PQ2,n-1,i}={\delta }_{n,1}(-Q_i-Im\left\{y_i\right\})-Im\left\{{\displaystyle \sum _{m=1}^{n-1}}{V}_i^{\ast }\left[m\right]\sum _{k\in \Omega }\sum _{p\in P}{Y}_{ik}{V}_k[n-m]\right\},i\in {\Omega }_{PQ}$$

(16b)

$$r_{PV2,n-1,i}=-{\displaystyle \sum _{m=1}^{n-1}}\sum _{p\in P}{V}_i^{\ast }\left[m\right]V_i[n-m]+(1+\alpha {(M_i-1)}^2)\left[n\right],i\in {\Omega }_{PV}$$

(16c)

where, $M_i$ is the target voltage magnitude for PV bus.

The power series were calculated for using the above equations and Viskovatov Pade approximant algorithm is used to determine the voltages and phase angles similar to the ones in [24,37].

3. Optimal Power Flow Model

As described in Section 1, OPF algorithms can be classified under three classes. There have been a lot of research on OPF and this can be seen in the vast array of work available in the literature.

Type C algorithms requires non-convex solvers to perform optimization problems. Non-convex solvers have been previously used to solve OPF problems but, they are used in the context of Class B algorithms. The authors in [40] have provided a method to plan reactive power flows optimally using generic algorithm as it provides optimum which is a global one. The proposed method is applied to two 51 and 224 bus networks. An OPF problem is solved using generic algorithm in [41] as a unified power flow controller to regulate branch voltages with respect to both angle and magnitude. It minimizes real power losses and security limits of power flows are maintained. Reactive power planning using hybrid genetic algorithm is presented in [42]. It uses genetic algorithm at the highest level and linear programming to get the optimum sequentially. This can be considered as a modified version of Class A. It uses genetic algorithm instead of just load flow to determine the initial converged solution to the OPF problem.

In [43], a feeder reconfiguration technique is presented. It uses genetic algorithm in the context of OPF to reduce losses in a distribution system. Switches are opened to determine the initial population. The authors in [44], have presented a hybrid evolutionary algorithm with multi-objective OPF. It is used to minimize losses, voltage, and power flow deviations and generator costs. In [45], optimal placement and sizing of capacitor banks in distributed networks using genetic algorithm is presented. The objective is to simultaneously improve the power quality and sizing of fixed capacitor banks.

In this paper, the OPF problem is formulated as follows,

$$\begin{array}{cc}\hfill \underset{x}{\mathrm{minimize}}& F(x,u)\hfill \\ \hfill \mathrm{subject}\mathrm{to}& \hspace{1em}\hspace{1em}u\in U\hfill \\ \hfill & \hfill G(x,u)=0,\\ \hfill & \hfill H(x,u)\le 0\end{array}$$

(17)

where,

x and u represent sets of state and input variables.

$F(x,u)$ is the objective function for the OPF problem. Typical objectives are total generator cost, loss minimization in network and in this paper, the objective function chosen is the three-phase unbalance minimization (see Section 4).

$G(x,u)$ and $H(x,u)$ represents the equality and inequality constraints of the OPF problem.

In the context of type C algorithms, accurate and reliable load flow is used as equality constraints. In this case, THELM is used.

Typical inequality constraints for a three-phase unbalanced distribution grids are enlisted below,

Limits on active power (kW) of a (generator) PV node:	$P_{Lo{w}_i}\le P_{P{V}_i}\le P_{Hig{h}_i}$
Limits on voltage (V (pu.)) of a PV or PQ node:	$|V_{Lo{w}_i}|\le |V_i|\le |V_{Hig{h}_i}|$
Limits on tap positions of a transformer:	$t_{Lo{w}_i}\le t_i\le t_{Hig{h}_i}$
Limits on phase shift angles of a transformer:	${\theta }_{Lo{w}_i}\le {\theta }_i\le {\theta }_{Hig{h}_i}$
Limits on shunt capacitances or reactances:	$s_{Lo{w}_i}\le s_i\le s_{Hig{h}_i}$
Limits on reactive power (kVAr) generation of a PV node:	$Q_{Lo{w}_i}\le Q_{P{V}_i}\le Q_{Hig{h}_i}$
Upper limits on active power flow in transmission lines or transformers:	$P_{i,j}\le P_{Hig{h}_{i,j}}$
Upper limits on MVA flows in lines or transformers:	$P_{i,j}^2+Q_{i,j}^2\le S_{Hig{h}_{i,j}}^2$
Upper limits on current magnitudes in lines or transformers:	$|I_{i,j}|\le |I_{Hig{h}_{i,j}}|$
Limits on voltage angles between nodes:	${\Theta }_{Lo{w}_i}\le {\Theta }_i-{\Theta }_j\le {\Theta }_{Hig{h}_i}$

In this paper, the non-convex solver used is a genetic algorithm, to minimize the objective function. Genetic algorithm is chosen due to its wide use in OPF techniques, ease of parallelizability to handle large networks and its probabilistic transition rule. The authors have used the method developed in [46]. Genetic algorithms of the kind, mixed integer non-linear non-convex is used to include all the constraints mentioned above. To accommodate THELM in genetic algorithm, it is included in the fitness function and penalty functions are used to include constraints.

4. Three-Phase Unbalance Minimization

As mentioned in Section 1, it is essential to manage the distribution network optimally and in a balanced fashion. Various methods have been presented in the literature to minimize three-phase unbalance. In [47], a method to minimize three-phase unbalance is presented. Reactive power compensation is performed using flexible AC transmission system (FACTS) devices to minimize the three-phase unbalance. It is applied to a simple study case of four bus system. This method does not provide optimal scheduling of loads and does not include all the buses in the network. It is applicable only to local grid where the FACTS devices are located. The authors in [48] have provided a method to minimizing network unbalance using phase swapping. A genetic and greedy algorithm is used to optimally swap the phases to generate a convenient solution, leading to a minimum number of swaps to minimize network unbalance. In [49], plug-in hybrid electric vehicles are used to minimize local three-phase unbalance. It does not include a grid perspective and is done only on the point of common coupling.

In this paper, OPF from Section 3 is applied to a real network in Austria. Figure 1 represents a real low-voltage distribution network. In this use case, the objective function is to minimize three-phase voltage unbalance which can be seen in Equation (18). This objective can be realized in multiple ways and in this paper, reference balanced voltages are used.

$$minimizeJ={\displaystyle \sum _{k\in \Omega }}{\displaystyle \sum _{p\in P}}{(real\left(V_{k,balanced}^p\right)-real\left(V_k^p\right))}^2+{(imag\left(V_{k,balanced}^p\right)-imag\left(V_k^p\right))}^2$$

(18)

where, $P\in phases(a,b,c)$ and B represented all the buses in the network. The voltages are represented in rectangular coordinate system with real part of voltage being the magnitude and imaginary part being the phase angle. Both real and imaginary value are considered because both phase and angle of voltages need to be balanced.

The controllable variables are per phase active (P) and reactive powers (Q) at buses 07, 15, 18, and 22. The single-phase loads are replaced with three-phase loads (see Figure 2).

For simplicity of representation, these three-phase loads are represented as single-phase loads with red coloring. This can be observed in Figure 1. P and Q on individual loads can be modulated by the OPF algorithm and can take values which are both positive and negative essentially, acting as a prosumer node.

5. Simulation Results

This section provides simulation results to the concepts presented in the previous sections. In Section 5.1, THELM described in Section 2 is validated against DIgSILENT PowerFactory Newton-Raphson algorithm. In Section 5.2, simulation results for three-phase unbalanced optimal power is presented with the three-phase unbalance minimization objective presented in Section 4.

5.1. Validation of THELM

THELM is benchmarked against load flow solver in an established power system analysis tool, DIgSILENT PowerFactory. Various simple networks are drawn with increased level of complexity (see Appendix A).

Voltages from 1000 random load flows by varying active and reactive power at load buses from ±10 kW and ±0.8 kVAr (which accounts for power factor 0.9) are generated and tabulated below.

Mann-Whitney U test is used to compare the sample means of voltage magnitudes from the two methods, since the samples are non-parametric in nature. It checks whether to accept or reject the null hypothesis. Mann-Whitney U test is similar to student’s T test but is suitable for non-parametric samples. A sample of 100 voltages from both THELM and Power Factory NR methods are used. Various statistical information and test results are tabulated in

Table 1. Benchmarking THELM against Power Factory NR method.

	Mean	Std	Min	25%	50%	75%	Max	Statistic	p Value
Test 00	0.00489	0.00350	1.14 × 10−5	0.00194	0.00430	0.00724	0.01591	2.38719	0.016977
Test 01	0.00572	0.00414	1.01 × 10−5	0.00237	0.00496	0.00846	0.01938	3.23993	0.001195
Test 02	0.00209	0.00027	0.00167	0.00186	0.00205	0.00228	0.00287	5.47075	4.48 × 10−8
Test 03	0.00015	0.00028	9.38 × 10−8	4.57 × 10−5	9.85 × 10−5	0.00018	0.00529	−4.69619	2.65 × 10−6
Test 04	0.00343	0.00246	1.08 × 10−5	0.00144	0.00296	0.00508	0.01125	−5.6259	1.84 × 10−8
Test 05	0.00178	0.00023	0.00142	0.00159	0.00175	0.00195	0.00241	−5.83359	5.42 × 10−9
Test 06	0.05586	0.01114	0.02464	0.04779	0.05540	0.06339	0.09113	5.69187	1.25 × 10−8

. Columns mean, standard deviation, min, 25%, 50%, 75%, max are calculated by taking the absolute difference between their respective voltages. All the data above is calculated by taking the average between various buses and phases. From the column statistic and p-value, it can be observed that their means are statistically insignificant, and the null hypothesis is accepted. The results from the test suggests that THELM produces results which are acceptable for load flow analysis with lower deviations from one another.

5.2. Simulation Results for Three-Phase Unbalanced Optimal Power Flow

Simulation is performed for the real grid detailed in Figure 1 using OPF algorithm described in Section 3. It is performed for one day from 2018-8-31 00:00:00 to 2018-9-01 00:00:00 with the sampling time of 15min (96 intervals). Load profiles are from smart-meter data, from real households and are acquired from all the buses in the network updating a database. Forecasted profiles are inputs to the OPF algorithm and optimal schedules are generated based on it. Load forecasting is performed for this time horizon using convolutional neural networks, using data until 2018-8-30 23:45:00. It is performed for one day (day ahead forecast) and more details able it is not provided since it is out of scope.

Load flow solution is non-causal in nature and to generate an optimal schedule for controllable buses, it must be run for all 96 intervals. OPF is performed using Class C algorithm presented in Section 3 for controllable buses described in Section 4. It can also be observed that the optimal schedules are generated for all the three phases and can take both positive and negative values. Real profile is recorded during day for uncontrollable loads at the buses.

Forecasted, optimal and real active and reactive power consumption profiles at one of the controllable buses (Bus 15) can be seen in Figure 3 and Figure 4 respectively. At Bus 15, all the loads are single phased (connected to phase C).

Active and reactive power for all the phases can be observed in Figure 3 and Figure 4.

Using the three schedules shown in Figure 3 and Figure 4, load flows are performed using THELM described in Section 2. Loads flows are performed for all intervals and are represented using box-plots.

Figure 5 describes the averaged objective function values based on Equation (18). It can be observed that the three-phase unbalance has been reduced from 0.879 for real and forecasted profiles to 0.529 for optimal profiles which accounts for 39% unbalance minimization based on the defined objective function (see Section 4).

From Figure 6, it can be observed that the voltages are indeed balanced, and the average values are close to balanced voltages. Additionally, the nature of the objective function used has also caused the voltages to cluster around 1 pu. since the balanced real part of the balanced voltage is exactly 1 pu.

6. Conclusions and Outlook

In this paper, a novel class of OPF algorithm is presented. It uses a novel three-phase unbalanced HELM presented in Section 2. Benchmarks are performed to test the performance of THELM and DIgSILENT Power Factory Newton-Raphson method described in Section 5.1. These benchmarks were performed on various test networks. Mann-Whitney U test was performed, and it was concluded that the results from both load flow methods were statistically indistinguishable and null hypothesis was accepted. Using THELM, optimal power flow method was developed using genetic algorithm in Section 3, describing the type C class of algorithms. The novel Class C algorithm provides various advantages over Class A and B OPF algorithms as discussed in Section 1. A use case with an objective function to minimize three-phase unbalance was applied to a real network in Austria in Section 5. The reason for choosing this objective is motivated by the requirements of the network operator and to handle the unbalance locally. It involves the generation of active and reactive power schedules for four controllable buses using smart-meter forecasts from other uncontrollable loads in the network (see Figure 1). Optimal schedules for these buses were generated and used to produce voltages using THELM and the results were described in Figure 6. It can be observed that the three-phase voltage unbalance has reduced up to 39% and the optimal average objective function values can be observed in Figure 5.

In future work, the scalability and replicability of the method needs to be analyzed. The method needs to be applied to various larger networks with large number of nodes. Simulation time and code optimization is not considered a priority for this study. To use this method in a real-time or near-real-time operation, the algorithm needs to be optimized. In this work, only three-phase unbalance minimization is used. OPF with various other objective functions need to be considered.