Multi-Time-Scale Optimal Scheduling in Active Distribution Network with Voltage Stability Constraints

: The uncertainty associated with loads and renewable-energy sources affects active distribution networks in terms of the operation and voltage stability on different time scales. To address this problem, a multi-time-scale voltage stability constrained optimal scheduling framework is proposed, which includes a day-ahead model with a coarse-grained time resolution and an intra-day model with a ﬁne-grained time resolution. The day-ahead economic-scheduling model maps out a scheme to operate different types of devices with the aim of minimizing the network losses. Following the scheme, the intra-day corrective-adjustment model based on model predictive control is proposed to regulate the ﬂexible devices, such as the energy storage systems and the photovoltaic converters. In particular, the proposed optimal scheduling framework embeds a voltage stability constraint which is constructed by using a novel index, deﬁned based on the Distﬂow model Jacobian. As the index at each bus is a linear function of the locally measurable power ﬂow variables, the proposed constraint does not introduce additional computational burdens. Simulation results demonstrate the necessity and effectiveness of the proposed multi-time-scale voltage stability constrained optimal scheduling model. The results also show that the variation trend of the proposed index is consistent with that of the commonly used voltage stability index.


Introduction
Active distribution networks (ADNs) are incorporating increasing demand-responsive loads, such as plug-in electric vehicles, smart buildings, and flexible domestic appliances. As a result, most ADNs are sometimes operating under stress and close to their loadability. This trend, together with the rapidly growing penetration of renewable-energy sources, introduces challenges to ADN operation and control. The ability to assess and maintain voltage stability in this operational context is crucial. Voltage instability, and ultimately voltage collapse, are the main factors in some well-known blackouts, e.g., the case on 12 July 2004, in Greece [1]. The high R/X ratio (0.5-7 [2]) of ADNs makes both the active power and the reactive power have a significant influence on the voltage amplitude and stability. To this end, this paper seeks to guarantee voltage stability in ADNs through a novel multi-time-scale voltage stability constrained optimal active and reactive power scheduling (hereinafter referred to as VSCOS for short).
In the early literature, the voltage stability of high-voltage transmission systems used to be of a greater concern than that of the distribution networks. Several competent methods, in terms of voltage stability analysis and assessment, have been developed. Typically, the continuation power flow proposed in [3] is used to map the P-V curves that provide a graphical tool for voltage stability analysis. The loadability point on the P-V curve corresponds to the nose point or the voltage collapse point. This contributes to an index, namely the load margin [4] (or voltage stability margin [5] or loadability margin [6]), defined as the distance between the current operating point and the loadability point along the power factor direction. An alternative option is based on the eigenvalue Despite the efforts on the voltage stability assessment and enhancement of the ADNs, the voltage stability indices and constraints in these studies either lead to nonlinear models or introduce a large number of variables for identifying the loadability point, resulting in great computational burdens. The loadability coefficient is calculated by a second-order cone program (SOCP) in [28], but the model should be used as a stability subproblem in distribution network optimization and thus a complicated bi-level problem needs to be solved. Based on the bus injection power flow model, reference [29] has proposed an index that is linear with respect to load increase. However, due to the nonlinear and nonconvex power flow model, the index is more suitable for stability assessment and monitoring rather than optimal-operation problems. The branch flow equation proposed in [30] (also referred to as the DistFlow model), thanks to its ability to be exactly relaxed to a convex SOCP problem, is widely used in the planning and optimal-operation problems of distribution networks. Under this power flow framework, reference [31] proposes a nonlinear approximate voltage stability index and verifies the performance and numerical robustness of the index in real-time stability monitoring and assessment. To the best of our knowledge, the linear indices and constraints applicable to the DistFlow-based optimal scheduling and planning problem of radial distribution networks have not been reported.
Furthermore, the errors in forecasting power injection and consumption are enlarged by demand-responsive loads and intermittent renewable-energy sources. Hence, the direct use of a day-ahead scheme is at risk of power imbalance, voltage instability, and voltage violation. A fine-grained time resolution can improve the situation, but this open-loop optimal scheduling considers none of the feedback of the state realization on operation and shows a modest effect. Although some day-ahead optimal scheduling approaches (e.g., [32,33]) consider in advance the uncertainty associated with the load demand and the renewable-energy generation of the next day, the impact of the uncertainty can only be mitigated, not eliminated. So, the intra-day adjustments are inevitable. Herein lies the motivation to perform an intra-day corrective adjustment based on model predictive control [34]. The intra-day model regulates the flexible devices according to the real-time operating state and thus requires the fast solvability of the optimal adjustment model. However, it is opposite to the computational burdens caused by the existing voltage stability constraints.
To address these problems, this paper proposes the VSCOS model. The framework of the model contains day-ahead (slow time scale) and intra-day (fast time scale) operation and coordinates the discrete controllers (i.e., the OLTCs and the capacitor banks) and the continuous controllers (i.e., the PV converters and the ESSs). The main contributions of this paper are concluded as follows.

1.
A voltage stability index is proposed based on the Jacobian of Distflow branch flow model with no assumption. The index is a linear function of the locally measurable power flow variables. The simulation shows that the index has the same change trend as the common index load margin over the time horizon. A linear voltage stability constraint is proposed by using the index and can be embedded in the optimization problem of the radial distribution networks without increasing the computational burdens.

2.
The VSCOS model is formulated as a multi-time-scale framework. The day-ahead model coordinates different types of devices to guarantee the voltage stability and minimize the network losses. To cope with the uncertainty of the loads and the renewable energy, the intra-day model reschedules the PV converters and the ESSs on a shorter fine-grained time scale.
The remainder of this paper is organized as follows. Section 2 introduces the distribution system model and proposes the voltage stability constraint. Section 3 provides the mathematical formulation of the multi-time-scale VSCOS model. In Section 4, the numerical results and analysis are provided. The paper is summarized in Section 5.

Radial Distribution Network Model
Let G = (N, E) be a directed tree representing a radial distribution network, where each node in N = {0, 1, . . . , n} represents a bus and each link in E represents a branch. Note that |E| = n. A directed edge in E is denoted by (i, j), where i is the parent of j. For each i ∈ N, let v i , p i , and q i denote the square of voltage amplitude, the active power and reactive power of the net load demand, respectively. The vectors v, p N , and q N are obtained by stacking v i , p i , and q i , respectively. For each (i, j) ∈ E, let r ij , x ij , l ij , p ij , and q ij denote the branch resistance, reactance, the square of branch current, and the sending-end active power and reactive power, respectively. The vectors r, x, l, p E , and q E are obtained by stacking r ij , x ij , l ij , p ij , and q ij , respectively. Node 0 in N corresponds to the substation bus with its voltage amplitude fixed (e.g., equal to 1 p.u. or 1.05 p.u.). Here, we assume that the notation [b] for an arbitrary vector b ∈ R n (R denotes the set of real numbers) represents the n × n matrix with elements of b on the diagonal and zeros everywhere else. The power flow in a radial ADN is formulated based on the DistFlow branch equations [30,31] in vector form, as follows: where Π ∈ R (n+1)×n and Γ ∈ R (n+1)×n are {0,1} coefficient matrices that denote the incidence of the tree G. We further write Equations (1)-(4) in the compact form: vector of parameters ξ = [p 1 , p 2 , · · · , p n , q 1 , q 2 , · · · , q n , v 0 ] T ∈ R 2n+1 .

Voltage Stability Constraint
In order to guarantee voltage stability, it is necessary to ensure a secure distance between the current operating point and the loadability point. The loadability point is obtained by the following optimization problem which is in line with the approaches proposed in [6,23,26,31]: max where f (ξ) is a linear function of ξ.
The optimal operation of ADNs generally takes benefit and cost as the objective function. To consider the voltage stability, the loadability point is required and usually obtained by solving model (6), resulting in a complex bi-level problem. An intuitive way to reformulate the model is to adopt the KKT conditions [35], and it results that the loadability point corresponds to the power flow Jacobi matrix J h (y), as given in (7), being singular, i.e., |J h (y)| = 0.
We will derive a voltage stability constraint based on the Jacobi matrix J h (y) in the remainder of this subsection. First, a new matrix J h (y) is obtained by removing the power balance equations associated with node 0 and the net power variables p 0 and q 0 from the Jacobi matrix J h (y). The matrix J h (y) is depicted in (8), and the relationship between its determinant J h (y) and the Jacobian |J h (y)| is shown in (9).
where A ∈ R n×n is a coefficient matrix for branch power flow, A = Γ 2 − Π 2 ; Γ 2 and Π 2 are obtained by removing the first row from Γ and Π, respectively. The determinant of J h (y) is calculated in (10): where the first two equalities follow from the direct application of the Schur complement; the third equality holds due to |−A| = 1 and A T = (−1) n ; the second Schur complement H ∈ R n×n is expressed explicitly in (11).
Therefore, we obtain (12) without making any assumptions. That is, the Jacobian |J h (y)| that reflects the voltage stability can be replaced exactly by the determinant |H|.
|J h (y)| = |H| (12) For i, j ∈ {1, 2, . . . , n}, let a ij represent the i-row and j-column element of H. Each a ij of H can be expressed explicitly as follows: For i = j, For i < j, For i > j, where the notation π(i) with respect to bus i represents the upstream bus of a branch (π(i), i) ∈ E; r 0π(i) and x 0π(i) are the sum of the resistances and the inductances of the branches connecting bus 0 to bus π(i), respectively. It can be observed that the off-diagonal elements given in (14) and (15) are non-positive.
If there is no power demand in the ADN, the power flows on branches are equal to zeros, and hence the voltage vector v = 1 (assume v 0 = 1 p.u.). In this case, the diagonal elements a ii = 1, and the off-diagonal elements a ij = 0 (i = j). Thus, the matrix H is strictly diagonal dominant with its determinant |H| = 1. With the power demand increasing, the value of each element in H decreases. At the loadability point, the determinant Energies 2021, 14, 7107 6 of 20 |H| = 0. Considering that each element in H is a continuous function of variables y, the first violation of strictly diagonal dominance in a row (or a column) of H must occur ahead of voltage instability. This feature provides the inspiration for establishing the voltage stability constraint by using diagonal dominance. It can also be found that all the elements of an arbitrary row i in the matrix H are linear with the variables corresponding to bus i, and the coefficients of the variables are related to the topology and branch parameters of the ADN. Thus, the diagonal dominance of rows is studied here and the voltage stability index H i is defined in (16).
As the value of the index H i relates only to the network parameters and local power flow variables, it can be regarded as a value reflecting the influence of bus i on the voltage stability level. Thus, the weak buses in the ADNs can be identified by H i . According to the above analysis, if H i ≥ 0 holds for every bus i, |H| > 0 and the voltage is stable. If ∃i ∈ N\{0}, H i = 0, the current operating point is still at a distance from the loadability point, indicating that H i ≥ 0 for every bus i is a sufficient condition for voltage stability. From the perspective of the ADN operation, ensuring row diagonal dominance is equivalent to leaving a certain stability margin for the ADN. Although the voltage stability index is linear, it does not provide a linear indication of the voltage stability level. The specific one-to-one match between H i and the voltage stability level is not studied in this paper because the purpose is to guarantee stability rather than to intuitively assess the stability level. It will be a part of future works.
As it is desired that the operating point keeps a secure distance from the loadability point, the voltage stability constraint is proposed to ensure that the index H i is greater than a non-negative threshold m, as in (17).
Recalling the power flow model and the entire process of derivation in this section, the advantages of constraint (17), together with the major limitations that could be addressed in future work, are summarized as follows.

•
First of all, the main advantage is its good linearity with respect to the decision variables. According to the definition of the index H i , the proposed voltage stability constraint of an arbitrary bus i ∈ N\{0} is a linear constraint and does not introduce any new variables. Therefore, it can be conveniently and directly integrated into the optimal scheduling and planning problems and cause no additional computational burdens.

•
Second, because all nodes in N\{0} are treated as PQ buses, the effect of the dynamic characteristics of the power injection and consumption on voltage stability cannot be considered in this model. Hence, the proposed constraint is essentially a steadystate voltage stability constraint. However, the aim of this paper is to ensure a safe voltage stability margin by coordinating various devices in the ADNs, so the proposed constraint is sufficiently effective. Nevertheless, it is important to point out that this constraint is not recommended when studying the short-term stability of the distribution networks hosting a large number of loads with dynamic voltage response characteristics. • Furthermore, the voltage stability constraint is applicable to balanced radial distribution networks. Although it is common that distribution networks are designed as loop topology and operated with radial network arrangements, extending the proposed constraint to all possible arrangements (i.e., radial, mesh, and loop) is a potential future improvement. As an increasing number of plug-in loads (especially electric vehicles) and behind-the-meter rooftop PVs are connected to the distribution networks through a single phase, it is promising to study the extension of the voltage stability constraint to unbalanced three-phase networks.

Multi-Time-Scale Optimal Scheduling Model
In this section, a multi-time-scale VSCOS model is proposed. The day-ahead model takes ∆t = 15 min as the time resolution. Due to the uncertainty associated with the distributed renewable generations and load demand, an intra-day operation model is developed to correct the flexible devices in the ADN every 5 min. The controllable device model and the objective function, as well as the constraints included in the VSCOS model, are provided, along with the notation used for the variables.

Controllable Device Model
The models of the controllable devices are given as follows.

On-Load Tap Changer
The action of the OLTC can be described by where E OLTC ⊆ E represents the set of branches with the OLTCs; the italic T represents the time horizon (otherwise the notation of transposition); k ij,t is the tap ratio of the OLTC in branch (i, j) at time span t; k 0 is the benchmark tap ratio; ∆k ij is the tap ratio step of the OLTC in branch (i, j) at time span t; K ij,t is the tap position of the OLTC in branch (i, j) at time span t; and K max ij and K min ij are the upper and lower tap position limits of the OLTC in branch (i, j). Moreover, the constraint (3) for the branches with the OLTCs should be modified as follows:

Energy Storage System
ESSs can charge from or discharge to the ADN, while satisfying the following operational constraints. Specifically, the temporal relationship between the energy and the power of the ESS is given in (21). The maximum and minimum allowable energy stored in the ESS is set by (22). The technical limits for the charging and discharging power are given in (23) and (24). Constraint (25) avoids simultaneous charging and discharging. Constraint (26) limits the state of charge at the end of the time horizon within a tolerance band of the initial condition. For ∀i ∈ N ESS , ∀t, e i,0 is the initial state of the ESS at bus i; and µ is a tolerance band for energy stored in the ESS at the end of the time horizon.

Capacitor Bank
The operational constraints for the capacitor banks (CBs) are shown in (27) and (28).
where N CB ⊆ N represents the set of buses equipped with the CBs; q C i,t is the reactive power of the CB at bus i; ∆q C i is the reactive power step of the CB at bus i; c i,t is the number of capacitors switched on at bus i; and c max i is the upper limit of the CB.

PV Converter
The converters of the PVs can provide reactive power subject to their capacities, as shown in (29). For ∀i ∈ N PV , ∀t, where N PV ⊆ N represents the set of buses with the PVs; s PV i is the capacity of the PV converter at bus i; p PV i,t is the PV active power output at bus i at time span t; and q PV i,t is the reactive power provided by the PV converter at bus i at time span t.

Day-Ahead Optimal Scheduling
Day-ahead optimal scheduling for the ADN aims to minimize the total branch losses. The devices, including the OLTCs, ESSs, CBs, and PVs, are scheduled using the predicted values of the load demand and PV outputs. The mathematical formulation of the day-ahead model is presented as follows.
subject to (18), (19), (21)- (29), and where n T denotes the number of time spans over the time horizon; the notation b 2 B for a vector b ∈ R n and a diagonal matrix B ∈ R n×n represents b T Bb; Θ, ∆, and Ω are n × n diagonal weight matrices; θ t , δ t , and ω t ∈ R n are vectors of slack variables with non-negative entries introduced to relax the voltage amplitude constraint and voltage stability constraint; p L i,t and q L i,t are the active and reactive load power at bus i at time span t; h O (y t , ξ t ) = 0 represents the power flow model of the ADN with the OLTCs; H t ∈ R n is the vector of the voltage stability indices at time span t, i.e., H t = [H 1,t , H 2,t , · · · , H n,t ] T ; m ∈ R n is the threshold vector for H t ; v max and v min are vectors of the upper and lower limits on the squared voltage amplitudes; and l max is a vector of the upper limits on the branch power flows. The objective function (30) minimizes the total branch losses over the time horizon as well as the penalties on violating the voltage amplitude limits and stability constraint. Equations (31) and (32) give the net active load power and the net reactive load power at bus i at time span t. Equation (33) represents the linearized Distflow model where the Equation (3) for the branches with the OLTCs is replaced by a linearized version of Equation (20). It is noteworthy that the non-linear Equations (4) and (20) make the model difficult to be solved by commercial solvers. To this end, we employ an exact linearization method [36] for Equation (20) and an exact second-order cone relaxation method [37] for Equation (4). Equation (33) is obtained by applying these methods to (5) and consists of linear constraints and second-order cone constraints. The constraints (34)-(36) set acceptable bounds on the voltage stability indices, the voltage amplitudes, and the current flows, respectively.

Intra-Day Corrective Adjustment
The day-ahead scheme cannot be used directly in practice due to the uncertainty of the loads and the PVs. We propose an intra-day optimal operation model in this subsection based on model predictive control. In the intra-day operation, the day-ahead scheme for discrete control devices (i.e., the OLTCs and CBs) is executed, and the scheme for the other devices is used as the reference trajectory. The intra-day optimization produces a corrective scheme that describes the adjustments of the ESS charging and discharging power and the PV converter reactive power at each time span over time interval t, t + n p ∆t . Only the optimal scheduling sequence for the first time span is executed. Then, the states are measured (or estimated) based on the realized load demand and the PV generations at time t + ∆t. This process is repeated every ∆t = 5 min for the intra-day operation to respond to the uncertainty. The mathematical formulation, whose objective is to track the day-ahead scheme as well as avoid penalties on violating the voltage-related constraints, is given as follows.
subject to (21)-(24), (29), (31)- (36), and where the state variables and Λ ∈ R (N ESS +1)×(N ESS +1) is a diagonal weight matrix; R t represents the reference trajectory corresponding to Y t at time span t; 0 denotes vector or matrix with all entries of 0; and I N ESS ∈ R N ESS ×N ESS and I N PV ∈ R N PV ×N PV are identity matrices. Note that the values of the matrices F X and F U are obtained according to (21) and the definition of the state variables X t and the control variables U t .

Case Study
In this section, we present numerical experiments on a modified 33-bus test system to demonstrate the effectiveness of the proposed multi-time-scale VSCOS model. All computations are carried out using MATLAB software with YALMIP toolbox and solved by GUROBI solver [38] on a computer with an Intel 2.90GHz CPU and 16GB RAM.

Test System
The modified 33-bus test system is shown in Figure 1. The network parameters and topology of the test system are set in line with the IEEE 33-bus system. The predicted daily profiles of the loads and the PV active power outputs are given in Figures 2 and 3.
Five ESSs are installed at buses 6, 15, 21, 24, and 30, with a 200 kW/1000 kWh capacity for each ESS. The maximum and minimum state of charge of the ESSs are 90% and 10%, respectively. The charging and discharging efficiency of the ESSs is 0.95. The initial state of charge is 40%. Three PVs are connected to buses 5, 19, and 24. The capacities of the PVs and their power converters are 600 kW and 600 kVA. This system has 3 OLTCs, and their locations are shown in Figure 1. The voltage control range of the OLTCs is between 0.95 p.u. and 1.05 p.u. (0.005 p.u. per step), and the initial voltage is 1.0 p.u. This system is equipped with a total of 4 CBs. The rated capacity of the CB at bus 26 is 400 kVar, comprising 8 capacitors (8 × 50 kVar = 400 kVar). The range of the other CBs (at buses 3, 9, and 16) are considered between 0 and 500 kVar with a discrete step size of 50 kVar, i.e., 11 switching positions. Additionally, the total number of control actions of the discrete devices is limited to no more than 8 over the time horizon.
The lower and upper bound of the voltage amplitude are 0.95 p.u. and 1.05 p.u., respectively. Without loss of generality, we assume that the voltage at bus 0 is fixed to be unity, i.e., v 0 = 1 p.u. The threshold for the stability constraint is 0.9. For the day-ahead model, the sample time span is selected as 15 min, and the time horizon n T is 24 h. The intra-day model is performed every 5 min, and the prediction horizon n p is 20 min.  locations are shown in Figure 1. The voltage control range of the OLTCs is between 0.95 p.u. and 1.05 p.u. (0.005 p.u. per step), and the initial voltage is 1.0 p.u. This system is equipped with a total of 4 CBs. The rated capacity of the CB at bus 26 is 400 kVar, comprising 8 capacitors (8 × 50 kVar = 400 kVar). The range of the other CBs (at buses 3, 9, and 16) are considered between 0 and 500 kVar with a discrete step size of 50 kVar, i.e., 11 switching positions. Additionally, the total number of control actions of the discrete devices is limited to no more than 8 over the time horizon.  and their power converters are 600 kW and 600 kVA. This system has 3 OLTCs, and their locations are shown in Figure 1. The voltage control range of the OLTCs is between 0.95 p.u. and 1.05 p.u. (0.005 p.u. per step), and the initial voltage is 1.0 p.u. This system is equipped with a total of 4 CBs. The rated capacity of the CB at bus 26 is 400 kVar, comprising 8 capacitors (8 × 50 kVar = 400 kVar). The range of the other CBs (at buses 3, 9, and 16) are considered between 0 and 500 kVar with a discrete step size of 50 kVar, i.e., 11 switching positions. Additionally, the total number of control actions of the discrete devices is limited to no more than 8 over the time horizon.

Simulation Results and Analysis
The day-ahead scheme is obtained by employing the predicted power profiles. The PV active power output could have high ramping characteristics where the ramping may go up to 90% of its rated capacity due to uncertain events such as the passing of fast-moving clouds [14]. However, such intermittent variations are difficult to predict and almost impossible to observe at a coarse time resolution. The intra-day model is optimized under two weather conditions, i.e., clear sky and cloudy sky. The simulation results are shown as follows. In some of the later figures, the results are marked with two types of lines, where the solid lines represent the schemes obtained by the VSCOS model, and the dotted lines of the same color correspond to the results of the model without considering the voltage stability. Figure 4 shows the tap position of the OLTCs. It can be seen that the turn ratio of OLTC1 keeps 0.955:1 over the time horizon. The tap position of OLTC2 varies between −2 and −3, while that of OLTC3 stays at −4 before 17:00 and decreases to −5 subsequently.

Day-ahead Scheme Obtained by the VSCOS Model
The reactive power of the CBs is shown in Figure 5. The CBs connected to Feeder1 provide a great deal of reactive power, and the feeder is overcompensated in order to achieve the acceptable voltage amplitudes. The CB at the end of the feeder (i.e., CB16) outputs more reactive power because the voltage amplitudes here are lower. The total reactive power injected by the CBs shows a significant rise from 9:15 to 11:00. This is related to the increase in reactive load power during this period.  The reactive power of the CBs is shown in Figure 5. The CBs connected to Feeder1 provide a great deal of reactive power, and the feeder is overcompensated in order to achieve the acceptable voltage amplitudes. The CB at the end of the feeder (i.e., CB16) outputs more reactive power because the voltage amplitudes here are lower. The total reactive power injected by the CBs shows a significant rise from 9:15 to 11:00. This is related to the increase in reactive load power during this period.   Figure 6 shows the reactive power of the PV converters. As can be seen, during the peak solar irradiation hours around noon (10:00-14:00), the PV converters can generate a small amount of reactive power limited by their capacities. At the other times, more reactive power is available. The PV reactive power outputs are adjusted frequently as the converters can be controlled flexibly and continuously. Due to the lack of reactive sources in Feeder2 and Feeder3, PV19 and PV24 inject a large amount of leading reactive power into the feeders. By contrast, PV5 regulates the power flow by generating both leading and lagging reactive power.   It is observed that the ESSs installed in Feeder1 and Feeder4 all charge before 8:00 and discharge from 15:00 to 21:15. As these two feeders have few active power sources, the ESSs tend to store a lot of energy to cope with the low voltage problem and guarantee the voltage stability via discharging during the peak load period. Feeder2 and Feeder3 are not as stressed as the others. ESS21 and ESS24 discharge in the early morning and then charge when the PV outputs rise between 7:15 and 14:45. Thus, the role of ESS21 and ESS24 is to improve the PV energy consumption during this period. Subsequently, the two ESSs also swing to discharge to support the voltage amplitudes and stability.  Figure 6 shows the reactive power of the PV converters. As can be seen, during the peak solar irradiation hours around noon (10:00-14:00), the PV converters can generate a small amount of reactive power limited by their capacities. At the other times, more reactive power is available. The PV reactive power outputs are adjusted frequently as the converters can be controlled flexibly and continuously. Due to the lack of reactive sources in Feeder2 and Feeder3, PV19 and PV24 inject a large amount of leading reactive power into the feeders. By contrast, PV5 regulates the power flow by generating both leading and lagging reactive power.   Figure 6 shows the reactive power of the PV converters. As can be seen, during the peak solar irradiation hours around noon (10:00-14:00), the PV converters can generate a small amount of reactive power limited by their capacities. At the other times, more reactive power is available. The PV reactive power outputs are adjusted frequently as the converters can be controlled flexibly and continuously. Due to the lack of reactive sources in Feeder2 and Feeder3, PV19 and PV24 inject a large amount of leading reactive power into the feeders. By contrast, PV5 regulates the power flow by generating both leading and lagging reactive power.   It is observed that the ESSs installed in Feeder1 and Feeder4 all charge before 8:00 and discharge from 15:00 to 21:15. As these two feeders have few active power sources, the ESSs tend to store a lot of energy to cope with the low voltage problem and guarantee the voltage stability via discharging during the peak load period. Feeder2 and Feeder3 are not as stressed as the others. ESS21 and ESS24 discharge in the early morning and then charge when the PV outputs rise between 7:15 and 14:45. Thus, the role of ESS21 and ESS24 is to improve the PV energy consumption during this period. Subsequently, the two ESSs also swing to discharge to support the voltage amplitudes and stability.

Effect of the Voltage Stability Constraint on Day-Ahead Scheme
The day-ahead scheme obtained by the proposed VSCOS model is regarded as the benchmark case and is compared with the solution without considering voltage stability. In terms of active power, compared with the benchmark case, the ESSs, excluding ESS21, generally start to discharge early in the afternoon, and their charging rates are high from 21:30 to 24:00, as shown in Figure 7. In terms of reactive power, it can be seen from Figure  5 that the reactive power compensation provided by the CBs on the over-compensated Feeder1 is significantly less than the benchmark case, while that on Feeder4 is reduced only in a few hours. As observed in Figure 6, the reactive power provided by the PV converters is similar to the benchmark case at most of the time spans. PV5 injects reactive power at close to its maximum value during the periods 1:00-4:30 and 22:00-22:30. Figure  4 shows that the OLTC scheme changes slightly compared with the benchmark case.
The model without the voltage stability constraint can find a solution with the total branch losses (the objective) no higher than our proposed method, but it results in unde- By tracking the H i and voltage of the test system over the time horizon, the weakest buses of every feeder are bus 17, bus 21, bus 23, and bus 32, respectively. The H i profiles of these buses are shown in Figure 8. The H i at bus 17 is significantly lower than that at the other buses, indicating that voltage instability is most likely to be associated with bus 17. The H i values are low during the heavy-load period and can be improved by the PV active power supply. The weak buses are usually at the end of the feeders and far away from the power sources, e.g., bus 17 and bus 32. The voltage stability can be enhanced if there are power sources near the end of the feeders, such as Feeder2 and Feeder3.

Effect of the Voltage Stability Constraint on Day-Ahead Scheme
The day-ahead scheme obtained by the proposed VSCOS model is regarded as the benchmark case and is compared with the solution without considering voltage stability. In terms of active power, compared with the benchmark case, the ESSs, excluding ESS21, generally start to discharge early in the afternoon, and their charging rates are high from 21:30 to 24:00, as shown in Figure 7. In terms of reactive power, it can be seen from Figure  5 that the reactive power compensation provided by the CBs on the over-compensated Feeder1 is significantly less than the benchmark case, while that on Feeder4 is reduced only in a few hours. As observed in Figure 6, the reactive power provided by the PV converters is similar to the benchmark case at most of the time spans. PV5 injects reactive power at close to its maximum value during the periods 1:00-4:30 and 22:00-22:30. Figure  4 shows that the OLTC scheme changes slightly compared with the benchmark case.
The model without the voltage stability constraint can find a solution with the total branch losses (the objective) no higher than our proposed method, but it results in undesired stability levels. As the comparison of these two cases shows, using the model with-

Effect of the Voltage Stability Constraint on Day-Ahead Scheme
The day-ahead scheme obtained by the proposed VSCOS model is regarded as the benchmark case and is compared with the solution without considering voltage stability. In terms of active power, compared with the benchmark case, the ESSs, excluding ESS21, generally start to discharge early in the afternoon, and their charging rates are high from 21:30 to 24:00, as shown in Figure 7. In terms of reactive power, it can be seen from Figure 5 that the reactive power compensation provided by the CBs on the over-compensated Feeder1 is significantly less than the benchmark case, while that on Feeder4 is reduced only in a few hours. As observed in Figure 6, the reactive power provided by the PV converters is similar to the benchmark case at most of the time spans. PV5 injects reactive power at close to its maximum value during the periods 1:00-4:30 and 22:00-22:30. Figure 4 shows that the OLTC scheme changes slightly compared with the benchmark case.
The model without the voltage stability constraint can find a solution with the total branch losses (the objective) no higher than our proposed method, but it results in undesired stability levels. As the comparison of these two cases shows, using the model without considering stability yields an objective function value 21.4 kWh lower than that of the benchmark case, but leads to the values of stability index H i below 0.9 during the period 22:00-22:30. However, it remains doubtful whether this situation is unique because the violation occurs during the off-peak period. Thus, we simulate two additional cases, raising and reducing all loads by 20%, and depicting the voltage stability level H i in both the cases at the weak bus 17 in Figure 9. It can be seen that there is almost no difference between the solutions of the two models if the loads are reduced by 20%. However, in the case where the loads are increased by 20%, the proposed VSCOS model can keep the stability level above the preset threshold by scheduling various devices, while the model without considering stability has an undesirable stability level during the peak-load hours between 19:15 and 20:30.
Energies 2021, 14, x FOR PEER REVIEW 15 of 21 benchmark case, but leads to the values of stability index i H below 0.9 during the period 22:00-22:30. However, it remains doubtful whether this situation is unique because the violation occurs during the off-peak period. Thus, we simulate two additional cases, raising and reducing all loads by 20%, and depicting the voltage stability level i H in both the cases at the weak bus 17 in Figure 9. It can be seen that there is almost no difference between the solutions of the two models if the loads are reduced by 20%. However, in the case where the loads are increased by 20%, the proposed VSCOS model can keep the stability level above the preset threshold by scheduling various devices, while the model without considering stability has an undesirable stability level during the peak-load hours between 19:15 and 20:30. In addition, it can be observed in Figure 6 that the PV converters sometimes provide a large amount of reactive power, for example from 1:45 to 4:30 in the case without the voltage stability constraint. Due to high power losses inside the converters, it is rare in practice, except in the following two possible scenarios: 1. The distribution network operator is the owner of the PVs or has full authority over the PVs. For the purposes of economic operation, the operator requires the PV converters to output a large amount of reactive power. The case study in this paper is consistent with this situation. 2. In a competitive market, the PV entities can be encouraged to output a large amount of reactive power by an effective compensation mechanism or market price signal. For the sake of profit, the active power generation can even be reduced to make way for reactive power. It is noteworthy that the PV converters may no longer adopt the maximum power point tracking control in this scenario.

Intra-Day Scheme under Cloudy Sky Conditions
In this scenario, ESS6 and ESS24 operate almost exactly according to the day-ahead scheme. The adjustments of the other ESSs are depicted in Figure 10. In addition, it can be observed in Figure 6 that the PV converters sometimes provide a large amount of reactive power, for example from 1:45 to 4:30 in the case without the voltage stability constraint. Due to high power losses inside the converters, it is rare in practice, except in the following two possible scenarios: 1.
The distribution network operator is the owner of the PVs or has full authority over the PVs. For the purposes of economic operation, the operator requires the PV converters to output a large amount of reactive power. The case study in this paper is consistent with this situation.

2.
In a competitive market, the PV entities can be encouraged to output a large amount of reactive power by an effective compensation mechanism or market price signal. For the sake of profit, the active power generation can even be reduced to make way for reactive power. It is noteworthy that the PV converters may no longer adopt the maximum power point tracking control in this scenario.

Intra-Day Scheme under Cloudy Sky Conditions
In this scenario, ESS6 and ESS24 operate almost exactly according to the day-ahead scheme. The adjustments of the other ESSs are depicted in Figure 10.
ESS6 is in the first half of Feeder1, where the PV power variations are not very significant with respect to the load level and can be addressed by the power provided upstream without rescheduling ESS6. The total load of Feeder3 accounts for 17-41% of the overall load, and hence the impact of the PV power fluctuations is also not significant. There is no need to adjust the scheme of ESS24 in Feeder3 to accommodate the PV outputs. In contrast, the peak load of Feeder2 is only 627.4 kW, which is close to the capacity of PV19. Thus, the intermittent PV power has a significant impact on this feeder and ESS21 should charge less energy between 8:00 and 16:00. It can also be seen that ESS15 and ESS30 discharge more between 7:30 and 13:45 to cover the decrease caused by the PV intermittency. In addition, the load deviation contributes to these adjustments as well. Figure 11 shows the real H i of the weakest bus 17. Due to the uncertainty, the direct utilization of the day-ahead scheme can lead to undesired power flow and H i profiles. For example, as shown in Figures 8 and 11, the real H i is worse than that calculated by the day-ahead model. A better H i profile can be achieved through the intra-day adjustments.  ESS6 is in the first half of Feeder1, where the PV power variations are not very significant with respect to the load level and can be addressed by the power provided upstream without rescheduling ESS6. The total load of Feeder3 accounts for 17-41% of the overall load, and hence the impact of the PV power fluctuations is also not significant. There is no need to adjust the scheme of ESS24 in Feeder3 to accommodate the PV outputs. In contrast, the peak load of Feeder2 is only 627.4 kW, which is close to the capacity of PV19. Thus, the intermittent PV power has a significant impact on this feeder and ESS21 should charge less energy between 8:00 and 16:00. It can also be seen that ESS15 and ESS30 discharge more between 7:30 and 13:45 to cover the decrease caused by the PV intermittency. In addition, the load deviation contributes to these adjustments as well. Figure 11 shows the real i H of the weakest bus 17. Due to the uncertainty, the direct utilization of the day-ahead scheme can lead to undesired power flow and i H profiles. For example, as shown in Figure 11 and Figure 8, the real i H is worse than that calculated by the day-ahead model. A better i H profile can be achieved through the intra-day adjustments.

Intra-Day Scheme under Clear Sky Conditions
The intra-day states of the ESSs under clear sky conditions are shown in Figure 10. Figure 12

Intra-Day Scheme under Clear Sky Conditions
The intra-day states of the ESSs under clear sky conditions are shown in Figure 10. Figure 12 depicts the real H i of bus 17 under clear sky conditions. Slight adjustments can be observed in this scenario, because the forecast errors of load and PV power are small. Similarly, the intra-day regulation achieves better H i profiles at most of the time spans. Sometimes the values of H i decrease (e.g., at some time spans from 9:45 to 14:30), but keep above the desired threshold.

Intra-Day Scheme under Clear Sky Conditions
The intra-day states of the ESSs under clear sky conditions are shown in Figure 10. Figure 12

Effectiveness of the Proposed Voltage Stability Index
To demonstrate the effectiveness of the proposed voltage stability index, a wellknown stability index, namely load margin, is introduced here to characterize the stability level. The model for load margin is consistent with (6) and the details can be found in [23]. The load margins for the day-ahead and intra-day scenarios are calculated where the dayahead scheme is employed here, as shown in Figure 13. Compared with Figure 8, as well as the day-ahead i H profiles in Figures 11 and 12, the trends of the load margins are almost identical to the profiles of i H at the weakest bus of the whole system, indicating the validity of our proposed index.

Effectiveness of the Proposed Voltage Stability Index
To demonstrate the effectiveness of the proposed voltage stability index, a well-known stability index, namely load margin, is introduced here to characterize the stability level. The model for load margin is consistent with (6) and the details can be found in [23]. The load margins for the day-ahead and intra-day scenarios are calculated where the day-ahead scheme is employed here, as shown in Figure 13. Compared with Figure 8, as well as the day-ahead H i profiles in Figures 11 and 12, the trends of the load margins are almost identical to the profiles of H i at the weakest bus of the whole system, indicating the validity of our proposed index.

Computational Time Comparison
A comparison of the computational time for the different models is shown in Tabl 1. No significant effect on the computational time is observed while incorporating the pro posed voltage stability constraints into either the day-ahead model or the intra-day model The longer time for solving the day-ahead model is due to a large number of integer var iables for the actions of the OLTCs and the CBs and the binary variables for the ESS charg ing/discharging states.

Computational Time Comparison
A comparison of the computational time for the different models is shown in Table 1. No significant effect on the computational time is observed while incorporating the proposed voltage stability constraints into either the day-ahead model or the intra-day model. The longer time for solving the day-ahead model is due to a large number of integer variables for the actions of the OLTCs and the CBs and the binary variables for the ESS charging/discharging states. In addition, the scalability of the proposed stability constraint is tested. We compare the computational efficiency of the model with and without the proposed voltage stability constraint on a larger-scale 141-node system. Its network and connectivity data can be found in [39]. For efficiency, only three types of devices are considered in this test system, i.e., PVs, CBs, and ESSs. See [40] for their detailed allocations. As the impact of the proposed constraint on the solution varies with the levels of demand, as stated in Section 4.2.2, the day-ahead scheme is solved 50 times in three cases with different loads, and the average solution time is shown in Table 2. In the table, each case is marked with its peak load value. It can be seen that in all the cases solving the stability-constrained model takes only less than 0.02 s longer, which is negligible. Therefore, the proposed voltage stability constraint is scalable. As only the solution time is the concern in this subsection, the other simulation results are uploaded to an online repository [40] and not shown here. The results lead to the same conclusion as the foregoing.

Effect of ESS on Voltage Stability
This subsection compares H i of the weakest bus in the test system with and without the ESSs, thereby verifying the effect of the ESSs on voltage stability. The comparisons are shown in Figure 14.  In the absence of the ESSs, the high PV outputs between 8:00 and 16:00 result in high values of i H , and during the subsequent heavy load period, the i H is significantly low compared to the case with the ESSs. Although sometimes the stability level in the case with the ESSs is slightly inferior to that without the ESSs, the result implies that the ESSs can help maintain the voltage stability at a higher and more balanced level over the time horizon. In the absence of the ESSs, the high PV outputs between 8:00 and 16:00 result in high values of H i , and during the subsequent heavy load period, the H i is significantly low compared to the case with the ESSs. Although sometimes the stability level in the case with the ESSs is slightly inferior to that without the ESSs, the result implies that the ESSs can help maintain the voltage stability at a higher and more balanced level over the time horizon.

Conclusions
This paper proposes a voltage stability constrained multi-time-scale optimal scheduling framework for ADNs. In this framework, the day-ahead model optimally regulates various types of devices in the ADN with the aim of minimizing the branch losses. The intra-day model takes the day-ahead scheme as a reference and adjusts the control actions of the ESSs and PV converters to cope with the uncertainty of the loads and renewable energy. To guarantee the voltage stability, a stability index is presented based on the Jacobian of the Distflow branch flow model and then a linear voltage stability constraint is further constructed. It is shown that the proposed index can track changes in the voltage stability level over the time horizon and the constraint outperforms existing ones in terms of weak bus identification, computational properties, and good scalability. The performance of employing the multi-time-scale model is also tested in different scenarios. It is verified through the simulations that the proposed multi-time-scale VSCOS model can effectively operate the ADN within a voltage stability region and at a desired distance from the loadability point. In particular, the ability of the ESSs to transfer energy helps achieve better and more balanced stability profiles.
The proposed VSCOS model is applicable to radial distribution networks and is especially effective in those with heavy load and abundant reactive power compensation. How to extend the proposed model to distribution networks with general arrangements and unbalanced power flow, and how to match the corresponding relationship between the proposed index and the voltage stability levels are suggested for future work.