Multi-period Two-stage Robust Optimization of Radial Distribution System with Cables Considering Time-of-use Price

Jian Zhang; Mingjian Cui; Yigang He; Fangxing Li

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

Multi-period Two-stage Robust Optimization of Radial Distribution System with Cables Considering Time-of-use Price PDF

- ORCID：
Jian Zhang
✉
- ORCID：
Mingjian Cui
✉
- ORCID：
Yigang He
✉
- ORCID：
Fangxing Li
✉

the School of Electrical Engineering and Automation, Hefei University of Technology, Hefei 230009, China； the University of Tennessee, Knoxville, TN 37996, USA； the School of Electrical Engineering and Automation, Wuhan University, Wuhan 430072, China

Updated：2023-01-25

DOI：10.35833/MPCE.2021.000283

OUTLINE

Abstract

In the existing multi-period robust optimization methods for the optimal power flow in radial distribution systems, the capability of distributed generators (DGs) to regulate the reactive power, the operation costs of the regulation equipment, and the current of the shunt capacitor of the cables are not considered. In this paper, a multi-period two-stage robust scheduling strategy that aims to minimize the total cost of the power supply is developed. This strategy considers the time-of-use price, the capability of the DGs to regulate the active and reactive power, the action costs of the regulation equipment, and the current of the shunt capacitors of the cables in a radial distribution system. Furthermore, the numbers of variables and constraints in the first-stage model remain constant during the iteration to enhance the computation efficiency. To solve the second-stage model, only the model of each period needs to be solved. Then, their objective values are accumulated, revealing that the computation rate using the proposed method is much higher than that of existing methods. The effectiveness of the proposed method is validated by actual 4-bus, IEEE 33-bus, and PG 69-bus distribution systems.

Keywords

Distribution system; robust optimization; mixed-integer second-order cone programming; cost of regulation equipment; coordinated optimization of active and reactive power

I. Introduction

ELECTRIC vehicles (EVs) have the potential to reduce fossil fuel dependence, environmental pollution, and greenhouse gas emissions. Therefore, EV ownership is expected to increase significantly in the next few years. A large number of EVs randomly connected to the power grid with uncoordinated or fast charging will aggravate the peak-valley difference in loads and deteriorate the safe and economic operation of distribution systems. Owing to the highly random charging times and demands for electricity, the increased presence of EVs presents a significant challenge to the optimal control of distribution systems.

Similarly, the renewable energy penetration has been increasing very rapidly, and the randomness, volatility, and anti-peak regulation of renewables pose serious threats to the real-time power balance of power grids [

1]. As the problem of renewable energy consumption in distribution systems has become increasingly urgent, multi-period rolling optimization can effectively solve the above problem [2].

Generally, the multi-period optimization of an active distribution system consists of centralized [

3], [4] and distributed [5], [6] methods. In the centralized optimization, the strategy is to relax a nonconvex model to a convex programming model using the radial operation structure of distribution systems [7]. Convex relaxation methods mainly include second-order cone programming (SOCP) [8] and semidefinite programming (SDP) [9]. The SOCP model is only applicable to balanced distribution systems. Although the SDP method is applicable to unbalanced distribution systems, its computation complexity sharply increases with the number of nodes. Considering that the method for solving the mixed-integer SOCP technique is also very mature, it is theoretically feasible to solve the optimal power flow (OPF) problem of active distribution systems with discrete variables.

In [

10]-[12], the precision, applicability, and feasibility of SOCP convex relaxation for the OPF in distribution systems are analyzed. In [13], it is pointed out that the traditional SOCP model and its convex relaxation conditions are no longer applicable when there are coaxial cables in the distribution system because the current of the shunt capacitor of the coaxial cable cannot be ignored. An OPF model of distribution systems is formulated considering the current of the shunt capacitor of the cables. A set of sufficient conditions for relaxing the model to an SOCP model that can be checked ex ante is developed.

The traditional deterministic optimization method may lead to voltage and current magnitudes outside their limits, considering the uncertainties in the intermittent distributed generators (DGs) and loads. Robust optimization is an effective way to hedge against the uncertainty [

14], [15]. In [16]-[18], a two-stage robust optimization model is formulated and solved using the conventional column and constraint generation (CCG) method. In our previous study [19], a two-stage robust optimization model that accounts for the travel distance constraints of energy storage systems (ESSs), switched capacitor reactors (SCRs), and on-load tap changers (OLTCs) is constructed with the minimum network loss as the objective function, and a fast solution method is proposed. However, when other types of objective functions are adopted, the method proposed in [19] may not be applicable. In [20], a robust day-ahead scheduling model for smart distribution systems considering demand response programs is presented. In [21], a multi-interval uncertainty-constrained robust dispatch model is proposed. The uncertainty budget is rationally divided according to the distribution probabilities to improve the over-conservativeness of traditional robust models. To address the min-max-min robust model with a mixed-integer recourse problem, a nested CCG method is adopted to quickly obtain the minimum operating cost in the worst-case scenario. In [22], a two-stage robust optimization model for planning the expansion of active distribution systems coupled with urban transportation networks is proposed. The interactions between the transportation networks and the active distribution system are explored. In [23], a novel tri-level robust planning-operation co-optimization model is proposed to determine the capacity, power, location, and scheduling strategy of distributed ESSs. In the literature, the CCG method is commonly used to solve robust optimization models.

At present, the action costs of the regulation equipment, the regulation potentials of reactive power of the DGs, and the current of the shunt capacitors of the cables are not considered in the existing robust optimization models. To this end, considering the action costs of ESSs, SCRs, and OLTCs, a robust optimization model based on the branch flow equations is developed for the active and reactive power coordinations of distribution systems with cables. A fast robust optimization method that iteratively solves on a cutting plane is proposed. The capability of the proposed method is validated by actual 4-bus, IEEE 33-bus, and PG 69-bus distribution systems.

Compared with our previous work in [

19], the main innovations of this paper are as follows.

1) The adjustment costs of the ESSs, SCRs, and OLTCs are taken into account.

2) Both the active and reactive power of all DGs are utilized, which can significantly increase regulation flexibility.

3) The distribution system model adopted is based on [

13], in which the shunt capacitance of the cables is considered. Moreover, the thermal loading constraint is different from that in [19].

4) The SOCP convex relaxation conditions can be checked ex ante using conditions C1-C5 in [

13].

5) The objective function of the robust optimization model consists of the first-stage variables, and the peak-valley price difference can be utilized with the time-of-use price.

6) A method is proposed to linearize the objective function of the second-stage model.

7) The model proposed in this paper contains more constraints of second-order cones, equalities, and inequalities than that in [

19]. As a result, the problem scale is much larger than that in [19]. Nevertheless, the proposed robust optimization method can efficiently solve the optimization problem.

The novelty and originality of this paper are as follows.

1) A two-stage multi-period mixed-integer second-order cone robust optimization model of a distribution system of cables considering the time-of-use price is developed on the basis of the branch flow equations. There are no dummy variables in the second-stage model.

2) In contrast to the CCG method, the numbers of optimization variables and constraints in the first-stage model remain constant and are less than those of the CCG method by approximately two orders.

3) For the second-stage multi-period model, the solution complexity of the second-stage model is greatly reduced compared with that of the CCG method.

4) Overall, the computation rate of the proposed method is significantly enhanced with a higher precision compared with those CCG method.

The remainder of this paper is organized as follows. In Section II, the model of a radial distribution system with cables is presented. A robust optimization model is developed in Section III. In Section IV, the solution method for the robust optimization model is presented. Finally, concluding remarks are summarized in Section V.

II. Models of Radial Distribution System with Cables

A. SOCP Model of Radial Distribution System with Cables

Figure 1 shows the radial distribution system with cables in which the susceptance $b_{l}$ cannot be ignored in general.

Fig. 1 Radial distribution system with cables.

Although the proposed method is mainly developed for distribution system with cables such as urban distribution system, the proposed method is still applicable to distribution systems with overhead lines. This is because the framework of the equivalent circuit for an overhead line is the same as that of a cable. Without loss of generality, we assume that only bus 1 is connected to the slack bus. When the line parameters in distribution systems meet conditions C1-C5 in [

13], the SOCP model of OPF for a radial distribution system can be formulated as follows on the basis of the branch flow equations [24]-[27]:

\underset{s, S, v, f, \hat{S}, \bar{v}, \bar{S}, \bar{f}}{m i n} \{\sum_{l \in 𝓁} (c^{e} (R (s_{l}), ℑ (s_{l}))) + c^{e} (P_{1}^{t})\}

(1)

s.t.

\begin{matrix} S_{l}^{t} = s_{l} + \sum_{m \in 𝓁} G_{l, m} S_{m}^{t} + z_{l} f_{l} - j (v_{u p (l)} + v_{l}) b_{l} & \forall l \in 𝓁 \end{matrix}

(2)

\begin{matrix} v_{l} = v_{u p (l)} - 2 R (z_{l}^{*} (S_{l}^{t} + j v_{u p (l)} b_{l})) + {|z_{l}|}^{2} f_{l} & \forall l \in 𝓁 \end{matrix}

(3)

\begin{matrix} f_{l} \geq \frac{{|S_{l}^{t} + j v_{u p (l)} b_{l}|}^{2}}{v_{u p (l)}} & \forall l \in 𝓁 \end{matrix}

(4)

\begin{matrix} {\hat{S}}_{l}^{t} = s_{l} + \sum_{m \in 𝓁} G_{l, m} {\hat{S}}_{m}^{t} - j ({\bar{v}}_{u p (l)} + {\bar{v}}_{l}) b_{l} & \forall l \in 𝓁 \end{matrix}

(5)

\begin{matrix} {\bar{v}}_{l} = {\bar{v}}_{u p (l)} - 2 R (z_{l}^{*} ({\hat{S}}_{l}^{t} + j {\bar{v}}_{u p (l)} b_{l})) & \forall l \in 𝓁 \end{matrix}

(6)

\begin{matrix} {\bar{S}}_{l}^{t} = s_{l} + \sum_{m \in 𝓁} G_{l, m} {\bar{S}}_{m}^{t} + z_{l} {\bar{f}}_{l} - j (v_{u p (l)} + v_{l}) b_{l} & \forall l \in 𝓁 \end{matrix}

(7)

\begin{array}{l} {\bar{f}}_{l} v_{l} \geq m a x \{{|{\hat{P}}_{l}^{b}|}^{2}, {|{\bar{P}}_{l}^{b}|}^{2}\} + m a x \{{|{\hat{Q}}_{l}^{b} - {\bar{v}}_{l} b_{l}|}^{2}, {|{\bar{Q}}_{l}^{b} - v_{l} b_{l}|}^{2}\} \\ \forall l \in 𝓁 \end{array}

(8)

\begin{array}{l} {\bar{f}}_{l} v_{u p (l)} \geq m a x \{{|{\hat{P}}_{l}^{t}|}^{2}, {|{\bar{P}}_{l}^{t}|}^{2}\} + m a x \{{|{\hat{Q}}_{l}^{t} + {\bar{v}}_{u p (l)} b_{l}|}^{2}, {|{\bar{Q}}_{l}^{t} + v_{u p (l)} b_{l}|}^{2}\} \\ \forall l \in 𝓁 \end{array}

(9)

\begin{matrix} {\bar{S}}_{l}^{b} = s_{l} + \sum_{m \in 𝓁} G_{l, m} {\bar{S}}_{m}^{t} & \forall l \in 𝓁 \end{matrix}

(10)

\begin{matrix} {\hat{S}}_{l}^{b} = s_{l} + \sum_{m \in 𝓁} G_{l, m} {\hat{S}}_{m}^{t} & \forall l \in 𝓁 \end{matrix}

(11)

\begin{matrix} v_{m i n} \leq v_{l} & \forall l \in 𝓁 \end{matrix}

(12)

\begin{matrix} {\bar{v}}_{l} \leq v_{m a x} & \forall l \in 𝓁 \end{matrix}

(13)

\begin{matrix} {|m a x \{|{\hat{P}}_{l}^{b}|, |{\bar{P}}_{l}^{b}|\} + j m a x \{|{\hat{Q}}_{l}^{b}|, |{\bar{Q}}_{l}^{b}|\}|}^{2} \leq v_{l} I_{l}^{m a x} & \forall l \in 𝓁 \end{matrix}

(14)

\begin{matrix} {|m a x \{|{\hat{P}}_{l}^{t}|, |{\bar{P}}_{l}^{t}|\} + j m a x \{|{\hat{Q}}_{l}^{t}|, |{\bar{Q}}_{l}^{t}|\}|}^{2} \leq v_{u p (l)} I_{l}^{m a x} & \forall l \in 𝓁 \end{matrix}

(15)

\begin{matrix} P_{l}^{t} \leq {\bar{P}}_{l}^{t} & \forall l \in 𝓁 \end{matrix}

(16)

\begin{matrix} Q_{l}^{t} \leq {\bar{Q}}_{l}^{t} & \forall l \in 𝓁 \end{matrix}

(17)

\begin{matrix} p_{l}^{m i n} \leq R (s_{l}) \leq p_{l}^{m a x} & \forall l \in 𝓁 \end{matrix}

(18)

\begin{matrix} q_{l}^{m i n} \leq ℑ (s_{l}) \leq q_{l}^{m a x} & \forall l \in 𝓁 \end{matrix}

(19)

Considering the regulation equipment, DGs, and different time intervals, (2) is formulated using:

\begin{matrix} P_{l, t}^{t} = p_{l, t} - P_{l, t}^{D G} - P_{l, t}^{E S S} + \sum_{m \in 𝓁} G_{l, m} P_{m, t}^{t} + r_{l} f_{l, t} & \forall l \in 𝓁 \end{matrix}

(20)

Q_{l, t}^{t} = q_{l, t} - Q_{l, t}^{D G} - Q_{l, t}^{S V C} - Q_{l, t}^{S C R} + \sum_{m \in 𝓁} G_{l, m} Q_{m, t}^{t} + x_{l} f_{l, t} - (v_{u p (l), t} + v_{l, t}) b_{l} \forall l \in 𝓁

(21)

Equation (3) can be formulated as:

\begin{matrix} v_{l, t} = v_{u p (l), t} - 2 r_{l} P_{l, t}^{t} - 2 x_{l} (Q_{l, t}^{t} + v_{u p (l), t} b_{l}) + {|z_{l}|}^{2} f_{l, t} & \forall l \in 𝓁 \end{matrix}

(22)

Equation (4) can be formulated as:

f_{l, t} v_{u p (l), t} \geq \begin{matrix} {|P_{l, t}^{t}|}^{2} + {|Q_{l, t}^{t} + v_{u p (l), t} b_{l}|}^{2} & \forall l \in 𝓁 \end{matrix}

(23)

Equation (5) can be formulated as:

\begin{matrix} {\hat{P}}_{l, t}^{t} = p_{l, t} - P_{l, t}^{D G} - P_{l, t}^{E S S} + \sum_{m \in 𝓁} G_{l, m} {\hat{P}}_{m, t}^{t} & \forall l \in 𝓁 \end{matrix}

(24)

\begin{array}{l} {\hat{Q}}_{l, t}^{t} = q_{l, t} - Q_{l, t}^{D G} - Q_{l, t}^{S V C} - Q_{l, t}^{S C R} + \sum_{m \in 𝓁} G_{l, m} {\hat{Q}}_{m, t}^{t} - ({\bar{v}}_{u p (l), t} + {\bar{v}}_{l, t}) b_{l} \\ \forall l \in 𝓁 \end{array}

(25)

Equation (6) can be formulated as:

\begin{matrix} {\bar{v}}_{l} = {\bar{v}}_{u p (l)} - 2 r_{l} {\hat{P}}_{l, t}^{t} - 2 x_{l} ({\hat{Q}}_{l, t}^{t} + {\bar{v}}_{u p (l), t} b_{l}) & \forall l \in 𝓁 \end{matrix}

(26)

Equation (7) can be formulated as:

\begin{matrix} {\bar{P}}_{l, t}^{t} = p_{l, t} - P_{l, t}^{D G} - P_{l, t}^{E S S} + \sum_{m \in 𝓁} G_{l, m} {\bar{P}}_{m, t}^{t} + r_{l} {\bar{f}}_{l, t} & \forall l \in 𝓁 \end{matrix}

(27)

{\bar{Q}}_{l, t}^{t} = q_{l, t} - Q_{l, t}^{D G} - Q_{l, t}^{S V C} - Q_{l, t}^{S C R} + \sum_{m \in 𝓁} G_{l, m} {\bar{Q}}_{m, t}^{t} + x_{l} {\bar{f}}_{l, t} - (v_{u p (l), t} + v_{l, t}) b_{l} \forall l \in 𝓁

(28)

Equation (8) can be formulated as:

\begin{matrix} \{\begin{array}{l} {\bar{f}}_{l, t} v_{l, t} \geq {|{\hat{P}}_{l, t}^{b}|}^{2} + {|{\hat{Q}}_{l, t}^{b} - {\bar{v}}_{l, t} b_{l}|}^{2} \\ {\bar{f}}_{l, t} v_{l, t} \geq {|{\hat{P}}_{l, t}^{b}|}^{2} + {|{\bar{Q}}_{l, t}^{b} - v_{l, t} b_{l}|}^{2} \\ {\bar{f}}_{l, t} v_{l, t} \geq {|{\bar{P}}_{l, t}^{b}|}^{2} + {|{\hat{Q}}_{l, t}^{b} - {\bar{v}}_{l, t} b_{l}|}^{2} \\ {\bar{f}}_{l, t} v_{l, t} \geq {|{\bar{P}}_{l, t}^{b}|}^{2} + {|{\bar{Q}}_{l, t}^{b} - v_{l, t} b_{l}|}^{2} \end{array} & \forall l \in 𝓁 \end{matrix}

(29)

Equation (9) can be formulated as:

\begin{matrix} \{\begin{array}{l} {\bar{f}}_{l, t} v_{u p (l), t} \geq {|{\hat{P}}_{l, t}^{t}|}^{2} + {|{\hat{Q}}_{l, t}^{t} + {\bar{v}}_{u p (l), t} b_{l}|}^{2} \\ {\bar{f}}_{l, t} v_{u p (l), t} \geq {|{\hat{P}}_{l, t}^{t}|}^{2} + {|{\bar{Q}}_{l, t}^{t} + v_{u p (l), t} b_{l}|}^{2} \\ {\bar{f}}_{l, t} v_{u p (l), t} \geq {|{\bar{P}}_{l, t}^{t}|}^{2} + {|{\hat{Q}}_{l, t}^{t} + {\bar{v}}_{u p (l), t} b_{l}|}^{2} \\ {\bar{f}}_{l, t} v_{u p (l), t} \geq {|{\bar{P}}_{l, t}^{t}|}^{2} + {|{\bar{Q}}_{l, t}^{t} + v_{u p (l), t} b_{l}|}^{2} \end{array} & \forall l \in 𝓁 \end{matrix}

(30)

Equation (10) can be formulated as:

\begin{matrix} {\bar{P}}_{l, t}^{b} = p_{l, t} - P_{l, t}^{D G} - P_{l, t}^{E S S} + \sum_{m \in 𝓁} G_{l, m} {\bar{P}}_{m, t}^{t} & \forall l \in 𝓁 \end{matrix}

(31)

\begin{matrix} {\bar{Q}}_{l, t}^{b} = q_{l, t} - Q_{l, t}^{D G} - Q_{l, t}^{S V C} - Q_{l, t}^{S C R} + \sum_{m \in 𝓁} G_{l, m} {\bar{Q}}_{m, t}^{t} & \forall l \in 𝓁 \end{matrix}

(32)

Equation (11) can be formulated as:

\begin{matrix} {\hat{P}}_{l, t}^{b} = p_{l, t} - P_{l, t}^{D G} - P_{l, t}^{E S S} + \sum_{m \in 𝓁} G_{l, m} {\hat{P}}_{m, t}^{t} & \forall l \in 𝓁 \end{matrix}

(33)

\begin{matrix} {\hat{Q}}_{l, t}^{b} = q_{l, t} - Q_{l, t}^{D G} - Q_{l, t}^{S V C} - Q_{l, t}^{S C R} + \sum_{m \in 𝓁} G_{l, m} {\hat{Q}}_{m, t}^{t} & \forall l \in 𝓁 \end{matrix}

(34)

Equation (14) can be formulated as:

\begin{matrix} \{\begin{array}{l} v_{l, t} I_{l}^{m a x} \geq {|{\hat{P}}_{l, t}^{b}|}^{2} + {|{\hat{Q}}_{l, t}^{b}|}^{2} \\ v_{l, t} I_{l}^{m a x} \geq {|{\hat{P}}_{l, t}^{b}|}^{2} + {|{\bar{Q}}_{l, t}^{b}|}^{2} \\ v_{l, t} I_{l}^{m a x} \geq {|{\bar{P}}_{l, t}^{b}|}^{2} + {|{\hat{Q}}_{l, t}^{b}|}^{2} \\ v_{l, t} I_{l}^{m a x} \geq {|{\bar{P}}_{l, t}^{b}|}^{2} + {|{\bar{Q}}_{l, t}^{b}|}^{2} \end{array} & \forall l \in 𝓁 \end{matrix}

(35)

Equation (15) can be formulated as:

\begin{matrix} \{\begin{array}{l} v_{u p (l), t} I_{l}^{m a x} \geq {|{\hat{P}}_{l, t}^{t}|}^{2} + {|{\hat{Q}}_{l, t}^{t}|}^{2} \\ v_{u p (l), t} I_{l}^{m a x} \geq {|{\hat{P}}_{l, t}^{t}|}^{2} + {|{\bar{Q}}_{l, t}^{t}|}^{2} \\ v_{u p (l), t} I_{l}^{m a x} \geq {|{\bar{P}}_{l, t}^{t}|}^{2} + {|{\hat{Q}}_{l, t}^{t}|}^{2} \\ v_{u p (l), t} I_{l}^{m a x} \geq {|{\bar{P}}_{l, t}^{t}|}^{2} + {|{\bar{Q}}_{l, t}^{t}|}^{2} \end{array} & \forall l \in 𝓁 \end{matrix}

(36)

B. Model of ESS

The power and energy constraints of the ESS are the same as those in [

19]. One major concern with the ESS is its limited lifecycle. Frequent charging or discharging may significantly affect the lifecycle of the ESS. To address this, the cycle limit for the ESS has been considered in this paper according to (20) in [19]. The cost of the ESS is given by:

C_{l}^{E S S} = \sum_{t = 1}^{T_{m a x}} Δ t [c_{l}^{E S S} (P_{l, t}^{c h} + P_{l, t}^{d i s}) + c^{e} (1 - η_{c h}) P_{l, t}^{c h} + c^{e} (\frac{1}{η_{d i s}} - 1) P_{l, t}^{d i s}]

(37)

C. Model of SVC

Q_{l}^{S V C, m i n} \leq Q_{l}^{S V C} \leq Q_{l}^{S V C, m a x}

(38)

D. Model of SCR

The operation constraints of the SCR are the same as those in [

19]. The action cost of the SCR

C_{l}^{S C R}

is given by:

C_{l}^{S C R} = c_{l}^{S C R} \sum_{t = 1}^{T_{m a x} - 1} X_{l, t}

(39)

E. Model of OLTC

The schematic of an OLTC is shown in Fig. 2.

Fig. 2 Schematic of an OLTC.

The branch impedance can be thought of as a cable with zero shunt capacitance. The operation constraints of the OLTC are the same as those in [

19]. The action cost of the OLTC

C^{O L T C}

is given by:

C^{O L T C} = c^{O L T C} \sum_{t = 1}^{T_{m a x} - 1} O_{t}

(40)

F. Model of DG

In this paper, only DGs with a full-capacity converter interface such as a photovoltaic power plant or permanent magnet synchronous generator (PMSG) wind turbine are considered. However, the proposed method can also be applied to distribution systems with other types of DGs such as doubly-fed induction generators. It is only necessary to slightly modify (41)-(43).

1)　Constraint on Power Factor

P_{l, t}^{D G} t a n θ_{l}^{m i n} \leq Q_{l, t}^{D G} \leq P_{l, t}^{D G} t a n θ_{l}^{m a x}

(41)

When the energy radiated by the sun and the wind speed are very low, a DG may be cut off. Further, the power factor of the DG should be in a reasonable range to enhance the efficiency. To address this, the reactive power is bounded as in (41).

2)　Constraint on Active Power

0 \leq P_{l, t}^{D G} \leq g_{l, t}

(42)

The injected active power of a DG is greater than zero and less than its predicted value.

3)　Constraint on Capacity

\sqrt[]{{(P_{l, t}^{D G})}^{2} + {(Q_{l, t}^{D G})}^{2}} \leq S_{l, t}^{D G}

(43)

The active power and reactive power of a DG are bounded by the capacity of the converter.

G. Objective Function

f_{o b j} = m i n \{\sum_{t = 1}^{T_{m a x}} c^{e} P_{1, t}^{t} Δ t + \sum_{l \in B^{E S S}} C_{l}^{E S S} + C^{O L T C} + \sum_{l \in B^{S C R}} C_{l}^{S C R}\}

(44)

The objective function is composed of the cost of purchasing electricity and the operation costs of the ESS, OLTC, and SCR.

III. Robust Optimization Model

A. Deterministic Model

Because the OLTC, SVR, and ESS cannot be frequently adjusted, two-stage optimization strategies are utilized in the proposed method. The first-stage variables include the charging power and discharging power, the energy stored in the ESS, the operation groups of the SCR, the tap position of the OLTC, and other discrete variables. All second-stage variables are continuous. They include the branch flow variables and their upper and lower bounds, the reactive power of the SVC, and the active and reactive powers of the DG. In every optimization time interval, the first-stage variables cannot be changed when they are set. However, the second-stage variables can be flexibly changed in response to actual operation conditions. However, the compensated reactive power of the SCR is set as a second-stage variable because it is proportional to the square of the voltage. The dummy variables associated with the second-stage variables, which are used to linearize the model, are also set as second-stage variables.

For a clear presentation, the deterministic optimization model can be written in a compact form as:

M a s t e r p r o b l e m (M P) : f_{m p} = \underset{{ψ_{t}}, {φ_{t}}}{m i n} f_{o b j}

(45)

s.t.

ψ_{t} \in Ψ_{t}

(46)

\sum_{t = 1}^{T_{m a x}} A_{t}^{T} ψ_{t} \leq b

(47)

B_{t}^{T} ψ_{t} \leq b_{1 t}

(48)

J_{t} φ_{t} \leq b_{2 t}

(49)

L_{t} ψ_{t} + N_{t} φ_{t} \leq b_{3 t}

(50)

D_{i, t}^{p} ψ_{t} + E_{i, t}^{p} φ_{t} = d_{t}^{p} i = 1,2, \dots, 5

(51)

D_{i, t}^{q} ψ_{t} + E_{i, t}^{q} φ_{t} = d_{t}^{q} i = 1,2, \dots, 5

(52)

H_{l, t} φ_{t} \leq g_{l, t} \forall l \in B^{D G}

(53)

‖K_{n t} φ_{t}‖ \leq h_{n t}^{T} φ_{t}

(54)

Equation (45) minimizes the total cost. Equation (46) represents a feasible set for $ψ_{t}$ . Equation (47) describes the coupling of $ψ_{t}$ between different intervals. Equation (48) describes the inequality for $ψ_{t}$ in each interval. Equation (49) summarizes the inequality for $φ_{t}$ in each interval. Equation (50) denotes the coupling between $ψ_{t}$ and $φ_{t}$ in each interval. Equations (51) and (52) represent the active and reactive power flows in each interval corresponding to (2), (5), (7), (10), and (11). Equation (53) represents the constraints on the active power of the DGs. Equation (54) represents all second-order cone constraints associated with all branches and DGs.

B. Robust Optimization Model

Because photovoltaics, wind power, and loads have large uncertainties, a robust optimization scheduling strategy is formulated as:

\underset{{ψ_{t}}}{m i n} \underset{\begin{array}{l} d_{t} \in D_{t} \\ g_{t} \in G_{t} \end{array}}{m a x} \underset{{φ_{t}}}{m i n} f_{o b j}

(55)

s.t.

(46)-(54)

D_{t} = \{(d_{t}^{p}, d_{t}^{q}) |d_{t}^{m i n} \leq d_{t}^{p} \leq d_{t}^{m a x}, d_{t}^{m i n} \circ t a n θ_{L} \leq d_{t}^{q} \leq d_{t}^{m a x} \circ t a n θ_{L}\}

(56)

G_{t} = \{g_{l, t} |g_{l, t}^{m i n} \leq g_{l, t} \leq g_{l, t}^{m a x}\} \forall l \in B^{D G}

(57)

In (55), the first-stage solutions of $ψ_{t}$ try to minimize the total costs for the worst-case scenario. The inner max-min bilevel solutions seek to determine the worst-case scenario.

In the uncertainty set of the loads and DGs in (56) and (57), the correlations between the wind, photovoltaics, and loads are not considered in this paper, which is the same as in [

3], [16], and [17]. However, these correlations can be considered by the addition of uncertainty budgets. If they are considered, the proposed method is still applicable, and the optimization results can be less conservative.

In the worst case of uncertainty, the system operates under more severe conditions compared with those of a normal OPF problem. However, the exactness of convex relaxation can be guaranteed if and only if the problem in (55)-(57) is feasible regardless of the worst case of uncertainty when the line parameters in distribution systems meet conditions C1-C5 in [

13]. Moreover, conditions C1-C5 can be checked ex ante and are mild. However, the exactness can be checked by calculating the gaps in convex relaxation after the SOCP problem is solved if conditions C1-C5 are not satisfied. If one of the gaps is larger than some tolerance, it indicates that inexactness occurs. Then, the method in [28] can be used to formulate a sequential SOCP problem. Although its computation rate is slower than that for an SOCP problem, a (near) global optimal solution can still be recovered with a relatively high solution efficiency.

IV. Solution Method for Robust Optimization Model

A. Subproblem (SP)

Equations (45)-(54) present the MP, while (58)-(64) present the SP by substituting the solutions for the first-stage variables into (20), (21), (24), (25), (27), (28), and (31)-(34).

S P : L (ψ_{1}^{*}, ψ_{2}^{*}, \dots, ψ_{T_{m a x}}^{*}) = \underset{\begin{array}{l} d_{t} \in D_{t} \\ g_{t} \in G_{t} \end{array}}{m a x} m i n \sum_{t = 1}^{T_{m a x}} c_{t}^{T} y_{t}

(58)

R_{t} y_{t} = b_{1 t} \forall t, π_{1 t}

(59)

M_{t} y_{t} \leq b_{2 t} \forall t, π_{2 t}

(60)

W_{i, t}^{p} φ_{t} = d_{t}^{p} + d_{0 t}^{p} i = 1,2, \dots, 5, \forall t, π_{3 t}^{(i)}

(61)

W_{i, t}^{q} φ_{t} = d_{t}^{q} i = 1,2, \dots, 5, \forall t, π_{4 t}^{(i)}

(62)

I_{l, t} y_{t} \leq g_{l, t} \forall l \in B^{D G}, \forall t, π_{5 t}^{(l)}

(63)

‖K_{n t} y_{t}‖ \leq h_{n t}^{T} y_{t} \forall t, (π_{6 n t}, π_{7 n t})

(64)

Formulas (58)-(64) can be transformed into (65)-(71) using dual theory.

S P : \sum_{t = 1}^{T_{m a x}} \underset{d_{t}, g_{t}, π_{1 t}, π_{2 t}, π_{3 t}, π_{4 t}, π_{5 t}, π_{6 n t}, π_{7 n t}}{m a x} \{b_{1 t}^{T} π_{1 t} + b_{2 t}^{T} π_{2 t} + {(d_{t}^{p} + d_{0 t}^{p})}^{T} \cdot \sum_{i = 1}^{5} π_{3 t}^{(i)} + {(d_{t}^{q})}^{T} \sum_{i = 1}^{5} π_{4 t}^{(i)} + \sum_{l \in B^{D G}} g_{l, t} π_{5 t}^{(l)}\}

(65)

s.t.

‖π_{6 n t}‖ \leq π_{7 n t} \forall t

(66)

\begin{array}{l} R_{t}^{T} π_{1 t} + M_{t}^{T} π_{2 t} + \sum_{i = 1}^{5} {(W_{i, t}^{p})}^{T} π_{3 t}^{(i)} + \sum_{i = 1}^{5} {(W_{i, t}^{q})}^{T} π_{4 t}^{(i)} + \\ \sum_{k = 1}^{n_{0}} I_{k, t}^{T} π_{5 t}^{(i)} + \sum_{n} (K_{n t}^{T} π_{6 n t} + h_{n t} π_{7 n t}) = c_{t} \forall t \end{array}

(67)

π_{2 t} \leq 0 \forall t

(68)

d_{t}^{m i n} \leq d_{t}^{p} \leq d_{t}^{m a x} \forall t

(69)

d_{t}^{m i n} \circ t a n θ_{L} \leq d_{t}^{q} \leq d_{t}^{m a x} \circ t a n θ_{L} \forall t

(70)

g_{l, t}^{m i n} \leq g_{l, t} \leq g_{l, t}^{m a x} \forall t, \forall l \in B^{D G}

(71)

The worst-case scenario must take place at the extreme point of each component of the uncertainty set. This is because the components of the uncertainty set are independent of each other [

3], [16], [17]. Let

π_{3 t} = \sum_{i = 1}^{5} π_{3 t}^{(i)}

and

γ_{t, n}^{p} = π_{3 t, n} δ_{t, n}

; then,

{(d_{t}^{p})}^{T} π_{3 t}

can be linearized as:

\{\begin{array}{l} {(d_{t}^{p})}^{T} π_{3 t} = \sum_{n} [d_{t, n}^{p, m i n} π_{3 t, n} + (d_{t, n}^{p, m a x} - d_{t, n}^{p, m i n}) γ_{t, n}^{p}] \\ - M δ_{t, n} \leq γ_{t, n}^{p} \leq M δ_{t, n} \\ - M (1 - δ_{t, n}) \leq γ_{t, n}^{p} - π_{3 t, n} \leq M (1 - δ_{t, n}) \\ δ_{t, n} \in \{0,1\} \end{array} \forall t, \forall n

(72)

Similarly, let $π_{4 t} = \sum_{i = 1}^{5} π_{4 t}^{(i)}$ and $γ_{t, n}^{q} = π_{4 t} δ_{t, n}$ ; then, ${(d_{t}^{q})}^{T} π_{4 t}$ can be linearized as:

\{\begin{array}{l} {(d_{t}^{q})}^{T} π_{4 t} = \sum_{n} [d_{t, n}^{p, m i n} t a n (θ_{L, n}) π_{4 t, n} + \\ (d_{t, n}^{p, m a x} - d_{t, n}^{p, m i n}) t a n (θ_{L, n}) γ_{t, n}^{q}] \\ - M δ_{t, n} \leq γ_{t, n}^{q} \leq M δ_{t, n} \\ - M (1 - δ_{t, n}) \leq γ_{t, n}^{q} - π_{4 t, n} \leq M (1 - δ_{t, n}) \end{array} \forall t, \forall n

(73)

Similarly, let $γ_{l, t}^{D G} = λ_{t} π_{5 t}^{(l)}$ for the same type of renewable DG such as a wind turbine, then, $g_{l, t} π_{5 t}^{(l)}$ can be linearized as:

\{\begin{array}{l} g_{l, t} π_{5 t}^{(l)} = g_{l, t}^{m i n} π_{5 t}^{(l)} + (g_{l, t}^{m a x} - g_{l, t}^{m i n}) γ_{l, t}^{D G} \\ - M λ_{t} \leq γ_{l, t}^{D G} \leq M λ_{t} \\ - M (1 - λ_{t}) \leq γ_{l, t}^{D G} - π_{5 t}^{(l)} \leq M (1 - λ_{t}) \\ λ_{t} \in \{0,1\} \end{array} \forall t, \forall l \in B^{D G}

(74)

It is worth noting that for the same type of DG, the uncertainties in their outputs are not independent of each other. The same statistical law must be followed. Therefore, $λ_{t}$ of each DG is identical during the same time period for the same type of DG.

B. Solving Steps

Step 1: let $k = 1$ , $U B = + \infty$ , and $L B = - \infty$ .

Step 2: obtain the optimal $(ψ_{t}^{*}, φ_{t}^{*}, f_{m p}^{*})$ by solving the MP. Compute $f_{φ}^{*} = \sum_{t = 1}^{T_{m a x}} (c^{e} P_{1, t}^{t} Δ t)$ , and update $L B = f_{φ}^{*}$ .

Step 3: fix ${ψ_{t}^{*}}$ . Solve the SP for the single-period model to obtain the worst-case scenario ${d_{t}^{*}}$ and ${g_{t}^{*}}$ . Calculate $f_{s p}^{*} (ψ_{t}^{*})$ . Compute $U B = f_{s p}^{*} (ψ_{t}^{*})$ . If $| L B - U B | < ε$ or $k > k_{m a x}$ , terminate the program, and output $ψ_{t}^{*}$ . Otherwise, ${d_{t}^{*}}$ and ${g_{t}^{*}}$ are substituted into the MP.

Step 4: update $k = k + 1$ , and go to Step 2.

Generally, $ε$ can be set to be 0.0001-0.01 (in units of thousands of dollars), while $k_{m a x}$ can be set to be 5-10.

C. Configurations of Three Test Systems

The actual 4-bus, IEEE 33-bus, and PG 69-bus distribution systems with uncertain wind power generation and loads are used as the simulation cases. The topology of the actual 4-bus distribution system is shown in Fig. 3, while the topologies of the IEEE 33-bus and PG 69-bus distribution systems are shown in Figs. 4 and 5, respectively.

Fig. 3 Topology of actual 4-bus distribution system.

Fig. 4 Topology of IEEE 33-bus distribution system.

Fig. 5 Topology of PG 69-bus distribution system.

All programs are developed using MATLAB R2018a. The mixed-integer SOCP toolbox of MOSEK 9.1.4 is used to solve the MP and SP. The computer is equipped with an Intel Xeon i5-E5640 CPU with four cores and eight threads running at 2.67 GHz and 40 GB of memory. To reduce memory usage, sparse matrices are compressed and stored.

The parameters of the actual 4-bus distribution system are the same as those in Table I of [

13], except that the length of both cables is 20 km. The shunt capacitances of the cables in the IEEE 33-bus and PG 69-bus distribution systems are listed in Tables I and II, respectively.

TABLE I Shunt Capacitances of IEEE 33-bus Distribution System

Branch	Capacitance (F)	Branch	Capacitance (F)	Branch	Capacitance (F)
1	0	12	2.0740×10^-7	23	5.1650×10^-7
2	7.8740×10^-8	13	1.9350×10^-6	24	1.1880×10^-6
3	4.2067×10^-7	14	1.1943×10^-6	25	1.1746×10^-6
4	3.1228×10^-7	15	8.8120×10^-7	26	1.7323×10^-7
5	3.2518×10^-7	16	9.1305×10^-7	27	2.4242×10^-7
6	1.1844×10^-6	17	2.8832×10^-6	28	1.5642×10^-6
7	1.0367×10^-6	18	9.6163×10^-7	29	1.1737×10^-6
8	3.9387×10^-7	19	2.6219×10^-7	30	4.3307×10^-7
9	1.2397×10^-6	20	2.2707×10^-6	31	1.6133×10^-6
10	1.2397×10^-6	21	8.0147×10^-7	32	6.0630×10^-7
11	1.0890×10^-7	22	1.5703×10^-6	33	8.9830×10^-7

TABLE II Shunt Capacitances of PG 69-bus Distribution System

Branch	Capacitance (F)	Branch	Capacitance (F)	Branch	Capacitance (F)
1	0	24	1.9182×10^-7	47	4.4061×10^-7
2	2.0104×10^-9	25	4.5987×10^-7	48	1.6854×10^-7
3	2.0104×10^-9	26	1.7105×10^-7	49	1.9635×10^-7
4	6.0311×10^-9	27	9.5828×10^-8	50	4.3307×10^-7
5	4.9254×10^-8	28	1.8093×10^-8	51	8.3096×10^-8
6	3.1228×10^-7	29	2.6219×10^-7	52	1.2364×10^-7
7	3.2518×10^-7	30	2.2030×10^-7	53	6.0630×10^-7
8	7.8740×10^-8	31	3.8867×10^-8	54	8.8825×10^-7
9	4.2050×10^-8	32	1.9434×10^-7	55	1.0236×10^-7
10	4.5351×10^-7	33	4.7177×10^-7	56	2.3454×10^-9
11	1.1576×10^-7	34	9.4588×10^-7	57	4.0945×10^-7
12	3.9387×10^-7	35	7.8287×10^-7	58	2.6805×10^-9
13	5.6961×10^-7	36	1.4073×10^-8	59	1.8093×10^-8
14	5.7798×10^-7	37	3.4897×10^-7	60	2.6219×10^-7
15	5.8569×10^-7	38	1.1880×10^-6	61	2.0606×10^-7
16	1.0890×10^-7	39	3.3691×10^-7	62	5.9474×10^-7
17	2.0740×10^-7	40	7.9242×10^-8	63	3.5182×10^-9
18	2.6805×10^-9	41	1.8663×10^-7	64	1.4255×10^-6
19	1.8144×10^-7	42	1.4843×10^-7	65	6.0698×10^-7
20	1.1660×10^-7	43	1.7323×10^-7	66	8.0080×10^-8
21	1.8914×10^-7	44	2.4242×10^-7	67	1.9434×10^-8
22	7.7064×10^-9	45	2.4007×10^-7	68	2.3002×10^-7
23	8.8122×10^-8	46	8.9412×10^-7	69	2.0104×10^-9

The base voltages of the three systems $V_{b}$ are chosen to be 24.9, 12.66, and 12.66 kV, respectively. $S_{b}$ is chosen to be 5, 10, and 10 MVA, respectively. There is one OLTC transformer connected to the root node for the three systems. The impedance of the transformer is $(0.02 + j 0.105)$ p.u.. The minimum and maximum turn ratios of the OLTC are 0.94 and 1.06, respectively. The step size of the turn ratio is 0.01. The voltage bound on each bus is [0.9, 1.1]p.u.. The root node is taken as the slack node whose voltage is fixed at 1.0 p.u.. The current limits for each branch are 120 A for the actual 4-bus system and 400 A for the IEEE 33-bus and PG 69-bus distribution systems. The value M in the Big M method is set to be 100. $ε$ is set to be 0.0001. The maximum number of iteration is set to be 7. Further, a number greater than the default value for the relative gap should be chosen for large-scale problems so that the program does not stall. The optimization period is 00:00-24:00 with 1-hour interval. The lifecycles of the ESS in the three systems are set to be 3.

There is one SCR connected to node 1 for the actual 4-bus system, one SCR connected to nodes 3 and 6 for the IEEE 33-bus system, and one SCR connected to nodes 19, 36, 41, 53, and 64 for the PG 69-bus system. The capacities of the SCRs are $[- 0.3,0.3]$ , $[- 0.6,0.6]$ , and $[- 0.6,0.6]$ Mvar, while the step sizes are 0.05, 0.1, and 0.1 Mvar for the three systems. The maximum travel distances of the OLTC and SCR are 24.

There is one SVC connected to node 3 for the actual 4-bus system, one SVC connected to node 18 for the IEEE 33-bus system, and one SVC connected to nodes 3 and 11 for the PG 69-bus system. The capacities of the SVCs in the three distribution systems are $[- 0.15,0.15]$ , $[- 0.5,0.5]$ , and $[- 0.5,0.5]$ Mvar.

There is one PMSG wind turbine connected to node 3 for the actual 4-bus system, one PMSG connected to nodes 13, 21, 24, and 31 for the IEEE 33-bus system, and one PMSG connected to nodes 19, 41, 54, 56, and 66 for the PG 69-bus system. The capacities of the PMSGs in the three distribution systems are 5, 0.4, and 0.3 MVA.

There is one ESS connected to node 1 for the actual 4-bus system, one ESS connected to nodes 17 and 33 for the IEEE 33-bus system, and one ESS connected to nodes 2 and 12 for the PG 69-bus system.

The capacity of the ESS in the actual 4-bus system is 1.5 MWh. The bound on the quantity of electric charge is [0.15, 1.5]MWh. The maximum charging power and discharging power are both 150 kW. The capacity of the ESS connected to node 17 in the IEEE 33-bus system and node 2 in the PG 69-bus system is 1.5 MWh. The bound on the quantity of electric charge is [0.15, 1.5]MWh. The maximum charging and discharging power are both 300 kW. The capacity of the ESS connected to node 33 in the IEEE 33-bus system and node 12 in the PG 69-bus system is 0.5 MWh. The bound on the quantity of electric charge is [0.05, 0.5]MWh. The maximum charging power and discharging power are both 100 kW.

The charging and discharging efficiencies of each ESS are 0.9. The maximum number of cycles of the ESS is set to be 3. The action costs of the OLTC, SCRs, and ESS, i.e., $c^{O L T C}$ , $c_{}^{S C R}$ , and $c_{}^{E S S}$ , are set to be 80, 40, and 50 $/MWh, respectively.

The electricity prices, normalized predicted loads, and wind power for each time period are listed in Table III. The maximum predicted load at nodes 2-4 for the actual 4-bus system is 1 MW. The power factor is 0.95. The maximum predicted power of the wind turbines in the actul 4-bus, IEEE 33-bus, and PG 69-bus systems is 2.5, 0.25, and 0.25 MW, respectively.

TABLE III Electricity Prices, Normalized Predicted Loads, and Wind Power

Time	Price($/MWh)	Normalized predicted load (%)	Wind power (%)	Time	Price ($/MWh)	Normalized predicted load (%)	Wind power (%)
1	50	65.8	82.7	13	78	80.0	9.2
2	38	63.2	68.7	14	85	75.3	1.3
3	39	62.1	85.3	15	100	83.2	2.0
4	40	62.6	94.6	16	82	84.2	0.0
5	46	62.9	100.0	17	70	84.7	3.9
6	45	63.6	91.2	18	115	90.5	9.7
7	145	70.5	89.1	19	160	100.0	36.2
8	150	75.3	79.8	20	200	95.8	45.9
9	64	77.9	75.4	21	220	93.7	36.4
10	60	84.2	48.2	22	210	89.5	43.7
11	64	85.3	29.0	23	60	80.0	46.5
12	75	84.7	21.2	24	40	72.1	33.7

D. Simulation Results

For the three systems, the optimization results for the turn ratio of the OLTC at different time periods and the prediction errors are all the same. For the actual 4-bus system, the optimization results for the compensated reactive power for the SCR during different time periods and the prediction errors are the same, all of which are -0.3 Mvar. For the IEEE 33-bus and PG 69-bus systems, the optimization results for the compensated reactive power for the SCR during different time periods and the prediction errors are the same, all of which are 0.6 Mvar.

For the actual 4-bus system, the charging power and discharging power during different time periods are shown in Fig. 6. When the electricity price is at its minimum during the time period of 02:00-04:00, the ESS is charged with the maximum power. When the electricity price is at its maximum during the time period of 20:00-22:00, the ESS is discharged with the maximum power.

Fig. 6 Charging power and discharging power of ESS of actual 4-bus system.

The charging and discharging states of the ESS are shown in Fig. 7. It is observed that the charging and discharging states are equal to 1 and 0, respectively, when the charging power is high. Moreover, the charging and discharging states are equal to 0 and 1, respectively, when the discharging power is high. The sum of the charging and discharging states is no more than 1.

Fig. 7 Charging and discharging states of ESS of actual 4-bus system.

The maximum gaps in conic relaxation are listed in Table IV. As can be observed, the gap in each case is less than 5×10^-6, which implies that conic relaxation is the same as the original nonconvex model. Therefore, the precision of the SOCP model is sufficiently high.

TABLE IV The Maximum Gaps in Conic Relaxation

$ζ$	The maximum gap
$ζ$	Actual 4-bus system	IEEE 33-bus system	PG 69-bus system
0.1	2.6482×10^-8	1.1511×10^-8	1.6303×10^-6
0.2	8.7915×10^-8	1.1682×10^-8	1.8231×10^-6
0.3	5.0155×10^-8	1.1737×10^-8	3.9148×10^-6
0.4	1.4872×10^-7	1.1835×10^-8	2.2381×10^-6
0.5	3.6849×10^-8	1.1859×10^-8	4.0481×10^-6
0.6	6.0642×10^-8	1.1916×10^-8	3.7460×10^-6

E. Comparison of Computational Performance of Proposed and CCG Methods

The total cost of actual 4-bus system, computational complexity of actual 4-bus system, total cost of IEEE 33-bus system, and computational complexity of IEEE 33-bus system [

19] are listed in Tables V-VIII, respectively.

TABLE V Total Cost of Actual 4-bus System

$ζ$	Total cost (10³ $)
	Improved CCG method		Proposed method
	MP	SP	MP	SP
0.1	4.0100	3.9470	4.0100	3.9470
0.2	4.8264	4.7635	4.8264	4.7635
0.3	5.6484	5.5854	5.6484	5.5854
0.4	6.4760	6.4130	6.4760	6.4130
0.5	7.3094	7.2464	7.3094	7.2464
0.6	8.1489	8.0859	8.1489	8.0859

As can be observed, the objective function values of the proposed method fit those of the improved CCG method very well for different prediction errors. Thus, the precision of the proposed method is relatively high. However, the computation rate of the proposed method is faster than that of the improved CCG method for all cases, except when the prediction error is 0.6 for the actual 4-bus system. Furthermore, the proposed method converges within only two iterations in all cases. Moreover, the values of the objective function progressively increase as the prediction error increases. That is, the tariff in the worst-case scenario increases when the uncertainties in the loads and DG outputs increase. This is obvious in practice. Therefore, the optimization results are consistent with the actual situation.

TABLE VI Computational Complexity of Actual 4-bus System

$ζ$	Improved CCG method		Proposed method
$ζ$	Iteration	Time (s)	Iteration	Time (s)
0.1	2	53.783	2	43.807
0.2	2	51.381	2	37.475
0.3	2	58.444	2	34.498
0.4	2	58.405	2	43.756
0.5	2	55.162	2	44.109
0.6	2	234.239	2	268.439

TABLE VII Total Cost of IEEE 33-bus System

$ζ$	Total costs (10³ $)
	Improved CCG method		Proposed method
	MP	SP	MP	SP
0.1	6.3069	6.1391	6.3069	6.1391
0.2	7.0422	6.8743	7.0422	6.8743
0.3	7.9154	7.7475	7.9154	7.7475
0.4	9.0912	8.9233	9.0912	8.9233
0.5	10.0154	9.8475	10.0156	9.8475
0.6	11.0535	10.8798	11.0536	10.8798

TABLE VIII Computational Complexity of IEEE 33-bus System

$ζ$	Improved CCG method		Proposed method
$ζ$	Iteration	Time (s)	Iteration	Time (s)
0.1	2	222.728	2	155.649
0.2	2	220.986	2	154.633
0.3	2	222.238	2	146.274
0.4	2	221.467	2	153.368
0.5	2	221.942	2	154.222
0.6	2	247.540	2	149.922

The worst-case scenario generated in the last iteration of the SP using the proposed method for actual 4-bus system is shown in Fig. 8 when $ζ$ is 0.2. As can be observed, the total load of the worst-case scenario is 1.2 times that of the predicted one, whereas the total wind power in the worst-case scenario is 0.8 times that of the predicted one.

Fig. 8 Worst-case scenario generated in the last iteration of SP using proposed method for actual 4-bus system.

Using the proposed method with the actual 4-bus system, the output active and reactive power of the PMSG wind turbine and the injected reactive power of the SVC in the last iteration of the MP are shown in Figs. 9-11, respectively, when the prediction error is 0.2.

Fig. 9 Active power of PMSG wind turbine in the last iteration of MP using proposed method for actual 4-bus system.

Fig. 10 Reactive power of PMSG wind turbine in the last iteration of MP using proposed method for actual 4-bus system.

Fig. 11 Injected reactive power of SVC in the last iteration of MP using proposed method for actual 4-bus system.

As can be observed, the output reactive power of the PMSG wind turbine and SVC is always negative. This is because the lengths of the cables are long, and the capacity of the PMSG wind turbine is large. As a result, the reactive power injected from the shunt capacitors of the cables and the active power from the PMSG wind turbine are large. Furthermore, the reactive power absorbed by the converter is at its maximum (minimum) when the output active power of the PMSG wind turbine is at its maximum (minimum). This is because the power factor angle of the PMSG wind turbine is set to be within $[- π / 4, π / 4]$ . Moreover, the reactive power absorbed by the SVC is low (high) when the absorbed reactive power of the PMSG wind turbine is high (low). The SVC cooperates with the converter of the PMSG wind turbine to regulate the voltage and reduce the losses.

Figure 12 shows the iterations of IEEE 33-bus distribution system when the prediction error is 0.2. As can be observed, the proposed method only requires two iterations to converge, where the upper bound remains constant and the lower bound is increased. The program stops until the gap between the upper and lower bounds is smaller than some given value.

Fig. 12 Iterations of IEEE 33-bus distribution system.

The maximum and minimum voltages and the maximum current of the proposed method for the IEEE 33-bus distribution system are shown in Figs. 13 and 14, respectively, where every load is multiplied by 1.7 and the prediction error is 0.2. As can be observed, all voltages and currents are within the rated ranges. However, the minimum voltage during the period of 18:00-20:00 when the load is the maximum is very close to the lower limit of 0.9 p.u.. Therefore, the voltage rather than the current is a key limiting factor for the safe operation of this distribution system.

Fig. 13 The maximum and minimum voltages of IEEE 33-bus distribution system.

Fig. 14 The maximum current of IEEE 33-bus distribution system.

The simulation results obtained with the proposed method and the CCG method in [

8] are summarized in Tables IX and X, respectively. Clearly, the objective function values for the proposed method are lower than those of the CCG method for different prediction errors. This is because the scale of the CCG method is much larger than that of the proposed method. For the same relative dual gap, the precision of the former is lower than that of the latter. The CCG method requires approximately 6095-8010 s and seven iterations before the program stops. However, less than 650 s, only two iterations are needed to reach the convergence using the proposed method. Furthermore, a large amount of memory can be saved.

TABLE IX Total Cost of PG 69-bus System

$ζ$	Total costs (10³ $)
	CCG method		Proposed method
	MP	SP	MP	SP
0.1	6.8862	6.7195	6.8394	6.6743
0.2	7.6412	7.4762	7.5821	7.4170
0.3	8.5639	8.3920	8.3356	8.1705
0.4	9.5894	9.4224	9.2259	9.0609
0.5	9.9721	9.8074	9.9080	9.7429
0.6	11.4549	11.2901	10.7125	10.5474

TABLE X Computational Complexity of PG 69-bus System

$ζ$	CCG method		Proposed method
$ζ$	Iteration	Time (s)	Iteration	Time (s)
0.1	7	6095.892	2	633.078
0.2	7	7730.409	2	638.009
0.3	7	7812.867	2	633.229
0.4	7	7945.767	2	620.362
0.5	7	7998.384	2	632.997
0.6	7	8010.235	2	647.081

Fast robust optimization is the main contribution of this paper. The reason why the proposed method is faster than the well-known CCG method is as follows. In contrast to the CCG method, the increase in the numbers of variables and constraints is not required to solve the first-stage model using the proposed method. Further, only a model of each single period needs to be simulated to solve the second-stage multi-period model. Consequently, the computation rate is significantly enhanced.

V. Conclusion

Generally, the shunt capacitance of a coaxial cable cannot be ignored. In this paper, a robust programming strategy for active and reactive power coordination based on the branch flow equations of a radial distribution system with cables is developed considering the action costs of regulation equipment and the regulation capability with DGs. The proposed method aims to find a robust optimal solution that can hedge against any possible realization with the uncertainties in the load, wind, or photovoltaic power outputs. Then, a fast solution method is formulated.

However, the computation rate is crucial for online rolling optimization of large-scale distribution systems. The time required to solve the MP in the CCG method becomes increasingly large during the iteration when many iterations are performed to reach the convergence for a large-scale distribution sytem. To address this issue, a fast robust optimization method is proposed in this paper. The numbers of constraints and variables for the MP remain constant during the iteration. Further, the SP only needs to be solved for each time period. Then, their objective function values are accumulated, and the worst-case scenarios of each time period are concatenated. Therefore, the solution complexity is significantly reduced. Consequently, the computation rate is much higher than that of the CCG method. The precision of the optimization results is also improved, and the amount of required computer storage space is reduced. Specifically, the simulation results of the PG 69-bus system indicate that the computation rate is enhanced by approximately one order of magnitude.

Whether the proposed method is valid for other types of uncertainty sets such as irregular and nonconvex uncertainty sets is a topic worthy of studying in the future. A comparison of the results with the practical hardware in loop models to validate the capabilities of the proposed method is another area of future work.

Nomenclature

Symbol	——	Definition
A.	——	Indices
$l$ , m	——	Index of branches other than the slack node and branches whose downstream branch is $l$
$t$	——	A subscript indicating time period
B.	——	Sets
$l$	——	Set of branches
$Ψ_{t}$	——	Set of discrete variables
$B^{E S S}$ , $B^{S C R}$ , $B^{D G}$	——	Sets of buses with energy storage system (ESS), switched capacitor reactor (SCR), and distributed generator (DG)
$D_{t}$	——	Set of active and reactive power of load
$G_{t}$	——	Set of active power of DG
C.	——	Parameters
$Δ t$	——	Scheduling interval
$η_{c h}$ , $η_{d i s}$	——	Charging and discharging efficiencies of ESS
$θ_{L, n}$	——	Power factor angle of load on node n
$θ_{L}$	——	Predicted power factor angle of load
$θ_{l}^{m i n}$ , $θ_{l}^{m a x}$	——	The minimum and maximum power factor angles of DG connected to branch l
$ζ$	——	Prediction error
$ε$	——	Convergence tolerance
$c^{e}$	——	Price for the main grid power
$c_{t}$	——	Cost matrix of the second-stage variables
$c^{E S S}$ , $c^{O L T C}$ , $c^{S C R}$	——	Action prices for ESS, on-load tap changers (OLTC), and SCR
$d_{t}^{m i n}$ , $d_{t}^{m a x}$	——	The minimum and maximum active and reactive power vectors of load during time period t
$d_{t, n}^{p, m i n}$ , $d_{t, n}^{p, m a x}$	——	The minimum and maximum active power vectors of load at node n during time period t
$d_{t}^{p}$ , $d_{t}^{q}$	——	Active and reactive power vectors of load during time period t
d_t	——	Active power vector of load
f_obj	——	Objective function of robust optimization problem
$f_{φ}^{*}$	——	Optimal cost of electricity
$f_{s p}^{*}$	——	Optimal objective function value of subproblem
G	——	Adjacency matrix of oriented graph of distribution system
$g_{l, t}^{m i n}$ , $g_{l, t}^{m a x}$	——	The minimum and maximum active power of DG on branch l during time period t
$G_{l, m}$	——	Adjacency matrix of the oriented graph of distribution system, i.e., $G_{l, m}$ is defined for $(l, m) \in 𝓁$ , and $G_{l, m} = 1$ if $l = u p (m)$ and 0; otherwise, diagonal elements are zero
$I_{l}^{m a x}$	——	Upper bound of current square on branch l
$k_{t}$	——	Transformation ratio
M	——	A big number
$n_{m a x}$	——	The maximum number of iterations
$p_{l}^{m i n}$ , $p_{l}^{m a x}$	——	Lower and upper bounds of active load connected to branch l
$P_{0}$ , $Q_{0}$	——	No-load active and reactive power losses of transformer
$q_{l}^{m i n}$ , $q_{l}^{m a x}$	——	Lower and upper bounds of reactive load connected to branch l
$Q_{l}^{S V C, m i n}$ , $Q_{l}^{S V C, m a x}$	——	The minimum and maximum reactive power for static var compensator (SVC) connected to branch l
$r_{l}$ , $x_{l}$ , $z_{l}$ , $b_{l}$	——	Resistance, reactance, impedance, and half-shunt susceptance on branch l
$r_{1}$ , $x_{1}$	——	Rsistance of transformer and transformer leakage reactance
$S_{l}^{D G N}$	——	Nominal capacity of DG on branch l
$T_{m a x}$	——	Total number of scheduling periods
$v_{m i n}$ , $v_{m a x}$	——	Lower and upper bounds of voltage magnitude square
$v_{0}$	——	Voltage square of substation
D.	——	Variables
$ψ_{t}$	——	The first-stage variable during time period t
$φ_{t}$	——	The second-stage variable during time period t
$π_{(\cdot)}^{(\cdot)}$ , $π_{(\cdot)}^{}$	——	The dual variables
$C_{l, t}^{}$	——	Capacitor of SCR during time period t
$g_{l, t}^{}$	——	Predicted active power vector of DG connected to branch l during time period t
g_t	——	Active power vector of DG
$O_{t}^{}$	——	OLTC tap travel distance during time period t
$P_{1}^{t}$	——	Active power injected from the root node to bulk power system
$p_{l}$ , $q_{l}$ , s_l	——	Active, reactive, and complex power loads of connected to branch l
$P_{l}^{t}$ , $Q_{l}^{t}$ , $S_{l}^{t}$	——	Active, reactive, and complex power injected into the top of branch l
$P_{l}^{b}$ , $Q_{l}^{b}$ , $S_{l}^{b}$	——	Active, reactive, and complex power from the bottom of branch l
$P_{m}^{b}$ , $Q_{m}^{b}$ , $S_{m}^{b}$	——	Active, reactive, and complex power from the bottom of branch m
$P_{m}^{t}$ , $Q_{m}^{t}$ , $S_{m}^{t}$	——	Active, reactive, and complex power injected into the top of branch m
$P_{l, t}^{b}$ , $Q_{l, t}^{b}$ , $S_{l, t}^{b}$	——	Active, reactive, and complex power from the bottom of branch l during time period t
$P_{l}^{c h}$ , $P_{l}^{d i s}$	——	Charging and discharging power for ESS on branch l
$P_{l, t}^{E S S}$	——	Charging minus discharging power for ESS on branch l during time period t
$P_{l, t}^{D G}$ , $Q_{l, t}^{D G}$ ,	——	Active, reactive, and rated power of DG with full capacity of converter connected to branch l during time period t
Sl,tDGN
$Q_{l}^{c}$	——	Reactive power injected into the center of the equivalent circuit of a cable on branch l
$Q_{l}^{S V C}$ , $Q_{l}^{S C R}$	——	Injected reactive power from SVC, SCR on branch l during time period t
s	——	Complex power vector of load
S	——	Complex power vector of branch flow
UB, LB	——	Upper and lower bounds of the robust optimization problems
v, f	——	Vectors of squared voltage and current magnitudes
$v_{u p (l)}^{b}$	——	Voltage magnitude squares of downstream branch l
$v_{l}^{}$ , $f_{l}^{}$	——	Voltage and current magnitude squares on branch l
$X_{t, t + 1}^{}$	——	Number change for SCR groups operating between time periods $t$ and $t + 1$
$y_{t}^{}$	——	Vector of second stage variables excluding dummy variables
E.	——	Operators
$R$	——	Real part
$ℑ$	——	Imaginary part
$*$	——	Conjugate operation
${(\cdot)}^{T}$	——	Transpose of a matrix
$\circ$	——	Hadamard product
F.	——	Symbols
$-$	——	Upper bound
^	——	Lower bound

References

J. Zhang, M. Cui, and Y. He, “Robustness and adaptability analysis for equivalent model of doubly fed induction generator wind farm using measured data,” Applied Energy, vol. 261, pp. 1-14, Mar. 2020. [Baidu Scholar]

J. Zhang, Y. He, M. Cui et al., “Primal dual interior point dynamic programming for coordinated charging of EVs,” Journal of Modern Power Systems and Clean Energy, vol. 5, no. 6, pp. 1004-1015, Nov. 2017. [Baidu Scholar]

H. Gao, J. Liu, and L. Wang, “Robust coordinated optimization of active and reactive power in active distribution systems,” IEEE Transactions on Smart Grid, vol. 9, no. 5, pp. 4436-4447, Sept. 2018. [Baidu Scholar]

S. Xie, Z. Hu, and J. Wang, “Two-stage robust optimization for expansion planning of active distribution systems coupled with urban transportation networks,” Applied Energy, vol. 261, pp. 1-17, Mar. 2020. [Baidu Scholar]

B. A. Robbins, H. Zhu, and A. D. Domínguez-García, “Optimal tap setting of voltage regulation transformers in unbalanced distribution systems,” IEEE Transactions on Power Systems, vol. 31, no. 1, pp. 256-267, Jan. 2016. [Baidu Scholar]

E. Dall’Anese, H. Zhu, and G. B. Giannakis, “Distributed optimal power flow for smart microgrids,” IEEE Transactions on Smart Grid, vol. 4, no. 3, pp. 1464-1475, Sept. 2013. [Baidu Scholar]

J. Zhang, M. Cui, H. Fang et al., “Two novel load balancing platforms using common DC buses,” IEEE Transactions on Sustainable Energy, vol. 9, no. 3, pp. 1099-1107, Jul. 2018. [Baidu Scholar]

J. Lavaei and S. H. Low, “Zero duality gap in optimal power flow problem,” IEEE Transactions on Power Systems, vol. 27, no. 1, pp. 92-107, Feb. 2012. [Baidu Scholar]

J. Lavaei, D. Tse, and B. Zhang, “Geometry of power flows and optimization in distribution networks,” IEEE Transactions on Power Systems, vol. 29, no. 2, pp. 572-583, Mar. 2014. [Baidu Scholar]

L. Gan, N. Li, U. Topcu et al., “Exact convex relaxation of optimal power flow in radial networks,” IEEE Transactions on Automatic Control, vol. 60, no. 1, pp. 72-87, Jan. 2015. [Baidu Scholar]

M. Farivar and S. H. Low, “Branch flow model: relaxations and convexification-part I,” IEEE Transactions on Power Systems, vol. 28, no. 3, pp. 2554-2564, Aug. 2013. [Baidu Scholar]

M. Farivar and S. H. Low, “Branch flow model: Relaxations and convexification-Part II,” IEEE Transactions on Power Systems, vol. 28, no. 3, pp. 2565-2572, Aug. 2013. [Baidu Scholar]

M. Nick, R. Cherkaoui, J. L. Boudec et al., “An exact convex formulation of the optimal power flow in radial distribution networks including transverse components,” IEEE Transactions on Automatic Control, vol. 63, no. 3, pp. 682-697, Mar. 2018. [Baidu Scholar]

Y. Xiang, J. Liu, and Y. Liu, “Robust energy management of microgrid with uncertain renewable generation and load,” IEEE Transactions on Smart Grid, vol. 7, no. 2, pp. 1034-1043, Mar. 2016. [Baidu Scholar]

C. Lee, C. Liu, S. Mehrotra et al., “Robust distribution network reconﬁguration,” IEEE Transactions on Smart Grid, vol. 6, no. 2, pp. 836-842, Mar. 2015. [Baidu Scholar]

T. Ding, C. Li, Y. Yang et al., “A two-stage robust optimization for centralized-optimal dispatch of photovoltaic inverters in active distribution networks,” IEEE Transactions on Sustainable Energy, vol. 8, no. 2, pp. 744-754, Apr. 2017. [Baidu Scholar]

T. Ding, S. Liu, W. Yuan et al., “A two-stage robust reactive power optimization considering uncertain wind power integration in active distribution networks,” IEEE Transactions on Sustainable Energy, vol. 7, no. 1, pp. 301-311, Jan. 2016. [Baidu Scholar]

B. Zeng and L. Zhao, “Solving two-stage robust optimization problems using a column-and-constraints generation method,” Operations Research Letters, vol. 41, no. 5, pp. 457-461, Sept. 2013. [Baidu Scholar]

J. Zhang, M. Cui, Y. He et al., “Fast robust optimization for partial distributed generators (DG) providing ancillary services,” Energies, vol. 14, p. 4911, Aug. 2021. [Baidu Scholar]

M. Mazidi, H. Monsef, and P. Siano, “Robust day-ahead scheduling of smart distribution networks considering demand response programs,” Applied Energy, vol. 178, pp. 929-942, Sept. 2016. [Baidu Scholar]

H. Qiu, W. Gu, J. Pan et al., “Multi-interval-uncertainty constrained robust dispatch for AC/DC hybrid microgrids with dynamic energy storage degradation,” Applied Energy, vol. 228, pp. 205-214, Oct. 2018. [Baidu Scholar]

B. Zhao, J. Ren, J. Chen et al., “Tri-level robust planning- operation co-optimization of distributed energy storage in distribution networks with high PV penetration,” Applied Energy, vol. 279, pp. 1-11, Dec. 2020. [Baidu Scholar]

A. Ehsana and Q. Yang, “State-of-the-art techniques for modelling of uncertainties in active distribution network planning: a review,” Applied Energy, vol. 239, pp. 1509-1523, Apr. 2019. [Baidu Scholar]

J. Zhang, M. Cui, H. Fang et al., “Two novel load-balancing platforms using common DC buses,” IEEE Transactions on Sustainable Energy, vol. 9, no. 3, pp. 1099-1107, Jul. 2018. [Baidu Scholar]

M. Cui and J. Wang, “Deeply hidden moving-target-defense for cybersecure unbalanced distribution systems considering voltage stability,” IEEE Transactions on Power Systems, vol. 36, no. 3, pp. 1961-1972, May 2021. [Baidu Scholar]

C. Chen, M. Cui, X. Fang et al., “Load altering attack-tolerant defense strategy for secondary frequency control system,” Applied Energy, vol. 280, p. 116015, Dec. 2020. [Baidu Scholar]

W. Wei, J. Wang, N. Li et al., “Optimal power flow of radial networks and its variations: a sequential convex optimization approach,” IEEE Transactions on Smart Grid, vol. 8, pp. 2974-2987, Nov. 2017. [Baidu Scholar]

Address:No.19 Chengxin Avenue, Jiangning District, Nanjing 211106, China

E-mail: mpce@alljournals.cn

Tel:86-25-81093060

Fax:86-25-81093040

Home

Introduction

Editorial Board

For Author

Call For Papers

APC

Sponsor & Publisher

Multi-period Two-stage Robust Optimization of Radial Distribution System with Cables Considering Time-of-use Price PDF

Abstract

Keywords

I. Introduction

II. Models of Radial Distribution System with Cables

A. SOCP Model of Radial Distribution System with Cables

B. Model of ESS

C. Model of SVC

D. Model of SCR

E. Model of OLTC

F. Model of DG

G. Objective Function

III. Robust Optimization Model

A. Deterministic Model

B. Robust Optimization Model

IV. Solution Method for Robust Optimization Model

A. Subproblem (SP)

B. Solving Steps

C. Configurations of Three Test Systems

D. Simulation Results

E. Comparison of Computational Performance of Proposed and CCG Methods

V. Conclusion

Nomenclature

References