Research on Stochastic Optimal Operation Strategy of Active Distribution Network Considering Intermittent Energy

Abstract

Active distribution networks characterized by high flexibility and controllability are an important development mode of future smart grids to be interconnected with large scale distributed generation sources including intermittent energies. However, the uncertainty of intermittent energy and the diversity of controllable devices make the optimal operation of distribution network a challenging issue. In this paper, we propose a stochastic optimal operation strategy for distribution networks with the objective function considering the operation state of the distribution network. Both distributed generations and flexible loads are taken into consideration in our strategy. The uncertainty of the intermittent energy is considered in this paper to obtain an optimized operation and an efficient utilization of intermittent energy under the worst scenario. Then, Benders decomposition is used in this paper to solve the two-stage max-min problem for stochastic optimal operation. Finally, we test the effectiveness of our strategy under different scenarios of the demonstration project of active distribution network located in Guizhou, China.

1. Introduction

In recent years, the access of intermittent energy, such as wind and solar, keeps increasing in distribution networks. Under such a circumstance, active distribution networks (ADNs) characterized by high distributed generation penetration and high control requirements have been widely accepted as an important development mode of future smart grids [1,2]. Compared with traditional distribution networks, active distribution networks have more controllable equipment, including energy storage systems, flexible loads, and other distributed generation sources. In order to ensure the optimal operation of active distribution networks, reasonable arrangement for the output of different types of controllable equipment is essential [3]. However, in active distribution network, the uncertainty of intermittent energy makes the traditional deterministic power system to a random-deterministic coupling power system. At the same time, distribution network operation concerns more about the economic benefit and the operating state, which brings great challenge to the power system optimization.

To realize optimal operation of an ADN considering the uncertainty of intermittent energy, two key issues need to be taken into consideration: the uncertainty of intermittent energy and the coordination between the risk and the economic benefit of the ADN.

Intermittent energy power forecasting is usually introduced to optimal operation to achieve a more precisely control. There are a variety of models and prediction methods such as artificial neural networks (ANN), Kalman filtering method and fuzzy logic method used in power forecasting to obtain accurate predictions [4]. However, the forecasting for intermittent energy influenced by complex factors related to the environment is still a challenging problem. Researchers have done a lot of works on optimal operation considering the uncertainty of intermittent energy including stochastic unit commitment (SUC) [5,6,7] or stabilizing the unpredicted power fluctuation by energy storage system and demand response [8]. Such as both [6,7], the uncertainty of wind power or load is divided into several scenarios which represent different wind power output probability. It is obvious that the accuracy of the solution mainly depends on the selection of uncertainty scenarios.

To deal with the coordination between the risk and the economic benefits in power system, most researchers focus on minimizing the operating cost [9,10]. Like [9], it takes the minimization of total UC schedule cost as main consideration, while considering the alleviation of wind power curtailment and the risk of inclusion of wind energy generation at the same time. Due to the different requirements of the system operators, some other objectives are also used such as minimizing the variance for the risk-averse case or even multi-objectives [11,12,13]. Unlike the transmission network concerns more about the operating cost, one of the priorities of ADN is to improve the economic benefits brought by different controllable devices. At the same time, the operational characteristics of ADN like the peak-valley difference and the network loss are also concerned. The reduction of peak-valley difference and network loss can reduce the reserve capacity of the transmission network and improve the efficiency of energy usage, which will also improve the economic efficiency of ADN. So, the objective function presented in this paper will take the peak-valley difference, the network loss and the economic benefits brought by different controllable devices into consideration. In order to avoid the risk of optimization due to intermittent energy prediction errors, the uncertainties of intermittent energies are also taken into account to ensure optimal operation of ADN under the worst scenario.

Due to the complexity of both controllable equipment and objective function, the optimal operation problem formulated as a mixed-integer linear programming (MILP) is non-deterministic polynomial-time hard (NP-hard). Several algorithms such as the hybrid genetic algorithm [14], differential evolution algorithm [15], discrete particle swarm optimization [11,16] and artificial bee colony [17] are used to solve the problem. However, the intelligent algorithms mentioned above have a common problem. They may not be able to obtain the optimal solution but a better one instead which will lead to a considerable uncertainty of the optimal result. Outer approximation methods are used in [18,19] to solve the mixed-integer nonlinear programming (MINLP) problem of unit commitment. Some researchers also have applied other methods such as Lagrangian relaxation (LR) and Benders decomposition (BD) to solve the problem [20,21]. This type of algorithm requires certain characteristics for constraints or objective functions, such as convex functions, etc. So it cannot be applied to solve all problems. But it is an effective way to simplify the problem-solving process.

In this paper, we mainly aimed to integrate the controllable resources in ADN to achieve a more efficient power system operation. As discussed before, the objective function for optimization will take the economic benefits of DGs, two key characteristics of ADN and the uncertainty of intermittent energies into consideration. While taking the constraint of power flow and so on into consideration, the optimization problem will form a complicated non-linear problem. To simplify the nonlinear programming problem proposed, we reformulated the objective function into a two-stage max-min problem. Benders decomposition is introduced to the approach to simplify the complexity of sub-problem. The algorithm can limit the uncertainty of intermittent energy output according to the requirement and find out the optimal controllable resource output.

The remaining part of this paper is organized as follows: Section 2 provides mathematical modeling of an optimal operation strategy. Section 3 introduces the solution approach based on Benders decomposition. The case study and optimization results are discussed in Section 4, combined with the active distribution network demonstration project located in Guizhou, China. Section 5 concludes with the main results of this paper.

2. Mathematical Formulation

In this section, we establish the optimization formulation for optimal operation of active distribution network considering multi-controllable resources. Combined with the conditions of Guizhou demonstration project, direct load control for flexible loads are also taken into account in this formulation as an adjustable resource for distribution network. In order to simplify the complexity of flexible load modeling, we will divide flexible load including Heating, Ventilation and Air Conditioning (HVAC) devices and electric vehicles into two categories, in accordance with their characteristics and unified it with the distributed generation model. This can simplify the optimization formulation.

In general, there are three kinds of distributed generations in a distribution network: power regulation generations, energy storage systems, and intermittent energy. Different types of distributed generations correspond to different constraints. In active distribution network, intermittent energy is generally desirable to achieve the maximum power consumption, so in this paper we only consider the uncertainty of intermittent energy. It won’t be taken as one of the controllable resources.

For power regulation generations such as thermal power generations and hydropower generations, the constrains mainly include the ramp-up rate, ramp-down rate, maximum power limit and minimum power limit in the form of (1)–(3):

$${\chi }_i^{t-1}-{\chi }_i^t\le {\alpha }_{down}$$

(1)

$${\chi }_i^t-{\chi }_i^{t-1}\le {\alpha }_{up}$$

(2)

$$L_{low,i}\le {\chi }_i^t\le L_{up,i}$$

(3)

In this paper, we lump flexible loads such as HVAC together with former power regulation generations as ‘power regulation devices’. For flexible load such as HVAC, the instantaneous power will keep changing when it turns on or off to maintain the room temperature. However, we only concern about the average power when taking optimization instead of instantaneous power. The average power for a HVAC when the set-temperature is ${\nu }_{set,i}(t)$ can be expressed as (4) [22]:

$$P_{eq,i}({\nu }_{set,i},t)=({\nu }_{set,i}(t)-{\nu }_{outside}(t))/{\eta }_i{R}_{eq,i}$$

(4)

$R_{eq,i}$ is the equivalent thermal resistance for $TC{L}_i$ (degrees Celsius/W). ${\nu }_{outside}(t)$ is the ambient temperature. ${\nu }_{set,i}(t)$ is the set temperature for $TC{L}_i$. ${\eta }_i$ is the conversion efficiency from electricity to heat. We assume that the set temperature for a customer can be adjusted from ${\nu }_{\min ,i}$ to ${\nu }_{\max ,i}$ when taking customers’ comfort into consideration. The adjustable power range can be expressed as (5) and (6):

$$L_{low,i}=({\nu }_{\min ,i}-{\nu }_{outside}(t))/{\eta }_i{R}_{eq,i}$$

(5)

$$L_{up,i}=({\nu }_{\max ,i}-{\nu }_{outside}(t))/{\eta }_i{R}_{eq,i}$$

(6)

So constraint for flexible load such as HVAC can be expressed as (7):

$$({\nu }_{\min ,i}-{\nu }_{outside}(t))/{\eta }_i{R}_{eq,i}\le {\chi }_i^t\le ({\nu }_{\max ,i}-{\nu }_{outside}(t))/{\eta }_i{R}_{eq,i}$$

(7)

Considering that the set temperature of HVAC can be changed instantaneously, Constraints (1) and (2) can be neglected or set ${\alpha }_{down}={\alpha }_{up}=M$. $M$ is a large number. Other flexible loads like brightness adjustable lighting can also be expressed as (7).

Apart from former ‘power regulation devices’, energy storage systems (ESSs) are another kind of distributed energy in ADNs. This type of equipment can realize the energy transfer between different time spans by charging and discharging. The characteristics of such equipment can be described by the index of state of charge (SOC). The constraints for SOC can be expressed as (8)–(10):

$$S_{soc,i}^{t+1}=S_{soc,i}^t-({\vartheta }_{i,eq}^t-{\vartheta }_i^r)/e_{\vartheta ,i}^r$$

(8)

$$0\le {\vartheta }_{i,eq}^t\le 2{\vartheta }_i^r$$

(9)

$$-\frac{S_{soc}^{\max }-S_{soc}^t}{{\kappa }_{step}}\cdot e_{\vartheta ,i}^r\le {\vartheta }_{i,eq}^t-{\vartheta }_i^r\le \frac{S_{soc}^t-S_{soc}^{\min }}{{\kappa }_{step}}\cdot e_{\vartheta ,i}^r$$

(10)

Here we introduce ${\vartheta }_{i,eq}^t$ into the formulation to make sure all variables are non-negative and ${\vartheta }_i^t={\vartheta }_{i,eq}^t-{\vartheta }_i^r$. It should be noticed that, $S_{soc,i}^t$ is usually used to express the SOC for energy storage device $i$ by the end of time period $t$. But to simplify the future solution process, while not affecting the variables, $S_{soc,i}^t$ in this paper is defined as the SOC before the start of time period $t$.

In this paper, we merge flexible load with energy storage characteristics such as water heater and electric vehicle with energy storage system as ‘energy storage devices’. For flexible load such as water heater, it is usually used once a day. And it is necessary to reach a specified temperature at a specified time. The state of household water heater can be expressed by temperature ${\nu }_i(t)$. The state of the water heater can be expressed as (11):

$$S_{soc,i}^t=1-({\nu }_i^E-{\nu }_i(t))/({\nu }_i^E-{\nu }_i^0)$$

(11)

${\nu }_i^E$ is the desired temperature for water heater $i$. ${\nu }_i^0$ is the initial temperature for water heater $i$. If we take the sensitivity of temperature to energy into account, the constraints for water heater can also be expressed as (8)–(10). As the water heater cannot release energy, constraint (9) need to be replaced by (12):

$${\vartheta }_i^r\le {\vartheta }_{i,eq}^t\le 2{\vartheta }_i^r$$

(12)

SOC limit for a water heater is $S_{soc}^{\max }=1$ when ${\nu }_i(t)={\nu }_i^E$ and $S_{soc}^{\min }=0$, when ${\nu }_i(t)={\nu }_i^0$. The rated capacity can be expressed as:

$$e_{\vartheta ,i}^r=({\nu }_i^E-{\nu }_i^0)e_{n0,i}$$

(13)

$e_{n0,i}$ is the energy required for the temperature of the water heater increased by one degree. Other flexible load like electric vehicle is similar to energy storage battery. So it is much easier to be reformulated as the mode of ESS. To reflect the customers’ requirement, an extra constraint will be added to origin constraints. The SOC of an electric vehicle $S_{soc,i}^{EV}(t)$ needs to arrive a specific value $S_{SOC,set}^{EV}$ at specific time $t_{set}^{EV}$ as expressed in (14).

$$S_{soc,i}^{EV}(t_{set}^{EV})=S_{SOC,set}^{EV}$$

(14)

The formula can also be meaningful for a BESS. We can assume that BESS should return to the initial state after a cycle of operation to support a continuous operation for next time span. The model of energy storage devices when electric vehicles are merged can add an extra constraint shown as (15).

$$S_{soc,i}^{t_{set}}=S_{soc}^{set}$$

(15)

Apart from former model merging, the uncertainty for the intermittent energy is also considered in this paper described as a solution space for specify scenario. The prediction for intermittent energy is $p_{\psi ,i}^t$. We use the normal distribution or Cauchy distribution to describe the prediction error for intermittent energy centered on the predicted value $p_{\psi ,i}^t$. In this paper, we set a power deviation to make sure the probability for intermittent energy output within the deviation range to be at least 0.95. We can assume the upper deviation of the confidence interval for intermittent energy $i$ is $p_{\psi +,i}^t$ and the lower one is $p_{\psi -,i}^t$. The output range for intermittent energy is $[p_{\psi ,i}^t-p_{\psi -,i}^t,p_{\psi ,i}^t+p_{\psi +,i}^t]$. To limit the conservatism degree for intermittent power output, we introduce the method described in [23] to employ an integer ${\varpi }_i$ as the ‘cardinality budget’. The budget is used to restrict the number of time intervals when the output of intermittent energy reaches the limitation. If ${\varpi }_i$ equals 12, there will be less than 12 time intervals when the output reaches the limitation and for other time intervals the output will be restricted to $p_{\psi ,i}^t$. When ${\varpi }_i\ge 8$, the optimization solution will be feasible with a probability greater than 0.95 [23]. In this paper, we assume ${\eta }_{t,i}^{+}$ and ${\eta }_{t,i}^{-}$ are two binary variables represent the status of deviation for $p_{\psi ,i}^t$. If ${\eta }_{t,i}^{+}=1$, the output of intermittent energy will be $p_{\psi ,i}^t+p_{\psi +,i}^t$. If ${\eta }_{t,i}^{-}=1$, the output of intermittent energy will be $p_{\psi ,i}^t-p_{\psi -,i}^t$, so the uncertainty set $\aleph$ to describe the intermittent energy output can be shown as (16):

$$\begin{array}{c}\aleph :=\left\{p_i^t=p_{\psi ,i}^t+{\eta }_{t,i}^{+}{p}_{\psi ,i}^{+}-{\eta }_{t,i}^{-}{p}_{\psi ,i}^{-}\right\}\\ s.t.{\displaystyle \sum _{t=1}^T({\eta }_{t,i}^{+}+{\eta }_{t,i}^{-})}\le {\varpi }_i,\forall i\in \mathrm{H},\forall t\in T\end{array}$$

(16)

Based on the aforementioned unified model and uncertainty set of intermittent energy, the optimization formulation for optimal operation of active distribution network can be described as follow:

$$\max f=\max {\displaystyle \sum _{t=1}^T\left({\epsilon }_1{f}_1\left({\chi }_i^t\right)+{\epsilon }_2{f}_2(S_{\lambda }^t)\right)}+\underset{p\in \aleph }{\min }\underset{\vartheta \in \Theta }{\max }({\epsilon }_3{f}_3({\vartheta }_{i,eq}^t)+{\epsilon }_4{f}_4(V_{\zeta }))$$

(17)

$$f_1({\chi }_i^t)={\displaystyle \sum _{i\in \Phi }{\lambda }_t\times {\chi }_i^t\times }{\kappa }_{step}$$

(18)

$$f_2(S_{\lambda }^t)=(S_{\lambda ,t}^0-S_{\lambda ,t}^{rslt})\times {\lambda }_t$$

(19)

$$f_3({\vartheta }_{i,eq}^t)={\displaystyle \sum _{t=1}^T\left({\displaystyle \sum _{i\in \Gamma }{\lambda }_t\times ({\vartheta }_{i,eq}^t-{\vartheta }_i^r)\times }{\kappa }_{step}\right)}$$

(20)

$$f_4(V_{\zeta })=(V_{\zeta }^0-V_{\zeta }^{rslt})\times \zeta$$

(21)

The former objection for optimal operation is consist by two parts. The first part includes the profit for power regulation devices (e.g., Equation (18)) and the profit for the reduction of network loss (e.g., Equation (19)). The second part includes the profit for energy storage devices (e.g., Equation (20)) and the profit for the reduction of peak-valley difference (e.g., Equation (21)).The objection is expressed by two parts because the energy storage devices mainly focuses on the effect of peak load shifting in this paper. The power regulation devices preferentially meet the network loss reduction requirement. Since the capacity of intermittent energy is relatively small in the demonstration project, the influence of its uncertainty on network loss is far less than that on peak-valley difference, so the uncertainty of intermittent energy is considered in the second part. The network loss and the peak-valley difference are two main indicators concerned by DSOs when assessing the operation of distribution network. The calculation for both indicators are described in (22)–(24):

$$\{\begin{array}{ll}{S}_{\lambda ,t}^0=f_{loss}(0,{\chi }_i^{Now},p_{\psi ,i}^t)& \forall t\in T(a)\\ S_{\lambda ,t}^{rslt}=f_{loss}(0,{\chi }_i^t,p_{\psi ,i}^t)& \forall t\in T(b)\end{array}$$

(22)

$$\{\begin{array}{ll}{V}_{\zeta }^0=\underset{t}{\max }\left({\sigma }_t^0\right)-\underset{t}{\min }\left({\sigma }_t^0\right)& \forall t\in T(a)\\ V_{\zeta }^{rslt}=\underset{t}{\max }({\sigma }_t^{rslt})-\underset{t}{min}({\sigma }_t^{rslt})& \forall t\in T(b)\end{array}$$

(23)

$$\{\begin{array}{ll}{\sigma }_t^0=D_t-{\displaystyle \sum _{i\in \Phi }{\chi }_i^{Now}}-{\displaystyle \sum _{i\in \mathrm{H}}{p}_{\psi ,i}^t},& \forall t\in T(a)\\ {\sigma }_t^{rslt}=D_t-{\displaystyle \sum _{i\in \Phi }{\chi }_i^t}-{\displaystyle \sum _{i\in \Gamma }({\vartheta }_{i,eq}^t-{\vartheta }_i^r)}-{\displaystyle \sum _{i\in \mathrm{H}}{p}_i^t}& \forall t\in T(b)\end{array}$$

(24)

The initial state before optimization is recorded when the output for power regulation devices maintain at current power ${\chi }_i^{Now}$ (for flexible loads such as HVAC, maintain at the user set value), the output for energy storage devices maintain at 0 (for flexible loads such as EV, the charging time will start immediately until target is reached), the outputs for intermittent energy maintain at $p_{\psi ,i}^t$. Correspondingly, the optimized state is recorded when the output for power regulation devices and energy storage devices are optimized, the output for intermittent energy is $p_i^t$ when taking uncertainty into account. This is to make sure we can get an optimized result when the intermittent energy output reaches the worst scenario. $f_{loss}(•)$ is the calculation formula for network loss, which won’t be discussed in detail in this paper. ${\epsilon }_1$ to ${\epsilon }_4$ are weights for each objective. We can adjust the weights to obtain a desired target combination. Then constraint for former objection is described as a set $\Omega$. Some constraints are aforementioned. We will rewrite it here for a more clearly expression.

$$\Omega (\chi ,\vartheta ,p)=\{(\chi ,\vartheta ,p,S):$$

$$(\chi ,\vartheta ,p)\in \Im$$

(25)

$${\chi }_i^{t-1}-{\chi }_i^t\le {\alpha }_{down,i}\forall t\in T,\forall i\in \Phi$$

(26)

$${\chi }_i^t-{\chi }_i^{t-1}\le {\alpha }_{up,i}\forall t\in T,\forall i\in \Phi$$

(27)

$$L_{low,i}\le {\chi }_i^t\le L_{up,i}\forall t\in T,\forall i\in \Phi$$

(28)

$$S_{soc,i}^{t+1}=S_{soc,i}^t-({\vartheta }_{i,eq}^t-{\vartheta }_i^r)/e_{\vartheta ,i}^r\forall t\in T,\forall i\in \Gamma$$

(29)

$$0\le {\vartheta }_{i,eq}^t\le 2{\vartheta }_i^r\forall t\in T,\forall i\in \Gamma$$

(30)

$$-\frac{S_{soc}^{\max }-S_{soc}^t}{{\kappa }_{step}}\cdot e_{\vartheta ,i}^r\le {\vartheta }_{i,eq}^t-{\vartheta }_i^r\le \frac{S_{soc}^t-S_{soc}^{\min }}{{\kappa }_{step}}\cdot e_{\vartheta ,i}^r\forall t\in T,\forall i\in \Gamma$$

(31)

$$S_{soc,i}^1=s_i^{ini}\forall i\in \Gamma$$

(32)

$$S_{soc,i}^{T+}=s_i^{end}\forall i\in \Gamma$$

(33)

$${\kappa }_{step}=1$$

(34)

$${\chi }_i^t,{\vartheta }_{k,eq}^t,S_{soc,i}^t,S_{soc,i}^{T+},p_l^t\ge 0\forall i\in \Phi ,\forall k\in \Gamma ,\forall l\in \mathrm{H},\forall t\in T\}$$

(35)

$$\Im :=\{P_f\left(t\right)=(\chi ,\vartheta ,p,t)$$

(36)

$$V_{b,\min }\le V_b^t\le V_{b,\max }\forall b\in B,\forall t\in T$$

(37)

$$I_l^t\le I_{l.Rated}\forall l\in L,\forall t\in T\}$$

(38)

Set $\Im$ described in (25) is the constraint for power flow, including power flow equality constraint (e.g., Equation (36)), voltage constraint for buses (e.g., Equation (37)), transmission capacity constraints (e.g., Equation (38)). All the three constraints are used to finalize the optimization results by more stringent constraints when ${\vartheta }_k^t=\pm {\vartheta }_k^r$ and $p_l^t=p_{\psi ,l}^t+p_{\psi +,i}^t\left|\right|p_l^t=p_{\psi ,l}^t-p_{\psi -,i}^t$. The constraints for power regulation devices are described in (26)–(28). The constraints for energy storage devices are described in (29)–(33). The complete representation for optimization formulation can be summarized as (39):

$$\begin{array}{l}\max f=\max {\displaystyle \sum _{t=1}^T({\epsilon }_1\times {\displaystyle \sum _{i\in \Phi }{\lambda }_t\times {\chi }_i^t\times }{\kappa }_{step}}+{\epsilon }_2\times (S_{\lambda ,t}^0-S_{\lambda ,t}^{rslt})\times {\lambda }_t\times {\kappa }_{step})\\ +\underset{p\in \aleph }{\min }\underset{i\in \Gamma }{\max }({\displaystyle \sum _{t=1}^T\left({\epsilon }_3\times {\displaystyle \sum _{i\in \Gamma }{\lambda }_t\times ({\vartheta }_{i,eq}^t-{\vartheta }_i^r)\times }{\kappa }_{step}\right)}+{\epsilon }_4\times ((\max ({\sigma }_t^0)-\min ({\sigma }_t^0))\\ -(\max ({\sigma }_t^{rslt})-min({\sigma }_t^{rslt})))\times \zeta )\\ s.t.\\ \mathrm{constraint}(16)\\ \mathrm{constraint}(24)\\ \mathrm{constraints}(26)\mathrm{to}(38)\end{array}$$

(39)

3. Solution Methodology

The formerly proposed equation it is a relatively complicated problem with several extremum calculations, thus the equation cannot be solved directly. The uncertainty for intermittent energy has increased the solving difficulty of the problem significantly. To decrease the solving complexity, Benders decomposition is introduced to this paper.

Benders decomposition is an effective method to solve complex programming problem [24]. An example of Benders decomposition can be expressed as (40):

$$\begin{array}{cc}\hfill \min & f(x)+c^{\mathrm{T}}y\hfill \\ \hfill s.t.& F(x)+Ay\le b\hfill \\ & x\in S,y\ge 0\hfill \end{array}$$

(40)

The Benders decomposition method will decompose the original optimization problem into several parts to be solved. By the process of iteration, the solution of reformulated problem will gradually approximate the optimal solution of the original problem. This method requires the separation of variables shown as $x$ and $y$, both variables can be vectors. $f(x)$ and $F(x)$ are the functions for variable $x$. $A$ is the coefficient matrix; $b$, $c$ is the coefficient vector. $S$ is the feasible set for variable $x$. Then the original problem (40) can be divided into two parts by Benders decomposition. Main problem is defined on the set $S$, it can be linear, non-linear or discrete. The problem is formulated when $y=\overline{y}$, where $\overline{y}$ is the solution of the sub-problem:

$$\begin{array}{cc}\hfill \min & f(x)+c^{\mathrm{T}}\overline{y}\hfill \\ \hfill s.t.& F(x)+A\overline{y}\le b\hfill \\ & x\in S\hfill \end{array}$$

(41)

The sub-problem should be a linear programming, shown as (42):

$$\begin{array}{cc}\hfill \min & f(\overline{x})+c^{\mathrm{T}}y\hfill \\ \hfill s.t.& Ay\le b-F(\overline{x})\hfill \\ & y\ge 0\hfill \end{array}$$

(42)

Firstly, the algorithm will solve the duration problem of sub-problem (42) based on a “possible solution” of the main problem (41). This “possible solution” is one of the feasible solutions for main problem. Then, additional constraints (feasibility cuts and optimal cuts) will be formed by solving the sub-problem and new constraints will be added to main problem. The formation of both cuts will be discussed in detail later. Furthermore, solve the main problem and bring the optimized results $\overline{x}$ back to sub-problem (42) again, and so on. Solve the two questions alternately until the optimal solution is obtained.

Based on the aforementioned algorithm, the main problem and the sub-problem for the optimization problem proposed in this paper are expressed as (43) and (44).

Main problem:

$$\begin{array}{l}\max f=\max {\displaystyle \sum _{t=1}^T({\epsilon }_1\times {\displaystyle \sum _{i\in \Phi }{\lambda }_t\times {\chi }_i^t\times }{\kappa }_{step}}+{\epsilon }_2\times (S_{\lambda ,t}^0-S_{\lambda ,t}^{rslt})\times {\lambda }_t\times {\kappa }_{step})+\iota \\ s.t.\\ \mathrm{constraints}(26)\mathrm{to}(28)\\ \mathrm{constraints}(34)\mathrm{to}(38)\end{array}$$

(43)

Sub-problem:

$$\begin{array}{l}\iota =\underset{p\in \aleph }{\min }\underset{i\in \Gamma }{\max }G({\vartheta }_{i,eq}^t)\\ =\underset{p\in \aleph }{\min }\underset{i\in \Gamma }{\max }({\displaystyle \sum _{t=1}^T\left({\epsilon }_3\times {\displaystyle \sum _{i\in \Gamma }{\lambda }_t\times ({\vartheta }_{i,eq}^t-{\vartheta }_i^r)\times }{\kappa }_{step}\right)}+{\epsilon }_4\times ((\max ({\sigma }_t^0)-\min ({\sigma }_t^0))\\ -(\max ({\sigma }_t^{rslt})-min({\sigma }_t^{rslt})))\times \zeta )\\ s.t.\\ \mathrm{constraint}(16)\\ \mathrm{constraint}(24)\\ \mathrm{constraints}(29)\mathrm{to}(35)\end{array}$$

(44)

The main problem takes the outputs of ‘power regulation devices’ as variables and the sub-problem takes the outputs of both energy storage devices and intermittent energies as variables.

Firstly, the sub-problem needs to be simplified to a single extremum value problem. In the sub-system $V_{\zeta }^0=\max \left({\sigma }_t^0\right)-\min \left({\sigma }_t^0\right)$ is determined by the predicted value, therefore belongs to a constant. To express $\max ({\sigma }_t^{rslt})-min({\sigma }_t^{rslt})$ by inequality, we introduce extra variable and constraints in (45)–(47):

$$\begin{array}{l}\max ({\sigma }_t^{rslt})-\min ({\sigma }_t^{rslt})\\ =\min (h_{up}-h_{low})\end{array}$$

(45)

$$h_{low}\le D_t-{\displaystyle \sum _{i\in \Phi }{\chi }_i^t}-{\displaystyle \sum _{i\in \Gamma }({\vartheta }_{i,eq}^t-{\vartheta }_i^r)}-{\displaystyle \sum _{i\in \mathrm{H}}{p}_i^t}\forall t\in T$$

(46)

$$h_{up}\ge D_t-{\displaystyle \sum _{i\in \Phi }{\chi }_i^t}-{\displaystyle \sum _{i\in \Gamma }({\vartheta }_{i,eq}^t-{\vartheta }_i^r)}-{\displaystyle \sum _{i\in \mathrm{H}}{p}_i^t}\forall t\in T$$

(47)

Therefore problem (44) can be described as (48):

$$\begin{array}{l}\underset{p\in \aleph }{\min }\underset{i\in \Gamma }{\max }\{{\displaystyle \sum _{t=1}^T\left({\epsilon }_3\times {\displaystyle \sum _{i\in \Gamma }{\lambda }_t\times ({\vartheta }_{i,eq}^t-{\vartheta }_i^r)\times }{\kappa }_{step}\right)}\\ +(V_{\zeta }^0-\min \{h_{up}-h_{low}\})\times \zeta \times {\epsilon }_4\}\\ s.t.\\ \mathrm{constraint}(16)\\ \mathrm{constraints}(29)\mathrm{to}(35)\\ \mathrm{constraints}(46)\mathrm{to}(47)\end{array}$$

(48)

Furthermore, (48) can be transformed as (49):

$$\begin{array}{l}\underset{p\in \aleph }{\min }\underset{i\in \Gamma }{\max }\{{\displaystyle \sum _{t=1}^T\left({\epsilon }_3\times {\displaystyle \sum _{i\in \Gamma }{\lambda }_t\times ({\vartheta }_{i,eq}^t-{\vartheta }_i^r)\times }{\kappa }_{step}\right)}+(V_{\zeta }^0-\min \{h_{up}-h_{low}\})\times \zeta \times {\epsilon }_4\}\\ =\underset{p\in \aleph }{\min }\underset{i\in \Gamma }{\max }({\displaystyle \sum _{t=1}^T\left({\displaystyle \sum _{i\in \Gamma }{\lambda }_t\times ({\vartheta }_{i,eq}^t-{\vartheta }_i^r)\times }{\kappa }_{step}\right)}\times {\epsilon }_3+V_{\zeta }^0\times \zeta \times {\epsilon }_4-h_{up}\times \zeta \times {\epsilon }_4+h_{low}\times \zeta \times {\epsilon }_4)\\ =\underset{p\in \aleph }{\min }\underset{i\in \Gamma }{\max }(V_{\zeta }^0\times \zeta \times {\epsilon }_4-{\displaystyle \sum _{t=1}^T\left({\lambda }_t\times {\displaystyle \sum _{i\in \Gamma }{\vartheta }_i^r}\times {\kappa }_{step}\right)\times {\epsilon }_3}+{\displaystyle \sum _{t=1}^T\left({\lambda }_t\times {\displaystyle \sum _{i\in \Gamma }{\vartheta }_{i,eq}^t}\times {\kappa }_{step}\right)}\times {\epsilon }_3\\ -h_{up}\times \zeta \times {\epsilon }_4+h_{low}\times \zeta \times {\epsilon }_4)\\ =\underset{p\in \aleph }{\min }\underset{i\in \Gamma }{\max }G’({\vartheta }_{i,eq}^t，h_{up}，h_{low})\end{array}$$

(49)

To solve the further min-max problem, we use the dual approach to translate the maximize problem of $G’({\vartheta }_{i,eq}^t，h_{up}，h_{low})$ to a minimize problem. The solution space is defined as $\Psi$:

$$\begin{array}{l}\underset{i\in \Gamma }{\max }G’({\vartheta }_{i,eq}^t,h_{up},h_{low})\\ =\underset{\Psi }{\min }F({\omega }_t^{+},{\omega }_t^{-},{\tau }_i^t,{\mu }_{+,i}^t,{\mu }_{-,i}^t,{\delta }_i^t,{\delta }_i^{begin},{\delta }_i^{end})\\ =V_{\zeta }^0\times \zeta \times {\epsilon }_4-{\displaystyle \sum _{t=1}^T\left({\lambda }_t\times {\displaystyle \sum _{i\in \Gamma }{\vartheta }_i^r}\right)\times {\epsilon }_3}+{\displaystyle \sum _{t=1}^T\left({\omega }_t^{+}\times (D_t-{\displaystyle \sum _{i\in \Phi }{\chi }_i^t}-{\displaystyle \sum _{i\in \mathrm{H}}{p}_i^t}+{\displaystyle \sum _{i\in \Gamma }{\vartheta }_i^r})\right)}\\ -{\displaystyle \sum _{t=1}^T\left({\omega }_t^{-}\times (D_t-{\displaystyle \sum _{i\in \Phi }{\chi }_i^t}-{\displaystyle \sum _{i\in \mathrm{H}}{p}_i^t}+{\displaystyle \sum _{i\in \Gamma }{\vartheta }_i^r})\right)}+{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i\in \Gamma }2{\tau }_i^t}{\vartheta }_i^r}+{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i\in \Gamma }{\mu }_{+,i}^t({\vartheta }_i^r-S_{soc,i}^{\min }{e}_{\vartheta ,i}^r)}}\\ +{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i\in \Gamma }{\mu }_{-,i}^t(S_{soc,i}^{\max }{e}_{\vartheta ,i}^r-{\vartheta }_i^r)}}+{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i\in \Gamma }{\delta }_i^t({\vartheta }_i^r/e_{\vartheta ,i}^r)}}+{\displaystyle \sum _{i\in \Gamma }{s}_i^{ini}{\delta }_i^{begin}}+{\displaystyle \sum _{i\in \Gamma }{s}_i^{end}{\delta }_i^{end}}\end{array}$$

(50)

$s.t.$

$${\displaystyle \sum _{t=1}^T{\omega }_t^{+}}\ge \zeta \times {\epsilon }_4$$

(51)

$${\displaystyle \sum _{t=1}^T{\omega }_t^{-}}\le \zeta \times {\epsilon }_4$$

(52)

$${\omega }_t^{+}-{\omega }_t^{-}+{\tau }_i^t+{\mu }_{+,i}^t-{\mu }_{-,i}^t+(1/e_{\vartheta ,i}^r){\delta }_i^t\ge {\lambda }_t\times {\epsilon }_3\forall i\in \Gamma \forall t\in T$$

(53)

$$-e_{\vartheta ,i}^r{\mu }_{+,i}^t+e_{\vartheta ,i}^r{\mu }_{-,i}^t+{\delta }_i^{t-1}-{\delta }_i^t\ge 0\forall i\in \Gamma \forall t\in [2,T]$$

(54)

$$-e_{\vartheta ,i}^r{\mu }_{+,i}^1+e_{\vartheta ,i}^r{\mu }_{-,i}^1-{\delta }_i^1+{\delta }_i^{begin}\ge 0$$

(55)

$${\delta }_i^T+{\delta }_i^{end}\ge 0$$

(56)

$${\omega }_t^{+},{\omega }_t^{-},{\tau }_i^t,{\mu }_{+,i}^t,{\mu }_{-,i}^t\ge 0\forall i\in \Gamma ,\forall t\in T$$

(57)

In the former equation, ${\omega }_t^{+}$ and ${\omega }_t^{-}$ correspond to the dualizing of constrains (46) and (47). ${\mu }_{+,i}^t$ and ${\mu }_{-,i}^t$ correspond to constraint (31). ${\tau }_i^t$ corresponds to constraint (30). ${\delta }_i^t$, ${\delta }_i^{begin}$, ${\delta }_i^{end}$ correspond to constraints (29) and (32) and (33). Constraints (51) to (56) correspond to initial variables $h_{low}$, $h_{up}$, ${\vartheta }_{i,eq}^t$ and $S_{soc,i}^t$. In the dual problem formation process, $p_i^t$ is considered as a known quantity, when taking the uncertainty set of intermittent energy into consideration, (50) is a bilinear problem. To linearize the equation, we first separate the bilinear part of problem as (58):

$$\begin{array}{l}\underset{p\in \aleph }{\min }\underset{\Psi }{\min }F({\omega }_t^{+},{\omega }_t^{-},{\tau }_i^t,{\mu }_{+,i}^t,{\mu }_{-,i}^t,{\delta }_i^t,{\delta }_i^{begin},{\delta }_i^{end})\\ =\min F({\omega }_t^{+},{\omega }_t^{-},{\tau }_i^t,{\mu }_{+,i}^t,{\mu }_{-,i}^t,{\delta }_i^t,{\delta }_i^{begin},{\delta }_i^{end},p_i)\\ =\min (L({\omega }_t^{+},{\omega }_t^{-},{\tau }_i^t,{\mu }_{+,i}^t,{\mu }_{-,i}^t,{\delta }_i^t,{\delta }_i^{begin},{\delta }_i^{end})-{\displaystyle \sum _{t=1}^T({\omega }_t^{+}\times {\displaystyle \sum _{i\in \mathrm{H}}{p}_i^t})}+{\displaystyle \sum _{t=1}^T({\omega }_t^{-}\times {\displaystyle \sum _{i\in \mathrm{H}}{p}_i^t})})\\ =\min (L({\omega }_t^{+},{\omega }_t^{-},{\tau }_i^t,{\mu }_{+,i}^t,{\mu }_{-,i}^t,{\delta }_i^t,{\delta }_i^{begin},{\delta }_i^{end})+{\displaystyle \sum _{t=1}^T(({\omega }_t^{-}-{\omega }_t^{+})\times {\displaystyle \sum _{i\in \mathrm{H}}{p}_i^t})})\\ =\min (L({\omega }_t^{+},{\omega }_t^{-},{\tau }_i^t,{\mu }_{+,i}^t,{\mu }_{-,i}^t,{\delta }_i^t,{\delta }_i^{begin},{\delta }_i^{end})+{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i\in \mathrm{H}}(({\omega }_t^{-}-{\omega }_t^{+})\times p_i^t)}})\end{array}$$

(58)

The bilinear pars of (58) can be reformulated by introduce some extra variables ${\overline{\omega }}_t$, ${\overline{\omega }}_t^{+}$, ${\overline{\omega }}_t^{-}$, ${\eta }_{t,i}^{+}$, and ${\eta }_{t,i}^{-}$. Among which ${\eta }_{t,i}^{+}$, and ${\eta }_{t,i}^{-}$ are binary variables:

$$\begin{array}{l}\min {\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i\in \mathrm{H}}\left(({\omega }_t^{-}-{\omega }_t^{+})\times p_i^t\right)}}\\ =\min {\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i\in \mathrm{H}}({\omega }_t^{-}-{\omega }_t^{+})\times (p_{\psi ,i}^t+{\eta }_{t,i}^{+}{p}_{\psi ,i}^{+}-{\eta }_{t,i}^{-}{p}_{\psi ,i}^{-})}}\\ =\min ({\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i\in \mathrm{H}}{\overline{\omega }}_t\times p_{\psi ,i}^t}}+{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i\in \mathrm{H}}{\overline{\omega }}_t^{+}\times p_{\psi ,i}^{+}}}+{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i\in \mathrm{H}}{\overline{\omega }}_t^{-}\times p_{\psi ,i}^{-}}})\\ =\min R({\overline{\omega }}_t,{\overline{\omega }}_t^{+},{\overline{\omega }}_t^{-})\end{array}$$

(59)

$${\overline{\omega }}_t={\omega }_t^{-}-{\omega }_t^{+}$$

(60)

$${\overline{\omega }}_t^{+}\ge -M{\eta }_{t,i}^{+}$$

(61)

$${\overline{\omega }}_t^{-}\ge -M{\eta }_{t,i}^{-}$$

(62)

$${\overline{\omega }}_t^{+}\ge {\overline{\omega }}_t-M(1-{\eta }_{t,i}^{+})$$

(63)

$${\overline{\omega }}_t^{-}\ge -{\overline{\omega }}_t-M(1-{\eta }_{t,i}^{-})$$

(64)

$${\displaystyle \sum _{t=1}^T({\eta }_{t,i}^{+}+{\eta }_{t,i}^{-})}\le {\varpi }_i$$

(65)

$${\eta }_{t,i}^{+},{\eta }_{t,i}^{-}\in \{0,1\}$$

(66)

Now we have already translate the former multi-extremum equation to a single minimum equlation. The optimization equation can be summarized as (67):

$$\begin{array}{l}Y({\chi }_i^t)=\underset{p\in \aleph }{\min }\underset{\Psi }{\min }F({\omega }_t^{+},{\omega }_t^{-},{\tau }_i^t,{\mu }_{+,i}^t,{\mu }_{-,i}^t,{\delta }_i^t,{\delta }_i^{begin},{\delta }_i^{end})\\ =\min (L({\omega }_t^{+},{\omega }_t^{-},{\tau }_i^t,{\mu }_{+,i}^t,{\mu }_{-,i}^t,{\delta }_i^t,{\delta }_i^{begin},{\delta }_i^{end})+R({\overline{\omega }}_t,{\overline{\omega }}_t^{+},{\overline{\omega }}_t^{-}))\\ s.t.\\ \mathrm{constraints}(51)\mathrm{to}(57)\\ \mathrm{constraints}(60)\mathrm{to}(66)\end{array}$$

(67)

When we use Benders decomposition algorithm to solve the problem, the sub-problem will formulate the feasibility-cuts and optimality-cuts. Take the problem shown in Equations (40)–(42) as example. The duration for sub-problem (42) can be expressed as:

$$\begin{array}{cc}\hfill \max & {(b-F(\overline{x}))}^{\mathrm{T}}u\hfill \\ \hfill s.t.& A^{\mathrm{T}}u\le c\hfill \\ & u\ge 0\hfill \end{array}$$

(68)

According to the dual principle, when problem (68) has an unbounded solution, the sub-problem (42) has no solution, which means original problem (40) also has no solution. In this situation, an extra cuts ${(b-F(\overline{x}))}^{\mathrm{T}}{u}^{fsb}\le 0$ will be added to the main problem as feasibility-cuts. $u^{fsb}$ is the corresponding solution when the objective function is unbounded, also known as extreme direction. When problem (68) has a bounded solution, in order to ensure the optimization result, an extra cuts $z\ge {(b-F(\overline{x}))}^{\mathrm{T}}{u}^{opt}$ will be added to the main problem as optimality-cuts. $u^{opt}$ is the solution when the objective function get extreme value, also known as optimal pole. For the main problem, $y$ is eliminated but not completely disappeared. The optimization results for the sub-problem will feedback to the main problem by the form of feasibility-cuts and optimality-cuts. Each iteration will form new cuts. These cuts will accumulate until the optimal solution is obtained. The express for main problem considering additional cuts can be expressed as:

$$\begin{array}{ccc}\hfill \min & f(x)+z\hfill & \\ \hfill s.t.& z\ge {(b-F(\overline{x}))}^{\mathrm{T}}{u}^{opt}\hfill & optimality-cuts\hfill \\ & {(b-F(\overline{x}))}^{\mathrm{T}}{u}^{fsb}\le 0\hfill & feasibility-cuts\hfill \\ & x\in S\hfill & \end{array}$$

(69)

Therefore the main problem proposed in this paper can be reformulated as:

$$\begin{array}{l}\max f=\max {\displaystyle \sum _{t=1}^T({\epsilon }_1\times {\displaystyle \sum _{i\in \Phi }{\lambda }_t\times {\chi }_i^t\times }{\kappa }_{step}}+{\epsilon }_2\times (S_{\lambda ,t}^0-S_{\lambda ,t}^{rslt})\times {\lambda }_t\times {\kappa }_{step})+\iota \\ s.t.\\ \mathrm{constraints}(26)\mathrm{to}(28)\\ \mathrm{constraints}(34)\mathrm{to}(38)\\ feasibility-cuts\\ optimality-cuts\end{array}$$

(70)

(1) Feasibility-cuts formulation:

We establish a relaxation function for the original sub-problem to check the feasibility, which is also called L-shaped method. Here, we don’t care about the objection, so the problem can be formulated as:

$$\begin{array}{l}\underset{i\in \Gamma }{\max }(-{\displaystyle \sum _{c=1}^2{\alpha }_{c+}^t}-{\displaystyle \sum _{c=3}^6{\displaystyle \sum _{i\in \Gamma }{\alpha }_{c+}^{i,t}}})-({\displaystyle \sum _{c=1}^2{\alpha }_{c-}^t}+{\displaystyle \sum _{c=3}^6{\displaystyle \sum _{i\in \Gamma }{\alpha }_{c-}^{i,t}}})\\ s.t.\\ h_{low}+{\alpha }_{1+}^t-{\alpha }_{1-}^t\le D_t-{\displaystyle \sum _{i\in \Phi }{\chi }_i^t}-{\displaystyle \sum _{i\in \Gamma }({\vartheta }_{i,eq}^t-{\vartheta }_i^r)}-{\displaystyle \sum _{i\in \mathrm{H}}{p}_i^t}\\ h_{up}+{\alpha }_{2+}^t-{\alpha }_{2-}^t\ge D_t-{\displaystyle \sum _{i\in \Phi }{\chi }_i^t}-{\displaystyle \sum _{i\in \Gamma }({\vartheta }_{i,eq}^t-{\vartheta }_i^r)}-{\displaystyle \sum _{i\in \mathrm{H}}{p}_i^t}\\ S_{soc,i}^{t+1}+{\alpha }_{3+}^{i,t}-{\alpha }_{3-}^{i,t}=S_{soc,i}^t-({\vartheta }_{i,eq}^t-{\vartheta }_i^r)/e_{\vartheta ,i}^r\\ {\vartheta }_{i,eq}^t+{\alpha }_{4+}^{i,t}-{\alpha }_{4-}^{i,t}\le 2{\vartheta }_i^r\\ {\vartheta }_{i,eq}^t-{\vartheta }_i^r+{\alpha }_{5+}^{i,t}-{\alpha }_{5-}^{i,t}\le \frac{S_{soc}^t-S_{soc}^{\min }}{{\kappa }_{step}}\cdot e_{\vartheta ,i}^r\\ {\vartheta }_i^r-{\vartheta }_{i,eq}^t+{\alpha }_{6+}^{i,t}-{\alpha }_{6-}^{i,t}\le \frac{S_{soc}^{\max }-S_{soc}^t}{{\kappa }_{step}}\cdot e_{\vartheta ,i}^r\\ {\alpha }_{c+}^t,{\alpha }_{c-}^t\ge 0\forall t\in T\\ {\alpha }_{c+}^{i,t},{\alpha }_{c-}^{i,t}\ge 0\forall t\in T,\forall i\in \Gamma \\ \mathrm{constraints}(32)\mathrm{to}(35)\end{array}$$

(71)

The dual form of the aforementioned problem can be expressed as (72):

$$\begin{array}{l}Y‘({\chi }_i^t)=\min (L({\omega }_t^{+},{\omega }_t^{-},{\tau }_i^t,{\mu }_{+,i}^t,{\mu }_{-,i}^t,{\delta }_i^t,{\delta }_i^{begin},{\delta }_i^{end})+R({\overline{\omega }}_t,{\overline{\omega }}_t^{+},{\overline{\omega }}_t^{-}))\\ s.t.\\ {\displaystyle \sum _{t=1}^T{\omega }_t^{+}}\ge 0\\ {\displaystyle \sum _{t=1}^T(-{\omega }_t^{-})}\ge 0\\ {\omega }_t^{+}-{\omega }_t^{-}+{\tau }_i^t+{\mu }_{+,i}^t-{\mu }_{-,i}^t+(1/e_{\vartheta ,i}^r){\delta }_i^t\ge 0\forall i\in \Gamma \forall t\in T\\ \mathrm{constraints}(54)\mathrm{to}(56)\\ \mathrm{constraints}(60)\mathrm{to}(66)\\ {\omega }_t^{+},{\omega }_t^{-},{\tau }_i^t,{\mu }_{+,i}^t,{\mu }_{-,i}^t\in [0,1]\forall i\in \Gamma ,\forall t\in T\end{array}$$

(72)

If $Y‘({\chi }_i^t)=0$, ${\chi }_i^t$ is a feasible solution for the check problem. If $Y‘({\chi }_i^t)<0$, the feasibility-cuts formulation can be expressed as:

$$Y‘({\chi }_i^t)\ge 0$$

(73)

(2) Optimality-cuts formulation:

As shown in (70), we use $\iota$ represent the optimized value for the sub-system. In the k-th iteration, we solve the main problem and introduce ${\chi }_i^t$ to the sub-problem (67). If the solution of the sub-problem is smaller than ${\iota }_i$, then this solution is not an optimal one, so the optimality-cuts can be expressed as follows:

$$Y({\chi }_i^t)\ge \iota$$

(74)

The computational algorithm framework is shown as Figure 1. The core steps include main problem solving, sub-problem solving and constraint set updating.

4. Demonstration Application for Optimal Operation

The demonstration application for active distribution network located in Guizhou, China. The strategy of stochastic optimal operation is verified based on the structure of the demonstration. The schematic diagram for the demonstration project is shown in Figure 2. The distributed generation and flexible load conditions are listed in

Table 1. Condition for Distributed Generation and Flexible Load in Demonstration Project.

Name	Types of Devices		Capacity of Devices
PV1	Photovoltaic	Intermittent energies	32.4 kW
PV2			91.8 kW
PV3			32.4 kW
PV4			91.8 kW
WTP1	Wind turbine		100 kW
WTP2			100 kW
BESS1	Battery energy storage system	Energy storage devices	100 kW/200 kWh
BESS2			100 kW/200 kWh
EV1	Electric vehicle		100 kW/200 kWh
EV2			10 kW/30 kWh
EV3			10 kW/30 kWh
EV4			10 kW/30 kWh
Hydro	Hydropower		10,000 kW
CCHP	Cooling-heating-power supply	Power regulation devices	500 kW
Flexild	HVAC

. All electric vehicles support vehicle-to-grid (V2G) mode.

The electricity price adopts the peak valley electricity price as Shanghai, the electricity price is ¥0.68 per kWh from 9:00 to 22:00, and electricity price is ¥0.33 per kWh for the rest of the time. $\zeta$ is set to ¥0.167 per kWh [25]. The controllable range for HVACs in the demonstration are shown in Figure 3.

The upper and lower deviations for intermittent energies are shown in Figure 4. It should be noticed that the deviation for photovoltaic during the night is 0, as it is obvious that the photovoltaic system cannot generate electricity during the night.

4.1. Basic Optimization Scenario

The optimization effect for stochastic operation strategy of ADN is firstly discussed under the basic optimization scenario. The weight selection for objective function (${\epsilon }_1$ to ${\epsilon }_4$) is set as

Table 2. Weight Selection for Optimization.

Name	Weight Selection
${\epsilon }_1$	0.005
${\epsilon }_2$	0.495
${\epsilon }_3$	0.25
${\epsilon }_4$	0.25

.

As the benefit brought by controllable devices is significantly higher than the network loss reduce benefit when ${\epsilon }_1$ is close to ${\epsilon }_2$, we select the weights by the sensitivity of distributed power output to network loss to ensure both the distributed power supply gains and loss benefits are considered simultaneously.

The uncertainty for intermittent energy ${\varpi }_i$ is set to 8.

Table 3. Intermittent energy uncertainty when ${\varpi }_i=8$.

Name	Intermittent Energy Uncertainty
	Time Period When ${\mathit{\eta }}_{\mathit{t},\mathit{i}}^{+}=1$	Time Period When ${\mathit{\eta }}_{\mathit{t},\mathit{i}}^{-}=1$
PV1	19	14,15,16,17
PV2	19,20	13,14,15,16,17,18
PV3	19	14,16,17
PV4	19,20	13,14,15,16,18
WTP1	19,20	13,14,15,16,17
WTP2	19,20	13,14,15,16,17

showed the worst case after the solution for intermittent energy. It can be seen that the worst case for intermittent tend to have a negative deviation (${\eta }_{t,i}^{+}=1$) when the load is relatively larger and have a positive deviation (${\eta }_{t,i}^{-}=1$) when the load is relatively smaller.

The optimization result for both energy storage devices and power regulation devices are shown in Figure 5 and Figure 6.

It should be noticed that CCHP is connected to the ShuiPei feeder. As the rated power for CCHP is relatively small compared to the load along the feeder, when the active power for CCHP is rising, the benefit for power regulation devices and the benefit for network loss are increasing at the same time, so the active power for CCHP rose to the upper limit as shown in Figure 6a.

The optimization effects are shown in Figure 7 and Figure 8. The initial state before optimization is recorded as (22a) and (23a). The total load after the optimization is shown as the red line in Figure 7 that has already taken the worst intermittent energy case into consideration. The peak-valley difference has reduced by 1531.36 kW (from 4218.08 kW to 2686.72 kW). The network loss has reduced by 591.689 kW (from 1538.78 kW to 947.091 kW).

Then the influence of both intermittent energy power output uncertainty and the objective function weight selection to the optimization are discussed under different scenarios.

4.2. Different Intermittent Energy Uncertainty

In this part of the paper we select two more representative scenarios where ${\varpi }_i=0$ and ${\varpi }_i=16$ to discuss the influence of different uncertainty for intermittent energy ${\varpi }_i$. ${\varpi }_i$ represents the number of time intervals when the intermittent energy output reaches the limits $p_{\psi ,i}^t-p_{\psi -,i}^t$ or $p_{\psi ,i}^t+p_{\psi +,i}^t$.

When ${\varpi }_i=0$ the active power for all intermittent energy will be strictly the predicted value, which means there is no ‘worst scenario’ for the optimization equlation. Figure 9 shows the different optimization results in ADN when ${\varpi }_i=0$ compared to the loads in ADN when ${\varpi }_i=8$. The peak-valley difference is reduced to 2561.18 kW when ${\varpi }_i=0$, which was better than ${\varpi }_i=8$.

The absolute difference of load is shown in Figure 10 for a more explicit analysis. It can be observed that the time periods when the load difference is relatively large are basically the same periods of intermittent energy prediction uncertainty selection periods.

As the active power for intermittent energy is small, the network loss wouldn’t be influenced that much. The total network loss when ${\varpi }_i=0$ is 949.29 kW, which is only a 0.23% difference compared to the network loss when ${\varpi }_i=8$.

When ${\varpi }_i=16$, the intermittent power supply uncertainty has more choices, which will lead to a worse scenario shown in

Table 4. Intermittent Energy Uncertainty when ${\varpi }_i=16$.

Name	Intermittent Energy Uncertainty
	Time Period when ${\mathit{\eta }}_{\mathit{t},\mathit{i}}^{+}=1$	Time Period when ${\mathit{\eta }}_{\mathit{t},\mathit{i}}^{-}=1$
PV1	7,10,11,12,19	8,9,13,14,15,16,17,18
PV2	7,10,12,19	8,9,13,14,15,16,17,18
PV3	7,10,12,19	8,9,13,14,15,16,17,18
PV4	7,10,12,19	8,9,13,14,15,16,17,18
WTP1	19,20	13,14,15,16,17
WTP2	19,20	13,14,15,16,17

.

Figure 11 compares the optimization result when ${\varpi }_i=16$ with the result when ${\varpi }_i=8$. The peak-valley difference increased to 2714.12 kW, which is the worst among all three scenarios, but the network loss is 945.74 kW, which has only a −0.14% difference compared to the network loss when ${\varpi }_i=8$.

4.3. Different Objective Function Weight Selection

In this part, we will discuss the influence of different objective function weight selections. In all scenarios of this part, ${\varpi }_i$ is set to 8 to exclude the impact of intermittent energy uncertainty.

The first scenario is set to neglect the benefit of all controllable devices including both energy storage devices and power regulation devices. This scenario happens when all controllable devices are owned by the user, so the benefit of controllable devices does not affect the DSO’s choice of optimization strategy. The weight selection for the objective function is set as in

Table 5. Weight Selection for Optimization.

Name	Weight Selection
${\epsilon }_1$	0
${\epsilon }_2$	0.5
${\epsilon }_3$	0
${\epsilon }_4$	0.5

.

We take the basic optimization scenario discussed in the first part of this section as the initial scenario. The optimization results for total load in ADN is shown in Figure 12.

The peak-valley difference in this scenario is reduced to 2648.93 kW, which is better than the initial scenario. It can be seen from Figure 12 that the total load in ADN in the first scenario is larger than the total load in the initial scenario. That’s because as the benefit brought by controllable devices was neglected, the objective function wouldn’t encourage the output of controllable devices. The network loss comparison is shown in Figure 13. The total network loss reduced from 947.091 kW in the initial scenario to 834.38 kW in this scenario.

The second scenario is set to be more biased in favor of the benefit brought by controllable devices. The weight selection for both ${\epsilon }_1$ and ${\epsilon }_3$ are set to be higher as in

Table 6. Weight Selection for Optimization.

Name	Weight Selection
${\epsilon }_1$	0.01
${\epsilon }_2$	0.49
${\epsilon }_3$	0.375
${\epsilon }_4$	0.125

.

The optimization results for both load and network loss in ADN are shown in Figure 14 and Figure 15.

As the objective function was set to encourage the benefit brought by controllable devices, the active power for controllable devices appears to rise. The total load in ADN decreased as shown in Figure 14. The peak-valley difference in this scenario was increased to 2711.80 kW. The network loss was also increased to 1279.946 kW. The active power for different distributed generations are shown in Figure 16 and Figure 17. The total active power is larger than in the initial scenario.

Based on the former result and discussion, different weight selection will lead to a different optimization effect. In practice, we can flexibly adjust the weight as needed in order to obtain an optimized result.

4.4. Optimization Results Analysis

The summary of optimization results for the former scenarios is shown in Figure 18. The optimal operation strategy proposed in this paper realizes an efficient usage of all kinds of controllable resources. It can be seen from the figure that both network losses and the peak-valley differences are improved in all scenarios compared to the original operation state.

The left part of Figure 18a,b show the influence of the intermittent energy uncertainties on the optimization result. The increase of intermittent energy uncertainties will lead to the decrease of the optimization effect, while increasing the compatibility for intermittent energy fluctuations in ADN at the same time.

The right part of Figure 18a,b shows the influence of different objective function weight selection on the optimization results. By adjusting the objective function weight, we can choose to obtain better power grid operation effect or to encourage more distributed power supply.

5. Conclusions

In this paper, we have developed an optimal operation strategy to realize an efficient usage of all kinds of controllable resources in active distribution networks including distributed generations and flexible loads. To decrease the formulation complexity, the models for all resources are summarized as three kinds of devices: intermittent energy, energy storage devices and power regulation devices. Uncertainty for intermittent energy is taken into account in this paper to ensure a high utilization of intermittent energy even under the worst scenario. The prediction error can be flexibly adjusted according to the requirements. Then, we establish an objective function taking two main indicators reflecting the operation state of distribution network into consideration, the peak-valley difference and the network loss. We can set the weight of optimization objection to realize different objective function selection. The mathematical equations are reformulated to a two-stage max-min problem. The Benders decomposition algorithm is used in this paper to solve the nonlinear programming problem. Finally, the robustness of the strategy is verified in the demonstration application project for ADN located in Guizhou, China. The simulation is performed under different simulation settings to ensure an optimal operation for ADN. However, in all demand side loads only the controllable loads are taken into consideration in this paper. In further research we will take both demand responses (such as price incentives) and direct load control as responsible resources to realize comprehensive optimization for both active control resources and passive response resources.