Real-Value Formulations of Unbalanced Optimal Power Flow
Geth and Ergun

mode = title]Real-Value Power-Voltage Formulations of, and Bounds for, Three-Wire Unbalanced Optimal Power Flow

Corresponding author<EMAIL_ADDRESS>
1]Frederik Geth

[1]Energy Systems, CSIRO Energy, 10 Murray Dwyer Circuit, Newcastle 2304 NSW, Australia

2,3]Hakan Ergun

[2]KU Leuven, Kasteelpark Arenberg 10, Heverlee, Belgium

[3]EnergyVille, Thor Park 8310, Belgium

Unbalanced optimal power flow refers to a class of optimization problems subject to the steady state physics of three-phase power grids with nonnegligible phase unbalance.
Significant progress on this problem has been made on the mathematical modelling side of unbalanced OPF, however there is a lack of information on implementation aspects as well as data sets for benchmarking. One of the key problems is the lack of definitions of current and voltage bounds across different classes of representations of the power flow equations.
Therefore, this tutorial-style paper summarizes the structural features of the unbalanced (optimal) power problem for three-phase systems.
The resulting nonlinear complex-value matrix formulations are presented for both the bus injection and branch flow formulation frameworks, which typically cannot be implemented as-is in optimization toolboxes.
Therefore, this paper  also  derives the equivalent real-value  formulations, and discusses challenges related to the implementation in optimization modeling toolboxes.
The derived formulations can be re-used easily for  continuous and discrete optimization problems in distribution networks for a variety of operational and planning problems.
Finally, bounds are derived for all variables involved, to further the development of benchmarks for unbalanced optimal power flow, where consensus on bound semantics is a pressing need. We believe benchmarks remain a cornerstone for the development and validation of scalable and reproducible optimization models and tools.
The soundness of the derivations is confirmed through numerical experiments, validated w.r.t. OpenDSS for IEEE test feeders with 3x3 impedance matrices.

mathematical optimization unbalanced power flow

§ NOMENCLATURE

This article depends on the definition of a  variety of scalar, vector and matrix parameters and variables related to grid buses and branches (see Tables 1-5). The core variables are current, voltage and power, whereas parameters are mainly impedance and variable bounds.
Fig. <ref> summarizes the variables and parameters defined in the fundamental $3\times3$ branch model for which the well known The $\Pi$-equivalent model used. With respect to balance networks, both series and shunt elements are represented by full complex-valued matrices including the mutual impedance coupling between the conductors. Using the $\Pi$-equivalent model, the branch current $\IijkSDP$ can be split into a series component $\IijskSDP$ and a shunt component $\IijshkSDP$, respectively.  All circuit element voltages are defined w.r.t. (local) ground voltage $\VikG=0\,V$.

Unbalanced $3\times3$ $\Pi$-model branch in voltage and current variables.

Table <ref> illustrates typography and mathematical notation used throughout;
Table <ref> defines sets and indices;
Table <ref> defines parameters;
Table <ref> defines typical engineering variables;
Table <ref> defines lifted variables. Black and red colors indicate real-valued variables and parameters, respectively. Blue and brown colors are used for complex-valued variables and parameters, respectively.

Typography and mathematical notation

[]L s

x    real scalar variable

𝐗    real vector or matrix variable

x    complex scalar variable

𝐗    complex vector or matrix variable

x    real scalar parameter

𝐗    real vector or matrix parameter

x    complex scalar parameter

𝐗    complex vector or matrix parameter

𝒳    set

𝐗^     transpose of   $\mathbf{X}$

𝐗^*    conjugate  of   $\textcolor{\complexcolor}{\mathbf{X}}$

𝐗^= (𝐗^*)^    conjugate transpose of  $\textcolor{\complexcolor}{\mathbf{X}}$

∘    element-wise multiplication

⊘    element-wise division

imaginary unit, satisfies $ \imagnumber^2 = -1$

ab     polar notation of complex number $ a \cdot e^{\imagnumber b}$

nm     set of real $n\times m$ matrices

nm       set of complex $n\times m$ matrices

n ⊂nn       set of Hermitian $n\times n$  matrices

(𝐗)     extract diagonal of $\textcolor{\complexcolor}{\mathbf{X}}$, $\diag: \setComplexMatrix{n}{n} \rightarrow \setComplexMatrix{n}{1} $

Sets and indices

Phases     $\indexPhases,\indexPhasesTwo \in \setPhases = \{\symPhaseA, \symPhaseB, \symPhaseC  \}  $

Branches    $\indexGridLines\in  \omegaGridLines $

Buses     $\indexGridNode \in  \omegaGridNodes $

Topology (forward)    $\indexGridLines \indexGridNode\indexGridNodeTwo \in  \setTopology \subseteq  \omegaGridLines \times \omegaGridNodes \times \omegaGridNodes $

Topology (reverse)      $\indexGridLines  \indexGridNode \indexGridNodeTwo \in  \setTopologyReverse = \{\indexGridLines \indexGridNodeTwo \indexGridNode  \, | \,  \indexGridLines \indexGridNode\indexGridNodeTwo \in  \setTopology \} $

Topology     $ \indexGridLines \indexGridNode\indexGridNodeTwo \in \setTopologyBoth = \setTopology \cup \setTopologyReverse  $

Bus pairs     $ \indexGridNode\indexGridNodeTwo \in  \setBuspairs  = \{ ij \, | \,  \indexGridLines \indexGridNode\indexGridNodeTwo \in  \setTopology   \} \subseteq \omegaGridNodes \times \omegaGridNodes $

Units     $\indexUnit \in  \setUnits $

Unit connectivity     $\indexUnit \indexGridNode \in  \setUnitTopology \subseteq  \setUnits \times \omegaGridNodes $

Shunt connectivity     $\indexShunt \indexGridNode \in  \setShuntTopology \subseteq  \setUnits \times \omegaGridNodes $

Bus voltage magnitude min./max. ($V$)     $\VikSDPmin, \VikSDPmax     \in \setRealMatrix{|\setPhases|}{1}$

Bus phase angle diff. min./max. (rad)    $\thetaiSDPmin, \thetaiSDPmax  \in \setRealMatrix{|\setPhases|}{1}$

Branch current rating ($A$)    $ \IijkSDPrated     \in \setRealMatrix{|\setPhases|}{1}$

Branch apparent power rating ($V\!A$)    $ \SijkSDPrated     \in \setRealMatrix{|\setPhases|}{1}$

Branch series  impedance  ($\Omega$)       $  \zijksSDP   \in \setComplexMatrix{|\setPhases|}{|\setPhases|} $

Branch series  admittance   ($S$)      $  \yijksSDP     \in \setComplexMatrix{|\setPhases|}{|\setPhases|} $

Branch from/to shunt  admittance  ($S$)       $  \yijkshSDP,\yjikshSDP   \in \setComplexMatrix{|\setPhases|}{|\setPhases|}$

Bus pair angle diff. min./max. (rad)    $\thetaijSDPmin, \thetaijSDPmax  \in \setRealMatrix{|\setPhases|}{1}$

Bus shunt  admittance ($S$)       $  \ybSDP   \in \setComplexMatrix{|\setPhases|}{|\setPhases|}$

Unit current rating ($A$)    $ \IuunitSDPrated     \in \setRealMatrix{|\setPhases|}{1}$

Unit active power bounds ($W$)    $ \PuunitSDPmin,\PuunitSDPmax      \in \setRealMatrix{|\setPhases|}{1}$

Unit reactive power bounds ($var$)    $ \QuunitSDPmin,\QuunitSDPmax       \in \setRealMatrix{|\setPhases|}{1}$

Optimization variables

Bus voltage  ($V$)       $  \VikSDP, \VjkSDP    \in \setComplexMatrix{|\setPhases|}{1}$

Branch current  ($A$)         $ \IijkSDP, \IjikSDP     \in \setComplexMatrix{|\setPhases|}{1} $

Branch series current    ($A$)    $  \IijskSDP,\IjiskSDP       \in \setComplexMatrix{|\setPhases|}{1} $

Branch shunt  current   ($A$)     $ \IijshkSDP, \IjishkSDP   \in \setComplexMatrix{|\setPhases|}{1} $

Branch power flow   ($W$)       $   \SijkSDP, \SjikSDP   \in \setComplexMatrix{|\setPhases|}{|\setPhases|}  $

Branch series power flow     ($W$)     $    \SijksSDP, \SjiksSDP   \in \setComplexMatrix{|\setPhases|}{|\setPhases|} $

Unit current   ($A$)   $  \IuunitSDP \in \setComplexMatrix{|\setPhases|}{1}$

Unit  power     ($W$)      $   \SuunitSDP   \in \setComplexMatrix{|\setPhases|}{|\setPhases|}$

Lifted optimization variables

\begin{equation*}
\begin{IEEEeqnarraybox}[
\IEEEeqnarraystrutmode
]{s L}
%Entity & \text{Symbol} \\
\hline
Bus voltage product  ($V^2$)  &    \VisqkSDP, \VjsqkSDP \in   \setHermitianMatrix{|\setPhases|} \\
Bus pair voltage product  ($V^2$)  &    \VijsqkSDP, \VjisqkSDP \in  \setComplexMatrix{|\setPhases|}{|\setPhases|} \\
Branch   current product    ($A^2$)   & \IijsqkSDP, \IjisqkSDP \in   \setHermitianMatrix{|\setPhases|} \\
Branch series  current product    ($A^2$)   & \IijsqskSDP \in   \setHermitianMatrix{|\setPhases|} \\
System voltage product  ($V^2$)  & \VsqkSDP\in  \setHermitianMatrix{|\setBuspairs||\setPhases|} \\
Bus pair  voltage product  ($V^2$) &  \VbpsqkSDP\in  \setHermitianMatrix{2|\setPhases|} \\
Branch voltage-current product   (mix) &  \VlinesqkSDP\in  \setHermitianMatrix{2|\setPhases|}\\
% Unit current  product ($A^2$) & \IusqkSDP  \in \setHermitianMatrix{|\setPhases|}\\
% Unit  voltage-current product   ($-$) & \VunitsqkSDPcurrent  \in \setHermitianMatrix{2|\setPhases|} \\
% Bus ground voltage-current product ($-$) &\MisqkSDP  \in \setHermitianMatrix{2} \\
\hline
\end{IEEEeqnarraybox}
\end{equation*}

§ INTRODUCTION

Driven by the increased rollout of distributed energy resources (DERs) such as PV, battery storage as well as electric vehicles, power distribution grids are facing a number of challenges associated with the large scale integration of these technologies.
In low-voltage grids specifically, one can observe phase unbalance readily, due to the presence of a significant amount of single-phase loads and DERs. Many of the current electric vehicles use single phase charging and small rooftop solar systems are connected via single phase inverters to the low voltage grid, further increasing the  unbalance. Furthermore, phase  unbalance can stem from insufficient conductor transposition in radial distribution networks.
Unbalance implies an underutilization of the grid's transfer capacity, as it leads to higher losses and to faster-than-expected congestion. Simulation techniques, i.e. deriving solutions to the power flow  (PF) equations in unbalanced networks has long been a topic of interest [1, 2, 3], and is used in deriving hosting capacity by means of scenario analysis.  The accuracy of modeling of low-voltage grids has  been studied in-depth by Urquhart [4].

Unbalanced optimal power flow (OPF) refers to the mathematical optimization of problems subject to the physics of unbalanced grids, and serves as the core for a variety of problem classes such as benchmarking centralized or distributed optimal control solutions [5], determination of  expansion options and  hosting capacity analysis considering control actions.

§.§ State of the Art on Unbalanced (O)PF

Similar to the balanced (positive sequence) modeling, one can distinguish between the bus injection model (BIM) and the branch flow model (BFM) formulations of the unbalanced PF equations. BIM forms eliminate all current variables, which leads to active and reactive power flows being expressed purely as a function of the voltage differences between connected buses. The series impedance is consequently used in admittance form, which makes it impossible to represent zero-impedance branches. Conversely, the BFM forms keep (a representation of) the current variable through the series impedance. In this case, series impedance is represented in impedance form, therefore the edge case of zero series impedance remains representable.

Most OPF problems are developed in the complex power-voltage variable space, instead of the current-voltage variable space common in PF solvers.
Table <ref> maps a number of published formulations to BIM/BFM categories and the variable spaces in which they are defined.
We refer to [6, 7]  for recent in-depth reviews of  mathematical formulations for the OPF problem.

A nonlinear programming (NLP) unbalanced OPF formulation for  branches is presented in  [8], in which the shunt impedances have been neglected. The unbalanced current-voltage form, with generation power dispatch constraints, is derived in [9]. The nonconvex power-lifted voltage form is presented in [10].
The first rank-constrained semi-definite programming (SDP) of unbalanced OPF, and its SDP relaxation, are proposed by Dall'Anese et al.  in [11].
Gan et al. developed the BFM variant, and  defined both BFM and BIM in a consistent notation in  [12], without the provision of the real-value form.
Extensions to these SDP relaxations have been proposed:
Zhao et al. develop models for delta-connected loads [13];
Bazrafshan et al. propose extensions for voltage regulators [14];
Usman et al. discuss neutral conductor modeling [15];
Claeys et al. present detailed transformer models for unbalanced OPF [16];
Claeys et al. also present convex relaxation of delta/wye ZIP unbalanced loads [17];
Vanin et al. explore further relaxation to SOC problems [18].

It is noted that the formulations so far published are not easily compared due to the inherent complexity of the notation, the variety of notations used, and lacking details of approximations applied during implementation in modeling software.
For instance, the branch shunt impedance is often neglected, assumed to be diagonal, or simply modeled as bus shunts.
This leads to inaccuracies in the power and current flow values and in the enforcement of the proper branch flow limits, e.g. when only the series current is bounded, not the total including the shunt current.
This in turn makes it hard to validate the feasibility and correctness of results on published data sets which include such components without modification. Examples include the IEEE PES distribution test feeders [19]; the IEEE123 bus system specifically contains non-diagonal branch shunt matrices.

§.§ Bound definitions and data

In an OPF problem, authors typically assume [20, 6, 7] the power flow is subject to a subset of:

* apparent power, upper bound;

* current magnitude, upper bound;

* voltage magnitude, lower and upper bound;

* voltage angle difference between adjacent buses, lower and upper bound.

These bounds can be generalized to the case with phase unbalance. Furthermore, in the context of phase unbalance, additional limits have been discussed [21]. The most obvious missing bound is that of phase angle differences at a certain bus.
Note that bounds represent a key aspect of optimization problems, however currently there are few unbalanced power flow data sets that include all of these parameters. This makes benchmarking optimality of different unbalanced OPF engines challenging.

§.§ Implementation of Complex-Value  Optimization Problems

Figure <ref> illustrates the typical sequence of formulating, implementing and solving optimization problems as commonly applied to engineering problems.
First, a mathematical model of the optimization problem is formulated.
The problem  is formulated so that it can be handled by the chosen modeling language and optimization solver.
Eventually, the modeling toolbox  translates the mathematical model into the appropriate form for the solver interface,   fills out parameter values, and dispatches the  optimization solver.
Finally, the mathematical solver uses one or a variety of different algorithms to solve the optimization problem, and returns a solution.

Formulation, implementation and solution of optimization problems.

In the literature, real-value formulations of the convex relaxations of unbalanced OPF are rarely discussed. CVX has been used as an optimization  modeling tool in [11]. It is one of the very few optimization packages with support for complex optimization variable primitives, as well as matrix variables. CVX automates the reformulation process from a complex-value matrix to a set of real-value scalar variables and constraints. Such functionality is also available in the JuMP extension [22], which provides JuMP [23] with the functionality to automate the real-value reformulation process.

Gilbert and Josz [24] illustrate with a limited-scope solver that there are potential advantages of solving SDP problems directly in the complex domain. Nevertheless, a mature implementation of such a solver has not been developed so far. Therefore, the interest in real-value forms remains high.

Unbalanced power flow formulations, key references and variable spaces

Formulation     Ref.     Variable space      Exact?     Section
NLP    Kirchhoff unbalanced 	    [1, 9]     $\IijkSDP, \VikSDP $		     by definition, discussion: [25]     <ref>

NLP    BIM unbalanced  scalar polar	     [8]    $\diag(\SijkSDP),\VikSDP $		     if $\VikSDPmin \neq \mathbf{0}$     <ref>

NLP    BIM unbalanced  scalar rect.	         $\diag(\SijkSDP),\VikSDP $		     if $\VikSDPmin \neq \mathbf{0}$      <ref>

SDP    BIM unbalanced   	   [11, 26, 27]     $\SijkSDP, \VisqkSDP  \succeq 0, \VijsqkSDP ,\VsqkSDP \succeq 0 $		      $\rank(\VsqkSDP)=1$    <ref>

SDP    BIM unbalanced  radial 	    [12]     $\SijkSDP, \VisqkSDP \succeq 0, \VijsqkSDP, \VbpsqkSDP \succeq 0 $		        $\rank(\VbpsqkSDP) = 1$, only  if radial    <ref>

SDP    BFM unbalanced     	    [12]      $ \SijkSDP, \VisqkSDP \succeq 0, \IijsqskSDP \succeq 0, \VlinesqkSDP \succeq 0 $		    $\rank(\VlinesqkSDP)=1$, only  if radial     <ref>

§.§ Scope, Contributions and Paper Structure

In summary, we observe  the following gaps in the literature:

* polar and rectangular real forms in the power-voltage variable space for  branches with shunts ($\Pi$-model);

* real form of the power - lifted voltage BIM and BFM SDP relaxations;

* definitions of current, power and voltage angle difference limits in both the power-voltage and power-lifted voltage variable spaces.

This work therefore provides consistent derivations of real-value formulations for variants of   unbalanced optimal power flow  in different variable spaces, to enable straight-forward  implementation in  optimization modeling toolboxes and  efficient comparison on published data sets.

First, the basic relationship between the different variables and parameters is defined in  <ref>, which builds the mathematical foundation of this paper; it also contains a generalized model to represent loads and generators in the system.
Section <ref> provides  Kirchhoff's and Ohm's laws in the multiconductor form.
Next, <ref> develops a formulation of nonlinear BFM unbalanced power flow and
<ref> the equivalent BIM.
Furthermore, <ref> derives the BIM and <ref> the BFM in the lifted variable space which is the basis various convex relaxations.

Moreover <ref> presents the feasible sets of the derived real-value formulations and discusses implementation aspects.
Finally, <ref> presents the conclusions.

§ NOTATION AND BASIC RELATIONSHIPS

§.§ Scalar and Matrix Variables

The  voltage $\VikSDP$ of a bus $i$ is a complex value, vector variable, encapsulating variables for each conductor:

+ .

Similarly,  the total current of a branch $l$ connecting a pair of buses $i$ and $j$ is a complex value vector variable,

The  sending-side current ($i \rightarrow j$) through the series element in the $\Pi$-model is defined as $\IijskSDP$, the sending-side shunt  current as $ \IijshkSDP$, respectively.

The complex power flow in branch $l$ from bus $\indexGridNode$ to $\indexGridNodeTwo$  depends on the bus voltage of the sending-side $\VikSDP$ and the conjugate transpose (indicated with superscript $\hermitiantranspose$) of the current  $\IijkSDP$,

= [

We observe that $ \diag(\SijkSDP  )=\VikSDP \circ (\IijkSDP)^*$,
where $\circ$ is the element-wise (Hadamard) product.
It is noted that $\rank(\SijkSDP) = 1$, as it is defined as the outer product of two vectors.

The off-diagonals relate to the diagonal elements according to,
/ = / = / = ()^*,

/ = / = / = ()^*,

/ = / = / = ()^* ,
which means that the off-diagonals of $\SijkSDP$ are scaled and rotated versions of the more easily interpretable diagonal elements.

§.§ Variable Bounds

Voltage magnitudes have minimum and maximum operational limits, which are specific to each bus and phase,

0 ≤

_  ≤

_ .

In this paper, we overload `$\leq$' for `$\geq$' for vectors and matrices to indicate element-wise inequality (conversely, `$\succeq$' is used as the symbol for matrix positive semidefiniteness).
Recognizing we obtain the magnitude squared by multiplying complex numbers with their own conjugates, the nodal voltage bounds can also be presented quadratically as:
0 ≤∘≤∘()^* ≤∘.

Note that the upper bound constraints are convex (space inside a circle), but the lower bound constraints are nonconvex (space outside of a circle).

Next,  apparent power limits
are defined using the absolute value of the diagonals of the branch flow matrix,
0 ≤

= ,

which is equivalent to second order cone (SOC) constraints,

()^2 + ()^2

()^2 + ()^2

()^2 + ()^2

(  )^2

(   )^2

(   )^2

and can succinctly be written as,
0≤() ∘()^* ≤∘.

The magnitudes of the branch current should stay below rated values and are bounded for the diagonals of the branch current matrix,

=  ,

which again can be written as a set of SOC constraints,
0 ≤∘()^* ≤∘.

A valid SOC representation of current magnitude limits  (<ref>) that does not require explicit current variables $\IijkSDP$ is useful for BIM forms.
Using the nodal voltage magnitudes, the branch current limits can directly be enforced on the power flows,

(  )^2

(   )^2

(   )^2

which can also be developed in matrix notation  as,
() ∘()^* ≤∘∘∘()^*.

Valid bounds on all elements of the power flow matrix $\SijkSDP = \PijkSDP + \imagnumber \QijkSDP$ are,
- ()^ ≤, ≤()^ ,

both for active and reactive power components, respectively.

The voltage angle differences between connected buses $i$ and $j$ are bounded,

-  ≤

For a relatively balanced voltage phasor, we expect $ \VikANangle -\VikBNangle \approx  \VikBNangle    -\VikCNangle  \approx  \VikCNangle   -\VikANangle  \approx \frac{2\pi}{3}$.
The voltage angle differences between phases on buses $i$ can be bounded to enforce angle balance relative to the expected 120 degrees,

The voltage phasor in  reference buses $i \in \setReferenceNodes \subset \omegaGridNodes$ is assumed fixed, e.g.,
=  =

§.§ Branch Impedance

The  circuit series impedance matrix  is defined as a full matrix with no assumption on the particular structure:

= +

where each element consists of a series resistive and reactive impedance. It is noted that in physical systems, we expect $\rijksSDP \succeq 0, \xijksSDP \succeq 0$.
The impedance matrix can be rewritten in the corresponding admittance form,
= ()^-1 =  + .

where $ {(\zijksSDP)^{-1}}$ is the matrix inverse of $\zijksSDP$. In case of missing conductors, e.g. single, or two-conductor connections,  $\zijksSDP$ is not invertible but it is valid to use the Moore-Penrose inverse instead.
The shunt admittances at the sending and receiving sides, respectively $\yijkshSDP, \yjikshSDP$, are defined,

= + =

Although  for typical distribution lines and cables the shunt admittances are diagonal and have equal values for both sides, in this generalized derivation they can be considered as a two different and full matrices. This also allows for re-use of the representation for other elements, such as  transformers.

§.§ Loads and Generators as Units

In this paper, we define units $\indexUnit \in \setUnits$ to generalize loads, generators and storage elements. The current flowing from the connected bus into the unit is

The unit current is bounded by the corresponding current rating $\IuunitSDPrated$ similarly to (<ref>).
Tuples of units and the buses they are connected to are defined in the connectivity set $\indexUnit \indexGridNode \in  \setUnitTopology $.
The active and reactive power consumed by a unit is defined $\SuunitSDP = \PuunitSDP + \imagnumber \QuunitSDP = \VikSDP (\IuunitSDP)^\hermitiantranspose$. We  define bounds on active/reactive power dispatch separately,
≤() ≤,

≤( ) ≤.

We do not provide any further detail on the specific modeling of delta or wye-connected units and refer to [13, 17].

§.§ Shunts

A shunt element, e.g. shunt capacitance, inductance or resistance, with index $\indexShunt$ has an admittance  $\ybSDP = \gbSDP + \imagnumber \bbSDP$.
Tuples of shunts and the bus they are connected to are defined in the connectivity set $\indexShunt\indexGridNode \in \setShuntTopology$ .
The current  from the bus to the shunt is,

and is bounded by the rated current  $\IhshuntSDPrated$.

§ APPLICATION OF CIRCUIT LAWS

This section illustrates how Kirchhoff's and Ohm's laws are used in the context of multi-conductor branches with matrix impedances for the representation of the circuit physics, providing the relationship between the voltage, current and power variables.

§.§ Branch Model

Ohm's law for branch $lij$ representing the voltage drop along the branch is formulated in matrix form using nodal voltages, the branch current and the impedance matrix,
=- .

Kirchhoff's current law (KCL) is used to split up the series and shunt (to ground) currents in the $\Pi$-section,
=   + ,
=   + ,

+ = 0.

Fig. <ref> defines the harmonized single-wire equivalent  vector/matrix variables in both natural and lifted (defined in upcoming sections)  variable spaces for a clear presentation of the basic relationships. The scalar representation of all variables and parameters is defined in Fig. <ref>.

Power flow and losses in unbalanced  $\Pi$-model branch

§.§ Shunts and Units at Nodes

The unit model can be adapted for a variety of optimization problems, e.g. optimal dispatch of distributed generation, optimal scheduling of storage and electric vehicle charging and optimal demand management. Such extensions essentially define feasible sets dependent on $\diag(\PuunitSDP), \diag(\QuunitSDP)$. Without loss of generality, in this paper, we focus on the feasible sets defined by (<ref>)-(<ref>).

Similarly, the shunt power $\SbSDP$ is defined as
= + =   ()^.

The shunt current relates to the voltage through Ohm's law,
=  .

Substituting (<ref>) into (<ref>), the shunt power consumption can be expressed  depending on the nodal voltage only,
=  ()^  ()^ .

§.§ Bus Model: Kirchhoff's Current Law

Kirchhoff's Current Law (KCL) is conventionally expressed in current variables, but can also be lifted to the complex power variable space.
KCL in current variables at each bus $i$ is,
∑_lij ∈ +  ∑_∈ +  ∑_ bi ∈ = 0.

KCL in the complex power variable space is obtained by taking the conjugate transpose of  (<ref>) and element-wise multiplying with $\VikSDP\neq 0$ on the left,
∑_lij ∈ () +  ∑_∈ ( ) +  ∑_ bi ∈ ( )=0,

which means that the diagonal elements of the apparent power matrices of the connected branch flows, units and shunt need to sum to zero for each bus. Commonly, the shunt power expression (<ref>) is substituted into this equation. The real-value equivalent forms are obtained,
∑_lij ∈( ) +  ∑_∈( )   +  ∑_ bi ∈ ( ) =0,

∑_lij ∈ () +  ∑_∈ ()  +  ∑_ bi ∈ ( ) =0.

§ UNBALANCED BRANCH FLOW MODEL

This section illustrates how the unbalanced BFM is derived.
First, we define a variable for the complex power flow $\SijksSDP$ through the series element of the $\Pi$-section,
= + = ( )^ .

§.§ Power Flow Model

Passive components in electrical circuits cause losses. In the BFM we distinguish between the losses associated with the from-side shunt, i.e. $\SijkshlossSDP$,  the series impedance, i.e. $\SijkslossSDP$, and the to-side shunt, i.e. $\SjikshlossSDP$.
The series losses  depend on the voltage drop over the impedance and the current flow through it,
= (- ) ()^.

The voltage drop itself is can be derived from Ohm's law,
(- )  =  .

Substituting (<ref>) into (<ref>) we obtain the series loss
=  ()^ = ()^ ,

which is symmetric because of (<ref>).
The shunt losses $\SijkshlossSDP$ are derived from the current through the shunt, i.e. $\IijshkSDP$, and the voltage at the shunt, i.e. $\VikSDP$,
=   ()^,

The shunt current relates to the voltage through Ohm's law,
=  .

Substituting (<ref>) into (<ref>), the shunt losses are,
=  ()^  ()^ ,
=   ( )^  ()^ .

As the the sum of the sending and receiving side power flows need to equal the branch losses, the branch loss balance (see Fig. <ref>) can be written using the different loss components (<ref>) and (<ref>),
+ = + +

= ( )^ ()^  + ()^ +  ()^ ()^ .

§ UNBALANCED BUS INJECTION MODEL

This section illustrates how the unbalanced BIM is obtained for which the branch flow current is substituted by nodal voltage varibales.

§.§ Power Flow Model

Therefore, we write (<ref>) in admittance form,

= (- ).

The sending end apparent power flow can be calculated using the sending end voltage and the sum of the series current and the sending end shunt current, respectively,

\begin{equation}
\SijkSDP = \VikSDP(   \IijskSDP + \IijshkSDP)^{\hermitiantranspose}. \label{eq_flow_from}
\end{equation}

By substituting (<ref>) and (<ref>) into (<ref>) we derive,

\begin{equation}
\SijkSDP = \VikSDP(        \VikSDP)^{\hermitiantranspose} ( \yijkshSDP       )^{\hermitiantranspose} +  \VikSDP  (\VikSDP - \VjkSDP)^{\hermitiantranspose} (\yijksSDP)^{\hermitiantranspose} \label{eq_bim_complex_power},
\end{equation}

which is the nonlinear complex matrix form of the BIM. Note that this constraint is defined $\forall \indexGridLines \indexGridNode\indexGridNodeTwo \in \setTopologyBoth$, i.e. for  both the sending end $\SijkSDP$ and the receiving end $\SjikSDP$, respectively.
Choosing a coordinate system for $\VikSDP$, e.g. polar or rectangular representations of the conductor voltages, one can derive the polar, respectively rectangular, unbalanced power flow formulations by expanding the diagonal entries in (<ref>). In the coming subsections we are deriving the real-valued formulation of the BIM in both the polar and rectangular coordinate systems.

§.§ Unbalanced BIM Polar Scalar Form

The active and reactive power flow in each conductor can be written in a paramterized way, using nodal voltages per conductor $\textcolor{\complexcolor}{U_{i,\indexPhases}}  = \VikpNmag \imagangle \VikpNangle$.
The expressions of the diagonal elements of $\SijkSDP$ can be parameterized, $\indexPhases,\indexPhasesTwo \in \setPhases = \{\symPhaseA, \symPhaseB, \symPhaseC  \} $ and the real-value formulation of active and reactive power is obtained using trigonometric functions,

\begin{multline}{}
\Pijkpp  =
\sum_{\indexPhasesTwo \in \setPhases} \VikpNmag
\VikppNmag \cos(\VikpNangle \!- \VikhNangle)  \left( \gijksph + \gijkshph \right)  \\
+  \sum_{\indexPhasesTwo \in \setPhases} \VikpNmag  \VikppNmag \sin(\VikpNangle \!- \VikhNangle) \left( \bijksph + \bijkshph \right)   \\
-     \sum_{\indexPhasesTwo \in\setPhases}  \VikpNmag \VjkppNmag\!   \cos(\VikpNangle \!- \VjkhNangle)\gijksph    \\
-    \sum_{\indexPhasesTwo \in \setPhases}  \VikpNmag \VjkppNmag\!   \sin(\VikpNangle \!- \VjkhNangle) \bijksph   , \label{eq_active_polar_pf}
\end{multline}

\begin{multline}{}
\Qijkpp =
-\sum_{\indexPhasesTwo \in \setPhases} \VikpNmag
\VikppNmag \cos(\VikpNangle \!- \VikhNangle)  \left( \bijksph + \bijkshph \right)  \\
+  \sum_{\indexPhasesTwo \in\setPhases} \VikpNmag  \VikppNmag \sin(\VikpNangle \!- \VikhNangle) \left( \gijksph + \gijkshph \right)  \\
+   \sum_{\indexPhasesTwo \in\setPhases}  \VikpNmag \VjkppNmag\!   \cos(\VikpNangle \!- \VjkhNangle)\bijksph    \\
-   \sum_{\indexPhasesTwo \in \setPhases}  \VikpNmag \VjkppNmag\!    \sin(\VikpNangle \!- \VjkhNangle) \gijksph   . \label{eq_reactive_polar_pf}
\end{multline}

Thus, the active and reactive power flow through each conductor is obtained using the mutual coupling of the nodal voltages and branch impedances.
The bus shunt expressions are,
∑_∈ cos(- )

+  ∑_∈ sin(- )  ,

-∑_∈ cos(- )

+  ∑_∈ sin(- ) .

The bus pair voltage angle difference constraint (<ref>) and phase angle difference constraint (<ref>) directly apply in this variable space.

§.§ Unbalanced BIM Rectangular Scalar Form

+   .

tan() ∘

+     ≤

-  ≤tan() ∘

+ .

tan( 2π/3 +  ) ∘

+     ≤

-  ≤tan( 2π/3 +  ) ∘

+ .

Using   $\textcolor{\complexcolor}{U_{i,\indexPhases}} =  \VikpNreal + \imagnumber\VikpNimag $, the  expressions for active power of the diagonal elements of $\SijkSDP$ can be parameterized, $\indexPhases,\indexPhasesTwo \in \setPhases$, as,

\begin{multline}{}
\Pijkpp =
\sum_{\indexPhasesTwo \in \setPhases} ( \VikpNreal   \VikppNreal + \VikpNimag   \VikppNimag)  \left( \gijksph + \gijkshph \right)  \\
+  \sum_{\indexPhasesTwo \in \setPhases} ( \VikpNimag   \VikppNreal - \VikpNreal   \VikppNimag)  \left( \bijksph + \bijkshph \right)   \\
-     \sum_{\indexPhasesTwo \in \setPhases}  ( \VikpNreal   \VjkppNreal + \VikpNimag   \VjkppNimag)\gijksph    \\
-     \sum_{\indexPhasesTwo \in \setPhases}  ( \VikpNimag   \VjkppNreal - \VikpNreal   \VjkppNimag) \bijksph   ,  \label{eq_active_rectangular_pf}
\end{multline}

\begin{multline}{}
\Qijkpp =
- \sum_{\indexPhasesTwo \in \setPhases} ( \VikpNreal   \VikppNreal + \VikpNimag   \VikppNimag)  \left( \bijksph + \bijkshph \right)  \\
+  \sum_{\indexPhasesTwo \in \setPhases} ( \VikpNimag   \VikppNreal - \VikpNreal   \VikppNimag)  \left( \gijksph + \gijkshph \right)   \\
+    \sum_{\indexPhasesTwo \in \setPhases}  ( \VikpNreal   \VjkppNreal + \VikpNimag   \VjkppNimag)\bijksph    \\
-   \sum_{\indexPhasesTwo \in \setPhases} ( \VikpNimag   \VjkppNreal - \VikpNreal   \VjkppNimag) \gijksph   .\label{eq_reactive_rectangular_pf}
\end{multline}

The bus shunt expressions are,
∑_∈ ( + )

+  ∑_∈ ( - )   ,

- ∑_∈ ( + )

+  ∑_∈ ( - )   .

In rectangular coordinates, the nodal voltage bounds as defined in (<ref>) become

( )^2

()^2 + ()^2

()^2 + ()^2

()^2 + ()^2

( )^2

(  )^2

(   )^2

which are nonconvex for strictly positive voltage lower bounds.

The voltage angle difference constraint between buses (<ref>) needs to be reformulated. We use the tangent identity (<ref>).

Because the angle difference is bounded by $[-\pi/2, \pi/2]$ for voltage stability, we can rewrite this as quadratic nonconvex constraint (<ref>),  where the tangent function is applied element-wise to the components of the vector.
The phase angle difference constraint (<ref>) is equivalent to,
2π/3 + ≤

-   ≤2π/3 +   .

For the tangent function to be increasing and invertible, we want to further restrict ourselves to,
π/2 ≤2π/3 + ≤2π/3 +   ≤π,

which implies,
- π/6 ≤≤≤π/3.

With this restriction we can derive the quadratic nonconvex constraint (<ref>) from (<ref>), using an identity similar to (<ref>).

§ LIFTING OF BUS INJECTION MODEL

This section defines new variables to represent products of voltages and currents for the BIM, to enable a lift-and-project approach. This approach is commonly used for SOC and SDP relaxations of the nonlinear power-voltage formulation.

§.§ Lifted Variables

§.§.§ Bus Voltage

We define an auxiliary variable for the  voltage products, $\VisqkSDP$, satisfying
=  + =     ()^ ,

≽0, () = 1.

Note that (<ref>) is a quadratic nonconvex constraint, and (<ref>) is the well-known rank-constrained SDP equivalent form.
We illustrate the structure of $\VisqkSDP =\VikSDP(\VikSDP)^{\hermitiantranspose} $ as a real-valued matrix in rectangular coordinates,

- 	     			0

-     -      0

Note that this representation requires 9 unique scalar variables (underlined) and that the diagonal is real-valued,
( )  =

One can define bounds on the matrix entries as,
∘≤( ) ≤∘,

- ()^ ≤, ≤()^ .

In this variable space, the phase voltage angle difference constraint (<ref>) becomes,
tan( 2π/3 + )  ∘

≤tan( 2π/3 + )  ∘

The reference bus phasor is fixed,
=  ()^.

§.§.§ Bus Voltage Cross Product

We define a variable $\VijsqkSDP$ for the cross-product of the voltages $\VikSDP$ and $\VjkSDP$ of the buses associated with a branch (i.e. all bus-pairs) as used in equations (<ref>-<ref>) and (<ref>-<ref>),

=   ()^  .

Although $\VijsqkSDP$ is rank-1 by construction, it is not Hermitian.
Note that this definition implies,
= ()^  .

One can define bounds on the matrix entries of $\VijsqkSDP$ as,
- ()^ ≤, ≤()^ .

§.§ Power Flow Model

Now power flow equation (<ref>) can be written using the lifted variables,
= ( )^ +  (-  ) ()^ .

The real-value equivalents are,
= ( + )^ + ( +)^

-  ()^   -  ()^ ,

= ( +)^ - ( +)^

-  ()^ +  ()^ .

Similarly, the bus shunt power with lifted variables is,
=  ()^,

and its real-value equivalents are,
= ( )^ + ( )^,

= ( )^ - ( )^   .

Using (<ref>), we lift the branch current limit (<ref>) to the $\VisqkSDP$ variable space as
() ∘()^* ≤∘∘(),

which is a convex SOC constraint.

§.§ Voltage Angle Difference Bound

We note the `tangent inequality' [28] can be extended to the three-phase case in the following way:

) = () ⊘  ()

The voltage angle difference bounds therefore are,
tan(  ) ∘()

≤tan(  ) ∘().

§.§ Rank-Constrained SDP Model

§.§.§ Meshed Grids

Note that (<ref>) can be generalized (to support meshed grids), by shaping the matrices into a block matrix $\VsqkSDP$,

⋮    ⋮    ⋱     ⋮

≽0, ( ) =1.

Note that $\VsqkSDP  \in \setHermitianMatrix{|\setPhases||\setBuspairs|}$ with $|\setBuspairs|$  the number of  unique bus pairs in the topology.
If there is no branch between buses $i$ and $j$, the corresponding $\VijsqkSDP = \mathbf{0} $.
The real-value equivalent SDP constraint [29], for $\VsqkSDP = \VsqkSDPreal + \imagnumber \VsqkSDPimag$, is,

-       ≽0,

-       ) =1.

Chordal relaxation   can be employed to replace $\VsqkSDP$  with a set of smaller matrices [26, 27].

§.§.§ Radial Grids

In case of radial grids, a decomposition of $\VsqkSDP$ is easily derived.
Only constraints of type (<ref>) are retained [27]:
∀ij ∈:  =

≽0, (  ) =1.

Note that $\VbpsqkSDP  \in \setHermitianMatrix{2|\setPhases|}$.
Note that $\VbpsqkSDP \succeq 0$ implies both  $\VisqkSDP \succeq 0$ and $\VjsqkSDP \succeq 0$.
The real-value equivalent form  [29] is,

()^           (-)^

-      -

()^    -       ()^        ≽0 .

() = 1

§ LIFTING  OF BRANCH FLOW MODEL

This section defines new variables to represent products of voltages and currents specific to the BFM and then details the lift-and-project approach taken.

§.§ Lifted Variables

§.§.§ Series Current

The auxiliary variable for the  current products, $\IijsqskSDP$ satisfies,

= ()^ = ()^  = +

≽0, () = 1

Which has the following  representation in scalar variables:

_ +

- 	      			0

- 	    -       0

This representation requires 9 unique variables (underlined).

§.§.§ Total Current

The total current $\IijkSDP$ is lifted as,
= + =        ()^

≽0 ,  () = 1

We start from the definition
$  \IijkSDP =   \IijskSDP + (\yijkshSDP)  \VikSDP $ and  multiply both sides with their conjugate transpose,
()^ = (+ ()  )(()^ +  ()^()^ ).

Substituting in the lifted variables we obtain,
=+ ()  ()^+ ()  + ()^ ()^ ,

proving the lifted total current variable is a linear combination of $\IijsqskSDP$, $\VisqkSDP$ and $\SijksSDP$. In the real domain this becomes,

= +  ()^  + ()^

+  ()^ - ()^

+  - +()^  - ()^  ,

= +  ()^  + ()^

()^ - ()^

+ + -  ()^  -  ()^ .

One can define current bounds (<ref>) on the matrix entries of $\IijsqkSDP$ as,
0≤() ≤∘,

- ()^ ≤, ≤()^ .

Note that (<ref>)-(<ref>) is equivalent to (<ref>), however the former is linear due to (<ref>) -(<ref>) and the latter is quadratic-convex (SOC).

There are no direct bounds on the auxiliary variable for the series current $\IijsqskSDP$, however we can derive valid bounds through the known bounds total current and voltage,

-( ()^ + ||  ()^ ||^ )

≤()^ + ||  ()^ ||^  ,

where $|\yijkshSDP|$ indicates the element-wise application of the absolute value operation to obtain the magnitude.

§.§ Power Flow Model

Given $\SlkslossSDP  =  \zijksSDP  \IijsqskSDP$,
$\SijkshlossSDP = \VisqkSDP (\yijkshSDP)^{\hermitiantranspose}$, and
$\SjikshlossSDP =  \VjsqkSDP (\yjikshSDP)^{\hermitiantranspose}$, the power balance (<ref>) of a branch becomes,

= ()^  + + ()^.

The equivalent active and reactive power expressions are,

+ =  ()^  +  ()^+ -

+()^  +  ()^ ,

+ =  ()^ -    ()^ +  +

+ ()^ -    ()^ .

Either $\SijkSDP$ or $\SijksSDP$ can be substituted out through,
= +    ()^,
= -    ()^ .

§.§ Ohm's Law

Ohm's law (<ref>) is reformulated  by multiplying both sides with their Hermitian adjoint,

()^ = (- ) (- )^  ,

=   ( )^ -  ()^ -  ()^
+ ()^  ()^         .
Substituting in the lifted variables into (<ref>), we obtain,

= -  ()^ -  ()^    + ()^      .

The equivalent real expressions are,
=  - ()^ - ()^
-()^ -()^

+ ()^-()^
+  ()^ + ()^,

=  -()^ + ()^
- ()^  + ()^

+ ()^ + ()^
- ()^ + ()^     .

These equations are symmetrical, so it is sufficient to generate the scalar constraints only on the upper triangle for $\VjsqkSDPimag$; upper triangle and diagonal for $\VjsqkSDPreal$ (generating all scalar constraints is redundant).

§.§ Voltage Angle Difference Bound

We can derive $\VijsqkSDP$ as a function of $\VisqkSDP$ and $\SijksSDP$,

= -   ()^ ,

which in the reals becomes,
= -   ()^ -   ()^ ,

= -   ()^ +   ()^ .

Therefore, to avoid introducing the variable $\VijsqkSDP$ in the BFM, we substitute this into (<ref>),
tan(  ) ∘(-   ()^ -   ()^)

≤(-   ()^ +   ()^)

≤tan(  ) ∘(-   ()^ -   ()^).

§.§ Rank-Constrained SDP Model

The product (<ref>)  formulated in the lifted variable space is,
∀lij ∈:  =

()^      ,

≽0, (  ) =1.

Note that $\VlinesqkSDP  \in \setHermitianMatrix{2|\setPhases|}$.
Furthermore $\VlinesqkSDP\succeq 0$ necessitates $\VisqkSDP \succeq 0$ and $\IusqkSDP \succeq0 $ but not $\VjsqkSDP \succeq 0$. So for leaf buses, an explicit $\VjsqkSDP \succeq 0$ is needed.
The real-value equivalent form is,

()^           -()^

-     -

()^     -      ()^      ≽0  ,

() = 1 .

§ IMPLEMENTATION OF REAL-VALUE FEASIBLE SETS

Table <ref> section summarizes the feasible sets in the real domain.
The  shunt    $\PbSDP, \QbSDP$  feasible sets can be substituted into KCL and is therefore not listed explicitly.

Formulation feasible sets: variable spaces, bounds, constraints and substitution options (NCNL = nonconvex nonlinear, NCQ = nonconvex quadratic, QCP = quadratically constrained programming)

2c$\SijkSDP-\VikSDP$ 	 	      	    2c$\SijkSDP-\VisqkSDP$

3-4  6-7
AC polar     AC rectangular        BIM lifted relax.     BFM lifted relax.

Overall complexity         NLP     QCP         SDP      SDP
Voltage variables     $\VikSDP$    $[ \VikpNmag \imagangle \VikpNangle]_{ \indexPhases \in \setPhases}$   (<ref>)	    $[ \VikpNreal + \imagnumber\VikpNimag]_{ \indexPhases \in \setPhases}$    (<ref>)        	      $\VisqkSDPreal, \VisqkSDPimag$ (<ref>)       $\VisqkSDPreal, \VisqkSDPimag$ (<ref>)  ,

Unit power variables    $\SuunitSDP$       $\PuunitSDP, \QuunitSDP $ (<ref>)(<ref>)	      $\PuunitSDP, \QuunitSDP $ (<ref>)(<ref>)     	        $\PuunitSDP, \QuunitSDP $ (<ref>)(<ref>)         $\PuunitSDP, \QuunitSDP $ (<ref>)(<ref>)

Power flow variables     $\SijkSDP$     $\PijkSDP, \QijkSDP$  (<ref>)	      $\PijkSDP, \QijkSDP$(<ref>)     	     $\PijkSDP, \QijkSDP$   (<ref>)      $\PijkSDP, \QijkSDP$  (<ref>)

Current variables     $\IijkSDP$     eliminated 	      eliminated      	      eliminated        $\IijsqkSDPreal, \IijsqkSDPimag  $

Additional variables     	    -      -         $\VijsqkSDPreal, \VijsqkSDPimag$  (<ref>)      $\IijsqskSDPreal, \IijsqskSDPimag$  (<ref>) (<ref>),
Voltage bounds     $\VikSDPmin, \VikSDPmax$ 	     LP (<ref>)     NCQ(<ref>)      	      LP(<ref>), (<ref>)     LP (<ref>), (<ref>),

Phase angle difference bounds      $\thetaiSDPmin, \thetaiSDPmax$	    LP (<ref>)      NCQ (<ref>)      	     LP (<ref>)      LP (<ref>)

Power flow bounds     $\SijkSDPrated$    SOC (<ref>)      SOC (<ref>)     	      SOC (<ref>)       SOC (<ref>)

Current bounds      $\IijkSDPrated$	    SOC (<ref>)      SOC(<ref>)      	     SOC (<ref>)      LP (<ref>), (<ref>)

Voltage angle difference bounds      $\thetaijSDPmin, \thetaijSDPmax$	     LP (<ref>)      NCQ (<ref>)      	      LP (<ref>)      LP (<ref>)

Kirchhoff's current law     	    LP (<ref>), (<ref>)      LP (<ref>), (<ref>)      	     LP (<ref>), (<ref>)       LP (<ref>), (<ref>)

Bus shunt equations (admittance)     	    NCNL (<ref>), (<ref>)     NCQ (<ref>), (<ref>)      	     LP (<ref>), (<ref>)     LP (<ref>), (<ref>)

Power flow equations  ($\Pi$-model)    	    NCNL (<ref>), (<ref>)     NCQ (<ref>), (<ref>)      	     LP (<ref>), (<ref>)     LP (<ref>), (<ref>)

Kirchhoff's voltage law     	    implicit      implicit      	     SDP (<ref>) or (<ref>), (<ref>)      implied by Ohm's law when radial

Ohm's law     	    implied       implicit      	      implicit       LP (<ref>), (<ref>)

Complex power definition          implicit       implicit      	      implicit       SDP (<ref>)

Reference bus fixed phasor      	    LP (<ref>)       LP (<ref>)        	     LP (<ref>)       LP (<ref>)

Substitution     	    -     -      	     -       $\IijsqkSDPreal, \IijsqkSDPimag \rightarrow \IijsqskSDPreal, \IijsqskSDPimag$ (<ref>), (<ref>)

$\PijkSDP, \QijkSDP \leftrightarrow \PijksSDP,   \QijksSDP$ (<ref>), (<ref>)

$ \VijsqkSDPreal, \VijsqkSDPimag \rightarrow  \VisqkSDPreal, \VisqkSDPimag, \PijksSDP, \QijksSDP$ (<ref>), (<ref>)

The AC polar and rectangular forms can be implemented in optimization toolboxes with support for nonlinear programming. The core power flow equations (<ref>), (<ref>) or (<ref>), (<ref>) are nonlinear, however only the rectangular form is quadratically representable.

Solving the lifted BIM and BFM forms requires dropping the rank constraints (<ref>) or (<ref>) to obtain a SDP formulation.
It is noted that further SOC relaxation can be performed [30, 18].
Furthermore,  for the tightness of the SOC relaxation of the SDP constraints, it is best to work in the complex domain as long as possible.
The process for obtaining the original (nonlifted) current and voltage variables is detailed in [31].

Implementations of the real-value formulations are available in the package PowerModelsDistribution.jl [32].
This package extends PowerModels.jl [23] and is built on top of JuMP  [23], a Julia package for mathematical programming.
In the implementation, where possible, the matrix forms are used, and the scalarization of the equations is performed by JuMP.
The discussed substitution and elimination of $\IijsqkSDPreal, \IijsqkSDPimag$, $\PijksSDP, \QijksSDP$ and $\VijsqkSDPreal, \VijsqkSDPimag$ in the BFM  are performed through JuMP expressions, which improves the readability of the mathematical model in the code. We note that nondimensionalization (per unit scaling) is performed in the implementation, in an effort to improve numerical conditioning.

Using PowerModelsDistribution.jl, one can read OpenDSS power flow case files and run power flow and optimal power flow case studies.
OpenDSS is used to validate the accuracy of the power flow formulations.
The worst relative voltage error for IEEE 13, 34, 123 and LV test case  [19] is  1.4E-7 [32].

§ CONCLUSIONS

In this paper, the electrical physics of three-phase grids are derived in complex-power matrix variables and different voltage variable spaces, i.e. polar, rectangular and lifted. Using the different variable spaces, the unbalanced power flow formulations are derived for  generic $\Pi$-model branches with asymmetric shunt impedance (full matrices, no structure assumed).
The real-value matrix and scalar equivalent formulations are derived, and implemented in an open-source software package.

The derived framework can easily be extended and applied to a variety of continuous and discrete optimization problems in distribution networks.
Expressions are presented to enforce voltage magnitude, angle, power and current bounds across the formulations in an exact manner, without needing to resort to auxiliary variables. Nevertheless, including auxiliary variables instead of eliminating them does not by default increase calculation time. Comparing these variants numerically requires the development of an extensive unbalanced OPF test case library, which includes values for the different bounds discussed, not just voltage magnitudes.
Future work includes adding more bounds, e.g. for sequence voltage magnitudes and unbalance metrics [21], supporting explicit representations of neutral conductors, and developing benchmarks for comparing optimality and scalability across formulations and solvers.

§ ACKNOWLEDGEMENT

Part of this work has been performed in the framework of the ADriaN project supported by Fluvius. Special thanks to  Reinhilde D'hulst, Carleton Coffrin, David Fobes, Sander Claeys, Thomas Brinsmead.

[1]
R. Berg, E. Hawkins, and W. Pleines, “Mechanized calculation of unbalanced
load flow on radial distribution circuits,” IEEE Trans. Power Appar.
Syst., vol. PAS-86, no. 4, pp. 415–421, 1967.

[2]
T. H. Chen, M. S. Chen, K. J. Hwang, P. Kotas, and E. A. Chebli,
“Distribution system power flow analysis - a rigid approach,” IEEE
Trans. Power Del., vol. 6, no. 3, pp. 1146–1152, 1991.

[3]
C. Cheng and D. Shirmohammadi, “A three-phase power flow method for real-time
distribution system analysis,” IEEE Trans. Power Syst., vol. 10,
no. 2, pp. 671–679, 1995.

[4]
A. J. Urquhart, “Accuracy of low voltage distribution network modelling,”
PhD Dissertation, Loughborough University, 2016.

[5]
N. Dkhili, J. Eynard, S. Thil, and S. Grieu, “A survey of modelling and smart
management tools for power grids with prolific distributed generation,”
Sustainable Energy, Grids and Networks, vol. 21, p. 100284, 2020.
[Online]. Available: <https://doi.org/10.1016/j.segan.2019.100284>

[6]
D. K. Molzahn and I. A. Hiskens, A survey of relaxations and
approximations of the power flow equations, 2019.

[7]
D. Bienstock, M. Escobar, C. Gentile, and L. Liberti, “Mathematical
Programming formulations for the alternating current optimal power flow
problem,” pp. 1–40, 2020. [Online]. Available:

[8]
B. Mahdad, T. Bouktir, and K. Srairi, “A three-phase power flow modelization:
a tool for optimal location and control of FACTS devices in unbalanced power
systems,” in IEEE IECON, 2006, pp. 2238–2243.

[9]
F. Geth, S. Claeys, and G. Deconinck, “Current-voltage formulation of the
unbalanced optimal power flow problem,” in 8th Workshop Modeling
Simulation Cyber-Physical Energy Syst., Sydney, Australia, 2020, pp. 1–6.

[10]
——, “Nonconvex lifted unbalanced branch flow model: derivation,
implementation and experiments,” Electric Power Systems Research,
vol. 189, no. December, p. 106558, 2020.

[11]
E. Dall'Anese, G. B. Giannakis, and B. F. Wollenberg, “Optimization of
unbalanced power distribution networks via semidefinite relaxation,” in
North American Power Symp., Champaign, IL, USA, 2012, pp. 1–6.

[12]
L. Gan and S. H. Low, “Convex relaxations and linear approximation for
optimal power flow in multiphase radial networks,” in Power Syst.
Comp. Conf., Wroclaw, Poland, 2014, pp. 1–9.

[13]
C. Zhao, E. Dall'Anese, and S. Low, “Convex relaxation of OPF in multiphase
radial networks with wye and delta connections,” in IREP, Espinho,
Portugal, 2017, pp. 1–7.

[14]
M. Bazrafshan, N. Gatsis, and H. Zhu, “Optimal tap selection of step-voltage
regulators in multi-phase distribution networks,” Power Syst. Comp.
Conf., pp. 1–7, 2018.

[15]
M. Usman, A. Cervi, M. Coppo, F. Bignucolo, and R. Turri, “Bus injection
relaxation based OPF in multi-phase neutral equipped distribution networks
embedding wye- and delta-connected loads and generators,” Int. J.
Elec. Power Energy Syst., vol. 114, no. Jan, p. 105394, 2020.

[16]
S. Claeys, G. Deconinck, and F. Geth, “Decomposition of n-winding
transformers for unbalanced optimal power flow,” IET Gener. Transm.
Distrib., 2020.

[17]
——, “Voltage-dependent load models in unbalanced optimal power flow using
power cones,” IEEE Trans. Smart Grid (prepublished), vol. -, no. -,
pp. –, 2021.

[18]
M. Vanin, H. Ergun, R. D'hulst, and D. Van Hertem, “Comparison of linear
and conic power flow formulations for unbalanced low voltage network
optimization,” Electric Power Systems Research, vol. 189, no.
September 2019, p. 106699, 2020. [Online]. Available:

[19]
B. A. Mather, L. R. de Araujo, S. Paudyal, T. E. McDermott, J. C. Fuller,
H. Zhu, W. Kersting, B. C. Pal, G. J. Shirek, L. F. Ochoa, J. L. R. Pereira,
R. C. Dugan, C.-W. Ten, K. P. Schneider, and S. Matthias, “Analytic
considerations and design basis for the IEEE distribution test feeders,”
IEEE Trans. Power Syst., vol. 33, no. 3, pp. 3181–3188, 2017.

[20]
C. Coffrin, H. L. Hijazi, and P. Van Hentenryck, “DistFlow Extensions for
AC Transmission Systems,” pp. 1–20, 2015. [Online]. Available:

[21]
K. Girigoudar, D. K. Molzahn, and L. A. Roald, “Analytical and empirical
comparisons of voltage unbalance definitions,” in North American
Power Symp., Wichita, Kansas, 2019, pp. 1–6.

[22]
J. Sliwak, M. Ruiz, M. F. Anjos, L. Létocart, and E. Traversi, “A Julia
module for polynomial optimization with complex variables applied to optimal
power flow,” 2019. [Online]. Available:

[23]
I. Dunning, J. Huchette, and M. Lubin, “JuMP: a modeling language for
mathematical optimization,” SIAM Review, vol. 59, no. 2, pp.
295–320, 2017.

[24]
J. C. Gilbert and C. Josz, “Plea for a semidefinite optimization solver in
complex numbers,” Inria, Paris, Tech. Rep., 2017. [Online]. Available:

[25]
A. Urquhart and M. Thomson, “Assumptions and approximations typically applied
in modelling LV networks with high penetrations of low carbon
technologies,” in Solar Integration Workshop.1em plus 0.5em
minus 0.4emLondon: Energynautics, 2013, p. 6.

[26]
Y. Liu, J. Li, L. Wu, and T. Ortmeyer, “Chordal relaxation based ACOPF for
unbalanced distribution systems with ders and voltage regulation devices,”
IEEE Trans. Power Syst., vol. 33, no. 1, pp. 970–984, 2018.

[27]
L. Gan and S. H. Low, “Chordal relaxation of OPF for multiphase radial
networks,” in Int. Symp. Circuits Syst., Melbourne VIC, 2014, pp.

[28]
R. A. Jabr, “A conic quadratic format for the load flow equations of meshed
networks,” IEEE Trans. Power Syst., vol. 22, no. 4, pp. 2285–2286,

[29]
M. Fazel, H. Hindi, and S. Boyd, “A rank minimization heuristic with
application to minimum order system approximation,” Proc. American
Control Conf., vol. 6, no. 2, pp. 4734–4739, 2001.

[30]
S. Kim, M. Kojima, and M. Yamashita, “Second order cone programming
relaxation of a positive semidefinite constraint,” Optimization
Methods Software, vol. 18, no. 5, pp. 535–541, 2003.

[31]
F. Geth and C. Coffrin, “Direct method to recover current and voltage in
multi-conductor optimal power flow models,” in IEEE PES General
Meeting, Atlanta, Georgia, 2019, pp. 1–5.

[32]
D. M. Fobes, C. Coffrin, F. Geth, and S. Claeys, “PowerModelsDistribution.
jl: an open-source framework for exploring distribution power flow
formulations,” Electric Power Systems Research, vol. 189, no.
December, p. 106664, 2020.