A Novel Linear Optimal Power Flow Model for Three-Phase Electrical Distribution Systems

This paper presents a new linear optimal power flow model for three-phase unbalanced electrical distribution systems considering binary variables. The proposed formulation is a mixed-integer linear programming problem, aiming at minimizing the operational costs of the network while guaranteeing operational constraints. Two new linearizations for branch current and nodal voltage magnitudes are introduced. The proposed branch current magnitude linearization provides a discretization of the Euclidean norm through a set of intersecting planes; while the bus voltage magnitude approximation uses a linear combination of the L1 and the L norm. Results were obtained for an unbalanced distribution system, in order to assess the accuracy of the linear formulation when compared to a nonlinear power flow with fixed power injections, showing errors of less than 4% for currents and 0.005% for voltages.

Abstract-This paper presents a new linear optimal power flow model for three-phase unbalanced electrical distribution systems considering binary variables. The proposed formulation is a mixed-integer linear programming problem, aiming at minimizing the operational costs of the network while guaranteeing operational constraints. Two new linearizations for branch current and nodal voltage magnitudes are introduced. The proposed branch current magnitude linearization provides a discretization of the Euclidean norm through a set of intersecting planes; while the bus voltage magnitude approximation uses a linear combination of the L1 and the L norm. Results were obtained for an unbalanced distribution system, in order to assess the accuracy of the linear formulation when compared to a nonlinear power flow with fixed power injections, showing errors of less than 4% for currents and 0.005% for voltages.
Index Terms-Linear optimal power flow, unbalanced, threephase, electrical distribution systems, mixed-integer linear programming The classical optimal power flow (OPF) is a nonlinear, non-convex, optimization problem, usually containing both continuous and discrete variables [1]. The solution of the OPF aims at defining the participation of controllable energy resources to the power balance, aiming at minimizing selected objectives while satisfying technical and physical operational limits. The inclusion of three-phase models increases the size and the complexity of the problem. Moreover, the addition of integer variables regarding the operating status of distributed generator (DG) units, energy storage system (ESS), capacitor banks, transformer's tap position, among others; increases the challenge of obtaining optimal solutions in polynomial time, since solving non-convex mixed-integer nonlinear programming (MINLP) models is NP-hard [2].
On the other hand, mixed-integer linear programming (MILP) models are convex, tractable, and easier to solve; which is why obtaining accurate approximated models has been important to the research community for more than 50 years [3].
A compendium of optimal power flow formulations applied to modern electrical distribution systems (EDS) is proposed in [4] including convex and non-convex approaches. Mathematically equivalent convex OPF formulations are available either for meshed transmission networks [5], or for radial distribution networks [6]. Convexification also includes linearization techniques, focused on obtaining approximated models which are good representations of the original nonlinear problem as in [7], where authors introduce a linear OPF considering voltage stability constraints or in [8], where a linearly constrained quadratic OPF is proposed. However, all the aforementioned works only consider single-phase or balanced configurations.
Three-phase OPF formulations for EDS have also been studied, as in [9] where a genetic algorithm is used for solving a non-convex MINLP. Reference [10] introduces an OPF for unbalanced distribution systems considering stochastic variables as a nonlinear programming model. Authors in [11] propose a MILP model considering droop-controlled generation units, however the formulation relies on the knowledge of previous operating conditions. A semidefinite formulation is proposed in [12], the formulation is convex, however, it does not consider integer variables.
This paper proposes a new MILP model for the OPF of unbalanced three-phase EDS. The proposed formulation accounts for dispatchable DG units only, although it can be easily extended to include ESS, renewable energy sources, under-load tap changers, etc. The objective function aims at minimizing the operational cost of the network, while respecting nodal voltage and branch current magnitude limits. The obtained linear model is based on two new linearization procedures for the Euclidean norm, present in voltage and current magnitudes. The linearizations are suitable for different operating conditions and have been tested under the particular characteristics of radial distribution networks. The new formulation is tested using an unbalanced distribution system of 25-buses using AMPL [13] and solved with CPLEX [14] under default settings. The results obtained with the new MILP model are compared with the results from a conventional nonlinear power flow, assessing the accuracy of the introduced approximations.

II. MINLP OPTIMAL POWER FLOW
The OPF problem can be modeled as a MINLP problem using (1)- (15), where italic characters represent state variables, e.g., i r jk,f , i i DG g,f , P DG g , etc. The objective function in (1) aims at minimizing the total operational cost of the EDS regarding the energy imported from the main grid and the energy injected by DG units. Constraints (2)- (14) are defined for each phase, i.e., f ∈ Ω f .
subject to: Real and imaginary current balances are defined in (2) and Similarly, voltages at each node are calculated using (4) and (5), fixing the voltage from the electrical substation . Three-phase active power through the substation is calculated in (6), which is linear since voltage at the substation is known. Similarly, (9) and (10) represent the three-phase active and reactive power from each DG unit. Expressions relating active and reactive power from loads are shown in (11) and (12), respectively. Operational limits regarding the capability curve of DG units are considered in (7) and (8), while nodal voltage and branch magnitude limits are considered in (14) and (13), respectively. Finally, (15) states the binary nature of variable µ g .

III. LINEARIZATION OF NONLINEAR CONSTRAINTS AND PROPOSED MILP MODEL
Note that nonlinearities in the MINLP model are related to the product of variables in power injections, as seen in (9) (12), and to the calculation of magnitudes in (13) and (14).

A. Linearization of power injections
Active and reactive power injections are of the form Linear expressions can be obtained using approximations based on estimated values for the real and imaginary parts of the nodal voltages {v r0 , v i0 }. The first-order approximations are: where y represents either loads (y = D) or generators (y = DG). The estimated values {v r0 , v i0 } can be selected using an initial three-phase load flow analysis, experiencebased values, or even a flat-start, as it was performed in [15].

B. Linearization of branch current magnitudes
Branch current magnitudes in (13) are handled using a new linear approximation to the Euclidean norm. Assume α = max{|i r kj,f |, |i i kj,f |} and β = min{|i r kj,f |, |i i kj,f |} for all f ∈ Ω f , then: Note that there are four combinations regarding the signs of α and β considering the absolute value, and two combinations relating the maximum and minimum values. These combinations result in the addition of 8 |Ω l | |Ω f | linear equations. This is done to consider the solution space where the currents are held, as seen in the contour of the constraints in Fig. 1a, where the 8 plains, which are product of the constraints, can be seen. Simulations showed acceptable errors for branch currents within a wide range, i.e., i r = i i = [−1000, 1000], as shown in Fig. 1b. It can be seen that the maximum asymptotic error is 4%, which implies that it does not increase further even for wider ranges.
v Fig. 2. Representation of voltage magnitudes, limits, and proposed range angle.

C. Linearization of nodal voltage magnitudes
Linearizing voltage magnitudes requires more accurate approximations, since even small errors can lead to severe voltage violations. Authors in [15], for example, used geometrical relationships to limit voltage magnitudes between an angle range. In this paper, it is proposed the use of a similar approach using a limited angle range, based on a linear approximation to the Euclidean norm (L 2 ). The approximation is a function of the L 1 norm and the L ∞ norm, as explained in [16]. Since L 1 ≥ L 2 ≥ L ∞ , it is possible to express L 2 as a linear combination of the other two, as: where L ∞ (x) = max{|x r |, |x i |} and L 1 (x) = |x r | + |x i |.
Hence, voltage magnitudes are approximated as: An angle range, namely ±Θ, is introduced to improve the quality of the approximation, as shown in Fig. 2. This is done considering that voltage angles remain between a small deviation from the nominal angles of each phase (θ f ) in EDS. For |Θ| < 15 • it can be shown that for each phase, Values for a f and b f are obtained using regression techniques. These coefficients depend on the value set for Θ as shown in Table I, where 10,000 randomly generated points were used for the regressions. The maximum percent error between the linear approximations and the exact norm is shown for each phase.
The maximum percent errors of the approximation as a function of the voltage angles using Θ = 2.5 • are shown for phase A and phase B in Fig. 3a and Fig. 3b, respectively. Phase C behaves as a mirror of phase B with negative angles. The contours of the voltage magnitude constraints are shown in Fig. 3c and Fig. 3d using Θ = 10 • . It can be seen that error  The maximum percent error obtained for different values of Θ is shown in Fig. 4. The figure was obtained after comparing the exact norm and the proposed approximation for 1,000 points within each range angle and selecting the biggest error among them. The errors show a quadratic behavior, as can be seen from the fitting equations defined as E (Θ) f ; with a steeper growing rate for phases B and C than for phase A. The selection of optimal values for Θ, i.e., producing the minimum error, depends on the expected maximum angle deviation presented in the system.

D. Proposed MILP model
The proposed OPF for unbalanced, three-phase EDS formulated as an MILP model is summarized in (21).

IV. TEST SYSTEM AND RESULTS
The proposed linearized model was tested on an unbalanced 25-bus test system, whose load data, topology, and lines' parameters can be found in [17]. Three DG units have been added to the original system at buses 13, 19, and 25. Parameters for DG units along with other system information can be found in    Fig. 4, it is expected that the maximum percent errors will be lower than 0.003% for phase A and lower than 0.010% for phases B and C.
Voltage magnitudes found with the nonlinear power flow are shown in Fig. 5a. It can be seen that the voltage magnitudes at buses 15 and 17 are the ones closer to the lower voltage limit. In fact, voltage at bus 15 for phase B is 4.75 × 10 −5 pu lower than the predefined limit. The percent error between the obtained voltage magnitudes using the proposed MILP model and the nonlinear power flow is shown in Fig. 5b. It can be seen that errors for all phases are lower than 0.006%. Specifically for bus 15, the error at phase B is approximately 0.005% which matches with the small voltage violation seen in Fig. 5a.
Current magnitudes from the nonlinear power flow are shown in Fig. 6a. It can be seen that all branch current magnitudes were held within the maximum current limit. The percent error between the obtained current magnitudes with the linear approximated model and the nonlinear power flow is shown in Fig. 6b. It can be seen that errors for all phases are lower than 4%. Most branch currents in Phase B show the biggest errors compared the errors at phases A and C; this is due to the particular angle of the current flowing through the branches at each phase, as can be inferred by Fig. 1.
Although the observed voltage violation is small (4.75 × 10 −5 pu), it could be managed by obtaining more conservative solutions, e.g., using V < 0.95 for the MILP model. Another way of solving this problem is reducing the error of the approximation for phases B and C by performing a linear rotation, for example. The last solution could be implemented in future work. The selection of the optimum value for Θ to reduce the error of the approximation demands a previous expectation of the maximum angle deviation in the system. This can be obtained by running a previous power flow or by historical data.
The approximation for branch current magnitudes presents a maximum percent error of 4%. Although not desirable, small overloads can be accepted during short periods of time in overhead distribution lines. On the other hand, reducing the approximation error is also possible by increasing the amount of intercepting planes. This would, however, increase the size of the problem.
A second test is performed, this time only fixing the binary variables µ g found with the MILP and allowing the nonlinear formulation to find the respective dispatches. Objective function costs are displayed in Table III for different range angles along with the percent error between the linear approximation and the nonlinear power flow. It can be seen that the error increases with the range angle, as expected. It can also be seen that the objective function for Θ = 10 • presents a different configuration, involving the commitment and dispatch of another DG unit. This last statement can be seen also from Fig. 7, where the power dispatch is shown for every DG unit depending on the range angle used. It can be inferred that the approximation errors are reflected in the units committed and the dispatches of these units.

V. CONCLUSIONS
A new linear OPF model for three-phase unbalanced EDS considering binary variables was proposed in this paper. The proposed formulation has been conceived as an MILP prob-  lem, aiming at minimizing the operational costs of the network while respecting operational constraints. Two new linearizations for branch current and nodal voltage magnitudes were introduced. The proposed branch current magnitude linearization provides a discretization of the Euclidean norm through a set of intersecting planes; while the bus voltage magnitude approximation uses a linear combination of the L 1 and the L ∞ norm. The maximum expected errors using both approaches were discussed, as well as the effect of the introduced angle range, Θ. Results were satisfactory in terms of the accuracy of the linear approximations when compared to a nonlinear power flow with fixed power injections. The effect of different values for Θ regarding the value of the objective function and the error with the nonlinear approach was also discussed. The proposed linearizations can be adapted to several other planning/operation problems in EDS, such as reconfiguration, energy management, Volt/Var control, expansion planning, among others.
Future work can be done towards including ESS, capacitor banks, voltage regulators, renewable energy sources, and the stochasticity of exogenous parameters.