# A Multistage Solution Approach for Dynamic Reactive Power Optimization Based on Interval Uncertainty

Authors: Xiaodong Shen; Yang Liu; Yan Liu
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

In order to solve the uncertainty and randomness of the output of the renewable energy resources and the load fluctuations in the reactive power optimization, this paper presents a novel approach focusing on dealing with the issues aforementioned in dynamic reactive power optimization (DRPO). The DRPO with large amounts of renewable resources can be presented by two determinate large-scale mixed integer nonlinear nonconvex programming problems using the theory of direct interval matching and the selection of the extreme value intervals. However, it has been admitted that the large-scale mixed integer nonlinear nonconvex programming is quite difficult to solve. Therefore, in order to simplify the solution, the heuristic search and variable correction approaches are employed to relax the nonconvex power flow equations to obtain a mixed integer quadratic programming model which can be solved using software packages such as CPLEX and GUROBI. The ultimate solution and the performance of the presented approach are compared to the traditional methods based on the evaluations using IEEE 14-, 118-, and 300-bus systems. The experimental results show the effectiveness of the presented approach, which potentially can be a significant tool in DRPO research.


---

## Body

## 1. Introduction

Reactive power optimization’s importance in enabling the operation safety in the power system has been proved. The traditional reactive power optimization researches consider a number of elements including the power balance constraints, the voltage magnitudes, the branch currents, the coordinated control of the reactive power compensators, and the transformer tap ratios under a given load level. However, along with the number and capacity of the renewable resources increase, new challenges, for example, the overvoltage and power loss, gradually appear. These issues significantly lead to the bidirection of power flow. As a result the traditional reactive power optimization model should be improved. From the mathematical point of view, the reactive power optimization can be formulated as the mixed-integer/nonlinear programming model, which is mainly categorized into the static reactive power optimization (SRPO) [1–3] and the dynamic reactive power optimization (DRPO) considering multiperiod coupling [4–9]. In research [3], an improved particle swarm optimization algorithm using eagle strategy (ESPSO) is proposed for solving reactive power optimization problem to minimize the power losses. The SRPO only focuses on single system state in certain time period without considering the potential significant change from successive operating points. For example, discrete controllable devices (e.g., capacitors/reactors and transformer taps) are controlled by their switching on/off, the daily operating times of which are limited by their service lifetime and existing manufacture techniques. As a result, researches pay more attention to DRPO which significantly contributes the online and practical operations of the power systems. A day-ahead voltage stability constrained dynamic optimal reactive power flow (VSC-DORPF) model is proposed in research [4], in which the discrete control variables and the time coupled constraints are handled by the proposed branching and pruning principles. Considering the intertemporal constraints on the operating times of switching discrete variables and the nonlinear power balance equations, DRPO can be presented by large-scale mixed-integer/nonlinear programming problems (mixed-integer nonlinear programming, MINLP).However, the MINLP based DRPO has difficulties in modeling and solving. The reasons mainly include the following: (1) DRPO contains the continuous and discrete variables as a typical mixed integer nonlinear nonconvex problem so that its global optimal solution cannot be obtained due to the nonconvex power flow equations [5–9]; (2) DRPO specially considers the intertemporal constraints for the operating times of adjusting discrete reactive power compensation devices. This point makes the solution become computational intensive and time consuming [6–10]. In order to facilitate the DRPO solution, a number of efforts have been done. Round-off approach [11–13] which is a simplified method in dealing with discrete variables is claimed as an underlying tool for serving the DRPO solution. The discrete variable is firstly relaxed into continuous variables in the optimization model. Then, the continuous variables obtained from the optimization are rounded off to the nearest integer value. Due to the numerical approximation, such approaches may introduce unnecessary overhead in the objective functions and even cause constraint violations. The round-off approach is particularly unsatisfactory for the maximum allowable daily operating times of the electric devices. Some other researchers presented that the artificial intelligence [14, 15] can deal with the discrete variables and the limitation of the number of actions in reactive power optimization well [10, 16]. However, the algorithms employed in their researches show low efficiency and unstable optimization results. Heuristic approaches are regarded as better choices to deal with the discrete variables and the limitation of the number of actions. A number of heuristic approaches segment the network loss or daily load curve according to the trend of the load changes. And then the approaches make the number of the segments equal to the upper limit of daily operating times for the control devices [17, 18]. However, how to determine the position of the segmentation point in this method lacks theoretical basis and cannot guarantee getting the optimal segmentation scheme. Moreover, research [19] puts forward a novel particle swarm optimization algorithm with quantum behavior (QPSO) to solve reactive power optimization in power system with distributed generation. Ionescu et al. [20] presented a mixed integer hybrid differential evolution method to handle multiperiod active power loss minimization. They formulated the task as a mixed-integer nonlinear programming (MINLP) problem, including constraints that specifically limit the number of switching actions between two successive anticipated system states. Also Rabiee et al. [21] presented a new voltage security constrained multiperiod ORPF approach which can determine the optimal settings for both continuous and discrete voltage control facilities in a coordinated manner. For a given horizon of time with definite time intervals, generalized Benders decomposition is adopted to convert RPO (reactive power optimization) into a simple mixed integer nonlinear programming and a relaxed nonlinear programming with fixed discrete variables.The traditional reactive power optimization model is based on certain power system conditions without considering the uncertainties. Additionally, a large number of renewable resources in modern power systems challenge the traditional planning and operations due to their stochastic features. Without considering the uncertainties introduced by the renewable energy especially the wind power, the performance of RPO will be significantly impacted. To solve the uncertainties of wind power output, a two-stage robust optimization model in [22] is presented to coordinate the discrete and continuous reactive power compensators and figure out a robust optimal solution that can hedge against any possible realization with the uncertainty of the wind power output. Generally, the research of the robust optimization model contains two-stage decisions. The first-stage decisions serve as the “here-and-now” decisions. And the second stage decisions are the “wait-and-see” decisions that can be adjusted after the first-stage decisions are determined and the wind power uncertainty is revealed. Rabih [23] presented a decentralized approach for controlling reactive power from a photovoltaic (PV) inverter through a linear decision rule that is in terms of the PV generated real power. The linear decision rules are computed using an affinely adjustable robust counterpart of an optimal power flow type formulation, with the PV real power generation specified in an uncertainty interval. In [24], an optimization scenario method is established to deal with the uncertainty of the output of the wind farm. The method decomposes the reactive power optimization into two phases, which are, respectively, solved by the interior point method and the genetic algorithm. The presented method can utilize the information of the probability distribution function of the output of the wind farm to effectively deal with the discrete variables.Motivated by the previous work, in order to solve the uncertainties of renewable resources in DRPO, this paper presents a multistage dynamic reactive power optimization method based on interval uncertainty (MISOCP), which belongs to the large-scale mixed integer nonlinear nonconvex programming. And then, MISOCP is decomposed into two determinist nonlinear programming subproblems using the theory of direct interval matching and selection of the extreme value intervals. Then, the heuristic search and the variable correction are employed to relax the nonconvex power flow equations of the determinist nonlinear programming subproblem to obtain the mixed integer quadratic programming model which can be finally solved.The rest of this paper is organized as follows: in Section2, a dynamic reactive power optimization model based on interval uncertainty is presented; Section 3 presents the solution of the model; Section 4 presents the solution of the MISOCP based relaxation formulation; Section 5 discusses the experimental results; Section 6 concludes the paper.


## 2. A Dynamic Reactive Power Optimization Model Based on Interval Uncertainty

### 2.1. Uncertain Variable Description

Renewable energy including photovoltaic (PV) and wind power is certainly affected by the light intensity and wind speed, which results in the strong uncertainty and uncontrollability of their output. Moreover the load forecasting in the power system also has the uncertainty of the prediction errors. In the RPO, the randomness is frequently employed to describe the uncertainty of the system, which requires a large number of historical data to determine the statistical law (probability characteristic parameters, probability distribution, etc.) of the system input uncertain variables. However, the data is difficult to achieve precisely from practical systems decisions. For solving the issue, interval number optimization is employed by this paper so that the upper and lower bounds of variables with less information can be achieved using the uncertainty of interval description variables. Based on the interval optimization method, the uncertainty of the wind power, PV output, and load forecasting can be represented using(1)PWT±t=PWT-t,PWT+tQWT±t=QWT-t,QWT+t(2)PPV±t=PPV-t,PPV+tQPV±t=QPV-t,QPV+t(3)PD±t=PD-t,PD+tQD±t=QD-t,QD+twhere PWT±t and QWT±t are the interval variables for the active and reactive wind power; PPV±t and QPV±t are the interval variables for the active and reactive PV power; PD±t and QD±t are the interval variables for the active and reactive load power; PWT-t, PWT+t, QWT-t, QWT+t, PPV-t, PPV+t, QPV-t, QPV+t, PD-t, PD+t, QD-t, and QD+t are the lower bound and upper bound of the wind active and reactive power, the PV active and reactive output, and the load active and reactive forecasting.


### 2.2. Dynamic Reactive Power Optimization Model Based on Interval Uncertainty

The MISOCP is an interval mixed-integer nonlinear optimization problem with multiple constraints. The objective is to minimize of the active power losses, which is given by(4)min∑t∈TPloss,tt=∑t∈T∑i,j∈Erij,tPij,t2+Qij,t2Ui,t2The MISOCP optimization is supposed to satisfy the equality constraint of the power flow equation, as shown by(5)Pg,i,t±+PWT,i,t±+PPV,i,t±-PD,i,t±=∑κ∈δiPiκ,k-∑j∈πiPij,k-rij,tPij,t2+Qij,t2Uj,t2,∀i∈B(6)Qg,i,t±+QWT,i,t±+QPV,i,t±-QD,i,t±=∑κ∈δiQiκ,k-∑j∈πiQij,k-xij,tPij,t2+Qij,t2Uj,t2+bs,i,tUi,t2,∀i∉Ω(7)Qg,i,t±+QWT,i,t±+QPV,i,t±-QD,i,t±+12Ui,t2Ci,t=∑κ∈δiQiκ,k-∑j∈πiQij,k-xij,tPij,t2+Qij,t2Uj,t2+bs,i,tUi,t2,∀i∈Ω(8)Ui,t2=Uj,t2-2rij,tPij,t+xij,tQij,t+rij,t2+xij,t2Pij,t2+Qij,t2Uj,t2,∀i,j∈ET(9)Ui,t2tsij,t2=Uj,t2-2rij,tPij,t+xij,tQij,t+rij,t2+xij,t2Pij,t2+Qij,t2Uj,t2,∀i,j∈EThe inequality constraints of the MISOCP optimization are given as shown by (10) to (16):Voltage and branch current constraints:(10)Ui,min,t≤Ui,t≤Ui,max,t,∀l∈B(11)Pij,t±2+Qij,t±2Vi,t2≤Iij,max,t,∀i,j∈ESwitched capacitor/reactors (SCRs) discrete constraints:(12)Ci,t=yi,tsi,t,∀i∈Ω(13)Ci,min,t≤Ci,t≤Ci,max,t,∀i∈Ω(14)∑t=0Tnyi,t+1-yi,t≤Yc,∀i∈ΩOn-load taps changers (OLTCs) discrete adjustment constraints:(15)tsij,min,t≤tsij,t≤tsij,max,t,∀i,j∈T(16)∑t=0Tntsij,t+1-tsij,t≤YTS,∀i,j∈Twhere t=0,1,⋯,Tn;Ploss,t represents the total transmission loss at time period t; B is the set of buses; E is the set of branches; T is the set of branches with transformers; Ω is the set of buses for shunt capacitors/reactors; πi represents the set of all parents of bus i; δi represents the set of all children of bus i; i,j∈E/T denotes i,j∈E, i,j∉E; Pij/Qij represents the active/reactive power flow from bus i to bus j; Pi/Qi represents the active/reactive of the generator at bus i; Ui represents the voltage magnitude of bus i; rij/xij represents the resistance/reactance of branch (i,j); bs,i represents the shunt susceptance from i to ground; Ui,max/Ui,min represent the upper/lower bound of voltage magnitude at bus i; Iij,max represents the current capacity limit of branch (i,j); Ci represents the value of shunt reactive power compensation of bus i; yi represents the optimal step of shunt capacitors/reactors at bus i; si represents the step size of shunt capacitors/reactors at bus j; Ci,max/Ci,min represents the upper/lower bound of reactive power compensation of bus i; YC are the specified maximum operational times of the SCRs capacitors; tsij,max/tsij,min represents the upper/lower bound of the tap ratio of the transformer branch (i,j); tsij represents the tap ratio of the transformer branch (i,j); and YTS are the specified maximum operational times of the tap changers.


## 2.1. Uncertain Variable Description

Renewable energy including photovoltaic (PV) and wind power is certainly affected by the light intensity and wind speed, which results in the strong uncertainty and uncontrollability of their output. Moreover the load forecasting in the power system also has the uncertainty of the prediction errors. In the RPO, the randomness is frequently employed to describe the uncertainty of the system, which requires a large number of historical data to determine the statistical law (probability characteristic parameters, probability distribution, etc.) of the system input uncertain variables. However, the data is difficult to achieve precisely from practical systems decisions. For solving the issue, interval number optimization is employed by this paper so that the upper and lower bounds of variables with less information can be achieved using the uncertainty of interval description variables. Based on the interval optimization method, the uncertainty of the wind power, PV output, and load forecasting can be represented using(1)PWT±t=PWT-t,PWT+tQWT±t=QWT-t,QWT+t(2)PPV±t=PPV-t,PPV+tQPV±t=QPV-t,QPV+t(3)PD±t=PD-t,PD+tQD±t=QD-t,QD+twhere PWT±t and QWT±t are the interval variables for the active and reactive wind power; PPV±t and QPV±t are the interval variables for the active and reactive PV power; PD±t and QD±t are the interval variables for the active and reactive load power; PWT-t, PWT+t, QWT-t, QWT+t, PPV-t, PPV+t, QPV-t, QPV+t, PD-t, PD+t, QD-t, and QD+t are the lower bound and upper bound of the wind active and reactive power, the PV active and reactive output, and the load active and reactive forecasting.


## 2.2. Dynamic Reactive Power Optimization Model Based on Interval Uncertainty

The MISOCP is an interval mixed-integer nonlinear optimization problem with multiple constraints. The objective is to minimize of the active power losses, which is given by(4)min∑t∈TPloss,tt=∑t∈T∑i,j∈Erij,tPij,t2+Qij,t2Ui,t2The MISOCP optimization is supposed to satisfy the equality constraint of the power flow equation, as shown by(5)Pg,i,t±+PWT,i,t±+PPV,i,t±-PD,i,t±=∑κ∈δiPiκ,k-∑j∈πiPij,k-rij,tPij,t2+Qij,t2Uj,t2,∀i∈B(6)Qg,i,t±+QWT,i,t±+QPV,i,t±-QD,i,t±=∑κ∈δiQiκ,k-∑j∈πiQij,k-xij,tPij,t2+Qij,t2Uj,t2+bs,i,tUi,t2,∀i∉Ω(7)Qg,i,t±+QWT,i,t±+QPV,i,t±-QD,i,t±+12Ui,t2Ci,t=∑κ∈δiQiκ,k-∑j∈πiQij,k-xij,tPij,t2+Qij,t2Uj,t2+bs,i,tUi,t2,∀i∈Ω(8)Ui,t2=Uj,t2-2rij,tPij,t+xij,tQij,t+rij,t2+xij,t2Pij,t2+Qij,t2Uj,t2,∀i,j∈ET(9)Ui,t2tsij,t2=Uj,t2-2rij,tPij,t+xij,tQij,t+rij,t2+xij,t2Pij,t2+Qij,t2Uj,t2,∀i,j∈EThe inequality constraints of the MISOCP optimization are given as shown by (10) to (16):Voltage and branch current constraints:(10)Ui,min,t≤Ui,t≤Ui,max,t,∀l∈B(11)Pij,t±2+Qij,t±2Vi,t2≤Iij,max,t,∀i,j∈ESwitched capacitor/reactors (SCRs) discrete constraints:(12)Ci,t=yi,tsi,t,∀i∈Ω(13)Ci,min,t≤Ci,t≤Ci,max,t,∀i∈Ω(14)∑t=0Tnyi,t+1-yi,t≤Yc,∀i∈ΩOn-load taps changers (OLTCs) discrete adjustment constraints:(15)tsij,min,t≤tsij,t≤tsij,max,t,∀i,j∈T(16)∑t=0Tntsij,t+1-tsij,t≤YTS,∀i,j∈Twhere t=0,1,⋯,Tn;Ploss,t represents the total transmission loss at time period t; B is the set of buses; E is the set of branches; T is the set of branches with transformers; Ω is the set of buses for shunt capacitors/reactors; πi represents the set of all parents of bus i; δi represents the set of all children of bus i; i,j∈E/T denotes i,j∈E, i,j∉E; Pij/Qij represents the active/reactive power flow from bus i to bus j; Pi/Qi represents the active/reactive of the generator at bus i; Ui represents the voltage magnitude of bus i; rij/xij represents the resistance/reactance of branch (i,j); bs,i represents the shunt susceptance from i to ground; Ui,max/Ui,min represent the upper/lower bound of voltage magnitude at bus i; Iij,max represents the current capacity limit of branch (i,j); Ci represents the value of shunt reactive power compensation of bus i; yi represents the optimal step of shunt capacitors/reactors at bus i; si represents the step size of shunt capacitors/reactors at bus j; Ci,max/Ci,min represents the upper/lower bound of reactive power compensation of bus i; YC are the specified maximum operational times of the SCRs capacitors; tsij,max/tsij,min represents the upper/lower bound of the tap ratio of the transformer branch (i,j); tsij represents the tap ratio of the transformer branch (i,j); and YTS are the specified maximum operational times of the tap changers.


## 3. Optimization Model Conversion Based on Interval Uncertainty

The interval number is defined as a set of random variables with the upper and the lower bounds:(17)x±=x-,x+=x∣x-≤x≤x+,x∈Rwhere x- and x+ are the upper and the lower bounds. When x-=x+, the interval degenerates to a real number.The general form of an interval optimization problem is(18)minf±x±s.t.g±x±≥0x±=x-,x+where x± is the vector of the variables which can be represented by(19)x±=U1±,⋯UNB±,PGslack±,Q1±⋯QNG±,C1±,⋯CNC±,ts1±,⋯,tsNT±U1±,⋯UNB± and PGslack± are the state variables; Q1±⋯QNG±, C1±,⋯CNC±, and ts1±,⋯,tsNT± are the control variables.According to the basic properties of the interval matching and the extreme interval selection, the following theorem can be obtained [25]: The necessary and sufficient condition for the interval function f±x±=f+x,f-x, to obtain the maximum (minimum) value at the point x∗ in the field G, is that the upper boundary functions f+x and lower boundary functions f-x obtain the maximum (minimum) value at this point. As a result, the problem of the interval planning can be converted into two regular deterministic optimization problems.Interval mixed-integer nonlinear optimization problem represented by (4) to (16) can be transformed into the deterministic nonlinear programming problems with the minimization of the upper and the lower boundary functions denoted by(20)minf-xs.t.g+x≤0(21)minf+xs.t.g+x≤0


## 4. The Multistage Mathematical Model of Dynamic Reactive Power Optimization

### 4.1. Relaxed Dynamic Reactive Power Optimization Model for Interior Points

Equations (20) and (21) are both deterministic nonlinear programming problems. Firstly, the discrete control variables such as switched capacitor/reactors and on-load taps changers are relaxed as continuous variables. And then the restrictions on the number of times of reactive power control equipment are relaxed, which means we do not consider the limit of the times of actions firstly. The optimal solution obtained by solving the relaxation model is the lower bound of the optimal solution of the original reactive power optimization model shown by (20) and (21). Relaxation model is as follows:(22)min∑t=0Tnfxt,x~t(23)s.t.htxt,x~t=0(24)gtxt,x~t≤0(25)xmin,t≤xt≤xmax,t(26)xmin,t≤x~t≤xmax,t(27)∑t=0Tn-1x~i,t+1-x~i,t≤Y,i=1,⋯,nEquation (27) contains the constraint of the absolute value of the variables, which does not satisfy the second-order continuous and differentiable properties. To overcome this problem, we present a method to transform absolute inequality constraints into general inequality constraints using(28)-δi,t≤x~i,t+1-x~i,t≤δi,t,i=1,⋯,n;t=0,⋯,Tn-10≤δi,t≤δmax,i,t,i=1,⋯,n;t=0,⋯,Tn-1∑t=0Tn-1δi,t≤Yi,i=1,⋯,nwhere δi,t is the auxiliary variable and δmax,i,t are the specified maximum operational times of the SCRs capacitors and the on-load tap changers at time t.


### 4.2. Loss in Increments as Small as the Target of Discrete Variables

The solution of the relaxation model isf-- and f-+, which are the lower bounds of the original problems (20) and (21). Based on the solution of the relaxation model, a model is established for each discrete device that minimizes the amount of collation error as the objective function with constraints of time coupling. The variables discretization model of problems (20) and (21) is essentially a mixed integer quadratic programming:(29)minf-x-f--2s.t.∑t=0Tn-1yi,t+1-yi,t≤Yyi∈Integer(30)minf+x-f-+2s.t.∑t=0Tn-1yi,t+1-yi,t≤Yyi∈Integer


### 4.3. Fine-Tuning for Continuous Variables

By solving (29) and (30), the optimal solution of all the discrete variables in the original optimization model represented by (20) and (21) can be achieved. And then, all the discrete variables in (20) and (21) can be fixed by an optimization solution obtained by (29) and (30), so that the original optimization model is transformed into a general nonlinear optimization model containing only continuous variables which can be solved quickly by interior point method.


## 4.1. Relaxed Dynamic Reactive Power Optimization Model for Interior Points

Equations (20) and (21) are both deterministic nonlinear programming problems. Firstly, the discrete control variables such as switched capacitor/reactors and on-load taps changers are relaxed as continuous variables. And then the restrictions on the number of times of reactive power control equipment are relaxed, which means we do not consider the limit of the times of actions firstly. The optimal solution obtained by solving the relaxation model is the lower bound of the optimal solution of the original reactive power optimization model shown by (20) and (21). Relaxation model is as follows:(22)min∑t=0Tnfxt,x~t(23)s.t.htxt,x~t=0(24)gtxt,x~t≤0(25)xmin,t≤xt≤xmax,t(26)xmin,t≤x~t≤xmax,t(27)∑t=0Tn-1x~i,t+1-x~i,t≤Y,i=1,⋯,nEquation (27) contains the constraint of the absolute value of the variables, which does not satisfy the second-order continuous and differentiable properties. To overcome this problem, we present a method to transform absolute inequality constraints into general inequality constraints using(28)-δi,t≤x~i,t+1-x~i,t≤δi,t,i=1,⋯,n;t=0,⋯,Tn-10≤δi,t≤δmax,i,t,i=1,⋯,n;t=0,⋯,Tn-1∑t=0Tn-1δi,t≤Yi,i=1,⋯,nwhere δi,t is the auxiliary variable and δmax,i,t are the specified maximum operational times of the SCRs capacitors and the on-load tap changers at time t.


## 4.2. Loss in Increments as Small as the Target of Discrete Variables

The solution of the relaxation model isf-- and f-+, which are the lower bounds of the original problems (20) and (21). Based on the solution of the relaxation model, a model is established for each discrete device that minimizes the amount of collation error as the objective function with constraints of time coupling. The variables discretization model of problems (20) and (21) is essentially a mixed integer quadratic programming:(29)minf-x-f--2s.t.∑t=0Tn-1yi,t+1-yi,t≤Yyi∈Integer(30)minf+x-f-+2s.t.∑t=0Tn-1yi,t+1-yi,t≤Yyi∈Integer


## 4.3. Fine-Tuning for Continuous Variables

By solving (29) and (30), the optimal solution of all the discrete variables in the original optimization model represented by (20) and (21) can be achieved. And then, all the discrete variables in (20) and (21) can be fixed by an optimization solution obtained by (29) and (30), so that the original optimization model is transformed into a general nonlinear optimization model containing only continuous variables which can be solved quickly by interior point method.


## 5. Experimental Result

### 5.1. Configuration of Three Test Systems

Three benchmark systems including IEEE 14-bus, 118-bus, and 300-bus system are employed to evaluate the presented approaches facilitating DPRO containing the renewable energy resources. The expanded IEEE 14-bus system has three transformer branches, eight load nodes, and five generators. The load data is calculated by the load of the IEEE 14-bus load multiplied by scale factor which is depicted in Figure1. The bounds of the loads and the generator output are set as [0.95,1.05] p.u.. The wind farm is connected to node 5, which has a total number of 40 wind turbines. Single rated power of wind turbine is 2.5MW. The total wind power output of the wind farm is within [20MW,80MW]. Photovoltaic power stations with a rated power of 50MW are connected to node 14, and the light intensity interval is given in [26]. The PV array output is within [29.9kW,36.6kW]. The parameters of the reactive power adjustment equipment in the IEEE 14-bus system are shown in Table 1.Table 1
Parameters of reactive power adjustment equipment (IEEE 14-bus).
Capacitor node upper bound/Mvar lower bound/Mvar Adjustment step /Mvar Number of actions 1 7 +30 - 30 2 5 2 13 +30 - 30 3 4 3 14 +30 - 30 3 4  Transformer Branch Lower bound/p.u. Upper bound/p.u. Adjustment step/p.u. Number of actions  1 4-7 0.925 1.075 0.025 4  2 4-9 0.925 1.075 0.025 5  3 5-6 0.925 1.075 0.025 5Figure 1
Typical daily load curve.The location, the active output, and the reactive power output of the renewable energy power generation in the IEEE 118 and 300 systems are shown in Table2. The parameters of the reactive power adjustment equipment in the IEEE 118-bus and 300-bus systems are shown in Tables 3 and 4.Table 2
Parameter of renewable energy sources (IEEE 118-bus and IEEE 300-bus).
Number Connected node Active power output interval/MW Reactive power output interval/Mvar IEEE 118 IEEE 300 1 6 8 [ 8.2,17.8 ] [ - 2.6,4.4 ] 2 11 20 [ 8.9,15.1 ] [ - 3.1,5.2 ] 3 19 28 [ 10.4,24.5 ] [ - 4.7,7.3 ] 4 24 99 [ 2.6,15.4 ] [ - 3.8,6.9 ] 5 39 147 [ 12.3,18.9 ] [ - 2.8,5.5 ] 6 69 176 [ 3.1,19.9 ] [ - 3.2,6.3 ] 7 85 213 [ 11.2,19.2 ] [ - 2.7,4.7 ] 8 95 222 [ 6.2,23.2 ] [ - 3.5,5.7 ] 9 107 236 [ 8.2,17.5 ] [ - 2.7,6.8 ] 10 111 239 [ 7.0,24.5 ] [ - 3.9,7.5 ]Table 3
Parameters of reactive power adjustment equipment (IEEE 118-bus).
(a) Capacitor node upper bound/Mvar lower bound/Mvar Adjustment step/Mvar Number of actions 1 2 +30 - 30 2 5 2 16 +30 - 30 3 3 3 33 +30 - 30 3 4 4 41 +30 - 30 2 3 5 43 +30 - 30 4 3 6 44 +30 - 30 6 4 7 48 +30 - 30 3 5 8 53 +30 - 30 4 4 9 67 +30 - 30 3 2 10 82 +30 - 30 5 5 11 88 +30 - 30 3 4 12 96 +30 - 30 5 3  (b) Transformer Branch upper bound/p.u. Upper bound/p.u. Adjustment step/p.u. Number of actions 1 8 - 5 0.925 1.075 0.025 4 2 26 - 25 0.925 1.075 0.025 3 3 37 - 38 0.925 1.075 0.025 4 4 63-59 0.925 1.075 0.025 4 5 65-66 0.925 1.075 0.025 3Table 4
Parameters of reactive power adjustment equipment (IEEE 300-bus).
(a) Capacitor node upper bound/Mvar lower bound/Mvar Adjustment step/Mvar Number of actions 1 7 +30 - 30 2 5 2 16 +30 - 30 3 3 3 35 +30 - 30 3 4 4 44 +30 - 30 2 3 5 48 +30 - 30 4 3 6 51 +30 - 30 6 4 7 58 +30 - 30 3 5 8 63 +30 - 30 4 4 9 67 +30 - 30 3 2 10 82 +30 - 30 5 5 11 88 +30 - 30 3 4 12 97 +30 - 30 5 3 13 136 +30 - 30 4 5 14 159 +30 - 30 3 3 15 164 +30 - 30 5 4 16 173 +30 - 30 3 4 17 278 +30 - 30 4 3 18 289 +30 - 30 6 4 19 7008 +30 - 30 4 2 20 7015 +30 - 30 4 3 21 7139 +30 - 30 2 4 22 9044 +30 - 30 5 3  (b) Transformer Branch upper bound/p.u. Upper bound/p.u. Adjustment step/p.u. Number of actions 1 9001-9006 0.925 1.075 0.025 4 2 9001-9012 0.925 1.075 0.025 3 3 9005-9051 0.925 1.075 0.025 4 4 9005-9052 0.925 1.075 0.025 4 5 9005-9055 0.925 1.075 0.025 3 6 3-1 0.925 1.075 0.025 4 7 3-2 0.925 1.075 0.025 3 8 7-5 0.925 1.075 0.025 3 9 7-6 0.925 1.075 0.025 4 10 10-11 0.925 1.075 0.025 3 11 12-11 0.925 1.075 0.025 3 12 21-20 0.925 1.075 0.025 3 13 24-23 0.925 1.075 0.025 4 14 36-35 0.925 1.075 0.025 4 15 160-121 0.925 1.075 0.025 3


### 5.2. Comparison of Solution Property between Interval Optimization and Deterministic Approach

In order to verify the effectiveness of the interval optimization algorithm in this paper, the optimization based on Monte Carlo (MC) is also conducted in terms of comparison. Assuming that the known interval variables are the uniformly distributed uncertain variables, the interval uncertainty in the dynamic multistage reactive power optimization is carried out using MC. The maximum and minimum values of the interval uncertainty in the calculation result are employed as the upper and lower bounds for the desired interval variable. Figure2 shows the IEEE 14-bus system network loss boundary along with the increasing number of the MC samples. It can be observed that the number of the iterations increases, and the upper and lower bounds of the network loss of the system tend to be stable. In our experiments the number of samples is set to be 5000 to make a compromise between calculation time and algorithm accuracy. The results of the MC method are used as a comparative standard.Figure 2
Number of sample points in MC method effect on the result of MISOCP.In order to evaluate the calculation error of the method in this paper, we define two interval indices including the interval meanEinter and the interval radius Rinter. For any interval [A-,A+], the mean and radius of the interval are defined as(31)Einter=A-+A+2(32)Rinter=A+-A-2The relative error between the MISOCP and MC method is defined as(33)εE=EinterMISOCP-EinterMCEinterMC×100%(34)εR=RinterMISOCP-RinterMCRinterMC×100%where EinterMISOCP and EinterMC are the interval mean of the total transmission losses obtained using the MISOCP and MC method and RinterMISOCP and RinterMC are the interval radii of the total transmission losses. The results of the interval radius and mean of each system are shown in Tables 5 and 6.Table 5
Relative error analysis of the interval radius of transmission losses in MISOCP.
Relative error of the interval radius of transmission losses/ % IEEE 14 IEEE 118 IEEE 300 0.64 7.85 7.16Table 6
Relative error analysis of the average of transmission losses in MC.
relative error of the interval radius of transmission losses/ % IEEE 14 IEEE 118 IEEE 300 0.00 0.04 0.05As can be seen from Tables5 and 6, the interval radius errors of the transmission loss do not exceed 8%; the interval mean error does not exceed 0.6%. Obviously, the result of the MISOCP method is highly close to that of the MC method. Based on the analysis mentioned in Tables 5 and 6, it can be seen that the results of MISOCP method and MC method are relatively close. However, the results of the MISOCP method are slightly conservative, which means that the overestimated phenomenon exists. The reasons are as follows: (1) The MC method tends to underestimate the bounds of interval variables in the case of interval variables with a mean of 0, which is presented in research [25]; (2) the MISOCP method does not deal with the correlation and the package effect between the interval variables, which causes its conservation and overestimation for the interval boundary.In this paper, we employ the interval matching and extreme interval selection to decompose the DRPO model considering renewable energy integration into two deterministic dynamic reactive power optimization models. And then our approach relaxes the discrete variables into the continuous variables and obtains the relaxation solution. The total system transmission losses without considering the renewable energy integration are 30.674 MW. The lower bound of the transmission loss is 32.346MW and the upper bound of the transmission loss is 34.672MW. Then, we can get the transmission loss interval results of[32.346MW,34.672MW]. Using the approaches presented in this paper, we can get the whole system node voltage amplitude distribution as shown in Table 7 which indicates that when the renewable energy is connected, the node voltages are basically within the voltage range using the system dynamic reactive power optimization.Table 7
Node interval voltage after optimization.
Number Node interval voltage after optimization Node 1 [ 1.076,1.076 ] Node 2 [ 1.041,1.041 ] Node 3 [ 1.032,1.032 ] Node 4 [ 1.015,1.043 ] Node 5 [ 1.013,1.046 ] Node 6 [ 1.046,1.046 ] Node 7 [ 1.024,1.044 ] Node 8 [ 1.072,1.072 ] Node 9 [ 1.037,1.049 ] Node 10 [ 1.019,1.036 ] Node 11 [ 1.043,1.056 ] Node 12 [ 1.031,1.051 ] Node 13 [ 1.026   1.047 ] Node 14 [ 1.016,1.038 ]We further analyze the impacts of the discrete device actions and the device adjustment steps of the optimization results. Obviously, the result of the discretization is closer to the relaxation solution, and the result is better. Based on this idea, the discretization results of the actual operation times of SCRs and OLTCs are shown in Figures3 and 4.Figure 3
Discretization of capacitor switching.Figure 4
Discretization of transformer tap adjustment.Figures3 and 4 indicate that when the number of the operations is constant, the large-step adjustment can achieve better results. When the number of the actions is small, the adjustment of the small steps is difficult to follow the system load fluctuations so that the adjustment result is less effective.


### 5.3. Comparison of Computational Performance between Interval Optimization and Deterministic Approach

Table8 shows the running time of the MC and MISOCP in calculating different systems. The MC method requires a large number of sampling points (5000 in this paper) for power flow calculation. However, the calculation of MISOCP method only needs to calculate the nonlinear interior programming problem twice, which performs efficiently. Table 8 also indicates that the efficiency of MISOCP method is more efficient than that of the MC method when the system size increases.Table 8
Comparison of calculation time of three test systems.
Time(s) IEEE 14 IEEE 118 IEEE 300 MC 213.7 2685.8 7231.4 MISOCP 5.45 21.38 43.73


## 5.1. Configuration of Three Test Systems

Three benchmark systems including IEEE 14-bus, 118-bus, and 300-bus system are employed to evaluate the presented approaches facilitating DPRO containing the renewable energy resources. The expanded IEEE 14-bus system has three transformer branches, eight load nodes, and five generators. The load data is calculated by the load of the IEEE 14-bus load multiplied by scale factor which is depicted in Figure1. The bounds of the loads and the generator output are set as [0.95,1.05] p.u.. The wind farm is connected to node 5, which has a total number of 40 wind turbines. Single rated power of wind turbine is 2.5MW. The total wind power output of the wind farm is within [20MW,80MW]. Photovoltaic power stations with a rated power of 50MW are connected to node 14, and the light intensity interval is given in [26]. The PV array output is within [29.9kW,36.6kW]. The parameters of the reactive power adjustment equipment in the IEEE 14-bus system are shown in Table 1.Table 1
Parameters of reactive power adjustment equipment (IEEE 14-bus).
Capacitor node upper bound/Mvar lower bound/Mvar Adjustment step /Mvar Number of actions 1 7 +30 - 30 2 5 2 13 +30 - 30 3 4 3 14 +30 - 30 3 4  Transformer Branch Lower bound/p.u. Upper bound/p.u. Adjustment step/p.u. Number of actions  1 4-7 0.925 1.075 0.025 4  2 4-9 0.925 1.075 0.025 5  3 5-6 0.925 1.075 0.025 5Figure 1
Typical daily load curve.The location, the active output, and the reactive power output of the renewable energy power generation in the IEEE 118 and 300 systems are shown in Table2. The parameters of the reactive power adjustment equipment in the IEEE 118-bus and 300-bus systems are shown in Tables 3 and 4.Table 2
Parameter of renewable energy sources (IEEE 118-bus and IEEE 300-bus).
Number Connected node Active power output interval/MW Reactive power output interval/Mvar IEEE 118 IEEE 300 1 6 8 [ 8.2,17.8 ] [ - 2.6,4.4 ] 2 11 20 [ 8.9,15.1 ] [ - 3.1,5.2 ] 3 19 28 [ 10.4,24.5 ] [ - 4.7,7.3 ] 4 24 99 [ 2.6,15.4 ] [ - 3.8,6.9 ] 5 39 147 [ 12.3,18.9 ] [ - 2.8,5.5 ] 6 69 176 [ 3.1,19.9 ] [ - 3.2,6.3 ] 7 85 213 [ 11.2,19.2 ] [ - 2.7,4.7 ] 8 95 222 [ 6.2,23.2 ] [ - 3.5,5.7 ] 9 107 236 [ 8.2,17.5 ] [ - 2.7,6.8 ] 10 111 239 [ 7.0,24.5 ] [ - 3.9,7.5 ]Table 3
Parameters of reactive power adjustment equipment (IEEE 118-bus).
(a) Capacitor node upper bound/Mvar lower bound/Mvar Adjustment step/Mvar Number of actions 1 2 +30 - 30 2 5 2 16 +30 - 30 3 3 3 33 +30 - 30 3 4 4 41 +30 - 30 2 3 5 43 +30 - 30 4 3 6 44 +30 - 30 6 4 7 48 +30 - 30 3 5 8 53 +30 - 30 4 4 9 67 +30 - 30 3 2 10 82 +30 - 30 5 5 11 88 +30 - 30 3 4 12 96 +30 - 30 5 3  (b) Transformer Branch upper bound/p.u. Upper bound/p.u. Adjustment step/p.u. Number of actions 1 8 - 5 0.925 1.075 0.025 4 2 26 - 25 0.925 1.075 0.025 3 3 37 - 38 0.925 1.075 0.025 4 4 63-59 0.925 1.075 0.025 4 5 65-66 0.925 1.075 0.025 3Table 4
Parameters of reactive power adjustment equipment (IEEE 300-bus).
(a) Capacitor node upper bound/Mvar lower bound/Mvar Adjustment step/Mvar Number of actions 1 7 +30 - 30 2 5 2 16 +30 - 30 3 3 3 35 +30 - 30 3 4 4 44 +30 - 30 2 3 5 48 +30 - 30 4 3 6 51 +30 - 30 6 4 7 58 +30 - 30 3 5 8 63 +30 - 30 4 4 9 67 +30 - 30 3 2 10 82 +30 - 30 5 5 11 88 +30 - 30 3 4 12 97 +30 - 30 5 3 13 136 +30 - 30 4 5 14 159 +30 - 30 3 3 15 164 +30 - 30 5 4 16 173 +30 - 30 3 4 17 278 +30 - 30 4 3 18 289 +30 - 30 6 4 19 7008 +30 - 30 4 2 20 7015 +30 - 30 4 3 21 7139 +30 - 30 2 4 22 9044 +30 - 30 5 3  (b) Transformer Branch upper bound/p.u. Upper bound/p.u. Adjustment step/p.u. Number of actions 1 9001-9006 0.925 1.075 0.025 4 2 9001-9012 0.925 1.075 0.025 3 3 9005-9051 0.925 1.075 0.025 4 4 9005-9052 0.925 1.075 0.025 4 5 9005-9055 0.925 1.075 0.025 3 6 3-1 0.925 1.075 0.025 4 7 3-2 0.925 1.075 0.025 3 8 7-5 0.925 1.075 0.025 3 9 7-6 0.925 1.075 0.025 4 10 10-11 0.925 1.075 0.025 3 11 12-11 0.925 1.075 0.025 3 12 21-20 0.925 1.075 0.025 3 13 24-23 0.925 1.075 0.025 4 14 36-35 0.925 1.075 0.025 4 15 160-121 0.925 1.075 0.025 3


## 5.2. Comparison of Solution Property between Interval Optimization and Deterministic Approach

In order to verify the effectiveness of the interval optimization algorithm in this paper, the optimization based on Monte Carlo (MC) is also conducted in terms of comparison. Assuming that the known interval variables are the uniformly distributed uncertain variables, the interval uncertainty in the dynamic multistage reactive power optimization is carried out using MC. The maximum and minimum values of the interval uncertainty in the calculation result are employed as the upper and lower bounds for the desired interval variable. Figure2 shows the IEEE 14-bus system network loss boundary along with the increasing number of the MC samples. It can be observed that the number of the iterations increases, and the upper and lower bounds of the network loss of the system tend to be stable. In our experiments the number of samples is set to be 5000 to make a compromise between calculation time and algorithm accuracy. The results of the MC method are used as a comparative standard.Figure 2
Number of sample points in MC method effect on the result of MISOCP.In order to evaluate the calculation error of the method in this paper, we define two interval indices including the interval meanEinter and the interval radius Rinter. For any interval [A-,A+], the mean and radius of the interval are defined as(31)Einter=A-+A+2(32)Rinter=A+-A-2The relative error between the MISOCP and MC method is defined as(33)εE=EinterMISOCP-EinterMCEinterMC×100%(34)εR=RinterMISOCP-RinterMCRinterMC×100%where EinterMISOCP and EinterMC are the interval mean of the total transmission losses obtained using the MISOCP and MC method and RinterMISOCP and RinterMC are the interval radii of the total transmission losses. The results of the interval radius and mean of each system are shown in Tables 5 and 6.Table 5
Relative error analysis of the interval radius of transmission losses in MISOCP.
Relative error of the interval radius of transmission losses/ % IEEE 14 IEEE 118 IEEE 300 0.64 7.85 7.16Table 6
Relative error analysis of the average of transmission losses in MC.
relative error of the interval radius of transmission losses/ % IEEE 14 IEEE 118 IEEE 300 0.00 0.04 0.05As can be seen from Tables5 and 6, the interval radius errors of the transmission loss do not exceed 8%; the interval mean error does not exceed 0.6%. Obviously, the result of the MISOCP method is highly close to that of the MC method. Based on the analysis mentioned in Tables 5 and 6, it can be seen that the results of MISOCP method and MC method are relatively close. However, the results of the MISOCP method are slightly conservative, which means that the overestimated phenomenon exists. The reasons are as follows: (1) The MC method tends to underestimate the bounds of interval variables in the case of interval variables with a mean of 0, which is presented in research [25]; (2) the MISOCP method does not deal with the correlation and the package effect between the interval variables, which causes its conservation and overestimation for the interval boundary.In this paper, we employ the interval matching and extreme interval selection to decompose the DRPO model considering renewable energy integration into two deterministic dynamic reactive power optimization models. And then our approach relaxes the discrete variables into the continuous variables and obtains the relaxation solution. The total system transmission losses without considering the renewable energy integration are 30.674 MW. The lower bound of the transmission loss is 32.346MW and the upper bound of the transmission loss is 34.672MW. Then, we can get the transmission loss interval results of[32.346MW,34.672MW]. Using the approaches presented in this paper, we can get the whole system node voltage amplitude distribution as shown in Table 7 which indicates that when the renewable energy is connected, the node voltages are basically within the voltage range using the system dynamic reactive power optimization.Table 7
Node interval voltage after optimization.
Number Node interval voltage after optimization Node 1 [ 1.076,1.076 ] Node 2 [ 1.041,1.041 ] Node 3 [ 1.032,1.032 ] Node 4 [ 1.015,1.043 ] Node 5 [ 1.013,1.046 ] Node 6 [ 1.046,1.046 ] Node 7 [ 1.024,1.044 ] Node 8 [ 1.072,1.072 ] Node 9 [ 1.037,1.049 ] Node 10 [ 1.019,1.036 ] Node 11 [ 1.043,1.056 ] Node 12 [ 1.031,1.051 ] Node 13 [ 1.026   1.047 ] Node 14 [ 1.016,1.038 ]We further analyze the impacts of the discrete device actions and the device adjustment steps of the optimization results. Obviously, the result of the discretization is closer to the relaxation solution, and the result is better. Based on this idea, the discretization results of the actual operation times of SCRs and OLTCs are shown in Figures3 and 4.Figure 3
Discretization of capacitor switching.Figure 4
Discretization of transformer tap adjustment.Figures3 and 4 indicate that when the number of the operations is constant, the large-step adjustment can achieve better results. When the number of the actions is small, the adjustment of the small steps is difficult to follow the system load fluctuations so that the adjustment result is less effective.


## 5.3. Comparison of Computational Performance between Interval Optimization and Deterministic Approach

Table8 shows the running time of the MC and MISOCP in calculating different systems. The MC method requires a large number of sampling points (5000 in this paper) for power flow calculation. However, the calculation of MISOCP method only needs to calculate the nonlinear interior programming problem twice, which performs efficiently. Table 8 also indicates that the efficiency of MISOCP method is more efficient than that of the MC method when the system size increases.Table 8
Comparison of calculation time of three test systems.
Time(s) IEEE 14 IEEE 118 IEEE 300 MC 213.7 2685.8 7231.4 MISOCP 5.45 21.38 43.73


## 6. Conclusion

In this paper, a multistage dynamic reactive power optimization method based on the interval uncertainty is presented to address the high penetration of uncertain renewable energy integrated into power system networks, which uses the selection theory of interval extremum to decompose the uncertainty problem into two deterministic same-type subproblems. To address the daily operating times of discrete control variables, intertemporal constraints and discrete control variables in the DRPO are relaxed to a large-scale nonlinear nonconvex programming. The experimental results show that the proposed method can solve the optimization with higher efficiency and less error.



---

*Source: 3854812-2018-10-18.xml*