**Contents**


## **About the Editors**

**Javier Reneses** is Senior Research Associate Professor in the Institute for Research in Technology (IIT) at the Universidad Pontificia Comillas in Madrid. He was a visiting scholar at the Berkeley Lab from 2016 to 2019. He has headed more than 80 research and consultancy projects and has participated in more than 150. He has worked and lectured extensively on the operation, planning, and regulation of energy systems, and particularly on medium-term operation of electricity markets, regulation of electric distribution business, tariff design, and natural gas markets. He has served as a consultant on these topics for governments, international institutions, industrial associations, and utilities in numerous countries. Dr. Reneses has published more than 100 papers and book chapters published in national and international journals and conference proceedings. He also has been teaching Electric Distribution Business at the Universidad Pontificia Comillas for the master's degree in the Electric Power Industry as well as Statistics and Energy and Sustainability at the School of Engineering for several years. He obtained his Electrical Engineering degree in 1996 from the School of Engineering (ICAI) of the Universidad Pontificia de Comillas, Madrid, Spain. He received his doctoral degree in May 2004 for the thesis "Mid-Term Operation of Electricity Markets". He obtained his degree in Mathematics in 2005 from the Universidad Nacional de Educacion a Distancia ´ (UNED), Madrid.

**Antonio Bello** is Research Assistant Professor at the Institute for Research in Technology (IIT) of ICAI School of Engineering at Comillas Pontifical University. Antonio holds an Industrial Engineering degree, an MSc degree in Electric Power Systems, and a Ph.D. degree in Electrical Engineering from the Comillas Pontifical University-ICAI, with a dissertation about probabilistic forecasting of electricity prices, which he defended in 2016. Throughout his career, he has participated in over 70 research and consultancy projects dealing with operation, simulation models, planning of electricity and gas markets, energy forecasting, and risk management support. Besides his research activities, he teaches time series analysis and forecasting in the energy sector. He has published numerous papers in national and international journals and conference proceedings on these topics. Antonio Bello has been a visiting researcher at London Business School (London, United Kingdom).

## *Article* **Expansion Planning Method of the Industrial Park Integrated Energy System Considering Regret Aversion**

#### **Haokai Xie 1,\*, Pu Zhao 2, Xudong Ji 2, Qun Lin <sup>2</sup> and Lianguang Liu <sup>1</sup>**


Received: 3 October 2019; Accepted: 23 October 2019; Published: 27 October 2019

**Abstract:** Industrial parks have various sources and conversion forms of energy. The many uncertainties in the planning of industrial park integrated energy systems (IPIES) pose a great risk of regret in planning schemes; thus, an expansion planning method for an IPIES, considering regret aversion, is proposed. Based on comprehensive regret value consisting of min–max regret aversion and the min average regret value, the method optimizes the comprehensive cost of the expansion planning scheme in IPIES under different natural gas price fluctuation scenarios, including costs of construction, operation and maintenance, and environmental protection. A multi-stage expansion planning scheme and typical daily operation plans under multiple natural gas price fluctuation scenarios of the IPIES in an economic and technological development zone in southeast China are used to demonstrate the validity of the method. The results show that, compared with a traditional planning method based on expectation, the proposed expansion planning method could reduce the maximum regret value by 14% on average, and greatly reduces the risk of decision-making regret by up to 18%. At the same time, the influence of natural gas price on expansion planning of the IPIES is discussed.

**Keywords:** industrial park integrated energy system; expansion planning; natural gas price uncertainty; regret aversion; min–max regret value

#### **1. Introduction**

The integrated energy system (IES) integrates energy production, conversion, storage, and consumption [1,2]. It is an important development trend of energy technology to achieve coordinated and complementary optimization of multiple energy sources in the future [3–5]. Industrial parks with intensive demand for electricity, steam, cold, and heat energy are typical application scenarios for IES. How to plan industrial park integrated energy systems (IPIES) is an important issue in current research [6,7]. In order to optimize the structure and capacity of the IES, domestic and foreign scholars have proposed various models, algorithms, and planning objectives. Geidl et al. [8] proposed the concept of an energy hub (EH) and established a planning model for an electric power and natural gas system with the objective function of minimum energy loss in the EH. Zhang et al. [9] considered a variety of combined heat and power (CHP) generation units, designing CHP units with the objective of system economics and environmental performance. Zhang et al. [10] decomposed the planning model into two aspects: Investment and operation feasibility of a power system or natural gas system. Zhao et al. [11] proposed a three-level collaborative global optimization method for a combined cooling, heating, and power (CCHP) system. These research and planning methods were mainly designed for a system plan in a specific state.

Based on EH, Zhou et al. [12] proposed a collaborative expansion optimization configuration method for a renewable power system and a natural gas system. Considering the topological constraints of grid and gas networks, a regional integrated energy system expansion planning model based on CCHP was proposed in [13], and the economic scheduling strategy of the system was analyzed by the scenario method. Considering the value chain of natural gas, the long-term, multi-regional, and multi-stage expansion planning of a gas-electric coupling system was studied in reference [14]. The authors of [15,16] proved that flexible expansion planning can better handle uncertainties, compared to traditional lowest-cost planning methods. Adaptation cost is defined as the additional investments required for a proposed plan, if an unexpected load growth happens. Qiu et al. [17] established a joint expansion planning model for natural gas networks and power grids with the objective of maximizing social benefits. The Monte Carlo simulation was applied to create scenarios that simulate random system characteristics in [18]. An extended planning model based on two-stage stochastic optimization was established to realize the construction planning of natural gas and power facilities under uncertainty of demand growth in [19].To deal with the uncertainty of renewable energy, the scenario method was used to deal with the wind power and energy storage and load, and the capacity of the internal device in the EH was configured in [20]. Robust optimization was used to obtain an optimized solution that is immune to the effects of all possible wind power realizations within the uncertainty interval in [21,22] Cesena et al. [23] proposed a unified planning and scheduling method to assess the flexibility of system investment and operation under long-term uncertainty. The proposed approach in [24] reflected real-options thinking borrowed from finance, and had been cast as a stochastic mixed-integer linear program. These integrated energy expansion planning studies did not consider the need for policymakers to avoid the risk of regret in planning schemes together with the lowest cost.

Regret is an emotion that affects decision-making behavior. When decision-makers face a variety of schemes, one scheme is selected and compared with other unselected schemes. When the uncertainty causes the actual situation to be different from the expected, resulting in the profit of the selected scheme being smaller than one or more of the unselected schemes, a feeling of regret in the decision-maker is generated [24,25]. Savage et al. [26] proposed a model based on minimum-maximization regret, which shows that decision-makers will choose a decision with a minimal regret value from the decision plan that maximizes regret. In [27], regret is considered as one of the objectives in a multi-objective optimization framework. By applying the min-max regret criterion, model obtained a solution that minimizes the worst-case regret over all possible scenarios while ensuring system robustness [28]. The theory of regret has been applied in the study of consumer behavior in economics, travel path planning, etc. [29,30]. In [31], the min-max regret criterion is considered for the unit commitment problem. Compared with the min-max cost criterion, it is concluded that min-max regret outperforms min-max cost for certain unit commitment problems. Min-max cost and min-max regret have been proposed to address wind power generation uncertainties in [32]. Both criteria provide good upper bounds for the total costs under scenarios contained in an uncertainty set. The min-max cost criterion provides a smaller upper bound while min-max regret has higher variability. Depending on the characteristics of uncertainty sets and the preference of decision-makers, different models outperform each other under different situations. In the IPIES, the prediction error caused by uncertainty will cause decision-makers to pay excessive construction and operation costs for the system, which will also cause the decision-makers regret. In this paper, there are many types of comprehensive energy sources and conversion forms in industrial parks, such as natural gas, electric, heat, cold, and steam. In view of the various sources and conversion forms of energy in industrial parks as well as the large amount of uncertainty posing a great risk of regret in planning schemes, an expansion planning method of the IPIES considering regret aversion is proposed. The method deals with uncertainty using the scenario method, sets the regret value as the main indicator, and the comprehensive regret value, which considers the min-max regret value with its distribution and average regret value together, as the objective function.

#### **2. IPIES Structure and Expansion Planning Method Model**

Referring to [8,11,14,23,33,34], the structure and energy flow of the IPIES studied in this paper is shown in Figure 1. The system is connected to a steam network provided by an external large thermal power plant. The IPIES includes four parts: (1) Supply part: Power grid, natural gas network, steam network, photovoltaic (PV); (2) conversion part: Micro turbine (MT), heat recovery device (HR), electric boiler (EB), gas boiler (GB), heat pump (HP), electric chiller (EC), heat exchanger (HC), absorption chiller (AC); (3) storage part: Battery (BAT), steam heat storage (HS), cold energy storage (CS); and 4) load part: Electrical loads, steam loads, heat loads, cold loads, gas loads.

**Figure 1.** Structure and energy flow of industrial park integrated energy system (IPIES).

Based on the structure and energy flow of the IPIES, the expansion planning method of the IPIES considering regret aversion proposed in this paper is shown in Figure 2. It includes three layers: The stage scenario analysis layer, the expansion planning layer, and the regret aversion layer.

**Figure 2.** Expansion planning method of IPIES considering regret aversion.

In the stage scenarios analysis layer, the method mainly studies the typical daily load scenarios of the system and the natural gas price fluctuation scenarios.

In the expansion planning layer, the method sets the capacity and the typical daily operating power corresponding to device capacity as the variables. The method establishes an expansion planning model by taking the net present value of the comprehensive cost, including costs of construction, operation and maintenance, and environmental protection, as the objective function.

In the regret aversion layer, the method firstly calculates the optimal alternative schemes under different natural gas price fluctuation scenarios with the lowest comprehensive cost. Then, the method sets a min–max regret value and the lowest average regret value between the final planning scheme and the optimal alternative schemes as the objective function to optimize the device capacity within the final multi-stage expansion planning scheme and the typical daily operation plans under multiple natural gas price fluctuation scenarios.

#### *2.1. Stage Scenarios Analysis Layer*

Firstly, the scenario analysis method in [35] is used to deal with the volatility and randomness of each energy load and photovoltaic (PV) unit output power. The number of typical day scenarios after the load scenes and the lighting scenes are reduced is M, and each typical day scenario, *m*, has *Dm* days in a whole year. In order to meet the normal operation of the system at all times, an additional daily limit scenario with a constant zero PV output power is added. In view of the growth of energy load during the planning stages, the paper refers to the multi-stage planning method of a power system [36], introduces the continuous load curve to describe the medium and long-term load growth expectations, and divides the load level in the planning period into several horizontal sections, though the simplification will affect the accuracy of the model to some extent. A typical day scenario's load characteristic curves at different load levels can be obtained by equal ratio changes.

In addition, considering the price of the system device will decrease during the planning stages with the development of science and technology, and will stabilize after the technology matures [37], a piecewise exponential function is used to represent the dynamic change in device prices [38]:

$$\mathcal{L}^{J,y} = \begin{cases} \mathcal{L}\_0^{J,k} (1 + \mathcal{g}\_k)^y & 1 < y < \mathcal{Y}\_k^c \\\ \mathcal{L}\_0^{J,k} (1 + \mathcal{g}\_k^c) & \mathcal{Y}\_k^c < y < \mathcal{Y}\_k^c \end{cases} \tag{1}$$

where *cI*,*k*,*<sup>y</sup>* is the construction price of the device *k* in year *y*; *gk* is the price correction coefficient of the device *k*; *g<sup>c</sup> <sup>k</sup>* is the critical price reduction factor; *<sup>Y</sup><sup>c</sup> <sup>k</sup>* is the time for the device *k* to reach the critical price; and *Y* is the operating period of the IPIES.

The electricity market is still not perfect in China, but electricity prices are relatively stable due to policy decisions; natural gas prices, however, will be affected by changes in global trade prices and domestic supply and demand factors, and there will be greater uncertainty in prices over time and space. In the past, natural gas was originally developed as a replacement for traditional fuel, and its pricing was linked to other energy such as oil and the oil-indexed gas imports in China accounted for the majority [39–41]. However, with the changes in the international natural gas supply and demand pattern and the continuous reform of China oil and gas market, natural gas, especially liquefied natural gas (LNG), is gradually becoming an independent energy product. In 2018, China LNG imports accounted for 60% of total natural gas imports [42]. LNG breaks the restrictions on natural gas transmission and trade between regions, greatly enhancing the transmission and impact range of natural gas prices. The 2019 Wholesale Gas Price Survey shows that Henry-Hub priced US LNG exports continued rising and there is more gas price convergence amongst countries since the global gas market and market-related pricing [43]. Figure 3 shows the China LNG ex-factory price index given by the Shanghai Petroleum and Natural Gas Exchange, reflecting the price trend of LNG in the domestic market [42]. The Shanghai Petroleum and Natural Gas Exchange, which was officially launched in 2015, opened the situation that China natural gas prices are determined by competition between supply and demand. It can be seen from Figure 3 that China's LNG spot price also has large fluctuations in different periods.

**Figure 3.** China liquefied natural gas ex-factory price index given by the Shanghai Petroleum and Natural Gas Exchange.

Changes in natural gas prices will directly affect the operating costs of the IPIES, which in turn will affect the economic operation of the system and the expansion of the plan; the greater fluctuation, the higher the impact on the plan. Seasonal or yearly consideration of natural gas price fluctuations during the planning cycle will lead to excessive calculation of the model. The method uses stage scenario analysis techniques to analyze the uncertainty of natural gas prices during expansion planning. The method takes the natural gas price in the initial planning stage as the benchmark price, and considers that there may be fluctuations, η, in the price of the subsequent stage compared to the benchmark price; that is, it may rise, fall, or remain unchanged from the benchmark price. In the resulting natural gas price fluctuation scenario set, S, *a* planning stages corresponds to 3*a*−<sup>1</sup> natural gas price fluctuation scenarios, and the scenario probability corresponding to the natural gas price fluctuation scenario, *s*, is π*<sup>s</sup>* .

#### *2.2. Expansion Planning Layer*

#### 2.2.1. Objective Function

The expansion planning layer takes the net present value of the comprehensive cost, *CCOM <sup>s</sup>* , of multi-stage planning in a natural gas price fluctuation scenario, *s* ∈ S, as the objective function, including the system construction cost, *C<sup>I</sup> <sup>s</sup>*, operation cost, *CO <sup>s</sup>* , maintenance cost, *C<sup>M</sup> <sup>s</sup>* , and environmental protection cost, *CENV <sup>s</sup>* :

$$\mathbf{C}\_s^{\rm COM} = \mathbf{C}\_s^l + \mathbf{C}\_s^O + \mathbf{C}\_s^M + \mathbf{C}\_s^{ENV},\tag{2}$$

$$\mathbf{C}\_{s}^{I} = \sum\_{a} \frac{1}{\left(1 + \lambda\right)^{T(a-1)}} \cdot \left[\mathbf{c}\_{\mathbf{a}}^{\mathrm{I,k}} \cdot \left(\mathbf{W}\_{\mathbf{s,a}}^{\mathrm{k}} - \mathbf{W}\_{\mathbf{s,a-1}}^{\mathrm{k}}\right)^{\mathrm{T}} + \mathbf{c}\_{a}^{I,T} \cdot \left(l\_{a}^{\mathrm{des}} - l\_{a-1}^{\mathrm{des}}\right)\right] \tag{3}$$

$$\mathbf{C}\_{s}^{O} = \sum\_{y} \delta\_{y} \cdot \sum\_{m} D\_{m} \cdot \sum\_{t} \left( \mathbf{c}\_{t}^{E} \cdot P\_{s,a,m,t}^{SYS} + \mathbf{c}\_{s,a}^{G} \cdot G\_{s,a,m,t}^{SYS} + \mathbf{c}^{S} \cdot S\_{s,a,m,t}^{SYS} \right) \tag{4}$$

$$\begin{aligned} \mathbf{C}\_{s}^{O} = \sum\_{\mathbf{y}} \delta\_{\mathbf{y}} \cdot \sum\_{m} D\_{m} \cdot \sum\_{t} \Big( \mathbf{c}^{M,RAT} \cdot \Big| \begin{vmatrix} P\_{s,\mu,m,t}^{RAT} \\ s,\mu,m,t \end{vmatrix} + \mathbf{c}^{MHS} \cdot \Big| \mathbf{S}\_{s,\mu,m,t}^{HS} \Big| + \mathbf{c}^{M,CS} \cdot \Big| \mathbf{C}\_{s,\mu,m,t}^{CS} \Big| \\ + \mathbf{c}^{M,\mathbf{k}1} \cdot \mathbf{Q}\_{\mathbf{s},\mathbf{a},\mathbf{m},t}^{\mathbf{k}1} \Big), \end{aligned} \tag{5}$$

$$\mathbf{C}\_{s}^{ENV} = \sum\_{y} \delta\_{y} \cdot \sum\_{m} D\_{m} \cdot \sum\_{t} \left( \gamma^{E} \cdot P\_{s,a,m,t}^{SYS} + \gamma^{G} \cdot G\_{s,a,m,t}^{SYS} + \gamma^{S} \cdot S\_{s,a,m,t}^{SYS} \right) \tag{6}$$

$$\delta\_{\mathcal{Y}} = \frac{1}{(1+\lambda)^{y-1}},\tag{7}$$

where *a* is the planning stage, and each stage has *T* years; *y* is the operation year of IPIES; *t* is 24 intraday hours; δ*<sup>y</sup>* is the year discount rate; λ is the annual discount rate; **c I,k <sup>a</sup>** is the unit construction cost matrix of each device, including BAT, HS, CS, PV, CHP, GB, EB, HP, EC, HC, and AC; **W<sup>k</sup> s,a** is the total capacity matrix of each device in stage *a* under scenario *s*; c *I*,*T <sup>a</sup>* is the power transmission capacity expansion cost; *I des <sup>a</sup>* is the 0–1 mark for the power transmission expansion status, 1 after expansion, 0 before expansion; *c<sup>E</sup> <sup>t</sup>* is the electricity price at time *<sup>t</sup>*; *<sup>c</sup><sup>G</sup> <sup>a</sup>* is the price for a unit kWh of energy natural gas in stage *a* under scenario *s*; *c<sup>S</sup>* is the price for a unit kW·h of energy steam, and the unit kW·h energy price can be calculated by the unit cubic meter price or the unit steaming price and the low calorific value of the energy; *PSYS s*,*a*,*m*,*t* , *GSYS s*,*a*,*m*,*t* , and *SSYS <sup>s</sup>*,*a*,*m*,*<sup>t</sup>* are the electricity, gas, and steam power, respectively, that the system interacts with in the external network at typical day *m*, time *t* in stage *a* under scenario *s*; *cM*,*BAT*, *cM*,*HS*, and *cM*,*CS* are the unit maintenance cost of BAT, HS, and CS, respectively; **cM,k**<sup>1</sup> is the unit maintenance cost matrix of the device except for energy storage in IPIES; *PBAT s*,*a*,*m*,*t* , *SHS s*,*a*,*m*,*t* , and *CCS s*,*a*,*m*,*t* are the power exchange of the energy storage device for BAT, HS, and CS, respectively; **Qk**<sup>1</sup> **s,a,m,t** is the operating power matrix of each device except for energy storage in IPIES; and γ*E*, γ*G*, and γ*<sup>S</sup>* are the environmental cost of emissions from unit electricity, gas, and steam power, respectively, and can be calculated from the environmental value of the pollutants discharged per kW·h of energy.

#### 2.2.2. Constraints

The constraints in the expansion planning layer include expansion constraints, load part constraints, supply part constraints, conversion part constraints and storage part constraints.

#### (1) Expansion constraints

Each stage can expand the capacity of each device in the IPIES or maintain the configuration of the previous stage. Decommissioning is required for the life of the device to expire:

$$\mathcal{W}\_{s,a}^{k} \ge \mathcal{W}\_{s,a-1}^{k} - \mathcal{W}\_{s,a}^{k,out} \, \prime \tag{8}$$

$$I\_a^{\text{des}} \ge I\_{a-1}^{\text{des}} \tag{9}$$

$$\mathcal{W}\_{s,a}^{k,out} = \mathcal{W}\_{s,a-n\_k}^{k} - \mathcal{W}\_{s,a-n\_k-1'}^{k} \tag{10}$$

where *W<sup>k</sup> <sup>s</sup>*,*<sup>a</sup>* is the capacity of the device *<sup>k</sup>* in stage *<sup>a</sup>* under scenario *<sup>s</sup>*; *<sup>W</sup>k*,*out <sup>s</sup>*,*<sup>a</sup>* is the capacity of the device *k* to be decommissioned in stage *a* under scenario *s*; and *nk* is the number of planned stages that device *k* can serve; when *a* ≤ 0, *W<sup>k</sup> <sup>s</sup>*,*<sup>a</sup>* is 0.

(2) Load part constraints:

$$P\_{s\mu,m,t}^{PV} + P\_{s\mu,m,t}^{CHP} + P\_{s\mu,m,t}^{SYS} = P\_{s\mu,m,t}^{LD} + P\_{s\mu,m,t}^{EB} + P\_{s\mu,m,t}^{EC} + P\_{s\mu,m,t}^{EH} + P\_{s\mu,m,t}^{BAT} \tag{11}$$

$$S\_{s,a,m,t}^{CHP} + S\_{s,a,m,t}^{EB} + S\_{s,a,m,t}^{GB} + S\_{s,a,m,t}^{SYS} = S\_{s,a,m,t}^{LD} + S\_{s,a,m,t}^{HC} + S\_{s,a,m,t}^{AC} + S\_{s,a,m,t}^{HS} \tag{12}$$

$$H\_{s,a,m,t}^{HC} + H\_{s,a,m,t}^{EH} = H\_{s,a,m,t}^{LD} \tag{13}$$

$$\mathbf{C}\_{s,\mathfrak{a},m,t}^{AC} + \mathbf{C}\_{s,\mathfrak{a},m,t}^{\mathrm{EC}} = \mathbf{C}\_{s,\mathfrak{a},m,t}^{LD} + \mathbf{C}\_{s,\mathfrak{a},m,t'}^{\mathrm{CS}} \tag{14}$$

$$\mathbf{G}\_{s\rho,m,t}^{SYS} = \mathbf{G}\_{s\rho,m,t}^{LD} + \mathbf{G}\_{s\rho,m,t}^{CHP} + \mathbf{G}\_{s\rho,m,t'}^{CB} \tag{15}$$

where *PPV s*,*a*,*m*,*t* , *PCHP s*,*a*,*m*,*t* , *PLD s*,*a*,*m*,*t* , *PEB s*,*a*,*m*,*t* , *PEC s*,*a*,*m*,*t* , and *PEH <sup>s</sup>*,*a*,*m*,*<sup>t</sup>* are the output power of the PV, CHP, the electric load power, and the electric power consumed by the EB, EC, and EH, respectively; *SCHP s*,*a*,*m*,*t* , *SEB s*,*a*,*m*,*t* , *SGB s*,*a*,*m*,*t* , *SLD s*,*a*,*m*,*t* , *SHC s*,*a*,*m*,*t* , and *SAC <sup>s</sup>*,*a*,*m*,*<sup>t</sup>* are the output steam power of CHP, EB, GB, and the steam load power, the steam power consumed by the HC, and AC, respectively; *HHC s*,*a*,*m*,*t* , *HEH s*,*a*,*m*,*t* , and *HLD <sup>s</sup>*,*a*,*m*,*<sup>t</sup>* are the output heat power of the HC, EH, and the heat load power, respectively; *CAC s*,*a*,*m*,*t* , *CEC s*,*a*,*m*,*t* , and *CLD <sup>s</sup>*,*a*,*m*,*<sup>t</sup>* are the output cold power of the AC, EC, and the cold load power, respectively; *GLD s*,*a*,*m*,*t* , *GCHP s*,*a*,*m*,*t* , and *GGB <sup>s</sup>*,*a*,*m*,*<sup>t</sup>* are

the natural gas load power, the natural gas power consumed by CHP, and GB, respectively. Considering that the accuracy requirements of the planning are not as high as the actual running, in order to improve the efficiency of the model solving, the transmission loss of power grid, gas network and steam network were neglected.

#### (3) Supply part constraints:

$$P\_{\rm min}^{SYS} \le P\_{s,a,m,t}^{SYS} \le P\_{\rm max}^{SYS} + P\_0 \cdot I\_a^{des} \,. \tag{16}$$

$$\mathbf{G}\_{\text{min}}^{SYS} \le \mathbf{G}\_{s,a,m,t}^{SYS} \le \mathbf{G}\_{\text{max}}^{SYS} \tag{17}$$

$$S\_{\rm min}^{SYS} \le S\_{s,a,m,t}^{SYS} \le S\_{\rm max'}^{SYS} \tag{18}$$

where *PSYS* max, *GSYS* max, and *SSYS* max are the upper limit of the interaction power between the system and the external electricity, gas, and steam networks, respectively; *P*<sup>0</sup> is the capacity for power transmission expansion; *PSYS* min, *<sup>G</sup>SYS* min, and *SSYS* min are the lower limit of the interaction power between the system and the external electricity, gas, and steam networks, respectively;

#### (4) Conversion part constraints

In order to simplify the analysis, the operating efficiency of each energy conversion device is constant, and the variable operating characteristics are neglected. The constraints of GB, EB, AC, EC, HP, and HC are uniformly stated as:

$$Q^{k1,out}\_{s,a,m,t} = \eta^{k1} \cdot Q^{k1,in}\_{s,a,m,t'} \tag{19}$$

$$
\varepsilon\_{\text{min}}^{k1} \cdot \mathcal{W}\_{s,a}^{k1} \le Q\_{s,a,m,t}^{k1,in} \le \mathcal{W}\_{s,a}^{k1} \tag{20}
$$

where *Qk*,*in <sup>s</sup>*,*a*,*m*,*<sup>t</sup>* and *Qk*1,*out <sup>s</sup>*,*a*,*m*,*<sup>t</sup>* are the input and output power, respectively, of the above device, *k*1, at typical day *m*, time *t* in stage *a* under scenario *s*; *Wk*<sup>1</sup> *<sup>s</sup>*,*<sup>a</sup>* is the total configuration capacity of the device in stage *a* under scenario *s*; η*k*<sup>1</sup> is the operating efficiency of the device; and ε*k*<sup>1</sup> min is the lowest power factor of the device.

The CHP is coupled to both electricity and steam with the following constraints:

$$P\_{s,a,m,t}^{CHP} = G\_{s,a,m,t}^{CHP} \cdot \eta\_{\mathcal{R}^t}^{CHP},\tag{21}$$

$$S\_{s,a,m,t}^{CHP} = G\_{s,a,m,t}^{CHP} \cdot \eta\_{\text{gph}}^{CHP} \, , \tag{22}$$

$$
\varepsilon\_{\text{min}}^{\text{CHIP}} \cdot \mathcal{W}\_{s,\mu}^{\text{CHIP}} \le P\_{s,a,\text{wt},t}^{\text{CHIP}} \le \mathcal{W}\_{s,a}^{\text{CHIP}}.\tag{23}
$$

#### (5) Storage part constraints

The three types of energy storage device, including BAT, HS, and CS, have similar operating characteristics:

$$\mathcal{W}\_{s,\mu,m,t}^{k\_2} = \mathcal{W}\_{s,\mu,m,t-1}^{k\_2} \left( 1 - \mu\_{\text{loss}}^{k\_2} \right) + \left( \eta\_{ch}^{k\_2} \cdot \max \left( P\_{s,\mu,m,t}^{k\_2}, 0 \right) + \frac{\min \left( P\_{s,\mu,m,t}^{k\_2}, 0 \right)}{\eta\_{dis}^{k\_2}} \right) \cdot \Delta t,\tag{24}$$

$$
\rho\_{\text{min}}^{k\_2} \cdot \mathcal{W}\_{s,a}^{k\_2} \le \mathcal{W}\_{s,m,a,t}^{k\_2} \le \rho\_{\text{max}}^{k\_2} \cdot \mathcal{W}\_{s,a'}^{k\_2} \tag{25}
$$

$$1 - P\_{\text{max}}^{k\_2} \le P\_{s,m,a,t}^{k\_2} \le P\_{\text{max}}^{k\_2} \tag{26}$$

$$\mathcal{W}^{k\_2}\_{s,a,m,24} = \mathcal{W}^{k\_2}\_{s,a,m,0'} \tag{27}$$

where *k*<sup>2</sup> is the type of energy storage device in the IPIES; *Wk*<sup>2</sup> *<sup>s</sup>*,*a*,*m*,*<sup>t</sup>* is the stored energy of energy storage device *k*<sup>2</sup> at typical day *m*, time *t* in stage *a* under scenario *s*; μ *k*2 *loss* is the self-consumption rate of the energy storage device, *k*2; η *k*2 *ch* and η *k*2 *dis* are the charging efficiency and discharging efficiency of the energy storage device *k*2, respectively; Δ*t* is the unit scheduling time; ϕ *k*2 max and ϕ *k*2 min are the upper and lower limit coefficients, respectively, of the energy storage device, *k*2,energy stored; *Pk*<sup>2</sup> max is the upper limit of the switching power of energy storage device *k*2, related to the converter device. In order to achieve continuous scheduling, constrained energy storage stores the same energy at the beginning and end of the day.

#### *2.3. Regret Aversion Layer*

#### 2.3.1. Optimal Alternative Schemes

The regret aversion layer can be divided into two parts: Optimal alternative schemes calculation and regret aversion optimization. The optimal alternative schemes part uses the expansion planning layer model, with the comprehensive cost in Equation (2) being lowest as the objective function, and Equations (3)–(27) as the constraint. The comprehensive costs *CCOM <sup>s</sup>* under all natural gas price fluctuation scenarios, *s* ∈ S, are calculated as well as the corresponding optimal alternative schemes, ωs, and operational plans, τω*<sup>s</sup> <sup>s</sup>* . The calculated *CCOM <sup>s</sup>* will then be substituted as a known parameter into the evasive optimization part.

#### 2.3.2. Regret Aversion Optimization

In response to the regret resulting from the fact that decision-makers did not choose a better expansion planning scheme, the additional aggregate comprehensive cost was chosen as the regret value: In view of the regret aversion optimization part, the decision-maker selects the additional comprehensive cost as the regret value:

$$\mathcal{C}\_{\sf s}^{REG}(\omega, \tau\_{\sf s}^{\sf av}) = \mathcal{C}\_{\sf s}^{COM}(\omega, \tau\_{\sf s}^{\sf av}) - \mathcal{C}\_{\sf s}^{COM}(\omega\_{\sf s}, \tau\_{\sf s}^{\sf av})\_{\sf s} \tag{28}$$

where *CREG <sup>s</sup>* (ω, τω *<sup>s</sup>* ) is the regret value of the expansion planning scheme, ω,under scenario *s*; τω *<sup>s</sup>* is the typical daily operation plans based on scheme ω under scenario *s*; *CCOM <sup>s</sup>* (ω, τω *<sup>s</sup>* ) is the comprehensive cost of the scheme, ω, and operation plans, τ<sup>ω</sup> *<sup>s</sup>* , which is obtained during the expansion planning layer; and *CCOM <sup>s</sup>* (ωs, <sup>τ</sup>ω<sup>s</sup> *<sup>s</sup>* ) is the lowest comprehensive cost under scenario *s*, which is obtained in the optimal alternative schemes part.

We used the minimum–maximum regret value under all natural gas price fluctuation scenarios considering the distribution of the scenarios as an objective to control the regret risk of the IPIES expansion planning scheme, ω, under different natural gas price fluctuation scenarios:

$$\underset{s \in \mathbb{S}}{\text{minimize}} \pi \tau^s \cdot \mathbb{C}\_s^{\text{REG}}(\omega, \tau\_s^{\omega}) . \tag{29}$$

Further, the minimum–maximum regret aversion objective can be considered together with the minimum objective of the average comprehensive cost. As the optimal comprehensive cost in each scenario is known, the average regret value is equivalent to the objective of the minimum expected comprehensive cost. The objective function of the expansion planning method of the IPIES is finally constructed as:

$$\min \qquad \mathsf{C}^{\mathsf{CRE}}(\omega, \mathsf{r}\_s^{\omega}) \quad , \tag{30}$$

$$\mathbb{C}^{\rm CRE}(\omega, \mathsf{r}\_s^{\omega}) = a \cdot \max\_{s \in \mathbf{S}} \pi^s \cdot \mathbb{C}\_s^{\rm REG}(\omega, \mathsf{r}\_s^{\omega}) + \beta \cdot \sum\_{s \in \mathbf{S}} \pi^s \cdot \mathbb{C}\_s^{\rm REG}(\omega, \mathsf{r}\_s^{\omega}), \tag{31}$$

where *CCRE*(ω, τ<sup>ω</sup> *<sup>s</sup>* ) is the comprehensive regret value; α and β is the weight coefficient of the minimum–maximum regret aversion objective and the minimum average regret objective, and α+β = 1. When α is taken as 1, Equation (30) is equivalent to Equation (29).

The proposed method of regret aversion optimization uses Equation (30) as the objective function and Equations (2)–(28) as the constraint condition to optimize the multi-stage expansion planning scheme of IPIES and the typical daily operation plans under multiple natural gas price fluctuation scenarios.

#### *2.4. Method Solution*

The expansion planning method of IPIES considering regret aversion proposed in this paper is a mixed integer nonlinear programming model. Considering the variables and constraints of the model, mathematical modeling was performed on the MATLAB platform through the YALMIP toolbox, and the commercial optimization solver GUROBI was used to solve the model. The model solution environment for this article was: Intel Core 2 Duo P8600 CPU; 6 GB memory; software version: MATLAB R2014A; YALMIP R20180612; GUROBI 8.1.

#### **3. Case Study**

#### *3.1. Basic Data*

Taking an economic and technological development zone that includes more than 3000 industrial enterprises in Zhejiang province in southeastern China as the case study, the planned total operating life of the IPIES is 15 years, and every 5 years is an expansion planning stage. The typical days of the industrial park can be divided into summer, winter, and ordinary days. The number of days in the whole year is 100, 60, and 200 days. There is cold load demand in summer and heat load demand in winter. There is no gas load demand. The garment industry and beverage processing industry in the park have steam demand. The typical daily energy load curve and PV output curve of the park are shown in Figures 4 and 5, respectively. In addition, consider an extreme scenario where PV units do not output power in a typical summer day. According to the high growth forecast of the park, the sustained load of the three extended planning stages is 150%, 200%, and 250% of the current load. The current electric, steam, and natural gas prices in the park are shown in Table 1. Considering the continuous advancement of China oil and gas marketization reform, referring to Figure 3, it is assumed that the future natural gas price may fluctuate by up to 50% compared with the current price, which has uniform distribution characteristics, corresponding to nine natural gas price fluctuation scenarios with the same probability: (S1) 100%-150%-150%; (S2) 100%-150%-100%; (S3) 100%-150%-50%; (S4) 100%-100%-150%; (S5) 100%-100%-100%; (S6) 100%-100%-50%; (S7) 100%-50%-150%; (S8) 100%-50%-100%; and (S9) 100%-50%-50%.

**Figure 4.** Typical daily electricity, steam, cold, and heat load curves.

**Figure 5.** Typical daily 350-kW photovoltaic output curve.



The approximate IPIES device operating parameters, unit construction costs, and operation and maintenance costs are shown in Table 2 [44–47]. The relevant parameters of the energy storage device are shown in Table 3 [33,44,47]. The upper limits of the interaction power between the park and the grid and steam network are related to the network capacity, and they are 280 and 192.25 MW (250 T/h), respectively. The system does not output energy to the external network; that is, the lower limit of interaction between each network is 0. The minimum power factor of each device is 0. The power system in the park can be expanded to 35 MW, and the expansion cost is 5 million yuan. The environmental protection costs of electricity, natural gas, and steam using kWh energy per unit in the park are 0.07, 0.01, and 0.04 yuan/kWh, respectively [44,45]. The annual discount rate of the park is 5%, the low heat value of steam is 769 kWh/T, and the low heat value of natural gas is 9.9 kWh/m3.

**Table 2.** Device parameters of IPIES.



**Table 3.** Parameters of energy storage.

#### *3.2. Results and Analysis*

The optimized IPIES expansion planning scheme, which takes α as 0.5, is shown in Table 4. The regret value and cost of the scheme under different natural gas price fluctuation scenarios are shown in Tables 5 and 6, respectively.

**Table 4.** The IPIES expansion planning scheme considering regret aversion.


**Table 5.** Regrets under multiple natural gas price fluctuation scenarios.


**Table 6.** Costs under multiple natural gas price fluctuation scenarios.


As seen from the expansion planning scheme in Table 4, with an increase in the load level of the park, the capacity of most energy conversion devices in the system expands. Due to the dynamic changes in the price of the system device, and the unpredictable fluctuations in the system's natural gas price from the second phase, the device capacity increase between stages 1 and 2 is greater than that between stages 2 and 3.

The PV in the system uses renewable energy and has excellent economic benefits with the most installed capacity. The CHP is the main natural gas drive device in the IPIES. At the current natural gas price, the CHP has a certain economic advantage over the grid peak hour electricity price and normal electricity price. However, in the case of natural gas price fluctuations, this situation will change. Figure 6 shows the system energy purchases of stage 2 under the scenarios of different natural gas price fluctuations. It can be seen that, when the price of natural gas rises by 50%, the purchase of natural gas in the system decreases, and the purchase of electric and steam increases. When the price of natural gas drops by 50%, the purchase of natural gas in the system increases, and the purchase of electric and steam declined, with the purchase of steam almost falling to zero.

**Figure 6.** Stage 2 energy purchase situation under different scenarios.

As the system is connected to an external steam network, the role of the GB in supplementing the electro-thermal ratio is replaced to some extent, with fewer configurations in the system. The EB works only at the peak of the PV output, making full use of the electric energy generated by the PV and making up for the lack of steam energy caused by the reduction of the power of the CHP. Heat storage stores energy when the load demand is low and discharges it when the load demand is high. The BAT can store energy in the electricity price valley, and the energy can be released at the peak and normal time to achieve the peaking and filling of the load and the economic improvement of the system.

The absorption capacity of the AC and HC in the IPIES is matched with the CHP. It can be seen from Figure 7a that, in the normal output power of the PV, the CHP reduces the output power to make full use of the renewable energy, and the cooling power in the system is mainly provided by the EC. As the PV output power decreases after 14:00, the output power of the CHP increases, and the output cooling power of the AC increases. When entering the electricity price valley at 22:00, the system energy supply turns to the grid and the output power of the EC rises again. In the case of unexpected failure of the photovoltaic unit, the daily load needs to be supplemented by the CHP. At the same time, the output of the EC will rise when the electricity price valley is between 11:00 and 13:00, and the output of the CHP and the AC will decrease.

Comparing Figure 7a,b, it can be seen that in order to ensure the user's energy demand, the system often needs to configure other backup devices to prevent the renewable sources from fluctuating or even zero output, causing the load shedding. Thus, although photovoltaic units have high economic benefits, there are certain restrictions on the permeability of renewable energy in the system. It is necessary to consider the extreme output scenarios of some renewable energy units in the planning process to optimize the capacity of the units and other related units to ensure the safe and economic operation of the system.

Tables 5 and 6 show the regret value and cost of the plan in different natural gas price scenarios, respectively. The planning scheme has a large regret value under scenario 1 and scenario 9, and the regret value of scenario 5 is the smallest. Comparing Table 4 and optimal alternative schemes' typical device planning in stage 3 under scenarios 1, 5, and 9 in Table 7, it can be seen that the high or low natural gas price in scenario 1 and scenario 9 leads to the significant difference in the configuration capacity of the equipment between the schemes, and further causes the actual planning scheme not to match the optimal alternative, resulting in an increase in regret value.

**Figure 7.** Summer typical day partial devices output power in stage 3 under scenario 5: (**a**) Typical device output power; (**b**) Typical device output power when the PV output is 0.



The cost difference between the different scenarios of the planning scheme is mainly due to the difference in operating costs. The overall energy consumption of scenario 1 is high, and the running cost is high. The overall energy consumption of scenario 9 is low, and the running cost is low. Because the energy demand of phase 3 is high, the high energy cost of phase 3 has a greater impact on the overall operating cost than phase 2, and the cost of scenario 7 is higher than that of scenario 3.

#### **4. Discussion**

To verify the effectiveness of the method, we considered three planning methods for comparison:

Case1: Expansion planning method that considers regret aversion proposed in this paper;

Case2: Expansion planning method based on the lowest expected cost; and

Case3: Expansion planning method that does not consider gas price fluctuations.

The regret value of the schemes obtained by different planning methods is shown in Figure 8. The typical device expansion planning of each scheme is shown in Table 8. The regret values under multiple scenarios with different planning methods are shown in Table 9.

**Figure 8.** Regret value under different planning methods.


**Table 8.** Typical device expansion planning under different planning methods.

**Table 9.** Regret value under multiple scenarios with different planning methods.


More PV units are deployed in the expansion plan with the lowest expected cost, showing the economic benefits of renewable energy in the system. However, in case 2 and case 3, the plan of natural gas-related devices, such as combined heat and power (CHP), is insufficient, resulting in a large increase in regret when the price of natural gas is low. Compared to the extended planning method without consideration of the fluctuation of gas prices, the proposed method reduces the maximum regret value by 17.8% and reduces the comprehensive regret value by 9.4%.

Compared with the extended planning method based on the lowest expected cost, the method proposed in this paper has a lower regret value. The maximum regret value in case 1 is effectively constrained by the objective function. By introducing more natural gas equipment such as CHP, plan scheme in case 1 has better performance in scenario 6 to 9, especially in the worst scenario 9, but worse performance in scenario 1 to 5. The proposed method reduces the maximum regret value by 13.8% and reduces the comprehensive regret value by 2.7%. Although the average regret value increases by 0.92%, the reduction in the comprehensive regret value indicates that the benefit of controlling the maximum regret value exceeds the control of the average regret value under the decision-maker's risk control requirement.

Further considering the influence of the minimum maximum regret aversion weight coefficient, which represents the risk control requirement of the decision maker, the reduction of comprehensive regret value between case 1 and case 2 can be calculated as:

$$\frac{\mathbf{C}^{\rm CRE}\Big(\boldsymbol{\omega}\_{2\prime}\boldsymbol{\tau}\_{s}^{\omega\_{2}}\big) - \mathbf{C}^{\rm CRE}\Big(\boldsymbol{\omega}\_{1\prime}\boldsymbol{\tau}\_{s}^{\omega\_{1}}\big)}{\mathbf{C}^{\rm CRE}\Big(\boldsymbol{\omega}\_{2\prime}\boldsymbol{\tau}\_{s}^{\omega\_{2}}\big)},\tag{32}$$

where ω1, ω<sup>2</sup> is the plan scheme in case 1 and case 2, respectively

The comprehensive regret reduction between case 1 and case 2 under different minimum–maximum regret aversion objective weights, α, are shown in Figure 9. It can be seen from Figure 9 that when the range of is changed from 0.1 to 0.9, the comprehensive regret reduction rises, which indicates that with the increase of the decision-makers' requirement for maximum regret risk control, the planning method proposed is better than traditional method based on the lowest expected cost, making the plan more adaptive when faced with uncertain natural gas prices. If the decision makers have low demand for risk control, the planning method proposed is similar to the traditional method but still provides a little reduction in the comprehensive regret value. It shows that in the industrial park integrated energy system expansion plan, due consideration is given to the regret aversion factor, which can effectively control the regret risk of system planning decisions, and make the plan more adaptive when faced with uncertain natural gas prices.

**Figure 9.** Comprehensive regret reduction between case1 and case 2 under different minimum–maximum regret aversion objective weights, α.

However, the model is relatively simple while the transmission loss of the power grid, gas network, and steam network were entirely neglected in the paper. The theory of how the regret value and expansion plan is affected by load growth expectation was also not put forward in this paper.

#### **5. Conclusions**

This paper proposed an expansion planning method for the industrial park integrated energy systems considering regret aversion. Based on the min–max regret aversion and the lowest average regret value, the method optimized the comprehensive cost of an expansion planning scheme in an IPIES under different natural gas price fluctuation scenarios, including costs of construction, operation and maintenance, and environmental protection. The example verifies the rationality and effectiveness of the proposed method. The optimized industrial park integrated energy system expansion plan greatly reduces the degree of decision-making regret and reduces the system cost compared with the traditional expansion plan, which does not consider natural gas price fluctuation. Compared with the expansion plan based on the lowest expected cost, it also effectively controls the system's decision-making regret risk. At the same time, the simulation results show that natural gas price fluctuations have a greater impact on system planning and operation.

With the deepening of the national power system reform, multi-regional integrated energy system collaborative planning and multi-subject integrated energy system planning and operation game theory will be the focus of future research.

**Author Contributions:** H.X. conceived the main idea and wrote the manuscript with guidance from Q.L., and L.L., who reviewed the work and gave helpful improvement suggestions. P.Z. and X.J. managed the project and provided case data.

**Funding:** This work was supported by the Science and Technology Project of the State Grid Corporation of China (No. SGZJWZ00FZJS1901007).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**




#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Financial Impacts of Net-Metered Distributed PV on a Prototypical Western Utility's Shareholders and Ratepayers**

#### **Peter Cappers 1,\*, Andrew Satchwell 1, Will Gorman <sup>1</sup> and Javier Reneses <sup>2</sup>**


Received: 5 November 2019; Accepted: 13 December 2019; Published: 16 December 2019 -

**Abstract:** Distributed solar photovoltaic (DPV) under net-energy metering with volumetric retail electricity pricing has raised concerns among utilities and regulators about adverse financial impacts for shareholders and ratepayers. Using a pro forma financial model, we estimate the financial impacts of different DPV deployment levels on a prototypical Western U.S. investor-owned utility under a varied set of operating conditions that would be expected to affect the value of DPV. Our results show that the financial impacts on shareholders and ratepayers increase as the level of DPV deployment increases, though the magnitude is small even at high DPV penetration levels. Even rather dramatic changes in DPV value result in modest changes to shareholder and ratepayer impacts, but the impacts on the former are greater than the latter (in percentage terms). The range of financial impacts are driven by differences in the amount of incremental capital investment that is deferred, as well as the amount of incremental distribution operating expenses that are incurred. While many of the impacts appear relatively small (on a percentage basis), they demonstrate how the magnitude of impacts depend critically on utility physical, financial, and operating characteristics.

**Keywords:** distributed solar PV; financial analysis; net-energy metering; investor-owned utility; earnings; return on equity; retail rates; ratepayer bills

#### **1. Introduction**

Residential solar power is rapidly expanding in the United States (U.S.). In 2018, there was a 7% increase in residential distributed solar photovoltaic (DPV) deployment [1]. Such large increases in the deployment of DPV in the U.S. over the previous 5–7 years has been attributed to significant declines in equipment costs [2], state and federal tax credits, and electric utility net-energy metering (NEM) compensation programs [3]. NEM is a billing mechanism that credits customers with distributed generation systems for any electricity they export to the grid [4]. Use of NEM in conjunction with volumetric retail electricity pricing (i.e., uniform compensation of generation in excess of consumption, regardless of its characteristics such as time of generation), however, has also raised concerns among utilities and regulators of higher retail electric power rates and shifting of costs from DPV to non-DPV customers [5]. While current amounts of DPV in many jurisdictions are small and thus any retail rate and cost-shifting concerns may be anticipatory in nature, NEM reforms have been proposed and, in certain cases, adopted by state public utility commissions [6]. Importantly, most reforms change the DPV system payback periods and have the potential to reduce distributed solar PV deployment [7]. Barbose et al. [8], for example, modeled effects of a reduction in NEM compensation for grid exports

from retail to wholesale electric power rates. They found that this reduction in NEM compensation would reduce residential 2050 solar PV deployment by approximately 20%.

Electric investor-owned utilities (IOUs), particularly those in the United States, are concerned about the effects of DPV on future earnings opportunities from deferred or avoided capital investments under existing regulatory and business models [9]. IOUs increase their earnings base by investing in capital, which may be growth-related (e.g., new distribution system investments and generating plants to serve increasing load) [10]. Thus, stagnant or declining sales as a result of DPV, as well as energy efficiency [11] and other forms of other forms of distributed energy resources (DERs), may reduce these future earnings opportunities [12]. Future growth in electric vehicle penetration, among other sources of electrification which increase electricity consumption, may counter the prevailing trend of declining sales.

Furthermore, the decrease in DPV and other forms of DER costs (e.g., battery storage) has led to increased financial pressure on the utility from customer self-supply [13]. Many utilities around the world typically allocate a significant portion of their fixed costs to volumetric energy charges. As a result, any reduction in electricity sales without a similar reduction in fixed costs erodes a utility's net revenues. Such impacts are especially disconcerting to utility shareholders due to the reduction in achieved earnings and return on equity (ROE) [14].

At the same time, utility regulators are increasingly concerned about possible increases in retail rates and cost-shifting from customers with DPV (i.e., participants) to non-DPV customers (i.e., non-participants) [15]. In instances where costs increase faster than sales, there is upward pressure on retail rates. In addition, customers who invest in DPV and can significantly reduce or even eliminate the volumetric portion of their bills via NEM may not adequately pay their full share of the utility's fixed costs, which places an increased cost burden on non-participating customers.

Regulators of any utility in such a situation must weigh utility and ratepayer concerns as they consider changes to NEM and retail rate design that directly affects DPV, battery storage, and other forms of DER. Ultimately, they must make a determination that they believe serves the broader public interest based on the information they have available to them.

While these concerns are qualitatively understood, there is a lack of empirical and quantitative analysis to bound the magnitude of these concerns and the efficacy of alternative utility regulatory and business models. Instead, most quantitative analyses evaluating DPV impacts on an electric utility exclusively focus on a simplified cost and/or benefit analysis without considering the financial implications on a utility [16]. Those studies which focus on DPV costs find incremental electric system costs which range from \$0–25/MWh [17–24]. On the other hand, those studies which focus on DPV benefits find electric system benefits which range from \$0–53/MWh [25–32]. In general, these analyses focus on system costs and benefits without considering the role of a regulated utility and their accompanying business practices.

Another subset of the literature evaluates how NEM and DPV adoption impacts utility ratepayers but notably does not incorporate a fully-integrated pro forma financial model. Poulilkkas [33] studied the effect of NEM on DPV adoption at one representative household while Christoforidis et al. [34] performed a similar analysis across a broader set of 31 customers. Neither study evaluated the impact on electric utility collected revenues or earnings. Eid et al. [35] and Picciariello et al. [36] expanded on this work by evaluating the cross-subsidization due to NEM between PV owners and non-owners. Furthermore, Johnson et al. [37] added an analysis of the DPV impact on utility load shape into the analysis of NEM cross-subsidization between PV owners and non-owners. However, none of the above quantitative studies: (1) calculate utility shareholder impacts, (2) take into account regulatory lag and other artifacts of utility regulation (e.g., test years), or (3) integrate a long-term analysis horizon that incorporates feedback effects between PV hourly loads and utility costs that accumulate over time to impact electricity rates.

Many of the concerns expressed by utilities and regulators, though, depend critically on the specific changes in costs and revenues that are a function of utility cost compositions (e.g., proportion of fixed versus variable costs), physical system attributes (e.g., hourly loads), and fixed cost recovery mechanisms, among others, which interact with PV adoption and utility regulatory considerations over time. Most of the prior studies reviewed above do not perform a comprehensive financial analysis using these key utility characteristics and do not incorporate a robust review of the costs and benefits of DPV in retail electricity rates.

To fill this research gap, this study quantifies the financial impacts of net-metered DPV on a prototypical Western U.S. IOU and identifies the key sensitivities driving lesser or greater magnitude of impacts. While we integrate the above cost-benefit studies into our financial modeling, we use a financial framework that better assesses the implications for a regulated IOU. Furthermore, we build on prior quantitative analysis of the financial impacts of net-metered PV [38,39] by assessing a wider range of sensitivities specific to the ability of DPV to avoid or defer utility costs (i.e., "DPV value"). Although the costs, revenue, and regulatory accounting assumptions are based on the U.S. context, we believe the results are generalizable and relevant for other utility circumstances around the world.

#### **2. Materials and Methods**

We quantify the shareholder and ratepayer impacts for a Western U.S., vertically-integrated IOU (i.e., that owns generation, transmission, and distribution assets) at two DPV deployment levels (i.e., 4%, and 8% of 2027 retail sales) representing the range of forecasted DPV deployment among Western states [40] using a proprietary pro forma financial model.

The FINancial impacts of Distributed Energy Resources (FINDER) pro forma financial model quantifies the utility's annual costs and revenues over a pre-defined analysis period. See Appendix A for a more detailed description of the pro forma financial model used in this analysis. The model performs all cost calculations at the total utility level but has the ability to allocate those costs to different rate classes in order to assess differential revenue impacts. Key outputs include achieved ROE and earnings, average all-in retail rates and customer bills, which can be used by utilities, policymakers, customer groups, and other stakeholders when assessing the impacts and implications of policy proposals and decisions.

Because the model derives customer class level ratepayer metrics, a more comprehensive bill impact analysis across different customer types, which assumes different hourly consumption profiles for customers and different DPV production profiles for participants, was not possible. A number of other studies have sought to quantify participant and/or non-participant bill impacts in substantial detail, but without any associated feedback effects on rates [41–43]. However, our estimates of the percentage change in average all-in retail rates can be used as a proxy for the percentage change in a non-participants' bill assuming there are no associated changes in their electricity consumption.

Results of our analysis using this pro forma financial model are compared against a case without DPV, incremental energy efficiency, or other DERs in order to isolate the DPV impacts. The DPV is installed linearly over ten years to reach its terminal deployment level (as a percentage of retail electricity sales) with impacts measured over 20 years to capture utility system economic end-effects (i.e., cost avoidance or deferral). We limit our analysis to 20 years despite DPV system lifetimes in excess of 20 years due to the effects of discounting costs and revenues, in addition to increasing uncertainty in utility cost and load forecasts. We also assess the sensitivity of impacts to different assumptions about the "DPV value" (i.e., ability of DPV to avoid or defer utility costs). Each of these different cases are discussed in more detail below.

#### *2.1. Prototypical Western Utility Characterization—Base Case*

We developed a 20-year cost and load forecast for a prototypical Western utility by starting with the 2014 general rate case filing of Public Service of Colorado. The general rate case filing was the most recent for the utility that included a cost-of-service study and provided reasonable starting year cost levels, starting year class-level retail sales, peak demand, and customer counts, and class-level cost allocation and rate design. Growth in retail sales, peak demand, and customers are based on Public Service of Colorado's 2016 integrated resource plan (IRP), which was the most recent one available. Growth in utility cost categories, specifically generation capital expenditures (CapEx), distribution CapEx, and operations and maintenance (O & M), are based on historical 5-year average annual growth rates among Western utilities derived from their FERC Form 1 filings. Last, the Western utility's hourly load shape is based on Public Service of Colorado's 2017 hourly load (reported in EIA Form 930) and we used a simple average of load in hours before and after missing values to derive a complete 8760-h load shape. Importantly, while many of the input assumptions were seeded with a single utility's data, the utility in this analysis is not meant to represent Public Service of Colorado specifically.

We refer to four key assumptions in the Western utility characterization when describing financial impacts. First, non-fuel costs (inclusive return of (i.e., depreciation) and on capital investments, fixed O & M, interest expense, and taxes) grow by 3.3% per year over the 20-year analysis period (2018–2037) (see Figure 1). Western utilities have seen median generation, transmission, and distribution CapEx costs increase by 6.4%, 3.6%, and 3.8% per year from 2012 to 2016, respectively (calculated based on utility FERC Form 1 data), and we assume similar, rounded CapEx cost escalations (i.e., 6%, 4%, and 4% annual growth in generation, transmission and distribution CapEx costs, respectively). Utility fuel and purchased power (energy and capacity) costs grow by 2.6% per year over the same 20-year analysis period. Our fuel and purchased power costs for non-renewable generation technologies are based on EIA fuel, heat rate, and variable O & M cost assumptions [44]. Wind and solar PPA costs are based on NREL's Annual Technology Baseline LCOEs in the "Wind TRG 4" and "Solar CF 20%" forecasts [45].

**Figure 1.** Forecasted Western utility costs (without DPV).

Second, the Western utility's retail sales grow by about 1.0% per year and annual peak demand grows by about 0.9% per year. Our load growth assumptions are based on a specific utility's IRP forecasts in order to match any incremental generation or purchased power, though the retail sales forecasts are higher than historical, median sales growth among Western utilities. From 2012 to 2016, Western utility median sales slightly declined by about 0.3% per year based on EIA-861 data.

Third, we assume no incremental generation additions in the Base case analysis, as Public Service of Colorado is forecasting PPAs to meet incremental load in its 2016 IRP. We make this assumption about no incremental generation additions in order to maintain consistency between our load and cost assumptions. Given that Western utilities averaged flat, or declining, load growth over the last 5-years, we believe our assumption is reasonable. Nevertheless, we consider the case of incremental generation additions in our sensitivity analysis. We also assume retirement of three generating units during the 20-year analysis period to maintain consistency with Public Service of Colorado's IRP loads and resources table. These are all input assumptions used to develop a pro forma revenue requirement and were not derived as part of a system planning module within FINDER. As such, the impacts of DPV

on utility capital costs are based on a coarser set of assumptions than might be possible with planning models representing the utility's generation dispatch and power flows on distribution feeders.

Fourth, we assume a flat retail rate design for all customer classes (as opposed to inclining block or time-of-use rates). Residential customer rates and bills collect 90% of revenues via volumetric energy rates with the remainder of revenues (10%) collected via a monthly, fixed customer charge. Commercial and industrial (C & I) customer rates and bills collect about 40% of revenues via an energy charge, 55% of revenues via a demand charge (based on class monthly non-coincident peak), and the remaining via a fixed, monthly customer charge.

#### *2.2. Alternative Assumptions in Utility Characterization—Sensitivity Cases*

We developed a set of cases to better understand the sensitivity of shareholder and ratepayer impacts from DPV to assumptions related to its capacity value and avoided generation, transmission, and distribution costs. These sensitivity cases involved modifying a number of parameters from the Base case (see Figure 2), based on ranges that exist in either the literature or Western utility historical data.

**Figure 2.** Modeled sensitivity cases among three key assumptions related to DPV value.

In our first sensitivity case, we change the assumed percent of transmission and distribution (T & D) CapEx that are growth-related. As previously discussed, we model two categories ofT&D CapEx: load growth-related and non-load growth-related. The Base Case assumes a portion (33%) of investments are growth-related CapEx and the addition of DPV reduces this growth-related CapEx proportional to reductions in annual peak demand. The 33% assumption is consistent with assumptions in prior studies [38,39,46]. This results in corresponding reductions in returns on ratebase, depreciation expenses, interest, and taxes For the sensitivity case, we bound the assumption with a low value of 10% growth-relatedT&D CapEx (i.e., Sensitivity Case 1.A.) and high value of 90% growth-related T&D CapEx (i.e., Sensitivity Case 1.B.). Appendix B shows the sensitivity case analysis results as change in DPV value, change in achieved earnings, and change in average all-in retail rates.

In our second sensitivity case, we change assumptions related to incremental CapEx. In one case, we increase distribution CapEx and O & M costs in conjunction with DPV to represent the possibility that integration costs for DPV could result in a net increase in distribution system expenditures (i.e., Sensitivity Case 2.A.). For the purposes of our study, system integration costs are borne by the utility (and all ratepayers via retail rates) and are different from interconnection costs that are paid by DPV owners. We assume incremental distribution O&M cost of \$15/kW-year installed DPV [47]

and incremental distribution CapEx of \$30/kW installed DPV [48]. Alternatively, we add incremental utility generation to meet future capacity needs motivated by the fact that some Western utility DPV value studies assumed deferred generation (i.e., Sensitivity Case 2.B.). As discussed in Section 2.1, the Base case utility characterization assumes the Western utility meets future capacity and energy through PPAs and short-term market purchases (as is consistent with the IRP data we used as the basis for our Western utility characterization). We base this sensitivity case on a simple loads and resources table and add mid-merit and natural gas generating plants in 100MW and 250MW increments to meet forecasted peak demand in the Base case without DPV. Capital and O & M costs of the incremental generation are based on EIA overnight capital cost estimates for generating plants and inflated at 2% per year.

In our third sensitivity case, we change the amount of DPV coincident with the Western utility's annual peak demand (i.e., capacity contribution to peak). The Base case assumes a 22% DPV capacity contribution to peak (discussed in the next section). We bound this assumption with a lower value of 12% (i.e., Sensitivity Case 3.A.) and higher value of 32% (i.e., Sensitivity Case 3.B.). The DPV capacity contribution to peak drives the reduction in annual system peak demand upon which capacity and T & D CapEx costs are based. Thus, an increase in DPV capacity contribution to peak would result in greater avoided capacity-driven costs at the same DPV deployment level.

#### *2.3. DPV Characterization*

Our pro forma financial model derives the impacts of DPV through several key static and dynamic interrelationships. DPV impacts utility billing determinants; specifically retail sales and peak demand, which has an effect on utility costs and subsequently retail rates. DPV reduces volumetric sales based on a direct relationship between the assumed annual DPV penetration, expressed as a percentage of annual sales on a customer-class basis, and the utility's class-level sales. The model derives reductions in the utility's peak demand through dynamic relationships of several variables that take into account: (1) the specific timing of DPV relative to the utility's hourly load, and (2) potential differences in alignment between when the DPV causes reductions in the utility's load and the utility system annual peak demand.

The timing of DPV production (savings) and the utility annual peak demand is a key driver in the analysis. DPV reduces energy only in hours when the PV system operates (i.e., during the daylight hours) and may also reduce utility system demand depending on whether the reduction in energy is coincident with the utility's system peak. Thus, the timing of DPV energy and demand savings in relation to customer class and utility peak demands (monthly and annual) drives projections of future costs, retail rates, and revenues collected on a volumetric basis through energy and demand charges.

To calculate the DPV shape, we relied on the National Renewable Energy Laboratory's (NREL) System Advisor Model (SAM) [49]. We simulated five solar profiles with typical meteorological year (TMY) weather data for different locations in Colorado's main metropolitan areas (i.e., Boulder, Aurora, Denver, Colorado Springs, and Pueblo). These simulations relied on PV Watts default assumptions (i.e., azimuth of 180 degrees, DC to AC ratio of 1.2, and tilt of 40 degrees). To estimate a single input profile for our pro forma financial model, we applied a population-weighted average solar output of the five metro area's solar shapes.

We further analyzed DPV's capacity contribution to peak load reduction by simulating DPV profiles using 2017 weather data. While the TMY weather data described above provides an ideal average profile, it does not provide an understanding of DPV's contribution to peak load reduction for our underlying load year of 2017. To determine this value, we simulated the additional solar profiles and sampled the capacity factor of our DPV simulation in the top-20 load hours of 2017 for Public Service of Colorado. The resulting capacity contribution to the top-20 load hours was 22%. DPV capacity factors are often calculated based on probabilistic simulation and modeling methods such as Effective Load Carrying Capability (ELCC), however, carrying out such a simulation is not the focus of this study; alternatively, we quantify DPV capacity contribution as the percentage of

the nameplate capacity that is available during top-20 load hours; this factor provides a simple and intuitive measure of the value provided by DPV in terms of capacity and can be represented in our pro forma financial model. This value became our Base Case capacity contribution to peak load reduction for DPV. We performed this analysis for a number of other Western utilities and found that this capacity value in the top-20 load hours ranged from 7% to 26%, which we use to inform our sensitivity cases discussed above.

The DPV portfolios reduce the Western utility's annual retail sales and peak demand. Retail sales grow by 1.0% per year in the case without PV, but the annual growth rate declines to 0.6%, and 0.2% in the 4% and 8% DPV deployment cases, respectively, from 2018 to 2027 as DPV systems are installed. Because the DPV penetration levels are specified in terms of a percent reduction of retail sales, they each reduce the Western utility's sales on a one-for-one basis. As shown in Figure 3, the reduction in retail sales increases as the DPV deployment level increases.

**Figure 3.** Forecasted Western utility annual retail sales without DPV and at increasing DPV deployment levels (4% and 8% of 2027 retail sales).

The impacts of DPV on the Western utility's annual peak demand depends on the timing and coincidence of DPV relative to the utility's annual peak demand. Figure 4 shows the forecasted annual peak demand for the Western utility from 2018 to 2027 and the coincident peak demand savings attributable to the DPV deployment cases. The prototypical Western utility modeled in this study has peak loads that occur in July generally between 2 p.m. and 6 p.m. prior to the addition of DPV systems. The coincident peak demand impact of DPV in our study is less than the retail sales impacts on a percentage basis (e.g., 0.8% per year reduction in retail sales and 0.6% per year reduction in peak demand in the 8% DPV deployment case) because the timing of maximum PV output does not coincide perfectly with the utility's annual peak demand. This is particularly the case for Public Service of Colorado that serves load near the Rocky Mountains, which results in lower DPV production in afternoon hours relative to other geographic locations due to the effect of mountain shadows [50]. We consider lower and higher contribution of DPV savings to peak in the DPV value sensitivities.

**Figure 4.** Forecasted Western utility annual peak demand without DPV and at increasing DPV deployment levels (4% and 8% of 2027 retail sales).

#### **3. Results**

#### *3.1. Sensitivity of DPV Value to Alternative Utility Assumptions*

Figure 5 shows the change in DPV value with Base and alternate assumptions of the proportion of T & D CapEx that is growth-related (i.e., Sensitivity Case 1.A. and 1.B.). A lower proportion of growth-related T & D CapEx results in a lower DPV value, and vice-versa, where the change in value occurs exclusively among non-fuel cost categories. In the 4% DPV deployment case, the DPV value ranges from \$51/MWh to \$57/MWh and, in the 8% DPV deployment case, the DPV value ranges from \$50/MWh to \$55/MWh. The modeled DPV value results are not particularly sensitive to this assumption, ranging from −4% to 6% relative to the Base Case assumption.

**Figure 5.** Sensitivity of DPV value to assumed proportion of growth-related CapEx.

Figure 6 shows the change in DPV value assuming incremental distribution and generation CapEx and distribution O & M costs to the Base Case (i.e., Sensitivity Case 2.A. and 2.B.). The incremental distribution CapEx and O & M costs (i.e., Sensitivity Case 2.A.) do not result in additional value, as the costs are added incrementally with the installed DPV and counteract many of the avoided costs. Thus, the DPV value declines significantly in both the 4% DPV and 8% DPV deployment cases. In fact, DPV value at 8% deployment and assuming incremental distribution CapEx and O & M is roughly half of the Base Case DPV value (\$27/MWh compared to \$52/MWh). The incremental utility generation assumption (i.e., Sensitivity Case 2.B.) reduces the avoided purchased capacity value, as would be expected due to the incremental generation installed in lieu of the capacity purchases. Also, as to be expected, the proportion of CapEx deferral value increases as the DPV defers or avoids some of the incremental generation. In total, however, the DPV value in the incremental utility generation case is about 10% lower than the Base Case assumption because the cost of utility-owned generation is lower relative to meeting the same capacity needs through PPAs and short-term market purchases. Thus, the incremental utility generation sensitivity case produces a lower total DPV value.

**Figure 6.** Sensitivity of DPV value to incremental non-fuel costs.

Figure 7 shows the change in DPV value with lower or higher assumed DPV capacity contribution to peak relative to the Base Case assumptions (i.e., Sensitivity Cases 3.A. and 3.B.). As expected, a lower capacity contribution to peak (i.e., Sensitivity Cases 3.A.) results in lower DPV value, and vice-versa, with the largest change in avoided capacity purchases. In the 4% DPV deployment case, the DPV value ranges from \$42/MWh to \$59/MWh and, in the 8% DPV deployment case, the DPV value ranges from \$40/MWh to \$54/MWh. The modeled DPV value results are quite sensitive to this assumption, ranging from −27% to 11% relative to the Base Case assumption. The results for all sensitivity cases show that DPV value is sensitive to alternate assumptions, but the degree depends on the specific assumption. For example, the assumed proportion of growth-related CapEx (i.e., Sensitivity Case 1.A. and 1.B.) has a small range whereas the DPV capacity contribution to peak (i.e., Sensitivity Case 3.A. and 3.B.) exhibits a much larger range of results.

**Figure 7.** Sensitivity of DPV value to assumed DPV capacity contribution to peak.

Importantly, we did not combine changes in more than one key assumption, which would likely drive greater change in DPV value than is observed in isolated cases (e.g., combine higher DPV capacity contribution to peak with higher CapEx deferral which would likely leader to greater DPV value).

#### *3.2. Sensitivity of Shareholder and Ratepayer Metrics to Alternative Utility Assumptions*

As shown in Figures 8 and 9, the impacts of DPV on shareholder earnings and ROE vary under these different assumptions related to the penetration of DPV and the value of DPV to the utility. In the Base Case, the after-tax earnings achieved by the Western utility decline as the DPV deployment level increases (1.5% reduction at 4% DPV and 3.1% reduction at 8% DPV) while its achieved average ROE declines as the DPV deployment level increases (1.6% reduction at 4% DPV and 3.2% reduction at 8% DPV) compared to the case without DPV. Adding incremental distribution CapEx and distribution O&M costs (i.e., Scenario 2.A.), or substituting PPAs with utility generation (i.e., Scenario 2.B.) alters the utility's non-fuel cost assumptions directly and, therefore, produce the most significant change in shareholder impacts. Across the range of sensitivity cases at 8% DPV, earnings impacts range from a 3.0% reduction (i.e., Scenario 1.B) to a 4.8% reduction (i.e., Scenario 2.A.), and ROE impacts range from a 2.8% reduction (i.e., Scenario 2.B) to a 5.4% reduction (i.e., Scenario 2.A) compared to the case without DPV (on a percentage, not absolute, basis).

Importantly, these percentage reductions are against small reductions in earnings and ROE on an absolute basis. For example, achieved earnings decline by \$123M (20-year present value) in the 8% DPV Base case out of a total earnings basis of \$3.96B (20-year present value). Even with DPV value assumptions driving the most impactful change in earnings that assume incremental distribution CapEx and O & M (i.e., Scenario 1.A.), the absolute reduction in earnings is \$190M (20-year present value). The same is true for achieved average ROE impacts that are a 25 basis-point reduction at 8% DPV in the Base Case. The reduction in achieved ROE assuming incremental distribution CapEx and O & M (i.e., Scenario 1.A.) is 42 basis points.

4% DPV 8% DPV

**Figure 9.** Sensitivity of Western utility achieved average ROE to assumptions related to DPV value.

As shown in Figures 10 and 11, the impacts of DPV on average all-in retail rates and customer bills vary under these different assumptions related to the penetration of DPV as well as the value of DPV. In the Base Case, the Western utility's all-in average retail rate increases as the DPV deployment level increases (1.1% increase at 4% DPV and a 2.4% increase at 8% DPV); however the reduction in sales associated with DPV exceeds the impact from the rate increase, resulting in a decline in total customer bills (1.8% reduction at 4% DPV and a 3.6% reduction at 8% DPV). Relative to the case without DPV, assumptions driving greater cost deferral (regardless of fuel and purchased power costs or CapEx-related costs), result in lower ratepayer impacts (i.e., lower average rate impacts and greater total customer bill reductions). Across the range of sensitivity cases at 8% DPV, average retail rate impacts range from a 2.3% increase (i.e., Scenario 1.B.) to a 3.4% increase (i.e., Scenario 2.A.) and total customer bill impacts range from a 2.6% reduction (i.e., Scenario 2.A.) to a 3.7% reduction (i.e., Scenario 1.B.) compared to the Base case without DPV. Importantly, these bill savings reflect the aggregate impact across all customers and do not reflect the distribution of customer bill impacts among participating and non-participating customers. However, for non-participants who are not assumed to change their electricity consumption across scenarios, the reported percentage changes in average all-in retail rates is a proxy for their associated bill impacts.

Like the shareholder impacts, though, the percentage reductions are against small changes in ratepayer metrics on an absolute basis. Specifically, average all-in retail rates increase by 0.22 cents/kWh in 8% DPV Base Case and by 0.32 cents/kWh with DPV value assumptions with the most impactful change in average all-in retail rates (i.e., Scenario 2.A.). These changes compared to an average all-in retail rate of 9.20 cents/kWh (20-year present value) without DPV in the Base case. Similarly, there is a \$1.31B decrease in *total* customer bills (out of ~\$36B basis) at 8% DPV. The total customer bill savings in the high CapEx deferral DPV value sensitivity (i.e., Scenario 1.B.) are \$1.35B at 8% DPV.

**Figure 10.** Sensitivity of Western utility average all-in retail rate to assumptions related to DPV value.

4% DPV 8% DPV

**Figure 11.** Sensitivity of Western utility total customer bills to assumptions related to DPV value.

#### **4. Discussion and Conclusions**

This analysis quantified the financial impacts of different levels of net-metered DPV deployment on a prototypical Western U.S. utility over 20 years and estimated changes in the utility's costs, revenues, achieved shareholder earnings and ROE, average all-in retail rates, and customer bills. We also quantified the sensitivity of results to different assumptions about the ability of DPV to avoid, defer, or increase utility fixed and variable costs.

Our analysis shows that DPV does in fact reduce utility achieved earnings, which occurs through two separate means. First, if DPV reduces utility revenues more than utility costs (as is likely under rate structures that recover the majority of utility costs via volumetric energy charges), then net revenues or earnings are likewise reduced (i.e., the "revenue erosion effect"). Second, and separately, DPV savings may also diminish future earnings opportunities by reducing the rate of growth or deferring capital investments (T & D CapEx in our Base Case assumptions, specifically) that would otherwise contribute to the utility's ratebase (i.e., the "lost earnings opportunity effect"). Although our analysis does illustrate that the financial impacts on shareholders increase as the level of DPV deployment increases, the magnitude is small even at high DPV penetration levels (e.g., 2 to 4% change in financial metrics at 8% DPV deployment).

Our analysis also shows that net-metered DPV increases average all-in retail rates, which occurs in two, interrelated ways. First, DPV affects the retail rates set within each general rate case (GRC) through the net result of reductions in the test-year utility costs and billing determinants used to establish rates. DPV generally tends to reduce utility sales more than costs and, as a result, average retail rates established through each GRC increase with the addition of DPV in order to ensure the utility's rates collect sufficient revenue to recover authorized costs. Second, DPV affects average all-in retail rates in the years between GRCs, though this effect is simply a mathematical artifact. Average all-in rates are, by definition, equal to total collected revenues divided by total retail sales in any given year. Retail sales (i.e., the denominator) are reduced due to the incremental DPV. Because

of these lower volumetric energy billing determinants, the revenues (i.e., the numerator) collected on an annual basis will also be reduced. However, the reduction in revenues are necessarily smaller than the reduction in retail sales, given that some portion of revenues are derived from fixed customer charges (which are unaffected by DPV) and demand charges (which are only marginally affected by DPV). Thus, DPV tends to increase average all-in retail rates in between GRCs as well. As with shareholder metrics, our analysis illustrates that ratepayer financial impacts increase as the level of DPV deployment increases, but the magnitude is small even at 8% DPV penetration levels.

We know that utilities in the Western U.S. are varied, and exist within dramatically different operating environments, face substantially different cost structures, and provide service to vastly different customer bases. To better understand the likely impacts of DPV on utilities in the West, and potentially extend the application of the results more broadly to other utilities around the world, this study also explicitly links different estimates of DPV value to shareholder and ratepayer impacts. Our analysis finds that even rather dramatic changes in DPV value result in modest changes to shareholder and ratepayer impacts. Also, the range of financial impacts under alternative DPV value assumptions are greater for shareholders than ratepayers on a percentage basis and driven by differences in the amount of incremental CapEx that is deferred, as well as the amount of incremental distribution O&M that is incurred. The sensitivity cases reflect the key drivers of our results, but are not a complete list of all the sources of uncertainty and variation in modeled assumptions. There are other utility characteristics that might also change the magnitude and, in more extreme instances, the direction of impacts (e.g., higher or lower assumed load growth, higher or lower proportion of revenues from fixed charges, current or future test years). See [46] for the results of a number of sensitivity cases beyond the value of DPV. To be sure, the shareholder and ratepayer impacts presume no change in the underlying regulatory model or ratemaking approaches. More fundamental changes in the way electric utilities price energy services and recover fixed and variable costs may suggest different impacts than reported herein (e.g., see [51]). As such, our purpose here is not to bound the full range of impacts, but rather to illustrate some key themes and considerations related specifically to DPV value.

It is worth noting two particular feedback effects that our pro forma financial model does not account for in the present study and that would have potential implications for our results. First, we do not represent the feedback effects between retail rate impacts and DPV adoption rates. An increase in retail rates will improve the economics of DPV investment to customers (i.e., lower payback time for PV system) which, all else being equal, would be expected to drive greater PV adoption and thus lead to increased reductions in the utility's future load. Though these effects have been found to be small [52]. Darghouth et al. [52] also addressed a separate feedback mechanism between increasing PV penetration and the timing of time-of-use (TOU) periods; their analysis found that greater PV penetration causes TOU peak periods to shift into the evening hours, which in turn dampens further adoption. Second, we do not represent the feedback effects of changes in retail rate designs or NEM alternatives on overall customer energy consumption (e.g., fixed customer charge may reduce financial incentive to invest in energy efficiency or net billing may encourage DPV system design to maximize exports), all else being equal. Instead, we construct an annual load and PV penetration forecast which is adhered to regardless of these factors. Incorporating such changes into the model will be pursued as a future research effort.

Most Western U.S. utilities, with the exception of some of those in California, currently have distributed generation deployments equivalent to less than 1% of annual retail sales. It will take them several years to see the kinds of impacts depicted in this analysis, at even the 4% penetration level, let alone the 8% level. Accordingly, policymakers and regulators likely have time to study and deliberate changes to NEM, as well as other policy and regulatory actions, before observing material financial impacts. While many of the impacts appear relatively small (on a percentage basis), they demonstrate how underlying ratemaking and regulatory mechanisms can change utility support

for or customer interest in DERs, and the magnitude of impacts depends critically on utility physical, financial, and operating characteristics.

**Author Contributions:** P.C. wrote the original draft manuscript, contributed to the project's conceptualization, performed a subset of the formal analysis, and reviewed and edited the final manuscript. A.S. provided supervision of the research activity, led the development of the project's conceptualization, performed a subset of the formal analysis, and reviewed and edited the manuscript. W.G. and J.R. contributed to the project's conceptualization, performed a subset of the formal analysis, and reviewed and edited the manuscript. All authors read and approved the final manuscript.

**Funding:** This work was supported by the U.S. Department of Energy Solar Energy Technologies Office under Lawrence Berkeley National Laboratory Contract No. DE-AC02-05CH11231. This manuscript has been authored by an author at Lawrence Berkeley National Laboratory under Contract No. DE-AC02-05CH11231 with the U.S. Department of Energy. The U.S. Government retains, and the publisher, by accepting the article for publication, acknowledges, that the U.S. Government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for U.S. Government purposes.

**Conflicts of Interest:** The authors acknowledge that one of them is also an editor for this journal.

#### **Appendix A FINDER Model Overview**

The FINDER pro forma financial model was developed to quantify the financial impacts on ratepayers, utilities, and shareholders from the deployment of DERs as well as the introduction of alternative regulatory/business models. The basic structure of the model is depicted in Figure A1. What follows is a more detailed description of the different modules within the model.

Utility costs are based on model inputs that characterize current and projected utility costs over the analysis period. The model represents major cost categories of the utility's physical, financial, and operating environment, including fuel and purchased power, operations and maintenance, and capital investments in generation and non-generation assets (i.e., transmission and distribution investments). Some costs are projected using stipulated first year values and compound annual growth rates (CAGRs); other costs are based on schedules of specific investments (e.g., generation expansion plans). The model calculates the utility's ratebase over the analysis period, accounting for increases due to additional capital investments as well as decreases due to depreciation of existing assets. The model also estimates interest payments for debt and returns for equity shareholders based on an authorized amount used to finance capital investments and includes taxes on earnings.

**Figure A1.** FINDER Model Overview.

The utility's collected revenues are based on retail rates that are set in periodic general rate cases (throughout the analysis period. By default, the model assumes that a GRC occurs at some specified frequency (e.g., every three years); the model also allows the utility to file a GRC that may be triggered by a significant capital investment (e.g., new power plant, proposal to install advanced metering infrastructure).

GRCs are used to establish new rates for customers based on the revenue requirement set in a test year, including an authorized ROE for capital investments, the test year billing determinants (i.e., retail sales, peak demand, and number of customers), and assumptions about how the test year revenue requirement is allocated to customer classes and among the billing determinants. The model allows for different types of test years (i.e., historical, current, and future test years). Many states allow the utility to file an adjustment to its historical test-year costs during a GRC (i.e., pro forma adjustment period) to update and correct them to reflect expectations about normal cost levels, however, our model uses unadjusted historic test year values for ratemaking purposes. The particular rate design of the utility consists of a combination of a volumetric energy charge (\$/kWh), volumetric demand charge (\$/kW), and fixed customer charge (\$/customer) for a particular customer class. Model inputs specify the relative share of different types of utility costs that are collected from each of these three rate components.

The rates established in a GRC are then applied to the actual billing determinants in future years to calculate utility collected revenue in those years. The model accounts for a period of regulatory lag whereby rates that are established in a GRC do not go into effect until some specified number of years after the GRC. In between general rate cases, certain costs are passed directly to customers through rate-riders (e.g., fuel-adjustment clause or FAC). The model derives an average all-in retail rate metric that reflects the average revenue collected per unit of sales at the utility or customer-class level and accounts for periodic setting of new rates, rate-riders, and delays in implementing new rates.

The financial performance of the utility is measured by achieved after-tax earnings and achieved after-tax ROE. We calculated the prototypical utilities' achieved after-tax ROE in each year as the current year's earnings divided by current year's outstanding equity (i.e., the equity portion of the ratebase). The FINDER Model does not take into account changes in financing costs that may result from underor over-recovery of costs, which may impact ROE. Achieved after-tax ROE may, and often does, differ from the utility's authorized ROE. The authorized ROE is typically established by regulators during a regulatory proceeding and used in a GRC to determine the amount of return that a utility may receive on its capital investments. Actual utility revenues and costs may—and nearly always do—differ from those in the test year, leading to achieved earnings, and hence achieved ROE, that deviates from the authorized level. In a GRC, utility rates are set such that the test-year revenue requirement produces earnings that are sufficient to reach the authorized after-tax ROE based on the test year costs and billing determinants. In general, achieved ROE will be less than authorized ROE if, between rate cases, utility costs grow faster than revenues. Conversely, achieved ROE will generally be greater than authorized ROE if utility costs grow slower than revenues between rate cases.

FINDER calculates the prototypical utilities' achieved after-tax earnings as collected revenues minus costs in each year. Achieved after-tax earnings can be different than the utility's authorized earnings, because the achieved earnings are based on actual profitability in a given year and the authorized earnings are set in the GRC revenue requirement, based on the authorized ROE. Technically, state regulators do not explicitly authorize earnings in a GRC; they authorize a ROE, which when applied to the undepreciated portion of a utility's share of equity financed ratebase produces a level of earnings. For our purposes in this report, we refer to that as authorized earnings.

Alternative regulatory mechanisms and rate structures can also be implemented in FINDER: decoupling mechanisms (i.e., sales based or revenue-per-customer), lost revenue adjustment mechanisms, and shareholder incentive mechanisms. Alternative rate structures (e.g., high fixed customer charge) are represented by changing the way utility revenues are collected among different billing determinants.

#### **Appendix B Detailed Sensitivity Case Results**

Table A1 shows the sensitivity case results expressed in percentage changes. The change in DPV value is relative to the DPV value in the Base case at the respective DPV deployment level (i.e., 4% or 8%). The change in achieved earnings and average all-in retail rates is relative to a case without DPV.


**Table A1.** Full sensitivity case results at 4% and 8% DPV.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Strategic-Agent Equilibria in the Operation of Natural Gas and Power Markets**

**Sheng Chen 1,\* and Antonio J. Conejo <sup>2</sup>**


Received: 20 January 2020; Accepted: 11 February 2020; Published: 17 February 2020

**Abstract:** We consider strategic gas/power producers and strategic gas/power consumers operating in both gas and power markets. We build a flexible multi-period complementarity model to characterize day-ahead equilibria in those markets. This model is an equilibrium program with equilibrium constraints that characterizes the market behavior of all market agents. Using a realistic case study, we analyze equilibria under perfect and oligopolistic competition. We also analyze equilibria under different levels of information disclosure regarding market outcomes. We study as well equilibria under different ownership schemes: no hybrid agent, some hybrid agents, and only hybrid agents. Finally, we derive policy recommendations for the regulators of both the gas and the power markets.

**Keywords:** natural-gas market; electricity market; equilibrium analysis

#### **1. Introduction**

Electric power systems and natural-gas systems are generally operated independently, with limited or no coordination [1]. This is the result of how these systems were created and have evolved over time. In fact, gas has not played a significant role as a primary fuel in electricity production until recently, and thus, gas and power system coordination has not been important until recently.

Due to the increasing availability of gas and its competitive price, during the last decade, an increasing number of combined cycle gas turbines (CCGTs) have been incorporated into the generation mix of many power systems. This has resulted in an increasingly strong interdependency between gas systems and power systems [2]. In fact, this interdependency can no longer be disregarded if the gas and power systems are to be operated efficiently [1].

However, tools to comprehend the effect of such interdependency are limited. Many of these tools adopt a centralized perspective, in which a single operator manages both the gas and power systems [3–14] , which is unrealistic. Representative references are briefly discussed below. Chen et al. [3] develop a unit commitment model that includes an enhanced second order conic gas flow model, where the interdependency between gas and power prices is investigated. Byeon and Van Hentenryck [4] introduce a unit commitment problem with gas network awareness, where bid-validity constraints are imposed on gas-fired units. He et al. [5] propose an integrated gas and power system operation model that considers demand response and uncertainty via distributionally robust optimization. He et al. [6] develop a decentralized operation model for multi-area gas and power systems. Chen et al. [7] develop a joint gas and power market model that addresses wind power uncertainty and gas system congestion. Ameli et al. [8] quantify the value of the flexibility of the gas system in accommodating intermittent renewable energy sources. Yang et al. [9]

propose a two-stage robust operation model that considers gas network dynamics and wind power uncertainty. Bai et al. [10] develop a robust scheduling model that considers N-1 contingencies of power transmission lines or gas pipelines. Zlotnik et al. [11] analyze the economic and security benefits of a coordinated scheduling of interdependent gas and power systems. Chen et al. [12] propose a two-stage robust day-ahead dispatch model for urban electric and gas systems. Antenucci and Sansavini [13] investigate the impacts of gas-system operational constraints on a stochastic unit commitment model with large renewable penetration. Ordoudis et al. [14] develop an integrated electricity and gas market-clearing model, in which the value of the gas system flexibility to accommodate high shares of renewables is discussed.

Complementarily, Ref. [15] proposes an equilibrium model of the type we propose in this paper, but for distribution systems and [16] describes an equilibrium model at the bulk level, but uses and heuristic solution approach.

We propose in this paper an equilibrium model that allows studying the interactions of both gas/power producers and gas/power consumers (referred generically to as agents) through both the gas and the power markets. This model expands the one reported in [17] as it considers a multi-period framework and carries out a comprehensive analysis. Each market agent (producer of consumer) is represented as a bi-level model (see Appendix A.3 of the Appendix) with an upper-level problem that pursues maximum profit (revenue minus cost or utility minus payment) for the agent (see Appendices A.3.1 and A.3.2, respectively of the Appendix), and two lower-level problems representing the clearing of the gas and the power markets (see Appendices A.1 and A.2 , respectively, of the Appendix). We then jointly consider the bi-levels problems of all the agents participating in the gas and power markets, and solve the resulting Equilibrium Problem with Equilibrium Constraints (EPEC) using a direct approach [18,19] that does not rely on heuristics.

We consider hybrid producers that own both gas and power production facilities as well as non-hybrid ones. Likewise we consider hybrid consumers that consume both gas and electricity and non-hybrid ones.

The study horizon that we consider for both the gas and the power markets is one day divided in a number of periods to capture inter-temporal effects, such as steep ramping requirements due to the variability of the production of renewable units.

The proposed model represents in detail the gas and power network, the latter using linear (dc) equations (see Appendix A.2 of the Appendix) and the former via second order conic equations (see Appendix A.1 of the Appendix).

We consider that gas/power producers and gas/power consumers are both strategic and seek to alter gas/power clearing prices to their respective benefits and analyze equilibria under three conditions, namely:


The equilibrium analysis reported in this paper is particularly relevant to the regulator, as it helps devising market adjustments and coordination rules to maximize social welfare in both the gas and power markets.

The contributions of this paper are twofold:


The rest of this paper is organized a follows. Section 2 describes in a generic manner the bi-level model of a strategic agent (producer or consumer), Section 3 describes the considered EPEC, Section 4 shows how to solve it, Section 5 provides an illustrative example, Section 6 and 7 describe and discuss results from two realistic test systems, and Section 8 draw conclusions. The Appendix provides detailed descriptions of the models considered and metrics used.

#### **2. Single-Agent Model**

A generic bi-level model to represent the profit-seeking behavior of a single strategic agent (producer or consumer) is provided below:

$$\max\_{\Xi^{(i)}} \quad \pi^{(i)}\left(\mathbf{x}\_{\mathcal{S}}^{(i)}, \mathbf{x}\_{\mathcal{E}}^{(i)}, \boldsymbol{\lambda}\_{\mathcal{S}}^{(i)}, \boldsymbol{\lambda}\_{\mathcal{E}}^{(i)}\right) \tag{1}$$

$$\text{s.t.} \qquad o\_{\mathcal{S}}^{(i)} \in \mathcal{O}\_{\mathcal{S}}^{(i)}, o\_{\mathcal{e}}^{(i)} \in \mathcal{O}\_{\mathcal{e}}^{(i)} \tag{2}$$

$$\min\_{\mathbf{x}\_{\mathcal{S}}} \qquad \qquad f\_{\mathcal{S}}(\mathbf{x}\_{\mathcal{S}}, a\_{\mathcal{S}}) \tag{3}$$

$$\text{s.t.}\qquad\qquad h\_{\mathcal{S}}(\mathbf{x}\_{\mathcal{S}}) = 0 \text{ : } \lambda\_{\mathcal{S}}\tag{4}$$

$$\mathcal{g}\_{\mathcal{S}}(\mathbf{x}\_{\mathcal{S}}, o\_{\mathcal{S}}) \le 0: \ \mu\_{\mathcal{S}} \tag{5}$$

$$\min\_{\mathbf{x}\_{\mathcal{E}}} \qquad \qquad f\_{\mathcal{E}}(\mathbf{x}\_{\mathcal{E}}, o\_{\mathcal{E}}) \tag{6}$$

$$\text{s.t.}\qquad\qquad h\_{\mathfrak{e}}(\mathbf{x}\_{\mathfrak{e}}) = 0 \;:\; \lambda\_{\mathfrak{e}}\tag{7}$$

$$g\_{\ell}(\mathfrak{x}\_{\ell}, o\_{\ell}) \leq 0 : \ \mu\_{\ell} \tag{8}$$

where <sup>Ξ</sup>(*i*) <sup>=</sup> {*<sup>o</sup>* (*i*) *<sup>g</sup>* , *o* (*i*) *<sup>e</sup>* }∪{*xg*, *xe*, *λg*, *μg*, *λe*, *μe*}.

The notation used is described below:

*<sup>π</sup>*(*i*) (·) is the profit of agent (*i*),

*xg* the vector of gas variables,

*x* (*i*) *<sup>g</sup>* the sub-vector (of vector *xg*) of gas variables that pertains to agent (*i*),

*λg*, *μ<sup>g</sup>* vectors of dual gas variables,

*λ*(*i*) *<sup>g</sup>* the sub-vector (of vector *λg*) of dual gas variables that pertains to agent (*i*),

*xe* the vector of power variables,

*x* (*i*) *<sup>e</sup>* the sub-vector (of vector *xe*) of power variables that pertains to agent (*i*),

*λe*, *μ<sup>e</sup>* vectors of dual power variables,

*λ*(*i*) *<sup>e</sup>* the sub-vector (of vector *λe*) of dual power variables that pertains to agent (*i*),

*og* the gas offer/bid vector,

*o* (*i*) *<sup>g</sup>* the gas offer/bid sub-vector (of vector *og*) pertaining to agent (*i*),

*oe* the power offer/bid vector,

*o* (*i*) *<sup>e</sup>* the power offer/bid sub-vector (of vector *oe*) pertaining to agent (*i*),

O(*i*) *<sup>g</sup>* the feasible set of gas offers/bids of agent (*i*), and

O(*i*) *<sup>e</sup>* the feasible set of power offers/bids of agent (*i*).

Upper-level problem (1) and (2) represents the profit of the agent (revenue minus cost for a producer and utility minus payment for a consumer), while lower-level problems (3)–(5) and (6)–(8) represent the clearing of the gas and power markets, respectively.

The detailed models of a strategic gas/power consumer and a strategic gas/power producer are provided in Appendices A.3.1 and A.3.2, respectively, of the Appendix. Detailed descriptions of the gas clearing model (3)–(5) and the power clearing model (6)–(8) are provided in Appendices A.1 and A.2, respectively, of the Appendix.

*Energies* **2020**, *13*, 868

Assuming that lower-level problems (3)–(5) and (6)–(8) are convex or have been convexified [21], we replace them with their corresponding Karush-Kuhn-Tucker (KKT) optimality conditions [18,19,22], rendering the Mathematical Program with Equilibrium Constraints (MPEC) below:

$$\max\_{\pi\_{\Xi^{(i)}}} \quad \pi^{(i)}\left(\mathbf{x}\_{\mathcal{S}}^{(i)}, \mathbf{x}\_{\varepsilon}^{(i)}, \lambda\_{\mathcal{S}'}, \lambda\_{\mathcal{E}}\right) \tag{9}$$

$$\text{s.t.} \qquad o\_{\mathcal{S}}^{(i)} \in \mathcal{O}\_{\mathcal{S}}^{(i)}, o\_{\varepsilon}^{(i)} \in \mathcal{O}\_{\varepsilon}^{(i)} \tag{10}$$

$$
\nabla\_{\mathbf{x}\_{\mathcal{S}}} f\_{\mathcal{S}}(\cdot) + \lambda\_{\underline{\mathcal{S}}}^{\top} \nabla\_{\mathbf{x}\_{\mathcal{S}}} h\_{\mathcal{S}}(\cdot) + \mu\_{\underline{\mathcal{S}}}^{\top} \nabla\_{\mathbf{x}\_{\mathcal{S}}} \mathbf{g}\_{\mathcal{S}}(\cdot), \quad h\_{\mathcal{S}}(\cdot) = 0, \quad 0 \le \mu\_{\mathcal{S}} \perp \mathbf{g}\_{\mathcal{S}}(\cdot) \le 0 \tag{11}
$$

$$
\nabla\_{\mathbf{x}\_{\varepsilon}} f\_{\varepsilon}(\cdot) + \lambda\_{\varepsilon}^{\top} \nabla\_{\mathbf{x}\_{\varepsilon}} h\_{\varepsilon}(\cdot) + \mu\_{\varepsilon}^{\top} \nabla\_{\mathbf{x}\_{\varepsilon}} g\_{\varepsilon}(\cdot), \quad h\_{\varepsilon}(\cdot) = 0, \quad 0 \le \mu\_{\varepsilon} \perp g\_{\varepsilon}(\cdot) \le 0,\tag{12}
$$

Since MPEC (9)–(12) might be complex to solve/transform and considering that the gas problem is formulated as a second order conic problem (SOCP) [21] and that the power problem is formulated as a linear programming problem, each of these problems can be replaced by its primal constraints, its dual constraints, and its strong duality equality. Thus, instead of considering (9)–(12), we consider:

$$\max\_{\mathbf{x}\_{\omega(i)}} \quad \pi^{(i)}\left(\mathbf{x}\_{\mathcal{S}}^{(i)}, \mathbf{x}\_{\mathfrak{e}}^{(i)}, \lambda\_{\mathcal{S}'}, \lambda\_{\mathfrak{e}}\right) \tag{13}$$

$$\text{s.t.} \qquad o\_{\mathcal{S}}^{(i)} \in \mathcal{O}\_{\mathcal{S}}^{(i)}, o\_{\varepsilon}^{(i)} \in \mathcal{O}\_{\varepsilon}^{(i)} \tag{14}$$

$$\text{primual-constraints}\_{\mathcal{K}'} \quad \text{dual-constraints}\_{\mathcal{K}} \quad \text{strong-duality-equality}\_{\mathcal{K}} \tag{15}$$

$$\text{\\_primal-constraints}\_{\text{\\_}} \quad \text{dual-constraints}\_{\text{\\_}} \quad \text{strong-duality-equality}\_{\text{\\_}} \tag{16}$$

Problem (13)–(16) is generally better behaved than problem (9)–(12), and the KKT optimality conditions of (13)–(16) (single agent optimality conditions) are easily obtained [22] and represented as:

$$\text{KKT}^{(i)}\tag{17}$$

Deriving KKT conditions is a relatively simple exercise. For example, the solver EMP (Extended Mathematical Programming), (https://www.gams.com/latest/docs/UG\_EMP.html) which is available in GAMS (General Algebraic Modeling System) (https://www.gams.com), derives KKT conditions automatically.

We note that since problem (13)–(16) is generally non-convex and its constraints might be non-regular, its optimality conditions as given by (17) identify points that might or might not be extrema.

#### **3. Multiple-Agent Model: EPEC**

To search for equilibria, we jointly solve (17) for all market agents, which constitutes an Equilibrium Problem with Equilibrium Constraints (EPEC) [19]. This is expressed as:

$$\begin{cases} \text{KKT}^{(i)} & \forall i \end{cases} \tag{18}$$

which is a system of nonlinear equalities and inequalities difficult to solve. How to solve EPEC (18) is addressed in Section 4 below.

We note that since the constraints of problem (13)–(16) might be non-regular, (18) identifies equilibria and other stationary points [23].

#### **4. EPEC Solution**

The auxiliary problem below can be used to solve (18), i.e., to search for equilibria:

$$\mathbf{\color{red}{\bf max}} \quad o(\cdot) \tag{19}$$

$$\text{s.t.} \qquad \text{KKT}^{(i)} \quad \forall i,\tag{20}$$

where *o*(·) is a suitable objective function. We consider in the example and case study (Section 5 and 6, respectively) three objective functions (19), namely:


The corresponding EPECs (19)–(20) are referred to as:


Since (19)–(20) is generally nonlinear and non-convex, its solution can be attempted via linearization or using global solvers, such as BARON [24].

Once potential equilibrium points (solutions of (19)–(20)) are found, a diagonalization algorithm [22] can be used to verify if these points are indeed equilibria.

#### **5. Illustrative Example**

For the sake of illustration, we consider in this section a simple example. We analyze a two-bus power system (bus is used to refer to a power-system node) and a two-node gas system (node is used to refer to a gas node), the topology of which is shown in Figure 1. The gas-fired power unit at power bus 2 receiving gas from gas node 2 couples the two systems.

We consider two hybrid agents:


For simplicity, we do not consider strategic bids by consumers in this example. In addition, we consider a perfect gas price information interchange between the gas market and the owner of gas-fired power unit 2 (Agent 2).

**Figure 1.** Example: two-bus power system and two-node gas system.

#### *5.1. Data*

The capacities of the two power units at buses 1 and 2 are 50 MW and 20 MW, respectively. The marginal production cost of the power unit at bus 1 is 18 \$/MWh. The non-fuel cost of the gas-fired unit at bus 2 is 1 \$/MWh, and its energy conversion coefficient associated with gas consumption is 0.0045 Mm3/MWh.

Regarding the two gas sources at nodes 1 and 2, their capacities are 0.5 Mm3/h and 0.7 Mm3/h, respectively, and their marginal production cost are 3000 \$/Mm3 and 3500 \$/Mm3, respectively.

The transmission capacity of the power transsmission line connecting buses 1 and 2 is 18 MW. The lower and upper gas pressure limits at gas nodes are 25 bar and 40 bar, respectively. We note that these gas nodal pressure bounds do not restrict the gas flows through the pipeline connecting nodes 1 and 2.

The baseline utility of the power demands at buses 1 and 2 are 30 \$/MWh and 35 \$/MWh, respectively. The baseline utility of the gas demands at buses 1 and 2 are 4000 \$/Mm3 and 4200 \$/Mm3, respectively. The marginal utility factors of both gas and power demands during time periods 1–8, 9–16, and 17–24 are 0.8, 1.0, and 1.2 relative to their baseline utilities, respectively.

Finally, Figure 2 depicts the 24-h total non-generation-related gas demand and the total electricity demand.

**Figure 2.** Example: non-generation-related gas demand and power demand.

#### *5.2. Results*

We considered two equilibrium models (19)–(20), whose objective functions were total producers' profit and social welfare of both markets, i.e., Max TPP EPEC, and Max SW EPEC, respectively. Table 1 summarizes the market equilibria obtained from the two models. We observed that the equilibrium model that maximized TPP yielded a lower SW but a higher TPP than the corresponding SW and TPP obtained from the equilibrium model that maximized SW. In addition, these two equilibrium models resulted in differences in the distribution of profits between the two production agents. Specifically, Agent 1 earned a higher profit from the model that maximized SW, while the model that maximized TPP was more beneficial for Agent 2.

**Table 1.** Example: profits and social welfare (\$ thousand) for two equilibrium models.


Additionally, we considered a gas-shortage case, where the capacity of gas-fired unit 2 was reduced to 10 MW. Table 2 provides results for the base case and the gas-shortage case obtained from the Max TPP EPEC. This table shows that the gas-shortage case resulted in a higher profit for Agent 1, earned from the electricity market. This is because power unit 1 accounted for an increased share of electricity supply. Additionally, the gas shortage resulted in lower profits for the two production agents earned from the gas market due to reduced generation-related gas demands.

These results show how the operation of the gas system impacts production agents' profits earned from both gas and power markets. In practice, gas-fired power producers should be aware of potential gas-system bottlenecks, which determine the availability and reliability of their fuel supply.


**Table 2.** Example: profits and social welfare (\$ thousand) for two cases (Max Total Producers' Profit (TPP) Equilibrium Problem with Equilibrium Constraints EPEC).

\* Agent 1/2 (E) and Agent 1/2 (G) represent Agent 1's/2's profits earned in the power and gas markets, respectively.

The EPEC model (19)–(20) was solved using BARON [24] under GAMS on a computer with a 2.1-GHZ Intel Core-i7 processor with 8 GB of memory. The solution time of any instance analyzed was below 190 seconds.

#### **6. Case Study**

This section examines a case study comprising the IEEE-57 bus system [25] and a tree-like 134-node Greek gas system (http://gaslib.zib.de/).

We consider (i) strategic offers/bids from both producers and consumers, (ii) disaggregated and aggregated gas price information, and (iii) diverse ownership of gas and power production units.

Taking into account the computational machinery used and for the sake of simplicity and tractability, we consider a time horizon of 3 h.

#### *6.1. Data*

The gas system consists of three gas sources, 45 demand nodes, 132 pipelines, and one gas compressor. The power system includes seven power units, being the units at buses 1, 2 and 3 gas-fired and connected to gas nodes 2, 8, and 15, respectively. This system includes 22 demand nodes and 80 transmission lines.

We consider three strategic agents, agents 1 and 2 being hybrid producers, and agent 3 a hybrid consumer. Specifically:


All power units and gas sources are owned by either by Agent 1 or 2 and submit strategic offers. However, a number of electricity/gas demands are not owned by Agent 3, and hence bid competitively.

#### *6.2. From Perfect to Oligopolistic Competition*

Table 3 summarizes the market equilibria obtained from the competitive model and three oligopolistic models:


**Table 3.** Case study: profits and social welfare (\$ thousand) under different equilibrium models.


\* In the competitive model, all agents are non-strategic and offer/bid at their marginal costs/utilities.

Figures 3 and 4 provide the load-weighted electricity and gas locational marginal prices (LMPs), respectively, obtained from the four models.

The results obtained allow the following conclusions:


The EPEC models that maximized SW, TPP, and TCP required approximately 1681 s, 4123 s, and 3124 s, respectively, of wall-clock time to solve.

**Figure 3.** Case study: load-weighted electricity locational marginal prices (LMPs) obtained from four equilibrium models.

**Figure 4.** Case study: load-weighted gas LMPs obtained from four equilibrium models.

#### *6.3. Aggregated Gas Prices*

This subsection investigates the impact of temporal/spatial aggregation of gas prices on the market equilibria reached.

Considering the Max SW EPEC, Table 4 summarizes the market equilibria obtained from perfect pricing, spatial averaging pricing, temporal averaging pricing, and combined spatial and temporal averaging pricing.

The spatial averaging pricing derived a single price per hour by performing a load-weighting average across nodes of all gas LMP that hour (see (A32) in the Appendix). Similarly, the temporal averaging pricing derived a single price per node by performing a load-weighting average across hours of all gas LMP in that node (see (A33) in the Appendix). Finally, the combined spatial and temporal averaging pricing did both, deriving a single gas price per day (see (A34) in the Appendix).

We observe from Table 4 that the imperfect-pricing cases resulted in lower SW. Specifically, both spatial averaging pricing and temporal averaging pricing models yielded a lower TPP and a slightly higher TCP. However, the combined averaging pricing model resulted in a loss of both TPP and TCP.


**Table 4.** Case study: profits and social welfare (\$ thousand) for a perfect pricing and three imperfect pricing cases (Max SW EPEC).

These results show that highly granular pricing practices are desirable to co-ordinate gas and power markets. This is so because such practices prevent loss of SW and increased profits of gas/power producers.

#### *6.4. Ownership Structure*

We investigate in this section the impact of ownership structure on market equilibria. This was done by considering three cases involving all hybrid agents, some hybrid agents, and no hybrid agent. The Max TPP EPEC was considered. Table 5 describes the three cases considered.

**Table 5.** Case study: ownership structure. A: power units at buses 1–3 and 12. B: power units at buses 6, 8, and 9. C: gas sources at nodes 1 and 20. D: gas source at node 80.


The resulting market equilibria are provided in Table 6. This table shows that the all hybrid agents' cases resulted in the highest TPP and the lowest SW. In comparison, the case of no hybrid agent resulted in the lowest TPP and the highest SW. These changes in TPP and SW are due to differences in the market power exercised by gas/power producers. In the all hybrid agents' cases, each agent accounted for a larger gas/power production capacity, and thus it could potentially exercise higher market power to its own profit, which, consequently, reduced the SW.

**Table 6.** Case study: profits and social welfare (\$ thousand) for different market ownership cases (Max TPP EPEC).


#### **7. Case Study 2**

This section summarizes numerical results from a realistic Belgian 24-node power system and 20-node gas system [17], the topology of which is shown in Figure 5. The power units at buses 2, 3, 6, 8, 16, 15, and 22 are gas-fired and connected to nodes 4, 3, 4, 4, 6, 11, and 13, respectively. We considered three strategic producers: agents 1, 2 and 3. Agent 1 owned power units in area 1 (see upper left-hand-side of Figure 5); agent 3 owned gas sources in area A (see upper right-hand-side of Figure 5); agent 2 owned power units in area 2 (see lower left-hand-side of Figure 5) and gas sources in

area B (see lower right-hand-side of Figure 5). A fourth strategic agent owned electricity demands at buses 7, 9, 23, and 24 and gas demands at nodes 10, 12, 19, and 20. We considered a time horizon of 6 h.

**Figure 5.** Case study 2: Belgian 24-node power system and 20-node gas system.

We investigated the impact of gas-pressure limits on the market equilibria reached. This was done by comparing the results obtained from two cases, in which the ranges of nodal gas pressures were between 30 bar and 70 bar and between 35 bar and 65 bar, respectively. Table 7 and Figure 6 summarize the equilibrium results obtained from the two cases. These results indicate that a strict gas-pressure limit resulted in 1) a lower TPP, TCP, and SW, 2) higher gas LMPs, and 3) lower profits of agents 1 and 2 obtained from the power market owing to increased fuel cost for gas-fired units.

**Table 7.** Case study 2: profits and social welfare (\$ thousand) for two sets of gas pressure limits (Max TPP EPEC).


\* Agent 2 (E) and Agent 2 (G) represent Agent 2's profits earned in the power and gas markets, respectively.

**Figure 6.** Case study 2: gas LMPs at the peak time period for two sets of gas pressure limits (Max TPP EPEC).

#### **8. Conclusions**

This paper proposes a multi-period EPEC model to analyze the interactions of both strategic power/gas producers and power/gas consumers that participate in power and gas markets. We investigate the impacts of (i) market power, (ii) aggregated gas prices and (iii) ownership structure of power/gas producers on the market equilibria reached. From the analysis carried out, the conclusions below are in order:


**Author Contributions:** Conceptualization, A.J.C. and S.C.; methodology, A.J.C. and S.C.; software, S.C. and A.J.C.; validation, S.C. and A.J.C.; formal analysis, A.J.C. and S.C.; investigation, S.C. and A.J.C.; data curation: S.C. and A.J.C.; writing-original draft preparation, A.J.C. and S.C.; writing–review and editing, A.J.C. and S.C.; visualization, S.C. and A.J.C.; supervision, A.J.C. and S.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** A. J. Conejo is partly supported by NSF Grant 1808169. S. Chen is partly supported by Fundamental Research Funds for the Central Universities under Grant 2019B05714.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Detailed Models**

#### *Appendix A.1. Gas Market Clearing*

An SOC (Second Order Conic) formulation of the gas operation problem (3)–(5) is:

$$\max\_{\boldsymbol{\Sigma}\_{\text{G}}^{\text{P}}} \sum\_{\boldsymbol{\varepsilon} \in \Lambda^{\text{GD}}, t \in T} \mathsf{C}\_{\text{et}}^{\text{GL}} F\_{\text{et}}^{\text{L}} + \sum\_{\boldsymbol{l} \in \mathbb{L}, \boldsymbol{\mu} \in \Lambda\_{\text{l}}^{\text{GL}}, t \in T} \mathsf{c}\_{\text{et}} F\_{\text{et}}^{\text{L}} + \sum\_{\boldsymbol{m} \in \mathbb{N}, \boldsymbol{\upmu} \in T} \left( \sum\_{\boldsymbol{v} \in \mathsf{Y}\_{\text{m}}^{\text{G}}} \gamma\_{\text{v}t}^{\text{G}} F\_{\text{vt}}^{\text{G}} - \sum\_{\boldsymbol{w} \in \mathsf{Y}\_{\text{m}}^{\text{S}}} \beta\_{\text{wt}} F\_{\text{wt}}^{\text{S}} \right) \tag{A1}$$

$$\text{s.t. } \sum\_{w \in \Psi\_m^S} F\_{\text{wt}}^S = \sum\_{k \in \mathbb{C}(m)} (1 + \theta\_k) F\_{\text{tf}}^C + \sum\_{c \in \Psi\_m^L} F\_{\text{tf}}^L + \sum\_{v \in \Psi\_m^G} F\_{\text{vt}}^G + \sum\_{n \in \mathbb{G}(m)} F\_{\text{mut}}; \forall m \in \mathbb{N}, t \in T \quad (u\_{\text{mt}}) \tag{A2}$$

$$\bar{F}\_{\rm mnt} = (F\_{\rm mnt} - F\_{\rm nmt}) / 2 ; \forall m, n \in \mathbb{N}, t \in T \tag{A3}$$

$$F\_{mnt} + F\_{mnt} = L\_{mnt} - L\_{mnt, t-1}; \forall m, n \in \mathbb{N}; t \in T \tag{A4}$$

$$L\_{m,n,t} = K\_{m,n} \cdot \left(\pi\_{m,t} + \pi\_{n,t}\right) / 2; \forall m, n \in \mathbb{N}; t \in T \tag{A5}$$

$$(\mathbb{F}\_{mnt}/\mathcal{W}\_{mn})^2 \le \Pi\_{mt}^2 - \Pi\_{nt}^2; \forall m \in \mathbb{N}, n \in \mathbb{G}(m), t \in T \tag{A6}$$

$$
\Pi\_{\rm kt}^{in} \rho\_{\mathbb{C},k}^{\rm min} \le \Pi\_{\rm kt}^{out} \le \Pi\_{\rm kt}^{in} \rho\_{\mathbb{C},k}^{\rm max}; \forall k \in \mathcal{K}, t \in T \tag{A7}
$$

$$0 \le F\_{kt}^C \le F\_k^{C, \max}; \forall k \in \mathcal{K}, t \in T \tag{A8}$$

$$0 \le F\_{wt}^{\mathbb{S}} \le \bar{F}\_w^{\mathbb{S}, \max}; \forall m \in \mathbb{N}, w \in \mathbb{Y}\_{m'}^{\mathbb{S}}, t \in T \tag{A9}$$

$$1 - F\_{\textit{w}}^{S,ramp} \le F\_{\textit{wt}}^{S} - F\_{\textit{w},t-1}^{S} \le F\_{\textit{w}}^{S,ramp}; \forall m \in \mathbb{N}, w \in \mathbf{Y}\_{\textit{m}}^{S}, t \in T \tag{A10}$$

$$0 \le F\_{et}^L \le F\_{et}^{L,\text{max}}; \forall e \in \Lambda^G, t \in T \tag{A1}$$

$$
\Pi\_m^{\min} \le \Pi\_{mt} \le \Pi\_m^{\max}; \forall m \in \mathbb{N}, t \in T \tag{A12}
$$

$$F\_{mnt} \ge 0; \forall m \in \mathbb{N}, n \in \mathbb{G}(m), t \in T \tag{A13}$$

$$0 \le F\_{vt}^G \le F\_v^{G, \max}; \forall m \in \mathbb{N}, \upsilon \in \Psi\_m^G, t \in T,\tag{A14}$$

where *e* is the index of gas demands in set Λ*GO*, *t* the index of operating periods in set *T*, *l* the set of agents in set L, *m* and *n* the indices of nodes in set N, *v* the index of power units, Ψ*<sup>G</sup> <sup>m</sup>* the set of gas-fired units connected to node *m*, Ψ*<sup>S</sup> <sup>m</sup>* the set of gas sources connected to node *m*, C(*m*) the set of gas compressors connected to node *m*, Ψ*<sup>L</sup> <sup>m</sup>* the set of gas demands connected to node m, G(*m*) the set of nodes that are connected directly to node m, *k* the index of gas compressors in the set K, and *w* the index of gas sources.

The parameters of the problem (A1)–(A14) are described below. *CGL et* is the marginal utility of demand *e* at time period *t*, *FC*,*max <sup>k</sup>* the gas transportation capacity of compressor *<sup>k</sup>*, *<sup>F</sup>S*,*max <sup>w</sup>* the production capacity of gas source *w*, *FL*,max *et* the quantity of gas demand *<sup>e</sup>* at time period *<sup>t</sup>*, *<sup>F</sup>G*,max *<sup>v</sup>* the maximum gas consumption of power unit *v*, *Kmn* the line-pack parameter of the pipeline connecting nodes *m* and *n*, *Wmn* the Weymouth constant of the pipeline connecting nodes *m* and *n*, *ρ*min *<sup>C</sup>*,*<sup>k</sup>* and *<sup>ρ</sup>*max *<sup>C</sup>*,*<sup>k</sup>* the minimum and maximum compression ratio of compressor *k*, *ϑ<sup>k</sup>* the conversion efficiency of gas compressor *k*, and Πmin *<sup>m</sup>* and Πmax *<sup>m</sup>* the minimum and maximum gas pressures of node *m*, respectively.

The variables of the problem (A1)–(A14) are as follows. *F<sup>L</sup> et* is the consumption of demand *e* in time period *t*, *F<sup>S</sup> wt* the production of gas source *w* in time period *t*, *F<sup>C</sup> kt* the gas flow through compressor *k* in time period *t*, *Fmnt* the gas flow through the pipeline connecting nodes *m* and *n* in time period *t*, *F*¯ *mnt* the average gas flow through the pipeline connecting nodes *m* and *n* in time period *t*, *Lmnt* the line-pack in pipeline connecting nodes *m* and *n* in time period *t*, *εet* the bid of demand *e* in time period *t*, *γ<sup>G</sup> vt* the bid of gas-fired power unit *v* in time period *t*, and *βwt* the offer of gas source *w* in time period *t*.

It should be noted that the variables *εet*, *γ<sup>G</sup> vt*, and *βwt* are fixed by upper-level problems and thus are constant for this problem. *umt* denotes the dual variable of (A2), and represents the gas LMP of node *m* at time period *t*. The variables of problem (A1)–(A14) are those in set ΞP <sup>G</sup> = *FL et*, *F<sup>S</sup> wt*, *F<sup>C</sup> kt*, *Fmnt*, *<sup>F</sup>*¯ *mnt*, *Lmnt* .

The objective function (A1) is the gas SW that incorporates strategic offers from gas producers and strategic bids from power producers and gas consumers. Constraints (A2) represent the gas nodal balances. Constraints (A3) calculate average gas flows through pipelines. Constraints (A4) give the relationship between hourly changes in flows and line-pack in pipelines. Constraints (A5) determine the hourly line-pack in each pipeline, which is considered to be linear with the average gas pressure at the two ends of the pipeline. Constraints (A6) relate the average gas flow with the change in squared gas pressures between the upstream and downstream nodes for each pipeline. (A6) represent an SOC formulation of an exact gas flow model [3]. Constraints (A7) enforce minimum and maximum gas pressure ratios of gas compressors. Constraints (A8) impose transportation limits on gas compressors. Constraints (A9) and (A10) impose production capacities and ramping limits on gas sources, respectively. Constraints (A11) limit the amount of gas demands served. Constraints (A12) enforce minimum and maximum gas pressures of each node. Constraints (A13) assume that the direction of gas flows are known a priori, which is generally reasonable in short-term operations [3]. Constraints (A14) limit the amount of generation-related demands.

#### *Appendix A.2. Power Market Clearing*

An LP (Linear Programming) formulation of the power operation problem (6)–(8) is:

$$\max\_{\mathsf{T}\_{\mathrm{E}}^{\mathrm{P}}} \sum\_{l \in \mathbb{L}, d \in \Lambda\_{l}^{\mathrm{EL}}, t \in T} \varsigma\_{d} t P\_{d\boldsymbol{t}}^{\mathrm{L}} + \sum\_{d \in \Lambda^{\mathrm{EO}}, t \in T} \mathbb{C}\_{d\boldsymbol{t}}^{\mathrm{EL}} P\_{d\boldsymbol{t}}^{\mathrm{L}} - \sum\_{v \in \Omega^{\mathrm{E}}, t \in T} a\_{v\boldsymbol{t}} P\_{v\boldsymbol{t}}^{\mathrm{G}} \tag{A15}$$

$$\text{s.t. } \sum\_{d \in \Theta\_i^D} P\_{\text{d}t}^L + \sum\_{j \in \mathbb{E}(i)} b\_{ij} \cdot (\delta\_{\text{it}} - \delta\_{\text{jt}}) = \sum\_{v \in \Theta\_i^G} P\_{\text{v}t}^G; \forall i \in \mathbb{B}, t \in T \quad (\lambda\_{\text{it}}) \tag{A16}$$

$$0 \le P\_{dt}^{L} \le P\_{dt}^{L,\text{max}}; \forall d \in \Lambda^E, t \in T \tag{A17}$$

$$b\_{ij} \cdot (\delta\_{it} - \delta\_{jt}) \le P\_{ij}^{\max}; \forall i \in \mathbb{B}, j \in \mathbb{E}(i), t \in T \tag{A18}$$

$$0 \le P\_{vt}^G \le P\_v^{G, \max}; \forall v \in \Omega^E, t \in T \tag{A19}$$

$$P\_v^{G,ramp} \le P\_{vt}^G - P\_{v, t-1}^G \le P\_v^{G,ramp}; \forall v \in \Omega^E, t \in T \tag{A20}$$

$$
\delta\_{REF,t} = 0; \forall t \in T,\tag{A21}
$$

where *d* is the index of electricity demands in set Λ*E*, *i* and *j* the indices of electric buses in set B, *v* the index of power units in set Ω*E*, *REF* the index of the reference bus. Λ*EL <sup>l</sup>* the set of strategic consumers owned by agent *l*, Λ*EO* the set of non-strategic consumers, Θ*<sup>D</sup> <sup>i</sup>* the set of electricity demands directly connected to bus *i*, Θ*<sup>G</sup> <sup>i</sup>* the set of power units directly connected to bus *i*, and E(*i*) the set of buses directly connected to bus *i*.

The parameters of the problem (A15)–(A21) are described below. *CEL dt* is the marginal utility of demand *d* in time period *t*, *bij* the susceptance of the line connecting buses *i* and *j*, *PL*,max *dt* the quantity of demand *d* in time period *t*, *P*max *ij* the transmission capacity of the line connecting buses *i* and *j*, *<sup>P</sup>G*,max *<sup>v</sup>* and *<sup>P</sup>G*,*ramp <sup>v</sup>* the capacity and ramping limit of power unit *<sup>v</sup>*, respectively.

The variables of the problem (A15)–(A21) are as follows. *P<sup>L</sup> dt* is the quantity of demand *d* served in time period *t*, *P<sup>G</sup> vt* the power production of unit *v* in time period *t*, *δit* the phase angle of bus *i* in time period *t*, *ςdt* the bid of demand *d* in time period *t*, and *αvt* the offer of power unit *v* in time period *t*.

It should be noted that variables *ςdt* and *αvt* are determined by upper-level problems and thus are constants for this problem. *λit* denotes the dual variable of (A16), and represents the electricity LMP of bus *i* in time period *t*. The variables of problem (A15)–(A21) are those in the set Ξ<sup>P</sup> <sup>E</sup> = *PL dt*, *<sup>δ</sup>it*, *<sup>P</sup><sup>G</sup> vt* .

The objective function (A15) is the power SW that considers the strategic offers of power producers and the strategic bids of power consumers. Constraints (A16) represent the active power balances, in which the DC power flow model is used. Constraints (A17) impose upper limits on power demands. Constraints (A18) enforce the transmission capacity of each line. Constraints (A19) and (A20) impose production capacities and ramping limits on power units, respectively. Constraints (A21) fix the phase angle of the reference bus to zero.

#### *Appendix A.3. Agent Model*

#### Appendix A.3.1. Strategic Consumer

The problem of a gas/power strategic consumer is:

$$\max\_{\Xi\_{\text{UC}}} \sum\_{d \in \Lambda\_l^{\text{EL}}, t \in T} (\mathsf{C}\_{dt}^{\text{EL}} - \lambda\_{i(d)}, t) P\_{dt}^{L} + \sum\_{e \in \Lambda\_l^{\text{GL}}, t \in T} (\mathsf{C}\_{et}^{\text{GL}} - \mathsf{u}\_{m(e), t}) F\_{et}^{L} \tag{A22}$$

$$\text{s.t. } \varepsilon\_{tt} \ge 0; \forall \varepsilon \in \Lambda\_l^{GL}, t \in T \tag{A23}$$

$$
\xi\_{dt} \ge 0; \forall d \in \Lambda\_l^{EL}, t \in T \tag{A24}
$$

$$(A1) - (A21),$$

where *i*(*d*) is the bus at which electricity demand *d* is located, *m*(*e*) the node at which gas demand *e* is located, and ΞUC = ΞP <sup>G</sup>, <sup>Ξ</sup><sup>E</sup> <sup>G</sup>,*εet*, *ςdt* .

The objective function (A22) is the profit of consumer *l*. Constraints (A23) and (A24) represent the non-negative bids of strategic gas consumer *e* and electricity consumer *d*, respectively. Constraints (A25) enforce market constraints.

#### Appendix A.3.2. Strategic Producer

The problem of a gas/power strategic producer is:

$$\begin{split} \max\_{\boldsymbol{\Sigma}\_{\text{UP}}} & \sum\_{\boldsymbol{v} \in \Omega\_{l}^{G} \cup \Omega\_{l}^{R}, t \in T} \lambda\_{i(\boldsymbol{v}), t} P\_{\text{vt}}^{G} - \sum\_{\boldsymbol{v} \in \Omega\_{l}^{R}, t \in T} \mathbb{C}\_{\boldsymbol{v}}^{G} P\_{\text{vt}}^{G} - \sum\_{\boldsymbol{v} \in \Omega\_{l}^{G}, t \in T} (\mathbb{C}\_{\boldsymbol{v}}^{\text{O}} + \eta\_{\boldsymbol{v}} u\_{\text{vt}}^{\text{GE}}) P\_{\text{vt}}^{G} \\ & + \sum\_{\boldsymbol{w} \in \Omega\_{l}^{S}, t \in T} (\mathbb{C}\_{\boldsymbol{w}}^{S} - u\_{\text{m}(\boldsymbol{w}), t}) P\_{\text{wt}}^{S} \end{split} \tag{A26}$$

$$\text{s.t. } u\_{\text{vt}} \ge 0; \forall v \in \Omega\_l^{\mathbb{C}} \cup \Omega\_l^{\mathbb{C}}, t \in T \tag{A27}$$

$$\{\beta\_{wt} \ge 0; \forall w \in \Omega\_l^S, t \in T\}\tag{A28}$$

$$\gamma\_{vt}^{\mathbb{C}} \ge 0; \forall v \in \Omega\_l^{\mathbb{C}}, t \in T \tag{A29}$$

$$(A1) - (A21). \tag{A30}$$

where *i*(*v*) is the bus at which power unit *v* is located, *uGE vt* the gas price information that the gas system sends to power unit *v*, and *m*(*w*) the node at which gas source *w* is located, and ΞUP = ΞP <sup>G</sup>, <sup>Ξ</sup><sup>E</sup> <sup>G</sup>, *<sup>α</sup>vt*, *<sup>β</sup>wt*, *<sup>γ</sup><sup>G</sup> vt* .

The objective function (A26) is the profit of each producer *l*. Constraints (A27) and (A28) represent non-negative offers of power producer *v* and gas producer *w*, respectively. Constraints (A29) represent non-negative bids of gas-fired power unit *v* in the gas market. Constraints (A30) enforce market constraints.

#### *Appendix A.4. Perfect and Imperfect Gas Price Disclosure*

We consider perfect and imperfect price coordination between gas and power markets. Specifically, we consider different levels of gas price granularity. In the perfect-pricing case, exact gas price information is exchanged as:

$$
\mu\_{vt}^{GE} = \mu\_{m(v),t'} \; \forall l \in \mathbb{L}, \upsilon \in \Omega\_l^G, t \in T. \tag{A31}
$$

where *m*(*v*) denotes the node at which gas-fired power unit *v* is located.

In the imperfect-pricing cases, we consider spatial, temporal, and combined averaging of gas LMPs, in which the information interchange on gas prices are given by (A32), (A33), and (A34), respectively.

$$\mu\_{vt}^{IN} = \sum\_{\boldsymbol{\varepsilon} \in \Lambda^G} F\_{\boldsymbol{\varepsilon}t}^L \mu\_{m(\boldsymbol{\varepsilon}),t} / \sum\_{\boldsymbol{\varepsilon} \in \Lambda^G} \mathbb{F}\_{\boldsymbol{\varepsilon}t}^L; \ \forall l \in \mathbb{L}, \boldsymbol{\upsilon} \in \Omega\_l^G, t \in T \tag{A32}$$

$$u\_{\rm vt}^{IN} = \frac{1}{|T|} \sum\_{t \in T} u\_{m(v), t}; \ \forall l \in \mathbb{L}, v \in \Omega\_l^G, t \in T \tag{A33}$$

$$\mu\_{\rm tr}^{IN} = \frac{1}{|T|} \sum\_{t \in T} \left( \sum\_{\boldsymbol{\sigma} \in \Lambda^G} F\_{\boldsymbol{\varepsilon}t}^{\boldsymbol{L}} \mu\_{\boldsymbol{m}(\boldsymbol{\varepsilon}),t} / \sum\_{\boldsymbol{\sigma} \in \Lambda^G} F\_{\boldsymbol{\varepsilon}t}^{\boldsymbol{L}} \right); \; \forall \boldsymbol{l} \in \mathbb{L}, \boldsymbol{\upsilon} \in \Omega\_{\boldsymbol{l}}^G, \boldsymbol{t} \in T. \tag{A34}$$

Additionally note that the true SW is given by:

$$\sum\_{d \in \Lambda^{\mathbb{E}}, t \in T} \mathsf{C}\_{dt}^{\mathbb{E}L} P\_{dt}^{L} - \sum\_{v \in \Omega^{\mathbb{E}}, t \in T} \mathsf{C}\_{vt}^{\mathbb{G}} P\_{vt}^{G} + \sum\_{e \in \Lambda^{\text{CO}}, t \in T} \mathsf{C}\_{et}^{GL} F\_{vt}^{L} - \sum\_{m \in \mathbb{N}, w \in \Psi\_{m}^{\mathbb{S}}, t \in T} \mathsf{C}\_{wt}^{\mathbb{S}} F\_{wt}^{\mathbb{S}}.\tag{A35}$$

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
