Coordinated Dispatch of Multi-Energy Microgrids and Distribution Network with a Flexible Structure

Abstract

With the increasing scale of multi-energy microgrids (MGs) and complicated operation modes, the coordinated operation of microgrids and the distribution network (DN) has posed great challenges. In this paper, a bi-level optimal coordinated dispatch framework of the DN and multi-energy MGs based on CCHP (combined cooling, heating, and power) is proposed. The first level studies the optimal operation of the DN with power interaction on tie lines between MGs considering the coupling relationship and constraints of the equipment and network. The network reconfiguration with limited control actions is considered to increase the flexibility of the topology and further improve the working state. For the second level, MGs receive orders from the DN and determine the optimal strategies of multi-energy devices to achieve optimized operation under the condition of satisfying the different types of load and requirement for the DN. To solve the optimal dispatch problem of both the DN and the multi-energy MGs considering the DN reconfiguration, a method combining particle swarm optimization algorithm (PSO) with mixed-integer linear programming (MILP) is proposed. Cases studied in an IEEE33-node DN with renewable power sources and grid-connected MGs validate that the proposed method is very effective in reducing the power loss and voltage offset of the DN while ensuring the benefits of the MGs.

1. Introduction

With the gradual depletion of traditional primary energy, various forms of renewable energy have been widely used [1]. As the penetration rate of renewable energy gradually increases, the distribution networks (DNs) transform into active networks, the characteristics of which are that control and dispatch are distributed, power flows are bidirectional [2]. Microgrids (MGs) are presented to connect to DNs in response to the growing power supply security issues, as well as considering the economics of energy, cleanliness and various ancillary services [3,4]. With the maturity of distributed power sources and energy storage technology, multi-MGs’ access to distribution networks gradually become the mainstream. As a controllable integrated unit of distributed power supply, energy storage device and user load, the MG can be flexibly controlled and managed and can run in grid-connected or island mode. It is of considerable significance to increase the economic and environmental benefits of the power system.

The traditional research of MGs mainly concentrates on the utilization of electric power. With the introduction of the concept of energy internet and the development of communication technology, the barriers between various forms of energy are being broken [5,6]. The multi-energy systems integrate multiple types of energy, such as natural gas, coal, petroleum, and renewable power sources, which are highly coupled with production, transmission, and management. Combined cooling, heating and power (CCHP) technology is widely used in small and medium-sized MGs, which leads to new challenges for optimized dispatch of MGs [7,8,9,10].

Much effort has been focused on the optimal dispatch of the MGs [11,12,13,14,15,16,17,18,19,20]. Wang Chengshan proposes a general modeling method using an energy bus to describe the injection of different forms of energy in [12]. To minimize fuel costs, Seon-Ju uses the direct search method to tackle the economic problem considering load and DG constraints [11]. Chunyang Liu develops a mixed-integer linear programming (MILP) dispatch model of MGs considering the lifetime of the battery under the weighted Wh throughput method [13]. For multi-energy microgrids, operating modes proceed following the electric load (FEL), following the thermal load (FTL) and optimal economic dispatch, etc. [14]. Li Lin proposes an economic dispatch method of CCHP-based multi-microgrids considering the power interaction [15]. Liu N proposes an energy sharing strategy, which let CHP units to break the limit of working modes [16]. Bui, V develops an internal trading method with barter trading, which can reduce the power trading with the external system [17]. For stochastic operation [18,19,20], Li Yong proposes a low carbon micro integrated electric power, natural gas, and heat delivery system, considering uncertainties of renewable power and demand under the scenario-based method in [18]. A hierarchical energy management strategy is presented in [19] to establish an efficient framework for MGs considering the uncertainty of power interaction among MGs. In [20], a joint-dispatch method of reserve and energy is presented for the study of CCHP MGs considering the thermal dynamics of buildings and forecast errors of renewable power. When considering demand response [21,22], Rakipour establishes an energy hub with multi-energy flow considering uncertainties of renewable power and different degrees of demand response of cooling and electric load to maximize income in [21]. In [22], a mechanism of bi-transaction mode between CCHP-based-MG and a load aggregator including pricing strategy and a two-stage optimization model is introduced to achieve integrated demand response.

As renewable power sources, MGs, distributed energy storage systems, and controllable loads are connected to the DN, power quality monitoring [23] and optimal dispatch [24] of the DN has become increasingly important. There are also many types of research on the coordinate dispatch of the DN and MGs [25,26,27,28,29]. The MGs can provide the DN auxiliary service [25] to improve the power quality and voltage stability of the DN. Moreover, the MGs have better ability in accommodation of renewable energy and some low-carbon distributed power sources, increasing environmental benefits, and having a positive impact on the DN operation [26]. As two stakeholders, MGs and DN should have their own interests. To maximize their benefits with coordination, in [27] DN is equivalent to a coupled multi-microgrid system, using a decentralized observable Markov decision process with a dual multiplier based on Lagrangian relaxation.

Some research not only considers the electric power in the MGs but also takes cooling and heating load into consideration. With the increasingly close coupling of various types of energy, the distribution and regulation of energy are more diversified, providing more choices for ensuring network operation security. CCHP-based MGs have greater flexibility in energy management due to its multi-energy coupling characteristics [28]. In [29], CCHP-based MGs and DN are regarded as two stakeholders, using the analytical target cascading algorithm (ATC) to achieve the optimization goals under the distributed controlling method. Nevertheless, the topology of the DN is not considered. Flexible topologies also have a crucial impact on reducing losses and improving the reliability of DN, which can be achieved through DN reconfiguration [30,31,32].

In summary, some of the main deficiencies of the above research can be generalized as follows:

(1)

Many literatures only consider the optimal operation of the MG itself, but do not consider the working status of the DN connected to the MG, which will have a bad impact on the entire distribution system;

(2)

Some other literature only considers the dispatch of electric loads when considering the optimal joint dispatching of the MG and the DN regardless of other types of loads;

(3)

A few works in the literature consider the joint dispatch of multi-energy MGs and the DN but do not consider the topology and power flow of the distribution network itself.

Thus, to better study the operation of the MGs-integrated DN, a bi-level optimal dispatch method of multi-energy MGs and DN is proposed in this paper. The framework of the DN with the integration of the distributed generators (DG) and the MGs is shown in Figure 1. A particle swarm optimization (PSO) algorithm with mixed-integer linear programming (MILP) is used to solve the optimization of the bi-level problem. The advantages of the proposed model are as follows:

(1)

Cooling and heating loads are taken into consideration, which will increase the complementary use of energy and more dispatch flexibility.

(2)

The load flow optimization problem of the DN under the integration of the multi-MGs is considered, the aim of which is to improve the operation conditions (minimize the power loss and voltage offset of the DN) and to make the model more accurate.

(3)

The DN reconfiguration is included in the proposed method to make further efforts to optimize the operation mode and the security of the DN under the integration of MGs and to improve control flexibility with limited control actions, which will also make the MGs sacrifice the minimum economic benefit under the framework of the coordinated dispatch method.

The rest of the paper is organized as follows. Section 2 presents the mathematical model of MGs and DN. The form of the problem and the objective function are introduced in Section 3. Section 4 presents the solving method. Numeral results are given and analyzed in Section 5. Finally, the conclusion is presented in Section 6.

2. Mathematical Model Formulation

2.1. Multi-Energy Microgrids

Figure 2 illustrates the underlying architecture of multi-energy MGs, whose connection relationships are illustrated in Figure 1. The devices in MGs include photovoltaic power stations, wind turbines, energy storage systems, gas turbines, gas boilers, waste heat boiler (WH), heat exchangers, absorption chillers, electric chillers, etc. [13,19]. The definitions of the symbols in Figure 2 are shown in

Table 1. Definitions of the symbols in Figure 2.

Symbol	Definition	Symbol	Definition
PV	Photovoltaic array	ESS	Energy storage system
WT	Wind turbine	EC	Electric chiller
GS	Gas station	AC	Absorption chiller
GB	Gas boiler	GT	Gas turbine
WH	Waste heat boiler	HX	Heat exchanger

.

2.1.1. Gas Turbine (GT)

The gas turbine (GT) burns natural gas to generate electricity and heat. The partial load characteristics of the efficiency of gas turbine power generation can be expressed as the following formula [10]:

$${\eta }_{GT}^E=(a\cdot PL{R}^3+b\cdot PL{R}^2+c\cdot PLR+d){\eta }_{GT}^{En}$$

(1)

where ${\eta }_{GT}^E$ is the electrical efficiency of GT. ${\eta }_{GT}^{En}$ is the rated electrical efficiency of the GT. PLR is the partial load rate (The ratio of actual load value to rated load value) of GT.$a,b,c,d$ are constants.

The mathematical model of gas consumption and residual heat produced by the exhaust gas of GT can be expressed as:

$$F_{GT}(t)=\frac{P_{GT}(t)}{R_G{\eta }_{GT}^E}$$

(2)

$$H_{GT}(t)=F_{GT}(t)R_G(1-{\eta }_{GT,l}^E-{\eta }_{GT}^E)$$

(3)

where F_GT(t) is the gas consumption at time slot t. P_GT(t) is the electrical power of the GT at time slot t. R_G is the calorific value of natural gas. η^E_GT,l is the loss factor of the GT. H_GT(t) is the waste heat generated by GT. The output electric power of the GT meets the operational constraints below:

$$u_{GT}(t)P_{GT}^{\min }\le P_{GT}(t)\le u_{GT}(t)P_{GT}^{\max }$$

(4)

where u_GT(t) is the binary variable of the working condition of GT. u_GT(t) = 1 means GT is operating at time slot t. u_GT(t) = 0 means GT is shutdown at time slot t. $P_{GT}^{\max }$, $P_{GT}^{\min }$ are maximal and minimal electrical power of GT respectively

2.1.2. Gas Boiler (GB)

A gas boiler (GB) is a gas-fired boiler directly supplying heating load by natural gas combustion. The mathematical model can be expressed as:

$$\begin{array}{l}{H}_{GB}(t)={\eta }_{GB}{F}_{GB}(t)R_G\\ H_{GB}^{\min }\le H_{GB}(t)\le H_{GB}^{\max }\end{array}$$

(5)

where H_GB(t) is the output heating power of GB at time slot t. F_GB(t) is the natural gas consumption of GB at time slot t. $H_{GB}^{\min }$, $H_{GB}^{\max }$ are minimal and maximal output heating power of GB respectively

2.1.3. Waste-Heat Boiler (WH)

The waste-heat boiler (WH) is used to recollect the waste heat to supply cooling and heating load and realize the cascade utilization of energy. The mathematical model of the waste heat boiler can be expressed as follows:

$$\begin{array}{l}{H}_{WH}(t)=H_{GT}(t){\eta }_{WH}\\ H_{WH}^{\min }\le H_{WH}(t)\le H_{WH}^{\max }\end{array}$$

(6)

where H_WH(t) is the output heating power of WH at time slot t. H_GT(t) is the heat generated by the combined supply equipment at time slot t. η_WH is the efficiency of WH. $H_{WH}^{\min },H_{WH}^{\max }$ are minimal and maximal output heating power of WH, respectively.

2.1.4. Heat-Exchange Devices (HX)

The heat-exchange (HX) devices convert the steam heat of the waste heat boiler into the heat energy required by the user, and its mathematical model is as follows:

$$\begin{array}{l}{H}_{HX}(t)=H_h(t){\eta }_{HX}\\ 0\le H_h(t)\le H_{WH}(t)\\ H_{HX}^{\min }\le H_{HX}(t)\le H_{HX}^{\max }\end{array}$$

(7)

where H_HX(t) is the output heat of HX at time slot t. H_h(t) is the heat gained from WH by HX. $H_{HX}^{\min },H_{HX}^{\max }$ are minimal and maximal output heating power of HX, respectively.

2.1.5. Absorption Chiller (AC)

The absorption chiller (AC) can convert steam absorbed from WH into cooling power. The refrigeration cycle of the AC is performed through various circulation processes by utilizing characteristics such as absorption of the solution and generation of refrigerant steam. The mathematical model can be expressed as follows:

$$\begin{array}{l}{H}_{AC}(t)=H_C(t)CO{P}_{AC}\\ 0\le H_C(t)\le H_{WH}(t)\\ H_{AC}^{\min }\le H_{AC}(t)\le H_{AC}^{\max }\end{array}$$

(8)

where H_AC(t) is the output cooling power of AC. COP_AC is the coefficient of performance of AC. H_C(t) is the heat gained from WH by AC. $H_{AC}^{\min },H_{AC}^{\max }$ are minimal and maximal output cooling power of AC, respectively.

2.1.6. Electric Chiller (EC)

The electric chiller (EC) converts the input electrical energy into a cooling power, and the conversion performance coefficient depends on the ratio of the input electrical power to the output cooling capacity. The mathematical model is as follows:

$$\begin{array}{l}{H}_{EC}(t)=CO{P}_{EC}{P}_{EC}(t)\\ P_{EC}^{\min }\le P_{EC}(t)\le P_{EC}^{\max }\end{array}$$

(9)

where H_EC(t) is the output cooling power of EC at time slot t. COP_EC is the coefficient of performance of EC. P_EC(t) is input electric power of EC at time slot t. $P_{EC}^{\min }$,$P_{EC}^{\max }$ are minimal and maximal output electric power of EC, respectively.

2.1.7. Energy Storage System (ESS)

The energy storage system (ESS) can be effectively used to balance the intermittent fluctuation of the output power of renewable power sources and ensure stable operation. The mathematical model of the ESS is modeled as follows [33]:

$$0\le u_t^{ch}+u_t^{dis}\le 1,\forall t$$

(10)

$$\begin{array}{l}{u}_t^{ch}{P}_{ch,\min }^{ESS}\le P_{ch}(t)\le u_t^{ch}{P}_{ch,\max }^{ESS}\\ u_t^{dis}{P}_{dis,\min }^{ESS}\le P_{dis}(t)\le u_t^{dis}{P}_{dis,\max }^{ESS}\end{array}$$

(11)

where $u_t^{ch}$ $u_t^{dis}$ are binary variables, which represent the working state of the ESS at time slot t. The constrains in Equation (11) mean that the ESS can’t charge and discharge at the same time. $P_{ch,\min }^{ESS},P_{ch,\max }^{ESS}$ are minimal and maximal charging power of ESS, respectively. $P_{dis,\min }^{ESS},P_{dis,\max }^{ESS}$ are minimal and maximal discharging power of ESS, respectively.

$$\begin{array}{c}{S}_{ESS}(t+1)=(1-\delta )S_{ESS}(t)+({\eta }^{ch}{P}_{ch}(t)-\frac{1}{{\eta }^{dis}}{P}_{dis}(t))\Delta t\\ S_{ESS}^{\min }\le S_{ESS}(t)\le S_{ESS}^{\max }\end{array}$$

(12)

where S_ESS(t) is the energy stored in the ESS at time slot t. $S_{ESS}^{\min }$,$S_{ESS}^{\max }$ are minimal and maximal energy stored in the ESS, respectively. Δt is the time scale of the operation (one hour in this paper).

2.1.8. Tie Lines

When the energy is not equal to the corresponding load demand. A multi-energy MG can buy or sell electricity to other MGs and the DN constraints to be met are shown as follows:

$$\begin{array}{l}{P}_{grid}^{\min }\le P_{grid}(t)\le P_{grid}^{\max }\\ P_{MG}^{\min }\le P_{MG}(t)\le P_{MG}^{\max }\end{array}$$

(13)

where P_grid(t) represents the electric power interaction between the MG and the DN. P_MG(t) represents the electric power bought or sold to other MGs.

2.2. Distribution Network (DN)

The DN in this paper is a distribution network with renewable energy sources. It can exchange power with the MGs. The mathematical model is shown in Equations (14)–(18):

$$\begin{array}{l}{\displaystyle \sum _{i\in \Gamma }{P}_{G,i}^t}-{\displaystyle \sum _{i\in \Omega }{P}_{L,i}^t}-{\displaystyle \sum _{i\in \Phi }{P}_{MG,i}^t}-\\ U_i^t{\displaystyle \sum _{i=1}^{N_{bus}-1}{U}_j^t(G_{ij}\cos {\theta }_{ij}^t+B_{ij}\sin {\theta }_{ij}^t)}=0\\ {\displaystyle \sum _{i\in \Gamma }{Q}_{G,i}^t}-{\displaystyle \sum _{i\in \Omega }{Q}_{L,i}^t}-{\displaystyle \sum _{i\in \Phi }{Q}_{MG,i}^t}-\\ U_i^t{\displaystyle \sum _{i=1}^{N_{bus}-1}{U}_j^t(G_{ij}\sin {\theta }_{ij}^t-B_{ij}\cos {\theta }_{ij}^t)}=0\end{array}$$

(14)

$$N_{operation}={\displaystyle \sum _{t=2}^T{\displaystyle \sum _{i=1}^{N_{br}}|s_{i,t}-s_{i,t-1}|}}$$

(15)

$$N_{operation}\le N_{\max }$$

(16)

$$U_i^{\min }\le U_i^t\le U_i^{\max }$$

(17)

$$S_{br,k}\le S_{br}^{\max }$$

(18)

where $P_{G,i}^t,P_{L,i}^t,P_{MG,i}^t$ are the active power of power sources, loads and the active power injected into MGs at bus i in DN at time slot t, respectively. Γ, Ω are the sets of nodes connected to the power sources and loads, respectively. Φ is the set of nodes connected to the MGs. U_i^t, U_j^t, U_min, U_max are the voltage of node i and node j at time slot t, the lower and upper boundaries of voltage, respectively. N_bus number of electric buses. G_ij, B_ij is conductance and susceptance of the branch between the node i and node j at time slot t, respectively. θ^t_ij is the phase angle difference between the ith and jth node at time slot t. $Q_{G,i}^t,Q_{L,i}^t,Q_{MG,i}^t$ are reactive power of power sources, loads and active power injected into MGs at bus i in the DN at time slot t, respectively. N_operation is the number of remote-controlled switches actions in a day. s_i,t is the state of the ith switch at time slot t. N_br is the number of remote-controlled switches in the DN. N_max is the upper limit of the switches’ actions. S_br,k, S_br^max represent the apparent power of branch k.

3. Objective Function Formulation

The logic diagram of the model is shown in Figure 3. The problem in this paper can be divided into two levels. The first level defines the inject power of MGs and the states of remote-controlled switches to minimize the power loss and the voltage offset. The model in the second level is performed to gain the minimum cost of the MGs. Each level has its objective functions and constraints. The decisions of the first level would determine the behavior of the second level while the second level would impact the decision of the first level. Thus the relation of two levels is established, the mathematical model of which can be expressed as follows.

$$\mathrm{Level} 1:\{\begin{array}{l}\underset{u}{\min }F(u,v)\\ g(u,v)\le 0\\ h(u,v)=0\end{array}\mathrm{Level} 2:\{\begin{array}{l}\underset{v}{\min }f(u,v)\\ J(u,v)\le 0\\ K(u,v)=0\end{array}$$

(19)

where F(u,v) and f(u,v) are the objective function of the first and second level, respectively. u, v are decision variables of the first and second level, respectively. g(u,v), J(u,v) are inequality constraints respectively. h(u,v) and K(u,v) are equality constraints respectively.

3.1. First Level

The goal of the first level is to minimize the power loss of the DN, as well as the voltage offset. The object function is shown in Equation (20).

$$\{\begin{array}{l}\min ({\displaystyle \sum _{t=1}^T{P}_{loss,t}}+\omega V_{off})\\ P_{loss,t}={\displaystyle \sum _{i=1}^{N_{br}}{s}_{i,t}{r}_i\frac{P_{i,t}^2+Q_{i,t}^2}{{(U_i^t)}^2}},V_{off}={\displaystyle \sum _{t=1}^T\sqrt{\frac{{\displaystyle \sum _{i=1}^{N_{bus}}{[([|U_i^t-1|/0.05]\times 10+|U_i^t-1|/0.05)\times 0.05]}^2}}{N_{bus}}}}\end{array}$$

(20)

$$\begin{array}{l}SW=[s_{1,t},s_{2,t},s_{3,t}\dots s_{N_{br},t}]\\ \Phi =[P_{MG,1}^t,P_{MG,2}^t,P_{MG,3}^t\dots P_{MG,N_{MG}}^t]\end{array}$$

(21)

$$N_f=N_{Bus}-1$$

(22)

The first term in Equation (20) represents the total network power loss of the DN. P_loss,t is the power loss of the DN at time slot t. The second term in Equation (20) represents the node voltage offset, which can be regarded as a penalty function. V_off is voltage offset in a day [26]. ω is a coefficient indicating the degree of importance. r_i is the resistance of the branch i. P_i,t, Q_i,t are the active power and reactive power of the end node of the branch i.

The controlled variables of the reconfiguration of the DN are shown in Equation (21), SW represents the set that controls the open/closed state of the switch on each branch, “0” represents the switch is open while “1” represents the switch is closed, the inequality constraint of SW is shown in Equation (15). Φ represents the set that controls power interaction between the MGs and the DN. The relevant constraints of Φ are shown in Equations (13) and (14). N_MG is the number of the MGs connected to the DN. The operating status of the DN can be adjusted by changing the values or states of the controlled variables above when the DN is running.

In order to facilitate the maintenance of the distribution network and reduce the short-circuit current when the distribution network fails, the topology of the distribution network is usually radial. The constraint is shown in Equation (22). N_f is the number of branches in the DN.

3.2. Second Level

The object function is shown in Equation (23). It is a mixed-integer linear programming (MILP) problem to search the optimal operating mode of the devices. The object is to minimize the operation cost of the MGs, which includes the cost of natural gas, cost of power interaction with the DN, cost of power interaction with the other MGs, and devices’ maintenance cost in the MGs.

$$\min C_{operation}={\displaystyle \sum _{i=1}^{N_{MG}}(C_{gas,i}+C_{grid,i}+C_{ex,i}+C_{mc,i}})$$

(23)

where C_operation represents the total operation cost of MGs. C_gas,i represents the cost of natural gas of ith MG. C_grid,i is the power purchase cost of ith MG. C_ex,i represents the power interaction cost with other MGs of ith MG. C_mc,i represents the device maintenance cost of ith MG. The details of each cost are shown below:

$$C_{gas,i}={\displaystyle \sum _{t=1}^T{c}_{gas,t}(F_{GT}(t)+F_{GB}(t))\Delta t}$$

(24)

$$C_{grid,i}={\displaystyle \sum _{t=1}^T{c}_{grid,t}{P}_{grid}(t)\Delta t}$$

(25)

$$C_{ex,i}={\displaystyle \sum _{t=1}^T{c}_{ex,t}{P}_{MG}(t)\Delta t}$$

(26)

$$C_{mc,i}={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{k=1}^{N_{dv}}{c}_k(t)P_k(t)}}\Delta t$$

(27)

where c_gas,t represents the gas price at time slot t. c_grid,t is electricity price at time slot t. c_ex,t is power interaction price between MGs at time slot t. c_k and P_k(t) are unit operation and maintenance cost and power generated or consumed by GT, GB, WH, HX, AC, EC, photovoltaic array (PV), wind turbine (WT), ESS. N_dv is the number of devices in the MG.

The multi-energy MGs contain a variety of constraints, including device output constraints and bus balance constraints. The device output constraints and power interaction constraints have been introduced in Section 2. As illustrated in Figure 2, the bus balance constraint in each MG can be written as:

• Electric bus balance

  $$\begin{array}{l}{P}_{PV}(t)+P_{WT}(T)+P_{GT}(t)+P_{grid}(t)+\\ P_{dis}(t)+P_{MG}(t)=P_{Load}(t)+P_{ch}(t)+P_{EC}(t)\end{array}$$

  (28)
• Cooling load bus balance

  $$P_{EC}(t)CO{P}_{EC}+H_{AC}=H_{cooling}(t)$$

  (29)
• Heating load bus balance

  $$H_{HX}(t)+H_{GB}(t)=H_{heating}(t)$$

  (30)
• Steam bus balance:

  $$H_C(t)+H_h(t)=H_{WH}(t)$$

  (31)

  where P_Load(t), H_cooling(t), H_heating(t) represent the electrical load, cooling load and heating load, respectively.

4. Solving Method

From the above, the problem is divided into two levels. The first level mainly solves the problem of DN reconfiguration and determination of the power injected from the MGs. The second level solves the optimal dispatch problem of interconnected MGs.

4.1. First Level

The problem can be solved by a particle swarm optimization (PSO) algorithm. PSO originates from the study of the predation behavior of birds. The basic idea of the particle swarm optimization algorithm is to find the optimal solution through cooperation and information sharing between individuals in the group. The variables are abstracted as particles (points) without mass and volume and extend to N-dimensional space. The position of particle k in N-dimensional space at ith iteration is represented as vector x_i,k, the flight speed at ith iteration is represented as v_i,k. Each particle has a fitness value, determined by the objective function. Pbest = (pbest₁, pbest₂ …, pbest_M) is the set of the best position of each particle. M is the number of particles. The best position the particles have ever found and the current positions can be seen as the particle’s own flight experience. In addition, each particle also knows that the best position gbest_i (the best value in Pbest). This can be regarded as the experience of particle companions. Particles update their own flight speed and position by comparing the fitness value fit(x_i,k) with the corresponding personal best fitness value in Pbest and gbest. The detailed procedures are shown as follows:

• Step 1: Initialization, acquire the system data, including line and bus parameters, generation information and states of switches. Generate a particle swarm with a sample number of M. Initialize the position x_0,k (including the vectors SW and Φ in Equation (21)) and speed v_0,k of each particle. Set pbest_k = x_0,k, evaluate the fitness value fit(x_0,k) of each particle using the function in Equation (20) by power flow calculation. Find the particles with the minimum fitness value of all particles and set its position to gbest₀.
• Step 2: For ith iteration. Update the position and speed of each particle using Equations (32) and (33), evaluate the fitness value fit(x_i,k) of each particle using the function in Equation (20).

  $$v_{i,k}=\{\begin{array}{l}w{v}_{i,k}+c_1\times ran{d}_1()\times (pbes{t}_k-x_{i,k})+\dots \\ \dots c_2\times ran{d}_2()\times (gbest-x_{i,k})\\ v_{\max }\hspace{1em}{v}_{i,k}>v_{\max }\\ v_{\min }\hspace{1em}{v}_{i,k}<v_{\min }\end{array}$$

  (32)

  $$x_{i,k}=\{\begin{array}{l}{x}_{i,k}+v_{i,k}\hspace{1em}{x}_{\min }\le x_{i,k}\le x_{\max }\\ x_{\max },v_{i,k}=-v_{i,k}\hspace{1em}{x}_{i,k}>x_{\max }\\ x_{\min },v_{i,k}=-v_{i,k}\hspace{1em}{x}_{i,k}<x_{\min }\end{array}$$

  (33)

  where rand₁(), rand₂() represents random numbers between 0 and 1. c₁, c₂ are learning factors, w is the inertia weight. If the particles are discrete variables like the state of the ith switch s_i,t in the vector SW, the following rules are required when performing speed updates:

  $$x_{i,k}=\{\begin{array}{l}1\hspace{1em}ran{d}_{\mathrm{i},k}()<S(v_{i,k})\\ 0\hspace{1em}ran{d}_{i,k}()\ge S(v_{i,k}))\end{array}$$

  (34)

  $$S(v_{i,k})=\frac{1}{1+e^{-v_{i,k}}}$$

  (35)

  where S(V_k) is a Sigmoid function, which is used to select the probability of discrete variable positions. Generate a random number rand_i,k() between 0 and 1 and compare it with S(v_i,k) to determine the value of the binary variable x_i,k.
• Step 3: For each particle, compare the fitness value fit(x_i,k) with fit(pbest_k), update pbest_k to the position with lower fitness value.
• Step 4: For all particles, select the particle who has the lowest fitness value. Update gbest_i to its position.
• Step 5: Judge whether the result meets the end condition. If the iteration number has reached the upper limit or the fitness value convergence, the calculation stops. Otherwise, return to Step 2.

4.2. Second Level

The second level problem in Equation (23) is essentially a 0–1 mixed-integer linear programming problem (MILP), whose general form is as follows:

$$\begin{array}{l}\min f(x)\\ s.t.A_{eq}=b_{eq}\\ Ax\le b\\ x_{Integer}=\{0,1\}\\ lb<x_{real}<ub\end{array}$$

(36)

In Equation (36), x is a variable matrix that changes with time and is updated once an hour. This problem can be solved by the MILP tool in MATLAB. The details are illustrated as follows:

• Step 1: Aquire the power interaction command vector Φ from the first level, and calculate the solution of the MGs according to the power command.
• Step 2: Judge whether there is a solution for the optimization of the second level. If the solution exists, feedback the information to the first level. If the solution does not exist, it needs to change the operation method (regardless of upper DN) and re-solve the first level. The flow chart of the solution process above is shown in Figure 4.

5. Simulation Result

5.1. Case Set Up

For the aim of verifying the validity of the proposed model, this paper established a test system including IEEE 33-bus distribution network (voltage reference value is 12.66 kV, 37 remote-controlled switches(K1-K37) with three wind farms, three photovoltaic power stations, and three multi-energy MGs, the upper limit of control actions is set to 24), as shown in Figure 5. The operation data of distributed generators and MGs is based on history. The test system was set up on Thinkpad T440p, 2.90 GHz with 8 cores and developed using MATLAB 2017b with toolbox CPLEX 12.9.

The voltage is subjected to the constraint of 0.93 and 1.07. Three wind farms are connected to node 8, 10, 30. Three photovoltaic power stations are connected to node 13, 16, 18 (forecast 24-h power generation is shown in Figure 6). Three microgrids are connected to the 22, 25, 33, respectively, the load flow of MGs is not considered in this paper [29].

MG1 is CCHP-based MG, while MG2 and MG3 are CHP-based MG and CCP-based MG, respectively. The parameters of devices in MGs are given in

Table 2. Device Parameters in MGs.

Parameters	Values	Parameters	Values
ηWH	0.8	$Q_{\max }^{HX}$	2000 kW
ηGB	0.9	$Q_{\max }^{AC}$	1500 kW
ηHX	0.9	$Q_{\max }^{EC}$	1000 kW
COPAC	1.2	$P_{ch,\max }^{ESS}$	100 kW
COPEC	4	$P_{dis,\max }^{ESS}$	100 kW
ηch/ηdis	0.96	$SO{C}_{\max }^{ESS}$	0.9
δ	0.02	$SO{C}_{\min }^{ESS}$	0.2
$P_{\max }^{GT}$	1500 kW	$P_{\max }^{grid}$	2000 kW
$Q_{\max }^{GB}$	2000 kW	$P_{\max }^{MG}$	1000 kW

. The price of natural gas per cubic meter is ¥2.2. The heat value of natural gas is 9.69 kWh/m³. The prices of power interaction between DN and MGs are shown in

Table 3. Purchasing price of MGs and DN.

Time Slot	Price between DN and MGs	Price between MGs
23:00–7:00	0.17	0.12
7:00–8:00 17:00–18:00	0.49	0.37
8:00–11:00 18:00–23:00	0.83	0.65

, which can show that the peak-to-valley price method is adopted. To protect the benefit of electric power enterprises, the price of electricity traded between MGs cannot be higher than that of trading between microgrid and DN.

To verify the rationality of the proposed method, there are three comparative cases to show the impact of the proposed method, given in

Table 4. Of Different Cases.

	MGs Optimal Dispatch	Operation of DN	DN Reconfiguration
Case1	√	×	×
Case2	√	√	×
Case3	√	√	√

.

• Case 1: Optimal dispatch of MGs. MGs purchase or sell power to DN and buy natural gas from the gas network without consideration of the DN [15].
• Case 2: Optimal dispatch of MGs and operation of DN (considering power loss) without consideration of the DN reconfiguration [29].
• Case 3: The method proposed in this paper. The bi-level optimal dispatch considering operations of MGs and DN considering reconfiguration, through controlling the switch on each tie line in DN, operating condition of devices and power interaction of the PCC (point of common coupling, the connection point between the MG and the DN).

5.2. Result Analysis

5.2.1. First Level

Table 5. Daily power loss and voltage offset.

Result	Case 1	Case 2	Case 3
Ploss/MW	2.5310	2.1543	1.8015
Voff/pu	12.6580	9.9851	9.0416

records the daily power loss and voltage offset of the distribution network. The daily power loss and voltage offset in Cases 2 and 3 are significantly lower than those in Case 1. Compared with Case 2, the daily power loss and voltage offset of Case 3 reduces more. It can be seen that under the framework of the coordinated dispatch of the MGs and the DN and, taking the topology of the DN into consideration, the power loss and voltage offset of the distribution network can be effectively reduced. Figure 7 illustrates the lowest voltage amplitude of all electrical buses at each time slot in a day. As can be seen from the Figure 7 that in all cases, the lowest voltage is within the limitation. The amplitude of the lowest voltage in Case 2 has the maximum fluctuation value and the voltage reaches lower bound at time slot 1, 6, 20, 22 and 24 compared with Case 1 and Case 3. Combined with data shown in Table 5, although the total voltage offset decreases, the minimum amplitude of voltage in Case 2 is at the risk of exceeding the lower limit at part of the time slots. In Case 3, through the DN reconfiguration, the voltage amplitude offset further reduces. At the same time, the fluctuation of the minimum value of all electrical bus voltage amplitudes reduces and does not fall to the lower limit at each time slot, leaving a certain safety margin. It can be seen that the reconfiguration of the DN can effectively restrain the voltage amplitude from exceeding the limit. which also proves the necessity of distribution network reconfiguration.

For Case 3, Figure 8 demonstrates the states of the remote-controlled switches in the DN during the day (red means closed while green means open). By adjusting the states of the branch switch at multiple time slots (1–2, 6–8, 11–12, 18–24), the topology of the DN is guaranteed to be radial and the operating state of the distribution network is optimized to reduce the DN network loss and voltage offset and also improve the condition of power interaction between the DN and the MGs. The power interaction on tie lines between the DN and the three MGs is shown in Figure 9. Figure 9a–c represents the power interaction between the DN and the three MGs in Cases 1–3, respectively. In Case 1, when only optimal dispatch of the MGs is considered, the power interaction between the MGs and the DN is relatively random, which will adversely affect the operation of the DN.

In Case 2, when the coordinated optimization of the MGs and the DN is considered but the reconfiguration problem of the DN is not considered, the power interaction between some MG1 and the DN continues to be large, which is unfavorable to the power flow distribution of the DN and will increase excessive operating cost of the MG.

In Case 3, through the reconfiguration of the DN, the power interaction distribution between the DN and each MG is more balanced, and the randomness of power interaction greatly reduces, which helps the MGs to better cooperate with DN operation, reducing power loss in the DN.

The trend of power interaction of MG2 and MG3 is similar to that in Case 3; the power purchased by the MGs to the DN shows a trend of increasing first and then decreasing, indicating that the MGs play a positive role in the renewable energy accommodation in the DN, which is affected by the reconfiguration of the DN.

5.2.2. Second Level

Figure 10, Figure 11 and Figure 12 present the load balance of MGs in Case 1. The MG1 makes full use of the cooling and heating power. There is no cold and heat power curtailment. In the beginning, the electricity price is low and load is light; to consume the energy generated by the renewable power sources, the MG1 tends to purchase electricity from the DN and the MG2, 3 with surplus power, and the GT is in the shutdown state. When the GT is not working, the GB and the EC are used to meet the cooling and heating load demand in the MG1, respectively. As the electricity price and the load increase, the GT begins to generate electricity, and the WH transforms the waste heat from the GT to the steam. AC and HX absorb the steam to satisfy the cooling and heating load. The GB and the EC supplement insufficient demand.

From 12:00 to 18:00, since all the steam transformed by the WH is to meet the demand of the heating load, the cooling load is fully satisfied by the EC. The battery is discharging during the peak electricity price period (9:00–11:00, 19:00–21:00). Similarly, the battery is charging in other periods (4:00–8:00,15:00–18:00), and the maximum charge-discharge power is up to 100 kW.

In MG2, as can be seen from Figure 12, renewable power sources generate more electricity than the MG needs, so MG2 tends to sell electricity to DN and other MGs with insufficient electricity. The GT only works at the peak electricity price period (9:00–11:00,19:00–23:00) to generate waste heat to meet the cooling load demand while selling excess electricity to DN and other MGs. The rest time period is satisfied by the EC. The power exchange only happens among MGs. The charging and discharging state of ESS is also related to electricity price and load level.

MG3 mainly sells electricity to other MGs and DN, and the operating state of GT and ESS (maximum charge-discharge power 100 kW) is similar to that of MG1. The waste heat generated by the GT is supplied to the HX through the WH, and the GT is used for heating when the GT is not working or the HX is unable to meet the heating load.

In Case 1, the MGs only consider their own operation object, regardless of the working condition of the DN. It is obvious that the total operating cost of the MGs is the lowest. However, the power loss and voltage offset of the distribution network have not been optimized in this case. Compared with Case 2 and Case 3 in Table 5, the above two values in Case 1 are higher, which is not permitted in the actual situation, since it will worsen the operating status of the DN and increase the failure rate. It is not beneficial to the global optimal operation of the whole system. The goal of Case 2 and Case 3 is to optimize the MGs’ dispatch under the consideration of the DN. The MGs should take the power exchange instructions into consideration not just considering their own economic benefits. In this case, although the operation cost of the MGs has increased, the power loss and voltage offset of the DN reduce significantly, which is very effective for the optimization of the whole system. The economic benefits of the MGs also need to be guaranteed, and the increase of the cost should not be too great. The reconfiguration of the DN can optimize the operation of the DN while increasing the operating cost of the microgrids as little as possible and control it within a reasonable range.

Table 6. MGs in different cases.

	Case 1	Case 2	Case 3
MG1	28,173.95	33,338.42	31,646.77
MG2	−2773.21	−1986.90	−3012.89
MG3	7558.35	9433.32	8802.69
Total	33,059.09	40,784.84	37,436.03

shows the cost of the three MGs in all cases.

As can be seen from Table 6, MG2 makes profits while MG1 and MG3 generate cost. Compared with Case 1, to consider the operation of the DN, in Case 2 and Case 3 the devices in the MGs adjust the output according to the command, as a result increasing the total system cost.

5.2.3. Summary

From the data above, we can see that the total power loss and voltage offset decline remarkably in Case 2 (17.49% and 21.12%) and Case 3 (28.82% and 28.58%) compared with Case 1. To satisfy the power demand of the MGs, the dispatch approach in Case1 will consume more fluctuations of power interaction, which causes the result of more power loss and voltage offset. Moreover, the dispatch method in Case 1 could not accommodate the power from the PVs and WTs in the DN well, which can be solved in the Case 2 and Case 3. Moreover, the curve of the power interaction in Case 3 (shown in Figure 9c) is more similar to the forecast power of the renewable power sources, which reflects the effectiveness of the DN reconfiguration.

The total cost of the MGs increases by 23.31% and 13.23%, respectively, in Case 2 and Case 3. Under the requirement of the DN, the devices in the MGs will adjust their working conditions. Although the operation cost of the MGs increase, the coordinate optimization method can effectively reduce the power loss, voltage offset, improve the power flow distribution and obtain a better global economic optimal value of the whole DN system. In addition, the network reconfiguration of the DN can realize further optimization.

6. Conclusions

In this paper, a bi-level optimization method is presented. The purpose of this method is to improve both the working condition of the DN and MGs. The first level is to minimize power loss and voltage offset in the DN based on the PSO algorithm containing discrete variables, considering reconfiguration of the DN and the power interaction between DN and MGs. The second level solves a MILP problem to control the operation modes of devices in CCHP-based MGs to minimize the total cost of MGs. The simulation results in all cases demonstrate that:

(1)

The power loss, voltage offset, and power interaction on tie lines have noticeable reductions based on the proposed method, reflecting that it is necessary to consider both the benefit of DN and MGs.

(2)

Compared with [25], the consideration of the structure of the DN can keep the voltage of partial nodes away from extreme conditions, while further reducing DN power loss and operating costs of the MGs.

(3)

The MGs provide support for the DN, and in the appropriate cases, some of the benefits can be sacrificed to meet the dispatch requirements under the framework proposed in this paper.

(4)

This paper only considers the electrical distribution network. Further research can focus on the effect of taking gas networks and thermal networks into consideration.