Predictive Control and Benefit Sharing in Multi-Energy Systems

This article proposes the design of a hierarchical control architecture capable of optimally coordinating multi-energy systems (MESs). A MES involves the synergetic operation of subsystems belonging to different energy domains (e.g., thermal, electrical, or gas), enhancing their energy efficiency and economic savings, at the price of significant control challenges. In fact, MESs imply an increased model complexity and the interaction of networked subsystems with largely different dynamics. This motivates the design of a multilayer control architecture where, at the upper level, a model predictive control (MPC) regulator relies on energy models of reduced order to coordinate power exchanges among MES subsystems, while, at the lower layer, decentralized MPC regulators locally control subsystems with different sampling rates, consistently with their local dynamics. On the other hand, the optimal MES regulation may imply additional costs to few subsystems, although the overall operational cost decreases. Thus, benefit sharing algorithms are also proposed, relying on game-theoretical methods, enabling to properly share the achieved economic benefit among subsystems, and guaranteeing that the MES operation is more convenient than the independent regulation for each single subsystem. The designed control strategy is tested on two different MES case studies, considering also the presence of referenced electrical and thermal networks, showing high versatility and enhanced performances.


Predictive Control and Benefit Sharing in Multi-Energy Systems
Alessio La Bella , Member, IEEE, Lorenzo Nigro, Graduate Student Member, IEEE, and Riccardo Scattolini Abstract-This article proposes the design of a hierarchical control architecture capable of optimally coordinating multienergy systems (MESs).A MES involves the synergetic operation of subsystems belonging to different energy domains (e.g., thermal, electrical, or gas), enhancing their energy efficiency and economic savings, at the price of significant control challenges.In fact, MESs imply an increased model complexity and the interaction of networked subsystems with largely different dynamics.This motivates the design of a multilayer control architecture where, at the upper level, a model predictive control (MPC) regulator relies on energy models of reduced order to coordinate power exchanges among MES subsystems, while, at the lower layer, decentralized MPC regulators locally control subsystems with different sampling rates, consistently with their local dynamics.On the other hand, the optimal MES regulation may imply additional costs to few subsystems, although the overall operational cost decreases.Thus, benefit sharing algorithms are also proposed, relying on game-theoretical methods, enabling to properly share the achieved economic benefit among subsystems, and guaranteeing that the MES operation is more convenient than the independent regulation for each single subsystem.The designed control strategy is tested on two different MES case studies, considering also the presence of referenced electrical and thermal networks, showing high versatility and enhanced performances.

I. INTRODUCTION
M ULTI-ENERGY systems (MESs) are today gaining much attention, being recognized as a key solution to achieve a reliable, resilient, and low-carbon energy system [1].The MES concept consists in the synergetic operation of subsystems belonging to different energy domains, which can exchange energy through proper conversion devices (e.g., interconnected electrical grids, district heating systems, and gas networks).The optimal coordination of these additional degrees of flexibility can significantly increase the efficiency and storage capacity of the energy system, which are crucial aspects to support a large-scale diffusion of intermittent renewable energy sources [2].For example, a surplus of renewable electrical energy can be stored in other forms, such as in a thermal or a gas storage tank through a proper conversion process, and used when necessary.The growing interest in MESs has been also pushed by the today commercial availability of different conversion devices, such as cogeneration systems (gas to electricity and heat), heat pumps (electricity to heat), hydrogen electrolysers (electricity to gas), and fuel cells (gas to electricity).
The benefits given by an integrated control of MESs come, however, at a price.In fact, MESs can be large-scale systems characterized by very different dynamic behaviors, such as "fast" electricity systems and "slow" thermal ones, and they can be governed by complex nonlinear equations (as for thermal and gas networks [3], [4]).Moreover, since energy systems of different domains may belong to different owners or operators, privacy issues can be of paramount importance.Finally, it is also worth noting that, despite the coordinated control of MESs can lead to consistent economic saving, this benefit could be a priori not equally shared among the involved energy systems or owners.In fact, for each single-energy system, it must be ensured that the coordinated MES optimization is more convenient than its autonomous regulation.

A. State of the Art
The problem of modeling and control of MESs has recently gained attention, given their inherent complexity and benefits, as discussed in [5].In [6], the energy hub concept is discussed, comprising the interconnection of different MES devices, and showing that applying model predictive control (MPC) strategies can actually lead to enhanced energy savings.A generalized modeling methodology for MESs has been proposed in [7], leading to the formulation of mixed-integer linear models suited for the design of MPC regulators.In [8], a MES is modeled using state-space linear systems, and an economic MPC is designed for its optimal regulation.Other optimization-based control approaches for MESs consider robust methods [9], [10], data-based approaches [11], [12], [13], or the presence of time-varying constraints [14].Nevertheless, all the mentioned contributions concern the design of centralized regulators, which could be an impractical solution for MESs due to computational issues and privacy concerns.
Motivated by this, distributed approaches have been also proposed.The design of a distributed MPC scheme for MESs is described in [15], relying on dual-decomposition techniques.In [16] and [17], distributed optimization algorithms for energy systems described by mixed-integer linear models are discussed.A distributed algorithm relying on the alternating direction method of multipliers (ADMM) is proposed in [18], assuming networked energy systems to be represented by models.Nevertheless, distributed control strategies can be critical in a multi-energy scenario.They often involve iterative optimization processes, which can be slow, and, moreover, are not guaranteed to converge to the optimal solution in the presence of nonconvex and nonlinear problems, typical in MESs.Because of this, a heuristic procedure based on the relaxation of nonlinear constraints is proposed in [19], able to recover at least a feasible solution if the designed distributed algorithm fails in optimizing the MES.Other control approaches for networked energy systems have relied on hierarchical architectures, enabling enhanced scalability and prompt system regulation [20], [21], [22], [23].In this context, a two-layer control architecture for large-scale networks is proposed in [24], exploiting a network partitioning procedure to identify the subsystems to control, but, however, focusing only on the electrical energy domain.
Concerning the distribution of the benefit achieved by the optimal control of interconnected energy subsystems, most literature has relied on game-theoretical approaches.A wellknown technique for this task is the Shapley value method, first proposed in [25], where the profit is distributed based on the effective contribution of each agent to the coalition [26].Profit redistribution via Shapley value is applied to a residential microgrid in [27], however highlighting that it can be intractable from the computational point of view, as each agent needs to know the achievable benefit by any possible coalition of agents [28].An alternative game-theoretical method to redistribute the benefit in cooperative fashion is the Nash bargaining problem, first proposed in [29].This benefit sharing technique is exploited in [30] and [31], considering the interconnection of electrical microgrids.Similar methods are applied in [32] and [33], proposing distributed bargaining algorithms, which lead to real-time trading frameworks among energy subsystems.
Nevertheless, note that most of the mentioned references consider simplified system models (e.g., mixed-integer linear), which are not suited to represent dynamical energy networks, such as thermal and gas ones, which are inherently nonlinear.On the other hand, the presented approaches use the same sampling period for controlling all devices, which can compromise performances in the presence of MESs, as they often comprehend subsystems with largely different dynamics and transients.Finally, note that the mentioned references do not investigate benefit sharing strategies applied to different multienergy domains, which usually belong to distinct owners.

B. Main Contribution
The mentioned challenges motivate the design of a novel hierarchical MPC architecture for MESs.In particular, we consider a MES including interconnected subsystems belonging to distinct energy domains (e.g., an electrical grid, a district heating system, and a gas network) and that can exchange energy through proper conversion devices (e.g., heat pumps, gas boilers, cogenerators, and so on).The proposed control architecture includes, at the high level, an MPC, which, exploiting properly energy models of reduced order, computes the optimal amount of energy to be stored in each MES subsystem and the power exchanges among them.The proposed reduced modeling is independent of the specific energy domain, and it enables to describe the subsystem energy dynamics at reduced complexity with respect to their detailed modeling, facilitating their multienergetic coordination.Power exchanges among subsystems computed at the high level are directly imposed on conversion devices, whereas the computed energy references are sent to a decentralized low-level layer of MPC, each one regulating a subsystem at the minimum cost.Each low-level MPC considers the detailed subsystem model, possibly nonlinear, and is executed with a sampling period consistent with the subsystem dynamical behavior.
Concerning the benefit sharing problem, a cooperative bargaining game is formulated, resulting in an agreement among subsystems on how the profit given by the multienergy interaction should be shared.In particular, the proposed bargaining game is supposed to be periodically executed, evaluating the benefit achieved by the designed control strategy during the past instants and then defining the optimal economic transactions among MES subsystems, so that the obtained benefit is properly distributed.
The main advantages and contributions of the proposed control strategy with respect to the literature are as follows.
1) At the high level, the definition of a novel reducedorder modeling of MES subsystems, leading to convex energetic models whose structure is independent of the specific energy domain, enabling the computationally efficient and prompt coordination of power exchanges among MES subsystems and of their stored energy.2) At the low level, the possibility of controlling each MES subsystem through a local MPC (L-MPC) regulator, considering the detailed system modeling, eventually nonlinear, tracking an optimized overall energy profile at the minimum cost.Moreover, each L-MPC regulator can regulate its subsystem with a specific sampling period, consistently with the dynamical characteristics of the corresponding energy domain.
3) The design of a benefit sharing algorithm, which, relying on bargaining games, enables each subsystem (or owner) to benefit from the MES operation.In particular, the proposed formulation guarantees that, for each subsystem, the multi-energy cooperation is always more convenient with respect to the independent regulation.Moreover, the presented approach enables subsystems to preserve their internal information and models, as they must periodically communicate just their total operational cost.The proposed control strategy is tested on two MES case studies.The first, more illustrative, has the purpose of clearly showing the validity of the proposed modeling and how the designed architecture actually solves the highlighted control challenges.The second case study is based on a more realistic MES scenario, where two subsystems are interconnected and jointly regulated, i.e., the IEEE 37-bus electrical grid [34] and the AROMA district heating network (DHN) [3], showing how the designed control architecture is effective also for the coordination of large-scale MES benchmarks.

C. Outline
This article is organized as follows.Section II introduces the MES modeling, the control objectives, and presents an illustrative example to motivate the addressed control challenges.The proposed two-layer control architecture is described in detail in Section III.Section IV illustrates the formulated bargaining game problem.Section V reports the analysis of the designed control architecture tested on two MES case studies.The final conclusions are gathered in Section VI.
Notation: Let R denote the set of real numbers, R ≥0 the set of positive or null real numbers including zero, and R >0 the set of strictly positive real numbers.Moreover, let N denote the set of natural numbers.Given a matrix A, A ∈ R n,m indicates that it has n rows and m columns, and all its entries are real numbers.Given a vector v, v ∈ R n indicates that it has n rows and one column, and all its entries are real numbers.Given a matrix A and a vector v, their transpose is denoted with A ′ and v ′ , respectively.The vector 1 n is used to define an n-dimensional vector where all entries are 1, i.e., 1 n = [1, . . ., 1] ′ ∈ R n .Similarly, the vector 0 n is used to define an n-dimensional vector where all entries are 0, i.e., 0 n = [0, . . ., 0] ′ ∈ R n .Given two vectors of variables x, y ∈ R n , the inequalities among the two, e.g., x > y, are intended elementwise.For a vector x ∈ R n , the vectors of the corresponding upper and lower bounds are x ∈ R n and x ∈ R n , respectively, with x > x.The operators max(x) and min(x) express the maximum and the minimum element of vector x, respectively.Considering a real variable a ∈ R ≥0 , b = ⌊a⌋ is the largest integer less than or equal to a, i.e., b ∈ N ∪ {0}.Given a sequence of variables a 1 , . . ., a n , and the set of their indexes N = {1, . . ., n}, the vector a = [a 1 , . . ., a n ] ′ is compactly written as a = {a i } ∀i∈N .Finally, given a set N , its cardinality is denoted with card(N ).

II. PROBLEM STATEMENT
Consider a MES S † composed by a set of M subsystems S 1 , . . ., S M , where each subsystem S i belongs to a specific energy domain (electrical, thermal, gas, and so on) and is assumed to be equipped with local generation and load units.Moreover, assume that power can be transmitted among some specific pairs of subsystems to optimize the overall energy behavior of S † .More formally, this structure can be described by a directed graph G(N , E), where the nodes N represent the subsystems, i.e., card(N ) = M, while the edges E ⊆ N × N model their energetic interconnections, i.e., (i, j) ∈ E if power can be transferred from S i to S j .For each S i , i ∈ N , consider the sets of output and input nodes N out i , N in i ⊆ N , defined as Each subsystem S i is modeled as follows: where x i ∈ R n i is the vector of local states and u i ∈ R m i is the vector of local inputs, both assumed to be bounded as follows: Moreover, in (1a), d i represents the vector of exogenous signals for S i , defined as follows: where d i ∈ R r i are local disturbances, while w i includes the power exchanges between S i and the other subsystems.
In particular, w i is structured as follows: where w out i (t) ∈ R o i and w in i ∈ R v i include the power transfers sent and received by S i , respectively.Expressing with p i j ∈ R the power produced by S i and transferred to S j ,1 it follows that: where the power transfers in w in i are weighted by the efficiency term η ji ∈ R, considering the conversion losses due to the power transfer from S j to S i , with j ∈ N in i .Now, assume that a set G i of controllable sources and a set L i of nondispatchable loads are present in each S i and denote with p g i ∈ R g i and p l i ∈ R l i their power injections and absorptions, respectively.These are described as the functions of the local variables, i.e., where, in particular, controllable power flows are the functions of local inputs, while the nondispatchable ones depend on local disturbances. 2oreover, introduce the variable e i ∈ R, representing the total energy stored in S i .Then, the following modeling assumption can be stated.
Assumption 1: Given a generic energy subsystem S i described by (1), the total stored energy e i ∈ R can be defined as a function of the local state, i.e., where ξ i : R n i → R.Moreover, its derivative with respect to time can be stated as follows: □ The formulated assumption is straightforward, as local states in an energy system are generally linked to the stored energy (e.g., batteries' state of charges, fluid temperatures, gas pressures, and so on), and the variation of the total energy is dictated by the power balance.As it will be shown, several energy systems of different domains can be represented by ( 1) and ( 2) and satisfy Assumption 1.However, the addressed multi-energy control problem and the proposed solution are firstly presented.

A. Control Problem and Proposed Solution
Considering the presented framework for networked MESs, the main control objectives can be synthesized as follows.
1) Objective 1: Locally control each S i to maintain its states in prescribed operating bounds, meanwhile minimizing local power production costs.2) Objective 2: Coordinate the power exchanges p i j among interconnected subsystems, possibly belonging to distinct energy domains, considering the energy availability, and the requirements of each S i , as well as the minimization of the operational cost of the whole networked MES S † .As discussed in Section I, pure centralized and distributed control approaches are not suited to accomplish these tasks in a multi-energy scenario, given the large-scale dimension of the problem, the different involved time scales, and the modeling complexity.This motivates the design of a novel hierarchical predictive control architecture, depicted in Fig. 1.In this structure, a lower layer made by decentralized L-MPC, one for each S i , runs at "small" sampling periods, assuming that the energy exchange terms p i j have been set.Each L-MPC i computes local inputs u i to satisfy local constraints, minimize internal operational costs, and maintain the total energy stored in S i around a given reference, denoted with e h i .This reference is computed by a coordination scheme at the higher layer, named high-level MPC (HL-MPC), which, exploiting a reduced energy model for each subsystem and according to an economic criterion, computes the optimal energy to be stored, i.e., e h i , and the optimal power exchanges among the interconnected subsystems, denoted with p h i j ∈ R, at "large" sampling periods (for the sake of clarity, all variables regulated by the HL-MPC are reported with the superscript h).As shown in the following, the proposed two-layer control architecture enables, on the one hand, to locally control each energy subsystem at a specific time rate, separately considering the dynamical model complexity of diverse energy domains, and, on the other hand, to optimally coordinate power exchanges among subsystems taking advantage of a reduced modeling approach.Before describing in detail the formulation of each control layer, an illustrative example of a simple MES is hereafter reported to better motivate the proposed system modeling and assumptions, as well as the designed control architecture.

B. Illustrative Example of a MES
Consider the MES S † depicted in Fig. 2, where three subsystems belonging to different energy domains are interconnected: a gas storage tank (subsystem S 1 ), a water thermal system (subsystem S 2 ), and an electrical circuit (subsystem S 3 ).
The gas storage tank has volume V g and contains gas at density ρ g (t).It is subjected to an input mass flow rate q in g (t), which is locally regulated, and to an output mass flow rate q out g (t), acting as a disturbance.Moreover, part of the stored gas energy is transferred to the water thermal system through a gas boiler absorbing the primary power p b g (t).The water thermal system includes a water vessel with mass m w and temperature T w (t), heated by a temperature-regulated thermal plate of area A and temperature T m (t).The water contained in the vessel is also heated by a thermal resistance, fed by the external electrical circuit with an electrical power p r e (t).A mass flow rate q w (t) of water is delivered by a pump to an heat exchanger, absorbing an unknown thermal power and returning output water at a temperature T r (t).In turn, this is heated by the gas boiler and injected again into the tank.
Finally, the electrical circuit is constituted by a locally regulated voltage source V (t) delivering the current I (t), an inductor L, a capacitor C, a time-varying resistor R(t), and a thermal resistance, automatically regulated to absorb the electrical power p r e (t).As evident from Fig. 2, the capacitor, the time-varying resistor, and the thermal resistance are connected in parallel, implying that they are subjected to the same voltage V c (t).
As suggested by Liu et al. [4], the gas temperature dynamics can be neglected given its low thermal capacity.Modeling the mass conservation in subsystem S 1 , it follows that: where c lhv g is the gas lower heating value.Setting , and w in 1 (t) = 0, it is clear that (5) belongs to the system class represented by (1).Moreover, the energy stored in the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 2. Schematic of an illustrative example of a MES S † , composed of a gas storage tank (subsystem S 1 ), a water thermal system (subsystem S 2 ), and an electrical circuit (subsystem S 3 ).
gas tank is and, consequently, the energy derivative follows: It is apparent that, defining the local injected and absorbed power as p g Considering the electrical system S 3 , and using the Kirchhoff's circuit laws, it holds that , and w in 3 (t) = 0. Thus, subsystem ( 8) is represented by a model with the structure of (1).The total energy of the electrical system is stored in the inductor and capacitor and is defined as follows: while its derivative is Indicating the injected and absorbed power as p g , respectively, it is evident that Assumption 1 holds.
Finally, modeling the thermal dynamics in the water vessel, assuming it to represent the dominant inertia in the water system, it follows that: where k wm is the conduction heat transfer coefficient and c w is the water specific heat capacity.The terms η gw and η ew represent the power conversion efficiencies from the gas and electrical systems to the water system, respectively.Setting ] ′ , also system (11) is represented by the model structure proposed in (1).Moreover, the total stored energy in the water vessel is its derivative is Defining the power injected by the thermal plate as p g 2 (t) = k wm A(T m (t) − T w (t)) and the one absorbed by the heat exchanger as p l 2 (t) = c w q w (t)(T w (t) − T r (t)), it is apparent that Assumption 1 holds.
Many other systems and energy domains can be formulated with the model class ( 1)-( 4) beyond the ones considered in the described example (e.g., mechanical, hydraulic systems, and so on), as well as more complex system architectures (e.g., electrical, thermal, or gas networks).The formulated example shows which are the possible interactions and control challenges in a MES, where both local control and coordination of power exchanges must be performed in a cost-effective and scalable manner.
Remark 1: It is worth noting that the multi-energy interaction may give many advantages from the energetic and economic point of view, but the achieved benefit could be not a priori equally shared among participating subsystems.In fact, considering the illustrative example in Fig. 2, transferring power to the water thermal system implies an additional cost both for the gas and electrical system, which may not want to afford without any revenue.This motivates the necessity of a benefit sharing algorithm, which Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
will be described in Section IV.The formulation of the proposed hierarchical control architecture is described in Section III.

III. HIERARCHICAL PREDICTIVE CONTROL OF MESS
Considering the hierarchical control architecture shown in Fig. 1, each L-MPC i is executed with sampling period τ i , while the HL-MPC runs with sampling period τ h ≥ τ i , ∀i ∈ N .For the sake of simplicity, τ h is chosen to be a common multiple of τ i , ∀i ∈ N , i.e., ∃H i ∈ N, so that τ h = H i τ i .
A discrete time index k i ∈ N is introduced for each subsystem S i , defining the time instants t = k i τ i when the L-MPC i is executed, while k h ∈ N is used for the HL-MPC layer, running at time instants A receding horizon strategy is adopted for the MPC control design, using N i prediction steps for L-MPC i and N h prediction steps for HL-MPC.This one, providing the optimal energy references to the L-MPC regulators, is designed with a longer prediction horizon with respect to each L-MPC i , precisely setting

A. High-Level Coordination of a MES
The HL-MPC layer considers a reduced energy model for each S i , formulating a model structure independent of the specific energy domain.Denote with p gh i (k h ), p lh i (k h ) ∈ R the total generated and absorbed power in S i at k h , respectively.Moreover, use e h i (k h ) ∈ R to indicate the total energy stored in S i at k h .Thus, (4) can be discretized with sampling period τ h as follows: enabling to obtain a dynamical model representing S i from the energetic point of view.
The main objective of HL-MPC is to optimize the energy stored in each S i and the power exchanges among S 1 , . . ., S M , considering the minimization of an overall production cost.To perform these tasks, the HL-MPC is assumed to have available the following aggregated information for each S i : 1) the bounds of the total storable energy ēh i , e h i ∈ R; 2) the required energy nominal level e ho i ∈ R; 3) the bounds of the total generated power pgh i , p gh i ∈ R; 4) the total generation cost function, i.e., c gh i : R g i → R ≥0 , expressing the total generation cost as a function of p gh i , i.e., c gh i ( p gh i ).Moreover, the HL-MPC receives at time instants t = k h τ h the following variables from each L-MPC i : 1) the total energy stored at k h in S i , denoted as 2) the expected total power absorption in S i , denoted as For the sake of clarity, the detailed definition of the mentioned variables will be given in Section III-B, where the L-MPC i problems are formulated.
Letting the prediction horizon used in HL-MPC defined by T h , the high-level control problem is formulated as follows: min 14) and The cost function (15a) comprises the overall cost for the produced power in each S i , the one to implement power exchanges among interconnected subsystems, and a terminal cost with coefficient λ i ∈ R ≥0 included to prevent excessive variations of e h i with respect to the nominal value e ho i .The power exchanges cost is expressed using the term c ph i j ∈ R ≥0 weighting the absolute value of p h i j , which is bounded between p h i j ≤ 0 and ph i j ≥ 0 in (15e), ∀(i, j) ∈ E. On the other hand, the variables of each S i are constrained by the total power and energy bounds (15b) and (15c), and by the initial energy state (15d).
The solution to (15) is the optimal sequence over T h of the total power production in each S i and of the power exchanges among S 1 , . . ., S M .The optimal energy to be stored in each S i , indicated with e h * i (k|k h ), k ∈ T h , can be obtained iterating ( 14) and is sent as a reference to each L-MPC i .The same holds for the optimal power exchanges computed at k h , indicated with p h * i j (k|k h ), k ∈ T h , where, in particular, the first step of the sequence is directly imposed by the HL-MPC according to the receding horizon approach, i.e., p i j (t Finally, the power production in each S i is reoptimized by the corresponding L-MPC i with sampling period τ i ≤ τ h , considering local constraints, costs, and disturbances.

B. L-MPC Regulation for Single Energy Subsystem
Each L-MPC i regulates S i considering its discretized dynamics, i.e., where f d i derives from (1) applying a suitable integration method.Letting the prediction horizon used in L-MPC i be T i = {k i , . . ., k i + N i − 1}, the local control problem formulation follows: s.t.∀k ∈ T i , to (16) and The L-MPC i cost function (17a) is constituted of two terms.The first term expresses the local power production cost, with c g i ∈ R g i ≥0 , including the cost of each single power generation in S i , and the second term minimizing the slack variable This slack variable σ i is introduced to ensure problem feasibility while imposing the HL-MPC energy reference through (17f)-(17i), where it is recalled that, given the definition of τ h , k h = ⌊(k i /H i )⌋, ∀i ∈ N .The optimal power exchanges are included as hard constraints in (17j) and (17k), as these are directly imposed by the HL-MPC.The L-MPC i regulates S i considering its discrete dynamics in (16) and the local states and inputs constraints in (17b) and (17c).Local power injections, defined in (17d), are bounded in (17e).Moreover, it is assumed that a prediction of local disturbances over T i is available; otherwise, they are considered constant.
Each L-MPC i enables the fine regulation of S i , considering local production costs and constraints, while tracking the HL-MPC references of energy and power.Precisely, the first step of the optimal input sequence computed by (17), i.e., Definition of Aggregated Variables of the Energy Subsystems: As mentioned in Section III-A, the HL-MPC requires some aggregated information regarding each S i , namely, the following hold.
1) The bounds of the total stored energy in S i , which are the functions of the local state ones.In particular, if the function ξ i (•) is monotonically increasing (as the ones of the example described in Section II-B), it follows that e h i = ξ i (x i ) and ēh i = ξ i ( xi ).
2) The nominal energy level e ho i in S i can be defined as a function of state nominal values, denoted as x o i , and implying that e ho i = ξ i (x o i ).
3) The bounds of total generated power in S i , which can be computed as the sum of the local bounds, i.e., pgh i = 1 ′ g i pg i and p gh i = 1 ′ g i p g i .4) The overall cost of the total generated power in S i , which can be derived through different methods.A simplistic solution is to exploit an average cost, as in [35], i.e., c gh i ( , with g i representing the number of generators in S i .On the other hand, a more precise convex quadratic approximation of c gh i ( p gh i ) can be computed, using the procedure reported [36].Finally, an alternative method exploiting the structure of problem (17) and leading to a piecewise definition of On the other hand, the HL-MPC requires also to measure specific variables from each L-MPC i at each t = k h τ h , defining the energy and power conditions of S i .These are the following.
1) The total energy stored, which can be derived from the local state variables as 2) The expected load power absorption in S i , denoted as p lh i (k) with k ∈ T h , estimated exploiting (2b) and assuming to have available predictions on the future states, inputs, and disturbances in S i , e.g., given by the optimal solution of ( 17).Thus, assuming to know future predictions of local load power absorptions in S i and indicating it with The proposed hierarchical control architecture can be implemented with a fully distributed scheme, e.g., similar to [37], avoiding to share even aggregated information of each S i to a central entity.In fact, the HL-MPC control problem (15), exploiting reduced energy models for S 1 , . . ., S M , is convex if the overall costs c gh i ( p gh i ) are defined with convex functions as well (e.g., using the method proposed in [36]).The convexity of (15) would enable its fully distributed computation with guaranteed optimality and convergence, e.g., using the algorithms proposed in [38] and [39].
Remark 3: The analysis of the recursive feasibility of ( 15) and ( 17) is a complex task due to model nonlinearities and disturbance uncertainties.From a practical point of view, the feasibility of ( 15) can be assumed provided that S 1 , . . ., S M are equipped with sufficient generating power to keep their total energy in the prescribed bounds, considering also the possibility of multi-energy power transfers.In particular, if subsystems are designed, such that p gh i ≤ p lh i ≤ pgh i ∀i ∈ N , it can be shown that a feasible solution to (15) always exist, because each S i can balance its own power absorption without exploiting power exchanges.
The feasibility analysis of ( 17) is more complex, mainly due to state constraints (17b).However, note that the HL-MPC computes the optimal power transfers considering total energy and power constraints for each S i , which bounds depend on the local states and power limits, as previously described.Thus, the feasibility of ( 17) can be assumed provided that local inputs are sufficiently sized to ensure that (17b) and (17e) are respected.It is also worth noting that, as usual in MPC approaches [40], the recursive feasibility can be guaranteed by the use of additional slack variables in the formulation of constraints (15b) at the high level and of (17b) and (17e) at the low level.

IV. BENEFIT SHARING IN MESS
The proposed hierarchical control architecture enables the optimized control of a MES S † , but it does not consider how the achieved benefit is shared among subsystems S 1 , . . ., S M .In fact, the multi-energy operation can imply an additional cost to few subsystems and an increased gain to others, even though the overall operational cost of S † decreases with respect to the independent regulation of each S i .This benefit sharing problem is here recast in the framework of cooperative game theory, 3 formulating it as a bargaining problem [29].
A bargaining problem, or game, involves a group of players having to reach an agreement about a set of feasible payoff solutions, with the common objective of maximizing their own utility.In particular, it involves M players and a pair ( , δ), where ⊆ R M identifies the set of feasible players payoff, whereas δ ∈ R M corresponds to a given disagreement outcome, i.e., the one that players achieve if they do not cooperate.It is assumed that ∃θ ∈ , so that θ > δ, meaning that there exists at least one cooperative feasible payoff solution exceeding the disagreement outcome for each subsystem.As described in detail in [41], there is a unique bargaining solution possessing specific common properties (e.g., the independence of equivalent utility representations).This is the one maximizing the weighted geometric mean of players utilities, 4 i.e., max where α i ∈ R >0 is a parameter denoting the bargaining power of the ith player, defined so that M i=1 α i = M.In particular, if α i = 1, ∀i ∈ {1, . . ., M}, the solution of ( 18) is named symmetric Nash bargaining solution; otherwise, it is named nonsymmetric Nash bargaining solution [42].Now, consider our case study and assume that the interaction over multiple energy vectors among subsystems S 1 , . . ., S M achieves an overall cost lower than the one achieved by their independent regulation.Moreover, suppose that subsystems agree with a time period τ b on how the achieved benefit must be shared.For the sake of simplicity, τ b is chosen, so that τ b = B i τ i , ∀i ∈ N , with B i ∈ N, and τ b = G h τ h , with G h = B i /H i > 1.Moreover, define with k b ∈ N the discrete time index used in the benefit sharing problem and with t = k b τ b the time instants where agreement must be 3 For the sake of simplicity, it is here assumed that all subsystems in S † take part in the multi-energy cooperation and share the achieved benefit.The scenario where some subsystems decide to not join the coalition S † is not here addressed, being the overall cooperation the most convenient scenario (for coalitional approaches, the reader is referred to [26] and [28]). 4For notational simplicity, the Mth power of the weighted geometric average is used in (18), since this is maximized by the same point of the weighted geometric average [41].
reached for the time period where p g * i (k|k) is the power generation solution implemented by the L-MPC i , i.e., the first step of the optimal sequence computed at t = kτ i by problem (17).On the other hand, as a convention, it is chosen in (19) that the costs to implement power transfers, if present, are applied to the receiving subsystems.It is worth noting that the cost incurred by each S i over the period [(k b −1)τ b , k b τ b ] could also be directly derived from local measurements, e.g., by substituting, in (19), the optimal generated and transferred power with the effectively measured ones.
To evaluate the benefit introduced by the multi-energy operation, assume that each S i can be independently regulated, i.e., local input and generation capability are sufficient to satisfy the internal demand without exploiting power exchanges.Thus, the cost that each S i would have incurred if independently regulated can be computed by solving, for each s.t.∀k ∈ T i , to ( 16), (17b)-(17e) and This corresponds to the L-MPC i problem in case each subsystem S i does not exchange energy with the others, aiming solely to the minimization of its internal cost.The optimal power generation computed by ( 20) is denoted with p g * * i (k|k).Therefore, the cost that each S i would have afforded over It is possible now to define the effective benefit that each S i has achieved because of the multi-energy interaction over the time period which is positive in case subsystem S i achieves a lower cost while exchanging energy with the others.A positive benefit is not guaranteed for each S i , while a total positive benefit is assumed to be achieved in S † , i.e., Thus, the total benefit J tot (k b ) must now be properly shared among S 1 , . . ., S M , so that all players take advantage from the multi-energy interaction.To do that, introduce π i j (k b ) ∈ R, ∀(i, j) ∈ E, expressing a payment transfer from S i to S j .Then, the following problem is solved at each t = k b τ b : max where, for the sake of clarity, the dependency on k b has been removed.The variable z i represents the net payment transfer of each S i , defined as the difference between the incoming and outgoing payment transfers from, or to, its neighbors, as dictated by (24c).Equation (24d) states that the sum of net payment z i and of the benefit J i must be positive for each S i ; otherwise, players would be not incentivized to cooperate.
It is worth noting that (24c) implies that meaning that the sum of the net payments of S 1 , . . ., S M is balanced, being these used to redistribute the achieved cooperation benefit among subsystems, i.e., z ∈ Z = {z ∈ R M : 1 M ′ z = 0}.Considering the bargaining problem formulation in (18), θ i = z i represents the payoff of each player S i , while δ i = − J i identifies its disagreement outcome, i.e., the cost loss (or gain) achieved if subsystems S 1 , . . ., S M do not cooperate [see (22)].
As mentioned before, α i in (24a) is a parameter expressing the bargaining power of each S i .If α i = 1, ∀i ∈ N , the bargaining problem is symmetric, meaning that the optimal solution of (24) will provide all players the same utility, i.e., Note that maximizing the utility of each S i in (24) actually corresponds to maximizing the cost reduction achieved by each subsystem with respect to its independent regulation.In fact, it holds that On the other hand, the total cooperation benefit could be shared according to different strategies exploiting the bargaining power parameters α i .For instance, α i can be defined proportionally to the operational cost of each S i , thus providing a higher cost reduction to subsystems characterized by higher expenses.
Remark 4: The benefit sharing problem ( 24) is assumed to be solved by a central entity, e.g., by an additional control layer or by the HL-MPC itself, as it requires as inputs the benefit achieved by each subsystem S i , i.e., J i .On the other hand, operations ( 19)-( 22) are separately solved by each S i , e.g., by the L-MPC i themselves.It is worth noting that the bargaining problem can be solved in a very reduced computational time and with a large sampling period, being a static problem, as it will be evident in Section V.
Remark 5: A different strategy that can be exploited for the benefit sharing problem relies on the Shapley value method [26].First of all, consider the set of all players S † = {S 1 , . . ., S M } and the achieved benefit by their joint regulation J tot (S † ) = M i=1 J i > 0, defined in (23) (in the following, the time index k b is removed for notation clarity).The idea of this technique is to redistribute the total achieved benefit J tot (S † ) among all subsystems S 1 , . . ., S M based on their effective contribution to the coalition.To do that, consider any possible coalition C l ⊆ S † , with l ∈ {1, . . ., 2 M − 1}, and denote with J tot (C l ) the benefit achieved by the joint regulation of C l .In particular, given a generic coalition C l , J tot (C l ) can be computed as in (23) by applying the proposed control strategy to the overall S † while adding the constraint p h i j (k) = 0, ∀S i , S j / ∈ C l to (15), i.e., imposing that power exchanges can occur just within C l .Now, suppose to have computed J tot (C l ) for any l ∈ {1, . . ., 2 M − 1}.Then, the Shapley value assigns the final utility to each S i , denoted as χ * i , through the following formula: where The final utility assigned to each S i through the Shapley value method is a fraction of the total benefit, implying that Considering the formulation of the final utility adopted in the bargaining problem, i.e., z * i + J i , the final net payment transfer for each S i using the Shapley value is z * i = χ * i − J i , with J i defined by (22).Conversely to the solution of the bargaining problem (24) with α i = 1 ∀i ∈ N , now, the benefit will be not distributed equally but based on the effective contribution of S i to the coalition, i.e., z * Moreover, while the proposed bargaining method implies to just solve (24), the Shapley method requires to solve 2 M − 1 optimization problems to compute (25), as J tot (C l ) for any possible coalition of subsystems is necessary, making this method often computationally intractable [27].

V. CASE STUDIES
The described hierarchical control architecture and benefit sharing algorithm have been tested on two MES case studies.First, the illustrative example described in Section II-B is considered, showing how the highlighted challenges are solved by the proposed control strategy.Then, a realistic MES benchmark is considered, where two energy networks referenced in the literature, i.e., an electrical and a thermal one, are effectively regulated by the proposed control architecture.

A. Control of the MES Illustrative Example
The MES illustrative example, depicted in Fig. 2, has been simulated in MATLAB numerically integrating the dynamical equations described in Section II-B and using the physical parameters reported in Table I.The system has been tested under the regulation of the proposed hierarchical control Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE I PHYSICAL PARAMETERS OF THE MES EXAMPLE TABLE II CONTROL PARAMETERS FOR THE MES EXAMPLE
architecture, implementing an L-MPC regulator for each subsystem, i.e., L-MPC 1 for the gas storage tank (subsystem S 1 ), L-MPC 2 for the water thermal system (subsystem S 2 ), and L-MPC 3 for the electrical circuit (subsystem S 3 ).The HL-MPC, placed at the upper layer, provides the energy set points to the L-MPC regulators and optimizes the power exchanges among the different subsystems, thus regulating the gas boiler (power transfer from S 1 to S 2 ) and the thermal resistance (power transfer from S 3 to S 2 ).The control problems ( 15) and ( 17) have been implemented in MATLAB and solved using CasAdi and the Interior-Point Optimizer (IPOPT) solver, which are particularly suited for optimizing nonlinear dynamical systems [43].The main control parameters are reported in Table II (please refer to Section II-B for the physical meaning of the subsystems variables).Please note that γ i = 10 3 and λ i = 10 3 , ∀i ∈ N .Moreover, note from Table II that each L-MPC i regulator runs with a specific sampling period and prediction horizon, chosen consistently with the local requirements of each subsystem.On the other hand, the HL-MPC layer is designed to operate with τ h = 15 min and N h = 16, optimizing the MES over a prediction horizon larger than the L-MPC regulators ones.The generation costs for each subsystem are depicted in Fig. 3(a), where c g 1 is the cost for power injected in the gas tank, c g 2 is the cost for the power generated by the thermal plate, and c g 3 is the cost for the electrical power generated by the voltage source.The costs of power interfaces are c ph 12 = 0.1 $/kW for the gas boiler and c ph 31 = 0.01 $/kW for the thermal resistance.Moreover, the total costs c gh i have been derived using the procedure described in Appendix A. The MES example benchmark is tested with time-varying disturbances, which normalized profiles are depicted in Fig. 3(b), since they are characterized by different units of measure.Fig. 4 reports a 24-h simulation of the considered case study, whose numerical results are reported as normalized with respect to their maximum bounds for the sake of clarity.It is evident from Fig. 4(a) that the HL-MPC optimal coordination involves a continuous power transfer to S 2 , principally from S 1 in the first half of the simulation and from S 3 in the second half.This is mainly due to the generation costs, being S 2 characterized always by the higher one [see Fig. 3(a)].While optimally coordinating power transfers among subsystems, the hierarchical control architecture ensures that all local constraints are respected.As an example, consider the time period between 16:00 and 20:00, characterized by an increase of the load water flow q w in S 2 [red dashed-dotted line in Fig. 3(b)].This implies an increment of the power transfer from S 3 and S 1 to S 2 in the same time period [see Fig. 4(a)], being local generation in S 2 still too costly.Nevertheless, S 1 is also characterized by a high output gas flowrate q g , and an extra power transfer to S 2 may cause to violate local state constraints in S 1 .Thus, the HL-MPC requires to increase the stored energy in S 1 between 12:00 and 16:00, as evident from Fig. 4(b), so that it can support S 2 by transferring the necessary power between 16:00 and 20:00 without violating the local states bounds.The regulation of the L-MPC controllers is depicted in Fig. 4(c), ensuring that local inputs and states are optimally controlled at their specific time rates and maintained within the desired bounds.Please note that there is a clear correspondence between the HL-MPC energy references and the local states, e.g., the variation of stored energy in S 1 directly implies a variation of the stored gas density [see Fig. 4(b) and (d)].
The benefit sharing problem discussed in Section IV has been solved through the bargaining method at the end of the overall simulation, i.e., using τ b = 24 h, and considering a symmetric bargaining, i.e., α i = 1, ∀i ∈ N .The values of the cost functions ( 19) and ( 21) are reported in Table III.As evident, the optimal MES coordination is economically convenient for the total MES S † , but it could be not convenient for each single subsystem.As evident, the water thermal system S 2 is the only one achieving a cost reduction with respect to its independent regulation (J c 2 < J nc 2 ), being it supported with the extra power produced by the other subsystems.In fact, solving the benefit sharing problem (24), its optimal solution implies a payment transfer z * i from S 2 to S 1 and S 3 , eventually resulting in a MES interaction economically convenient for all subsystems.As clear from Table III, the effective final cost, i.e., J c i −z * i , is lower than the one achieved by the independent regulation, i.e., J nc i , for each subsystem belonging to S † .For the sake of comparison, the benefit sharing problem is solved also through the Shapley value method, described in Remark 5.In particular, calculating the achieved benefit by each possible coalition (see Remark 5), it is obtained that:  III, it is evident that the Shapley method assigns to the water thermal system S 2 a lower final effective cost with respect to the bargaining method.This is motivated by the fact that the Shapley method redistributes the benefit based on the contribution of each S i to the coalition, and S 2 is the element that actually enables the power exchanges among the different subsystems and, consequently, leads to cost saving.Note that, given the benefit problem formulation, also the Shapley value method ensures that the final cost of each S i , i.e., J c i − z * i , is lower than the one achieved by the independent regulation, i.e., J nc i .Finally, the performances achieved from a centralized controller running with τ centr = 15 min are also reported, computed as described in Appendix B. The obtained optimal cost functions are J centr 1 = 473.14,J centr 2 = 30.59,and J centr 3 = 2356.6,whereas the total cost is ∀i J centr i = 2860.3.Comparing these values with J c i in Table III, it is evident that the suboptimality introduced by the proposed hierarchical strategy is quite limited, i.e., ( ∀i

B. Control of a Realistic Networked MES Benchmark
The proposed hierarchical control architecture has been also tested considering the realistic multi-energy scenario depicted in Fig. 5, where two energy networks are jointly coordinated, i.e., the AROMA DHN [3] and the IEEE 37-bus electrical network [44].According to the modeling framework presented in Section II, we denote the IEEE 37-bus electrical subsystem as S 1 and the AROMA DHN subsystem as S 2 .The two subsystems are energetically connected through a heat pump, which converts electrical power into thermal one, transferring energy from S 1 to S 2 , as clear from Fig. 5.As a consequence, S † can be represented by a graph G(N , E) with N = {S 1 , S 2 } and E = (S 1 , S 2 ).
The detailed models of S 1 and S 2 are not here reported in detail for the sake of compactness, being them relying on wellknown physical equations for electrical and thermohydraulics network systems [3], [18].Nevertheless, the reader is referred to [18] and [45,Ch. 1] for the detailed modeling of S 1 and S 2 .Hereafter, a general description of each S i in S † and of the involved main variables is given.
The IEEE 37-bus electrical subsystem S 1 is equipped with three storage batteries (nodes E5, E21, and E30), two photovoltaic generators (nodes E15 and E32), 12 loads (nodes E6-E9, E16-E19, and E33-E36), and a converterinterfaced connection with the main utility (node E26), as depicted in Fig. 5.The vector of local states x 1 (t) comprises the state of charges of batteries, which dynamics is described through the battery integrator model [46].The local control inputs are the active power absorbed (delivered) from (to) the main utility, i.e., p g 1 (t) = u 1 (t) (defined to be positive if absorbed from the main utility), and the active power output of batteries in nodes E5 and E30, i.e., p g 2 (t) = u 2 (t) and p g 3 (t) = u 3 (t) (defined to be positive if injected into the grid).In fact, the battery in node E21 is not power-controlled, but it is regulated to maintain a fixed voltage at its node (similar to [46]), implying that its output power automatically ensures the overall power balance.The load power absorption and the photovoltaic power production are local disturbances, and so compacted in d 1 (t).The electrical power absorbed by the heat pump (node E11) is w 1 (t) = p 12 (t), whereas the total energy in S 1 is the sum of the ones stored in the batteries, implying that e 1 (t) linearly depends on x 1 (t).On the other hand, the electrical model of S 1 is nonlinear, being it based on the power-flow equations [18].
The AROMA DHN subsystem S 2 is equipped with a thermal boiler, a thermal storage (with a constant water mass), Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 5. MES realistic case study: IEEE 37-bus electrical network [44] (subsystem S 1 ) and AROMA DHN [3] (subsystem S 2 ).and eight thermal loads (see Fig. 5).In particular, the DHN is constituted of two parallel layers: an upper network, which supplies hot water to the loads, and a lower returning layer where the output load water flows back toward the heating stages.The DHN and the thermal storage are modeled through the finite-volume method, where thermal pipes and the thermal storage are partitioned in finite-volume sections parallel to the water flow direction, and each finite-volume section is described by a local lumped temperature model [3], [47].On the other hand, the thermal boiler and loads are described through static models, given their smaller time constants with respect to the thermal network ones [48].Considering the model of S 2 , the local state x 2 (t) = [x 21 (t) ′ x 22 (t) ′ ] ′ comprehends the water temperature in the sections of the heating supply network, included in x 21 (t), and in the ones of the thermal storage, included in x 22 (t).Moreover, as common in practice [48], the thermal boiler is temperature-regulated, and the inlet/outlet water flowrate of the thermal storage is regulated through a local pump.Thus, the local inputs u 2 (t) = [u 21 (t)u 22 (t)] ′ are the output temperature of the thermal boiler, denoted with u 21 (t), and the water flowrate for the thermal storage, denoted as u 22 (t) (defined as positive if hot water from the upper network layer is injected in the thermal storage).The local power generation p g 2 (t) includes the power injected by the thermal boiler, whereas the vector of disturbances d 2 (t) is the power absorptions of thermal loads.The energy e 2 (t) stored in S 2 is a linear function of x 2 (t), i.e., of the temperatures of pipes and storage sections, as in (12) for the MES illustrative example.Nevertheless, the model of S 2 is nonlinear, mainly due to transport dynamics in thermohydraulic networks and their nonlinear dependence on the water flowrate.
The physical parameters of the IEEE 37-bus electrical network and of the AROMA DHN, as well as their topologies, are reported in detail in [3] and [44], whereas additional parameters regarding generation and storage units are available in Table IV.The power profiles of the electric loads and of the photovoltaic generators are depicted in Fig. 6(a) and (b), respectively, whereas the power absorbed by thermal loads is reported in Fig. 6(c).The cost of the electrical power exchanged with the main utility and the one of the gas power consumed by the thermal boiler are reported in Fig. 6(d).
The main design parameters adopted for the numerical tests are reported in Table V, considering also that γ i = 10 3 and λ i = 10, ∀i ∈ N .As clear from Table V, also in this case, the L-MPC regulators operate with different sampling periods and prediction horizons, given the different involved time constants in electrical and thermal systems.On the other hand, the HL-MPC layer is executed at a slower time rate, having a sampling period of τ h = 1 h and a prediction horizon of N h = 24 steps.Also, for the considered case study, control problems have been implemented and solved using CasAdi and the IPOPT solver in MATLAB.The overall MES benchmark is simulated using Simscape, a MATLAB library, which enables the cosimulation of different physical domains (e.g., electrical, thermal liquid, and so on).Fig. 7 reports the numerical results obtained by running the proposed hierarchical architecture on the considered MES benchmark for one day.Considering the electrical energy cost depicted in Fig. 6(d), it is clear from Fig. 7(b) that the HL-MPC coordinates S 1 and S 2 so as to increase the stored energy when the cost is lower (i.e., before 7:00 and between 9:00 and 16:00) and to use the stored energy when the cost is higher (i.e., around 8:00 and 19:00).Coherently, the HL-MPC also regulates the heat pump, so that its maximum power is transferred from S 1 to S 2 when the electrical energy

TABLE V CONTROL DESIGN PARAMETERS FOR THE MES REALISTIC BENCHMARK
cost is lower, as shown in Fig. 7(a).In turn, the L-MPC controllers regulate their own subsystems to track the energy references defined by the HL-MPC, while optimizing their local resources.In particular, L-MPC 1 modulates the power exchanged with the main utility, depicted in Fig. 7(c), so that the total energy stored in batteries in S 1 tracks the HL-MPC reference, as evident from the states of charge profiles depicted in Fig. 7(d), while ensuring the optimal electrical power supply to the heat pump.The L-MPC 2 regulates the output temperature of the thermal boiler, depicted in Fig. 7(e), and the water flowrate of the thermal storage, depicted in Fig. 7(g), so as to store the energy transferred from the electrical system, tracking the energy reference dictated at the HL-MPC layer.In particular, it is worth noticing that most power produced by the photovoltaic generators is not sold to the main utility by S 1 , but it is transferred to S 2 between 8:00 and 16:00, whereas L-MPC 2 , in turn, acts on the local inputs so as to increase the water temperature at each node of the heating supply network [see Fig. 7(f)] and at the different sections of the thermal storage [see Fig. 7(h)].The stored thermal energy is later used, i.e., between 16:00 and 20:00, while keeping the output temperature of the thermal boiler at its lower bound so as to minimize the local gas consumption.
The benefit sharing problem discussed in Section IV is solved for the considered scenario at the end of the day, using τ b = 24 h, and considering a symmetric bargaining, i.e., α i = 1, ∀i ∈ N .The obtained results are reported in Table VI.It is evident that the total operative cost of S † significant decreases when the proposed hierarchical control architecture is applied (i.e., ∀i J c i < ∀i J nc i ).This operation is, in principle, not convenient for all subsystems, as the operative cost of S 1 increases (it delivers additional energy to S 2 instead of selling it to the main utility), whereas the operative cost of S 2 significantly decreases.This issue is solved by running the benefit sharing problem formulated in (18), achieving an effective final cost for each subsystem, i.e., (J c i − z * i ), much lower than the one achievable by their independent regulation.The described case study witnesses that the designed control strategy is effective in enabling the optimal and efficient operation of MES, also when applied to large-scale energy networks.The control problem has been solved on a laptop with an Intel i7-11850H processor, recording an average computational time of 0.03 s for solving the HL-MPC coordination problem, 3.7 s for solving the L-MPC 1 control problem, and 0.35 s for solving the L-MPC 2 one.On the other hand, the computational time for solving the benefit sharing problem in ( 18) is 0.37 s, whereas the one to solve (20) is, for each subsystem, similar to the one of the corresponding L-MPC regulator.

VI. CONCLUSION
This article describes the design of a novel hierarchical predictive control architecture for the coordination of MESs, i.e., interconnected energy systems of different energy domains (e.g., thermal, electrical, or gas).The proposed control strategy involves, at the high level, an energy modeling of reduced order for each MES subsystem, which enables the computational-effective and optimal coordination of power exchanges among MES subsystems.Moreover, at the high level, an energy reference for each MES subsystem is computed, which is then tracked by a low-level MPC regulator regulating local subsystem inputs.On the other hand, a benefit sharing algorithm is designed, with the goal of properly sharing the economic saving achieved by the multi-energy interaction among MES subsystems.As witnessed by the numerical results, the formulated control can be easily applied to different MES benchmarks, being capable of dealing with nonlinear large-scale systems, as most model complexity is addressed by the local low-level regulators and not by the high-level one.The approach described in this article allows for many extensions, currently underway.Among them, the following hold: 1) the design of a distributed control algorithm for the high-level coordination; 2) the formulation of robust and/or stochastic MPC algorithms to deal with possible uncertainties in the prediction of disturbances; 3) the investigation of a more structured game-theoretical framework for the benefit sharing, where MES subsystems may directly compute the optimal energy exchanges through negotiation 8. Cost of total produced power in S i with g i = 3 local power sources.
strategies; and 4) the investigation of more advanced coordination strategies among multi-energy subsystems, e.g., involving demand-response services.

APPENDIX A PIECEWISE DEFINITION OF TOTAL GENERATION COST
It is here presented an alternative method, other than the ones discussed in Section III-B, to approximate the total produced power cost c gh i ( p gh i ), used by HL-MPC in (15), using the local generation costs c g i exploited by L-MPC i in (17a).For the sake of readability and without any loss of generality, this procedure is here described considering p g i = 0 g i .First of all, number power generators in S i in increasing order with respect to their costs, assumed to be different among each other, i.e., c g i1 < c g i2 < • • • < c g ig i .Then, assuming local state and input constraints (17b) and (17c) to be respected, each L-MPC i will use local power generators in increasing order with respect to their index, until their power bounds are reached, minimizing the operational cost.Thus, the overall cost of the total generation in each S i can be modeled at the HL-MPC level as follows: The idea behind this cost formulation can be better appreciated considering Fig. 8, where the function c gh i ( p gh i ) is depicted for a system with g i = 3 power sources.As evident, the overall cost represents, at the HL-MPC level, the fact that generators in S i would be exploited by the L-MPC i regulator in sequence with respect to their cost.

APPENDIX B CENTRALIZED CONTROL PROBLEM
Hereafter, the centralized control problem for the coordination of S † is formulated.Consider a sampling period τ centr and a prediction horizon T centr = {k, . . ., k + N centr − 1}, with N centr > 0. Thus, at each k, a centralized predictive regulator where, as a convention, the cost of power exchanges is applied to the receiving subsystems [as in (19)].

Fig. 1 .
Fig. 1.Schematic of the proposed hierarchical control architecture for MESs.The variables transmitted along the arrows going from each L-MPC to the HL-MPC are not here depicted for the sake of clarity.
}) = 483.4$,and J tot ({S 1 , S 2 , S 3 }) = 652.4$.First of all, these results witness that the cooperation among all subsystems is the scenario that achieves the maximum total benefit.Applying the Shapley value method discussed in Remark 5, the final outcome is z * 1 = 190.48,z * 2 = −426.95,and z * 3 = 235.57,leading to the following effective final cost for each subsystem: J c 1 − z * 1 = 284.8,J c 2 − z * 2 = 462.2,and J c 2 − z * 3 = 2131.5.Comparing these results with the one reported in Table

Fig. 6 .
Fig. 6.(a) Active power absorbed by loads in nodes E6-E9 (blue dashed line), load in nodes E16-E19 (red dashed-dot line), and load in nodes E33-E36 (green solid line).(b) Active power produced by photovoltaic generators in nodes E15 and E32.(c) Thermal power absorbed by loads in nodes T 1, T 2, and T 6 (blue dotted line), by loads in nodes T 3, T 4, and T 7 (red dashed-dotted line), and by loads in nodes T 5 and T 8 (yellow dashed line).(d) Cost of consumed gas power by the thermal boiler (red solid line) and cost of exchanged electrical energy with the main utility (green dotted line).

Fig. 7 .
Fig. 7. (a) Electrical power absorbed by the heat pump.(b) Energy references for S 1 and S 2 defined by the HL-MPC.(c) Electrical power exchanged with the main utility in S 1 .(d) States of charge of batteries in S 1 : node E5 (red solid line), node E21 (blue dotted line), and node E30 (green dashed line); (e) Output temperature of the thermal boiler in S 2 .(f) Temperature at each node of the supply network in S 2 .(g) Water flowrate of the thermal storage in S 2 .(h) Temperature at each layer of the thermal storage in S 2 .
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.shouldsolve the following problem:min u ∀i (•) ∀k∈T i .∀k ∈ T i , to(16), (17b)-(17e) ∀i = 1, . . ., M p i j ≤ p i j (k) ≤ pi j ∀(i, j) ∈ E. , respectively.Similar to Section IV, considering a sampling period τ b > 0 and the generic time instant t = k b τ b , with k b ∈ N, the operational cost of each S i over [(k b − 1)τ b , k b τ b ] is

TABLE III BARGAINING
RESULTS FOR THE MES EXAMPLE

TABLE VI BARGAINING
RESULTS FOR THE MES REALISTIC EXAMPLE