A Two-Stage Solution to the Bi-Objective Optimal Voltage Regulation Problem

In this paper, a two-stage approach is introduced to efficiently solve the optimal voltage regulation problem in radial medium-voltage (MV) distribution networks. The proposed method uses the available reactive power of distributed generation units and the on-load tap changer (OLTC) of the high-voltage/MV transformer to regulate network voltages, while also minimizing two conflicting objectives, namely the network energy losses and the frequency of tap changes. This bi-objective optimal voltage regulation problem is addressed in two distinct stages. In the first stage, the network is linearized and a simplified optimal voltage regulation problem is solved to determine the candidate OLTC operating plans (COOPs). For each COOP, a new reactive power allocation method is employed in the second stage to regulate the network voltages and minimize network losses. This method follows a rule-based approach, also allowing its implementation under real-field conditions. The final outcome of this two-stage process is the Pareto-front which can be used by the distribution system operators to select the most preferable solution for the real-time network operation. The proposed methodology is characterized by reduced computational complexity compared to the application of conventional optimization techniques. Time-domain and time-series simulations on a radial MV distribution network are employed to thoroughly evaluate the performance of the proposed method.

DSOs can easily handle this problem using all the available network elements [3], their uncoordinated operation may increase the stress on the grid by means of both excessive operation of the network devices and high network losses [4]. Thus, to improve the overall network performance, various optimization objectives can be considered, such as loss minimization. The problem is then transformed into an optimal voltage regulation (OVR) problem, which is a newly introduced research topic for distribution networks, characterized by increased complexity [5].
In the literature, several optimization-based methods have been proposed to address the OVR problem in medium-voltage (MV) distribution networks. The vast majority of these strategies use the reactive power capability of DG units in combination with the basic network elements, i.e., the on-load tap changer (OLTC) of the high-/medium-voltage (HV/MV) transformer, feeder capacitors, and series voltage regulators (SVRs), to regulate the network voltages in an optimal way. More specifically, the authors in [6] propose a distributed method to optimally reallocate the reactive power among the DG units in order to minimize network losses. An additional operation constraint is introduced in [7] to limit the daily number of the tap changes. Furthermore, an offline coordination of the Q(V ) droop characteristics is proposed in [8] to minimize the overall reactive power of the DG units. Nevertheless, the above-mentioned methods attempt to solve a single-objective optimization problem, focusing on the coordinated operation of the DG units and neglecting possible interactions with other network devices, such as the OLTC [9].
To deal with this problem, multiple optimization objectives can be considered within the OVR concept. Among a wide range of optimization objectives, the most well-established for distribution networks include loss minimization [10], voltage profile improvement [11], and OLTC switching operation minimization [12]. However, the objective of voltage profile improvement can be disregarded, since it does not provide any added value compared to the simple use of voltage limits constraints that is already foreseen in the OVR concept. Consequently, the problem is simplified to a bi-objective optimal voltage regulation problem. Nevertheless, conflicts between the remaining objectives may occur in case of highly variable generated power due to the intermittent nature of DG units consisting mainly of renewable energy sources [9]. An attempt to address these conflicts is proposed in [13] and [14] by developing a coordinated voltage regulation method, without, however, ensuring the optimal network operation.
The proposed solutions to the bi-objective OVR problem can be classified into two main categories. In the first category, conventional optimization techniques are employed to solve the OVR problem [15]. Nevertheless, these methods present an increased computational complexity, limiting their applicability in real-time operation of distribution networks. To overcome this, a real-time optimal power flow algorithm based on quasi-Newton methods has been developed in [16], operating on a fast timescale but with suboptimal solutions. The authors in [17] propose a linearized approach based on Taylor series, where the coupling between network voltages and reactive power is taken into account, thus presenting an improved performance compared to the well-known dc optimal power flow (DC-OPF) method. However, inaccuracies and possible infeasible solutions cannot be sufficiently prevented. Sophisticated techniques have recently been developed to reduce the computational burden and improve the performance of the conventional optimization techniques [18], [19]. In particular, the conic relaxation method is employed in [18], while an enhancement of this method is introduced in [19], where the quadratic constraints are replaced with piecewise linear expressions. Although these techniques are valid, the above-mentioned assumptions may result in infeasible solutions during the network operation.
The second category includes the use of metaheuristic techniques to solve the bi-objective OVR problem. More specifically, a particle swarm optimization and a teaching-learning method are proposed in [12] and [20], respectively, whereas a differential evolution algorithm using the fuzzy set theory has been developed in [21]. Nevertheless, these methods present slow convergence rates and may lead to suboptimal solutions.
The common drawback of the developed methods lies in the use of generation and consumption forecasts to determine the operation settings of the participating network elements in the OVR process. As a result, possible security and reliability issues cannot be effectively tackled in case of forecast errors. Furthermore, a clear overview of the trade-off among the conflicting objectives must be available to the DSOs in order to decide the best compromise during the network operation. This is attained by obtaining the Pareto-front solutions of the OVR problem [20], which is, however, considered a computationally intensive process.
In this paper, a two-stage solution is proposed to tackle the bi-objective OVR problem in radial MV networks, using the available reactive power of DG units and the OLTC of the HV/MV transformer. The proposed methodology consists of two distinct stages. In the first stage, the sensitivity theory is employed to linearize the network operation. Afterward, the linearized OVR problem is solved multiple times to obtain the candidate OLTC operating plans (COOPs). In the second stage, an optimal reactive power allocation method is implemented for each COOP, following a rule-based approach. Scope of this method is to regulate the network voltages within permissible limits and minimize network losses. Finally, the Pareto-front is obtained which provides a clear overview of the trade-off between the conflicting objectives. By exploiting the Pareto-front, DSOs can select the most preferable OLTC schedule that will be applied in the real-time network operation along with the second stage of the proposed methodology.

II. THEORETICAL FRAMEWORK
The bi-objective OVR problem is considered as a mixedinteger nonlinear optimization problem that can be formulated as follows: where (1a) denotes the objective function, while (1b) and (1c) are the equality and inequality constraints, respectively. x is the vector of dependent variables containing the voltages of the distribution network, u c stands for the vector of the continuous control variables, and u d is the vector of the discrete control variables. In this paper, the reactive power of the DG units is treated as a continuous variable, whereas the tap position of the OLTC is a discrete control variable.
The objective function consists of the linear weighted combination of two conflicting objectives, i.e., the minimization of network energy losses (E loss ) and of the number of tap changes (T aps), calculated as follows: where w 1 and w 2 are two weight coefficients. Although costweighting-based coefficients can be used in (2) to scalarize the bi-objective optimization problem into a single-objective one, no unique solution exists due to existence of conflicting objectives [22]. To address this issue, the Pareto-front solutions of the bi-objective optimization problem are derived by varying the weight coefficients w 1 and w 2 . The analytical expressions of the optimization objectives are presented in the following: where tap t is the discrete control variable denoting the OLTC position during time instant t, whereas T is the time horizon set. The network energy losses are calculated using (4), which is the exact loss formula as presented in [23]. N is the set of network nodes, while P t i and Q t i stand for the net injected active and reactive power at the i-th node during time instant t, respectively. Furthermore, α t ij and β t ij are two coefficients calculated by: where V t i and θ t i denote the magnitude and angle of the voltage at i-th node during time instant t, while R ij is the real part of the ij-th element of the Z-matrix, which is the inverse of the network admittance matrix.
The equality constraints of (1b) describe the power flow equations and the OLTC operation. In particular, the power flow equations are expressed mathematically by [17]: where G ij and B ij are the real and imaginary part of the ij-th element of the network admittance matrix. Considering slack bus during time instant t, the voltage angle is equal to zero, while the voltage magnitude (V t sb ) is calculated according to (9) to model the OLTC operation.
Here, V t hv denotes the voltage magnitude of the HV network at time instant t, m stands for the voltage transformation ratio, and δ is the variation of the transformation ratio per tap position change. Furthermore, the inequality constraints of (10) are used to maintain the network voltages within permissible limits.
Here, V min and V max are the minimum and maximum permissible voltage limits, respectively. These limits are determined by the DSOs in compliance with national or international regulations, e.g., Standard EN 50160 [24]. It is worth noticing that (10) indirectly imposes limits on the voltage stability margin of the distribution network, since there exists a strong dependence of the voltage deviation on the network stability, as presented in [25]. Finally, (11) and (12) are introduced to model the reactive power capability of the DG units and the operation range of the OLTC, as follows: where N dg is the set of network nodes with DG units and D is the set of the available tap positions. Q t dg,i is a continuous control variable denoting the reactive power of the DG unit located at the i-th node during time instant t, while Q t min,i and Q t max,i are the corresponding permissible limits. These limits are timevariant and are calculated using (13) and (14), as presented in [26].
Here, S i denotes the rated apparent power of the DG unit located at i-th node, whereas P t dg,i is the corresponding injected power at time instant t.
Solving the optimization problem of (2)- (14) presents several difficulties which can be classified into the following categories based on whether metaheuristics or conventional optimization techniques are used: r Conventional optimization techniques: To obtain the Pareto-front solutions, multiple solutions of the optimization problem of (2)- (14) are required, assuming different combinations of the weight coefficients in (2). Nevertheless, this is a time-consuming process, preventing its implementation in real-field conditions. r Metaheuristic techniques: The above-mentioned problem can be effectively addressed using well-established metaheuristic techniques with an inherent ability of obtaining the Pareto-front solutions on a single simulation run, e.g., the non-dominated sorting genetic algorithm II (NSGA-II) or the group search optimizer with multiple producers (GSOMP) method [27]. However, these methods are characterized by slow convergence rates and may lead to suboptimal solutions.
r Common drawbacks: The complexity of the optimization problem is strongly related to the network characteristics and the time step of the performed analysis. The nonlinear behavior and the extended size of distribution networks are the most important factors that increase the problem complexity. Moreover, small time steps are required to examine the network performance under highly variable generated power, resulting in increased computational burden. Additionally, the use of forecast data in the optimization problem of (2)- (14) to determine the operation settings of the participating network elements cannot always guarantee the real-time safe and optimal network operation. For example, in case of forecast errors, the operation settings are miscalculated leading to suboptimal solutions during the real-time network operation, whereas possible voltage violations may also occur. To overcome the above-mentioned drawbacks, a new method is proposed that solves the optimization problem of (2)- (14) in two phases, namely the planning and the real-time operation phase. Scope of the planning phase is to determine the OLTC schedule that will be applied during the real-time operation phase. This is attained by solving the optimization problem of (2)- (14) in two successive stages. In the first stage, a linearized version of the OVR problem is solved to obtain the COOPs. For each COOP, a rule-based reactive power allocation method is applied to regulate network voltages and minimize energy losses, constituting the second stage of the proposed solution. Afterward, the Pareto-front is obtained and thus DSOs have a clear overview of the trade-off between the optimization objectives to decide the most preferable OLTC schedule. The time frame of the planning phase is determined by the DSO and can vary from few hours up to a day, depending on the availability of forecasts. Finally, in the real-time operation phase, the selected OLTC schedule is applied in conjunction with the rule-based reactive power allocation method used in the second stage of the planning phase.

III. FIRST STAGE
Scope of the first stage of the planning phase is to determine the COOPs. This is attained by employing the sensitivity theory to reduce the computational complexity of the bi-objective OVR problem by linearizing the network operation. Therefore, it is transformed into a mixed-integer linear optimization problem by replacing (3)-(9) of the conventional formulation with (15)- (20). More specifically, (15)- (18) are employed to linearize (3) as follows [18]: where T C t real is the real value of tap changes at time instant t calculated using (16). T C t abs is the absolute value of T C t real calculated according to (17) and (18) using two auxiliary nonnegative variables, namely T C t + and T C t − . The reason behind the use of these auxiliary variables lies in the following rationale: Every real number x real can be mathematically expressed as the difference of two non-negative numbers, i.e., x + − x − . Among the infinite combinations, the absolute value x abs corresponds to the pair of the non-negative numbers that minimizes the sum x + + x − . Thus, due to the use of (15) as a minimization objective, the pair of the non-negative numbers and thus the absolute value are correctly calculated. It is worth noticing that the use of these auxiliary variables is common practice in linear optimization theory [28].
Furthermore, the network energy energy losses are estimated using (19).
Eq. (19) is a linear approximation of the network energy losses with respect to the control variables, i.e., the OLTC position and the reactive power of the DG units. P t loss,pf denotes the network losses at time instant t multiplied by a linear function of tap t to model the dependence on the OLTC position. Further details regarding the derivation of this linear term are provided in the Appendix. Moreover, l t i stands for the sensitivity of network losses against the reactive power of the DG unit at i-th node during time instant t. Finally, the network voltages are calculated as follows: pf ,i is the voltage at the i-th node during time instant t multiplied by a linear term to model the dependence on the OLTC position, as analytically explained in the Appendix. Furthermore, s t ij is the ij-element of the sensitivity sub-matrix during time instant t, used to quantify the voltage variations at i-th node with respect to reactive power variations of the DG unit at j-th node.
The process to obtain the Pareto-front solutions is depicted in Fig. 1 by means of a flowchart. Initially, the generation and consumption forecasts as well as the network characteristics are determined and a power flow analysis is performed to calculate the matrices of the network voltages (V pf ) and losses (P loss,pf ). These values are calculated for each time instant t, assuming that the control variables, i.e., the reactive power of the DG units and the tap position, are equal to zero. The size of V pf is N × T , while the size of P loss,pf is 1 × T . Next, the sensitivity matrices (l, s) are calculated using the results obtained from the power flow analysis. More specifically, the size of sensitivity matrix l is N × T and can be analytically calculated following the analysis presented in [23]. Additionally, the size of sensitivity matrix s is N × N × T and is calculated for each time instant t by inverting the power flow Jacobian matrix [8]. Afterward, the linearized bi-objective OVR problem is solved multiple times for different combinations of the weight coefficients. Among a variety of well-established solvers that can handle mixed-integer linear optimization problems, e.g., CPLEX, Gurobi, MOSEK, etc. [7], CPLEX is employed in this paper to solve the linearized OVR problem. Thus, the Pareto-front solutions of the linearized OVR problem and the corresponding COOPs are derived resulting in reduced execution time compared to the initial nonlinear optimization problem.

IV. SECOND STAGE
As mentioned in Section II, scope of the planning phase is to determine the OLTC schedule that will be applied to the real-time operation phase. The decision regarding the most preferable OLTC schedule is made by the DSO, after evaluating the impact of each COOP on the network energy losses. Although the Pareto-front solutions of the linearized OVR problem can be used in the evaluation process, the corresponding losses are approximately calculated using (19). Thus, to accurately calculate the energy losses for each COOP, a coordinated reactive power control (CRPC) strategy is proposed in the second stage of the planning phase, assuming the full, nonlinear representation of the network. The proposed method follows a rule-based approach to optimally reallocate the reactive power among the DG units, targeting at the regulation of network voltages with minimum energy losses. Furthermore, the proposed CRPC strategy is applied in conjunction with the selected OLTC schedule during the real-time operation phase.
In the next subsections, the theoretical background and the actual implementation of the proposed method are analytically presented to address both overvoltages and undervoltages.

A. Overvoltage Mitigation
The first part of the CRPC method is related to the optimal dispatch of the consumed reactive power among the DG units to tackle overvoltages while also minimizing network losses. The theoretical background of the optimal dispatch strategy is based on the following concept: During time instant t, the DG unit connected to the i-th node will start to absorb reactive power only when the following conditions are met: 1) V t i exceeds V max and 2) V t i is the maximum network voltage. This concept has been mathematically proved in [29] but the therein developed implementation is applicable to dedicated MV feeders where only DG units are connected.
To overcome this limitation, a new control scheme is proposed. Its distinct feature involves the integration of time delays on the DG units to prioritize their response towards overvoltage mitigation. The introduction of the time delays is justified by the following analysis: Assuming a given time instant t, the reactive power absorption is initially undertaken by the DG unit with the maximum network voltage. Nevertheless, in the updated voltage profile, the maximum network voltage may be shifted to a different DG unit. In such a case, the reactive power absorption must be undertaken by the new DG unit, as imposed by the above-mentioned concept. Thus, time delays are introduced to determine the activation sequence of the DG units, whose location also coincides with the location of the maximum network voltages during the overvoltage mitigation process.
The allocation of the time delays among the DG units depends on the loading condition of the network and is carried out by a central controller, which is usually located at the DSO level. More specifically, the central controller constantly monitors the voltage at the point of common coupling (PCC) of each DG unit. Based on the network topology and the acquired voltages, the following time delay allocation process is implemented: Initially, the central controller categorizes the DG units based on whether they will participate in the overvoltage mitigation process or not. This is attained by comparing the PCC voltage of each DG unit with the corresponding voltages of the upstream DG units. A lower PCC voltage precludes the possibility that the maximum network voltage will occur at this DG unit. Thus, this DG unit will not participate in the overvoltage mitigation process, while the PCC voltages of the remaining DG units will form an increasing voltage profile. Afterward, the time delays of the participating DG units are determined by the central controller. Considering radial tree-like MV networks, the following twofold condition must be met: r In every tree path, the time delay allocation process follows the activation sequence of the participating DG units, as foreseen in the above concept. Thus, the smallest and the highest time delays are applied to the DG units with the maximum and minimum PCC voltages, respectively.
r In case of DG units with a common upstream node, the same time delays must be applied in order to be activated simultaneously and proportionally share their reactive power consumption during the overvoltage mitigation process. Finally, the time delays are forwarded to the participating DG units. This process is carried out at discrete time intervals, which are DSO-defined and may vary from seconds to minutes.
The control logic applied to the DG units is presented in Fig. 2 and consists of two operation modes. To avoid oscillations and repeated activation-deactivation cycles between these modes, a small voltage deadband (db) is introduced. The first operation mode is related to the proper increase of the reactive power consumption of the DG units to mitigate overvoltages, while also taking into account their reactive power capability limits as defined in (13). More specifically, initially, the DG unit acquires the PCC voltage (V t i ) and checks whether an overvoltage occurs. If so, after a predefined time delay (d t up,i ), the DG unit starts absorbing reactive power by employing a proportional-integral (PI) controller to zero the difference between V t i and the target voltage (V max − 0.5db). This process continues till the PCC voltage is effectively regulated or the equivalent network power factor (ENPF) seen downstream from the i-th node (spf t i ) is the minimum one (spf min,i ). The latter condition introduces an upper limit to the reactive power flowing through the upstream branch calculated using (21).
Here, subscript u denotes the upstream node of the i-th node. As shown in (21), spf t min,i depends only the line characteristics and can be obtained from the LinDistFlow equations of [30], assuming zero voltage drop between two consecutive nodes. The zero voltage drop indicates that there exist upstream nodes with voltages equal to or greater than of the i-th node. Thus, this limit is introduced to avoid the excessive and unnecessary reactive power consumption. Finally, the DG unit switches to a constant ENPF operation, leading to a constant voltage drop on the upstream branch.
The second operation mode involves the reverse process of decreasing the reactive power consumption. This mode is activated when V t i falls below the voltage threshold (V max − db). After a predefined time delay (d t down,i ), the reactive power of the DG unit is reduced using a PI controller until zero is reached or the voltage regulation is accomplished. In the latter case, the constant spf t i operation is activated. It is worth mentioning that the time delays related to this operation mode are determined following a similar rational as presented above. In particular, in every tree path, the smallest and the highest time delays are applied to the participating DG units with the minimum and maximum PCC voltages, respectively.

B. Undervoltage Regulation
During the overvoltage mitigation process, the reactive power absorbed by the DG units is provided by the HV grid through the HV/MV transformer. As a result, the voltage drop on the transformer is increased, which, in turn, reduces the network voltages. Thus, there is a possibility of undervoltages, especially in passive feeders where only loads are connected. Scope of the second part of the CRPC method is to effectively address these undervoltages. This is achieved by exploiting the reactive power of the DG units that do not participate in the overvoltage mitigation process in order to provide locally reactive power support.
The proposed reactive power support strategy follows a similar approach to the overvoltage mitigation control scheme. In particular, every DG unit acquires the minimum network voltage (V t end ), which usually occurs at the end nodes of the feeders and especially of passive feeders. This value is forwarded to the DG units by the central controller which constantly monitors the end node voltages of these feeders. In case of undervoltages, the DG unit connected to the i-th node starts to produce reactive power following a predefined time delay (c t up,i ). The time delays are introduced to prioritize the response of the DG units based on the impact that may have on the corresponding PCC voltage. This impact can be quantified using the voltage  sensitivity sub-matrix s. The smallest time delay is applied to the DG unit with the minimum voltage sensitivity coefficient, i.e., located close to the HV/MV transformer, indicating that it can provide reactive power support with the minimum impact on the network voltages. This process continues till the maximum reactive power limit of the DG unit, as expressed by (14), is reached or the undervoltage is effectively addressed. In the latter case, the DG unit switches to constant power factor operation. A similar rational is implemented on the reverse process of properly reducing the reactive power production, where the smallest and the highest time delays (c t down,i ) are applied to the DG units with the maximum and minimum voltage sensitivity coefficients, respectively.

V. NUMERICAL RESULTS
The performance of the proposed methodology is assessed on the 20 kV radial MV distribution network of Fig. 3. Further details regarding the technical characteristics of the network can be found in [31]. Although different DG technologies can be considered, only PV units are assumed in the simulated scenarios without affecting the performance and the general applicability of the proposed methodology. This is justified by the fact that the different technical features of these technologies mainly affect the generation profile and the reactive power capability limits, which, however, have been already taken into account in the proposed methodology presented in Sections III and IV. PV units and loads are numbered according to their connection node, whereas their rated power is presented in Table I. The  TABLE II  VOLTAGE SENSITIVITY COEFFICIENTS OF NODE 25 WITH RESPECT TO  REACTIVE POWER VARIATIONS OF PV NODES negative sign indicates active power consumption. The power factor of all loads is considered equal to 0.95 lagging, while the nominal power factor of the grid-interfaced inverters of the PV units is considered equal to 0.85. The permissible voltage limits are equal to ± 5 % of the nominal value. Considering the CRPC method, PV3 and PV4 have been selected to address undervoltages, since they present the smallest voltage sensitivity coefficients against reactive power variations, as presented in Table II. These coefficients have been calculated assuming the rated generation and consumption values of Table I. The remaining PV units participate in the overvoltage mitigation process.
In the next subsections, dynamic simulations are performed to validate the CRPC method during the real-time operation phase. Furthermore, time-series simulations are employed to evaluate the overall long-term performance of the proposed two-stage method of the planning phase.

A. Dynamic Simulations
The examined network of Fig. 3 is modeled with the PSIM software that is widely used for time-domain simulations. The simulation period is equal to 3.2 s. The voltage of the HV busbar is kept constant to 1.05 p.u. and the tap position of the HV/MV transformer is equal to zero. The loads connected between nodes 27 and 36 operate at 50 % of the nominal capacity, whereas the remaining ones operate under nominal conditions. Furthermore, the active power output of some indicative PV units is depicted in Fig. 4a, while the remaining ones operate at their rated power. The voltage magnitude at the end-nodes of the passive feeder is forwarded to PV3 and PV4 that participate in the undervoltage mitigation process, assuming a communication delay of 0.03 s. Finally, the voltage deadband db is set to 0.002 p.u., while the time delay between two sequential activations of the overvoltage mitigation process is equal to 0.2 s. The same time delay is also applied between two sequential activations of the undervoltage regulation process.
The reactive power of the PV units is presented in Fig. 4b, whereas the voltage magnitude of several network nodes is depicted in Fig. 5.
The first part of the proposed CRPC method, i.e., the overvoltage mitigation, is activated at 0.1 s. The time delay allocation has been implemented according to the network voltages prior to 0.1 s, following the approach presented in Section IV-A. For example, PV25 will instantly react on overvoltages, PV23 and PV24 after a time delay of 0.2 s, PV15 and PV20 after 0.4 s, and so on. Considering the second operation mode of properly reducing the reactive power absorption, PV14 will instantly At 0.1 s, PV25 starts reducing the PCC voltage by absorbing reactive power till the minimum ENPF is reached at 0.13 s, which indirectly denotes the equality of the PCC voltages between PV24 and PV25. After a time delay of 0.2 s, at 0.3 s, the reactive power absorption process is undertaken by PV23 and PV24 and the corresponding PCC voltages, as well as the voltages at the load nodes of the active feeder, are regulated at 0.48 s, as shown in Figs. 5a and 5b. Nevertheless, although the absorbed reactive power reduces the network voltages bellow the maximum permissible limit, an undervoltage violation occurs at node 43 of the passive feeder according to Fig. 5c.
To address this problem, the second part of the proposed CRPC method is activated at 1 s. According to the analysis presented in Section IV-B, PV3 will instantly react on undervoltages, while the undervoltage regulation control scheme of PV4 will be activated after a time delay of 0.2 s. As a result, at 1 s, PV3 starts to produce reactive power till the minimum power factor is reached at 1.18 s and the undervoltage at node 43 is successfully addressed. Nevertheless, at 1.15 s, the reactive power production of PV3 leads to overvoltages at the PCC of PV units, as depicted in Fig. 5a. Thus, the overvoltage mitigation process of PV23 and PV24 is reactivated after a time delay of 0.2 s, i.e., at 1.35 s. Once again, the reactive power absorption of PV23 and PV24 results in an undervoltage violation. Therefore, after a time delay of 0.2 s, PV4 starts to produce reactive power and the network voltages are finally regulated at 1.8 s, as shown in Fig. 5.
The reverse process is activated at 2 s by reducing the active power injection of PV14 and PV20, as depicted in Fig. 4a, which, in turn, reduces the PCC voltages of the PV units below the voltage threshold of V max − db. After a time delay of 0.4 s, PV23 and PV24 start to reduce the absorbed reactive power. This action increases the voltage at node 43 of the passive feeder above the voltage threshold of V min + db at 2.5 s. PV4 reacts instantly by reducing the reactive power production. This process continues till the network voltages are finally regulated, as shown in Fig. 5.

B. Long-Term Evaluation
In this section, the long-term performance of the proposed two-stage method of the planning phase is assessed by performing time-series simulations. The simulation period is one day with one minute time interval. Normalized generation and consumption profiles, similar to those presented in [31], are arbitrary distributed to PV units and loads, respectively. Considering PV units, the following process is adopted to derive the generation profiles: Initially, active power measurements from existing PV installations are acquired for a cloudy day with one minute resolution. A cloudy day has been selected instead of a sunny day to evaluate the performance of the proposed method under highly variable generation. Afterward, the generation profiles are normalized by dividing the acquired measurements by the rated power of the corresponding PV installation.
The coefficients of the linear terms of (19) and (20) that model the dependence of the OLTC on the network voltages and losses are presented in Table III. It can be observed that the linear  models can effectively replace the non-linear terms, as verified by either the small root mean square error (RMSE) or the high coefficient of determination (R 2 ).
Considering the first stage of the proposed method, the linearized model as expressed by (10)-(20) is solved 100 times, assuming different combinations of the weight coefficients. More specifically, w 1 varies from 0 to 1 with a step of 0.01, whereas w 2 is the complement of w 1 . The linearized OVR problem is solved using the CPLEX solver in GAMS. The average solution time is equal to 12 s, whereas the solution time of the nonlinear problem, which is solved using the BONMIN solver in GAMS, exceeds four orders of magnitude. This can be regarded as an important shortcoming towards the adoption of a nonlinear model for the operation of the distribution network under real-time conditions.
The Pareto-front of the linearized OVR problem and four daily COOPs are presented in Figs. 6 and 7, respectively. Concerning Fig. 6, the percentage increase of the energy losses is calculated with respect to the Pareto-front solution of 26 daily tap changes.
The results of the second stage are summarized in Table IV, where the proposed method is compared to an optimized solution. More specifically, for each COOP, an optimization problem is solved aiming to minimize network losses by optimally  reallocating the reactive power among the PV units. The optimization problem is solved in GAMS using the BONMIN solver. It can be observed that the proposed rule-based reactive power allocation method presents a similar performance against the use of optimization techniques.
Finally, a Monte Carlo (MC) analysis is performed to investigate the optimality of the calculated COOPs. In particular, for each Pareto-front solution, 100 unique OLTC operating plans are randomly created. For each OLTC operating plan, an optimization problem is solved aiming to minimize network losses by optimally distributing the reactive power among the PV units. The results of this analysis are presented in Fig. 8, where the percentage increase of the energy losses is calculated with respect to the case the optimization problem of (2)-(14) is solved, assuming w 1 and w 2 are equal to 0 and 1, respectively. It can be observed that the proposed two-stage solution presents a similar performance compared to the optimized solution and the calculated COOPs are near-optimum ones.

VI. CONCLUSIONS
In this paper, a two-stage solution to the bi-objective OVR problem is proposed. The developed method is characterized by reduced computational complexity, compared to the use of conventional optimization techniques, thus allowing its adoption for the real-time network operation.
The proposed method is evaluated by performing timedomain and time-series simulations. The former clearly shows the quick response and the ability of the proposed methodology to regulate network voltages in real-time, while also minimizing network losses. The latter highlights its near-optimal performance, compared to the use of non-linear optimization methods. Therefore, it can be considered as a valuable real-time control scheme for DSOs towards the enhanced operation of distribution networks.

APPENDIX
Following the analysis in [32] where the simplified DistFlow equations are presented and using (9), the network losses and voltages can be approximately calculated as follows: where P t ij and Q t ij are the active and reactive power flowing between the i-th and j-th nodes, respectively, while V t nt,i denotes the voltage at i-th node assuming zero tap position. In both expressions, the nonlinear terms that model the OLTC operation increase the computational complexity of the optimization problem. Nevertheless, due to the limited range of the available tap positions, these nonlinear terms approximate a linear behavior. Thus, to overcome this drawback, the nonlinear terms are replaced with linear ones as shown in (19) and (20), where coefficients a, b, c, and d can be calculated by simply applying curve fitting techniques. These coefficients depend on the characteristics of the OLTC, namely the available tap positions and the parameter δ. Therefore, for a given OLTC configuration, they are considered as constant values in the linearized OVR problem.