Economic–Energy–Environmental Optimization of a Multi-Energy System in a University District

Abstract

The integration of energy generation and consumption is one of the most effective ways to reduce energy-system-related waste, costs, and emissions in cities. This paper considers a university district consisting of 32 buildings where electrical demand is currently met by the national grid, and 31% of thermal demand is supplied by a centralized heating station through a district heating network; the remainder is covered by small, dedicated boilers. Starting from the present system, the goal is to identify “retrofit” design solutions to reduce cost, environmental impact, and the primary energy consumption of the district. To this end, three new configurations of the multi-energy system (MES) of the district are proposed considering (i) the installation of new energy conversion and storage units, (ii) the enlargement of the existing district heating network, and (iii) the inclusion of new branches of the electrical and heating network. The configurations differ in increasing levels of integration through the energy networks. The results show that the installation of cogeneration engines leads to significant benefits in both economic (up to −12.3% of total annual costs) and energy (up to −10.2% of the primary energy consumption) terms; these benefits increase as the level of integration increases. On the other hand, the limited availability of space for photovoltaics results in increased CO₂ emissions when only total cost minimization is considered. However, by accepting a cost increase of 8.4% over the least expensive solution, a significant reduction in CO₂ (−23.9%) can be achieved while still keeping total costs lower than the existing MES.

1. Introduction

The urgent need to reduce carbon emissions is driving the energy sector toward harnessing distributed resources, integrating generation and consumption at the local level, and exploiting synergies between end-use sectors [1]. In fact, more interconnected energy infrastructures allow increasing efficiency, reducing primary energy consumption and increasing the penetration of renewable sources [2]. Accordingly, multi-energy systems are emerging as promising configurations due to their holistic approach in decarbonizing power, thermal, and other sectors while minimizing economic impacts [3]. These benefits can be promoted by the increased flexibility deriving from the active role of end users in energy management, including active participation in demand-side management programs [4].

In general terms, a multi-energy system (MES) is an energy system of any geographical extension, from a single building to a national system, that combines multiple energy vectors (as, for instance, electricity, heat, and fuels) and fulfils the end users’ demand for different forms of energy [5]. The design and operation of an MES is a challenging task that requires proper optimization methods, due to the necessity of contextually considering energy conversion, storage, delivery, and utilization [6]. The goal is to properly evaluate the number, type, size, and management over time of the energy conversion and storage units meeting the demand of the end users, as well as the type, layout, and capacity of the energy networks needed for interconnection. This can be formulated, in general, as a mixed-integer non-linear programming (MINLP) problem due to the coexistence of both binary and continuous decision variables (as the existence of a component and its size, respectively) and the presence of non-linearities (e.g., off-design maps) [7]. The large number of decision variables involved can easily make the problem unsolvable with currently available computational technologies [8]. However, linearization approaches that transform the problem into a mixed-integer linear programming (MILP) one allow simplifying the optimization process, thereby drastically reducing computational requirements [9]. Rech [10] proposed a rigorous methodology for modeling the components of an energy system within an MILP problem. This methodology is used in this work as well.

The aforementioned design and operation optimization of an MES is an inherently integrated problem that requires being solved in “one shot”. However, the common approach in the literature is to consider the optimization of energy conversion and storage systems separately from the optimization of energy networks. This may lead to solutions that are not optimal for the system as a whole. For instance, in a cost minimization problem, the optimal solution for energy conversion may not be the global optimum because of the failure to consider the interconnection costs with end users.

Few examples are provided of works dealing with the optimization of energy conversion and storage units while neglecting the design of the networks. Dal Cin et al. [11] solved the multi-objective design and operation optimization of an MES serving the electricity and heating demand of a renewable energy community, with the goal of minimizing life-cycle costs and CO₂ emissions. They considered photovoltaic (PV) panels, cogeneration units, and storage systems (both electrical and thermal) but did not model the energy networks. Rech et al. [12] minimized the total cost of an MES providing electricity and heat to a small mountain town by means of renewable plants. Heat is distributed to the users through a district heating network (DHN) that is assumed ideal and modeled as a “black box” connecting generators, storage systems, and consumers. Wirtz [13] proposed a web tool for the optimal conceptual design of the generation mix in a district MES, to define the optimal size of the energy conversion and storage plants available in an energy hub that fulfils electricity, heating, and cooling demands. Heating and cooling networks can also be considered and optimized in terms of operation but not in terms of design (i.e., the capacity of a line is assumed to be sufficient to transport all the energy required). Mashayekh et al. [14] developed an optimization method for the design of the energy conversion and storage plants supplying a multi-energy microgrid that involves both electricity and heating networks. The optimization addresses the generation mix selection and sizing, the resource siting and allocation, and the operation scheduling. However, the design of the networks is provided as input and cannot be improved.

On the other hand, works that optimize the topology and size (i.e., rated power that can flow in each branch) of energy networks usually consider energy conversion systems as input data, thereby neglecting their design. For instance, Pizzolato et al. [15] optimized the topology and capacity of a DHN in Turin, Italy, and analyzed the inclusion of network loops to increase the reliability of the system. They considered an existing heat generation facility to provide the required thermal energy. Röder et al. [16] optimized the design of a DHN serving a mixed residential–commercial district in Germany by searching for the cost-optimal solution with the inclusion of distributed thermal storage fixing in advance of the position, type, and size of the energy conversion plants and storage systems.

Works in the literature that considered the design and operation optimization of both the energy conversion and storage systems and energy networks of an MES as a single problem applied many simplifications to ensure the achievement of a solution in a reasonable time. For instance, Keirstead et al. [17] optimized the design and operation of the electricity, heating, and gas networks, as well as conversion systems (mainly cogeneration units), of an urban MES. The one-day operation of the MES is aggregated into two representative time intervals that allow solving the problem in an acceptable time (the yearly operation, composed of few representative days, is still described by a small number of time intervals). The drawback is that the coarse representation of the simulated time frame is not sufficient to properly model the short-term variability of time-dependent quantities (e.g., availability of renewable sources, changes in energy demand). Sidnell et al. [18] optimized the design and operation of a neighborhood MES integrating electricity and heating networks with renewable plants (PV) and cogeneration units. They considered a simplified network representation in which network lines are modeled as linear segments that connect adjacent buildings, regardless of the real geographical layout of the MES. Lerbinger et al. [19] provided a more detailed model of the DHN in an MES by constraining heating network lines to follow only predetermined paths (such as the streets of a neighborhood), thus avoiding unfeasible connections between nodes. Morvaj et al. [20] highlighted that a “street-following approach” for the optimal design of energy networks guarantees more realistic solutions than a “green-field approach”. Dal Cin et al. [21] adopted the street-following approach to optimize the design and operation of a district MES including both electricity and heating networks, as well as PV plants, cogeneration engines, heat pumps, and both thermal and electrical storage systems. Moreover, they made use of typical days obtained by K-medoids clustering to reduce computational complexity while ensuring a sufficient accuracy of the simulated time frame. It should be emphasized that in the mentioned paper, the authors also introduced the concept of “retrofit design”, which aims at adding new components or additional installed capacities to an existing system in order to improve its initial layout. Finally, the same authors presented in [22] a general method for the integrated optimization of an MES considered as a whole, i.e., by considering energy conversion units, storage units, and the energy network in the same synthesis, design, and operation (SDO) optimization problem. However, the latter works [21,22] focus on the optimization methodology, the validity of which is demonstrated through hypothetical case studies. They lack, instead, implementations in real-case studies that further strengthen the potential of the proposed method in real applications.

Among the several studies on MESs available in the literature, this paper focuses on those related to university districts, where proposed MES systems have proven to be by far superior solutions in terms of reducing cost and environmental impact compared to the conventional systems based on boiler and electricity withdrawal from the grid. Moreover, the analysis of this type of MES is interesting because the associated studies are often based on measured demand data rather than simulated or standard consumption profiles, thus providing realistic solutions to practical problems. Gabrielli et al. [23] optimized the MES operation of a university facility in Zurich, Switzerland, by focusing on geothermal storage. It turned out that CO₂ emissions can be reduced by 87% compared to a conventional system based on centralized heating and cooling. Martelli et al. [24] studied the subsidies and taxes required to achieve specified decarbonization targets in MESs. Considering a university campus in Parma, Italy, and optimizing both the design and operation of the associated MES, they found that a carbon tax of about 160 EUR/ton favors the installation of PV systems, thereby reducing CO₂ emissions by 25%. Venturini et al. [25] also included a life-cycle assessment in the design and operation optimization problem. The analysis of a university facility in Parma, Italy, showed that the optimal sizing and operation of cogeneration systems results in primary energy savings of about 15% and total cost reductions of about 12% over the life cycle of the MES. Testi et al. [26] developed a multi-objective stochastic optimization methodology to evaluate the integrated optimal sizing and operation of an MES under uncertainties in climate, space occupancy, energy loads, and fuel costs. They considered a university campus in Trieste, Italy, as a test case by focusing on the integration of heat pumps, which can increase the share of renewable energy, reduce the operating cost of the system, and moderate the investment risk. Dos Santos et al. [27] focused on the optimal design of the electrical microgrid serving a university campus in São Paulo, Brazil. The authors of the study optimized the size and location of distributed energy conversion and storage units (mainly photovoltaic plants and battery storage systems) and also the capacity of the cables constituting the microgrid. Finally, Comodi et al. [28] modeled the MES of a university facility in Singapore, with the goal of defining the optimal mix of energy conversion and storage systems and energy network infrastructure needed to meet electricity and cooling demands. They found that a district cooling network can reduce the total installed capacity of electric chillers, which results in capital cost savings. It is worth noting that the aforementioned design optimization studies modeled the MES of the university districts “from scratch” regardless of the pre-existing system configurations, thus achieving optimal layouts that are substantially alternative to the available ones. How to apply “retrofit design” optimization to this type of MES in order to improve existing system configurations (i.e., by starting with the components already available and replacing or resizing them) requires further research.

This paper considers a university district in Padova, northern Italy, and explores new and “smart” solutions to reduce the cost, primary energy consumption, and environmental impact of meeting its electricity and heating demand. This study takes the existing system configuration as a starting point and applies the retrofit design approach proposed in [21] for the optimization of the associated MES. In the existing system layout, the heating and summer-cooling demands of buildings are mostly covered by autonomous small boilers fired by natural gas and compression refrigeration systems, respectively. The buildings belonging to a part of the district are already interconnected through a district heating network (DHN) and a local electrical distribution network (EDN). A centralized thermal power generation system covers via the DHN the heat demand of this part, corresponding to 31% of the total heating demand of the district. The goal of this study is to identify the capacity and operation of the new plants to be installed, including photovoltaic (PV) panels, air–water heat pumps (HPs), and gas-fired cogeneration internal combustion engines (CHP ICEs), as well as the additional capacity and possible expansion of the available DHN and EDN. More specifically, three new energy system layouts are proposed and compared to determine the best solution for the district in terms of cost, primary energy consumption, and environmental impact. The first layout maintains the same structure and capacity of the existing DHN and EDN. The second and third layouts consider instead the expansion of the DHN and EDN. The optimization approach proposed in [21] is then applied to the three layouts to find the most cost-effective design and operation of the district MES (Cost Minimization scenario, CM). It is worth emphasizing that the optimization problem is based on electricity and natural gas consumption data measured over multiple years, from which annual electricity and heating demands are derived with an hourly resolution. This allows the optimization method proposed in [21] to be tested on a real-case study (conversely, the original study made use of standardized demand curves). Two additional scenarios are finally identified by imposing target reduction values on carbon emissions (Low Carbon Emissions scenario, LCE) and primary energy consumption (Primary Energy Saving scenario, PES) as secondary objectives (“epsilon-constraints”) in the optimization procedure. The goal is to identify the suitable trade-off between cost, energy efficiency, and environmental impact.

What is new compared to the literature on MESs in university districts is the implementation of the retrofit design problem to improve the economic and environmental performance of the existing system and its components. This makes it possible to take into account the energy infrastructure already available and to decide on new interventions in the current system. This approach contrasts with that commonly used in the literature, which is based on the design “from scratch” of new system configurations regardless of the components already available, thus providing MES configurations that are completely alternative to the existing ones. In addition, through the consideration of a real-case study based on measured data, this paper provides a practical validation to the optimization methodology proposed in [22], which has been used for “retrofit design” in other contexts but has been applied so far only to hypothetical test cases.

The rest of this paper is organized as follows. Section 2 describes in detail the considered university district with the existing energy infrastructure and outlines the three new MES configurations proposed to improve the existing configuration. Section 3 concerns the methodology and focuses on the input data to be provided to the optimization problem, as well as on the mathematical formulation of the optimization problem. Section 4 shows the results in terms of installed capacities, investment and operational costs, primary energy consumption, and CO₂ emissions for the three new MES configurations in the considered scenarios (CM, LCE, and PES). Section 5 discusses the results in order to determine the best solution for the university district among those proposed. Finally, Section 6 draws conclusions.

2. System Description

This section describes the university district and the proposed configurations considered in the optimization problem in order to explore new solutions to reduce the costs, primary energy consumption, and environmental impact of the district.

2.1. Reference Case

The multi-energy system (MES) considered as a case study is a university district located in the northeast of Padova, Italy. The size of the district area is 422,000 m², consisting of 32 university buildings placed in a residential neighborhood. The entire district, shown in Figure 1a, is crossed by a river that separates the north and south areas.

Currently, the yearly electrical demand of 22.5 GWh is covered by withdrawing electricity from 14 PODs connected to the national grid, all referring to the same primary substation. The yearly heat demand of 12.6 GWh is covered only by natural-gas-fired boilers. The north area accounts for 50% of the electrical demand and 40% of the heating demand, fulfilled mostly by a single POD and a centralized boiler, exploiting the existing local EDN and DHN (see blue/red line connecting node 6 to 9/8 in Figure 1b). The energy system of the south area presents a more fragmented situation; the electrical energy is withdrawn by 11 PODs, and each building is equipped with small boilers.

Hourly data of the electrical consumption are available for each POD, starting from 2019 to 2022. For the heat demand, the only data available are the cumulative monthly gas consumptions measured at each of the 21 gas delivery points (GDPs) during 2019 and 2020.

The MES outlined in Figure 1b is modeled as a multi-nodal system composed of N = 20 nodes subdivided into two types, i.e., “aggregation” and “connection” nodes. The former represent energy withdrawal points associated with aggregated buildings, whereas the latter shape both existing and new energy networks by identifying key points where network branches divide or join.

Table 1. Characterization of the nodes of the starting MES layout.

Node	Type	POD	Installed GB Capacity [kWth]
0	Aggregation node	✓	2427
1	Aggregation node	✓	2274
2	Aggregation node	✓	540
3	Aggregation node	✓	750
4	Aggregation node	✓	1291
5	Aggregation node	✓	1368
6	Connection node	-	-
7	Aggregation node	✓	9000
8	Connection node	-	-
9	Aggregation node	-	220
10	Aggregation node	✓	1980
11	Aggregation node	-	380
12	Aggregation node	✓	313
13	Connection node	-	-
14	Connection node	-	-
15	Connection node	-	-
16	Connection node	-	-
17	Connection node	-	-
18	Connection node	-	-
19	Connection node	-	-

provides a comprehensive description of each node in the MES, including the existing installed capacity of gas boilers (GBs). The current total capacity of GBs, most of which are condensing boilers, is 20.5 MW (including backup units). “Aggregation” nodes are nodes 0 to 5, node 7, and nodes 9 to 12. The criteria adopted for the aggregation of buildings are as follows:

• Sharing: buildings connected to the same POD, conversion unit, or pre-existing distribution network;
• Proximity: nearby buildings, not separated by public streets;
• Space availability: availability of existing indoor rooms or usable outdoor spaces to be used for technical rooms.

“Connection“ nodes are all the remaining nodes (6, 8, and 13 to 20). Their location is chosen to follow the paths of the existing network branches (nodes 6 and 8 identify the existing EDN and DHN together with aggregation nodes 7 and 9, see Figure 1b) and to draw possible new network branches, considering all constraints associated with their installation in accordance with a street-following approach.

Despite its central location (1 km from the city center), the district presents a high spatial density of energy consumption compared with residential neighborhoods; thus, it represents a good test case for the study of energy intervention strategies, combining the advantage of working on a relatively small area with the opportunity to affect a substantial amount of consumptions and related emissions. However, the historical-interest restriction on many of the buildings, the limited roof availability for PV systems, and the absence of other renewable sources strongly limit the renewable share and possible range of interventions.

2.2. Proposed Configurations

Three different MES configurations are considered in the optimization problem. Given the limitations of the case study in terms of renewable energy share, the proposed configurations consider the extension of energy networks that may favor centralized systems over autonomous plants and lead to reduced costs while achieving a higher average conversion efficiency. An increasing number of potential new network branches characterizes the three MES configurations. This may lead to an increasing integration of the MES, intended as the degree of interconnection between nodes through the energy networks and the consequent flexibility in energy distribution:

• Existing networks: only the pre-existing DHN and EDN are considered (Figure 1b);
• Medium integration: the north, south-east, and south-west areas constitute three independent DHNs and EDNs, each with a sufficient availability of technical space to accommodate all types of energy conversion and storage units (Figure 2a);
• High integration: a single DHN and EDN extend throughout the district (Figure 2b).

The design of each configuration starts with the current MES assets (i.e., existing boilers, PODs, DHN, and EDN) and evaluates the installation of new energy conversion and storage units and additional network infrastructure. In particular, the proposed solutions combine the installation of PV systems, where possible, and HPs with an effective use of the available fossil resources (natural gas). To this end, integration is promoted with the installation of centralized combined heat and power internal combustion engines (CHP ICEs). The possible installation of electric and thermal energy storage (EES and TES) units is also considered.

The configurations are all characterized by the same dataset, reported in

Table 2. Candidate energy conversion and storage units to be included in the optimal MES configuration and their maximum capacities.

Node	POD	Max GB Capacity [kWth]	Max PV Capacity [kWel]	Max CHP ICE Capacity [kWel]	Max HP Capacity [kWth]	Max EES Capacity [kWh]	Max TES Capacity [kWh]
0	✓	3000	39	-	2500	250	-
1	✓	3000	-	4000	2500	10,000	20,000
2	✓	1000	-	-	600	-	-
3	✓	1000	58	-	800	250	-
4	✓	2000	-	-	1500	-	-
5	✓	2000	-	-	1500	-	-
7	✓	9000	272	4000	4000	10,000	20,000
9	-	500	-	-	300	-	-
10	✓	2500	-	4000	2000	10,000	20,000
11	-	500	116	-	600	600	-
12	✓	500	-	-	600	-	-

, which defines the type and maximum size of the energy conversion and storage units that can be installed at the aggregation nodes. The possible locations and maximum sizes of new installations have been identified by considering the available roof area and technical rooms for PV systems (based on a feasibility study initiated by the University of Padova) and CHP ICEs, respectively; the size of existing GBs for HPs (which can be installed at any node associated with a heating demand); the installation of PV systems for small EESs; and the installation of ICEs for large EESs and TESs. In this regard, the maximum sizes of the EESs and TESs at nodes 1, 7, and 10 are intended as sufficiently large thresholds for exploring the optimization results.

Finally, it is worth noting that the proposed optimization of the design considers not only the installation of new units and network branches but also the retrofit of the existing ones, i.e., the addition of new capacity to units and networks initially available in the reference case.

3. Methods

Figure 3 provides an overview of the procedure used to optimize the design and operation of the MES configurations in Section 2.2.

The “design choices” are already discussed for each configuration in Section 2.2; the “case study” and “techno-economic data” steps are described in Section 3.1. Finally, Section 3.2 describes the mathematical formulation of the “MILP design and operation optimization problem” and the “criteria” adopted for each scenario.

3.1. Input Data

The optimization problem requires as input data the following time series with an hourly resolution:

• electricity demand of each node associated with an electrical load;
• heating demand of each node associated with a thermal load;
• global solar irradiance, required to calculate the PV generation (source PVGIS [29]);
• ambient temperature, required to assess the coefficient of performance (COP) of HPs (source PVGIS [29]).

The yearly time series of electricity demand has been obtained by averaging the hourly data available for the years 2019 and 2022, neglecting 2020 and 2021 due to the bias introduced by the COVID-19 pandemic. For the heating demand, only the cumulated monthly gas consumptions of each GDP were available, requiring some assumptions in order to build up the hourly time series needed to solve the optimization problem. A simple approach has been adopted, assuming a flat demand curve during the scheduled operating hours of the gas boilers on weekdays, Saturdays, and holidays, for each month of the year. The resulting curves are then normalized with respect to the measured monthly consumption, ensuring that the integral of the monthly demand equals the measured data. Finally, data from PODs and GDPs are stacked into node demands, according to the criteria adopted for the aggregation presented in Section 2 (energy demands of the aggregated buildings are summed hour by hour).

To reduce the computational effort required to solve the optimization problem, the yearly demand time series data are aggregated into K = 36 typical days, considering three days for each month of the year, representative of the average weekday, Saturday, and Sunday/holiday, respectively. Each typical day is defined as H = 24 hourly time steps and is associated with a weight that is equal to the number of represented days in the year (the sum of all weights is 365). Figure 4 shows the energy demand curves of two typical days at one of the aggregation nodes. The shape of the demand profile is homogeneous throughout the district and characterized by the simultaneity of the electrical and thermal demands during the heating period. The choice of 36 typical days allows accounting for both weekly and annual energy demand seasonality. The weekly seasonality is related to the different occupancy of buildings during weekdays, Saturdays, and Sunday/holidays, while seasonal seasonality is mainly related to different environmental conditions (e.g., ambient temperature and daylight hours). For instance, during the heating period (October 15–April 15), the existing heating system (gas boilers) operates from 7 to 19 on weekdays and from 7 to 17 on Saturdays from November to March; from 7 to 17 on weekdays and from 8 to 14 on Saturdays in October and April; and is always off during holidays. Seasonality can be observed in the electric and thermal energy balances (Figure 5), for each day of the year, in the reference case. The heat balance shows clearly the periods when the heating system is in operation (winter, weekdays, and Saturday, daytime); the electric balance shows that the electrical demand during summer is the highest due to the use of compression refrigeration systems. It is worth recalling that in the reference case, the entire heating demand is covered by gas boilers, while the entire electricity demand is met by the national electric grid.

In addition, the data related to the system topology must be provided. These data include the Cartesian position of each node with respect to two reference axes, the position and capacity of the available network branches (identified by the extreme nodes), and the position and capacity of each available conversion unit along with those new units that can be installed. To this end, the considered technologies are as follows:

• Photovoltaic (PV);
• Combined heat and power gas-fired internal combustion engine (CHP ICE);
• Gas-fired boiler (GB);
• Air–water heat pump (HP);
• Thermal energy storage (TES) based on hot water tanks;
• Electrical energy storage (EES) based on lithium batteries.

The other information required to solve the model is the techno-economic data of the energy conversion and storage technologies (e.g., lifetime of the units, investment costs), which are summarized in

Table 3. Linearized investment cost and lifetime of the energy conversion and storage units and distribution networks.

Technology	Investment Cost (cinv)		Lifetime [Year]
	cinv,var	cinv,fix
PV	1250 EUR/kW	-	20
CHP ICE	1398 EUR/kWel	25.8 kEUR	20
GB	65 EUR/kW	1.6 kEUR	20
HP	343 EUR/kWth	6.3 kEUR	20
EES	880 EUR/kWh	3.5 kEUR	20
TES	244 EUR/kWh	1.0 kEUR	20
DHN	0.2 EUR/kWth/m	103 EUR/m	40
EDN	0.006 EUR/kWel/m	34 EUR/m	40

.

Finally,

Table 4. Specific cost and emissions of the considered energy carriers.

Carrier	Cost [EUR/MWh]	Emission Factor [kgCO2/MWh]
Electricity from the grid	259	271
Natural gas	101	197

summarizes the specific costs of each energy carrier (i.e., electricity and natural gas), which are evaluated by averaging the university bills collected from October 2022 to September 2023, and their CO₂ emission factors are also provided, which consider the Italian average generation mix of the electricity withdrawn from the national grid.

3.2. Optimization Problem

The design and operation optimization of the MES configurations in Section 2.2 is formulated as a mixed-integer linear programming (MILP) problem:

$$\begin{array}{l}\mathrm{F}\mathrm{i}\mathrm{n}\mathrm{d}{\mathbf{x}}_D^{\mathbf{*}}\mathrm{and}{\mathbf{x}}_O^{\mathbf{*}}\left(t\right)\in {\mathbb{R}}^{\mathrm{n}}\mathrm{or}{\mathbb{I}}^{\mathrm{m}}\mathrm{t}\mathrm{h}\mathrm{a}\mathrm{t}\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}f{\left({\mathbf{x}}_D,{\mathbf{x}}_O\left(t\right)\right)}_Y={\mathit{a}}^T{\mathbf{x}}_D+{\mathit{b}}^T{\mathbf{x}}_O\left(t\right)\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\mathit{g}\left({\mathbf{x}}_D,{\mathbf{x}}_O\left(t\right)\right)=\mathit{C}{\mathbf{x}}_D+\mathit{D}{\mathbf{x}}_O\left(t\right)\ge \mathit{e}\end{array}$$

(1)

where $f$ is the linear objective function; ${\mathbf{x}}_D$ and ${\mathbf{x}}_{\mathrm{O}}\left(t\right)$ are the decision variables associated with the design (constant in the whole period $Y$) and operation (time-varying) of the MES, respectively; and $\mathit{g}$ represents the linear equality/inequality constraints that make up the model of the MES. The variables and equations involved in the problem in Equation (1) are discussed in detail in the following.

The relationships $\mathit{g}$ include (i) the electric and thermal balances of the MES, (ii) the characteristic equations describing the performance of the energy conversion and storage units, and (iii) the constraints related to the capacity of the energy networks lines.

The electric and thermal balances (relationships (i)) are defined for each node of the MES (see Figure 1a) and written in accordance with [22]. Each balance, as shown in Equation (2), imposes that the sum of all electric/thermal power flows that enter, are produced, or discharged in the node must be equal to the sum of the electric/thermal power flows that exit, are consumed, or charged in the node.

$$\left(F_{imp}\left(t\right)\right)+\sum _i{P}_i\left(t\right)+\sum _j{P}_j\left(t\right)+\sum _z{T}_z\left(t\right)·\left(1-{\lambda }_z\left(t\right)·l_z\right)=\left(P_{exp}\left(t\right)\right)+F_{dem}\left(t\right)+\sum _{i^{\prime }}{F}_{i^{\prime }}\left(t\right)+\sum _{j^{\prime }}{F}_{j^{\prime }}\left(t\right)+\sum _{z^{\prime }}{T}_{z^{\prime }}\left(t\right)$$

(2)

where $F_{imp}\left(t\right)$ and $P_{exp}\left(t\right)$ are included in the electric balances only and represent the electric power flows that are imported and exported through the node POD, respectively; $P_i\left(t\right)$ is the power output of each energy conversion unit $i$ in the node; $P_j\left(t\right)$ is the power flow discharged from each energy storage unit $j$ in the node; $T_z\left(t\right)$ and $T_{z^{\prime }}\left(t\right)$ are the power flows entering and exiting the node through the local EDN/DHN, respectively; $F_{dem}\left(t\right)$ is the power flow associated with the electrical/thermal demand; $F_{i^{\prime }}\left(t\right)$ is the power input of each energy conversion unit $i^{\prime }$ in the node; and $F_{j^{\prime }}\left(t\right)$ is the power flow charged in each energy storage unit $j^{\prime }$ in the node. In Equation (2), the subscripts $i$ and $j$ refer exclusively to units that produce or store the considered energy carrier (electricity or heat) in the node, while the subscripts $i^{\prime }$ and $j^{\prime }$ refer to units that consume or store the considered energy carrier in the node. Note that the power flow entering a node through the local EDN/DHN equals the power flow leaving the node at the other end of the branch, minus the losses (${\lambda }_z\left(t\right)$ is the network-specific loss per unit of network length $l_z$).

Relationships (ii) are written in accordance with [10,22], subdividing the energy conversion units into dispatchable and non-dispatchable ones.

The input power (fuel of the unit, $F_i\left(t\right)$) of a dispatchable conversion unit $i$ (i.e., CHP ICE, GB, HP) varies linearly with the power output (product of the unit, $P_i\left(t\right)$), as given in Equation (3), where $k_i\left(t\right)$ is a correction factor depending on the input data (e.g., ambient temperature in the case of HP), $C_i$ and $D_i$ are constant coefficients linearizing the off-design performance map of the unit, and ${\delta }_i\left(t\right)$ is a binary variable describing the on/off status of the unit.

$$F_i\left(t\right)=k_i\left(t\right)·\left(C_i·P_i\left(t\right)+D_i·{\delta }_i\left(t\right)\right)$$

(3)

The power output is upper and lower bounded according to the unit rated capacity $P_i^{max}$, as shown in Equation (4), which includes the auxiliary variable $M_i\left(t\right)$ to avoid non-linearities [30].

$$\begin{array}{c}\begin{array}{c}{P}_i\left(t\right)\le P_i^{max}·{\beta }_i\\ M_i\left(t\right)\le P_i^{max}·{\delta }_i\left(t\right)\end{array}\\ P_i\left(t\right)-M_i\left(t\right)\le \left(1-{\delta }_i\left(t\right)\right)·P_i^{max}\\ m_i·P_i^{max}\le P_i\left(t\right)\le P_i^{max}\end{array}$$

(4)

In Equation (4), ${\beta }_i$ is a binary variable (constant in the total period) describing the existence of the unit $i$, and $m_i$ is the minimum load of the unit $i$ referred to its capacity $P_i^{max}$. Note that if ${\beta }_i=0$, the power output $P_i\left(t\right)$ is zero for all the time steps in the total period (i.e., the unit $i$ is excluded from the MES), whereas when ${\beta }_i=1$, the power output $P_i\left(t\right)$ can vary between the minimum ($k_i·P_i^{max}$) and maximum ($P_i^{max}$) load when ${\delta }_i\left(t\right)=1$ (the unit $i$ is on) or is equal to zero when ${\delta }_i\left(t\right)=0$ (the unit $i$ is off).

For CHP ICEs, only the additional linear relationship in Equation (5) is added to calculate the available thermal power that can be recovered (second product of the unit, $Q_i$) as a function of the unit load ($P_i$)

$$Q_i\left(t\right)\le C_{2,i}·P_i\left(t\right)+D_{2,i}·{\delta }_i\left(t\right)$$

(5)

The power output of PV is calculated using Equation (6), where the correction factor $k_{PV}\left(t\right)$ is a function of the solar irradiance in the time step $t$.

$$P_{PV}\left(t\right)=k_{PV}\left(t\right)·P_{PV}^{max}$$

(6)

The state of charge of a storage unit $j$ (TES or EES) is modeled as the sum of the contribution from an intra-day state of charge (periodic within each typical day) and an inter-day state of charge (periodic within the total period of one year) [31].

The intra-day state of charge ($E_j\left(t\right)$) accounts for the daily storage and is calculated using the dynamic energy balance in Equation (7), where $k_j\left(t\right)$ is a correction factor depending on the input data and the relative coefficient of self-discharge losses of the unit, and ${\eta }_{char,j}$ and ${\eta }_{disc,j}$ are the charging and discharging efficiencies, respectively.

$$E_j\left(t\right)=k_j\left(t\right)·E_j\left(t-1\right)+{\eta }_{char,j}·F_j\left(t\right)-{\displaystyle \frac{1}{{\eta }_{disc,j}}}·P_j\left(t\right)$$

(7)

$E_j\left(t\right)$ is upper and lower bounded according to the maximum capacity of the storage unit $E_j^{max}$, as shown in Equation (8). Moreover, it is assumed that each storage unit is completely empty at the beginning of each typical day (second relationship in Equation (8)).

$$\begin{array}{c}{m}_j·E_j^{max}\le E_j\left(t\right)\le E_j^{max}\\ E_j\left(t=0\right)=0\end{array}$$

(8)

The inter-day state of charge accounts for seasonal storage and is calculated as proposed in [31].

The constraints at point (iii) above related to the capacity of the energy networks (EDN and DHN) are taken from [22]. The capacity $T_z^{max}$ of each energy network branch $z$ to carry electricity (heat) between two connected nodes (see Figure 1 and Figure 2) is upper bounded by a “big-M constraint”, where $T_z\left(t\right)$ is the power flow in the branch $z$, and $M_z$ is a “large enough” value, as given in Equation (9).

$$T_z^{max}\le M_z$$

(9)

Finally, Equation (10) constrains the power flow $T_z\left(t\right)$ in each network branch to be lower than the branch capacity.

$$T_z\left(t\right)\le T_z^{max}$$

(10)

For each dispatchable conversion unit $i$, the design decision variables (${\mathbf{x}}_D$) include the unit rated capacity $P_i^{max}$ and the binary variable ${\beta }_i$ in Equation (4), the latter describing the existence of the unit $i$. For PV and TES/EES, the only design decision variable is the maximum unit capacity ($P_{PV}^{max}$ and $E_j^{max}$), as well as for the energy network branches ($T_z^{max}$). A value of $P_{PV}^{max}$, $E_j^{max}$, or $T_z^{max}$ equal to zero corresponds to the exclusion of the specific unit/network branch from the optimal MES design.

For all dispatchable energy conversion units, storage units, and energy network branches, the operation decision variables (${\mathbf{x}}_O\left(t\right)$) include a power flow associated with unit/branch load ($P_i\left(t\right)$, $F_j\left(t\right)$ and $P_j\left(t\right)$, $T_z\left(t\right)$) for each of the $864$ hourly time steps (24 h $\times$ 36 typical days). In addition to these, all binary variables ${\delta }_i\left(t\right)$ describing the on/off status of the dispatchable energy conversion units (Equation (3)), and the thermal power output $Q_i\left(t\right)$ for CHP ICEs only (Equation (5)) are included among the operation decision variables.

The objective function to be minimized ($f$ in Equation (1)) is the annual cost of the MES, i.e., the annual levelized cost of investments $C_{inv}$ plus the annual operating costs $C_{oper}$:

$$f{\left({\mathbf{x}}_D,{\mathbf{x}}_O\left(t\right)\right)}_Y=C_{inv}+C_{oper}={\sum }_{i,j,z}{c}_{in{v}_{i,j,z}}+C_{oper}$$

(11)

$C_{inv}$ in Equation (11) is obtained as the sum of the annual levelized cost of investments of each unit ($c_{in{v}_{i,j,z}}$), which are calculated according to the following linear equations:

$$\begin{array}{ll}{c}_{in{v}_i}=\left(c_{inv,va{r}_i}·P_i^{max}+c_{inv,fi{x}_i}\right)·F_D& \mathrm{for}\mathrm{energy}\mathrm{conversion}\mathrm{units}\\ c_{in{v}_j}=\left(c_{inv,va{r}_j}·E_j^{max}+c_{inv,fi{x}_j}\right)·F_D& \mathrm{for}\mathrm{storage}\mathrm{units}\\ c_{in{v}_z}=\left(c_{inv,va{r}_z}·T_z^{max}+c_{inv,fi{x}_z}\right)·F_D& \mathrm{for}\mathrm{energy}\mathrm{network}\mathrm{branches}\end{array}$$

(12)

where $c_{inv,va{r}_{i,j,z}}$ and $c_{inv,fi{x}_{i,j,z}}$ are the variable and fixed linearization coefficients, respectively; $P_i^{max},E_j^{max}$, and $T_z^{max}$ are the unit capacities; and $F_D$ is the discount factor based on a 5% interest rate and the expected lifetime of the specific technology.

Finally, $C_{oper}$ is calculated as the sum of expenditures for the imported energy (gas and electricity, subscript $imp$) and revenues for the exported electricity (subscript $exp$):

$$C_{oper}={\sum }_Y\left(F_{imp}\left(t\right)·c_{el,imp}\left(t\right)+G_{imp}\left(t\right)·c_{gas}\left(t\right)-P_{exp}\left(t\right)·c_{el,exp}\left(t\right)\right)·h$$

(13)

where $F_{imp}\left(t\right)$, $G_{imp}\left(t\right)$, and $P_{exp}\left(t\right)$ are the power flows associated with the imported electricity, imported natural gas, and exported electricity in each hourly time step, respectively; $c_{el,imp}$, $c_{gas}$, and $c_{el,exp}$ are their specific costs/prices; and $h$ is the considered time step size (1 h).

Optimization results also include the assessment of the annual carbon emissions ($f\prime {\left({\mathbf{x}}_D,{\mathbf{x}}_O\left(t\right)\right)}_Y$ in Equation (14)) and the annual primary energy consumption of the MES ($f″{\left({\mathbf{x}}_D,{\mathbf{x}}_O\left(t\right)\right)}_Y$ in Equation (15)) due to the system operation.

$$f\prime {\left({\mathbf{x}}_D,{\mathbf{x}}_O\left(t\right)\right)}_Y={\sum }_Y\left(\left(F_{imp}\left(t\right)-P_{exp}\left(t\right)\right)·e_{el}+G_{imp}\left(t\right)·e_{gas}\right)·h$$

(14)

where $e_{el}$ and $e_{gas}$ are the average emission factors related to the electricity produced by the Italian electric system and the combustion of the natural gas, respectively.

$$f″{\left({\mathbf{x}}_D,{\mathbf{x}}_O\left(t\right)\right)}_Y={\sum }_Y\left({\displaystyle \frac{\left(F_{imp}\left(t\right)-P_{exp}\left(t\right)\right)}{{\eta }_{el}}}+G_{imp}\left(t\right)+P_{PV}(t)\right)·h$$

(15)

where ${\eta }_{el}$ is the average primary energy-to-electricity conversion efficiency of the Italian electricity system, and $P_{PV}\left(t\right)$ is the electric power produced by the PV. Note that in Equation (15), the contribution of photovoltaics to primary energy consumption is equal to the electricity produced by this technology, as it is proposed in [32].

The solution of the optimization problem in Equation (1) with the only objective of minimizing the cost of the MES (Equation (11)) falls under what is called the Cost Minimization (CM) scenario.

The fixed objective values of the reduction in carbon emissions (Equation (14)) and primary energy consumption (Equation (15)) compared to the CM scenario are imposed in specific additional runs of the optimization problem (Low Carbon Emissions—LCE scenario and Primary Energy Saving—PES scenario, respectively) to identify the optimal trade-off between cost, energy efficiency, and environmental impact.

4. Results

The optimization problem of the design and operation of the MES has been solved for each proposed configuration in the CM scenario. The following sections present the results considering the reference case as a term of comparison. The results are then compared with the other scenarios (LCE and PES) in which additional constraints are included to reduce both the primary energy consumption and CO₂ emissions.

All optimization runs were performed using the Gurobi (version 11.0.2) Python API [33] optimization software library on a standard laptop equipped with a 9th Gen i7 processor and 16 GB RAM. Execution times varied widely, ranging from 30 min to 4 days depending on the complexity of the configuration (the optimization runs involving configuration C are the most time-consuming), with tolerances in the solution on the order of 1%.

4.1. Reference Case

The reference case represents the current operation of the district. It considers all the electricity demand met by withdrawing electricity from the national grid and all the heating demand covered by gas-fired boilers.

The investment costs of the existing units and networks are not included because the capital cost of all the existing facilities is assumed to be already depreciated (all considered existing facilities are dated), thus the yearly cost of 7029.8 kEUR is only composed of the operation contribution. The electrical energy withdrawn by the national grid is 22,478 MWh/year, while the annual gas consumption is 11,960.5 MWh, resulting in a total primary energy consumption of 59,183.2 MWh/year and 8447.8 t/year of CO₂ emissions.

In the following, the electricity and heating demands of the reference case are named “reference electricity demand” and “reference heating demand”, respectively, to provide homogeneous terms of comparison for the proposed new configurations. In fact, in each of the proposed configurations, the total electricity demand also includes the electricity used by HPs and the EDN losses associated with the specific configuration, while the total heating demand also includes the DHN losses associated with the specific configuration.

4.2. Configuration A: Existing Networks

The first configuration considers the design and operation optimization of the university district MES without modifications to the topology of the pre-existing DHN and EDN.

Figure 6 shows the resulting position and size of the energy conversion and storage units. The results (

Table 5. Results of the design and operation optimization of the three configurations in the CM scenario.

	Configurations			Reference
	A	B	C
Total PV capacity [kWp]	485	485	485	0
Total CHP ICE capacity [kWel]	1705	2643	3723	0
Total GB capacity [kWth]	20,543	20,543	20,543	20,543
Total HP capacity [kWth]	2992	2448	3670	0
Overall yearly cost [kEUR/year]	6609.2	6367.9	6163.3	7029.8
Overall yearly cost [%ref]	−6.0	−9.4	−12.3	-
Operation costs [kEUR/year]	6213.3	5862.3	5467.8	7029.8
Investment costs [kEUR/year]	395.9	505.7	695.5	0.0
Imported electricity [MWh/year]	13,079.9	4371.6	1059.0	22,478.0
Exported electricity [MWh/year]	30.0	53.2	23.5	0.0
Ren 1 electricity generation [MWh/year]	704.0	704.0	704.0	0.0
NRen 2 electricity generation [MWh/year]	9835.5	18,337.6	21,936.7	0.0
HP heat generation [MWh/year]	4465.3	3487.4	4602.5	0.0
GB heat generation [MWh/year]	5125.3	4608.3	2557.3	12,594.5
CHP heat generation [MWh/year]	3027.7	4597.5	5582.2	0.0
CHP factor [%]	30.8	25.1	25.4	-
Usage of natural gas [MWh/year]	27,464.7	45,850.8	50,202.7	11,960.5
HP usage of electricity [MWh/year]	1117.0	879.2	1172.8	0.0
PE 3 consumption [MWh/year]	55,647.4	55,738.8	53,131.4	59,183.2
PE 3 consumption [%ref]	−6.0	−5.8	−10.2	-
CO2 emissions [t/year]	8955.2	10,217.3	10,176.9	8447.8
CO2 emissions [%ref]	6.0	20.9	20.5	-

¹ Ren: renewable; ² NRen: non-renewable; ³ PE: primary energy.

, column A) suggest that even if gas-fired boilers can completely cover the heating demand, it is convenient to install CHP ICEs and HPs to minimize the costs. ICEs and HPs contribute to meet the heating demand with 3027.7 MWh/year (24% of the reference heating demand) and 4465.3 MWh/year (35.5% of the reference heating demand), respectively. The installed PV capacity of 485 kW_p fully exploits the limited availability of roof area and contributes to the overall electrical generation with 704 MWh/year (3.1% of the reference electricity demand). The electricity generation from CHP ICEs instead is 9835.3 MWh/year and covers 43.8% of the reference electricity demand. No EES units are included in the solution due to the very limited PV capacity. In addition, the shape of the heating demand is well met by the installed CHP ICE and HP units, with no need for TES, while its simultaneity with the electrical demand favors cogeneration. Finally, the design optimization shows that an increase in the capacity of the existing DHN and EDN is not convenient (i.e., the existing energy networks’ capacity is sufficient to distribute the flows produced by existing and new units in the optimal operation).

The resulting overall cost, composed of the operation and investment costs, is 6609.2 kEUR/year, 6% less than the reference case also considering the amortization of investment costs for the new units, which is equal to 395.9 kEUR/year. This is because the operation costs in the optimized configuration A are 6213.3 kEUR/year, 11.6% less than the reference case. The primary energy consumption is 55,647.4 MWh/year, 6% less than the reference case, while the annual CO₂ emissions value of 9282.5 t is 6% higher than the reference case.

Higher emissions are due to the relatively low CHP factor, i.e., the ratio between the recovered waste heat and the electricity produced by the ICEs, which is equal to 30.8%. Table 5 summarizes the results.

4.3. Configuration B: Medium Integration

The second configuration allows a partial integration of the university district according to the possible network paths shown in Figure 2a. Figure 7 shows the resulting position and size of the energy conversion and storage units, and the topology of the DHNs and EDNs of the three distinct areas in which the system is divided.

The CHP ICEs’ capacity of 2643 kW_el (Table 5, column B) contributes 18,337.6 MWh/year (81.6% of the reference electricity demand) and 4597.5 MWh/year (36.5% of the reference heating demand) to the electricity and heating production, respectively, with a CHP factor of 25.1%. The optimal HPs’ capacity is concentrated into six nodes with an overall heat production of 3487.4 MWh/year (27.7% of the reference heating demand). The remaining 36.3% of the overall heat production, which includes also the DHN losses, is covered by the existing boilers, the capacity of which is not increased. The installed PV capacity is the maximum available. PV and CHP ICEs together cover 84.7% of the reference electricity demand. No energy storage capacities are included in the optimal solution.

All nodes are connected to the EDN, while the results exclude nodes 3 and 9 from the DHN. As expected, the highest capacities of the DHN and EDN branches are concentrated in the nodes where the CHP ICEs are installed and in adjacent nodes.

Finally, the results of the optimization problem show an overall cost of 6367.9 kEUR/year (9.4% less than the reference case), with operation costs of 5862.3 kEUR/year (16.6% less than the reference case) and investment costs of 505.7 kEUR/year. The resulting primary energy consumption is 55,738.8 MWh/year, 5.8% lower than the reference case, and the annual CO₂ emissions are 10,217.3 t, 20.9% higher than the reference. Both results are mainly influenced by the large amount of electrical energy generated by the CHP ICEs without taking advantage of the total available waste heat.

4.4. Configuration C: High Integration

The third configuration represents the highest integrated solution, with the whole university district served by the DHN and EDN. Figure 8 shows the topology of the networks and the resulting position and size of the energy conversion and storage units.

CHP ICEs are installed at nodes 1 and 7, with an overall capacity of 3723 kW_el (Table 5, column C), while the total capacity of HPs is 3670 kW_th. Both capacities are the highest among the three considered configurations. CHP ICEs and HPs contribute 5582.2 MWh/year and 4602.5 MWh/year, respectively, to cover 80.9% of the reference heating demand. At 21,936.7 MWh/year, the electricity generation of ICEs covers 97.6% of the reference electricity demand; thus, in this configuration, the autonomous electricity generation of CHP ICEs and PV completely fulfills the reference electricity demand. Once again, the optimal solution excludes electrical and thermal energy storage units.

Figure 9 shows the progressive expansion of both the DHN and EDN across the three layouts resulting from the optimization in the CM scenario. In configuration C, the resulting capacities of branches 6–16, 15–16, and 16–17, laid out for both the DHN and EDN, allow the three main areas of the MES to be effectively interconnected. Only node 3 and branch 18–19 are excluded by the resulting DHN, while all the available links are installed in the EDN.

The resulting overall cost of 6163.3 kEUR/year is the lowest among the three configurations, 12.3% less than the reference case, with operation costs of 5467.8 kEUR/year, 22.2% less than the reference. In contrast, the investment costs of 695.5 kEUR/year are the highest among the proposed configurations, due to the higher capacity of CHP ICEs and HPs and the greater extent of the distribution networks. Finally, this configuration also achieves the best result in terms of primary energy savings, with a consumption of 53,131.4 MWh/year, 10.2% less than the reference case, while the annual CO₂ emissions are 10,176.9 t, about the same as configuration B.

4.5. Low Carbon Emissions and Primary Energy Saving Scenarios

In the CM scenario, an increased network integration (i.e., greater interconnection between nodes) promotes the energy generation from CHP ICEs, leading the proposed configurations to perform better than the reference case, both in terms of overall costs (up to −12.3%) and primary energy consumption (up to −10.2%). On the other hand, the installation of CHP ICEs, particularly with low CHP factors, increases the CO₂ emissions of the system, due to the lower share of the renewables of this “local” system compared to that of the national grid.

Starting from these considerations, the Cost Minimization (CM) scenario is compared with two different scenarios: Low Carbon Emissions (LCE) and Primary Energy Saving (PES). Cost minimization remains the objective in both the LCE and PES scenarios, but the former also considers a cap of −17% (referred to the reference case) on the maximum CO₂ emissions, while the latter imposes a cap of −11% (referred to the reference case) on the maximum primary energy consumption and constrains the minimum hourly utilization of CHP ICE waste heat to be equal to 10% to foster the CHP operation of ICEs. The caps are obtained by minimizing the CO₂ emissions and primary energy consumption of configuration A in two separate optimization problems.

Table 6. Results of the optimization of the proposed configurations A, B, and C in the three scenarios: Cost Minimization (CM), Low Carbon Emissions (LCE), and Primary Energy Saving (PES).

	CM Scenario			LCE Scenario			PES Scenario			Reference
	A	B	C	A	B	C	A	B	C
Total PV capacity [kWp]	485	485	485	485	485	485	485	485	485	0
Total CHP ICE capacity [kWel]	1705	2643	3723	387	0	0	2334	3318	3268	0
Total GB capacity [kWth]	20,543	20,543	20,543	20,543	20,543	20,543	20,543	20,543	20,543	20,543
Total HP capacity [kWth]	2992	2448	3670	9716	8136	7874	6682	3430	2864	0
Overall yearly cost [kEUR/year]	6609.2	6367.9	6163.3	6952.4	6904.9	6908.2	6819.8	6733.2	6683.9	7029.8
Overall yearly cost [%ref]	−6.0	−9.4	−12.3	−1.1	−1.8	−1.7	−3.0	−4.2	−4.9	-
Operation costs [kEUR/year]	6213.3	5862.3	5467.8	6478.0	6531.4	6533.7	6207.6	6099.3	6077.8	7029.8
Investment costs [kEUR/year]	395.9	505.7	695.5	474.4	373.6	374.5	612.2	633.9	606.1	0.0
Imported electricity [MWh/year]	13,079.9	4371.6	1059.0	23,993.0	24,436.2	24,443.9	20,070.2	17,447.8	16,938.9	22,478.0
Exported electricity [MWh/year]	30.0	53.2	23.5	0.0	0.0	0.0	14.6	26.4	1.7	0.0
Ren. electricity generation [MWh/year]	704.0	704.0	704.0	704.0	704.0	704.0	704.0	704.0	704.0	0.0
NRen. electricity generation [MWh/year]	9835.5	18,337.6	21,936.7	480.5	0.0	0.0	3672.6	5462.0	5805.4	0.0
HP heat generation [MWh/year]	4465.3	3487.4	4602.5	10,550.9	10,611.8	10,670.5	7685.0	4390.8	3826.3	0.0
GB heat generation [MWh/year]	5125.3	4608.3	2557.3	1480	2081.5	2070.4	1480.9	3287.7	3858.9	12,594.5
CHP heat generation [MWh/year]	3027.7	4597.5	5582.2	587.4	0.0	0.0	3452.4	5016.0	5056.2	0.0
CHP factor [%]	30.8	25.1	25.4	122.2	-	-	94.0	91.8	87.1	-
Usage of natural gas [MWh/year]	27,464.7	45,850.8	50202.7	2586.5	1976.7	1966.2	9804.7	15,313.9	16,383.2	11,960.5
HP usage of electricity [MWh/year]	1117.0	879.2	1172.8	2705.0	2667.8	2675.4	1959.8	1112.4	969.1	0.0
PE consumption [MWh/year]	55,647.4	55,738.8	53,131.4	53,695.9	54,017.3	54,022.9	52,673.0	52,673.0	52,673.0	59,183.2
PE consumption [%ref]	−6.0	−5.8	−10.2	−9.3	−8.7	−8.7	−11.0 *	−11.0 *	−11.0 *	-
CO2 emissions [t/year]	8955.2	10,217.3	10,176.9	7011.6	7011.6	7011.6	7370.6	7745.2	7817.9	8447.8
CO2 emissions [%ref]	6.0	20.9	20.5	−17.0 *	−17.0 *	−17.0 *	−12.8	−8.3	−7.5	-

* Fixed objective values in the optimization procedure.

shows the results of the three scenarios.

In the LCE scenario, the installed capacity of CHP ICEs is reduced to zero (a marginal CHP ICE capacity is maintained only in configuration A) in favor of the HPs, whose capacity in configuration A is more than tripled compared to the CM scenario. Again, in this scenario, as well as in the PES scenario, PV capacity saturates roof availability (485 kW_p) while the installation of EESs and TESs is not beneficial. It is worth noting that given the limitations of PV availability and the absence of other renewable sources, a large reduction in CO₂ emissions can only be achieved by using the renewable share of the national electric generation mix for heat generation as well, thus employing HPs. Installing HPs reduces investment costs but increases operation costs, so that the overall yearly costs are the highest among the three scenarios, although lower than the reference case (up to −1.8%). HPs cover 84% of the reference heating demand in all three configurations of the LCE scenario, which results in a lower primary energy consumption compared to the reference case (between −9.3% and −8.7%). Regarding CO₂ emissions, all the configurations reach the −17% cap imposed by the scenario.

The constraints imposed on the PES scenario result in high CHP factors: 94%, 91.8%, and 87.1% in configurations A, B, and C, respectively. The installed capacity of HPs is higher in configuration A (6682 kW_th), where a reduced network integration limits the amount of cogenerated heat. In fact, CHP ICE heat generation is higher in configurations B and C (up to 46.5% more in configuration C than in configuration A). CHP ICEs’ operation results in lower overall yearly costs compared to the LCE scenario: 3%, 4.2%, and 4.9% less than the reference case in configurations A, B, and C, respectively. Primary energy consumption is 52,673 MWh/year in the three configurations, reaching the −11% cap imposed by the scenario. Finally, regarding CO₂ emissions, the best result within the PES scenario is obtained in configuration A (12.8% less than the reference case), due to the higher heat generation of HPs. In any case, compared with the CM scenario, the high CHP factors in the PES scenario entail a strong reduction in CO₂ emissions in all three configurations.

Finally, Figure 10 shows a comparison of the energy balances for configuration C of the district MES for the three considered scenarios. In the CM scenario (Figure 10a), the economic convenience of electricity generation from the CHP ICEs leads these conversion units in operation throughout the entire year to cover 97.6% of the electricity demand (green area in Figure 10a). In contrast, the operating conditions required to achieve the CO₂ emissions target of the LCE scenario (Figure 10b) result in a completely different solution, in which all the electricity is withdrawn from the national grid (except for the 3.1% produced by PV), and HPs cover 84.7% of the reference heating demand. The additional constraint on the minimum hourly utilization of CHP ICE waste heat (10%, PES scenario) leads to a CHP factor of 87.1%, and therefore, ICEs’ operation is limited to the heating period.

5. Discussion

In terms of overall annual costs and primary energy consumption, all the proposed configurations in the three different scenarios perform better than the reference case. A higher integration (i.e., greater interconnection between nodes through energy networks) is beneficial from the economic point of view, with configuration C achieving the best result in both the PES scenario and, in particular, the CM scenario, with an overall yearly cost 12.3% lower than the reference case. In the CM scenario, network integration is also associated with lower primary energy consumptions (10.2% less than the reference case in configuration C). Both trends are related to the increased use of CHP ICEs promoted by a greater extension of the DHN and EDN (note that CHP ICEs can only be installed at nodes 1, 7, and 10), which allows for a better match with the district’s energy demands. In particular, the reduction in primary energy consumption is associated with the increase in heat cogeneration: 84% more in configuration C than in configuration A (CM scenario), with only a 3% increase in HPs’ heat generation.

Network integration and CHP are also associated with a cost-effective reduction in primary energy consumption. This is evident in the PES scenario, where for the same primary energy consumption, configuration C achieves an additional 1.9% reduction in overall yearly cost compared to configuration A, reducing both operation and investment costs.

On the other hand, when the goal is a strong reduction in CO₂ emissions, the results of the LCE scenario show that the installed capacity of CHP ICEs is reduced to zero in favor of HPs. In this scenario, DHN integration in configuration B allows the reduction of the total HPs’ installed capacity (16.3% less than configuration A), resulting in an overall yearly cost reduction of 1.1% in configuration A and 1.8% in configuration B (referred to the reference case). Configuration C, however, does not benefit further from the extended network topology.

Nevertheless, a moderate reduction in CO₂ emissions is still feasible in a hybrid context of CHP ICEs and HPs, as shown by the three configurations in the PES scenario, which achieve 12.8%, 8.3%, and 7.5% reductions in CO₂ emissions compared with the reference case, respectively.

Summarizing, in the context of this case study and considering the limited roof availability for PV systems and the absence of other renewable sources, the extension of the DHN and EDN combined with the introduction of CHP ICE units is effective in terms of both overall yearly cost reduction and primary energy saving. The use of CHP ICEs, however, limits the reduction in CO₂ emissions otherwise achievable by installing only HP units. It is also true that the penalization of CHP in this sense is strongly related to the comparison with the electricity supply from the national grid coupled with heating electrification, which can rely on a much larger renewable share: 42% in 2020 (source IEA [34]) on a national basis versus 3.1% locally.

The use of gas boilers in the three different configurations deserves further consideration. Although design optimization never affects their size, in all the proposed solutions, a part of the heating demand is met by GBs (from 11.8% to 40.7% of the reference heating demand). The cost-effectiveness of heat generation from GBs is partly biased by the high installed capacity, which favors this solution over HPs; moreover, the absence of a minimum load threshold in the GB model (unlike CHP ICEs and HPs) makes the use of boilers necessary to fulfill the heating demand at very low load conditions. This is evidenced by the results of configuration A in the LCE scenario, where despite the constraint of minimizing CO₂ emissions, HPs’ heat generation cannot completely replace GBs’ heat generation (11.7% of total heat generation). In this sense, the optimization process would benefit from the implementation of a database of commercial energy conversion units, in order to consider the installation of units of discrete sizes (instead of considering size as a continuous variable), thus improving the prediction of the part-load performance of the MES.

It is worth noting that the obtained results can be easily extended also to other districts of different sizes and similar characteristics in terms of the limited availability of renewable sources (urban context) and use of the buildings (which results in an almost simultaneous demand for electricity and heating). In fact, in these cases, the inclusion of district CHP units and networks typically leads to reduced costs and primary energy consumption, as demonstrated in [25,28], and a trade-off between cost reduction and the massive installation of heat pumps must be found to limit emissions of CO₂.

Finally, cost assumptions on technologies and energy carriers play a crucial role in the optimization problem, driving the design choices and thus influencing the economic and environmental outcomes. In this regard, the data available from the bills used to calculate the average energy costs are affected by the increase in energy costs that occurred in Europe in winter 2022–2023 due to the Russo–Ukrainian war. To this end, the subsequent reduction in energy costs registered globally could lead to favorable economic conditions.

In conclusion, it is important to emphasize that the retrofit design and operation optimization approach applied in this paper is a completely general tool with broad applicability. While this study focuses on a specific case in terms of renewable energy availability, energy demands, and economic context, the approach can be applied to other districts or systems characterized by very different technical (e.g., high availability of different renewable sources) and economic constraints, yielding variable results tailored to each specific case.

6. Conclusions

Identifying effective retrofit solutions to improve existing energy systems is a key challenge to reduce energy-related costs and environmental impacts in cities. To this end, design and operation optimization tools that consider both energy conversion and storage units and energy distribution networks play a crucial role in choosing the best investment strategy.

This paper investigates the case study of a university district in Padova, proposing three new configurations of the multi-energy system (MES) with increasing possibilities for DHN and EDN interconnection and optimizes the design and operation of additional energy conversion units and energy distribution networks applying a new MILP optimization approach suggested by the same group of authors. Three retrofit design and operation optimizations are solved for each MES configuration, considering Cost Minimization (CM), Low Carbon Emissions (LCE), and Primary Energy Saving (PES) scenarios. Results are compared to search for a suitable trade-off between the economic, environmental, and energy efficiency objectives.

All the optimized configurations in the three different scenarios (CM, LCE, and PES) perform better than the reference case in terms of overall cost reduction and primary energy saving. To this end, the results show that DHN and EDN expansion combined with CHP installation are effective intervention strategies. In fact, the most integrated configuration (C), in both the CM and PES scenarios, has the highest CHP penetration (43.8% and 39.7% of cogenerated heat, respectively) and achieves the best economic results, with overall cost reductions of 12.3% and 4.9%, respectively, compared to the reference case. In addition, configuration C also achieves the best result in terms of primary energy consumption within the CM scenario, 10.2% less than the reference case. On the other hand, the operation of ICEs at low CHP factors results in high CO₂ emissions (up to 20.9% more than the reference case) in the CM scenario, and the LCE scenario shows that ICE installation is reduced to zero if a strong reduction in CO₂ emissions is required (17% less than the reference case), albeit with a negative impact on total costs that still remain lower than in the reference case (up to −1.8%). Nevertheless, a moderate reduction in CO₂ emissions is still feasible in hybrid solutions combining ICEs (operating at high CHP factors) and HPs, as evidenced by the 8.3% reduction in CO₂ emissions and the 39.5% of cogenerated heat achieved by configuration B in the PES scenario.

In conclusion, current needs for drastic emission reduction and consumption containment make the environmental objective of any intervention imperative. To this end, MES integration through expanded energy networks and cogenerating solutions that emerge from the PES scenario represent cost-effective intervention strategies capable of mitigating the environmental impact of the district while preserving the economic sustainability of the investment. Further reduction in CO₂ emissions is also possible at the expense of overall annual cost and primary energy consumption.

Future developments of this work will consider the addition of the “synthesis” problem to the current design and operation optimization problem, i.e., the addition of the variable associated with MES topology among the decision variables of the optimization problem. Moreover, the performance prediction of the energy conversion and storage units will be improved by implementing a database of commercial units of known sizes. Finally, a cost sensitivity analysis will be conducted to determine the influence of cost assumptions.