Multi-Energy System Optimization of Costs Versus Carbon Dioxide Emissions for Flexibility. A Case Study in Italy

1  Introduction

Current energy systems are increasingly complex due to the need to manage multiple energy vectors, the growing integration of non-programmable renewable energy sources, and the increased focus on reducing carbon dioxide emissions (CO2) while still trying to reduce costs. This complexity is the opportunity explored in this paper, with the additional idea of handling the complexity systematically and automatically.

The paper focuses on implementing a multi-energy system (MES) for managing the production of electricity, steam, and cold, integrated with photovoltaic (PV), a battery energy storage system (BESS), cogenerators (CCHP), and boilers. Since this research is part of the results of the FLEXIndustries Project, from which we obtain the requirements and case studies, the evaluation was carried out in a pharmaceutical industry established in the North of Italy. The FlexIndustries project aims to contribute to the transition towards more efficient energy sources and improved production processes to enable energy flexibility.

This paper proposes a multi-objective and multi-energy optimization system that implements a Linear Programming algorithm to optimize costs vs carbon dioxide emissions on an hourly basis over three months. The evaluated industry provides historical data on electric energy and steam demand, while data on prices (electric energy and natural gas) are taken from historical hourly energy prices in Italy. The implementation enables automatic optimization of the defined multi-energy system, providing energy supply alternatives for evaluating cost reductions or reductions in carbon dioxide emissions. The optimization was conducted using a Linear Programming (LP) formulation to achieve fast, reliable solution times and strong optimality guarantees, which are essential in such decision-support systems. Nevertheless, we recognize that a more detailed MES formulation may introduce nonlinearities and nonconvexities, and metaheuristic algorithms with adaptive parameter-control approaches may be needed to handle these extended problem variants.

The results show that the optimization system allows a Pareto front comparison to choose between reducing costs and reducing carbon dioxide emissions, with the values expressed as percentage differences from the actual values. The cost differences between the optimized solutions and the real cost at market prices range from −9% to 42%. Then, the differences between the carbon dioxide emissions of the optimized solutions and the actual carbon dioxide emissions range from −61% to 16%. Also, the paper examined in greater depth the flexibility enabled by the energy system’s storage and cogeneration systems, as well as the system’s overall management.

2  State of the Art

Related research has demonstrated interest in addressing multi-energy systems to optimize processes and reduce costs, consumption, and CO2 emissions. However, existing literature typically focuses on a single technology, energy vector, or evaluation objective for modelling energy systems.

The closest research publications on multi-objective optimization are implemented for an Energy Community in Italy [1] and another for a public building in Spain [2]. The former research [1] did not consider BESS technologies, and the data were not obtained directly from the buildings; part of their simulation data is based on research in the literature, technical reports, webpages, and contact with the building administration. The latter study [2] focuses on a yearly optimization to compare different configurations of the polygeneration system for building energy system retrofit.

Therefore, to the best of the author’s knowledge, we did not find a MES optimization framework that simultaneously addresses the trade-off between energy procurement cost and carbon dioxide emissions with an integrated energy system comprising PV, CCHP, boilers, and BESS, using as energy vectors electricity, heat, and cold. Additionally, this research collaborates on the implementation of a holistic decision support system, using the optimization system with real data, the energy demands from the company data collections, and the energy prices from historical data of the current Italian energy market. Our approach to support energy production with an optimization system gives this research the unique objective of providing energy flexibility information to the energy managers.

The following subsections briefly elaborate on the state-of-the-art in the corresponding research areas of this paper, namely MES Process Optimization and the optimization algorithm used to address it. The Multi-Energy System Process Optimization subsection shows a literature review of MES in buildings, an optimal design approach, a model for investment planning, and optimization concentrating on different technologies. The linear programming subsection presents related research that implements similar energy systems, describes the pros and cons of LP, and opens the discussion to possible extensions of this research.

2.1 Multi-Energy System Process Optimization

The simultaneous production of different energy vectors in multi-energy systems is a promising approach for enhancing energy efficiency and facilitating the development of distributed energy systems. Also, system optimization is a viable way to exploit polygeneration systems [3]. This opens the opportunity for this research paper by focusing on operating and coordinating multiple technologies for the management of Multi-Energy Systems, opening the model analysis to different sections of the energy systems, namely, energy production, energy consumption, energy storage, and energy transportation. In recent years, considerable effort has been invested in aggregating and integrating energy vectors, including Combined Heat and Power (CHP) generation, photovoltaic systems, absorption chillers, and other technologies. This effort demonstrates an intention toward energy efficiency and enables possibilities for energy flexibility. For instance, a multiobjective (economic and environmental criteria) optimization model for the optimal design of a methanol/electricity polygeneration plant has been proposed in [4].

Energy system optimization algorithms can fall into two broad categories: model-based and data-driven. The traditional model-based approaches rely on precise representations of the energy system’s behavior using physics or engineering models, along with the corresponding technical limits. On the other hand, data-driven approaches rely on the industry’s 4.0 data connectivity in the energy system to enable dynamic, automated evaluations. In both cases, a model is needed to represent the energy system, where traditional approaches require precise information about machines and production capacities, and data-driven approaches could use data streams to learn the behavior of machines and energy systems, enhancing operational efficiency and reducing costs. However, to achieve data-driven advantages, the implementation has its additional costs, and the operation is resource-intensive (sensors, technology services, etc.).

The increasing need for MES creates new opportunities for energy system optimization algorithms to address the rising complexity. Different solutions can target varying levels of complexity, as shown in [5], ranging from a highly simplified case with simple operation optimization, storage, constant efficiency, and no technical limitations, to models that incorporate uncertainties and flexibility measures. An essential consideration in evaluating the level of complexity is that increasing the model’s complexity also increases the level of detail in the shared data. Thus, given this critical aspect, sharing details of the company’s business processes is sometimes limited to avoid disclosing sensitive information. Furthermore, multiobjective optimization is needed when considering the increasing control over carbon dioxide emissions, and power generation is a colossal source of global greenhouse gas (GHG) emissions, being the cause for the release of more than 7.7 billion tons of carbon dioxide (CO2) annually, accounting for 37.5% of the total annual carbon dioxide emissions [5].

Moreover, in [6], MES is analyzed and optimized with a combined cooling-heat-power (CCHP) system with a gas turbine and an air-source heat pump, highlighting the trade-offs among operating cost, energy efficiency, and load matching under realistic dynamic demand patterns. In [7], a set of smart railway stations was modeled as networked microgrids that can exchange energy to minimize power drawn from the utility grid, explicitly accounting for PV generation, on-site energy storage systems (ESS), and regenerative braking energy from trains. The work addresses uncertainty through a scenario-based formulation and formulates scheduling as a Mixed-Integer Linear Programming (MILP), showing that coordinated operation can reduce both grid imports and station operating costs. Finally, in [8], multi-energy configurations are increasingly moving towards emerging technologies. Specifically, an integrated energy system (IES) coupling wind power, photovoltaics, and biogas generation with carbon capture and storage (CCS) and power-to-gas (P2G) was proposed and solved using CPLEX, demonstrating the significance of a collaborative scheduling model for enhancing local energy utilization and low-carbon operations.

2.2 Linear Programming

The short computational times characteristic of LP problems can be particularly advantageous when the optimization process needs to be executed frequently, as in approximating the Pareto curve in multi-objective optimization [9]. LP problems are relaxed representations of MILP problems, in which the constraints are not integer-valued.

MILP is a state-of-the-art approach to tackling the optimization problem of polygeneration systems. The guarantee of global optimality for linear problems and the effectiveness of available commercial solvers make MILP very attractive and widely used for optimization problems in multi-energy systems [10].

MILP formulations have some drawbacks, such as the inability to model nonlinear effects, the need to consider all periods simultaneously, and the risk of high-dimensional problems [10]. However, many techniques have been proposed to address such limitations, such as piecewise linearization methods [11], rolling horizon approaches [12], and dimensionality reduction [13].

When analyzing energy systems optimization, a common characterization is the difference between solving the synthesis and design problem and adding the operation problem. For instance, in [14] they showed that when the optimization target is synthesis, design, and operation, heuristic methods may not be very efficient if used alone in comparison with deterministic approaches. Building on that insight, in [15] the authors applied linear programming to design optimal sizing and operational strategies for an energy system using trigeneration. Poland also evaluated the distributed energy system to optimize energy device handling efficiency using linear programming methods, confirming that an integrated synthesis-operation framework improves overall system performance.

A case study analysis investigated the multi-objective optimization model for the integration of polygeneration systems in energy communities, developing an MILP algorithm as the proposed solution in [1]. Recent works, such as [14], used MILP models for multi-objective optimization in Energy Community scenarios involving Heating and Cooling Networks. This analysis aimed to reduce carbon emissions while accounting for economic consequences. Similarly, in [2] an MILP model for optimizing a polygeneration system that considers real building demand, battery capabilities, and multi-objective considerations (cost and CO2 emissions) was proposed. The research summarizes recent studies on polygeneration systems for buildings, second-life batteries in stationary applications, and building retrofitting to enhance energy efficiency.

Also, a multiobjective optimization model using fuzzy linear programming (FLP) has been developed to seamlessly integrate a biomass-based polygeneration system (BBPS) in an iron plant while obtaining negative carbon emissions [16]. The study highlighted the influence of biomass prices and the trade-off between costs and carbon emissions.

3  Energy Systems Optimization

The methodology follows a comparative optimization framework that evaluates the potential of an integrated energy system by contrasting historical performance against an optimized scenario. The input data are real-world historical consumption for electricity, steam, and cooling, photovoltaics production, and Italian day-ahead market energy prices. These inputs are fed into the energy model, which simulates the conversion of energy among the PV, CCHP, BESS, and boilers. Using a multi-objective optimization algorithm, the system performs a dispatch analysis, determining the hourly operational schedules to minimize both total costs and CO2 emissions. Finally, the results are visualized as an aggregated Pareto Front, allowing for a direct comparison between the actual historical baseline and the theoretical maximum benefits achievable through the energy system optimization.

The energy system design uses as primary energy inputs electricity and natural gas. To support the production processes, these inputs are converted into secondary energy vectors, including head carriers (steam and hot water) and cold carriers (cold water). Additionally, but of greater strategic importance, is the commitment to integrate the generation of energy carriers from non-programmable renewable sources to support decarbonization, such as photovoltaics.

Another critical design factor is energy flexibility, the ability of energy systems to produce or consume energy using more than one alternative, thereby enabling adjustments to production and consumption based on supply and demand. This research focuses on energy production flexibility. The final energy system design layout is shown in Fig. 1. The external energy resources are electricity (EE) and natural gas (NG) network providers. Photovoltaics (PV), Trigenerators (CCHP), and Boilers are the main energy production utilities. The only Storage utility is the Battery Storage System (BESS). The Demands are Electric Energy (EE), Steam (ST), and Cold (CO).

Figure 1: General energy system model for process optimization.

The energy system considers the design of two utilities that allow energy flexibility: the BESS and the CCHP. The BESS allows time-shifting of energy production (from low-cost to higher-cost periods). The CCHP allows for choosing EE production over importing from the electricity network and ST production over producing from the Boilers. For the implementation of such systems, it is essential to establish an appropriate level of automation for the data required for these data-driven optimization processes. This is important for the correct deployment of the energy system model within the project; however, it is less relevant to this paper’s focus. Also, the general energy design allows multipliers (N1, N2, N3, N4, N5, and N6) to adjust the configuration to different scenarios. The only ones that should be higher than zero are N4 and N5. This generic energy system includes services to manage configurations, inputs, and outputs to ensure correct implementation and integration within the FlexIndustries Project.

3.1 Case Study

The case study is a company that produces and distributes active pharmaceutical ingredients (APIs) for human and animal health, based in Rovereto, Italy. The company focuses on production based on chemical synthesis or fermentation processes, which require a reliable supply of electricity, steam, and cold.

The energy system utilities and energy networks have the following details:

•   EE_Net: one connection for buying and selling electricity. The network limit is set equal to 20 MW.

•   NG_Net: one connection for buying natural gas.

•   PV System: one PV system of 500 kWp that produces about 550 MWh/y of electricity.

•   CCHP: two CCHPs with an electricity production capacity of 4.5 MWel each.

•   Boiler: three Boilers with a total Steam production capacity of 40 t/h. The Boilers’ natural gas flow rate is 3000 Sm3/h.

•   BESS System: one BESS system with a storage capacity of 0.4 MWh and a C_Rate of 0.66.

The data analysis was conducted in July, August, and September 2024. The energy demands used in the simulation are real hourly historical data for the utilities’ electricity and steam usage. The wholesale price from the day-ahead market of GME (https://www.mercatoelettrico.org/it-it/) was used. Specifically, historical data for electricity and natural gas prices for the Northern Italy market zone were used in this study. The simulation of the optimization algorithm used real values for electricity and steam, while, for simplicity, cold values were assigned a value of zero. These are important assumptions because the forecasting was considered perfect, which does not apply to real scenarios; therefore, these results should be considered ideal.

3.2 Implementation

We implemented the energy system in the case study using a Linear Programming algorithm. The optimization objectives are two: Costs (TotalCosts Eq. (1)) and CO2 Emissions (TotalCO2eq Eq. (2)). The default optimization horizon is 24 h, and the optimization step is hours. This representation allows different optimization horizons (from 1 to 720) to allow a generic design, where the minimum optimization is one step, and the maximum is 30 days (720 h). The proposed optimization framework solves sequentially two optimization problems characterized by two distinct objective functions: Costs (TotalCosts Eq. (1)) and CO2 Emissions (TotalCO2eq Eq. (2)). The time step of the optimization framework is 1 h. The optimization horizon can be customized from 1 h to 30 days (720 h). Given the previously mentioned assumptions on the perfectly forecasted input data, the optimization horizon is set to 24 h.

We also implement additional services to manage the coefficients used to manage the linear representation. These services help to simplify the deployment in the real scenario and enable the possibility of implementing external services to update the coefficients. The following are the equations for the optimization objectives, Costs, and CO2 Emissions over the optimization horizon N∈{1,…,720}.

Formulation of Eq. (1), minimize:

TotalCosts=∑i=0N{EEImportCosts,i+EEExportCosts,i+NGImportCosts,i}(1)

Subject to Operational Constraints and Emission Constraints.

Formulation of Eq. (2), minimize:

TotalCO2eq=∑i=0N{EENeti∗εEENet+EEPVi∗εPV+NGImporti∗εNG+EEBessi∗εBESS}(2)

Subject to Operational Constraints.

Auxiliary functions Eqs. (3)–(6) extend the explanation for Eq. (1):

XNG_Net_INi=(XNG_Trigeni+XNG_Boileri)(3)

EEImportCosts,i=XEE_Net_INi∗cEE_ImportCosti(4)

EEExportCosts,i=XEE_Net_OUTi∗cEE_ExportCosti(5)

NGImportCosts,i=XNG_Net_INi∗cNG_ImportCosti(6)

With respect to the decision variables X. The c variables are cost-related parameters, and the script epsilon variables are emission-factor parameters. The Emission Constraint refers to the emission limit established by the epsilon constraint algorithm. The following Tables 1–4 describe the Operational Constraints implemented in the optimization model.

Some auxiliary variables and equations are shown in Table 1, and the mathematical limits of the main variables are shown in Table 2. The optimization constraints are shown in Table 3. The coefficients used for the optimization to maintain the linear representation are shown in Table 4. For the management of the coefficients in Table 4, we made them configurable with software services to allow adjusting these values upon request, leaving this estimation to the energy manager. For the current study, the average over the historical dataset is adopted and fixed.

This optimization representation can be implemented with different available libraries and programming languages, and we decided to use the Pyomo library of Python. This library has an optimization algorithm based on [17], and only supports mono-objective optimization. Therefore, the multi-objective requirement was implemented using the Epsilon Constraint method. This algorithm uses the other objective functions as constraints of the optimization algorithms. The selected solver was the GNU Linear Programming Kit (GLPK).

The algorithm implementation of the epsilon constraint consists of calculating the extreme optimizations by considering only costs and only emissions. The first step is to calculate the mono-objective optimization for cost, where the emissions objective function is deactivated but calculated. The second step is to calculate the mono-objective optimization for emissions, where the cost function is deactivated from the objective function but calculated (thus considered only as an auxiliary function). Then, using these two extremes for emission values, we create N evenly spaced numbers between the calculated values of emissions (where N is the configured pareto size). Later we use the list of evenly spaced numbers of emissions as epsilon values in the Emission Constraints of formulation (1), optimizing for cost with the epsilon assigned, building the pareto front. Later, if needed, we have a configuration of the position in the pareto front that will determine the actual solution considered from the pareto front.

The Linear Programming implementation follows the design presented in Fig. 1 and implements the following definitions:

•   Decision Variables: XEE_Net_OUT,XEE_Net_IN,XNG_Trigen,XNG_Boiler,XEE_BESS_IN,XEE_BESS_OUT.

•   Objectives functions: TotalCosts (Euros) and TotalCO2eq (tons of CO2).

•   Demand Electricity Constraints: Since the electrical network allows In and Out, the constraints are equality constraints. They are about the dv_length constraints.

•   Demand Steam Constraints: Configurations to establish production equal to or greater than the demand. They are about the dv_length constraints.

•   BESS Constraints: BESS constraints are about four sets of constraints (top capacity, discharge limit, equilibrium, and EE_BESS_IN inside EE production) multiplied by dv_length.

•   Demand Cold Constraints: This set of constraints is related to the Cooling produced by the CCHPs concerning the CO demand. These configurations establish production equal to or greater than demand.

4  Results

The simulation for the optimization considers all the days of the three analyzed months, calculating 11 solutions for the pareto front. Then, we average the Costs and the CO2 emissions for all the days and solution types (optimal toward cost or CO2 emissions).

The analysis shows a pareto front (cost (Euros) vs. CO2 emissions (tons of CO2)) with multiple solutions for every analyzed day, as shown in Fig. 2. The pareto front of that day showed solutions that vary by about 11,000 Euros between the extreme solutions. Also, the CO2 emissions vary by about 85 t of CO2. These results are not normalized, but we implemented a post-processing configuration to change the values to percentages in relation to the maximum value.

Figure 2: Pareto front’s example for a cost-emission optimization of 24 h.

Analyzing a single optimization day in September, the dispatch profiles of the energy utilities supplying electricity demand are shown in Figs. 3 and 4, while those supplying steam demand are reported in Figs. 5 and 6. Two solutions are compared: the cost-optimal solution (Solution-10, Figs. 3 and 5) and the balanced solution (Solution-05, Figs. 4 and 6). The different operational strategies adopted by the two solutions can be interpreted by referring to the market price signals of the respective sectors, shown as red dotted lines with the corresponding scale on the right-hand side of each figure.

Figure 3: Electricity (EE) comparison between estimation and production of a solution for optimal Cost (Solution-10).

Figure 4: Electricity (EE) comparison between estimation and production of a balanced solution (Solution-05).

Figure 5: Steam (ST) comparison between estimation and production of a solution for optimal Cost (Solution-10).

Figure 6: Steam (ST) comparison between estimation and production of a balanced solution (Solution-05).

In the electricity sector, it can be observed that the balanced solution, which is more conservative in terms of emissions, strategically increases electricity imports from the grid during the most economically convenient hours. In the case analyzed, the central hours of the day are characterized by very low electricity prices, which induces both solutions to reduce the output of the CCHP unit. However, in the balanced solution (Solution-05), the nighttime hours are further exploited to increase grid imports. In fact, the site benefits from an electricity supply characterized by a low emission factor. Due to the coupling between the electricity and heat sectors, the partial load operation of the CCHP results in increased production in the coupled thermal sector.

Regarding the execution time, a Pareto front result of size 11, using the LP implementation, runs in under 10 s on a personal computer (without considering network transmission time).

The electricity production shown in Fig. 3 illustrates the five charge-discharge cycles of the BESS to exploit electricity costs and the distribution of renewable energy production throughout the day.

Finally, the Cold energy vector was designed in the energy system to have a more realistic view of the operational situation for the utilities. But, for simplicity, we assume the cold demand values are zero.

5  Discussion

Energy flexibility is the ability of energy systems to produce or consume energy using multiple alternatives, enabling adjustments to production and consumption in response to supply and demand. Renewable energy benefits from flexible energy systems. Shifting loads from one period of time to another is the dominant flexibility type in 60% of applications [18].

The energy system optimization implemented considers energy production shifting to provide flexibility through two energy devices: BESS and CCHP. BESS allows the storage of electricity during low-cost periods and the discharge of electricity during high-cost periods, enabling energy production and consumption to shift. As a non-variable generator, a CCHP system enables high-efficiency electricity production while simultaneously integrating heating and cooling. Combining multiple energy vectors, electricity, steam, and cold, offers operational flexibility to reduce reliance on grid electricity and steam from boilers. This allows the system to optimize energy use by choosing when to generate electricity, steam, and cooling from natural gas instead of consuming them directly from other sources.

The solutions are described as the amount of primary energy input (electricity and natural gas) used by the utilities. However, as well as Figs. 7–9, the optimization has a service to post-process the solution into percentage utilization of the utilities, based on their capacities. This is important to mention because the way we suggest visualizing the optimized solutions by the energy managers is by percentages of usage, facilitating the data visualization in the real environment. We also take this percentage visualization to the figures discussed in this section, to facilitate the visualization of the available flexibility in terms of percentages of cost and CO2 emissions.

Figure 7: Aggregated results for July 2024 grouped by solutions.

Figure 8: Aggregated results for August 2024 grouped by solutions.

Figure 9: Aggregated results for September 2024 grouped by solutions.

We aggregated the solutions based on their inclination towards an objective, calculating the average values of Cost and CO2 emissions for the solutions from three months: July, August, and September (Figs. 7–9, respectively). The inclination towards an objective is described with Solution numbers 0 to 10, where 0 is fully tilted to CO2 emissions, and 10 is fully tilted to Cost. The Figs. 7–9 are expressed as percentages, representing the ratio between the optimized values obtained from the linear programming (LP) model and the actual cost and CO2 emissions, derived from the actual scheduling of utilities, with the same price cost calculation method for electricity and natural gas. It is important to highlight that the origin points of Figs. 7–9 are the baseline values for each month, meaning the average monthly actual cost and CO2 emission. Presenting the results as percentages allows us to clearly illustrate the extent of potential improvements, showing how much costs and emissions could be reduced through the proposed optimization algorithm compared to real operational data.

Importantly, the solutions also capture the flexibility of industrial players responding to market-driven signals, enabling them to identify and capitalize on opportunities arising in the energy market while simultaneously meeting the energy requirements necessary to satisfy their production objectives. This capability highlights the practical value of the optimization approach, as it reflects both economic benefits and operational feasibility within dynamic market conditions.

The results for July allow a shift in Cost between −8% and 35%, and a shift in CO2 emissions between 16% and −52%. The results for August allow a change in Cost between −9% and 37%, and a shift in CO2 emissions between 14% and −60%. The results for September allow a change in Cost between −6% and 42%, and a shift in CO2 emissions between 10% and −61%. Finally, it is essential to note that the Pareto front passes lower than the origin point, which is the actual cost and CO2 emission, which could be explained by energy production plans that did not take advantage of production capabilities, or the energy system model did not consider relevant production limitations.

Moreover, the Figs. 7–9 show about 2%–5% distance between the approximated pareto front and the origin of the axes. It could be that the actual solution is a suboptimal production decision because of a lack of consideration of cost signals or production capabilities. Another hypothesis could be that the proposed model has a modelling error of about 2–5%, based on relevant operational constraints not considered in the energy system model. However, to better understand this situation, future work is needed on the implementations, analysis, and testing of the actual production planning and the operational constraints.

Forecasting demand, renewable energy production, and costs would likely lower optimization quality due to the well-known challenges of managing production uncertainties. For example, in real-world deployments, forecast errors can reduce realized gains, primarily by shifting or weakening price-arbitrage opportunities when price peaks are mispredicted, and by causing over- or under-charging decisions when load and PV are misestimated. To contextualize the perfect-foresight results without diminishing the primary results and purpose of this paper, the forecasting models used in the FlexIndustries project, trained on real operational data, yield errors of approximately 6%–24% MAPE (Mean Absolute Percentage Error) depending on the stream for demand, 6.6%–21.9% nMAE (normalized Mean Absolute Error) for PV with weather inputs consistently improving accuracy in various locations, and 3.6%–7.0% nMAE across many countries for prices. In practice, the impact of forecast uncertainty can be mitigated by executing the MES in a rolling-horizon manner with periodic re-optimization using updated measurements and refreshed forecasts, which limits performance degradation relative to the ideal Pareto front. A complete analysis of MES, along with the forecasting models, will be addressed in an extended version of this paper.

Furthermore, energy flexibility can be improved in the following cases:

1.   If we add a third CCHP, the electrical production of the CCHPs can be higher than the electrical demand. Then, the optimization has room to choose between generating electricity and Steam with CCHP or importing electricity from the network and using the Boilers for steam generation, thereby expanding the potential flexibility offered to the grid while sustaining efficient generation.

2.   Using the BESS, it is possible to choose whether to store cheap energy or to discharge it when it is expensive. Increasing the capacity of the BESS, along with a variable Renewable Energy Source (RES), such as photovoltaics, can improve the cost and emission savings.

The model developed, and the scenario analysis can suggest technological upgrades to the production utilities to achieve greater energy flexibility, with better cost and CO2 emissions control. The multiobjective optimization enabled precise control and informed decision-making regarding CO2 emissions, representing a significant contribution to enhancing industrial participation in energy markets and supporting broader decarbonization goals.

6  Conclusion

Process and Energy planning is essential to significantly improve the efficiency, sustainability, and overall performance of various industries. To ensure proper planning, it is vital to understand the level of systematization of operational data for energy utilities within the energy system. This includes installing sensors and managing data with an Energy Management System (EMS) to enable data-driven analysis.

This paper presents the design and development of an energy model that accounts for CCHPs, PV, Boilers, and BESS. The model considers the usage of Electrical and Natural Gas Networks and includes the benefits of integrating and managing Renewable Energy Sources.

The contributions of this paper are as follows:

•   Design and development of an LP-based multi-energy and multi-objective energy system that accounts for integrated operations of CCHPs, PV, Boilers, and BESS.

•   It demonstrates, in a real-world use case, the advantage of energy flexibility, showing results that reduce costs considering environmental impact (CO2 emissions).

•   It demonstrates the advantages of the capacity enhancements for some energy production technologies (BESS and CCHP), increasing energy flexibility.

The results presented in this paper compare the actual historical baseline with the theoretical maximum benefits achievable through the proposed energy system optimization. This gives a good idea of the potential benefits of implementing this type of energy system optimization in real scenarios.

Additionally, the flexibility analysis shows that energy utilities such as BESS and CCHPs enabled the implementation and evaluation of energy flexibility, helping energy managers better understand improvements to the energy system and possibly impacting demand planning.

If we want to use this paper’s implementation in a stochastic approach, considering PV production forecast, demand forecast (for EE, ST, and CO), and energy prices forecasting (for EE (import and export) and NG), it is essential to note that the quality of the optimization results depends on the accuracy of the forecasts.

The usage of actual forecasting values and algorithms will be examined comprehensively in future work as part of a planned effort to address the well-known challenges of managing uncertainty in decision-making related to future events, as well as different optimizers, such as metaheuristics. This could strengthen the robustness of the system’s decisions.

Another interesting future work is the extension of the solution with metaheuristics optimizations, allowing nonlinear behavior, nonconvex efficiency curves, and uncertainty-aware constraints. Even though with the dimensions of this paper’s energy system model, metaheuristics may not give any advantages. In this direction, we can consider the implementation of more detailed operational constraints to allow a comparison of the results to know if the improvements are significant enough with the additional complexity. When trying to move into functional optimization services for a real scenario, experimentation between different levels of complexity could provide good insights. Further extension of this work is to analyze a full year of data to understand the sensitivity to seasonality and recognize common patterns.

Finally, it could be interesting to adjust the cost calculation of energy to the company’s actual behavior by combining participation in the different energy markets (forward markets, futures markets, etc.). These changes will probably suggest changing the optimization horizons to longer periods and can be combined with other future work.

Acknowledgement: Not applicable.

Funding Statement: This research has received funding from the European Union’s Horizon Europe research and innovation program under the project FLEXIndustries, Grant Agreement No. 101058453. The authors would like to acknowledge the European Commission for its support.

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, Marcelo Dario Rodas Britez; methodology, Marcelo Dario Rodas Britez and Francesco Ghionda; software, Marcelo Dario Rodas Britez; validation, Marcelo Dario Rodas Britez and Francesco Ghionda; formal analysis, Marcelo Dario Rodas Britez and Francesco Ghionda; investigation, Marcelo Dario Rodas Britez, Francesco Ghionda and Vasileios Tatsis; resources, Marcelo Dario Rodas Britez; data curation, Marcelo Dario Rodas Britez and Francesco Ghionda; writing—original draft preparation, Marcelo Dario Rodas Britez and Francesco Ghionda; writing—review and editing, Marcelo Dario Rodas Britez, Francesco Ghionda, Vasileios Tatsis and Dimosthenis Ioannidis; visualization, Marcelo Dario Rodas Britez and Francesco Ghionda; project administration, Dimosthenis Ioannidis. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The original contributions presented in the study are included in the article. Further inquiries can be directed to the Corresponding Author, MRB, at reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

Abbreviations

References

1. De Souza RJ, Reini M, Serra LM, Lozano MA, Nadalon E, Casisi M. Multi-objective optimization of an energy community powered by a distributed polygeneration system. Energies. 2024;17(13):3085. doi:10.3390/en17133085. [Google Scholar] [CrossRef]

2. Pinto ES, Amante B. Polygeneration system optimization for building energy system retrofit: a case of study for TR5 building of UPC-Terrassa. Energy Build. 2022;273(5):112375. doi:10.1016/j.enbuild.2022.112375. [Google Scholar] [CrossRef]

3. Rong A, Su Y. Polygeneration systems in buildings: a survey on optimization approaches. Energy Build. 2017;151(1):439–54. doi:10.1016/j.enbuild.2017.06.077. [Google Scholar] [CrossRef]

4. Liu P, Pistikopoulos EN, Li Z. A multi-objective optimization approach to polygeneration energy systems design. AlChE J. 2010;56(5):1218–34. doi:10.1002/aic.12058. [Google Scholar] [CrossRef]

5. Liu P, Gerogiorgis DI, Pistikopoulos EN. Modeling and optimization of polygeneration energy systems. Catal Today. 2007;127(1–4):347–59. doi:10.1016/j.cattod.2007.05.024. [Google Scholar] [CrossRef]

6. Tian X, Sun J, Xu T, Cui M, Wang X, Guo J, et al. Dynamic simulation and performance analysis on multi-energy coupled CCHP system. Energy Eng. 2022;119(2):723–37. doi:10.32604/ee.2022.015982. [Google Scholar] [CrossRef]

7. Akbari S, Saeed Fazel S, Hashemi-Dezaki H. Energy management of networked smart railway stations considering regenerative braking, energy storage system, and photovoltaic units. Energy Eng. 2023;120(1):69–86. doi:10.32604/ee.2022.024121. [Google Scholar] [CrossRef]

8. Xu Y, Liu J, Sun T, Kang H, Hao X. Research on optimal scheduling of integrated energy systems with wind-photovoltaic-biogas-storage considering carbon capture systems and power-to-gas coordination. Energy Eng. 2025;122(8):3155–76. doi:10.32604/ee.2025.065753. [Google Scholar] [CrossRef]

9. Ren H, Zhou W, Nakagami K, Gao W, Wu Q. Multi-objective optimization for the operation of distributed energy systems considering economic and environmental aspects. Appl Energy. 2010;87(12):3642–51. doi:10.1016/j.apenergy.2010.06.013. [Google Scholar] [CrossRef]

10. Urbanucci L. Limits and potentials of mixed integer linear programming methods for optimization of polygeneration energy systems. Energy Proc. 2018;148(52):1199–205. doi:10.1016/j.egypro.2018.08.021. [Google Scholar] [CrossRef]

11. Lin MH, Carlsson JG, Ge D, Shi J, Tsai JF. A review of piecewise linearization methods. Math Probl Eng. 2013;2013(2):101376–8. doi:10.1155/2013/101376. [Google Scholar] [CrossRef]

12. Kämper A, Geers P, Leenders L, Bardow A. Adaptive rolling horizon for operational optimization of multi-energy systems. In: Proceedings of the 34th International Conference on Efficiency, Cost, Optimization, Simulation and Environmental Impact of Energy Systems (ECOS 2021); 2021 Jun 27–Jul 2; Taormina, Italy. p. 116–27. doi:10.52202/062738-0011. [Google Scholar] [CrossRef]

13. Teichgraeber H, Brandt AR. Time-series aggregation for the optimization of energy systems: goals, challenges, approaches, and opportunities. Renew Sustain Energy Rev. 2022;157(11):111984. doi:10.1016/j.rser.2021.111984. [Google Scholar] [CrossRef]

14. Abdmouleh Z, Gastli A, Ben-Brahim L, Haouari M, Al-Emadi NA. Review of optimization techniques applied for the integration of distributed generation from renewable energy sources. Renew Energy. 2017;113(3):266–80. doi:10.1016/j.renene.2017.05.087. [Google Scholar] [CrossRef]

15. Rabe M, Bilan Y, Widera K, Vasa L. Application of the linear programming method in the construction of a mathematical model of optimization distributed energy. Energies. 2022;15(5):1872. doi:10.3390/en15051872. [Google Scholar] [CrossRef]

16. Ubando AT, Chen WH, Tan RR, Naqvi SR. Optimal integration of a biomass-based polygeneration system in an iron production plant for negative carbon emissions. Int J Energy Res. 2020;44(12):9350–66. doi:10.1002/er.4902. [Google Scholar] [CrossRef]

17. Hart WE, Laird CD, Watson JP, Woodruff DL, Hackebeil GA, Nicholson BL, et al. Pyomo—optimization modeling in python. Cham, Switzerland: Springer; 2017. doi:10.1007/978-3-319-58821-6. [Google Scholar] [CrossRef]

18. Li H, Wang Z, Hong T, Piette MA. Energy flexibility of residential buildings: a systematic review of characterization and quantification methods and applications. Adv Appl Energy. 2021;3:100054. doi:10.1016/j.adapen.2021.100054. [Google Scholar] [CrossRef]