A Comprehensive Multi-Period Optimal Power Flow Framework for Smart LV Networks

This paper presents an extensive multi-period optimal power flow framework, with new modelling elements, for smart LV distribution systems that rely on residential flexibility for combating operational issues. A detailed performance assessment of different setups is performed, including: ZIP flexible loads (FLs), varying degrees of controllability of conventional residential devices, such as electric vehicles (EVs) or photovoltaics (PVs), by the distribution system operator (DSO) (adhering to customer-dependent restrictions) and full exploitation of the capabilities offered by state-of-the-art inverter technologies. A comprehensive model-dependent impact assessment is performed, including phase imbalances, neutral and ground wires and load dependencies. The de-congestion potential of common residential devices is highlighted, analyzing capabilities such as active power redistribution, reactive power support and phase balancing. Said potential is explored on setups where the DSO can make only partial adjustments on customer profiles, rather than (as is common) deciding on the full profiles. The extensive analysis can be used by DSOs and researchers alike to make informed decisions on the required levels of modelling detail, the connected devices and the degrees of controlability. The formulation is computationally efficient, scaling well to medium-size systems, and can serve as an excellent basis for building more tractable or more targeted approaches.

EV original active charge, phase z, period t, p.u. R ij,f θ Resistance of branch ij, between phases f , θ, p.u. X ij,f θ Reactance of branch ij, between phases f , θ, p.u.
Active power import/export, phase z, period t, p.u. S inv p,z,t PV inverter apparent power, phase z, period t, p.u. P inj p,z,t PV active grid injection, phase z, period t, p.u. Q P V p,z,t PV reactive power, phase z, period t, p.u. P Oc e,z,t EV active "overcharge", phase z, period t, p.u. P Uc e,z,t EV active "undercharge", phase z, period t, p.u. P ds e,z,t EV active discharge, phase z, period t, p.u. P D i,z,t Active load demand, bus i, phase z, period t, p.u. Q D i,z,t Reactive load demand, bus i, phase z, period t, p.u. Q I/E z,t Reactive power import/export, phase z, period t, p.u. u i,f,t Voltage magnitude, bus i, phase f , period t, p.u. σ up i,z,t Overvoltage violation, bus i, phase z, period t, p.u. σ down i,z,t Undervoltage violation, bus i, phase z, period t, p.u. σ ij,z,t Thermal violation, branch ij, phase z, period t, p.u.

I. INTRODUCTION A. Motivation
In adhering with the smart grid vision, distribution systems are transforming into evermore active systems, characterized by high shares of distributed energy resources (DERs), and high degrees of operational controllability by distribution system operators (DSOs) [1]. The proper management of such distribution systems is crucial; having been designed under the (now) archaic philosophy of fit-and-forget, they are usually illequipped to handle the uncoordinated, large-scale integration of DERs [2]. Especially in LV networks, which are inherently unbalanced, operational issues are usually more prevalent and of elevated severity. For the highest penetration levels, the stress inflicted on such systems can result in harmful voltage spikes or dips and damaging thermal loading of the distribution equipment, all of which are difficult to effectively contain [3].
To fully understand the benefits of residential flexibility resources (FRs), detailed models of the various devices and the distribution systems themselves are needed. Points of interest include the impact of the load modelling detail on the operational profile, the behavior of the neutral and ground voltages (for protection studies) and the interactions between differently loaded phases. However, most research works utilize simplified load and network models, or/and convex relaxations, targeting scalability rather than accuracy.

B. Literature review
Works that opt for solving exact (i.e., non-relaxed) formulations usually make non-generic, case-specific simplifications, such as ignoring the neutral wire or assuming small load imbalances to name a few. The seminal paper on multi-period optimal power flow (MP-OPF) for active distribution systems, [4], employed the single-phase (1Φ) network representation and the constant P/Q load model. While subsequent papers have since introduced more advanced models for both MV and LV systems, see [5], most papers prioritize the solution technique rather than the model, assuming that non-generic, case-specific simplifications are always expected to hold.
The MP-OPF problems between the MV and LV level share some conceptual similarities, such as the radial network structure, the unbalanced system conditions or the non-negligible impact of line resistance/reactance. In MV systems, the DSO has various resources at its disposal, such as capacitors banks, network switches (reconfiguration), distributed generators and tap changers, see [6]- [8]. However, the situation is very different in LV systems, where the DSO has far less controllability and equipment available. System management is achieved primarily though electric vehicles (EVs) and photovoltaics (PVs), and rarely through energy storage (ES) systems.
In terms of MP-OPF features, the authors of [9] eliminate the neutral phase through Kron reduction and calculate an optimal control strategy through an iterative approach (constant P/Q load). The technique is also used in [10], where the current-based formulation is employed instead (simultaneous study of MV and LV network). The work [11] proposes a current-mismatch MP-OPF to optimize a feeder's operation, ignoring the grounding and assuming full controlability of residential ES systems by the DSO (constant P/Q load).
In [12], a local EV charging strategy is applied with respect to the PV's operation, though assuming that the EV is available at all times (constant P/Q load, 3-phase network). A multi-period approach is proposed in [13] for designing the entire charging strategy of some EVs, subject to to dynamic electricity pricing (current mismatch formulation, 3phase network). The authors of [14] combine central and local control strategies for managing distribution systems through PV utilization (constant P/Q load, neutral consideration).
Approaches based on PVs and EVs that are more specialized have also been proposed. For example, in [15], the 3stage, centralized, single-period PV curtailment and reactive management problem is solved for 24 consecutive hours for a four-wire LV distribution system (constant P/Q load). The authors of [16] employ three-phase (3Φ) inverters for phase balancing in pure 3Φ networks (constant P/Q load). The usage of ES is more rarely tackled, due to their low penetrations in LV networks and the undesirable level DSO involvement in household equipment management. In [17], the centralized control of ES in unbalanced networks is studied, though their temporal constraints are ignored (neutral consideration, constant P/Q load). The authors of [18] present a highly detailed ES model specifically for LV systems and propose a local area control strategy involving several network agents (Kron reduction, constant P/Q load).
While the constant P/Q load model is the one most commonly employed, more intricate models have occasionally been explored. Explicit ZIP coefficients for commercial, residential and industrial loads have been proposed for conservation of voltage reduction (CVR) studies in unbalanced distribution systems [19], though this customer variety is only found in MV systems. A common ZIP load structure has even been employed in studies following model-free approaches of distribution system optimization, see [20] for example (phasor regulation strategy), with the work assuming the simplest form of an unbalanced system. Comparisons between different load models have been performed, though these are either on a purely technical level, i.e., behavior of single, standalone load [21], or comparisons between the standard ZIP model and different approximations of it [22]. Comprehensive studies of all possible ZIP structures as they related to optimizing distribution system behavior are lacking from the literature.
In terms of solution techniques, convex relaxations are often employed for multi-phase networks. The second-order cone programming (SOCP) and semi-definite programming (SDP) relaxations are popular choices that have spawned different variations (based on types of connection and imbalance), though the neutral wire and ground are rarely included [23], [24]. The branch flow formulation is also commonly employed in radial systems [25]. While extendable to unbalanced systems, it usually includes power flow linearizations or assumptions of very small imbalances [26], [27]. However, most relaxations hold very rarely for realistic power systems [28].
Approximations for specialized versions of the unbalanced MP-OPF have also been proposed. The authors of [29] employ the SDP relaxation for distribution systems with neutral cables and fixed ZIP loads. The work [30] proposes linear and quadratic simplifications of an MI optimization problem that addresses different operating conditions between MV and LV systems (exponential ZIP load). Aside from the simplifications, these formulations are not entirely practical in representing the behavior of realistic distribution networks.
In recapitulating, all proposed approaches (simplified or not) do produce promising results. However, if the solved problem is a relaxed one, the original system is largely simplified and the results are rarely feasible/reliable. On the other hand, in exact approaches, the neutral and ground cables are usually ignored, while the effects of the mutual couplings and the load types are not always properly represented.

C. Contributions & Paper structure
To focus on scalability, most works sacrifice some accuracy, leaving the behavior of modern (vast flexibility options array) and realistic LV systems in MP settings insufficiently addressed. This work develops a comprehensive, versatile and easily reproducible MP-OPF tool for unbalanced distribution systems and realistic device models.
This work draws inspiration from and significantly extends several past works. The current-based formulation is a combi-nation of [9], [10], [31], extended to include the interactions with the neutral and ground wires and the rotation of each phase to a common reference plane. The load modelling is adapted from [5], also accounting for partial load flexibility. The PV modelling is adapted from [5], [9], also including limited curtailment capability and a realistic production capability curve (PCC). The modelling of EV flexibility is original. For the phase balancing through 3Φ inverters, this work extends [16], in accounting for the re-scheduling of the user-driven charging profile. The proposed approach for remunerating customer participation is also a novel addition.
On top of the novel modelling elements, the more conceptual contributions of the paper can be summarized as follows: • It provides an extensive analysis of a comprehensive set of available modelling decisions (many often disregarded) for the optimal management of LV distribution systems. • It construct a generic MP-OPF model that can provide valuable information on how to unlock the full flexibility potential of common LV networks. The framework is easily reproducible, adaptable to each researcher's needs and ideal as a basis for more sophisticated developments. • As the topic is underaddressed, the paper proposes a first crude DSO/customer collaboration framework through which the DSO can utilize residential devices to achieve better system management. The importance of the framework is paramount, given the DSO's traditionally minimal involvement in managing LV systems. The great advantage of the developed tool is its adaptability to the specific needs of each problem (load/network/device model). It can be used under a vast array of optimization setups, including active power redistribution, reactive power management schemes or phase balancing. This offers great insight to DSOs, which can explore the behavior of often neglected system parts in a reliable manner, and unlock the network's full potential for residential flexibility utilization. Informed decisions can be made on which modelling elements are required and which can be reliably ignored.
The remainder of this paper is structured as follows. The problem formulation and the main assumptions are presented in Section II. The case study is extensively analyzed in Section III. Conclusions are drawn in Section IV.

A. Problem assumptions
For the sake of clarity, we lay out the main problem assumptions. This is day-ahead (planning), multi-period (24-hour horizon, hourly resolution), centralized control optimization problem, where the DSO has partial controllability of the available flexibility resources (FRs). The work is contained to the deterministic setting, assuming the DSO uses a mostlikely-to-occur forecast scenario for its planning (real-time deviations are addressed on-the-spot, though this is out of scope). We assume the existence of a digital platform through which customers inform the DSO of their ideal controllable device schedules, initially designed through a rule-based approach or by a sophisticated software (e.g., energy management system). Each customer may or may not receive new set-points for their controllable devices. The expected difference (not the realtime deviations) between the original customer profile and the designed, post-request profile, is the basis for the customer's remuneration. All the necessary software is pre-installed. The generic formulation is applicable for most setups.
This work covers radial, LV distribution systems. Given the proposed framework's generic structure, it can be used to model 4-wire, multi-grounded systems (found in North America, Europe), 3-wire, grounded or ungrounded systems (found in Europe, UK), single-wire with earth return (found in Australia), and cases of highly specialized grounding (high resistance/reactance, Petersen coils) [32]. More intricate configurations, e.g., 5-wire systems, are out of scope.
This work fits within the generic MP-OPF formulation originally developed for 1Φ networks representations in [4]. The proposed formulation employs the rectangular coordinates formulation for voltages and currents (1)-(2):

B. Objective function
The DSO plans an hourly cost-optimal control strategy for the various FRs, to maintain acceptable conditions across its system. The objective (3) is composed of the following costs: import/export (4), PV utilization (5)-(6), FL alteration (7), EV profile re-design (8) and limit violation (9): C P P V = c P P V · t,z,p pf p,z,t · (S gen p,z,t − S inv p,z,t ) (5) t,z,l (P Od,1 l,z,t + P Od,2 l,z,t + P Ud,1 l,z,t + P Ud,2 l,z,t ) (7) Active power imported/exported from/to the MV level incurs a proportional cost for the DSO (4). PV active power curtailment is highly priced, while the cost for utilizing reactive capabilities is tied to the active power that is not injected due to producing reactive power instead (5)- (6). Increasing or decreasing the customer-forecasted load demand has a proportionally associated cost (7). Deviations from an EV's customer-desired profile, or discharging, are proportionally priced (8). Technical limits violations carry a high penalty (9).
Each cost is chosen based on the desired activation priority order (PO). For example, the less desirable PV curtailment has higher associated costs than load profile alteration. The applied costs (Table I) simply reflect a certain PO. However, the framework is applicable with any (more realistic) costs.

C. Device modelling
For all subsequent equations, for any device a (load, PV, EV), the symbol u a,z,t refers to the voltage difference that concerns device a, originally connected to phase z, at time period t: for 4-wire systems, the difference is between phase and neutral (Wye) or phase (delta); for 3-wire systems, it is between phase and ground (Wye) or phase (delta).
1) Customer loads: The active and reactive demand of a load is based on fixed impedance, current and consumption elements, i.e., the ZIP model (10)-(11), represented by their respective consumption percentages (12). A power factor (pf ) of 0.95 is assumed for all loads. However, the user may define any pf they deem appropriate. The nominal active and reactive powers, P 0 l,z,t , Q 0 l,z,t , can thus be interrelated through (13), where Ω represents the type of pf (1 for lagging, -1 for leading). The load currents are calculated based on (14)- (15). The following hold ∀l ∈ L, ∀z ∈ Z, ∀t ∈ T : A customer load may be flexible (FL) in some ways, in which case its nominal power can be generalized, as shown in (16). Highly accurate models split the FL into 4 parts: fixed consumption, P f x l,z,t , fully flexible consumption, P f c l,z,t , energy (E) shiftable consumption, P es l,z,t and time (T) shiftable consumption P ts l,z,t . The first two parts are most commonly considered in the vast majority of cases. The fully flexible consumption is complemented by the variables P Od l,z,t , P Ud l,z,t , i.e., "overdemand" and "underdemand", respectively. "Overdemand" means higher consumption than originally forecasted U: Upward, D: Downward and vice versa for "underdemand". The two variables are limited by a maximum alteration, M F L , and designed to not simultaneously be positive (sub-optimal, see [4]). While they are less common, the last two parts are also added for the sake of model comprehensiveness. The third part of the FL is energy-shiftable, i.e., in-day alterations are allowed (similarly to the second part), but their net sum must be zero (fixed daily energy consumption), (18). The final part of the FL is time-shiftable, i.e., it can be moved through time, though unaltered. This is modelled through the use of binary variables (δ l,t ), see (19) which dictates that the fixed part, P SL l,z,t , must be active exactly for A time periods (akin to cycle time). The resulting MINLP problem is generally intractable and requires the use of approximations or heuristics to effectively manage, see [33]. Common: E-shiftable: 2) Photovoltaics: The apparent power generated by the applied solar irradiation, S gen p,z,t , can be curtailed up to a certain percentage, M P V (20), resulting into the apparent power at the inverter level, S inv p,z,t . The PV's power factor (pf) is partially flexible, allowing for the utilization of a PV's reactive power (21)- (22). Said flexibility is limited by the PV's traditional PCC, as dictated by (23). The injected PV currents are calculated through (24)- (25). When the PV is not connected to a particular phase, all corresponding variables are equal to zero. The following hold ∀p ∈ P, ∀z ∈ Z, ∀t ∈ T : The above hold for 1Φ PV inverters. However, using 3Φ inverters instead allows for redistributing a PV's output among the three phases [16]. For simplicity, we ignore the PCC (unity pf ) when phase balancing is possible; nonetheless, per-phase formulations of (21)-(23) can be defined. The total injection among the three phases must be equal to the original 1Φ injection (26). There is a limit to the redistribution of active power, based on the rate of each 1Φ inverter, P P V,rate p , (27). An individual inverter may consume active power, hence the possible negative values for 1Φ PV injection. This can be fully exploited at night, when PV production is zero [16].
3) Electric vehicles: An EV can "overcharge" or "undercharge" with respect to its originally forecasted profile (similar logic as in FLs), or discharge if the V2G capability is available. Since V2G is a much more sought-after service, its cost is much higher than that of profile rescheduling. This design ensures that the "undercharge" capability (i.e., charging modulation) will be first fully utilized before before discharging can actually be activated (see also [4]). The EV can only interact with the grid when the resident is home (28). The EV must "conclude" the optimization horizon under the same status that its owner originally intended (29)-(30), where n c , n d are the charging and discharging efficiencies. Constraint (32) dictates the technical limits of the EV's charging/discharging behaviour. The binary parameter ξ EV represents whether the EV can only charge (ξ EV = 0) or if it has the V2G capability (ξ EV = 1), in which case it also dabbles as a traditional energy storage system. Constraint (33) ensures that the EV does not charge/discharge above/below certain energy limits, E cap e , E min e . The EV currents are calculated by (34)- (35). EVs are assumed to operate with a fixed pf . Advanced EVs may even be capable of employing 4-quadrant control [34], though such instances are rare; (21)-(23) could be easily adapted and employed if need be. When the EV is not connected to a particular phase, all corresponding variables are equal to zero. The following hold ∀e ∈ E, ∀z ∈ Z, ∀t ∈ T : P Oc e,z,t = P Uc e,z,t = P ds e,z,t = 0 ∀t ∈ T nc e (28) P demand e,z,t = P C e,z,t + P Oc e,z,t − P Uc e,z,t − ξ EV · P ds e,z,t (P Oc e,z,t , P Uc e,z,t , P ds e,z,t ) ≥ 0 (30) P demand e,z,t = u re e,z,t i EV,re e,z,t + u im e,z,t i EV,im e,z,t Q EV e,z,t = u im e,z,t i EV,re e,z,t − u re e,z,t i EV,im e,z,t EVs may also be equipped with 3Φ inverters [16]. Constraint (32) is modified, due to the limit on the redistribution of active power, based on the rate of each 1Φ inverter, P EV,rate e , (36). While power can be injected in a phase, the EV type is essentially the same (net charger); pure net discharge (V2G) is not possible, as this would involve a different (and far more expensive) kind of EV. In addition, the variables and parameters P C e,t , P Oc e,t , P Uc e,t are now only considered for the EV as a whole (no z index). (29) is thus modified as (37):

4) Stand-alone battery:
The EV model is easily adaptable to a stand-alone battery model. This is achievable through the following modifications, where P D e,z,t is the original, userdriven discharging profile, and P Ods e,z,t , P Uds e,z,t are the "overdischarge" and "under-discharge" variables, respectively:

D. Low voltage network technical constraints
The distribution line model depicted in Fig. 1 is adopted. Each "phase" f ∈ F is characterized by its self-impedance, Z f f and its mutual coupling with other "phases" θ ∈ F : f = θ, Z f θ . A grounding impedance, Z gr , may also be present. The neutral and ground currents are calculated using current dividers. As is common for LV networks, devices are assumed to be wye-connected, though the adaptation to delta connections is also presented. For neighboring nodes, currents are assumed equal at origin and destination (small line shunts [35]).
For the LV network, we have the current injection balance (42)-(43), where (42) concerns wye-connected devices and (43) details the necessary modification to address delta connections (the superscript demand refers to any considered device, e.g., PV, E, FL, and z, z * , z * * represent different phases). The modification is necessary, since employing the conventional Y → ∆ transformations would only lead to an approximated MP-OPF solution. We also have the matching of the injections of the three phases with the neutral phase (44), Fig. 1: 3Φ, four-wire, multi-grounded distribution line [35] where binary parameter φ ∆ represents the existence of a delta connection, the constraints that guarantee the current balance at common ground points (45), the voltage drop across lines (46)-(47), the voltage and branch currents technical limits (48)-(49). The above hold both for the real and imaginary parts. Note that all slack variables are positive (50). No constraints are enforced for the neutral and ground "phases" since they are completely dependent on phases a, b, c. The balances between the currents at the neutral-ground connection are enforced by the current dividers; no additional constraints are necessary to capture this behavior. The following hold ∀i, j ∈ I : i = j, ∀z ∈ Z, ∀f ∈ F , ∀t ∈ T : At the feeder head, voltages are defined as {1∠0 • , 1∠ − 120 • , 1∠120 • } → {1 + 0j, −0.5 − 0.87j, −0.5 + 0.87j} for phases a, b, c, respectively. While constraint (48) is perfectly valid, a computational inefficiency stems from the fact that the phases are rotated with respect to each other; the rectangular modelling of voltages causes constraint (48) to be non-convex. Assuming that the voltage angles (per phase) diverge only slightly from their reference value, a convex reformulation is employed, similarly to [9]. This allows overcoming one of the several non-convexities, while making the evaluation of the results far easier and more intuitive. The three phases are • } so that they lie close to the reference axis 0 • , and the same feasible space is defined for each [9]. A visual representation is presented in Fig. 2. Constraint (48) is thus re-defined as: A crucial point not addressed in [9] (as the paper assumed a perfectly grounded system and consequently used the Kron reduction to remove the neutral cable) is that the phase rotation affects the interactions between phases a, b, c and the neutral wire (44). These constraints must be updated to account for the voltage rotation performed in (51). As such, the calculated current for each phase is counter-rotated by ROT . If we define the total current injections at node bus i, phase z, period t, (52), then (44) are updated as (53)-(54):

F. Remarks
The problem at hand is an MP-OPF. The conventional setup (1Φ inverters, no voltage shifting) is comprised of the rectangular coordinates modelling and the various costs (1)-(9), the modelling of FLs (10) The physical problem modelled is an NLP, regardless of the modelling choices concerning the network and the residential devices. The nonlinearity stems from the ZIP load models, the PCC of PVs and the relations between power and voltage/current. The complex formulation, including all the peculiarities of MP-OPF in multi-phase LV networks, is not directly amenable to exact convex formulations (particularly SDP/SOCP). The authors prefer adopting an NLP formulation and solver, in order to obtain a feasible and at least local optimal solution, as compared to the (high) risk of obtaining physically meaningless solution from a relaxed problem [28].

A. Simulation environment
The proposed formulation is applied on the 18-node, modified CIGRE LV benchmark distribution network [36], depicted in Fig. 3. The original network dataset does not include data regarding the neutral wire and the ground; artificial values are used, based on observations regarding the relationship between phases a, b, c and the missing values, as derived from [35]. When grounding is included, a value of R g = 1Ω (adjusted in p.u.) is considered for all customer nodes, though other kinds of grounding are also possible (e.g., grounding impedance or Petersen coil). All other nodes are assumed to be ungrounded. Base power and (3Φ) voltage are chosen as 1kW and 240V , respectively. The Z/I/P percentages are randomly assigned.
The connected distributed devices throughout the network are FLs, EVs and PVs (connection points shown in Fig. 3). The characteristics of each device are available in Table II. Their actual 24-hour profiles are available in [3]. While the framework can accommodate any residential device, in order to better examine the impact of different flexibility setups, devices with high power rates are purposefully selected. All simulations are performed on a PC of 2.7-GHz and 8-GB RAM, using the general purpose NLP solver IPOPT [37], with default settings, through GAMS [38]. The framework was previously evaluated (via power flow comparisons) on a number of MV and LV networks. The errors (in p.u.) on complex voltages varied between 10 −4 and 10 −6 .

B. Modelling and operational scenarios
Apart from proposing a comprehensive MP-OPF framework for realistic distribution systems, this work also provides infor-   Table V. Focusing on smart distribution feeders with high penetration of potential FRs, the main purpose of this work is reflected in Cases 1 and 2. That purpose is to provide a comprehensive understanding to DSOs with regards to how their system realistically behaves, how the accuracy of the solution is affected with cumulative simplifications/approximations, and how practical each level of modelling detail actually is. The inclusion of Case 3 is done in order to demonstrate to DSOs how realistic distribution systems would behave if the full potential of the various interconnected FRs were unlocked.        Table VI. It becomes immediately apparent that the (often neglected) load modelling choices affect the solution in a non-negligible way (though solution times and total violations are comparable). When the Z-component is dominant, the network's voltage profile tends to rise. Contrariwise, when the P-component is dominant, the voltage profile tends to drop. The I-component appears to have the smallest impact on the voltage profile. Similar observations are much more prevalent in nodes that are further away from the substation.
A higher average demand does not necessarily correspond to a reduced voltage profile. When there are binding interactions between the load model and the voltage profile, it is not always straightforward to estimate how the two will be affected. Formulations that include the P-component (higher average demand) such as ZP, IP, P, tend to show pictures of higher stress for the system, requiring more reactionary measures in response (i.e., higher operational costs). Formulations that can produce a more "malleable" demand such as ZI, Z, I, have the opposite effect (lower costs on average). The smallest, almost negligible difference in the objective is observed under the ZP load model; the I-aspect of loads could perhaps safely be dropped from simplified formulations. Nevertheless, less conservative load models (assuming there is sufficient data to support them) can bring about operational benefits (more accurate coordination of flexibility options). Do note that while the differences (not necessarily errors) in objective values are not very high (on average about 3.3%), the cumulative impact of the load modelling decisions are certainly non-negligible.
2) Case 2: Results for various levels of network modelling accuracy are presented in Table VII (assuming N 1 is the ideal objective). Simpler network models produce inaccuracies in the voltage profile (elevated or reduced levels), leading to significant overestimations or underestimations of operational issues (ignoring imbalances, mutual coupling or both produces an average error of 29.25%). Models N 4 , N 5 have unrealistically low violation costs, despite exhibiting vastly different voltage profiles. While simpler models have benefits (speedwise), they could be inappropriate for many applications.
When the neutral wire is explicitly modelled the error drops to approximately 1.25%. The large errors observed are not unexpected. As was observed in [17], load imbalances of about 50% can lead to the corresponding neutral wire being responsible for approximately 20% of system losses. In our case, some nodes have 100% imbalance, with the neutral wire carrying very high currents. While the neutral voltages are very low (almost always below 0.1 p.u.), the coupling effects of the large neutral currents create several issues. The Kron reduction is acceptable only on systems with perfectly grounded neutral wires; this is not our case, hence we observe large errors. Depending on the level of imbalance, the modelling aspect that could be reliably dropped is the grounding.
3) Case 3: Results for various levels of FR controllability, using both 1Φ and 3Φ inverters, are presented in Table VIII. Expectedly (for the 1Φ case), as more flexibility is added to the DSO's "toolkit" the operational conditions improve. In our case, increased PV controllability is followed by increased FL controllability, finally adding EV controllability into the mix. The voltage profile remains relatively steady, though the upper and lower values are improved, resulting in reduced violations costs. In addition, the available residential flexibility also reduces the import costs, as an increasing amount of insystem energy is re-distributed, decreased or converted (e.g., PVs changing their pf to avoid curtailment).
For the 3Φ inverters case the IE price is increased 10 times in order to motivate system self-sufficiency, discourage (as much as possible) interactions with the MV level and maximize the utility of residential FRs. Even without employing any flexibility, the usage of 3Φ inverters improves the voltage profile; even for S 1 , we observe total violations of almost zero. The usage of 3Φ inverters offers tremendous amounts of 3Φ case: PV pf = 1, IE price increased tenfold to motivate system self-sufficiency. *Lower voltage limit set to 0.95 p.u.
flexibility to the operator (despite the fact that we have also removed the reactive capabilities of PVs), without the need to actually engage in profile alterations. A thing of note is that when the connected devices are able to distribute their profiles between the three phases, the neutral wires carry almost no current and as such most of the original active losses are no longer observed; had the original IE cost been kept then the objective function value (as compared to the 1Φ inverters case) would have been reduced by about 55%.

D. Assessing scalability 1) System vs customer base expansion:
Similarly to most NLP problems, full scalability for very large problems is not guaranteed. Nonetheless, the formulation is in general computationally efficient for small/medium-size systems and sometimes even for large systems (see Section III.D.2). For the 18-node system, despite the already large problem size and the variety of available flexibility options, a solution is achieved (on average) in less than 30 seconds.
The proposed formulation is applied on a partial reallife distribution feeder, composed of 180 nodes hosting 23 residential customers (13% customer-to-node ratio, see [3] for original), of which 5 own EVs, 8 own PVs and 10 own both (sizes of Table II). The system is assumed perfectly grounded and as such, the neutral wire can be reliably reduced. Five  Fig. 4 and Table IX.
The results illustrate some interesting points. While the problems sizes are comparable 1 , the 3Φ setup has the superior performance. After increasing the degree of stressing of the system per step (indicated by the increasing objective value), an optimal solution without violations is always achieved, with the solution speed increasing, but not significantly. Contrariwise, the 1Φ setup cannot cope efficiently with high levels of stress, requiring more time to reach a solution (on average 270% slower). In fact, to avoid infeasibility for the H case, much more residential flexibility is "activated" than would be reasonably expected. The 3Φ setup provides far better solutions (on average 54.4% better), due to not directly utilizing flexibility per se but rather by alleviating the negative effects of load imbalances through consumption redistribution amongst phases. The above are indicative of the superiority of "investing" in flexibility that is characterized by quality (phase balancing) instead of quantity (profile alteration).
The solution is calculated faster, owing to better tailored initialization and elimination of variables. Concerning the purely technical aspects, the 1Φ case has higher costs, owing to to the extensive re-shaping of FL/EV profiles, as well as engaging in extensive PV curtailment and reactive power injection. The 3Φ case makes extensive use of the EV/PV inverters (almost all are utilized), avoiding the need to resort to flexibility procurement almost entirely. As such, the feeder is managed much more cheaply and efficiently. In both cases, however, the DSO has access to a very large pool of residential flexibility; the subsequent technical issues are minor and the feeder is (nearly) always operated within acceptable conditions.
2) Further expansions and limitations: The problem size expansion can follow two directions: system expansion (nodes) and customer base expansion (controllable FRs). As such, we must explore the framework's performance with respect to two kinds of expansion. Four additional feeders (originals available in [3]) of 200, 400, 600 and 800 nodes, respectively, were examined (using the 3Φ inverter configuration for PVs/EVs), hosting different sizes of customer bases. A single simulation is performed per case. The customer distribution (node/phase) is random. PF-based initialization was again employed. Due to technical issues encountered (solver crashing without clear reason), the commercial solver KNITRO [39] is employed instead of IPOPT. The results are presented in Table XI.
As is obvious, some of the examined setups results in very highly loaded systems, requiring the "activation" of high percentages of FRs, thus stressing the solution process itself. As expected, increasing the system size (nodes) subse- Int: Solution manually interrupted after 10,000 seconds quently increases the solution time. However, if the number of controllable elements remains low (lower system stress, less variables), a locally optimal solution is (generally) achieved within acceptable time-frames for day-ahead settings. For very high penetrations of controllable elements the system is much more stressed, increasing the solution time. In fact, for larger systems hosting unrealistically high numbers of FRs, no solution is returned even after 10,000 seconds. However, it appears that the negative impact on the solution time stems on a larger part from the number of controllable elements, rather than from the size of the system (very large system with no controllable elements are simulated very fast). For most LV distribution systems (10-30% customer-to-node ratio), a solution can be calculated within reasonable time-frames.
On a final note, the authors wish to re-stress that despite the several pros of the framework in and of itself, no claim is made on the quality of the chosen modelling/simulation tools. However, the exact actions to be taken for fine-tuning the solution process are out of paper scope.

E. Additional cases of interest
The results so far demonstrated the capabilities of the framework for the most common applications of optimization in LV distribution systems. However, for the sake of comprehensiveness, we re-turn our focus to the 18-node system and examine the impact (on the solution) of five more intricate modelling choices: 1) Specialized load models: While the traditional flexible load model (17) is the one most commonly used, it is still useful to perform an impact analysis for the expanded (generic) FL model. For that purpose, we consider FLs as the only controllable elements. We examine three different FL models, where 50% of the load is fixed and the remainder 50% is either fully flexible, energy-shiftable or time-shiftable. A fourth model (most generic) is also examined, where 25% of the load is fixed and the remainder 75% is equally distributed to the first three models. Results are presented in Fig. 5.   5 illustrates through various metrics the percentage difference between using the most generic FL model and different versions of it. Obviously, more refined models translate to higher FL costs, albeit to somewhat lower objective function values and total technical violations. Specifically, the first three load models result in 20%, 19% and 5% more violations, respectively, than their combination, while the total FL cost is 76%, 47% and 22% lower, respectively. However, the objective function values are practically the same for all FL models. Whatever small improvement is achieved comes at a cost of drastically higher solution times, with the first three models reaching an optimal solution 72%, 61% and 23% faster, respectively. This is an important reason why the first FL model is the one most commonly employed. Reducing the solution time would require specialized approximation (a necessity for much larger networks), which are out of paper scope. The reader is referred to [33] for further details.
2) Modelling of multi-phase lines: While the presented results focused on the most commonly employed network models, i.e., either all-1Φ or all-3Φ, the proposed framework is fully compatible with multi-phase networks as well. An example of the framework's compatibility is the solution of the MP-OPF problem on the CIGRE LV feeder, modified to an arbitrary multi-phase system, see Fig. 6. As can be seen, while the main body of the feeder is 3Φ, the lines directly leading to buildings are not. Specifically, the more active buildings (hosting FLs, EV, PV) are connected to the main feeder through 2Φ lines, while 1Φ lines are used for less active buildings. All devices are accordingly modified so that their phases correspond to these of their respective lines. Modelling a multi-phase line is simply achieved by disregarding some of the parameters R l,f θ , X l,f θ , i.e., omitting them from the formulation after network topology processing. The phase allocation per line is completely random.
The random allocation of phases and solution of the MP-OPF was performed a number of times to validate the framework's capabilities. The main observation was that the MP-OPF always reached an optimal solution without issue. It is however worth stating that, on average, the solution time was slightly higher than that for pure 3Φ networks. This was expected, since when the number of phases is reduced the  power is not distributed as well, resulting in higher thermal stressing for some lines. As previously discussed, when the system is more stressed, i.e., requires the "activation" of additional slack variables, the solution time is generally higher. Nonetheless, this is a natural outcome of the problem setup and is not attributed to some weakness in the framework.
3) Validation of proposed stand-alone battery model: The proposed stand-alone battery model is validated on the 18node feeder. The battery, originally based on the proposed EV model, has a standard user-driven profile, sometimes charging (P C ) and sometimes discharging (P D ). Contrary to common battery models, e.g., the one presented in [4], instead of "penalizing" all charging and discharging activities, only the deviations from the user-driven profile are remunerated. As can be seen in Fig. 7, which presents the behavior of a battery at node 17, the battery originally charges at noon and discharges during the night, heuristically designed to do so by the user. Post-optimization, the battery charges more heavily at noon to counter the high PV production (avoid overvoltages), discharges more heavily at night to serve the high EV demand (avoid undervoltages). The extra stored energy is distributed to early morning hours. The simulation time remains small (comparable to what would be expected from a traditional battery model), and the battery behaves within the lines of "normalcy". Even after optimization, the battery exchanges the exact same energy amount as originally designed, a strong point of the proposed model.

4) Wye-connections vs delta-connections:
At this point it is still useful to examine how the solution is affected by the device connection type (remember that the framework is fully compatible with both connection types). For that reason, we examine three cases: the standard case, where only Y connections are assumed, the uncommon case, where only delta connections are assumed, and an intermediate case, where both connection types are assumed. All three cases are simulated through the proposed framework, and results are presented in Fig. 8. Delta-connected devices generally lead to higher currents (total demand), resulting in slightly higher objective function values. Systems with more delta connections may also require more time to solve, since not only do the constraints become more intertwined, but the increased demand also leads to higher system stress, which, as

5) Evaluation of the V2G capability:
Though not as common, some EVs also have the V2G capability. This capability has limited impact during weekdays: since the EV is not present for a big part of the day (owner is away), this cripples the battery's potential for providing flexibility, as there are no clear opportunities for it to charge and serve any possible need for self-consumption during nighttime. On the other hand, the consumption re-distribution capability has limited potential during weekends, where the limited need for any charging severely constrains the potential of 3Φ inverters. The aforementioned observations are confirmed by simulations with both inverter configurations. The best case scenario would be to consider an EV with a 3Φ inverter configuration and the V2G capability, thus this setup is rare itself.
To make a fair comparison, we examine an idealized setup where the EV must charge its usual amount, but without temporal restrictions, i.e., the EV may interact with the grid at all times. We compare the performance of the V2G capability (1Φ configuration) against the consumption re-distribution capability (3Φ configuration). We also set the exporting price to 20 e/kW to motivate system self-sufficiency (see also Case 3). As can be seen in Fig. 9, the 3Φ inverter configuration beats the V2G capability by a clear margin. Under the former, the average per-phase voltage level is generally more elevated and fluctuates more closely to unity. In addition, the remuneration cost with the 3Φ inverter configuration is significantly lower. When the V2G capability is employed, the phase to which the EV is connected must engage in drastic action (charging alteration and discharging), utilizing two different flexibility services (thus driving up the cost) and only partially managing the problem. With the 3Φ inverter configuration, the EV must only alter its charging profile. By constantly balancing its operation between the three phases, the flexibility cost is kept No single phase disproportionately affects the others, as the 3Φ inverter configuration ensures that the interactions between phases is more balanced and more harmonious.

IV. CONCLUSIONS
The authors have constructed a versatile MP-OPF framework that serves two key functions: it proposes and compares the performance of state-of-the-art device models for unlocking the flexibility potential of smart distribution grids and it provides up-to-date guidelines with respect to how each of the most commonly employed (or ignored) modelling choices (concerning the loads, the network and the degrees of controlability) affect the quality and reliability of the solution. The reported results can guide researchers into picking proper equipment models depending on their respective needs. The formulation also scales well for larger systems (under proper conditions). As such, it can serve as a solid basis for approaches aiming specifically at scalability; its results can also be safely contemplated by the DSO for hours-ahead use.
The impact of different versions of common FR devices was analyzed, based on novel and realistic models. Namely, the authors proposed flexible ZIP load models (profile alteration affects Z, I, P components), both 1Φ and 3Φ (balancing) versions of PVs with realistic associated costs, reactive capabilities and PCC curves and both 1Φ and 3Φ (balancing) versions of EVs with the added novelty of building on the original customer-desired profile, rather than determining an original profile altogether. The proposed models can be used in devising new approaches for unlocking the full potential of FRs to manage violated constraints in distribution systems.
In the future, the authors plan to extend their MP-OPF framework to to address the issues of uncertainty, an area were, due to the huge computational challenge, the research is scarce and often limited to small systems.