Assessing the Effects of Uncertain Energy and Carbon Prices on the Operational Patterns and Economic Results of CHP Systems

: In the power and heat sectors, the uncertainty of energy and carbon prices plays a decisive role in the rationale for decommissioning/repurposing coal-ﬁred CHP (combined heat and power) systems and on investment decisions of energy storage units. Therefore, there is a growing need for advanced methods that incorporate the stochastic disturbances of energy and carbon emission prices into the optimization process of an energy system. In this context, this paper proposes an integrated method for investigating the effects of uncertain energy and carbon prices on the operational patterns and ﬁnancial results of CHP systems with thermal energy storage units. The approach combines mathematical programming and Monte Carlo simulation. The computational process generates feasible solutions for proﬁt maximization considering the technical constraints of the CHP system and the variation of energy and carbon emission prices. Four scenarios are established to compare the operational patterns and economic performance of a CHP system in 2020 and 2030. Results show that in 2020, there is an 80% probability that the system’s annual proﬁt will be less than or equal to € 30.98 M. However, at the same probability level, the annual proﬁt in 2030 could fall below € 11.88 M. Furthermore, the scenarios indicate that the incorporation of a thermal energy storage unit leads to higher expected proﬁts ( € 0.74 M in 2020 and € 0.71 M in 2030). This research shows that coal-ﬁred CHP plant operators will face costly risks and potentially greater challenges in the upcoming years with the increasing regulatory and ﬁnancial pressure on CO 2 emissions and the EU’s plan of phasing out fossil fuels from electricity and heat generation.


Introduction
The ability to meet the demand for electricity and heat in energy systems is an important component of energy security. Consequently, energy security has become one of the main pillars of the European Union's energy policies. However, the achievement of this policy goal must also take into account different aspects, such as local market conditions and environmental issues. As a result, mechanisms for increasing the efficiency of primary energy use have gained significant attention in recent years, especially mechanisms that support the development of highly efficient cogeneration systems. In the EU heating sector, the deployment of new combined heat and power plants (CHPs) and the modernization of existing CHP installations have become vital for maintaining energy security.
However, with the rising CO 2 emission allowance prices, the uncertain situation on the fuel market, and the high volatility of the wholesale electricity prices, CHP plant owners and operators now face major financial and operational challenges; this is particularly true for systems with low operational flexibility and powered by fossil fuels. In this regard, several European countries have implemented new financial support schemes targeted to mitigate the inherent uncertainties of contemporary energy markets. An example of such schemes is the capacity market introduced in the United Kingdom and Poland or ancillary prices, and carbon mission prices with a Monte Carlo method. More recently, the authors of Ref. [14] developed a Monte Carlo simulation and a multi-objective optimization criterion to investigate the influence of uncertainties on the optimal size and annual costs of a CHP system. The study focused on evaluating three different operational strategies of the energy system. In Ref. [15], the optimal design and operational planning of a residential and urban energy network using a Monte Carlo-based framework was investigated. The analysis explored the effects of uncertainty in heat demand with two probability distribution types. In Ref. [16], the author proposed a mixed-integer linear programming model integrated with a Monte Carlo simulation method to investigate the uncertainties of electrical and thermal demand as well as the intermittency of PV arrays and wind turbined on the operation and sensitivity of a microgrid with CHP units. In the abovementioned studies, the effects of possible future carbon emission and coal and electricity prices on the operational and financial behavior of CHP systems were ignored.
Taking into account the above-described circumstances of EU energy markets and the limited literature on the study of long-term uncertainties associated with energy and carbon prices, the main purpose of this article is to investigate the impact of uncertain parameters on the financial and operational patterns of large-scale coal-fired CHP systems coupled with thermal energy storage units. To achieve this goal, this paper develops a computational framework based on a mathematical programming method and a Monte Carlo technique. In this regard, this paper contributes to the literature by: (i) Proposing a computational approach based on mathematical programming and a Monte Carlo technique to facilitate understanding and evaluation of stochastic parameters that affect the dynamic behavior of CHP systems. (ii) Exploring the impact of observed and projected energy and EUA prices (for 2020 and 2030) on the financial and operational patterns of a stand-alone CHP system and an integrated CHP-TES system.
With this scope in mind, the remainder of this paper is organized as follows. Section 2 details the method developed to investigate the effects of stochastic energy and carbon prices on the economic performance and operational patterns of CHP systems with thermal energy storage units. Section 3 describes the case study, data, and research scenarios. Section 4 presents the results obtained from the application of the method. Section 5 concludes.

Materials and Methods
This work proposes a combined method to investigate the effects of uncertain energy and carbon prices on the optimal operation of a coal-fired CHP system with thermal energy storage. The approach combines mathematical programming and a Monte Carlo computational technique as a way to incorporate the stochastic disturbances of energy and carbon emission prices into the optimization process of the energy system. A flowchart of the proposed method is shown in Figure 1.

Stochastic Simulation
In the power and heat sectors, the uncertainty of energy and carbon prices plays a decisive role in the rationale for decommissioning/repurposing coal-fired CHP systems and on investment decisions of energy storage units. Therefore, there is a growing need for combined methods that consider the stochasticity of specific parameters and facilitate understanding the financial implications of future electricity, coal, and carbon emission prices.
In this context, the approach proposed in this work considers the uncertainty of sensitive parameters through the application of a Monte Carlo method. Monte Carlo is a stateof-the-art computational technique for investigating the effects of uncertain parameters on the behavior of energy systems [14]. Thus, this method is applied in the present study to generate scenarios, taking into consideration a set of inputs chosen from random samples drawn from independent continuous probability distributions.
Unlike previous works that account for uncertainties as simple sensitivity analysis with discrete scenarios (case-based scenarios representing possible changes (e.g., low, mid, and high), this study incorporates a significant larger proportion of scenarios with stochastic inputs contained within minimum and maximum projected values of energy and carbon emission prices. The projected prices adopted as parameters for the distributions are based on the results of global transition scenarios generated by the International Energy Agency [17] and European Scenarios published by ENTSO-E [18].
The process employed for modeling and generating the Monte Carlo scenarios follows the guidelines outlined by [14,19]. It involves multiple steps that can be summarized as follows:

•
Identification of model parameters that are considered uncertain and data collection from various information sources-the data consist of historical and projected values for the uncertain model inputs (i.e., electricity, coal, and carbon emission prices). • Data analysis and selection of probability density functions for each model input.

•
Generation of random samples from probability density functions. The random samples are then converted into a set of possible inputs or Monte Carlo scenarios (coal, electricity, and carbon emission prices), which are subsequently transferred to the deterministic model. In this study, the deterministic model is a mixed-integer linear programming model.

Stochastic Simulation
In the power and heat sectors, the uncertainty of energy and carbon prices plays a decisive role in the rationale for decommissioning/repurposing coal-fired CHP systems and on investment decisions of energy storage units. Therefore, there is a growing need for combined methods that consider the stochasticity of specific parameters and facilitate understanding the financial implications of future electricity, coal, and carbon emission prices.
In this context, the approach proposed in this work considers the uncertainty of sensitive parameters through the application of a Monte Carlo method. Monte Carlo is a state-of-the-art computational technique for investigating the effects of uncertain parameters on the behavior of energy systems [14]. Thus, this method is applied in the present study to generate scenarios, taking into consideration a set of inputs chosen from random samples drawn from independent continuous probability distributions.
Unlike previous works that account for uncertainties as simple sensitivity analysis with discrete scenarios (case-based scenarios representing possible changes (e.g., low, mid, and high), this study incorporates a significant larger proportion of scenarios with stochastic inputs contained within minimum and maximum projected values of energy and carbon emission prices. The projected prices adopted as parameters for the distributions are based on the results of global transition scenarios generated by the International Energy Agency [17] and European Scenarios published by ENTSO-E [18].
The process employed for modeling and generating the Monte Carlo scenarios follows the guidelines outlined by [14,19]. It involves multiple steps that can be summarized as follows:

•
Identification of model parameters that are considered uncertain and data collection from various information sources-the data consist of historical and projected values for the uncertain model inputs (i.e., electricity, coal, and carbon emission prices). • Data analysis and selection of probability density functions for each model input. It is worth highlighting that the approach proposed in this study is particularly valuable for CHP-plant and district heating network operators, since it enables the testing of new operation strategies and the comprehensive analysis of the system's dynamic behavior considering multiple sources of uncertainty.

Energy System Model
As described in Section 1, this work expands the scope of our previous study [20] by (a) modeling and incorporating energy and electricity prices uncertainties and (b) investigating the economic and operational impacts of future energy and carbon prices on the operation of a coal-fired CHP system with thermal energy storage. Furthermore, the model used in the present study is an adapted version of the mixed-integer linear programming approach described in [21,22]. In our previous work, the feasible operating region of the combined heat and power plant equipped with an extraction-condensing steam turbine was modeled as a combination of convex points that represented the hourly thermal and electrical power generation. Using binary variables to represent the bounded polyhedron can be a major drawback, since it considerably increases the difficulty of finding an optimal solution. Therefore, in the present work, the model has been improved in several points. First, the modeling of the feasible operating region is based on a linear formulation that does not require binary variables and leads to shorter computation times. Second, the mathematical formulation of the storage thermal energy losses to the surroundings is modeled using a constant loss factor and predefined charge/discharge limits, which considerably simplify the problem. Consequently, the proposed formulation is much more efficient and allows the model to be used in an integrated stochastic approach based on a Monte Carlo method.
In addition, modeling the uncertainty of the future energy and carbon emission prices with a Monte Carlo technique allows to investigate the dynamics of the objective function and decision variables at the temporal resolution defined in the deterministic programming model. Generally, optimization models implemented in a Monte Carlo framework are markedly reduced and simplified to overcome the computational expense encountered during the iterative process. One standard method to reduce the computational burden is to find a subset of representative operating periods (e.g., daily, weekly, seasonal). However, a drawback of this method is that the results may lose their chronology and the thermal cycling of the storage units is not well represented [23]. Consequently, this study differs from previous ones in two aspects. First, the system CHP system is optimized for one year at hourly intervals, therefore retaining the chronology. Second, the short-term thermal variation of the energy storage modeled for one year makes it possible to quantitatively measure the contribution of the TES in reducing the mismatch of energy production and demand.
For the sake of brevity, the following subsections describe only some important features of the constraints and model. The equations of the optimization model can be grouped into four categories: objective function, thermal power balance, combined heat and power plant, auxiliary boilers, and thermal energy storage. Appendix A provides a description of the sets, parameters and variables used in the model.

Objective Function
The objective function concerns the maximization of the total system profit (P r ). It is defined as the sum of the hourly revenues from heat and electricity sales (R tot ) minus the sum of hourly costs (C tot ). The objective function reflects the operation strategy of a combined heat and power system exposed to a liberalized electricity market, and that can benefit from the possibility of selling electricity to the grid. The total revenue (R tot ) consists of the income collected from electricity and heat sales. In Equation (2), p EC t indicates the electrical power output of the CHP plant, C Electricity t stands for the electricity price at hour t, Q Demand t is the heat demand at hour t, and C Heat t is the heat price at time t: Equation (3) describes the total costs of the system (C tot ), which can be decomposed into fuel costs, emission costs, variable operating costs, and start-up costs. Please note that formulation of Equation (3) excludes capital and non-variable costs, since the model is constructed and intended for the generation scheduling of a CHP system: where f EC t stands for the fuel consumption of the CHP in hour h, C EC f is the fuel cost of the CHP, q EC t is the thermal power output of the CHP, E EC p is the emission factor of pollutant p related to the power output of the CHP, and C p is the emission cost of pollutant p. Variable O&M costs related to the thermal and electrical output of the CHP are described by C EC−Q VOM and C EC−E VOM . Moreover, τ B t,b is the slack variable used for modeling the heat-only boiler fuel consumption and C B f is the fuel cost of the heat-only boiler in hour h. The emission factor of pollutant p related to the power output of the heat-only boiler (HOB) is described by E B p and the variable O&M cost of the HOB by C B VOM . Lastly, the z B t,b is the binary variable used to model the start-up trajectory of the heat-only boiler and C B SU stands for the start-up costs of the HOB.

Thermal Power Balance
Equation (4) reflects the overall thermal power balance of the system. In the model, we assume that the thermal demand must be satisfied by the input sources (i.e., CHP plant, heat-only boilers, and discharged energy from the thermal energy storage Q TES,dis t ) and balanced by the output sources (i.e., the energy directed to the thermal energy storage q TES,chr t ). Q Demand t stands for the hourly heat demand of the district heating network:

Combined Heat and Power Plant
The operation of the CHP plant is modeled by Equations (5)- (13). Because the CHP plant is assumed to be equipped with an extraction-condensing steam turbine, the feasible operation zone (FOZ) is constructed as a two-dimensional region that describes the production possibility sets p EC t and q EC t . The production possibility sets are dependent on the power-to-loss ratio (β) and power-to-heat ratio (σ) [24,25]. Equations (5)-(9) express the FOZ considering the maximum thermal and electrical power output of the steam turbine. The fuel consumption of the pulverized coal-fired boilers supplying high-pressure steam to the turbine is given by Equations (10) and (11). Although condensing-extraction steam turbines provide high levels of operational flexibility, maximum ramping rates must be considered for the safety and stability of the system, as shown in Equations (12) and (13) (13) where P EC,Max t and P EC,Min t describe the maximum and minimum electrical power output of the CHP, the parameters m EC,Max , m EC,Min , b EC, Max , and b EC, Min represent the slopes and intercepts of the linear functions adopted to model the fuel consumption of the steam boilers, and ∆R Max and ∆R Min are the maximum and minimum rate of change, respectively, of the combined power output of the CHP system.

Auxiliary Boilers
The operational strategy and technical limitations of the auxiliary boilers are described by Equations (14)- (19). The maximum power output and logical state of the boilers are given in Equations (14)- (16). The auxiliary boilers can enter operation when the thermal demand exceeds the thermal power output of the CHP system. The fuel consumption of the boilers is linearized considering the approach proposed in [26], as expressed in Equation (19): where Q B,Max is the maximum power output of the heat-only boilers, Q EC,Max is the maximum power output of the CHP system, the commitment status of the heat-only boiler is given by the binary variable u B t,b , and the fuel consumption of the boilers is described by the slack variable τ B t,b and coefficients f b B j , f m B j .

Thermal Energy Storage
Equations (20)-(27) describe the operation of the thermal energy storage in each time period. The energy balance of the thermal energy storage is given in Equations (20) and (21). The mathematical formulation of the thermal energy storage unit considers a predefined capacity Q TES Cap and C-factors C f MAX chr , C f MAX dis (i.e., upper and lower bounds of charge and discharge rates). As mentioned in Section 2, the thermal energy losses to the surroundings are modeled using constant loss factors (η s ), which simplifies the optimization problem: where the thermal energy storage charge and discharge efficiencies are described by the parameters η chr and η dis .

Case Study
The combined heat and power plant examined in this study consists of two pulverized coal-fired boilers supplying high-pressure steam to an extraction-condensing steam turbine (120 MW e and 205 MW t ). Additionally, two auxiliary coal-fired boilers (80 MW each, commonly referred to as "heat-only boilers") are installed to satisfy the thermal demand at peak load hours or during scheduled/unscheduled downtimes and unexpected breakdowns. The CHP-HOB system covers the thermal demand for space heating and domestic hot water of customers in a district heating network. A boxplot of the hourly heat load for different months of the year is depicted in Figure 2.

Thermal Energy Storage
Equations (20)-(27) describe the operation of the thermal energy storage in each time period. The energy balance of the thermal energy storage is given in Equations (20) and (21). The mathematical formulation of the thermal energy storage unit considers a predefined capacity and C-factors , (i.e., upper and lower bounds of charge and discharge rates). As mentioned in Section 2, the thermal energy losses to the surroundings are modeled using constant loss factors , which simplifies the optimization problem: where the thermal energy storage charge and discharge efficiencies are described by the parameters and .

Case Study
The combined heat and power plant examined in this study consists of two pulverized coal-fired boilers supplying high-pressure steam to an extraction-condensing steam turbine (120 MWe and 205 MWt). Additionally, two auxiliary coal-fired boilers (80 MW each, commonly referred to as "heat-only boilers") are installed to satisfy the thermal demand at peak load hours or during scheduled/unscheduled downtimes and unexpected breakdowns. The CHP-HOB system covers the thermal demand for space heating and domestic hot water of customers in a district heating network. A boxplot of the hourly heat load for different months of the year is depicted in Figure 2. Besides the heat-only boilers, the CHP plant is connected to a tank thermal energy storage unit (TES). This type of storage technology is generally designed for short-term applications. Consequently, in this study, the storage capacity is limited to one-week thermal cycles; in other words, the TES must charge and discharge to the initial conditions of the storage capacity within 168 h of operation. The total storage capacity of the tank is 400 MWh. Because the CHP plant operators aim to maximize the expected profit from the electricity traded at the day-ahead power market, the storage unit is used to avoid thermal load imbalances and acts as a thermal energy buffer. Figure 3 shows the schematic of the CHP system. Besides the heat-only boilers, the CHP plant is connected to a tank thermal energy storage unit (TES). This type of storage technology is generally designed for short-term applications. Consequently, in this study, the storage capacity is limited to one-week thermal cycles; in other words, the TES must charge and discharge to the initial conditions of the storage capacity within 168 h of operation. The total storage capacity of the tank is 400 MWh. Because the CHP plant operators aim to maximize the expected profit from the electricity traded at the day-ahead power market, the storage unit is used to avoid thermal load imbalances and acts as a thermal energy buffer. Figure 3 shows the schematic of the CHP system.

Data and Research Scenarios
Four scenarios are designed to investigate the economic effects of uncertain energy and carbon emission prices on the financial performance of a coal-fired combined heat and power system. Scenario I and Scenario II aim to illustrate the effects of the commodity prices observed in 2020 (1 January 2020 to 31 December 2020) on the dynamic performance of the CHP system configuration. Scenario I considers the case of the stand-alone CHP-HOB system, while Scenario II describes the same system integrated with a 400-MWh tank thermal energy storage unit.
As future electricity, coal, and carbon emission allowance prices are uncertain, the present study investigates two additional scenarios. Scenario III and Scenario IV explore the impact of energy and EUA prices in 2030 on the financial and operational performance of the stand-alone CHP system and integrated CHP-HOB-TES system. Possible future energy and carbon emission prices reported by the International Energy Agency [17] and ENTSO-E [18] are used to define the bound and parameters that describe the probability distributions. Table 1 summarizes the scenarios under which the CHP system is evaluated.

Data and Research Scenarios
Four scenarios are designed to investigate the economic effects of uncertain energy and carbon emission prices on the financial performance of a coal-fired combined heat and power system. Scenario I and Scenario II aim to illustrate the effects of the commodity prices observed in 2020 (1 January 2020 to 31 December 2020) on the dynamic performance of the CHP system configuration. Scenario I considers the case of the stand-alone CHP-HOB system, while Scenario II describes the same system integrated with a 400-MWh tank thermal energy storage unit.
As future electricity, coal, and carbon emission allowance prices are uncertain, the present study investigates two additional scenarios. Scenario III and Scenario IV explore the impact of energy and EUA prices in 2030 on the financial and operational performance of the stand-alone CHP system and integrated CHP-HOB-TES system. Possible future energy and carbon emission prices reported by the International Energy Agency [17] and ENTSO-E [18] are used to define the bound and parameters that describe the probability distributions. Table 1 summarizes the scenarios under which the CHP system is evaluated. In the present set of scenarios, a special case of Beta distribution was used to model the potential ranges of future energy and EUA prices. This type of continuous distributioncommonly referred to as PERT-uses three estimators: minimum, maximum, and most likely value. This distribution was selected since it has been extensively used in several fields to model projections and expectations of commodity prices [11,19]. Moreover, it constructs a smooth curve and emphasizes the most likely value, which is of outstanding importance in exploratory studies (research that accounts for possible versions of the future). The data used to estimate the distribution parameters were collected from multiple information sources, and are summarized in Table 2. Data for coal prices in 2020 were obtained from the Polish Steam Coal market index [27]. Coal price projections for 2030 were extracted from the most recent World Energy Outlook (WEO 2021) [17]. In the case of electricity day-ahead prices for 2020, the data were taken from the Polish Power Exchange and the Transmission System Operator [28]. Electricity price projections for 2030 were obtained from the most recent ENTSO-E Ten-Year Network Development Plan (TYNDP 2022) [18]. Prices of CO 2 emission allowances for 2020 were extracted from the Quandl [29], whereas price projections for 2030 were taken from the WEO 2021 [17]. It is assumed that heating prices are regulated with a single tariff for the one-year period in all scenarios. The tariff for 2020 was estimated from heat prices reported by the Energy Regulatory Office of Poland [30]. Heat prices for 2030 were projected by extending the upward trend in local municipal heat prices observed in recent years.    In the present set of scenarios, a special case of Beta distribution was used to model the potential ranges of future energy and EUA prices. This type of continuous distribution-commonly referred to as PERT-uses three estimators: minimum, maximum, and most likely value. This distribution was selected since it has been extensively used in several fields to model projections and expectations of commodity prices [11,19]. Moreover, it constructs a smooth curve and emphasizes the most likely value, which is of outstanding importance in exploratory studies (research that accounts for possible versions of the future).
The data used to estimate the distribution parameters were collected from multiple information sources, and are summarized in Table 2. Data for coal prices in 2020 were obtained from the Polish Steam Coal market index [27]. Coal price projections for 2030 were extracted from the most recent World Energy Outlook (WEO 2021) [17]. In the case of electricity day-ahead prices for 2020, the data were taken from the Polish Power Exchange and the Transmission System Operator [28]. Electricity price projections for 2030 were obtained from the most recent ENTSO-E Ten-Year Network Development Plan (TYNDP 2022) [18]. Prices of CO2 emission allowances for 2020 were extracted from the Quandl [29], whereas price projections for 2030 were taken from the WEO 2021 [17]. It is assumed that heating prices are regulated with a single tariff for the one-year period in all scenarios. The tariff for 2020 was estimated from heat prices reported by the Energy Regulatory Office of Poland [30]. Heat prices for 2030 were projected by extending the upward trend in local municipal heat prices observed in recent years.

Model Implementation
The proposed computational framework was coded in the General Algebraic Modeling Systems (GAMS) and soft-linked to MATLAB. Random samples from PERT distributions for each input parameter were generated with MATLAB using the Statistics and Machine learning toolbox. The mixed-integer linear programming model was implemented in GAMS and solved using CPLEX 20.1 in a desktop computer with 46 GB of RAM and an Intel Core i7-8086 4.0 GHz.

Result
As described in Section 2, the method involved solving the optimization model number of times until the average value and standard deviation of the expected profit stabilized. For each of the Monte Carlo scenarios, the optimization model was run for a full year. Each scenario consisted of a set of possible coal, electricity, and carbon emission prices. Figure 6 shows the results of the Monte Carlo simulation for Scenarios I and III after 1200 replications. From the figure, it can be observed that the average value and standard deviation stabilized at approximately 1000 sample sets. It is worth noting that the addition of the energy storage in Scenarios II and IV increased the solution times significantly. This was mainly due to the relatively higher number of binary variables used for modeling the behavior of the thermal energy storage unit.

Model Implementation
The proposed computational framework was coded in the General Algebraic Modeling Systems (GAMS) and soft-linked to MATLAB. Random samples from PERT distributions for each input parameter were generated with MATLAB using the Statistics and Machine learning toolbox. The mixed-integer linear programming model was implemented in GAMS and solved using CPLEX 20.1 in a desktop computer with 46 GB of RAM and an Intel Core i7-8086 4.0 GHz.

Result
As described in Section 2, the method involved solving the optimization model N number of times until the average value and standard deviation of the expected profit stabilized. For each of the Monte Carlo scenarios, the optimization model was run for a full year. Each scenario consisted of a set of possible coal, electricity, and carbon emission prices. Figure 6 shows the results of the Monte Carlo simulation for Scenarios I and III after 1200 replications. From the figure, it can be observed that the average value and standard deviation stabilized at approximately 1000 sample sets.

Model Implementation
The proposed computational framework was coded in the General Algebraic Modeling Systems (GAMS) and soft-linked to MATLAB. Random samples from PERT distributions for each input parameter were generated with MATLAB using the Statistics and Machine learning toolbox. The mixed-integer linear programming model was implemented in GAMS and solved using CPLEX 20.1 in a desktop computer with 46 GB of RAM and an Intel Core i7-8086 4.0 GHz.

Result
As described in Section 2, the method involved solving the optimization model number of times until the average value and standard deviation of the expected profit stabilized. For each of the Monte Carlo scenarios, the optimization model was run for a full year. Each scenario consisted of a set of possible coal, electricity, and carbon emission prices. Figure 6 shows the results of the Monte Carlo simulation for Scenarios I and III after 1200 replications. From the figure, it can be observed that the average value and standard deviation stabilized at approximately 1000 sample sets. It is worth noting that the addition of the energy storage in Scenarios II and IV increased the solution times significantly. This was mainly due to the relatively higher number of binary variables used for modeling the behavior of the thermal energy storage unit. It is worth noting that the addition of the energy storage in Scenarios II and IV increased the solution times significantly. This was mainly due to the relatively higher number of binary variables used for modeling the behavior of the thermal energy storage unit. Table 3 summarizes the results obtained after 1200 iterations. Additionally, it shows the fluctuations that occurred between the first two hundred iterations and the evolution of the mean profit as the total number of iterations increases. Note: SD stands for standard deviation; t represents computation time.
The iterative computational process generates feasible solutions for profit maximization considering the technical constraints of the CHP system and the variation of energy and carbon emission prices. The solutions can be described in the form of histograms and used to assess the variation in profit and the financial contribution of the thermal energy storage unit. Figure 7 provides a visual comparison of the four scenarios examined in this study. The results showed that integrating a short-term thermal energy storage unit increased the profitability of the system and helped reduce the risk associated with fluctuating energy prices. This can be observed in Figure 7b Table 3 summarizes the results obtained after 1200 iterations. Additionally, it shows the fluctuations that occurred between the first two hundred iterations and the evolution of the mean profit as the total number of iterations increases. The iterative computational process generates feasible solutions for profit maximization considering the technical constraints of the CHP system and the variation of energy and carbon emission prices. The solutions can be described in the form of histograms and used to assess the variation in profit and the financial contribution of the thermal energy storage unit. Figure 7 provides a visual comparison of the four scenarios examined in this study. The results showed that integrating a short-term thermal energy storage unit increased the profitability of the system and helped reduce the risk associated with fluctuating energy prices. This can be observed in Figure 7b Furthermore, based on the outcomes of Scenario II (2020, CHP-HOB-TES) and Scenario IV (2030, CHP-HOB-TES), it can be noticed that in 2030 coal-powered CHP systems will face the risk of very low returns. This risk is mainly triggered by the continued upward movement in EUA prices.
The distributions of the different outcomes also allowed to estimate the probability of the potential annual profits. For the case of the stand-alone CHP system operating in Furthermore, based on the outcomes of Scenario II (2020, CHP-HOB-TES) and Scenario IV (2030, CHP-HOB-TES), it can be noticed that in 2030 coal-powered CHP systems will face the risk of very low returns. This risk is mainly triggered by the continued upward movement in EUA prices.
The distributions of the different outcomes also allowed to estimate the probability of the potential annual profits. For the case of the stand-alone CHP system operating in the market scenario of 2020, the results showed that there is an 80% probability that the annual profit will be less than or equal to €30.98 M. On the other hand, with the installation of a tank thermal energy storage unit, the cumulative probability of 80% was at €31.72 M. Based on the cumulative distribution functions (CDFs) of Scenarios III and IV (2030), it can be stated that there is a 95% probability that the annual profit of the stand-alone and the integrated CHP-HOB-TES system will be below €14.64 M. Furthermore, the analysis of these two scenarios showed that the thermal energy storage increases the chances of receiving additional profits. There is an 80% probability that in 2030, the annual profit of the stand-alone CHP system will be less than or equal to €11.88 M, while for the integrated CHP-HOB-TES system, the profit may be less than or equal to €12.6 M. The findings above are particularly important for potential investors in new cogeneration systems and thermal energy storage units, since they offer valuable insights into the economic consequences of integrating the two technologies.
During each iteration, the optimization model solves the coal-fired CHP system's operational planning problem, taking into account the scenarios drawn by the Monte Carlo procedure. This computational process allows monitoring and collecting information about the system's economic performance in each hour of the simulated year. Figure 8 illustrates the hourly generation costs and revenues from electricity sales for one week in 2030. The stochastic simulated time series capture the variability in generation costs and revenues from electricity sales of the stand-alone and the integrated CHP system. The large variation envelope in generation costs indicated that coal and carbon emission prices have a more significant impact on the optimal behavior of the system as compared to the variation in electricity prices. The average annual results obtained from the Monte Carlo procedure for 2020 and 2030 are illustrated in Figure 9. The figure breaks down the estimates into three separate components: revenues, costs, and profits. The total annual revenue of the stand-alone CHP increased from nearly €72 M in 2020 to €102 M in 2030, representing a rise in revenue of approximately 40%. However, because of the projected increase in carbon emission allowances prices, the average annual generation costs for the same CHP system configuration nearly doubled. Despite the higher revenues in 2030, the substantial increase in generation costs resulted in a drop in expected profits from €24 M to €9.15 M, or approxi- The average annual results obtained from the Monte Carlo procedure for 2020 and 2030 are illustrated in Figure 9. The figure breaks down the estimates into three separate components: revenues, costs, and profits. The total annual revenue of the stand-alone CHP increased from nearly €72 M in 2020 to €102 M in 2030, representing a rise in revenue of approximately 40%. However, because of the projected increase in carbon emission allowances prices, the average annual generation costs for the same CHP system configuration nearly doubled. Despite the higher revenues in 2030, the substantial increase in generation costs resulted in a drop in expected profits from €24 M to €9.15 M, or approximately 62%. Similar variations were observed for the scenarios that incorporate a tank thermal energy storage unit. These findings indicate that coal-fired CHP plant operators will face costly risks and potentially greater challenges in the upcoming years with the increasing regulatory and financial pressure on CO 2 emissions and the EU's plan of phasing out coal and other fossil fuels from electricity and heat generation. The average annual results obtained from the Monte Carlo procedure for 2020 and 2030 are illustrated in Figure 9. The figure breaks down the estimates into three separate components: revenues, costs, and profits. The total annual revenue of the stand-alone CHP increased from nearly €72 M in 2020 to €102 M in 2030, representing a rise in revenue of approximately 40%. However, because of the projected increase in carbon emission allowances prices, the average annual generation costs for the same CHP system configuration nearly doubled. Despite the higher revenues in 2030, the substantial increase in generation costs resulted in a drop in expected profits from €24 M to €9.15 M, or approximately 62%. Similar variations were observed for the scenarios that incorporate a tank thermal energy storage unit. These findings indicate that coal-fired CHP plant operators will face costly risks and potentially greater challenges in the upcoming years with the increasing regulatory and financial pressure on CO2 emissions and the EU's plan of phasing out coal and other fossil fuels from electricity and heat generation.

Conclusions
The main goal of the study was to investigate the effects of uncertain energy and carbon prices on the operation and financial performance of CHP systems with thermal energy storage. This objective was achieved by developing a stochastic approach composed of a mathematical programming method and a Monte Carlo technique. The approach was designed to deal with the uncertainty of fluctuating energy and carbon prices and assess the financial contribution of thermal energy storage. The proposed computational framework was coded in the General Algebraic Modeling Systems (GAMS) and soft-linked to MATLAB. The stochastic approach was applied to generate scenarios taking into consideration a set of inputs chosen from random samples drawn from independent continuous probability distributions.
In the study, four scenarios were investigated. Scenario I and Scenario II aimed to illustrate the effects of the commodity prices observed in 2020. Scenario III and Scenario IV explored the impact of energy and EUA prices in 2030. From the results, the following conclusions can be derived:

1.
The iterative computational process generates feasible solutions for profit maximization considering the technical constraints of the CHP system and the variation of energy and carbon emission prices.

2.
The distributions of the different outcomes allowed to estimate the probability of the potential annual profits. For the case of the stand-alone CHP system operating in the market scenario of 2020, there is an 80% probability that the annual profit will be less than or equal to €30.98 M.

3.
There is an 80% probability that in 2030 the annual profit of the stand-alone CHP system will be less than or equal to €11.88 M, while for the integrated CHP-HOB-TES system, the profit may be less than or equal to €12.6 M.

4.
Integrating a short-term thermal energy storage unit increased the profitability of the system and helped reduce the risk associated with fluctuating energy prices. Profit increased on average by €0.74 M (with the implementation of a TES) in 2020 and €0.71 M in 2030.
In the coming years, the operational patterns and economic results of CHP systems will be significantly affected by new electricity and heat consumption patterns and market changes. Moreover, further research challenges will arise because of the increasing penetration of renewable generation and large-scale electrical energy storage deployment. The issues mentioned above will require comprehensive models that consider multiple interdependent sources of uncertainty (e.g., short-term economic factors, environmental and operational aspects of renewable power technologies, power and heat consumption patterns, and thermal comfort levels, among others). In this regard, an important avenue for future research is the incorporation of wind and solar systems along with their high degree of uncertainty (wind speed and solar irradiation) into the proposed Monte Carlo-based method. Another direction for future research is the integration of new computational techniques such as deep learning (neural networks) to reduce the computational complexity of the Monte Carlo approach and the mixed-integer linear programming model. Author Contributions: Conceptualization, P.B., P.K. and J.K.; data curation, P.B.; methodology, P.B. and P.K.; software, P.B.; validation, P.B., P.K. and J.K.; formal analysis, P.B. and P.K.; writing-original draft preparation, P.B. and P.K.; writing-review and editing, P.B. and J.K.; visualization, P.B. and P.K.; supervision, J.K. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.  Table A2. Parameters of the mathematical model.  Slack variable-heat-only boiler fuel consumption, (Mg)

Variables Description
Binary variables Commitment status of heat-only boiler Start-up of heat-only boiler