An Enhanced Second-Order Cone Programming-Based Evaluation Method on Maximum Hosting Capacity of Solar Energy in Distribution Systems with Integrated Energy ^†

Abstract

In order to adjust to the change of the large-scale deployment of photovoltaic (PV) power generation and fully exploit the potentialities of an integrated energy distribution system (IEDS) in solar energy accommodation, an evaluation method on maximum hosting capacity of solar energy in IEDS based on convex relaxation optimization algorithm is proposed in this paper. Firstly, an evaluation model of maximum hosting capacity of solar energy for IEDS considering the electrical-thermal comprehensive utilization of solar energy is proposed, in which the maximization of PV capacity and solar collector (SC) capacity are fully considered. Secondly, IEDS’s potential in electricity, heat, and gas energy coordinated optimization is fully exploited to enhance the hosting capacity of solar energy in which the electric distribution network, heating network, and natural gas network constraints are fully modeled. Then, an enhanced second-order cone programming (SOCP)-based method is employed to solve the proposed maximum hosting capacity model. Through SOCP relaxation and linearization, the original nonconvex nonlinear programming model is converted into the mixed-integer second-order cone programming model. Meanwhile, to ensure the exactness of SOCP relaxation and improve the computation efficiency, increasingly tight linear cuts of distribution system and natural gas system are added to the SOCP relaxation. Finally, an example is given to verify the effectiveness of the proposed method. The analysis results show that the maximum hosting capacity of solar energy can be improved significantly by realizing the coordination of an integrated multi-energy system and the optimal utilization of electricity, heat, and gas energy. By applying SOCP relaxation, linearization, and adding increasingly tight linear cuts of distribution system and natural gas system to the SOCP relaxation, the proposed model can be solved accurately and efficiently.

1. Introduction

Solar energy has been widely deployed in the world with its huge potential in reducing carbon emissions. There are mainly two kinds of technologies of harvesting solar energy: photovoltaic (PV) and solar collector (SC) [1,2]. However, the output power of PV or SC from solar energy is intermittent and susceptible to the meteorological conditions, which will pose new challenges to the operation of distribution network (DN). For example, the high penetration of distributed PVs will inevitably increase the probability of current spillage and voltage violation, thus limiting the PV integration capacity in distribution networks [3]. It is of great significance to reasonably evaluate the maximum hosting capacity of solar energy for fully efficient utilization of solar energy resources and ensuring the safe operation of the energy systems at the same time.

Nowadays, many references have carried out research on the distributed renewable energy sources hosting capacity problems. Most of the research findings primarily concentrate on hosting capacity of three ways, as shown in

Table 1. Primarily concentrating on hosting capacity of three ways.

Ways	Features
MG [4,5,6]	The integration and coordination of distributed generators, distributed storage systems, controllable loads, etc.
ADN [7,8,9,10,11]	The network reconfiguration, the power factor control strategy, the flexible interconnection technology, etc.
IEDS [12,13,14,15]	The optimal utilization of multiple types of energies.

: (1) based on microgrid (MG) [4,5,6], (2) based on active distribution network (ADN) [7,8,9,10,11], and (3) based on integrated energy distribution system (IEDS) [12,13,14,15].

Potential solutions can be adopted to improve the distributed renewable energy sources hosting capacity through MG to realize the integration and coordination of distributed generators, distributed storage systems, and controllable loads [4,5,6]. Based on ADN, the network reconfiguration, the power factor control strategy, the reactive power compensation devices, the flexible interconnection technology, etc., can be used to improve the distributed renewable energy sources hosting capacity [7,8,9,10,11]. As the integration and coordination of a multi-energy system can realize optimal utilization of multiple types of energies, for example, electricity, heat, and gas, and can provide a “buffer” to accommodate more distributed renewable energy, IEDS has been considered to be able to provide more effective ways to enhance the distributed renewable energy sources hosting capacity [12,13,14,15].

The hosting capacity of solar energy in an electrical–thermal integrated energy distribution system can be improved by installing electrical boiler (EB), gas boiler (GB), SC, and heat storage tanks [2,16,17,18,19]. Not only the hosting capacity of solar energy in the form of heat directly, but also the hosting capacity of PV can be improved by the installation of SC [2,16,17,18]. The hosting capacity of solar energy in an electrical–natural gas integrated energy distribution system can be improved by the power to gas (P2G) technology which can eliminate the surplus power generation on a large scale by using the electricity coming from renewable energy to generate natural gas [20,21]. As IEDS are characterized by multi-energy complementarity and coordinated utilization, it is greatly helpful to improve the hosting capacity of solar energy by fully utilizing the potential of distribution system, natural gas system, and heating system.

Generally, the maximum hosting capacity problem of solar energy in IEDS is a highly nonconvex nonlinear programming (NLP) problem which may take a lot of time to solve due to multiple locally optimal points [22,23,24], when considering the electrical power flow equations, the pipeline of natural gas network (NGN), and heating network (HN) equations. Aiming at the convexification problems in IEDS, researchers have proposed many methods to solve the electrical power flow, natural gas flow, and thermal power flow problem [25,26,27,28,29,30]. For example, the second-order cone programming (SOCP) relaxation method was used to solve the optimal electrical power flow problem [25,26]. The piecewise linearization or SOCP relaxation method was used to solve the optimal natural gas flow problem [27,28], and the linearization method was used to convert the original nonlinear heating network model into a linear heating network model [29,30]. Based on the above methods, the original nonconvex NLP model can be converted into a mixed-integer second-order cone programming (MISOCP) model that can be solved easily by global optimization solvers to obtain the globally optimal solution and reduce the computation time [31]. However, the exactness of SOCP relaxation is always related to the selected objective function. In order to ensure the exactness of SOCP relaxation for various objective functions, increasingly tight linear cuts can be added to the SOCP relaxation problem [32].

In order to make full use of the potential of distribution system, natural gas system, and heating system to utilize more solar energies, an enhanced SOCP-based method is developed to evaluate the maximum hosting capacity of solar energy in electrical–natural gas–thermal IEDS. The major contributions of this paper are summarized as follows:

(1)

Model: An optimization model of the maximum hosting capacity evaluation of solar energy in IEDS is proposed, in which the maximization of PV capacity and solar collector (SC) capacity are fully considered.

(2)

Mechanism: IEDS’s potential in multi-energy coordinated optimization is fully exploited to enhance the hosting capacity of solar energy in which the electric distribution network, heating network, and natural gas network constraints are fully modeled.

(3)

Method: An enhanced-SOCP-based solving method is developed to solve the proposed maximum hosting capacity model, which can output a satisfactory solution and reduce the computation time.

The remainder of this paper is structured as follows: Section 2 introduces the optimization model for the maximum hosting capacity evaluation of solar energy in an IEDS with PV and SC. In this section, the distribution network, the heating network, and the natural gas network constraints are fully considered to solve the problem of current spillage and voltage violation, and, meanwhile, to ensure the safe operation of the energy systems. Section 3 develops the solution methods. By applying SOCP relaxation for the distribution and natural gas system model and linearization for the heating system model, the original NLP model of the maximum hosting capacity is converted into an MISOCP model, which can be solved efficiently to obtain the globally optimal solution of a large problem. In order to ensure the exactness of SOCP relaxation and improve the computation efficiency, increasingly tight linear cuts of distribution system and natural gas system are added to the SOCP relaxation. Section 4 presents the optimization results of two cases to verify the effectiveness of the proposed model and method. Section 5 concludes the paper.

2. Model of Maximum Hosting Capacity

2.1. Objective Function

The objective function of the maximum hosting capacity model is composed of four components: maximized PV capacity, maximized SC capacity, minimized electric power loss, and minimized heating power loss [2], as shown in (1):

$$\max f={\varphi }_1{f}_{PV}+{\varphi }_2{f}_{SC}-f_{loss1}-f_{loss2}$$

(1)

where f_PV is the output power of the PV, f_SC is the output power of the SC, f_loss1 is the distribution system power loss, and f_loss2 is the heating system power loss. The quantities ϕ₁ and ϕ₂ are the coefficient of the output power of the PV and the coefficient of the output power of the SC, respectively. The distribution system power losses f_loss1 and the heating system power losses f_loss2 are formulated in (2)–(5), respectively.

$$f_{PV}={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i=1}^N{\eta }_{\mathrm{PV}-\mathrm{P}}{\eta }_{i,t}^{\mathrm{SE}}S{E}_i^{PV}}}$$

(2)

where T is the total periods of time horizon; N is the total number of nodes of the DN; ${\eta }_{i,t}^{\mathrm{SE}}$ is the actual irradiance intensity of the system at node i; ${\eta }_{\mathrm{PV}-\mathrm{P}}$ is the efficiency of PV; $S{E}_i^{PV}$ is the alternative installation area of PV at node i.

The solar collector (SC) is a device that can collect solar radiation from the sun to heat the water for users. The evacuated tube solar collector [18,19] is considered in this paper.

$$f_{SC}={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i=1}^N{\eta }_{\mathrm{SC}}{\eta }_{i,t}^{\mathrm{SE}}S{E}_i^{SC}}}$$

(3)

where ${\eta }_{\mathrm{SC}}$ is the efficiency of SC; $S{E}_i^{\mathrm{SC}}$ is the alternative installation area of SC at node i.

$$f_{loss1}={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{ij\in {\mathsf{\Omega }}_{\mathrm{b}}}^{}{r}_{ij}{i}_{ij,t}^{}}}$$

(4)

where Ω_b is the set of all branches of the DN; r_ij is the resistance of branch ij; i_ij,t is the current magnitude square of branch ij at time t.

$$f_{loss2}={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{ij\in {\mathsf{\Omega }}_C}^{}{\eta }_{ij}(H_{ij,t}^{}+H_{ji,t}^{})}}$$

(5)

where Ω_c is the set of all pipelines of the HN; H_ij,t and H_ji,t are heat power of the pipeline from node i to node j and heat power of the pipeline from node j to node i in the HN at time t, respectively; η_ij is heat loss ratio of pipeline ij in the HN.

In addition, the total alternative installation area of PV and SC in the system is limited, constrained by the geographic location.

2.2. Constraints of the Distribution Network

In this paper, the distflow branch model [33] is used to model the distribution network. The electric distribution network constraints include the active and reactive power balance constraints, the Ohm’s law constraints, the constraints of the relationship between current, voltage, and power, the security constraints of DN, etc.

2.3. Constraints of the Heating Network

This paper adopts the heating network model proposed in [29], which is formulated as follows:

$$\left\{\begin{array}{l}{H}_{i,t}+{\displaystyle \sum _{j\in i}{\eta }_{ji}{H}_{ji,t}}={\displaystyle \sum _{j\in i}{H}_{ij,t}}\\ {\eta }_{ji}=1-\delta l_{ji}\end{array}\right.$$

(6)

$$\left\{\begin{array}{l}0\le H_{ij,t}\le u_{ij,t}\overline{H_{ij}}\\ 0\le H_{ji,t}\le u_{ji,t}\overline{H_{ji}}\\ u_{ij,t}+u_{ji,t}=1\end{array}\right.$$

(7)

where j is the set of nodes which can be directly connected to node i in the HN; H_i,t is the heat power injected into node i at time t; δ is the heat loss ratio per unit length in the HN; l_ji is the length of the pipeline ij in the HN; u_ij,t is equal to 1, and u_ji,t is equal to 0 if the heat power direction of the pipeline ij is from node i to node j in the HN at time t. $\overline{H_{ij}}$ is the maximum heat power of the pipeline ij in the HN, whose detailed expression can be referred to in [29].

2.4. Constraints of the Natural Gas Network

This paper introduces a 0-1 integer variable to represent the pipeline flow direction. It is an improved model of the natural gas network model proposed in [23], which is formulated as follows:

$$p_{l,i,t}-p_{l,j,t}=s_{l,ij,t}{F}_{l,ij}{q}_{l,ij,t}^2$$

(8)

$${\displaystyle \sum _{j\in i}{q}_{l,ij,t}}+q_{i,t}=0$$

(9)

$$s_{l,ij,t}=\left\{\begin{array}{l}1\\ 0\\ -1\end{array}\right.\begin{array}{cc}& \begin{array}{l}(p_{l,i,t}-p_{l,j,t})>0\\ (p_{l,i,t}-p_{l,j,t})=0\\ (p_{l,i,t}-p_{l,j,t})<0\end{array}\end{array}$$

(10)

$$p^{\min }\le p_{l,i,t}\le p^{\max }$$

(11)

where p_l,i,t and p_l,j,t are the begin node pressure of node i and the end node pressure of node j of gas pipeline l at time t, respectively; s_l,ij,t is the 0-1 integer variable to represent the pipeline flow direction; q_l,ij,t is the gas pipeline l volume flow at time t; q_i,t is the volume flow injected into node i in the NGN at time t; F_l,ij is the resistance coefficient of gas pipeline l, whose detailed expression can be referred to in [23]; p^max and p^min are the allowable maximum pressure and minimum pressure in the NGN.

2.5. Constraints of the Energy Station

This paper proposes a standardized matrix modeling method of the energy station proposed in [34], which is formulated as follows:

$$\left\{\begin{array}{l}\mathit{C}=\left[{\mathit{C}}_1,\hspace{0.33em}{\mathit{C}}_2,\hspace{0.33em}\cdots ,\hspace{0.33em}{\mathit{C}}_N\right]\\ \mathit{C}\mathit{I}=\mathit{L}\\ 0\le \mathit{I}\le {\mathit{I}}_{\max }\end{array}\right.$$

(12)

where C is the energy station energy conversion matrix, C_n is the energy converters’ conversion vector; I is the energy input power vector of the energy converters; L is the load vector of the energy station; I_max is the capacity vector of the energy converters.

The energy relationships between the energy station and each subsystem are shown in (13)–(16):

$$P_{j,t}^{line}=-P_{i,t}$$

(13)

$$Q_{j,t}^{line}=-Q_{i,t}$$

(14)

$$G_{j,t}^{line}=\frac{-{\mathrm{G}}_{\mathrm{CV}}{q}_{i,t}}{3.6}$$

(15)

$$H_{j,t}^{line}=-H_{i,t}$$

(16)

where P_i,t and Q_i,t are the active and reactive power injected into node I in the DN, respectively; $P_{j,t}^{line}$, $Q_{j,t}^{line}$, $G_{j,t}^{line}$, and $H_{j,t}^{line}$ are the active power, the reactive power, the gas power, and the heat power injected into the jth energy station, respectively; G_CV is the gross calorific value of gas.

3. Solution Methodology

By using SOCP relaxation and linearization [22], the DN model can be converted into a second-order cone model. Following the same path, in order to apply convex relaxations to the NGN model, the nonlinear constraint (8) is preprocessed to facilitate the convexification by adding the auxiliary variables $x_{l,t}$, $y_{l,t}$, M, ${\overline{v}}_{l,t}$, and ${\underset{¯}{v}}_{l,t}$. The detailed formulas are as follows:

$$\left\{\begin{array}{l}{p}_{l.i.t}\le {\overline{v}}_{l,t}\\ p_{l.j.t}\le {\overline{v}}_{l,t}\\ p_{l.i.t}\ge {\overline{v}}_{l,t}-M·(1-x_{l,t})\\ p_{l.j.t}\ge {\overline{v}}_{l,t}-M·(1-y_{l,t})\end{array}\right.$$

(17)

$$\left\{\begin{array}{l}{p}_{l.i.t}\ge {\underset{¯}{v}}_{l,t}\\ p_{l.j.t}\ge {\underset{¯}{v}}_{l,t}\\ p_{l.i.t}\le {\underset{¯}{v}}_{l,t}-M·(1-y_{l,t})\\ p_{l.j.t}\le {\underset{¯}{v}}_{l,t}-M·(1-x_{l,t})\end{array}\right.$$

(18)

$$\left\{\begin{array}{l}-M·y_{l,t}\le q_{l,ij,t}\le M·(1-y_{l,t})\\ -M·y_{l,t}\le p_{l,i,t}-p_{l,j,t}\le M·(1-y_{l,t})\end{array}\right.$$

(19)

$$\left\{\begin{array}{l}{x}_{l,t}+y_{l,t}\ge 1\\ x_{l,t},y_{l,t}\in \left\{0,1\right\}\end{array}\right.$$

(20)

where $y_{l,t}$ is equal to 1 if $p_{l.i.t}$ is less than or equal to $p_{l.j.t}$ in the NGN at time t; ${\overline{v}}_{l,t}$ is the larger of $p_{l.i.t}$ and $p_{l,j,t}$, ${\underset{¯}{v}}_{l,t}$ is the smaller of $p_{l.i.t}$ and $p_{l,j,t}$; M is an arbitrarily large positive number that is not infinite. Then, defining new variables $F_{l,ij}^{\prime }=\sqrt{F_{l,ij}}$, the natural gas flow (8) can be converted into (21):

$$({\overline{v}}_{l,t}-{\underset{¯}{v}}_{l,t})={(F_{l,ij}^{\prime }{q}_{l,ij,t}^{})}^2$$

(21)

Then, (21) can be further relaxed to the following second-order cone constraint (22):

$${‖{\left[1-({\overline{v}}_{l,t}-{\underset{¯}{v}}_{l,t})\hspace{1em}\hspace{1em}2F_{l,ij}^{\prime }{q}_{l,ij,t}^{}\right]}^{\mathrm{T}}‖}_2\le 1+({\overline{v}}_{l,t}-{\underset{¯}{v}}_{l,t})$$

(22)

However, the exactness of SOCP relaxation is greatly related to the selected objective function. Two indexes are defined to quantify the relaxation deviation. The maximum SOCP relaxation deviation of the distribution system is defined in [32], as shown in (23):

$$r{\mathrm{Gap}}_{\mathrm{DN}}(x)={\displaystyle \sum _{ij\in {\mathsf{\Omega }}_b}{r}_{ij}}\left(i_{t,ij,k}{u}_{t,i,k}-{\left(P_{t,ij,k}\right)}^2-{\left(Q_{t,ij,k}\right)}^2\right)$$

(23)

where P_t,ij,k, Q_t,ij,k, and i_t,ij,k are the active power flow, the reactive power flow, and the current magnitude square of branch ij in the kth iteration, respectively; u_t,i,k is the voltage magnitude square of node i in the kth iteration.

The maximum SOCP relaxation deviation of the natural gas system is defined as follows:

$$r{\mathrm{Gap}}_{\mathrm{NGN}}(x)={\displaystyle \sum _{ij\in {\mathsf{\Omega }}_n}({\overline{v}}_{l,t,k}-{\underset{¯}{v}}_{l,t,k})-{(F_{l,ij}^{\prime }{q}_{l,ij,t,k}^{})}^2}$$

(24)

where Ω_n is the set of all pipelines in the NGN; ${\overline{v}}_{l,t,k}$ is the larger of $p_{l.i.t}$ and $p_{l,j,t}$ in the kth iteration; ${\underset{¯}{v}}_{l,t,k}$ is the smaller of $p_{l.i.t}$ and $p_{l,j,t}$ in the kth iteration; q_l,ij,t,k is the gas pipeline l volume flow in the kth iteration.

In order to ensure the accuracy of SOCP relaxation, increasingly tight linear cuts of distribution system and natural gas system can be expressed in (25) and (26).

$${\displaystyle \sum _{ij\in {\mathsf{\Omega }}_b}{r}_{ij}}{i}_{t,ij}\le {\displaystyle \sum _{ij\in {\mathsf{\Omega }}_b}{r}_{ij}}\frac{{\left(P_{t,ij,k-1}\right)}^2+{\left(Q_{t,ij,k-1}\right)}^2}{u_{t,i,k-1}}$$

(25)

$${\displaystyle \sum _{ij\in {\mathsf{\Omega }}_n}({\overline{v}}_{l,t}-{\underset{¯}{v}}_{l,t})}\le {\displaystyle \sum _{ij\in {\mathsf{\Omega }}_n}{(F_{l,ij}^{\prime }{q}_{l,ij,t,k-1}^{})}^2}$$

(26)

By now, through SOCP relaxation and linearization, the maximum hosting capacity model of solar energy in IEDS with PV and SC is reformulated as the MISOCP model.

The enhanced SOCP-based method for evaluating the maximum hosting capacity of solar energy in this paper is shown in Figure 1. The specific operation process includes nine steps:

①

Basic data inputting;

②

Initialization parameters setting;

③

Check whether k is fewer than or equal to k_max. If so, continue to step 4. Otherwise, terminate the process;

④

Model constructing;

⑤

Model converting;

⑥

Model solving;

⑦

Check whether $r{\mathrm{Gap}}_{\mathrm{DN}}\le {\epsilon }_1 \&\& r{\mathrm{Gap}}_{\mathrm{GN}}\le {\epsilon }_2$. If so, move to step 9. Otherwise, continue to step 8;

⑧

Cutting planes adding and move to step 3;

⑨

Results outputting.

4. Case Study

4.1. Case Introduction

The structure diagram of IEDS is shown in Figure 2, including the modified 11-node natural gas network [23], the modified IEEE 33-node distribution network [33], the modified 32-node heating network [35], and three energy stations. The No. 1 energy station is coupled with No. 10 (E10) of the distribution network, No. 2 (G2) of the natural gas network, and No. 1 (H1) of the heating network. The No. 2 energy station is coupled with No. 24 node (E24) of the distribution network, No. 6 node (G6) of the natural gas network, and No. 31 node (H31) of the heating network. The No. 3 energy station is coupled with No. 31 node (E31) of the distribution network, No. 7 node (G7) of the natural gas network, and No. 32 node (H32) of the heating network. The detailed parameters of NGN, DN, and HN are provided in [23,33,35,36].

The structure diagram of the energy station is shown in Figure 3, including PV, SC, combined heat and power (CHP) units, GB, EB, and P2G.

The assumptions on operation conditions of the integrated energy distribution system [2,33,35,37] are as follows: The maximum installation area of PV and SC at each energy station is 15,000 m². The predicted value of irradiance intensity is taken as 700 W/m² [37]. The maximum allowable branch current is 250 A. The allowable range of the DN voltage is 0.9–1.1 p.u. and the allowable pressure of the NGN is 35–75 mbar. We assume that the maximum acceptable water velocity in pipelines is 2 m/s, the temperature difference between water at the inlet and outlet of the pipe is 25 °C, and the heat loss ratio per unit length is 0.15/km. The predefined precision with regard to the SOCP relaxation deviation of the DN and the NGN are set to 1 × 10⁻⁶ and 1 × 10⁻², respectively. The G_CV of natural gas is 41.04 MJ/m³.

The prediction curves of typical daily solar irradiance and load data are shown in Figure 4.

Table 2. Parameters of the energy converters.

Converter	Capacity/kW	Efficiency
CHP	$G_{i,t}^{CHP}=300/0.3$	${\eta }_{\mathrm{CHP}-\mathrm{P}}=0.3\hspace{0.33em}{\eta }_{\mathrm{CHP}-\mathrm{Q}}=0\hspace{0.33em}{\eta }_{\mathrm{CHP}-\mathrm{H}}=0.39$
EB	$P_{i,t}^{EB}=200$	${\eta }_{\mathrm{EB}}=0.95$
GB	$G_{i,t}^{GB}=300/0.85$	${\eta }_{\mathrm{GB}}=0.85$
P2G	$P_{i,t}^{P2G}=200$	${\eta }_{\mathrm{P}2\mathrm{G}}=0.7$
PV	${\eta }_{i,t}^{\mathrm{SE}}S{E}_i^{PV}=1873.5$	${\eta }_{\mathrm{PV}-\mathrm{P}}=0.175\hspace{0.33em}{\eta }_{\mathrm{PV}-\mathrm{Q}}=0$
SC	${\eta }_{i,t}^{\mathrm{SE}}S{E}_i^{SC}=5250$	${\eta }_{\mathrm{SC}}=0.5$

lists the parameters of the energy converters of the energy station.

Table 3. Nodal energy demand and source pressure.

Node Number	Energy Demand (kJ/s)	Pressure (mbar)
1 (Source Node)	0	75
2	1250	/
3	1100	/
4	1000	/
5	1300	/
6	900	/
7	250	/
8	1175	/
9	275	/
10	237.5	/
11	175	/

and

Table 4. Network pipe data.

Branch	From–To	Pipe Length (m)	Pipe Diameter (mm)
1	1–2	50	160
2	2–3	500	160
3	2–4	500	110
4	2–5	500	110
5	3–6	600	110
6	3–7	600	110
7	3–8	500	110
8	7–9	200	80
9	9–10	200	80
10	10–11	200	80

list the parameters of the NGN.

The energy station mode of IEDS in this paper is formulated as follows:

$$\mathit{C}=\left[\begin{array}{cccccccccc}{\eta }_{\mathrm{CHP}-\mathrm{P}}& -1& 0& -1& {\eta }_{\mathrm{PV}-\mathrm{P}}& 0& 1& 0& 0& 0\\ {\eta }_{\mathrm{CHP}-\mathrm{Q}}& 0& 0& 0& {\eta }_{\mathrm{PV}-\mathrm{Q}}& 0& 0& 1& 0& 0\\ -1& 0& -1& {\eta }_{\mathrm{P}2\mathrm{G}}& 0& 0& 0& 0& 1& 0\\ {\eta }_{\mathrm{CHP}-\mathrm{H}}& {\eta }_{\mathrm{EB}}& {\eta }_{\mathrm{GB}}& 0& 0& {\eta }_{\mathrm{SC}}& 0& 0& 0& 1\end{array}\right]$$

(27)

$$\mathit{I}={\left[G_{i,t}^{CHP}\hspace{0.33em}{P}_{i,t}^{EB}\hspace{0.33em}{G}_{i,t}^{GB}\hspace{0.33em}{P}_{i,t}^{P2G}\hspace{0.33em}{\eta }_{i,t}^{\mathrm{SE}}S{E}_i^{PV}\hspace{0.33em}{\eta }_{i,t}^{\mathrm{SE}}S{E}_i^{SC}\hspace{0.33em}{P}_{i,t}^{line}\hspace{0.33em}{Q}_{i,t}^{line}\hspace{0.33em}{G}_{i,t}^{line}\hspace{0.33em}{H}_{i,t}^{line}\right]}^T$$

(28)

$$\mathit{L}=\left[P_{i,t}^{LOAD}\hspace{0.33em}{Q}_{i,t}^{LOAD}\hspace{0.33em}{G}_{i,t}^{LOAD}\hspace{0.33em}{H}_{i,t}^{LOAD}\right]$$

(29)

The proposed method was implemented in the YALMIP [38] optimization toolbox (version 20200930) using MATLAB R2020a and solved by IBM ILOG CPLEX 12.6. The numerical experiments were performed on a computer with an Intel CORE CPU i7-8750H processor running at 2.20 GHz and 16 GB of RAM.

4.2. The Single-Period Case (Case 1)

In the single-period case (named Case 1), the basic data are shown in Section 4.1, and power cannot be sent back to the upstream power grid. The quantity ϕ₁ is set to 1, the same as the value of the quantity ϕ₂. Based on the above data, five scenarios are set as follows:

• Scenario I: Only PV are considered based on the DN.
• Scenario II: PV, CHP, GB, and EB are considered based on the IEDS.
• Scenario III: Based on Scenario II, P2G is considered.
• Scenario IV: Based on Scenario II, SC is considered.
• Scenario V: Based on Scenario II, SC and P2G are considered.

The accuracy and efficiency of the proposed method are verified as follows.

• Step 1: Input the basic network data and parameters of the devices.
• Step 2: The predefined precision about the SOCP relaxation deviation of the DN and the NGN are set to 1 × 10⁻⁶ and 1 × 10⁻², respectively.
• Step 3: The maximum number of iteration steps is set to 30. Initialize the iteration step k = 1. Check whether k is fewer than or equal to 30. If so, proceed to Step 4. Otherwise, the process terminates.
• Step 4: Build the optimization model for the maximum hosting capacity evaluation of solar energy.
• Step 5: Convert this model into an MISOCP model through SOCP relaxation and linearization.
• Step 6: Solve the MISOCP model to obtain the maximum relaxation deviation of the DN and the NGN.
• Step 7: Check whether $r{\mathrm{Gap}}_{\mathrm{DN}}\le 1\times {10}^{-6} \&\& r{\mathrm{Gap}}_{\mathrm{NGN}}\le 1\times {10}^{-2}$. If so, move to Step 9. Otherwise, continue to Step 8.
• Step 8: Update k = k + 1. Add the cutting plane constraint (14) and (15), and return to Step 3.
• Step 9: Output the optimization results and end the solving process.

The hosting capacities of solar energy in five scenarios are listed in

Table 5. The hosting capacity of solar energy in five scenarios.

Scenario	PV Capacity/MW	SC Capacity/MW
I	3.795	0.000
II	3.739	0.000
III	4.341	0.000
IV	4.379	1.810
V	4.854	1.882

.

Based on the comparison of Scenario II and Scenario III in Table 5, the hosting capacity of PV is increased by 16.10% relative to that of Scenario II because of the utilization of P2G.

Based on the comparison of Scenario II and Scenario IV in Table 5, the hosting capacity of PV is increased by 17.12% relative to that of Scenario II, and the hosting capacity of SC is increased from 0 MW to 1.810 MW because of the utilization of SC.

Based on the comparison of Scenario IV and Scenario V in Table 5, the hosting capacities of PV and SC are also increased because of the utilization of P2G.

Figure 5 shows the optimal dispatch results of electrical power in five scenarios. Considering that electrical power cannot be sent back to the upstream power grid, the sum of the electrical output power of the PV and CHP unit is exactly equal to the sum of the electrical input power of EB, the electrical input power of P2G, the distribution system active power load, and the distribution system active power losses.

Figure 6 shows the optimal dispatch results of thermal power in four scenarios. Considering that thermal power can be produced only by the SC, CHP, EB, and GB, the sum of the thermal output power of the SC, CHP, EB, and GB is exactly equal to the sum of the heating system thermal power load and thermal power losses.

As shown in Figure 5 and Figure 6, based on the comparison of Scenario III and Scenario V, the utilization of SC can reduce the thermal output power of CHP and GB in Scenario V, and the electrical output power of CHP in Scenario V is also reduced, which provides a “buffer” to accommodate more solar energy. Based on the comparison of Scenario IV and Scenario V, the utilization of P2G can increase the electrical input power in Scenario V, and the thermal output power of EB in Scenario V is also reduced, which also provides a “buffer” to accommodate more solar energy. Therefore, the optimal utilization of multiple energy in Case 1 can effectively increase the hosting capacity of PV and SC.

The impact of distribution system and heating system on nodal pressure across the gas system for Case 1 is shown in Figure 7. Due to the distributed injection of P2G, the nodal pressure of Scenario III increases compared to that of Scenario II. Due to the utilization of SC and the decrease of the natural gas flow demand of the CHP and GB, the nodal pressure of Scenario IV increases compared to that of Scenario II. Therefore, the utilization of P2G and SC in Case 1 can efficiently support the pressure management of the network.

The objective function consists of two main parts: the output power of the PV and the output power of the SC, whose weight could influence the results of the proposed model. The optimization results of scenarios IV and V considering the influence of the quantity ϕ₁ and the quantity ϕ₂ in Case 1 are listed in

Table 6. Optimization results of Scenario IV (Case 1).

Capacity/MW	${\mathit{\varphi }}_1 = 0.2,\hspace{0.33em}{\mathit{\varphi }}_2 = 0.8$	${\mathit{\varphi }}_1 = 0.4,\hspace{0.33em}{\mathit{\varphi }}_2 = 0.6$	${\mathit{\varphi }}_1 = 0.5,\hspace{0.33em}{\mathit{\varphi }}_2 = 0.5$	${\mathit{\varphi }}_1 = 0.6,\hspace{0.33em}{\mathit{\varphi }}_2 = 0.4$	${\mathit{\varphi }}_1 = 0.8,\hspace{0.33em}{\mathit{\varphi }}_2 = 0.2$
Total PV	3.796	3.796	4.398	4.398	4.398
Total SC	2.313	2.313	1.743	1.743	1.743
Total PV + SC	6.109	6.109	6.141	6.141	6.141

and

Table 7. Optimization results of Scenario V (Case 1).

Capacity/MW	${\mathit{\varphi }}_1 = 0.2,\hspace{0.33em}{\mathit{\varphi }}_2 = 0.8$	${\mathit{\varphi }}_1 = 0.4,\hspace{0.33em}{\mathit{\varphi }}_2 = 0.6$	${\mathit{\varphi }}_1 = 0.5,\hspace{0.33em}{\mathit{\varphi }}_2 = 0.5$	${\mathit{\varphi }}_1 = 0.6,\hspace{0.33em}{\mathit{\varphi }}_2 = 0.4$	${\mathit{\varphi }}_1 = 0.8,\hspace{0.33em}{\mathit{\varphi }}_2 = 0.2$
Total PV	4.399	4.399	4.854	4.854	5.000
Total SC	2.313	2.313	1.881	1.881	1.465
Total PV + SC	6.712	6.712	6.735	6.735	6.465

, respectively. Because the total alternative installation area of PV and SC in the system is limited, the results of the proposed model are affected by the quantity ϕ₁ and the quantity ϕ₂. From Table 6 and Table 7, we can see that as the quantity ϕ₁ becomes larger, the hosting capacity of PV is increased and the hosting capacity of SC is decreased. Because the quantity ϕ₁ becomes larger, more space can be provided for the installation area of PV to maximize the objective function.

4.3. The 24 h Period Case (Case 2)

In Case 2, the basic data are shown in Section 4.1, and power cannot be sent back to the upstream power grid. The quantity ϕ₁ is set to 1 and the quantity ϕ₂ is set to 1 in Case 2. Only PV considered based on the DN is set as Scenario VI. PV, SC, CHP, GB, EB, and P2G considered based on the IEDS is set as Scenario VII. Because the utilization of P2G and SC can support the pressure management of the NGN, the allowable pressure of the NGN is set to 50–75 mbar.

The hosting capacities of solar energy in Scenario VI and VII are listed in

Table 8. The hosting capacity of solar energy in Scenario VI and VII.

Capacity/MW	Scenario VI	Scenario VII
Total PV	4.684	4.911
Total SC	0.000	1.718

.

From Table 8, we can see that in Scenario VII, the hosting capacity of PV is increased by 4.85% relative to the result of Scenario VI, and the hosting capacity of SC is increased from 0 MW to 1.718 MW. Compared with Scenario VI, Scenario VII realizes the integration and coordination of the distribution system, heating system, and gas system, which can effectively accommodate more solar energy.

The optimal dispatch results of electrical power, thermal power, and gas power are shown in Figure 8.

In Figure 8, during the period of 12:00 to 15:00, the source of the electrical power is mainly the PV, and the source of the thermal power is mainly the EB and the SC. At night, the output of CHP is high, and the natural gas flow demand is increased. In the period of high irradiance intensity, the CHP and the GB maintain a state of zero output and the electrical input power of the EB and the P2G is high, which provides extra space to accommodate more solar energy.

The minimum pressure of Scenario VII in each time period is shown in Figure 9, where the minimum pressure of Scenario VII is greater than or equal to 50 mbar through the utilization of P2G and SC and the IEDS operation optimization. During the period of 13:00 to 15:00, the minimum pressure can be increased because the lower gas load and the thermal output power of SC reduced the natural gas flow demand of the CHP, and, thus, reduced the natural gas flow from main supply source.

As shown in Figure 10, the distribution system active power loss over a day of scenarios VI and VII is reduced from 1106.21 kWh to 696.95 kWh, which is a considerable improvement of economic efficiency.

The minimum distribution system voltages in each time period are shown in Figure 11, where a flat voltage profile was attained through the IEDS operation optimization under Scenario VII. The minimum distribution system voltages in Scenario VII are greater than or equal to 0.95 p.u.

Compared with Scenario VI, Scenario VII realizes optimal utilization of multiple energy, which can reduce distribution system active power losses and the minimum distribution system voltage deviation from the nominal value (1 p.u.).

4.4. Algorithm Validation

To verify the exactness of SOCP relaxation, the SOCP relaxation deviations of different scenarios in each iteration (Case 1) are shown in Figure 12, and the relaxation deviations of distribution system and natural gas system in each time period are shown in Figure 13. All relaxation deviation values of distribution system and natural gas system in Case 1 and Case 2 are smaller than 1 × 10⁻⁶ and 1 × 10⁻², respectively.

The relaxation deviation results in Figure 12 and Figure 13 show that the enhanced SOCP-based approach can calculate the maximum hosting capacity of solar energy in IEDS with acceptable accuracy.

BONMIN is an experimental open-source C++ code for solving general (mixed-integer nonlinear programming) MINLP problems [39]. For that reason, to further verify the accuracy and efficiency of the proposed method, the BONMIN solver is much more suitable for the initial maximum hosting capacity evaluation problem which can be solved using MINLP.

Table 9. Optimization results of the proposed method and BONMIN for Case 1.

Scenario	The Objective Function (MWh)		Time (s)
	BONMIN	Proposed Method	BONMIN	Proposed Method
I	3.715	3.715	0.218	0.244
II	3.490	3.500	6.138	0.581
III	3.940	4.105	5.588	0.494
IV	5.909	5.908	15.223	2.375
V	6.502	6.502	6.718	0.165

shows that compared with the BONMIN package, the proposed method can obtain an accurate solution and can greatly improve the computation efficiency. Because of the convex relaxation for the original problem, the proposed method has the advantages that the computation time will not increase greatly with the increase of problem scale caused by a larger system and the computation performance is better. However, the increase of problem scale tends to cause “the curse of dimensionality” when using the BONMIN package, which may greatly increase the computation time. The maximum hosting capacity evaluation method of solar energy in an integrated energy distribution system based on enhanced SOCP can output a satisfactory solution and reduce the computation time.

5. Conclusions

This paper proposes an evaluation method on maximum hosting capacity of solar energy in an integrated energy distribution system based on enhanced SOCP. When comparisons are made between different scenarios, the following conclusions are drawn:

(1)

The optimization results show that the maximum hosting capacity of solar energy is improved significantly by realizing the coordination of integrated multi-energy system and the optimal utilization of electricity, heat, and gas energy. With the utilization of gas energy (P2G, etc.), the hosting capacity of PV increases from 3.795 MW in Scenario II to 4.341 MW. With the utilization of gas energy (SC, etc.), the hosting capacity of PV increases from 3.795 MW in Scenario II to 4.379 MW and the hosting capacity of SC is increased from 0 MW to 1.810 MW.

(2)

Meanwhile, the distribution system power losses and the voltage fluctuations are effectively decreased with the optimal utilization of multiple energy. The distribution system active power loss over a day reduced from 1106.21 kWh in Scenario VI to 696.95 kWh, and a flat voltage profile was attained through the IEDS operation optimization.

(3)

By applying SOCP relaxation, linearization, and adding increasingly tight linear cuts of distribution system and natural gas system to the SOCP relaxation, the proposed model can be solved accurately and efficiently.

Several notable issues are worth further research. In future work, the influence of multiple energy storage (electrical energy storage system, thermal energy storage system, cooling energy storage system) on an integrated energy distribution system should be considered to improve the maximum hosting capacity of solar energy.