Probabilistic energy flow for multi-carrier energy systems

Renewable and Sustainable Energy Reviews 94 (2018) 989–997

Contents lists available at ScienceDirect

Renewable and Sustainable Energy Reviews journal homepage: www.elsevier.com/locate/rser

Probabilistic energy flow for multi-carrier energy systems Hosein Khorsand, Ali Reza Seifi

⁎

T

School of Electrical and Computer Engineering, Shiraz University, Shiraz, Iran

A R T I C LE I N FO

A B S T R A C T

Keywords: Multi-carrier energy system Energy hub Probabilistic Monte Carlo simulation Point estimate method

This paper investigates energy flow of Multi-Carrier Energy Systems using Point Estimate Method and Monte Carlo simulation considering uncertainties, which may happen in the Electrical, Natural Gas, and District Heating Networks all together. These uncertainties could occur from the probabilistic behavior of loads or unforeseen faults. An innovative probabilistic energy flow of Multi-Carrier Energy Systems based on Point Estimate Method has been presented in this paper. Results for two different case studies are investigated and compared against those achieved from the Monte Carlo simulation. The results prove that the presented point estimate schemes have precise results, smaller computational burden, and time, comparing than Monte Carlo simulation method.

1. Introduction Nowadays, by increasing the importance of utilization of cogeneration plants, which makes a strong coupling between the network industries such as Electrical, Natural Gas, and District Heating Networks (DHNs), the uncertainties in energy system, such as load variation, and faults of lines have been in great attention [1]. Therefore, the use of techniques able to account for uncertainty are requested to control and minimize the risks related to design and operation [2]. In this way, probabilistic methods are more appropriate for this purpose. Several techniques such as Monte Carlo Simulation (MCS) method, analytical methods, and approximate methods have been studied to overcome problems under uncertainty [3,4]. MCS generates random values for uncertain input variables, and these values are applied to simulate a deterministic problem in each simulation [5]. This technique has been widely utilized in electrical network analysis to model uncertainty. The main drawback of the MCS is the excessive number of computation required to achieve convergence [6]. Thus to decrease the computational burden in solving Probabilistic Load-Flow (PLF) problem, several analytical approaches and approximate methods were suggested to simulate the load-flow solution distributions [7]. Analytical methods are computationally adequate, but they need some mathematical assumptions in order to simplify the problem. These techniques usually simplify the traditional PLF formulations and calculate the convolution of different probabilistic variables through their linear relationships [8]. Besides, the convolution techniques achieve a mathematical description of the behavior of output random variables

⁎

Corresponding author. E-mail addresses: [email protected] (H. Khorsand), seifi@shirazu.ac.ir (A.R. Seifi).

https://doi.org/10.1016/j.rser.2018.07.008 Received 19 October 2017; Received in revised form 27 January 2018; Accepted 1 July 2018

1364-0321/ © 2018 Elsevier Ltd. All rights reserved.

(RVs). Approximate methods such as Point Estimate Method (PEM) use deterministic process in order to simulate probabilistic problems by applying only some statistical moments of probability functions such as mean, variance, skewness, and kurtosis coefficients [9]. Hence, PEM overcomes the difficulties associated without perfect knowledge of the probability functions of stochastic variables [10]. The PEM in this approach is somehow better than MCS method because it has the less computational burden and consumed time since a smaller volume of data is requested [11,12]. This paper presented the point estimate method to solve the probabilistic energy flow of Multi-carrier Energy Systems (MES). Various energy carriers such as electrical, natural gas, and DHN could be considered in an energy network. Electrical grids are the most favorite energy carriers employed in many countries [13]. Usually, an electrical grid is made up of electrical generation plants, which generate electric power by consuming suitable fuel such as gas, oil, and coal. Electrical loads are supplied by the generated electric power through transmission lines [14,15]. Natural gas is getting one of the most employed sources of energy in the world since the significance of its offer, market potential, adjustability, and smaller environmental issues in comparison to other effective fossil sources [16]. Natural gas supplies many commercial and residential consumers over the world by means of a vast and complex grid. A typical natural gas network is made up of one or more gas resources, numerous loads, pipelines, compressors, and other elements such as valves or regulators [17].

Renewable and Sustainable Energy Reviews 94 (2018) 989–997

H. Khorsand, A.R. Seifi

Nomenclature

B TR H C D L τ η N a, b, c, d,

Subscript

E G H S r g

electric natural gas heat supply pipeline return pipeline ground

susceptance of transmission lines or shunt capacitances (s) tap ratio of tap-transformers compression ratio of compressors constant of natural gas pipelines diameter of natural gas pipelines (m) total length of gas/heat pipelines (m) pump head of the district heating network(m) efficiency of units total number of units e heat rate coefficients of generators/characteristic coefficients of CHPs and boilers

Superscript Variables/parameters of the Point Estimate Method hub gen chp boil gs shc dem comp pump fuel line bus min max

energy hub electric generator combined heat and power plant boiler natural gas source shunt capacitance demand compressor circulating pump fuel transmission line/pipeline electrical/gas/heat bus minimum limit of a variable maximum limit of a variable

Pl, k ωl, k μ σ ζ k Z (l, k ) Mj (pl ) fpl λ l, j m ε RV εMCS Constants

Variables/parameters of the MES model

P Q S f ṁ V θ π T G

Kth value of random variable Pl weighting factor mean value standard deviation standard location total number of random variables output random variables jth central moment of the random variable Pl PDF of random variable Pl coefficient of skewness number of input random variables error index of the random variable relative error of Monte Carlo method

ε z δ c ρ ψ, ξ

electric/gas/heat power (MW) reactive power (MVar) apparent power (MVA) natural gas flow (m3/day) mass flow rate of water (kg/s) voltage magnitude (p.u.) voltage angle (o) natural gas pressure (bar) temperature (°C) conductance of transmission lines (s)

g

The International Energy Agency (IEA) declared that due to all energy consumption, heating, cooling estimated for about 46% of the entire global energy consumes in 2012, and about half of end use of energy in Europe is heat [18]. For that reason, DHNs have considerable ability in satisfying consumer's necessity. In addition, by using such a network, the environmental issues such as carbon dioxide emissions have been decreased sufficiently [19]. The effort of a DHN is to yield the consumers demand heat for its loads by applying sufficient heat production plants and a network of pipelines to transfer this produced heat power to the consumers. Hence, Shabanpour utilized a proper set of state-variables to optimize energy flow of electricity, natural gas, and DHN simultaneously [20]. Combined Heat and Power (CHP) plant is playing an indispensable role in producing the district heat by consuming natural gas. Thus, the entire fuel consumption can be diminished sufficiently. Shabanpour proposed the modified teaching-learning based optimization method, which analyzed energy flow optimization in MES. The presented method, unlike conventional techniques, did not use any extra variable such as dispatch factors or dummy variables [21]. Geidl presented a general optimization technique, which studied energy flow optimization in multi-carrier energy systems [22]. He considered the possibility of conversion between electricity, natural gas, and DHN.

absolute rugosity of natural gas pipelines (0.05 mm) natural gas compressibility factor (z = 0.8) density of natural gas relative to air (δ = 0.6106 ) specific heat capacity of water (cp = 4182 J/kg K) heat transition coefficient (U = 0.455 W/m K) constants of compressors (ψ = 0.167 for turbo-compressor, ψ = 0.157 for a moto-compressor, and ξ = 0.236 for both types) standard gravity constant (g = 9.81 m/s2)

Fig. 1. An example for an energy hub.

Shabanpour presented the new approach based on teachinglearning algorithm, which analyzed multi-objective optimal energy flow problem taking into account the total fuel cost and total emission of the MES [23]. The proposed algorithm applies a self-adapting wavelet mutation technique. In addition, a fuzzy clustering approach is presented to reduce repository size besides a smart population selection for the next iteration. Mohammadi presented a new point estimate scheme to solve the probabilistic power flow of electrical networks while he considered load variation and faults of lines [24]. Morales presented point estimate method to solve the probabilistic 990

Renewable and Sustainable Energy Reviews 94 (2018) 989–997

H. Khorsand, A.R. Seifi

Table 1 Configuration of the MES-Case1.

From stotage Node Compressor

From stotage Zeebrugge

1

5

Dudzele 4 Brugge 3 Zomergem

Heat node

Generator1 Generator2 Compressor1 Compressor2 Slack boiler and pump CHP, boiler, pump CHP, boiler, pump CHP, boiler, pump CHP, boiler, pump

1 2 4 5 12 6 13 14 8

12 19 9 18 3 15 7 6 10

– – – – 1 4 9 10 13

2. System modeling 2.1. Energy hub concept

6 7 Antwerpen

Norvegian gas

Gent Berneau Anderlues

Péronnes Mons 14

15

Warnand-Dreye

13

9

An energy hub can be illustrated as a unit that interfaces with consumers, producers, storage devices, electrical transformers, power electronic devices, gas compressors, heat exchangers, or boilers, and other equipment. An example of an energy hub is shown in Fig. 1. It contains a CHP unit, which consumes natural gas and generates electricity and heat. In addition, a boiler is applied to produce heat from natural gas. An electric pump circulates water for heating purposes. Either unit of an energy hub may have constant efficiencies, which make a linear relationship between the input and output of the energy hub, or it can indicate a dependency that yields a nonlinear relationship. Any energy hub is linked to various sub-networks such as electricity, natural gas, and DHN. The energy flow of MES used in this paper can be derived from (1−9) which is fully described in [30] by the authors. From a system point of view, the rest of electricity and heat demand of an energy hub is provided by the electrical and DHN, respectively. Therefore, if the produced electricity or heat power within the hub is larger than the required power for its loads the surplus power can be pass to the respective sub-network. Therefore, the energy flow of the energy hub can be computed as follows:

8 Voeren

11 10 Liège 12 17 Wanze

From stotage

Namur

16 Blaregnies

Sinsin

18

To France

Arlon

19 20 To Luxemburg

Fig. 3. Schematic diagram of the natural gas sub-network-Case1&2.

30km

10

25 km

m

km

20

4

25 km

8

10k

20km

m

5

m

6

5k

7

12

13 10 km

3

11

5km

30k

20km

5k

9

m

1

5km

Gas node

Loenhout

2

2

Electrical bus

Chun-Lien Su Indicates how to estimate the corresponding uncertainty in the probabilistic power flow solution. Likewise, the presented approach may be employed to solve the deterministic power flow problems [26]. Many studies have been done on the probabilistic power flow in particular for electricity networks [27–29]. However, to our knowledge, no general probabilistic energy flow technique regarding the operation of Multi-carrier Energy Systems has been published. Therefore, in this paper, the point estimate method is applied to solve the probabilistic energy flow of Multi-carrier Energy Systems.

Fig. 2. Schematic diagram of the electrical sub-network-Case1.

Algerian gas

Unit

PEhub + PEchp − PEpump = PEdem

(1a)

PGhub = PGchp + PGboil

(1b)

PHhub + PHchp + PHboil = PHdem

(1c)

2.2. Equality constraints

14

Equality constraints indicate the energy flow for each sub-network and make the energy flow feasibility in these systems. The equality constraints for the electrical sub-network can be calculated as follows:

Fig. 4. Schematic diagram of the district heating sub-network-Case1&2.

power flow of electrical networks [25]. He utilized Hong's point estimate schemes to prepare a proper solution for the tradeoff between the accuracy of the results and the efficiency of the computational burden for large-scale power system problems.

bus

dem PEgen . i = PE . i +

NE

∑ j=1

991

line Vi Vj (GEline . ij cos θij + BE . ij sin θij )

(2a)

Renewable and Sustainable Energy Reviews 94 (2018) 989–997

H. Khorsand, A.R. Seifi

Table 2 Nodal probabilistic data for natural gas and DHN.

Gas Network ( fGdem ,i )

dem Heat Network (PH ,i )

Node number

Distribution Type

3

Normal: Mean value = 4,standard deviation = 0.3

7 10 12 15 3

Normal: Mean value = 5, standard deviation = 0.1 Uniform: Upper band = 6.3,Lower band = 6 Normal: Mean value = 2, standard deviation = 1 Normal: Mean value = 6, standard deviation = 0.2 Normal: Mean value = 10, standard deviation = 0.3

4 7 9 13 14

Uniform: Upper band = 11,Lower band = 10 Normal: Mean value = 15, standard deviation = 0.1 Uniform: Upper band = 11,Lower band = 10 Uniform: Upper band = 26,Lower band = 25 Uniform: Upper band = 16,Lower band = 15

Table 3 Results of buses voltage amplitude (pu)-Case1. Electrical bus

Table 6 Results of mass flow rate of water (kg/s)-Case1. Heat line i-j

Vj MCS

1 3 4 7 9 10 14

2 m PEM

2 m+ 1 PEM

MCS

μ

σ

μ

σ

μ

σ

1.0600 0.9641 0.9823 0.9964 0.9914 0.9903 1.0000

5.4184e-14 3.2784e-04 5.4592e-04 3.0447e-04 4.6355e-04 7.2564e-04 8.7069e-16

1.0600 0.9471 0.9564 0.9810 0.9741 0.9643 1.0000

0 0.0163 0.0247 0.0148 0.0165 0.0249 0

1.0600 0.9595 0.9753 0.9922 0.9867 0.9834 1.0000

0 0.0005 0.0008 0.0004 0.0006 0.0009 0

1–2 2–3 5–6 6–8 5–9 10–11 13–14

Table 4 Results of buses voltage angle (o )-Case1. Electrical bus

1 3 4 7 9 10 14

Heat node 2 m PEM

2 m+ 1 PEM

σ

μ

σ

μ

σ

0 − 11.7904 − 9.5762 − 12.8830 − 14.6356 − 14.8647 − 16.5542

0 0.0612 0.0803 0.1380 0.1685 0.1836 0.2536

0 − 14.7575 − 13.1748 − 18.9575 − 22.0110 − 22.6627 − 26.9801

0 2.8427 3.4474 5.8190 7.0651 7.4701 9.9850

0 − 12.5737 − 10.5283 − 14.5032 − 16.6061 − 16.9508 − 19.3579

0 0.0822 0.1084 0.1871 0.2284 0.2486 0.3462

1 4 7 10 12 14

1 3 6 9 10 13 16 18

σ

μ

σ

μ

σ

41.5022 36.4668 86.3311 34.6050 − 36.5559 42.9495 21.4712

1.2219 1.0471 0.3473 0.0304 0.7162 0.7521 0.6158

87.8356 57.2234 94.9156 35.2078 − 34.3741 59.0654 23.1316

44.3118 19.8348 8.2085 0.5732 1.8022 15.3785 1.3405

56.4725 41.9533 88.7534 34.7881 − 35.0519 48.5893 21.9716

1.6418 1.4931 0.4853 0.0281 1.0533 1.0042 0.8529

Ts . j 2 m PEM

2 m+ 1 PEM

μ

σ

μ

σ

μ

σ

110 113.3907 109.3496 114.8234 110.8815 108.6770

0 0.1090 0.0611 0.0351 0.0212 0.0824

110 111.1989 108.2090 113.9305 111.1376 109.5576

0 2.0845 1.0839 0.8526 0.2425 0.8343

110 112.7448 109.0021 114.5241 110.9492 108.9622

0 0.0902 0.0569 0.0442 0.0225 0.1150

Table 8 The return temperature (oC ) at the end of pipelines-Case1. Heat node

πj MCS

2 m+ 1 PEM

μ

MCS

μ

Table 5 Results of natural gas pressures (bar) -Case1. Gas node

2 m PEM

Table 7 The supply temperature (oC ) at the end of pipelines-Case1.

θj MCS

line ṁ H . ij

2 m PEM

Tr . j MCS

2 m+ 1 PEM

μ

σ

μ

σ

μ

σ

56 55.8300 52.5287 49.4288 57.4980 53.3097 50.2506 45.8981

0 0.0212 0.2488 0.5561 0.6885 0.5242 0.5297 0.8609

56 55.1130 44.3005 34.3758 38.8483 38.2066 34.3419 22.6260

0 0.6867 7.8806 14.4033 17.8451 14.4643 15.2392 22.2697

56 55.7058 51.2058 47.4322 55.0231 51.2791 48.0161 42.8517

0 0.0231 0.2489 0.5225 0.6470 0.4851 0.5000 0.8073

1 4 7 10 12 14

2 m PEM

2 m+ 1 PEM

μ

σ

μ

σ

μ

σ

37.1650 39.2925 40 38.4303 40 40

0.0491 0.0082 0 0.0140 0 0

38.9815 39.4729 40 38.7324 40 40

1.7381 0.1721 0 0.2884 0 0

37.6672 39.3526 40 38.5344 40 40

0.0456 0.0100 0 0.0183 0 0

Table 9 The relative error (%) of different iteration of MCS for Case 1 test system. bus

QEgen .i

=

QEdem .i

NE

+

∑ j=1

line Vi Vj (GEline . ij sin θij − BE . ij cos θij )

(2b)

The Eqs. (2a) and (2b) convince the balance of generation value 992

NMCS Value

100

500

1000

3000

5000

7000

10,000

μ σ

0.1105 8.1154

0.0561 5.8194

0.0387 2.0379

0.0218 1.3203

0.0167 1.0120

0.0141 0.8219

0.0114 0.6901

Renewable and Sustainable Energy Reviews 94 (2018) 989–997

H. Khorsand, A.R. Seifi

Table 10 The average error indices (%) of PEM for Case 1 test system. Case 1

2 m PEM

Table 12 Results of buses voltage amplitude (pu)-Case2.

2 m+ 1 PEM

Electrical bus

Vj

εμRV

εσRV

εμRV

εσRV

MCS

Vj

1.4109

2256.6216

0.3787

64.3951

μ

σ

μ

σ

μ

σ (×10−3)

θj πj

46.0732

3796.1946

12.2764

32.8800

28.0376

2412.5113

6.2091

5.1266

line ṁ H . ij

77.4126

1831.4020

27.2476

33.5974

Ts . j Tr . j

1.2161

1832.5473

0.3468

12.2546

0.6650

1413.6991

0.1981

16.4167

1.0600 1.0160 1.0060 1.0485 1.0528 1.0493 1.0304 1.0373 1.0029 0.9933

1.5300e-13 5.89e-05 4.79e-05 2.22e-04 2.13e-04 2.4e-04 3.88e-04 2.54e-04 4.92e-04 4.98e-04

1.0600 1.0053 0.9998 1.0188 1.0232 1.0192 0.9827 1.0027 0.9485 0.9428

0 0.0104 0.0060 0.0290 0.0287 0.0294 0.0464 0.0337 0.0527 0.0491

1.0600 1.0136 1.0046 1.0419 1.0462 1.0426 1.0197 1.0296 0.9912 0.9821

0 0.0719 0.0614 0.2930 0.2809 0.3187 0.5228 0.3349 0.6805 0.7033

1 4 7 10 14 16 19 22 26 30

2 m PEM

2 m+ 1 PEM

Table 13 Results of buses voltage angle (o )-Case2. Electrical bus

θj MCS

1 4 7 10 14 16 19 22 26 30

2 m PEM

2 m+ 1 PEM

μ

σ

μ

σ

μ

σ

0 − 7.9825 − 11.2872 − 12.9745 − 12.3972 − 12.5737 − 13.8027 − 13.3563 − 14.0724 − 15.5732

0 0.0149 0.0162 0.0309 0.0324 0.0326 0.0393 0.0322 0.0385 0.0369

0 − 10.0451 − 13.5533 − 17.0655 − 16.6288 − 16.7332 − 18.3113 − 17.5615 − 18.6692 − 20.0674

0 2.0068 2.2049 3.9805 4.1175 4.0474 4.3870 4.0912 4.4398 4.3684

0 − 8.4384 − 11.7878 − 13.8779 − 13.3352 − 13.4947 − 14.8027 − 14.2830 − 15.0341 − 16.5562

0 0.0191 0.0208 0.0397 0.0419 0.0422 0.0510 0.0416 0.0519 0.0504

Table 14 Results of Natural Gas pressures (bar) -Case2. Gas node

πj MCS

1 3 6 9 10 13 16 18

Fig. 5. Schematic diagram of the electrical sub-network-Case2. Table 11 Configuration of the MES-Case2. Unit

Electrical bus

Gas node

Heat node

Generator1 Generator2 Compressor1 Compressor2 Slack boiler and pump CHP, boiler, pump CHP, boiler, pump CHP, boiler, pump CHP, boiler, pump

1 2 5 3 30 14 23 17 7

12 19 9 18 3 15 7 6 10

– – – – 1 4 9 10 13

2 m PEM

2 m+ 1 PEM

μ

σ

μ

σ

μ

σ

56 55.8306 52.5340 49.4706 57.5498 53.3229 50.2634 47.0784

0 0.0212 0.2457 0.5524 0.6838 0.5204 0.5262 0.8342

56 54.6718 39.6428 25.6397 28.0263 29.4654 25.2075 11.1504

0 1.1279 12.5487 23.1845 28.7231 23.2245 24.3941 34.9562

56 55.5465 49.4591 44.0926 50.8661 47.4461 43.9498 38.7547

0 0.0357 0.3807 0.8055 0.9995 0.8078 0.8369 1.2452

Table 15 Results of mass flow rate of water (kg/s) -Case2. Heat line i-j

line ṁ H . ij

MCS

1–2 2–3 5–6 6–8 5–9 10–11 13–14

with the load demand plus losses in this sub-network. It is worthwhile to note that the slack generator generates the difference of the total generated power from the total demand plus losses [31]. In the natural gas sub-network, the equality constraint can be

993

2 m PEM

2 m+ 1 PEM

μ

σ

μ

σ

μ

σ

41.5021 36.4744 86.3212 34.6049 − 36.5512 42.9515 21.4686

1.2207 1.0448 0.3491 0.0308 0.7106 0.7382 0.6148

114.3540 66.2047 100.0080 35.6251 − 32.2100 69.0029 23.8555

70.9045 28.9047 13.3053 0.9911 4.1172 25.3555 2.1596

60.6711 43.3463 89.3825 34.8267 − 34.3808 50.4761 22.0834

1.6638 1.4868 0.4911 0.0265 1.0767 1.0508 0.8881

Renewable and Sustainable Energy Reviews 94 (2018) 989–997

H. Khorsand, A.R. Seifi

Table 16 The supply temperature (oC ) at the end of pipelines-Case2. Heat node

Ts . j MCS

1 4 7 10 12 14

line

⎛ 3.7DG . ij ⎞ ⎤ −7 ⎡ CGline . ij = 3.0996 × 10 ⎢2 log ⎜ ⎟⎥ ε ⎝ ⎠⎦ ⎣

2 m PEM

σ

μ

σ

μ

σ

110 113.3917 109.3498 114.8234 110.8816 108.6769

0 0.1109 0.0620 0.0344 0.0209 0.0805

110 109.6719 107.4111 113.3850 111.3025 110.0716

0 3.6129 1.8830 1.4000 0.4085 1.3535

110 112.6110 108.9302 114.4335 110.9633 109.0449

0 0.0830 0.0538 0.0456 0.0217 0.1215

1 4 7 10 12 14

2 m PEM

2 m+ 1 PEM

μ

σ

μ

σ

μ

σ

37.1651 39.2924 40 38.4303 40 40

0.0489 0.0083 0 0.0137 0 0

40.0118 39.5949 40 38.9182 40 40

2.7715 0.2941 0 0.4749 0 0

37.7736 39.3674 40 38.5676 40 40

0.0421 0.0096 0 0.0189 0 0

100

500

1000

3000

5000

7000

10,000

μ σ

0.0765 8.0265

0.0245 3.3523

0.0191 2.4213

0.0135 1.3813

0.0107 1.0218

0.0093 0.8914

0.0084 0.7371

∑

PHline . ij

(4b)

line ⎛ ρLH . ij ⎞ Ts . j = (Ts . i − Tg )exp ⎜− line ⎟ + Tg ⎝ cṁ H . ij ⎠

(4c)

line ⎛ ρLH . ij ⎞ Tr . j = (Tr . i − Tg )exp ⎜− line ⎟ + Tg ⎝ cṁ H . ij ⎠

(4d)

chp chp chp chp . max PECHP PH . i + bichp Tschp ≤ PHchp . i = {a i . i + ci , d1. i PH . i .i chp . max chp chp . max ≤ PHchp ai PH . i + bichp Tschp − y1. i, d2. i PHchp ≤ PHchp .i . i + ci .i .i chp . max chp chp . min ≤ d1. i PHchp ai PH . i + bichp Tschp − y1. i − y2. i , PHchp .i . i + ci .i

Table 19 The average error indices (%) of PEM for Case 2 test system. 2 m PEM

2 m+ 1 PEM

εμRV

εσRV

εμRV

εσRV

Vj

2.4775

10,543.403

0.549

45.8295

θj πj

28.5357

12,216.6111

6.2804

28.7453

42.8251

4021.4251

15.9798

57.9458

line ṁ H . ij

126.532

3086.7213

35.5985

37.7505

Ts . j Tr . j

1.974

3113.2037

0.4219

16.6943

1.041

2332.9682

0.2466

17.7272

(5b)

. max y2. i = (d2. i PHchp − PHchp .i . i ) e2. i

(5c)

The equation (5) regards the part-load performance of the CHP. The electrical efficiency of CHPs is reduced at partial load which is fully described in [37,38] by the authors. Consequently, the gaseous input power of the ith CHP, boiler and compressor can be obtained by (6), (7) and (8) respectively where ηtchp is the total efficiency of the CHP and .i aiboil , biboil are constants related to the part-load performance of the ith boiler [39,40].

Table 20 CPU Time (in seconds). Monte Carlo

2m PEM

2m+ 1 PEM

Case1 Case2

354.95 651.56

54.20 89.53

40.22 86.70

j=1

(7)

ξ ⎞ ⎛ ⎛ πj ⎞ PEcomp = ψfGcomp − 1⎟ .i . ij ⎜ π i⎠ ⎝ ⎠ ⎝

(8)

(3a)

PEpump = .i 2 2 fGline = CGline . ij SG . ij (πi , πj ) SG . ij (πi, πj ). (πi − π j ) . ij

(6)

⎟

Where πj / πi = H comp is the compressor's compression ratio. The consumed electrical power of the ith circulating pump, which creates and holds pressure difference between supplies and returns pipelines, can be achieved by:

line

fGline . ij

ηtchp .i

⎜

NG

∑

chp PEchp . i + PH . i

boil boil . max PHboil PH . i . i + ai biboil

PGboil .i =

derived using the following equation:

fGgs. i = fGdem + .i

(5a)

. max y1. i = (d1. i PHchp − PHchp .i . i ) e1. i

PGchp .i =

CPU Time

(4a)

line PHline . ij = Cṁ H . ij (Ts . i − Tr . i )

chp . max ≤ PHchp . i ≤ d2. i PH . i

Case 2

(3c)

The Eq. (4a) expresses the balance of heat generation value with the load demand plus the outgoing power from it [34–36]. The heat power of the pipeline connected between node i and j can be expressed by (4b). (4c) and (4d) can derive the supply and return temperature at the end of this pipeline, respectively. Where the ground temperature at the pipeline's location is 10 oC. For a CHP unit, the produced electrical power can be achieved using the generated heat power of the unit and its temperature:

Table 18 The relative error (%) of different iteration of MCS for Case 2 test system. NMCS Value

NH

j=1

Tr . j MCS

zTδLGline . ij

line

boil dem PHchp . i + PH . i = PH . i +

Table 17 The return temperature (oC ) at the end of pipelines-Case2. Heat node

5

The Eq. (3a) defines the power balance between the input and output of each node of the gas network. The gas flow of each pipeline depends on the pressures at both ends of the pipeline and can be calculated by (3b) where SG . ij (πi , πj ) is a sign function which is + 1 if πi ≥ πj and −1 elsewhere [32]. Cij is the pipeline constant that can be obtained by (3c) [33].

2 m+ 1 PEM

μ

(DGline . ij )

(3b)

pump ṁ Hpump . i gτ pump ηi

(9)

More information about computing energy flows in a natural gas network and a DHN could be found in [20]. 994

Renewable and Sustainable Energy Reviews 94 (2018) 989–997

H. Khorsand, A.R. Seifi

3. Proposed methods

The standard location ζl, k and the weight ωl, k are achieved by solving the system of equations as below [25]:

3.1. Monte Carlo simulation method

K

∑ The complication of the analytical equations between the conversion system parameters, and especially, their nonlinearity, are the main reason for using the Monte Carlo simulation [28]. MSC is a direct technique for the probabilistic energy flow problem. It is to be noted that, unlike the deterministic approach, running Monte Carlo simulation does not need any additional calculations; it requires only the updating of the already generated system equations. In this approach, a wide variety of random load characteristics in Multi-Carrier Energy Systems is simulated until the resulting statistics agree with available field measurements. Therefore, variables are obtained for each randomly generated set of input variables as follow [29]. Let Z be a function of n independent RV as below:

Z (k ) = f (x1k , x 2k , ... ,x nk )

1 k

Var (Z ) =

1 k

λ l, j =

εMCS

NMCS ng

∑ ∑ i=1

j=2

Zi j

∫−∞

(pl − μpl ) j fpl dpl

(15b)

(

m

μj = E [Z j] ≅

)

K

∑ ∑

ωl, k (z (l, k )) j

(16)

l=1 k=1

where the Z (l, k ) values are employed to estimate the E (Z ) and E (Z j ) as below:

E (Z ) ≅ E (Z ) + ωl, k Z (l, k )

(17a)

E (Z j ) ≅ E (Z j ) + ωl, k (Z (l, j )) j

(17b)

× 100[%] (12)

3.2.1. 2m scheme (K = 2) This scheme yields from solving (14) for K = 2 to obtain the standard locations and weights of a random variable Pl . Hence, we use two concentrations for each input random variable to achieve ζl, k and ωl, k using the following equations:

j

where Z i is the mean random variable of Monte Carlo simulation at bus/node ith for the test group jth, given a number of simulations NMCS. Hence, the Eq. (12) applied the MCS to perform a given number of simulations NMCS (e.g. 7000 simulations) several times ng (e.g. 50 times). In this approach, the sensitivity of the number of iteration for presented MCS has been evaluated in Multi-carrier Energy Systems. The Probabilistic Energy Flow (PEF) has been solved for different iterations of MCS (100, 500, 1000, 3000, 5000, 7000, 10,000). It should be noted that by utilizing the Eq. (12), adequate results have been achieved. The Mean and variance have been used as Indices to evaluate the MCS accuracy.

ζl,1 =

λl,3 + 2

ωl,1 = −

2

λl,3 ⎞ m+⎛ ⎝ 2 ⎠ ⎜

ζl,2 =

⎟

ζ1,2 1 m ζl,1 − ζ1,2

ωl,2 =

λ1,3 − 2

2

λ1,3 ⎞ m+⎛ ⎝ 2 ⎠ ⎜

⎟

ζl,2 1 m ζl,1 − ζl,2

(18a)

(18b)

Moreover, the mean and standard deviation of Pl , the locations Pl,1 and Pl,2 can be derived by (13). The 2 m scheme always yields real value solutions for the concentrations. In addition, it has important benefits such as its simplicity and its low computational burden.

3.2. Point estimate method Point estimate method prepared by the first few central moments of a problem input random variable on K points for each variable, named concentration. The Kth concentration of a random variable Pl can be described as (Pl, k , ωl, k ) . The location Pl, k is the Kth value of the variable Pl at which the function F is evaluated. The weight ωl, k is a weighting factor, which accounts for the relative significance of this evaluation in the output random variables. According to this scheme, the function F, which relates input and output variables, has to be solved K times for each input random variable Pl . The location Pl, k can be derived by:

Pl, k = μpl + ζl, k σpl

(15a) ∞

of output random variables Z (l, k ) . Therefore, according to the weighting factors ωl, k and the Z (l, k ) values, the jth moment of the output random variables can be achieved using the following equation:

(11b)

Zi j−1 − Zi j

(σpl ) j

estimated at the points μ p1 , μ p2 , ... ,μ pl, k , ... ,μ pm affording the vector

For a precise responding, many simulations should be regarded. Since the calculation time of the Monte Carlo simulation in MES is longer as the number of simulations enlarges, it is significant to determine the accuracy of the results. This can be calculated by means of the relative error of Monte Carlo method for a given number of simulations [41]:

1 = (ng − 1). NMCS

(14b)

where λl,3 and λl,4 are known as the coefficient of skewness and kurtosis of Pl respectively. Now by making all the concentrations (Pl, k , ωl, k ) , the function F is

(11a)

k=1

j = 1, ... , 2k − 1

Mj (pl )

Mj (pl ) =

n

∑ (Z (k) )2 − (E (Z ))2

ωl, k (ζl, k ) j = λl, j

The jth central moment of the random variable Pl with PDF of fpl is calculated as follow:

n k=1

(14a)

k=1

(10)

∑ Z (k)

1 m

K

∑

where the mean and variance of outputs are calculated as:

E (Z ) =

ωl, k =

k=1

1. 2m+ 1 scheme (k = 3, ζl,3 = 0 ) The 2m+ 1 scheme utilizes one more assessment of function F than the 2m scheme. Hence by solving equation (14) for K= 3 and (ζl,3 = 0 ), the standard location ζl, k and the weight ωl, k can be achieved by:

(13)

ζl, k =

λl,3 3 + (−1)3 − k λ1,4 − λl2,3 2 4

ωl, k =

(−1)3 − k ζl, k (ζl,1 − ζ1,2)

,

ω1,3 =

k = 1, 2

ζ1,3 = 0

1 1 − m λ1,4 − λl2,3

(19a)

k = 1, 2 (19b)

It should be noted that in (19a) the standard location values of the 2 m+ 1 scheme do not relate to the number of input random variables (m), as do the K×m type schemes. In addition, the 2 m+ 1 scheme is

where ζl, k is the standard location, and the mean value and standard deviation of the input random variable Pl are μpl and σpl respectively. 995

Renewable and Sustainable Energy Reviews 94 (2018) 989–997

H. Khorsand, A.R. Seifi

electrical sub-network. The utilized natural gas sub-network and DHN are the same as case 1. Table 11 indicates the configurations and connections between units of Multi-carrier Energy Systems. The probabilistic data of electrical subnetwork such as load demand can be determined as follows. The active and reactive power of the load buses are modeled as normal distributions, whose means equal the base case data, and whose standard deviations are set arbitrarily as follows: 7% for bus 2,3,4,5,7,8, 4% for bus 10,12,14, 9% for bus 15 to bus 21, and bus 23,24, 5% for bus 26,29,30. The precision and efficiency of the proposed point estimate schemes are evaluated by comparing the results achieved with those derived from the Monte Carlo simulation with 7000 samples. Results have been shown in Table 18. It can be seen that 7000 are the appropriate number and raising the number of iteration only causes the simulation more time-consuming. In this case study with 53 input random variables, the number of estimations in the 2m or 2m+ 1 PEM scheme required to achieve a solution is 106 or 107 respectively. Tables 12–17 show the mean and standard deviation results of the case study 2 for a selected set of output random variables, which have been regarded as a sample of general results, achieved from different point estimate schemes. It can be mentioned that the 2m+ 1 scheme prepares good results compared with the Monte Carlo simulation, both for the mean and the standard deviation. The 2 m scheme mean results also produce adequate values, except for node 9&18 of Natural Gas pressure and mass flow rate of water in line(1−2) and line(2−3); but the differences increase for the standard deviation values, in comparison with Monte Carlo simulation results.(Table 18) Tables 10 and 19 show the average error indices relating to the estimations of the 2m and 2m+ 1 schemes presented in this paper for the Cases 1&2 test systems, respectively. The average error indices of the 2m+ 1 scheme are much smaller than those corresponding to the estimations of the 2 m scheme for both cases. For example, with taking a look at mean values achieved by the 2m+ 1 scheme in Case 2, the largest average error corresponds to the mass flow rate of water (ṁ Hline . ij ) and does not exceed 36%. Although, in the case of the 2m scheme, the highest average error corresponds to (ṁ Hline . ij ) and is about 127%. This means that the 2m+ 1 scheme yields better results than 2 m scheme. It can also be noted that the standard deviation average errors are greater than the mean ones in almost all the cases. This difference is especially remarkable for the case of the 2 m scheme applied to both case studies. Note, by comparing Tables 10 and 19, that the standard deviation average errors of the 2 m scheme increase when the total number of input random variables goes from 33 in Case 1–53 in Case 2 test system. As an example, in Case 1 the average error incurred by this scheme for the standard deviation of the voltage angle is 3796.19%, while this error reaches a value of 12,216.61% in Case 2. These simulations have been conducted on an Intel 2.4 GHz Core 2 Duo PC with 2 GB of RAM. Table 20 presents the CPU time required to compute the output variables for each scheme considered in this paper. It should be noted that the 2 m+ 1 scheme is faster than the Monte Carlo simulation and 2 m scheme.

more precise than the 2 m scheme since it uses the kurtosis λl,4 of the input random variables [26]. Equation (20) describes a general overview of the overall performance of the proposed method, for each RV:

εμRV =

εσRV =

RV RV μMCS − μPEM RV μMCS

× 100[%] (20a)

RV RV σMCS − σPEM × 100[%] RV σMCS

RV μMCS

(20b)

RV σMCS

and are the mean and standard deviation obtained by Where Monte Carlo simulation and, hence, taken as reference values. RV RV and σPEM are the mean and standard deviation estiSimilarly, μPEM mated by a proposed point estimate method. Moreover, the average error indices εμRV and εσRV are defined as the mean values of εμRV and εσRV respectively for each set of random variables, where RV may refer to voltages Vj , angles θj , Natural Gas pressures πj , Mass flow rate of water ṁ Hline . ij , Supply temperature at the end of pipelines Ts . j , or Return temperature at the end of pipelines Tr . j . 4. Simulations and results A Multi-carrier Energy System with an electrical, a natural gas, and a DHN is presented to illustrate the overall performance of the point estimate method presented in this paper. Hence, the probabilistic energy flow problem is solved for two case studies as follows. 4.1. Case 1 Schematic of this case study is illustrated in Figs. 2–4. For the electrical sub-network, the standard IEEE 14-bus system is used [20]. It has two generators that are assumed to consume natural gas in this paper. For the natural gas sub-network, the Belgium natural gas transmission system is applied (Fig. 3). It has 20-nodes, 24 pipelines, 6 gas sources, and 2 compressors driven by the electrical sub-network. Fig. 4 indicates the DHN with 14 nodes, which has a boiler as a slack source of the grid at node 1. This paper investigates four energy hubs in the Multi-Carrier Energy Systems. The Schematic of these hubs is shown in Fig. 1. Their parameters are reported in [20] by the authors. Configurations and connections between units of these sub-networks are shown in Table 1. In the MES network, these parameters were supposed to be stochastic with particular PDFs. The probabilistic data of electrical sub-network can be found in [42]. Table 2 presents nodal probabilistic data for natural gas and DHN of the test system. MCS with 5000 iterations has been applied to the MES network. It can be seen from Table 9 that 5000 are the proper number for Case 1 and increasing the number of iteration only causes the simulation more time-consuming. The outcomes of the simulations have been utilized as a reference for comparison with point estimate method. In the PEM with 33 input random variables, the number of estimations in the 2 m or 2m+ 1 PEM scheme required to obtain a solution is 66 or 67 respectively. Therefore, to illustrate the overall performance of the point estimate method a selection of typical value for the mean and standard deviation of the energy flow of Multi-carrier Energy Systems has been shown in Tables 3–8 via three different methods. (Tables 9 and 10)

5. Conclusion Point Estimate Method uses deterministic procedures in order to solve probabilistic problems by applying only some statistical moments of probability functions such as mean, variance, skewness, and kurtosis coefficients. Hence, PEM overcomes the difficulties associated with the lack of perfect knowledge of the probability functions of stochastic variables. So, good results can be achieved to handle stochastic problems in Multi-carrier Energy systems while keeping lower the computational burden involved. Two different point estimate schemes considered in this paper have been tested on the probabilistic energy flow problem of electrical, natural gas, and DHN simultaneously. Results over two different case

4.2. Case 2 To further analyze the probabilistic energy flow of Multi-carrier Energy Systems, a much larger network is applied for the second case study. For the electrical sub-network, the standard IEEE 30-bus system is used [43]. It has six generators, four tap-transformers, and two shunt capacitors. It is supposed that generators of bus 1 and 2 are gas-powered units while others are coal-powered generators. Fig. 5 indicates the 996

Renewable and Sustainable Energy Reviews 94 (2018) 989–997

H. Khorsand, A.R. Seifi

studies have been presented and compared against those obtained from the Monte Carlo simulation. It can be noticed that the PEM provides good results compared with the Monte Carlo simulation, both for the mean and the standard deviation. In addition, the average error indices of the 2m+ 1 scheme are much smaller than those obtained from the 2 m scheme for both case studies. Furthermore, by increasing the number of input random variables, the 2m+ 1 scheme yields adequate results while the 2 m scheme accuracy decreases. Hence, it seems that the 2 m scheme is not sufficient for large-scale system problems. Therefore, the 2m+ 1 scheme in this approach is somehow better than the MCS and 2m scheme since it has the less computational burden and consumed time due to a smaller volume of data is requested.

production and energy storage. Appl Energy 2015;159:401–21. [19] Rong A, Figueira JR, Lahdelma R. An efficient algorithm for bi-objective combined heat and power production planning under the emission trading scheme. Energy Convers Manag 2014;88:525–34. [20] Shabanpour-Haghighi A, Seifi AR. An integrated steady-state operation assessment of electrical, natural gas, and district heating networks. Power Syst IEEE Trans 2016;31:3636–47. [21] Shabanpour-Haghighi A, Seifi AR. Energy flow optimization in multi-carrier systems. Ind Inform IEEE Trans 2015;11:1067–77. [22] Geidl M, Andersson G Optimal power dispatch and conversion in systems with multiple energy carriers. In: Proceedings of the 15th power systems computation conference. PSCC; 2005. p. 1–7. [23] Shabanpour-Haghighi A, Seifi AR, Niknam T. A modified teaching–learning based optimization for multi-objective optimal power flow problem. Energy Convers Manag 2014;77:597–607. [24] Mohammadi M, Shayegani A, Adaminejad H. A new approach of point estimate method for probabilistic load flow. Electr Power Energy Syst 2013;51:54–60. [25] Morales M, Pérez-Ruiz J. Point estimate schemes to solve the probabilistic power flow. Power Syst IEEE Trans 2007;22:1594–601. [26] Su C. Probabilistic load-flow computation using point estimate method. Power Syst IEEE Trans 2005;20:1843–51. [27] Caramia P, Carpinelli G, Varilone P. Point estimate schemes for probabilistic threephase load flow. Electr Power Syst Res 2010;80:168–75. [28] Schellenberg A, Rosehart W, Aguado J. Cumulant-based probabilistic optimal power flow (P-OPF) with gaussian and gamma distributions. Power Syst IEEE Trans 2005;20:773–81. [29] Xiao Q. Comparing three methods for solving probabilistic optimal power flow. Electr Power Syst Res 2015;124:92–9. [30] Shabanpour-Haghighi A, Seifi AR. Simultaneous integrated optimal energy flow of electricity, gas, and heat. Energy Convers Manag 2015;101:579–91. [31] Chicco G, Mancarella P. A unified model for energy and environmental performance assessment of natural gas-fueled poly-generation systems. Energy Convers Manag 2008;49(8):2069–77. [32] Marangon Lima JW, Zambroni de Souza AC Modeling the Integrated Natural Gas and Electricity Optimal Power Flow. In: Proceedings of the power engineering society general meeting. IEEE; 2007. p. 1–7. [33] Shabanpour-Haghighi A, Seifi AR. Multi-objective operation management of a multi-carrier energy system. Energy 2015;88:430–42. [34] Hopper G, Robison M, Succar S, Vidashttps H New england sets the pace for electricity and natural gas integration. 〈http://www.icf.com/perspectives/whitepapers/2015/new-england-sets-pace-electricity-natural-gas-integration/〉. [35] Genon G, Torchio MF, Poggio A, Poggio M. Energy and environmental assessment of small district heating systems: global and local effects in two case-studies. Energy Convers Manag 2009;50(3):522–9. [36] Benonysson A, Bohm B, Ravn HF. Operational optimization in a district heating system. Energy Convers Manag 1995;36(5):297–314. [37] Dall'Anese E, Mancarella P, Monti A. Integrated optimization and control of multienergy systems. IEEE Power Energy Mag 2017;15:43–52. [38] Sundberg G, Henning D. Investments in combined heat and power plants: influence of fuel price on cost minimized operation. Energy Convers Manag 2002;43:639–50. [39] Mancarella P. Cogeneration systems with electric heat pumps: energy shifting properties and equivalent plant modelling. Energy Convers Manag 2009;50(8):1991–9. [40] Nesheim SJ, Ertesvag IS. Efficiencies and indicators defined to promote combined heat and power. Energy Convers Manag 2007;48:1004–15. [41] Ruiz-Rodriguez F, Herna andndez J, Jurado F. Probabilistic load flow for photovoltaic distributed generation using the Cornish–Fisher expansion. Electr Power Syst Res 2012;89:129–38. [42] Allan RN, Al-Shakarchi MRG. Probabilistic techniques in a.c. load flow analysis. Proc IEEE 1977;124:154–60. [43] IEEE 30-bus test system. 〈http://ee.washington.edu/research/pstca/pf30/pg_ tca30bus.htm/〉.

References [1] Geidl M, Andersson G. Optimal power flow of multiple energy carriers. Power Syst IEEE Trans 2007;22:145–55. [2] El-Khattam W, Hegazy YG, Salama MMA. Investigating distributed generation systems performance using Monte Carlo simulation. Power Syst IEEE Trans 2006;21:524–32. [3] Chen P, Chen Z, Bak-Jensen B Probabilistic load flow: a review. In: Proceedings of the third international conference on electric utility deregulation and restructuring and power technologies. DRPT, IEEE; 2008. p. 1586–91. [4] Ruiz-Rodriguez F, Herna andndez J, Jurado F. Probabilistic load flow for radial distribution networks with photovoltaic generators. Renew Power Gener IET 2012;6:110–21. [5] Verbic G, Cañizares CA. Probabilistic optimal power flow in electricity markets based on a two point estimate method. Power Syst IEEE Trans 2006;21:1883–93. [6] Soleimanpour N, Mohammadi M. Probabilistic load flow by using nonparametric density estimators. Power Syst IEEE Trans 2013;28(4):3747–55. [7] Carpinelli G, Di Vito V, Varilone P. Multi-linear Monte Carlo simulation for probabilistic three-phase load flow. Eur Trans Electr Power 2007;17:1–19. [8] Li X, Li Y, Zhang S. Analysis of probabilistic optimal power flow taking account of the variation of load power. Power Syst IEEE Trans 2008;23:992–9. [9] Tamtum A, Schellenberg A, Rosehart W. Enhancements to the cumulant method for probabilistic optimal power flow studies. Power Syst IEEE Trans 2009;24:1739–46. [10] Chen C, Wu W, Zhang B, Sun H. Correlated probabilistic load flow using a point estimate method with Nataf transformation. Electr Power Energy Syst 2015;65:325–33. [11] Delgado C, Domínguez-Navarro JA. Point estimate method for probabilistic load flow of an unbalanced power distribution system with correlated wind and solar sources. Electr Power Energy Syst 2014;61:267–78. [12] Abdullah MA, Agalgaonkar AP, Muttaqi KM. Probabilistic load flow incorporating correlation between timevarying electricity demand and renewable power generation. Renew Energy 2013;55:532–43. [13] Mancarella P. MES (multi-energy systems): an overview of concepts and evaluation models. Energy 2014;65:1–17. [14] Geidl M, Andersson G. Operational and structural optimization of multi-carrier energy systems. Electr Power Eur Trans 2006;16(5):463–77. [15] Moeini-Aghtaie M, Abbaspour A, Fotuhi-Firuzabad M, Hajipour E. A decomposed solution to multiple-energy carriers optimal power flow. Power Syst IEEE Trans 2013;29:707–16. [16] Woldeyohannes A, Abd Majid MA. Simulation model for natural gas transmission pipeline network system. Simul Model Pract Theory 2011;19:196–212. [17] Martínez-Mares A, Fuerte-Esquivel C. A unified gas and power flow analysis in natural gas and electricity coupled networks. Power Syst IEEE Trans 2012;27(4):2156–66. [18] Haichao Wanga C, Wusong Yin b, Abdollahi E, Lahdelma R, Jiao W. Modelling and optimization of CHP based district heating system with renewable energy

997