Fast Solving Method for Two-Stage Multi-Period Robust Optimization of Active and Reactive Power Coordination in Active Distribution Networks

This paper considers constraints on the cycle number of energy storage systems (ESSs), travel distances of switchable capacitors reactors (SCRs), on load tap changers (OLTCs) and step voltage regulators (SVRs). The characteristics and operational constraints of these equipments are properly formulated with binaries and big-M method, obtaining desirable linear descriptions. The branch flow equations of distribution systems for active and reactive power coordination are constructed. Considering physical constraints, a multi-period mixed-integer deterministic second-order conic programming (SOCP) to minimize network losses is formulated. Then, considering uncertainties of loads and active power of renewable distributed generators (DGs), a two-stage robust model is developed. A directly iteratively solving method between the first and second stage models is proposed. In contrast to the traditional column-and-constraint generation (CCG) method, increases to variables and constraints are not needed to solve the first stage model. To solve the second stage multi-period model, only the model of each single period needs to be solved first. Then their objective function values are accumulated, and the worst scenarios of each period are concatenated. As a result, the solving complexity is greatly reduced. The capabilities of the proposed method are validated by three simulation cases. It is found that the computational rate using the proposed method is significantly enhanced with much less computer storage. Specifically, the simulation results of 4, 33 and 69-bus systems indicate that the computing rate of the proposed method is about 28 times higher than CCG method when 8 iterations are performed. Meanwhile, the precision of optimization results is also improved.


I. INTRODUCTION
A large number of DGs and electric vehicles (EVs) connected to distribution network bring about issues of frequent voltage fluctuation and violation, power congestion, high power loss, abandoning wind and solar power etc [1], [2]. The capability of distribution network to accommodate DGs and EVs may be The associate editor coordinating the review of this manuscript and approving it for publication was Siqi Bu . greatly restricted [3]. On the other hand, the installed ESSs, static VAR compensators (SVCs) and generators (SVGs), SCRs, OLTCs, SVRs and sophisticated communication system provide effective tools to regulate distribution system [4], [5]. The traditional passive distribution network is gradually evolving into active distribution network [6], [7]. However, as wind, photovoltaic power and EVs have strong randomness and intermittency, they make the voltage violation has characteristics of randomness and volatility. If the traditional autonomous control mode is still adopted for these controllable equipments, the interaction between them and DGs may result in frequent actions of OLTCs, SVRs and SCRs, speeding up wear and pushing up operation and maintenance costs [8]. As an outcome, the investments in grid resources are wasted and security of distribution grid is deteriorated.
To conquer these challenges, var/voltage optimization is usually used to regulate voltages magnitudes and reduce network losses [9]. However, the decoupling between active power and voltages magnitudes in distribution network does not exist anymore. The effect of active power on voltages magnitudes is significant in distribution network for high r/x of lines. In low voltage distribution network or distribution network with many cables, as the resistances of cables and low voltage lines are far greater than reactance, the effect of active power on voltages magnitudes is much larger than reactive power. In distribution network with high penetration of renewable DGs and EVs, sometimes it is impossible to solve the voltages magnitudes overrun issues by only regulating reactive power. Therefore, a coordinated optimization of active and reactive power must be conducted especially when there are ESSs connected to distribution system. The multiperiod scheduling of active and reactive power coordination considering day-ahead data and characteristics of equipments from a prospective vision is an effective tool to ensure safe and economic operation of distribution network [10].
Traditionally, the coordinated optimization of active and reactive power could be formulated as a mixed-integer nonlinear programming. A lot of techniques are proposed to solve it. In the light of potential economic benefits, there are growing interests on obtaining global optimal solution. In recent years, important progresses on convex relaxation methods are released. Semidefinite programming (SDP) was introduced in [11] to solve the optimal power flows (OPF) for mesh balanced and unbalanced networks based on bus injection model. Compared to SDP, the SOCP based on bus injection model as well as branch flow model has attracted more attentions in radial balanced distribution networks for its much higher solving efficiency [12], [13], [14], [15]. While SDP is widely used in three-phase balanced and unbalanced distribution systems, most of literatures formulate SOCP of OPF in balanced radical distribution systems. Recently, reference [16] first formulates a SOCP based on bus injection model in unbalanced distribution systems. The solving efficiency and exactness are also analyzed. Moreover, some sufficient conditions are given in [17], [18], [19], and [20] for radial networks. For the active distribution system with high penetration of DGs, the exactness of the SOCP relaxation is verified in [20] using four IEEE benchmark and practical systems. But only under some mild conditions the relaxation is exact (e.g., the objective function is convex and monotonically increasing with each active power injection and the original OPF is feasible).
In [21], it is pointed out that the traditional SOCP and the conditions for exactness are no longer applicable when there are coaxial cables in the distribution network since the shunt capacitors of coaxial cable cannot be ignored. The OPF model of distribution networks with cables is formulated. A set of sufficient conditions for relaxing the model to SOCP exactly which can be checked ex ante are developed. In [22], a sequential SOCP is developed when the original SOCP is not exact. Simulation cases indicate the proposed method outperforms commercial nonlinear solvers in terms of robustness and efficiency without compromising precision. Considering that the solving efficiency of mixed-integer SOCP is also very high, it is feasible to formulate a mixed -integer SOCP for the OPF of distribution network with discrete regulating equipments.
To address the uncertainty of loads and renewable DGs, two major methods are widely applied, i.e., stochastic and robust optimization. The robust optimization methods have advantages in terms of full solution robustness and high computing efficiency [23]. Thus, to guarantee operating robustness against uncertainty, the robust optimization methods are preferred. At present, the CCG algorithm is widely used in solving robust model [7], [9], [10], [24], [25]. In [1], a robust day-ahead scheduling model of smart distribution networks considering demand response programs is presented. In [2], a multi-interval-uncertainty constrained robust dispatch model is proposed. The uncertainty budget is rationally divided according to distribution probabilities to improve conservativeness of traditional robust models. To address the min-max-min robust model with a mixed-integer recourse problem, a nested CCG algorithm is adopted to quickly obtain the minimum operating cost in the worst-case scenario. In [3], a two-stage robust optimization for expansion planning of active distribution systems coupled with urban transportation networks is proposed. The interactions between transportation and active distribution system are explored. In [4], a novel tri-level robust planning-operation co-optimization model to determine the capacity, power, location and scheduling strategy of distributed ESSs is proposed. In the above four literatures, the CCG algorithm is used to solve the robust optimization model. However, to solve the first stage model using CCG method, more and more variables and constraints are generated during iterations. Therefore, although CCG algorithm has been validated by many literatures, the computational rate of the CCG algorithm may be very slow when the number of iterations becomes large. To fill this gap, we propose a new method to quickly solve the two-stage robust optimization problem.
According to the above literature review, it is found that fast robust optimization method is still needed to develop. As far as the knowledge of the authors, there are no works that modified the traditional CCG method in terms of computing rate. We are the first team that has proposed a method modifying the traditional CCG algorithm in terms of enhancing computing rate greatly. In this paper, a two-stage multi-period mixed-integer second-order cone robust model for coordinated active and reactive power optimization of distribution VOLUME 11, 2023 network is formulated based on branch flow equations. A fast solving method is proposed. Compared to existing works, the novelty and contributions of this paper are summarized as follows. 1) We proposed a method to remove the absolute value symbol |·|. To this end, the characteristics and operational constraints of OLTCs, SVRs, SCRs were properly formulated with binaries and big-M method, obtaining desirable linear descriptions. 2) In contrast to CCG method, increases to variables and constraints are not needed to solve the first stage model using our proposed method. As a result, the number of optimization variables and constraints of the first stage model of the proposed method are much less than CCG method. 3) To solve the second stage multi-period model, only the model of each single period needs to be solved first. Then their objective function values are accumulated and the worst scenarios of each period are concatenated. As an outcome, solving complexity of the second stage model is greatly reduced. 4) Overall, the computational rate of the proposed method is much faster than CCG method with higher precision for optimization results and less computer storage.
The structure of this paper is as follows. In section II, the model for dynamic optimization of active and reactive coordination in active distribution network is formulated. In section III, a two-stage robust model is developed. In section IV, the fast solving method is proposed. In section V, simulation cases are performed and analyzed. Section VI concludes the paper.

A. BRANCH FLOW MODEL OF RADIAL DISTRIBUTION NETWORK
This paper focuses on the operation problems of distribution systems. We assume that the connecting locations of regulating equipments are already known. As shown in Fig. 1, the radial distribution network can be described by branch flow equations. For node j, the power equations are given by: where P ij,t and Q ij,t are the active and reactive power entering the top of line ij at time period t, respectively. r ij and x ij are the resistance and reactance of line ij, respectively. P DG j,t and P ESS j,t are the injected active power at node j time period t from DG and ESS respectively. Q SVC SCR, respectively. P Load j,t and Q Load j,t are the active and reactive power of loads at node j, time period t, respectively. I ij,t is the current magnitude of line ij at time period t. δ(j) is the set of nodes whose upstream node is j. (j) is the set of nodes whose immediately downstream node is j. b s,j is the sum of the susceptance of the line connected to node j. U j,t is the voltages magnitude of node j at time period t. The voltage magnitude drop across branch ij without any transformer is given by: where E is the set of lines. TR is the set of branches with transformers.
As shown in Fig.2, the square of voltage magnitude drop across branch ij with a transformer is given by: where t ij,t is the transformer turn ratio at time period t.
Given l ij,t = I 2 ij,t and u i,t = U 2 i,t , Eqs. (1) and (2) can be converted into: Eq. (3) can be converted into: A new variable is defined as u Tj,t = t 2 ij,t u j,t . Eq. (4) can be converted into: (8) For branch ij, the relationship between squared voltage and current magnitudes, active power, and reactive power is given by: 30210 VOLUME 11, 2023 Eq. (9) is nonconvex, conic relaxation can be performed by relaxing the quadratic equality into inequality as: Thus, in Eqs. (5)∼(8) and inequality (10), only u Tj,t = t 2 ij,t u j,t is non-convex, it will be linearized in the model of SVR and OLTC.

B. MODEL OF ESS
To mitigate fluctuations of photovoltaic and wind power, an ESS may be installed in the active distribution network. Its mathematical model is formulated as follows.
where P disch arge j,t and P dischargemax j are the discharge power of ESS at node j, time period t and its maximum, respectively. P are the discharge and charge states of ESS at node j, time period t, respectively.

2) CAPACITY CONSTRAINTS
where E ESS j,t , E ESSmin j and E ESSmax j are the stored electric energy at the beginning of time period t and its lower and upper limits, respectively. η charge and η discharge are the charge and discharge efficiency. T max is the scheduling interval. t is the total number of scheduling periods. In Eq. (17), to ensure that the ESS has the same regulation performance in the new scheduling cycle, the electric energy stored in ESS at the initial time of a dispatching cycle E ESS j,t is equal to that of the next cycle E ESS j,T max .

3) CHARGE AND DISCHARGE STATE CONSTRAINT
At any time period, the ESS can be in one of three states. That is, charging, discharging, or neither charging nor discharging. Thus, the states of charge and discharge must meet the following constraint: (20) where N ESS is the upper cycle limit of the ESS during the scheduling periods. E j is the electrical energy stored in the ESS when it is fully charged.

C. MODEL OF SVC AND SVG
where Q SVCmin j and Q SVCmax j are the upper and lower limits of reactive power injected from the SVC connected to node j at time period t, respectively.
where Q SVGmax j and Q SVGmin j are the upper and lower limits of reactive power injected from the SVG connected to node j at time period t, respectively.

D. MODEL OF SCR 1) LINEARIZATION OF INJECTING REACTIVE POWER FROM SCR
The reactive power injected from the SCR connected to node j at time period t is: where C j,t is the capacitive admittance of SCR. Eq. (22) can be linearized as follows. Let the total number of SCR groups be and C max j are the minimum and maximum capacitive admittance of SCR. s j is the step size of capacitive admittance. Let v j + 1 be the length of the binary number (an integer number built using 0 and 1) of SCR groups that are put into operations. Thus, the integer v j can be determined as: The capacitive admittance of SCR can be expressed as: The reactive power injected from SCR is given by: Define σ jv j ,t = z jv j ,t u j,t . By introducing dummy variables σ j0,t , σ j1,t . . . σ jv j ,t , Eq. (27) can be replaced with Eq. (28). Considering z jv j ,t is a binary and u j,t is a positive number, the bilinear expression σ jv j ,t = z jv j ,t u j,t can be linearized with inequalities (29) and (30) using the big-M method.
where M is a sufficiently big number. Its value depends on the range of the value taken by the substituted variable.

2) LINEARIZATION OF SCR TRAVEL DISTANCE CONSTRAINT
Between time period t and t + 1, the number of changes for SCR groups put into operation is given by: The absolute value sign in Eq. (31) is very hard to process and can be removed by introducing a binary dummy variable ρ j,t . Thus, Eq. (31) can be equivalently replaced with: Define is negative which implies v j n=0 2 n b z jn,t+1 − z jn,t is also negative, then the binary ρ j,t must equal to 0 since X j,t ≥ 0 and X j, can be linearized with inequalities (36) and (37) using the big-M method.
The travel distance constraint in a scheduling cycle is given by: where X max j is the upper limit of the total travel distance of SCR in a scheduling cycle.
The variable u Tj,t in Eq. (8) can be linearized as follows. The turn ratio of the SVR at time period t can be expressed as [26]: where t min ij and t max ij are the minimum and maximum turn ratio. t ij is the step size of the turn ratio. The discrete variable T ij,t is the tap position at time period t. K ij is the highest tap position.
Let LB be represented by a binary number. Eqs. (39) and (41) can be expressed as [26]: where N ij is the length of the binary number representing T ij,t . It can be determined similarly to (23). By multiplying the left and right sides of Eq. (42) by u j,t and defining new variables m ij,t = t min ij u j,t and x ij,n,t = λ ij,n,t u j,t , it yields [26]: Considering λ ij,n,t is a binary and u j,t is a positive number, the bilinear expression x ij,n,t = λ ij,n,t u j,t can be linearized with inequalities (46) and (47) using the big-M method [26].
Similarly, by multiplying the left and right sides of Eq. (42) by m ij,t and defining new variable y ij,n,t = λ ij,n,t m ij,t , and remember that u Tj,t = t 2 ij,t u j,t and m ij,t = t ij,t u j,t , it yields [26]: Considering λ ij,n,t is a binary and m ij,t is a positive number, and using the big-M method, the bilinear expression λ ij,n,t = λ ij,n,t m ij,t can be linearized with the inequalities [26]: Thus, the linearization of u Tj,t is completed.

2) LINEARIZATION OF TRAVEL DISTANCE CONSTRAINT FOR THE SVR TAP
The tap position difference of the SVR between adjacent time periods t and t + 1 is given by: The absolute value function in Eq. (51) can be removed by introducing a dummy binary τ ij,t , given by: By defining γ ij,t = τ ij,t N ij n=0 2 n+1 λ ij,n,t+1 − λ ij,n,t and using the big-M method, Eq. (53) can be linearized as: The total travel distance constraint on the tap of the SVR in a scheduling cycle is: where W max ij is the upper limit of the total tap travel distance in a scheduling cycle.

F. MODEL OF OLTC 1) LINEARIZATION OF TURN RATIO SQUARE FOR OLTC
Since the primary side voltage magnitude of OLTC is a fixed value at each time interval, the model of OLTC is different from SVR. Less variables and constraints than SVR are needed to formulate a linear model of OLTC. The turn ratio of the OLTC at the root node is expressed as: where k min and k max are the minimum and maximum turn ratio. k is the step size of the turn ratio. The discrete variable T t is the tap position at time period t. K is the highest tap position.
Let T t be represented by a binary number. Eqs. (59) and (61) can be expressed as: where L is the length of the binary number representing the tap position. It can be determined similarly to (23). The square of secondary side voltage magnitude for the OLTC can be expressed as: where r ij is the square of the root node voltage magnitude. By defining Y n,t = k t β n,t , Eq. (65) can be expressed as: Thus, the linearization of u Tj,t = u 0 k 2 t is completed.

2) LINEARIZATION OF TRAVEL DISTANCE CONSTRAINTS FOR OLTC
The tap position difference of the OLTC between adjacent time periods t and t + 1 is: 2 n β n,t , t = 1, 2, · · · T max − 1 (69) VOLUME 11, 2023 The absolute value function can be removed by introducing a dummy variable µ t . Eq. (69) can be expressed as: By defining α t = µ t L n=0 2 n+1 β n,t+1 − β n,t , Eq. (71) can be linearied as: The total travel distance constraint of tap position in a scheduling cycle is: where O max is the upper limit of the total tap travel distance in a scheduling cycle.
where U max j and U min j are the upper and lower limits of the voltage magnitude at node j, respectively.

H. OBJECTIVE FUNCTION
At present, the utilities are willing to achieve multiple objectives with a single solution. Nevertheless, network losses minimization can ensure the exactness of the relaxed SOCP when DGs penetration is not very high. However, if other objectives are added, voltages magnitudes deviation minimization for example, the exactness of the relaxed SOCP may not be guaranteed any more. To this end, the objective function is to minimize total network losses given by: DSTATCOM functionality of DGs and active power curtailment are not considered in this paper. However, if they are considered, the objective function must be modified as total network losses minus total DGs active power output. One second-order conic constraint of active and reactive power and one box constraint for active power for each DG should be added. The model of maximizing hosting capacity of uncertain photovoltaic by coordinated management of OLTC, reactive power sources and stochastic EVs can be formulated by modifying the objective function as total network losses minus total photovoltaic active power too. The models of smart charging of EVs must be added. The proposed fast solving method can be still applicable.
It is worthy of pointing out that Eqs. (11)∼(20) related to ESSs are self sufficient. The models of ESSs, other electric equipments are involved in the branch flow equations of distribution network as a whole. All the physical constraints related to ESSs, other equipments and distribution network are taken into account in the proposed mixed integer SOCP.
Although some of those models in Section II can be found in other papers. Their origins are indicated with references. Nevertheless, we propose a method to model OLTCs with less dummy variables than SVRs. We consider constraints on the cycle number of ESSs, travel distances of SCRs, OLTCs and SVRs. These constraints and the model of OLTCs are properly formulated with binaries and big-M method, obtaining desirable linear descriptions. They cannot be found in other papers.

III. ROBUST OPTIMIZATION MODEL A. DETERMINISTIC OPTIMIZATION MODEL
A two-stage optimization method is adopted in this paper. The discrete variables associated with OLTCs, SVRs, SCRs and the discrete as well as continuous variables associated with ESSs are the first stage variables. The first stage variables denoted by ψ t cannot be adjusted after the uncertainties are revealed. They are regarded as the ''here and after'' decisions. The continuous variables, such as reactive power of the SVC, branch power, currents, nodes voltages magnitudes are the second stage variables. The second stage variables denoted by ϕ t can be adaptively changed in the real-time operation when the uncertainties are revealed. Branch power flows, currents and nodes voltages magnitudes are the state variables and the reactive power of SVC are the ''wait and see'' decisions.
For a better illustration, the optimization model can be formulated in a compact form, given by: It is worthy of pointing out that the multi-period correlation coupling constraints only involve ESSs. The optimization variables related to ESSs have been set as the optimization variables in the first stage. There is no coupling relationship between the second stage optimization variables in different periods. Therefore, the second stage multi-period model only needs to be solved for each period. Therefore, the solution complexity can be greatly reduced.

B. ROBUST OPTIMIZATION MODEL
As intermittent DGs and loads have strong randomness, it is challenging to adjust all the voltages or currents magnitudes within rated ranges if the deterministic optimization model is adopted. Thus, a robust optimization model is constructed as: In the uncertainty sets (83c) and (83d), the lower and upper bounds d min t , d max t , g min t , g max t for load demand and output power of DGs are obtained by the interval prediction. These bounds limit the variation range of each uncertain variable d t and g t . The physical meaning of Eq. (83a) is that the first stage solution minimizes network losses in the worst scenario. The inner max-min model is to find out the worst scenario with the largest losses.
As the same as [7], [9], and [10], the temporal and spatial correlation of load, wind power and photovoltaic output uncertainties are not taken into account in (83c) and (83d). However, the degree of uncertainty can be reasonably set by using uncertainty adjustment parameters to take it into account. When it is considered, the method proposed in this paper is still applicable. The conservatism of the optimization results can be reduced at the cost of expanding the scale of the optimization problem.
Compared with the conventional deterministic optimal power flow problem, after considering the uncertainties, the distribution network operates in a worse scenario. However, if and only if the robust optimization problem represented by equation (82) is feasible, and the distribution network line parameters meet the conditions C1 and C2 in [20], despite degree of uncertainties, the exactness of the SOCP can be guaranteed. In addition, C1 and C2 are mild conditions, easy to meet, and can be checked in advance. However, if the distribution network line parameters do not meet the conditions C1 or C2, the convex relaxation gap of each branch can be calculated after solving the SOCP. If the convex relaxation gaps of all branches are less than tolerance, then the convex relaxation is exact. If the convex relaxation gap of any branch is greater than tolerance, the convex relaxation is inexact. To this end, a sequence SOCP can be constructed using the method in [22]. Although its solving efficiency is lower than SOCP, the (near) global optimal solution can still be recovered, and the solving efficiency is still high.
It is worthy of pointing out that although in [7], [9], and [10], a two-stage robust SOCP is constructed, the strong duality of SOCP problems usually requires certain conditions to meet [27], [28]. If the strong duality condition is not satisfied, add the allowable load shedding (set the upper limit of the load shedding sufficiently large) in the constraint conditions, and add the penalty term for the load shedding (set the penalty coefficient to a large positive number) in the objective function, then the strong duality of SOCP can be met on the premise of ensuring the accuracy of the model [27], [28].

IV. SOLVING METHOD A. MASTER AND SUB-PROBLEM OF PROPOSED METHOD
The master problem is shown in (82a)∼(82g). The objective of the master problem is to find the optimal values of the first stage variables given the worst scenario generated by the sub-problem. It gives a lower bound of the original problem given by (83a)∼(83d). The optimal values of the first stage variables such as the capacitive admittance of the SCR, the turn ratios of the OLTC and SVR, and the injected power of the ESS, are substituted into the power flow Eqs. (5)∼(8) and conic relaxation (10). As an outcome, the sub-problem can C. PHYSICAL MEANING OF THE PROPOSED SOLVING METHOD Similar to CCG method, given an initial scenario, the first stage variables are optimized by solving the MP. Then their setting values are passed to the SP. After that, the worst scenario is generated by solving the SP. Then they are passed to the MP. After that, the first stage variables are optimized again by solving the MP. This process is repeated until the end condition is met. During iterations, the scenario generated by the SP becomes worse and worse. Moreover, the objective function value of the MP becomes larger and larger. Meanwhile, the optimized values of the first stage variables become better and better. As a result, the objective function value of the SP becomes smaller and smaller. Finally, when the gap between the objective function values of SP and MP is less than the tolerance or the maximum iterations are reached, the program is terminated.

V. SIMULATION CASES
To verify the high computing efficiency, three test cases including 4 and 33, 69-bus distribution systems available from [30] are analyzed using the proposed method considering the uncertainty of loads and wind power. All the programs are developed in MATLAB R2018a. The MOSEK package is used to solve the MP and SP. The hardware configuration of the computer is Intel dual core i3-4150M with 3.5GHz CPU, and 32 GB memory. All the figures are plotted with VISIO or MATLAB. In order to make the program run smoothly, select the appropriate relative dual gaps, which could be much larger than the default value 0.0001 for large systems when necessary.
The maximal active power of each DG is set to 1, 0.2, 0.2 MW in the 4, 33 and 69-bus systems respectively. The total capacities of DGs are 1 × 1 MW = 1MW, 4 × 0.2 MW = 0.8MW and 5 × 0.25 MW = 1.25 MW in the 4, 33 and 69-bus systems respectively. Voltages magnitudes boosting problem is not severe. The regulating equipments can mitigate it using the proposed method. The simulation conditions for each case are selected mostly according to past literatures and the rated capacities of distribution systems.
A. CASE 1 1) SIMULATION CONDITIONS A 4-bus system shown in Fig. 3 is used to test the precision of the proposed method. The resistance and inductance of all three lines are set to 0.193 ohm and 0.38 mH respectively. There is not any OLTC installed in this distribution system. Other simulation conditions are set as follows.
a) The capacity of the ESS at bus 2 is 0.5 MWh. Both the charging and discharging efficiency of the ESS are set to 0.9. Both the maximum charging and discharging power are set as 50 kW. The minimum and maximum stored energies are 0.05 MWh and 0.5 MWh, respectively. Since battery life is greatly affected by cycle number, the maximum number of cycles from 0:00 to 24:00 is set to 3. b) The capacity of SVC at bus 4 is [−0.3 0.3] MVar. c) The minimum and maximum turn ratios of the SVR connected between buses 2 and 3 are 0.94 and 1.06, respectively. The step size is 0.01. d) The maximum and minimum compensation reactive power of the SCR at bus 3 are 0.6 and 0 MVar, respectively. The step size is 0.05 MVar. e) The maximum travel distance of the SVR and SCR is set as 24. f) The prediction curve of active power for load at buss 2 to 4 is shown in Fig. 4. The maximum load is 1 MW. The power factor is 0.95.
g) The prediction curve of wind generator (WG) at bus 4 is shown in Fig. 5. The maximum power is 1MW. The power factor is 1.0.   h) The base voltage magnitude and power are chosen to be 24.9 kV and 1 MVA, respectively. i) The root node is the slack node. Its voltage magnitude is fixed to 1.0 p.u.. j) The upper and lower limits of voltages magnitudes are 1.1 and 0.9 p.u., respectively. k) The upper current limit of each line is 120 A. l) The optimization time periods are from 0:00 to 24:00.
The time interval of optimization is set to be 1 hour. m) The optimization tolerance ε is set to be 0.0002. n) The big M in the MP and SP are set to be 10,000 and 100, respectively. o) The number of maximum iterations n max is set to be 8. p) The uncertainty interval is given by [ where P f is the forecasted load or wind power output and α reflects the confidence interval [9].

2) SIMULATION RESULTS
Using the proposed and improved CCG method (the SP is modeled and solved with the proposed method in this paper), the optimization results of reactive power compensation for VOLUME 11, 2023 the SCR at each time period under different α are exactly the same, 0.6 MVar. Moreover, the optimization results of the turn ratio for the SVR at each time period under different α are the same, 1.06. The total network losses under different α are shown in Table 1. It can be seen that, for each α, corresponding losses are almost the same. When α is 0.2, the optimal charging/discharging power and state of charge of ESSs using the proposed and improved CCG methods are shown in Fig. 6 and Fig. 7 respectively. It can be seen that, the two curves nearly overlap. Moreover, the ESS is charged when the loads are low and the wind power is high, and discharged when the loads are high and the wind power is low. To this end, it concludes the precision of the proposed method is high enough.  The iterations and computational time of the proposed and improved CCG methods under different α are shown in Table 2. It can be seen that, due to the small scale of the simulation system, the computational time of both methods is tens of seconds. However, the iterations of the proposed method are less than the improved CCG method.
The maximum and minimum voltages and maximum current using the proposed method and CCG algorithm for the 4-bus system when prediction error α is 0.2 are shown in Fig. 8 and Fig. 9, respectively. As can be seen, all the voltages and currents are within the specified ranges. However, the maximum current during time period 19 when load is maximum while output of DG is low is close to the upper limit 120 A. Therefore, the current rather than voltage is really a key limiting factor for safe operation of this distribution system.

B. CASE 2 1) SIMULATION CONDITIONS
In this section, the IEEE 33-bus distribution network shown in Fig. 10 is used to test the performance of the proposed method. For the sake of fairness, the selection of relative dual gap is the same as that of CCG algorithm. In order to make the program run smoothly, select the appropriate relative dual gaps, which could be much larger than the default value 0.0001. Some simulation conditions are set as follows.
a) Two ESS are connected to nodes 17 and 33 whose settings are shown in Tab.1 of [7]. b) One SVC connected to node 16, whose capacity is [-0.5 0.5] MVar. c) One OLTC is connected to the root node. The maximum and minimum turn ratios are 1.06, 0.94 respectively. The step size is 0.01. d) Two SCRs are connected to nodes 3, 6. The maximum and minimum reactive power are 0.6, -0.6 MVar respectively. The step size is 0.01 MVar. e) The prediction curve of active power for load at nodes 2 to 33 is shown in Fig. 4. f) Four wind turbines are connected to nodes 13, 21, 24, 31. The prediction curve of wind power is shown in Fig. 5. The maximum power is 0.2 MW. The power factor is 1. g) The base voltage magnitude and power are chosen to be 12.66 kV and 10 MVA respectively. h) The upper current limit of each line is 400 A. i) Optimization tolerance 0 ≤ T t ≤ K is set to be 0.001. j) The big M in the master and sub problems is set to be 10,000. k) There is not any SVR in this system. l) The tolerances of relative gaps are properly chosen such that the programs can run fluently. For fairness of comparison, they are set the same for the proposed and CCG methods. Other simulation conditions are set similarly to case 1.

2) SIMULATION RESULTS
The total network losses, iterations and computational time under different α are shown in Table 3 and Table 4, respectively. It can be seen that the network losses of the proposed method are less than or the same as the improved CCG method for each α. Moreover, the computing time of the proposed method is less than the CCG method, especially when α is 0.1.  When α is 0.6, the charging/discharging power and state of charge of ESS1 and ESS2 with the two methods are shown in Fig. 11 ∼ Fig. 14, respectively. It can be seen that curves are nearly overlapped. Further, both ESSs are charged when the loads are low and the wind power is high, and discharged when the loads are high and the wind power is low. Therefore, the precision of the proposed method is very high.  The maximum and minimum voltages and maximum current using the proposed method and CCG algorithm for the 33-bus system when prediction error α is 0.6 are shown in Fig. 15 and Fig. 16, respectively. As can be seen, all the voltages and currents are within the specified ranges. However, the maximum current and minimum voltage during period 19 when load is maximum while the output of DG is low is close to the upper limit 400 A and lower limit 0.9 p.u. respectively. Therefore, both current and voltage are key limiting factors for safe operation of this distribution system when the prediction error α is high.

C. CASE 3 1) SIMULATION CONDITIONS
In this section, the modified PG69-bus distribution network shown in Fig.17 is used to further verify the capability of the proposed method. For the sake of fairness, the selection of relative dual gap is the same as that of CCG algorithm. In order to make the program run smoothly, select the appropriate relative dual gaps, which could be much larger than the default value 0.0001. The connecting locations of regulating equipments are shown in Table 5. The maximum power of each WT is 0.25 MW. Other specifications of WTs are the same as f) of case 2. Other simulation conditions are similar to case 2. They are not listed anymore to avoid repetition. Since the developed mixed integer SOCP is a 24-hours multiperiod model, there are tens of thousands of variables and constraints in it. Therefore, it can be considered as a large system. In [7], the 69-bus system is also used as the biggest simulation system.

2) SIMULATION RESULTS
The tolerance ε is set to 0.01. The network losses and the iterations as well as computational time using the proposed method and the improved CCG method are shown in Table 6 and Table 7, respectively. It can be seen that the network losses for each α using the proposed method are lower than the improved CCG method. Further, the computational rate    of the proposed method is about 3.5 times faster than the improved CCG method.
To enhance the optimization precision, the tolerance ε is set to 0.0001. The optimization results are shown in Table 8 and  Table 9. It can be seen that, the corresponding network losses of the two methods are close. Nevertheless, the computational rate of the proposed method is about 28 times faster than the improved CCG method. This is because in the CCG method, a large number of new dummy variables and constraints are generated in MP during iteration. Therefore, the computational rate of the CCG method becomes very low when lots of iterations are performed to reach convergence for large scale distribution network. Nevertheless, in our proposed method, the number of constraints and variables for master problem remains constant during iteration. Further, the sub-problem only needs to be solved for each time period. As a result, the computational rate of the proposed method is much higher than the CCG method.

VI. CONCLUSION
In the CCG method, a large number of new dummy variables and constraints are generated in the model of master problem during iteration. Therefore, the computational rate of the CCG method becomes very low when lots of iterations are performed to reach convergence for large scale distribution network. In this paper, a fast robust optimization method is proposed, in which the number of constraints and variables for master problem remains constant during iteration. Further, the sub-problem only needs to be solved for each time period. Then their objective function values are accumulated, and the worst scenarios of each time period are concatenated. As a result, the computational rate of the proposed method is much higher than the CCG method. The precision of optimization results with the proposed method is also improved. Specifically, the simulation results of the 69-bus system indicate that the computing rate of the proposed method is about 28 times higher than CCG method when 8 iterations are performed. Therefore, the proposed method can be applied to large scale distribution systems that have thousands of nodes.
It is worthy of pointing out that, since the developed mixed integer SOCP is a 24-hours multi-period model, there are tens of thousands of variables and constraints for the 69-bus system. In [7], the 69-bus system is also used as the biggest simulation system. Therefore, the 69-bus system is large enough to validate the capabilities of the proposed method. To this end, the above conclusions are self-sufficient.
Whether the proposed method is valid for other types of uncertainty sets, irregular and nonconvex uncertainty sets for example, is also a topic worthy of study in future. Another future work is to compare the results with practical/hardware in loop models to validate the capabilities of the proposed method. Further, in future work, DSTATCOM functionality, active power curtailment of DGs and EVs requirements will be considered. Moreover, Enhancing hosting capacity of intermittent wind turbine systems using bi-level optimization considering OLTCs, SVRs, SCRs, SVCs, SVGs, ESSs and EVs charging stations will be considered too.