Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Enhancing Trust in Transactive Energy with Individually Linkable Pseudonymous Trading Using Smart Contracts
Previous Article in Journal
Ambient Temperature Effects on Energy Consumption and CO2 Emissions of a Plug-in Hybrid Electric Vehicle
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Optimal Operation of an Industrial Microgrid within a Renewable Energy Community: A Case Study of a Greentech Company

1
Department of Electrical, Electronics, and Telecommunication Engineering and Naval Architecture, University of Genoa, Via Opera Pia 11a, 16145 Genoa, Italy
2
Department of Energy Technologies and Renewable Sources, Italian National Agency for New Technologies, Energy and Sustainable Economic Development, 00196 Rome, Italy
*
Author to whom correspondence should be addressed.
Energies 2024, 17(14), 3567; https://doi.org/10.3390/en17143567
Submission received: 14 June 2024 / Revised: 5 July 2024 / Accepted: 16 July 2024 / Published: 20 July 2024
(This article belongs to the Section A1: Smart Grids and Microgrids)

Abstract

:
The integration of renewable energy sources in the European power system is one of the main goals set by the European Union. In order to ease this integration, in recent years, Renewable Energy Communities (RECs) have been introduced that aim to increase the exploitation of renewable energy at the local level. This paper presents an Energy Management System (EMS) for an industrial microgrid owned and operated by a greentech company located in the north of Italy. The company is a member of an REC. The microgrid is made of interconnected busbars, integrating photovoltaic power plants, a fleet of electric vehicles, including company cars and delivery trucks supporting Vehicle-to-Grid (V2G), dedicated charging stations, and a centralized battery energy storage system. The industrial site includes two warehouses, an office building, and a connection to the external medium-voltage network. The EMS is designed to optimize the operation of the microgrid and minimize the operating costs related to the sale and purchase of energy from the external network. Furthermore, as the company is a member of an REC, the EMS must try to follow a desired power exchange profile with the grid, suggested by the REC manager, with the purpose of maximizing the energy that is shared within the community and incentivized. The results demonstrate that, when minimizing only costs, local self-consumption is favored, leading to a Self-Sufficiency Rate (SSR) of 65.37%. On the other hand, when only the adherence to the REC manager’s desired power exchange profile is considered in the objective function, the SSR decreases to 56.43%, net operating costs increase, and the energy shared within the REC is maximized.

1. Introduction

The Green Deal Industrial Plan strengthens the competitiveness of Europe’s net-zero industry and accelerates the transition to climate neutrality by creating a more favorable environment for increasing European Union (EU) production capacity for the net-zero technologies and products needed to meet Europe’s ambitious climate goals [1]. This approach intends to make factories greener by the implementation of a sustainable energy supply to warehouses, office buildings, and industrial sheds. The industry sector, as one of the highest energy consumers in the EU, accounted for 25.1% of final energy consumption in 2022. Within the industrial sector, the highest energy consumers in the EU in 2022 were the chemical and petrochemical industry, the non-metallic minerals industry, and the paper, pulp, and printing industry [2]. Electricity (33.3%) and natural gas (31.2%) accounted for nearly two-thirds of the final energy consumption in the EU’s industry sector in 2022 [2]. In addition, the building sector is one of the most energy-consuming sectors, representing, at the EU level, 40% of the final energy consumption and 36% of Greenhouse Gas (GHG) emissions. Electricity consumption represents around 35% of buildings’ energy needs [3,4]. To reach the aforesaid challenging goals, EU countries have put in place several measures in recent years to reduce energy consumption in buildings and increase the exploitation of Renewable Energy Sources (RESs) with the so-called “Fit for 55” package that aims to reduce GHG emissions by at least 55% by 2030, with respect to 1990, and increase the share of renewable energy in the market by up to 40%.
In the industrial sector, buildings, warehouses, and halls constitute large energy consumers with a high demand for electricity, heating, and cooling to guarantee both an adequate level of comfort for workers and the proper conservation of goods (e.g., in the case of warehouses for food products). The installation of RES plants close to the aforesaid facilities can reduce their carbon footprint. Photovoltaic (PV) modules, either installed on the roofs, integrated with the facades of buildings, or adopted in parking lots, as well as small wind turbines and hydro plants, can also provide sustainable energy to industrial sites. On the other hand, to cover the thermal needs (heating, steam, hot water) of industrial facilities in a sustainable way, thermal solar collectors, geothermal heat pumps, and boilers fed by renewable sources are usually adopted, together with Combined Heat and Power (CHP) units [5] comprising microturbines and internal combustion engines fed by natural gas and biogas, and also by hydrogen produced by electrolyzers directly fed by RES plants. When cooling energy is also needed, Combined Cooling, Heat, and Power (CCHP) plants [6] are typically installed by coupling one or more CHP units with absorption chillers. Nowadays, several industrial sites are converting to green sites by adopting many of the above technologies and transforming their buildings into real prosumers, as highlighted in [7], where a multi-energy industry facility offering flexibility services in the electricity markets is analyzed. As reported in the scientific literature, and as illustrated in Figure 1, more and more papers focus on “green warehouses” and “sustainable warehouses”, proposing mathematical models and tools to optimize and simulate the design and operation of warehouses.
Battery Energy Storage Systems (BESSs) are usually installed at green industrial sites in order to mitigate the stochasticity of RES sources and loads, and to exploit the variability of electricity prices. RESs and BESSs can also provide ancillary services to the network and jointly work with the electric mobility infrastructure, as discussed in [8], where the example of a beverage company in northern Germany with electrical industrial forklifts, power-to-heat plants, wind turbines, PV systems, and energy storage devices (thermal, electrical, and hydrogen) is described and analyzed in detail. Moreover, RES generation and BESS systems can be exploited not only for the provision of active power but also for the exchange of reactive power to satisfy the reactive power demand of loads and to provide voltage support, even at night [9]. All these technologies and best practices are now commonly used in new logistics parks, as reported in the literature. For instance, the interaction between green logistics practices, the circular economy, and Industry 4.0 technologies is explored in [10], while a detailed review of key green practices in warehouses is proposed in [11]. In [12], the advantages of installing Microgrids (MGs) with wind, solar and battery technologies in sustainable factories and warehouses are illustrated, while in [13] the authors report on a new methodology to optimally locate and size wind turbines, solar power systems, and batteries in industrial facilities, with the aim of reducing net supply-chain costs. A similar point of view is adopted by the authors in [14], where the exploitation of the local production of energy from wind and solar sources is seen as a way to decarbonize manufacturing, transportation, and warehousing operations. On the other hand, in [15], a multi-objective optimization model is defined in order to minimize the global costs of an industrial site acting as a prosumer within the electricity market. As mentioned above, it is important to note that many industrial prosumers have MGs within their sites. In [16], the concept of industrial MGs is defined and real-world applications are described, with the aim of illustrating the main power control and energy management strategies in this context. In [17] an industrial MG with wind turbines, CHP units and BESS systems is analyzed by also including the possible application of demand response strategies.
While MGs are purely physical systems and can be installed in industrial contexts to provide a local and reliable energy supply, it should not be forgotten that green industrial sites can now also be an active part of virtual power plants and energy communities. In the first case, we deal with aggregated systems of energy assets remotely and automatically optimized to dispatch services for distribution or wholesale markets [18], while in the second case we mainly refer to Renewable Energy Communities (RECs) [19] as aggregated legal entities made of members (natural persons, small and medium enterprises, local authorities, etc.) that decide to jointly share energy locally produced by an RES in order to obtain environmental, economic and social community benefits, rather than final ones [20,21]. Many examples of energy communities can be found in the residential and tertiary sectors [22,23], but an excellent solution to maximize shared energy is also to bring together residential and industrial members. The latter is usually characterized by peaks in electricity demand at central times of the day, coinciding with peaks in PV production. Indeed, in [24] the authors prove that an REC can achieve economic profitability, including the industrial demand of a desalination unit, while the study reported in [25] highlights that the industrial sector presents a good opportunity to install RESs and provide renewable energy to closed city members. On the other hand, demonstrates how the aggregation of different types of distributed storage systems in RECs is a way to maximize the self-consumption of the community and provide ancillary services to the electric power system [26]. In [27], an optimal planning approach for renewable energy communities is proposed. A similar approach is presented in [28], where demand response with Electric Vehicles (EVs) is also considered. In [29], the impact of RECs on the grid is investigated.
Several considerations can also be made in relation to environmental issues at industrial sites, where not only office buildings, warehouses and halls are responsible for environmental emissions. A significant amount of emissions is due to transportation, i.e., the means of transport used to move goods within industrial sites, the vehicles (vans and trucks) used to deliver freight or transport goods to other destinations, and the vehicles used to load and unload trains and ships. As set out in the European Commission’s latest regulations, a gradual and massive reduction in CO2 emitted by the transport sector is planned. Regarding light-duty vehicles, the regulations indicate a 50% reduction by 2030, with zero emissions by 2035, while for heavy-duty vehicles, a progressive reduction, of 45% by 2030, 65% by 2035, and 90% by 2040, is envisaged [30]. As reported by the Italian Association MOTUS-E, nowadays, in many EU countries, light commercial EVs constitute less than 10% of annual new registrations, and even smaller percentages are found for electric trucks [30]. The proliferation of the EV infrastructure, coupled with RES plants used to directly charge them, will have an impact on public utilities and buildings [31], and also on the industrial sector, where many companies still use highly polluting vehicles, thus improving air quality in urban areas. Electric mobility infrastructures have a significant impact on the electric power system, considering the current deployment of quick, fast, and ultra-fast charging stations, which can lead to maximum charging power values of around 43 kW in alternating current, and greater than 300 kW in direct current systems, for each station. The increase in charging power allows for a reduction in charging times, which is particularly important in industrial applications, such as logistics, where EVs need to be parked for short periods during the day and therefore require fast charging; at the same time, they can also be charged at low/medium power at depots, typically during overnight charging. In any case, an increase in the number of EVs will mean an increased deployment of charging infrastructure, both in terms of the number of chargers and the total power supplied, which will have a significant impact on electricity distribution networks. In order to maximize the potential of electricity grids without creating congestion and inefficiencies, it is necessary to define intelligent EV smart-charging strategies and promote an optimized integration of e-mobility infrastructure and RES systems [32]. This implies the modulation of the charging power as a function of different factors, such as electricity prices, the availability of RES, current operating conditions, and the vulnerabilities of power networks. Furthermore, by implementing Vehicle-to-Everything (V2X) technologies, namely, Vehicle-to-Grid (V2G) and Vehicle-to-Building (V2B), EVs can truly take an active role in the electricity system and behave like BESSs [33]. They can be charged when electricity purchase prices are lower, the grid is not overloaded, or when there is a surplus of energy available from an RES; conversely, they can supply electricity to the grid (V2G) or to a building (V2B) when electricity purchase prices increase, or when the grid needs additional power to cover the decrease in RES production or the lack of energy supply from conventional power plants. This helps to reduce the variability of energy production by RESs and to increase the self-consumption of renewable energy, which is also beneficial from an economic point of view. The importance of these recent innovations in the field of electric mobility can also be seen in the literature, where an increasing number of scientific articles on the subject can be found, as shown in Figure 2.
There are several papers in the scientific literature dealing with the optimal design and operation of MGs, nanogrids, energy communities, and prosumer buildings, but only a few of them refer to industrial applications, and, in particular, to logistics. These types of decision problems are typically formalized by mathematical programming models, usually Mixed Integer Linear Programming (MILP) models with different time horizons, an equivalent year or several years, for optimal design problems, and a day or a week for optimal operation problems, the latter being typically used to define an Energy Management System (EMS). There are also some papers that combine optimal design with optimal management, where the main objective is to optimally size the energy system, but that provide, at the same time, some suggestions on how to operate it on a daily basis for some typical days. In general, MILP models are considered to be a good compromise for representing real system behavior well, albeit with some approximations in the representation of the physical relationships describing the different technologies, without greatly increasing the computational effort. In [34], a day-ahead EMS for the optimal operation of an MG with renewable generation and CHP units is presented, taking into account reactive power flows and voltage phenomena. In [35], a methodology to optimally design a multi-vector energy hub created to supply electricity and heat to a set of buildings within a sustainable district acting as a local energy community is described, while in [36] the authors describe a stochastic optimization model to optimally manage the energy infrastructure of the Savona Campus, seen as a local energy community, while considering innovative bidding strategies within the electricity markets. The impact of electric mobility on the optimal design of renewable energy collective self-consumers is investigated in [37], whereas in [38], a methodology is proposed to optimally design a local energy community that consumes energy from an MG with RES and storage systems. The issue of disordered EV charging demand and its impact on distribution networks is investigated in [39], where the authors propose a bi-level optimal dispatching scheme for a community charging station, considering time-of-use and real-time pricing mechanisms. The integration of EVs with power-to-gas technologies is investigated in [40], where a scheduling model based on demand response is proposed. A smart energy infrastructure integrating different types of batteries and EV charging stations is analyzed in [41], and, similarly, in [42], the focus is on the integration of V2G with stationary batteries. In [43], several EMS and advanced energy control strategies are proposed to investigate the potential of the interaction between buildings and EVs, while in [44], the focus is on EV smart-charging strategies in facilities characterized by the integration of RES and electric mobility. With regard to V2B applications, it is worth mentioning the study reported in [45], where the benefits of EVs acting as storage systems in buildings operating in island mode are evaluated, and the analyses carried out in [46,47] that investigated how V2B can improve the energy and economic key performance indicators of a prosumer building.
The research described in the present paper contributes to the growing field of REC optimization and provides a practical example of integrating advanced energy management technologies in the industrial sector. In particular, an EMS for an industrial MG, owned and operated by a greentech company located in the north of Italy, is presented. The peculiarity is that the site where the MG is installed is the only industrial member of an REC whose other members are only residential users, some of them just consumers and some of them producers. The Low-Voltage (LV) MG has a ring topology with only renewable power plants, namely PV modules, a centralized BESS, and some charging stations for a fleet of EVs, including company cars and delivery trucks able to exploit V2G. As far as the electrical demand of the site, the loads are represented by two warehouses, an office building and the EVs when connected to the dedicated chargers. The site has a Medium-Voltage (MV) connection to the public network with which it can exchange energy in both directions (absorption and injection). The EMS has been conceived in order to optimize the daily operation of the MG, trying to minimize a multi-objective function. On the one hand, it must minimize the net operating costs related to the sale and purchase of energy from the external network; on the other hand, the EMS should follow a desired power exchange profile with the external network, provided by the REC manager in order to maximize the energy shared within the REC, which is incentivized. The aforesaid profile can be the result of another EMS, the one used by the REC manager not described here, which aims to maximize the incentives for the REC members. A further objective of the industrial MG EMS is to limit the inductive or capacitive reactive power exchange with the network in particular time periods, in order to avoid overvoltages or undervoltages in the distribution system; hence, reactive power exchange penalties are included in net operating costs. The multi-objective function that characterizes the EMS in this paper has been defined precisely to take into account two aspects: achieving the objectives of the industrial user, who typically wants to minimize his energy bill, and maximizing the REC incentives. Thus, the EMS was developed to assess whether or not these objectives are in conflict.
A strong novelty of the paper is the introduction of a multi-objective function that considers the trade-off between the interest of the REC to maximize the shared energy and the interest of the single REC member that, in principle, is mainly the minimization of its own net operating costs. In addition, the case study is innovative, considering an industrial user member of an REC, not only exploiting common RESs and BESS but also the V2G application of electric trucks.
The results of the EMS application demonstrate significant improvements in energy efficiency and cost savings for the industrial user, highlighting the potential benefits for the inclusion of similar members within RECs. The EMS takes into account also the capability curves of inverters, that are inherently non-linear: therefore, the EMS is formalized as a Mixed Integer Quadratically Constrained Quadratic Programming (MIQCQP) problem.
The paper is organized as follows. In Section 2, the MILP mathematical models are described by reporting the multi-objective function and all the constraints and decision variables used to model the MG. Section 3 presents all the relevant technical and economic data concerning the installed technologies and the operation parameters. In the same section, the main results of the analysis are discussed for the two main scenarios, i.e., the one where we only minimize the operating costs of the MG and the one where only the deviation of the power exchange with the external network from the reference profile supplied by the REC manager is minimized. Finally, concluding remarks are given in Section 4, where future developments of the work are also envisaged.

2. Materials and Methods

This section provides a detailed insight into the mathematical model of the EMS. The EMS has been modelled as a MIQCQP problem, with continuous and binary decision variables and with linear and quadratic constraints and a quadratic objective function. For each technology installed in the MG, the operational constraints are presented and commented upon.

2.1. Methodology Overview

The flowchart representing the proposed methodology is presented in Figure 3. The manager of the REC is in charge of defining the optimal power exchange profile of the REC members in order to maximize the shared energy within the REC that is, in turn, incentivized. The REC manager provides reference signals to the members, representing a reference active power to possibly be exchanged with the local distribution network. Each member is supposed to be equipped with an EMS (with the exception of pure consumer members) that must define the optimal scheduling of generation units and storage systems locally installed, along with the EV charging infrastructure. Figure 3 highlights the main inputs and outputs of the EMS of the i-th member, that is supposed to be equipped with PV, BESS, and EV charging stations, with the possibility of exploiting V2G. The active power that is actually exchanged by each member with the local distribution network is the one that is accounted for in the assessment of the shared energy.

2.2. Photovoltaic System

The decision variables associated with the PV systems are both continuous and binary ones. Regarding active power, decision variables are continuous and are, namely, the actual active power production, P b , t P V , and the curtailed active power, P b , t P V , c u r t , with b = 1 B and t = 1 T , where B is the number of busses of the MG and T is the number of time intervals into which the optimization horizon is divided. Concerning reactive power, the continuous decision variables are the injected inductive reactive power Q b , t P V , i n d and the injected capacitive reactive power Q b , t P V , c a p ; the binary decision variables are y b , t P V , i n d and y b , t P V , c a p , respectively, equal to 1 when the PV inverter at bus b injects inductive or capacitive reactive power at the time interval t and 0 otherwise.
The only input data are the available active power, P b , t P V , a v and the size of the PV inverter, A b P V . The available active power varies depending on the bus since, even if the MG is located in a limited geographical area, PV panels may have different orientations, tilt angles, and shading.
The following continuous decision variables are considered to be positive:
P b , t P V 0           b = 1 B , t = 1 T
P t P V , c u r t 0           b = 1 B , t = 1 T
Q t P V , i n d 0           b = 1 B , t = 1 T
Q t P V , c a p 0           b = 1 B , t = 1 T
The active power production is given by the difference between the available power and the curtailed quantity, as follows:
P b , t P V = P b , t P V , a v P b , t P V , c u r t           b = 1 B , t = 1 T
The curtailed power can be, at most, equal to the available one, as follows:
P b . t P V , c u r t P b , t P V , a v           b = 1 B , t = 1 T
The relation between active and reactive power, also known as the capability curve of the inverter, is given by Italian technical standards (see CEI 0-16). The capability curve for PV inverters with a rated power lower than 400 [kW] is the one depicted in Figure 4. Since it is formed by a circular sector, it can be described with linear and quadratic constraints, as follows:
Q b , t P V , i n d 0.436 · A b P V · y b , t P V , i n d           b = 1 B , t = 1 T
Q b , t P V , c a p 0.436 · A b P V · y b , t P V , c a p           b = 1 B , t = 1 T
y b , t P V , i n d + y b , t P V , c a p 1           b = 1 B , t = 1 T
( P b , t P V ) 2 + ( Q b , t P V , i n d ) 2 ( A b P V ) 2           b = 1 B , t = 1 T
( P b , t P V ) 2 + ( Q b , t P V , c a p ) 2 ( A b P V ) 2           b = 1 B , t = 1 T

2.3. Battery Energy Storage System

The decision variables associated with the BESS are both continuous and binary. Regarding active power, the continuous decision variables are the charging and discharging power, P b , t B E S S , c h and P b , t B E S S , d c h , while the binary decision variables are x b , t B E S S , c h and x b , t B E S S , d c h , respectively, equal to 1 when the BESS at bus b is in charging or discharging mode at the time interval t and 0 otherwise. Another continuous decision variable is the State of Charge (SOC) of the BESS, S O C b , t B E S S . Concerning reactive power, the continuous decision variables are the injected inductive reactive power Q b , t B E S S , i n d and the injected capacitive reactive power Q b , t B E S S , c a p ; the binary decision variables are y b , t B E S S , i n d and y b , t B E S S , c a p , respectively, equal to 1 when the BESS inverter at bus b injects inductive or capacitive reactive power at the time interval t and 0 otherwise.
The input data were as follows: the BESS inverter rating, A b B E S S ; the maximum charging and discharging power, P b B E S S , c h , m a x and P b B E S S , d c h , m a x ; the charging and discharging efficiencies η b B E S S , c h and η b B E S S , d c h ; the self-discharge rate, γ b B E S S ; the capacity of the BESS, C b B E S S ; and the minimum, maximum, and initial SOC, respectively, indicated by S O C b B E S S , m i n , S O C b B E S S , m a x and S O C b B E S S , i n i .
The following continuous decision variables are considered positive:
P b , t B E S S , c h 0           b = 1 B , t = 1 T
P b , t B E S S , d c h 0           b = 1 B , t = 1 T
Q b , t B E S S , i n d 0           b = 1 B , t = 1 T
Q b , t B E S S , c a p 0           b = 1 B , t = 1 T
S O C b , t B E S S 0           b = 1 B , t = 1 T
The charging and discharging powers must be limited by the maximum values provided by the manufacturer, as follows:
P b , t B E S S , c h P b , t B E S S , c h , m a x · x b , t B E S S , c h           b = 1 B , t = 1 T
P b , t B E S S , d c h P b , t B E S S , d c h , m a x · x b , t B E S S , d c h           b = 1 B , t = 1 T
The following simultaneous charging and discharging of the BESS are forbidden:
x b , t B E S S , c h + x b , t B E S S , d c h 1           b = 1 B , t = 1 T
Regarding reactive power, Standard CEI 0-16 states that the capability curve of BESS inverters is circular, or may possibly be horizontally sectioned if the maximum charging and discharging powers are lower than the inverter rating or if they are considered to be dependent on SOC, cases that are not taken into account here. The standard capability curve is shown in Figure 5. The relevant constraints are as follows:
Q b , t B E S S , i n d A b B E S S · y b , t B E S S , i n d           b = 1 B , t = 1 T
Q b , t B E S S , c a p A b B E S S · y b , t B E S S , c a p           b = 1 B , t = 1 T
y b , t B E S S , i n d + y b , t B E S S , c a p 1           b = 1 B , t = 1 T
( P b , t B E S S , c h ) 2 + ( Q b , t B E S S , i n d ) 2 ( A b B E S S ) 2           b = 1 B , t = 1 T
( P b , t B E S S , c h ) 2 + ( Q b , t B E S S , c a p ) 2 ( A b B E S S ) 2           b = 1 B , t = 1 T
( P b , t B E S S , d c h ) 2 + ( Q b , t B E S S , i n d ) 2 ( A b B E S S ) 2           b = 1 B , t = 1 T
( P b , t B E S S , d c h ) 2 + ( Q b , t B E S S , c a p ) 2 ( A b B E S S ) 2           b = 1 B , t = 1 T
As the decision variables are positive, four sets of constraints (from (23) to (26)) are needed to model the capability curve of the inverter. Constraints (20) to (22) forbid the simultaneous injection of inductive and capacitive reactive power by the BESS inverter.
The energy content of the BESS follows the energy balance and must be limited between a minimum and a maximum value defined by the manufacturer. In addition, due to the structure of the discretized energy balance, the initial SOC must be set in accordance with the input data in the following formulas:
S O C b , t + 1 B E S S = ( 1 γ b B E S S ) · S O C b , t B E S S + Δ C b B E S S · ( η b B E S S , c h · P b , t B E S S , c h P b , t B E S S , d c h η b B E S S , d c h ) b = 1 B , t = 1 T 1
S O C b B E S S , m i n S O C b , t B E S S S O C b B E S S , m a x           b = 1 B , t = 1 T
S O C b , 1 B E S S = S O C b B E S S , i n i           b = 1 B
where Δ is the duration of the time interval.

2.4. Electric Vehicles

The decision variables associated with the EVs, and the related charging stations are both continuous and binary. It is supposed that each EV has one charging station available; therefore, no waiting time is modelled. Regarding active power exchange, it is assumed that some EVs are enabled to V2G and are therefore able to exchange a bidirectional power flow with the MG. The continuous decision variables are the charging and discharging power, P b , v , t E V , c h and P b , v , t E V , d c h , while the binary decision variables are x b , v , t E V , c h and x b , v , t E V , d c h , with v = 1 N b E V , where N b E V is the number of EVs that can be connected at the bus b. x b , v , t E V , c h and x b , v , t E V , d c h are, respectively, equal to 1 when the v-th EV is charging or discharging at bus b at the time interval t and 0 otherwise. Another continuous decision variable is the SOC of the EVs, S O C b , v , t E V . Concerning reactive power, the continuous decision variables are the injected inductive reactive power Q b , v , t E V , i n d and the injected capacitive reactive power Q b , v , t E V , c a p ; the binary decision variables are y b , v , t E V , i n d and y b , v , t E V , c a p , respectively, equal to 1 when the charger of the v-th EV installed at bus b injects inductive or capacitive reactive power at the time interval t and 0 otherwise.
Input data involving EV chargers is the rating of the EV charger converters, A b , v E V , the maximum charging and discharging power, P b , v E V , c h , m a x and P b , v E V , d c h , m a x and the charging and discharging efficiencies, η b , v E V , c h and η b , v E V , d c h . Regarding the availability of the EVs at the charging stations, input data are the sets of arrival, departure and stopover times for each EV, respectively, τ b , v a r r , τ b , v d e p and τ b , v s t o p , and a binary coefficient ε b , v , t E V , equal to 1 if the v-th EV is present at bus b at time t or 0 otherwise (i.e., equal to one for each t belonging to τ s t o p ). Regarding the energy content of the EV battery, the input data are the capacity of the battery C b , v E V , the minimum SOC of the EV, S O C b , v E V , m i n , the expected minimum and maximum SOC at arrival time, S O C b , v E V , a r r , m i n and S O C b , v E V , a r r , m a x , and the expected minimum SOC at the departure time, S O C b , v E V , d e p , m i n .
The following continuous decision variables are considered positive:
P b , v , t E V , c h 0           b = 1 B , v = 1 N b E V , t = 1 T
P b , v , t E V , d c h 0           b = 1 B , v = 1 N b E V , t = 1 T
Q b , v , t E V , i n d 0           b = 1 B , v = 1 N b E V , t = 1 T
Q b , v , t E V , c a p 0           b = 1 B , v = 1 N b E V , t = 1 T
S O C b , v , t E V 0           b = 1 B , v = 1 N b E V , t = 1 T
Charging and discharging powers must be limited by the maximum values, which are evaluated as the minimum between the maximum power that can be exchanged by the charging equipment and the one that can be exchanged by the EV, as follows:
P b , v , t E V , c h P b , v E V , c h , m a x · x b , v , t E V , c h           b = 1 B , v = 1 N b E V , t = 1 T
P b , v , t E V , d c h P b , v E V , d c h , m a x · x b , v , t E V , d c h           b = 1 B , v = 1 N b E V , t = 1 T
Charging and discharging of the EVs are possible only when the EV is present at the facility, as follows; in any case, the simultaneous charging and discharging of the EVs is forbidden:
x b , v , t B E S S , c h + x b , v , t B E S S , d c h ε b , v , t E V           b = 1 B , v = 1 N b E V , t = 1 T
The Italian standards do not provide specific references for the exchange of reactive power by dedicated inverters for EV chargers; therefore, a circular capability curve is applied as follows:
Q b , v , t E V , i n d A b , v E V · y b , v , t E V , i n d           b = 1 B , v = 1 N b E V , t = 1 T
Q b , v , t E V , c a p A b , v E V · y b , v , t E V , c a p           b = 1 B , v = 1 N b E V , t = 1 T
y b , v , t E V , i n d + y b , v , t E V , c a p ε b , v , t E V           b = 1 B , v = 1 N b E V , t = 1 T
( P b , v , t E V , c h ) 2 + ( Q b , v , t E V , i n d ) 2 ( A b , v E V ) 2           b = 1 B , v = 1 N b E V , t = 1 T
( P b , v , t E V , c h ) 2 + ( Q b , v , t E V , c a p ) 2 ( A b , v E V ) 2           b = 1 B , v = 1 N b E V , t = 1 T
( P b , v , t E V , d c h ) 2 + ( Q b , v , t E V , i n d ) 2 ( A b , v E V ) 2           b = 1 B , v = 1 N b E V , t = 1 T
( P b , v , t E V , d c h ) 2 + ( Q b , v , t E V , c a p ) 2 ( A b , v E V ) 2           b = 1 B , v = 1 N b E V , t = 1 T
EV chargers can exchange reactive power only when the EVs are present at the facility, as stated by (40). The same constraints state that inductive and capacitive reactive power cannot be injected simultaneously.
Since the capability curve is fully circular, four quadratic sets of constraints are needed to model it (from (41) to (44)).
Finally, the SOC of the EV batteries is described by the energy balance. The energy balance is defined for the stopover time intervals, when the v-th EV is present at the facility, as follows:
S O C b , v , t + 1 E V = S O C b , v , t E V + Δ C b , v E V · ( η b , v E V , c h · P b , v , t E V , c h P b , v , t E V , d c h η b , v E V , d c h ) b = 1 B , v = 1 N b E V , t τ b , v s t o p
The EV SOC must always lie between a minimum and a maximum value, provided by the manufacturer, as follows:
S O C b , v E V , m i n S O C b , v , t E V S O C b , v E V , m a x           b = 1 B , v = 1 N b E V , t = 1 T
At the arrival times, the EV SOC should be within an expected range, as follows:
S O C b , v E V , a r r , m i n S O C b , v , t E V S O C b , v E V , a r r , m a x           b = 1 B , v = 1 N b E V , t τ b , v a r r
At the departure times, the EVs must have a proper SOC in order to be able to perform their daily activities, as follows:
S O C b , v E V , d e p , m i n S O C b , v , t E V           b = 1 B , v = 1 N b E V , t τ b , v d e p

2.5. Connection to the Distribution Network

The MG is supposed to be connected to the local MV distribution network by means of a step-down transformer. The decision variables are both continuous and binary. Regarding active power, the continuous decision variables are the active power withdrawn from the distribution network and the active power injected into the distribution network, respectively, P t g r i d , b and P t g r i d , s , while the binary decision variables are x t g r i d , b and x t g r i d , s , respectively, equal to 1 when active power is sold to the external network or when it is withdrawn from the distribution network. Regarding reactive power, the continuous decision variables are the inductive and capacitive reactive power provided by the distribution network to the MG, Q t g r i d ,   i n d and Q t g r i d , c a p , while the binary decision variables are y t g r i d , i n d and y t g r i d , c a p , respectively, equal to 1 if inductive or capacitive reactive power is provided by the external network and 0 otherwise.
The data input is the size of the connection transformer, A g r i d .
The following continuous decision variables are considered positive:
P t g r i d , b 0           t = 1 T
P t g r i d , b 0           t = 1 T
Q t g r i d ,   i n d 0           t = 1 T
Q t g r i d , c a p 0           t = 1 T
The active power exchange is limited by the rating of the transformer, as follows:
P t g r i d , b A g r i d · x t g r i d , b           t = 1 T
P t g r i d , s A g r i d · x t g r i d , s           t = 1 T
The active power exchange cannot be bidirectional in the same time interval, as follows:
x t g r i d , b + x t g r i d , s 1           t = 1 T
A circular capability curve is considered for the transformer, modelled through four sets of quadratic constraints (the simultaneous withdrawal of inductive and capacitive reactive power from the external network is, again, forbidden), as follows:
Q t g r i d ,   i n d A g r i d · y t g r i d , i n d           t = 1 T
Q t g r i d ,   c a p A g r i d · y t g r i d , c a p           t = 1 T
y t g r i d , i n d + y t g r i d , c a p 1           t = 1 T
( P t g r i d , b ) 2 + ( Q t g r i d ,   i n d ) 2 ( A g r i d ) 2           t = 1 T
( P t g r i d , b ) 2 + ( Q t g r i d ,   c a p ) 2 ( A g r i d ) 2           t = 1 T
( P t g r i d , s ) 2 + ( Q t g r i d ,   i n d ) 2 ( A g r i d ) 2           t = 1 T
( P t g r i d , s ) 2 + ( Q t g r i d ,   c a p ) 2 ( A g r i d ) 2           t = 1 T

2.6. Optimal Power Flow

In order to model the active and reactive power flows over the MG internal lines and in order to evaluate the effect of those flows on bus voltages, OPF equations are implemented within the model of the EMS. Since classical OPF equations are inherently non-linear, a linearization like the one proposed in [34] has been applied. The formulas are as follows:
p b , k , t = μ b , k · ( v b , t v k , t ) + ξ b , k · ( δ b , t δ k , t )           b , k = 1 B , t = 1 T
q b , k , t = ξ b , k · ( v b , t v k , t ) μ b , k · ( δ b , t δ k , t )           b , k = 1 B , t = 1 T
where p b , k , t and q b , k , t are the per-unit active and reactive power flows over the link between busses b and k, v b , t and v k , t are the per-unit voltage modules at the sending and receiving end of the links, and δ b , t and δ k , t are the voltage phases at the sending and receiving end of the links. p b , k , t , q b , k , t , v b , t and δ b , t are the continuous decision variables.
μ b , k and ξ b , k are obtained as follows:
μ b , k = r b , k r b , k 2 + x b , k 2           b , k = 1 B
ξ b , k = x b , k r b , k 2 + x b , k 2           b , k = 1 B
where r b , k   and x b , k   are the per-unit resistance and reactance of the link between busses b and k.
Voltage modules must be limited between a minimum and a maximum value, as follows:
v min v b , t v max           b = 1 B ,   t = 1 T
Node 1 is considered a slack node; therefore, it is set as the following phase reference node:
v 1 , t = 1           t = 1 T
δ 1 , t = 0           t = 1 T

2.7. Nodal Power Balances

Active and reactive power balances must be satisfied at each time instant at each node.
For node 1, the active and reactive power balances are as follows:
P t g r i d , b P t g r i d , s = p 1 , 2 , t · S b a s e           t = 1 T
Q t g r i d , i n d Q t g r i d , c a p = q 1 , 2 , t · S b a s e           t = 1 T
where S b a s e is the system power base.
For the other nodes, the active and reactive power balances are as follows:
k = 1 B p b , k , t · S b a s e + D b , t e l , p + v = 1 N b E V P b , v , t E V , c h + P b , t B E S S , c h = P b , t P V + v = 1 N b E V P b , v , t E V , d c h + P b , t B E S S , d c h b = 1 B , t = 1 T
k = 1 B q b , k , t · S b a s e + D b , t e l , r + v = 1 N b E V Q b , v , t E V , c a p + Q b , t B E S S , c a p + Q b , t P V , c a p = Q b , t P V , i n d + v = 1 N b E V Q b , v , t E V , i n d + Q b , t B E S S , i n d b = 1 B , t = 1 T
where D b , t e l , p and D b , t e l , r are the active and the reactive power demands of the bus b at time t.

2.8. Renewable Energy Community

Even though the Italian legislation about RECs has still not provided clear rules about the role of BESS and V2G-enabled EVs within the incentive calculation mechanism on shared energy, the authors decided to include the following additional constraints:
x t g r i d , s + x b , t B E S S , d c h 1           b = 1 B , t = 1 T
x t g r i d , s + x b , v , t E V , d c h   1         b = 1 B , v = 1 N b E V , t = 1 T
These constraints prevent the discharging of BESS and EVs when active power is injected into the external network, since up to this moment it is not clear whether the discharged active energy would be taken into account or not for the calculation of the REC shared energy, as the stored energy could have been previously withdrawn from the external grid, rather than being locally generated by PV systems.

2.9. Objective Function

The objective function is made of two terms. The first term is related to the minimization of the net operating costs of the MG, which consider the purchase cost of active energy from the external network, the sale revenues derived from the sale of active energy to the external network, and the penalties on reactive power exchange with the external network. During the central hours of the day, the Distribution System Operator (DSO) sets a penalty for the withdrawal of inductive reactive power from the network, while during off-peak hours, it sets a penalty for the injection of inductive reactive power into the network. This is done in order to avoid, respectively, undervoltages during peak hours and overvoltages during off-peak hours. The second term aims to represent the behavior of the industrial MG within an REC. These are the input signals about the active power that each user must exchange with the external network and, through it, with the other members of the REC provided by the REC manager. These signals should be able to maximize the shared energy in the REC and, therefore, the relevant incentives, in the following formulas:
O b j = α · ( Δ · t = 1 T ( p t e l · P t g r i d , b r t e l · P t g r i d , s + σ t i n d · Q t g r i d , i n d + σ t c a p · · Q t g r i d , c a p ) ) + β · ( 1 α ) · · ( Δ · t = 1 T ( ( P t g r i d , b P t g r i d , b * ) 2 + ( P t g r i d , s P t g r i d , s * ) 2 ) )
where p t e l and r t e l are the unitary purchase costs and unitary sale revenues in [€/kWh], P t g r i d , b * and P t g r i d , s * represent reference signals of purchased and sold active power that the aggregator of the REC sent to each user; α and β are two weight coefficients; and σ t i n d and σ t c a p represent, respectively, the unitary penalties, in [€/kVArh], related to the withdrawal of inductive reactive power during the central hours of the day and the injection of inductive reactive power (dually, the withdrawal of capacitive reactive power) during off-peak hours.

3. Results

This section first shows and discusses the main input data of the optimization problem. Then, the optimization results are presented and discussed in detail. A sensitivity analysis performed varying the weight coefficient α of the objective function is also presented.

3.1. Input Data Description

This subsection describes the topology of the MG and the main input data of the optimization problem. The optimization problem has been set up as a MIQCQP problem, with Δ equal to 1 [h] and T equal to 168 [-].

Microgrid Topology

The MG is composed of six busses (B equal to 6 [-]) in ring configuration, as depicted in Figure 6. The lengths of the lines are as follows:
  • L 2 , 3 = 50   m ;
  • L 3 , 4 = 150   m ;
  • L 4 , 5 = 20   m ;
  • L 5 , 6 = 50   m ;
  • L 6 , 2 = 150   m .
As the lines’ lengths are in the order of hundreds of meters, the approximation introduced linearizing the OPF equations is acceptable [34].
The MG should be made of two sets of three-phase cables, leading to equivalent resistances and reactances, r b , k and x b , k , respectively, equal to 0.08 [Ω/km] and 0.04 [Ω/km].
Bus 1 represents the external distribution network at a rated voltage of 15 [kV] (medium voltage).
The connection transformer rated data are as follows:
  • Transformer rating A   g r i d equal to 800 [kVA];
  • Transformation ratio equal to 15/0.4 [kV]/[kV];
  • Short-circuit voltage equal to 6 [%].
Bus 3 is represented by a building hosting the offices of the company. The building is equipped with a PV plant on the roof and two EV chargers; each charger is dedicated to a Renault Zoe at the disposal of the company’s managers.
The PV input data for bus 3 are as follows:
  • Tilt angle equal to 30°;
  • Azimuth angle equal to 0° (facing south);
  • Inverter rating A 3 P V equal to 10 [kVA];
  • Maximum available power equal to 7.75 [kW].
Load input data for bus 3 are as follows:
  • The maximum active power load is equal to 18.74 [kW];
  • The maximum reactive power load: 13.10 [kVAr];
  • Power factor between 0.8 and 0.995 (lagging).
The Renault Zoe input data for bus 3 are as follows:
  • The size of the converter A 3 , v E V is equal to 22 [kVA];
  • The maximum charging power P 3 , v E V , c h , m a x is equal to 22 [kW] in AC;
  • The maximum discharging power P 3 , v E V , d c h , m a x is equal to 0, since the Renault Zoes are not enabled to V2G;
  • The charging efficiency η 3 , v E V , c h is equal to 95 [%];
  • η 3 , v E V , d c h is not applicable since V2G is not enabled;
  • The battery capacity C 3 , v E V is equal to 52 [kWh];
  • The minimum SOC S O C 3 , v E V , m i n is equal to 10 [%];
  • The expected minimum SOC at arrival time S O C 3 , v E V , a r r , m i n is equal to 20 [%];
  • The expected maximum SOC at arrival time S O C 3 , v E V , a r r , m a x is equal to 40 [%];
  • The expected minimum SOC at the departure time S O C 3 , v E V , d e p , m i n is equal to 80 [%].
The Renault Zoe EVs should arrive at the office building at 9 a.m. and leave at 6 p.m. on weekdays. In addition, one of the two Zoes is present at the facility from 9 a.m. to 6 p.m. on Saturday.
At bus 4, a centralized BESS is installed. The relevant rated data are as follows:
  • BESS inverter rating A 4 B E S S equal to 180 [kVA];
  • Maximum charging and discharging power P 4 B E S S , c h , m a x and P 4 B E S S , d c h , m a x equal to 180 [kW];
  • Charging and discharging efficiencies η 4 B E S S , c h and η 4 B E S S , d c h are equal to 95 [%];
  • Self-discharge rate γ 4 B E S S is equal to 0, since on a short-term time horizon the self-discharge contribution can be neglected;
  • Capacity of the BESS C 4 B E S S equal to 250 [kWh];
  • Minimum SOC S O C 4 B E S S , m i n equal to 8 [%];
  • Maximum SOC S O C 4 B E S S , m a x equal to 98 [%];
  • Initial SOC S O C 4 B E S S , i n i equal to 50 [%].
Bus 5 and 6 represent two industrial warehouses. On the rooftops of the buildings, two PV systems are installed. At each building, one charging station for two DAF XB Electric trucks is installed.
Regarding bus 5, the input data for the PV system are as follows:
  • Tilt angle equal to 20°;
  • Azimuth angle equal to −30°;
  • Inverter rating A 5 P V equal to 100 [kVA];
  • Maximum available power equal to 74.25 [kW].
The load input data for bus 5 are as follows:
  • The maximum active power load is equal to 65.66 [kW];
  • The maximum reactive power load: 87.98 [kVAr];
  • Power factor between 0.57 and 0.88 (lagging).
The DAF XB Electric input data for bus 5 are as follows:
  • The size of the inverter A 5 , 1 E V is equal to 150 [kVA];
  • The maximum charging power P 5 , 1 E V , c h , m a x is equal to 150 [kW] in DC;
  • The maximum discharging power P 5 , 1 E V , d c h , m a x is equal to 150 [kW], since the EV truck should be enabled to V2G;
  • The charging and discharging efficiencies η 5 , 1 E V , c h and η 5 , 1 E V , d c h are equal to 95 [%];
  • The battery capacity C 5 , 1 E V is equal to 315 [kWh];
  • The minimum SOC S O C 5 , v E V , m i n is equal to 10 [%];
  • The expected minimum SOC at arrival time S O C 5 , 1 E V , a r r , m i n is equal to 20 [%];
  • The expected maximum SOC at arrival time S O C 5 , 1 E V , a r r , m a x is equal to 40 [%];
  • The expected minimum SOC at the departure time S O C 5 , 1 E V , d e p , m i n is equal to 80 [%].
The truck is not available at the facility, due to freight delivery, from 8 AM to 5 PM during working days and Saturday; on Sunday it is connected at the facility for the whole day.
Regarding bus 6, the input data for the PV system are as follows:
  • Tilt angle equal to 35°;
  • Azimuth angle equal to −10°;
  • Inverter rating A 6 P V equal to 100 [kVA];
  • Maximum available power equal to 75.58 [kW].
The load input data for bus 6 are as follows:
  • The maximum active power load is equal to 76.12 [kW];
  • The maximum reactive power load: 88.53 [kVAr];
  • Power factor between 0.57 and 0.88 (lagging).
The rated data for the DAF XB Electric connecting at bus 6 are the same as the one that connects at bus 5. The only difference is the availability of the truck at the facility. It was assumed that this truck is used for long-distance freight delivery. Its availability at the facility is described in Figure 7, where 1 means that the truck is connected to the charging station and 0 means that it is not present at the facility.
Regarding economic quantities, the unitary purchase costs p t e l have been considered equal to the zonal market clearing price, increased to 0.12 [€/kWh], while the unitary revenues r t e l have been considered equal to the zonal market clearing price. σ t i n d is set equal to 0.00408 [€/kVArh] during F1 and F2 Italian time bands, while σ t c a p is set equal to 0.00408 [€/kVArh] during the F3 Italian time band. Time bands are defined by the Italian energy authority, ARERA, as follows: F1 is valid on weekdays from 8 a.m. to 7 p.m.; F2 is valid on weekdays from 7 a.m. to 8 a.m. and from 7 p.m. to 11 p.m. and on Saturdays from 7 a.m. to 11 p.m.; and F3 is valid in the remaining time periods and on Sundays and Holidays.
The reference signals provided by the central aggregator of the REC to the industrial user represent the active power that the industrial MG should exchange with the local distribution network in order to maximize the energy that is shared within the REC. The REC is assumed to be composed both of passive users, i.e., users characterized only by active power absorption, like residential users without PV systems, and prosumers, i.e., houses with PV and BESS systems. The reference signals have been developed assuming that the REC aggregator forecasts the load and PV production profiles of the industrial MG without perfectly forecasting it; thus, the forecasting uncertainty is taken into account. The reference profiles for purchased and sold active power are reported, respectively, in Figure 8 and Figure 9.

3.2. Optimization Results

This subsection discusses the results of the optimization model. The EMS was implemented in a Matlab R2022B/Yalmip R20230622 [48] environment and solved with a Gurobi solver [49].
Two main scenarios are analyzed: a scenario in which the EMS must minimize only the economic part of the objective function (i.e., α = 1 ) and a scenario in which the EMS must only minimize the difference between reference signals provided by the REC aggregator and the actual power exchanged with the network (i.e., α = 0 ). Then, intermediate combinations of the objective function weight coefficients are considered in order to build the Pareto front and to evaluate their impact on the facility operations.
Coefficient β has been chosen such that the two terms of the objective function are in the same order of magnitude.

3.2.1. Minimization of the Operating Costs of the Industrial MG

In this scenario, only the operating costs of the facility are minimized, disregarding the maximization of the shared energy within the REC; substantially, the industrial company behaves as though it is not a member of an REC. The resolution time is around 30 s.
The optimal active and reactive power dispatchments of the units are shown in Figure 10 and Figure 11.
In Figure 10, produced active power is shown on the positive y-axis, while absorbed active power is shown on the negative y-axis. PV active power production follows the trend of the available power. On the fifth day, PV production decreases suddenly in the central hours of the day, due to some clouds. On working days, active power is purchased from the network almost for the whole day, due to the high active power demand of the industry. Instead, on Saturday and Sunday, the coupling of the lower active power demand and V2G allows for the purchase of power from the external network to be avoided. The BESS and the EVs are charged during low-price time intervals (e.g., off-peak hours), thereby reducing the overall operating costs. Even during weekdays, in some intervals, V2G and BESS discharging are able to satisfy the industrial demand, either coupled with PV production during peak hours or on their own during off-peak hours. In addition, the EVs can be charged in order to avoid PV curtailment. On Saturday and Sunday, due to the low demand and high PV production, power is injected into the external network.
In Figure 11, injected inductive reactive power is shown on the positive y-axis, while injected capacitive reactive power (or, dually, absorbed inductive reactive power) is shown on the negative y-axis. As can be seen on the graph, the inductive reactive power demand of the loads is always satisfied by the PV and BESS inverters in order to avoid withdrawing it from the external network during the day and allowing for it to be injected into the external network during the night, preventing, respectively, undervoltages and overvoltages in the distribution network.
Figure 12 shows the trend of the voltage amplitude at each bus. Node 1 is the slack bus; therefore, the voltage module is constant at 1 [p.u.]. At the other busses, some overvoltage trends can be observed during the central hours of the day, when PV systems are injecting active and inductive reactive power at the nodes. Undervoltage peaks are evident in correspondence with the charging of the EV trucks, due to the high active power absorbed. Finally, on Saturday and Sunday, overvoltage trends are more evident due to the fact that the active power demand is low; therefore, the impact of PV active power production on the voltage is more significant.
Figure 13 and Figure 14 show the trend of the SOC of the four EVs owned by the company. The two Renault Zoes are not V2G-enabled; they are only charged during the day, applying a smart-charging approach aimed at minimizing operational costs. At the departure time, the two EVs have an SOC equal to 0.8 in order to ensure that the user is able to return home. DAF electric trucks are V2G-enabled; therefore, they can establish a bidirectional power flow with the industrial MG. “DAF XB Electric 1” is the truck connected to Bus 5, while “DAF XB Electric 2” is the truck connected to Bus 6. Both are discharged by the EMS, especially during off-peak periods, namely, during the night and during non-working days, in order to limit or avoid the purchase of power from the external network. In addition, they are used as mobile storage systems, being charged during low electricity price periods.

3.2.2. Maximization of the Shared Energy

In this scenario, only the difference between reference signals of purchased and sold energy and actual power exchanges is minimized. This is the equivalent of maximizing the shared energy within the REC. The resolution time is around 35 s.
The optimal active and reactive power dispatchments of the units are shown in Figure 15 and Figure 16.
In Figure 15, the produced active power is shown on the positive y-axis, while the absorbed active power is shown on the negative y-axis. As with the first scenario, PV active power production follows the trend of the available power. During working days, active power is purchased from the network almost for the whole day, due to the high active power demand of the industry; the only exceptions are the central hours, when the PV production reaches its peak. The same applies on Saturday and Sunday. In this scenario, V2G is not exploited at the weekend in order to adhere, as much as possible, to the reference profile of purchased active power provided by the REC aggregator. In addition, the charging of the EVs is distributed over the entire time horizon, thereby avoiding large peaks of power absorption that would limit costs but that, on the other hand, would decrease the adherence to the reference profile and, therefore, the shared energy. Active power is sold to the local distribution network basically only during the weekend, with some exceptions during the week.
In Figure 16, injected inductive reactive power is shown on the positive y-axis, while injected capacitive reactive power (or, dually, absorbed inductive reactive power) is shown on the negative y-axis. As it is evident from the graph, the inductive reactive power demand of the loads is always satisfied by the PV inverters and by the external network: in fact, in this scenario that aims only at maximizing the shared energy, the reactive power exchange with the network is not penalized.
Figure 17 shows the trend of voltage amplitude at each bus. Node 1 is the slack bus; therefore, the voltage module is constant at 1 [p.u.]. At the other busses, during weekdays, undervoltage behavior is evident, due to the local active power absorption and reactive power exchanges. Finally, on Saturday and Sunday, overvoltage trends are more evident due to the fact that the active power demand is low; therefore, the impact of PV active power production on the voltage is more significant. Globally, voltage values are closer to 1 [p.u.] than in the previous scenario.
Figure 18 and Figure 19 present the trend of the SOC of the four EVs owned by the company. The two Renault Zoes are not V2B-enabled; therefore, their profiles are almost identical to those of the first scenario, except that the charging power has been set to minimize the deviation of the actual power exchange from the reference desired one. For the same reason, DAF electric trucks, despite being V2G-enabled, are mainly used in smart-charging mode during weekdays to allow for the maximization of the shared energy. Instead, V2G is used at the weekend in order to follow the desired profile of active power exchange with the network.

3.2.3. Discussion of Results

The optimal values of the two terms of the objective function are reported in Table 1, varying the weight coefficient α. “Obj1” represents the value of the net operating costs of the industrial MG, while “Obj2” represents the value of the overall deviation of the actual active power exchange with the distribution network from the reference one provided by the REC aggregator. The values presented in the table are expressed in percentages with respect to their maximum values, namely 915.33 [€] for Obj1 (obtained when the EMS minimizes only Obj2) and 696,608 [kWh]2 for Obj2 (obtained when the EMS minimizes only Obj1). The relevant Pareto front is presented in Figure 20. It is evident that the minimization of only the net operating costs leads to a huge increase in the deviation from the REC aggregator’s active power exchange signal, thereby penalizing the energy that is shared within the REC.
In order to assess the operational efficiency of the MG under multi-objective operation, the Self-Sufficiency Rate (SSR) of the MG is calculated as follows:
S S R   [ % ] = E P V + E E V , d c h + E B E S S , d c h E g r i d , s D e l + E E V , c h + E B E S S , c h · 100
The SSR describes the share of the total energy demand of the facility, related to the satisfaction of the buildings’ loads and to the charging of the EVs and the BESS, covered by local resources, namely, the PV plants, the BESS and the EV trucks in V2G mode.
The SSR values are reported in Table 1 for the different values of α . The higher value of the SSR appears as α = 1 , i.e., the scenario in which the net operating costs are minimized. The reason for this is that, in order to limit the withdrawal of active power from the network, the EMS tends to store PV surplus production instead of selling it at a price lower than the purchase price. On the other hand, when α = 0 , i.e., the scenario in which the REC member tries to follow the reference signal of the REC manager, the PV surplus production is sold to cover the requests for power from the other REC members.
To conclude, a further analysis was carried out considering the operation of the MG with the BESS out of service. In this configuration, the net operating costs increase by an average of 3% for all the scenarios, while the deviation from the reference signal provided by the REC manager increases by around 40% for all the scenarios, apart from the scenario where α = 1 . In that scenario, the net operating costs increase by 3.7%, similar to the increases in the other scenarios, but the deviation from the REC manager’s signal decreases by 38% with respect to the configuration with the BESS in service. The reason for this is the fact that, with the BESS in service, in order to minimize the net operating costs, the EMS opts for the maximization of local self-consumption, thereby storing PV surplus production when possible. This is in contrast with the aim of the REC, which is the maximization of energy sharing among REC members. Instead, in the absence of the BESS, the PV surplus production is sold, allowing for it to be shared within the REC.

4. Conclusions

This paper developed and implemented an EMS for an industrial MG owned by a greentech company in northern Italy. This company is a member of an REC, together with other residential users (either consumers or prosumers). The EMS was designed with a multi-objective function to optimize the daily operation of the MG, focusing on minimizing net operating costs and reactive power exchange with the local distribution network and on maximizing adherence to a desired power exchange profile with the external grid, in order to maximize shared energy incentives within the REC. It was assumed that the REC aggregator, through some forecasting about the power demand and renewable production of its users, is able to define an optimal dispatching profile of the power that each user must exchange with the local distribution network. This profile, given as an input to the users’ own EMSs, ensures maximization of the energy that is both shared within the REC and also incentivized.
The industrial MG, modelled as a multi-busbar network, includes PV modules, a centralized BESS, four EV charging stations for two Renault Zoe electric cars, two DAF XB Electric trucks, and the industrial loads related to the warehouses’ activities. The inverters’ capability curves were modelled in accordance with Italian technical standards. Reactive power flows and voltage phenomena were included in the study and were modelled through OPF equations.
The results of the EMS were analyzed for the following two scenarios where: (1) the minimization of operating costs was considered; and (2) where only the minimization of the deviation of the actual active power exchange with the local distribution network from the reference values provided by the REC aggregator was considered. Finally, the Pareto front was built to assess the impact of the objective function’s weight coefficient on the operation of the facility. Regarding the MG operating costs, an increase of 18% was evident when considering the maximization of the energy shared within the REC, if compared to the scenario in which only the operating costs are minimized. On the other hand, in the former scenario, the weekly cumulative quadratic difference between the reference active power exchange and the actual power exchange was 67 times larger than the scenario in which the maximum adherence was granted. Nevertheless, for the industrial user it may not be convenient to follow the reference signal provided by the REC aggregator, given that the shared energy is remunerated at a low tariff and the incentive must be shared with other members.
Future developments of the EMS could involve the application of a multi-level optimization, developing separate EMSs for the REC aggregator and the different typologies of REC members in order to analyze the operation of the entire REC.

Author Contributions

Conceptualization, S.B., M.F. and T.R.; methodology, S.B., M.F. and T.R.; software, M.F. and T.R.; data curation, T.R.; writing—original draft preparation, M.F. and T.R.; writing—review and editing, M.C.; supervision, S.B. and F.D.; project administration, M.C.; funding acquisition, S.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by ENEA, grant Project 1.7 “Technologies for the efficient penetration of the electric vector in the final uses” within the “Electrical System Research” Program Agreement 22-24 between ENEA and the Ministry of the Environment (PTR 22-24).

Data Availability Statement

Data cannot be disclosed due to confidentiality reasons.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. The Green Deal Industrial Plan. Available online: https://commission.europa.eu/strategy-and-policy/priorities-2019-2024/european-green-deal/green-deal-industrial-plan_en (accessed on 12 June 2024).
  2. Eurostat Statistics Explained. Available online: https://ec.europa.eu/eurostat/statistics-explained/index.php?title=Main_Page (accessed on 10 June 2024).
  3. IEA. Buildings; IEA: Paris, France, 2022. [Google Scholar]
  4. Renovation Wave: Creating Green Buildings for the Future. Available online: https://europa.eu/!bG4m7p (accessed on 10 June 2024).
  5. Zhang, D.; Evangelisti, S.; Lettieri, P.; Papageorgiou, L.G. Optimal design of CHP-based microgrids: Multiobjective optimisation and life cycle assessment. Energy 2015, 85, 181–193. [Google Scholar] [CrossRef]
  6. Yang, X.; Leng, Z.; Xu, S.; Yang, C.; Yang, L.; Liu, K.; Song, Y.; Zhang, L. Multi-objective optimal scheduling for CCHP microgrids considering peak-load reduction by augmented ε-constraint method. Renew. Energy 2021, 172, 408–423. [Google Scholar] [CrossRef]
  7. Kostelac, M.; Pavić, I.; Zhang, N.; Capuder, T. Uncertainty modelling of an industry facility as a multi-energy demand response provider. Appl. Energy 2022, 307, 118215. [Google Scholar] [CrossRef]
  8. Fleschutz, M.; Bohlayer, M.; Braun, M.; Murphy, M.D. From prosumer to flexumer: Case study on the value of flexibility in decarbonizing the multi-energy system of a manufacturing company. Appl. Energy 2023, 347, 121430. [Google Scholar] [CrossRef]
  9. Lavi, Y.; Apt, J. Using PV inverters for voltage support at night can lower grid costs. Energy Rep. 2022, 8, 6347–6354. [Google Scholar] [CrossRef]
  10. Sharma, M.; Luthra, S.; Joshi, S.; Kumar, A.; Jain, A. Green logistics driven circular practices adoption in industry 4.0 Era: A moderating effect of institution pressure and supply chain flexibility. J. Clean. Prod. 2023, 383, 135284. [Google Scholar] [CrossRef]
  11. Oloruntobi, O.; Mokhtar, K.; Mohd Rozar, N.; Gohari, A.; Asif, S.; Chuah, L.F. Effective technologies and practices for reducing pollution in warehouses—A review. Clean. Eng. Technol. 2023, 13, 100622. [Google Scholar] [CrossRef]
  12. Islam, S.R.; Novoa, C.; Jin, T. Multi-facility aggregate production planning with prosumer microgrid: A two-stage stochastic program. J. Clean. Prod. 2022, 367, 132911. [Google Scholar] [CrossRef]
  13. Li, H.; Jin, T.; Novoa, C. Facility and microgrid location-allocation for integrated supply chain and transactive energy operations. Appl. Math. Model. 2023, 119, 119–136. [Google Scholar] [CrossRef]
  14. Pham, A.; Jin, T.; Novoa, C.; Qin, J. A multi-site production and microgrid planning model for net-zero energy operations. Int. J. Prod. Econ. 2019, 218, 260–274. [Google Scholar] [CrossRef]
  15. Perković, L.; Mikulčić, H.; Duić, N. Multi-objective optimization of a simplified factory model acting as a prosumer on the electricity market. J. Clean. Prod. 2017, 167, 1438–1449. [Google Scholar] [CrossRef]
  16. Polleux, L.; Guerassimoff, G.; Marmorat, J.-P.; Sandoval-Moreno, J.; Schuhler, T. An overview of the challenges of solar power integration in isolated industrial microgrids with reliability constraints. Renew. Sustain. Energy Rev. 2022, 155, 111955. [Google Scholar] [CrossRef]
  17. Blake, S.T.; O’Sullivan, D.T.J. Optimization of Distributed Energy Resources in an Industrial Microgrid. Procedia CIRP 2018, 67, 104–109. [Google Scholar] [CrossRef]
  18. Viana, M.S.; Ramos, D.S.; Manassero Junior, G.; Udaeta, M.E.M. Analysis of the Implementation of Virtual Power Plants and Their Impacts on Electrical Systems. Energies 2023, 16, 7682. [Google Scholar] [CrossRef]
  19. López, I.; Goitia-Zabaleta, N.; Milo, A.; Gómez-Cornejo, J.; Aranzabal, I.; Gaztañaga, H.; Fernandez, E. European energy communities: Characteristics, trends, business models and legal framework. Renew. Sustain. Energy Rev. 2024, 197, 114403. [Google Scholar] [CrossRef]
  20. Marchetti, B.; Vitali, M.; Biancini, G. Renewable Energy Proliferation and the New Local Energy Community Paradigm: Analysis of a Case Study in Italy. Energies 2024, 17, 1599. [Google Scholar] [CrossRef]
  21. Ahmed, S.; Ali, A.; D’Angola, A. A Review of Renewable Energy Communities: Concepts, Scope, Progress, Challenges, and Recommendations. Sustainability 2024, 16, 1749. [Google Scholar] [CrossRef]
  22. Battaglia, V.; Vanoli, L.; Zagni, M. Economic benefits of Renewable energy communities in smart districts: A comparative analysis of incentive schemes for NZEBs. Energy Build. 2024, 305, 113911. [Google Scholar] [CrossRef]
  23. Aruta, G.; Ascione, F.; Bianco, N.; Iovane, T.; Mastellone, M.; Maria Mauro, G. Optimizing the energy transition of social housing to renewable nearly zero-energy community: The goal of sustainability. Energy Build. 2023, 282, 112798. [Google Scholar] [CrossRef]
  24. Corsini, A.; Delibra, G.; Pizzuti, I.; Tajalli-Ardekani, E. Challenges of renewable energy communities on small Mediterranean islands: A case study on Ponza island. Renew. Energy 2023, 215, 118986. [Google Scholar] [CrossRef]
  25. Caballero, V.; Briones, A.; Coca-Ortegón, A.; Pérez, A.; Barrios, B.; de la Mano, M. Analysis and simulation of an Urban-Industrial Sustainable Energy Community: A use case in San Juan de Mozarrifar using photovoltaic energy. Energy Rep. 2023, 9, 1589–1605. [Google Scholar] [CrossRef]
  26. Brusco, G.; Menniti, D.; Pinnarelli, A.; Sorrentino, N. Renewable Energy Community with distributed storage optimization to provide energy sharing and additional ancillary services. Sustain. Energy Grids Netw. 2023, 36, 101173. [Google Scholar] [CrossRef]
  27. Cosic, A.; Stadler, M.; Mansoor, M.; Zellinger, M. Mixed-integer linear programming based optimization strategies for renewable energy communities. Energy 2021, 237, 121559. [Google Scholar] [CrossRef]
  28. Mariuzzo, I.; Fioriti, D.; Guerrazzi, E.; Thomopulos, D.; Raugi, M. Multi-objective planning method for renewable energy communities with economic, environmental and social goals. Int. J. Electr. Power Energy Syst. 2023, 153, 109331. [Google Scholar] [CrossRef]
  29. Weckesser, T.; Dominković, D.F.; Blomgren, E.M.V.; Schledorn, A.; Madsen, H. Renewable Energy Communities: Optimal sizing and distribution grid impact of photovoltaics and battery storage. Appl. Energy 2021, 301, 117408. [Google Scholar] [CrossRef]
  30. MOTUS-E Italian Association. Available online: https://www.motus-e.org/ (accessed on 12 June 2024).
  31. Fiaschi, D.; Bandinelli, R.; Conti, S. A case study for energy issues of public buildings and utilities in a small municipality: Investigation of possible improvements and integration with renewables. Appl. Energy 2012, 97, 101–114. [Google Scholar] [CrossRef]
  32. Sadeghian, O.; Oshnoei, A.; Mohammadi-Ivatloo, B.; Vahidinasab, V.; Anvari-Moghaddam, A. A comprehensive review on electric vehicles smart charging: Solutions, strategies, technologies, and challenges. J. Energy Storage 2022, 54, 105241. [Google Scholar] [CrossRef]
  33. Yu, H.; Niu, S.; Shang, Y.; Shao, Z.; Jia, Y.; Jian, L. Electric vehicles integration and vehicle-to-grid operation in active distribution grids: A comprehensive review on power architectures, grid connection standards and typical applications. Renew. Sustain. Energy Rev. 2022, 168, 112812. [Google Scholar] [CrossRef]
  34. Fresia, M.; Bordo, L.; Delfino, F.; Bracco, S. Optimal day-ahead active and reactive power management for microgrids with high penetration of renewables. Energy Convers. Manag. X 2024, 23, 100598. [Google Scholar] [CrossRef]
  35. Piazza, G.; Delfino, F.; Bergero, S.; Di Somma, M.; Graditi, G.; Bracco, S. Economic and environmental optimal design of a multi-vector energy hub feeding a Local Energy Community. Appl. Energy 2023, 347, 121259. [Google Scholar] [CrossRef]
  36. Di Somma, M.; Buonanno, A.; Caliano, M.; Graditi, G.; Piazza, G.; Bracco, S.; Delfino, F. Stochastic Operation Optimization of the Smart Savona Campus as an Integrated Local Energy Community Considering Energy Costs and Carbon Emissions. Energies 2022, 15, 8418. [Google Scholar] [CrossRef]
  37. Piazza, G.; Bracco, S.; Delfino, F.; Somma, M.D.; Graditi, G. Impact of electric mobility on the design of renewable energy collective self-consumers. Sustain. Energy Grids Netw. 2023, 33, 100963. [Google Scholar] [CrossRef]
  38. Piazza, G.; Bracco, S.; Delfino, F.; Siri, S. Optimal design of electric mobility services for a Local Energy Community. Sustain. Energy Grids Netw. 2021, 26, 100440. [Google Scholar] [CrossRef]
  39. Li, Y.; Han, M.; Yang, Z.; Li, G. Coordinating Flexible Demand Response and Renewable Uncertainties for Scheduling of Community Integrated Energy Systems With an Electric Vehicle Charging Station: A Bi-Level Approach. IEEE Trans. Sustain. Energy 2021, 12, 2321–2331. [Google Scholar] [CrossRef]
  40. Li, Y.; Wang, B.; Yang, Z.; Li, J.; Li, G. Optimal Scheduling of Integrated Demand Response-Enabled Community-Integrated Energy Systems in Uncertain Environments. IEEE Trans. Ind. Appl. 2022, 58, 2640–2651. [Google Scholar] [CrossRef]
  41. Bracco, S.; Delfino, F.; Longo, M.; Siri, S. Electric Vehicles and Storage Systems Integrated within a Sustainable Urban District Fed by Solar Energy. J. Adv. Transp. 2019, 2019, 9572746. [Google Scholar] [CrossRef]
  42. Bracco, S.; Cancemi, C.; Causa, F.; Longo, M.; Siri, S. Optimization model for the design of a smart energy infrastructure with electric mobility. IFAC-Pap. 2018, 51, 200–205. [Google Scholar] [CrossRef]
  43. Zhou, Y.; Cao, S.; Hensen, J.L.; Lund, P.D. Energy integration and interaction between buildings and vehicles: A state-of-the-art review. Renew. Sustain. Energy Rev. 2019, 114, 109337. [Google Scholar] [CrossRef]
  44. Zheng, S.; Huang, G.; Lai, A.C. Coordinated energy management for commercial prosumers integrated with distributed stationary storages and EV fleets. Energy Build. 2023, 282, 112773. [Google Scholar] [CrossRef]
  45. He, Z.; Khazaei, J.; Freihaut, J. Optimal integration of Vehicle to Building (V2B) and Building to Vehicle (B2V) technologies for commercial buildings. Sustain. Energy Grids Netw. 2022, 32, 100921. [Google Scholar] [CrossRef]
  46. Bracco, S.; Fresia, M. Energy Management System for the Optimal Operation of a Grid-Connected Building with Renewables and an Electric Delivery Vehicle. In Proceedings of the IEEE EUROCON 2023-20th International Conference on Smart Technologies, Turin, Italy, 6–8 July 2023; pp. 472–477. [Google Scholar]
  47. Fresia, M.; Bracco, S. Electric Vehicle Fleet Management for a Prosumer Building with Renewable Generation. Energies 2023, 16, 7213. [Google Scholar] [CrossRef]
  48. Lofberg, J. YALMIP: A toolbox for modeling and optimization in MATLAB. In Proceedings of the 2004 IEEE International Conference on Robotics and Automation (IEEE Cat. No. 04CH37508), New Orleans, LA, USA, 26 April–1 May 2004; pp. 284–289. [Google Scholar]
  49. Gurobi Optimization. Gurobi Optimizer Reference Manual. Available online: https://www.gurobi.com (accessed on 10 June 2024).
Figure 1. Occurrence of “Green warehouse” and “Sustainable warehouse” terms in the Scopus database.
Figure 1. Occurrence of “Green warehouse” and “Sustainable warehouse” terms in the Scopus database.
Energies 17 03567 g001
Figure 2. Occurrence of EV smart-charging and V2X technology topics in the Scopus database.
Figure 2. Occurrence of EV smart-charging and V2X technology topics in the Scopus database.
Energies 17 03567 g002
Figure 3. Flowchart of the methodology.
Figure 3. Flowchart of the methodology.
Energies 17 03567 g003
Figure 4. PV inverter capability curve.
Figure 4. PV inverter capability curve.
Energies 17 03567 g004
Figure 5. BESS inverter capability curve.
Figure 5. BESS inverter capability curve.
Energies 17 03567 g005
Figure 6. Industrial MG topology.
Figure 6. Industrial MG topology.
Energies 17 03567 g006
Figure 7. DAF XB Electric availability at Bus 6. Green when available, red when not available.
Figure 7. DAF XB Electric availability at Bus 6. Green when available, red when not available.
Energies 17 03567 g007
Figure 8. Reference profile for active power purchased from the external network.
Figure 8. Reference profile for active power purchased from the external network.
Energies 17 03567 g008
Figure 9. Reference profile for active power sold to the external network.
Figure 9. Reference profile for active power sold to the external network.
Energies 17 03567 g009
Figure 10. EMS active power results—minimization of only operating costs.
Figure 10. EMS active power results—minimization of only operating costs.
Energies 17 03567 g010
Figure 11. EMS reactive power results—minimization of only operating costs.
Figure 11. EMS reactive power results—minimization of only operating costs.
Energies 17 03567 g011
Figure 12. EMS voltage magnitude results—minimization of only operating costs.
Figure 12. EMS voltage magnitude results—minimization of only operating costs.
Energies 17 03567 g012
Figure 13. SOC of the two Renault Zoes—minimization of only operating costs.
Figure 13. SOC of the two Renault Zoes—minimization of only operating costs.
Energies 17 03567 g013
Figure 14. SOC of the two DAF XB Electric trucks—minimization of only operating costs.
Figure 14. SOC of the two DAF XB Electric trucks—minimization of only operating costs.
Energies 17 03567 g014
Figure 15. EMS active power results—maximization of shared energy.
Figure 15. EMS active power results—maximization of shared energy.
Energies 17 03567 g015
Figure 16. EMS reactive power results—maximization of shared energy.
Figure 16. EMS reactive power results—maximization of shared energy.
Energies 17 03567 g016
Figure 17. EMS voltage magnitude results—msaximization of shared energy.
Figure 17. EMS voltage magnitude results—msaximization of shared energy.
Energies 17 03567 g017
Figure 18. SOC of the two Renault Zoe—maximization of shared energy.
Figure 18. SOC of the two Renault Zoe—maximization of shared energy.
Energies 17 03567 g018
Figure 19. SOC of the two DAF XB Electric trucks—maximization of shared energy.
Figure 19. SOC of the two DAF XB Electric trucks—maximization of shared energy.
Energies 17 03567 g019
Figure 20. Pareto front.
Figure 20. Pareto front.
Energies 17 03567 g020
Table 1. Optimal results varying the weight coefficient.
Table 1. Optimal results varying the weight coefficient.
ScenarioObj1 [%]Obj2 [%]SSR [%]
α = 0 100 1.48756.43
α = 0.1 99.841.48856.45
α = 0.2 99.651.48956.46
α = 0.3 99.391.48956.51
α = 0.4 99.151.49156.81
α = 0.5 98.871.49456.86
α = 0.6 98.501.50057.30
α = 0.7 97.041.53857.41
α = 0.8 95.031.62157.81
α = 0.9 93.481.72458.06
α = 1 84.8610065.37
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Fresia, M.; Robbiano, T.; Caliano, M.; Delfino, F.; Bracco, S. Optimal Operation of an Industrial Microgrid within a Renewable Energy Community: A Case Study of a Greentech Company. Energies 2024, 17, 3567. https://doi.org/10.3390/en17143567

AMA Style

Fresia M, Robbiano T, Caliano M, Delfino F, Bracco S. Optimal Operation of an Industrial Microgrid within a Renewable Energy Community: A Case Study of a Greentech Company. Energies. 2024; 17(14):3567. https://doi.org/10.3390/en17143567

Chicago/Turabian Style

Fresia, Matteo, Tommaso Robbiano, Martina Caliano, Federico Delfino, and Stefano Bracco. 2024. "Optimal Operation of an Industrial Microgrid within a Renewable Energy Community: A Case Study of a Greentech Company" Energies 17, no. 14: 3567. https://doi.org/10.3390/en17143567

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop