The power flow equation of radial distribution networks is formulated as follows (Farivar and Low, 2013a): P j t = ∑ l ∈ m j P j l t − ∑ l ∈ n j P i j t − R i j I i j t 2 Q j t = ∑ l ∈ m j Q j l t − ∑ l ∈ n j Q i j t − X i j I i j t 2 , (1) where P j t and Q j t denote the active and reactive powers, respectively, injected into node j at time t ; P j l t and Q j l t are the active and reactive powers, respectively, flowing from node j to node l at time t ; R i j and X i j are the resistance and reactance, respectively, of line i j ; I i j t is the current flowing through line i j at time t ; n j is the set of sending nodes with node j as the end node; and m j is the set of the end nodes with node j as the sending node.

In addition, the relationship between the current I i j t flowing through line i j and the node voltage is given as follows: U i t 2 − U j t 2 = 2 R i j P i j t + X i j Q i j t − I i j t 2 R i j 2 − I i j t 2 X i j 2 , (2) I i j t 2 U i t 2 = P i j t 2 + Q i j t 2 , (3) where U i t is the voltage of node i at time t .

In the process of transmitting power in the distribution network, reactive power compensation can effectively improve the power factor and reduce line losses. Common reactive power compensation devices include continuous-type Static Var Generators (SVG) and discrete-type capacitors (CP), where CP satisfies Q C P t = Q C P , U ∑ m = 1 M θ C P , m t , (4) θ C P , m t > = θ C P , m + 1 t (5) θ C P _ a d d t + θ C P _ d e l t < = 1 (6) ∑ t = 1 T θ C P _ a d d t + θ C P _ d e l t < = θ max (7) ∑ m = 1 M θ C P , m t + 1 − ∑ m = 1 M θ C P , m t ≤ γ θ C P _ a d d t − θ C P _ d e l t , (8) ∑ m = 1 M θ C P , m t + 1 − ∑ m = 1 M θ C P , m t ≥ θ C P _ a d d t − γ θ C P _ d e l t , (9) where Q C P t is the compensated reactive power of the CP at time t ; Q C P , U is the reactive power output corresponding to each step of the CP; θ C P , m t , θ C P _ a d d t , and θ C P _ d e l t denote the 0–1 variables at time t ; γ is the maximum number of steps that can be adjusted; and θ max is the maximum number of adjustments per day.

It should be pointed out that in Eq. 4, a 0–1 variable is set for each grade, and the sum of each grade is the current grade. Equation 5 limits the gears to incremental steps, and it is not possible to have a high gear of 1 and a low gear of 0. Equation 6 ensures that the increase and decrease of gears cannot occur simultaneously, and Eq. 7 shows that the number of gear adjustments per day is limited. Equations 8, 9 together restrict the upward or downward adjustment of gears to the maximum adjustment range of gears.

For the SVG, the model is simpler than the CP and satisfies Q S V G min ≤ Q S V G t ≤ Q S V G max , (10) where Q S V G min and Q S V G max are the lower and upper limits of the regulation range, respectively, and Q S V G t is the reactive power of the SVG at time t .

The input and output energy of a power system should satisfy the law of energy conservation, so the following relationship should be ensured: Q j t + Q L , j t = Q S V G , j t + Q C P , j t + Q G , j t , P j t + P L , j t = P G , j t (11) where Q L , j t and P L , j t are the reactive and active loads of node j at time t , respectively. P G , j t and Q G , j t are the active and reactive powers of the generator of node j at time t , respectively.

The gas system is mainly composed of a gas source, gas transmission pipelines, a compressor, and loads, and the relationship between its flow and pressure is formulated as follows: f q t = M q ν m n ν m n p m t 2 − p n t 2 (12) M q = ς 18.73 T 0 D q 8 3 p 0 G L q T q a Z q a (13) ν m n = 1 − 1   p m t − p n t ≥ 0 p m t − p n t < 0 (14) where f q t is the flow of the pipe q at time t , p m t is the pressure of node m at time t , ς represents the efficiency factor of the pipe q , T 0 is the standard temperature, D q denotes the inner diameter of the pipe q , p 0 is the standard pressure, G is the gas relative density, L q is the length of the pipe q , T q a is the average temperature of the pipe q , and Z q a is the average compression factor of the pipe q .

There is a pressure loss during gas transmission, which needs to be replenished by a compressor. The gas flow to be consumed is formulated as σ u t = α u + β u P c o m t + γ u P c o m t 2 (15) P c o m t = G u f c o m , u ( p m t p n t Z u − 1 (16) where σ u t is the gas flow consumed by the compressor u at time t ; α u , β u , and γ u are the energy consumption coefficients of the compressor u ; P c o m t denotes the gas energy consumed by the compressor u at time t ; G u is the coefficient related to temperature and efficiency; f c o m , u represents the flow through the compressor u at time t ; and Z u is the other compression factor.

Similar to the law of conservation of energy in electric power systems, the input and output of gas at each node in a gas system should be equal, which satisfies A + U f + ϖ − T τ = 0 , (17) where f denotes the branch flow vector, ϖ is the node gas injection flow vector, τ is the compressor consumption flow vector, A is the branch node association matrix, U represents the unit node association matrix, and T is the compressor node association matrix.

The power flow model and natural gas model is the non-convex and non-linear. In order to facilitate the subsequent solution, second-order cone scaling (Farivar and Low, 2013a; Farivar and Low, 2013b) is used to convert the distribution network and natural gas system models into the convex model.

The operation of the upper-layer ESN is optimized by minimizing net costs as follows: C u p = ∑ t = 1 T a + b P G t + c ( P G t 2 + h G G t − ∑ t = 1 T p r i c e e t P g r i d t + p r i c e g t G g a s t P G t = ∑ k = 1 N P G , k t , G G t = ∑ k = 1 M G G , k t o b j = min C u p , (18) where a , b , and c denote the generators’ cost coefficients; h is the cost coefficient of gas supply; G G t is the gas supply quantity of the gas source at time t ; p r i c e e t and p r i c e g t represent the purchase prices of electricity and gas of the lower-layer DES at time t , respectively; P g r i d t and P g a s t are the purchased power and gas power of the lower-layer DES from the upper-layer ESN at time t , respectively; and T is the dispatch period. P G , k t and G G , k t are the k-th generator and air source supply, respectively. N and M are the total number of generators and gas sources, respectively.

In the upper-layer ESN, upper- and lower-limit constraints should be added to the individual equipment output, voltage value, gas flow value, node pressure, and gas supply from the gas source, which are expressed uniformly by Θ min ≤ Θ t ≤ Θ max , (19) where Θ min and Θ max are the lower and upper limits of the variable Θ t , respectively. Θ t represents the variables such as equipment output, voltage value, gas flow value, node pressure, and gas supply from the gas source.

Energy storage equipment is modeled as E i t = E i t − 1 1 − δ i + υ i C i t η i , c h a r − 1 − υ i D i t η i , d i s (24) where E i t is the reserves of energy storage at time t ; i is the type of energy storage device, which includes electricity and heat storage; δ i is the self-loss rate of the energy storage device; C i t and D i t are the charging and discharging powers of the energy storage device, respectively; η i , c h a r and η i , d i s represent the charging and discharging efficiencies of the energy storage device, respectively; and υ i is a 0–1 variable to represent the storage mode: 0 if discharging, and 1 if charging.

There are also various constraints on the energy storage device, for example, charging and discharging cannot be done simultaneously and the total energy storage remains unchanged in a scheduling cycle, in addition to upper and lower limits on variable parameters, which are formulated as C i t ≤ υ i t C i max D i t ≤ 1 − υ i D i max E i 1 = E i 24 E i t ∈ E i min , E i max (25)

PV and various load data are mainly derived by forecasting; however, there are always various random factors that affect the accuracy of forecasting in reality. Random factors cannot be avoided, and if they are not considered when developing the operation scheme, they may lead to inaccurate operation schemes and even cause operational accidents. To equip the system with a certain anti-interference ability and to ensure the system’s safe and reliable operation in different environments, it is necessary to deal with those uncertain parameters. The uncertainty U model of PV output and load is given by U = u = u P V t , u i , L t ∈ R T × 4 , t = 1,2 , ⋯ , T u P V t = u ^ P V t − B P V t Δ u P V max ⁡ , t + D P V t Δ u P V max ⁡ , t ∑ t = 1 T B P V t + D P V t ≤ Γ P V u i , L t = u ^ i , L t − B i , L t Δ u i , L max ⁡ , t + D i , L t Δ u i , L max ⁡ , t ∑ t = 1 T B i , L t + D i , L t ≤ Γ i , L (26) where u P V t and u i , L t are the uncertainty variables of PV output and load power at time t after introducing uncertainty, respectively. u ^ P V t and u ^ i , L t are the expected values of PV output and load power at time t , respectively. Δ u P V max ⁡ , t and Δ u i , L max ⁡ , t denote the maximum allowable deviations of PV output and load power at time t , respectively; i represents the energy types including electricity, heat, and cooling; and B P V t , D P V t , B i , L t , and D i , L t are 0–1 variables, while Γ P V and Γ i , L are the uncertainty conservativeness adjustment parameters.

The lower-layer DES also needs to satisfy the energy conservation law, and the balance between the input and output quantities of each type of energy should be maintained as P g r i d t + u P V t + P G T t + D e t = u e , L t + C e t + P E C t H G T t + H G B t + D h t = u h , L t + C h t + P A C t C E C t + C A C t = u c , L t G g a s t = G G T t + G G B t , (27) where D e t , C e t , D h t , and C h t are the variables of charging and discharging powers corresponding to the electric and thermal energies, respectively, while u e , L t , u h , L t , and u c , L t are the variables of load uncertainty corresponding to the electricity, heat, and cooling, respectively.

Like the upper layer, there are also some limit constraints on the output of each device in the lower layer, which are omitted here due to the length of the paper.

After presenting the mathematical models of both layers, the following objective function can be formulated: C d n = ∑ t = 1 T p r i c e e t P g r i d t + p r i c e g t G g a s t o b j = min x max u ∈ U min y ∈ Ω x , u C d n x = υ i t , C i t , D i t , E i t y = P g r i d t , G g a s t , G G T t , P G T t , H G T t , G G B t , H G B t , H A C t , C A C t , P E C t , C E C t . (28)

To facilitate subsequent simplification and solution, the lower-layer DES model is written in a matrix compact form: min x max u ∈ U min y ∈ Ω x , u C T y s . t . D y ≥ d L y + T x = 0 I u y = u P x = p F x ≥ f . (29)

The objective function can be split into two stages: the first stage is to optimize the charging and discharging behavior of the power and heat storage equipment, called the main problem (MP), and the second stage is to search for the optimal equipment output and energy purchased under the worst uncertain parameters, called the sub-problem (SP). Uncertainty is expressed by the box uncertainty set. The sub-problem is to solve the worst scenario in the interval according to the conservatism parameter and complete the two-stage robust optimization together with the main problem. The main problem and sub-problem of the objective function are described as follows: M P : min x σ s . t . D y ≥ d L y + T x = 0 I u y = u P x = p F x ≥ f σ ≥ C y , (30) S P : Q x = max u ∈ U min y ∈ Ω x , u C T y s . t . D y ≥ d L y + T x = 0 I u y = u . (31)

The sub-problem is a max–min type problem, which is not easy to solve. By the strong dual principle, the inner min problem can be converted into a max problem and merged with the outer one to obtain the following form of the sub-problem (Ding et al., 2016): S P : Q x = max u ∈ U , ω d T ω 1 − T x T ω 2 + u T ω 3 s . t . D T ω 1 + L T ω 2 + I u T ω 3 ≤ C ω 1 ≥ 0 , ω 3 ≥ 0 , (32) where ω 1 , ω 2 , and ω 3 are dual variables.

In Eq. 32, u T ω 3 is a bilinear variable, which is difficult to solve. Therefore, it is linearized by using the big M method: S P : Q x = max u ∈ U , ω d T ω 1 − T x T ω 2 + u ^ T ω 3 + Δ u T B ′ s . t . D T ω 1 + L T ω 2 + I u T ω 3 ≤ C 0 ≤ B ′ ≤ M B I G B ω 3 − M B I G 1 − B ≤ B ′ ≤ ω 3 ∑ t = 1 T B P V t + D P V t ≤ Γ P V ∑ t = 1 T B i , L t + D i , L t ≤ Γ i , L ω 1 ≥ 0 , ω 3 ≥ 0 , (33) where u ^ = u ^ P V t , u ^ i , L t , Δ u = Δ u P V max ⁡ , t , Δ u i , L max ⁡ , t , B ′ is the auxiliary variable, and M B I G is a large positive real number, which is taken as 10⁹ in this study.

In this study, the RIES is divided into two layers. The upper layer is a traditional ESN model, while the lower-layer DES is a complex two-stage robust optimization model. The CCG algorithm is used to obtain the optimal solution to the lower layer. The process is described as follows:

1) The initial uncertainty parameter u 1 * , the iteration convergence threshold ε , the upper and lower bounds for the iterative solution as L B = − ∞ , U B = + ∞ , and the number of iterations as k = 1 are set.

It is worth pointing out that the sum of the objective functions for the upper and lower layers is the complete objective function of RIES optimization, as the two layers are not separate. However, if the upper and lower layers are structured as one model, its complexity will greatly increase, and it will be more difficult to solve it. For this reason, the ADMM algorithm is used to maintain the correlation between the upper and lower layers while solving the models of the upper and lower layers independently. The flow of the ADMM algorithm is described as follows:

Step 1:

The augmented Lagrange function is established for the upper and lower layers, where it is formulated for the upper layer as min L u p = C u p + λ 1 , u p t P g r i d , u p t − P g r i d t + λ 2 , u p t G g a s , u p t − G g a s t + η 2 P g r i d , u p t − P g r i d t 2 + η 2 G g a s , u p t − G g a s t 2 . (35)

For the lower layer, it is min L d n = C d n + λ 1 , d n t P g r i d , d n t − P g r i d t + λ 2 , d n t G g a s , d n t − G g a s t + η 2 P g r i d , d n t − P g r i d t 2 + η 2 G g a s , d n t − G g a s t 2 (36) where L u p and L d n are the augmented Lagrange function of the upper and lower layers, respectively; P g r i d t and G g a s t are the global variables; P g r i d , u p t , G g a s , u p t , P g r i d , d n t , and G g a s , d n t represent the local variables of purchased electricity and gas in the upper and lower layers; η is the coefficient of penalty term; and λ 1 , u p t , λ 2 , u p t , λ 1 , d n t , and λ 2 , d n t are Lagrange multipliers.

The lower layer is a two-stage robust optimization model: MP and SP. The MP can solve the electricity and gas purchase, so only the objective function of the MP is changed to the augmented Lagrange function: min L d n = σ + λ 1 , d n t P g r i d , d n t − P g r i d t + λ 2 , d n t G g a s , d n t − G g a s t + η 2 P g r i d , d n t − P g r i d t 2 + η 2 G g a s , d n t − G g a s t 2 (37)

Step 2:

Initialization: We assume P g r i d , u p , 0 t = 0 , G g a s , u p , 0 t = 0 , λ 1 , u p t = λ 0 , λ 2 , u p t = λ 0 , λ 1 , d n t = λ 0 , λ 2 , d n t = λ 0 , and λ 0 as the initial values and the number of iterations as r = 0 .

Step 3:

Let P g r i d t = P g r i d , u p , r t and G g a s t = G g a s , u p , r t are substituted into the lower-layer problem to find P g r i d , d n , r + 1 t and G g a s , d n , r + 1 t .

Step 4:

Let P g r i d t = P g r i d , d n , r + 1 t and G g a s t = G g a s , d n , r + 1 t are substituted into the upper-layer problem to find P g r i d , u p , r + 1 t and G g a s , u p , r + 1 t .

Step 5:

The Lagrange multipliers are updated according to the following equations: λ 1 , d n , r + 1 t = λ 1 , d n , r t + η P g r i d , d n , r + 1 t − P g r i d , u p , r + 1 t λ 2 , d n , r + 1 t = λ 2 , d n , r t + η G g a s , d n , r + 1 t − G g a s , u p , r + 1 t λ 1 , u p , r + 1 t = λ 1 , u p , r t + η P g r i d , u p , r + 1 t − P g r i d , d n , r + 1 t λ 2 , u p , r + 1 t = λ 2 , u p , r t + η G g a s , u p , r + 1 t − G g a s , d n , r + 1 t (38)

Step 6:

The maximum deviation value is calculated. When the maximum deviation is less than the convergence threshold, the iteration is terminated and the results are output; otherwise, we return to step 3. The maximum deviation is calculated by φ = max P g r i d , d n t − P g r i d , u p t , G g a s , d n t − G g a s , u p t , (39) where φ denotes the maximum deviation.

It should be remarked that as the iteration proceeds, the Lagrange multipliers gradually become larger. Consequently, the proportion of the difference between the local and global variables in the objective function increases. So, the minimum solution will be made as small as possible until it is close to 0. Because the global variable is constantly taken from the local variable in the upper and lower layers, the difference between the corresponding local variable in the upper and lower layers will become smaller as the iteration proceeds, until it is close to 0, which is the principle of the ADMM algorithm. Because the expected value of the predicted value is not much different from the final solution value, the expectation of the predicted value is used as the initial value of the parameter of the CCG algorithm. The flowchart of the proposed solving method is given in Figure 1 for the overall optimization process.

The arithmetic model selected as a test example is shown in Figure 2. The upper-layer ESN includes a modified IEEE 33-node power system and a Belgian 20-node natural gas system, and the lower-layer DES model is a typical campus energy supply system. It should be pointed out that the method proposed in this study can be applied to different test examples, and only the relevant parameters need to be changed when applying to other test examples.

The forecast curves for each type of load in the lower-layer DES and PV are shown in Figure 3. Without loss of generality, the fluctuation range of the load is 0.1 of the predicted value, and the fluctuation range of PV is 0.15 of the predicted value. Figure 4 shows the time-of-use electricity price, and the gas price is set as 0.35 CNY/kW. The gas turbine’s maximum capacity is 1,400 kW, the maximum capacity of the gas boiler is 800 kW, the capacity of the electric chiller is 500 kW, the absorption chiller’s maximum capacity is 300 kW, the maximum storage capacity of electricity storage equipment and heat storage equipment is 1,500 kW h, and the minimum storage capacity is 200 kW h. The maximum value of charging and discharging power is 600 kW. Other parameters are selected according to references (Ding et al., 2017; Miao et al., 2020).