Risk Assessment of Industrial Energy Hubs and Peer-to-Peer Heat and Power Transaction in the Presence of Electric Vehicles

Abstract

The peer-to-peer (P2P) strategy as a new trading scheme has recently gained attention in local electricity markets. This is a practical framework to enhance the flexibility and reliability of energy hubs, specifically for industrial prosumers dealing with high energy costs. In this paper, a Norwegian industrial site with multi-energy hubs (MEHs) is considered, in which they are equipped with various energy sources, namely wind turbines (WT), photovoltaic (PV) systems, combined heat and power (CHP) units (convex and non-convex types), plug-in electric vehicles (EVs), and load-shifting flexibility. The objective is to evaluate the importance of P2P energy transaction with on-site flexibility resources for the industrial site. Regarding the substantial peak power charge in the case of grid power usage, this study analyzes the effects of P2P energy transaction under uncertain parameters. The uncertainties of electricity price, heat and power demands, and renewable generations (WT and PV) are challenges for industrial MEHs. Thus, a stochastically based optimization approach called downside risk constraint (DRC) is applied for risk assessment under the risk-averse and risk-neutral modes. According to the results, applying the DRC approach increased by 35% the operation cost (risk-averse mode) to achieve a zero-based risk level. However, the conservative behavior of the decision maker secures the system from financial losses despite a growth in the operation cost.

1. Introduction

In the power market, the role of local energy systems including energy storage units, wind farms, distributed energy resources, and solar photovoltaic (PV) systems has become significant [1]. The development of smart grid facilities along with the energy and prosumers communities has further accelerated the trend [2]. Hence, consumer-centric types of energy systems with modern market designs are required for future power systems to adapt to local energy systems and buildings for the management of DERs. An emerging method is to develop smaller units and collect them in energy hubs [3]. Therefore, due to the differences in the trading price of energy as well as losses, sharing DERs on a local scale can be quite effective [4]. The peer-to-peer (P2P) strategy is an emerging alternative that encourages neighbors in a community to share excess energy to manage peak power demands [5], through which both consumers and prosumers can be supplied. In this regard, the effectiveness of the DERs, self-consumption, local energy balance, and grid operation flexibility can be strengthened [6]. The P2P concept can promote not only the flexibility of a system (storage unit, demand response program) but also energy transaction based on local prices [7]. Hence, the grid utility tariff is a market-based feature that potentially influences P2P trading. For instance, the peak demand charge is a promising solution in Norway by which consumers (commercial and industrial) are incentivized to reduce their power demand, as they are already subject to the peak demand charge [8]. Due to the high level of energy, large consumers can supply a major part of their energy consumption using a distribution network [9].

1.1. Research Review

As P2P energy transaction is a new concept, there is no agreement on a pricing scheme or market design aiming to support the development of local markets. In this field, recent research has mainly focused on several perspectives. Firstly, the role of storage units in coordinating local resources is investigated in [10] to create more balance in the system functioning. Secondly, a bidding mechanism for local energy trading is developed in [11]. Then, while the digitalization of the system components is carried out in [12], the requirement of computational properties as well as the coordination algorithm is taken into account for a P2P structure in [13]. Many of the research works related to P2P power transaction are rarely considered for real-life projects, including the Enerchain [14], Brooklyn Microgrid [15], and Sonnen Community [16]. A further step is taken in [17] by considering possible market frameworks for the consumers and prosumers of a community participating in the electricity market. The development of smart grid technologies will facilitate the establishment of local P2P energy trading with consumer-based electricity markets having access to the wholesale electricity market [18]. In this regard, blockchain technology can be taken into account to create an affordable and secure platform for energy transaction [19].

Several pieces of research have also been conducted on the management of multi-energy hubs. In [20], a two-stage stochastic programming approach is used to manage several energy hubs, in which the reliability aspect of the system is investigated. In [21], the authors introduced a new method for the optimal design and operation of several energy hubs concerning the cost of cables and operation costs using a two-stage stochastic optimization method. In [22], researchers increased the operational flexibility of several microgrids by using P2P power transaction. In addition to the incrementing of renewable penetration, the emission cost of the system is intended to be minimized. Random EV charging as an uncertainty is investigated in [23] for the energy management of several EHs to minimize the operation cost along with the power losses and greenhouse gas emissions in a multi-objective study. In [24], the authors propose a software-defined grid system to facilitate energy sharing in an MEH using a transactive energy framework. The result of this implementation is the reduction of the overall cost of the EHs. In [25], the chance constraint method is applied for the optimal day-ahead planning of an MEH, considering environmental constraints.

In a P2P strategy, thermal energy generated by CHP units can improve not only the energy flexibility but the proficiency of the combined heat and power unit. In [26], the transaction of power and heat among smart energy hubs is investigated in a two-stage process. While energy trading is considered in the first stage, financial issues related to power and heat energy are evaluated in the second stage to ensure the privacy of the EHs. In [27], a genetic algorithm is used to analyze P2P multi-energy sharing to increase the local energy balance and minimize the energy cost of the proposed system. In [28], a new P2P trading model is taken into account among FC-CHP systems to enhance the resiliency and self-sufficiency of local energy systems by the alternating direction method of multipliers (ADMM) algorithm. In [29], the effectiveness of a novel P2P energy transaction of multiple energy hubs is validated by a cooperative game to establish a proper payoff allocation.

In the existing literature, the importance of P2P power trading among residential prosumers has been mainly considered to evaluate the operational and economic aspects of systems. However, a few of them focused on the importance of P2P heat and power trading among industrial energy hubs and the functions of system components under uncertainties.

1.2. Novelties and Contributions

In this paper, a framework of a Norwegian industrial multi-energy hub (MEH) is developed, in which the hubs are equipped with various energy resources to supply their power and heat demands and share their excess energy with other peers and the electricity network. The CHP unit, solar photovoltaic (PV) system, wind turbine (WT), and plug-in electric vehicles (PEVs) are the generation facilities being used in the industrial hubs. Also, load shifting is a flexibility asset that is considered for two hubs, upon which the decision-maker can shift a part of the energy demand from peak demand to valley demand. However, the performance of stochastic renewable generations as well as the thermal and electrical demands with market electricity price in the energy hubs (EHs) have potential effects on the optimal function of the system. In this concern, a scenario-based stochastic optimization procedure called the downside risk constraint (DRC) approach is applied to investigate the functions of the system components under uncertainties. The risk-averse ($\lambda =1$) and risk-neutral ($\lambda =0$) modes are used to forecast the impact of uncertain parameters, concisely. The contribution of this paper can be summarized below:

❖

Techno-economic analysis of an industrial MEH with P2P heat and power transaction.

❖

Risk analysis of an industrial MEH with the downside risk constraint method (DRC) as a stochastic optimization procedure.

❖

Load-shifting flexibility asset and distributed energy resources, namely WTs, PVs, convex and non-convex CHP units, and plug-in electric vehicles (PEVs) to support energy demands.

In order to achieve a better perception, the overall structure of the proposed system is represented in Figure 1.

1.3. Paper Structure

This paper is organized as follows: Section 2 expresses the objective function as well as the mathematical modeling of the MEHs. The constraints related to the DRC method are presented in Section 3. The study case along with the input data and the simulation results of the MEH are shown in Section 4. Finally, the conclusion of this paper is discussed in Section 5.

2. Mathematical Formulation

The mathematical model of the MEH is a linear optimization model representing the P2P interaction of the EHs and the operational decisions related to the DERs and flexibility assets. The objective function is to minimize the operational cost and analyze the risk level in the industrial EHs that are equipped with the CHP units and other distribution resources.

2.1. Objective Function

Based on the market features, the industrial EHs have the opportunity to trade their excess power through local P2P transaction. However, the operational cost arises in the cases of grid power consumption, load-shifting practice, and importing electricity from other peers. As there is a capability of trading power with the network and other peers, the units can obtain an income from grid feed-in and exporting electricity to a peer, which affects the optimal operation of the individual hubs. As shown in Equation (1), the cost of trading power with the power network is given in the first three terms. While the fourth term indicates the operational cost of the CHP units, the load-shifting cost is shown in the fifth term. Finally, the power and thermal energy transaction among the industrial peers is demonstrated in the last four terms.

$$\underset{\begin{array}{l}\forall t\in T\\ \forall b\in B\end{array}}{\min }C=\left\{{\displaystyle \sum _b^B\left(\begin{array}{l}{\displaystyle \sum _t^T\left[\left(C_{g,eng}+C_{(t)}^{g,SP}\right)\cdot P_{(t,h)}^{g,buy}\right]}+{\displaystyle \sum _m^M\left[C_{g,fix}+C_{g,peak}\cdot P_{(t,h)}^{g,peak}\right]}-\\ {\displaystyle \sum _t^T\left[C_{(t)}^{feed-in}\cdot P_{(t,h)}^{g,sell}\right]}+{\displaystyle \sum _t^T{\displaystyle \sum _h^H{\displaystyle \sum _c^C\left[a\_chp\times P_{t,h,c}^{chp}+b\_chp+c\_chp\times T_{t,h,c}^{chp}\right]}}}+\\ {\displaystyle \sum _t^T\left[C_{(h)}^{LS}\cdot P_{(t,h)}^{ls,sh}\right]}+{\displaystyle \sum _t^T\left[C_{(t,h)}^{P,P2P}\cdot P_{(t,h)}^{imp}\cdot \frac{1}{{\psi }_{P,P2P}}\right]}-{\displaystyle \sum _t^T\left[{\displaystyle \sum _{p\ne h}^H{C}_{(t,p)}^{P,P2P}\cdot P_{(t,h\to p)}^{\exp ,p}}\right]}+\\ {\displaystyle \sum _t^T\left[C_{(t,b)}^{Th,P2P}\cdot T_{(t,h)}^{imp}\cdot \frac{1}{{\psi }_{Th,P2P}}\right]}-{\displaystyle \sum _t^T\left[{\displaystyle \sum _{p\ne h}^H{C}_{(t,p)}^{Th,P2P}\cdot T_{(t,h\to p)}^{\exp ,p}}\right]}\end{array}\right)}\right\}$$

(1)

The overall cost of power network electricity for the Norwegian industrial MEHs is presented in Equation (2). As the local network company is responsible for determining the utility tariff, the utility tariff system in Norway differs from the flat rate tariffs towards time-of-use pricing, which has a peak-demand cost [30]. As shown in the equation, the first and second parenthetical terms represent buying electrical energy at a spot price and the price of energy during the peak hours, respectively. Moreover, in the second term, feed-in energy cost is calculated and considered for selling power to the power grid.

$$C_{(h)}^{g.tot}={\displaystyle \sum _{t\in T}\left(C_{(t)}^{g,SP}\cdot P_{(t,h)}^{g.buy}+C_{g,eng}\cdot P_{(t,h)}^{g.buy}\right)}+\left(C_{g,fix}+C_{g,peak}\cdot P_{(h)}^{g.peak}\right)-{\displaystyle \sum _{t\in T}\left(C_{(t)}^{feed-in}{P}_{(t,h)}^{g.sell}\right)}$$

(2)

2.2. Energy Balance

P2P power transaction between the industrial hubs will affect the balance of the system. In this regard, the constraint of power balance is shown in Equation (3), in which the total amount of electrical demand is equal to the overall level of power generated by the DERs. Also, thermal energy transaction is considered to increase the energy flexibility, though the related constraint is given in Equation (4). As mentioned in the equation, the overall amount of heat energy produced by CHP units along with the imported power from other peers, bought power from the main grid, and discharged energy of PEVs must be equal to the overall level of electrical demand. Also, the thermal energy produced by the CHP unit and imported from other peers must be equivalent to the demand and the exported heat energy. Finally, the limitations of trading power with a power network are shown in Equations (5) and (6).

$$\begin{array}{l}{\displaystyle \sum _{c=1}^C{P}_{(t,h,c)}^{chp}}+P_{(t,h)}^{dem}+P_{(t,h)}^{g.sell}+P_{(t,h)}^{\exp }+P_{(t,h)}^{ev,ch}+P_{(t,h)}^{ls,dem}+\\ P_{(t,h)}^{curtail}=P_{(t,h)}^{DER}+P_{(t,h)}^{g.buy}+P_{(t,h)}^{imp}+P_{(t,h)}^{ev.dch}+P_{(t,h)}^{ls.sh}\end{array}$$

(3)

$${\displaystyle \sum _c^C\left(T_{t,h,c}^{chp}\right)}+T_{t,h}^{imp}=T_{t,h}^{dem}+T_{t,h}^{exp}$$

(4)

$$0\le P_{(t,h)}^{g.buy}\le P_{(h)}^{g.peak}$$

(5)

$$0\le P_{(t,h)}^{g.sell}\le P_{feed-in}^{\max }$$

(6)

2.3. Constraints of P2P Energy Transaction

As the interconnected industrial hubs have the opportunity to trade energy among each other, specific market mechanisms defined by local markets can secure the trading process [31]. As the proposed P2P energy transaction methodology is a general model, the implementation of the system on different decentralized platforms like blockchains can be carried out according to [32,33]. The total sums of exported and imported power by hub h are shown in Equations (7) and (8), respectively. As indicated in Equation (9), the imported electricity from a peer must be equal to the electricity exported from the peer to the EH, which involves power losses (${\psi }_{P2P}$). Also, the total level of traded power among the EHs is given by Equation (10). The mentioned constraints are also indicated for thermal energy transaction between Equations (11) and (14)

$$P_{(t,h)}^{\exp }={\displaystyle \sum _{p\ne h}{P}_{(t,h\to p)}^{\exp ,p}}$$

(7)

$$P_{(t,h)}^{imp}={\displaystyle \sum _{p\ne b}{P}_{(t,h\leftarrow p)}^{imp,p}}$$

(8)

$$P_{(t,h\leftarrow p)}^{imp,p}=P_{(t,h\to p)}^{\exp ,p}\times {\psi }_{P,P2P},\forall p\ne h$$

(9)

$${\displaystyle \sum _h^H{P}_{(t,h)}^{\exp }\times {\psi }_{P,P2P}}={\displaystyle \sum _h^H{P}_{(t,h)}^{imp}}$$

(10)

$$T_{(t,h)}^{\exp }={\displaystyle \sum _{p\ne h}{T}_{(t,h\to p)}^{\exp ,p}}$$

(11)

$$T_{(t,h)}^{imp}={\displaystyle \sum _{p\ne b}{T}_{(t,h\leftarrow p)}^{imp,p}}$$

(12)

$$T_{(t,h\leftarrow p)}^{imp,p}=T_{(t,h\to p)}^{\exp ,p}\times {\psi }_{Th,P2P},\forall p\ne h$$

(13)

$${\displaystyle \sum _h^H{T}_{(t,h)}^{\exp }\times {\psi }_{Th,P2P}}={\displaystyle \sum _h^H{T}_{(t,h)}^{imp}}$$

(14)

2.4. Load-Shifting Constraints

Load shifting as a crucial strategy can help industrial hubs to reduce their operational cost by running a production process in the low–peak interval instead of the peak demand period. However, this strategy imposes productivity losses and labor rescheduling costs, both of which are considered penalty costs in the objective function. The mathematical modeling of load shifting is regarded in the form of a storage unit with 10% capacity to ease the computational burden [34], which is defined in the form of hourly rescheduled load in Equation (15). In this regard, while the storage balance for each EH is shown in Equation (16), their energy level with the maximum power shift is indicated in Equation (17).

$$0\le P_{(t,h)}^{ls,sh},P_{(t,h)}^{ls,dem}\le 0.1\times P_{(h)}^{g,peak}$$

(15)

$$E_{(t,h)}^{ls}=E_{(t-1,h)}^{ls}+P_{(t,h)}^{ls,sh}-P_{(t,h)}^{ls,dem}$$

(16)

$$0\le E_{(t,h)}^{ls}\le 0.4\times P_{(h)}^{g.peak}$$

(17)

2.5. Electric Vehicle (EV) Constraints

Vehicle-to-grid (V2G) is an on-site flexibility option providing a bi-directional use of electricity for EHs as a fast-responding storage unit that can be used for spinning reserve and peak shaving. As the industrial units hold a large number of employees, considering V2G technology for the parking lots can be an alternative flexibility asset. As given in Equation (18), the EV parking lot is a storage unit balancing the energy consumption of EHs. Equation (19) shows the limitation of the charging/discharging process by the nominal capacity of the EV charger. Finally, the start and end of a workday are limited by Equation (20) and Equation (21), respectively [35].

$$E_{(t,h)}^{ev}=E_{(t-1,h)}^{ev}+{\eta }_{ev,ch}\cdot P_{(t,h)}^{ev,ch}-\frac{1}{{\eta }_{ev,ch}}\cdot P_{(t,h)}^{ev,dch}$$

(18)

$$0\le P_{(t,h)}^{ev,ch},P_{(t,h)}^{ev,dch}\le P_{ev,ch\arg er}^{nom}\cdot E{V}_{num}$$

(19)

$$E_{(d_{start}(t),h)}^{ev}=E_{ev}^{nom}\cdot E{V}_{num}\cdot E_{start}{d}_{start}(t)\in T$$

(20)

$$E_{(d_{end}(t),h)}^{ev}\ge E_{ev}^{nom}\cdot E{V}_{num}\cdot E_{end}{d}_{end}(t)\in T$$

(21)

2.6. Constraints of CHP Units

Due to the high proficiency of CHP units, they can be used to supply both power and thermal energy in the EHs. The concept of the feasible operation region (FOR), which is shown in Figure 2, is used to model the constraints of the CHP units. The installed CHP units are convex and non-convex; hence, the start-up and shutdown limitations are given in Equations (22) and (23), respectively. While the convex unit is modeled by Equations (24)–(28), the model of the non-convex unit is formulated by Equations (29)–(34).

$$S{U}_{t,h,s,c}^{chp}\ge C_{h,c}^{SU}\left[V_{t,h,s,c}^{chp}-V_{t-1,h,s,c}^{chp}\right]$$

(22)

$$S{D}_{t,h,s,c}^{chp}\ge C_{h,c}^{SHD}\left[V_{t-1,h,s,c}^{chp}-V_{t,h,s,c}^{chp}\right]$$

(23)

$$P_{t,h,s,c}^{chp}-P_{h,c}^{chp,A}-\frac{P_{h,c}^{chp,A}-P_{h,c}^{chp,B}}{T_{h,c}^{chp,A}-T_{h,c}^{chp,B}}\left[T_{t,h,s,c}^{chp}-T_{h,c}^{chp,A}\right]\le 0$$

(24)

$$P_{t,h,c}^{chp}-P_{h,c}^{chp,B}-\frac{P_{h,c}^{chp,B}-P_{h,c}^{chp,C}}{T_{h,c}^{chp,B}-T_{h,c}^{chp,C}}\left[T_{t,h,c}^{chp}-T_{h,c}^{chp,B}\right]\ge -M$$

(25)

$$P_{t,s,c}^{chp}-P_{h,c}^{chp,C}-\frac{P_{h,c}^{chp,C}-P_{h,c}^{chp,D}}{T_{h,c}^{chp,C}-T_{h,c}^{chp,D}}\left[T_{t,s,c}^{chp}-T_{h,c}^{chp,C}\right]\ge -M$$

(26)

$$0\le P_{t,h,c}^{chp}\le P_{h,c}^{chp,A}\cdot V_{t,h,c}^{chp}$$

(27)

$$0\le T_{t,h,c}^{chp}\le T_{h,c}^{chp,B}\cdot V_{t,h,c}^{chp}$$

(28)

$$P_{t,h,c}^{chp}-P_{h,c}^{chp,B}-\frac{P_{h,c}^{chp,B}-P_{h,c}^{chp,C}}{T_{h,c}^{chp,B}-T_{h,c}^{chp,C}}\left[T_{t,h,c}^{chp}-T_{h,c}^{chp,B}\right]\le 0$$

(29)

$$P_{t,h,c}^{chp}-P_{h,c}^{chp,C}-\frac{P_{h,c}^{chp,C}-P_{h,c}^{chp,D}}{T_{h,c}^{chp,C}-T_{h,c}^{chp,D}}\left[T_{t,h,c}^{chp}-T_{h,c}^{chp,C}\right]\ge 0$$

(30)

$$P_{t,h,c}^{chp}-P_{h,c}^{chp,D}-\frac{P_{h,c}^{chp,D}-P_{h,c}^{chp,E}}{T_{h,c}^{chp,D}-T_{h,c}^{chp,E}}\left[T_{t,h,s,c}^{chp}-T_{h,c}^{chp,D}\right]\ge -\left[1-X_t^a\right]M$$

(31)

$$P_{t,h,c}^{chp}-P_{h,c}^{chp,D}-\frac{P_{h,c}^{chp,D}-P_{h,c}^{chp,F}}{T_{h,c}^{chp,D}-T_{h,c}^{chp,F}}\left[T_{t,h,c}^{chp}-T_{h,c}^{chp,D}\right]\ge -\left[1-X_t^b\right]M$$

(32)

$$0\le P_{t,h,c}^{chp}\le P_{h,c}^{chp,A}\cdot V_{t,h,c}^{chp}$$

(33)

$$0\le T_{t,h,c}^{chp}\le T_{h,c}^{chp,C}\cdot V_{t,h,c}^{chp}$$

(34)

2.7. Constraints of Renewable Energy Sources

Wind turbines and solar photovoltaic (PV) systems are renewable resources that are considered in this model. The generated available power of wind turbines and PV systems are modeled in Equation (35) and Equation (36), respectively. As shown in the equations, while wind speed potentially affects the power generation of WTs, solar radiation influences the PV system function.

$$P_{t,h}^{wt}=\{\begin{array}{ll}0& V_t<V_{cut-in},V_t\ge V_{cut-out}\\ P_{rated}^{wt}\times {\left(\frac{V_t-V_{cut-in}}{V_{rated}-V_{cut-in}}\right)}^3& V_{cut-in}\le V_t\le V_{rated}\\ P_{rated}^{wt}& V_{rated}\le V_t\le V_{cut-out}\end{array}$$

(35)

$$P_{t,h}^{pv}=\{\begin{array}{ll}{P}_{rated}^{pv}\left(\frac{I_t^2}{I_{std}{I}_C}\right)& I_t\le I_C\\ P_{rated}^{pv}\left(\frac{I_t^2}{I_{std}}\right)& I_C\le I_t\end{array}$$

(36)

3. Downside Risk Constraint (DRC) Method

The downside risk constraint (DRC) method is applied to control the risk of financial losses. Regarding decision variables, convexity is the most important feature of the DRC method. Compared to other risk measures, the DRC method has more robustness [36], with substantial advantages from the risk-management point of view. As the DRC has a convex function with a set of minimum points, it simplifies the optimization and control of the uncertainty through mathematical programming. As shown in Equation (37), when the cost of the scenario is less than the expected cost, the risk level is zero; otherwise, the level of risk can be obtained from the difference between the scenario cost and the expected cost. Also, the expected risk of the DRC method is given in Equation (38), in which the operator aims to achieve a small value for $DRC(C_0)$.

$$RISK(\omega )=\{\begin{array}{lll}\mathrm{Cos}t(\omega )-C_0& if& \mathrm{Cos}t(\omega )>C_0\\ 0& if& \mathrm{Cos}t(\omega )\le C_0\end{array}$$

(37)

$$DRC(C_0)=E\left[RISK\left(\omega \right)\right]={\displaystyle \sum _{\omega ={\Omega }}\pi \left(\omega \right)\cdot RISK\left(\omega \right)}$$

(38)

Equation (39) indicates that the scenario costs have more value than the expected cost. In this equation, the term $P\left(\omega |\mathrm{Cos}t\left(\omega \right)\ge C_0\right)$ shows the probability of a cost that is higher than the target cost. If the operator is not satisfied with the obtained risk level, a risk constraint like Equation (40) can be added to the main formulation as below:

$$DRC(C_0)=C_0-\frac{1}{P\left(\omega |\mathrm{Cos}t\left(\omega \right)\ge C_0\right)}\times {\displaystyle \sum _{\omega =1}^{{\Omega }}\left\{\pi \left(\omega \right)\cdot \max \left[\left(\mathrm{Cos}t\left(\omega \right)-C_0\right),0\right]\right\}}$$

(39)

$$DRC(C_0)\le DR{C}_0$$

(40)

where the term $DR{C}_0$ indicates the tolerance of the downside risk constraint. The flowchart of the proposed method with the strategy of the system operator is demonstrated in Figure 3.

As shown, the initial data related to the uncertain parameters are provided based on the historic data, upon which the scenario generation of uncertainties is performed by applying the Monte Carlo method for the wind turbine (WT), PV, electricity price, power, and thermal energy demands. In order to ease the computational burden, the Kantorovich method is used to achieve an optimal allocation of scenarios being used in the DRC programming.

4. Study Case

In this paper, a Norwegian industrial area with five energy hubs is considered, the prosumers of which are connected to each other to become more flexible by applying a P2P power and heat transaction strategy. These industrial units are related to food-processing industries, mechanical workshops, and manufacturing factories. In this regard, supporting their energy demand due to their high power consumption is a big challenge. The EHs are equipped with several energy sources, including PV systems and wind turbines as renewables, as well as EVs and CHP units. From a flexibility perspective, the load-shifting strategy increases the productivity of the energy hubs significantly. Intuitively speaking, load shifting means that the industrial building is willing to move demand from the peak demand period to valley demand, making a production process run at a later time. However, this strategy has rescheduling costs, namely overtime pay for laborers, rescheduling, and productivity losses. Moreover, the interconnected EHs have a connection with the main power grid, as the installed sources may not be able to meet the power demand of prosumers on their own. The industrial hubs are located in different places, and they are completely different in terms of load consumption and size [37]. As a final remark, the simulation of the proposed system is carried out through the GAMS software as a mixed-integer linear programming model with the CPLEX solver.

4.1. Input Data

This subsection presents the input data related to the generated scenarios for the uncertain parameters. The market electricity price, power demand, thermal energy demand, and renewable generations are taken into account as uncertainties. In this regard,

Table 1. The electricity price generated in 10 scenarios (Nok).

	SC = 1	SC = 2	SC = 3	SC = 4	SC = 5	SC = 6	SC = 7	SC = 8	SC = 9	SC = 10
t = 1	0.2125	0.2762	0.4314	0.3241	0.3631	0.3883	0.2539	0.3572	0.3039	0.2153
t = 2	0.3486	0.2504	0.2775	0.2111	0.2695	0.3523	0.4421	0.2276	0.3448	0.3589
t = 3	0.3888	0.4054	0.2231	0.3130	0.4031	0.3503	0.4162	0.4021	0.4889	0.1063
t = 4	0.2782	0.2776	0.1954	0.2638	0.4112	0.3078	0.3762	0.4273	0.2813	0.2723
t = 5	0.2078	0.2883	0.2545	0.2928	0.2923	0.2534	0.2211	0.2805	0.3109	0.2404
t = 6	0.3793	0.3360	0.3259	0.4281	0.5496	0.4225	0.4046	0.3490	0.3599	0.3993
t = 7	0.3250	0.2364	0.4177	0.3110	0.4144	0.2988	0.2590	0.2448	0.4037	0.5615
t = 8	0.403	0.4308	0.3975	0.4058	0.6106	0.3087	0.3684	0.3191	0.4522	0.4082
t = 9	0.359	0.4635	0.2927	0.4026	0.4105	0.5448	0.4436	0.3942	0.2596	0.2471
t = 10	0.4481	0.3221	0.3669	0.3056	0.4082	0.4127	0.5607	0.5679	0.4941	0.4963
t = 11	0.2697	0.5027	0.2914	0.2287	0.6829	0.4115	0.3180	0.3819	0.3877	0.1839
t = 12	0.3095	0.4268	0.3803	0.2202	0.3306	0.3628	0.3846	0.2187	0.4010	0.5211
t = 13	0.2244	0.4139	0.4983	0.4084	0.4279	0.3287	0.3557	0.2662	0.3616	0.3214
t = 14	0.3112	0.5073	0.3616	0.4311	0.4074	0.2780	0.1965	0.2389	0.3982	0.3444
t = 15	0.4501	0.3423	0.2532	0.3879	0.4824	0.4033	0.3606	0.3384	0.4982	0.4647
t = 16	0.3398	0.2525	0.2937	0.4124	0.2767	0.3438	0.3690	0.3179	0.3212	0.3591
t = 17	0.3301	0.2543	0.4533	0.2948	0.4819	0.3499	0.3534	0.5068	0.4078	0.3443
t = 18	0.3051	0.3345	0.4230	0.4660	0.3892	0.3372	0.4168	0.4491	0.3440	0.5233
t = 19	0.3493	0.4706	0.4773	0.3897	0.3761	0.4688	0.2110	0.3875	0.3909	0.4194
t = 20	0.2996	0.2491	0.1309	0.4372	0.4885	0.3420	0.6137	0.3232	0.2586	0.370
t = 21	0.3960	0.4517	0.3976	0.3122	0.2389	0.3292	0.2543	0.4556	0.3758	0.2071
t = 22	0.4198	0.2960	0.5145	0.3897	0.4067	0.3089	0.4350	0.3909	0.4055	0.3539
t = 23	0.4451	0.1496	0.3993	0.4839	0.3651	0.4388	0.3263	0.4000	0.3274	0.4932
t = 24	0.4254	0.5116	0.3767	0.3750	0.4105	0.5185	0.1532	0.2211	0.5013	0.5609

is given to show the price of electricity over the period of a 24 h scheduling horizon in 10 scenarios [37]. For each energy EH, 10 scenarios of power and heat demand are generated, which are shown in Figure 4 and Figure 5, respectively. In fact, the scenario generation is carried out based on the Monte Carlo method, in which a probability distribution is used to generate 100 scenarios for each uncertainty [38]. Whereas a large number of scenarios results in a computational burden, the Kantorovich procedure is applied to select 10 scenarios with high probability [39]. Based on the figures, hub5 has more energy demand compared to the other hubs. Meanwhile, in Figure 6 the PV function is demonstrated in hub1, hub4, and hub5, Figure 7 shows the function of WTs in hub2 and hub3. Also,

Table 2. Parameters related to the operation region of the CHP unit.

CHP Units	a ($/kW2)	b ($/kW2)	c ($)	d ($/kWth2)	e ($/kWth)	f ($/kW.kWth)	Feasible Region Coordinates
CHP 1	0.0345	44.5	26.5	0.03	4.2	0.031	[1.258 0], [1.258 0.324], [1.102 1.356], [0.4 0.75], [0.44 0.159], [0.44 0]
CHP 2	0.0435	56	12.5	0.027	0.6	0.011	[2.47 0], [2.15 1.8], [0.81 1.048], [0.988 0]

gives related data about the convex and non-convex CHP units [40].

4.2. Numerical Results

In this section, the levels of operational cost and risk are shown for different iterations. In Figure 8, the simulation result is obtained in 11 iterations, in which the lambda value changes from zero ($\lambda =0$) to one ($\lambda =1$). In this case, $\lambda =0$ and $\lambda =1$ imply the risk-neutral and risk-averse modes, respectively. Despite an increment in operational cost, the existing risk level decreases as we approach a high value of $\lambda$. In order to ease the understanding of this concept, the amount of risk in each iteration is shown in Figure 9. As shown, In the risk-neutral mode, the amount of risk is at its maximum value, but in the risk-averse mode, the operational cost is at its maximum amount (4.5 NOK), while the amount of risk for the system operator is flattened. The numerical results of the operational cost and risk level are demonstrated in

Table 3. Operational cost for each scenario and iteration (NOK).

n1	SC = 1	SC = 2	SC = 3	SC = 4	SC = 5	SC = 6	SC = 7	SC = 8	SC = 9	SC = 10
λ = 0	223,296	663,866	478,104	419,689	252,088	586,506	377,805	366,694	560,148	560,861
λ = 0.1	239,965	663,866	478,104	441,956	279,744	586,506	411,143	387,756	560,148	560,861
λ = 0.2	295,246	663,867	478,105	454,059	301,644	586,506	416,584	414,024	560,148	560,861
λ = 0.3	282,311	663,867	486,978	478,361	339,428	586,506	448,238	463,084	560,148	560,861
λ = 0.4	283,285	663,867	502,102	492,357	406,741	586,506	470,921	494,234	560,148	560,861
λ = 0.5	322,415	663,867	517,226	477,362	484,446	586,506	482,204	517,226	560,148	560,861
λ = 0.6	381,447	663,867	532,350	515,415	532,350	586,506	485,118	505,438	560,148	560,861
λ = 0.7	547,474	663,867	547,474	547,474	464,808	586,506	489,965	506,162	560,148	560,861
λ = 0.8	514,006	663,867	564,691	564,691	564,691	586,506	494,385	564,691	564,691	564,691
λ = 0.9	603,371	663,867	603,371	603,371	598,991	603,371	603,371	547,256	603,371	603,371
λ = 1	663,867	663,867	663,867	663,867	663,867	663,867	663,867	663,867	663,867	663,867

and

Table 4. Risk level for each scenario and iteration (NOK).

C	SC = 1	SC = 2	SC = 3	SC = 4	SC = 5	SC = 6	SC = 7	SC = 8	SC = 9	SC = 10
λ = 0	0	214,960	29,198	0	0	137,600	0	0	111,242	111,955
λ = 0.1	0	202,861	17,099	0	0	125,501	0	0	99,143	99,856
λ = 0.2	0	190,762	5000	0	0	113,402	0	0	87,044	87,757
λ = 0.3	0	176,888	0	0	0	99,528	0	0	73,170	73,883
λ = 0.4	0	161,764	0	0	0	84,404	0	0	58,046	58,759
λ = 0.5	0	146,641	0	0	0	69,280	0	0	42,922	43,635
λ = 0.6	0	131,517	0	0	0	54,156	0	0	27,798	28,511
λ = 0.7	0	116,393	0	0	0	39,032	0	0	12,674	13,387
λ = 0.8	0	99,176	0	0	0	21,815	0	0	0	0
λ = 0.9	0	60,495	0	0	0	0	0	0	0	0
λ = 1	0	0	0	0	0	0	0	0	0	0

, respectively. According to Table 3, the operational cost of the system for different lambdas is shown. When $\lambda =0$ (risk-neutral), a low level of cost is obtained for all scenarios. However, in the second scenario, which is the worst-case scenario, a high level of cost (663866 NOK) is obtained because there is a probability of scenarios with maximum financial losses.

Figure 10 and Figure 11 show the amounts of electrical and thermal energy exchanged among the EHs, respectively. According to Figure 10, all energy hubs have exported a large proportion of their generated power to meet the power demand of hub5. Because the load demand in hub5 is higher than in the other hubs (with a maximum value of 3000 kW), each of them has exported an amount of power between 1500 kW and 3000 kW to hub5. Also, the same trend is achieved for thermal energy transaction, as shown in Figure 11. As the evaluation of the energy transaction is carried out in both the risk-averse and risk-neutral modes, the level of transacted energy in the risk-averse mode is less than that in the risk-neutral mode due to the conservative behavior of the system operator.

Figure 12 shows the amount of electrical energy exchanged between the network and the EHs that are connected to the power network. Although the first, third, and fifth hubs bought electricity from the network in the risk-neutral mode, all hubs sold a portion of their energy to the grid which was between 400 kW and 800 kW. As shown in the figure, the EHs sold 70% more electrical energy in the risk-neutral mode compared to the risk-averse mdoe.

Figure 13 and Figure 14 demonstrate the power and heat energy produced by the CHP units in the EHs, respectively. By comparing the figures, it can be deduced that the second CHP unit (non-convex) generates more power, but the first CHP unit (convex) produces more thermal energy. The functioning of the CHP units is considerable in hub5 (CHP1 (1600 kW in risk-neutral mode and 1250 kW in risk-averse mode) and CHP2 (3400 kW in risk-neutral mode and 2600 kW in risk-averse mode)).

Figure 15 shows the charging and discharging statements of the EVs. Because we have considered the presence of EVs for 8 h in hub5, the functioning is obtained from 8 to 16. Due to the high peak demand, a significant amount of power has been delivered to hub5—about 20 kW. Finally, the function of the load-shifting strategy is demonstrated in Figure 16 for hub1 and hub3. This strategy is quite helpful in terms of energy cost reduction by moving the demand from peak demand to other time intervals. Based on this strategy, 20 kW and 15 kW of power loads in hub1 and hub3, respectively, are moved out of peak demand to reduce the energy cost of the hubs.

5. Results Validation

In this section, the generation of scenarios is carried out for different data representing a risky condition to validate the obtained results. Under such a condition, the standard deviation, as well as the fluctuation level, are a somewhat larger compared to the previous model. However, there is a slight difference when there is a comparison between uncertain parameters. In this regard, Figure 17, Figure 18, Figure 19, Figure 20 and Figure 21 indicate the generated scenarios for PV systems and WTs, market price, power demand, and thermal energy demand, respectively. As shown in

Table 5. Obtained scenario costs versus different levels of lambda.

	SC = 1	SC = 2	SC = 3	SC = 4	SC = 5	SC = 6	SC = 7	SC = 8	SC = 9	SC = 10
λ = 0	223,256	550,148	477,104	418,689	251,083	576,426	374,801	365,674	633,806	560,661
λ = 0.1	249,265	564,361	478,604	448,456	279,244	586,402	411,243	387,550	663,840	560,661
λ = 0.2	295,416	564,361	478,605	453,049	301,534	586,402	416,514	413,064	663,852	560,661
λ = 0.3	282,231	564,361	486,278	476,351	338,328	586,402	448,248	463,124	663,852	560,661
λ = 0.4	283,755	564,361	502,001	492,357	405,541	586,402	470,451	494,144	663,852	560,661
λ = 0.5	322,345	564,361	517,026	476,342	484,356	586,402	482,164	517,247	663,852	560,661
λ = 0.6	381,627	564,361	532,140	518,412	532,530	586,402	485,248	505,531	663,852	560,661
λ = 0.7	547,284	564,361	547,974	547,671	464,758	586,402	489,855	505,212	663,852	560,661
λ = 0.8	514,126	564,361	564,792	564,664	564,521	586,402	494,245	564,951	663,852	565,691
λ = 0.9	603,321	603,451	603,470	603,571	598,691	603,371	603,191	547,426	663,852	602,371
λ = 1	663,852	663,852	663,852	663,852	663,852	663,852	663,852	663,852	663,852	663,852

, the scenario costs are obtained for different lambda values, for which the operating costs were led to a certain number (663852 NOK) in the last iteration. In the previous section, however, the obtained value when $\lambda =1$ is 663867 NOK, representing the fact that in all uncertain situations, the scenario cost is quite close to our results.

6. Conclusions

In this paper, the value of P2P heat and power trading in combination with different resources of on-site flexibility is investigated for a Norwegian industrial site, where industrial units are considered as EHs and equipped with energy sources including renewables (WTs and PV systems), CHP units (convex and non-convex), and EVs. Due to the presence of uncertain parameters that greatly reduce the flexibility of the system, the downside risk constraint (DRC) method is applied to evaluate the flexibility of the system under risk-averse and risk-neutral modes. In comparison with the risk-neutral mode, the operator acts more conservatively in the risk-averse mode. For instance, as hub5 has a significant level of electrical demand, more power is imported from the network (4500 kW), and a small amount of value is sold to the grid (100 kW) based on the conservative behavior of the decision maker. Also, by increasing the system’s operational cost, the amount of risk is set at zero. In this concern, the operational cost to achieve a zero-risk condition is increased by nearly 36%. Because the CHP units produce power and heat simultaneously, thermal energy is exchanged among the EHs to meet their thermal loads. Also, the consideration of the load-shifting strategy in the first and third hubs resulted in a significant electrical load reduction in the risk-neutral mode (20 kW in hub1 and 15 kW in hub3). By and large, although the operational cost rises in the risk-averse mode, the decision maker becomes capable enough to face uncertain parameters.

In future research studies, the game-theoretic modeling of energy hubs could be examined based on free competition. Also, applying the ADMM algorithm to analyze the P2P heat and power transaction in a decentralized mode could be another research direction for studying the power flow among the energy hubs.