**1. Introduction**

Cogeneration of useful heat and electrical power for an urban district or a cluster of buildings is a technically mature, environmentally-friendly and cost-effective solution, supported by the European Union Directive [1] on energy efficiency, together with the use of renewable energy sources (RES). Indeed, the Directive [2] on energy performance of buildings indicates four high-efficiency technologies, whose feasibility should be evaluated prior to construction of any new building: (a) decentralized energy supply systems based on RES; (b) cogeneration; (c) district or block heating or cooling; and (d) heat pumps.

In view of this, several different configurations of distributed energy systems (DES) have been investigated in recent years, mostly focusing on cogeneration units and RES technologies. Pagliarini and Ranieri [3] studied the effectiveness of thermal storage coupled with a cogeneration engine to satisfy the energy requirements of a university campus, and stressed the importance of the sizing of the storage. Bracco et al. [4] dealt with the topic of distributed generation, by presenting the University of Genoa polygeneration microgrid, which is based on RES and cogeneration units. A tool for the optimal integrated design and operation of a trigeneration system serving a cluster of buildings was proposed by Piacentino et al. [5,6]. The optimal design and operation of a hybrid renewable energy system based on an internal combustion engine and photovoltaic panels was investigated by Destro et al. [7].

In addition to the synthesis and design problems, the optimal operational strategy of energy microgrids has also received considerable attention. Indeed, the adoption of smart control techniques can significantly improve the economic and environmental performances of those systems [8]. For example, Roldán-Blay et al. [9] developed an algorithm for the optimal management of a RES-based electric microgrid. Similarly, Phan et al. [10] investigated schedule strategies to minimize the operating cost of a building energy system with photovoltaic panels and a wind micro-turbine. Asaleye et al. [11] proposed a decision-support tool that identifies the optimal operation of renewable energy microgrids by considering forecast of climate variables.

All the above-mentioned works show the energy, environmental and economic effectiveness of cogeneration systems and renewable energy technologies in smart energy grids. Nevertheless, some unresolved issues remain. Indeed, the intermittent and stochastic nature of RES limits their use, and cogeneration units are characterized by a fixed power-to-heat ratio, thus, they fail to match both fluctuating thermal and electric demands. In this context, electrically driven heat pumps may represent an interesting solution due to their ability to shift thermal loads into electric loads. Moreover, heat pumps are a mature and efficient technology, and they are especially suited to the implementation of smart control strategies [12].

For those reasons, the present work discusses a novel configuration for smart multi-energy microgrids, which consists of distributed energy units and a centralized cogeneration unit feeding a micro-district heating network. Specifically, we investigate the benefits of integrating reversible heat pumps for heating and cooling purposes at the building level. The heat pumps represent an interconnection between the electricity and heating networks, therefore, they can be used to increase the operational flexibility of the microgrid and support the integration of renewable energy technologies, i.e., wind turbine, photovoltaic panels, and solar thermal collectors.

The paper is structured as follows. Section 2 presents the design and modeling methodology of the smart multi-energy microgrid. Sections 3 and 4 present the case study and the optimization problem and methodology, which are used in Section 5 to compare the proposed configuration using distributed heat pumps, with a more conventional solution that employs a centralized CHP (Combined Heat and Power) system and natural-gas boilers. Finally, Section 6 presents the concluding remarks.

## **2. Energy System Overview and Modeling**

In this work, we refer to multi-energy microgrids of medium dimensions with different buildings and loads, using an integrated thermal and electrical energy production system fed by traditional and renewable sources to concurrently satisfy various services (heating, cooling, electrical energy and domestic hot water). The considered generators are: (i) a CHP consisting in an internal combustion engine (ICE); (ii) natural gas boilers; (iii) heat pumps and chillers; (iv) solar thermal collectors; (v) wind turbines; and (vi) photovoltaic modules. Thermal storage is also considered. Figure 1 shows a simplified classification scheme of the reference energy system.

As is well-known, the traditional design approach based on a separate analysis of each component represents a suboptimal design method for multi-energy systems [13]. The so-called simulation-based optimization methods are the most recognized procedures to investigate the best synthesis, sizing and control of integrated system through the simulation of the operative performances. Therefore, in the following sub-section we present the operative dynamic model of each block listed in Figure 1.

**Figure 1.** Schematic of the energy system.

#### *Modeling of the System Components*

The components models must be based on a proper trade-off between the accuracy of the results and computational effort. The latter feature is essential to allow their employment within the optimization procedure to identify the most efficient design and operation of the smart multi-energy microgrid.

The ICE is modeled through performance curves taken from [14], which provide the thermal and electric efficiency, *ηel* and *ηth*, respectively, as a function of the engine load factor *LICE*. The electric (*EICE*) and thermal (*QICE*) output are evaluated based on the load factor, *LICE*, and the ICE nominal electric power capacity (*Enom ICE* ), as follows:

$$\begin{cases} \begin{array}{c} E\_{ICE} = E\_{ICE}^{nom} \cdot L\_{ICE} \\ Q\_{ICE} = \left[ \eta\_{th}(L\_{ICE}) / \eta\_{el}(L\_{ICE}) \right] \cdot E\_{ICE}^{nom} \cdot L\_{ICE} \end{array} \tag{1} \end{cases} \tag{1}$$

The boilers can be modeled with a constant efficiency (*ηB*) over their whole operating range. Photovoltaic panels are simulated through the model provided by [15], which considers the PV performance as a function of the solar irradiance, PV characteristics and cell array temperature:

$$E\_{PV} = n\_{PV} S\_{PV} \eta\_{PV} \eta\_{inv} I\_{sol,PV} \tag{2}$$

$$\eta\_{PV} = \eta\_{PV,ref} \left[ 1 - \beta\_{T,PV} \left( T\_{PV} - T\_{PV,ref} \right) \right] \tag{3}$$

$$T\_{PV} = T\_{\text{ext}} + (219 + 819K\_l) \frac{\text{NOCT} - 20}{200} \tag{4}$$

The thermal performances of the ST are evaluated through the classical model illustrated in [16], based on the characteristics of the panel in terms of transmittance and absorptance factors for normal irradiance (< *τα* >*n*), removal factor (*Fr*) and frontal losses (*Ul*). The equations read:

$$\eta\_{ST} = F\_r(\tau a)\_n \left[ 1 - b\_0 \left( \frac{1}{\cos \theta} - 1 \right) \right] - \frac{F\_r l l\_L (T\_{ST,in} - T\_{ext})}{I\_{sol, ST}} \tag{5}$$

$$Q\_{ST} = n\_{ST} S\_{ST} \eta\_{ST} I\_{sol, ST} \tag{6}$$

The electrical power generated by the wind turbine varies as the cube of the wind speed, between a cut-in speed and a nominal speed; the latter corresponds to the nominal electrical power generated by the wind turbine. The nominal electrical power is generated between the nominal wind speed and a cut-out speed. The equations, in accordance with [17], read:

$$E\_{WT} = \begin{cases} 0, & w \langle w\_{cut-in} \text{ or } w \rangle w\_{cut-out} \\ & kw^3, & w\_{cut-in} \le w \le w\_{nom} \\ kw^3\_{nom}, & w\_{nom} \le w \le w\_{cut-out} \end{cases} \tag{7}$$

where *k* is a coefficient that depends on the characteristic curve of the generator.

The internal energy variation in the thermal storage is calculated considering the thermal fluxes provided to the water volume from all the connected generators and the heat delivered to the load:

$$V\_{TS} \rho\_W c\_W \Delta T\_{TS} = \sum\_i Q\_{TS, in, i} - \sum\_j Q\_{TS, out, j} - Q\_{TS, ls} \tag{8}$$

where the heat losses of the storage tank are evaluated as

$$Q\_{TS,ls} = UA\_{TS}(T\_{TS} - T\_{ext,TS}) \tag{9}$$

Reversible heat pumps and chiller performance are evaluated by means of the so-called second-law efficiency [18]. The method reads:

$$\text{COP} = \eta\_H^{II} \cdot \text{COP}\_{\text{id}} = \eta\_H^{II} \frac{T\_{\text{cond}}}{T\_{\text{cond}} - T\_{\text{eva}}} \tag{10}$$

$$EER = \eta\_{\mathbb{C}}^{II} \cdot EER\_{id} = \eta\_{H}^{II} \frac{T\_{cva}}{T\_{cond} - T\_{cva}} \tag{11}$$

where *COPid* and *EERid* are the coefficients of performance of a reversed Carnot cycle operating between the source and sink temperatures. According to manufacturers, both *ηI I <sup>H</sup>* and *<sup>η</sup>I I <sup>C</sup>* can be assumed as constant.

The generators are connected to the thermal storages and the buildings through a district heating network (DHN), whose heat losses can be modeled as follows:

$$Q\_{DHN,ls} = \mathcal{U}\_{DHN} L\_{DHN} \left( T\_{\text{avg},DHN} - T\_{\text{grround}} \right) \tag{12}$$

Finally, the heating/cooling loads of the buildings can be evaluated through a model that correlates the sol-air temperature [19] with the energy load of the building, based on the standard EN 15306 [20]. This model is further improved by considering the effect of the building thermo-physical properties in shifting the influence of the external climate on the heating/cooling load.

$$Q\_{th,H/\mathbb{C}} = P\_{H/\mathbb{C}} \left( 1 - \frac{\overline{T\_{\text{ext}}^\*} - T\_{\text{des},H/\mathbb{C}}}{T\_{off,H/\mathbb{C}} - T\_{\text{des},H/\mathbb{C}}^\*} \right) \tag{13}$$

$$\overline{T^\*\_{ext}}(t) = \frac{1}{\overline{\Phi}} \sum\_{i=0}^{\overline{\Phi}} T^\*\_{ext} \left( t - \overline{\Phi} + i \right) \tag{14}$$

$$\overline{\Phi} = \sum\_{i} \frac{(lIA)\_{i} \phi\_{i}}{[\sum\_{i} (lIA)\_{i} + H\_{\text{ref}}]} \tag{15}$$

$$T\_{ext}^\*(t) = T\_{ext}(t) + \frac{\alpha\_S}{h\_\mathcal{e}} I\_{sol}(t) \tag{16}$$

Further details on Equations (13)–(16) can be found in [21]. This model represents a good trade-off between simplified models (e.g., the energy signature method [22]), which simply correlates external temperature and heating/cooling load, and dynamic models (e.g., TRNSYS or EnergyPlus), which include the building inertia characteristics, solar radiation, and internal loads, providing more accurate results, but requiring a detailed knowledge of the building envelope and heat gain profiles.

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

In this work, we refer to an integrated energy system serving a hypothetical campus, located in Trieste, Italy. This city has a favorable climate, where RES (solar thermal, photovoltaic modules, wind turbines) can provide a significant amount of energy. The Italian Thermotechnical Committee (CTI) provides hourly profiles of external temperature, global solar irradiance on the horizontal plane and wind speed [23]. The monthly-average values of the external temperature and irradiance on horizontal plane are reported in Figure 2.

**Figure 2.** Average monthly temperature and daily irradiance on horizontal plane.

The campus is located far enough from the city to avoid airflow obstructions and shading, which would reduce the renewable energy share. As a whole, the 1000-student campus occupies a surface of <sup>1</sup> × 0.5 km2 and includes five dormitories, a dining hall, a gym, a students' center with classrooms and administrative offices. A schematic representation of the campus is shown in Figure 3.

**Figure 3.** Scale representation of the campus: buildings, district heating, and generation systems.

All the buildings have a similar structure in terms of thermal transmittance of the walls, roofs, floors, and windows (0.29 W/m2·K, 0.19 W/m2·K, 0.26 W/m2·K, 1.80 W/m2·K, respectively). Specific profiles of internal gains and electrical energy requirements were chosen for each building type, according to their use, periods of presence, and appliances. The internal gains and electrical energy requirements have lower values during weekends and holidays. The buildings have different terminal units that require two supply temperature levels (i.e., low and medium). The power peak and total energy for the heating, cooling, domestic hot water (DHW) and electrical energy services are reported in Table 1.

**Table 1.** Peak values and energy needs for the four services of the campus and for the five types of buildings.


In this work, we compare a "centralized" configuration (see Figure 4a), in which the DHN is fed by a cogeneration unit and a centralized gas boiler, with a "distributed" one (see Figure 4b), in which reversible heat pumps are installed in the buildings. In both configurations, the sizing of the solar and wind generators remains the same, while an optimization analysis is performed for the CHP and heat pumps, together with the optimal control strategy (see Section 4).

**Figure 4.** Analyzed configuration for the reference multi-energy microgrid: (**a**) centralized, (**b**) distributed.

The generators concur to satisfy the heating/cooling/DHW/electrical energy requirements of the campus with the following strategy:


Table 2 reports all the parameters of the analyzed smart multi-energy microgrid.


**Table 2.** Parameters used in the analysis.

#### **4. Optimization Problem and Methodology**

In this work, we aim to compare the centralized and the distributed configurations defined in Section 3, from an economical point of view. The selected performance index is the annual total cost, *TC*, defined as:

$$TC = \frac{INV}{t\_{life}} + O\&M + OC\tag{17}$$

where *INV*/*tlif e*, *O*&*M*, and *OC* are the yearly capital, maintenance, and energy operational costs, respectively, and *tlif e* is the considered lifetime of the microgrid (i.e., 20 years). Since we are making a comparative analysis, *TC* only includes the costs that differ in the two configurations, i.e., the purchased equipment cost (PEC) for the ICE unit, reversible HPs, chillers, and boilers, the associated operations and maintenance costs (O&M), and the net cost of the energy purchased from the gas and power grids. We do not consider the sizing of the RES technologies, which are assumed to have the same design and energy production in both the distributed and centralized configurations; therefore, they are not included in the economic analysis and optimization process. The terms in Equation (17) read:

$$INV = PEC\_{ICE} + PEC\_{HP} + PEC\_{C} + PEC\_{B} \tag{18}$$

$$\text{O}\&M = \text{O}\&M\_{ICE} + \text{O}\&M\_{HP} + \text{O}\&M\_C + \text{O}\&M\_B\tag{19}$$

$$\text{OC} = \sum\_{i}^{8760} \text{OC}^{i} = \sum\_{i=1}^{8760} \left( c\_F^i F\_B^i + c\_F^i F\_{ICE}^i + c\_{el,P}^i E\_P^i - c\_{el,S}^i E\_S^i \right) \tag{20}$$

The timestep length adopted for the energy system simulation is one hour. The investment and maintenance cost functions are presented in Table 3. The costs of HPs, chillers and boilers were obtained through a linear regression of actual manufacturers' data.



According to the operating strategy described in Section 3, both the sizing and control optimization can be written as a function of the ICE nominal electrical capacity, *Enom ICE* , and load factor profile, *<sup>L</sup><sup>i</sup> ICE*.


### *4.1. Operational Optimization Problem*

The operational optimization problem consists in identifying the scheduling of the generators that meets the energy demands at minimum cost (i.e., cost for purchasing electricity from the grid, income for selling electricity to the grid, cost of natural gas). The operational optimization problem is therefore defined as the minimization of the total annual energy cost.

$$\min\{\text{TC}\} = \min\left\{\sum\_{i=1}^{8760} c\_{\text{F}}^{i} F\_{\text{B}}^{i} + \sum\_{i=1}^{8760} c\_{\text{F}}^{i} F\_{\text{ICE}}^{i} + \sum\_{i=1}^{8760} + c\_{\text{el,P}}^{i} E\_{\text{P}}^{i} - \sum\_{i=1}^{8760} c\_{\text{el,S}}^{i} E\_{\text{S}}^{i} \right\} \tag{21}$$

where *i* = 1, ... , 8760 timesteps, *cF* is the fuel price (0.04 €/kWh), and *cel*,*<sup>P</sup>* and *cel*,*<sup>S</sup>* are the prices of purchased and sold electricity (0.18 and 0.04 €/kWh), respectively. Consequently, the following decision variables are considered:

$$E^i\_{\rm ICE}, Q^i\_{\rm ICE}, Q^i\_{\rm B,j}, Q^i\_{\rm HP,IT,k}, Q^i\_{\rm HP,MT,k'}, C^i\_{\rm C,LT,k'} C^i\_{\rm C,MT,k'} E^i\_{\rm S'} E^i\_{\rm P}. \tag{22}$$

and demand constraints and balance equations and inequalities are defined as follows

$$Q\_{\rm DHN,HT,j}^{l} + Q\_{\rm B,j}^{l} - Q\_{\rm HT,netD,j}^{l} = 0 \tag{23}$$

$$Q^i\_{DHN,MT,k} + Q^i\_{HP,MT,k} - Q^i\_{MT,D,k} = 0 \tag{24}$$

$$Q^i\_{DHN,LT,k} + Q^i\_{HP,LT,k} - Q^i\_{LT,D,k} = 0 \tag{25}$$

$$Q^i\_{ICE} - Q^i\_{DHN,ls} - \sum\_j Q^i\_{DHN,HT,j} - \sum\_k Q^i\_{DHN,MT,k} - \sum\_k Q^i\_{DHN,LT,k} \ge 0 \tag{26}$$

$$\mathsf{C}\_{\mathsf{C},LT,k}^{i} - \mathsf{C}\_{LT,D,k}^{i} = 0 \tag{27}$$

$$\mathbf{C}\_{\mathbf{C},MT,k}^{i} - \mathbf{C}\_{MT,D,k}^{i} = 0 \tag{28}$$

$$E\_P^i - E\_\mathcal{S}^i - \sum\_k E\_{HP,k}^i - \sum\_k E\_{\mathcal{C},k}^i + E\_{ICE}^i + E\_{PV}^i + E\_{WIND}^i - E\_D^i - E\_{\text{aux}}^i = 0 \tag{29}$$

To solve the optimal operation problem, an ad-hoc dispatch strategy algorithm has been developed, based on the following considerations:


Indeed, in the energy system under investigation, the overall optimum coincides with the sum of optimums of every single timestep, since the behavior of the TSs (Thermal Storages) linked to the solar thermal is independent from the operational control. For this reason, as already shown in [29], the overall operational problem can be split into 8760 subproblems, one for each timestep, and the problem can be considered "static". Moreover, the DHN must always satisfy with higher priority the medium-temperature heat demand *Q<sup>i</sup> MT*,*D*, as opposed to the low-temperature heat demand *<sup>Q</sup><sup>i</sup> LT*,*D*, since the HPs operate with higher COP at lower supply temperatures. Therefore, only the three following combinations must be evaluated:


For each possible dispatch priority, once the *LICE* is set and the amount of electricity and heat produced by the ICE is defined, the thermal losses and the net amount of heat available at the DHN are known from Equations (23) and (24). Then, Equations (25)–(27) state that the boiler and the HP production must meet the remaining heat demand, if any. Furthermore, Equations (27) and (28) require that the electrical chillers or reversible heat pumps meet the chilled water demand, and Equation (29) defines the electrical energy exchange with the grid. Therefore, as mentioned, the nine decision variables are bound to each other and the problem is conveniently reduced to finding the dispatch priority order and the optimal ICE load factor *LICE* that minimize the cost of energy at each timestep. Eleven discrete values of *LICE* have been considered, and an exhaustive search algorithm was adopted to identify the optimal solution, among all the possible combinations. This allows the development of a low computational-cost algorithm, compatible with the need for a quick response for real-time implementation [30] and further advanced analyses (e.g., optimal design and uncertainty analysis, as in [29]).

#### **5. Results and Discussion**

Figure 5 shows the total annual costs depending on the nominal electrical capacity of the CHP unit. We note that the optimal CHP size is practically the same for both of the configurations (75 kWel), with a total cost reduction of about 8% for the distributed solution. Economic details are presented in Table 4. The size of all generators are shown in Table 5.

(**a**) Distributed configuration. (**b**) Centralized configuration.

**Figure 5.** Installation, energy, and maintenance annual costs using the optimal control strategy.


**Table 4.** Annual costs of the optimal distributed and centralized configurations, k€/yr.

**Table 5.** Optimal sizes for the distributed and centralized configurations, kW.


Figure 6 shows some examples of the heat and power profiles resulting from the optimization procedure in three weeks of the year for the distributed configuration. In addition to the areas and lines explained in the chart legend, we specify that the white area under the blue curve in the thermal plots represents the thermal overproduction by the CHP unit; the white area under the green, orange, or yellow lines in the electricity plots quantifies the electrical energy sold to the grid.

(**a**)1st week of January (**b**)1st week of January
