A Fluid Dynamic Approach to Model and Optimize Energy Flows in Networked Systems

: In this paper, attention is focused on the analysis and optimization of energy flows in networked systems via a fluid-dynamic approach. Considering the real case of an energy hub, the proposed model deals with conservation laws on arcs and linear programming problems at nodes. Optimization of the energy flows is accomplished by considering a cost functional, which estimates a term proportional to the kinetic energy of the overall system in consideration. As the real optimization issue deals with an integral formulation for which precise solutions have to be studied through variational methods, a decentralized approach is considered. First, the functional is optimized for a simple network having a unique node, with an incoming arc and two outgoing ones. The optimization deals with distribution coefficients, and explicit solutions are found. Then, global optimization is obtained via the local optimal parameters at the various nodes of the real system. The obtained results prove the correctness of the proposed approach and show the evident advantages of optimization procedures dealing with variational approaches.


Introduction
In the context of energy sectors, primary importance is given to the problems of energy conversion and management ( [1,2]), as well as possible interactions between renewable sources and the environment (see [3] for power planning, [4] for control issues, and [5] for a survey).This especially occurs in cases of multi-generation systems, which are useful for producing electricity as well as hydrogen, heat, and cooling power, with the consequent advantages of high efficiency and reduced CO 2 emissions.
In this direction, scientific communities are focusing on possible models for energy networked systems, with emphasis on optimization problems that arise naturally in daily situations (see [6,7] for possible applications as well as [8] for a complete version).In particular, energy hubs are a class of multi-generation system where multiple energy carriers are converted, stored, and dissipated [9,10].
In this paper, following the ideas proposed in [11,12], an energy hub designed in Waterloo, Ontario, Canada is modeled, and the energy flows at its nodes are optimized.Indeed, this problem was already considered in [11], where the authors provided a controloriented methodology based on a mixed-integer dynamic model and optimal scheduling (see also [13], where the authors accepted suboptimal solutions for some studied data and ensured that the solution remained feasible for data changes), which is robust to uncertainties in specific scenarios.In this case, the main discussion is different, as we want to guarantee robust optimization focusing on one fundamental issue: a description of the spacetime behavior of the energy flows in order to explicitly mitigate the power fluctuations that could affect correct and reliable energy system operation.
In order to achieve this aim, we deal with the conservation of energy flows, a situation that should be considered in most cases involving multi-generation systems.This requirement leads to a continuous model that foresees conservation laws (i.e., partial differential equations (PDEs); see [14]).In particular, we use an approach taken from road traffic and suggested by the Lighthill-Whitham-Richards (LWR) equation ( [15,16]) enriched by junction traffic assignments so as to model with real networks [17][18][19].The outcome is a model that reproduces the main features of car traffic (queue formation and backward propagation) and that, in the case of energy systems, offers the possibility of analyzing the energy flows in each part of the network and for any time.In spite of the apparent simplicity of the LWR equation, this network approach for energy systems represents, to the best knowledge of the authors, a good compromise for focusing on scenarios that are not always in steady states.Indeed, the superiority of the LWR model is also due to the results of existence and uniqueness of solutions for large networks, guaranteeing a solid analytical theory for numerical approximation and optimization problems (similar drawbacks are also described for different types of PDEs in [20][21][22][23] and for energy issues in [24]).For instance, in [25], efficient numerical algorithms are described for treating complex networks in acceptable computational times.On the other hand, for the particular case of supply chains, the authors of [26] provided a systematic presentation of continuous models and their possible numerical computation.This is the main issue for addressing optimization problems, which otherwise would be computationally too expensive.For similar results, we have ref.[27], where urban traffic networks were simulated by considering road sections of various capacities and lengths, and ref. [28], where the authors focused on an optimization strategy for a continuous model of supply chains with queues of unworked materials.Based on this last work, the authors of [29,30] both dealt with an optimal control problem for the minimization of queues of goods and the quadratic difference between a desired outflow and the real one.The former describes an analytical strategy where the input flow is a piecewise constant function, which is assumed as control.The latter aims at optimization via genetic algorithms.Finally, the authors of [31] considered optimization and control problems through machine learning techniques within the context of chemical processes.
Finally, the choice of the LWR approach for energy systems is motivated by three main reasons:

•
Modeling: The LWR model for networks allows the reproduction of features dealing with the spacetime behaviors of energy flows.The same does not happen for richer fluid dynamic models (see, for instance, [32], where the authors dealt with a system of PDEs that described the density of the cars, generalized momentum, which is itself a function of the density, and pressure) and static models, which consider only steady states.

•
Analysis: As for the theory of networks, the LWR model has fundamental and detailed results.To the best knowledge of the authors, there are no similar and complete theoretical developments for other models, especially those of the fluid dynamic type.

•
Numerics and optimization: The robust theory of LWR also gives rise to fast numerical algorithms, (an example is in [33], where the authors proposed simulations of some classical topologies of networks using kinetic schemes) which allow considering complicated optimization strategies.
Following the approach just described, an energy hub is modeled with a finite set of arcs that meet at some nodes (i.e., junctions).The spacetime evolutions of the energy flows within the arcs are found using conservation laws.The dynamics at the junctions is solved via linear programming problems that consider the maximization of the throughflows under constraints that foresee the bounded incoming and outgoing flows and the distribution coefficients that determine how the flows on the incoming arcs are distributed to outgoing ones.
Once assuming the model for energy flows, we consider a cost functional to estimate a term proportional to the kinetic energy in the energy hub.In particular, we want to maximize the functional with respect to (w.r.t.) the distribution coefficients at the nodes.Unfortunately, as it is difficult to foresee a priori an exact evolution of energy flows in the hub, analytical optimization of the cost functional is not possible.For this reason, we adopt a strategy that consists of the following three steps (a similar technique was used in [34]): 1.
For network topologies with a unique node and all initial data, we compute the optimal distribution coefficients.Then, we consider the asymptotic solution, assuming there are infinite-length arcs to avoid boundary data effects.2.
For generic networks with complex topologies, we use the (locally) optimal distribution coefficients at each node, updating the values of the parameters at each time instant through the actual flows on the arcs near the junction.

3.
Using simulations, we verify the performances obtained by the (locally) optimal distribution coefficients through comparisons with randomly chosen parameters.
Notice that the first step is non-trivial even for simple nodes.In fact, we deal with a hybrid problem, as continuous flows are influenced by discrete variables such as the distribution coefficients at the junctions.For this reason, and also considering the topology of the energy hub under discussion, we focus on the particular case of nodes of the 1 × 2 type (i.e., one incoming arc and two outgoing ones).The second step is carried out for the energy hub described in [11], while step three considers two different types of distribution parameters: (locally) optimal, according to the step one, and random (i.e., the coefficients are chosen randomly at the beginning of the simulation and then are kept constant).In particular, numerical approximations for the simulations are obtained via some of the methods for PDEs described in [35][36][37].
Considering that the optimization approach considers a local optimization in an asymptotic state, the obtained results are quite good.Using the optimal parameters, the behavior of the cost functional is better than the ones achieved using random distribution coefficients.Indeed, the approach is also quite robust, as indicated by further analysis of the asymptotic values of the cost functional in random cases and in the optimal situation.This also indicates that the followed approach is suitable for energy hub control as well as the possible scheduling of resources over a long time interval.
In short, the main content of this paper is indicated by the keywords listed as follows.
• Energy flows: These represent the main topic of this research and are modeled via a fluid dynamic approach dealing with conservation laws, namely PDEs of the hyperbolic type.• Optimization: The performances of a real energy system are optimized through a cost functional, which estimates a term proportional to the overall kinetic energy.The cost functional is minimized by using a decentralized approach to find the distribution coefficients at the nodes of the network.• Simulation: The performances of the real energy system are tested via simulations.In particular, the cost functional is studied by considering either optimal distribution coefficients or random ones.Suitable comparisons are made for the different behaviors.
Notice that a possible empirical validation of the proposed model for energy flows is under investigation, with emphasis on other real case studies.Various examples of validations for fluid dynamic models include ref. [38], where the authors discuss validation for a traffic model by considering data gathered in a part of the urban network in Rome (Italy), and refs.[39,40], which discuss possible validation, together with an analysis of emissions and pollutants, for fluid dynamic models of the second order for traffic networks.
This paper has the following structure.Section 2 presents a model for energy flows on networks.Section 3 is devoted to the optimization results for energy networks via the cost functional described above.Section 4 provides an example of a real energy hub designed in Waterloo, Ontario, Canada.Section 4.2 presents some simulations of energy transitions over the energy hub under discussion for two possible choices of distribution coefficients: optimal and random.Our conclusions (Section 5) end the paper.

Theoretical Foundations
An energy network is described by a couple (I, J ), where I ={I n } n=1,...,N is the set of arcs I n , n = 1, . . ., N, with each one represented by [a n , b n ] ⊂ R, and J ={j k } k=1,...,J is the collection of nodes j k , k = 1, . . ., J.
Each arc I n ∈ I, n = 1, . . ., N, is described by the following: where ρ n max is the maximal allowed density for arc , where v n max indicates the maximal velocity for particles traveling on arc The three above quantities have the following interpretation: ρ n is the energy density at time t at the point x of arc I i , v n is the average velocity of each energy particle, and f n (ρ n ) is the flux associated with ρ n .Notice that, as we are dealing with macroscopic quantities, particles are assumed to be of various type.Such an assumption is essential for the case study, where different quantities deal with gas, heat, electricity, and so on.Moreover, the proposed model, dealing with average quantities, considers all the energy phenomena as an average.Hence, there is not a clear distinction (which is typical of microscopic models) between active and reactive powers, which are indeed important for modeling power losses.Focusing only on conservative aspects, the proposed model deals mainly with active powers.The aspects for a more suitable presentation of losses and reactive power are under investigation.
For I n ∈ I, n = 1, . . ., N, the evolution of ρ n (t, x) follows the Lighthill-Whitham-Richards (LWR) model ( [15,16]), which is expressed by the conservation law Without loss of generality, by choosing The evolution at a node j k ∈ J , k = 1, . . ., J, obeys Riemann problems (RPs) (i.e., Cauchy problems that have constant initial data for incoming and outgoing arcs).
We fix a node j k of the r × s type, namely r incoming arcs I k φ , φ = 1, . . ., r and s outgoing ones Definition 1.For the node j k , a Riemann solver (RS) is a function is a weak solution [14] to Equation (1) with the initial datum ρ k i,0 and boundary condition ρ k i such that (C1) RS RS ρ k 0 = RS ρ k 0 and (C2) on arc I k φ , φ = 1, . . ., r, the wave ρ k φ,0 , ρ k φ has negative speed, and on arc I k ψ , ψ = r + 1, . . ., r + s, the wave ρ k ψ , ρ k ψ,0 has positive speed.
Intuitively, for the assigned initial datum ρ k 0 at j k , an RS associates a vector 1) with the initial datum ρ k i,0 and boundary condition ρ k i .More precisely, for φ ∈ {1, . . ., r} and ψ ∈ {r + 1, . . ., r + s}, the solution consists of the wave ρ k φ,0 , ρ k φ on I k φ and the wave ρ k ψ , ρ k ψ,0 on I k ψ .If r ≤ s, then a possible RS at the node j k is constructed using the following rules (see [18]): (A) The traffic of particles distributes at j k according to some coefficients, collected in a The φth column of A j k represents the percentages of particles that, from I k φ , distribute to the outgoing arcs; (B) The flux through j k is maximized w.r.t.rule (A).If r > s, then aside from rules (A) and (B), a further criterion is needed.For instance, if j k is of the r × 1 type, then a possible rule is the following: (C r×1 ) Not all particles enter the outgoing arc.Assume that Q is the quantity that can.
Then, p k φ Q particles come from I k φ to cross the arc junction, where 0 Remark 1.Notice that if r = 1 and s = 2 (i.e., the junction j k is of the 1 × 2 type), then the matrix A j k only has the parameters α k := α From rules (A), (B) and (C) (if this last one is necessary), we find that for a junction j k of the r × s type, with initial datum ρ k 0 and the flux function in Equation ( 2), the solution ρ k to the RP at j k is as follows.Consider the function ω : [0, ρ max ] → [0, ρ max ], which for each arc I k i , i = 1, . . ., r + s satisfies the following properties: Then, for I k φ , φ = 1, . . ., r, we have the following: φ and ρ k ψ , we find the maximal flux values on I k φ , φ = 1, . . ., r and I k ψ , ψ = r + 1, . . ., r + s.In other words, we have Remark 2. Notice that such an approach allows defining solutions to the Cauchy problem on the network (I, J ) via a wave-front tracking algorithm (see [17,18]).Refer to Appendix A for details about the construction of the solutions.

Energy Optimization
Consider an energy network (I, J ) as described in Section 2. For optimization of the network performances, we define the cost functional which represents a term proportional to the kinetic energy on the whole network.
Considering bounded ρ n (t, x), n = 1, . . ., N, the aim is to maximize E (t) w.r.t. the distribution coefficients of the matrices A j k ∀ j k ∈ J .
As the solution of such an optimization control problem involves spacetime variables, and the optimization itself refers only to 1 × 2 nodes for the energy hub described in Section 4, we consider an approach defined by the following steps: 1.
Consider a node j k ∈ J of the 1 × 2 type (one incoming arc (I k 1 ) and two outgoing arcs (I k 2 and I k 3 )) for which only one distribution coefficient α k is considered (see Remark 1).Assuming an initial datum ρ k 1,0 , ρ k 2,0 , ρ k 3,0 at j k , fix the local cost functional to

2.
For a time horizon [0, T], with T being quite large, assume that the traffic distribution coefficient α k is the control, and maximize E j k (T) w.r.t.α k .

3.
Construct the optimal solution to the overall network by localization (i.e., by using the single optimization solutions at each node j k ∈ J of the 1 × 2 type).
For step 2, assume the following conditions: In addition, define We obtain the following: Theorem 1.Consider a node j k ∈ J of the 1 × 2 type, and assume that T is sufficiently large.
Then, E j k (T) is maximized for the value α opt k (for some cases, the optimal control does not exist, and it is approximated through a positive and small constant ε): if T 2 is satisfied; Proof.Assume that v max = ρ max = 1.Other cases for which v max ̸ = 1 and/or ρ max ̸ = 1 lead to the same results.Fix a node j k ∈ J of the 1 × 2 type, where T is quite large.
Considering the solution to the RP at j k , E j k (T) is written as follows: where coefficients s k 1 and s k ψ , ψ = 2, 3, are For simplicity, from now on, we will drop the dependence on j k and T from E .As the solution to the RP at j k depends on α k , we have various cases.Here, for the sake of brevity, we consider only two of them, as the proofs for the other cases are similar.
Assume that γ k,max and we have to maximize where Consider now the case where γ k,max . We then find the following: where

31
, then where It is possible to verify that E is an increasing function.If 1 − g k,max 31 < α k < 1, then ∂E ∂α k is the expression in Equation (6).Hence, we conclude the following: , E does not have an optimal value, then α opt k is chosen for + ε, where ε is a positive and small constant.

Application Deployment 4.1. Energy Hub Operation Scheduling
To assess the effectiveness of the proposed methodology in the task of solving real operation problems in complex networked systems, the problem of optimal energy flow management of a realistic energy hub is considered here.This is a relevant problem in modern energy systems, where the increasing interdependencies between heterogeneous energy infrastructures is introducing new and more complex vulnerabilities, requiring effective modeling and optimization tools aimed at improving the accuracy and robustness of coordinated control actions.In this context, large-scale deployment for the energy hub could have a strategic role, since this allows improving the energy network's flexibility by providing reliable energy services, such as electricity and heating, by exploiting different combinations of the energy carriers available at the hub inputs.To this aim, a system designed in Waterloo, Ontario, Canada for the supply of commercial load was analyzed.More details on this system, which is schematically depicted in Figure 1, can be found in [11].In Figure 1, it is worth noting that the input power flows were electricity P E and natural gas P G , while the output power flows were electricity L E and heat L H . Notice that P E is split into β 1 P E and β 2 P E , while P G is divided into β 3 P G and β 4 P G .To be precise, β 1 P E and β 2 P E are the inputs of the hydrogen production plant (HPP) and the transformer (T), respectively, and β 3 P G is the input of the combined heat power (CHP), while β 4 P G is the input of the furnace (F).The system works as follows:

•
The HPP subsystem, characterized by the electric-hydrogen and heat-hydrogen efficiences η HPP Using the hydrogen-electricity and hydrogen-heat efficiencies η FC 1 and η FC 2 , respectively, the FC subsystem transforms a part of the hydrogen β 1 P E η HPP 1 η FC 1 , into electricity and another part The T subsystem, due to its efficiency η T , has the electricity power flow β 2 P E η T as its output.

•
The CHP subsystem, considering the gas-electric and gas-heat efficiencies η CHP Finally, the F subsystem, characterized by its efficiency η F , has an output As for the outputs of the energy hub, we simply find that Notice that Equations ( 7) and ( 8) are obtained while considering some losses that depend on the coefficients η HPP , and η F .Such parameters are usually fixed, and for the topology described in Figure 1, they are the following: and 0 < β 2 < 1, and hence there is the possibility of choosing how to redistribute incoming flows over the energy hub.
The definition of optimal energy hub operation strategies could follow different directions.An example is described in [11].In our case, for fixed incoming flows P E and P G , we want to find the optimal coefficients β 1 and β 2 (for nodes with one incoming arc and two outgoing ones) such that the global energy conversion losses are minimized.Obviously, different and more general operation criteria could be defined without affecting the effectiveness of the proposed methodology.
In order to achieve this aim, we chose to model the energy hub using an approach dealing with conservation laws on networks (an exhaustive overview is in [18]), considering that eventual losses over the system are considered "outputs" different from L E and L H .
For simplicity of discussion, from now on, we indicate a node j k and its distribution matrix A j k simply by k and A k , respectively.The same applies to an arc I m , named simply m.
Notice that arcs 1 and 14 are the inputs of the systems, while the outputs are arcs 13 and 23, indicated by OUT 2 and OUT 1, respectively.The performances of the network were evaluated through the cost functional E (t), whose evolution is deeply influenced by the distribution parameters.Indeed, due to the real characteristics of the energy hub, for nodes 2, 3, 4, 8, and 9, the distribution matrices A J , J ∈ {2, 3, 4, 8, 9} assume the form Moreover, while always referring to measures on the real hub, the priority parameters were all chosen to be 0.5 for the incoming arcs of nodes 5, 6, 10, and 11.Hence, no control was considered for junctions 2, 3, 4, 5, 6, 8, 9, 10, and 11, and the optimization of E (t) dealt only with the distribution coefficients at nodes 1 and 7 of the 1 × 2 type.
As for the numerical construction of E (t), a suitable approximation for the densities ρ i (t, x), i = 1, . . ., 23 is necessary, and their evolution is ruled by Equation (1).
In this paper, we apply the Godunov scheme (see [35][36][37]), using a numerical grid with constant space and time sizes ∆x = 0.0125 and ∆t = 0.5∆x, respectively (see Section 4.4 for details about the computational cost).The network of Figure 2 was simulated in the following conditions: a time interval [0, T] for the simulation, with T = 150 min; empty arcs when the simulation started (t = 0); and boundary data of the Dirichlet type equal to 0.3 for arcs 1 and 14, while for arcs 5, 7, 12, 13, 17, 19, 21, and 23, we chose Dirichlet boundary data equal to 0.9.
Notice that the typical maximal values for the inputs of our hub were 15 MWh and 20 MWh for arcs 1 and 14, respectively.In our case, by associating ρ max = 1 with the quantity 15 MWh, we simply found that the boundary data 0.3 and 0.9 corresponded to 4.5 MWh and 13.5 MWh, respectively.
Two different choices for the distribution parameters were assumed for nodes 1 and 7: • Optimal case: The parameters that optimize the asymptotic behavior of E (t) locally (i.e., distribution coefficients that refer to Theorem 1 for junctions 1 and 7).Such a type of simulation is useful for testing the global performance, starting from analytical results that consider only part of the nodes of the network.

•
Random case: The parameters at nodes 1 and 7 are chosen in a random way at t = 0 and then are kept constant in [0, T].A random simulation allows comparisons with network performances obtained via the local optimal distribution coefficients.

Results Discussion
In Figures 3 and 4, we show some simulation results for the energy hub.More precisely, the values of E (t), computed on the whole network, are represented as a function of time.
In particular, the behavior of E (t) in the optimal case was compared with the ones obtained via 10 different random simulation studies.As for these last cases, Figure 3 shows the first five, while Figure 4 shows the remaining ones.The optimization algorithm for the 1 × 2 nodes, which are of the local type, can be applied to the complex topology of the energy hub without compromising the possibility of global optimization.Such a situation is evident in the optimal case for E (t), which was compared to the behaviors in 10 different random cases.In Figures 3 and 4, it is shown that the optimal case was always higher than the random ones when t tended toward infinity.This is exactly what Theorem 1 establishes; for long times, the choice of optimal coefficients provides better performances for the energy hub.For short times, it might happen that a random simulation has values for the cost functional higher w.r.t. the optimal case.Indeed, in order to test the validity of the proposed approach, 100 random cases were simulated and compared to the optimal behavior for E (t).Table 1 reports the value of E (t) in the optimal configuration at T = 150 (i.e., OPTconf T), with the average value (RAND T) of random simulations set to T = 150.Notice that OPTconf T was, as expected, higher than RAND T, and the global optimization of the local type had strong robustness.Such a result is also represented in Figure 5, where a histogram reports the values of E (t) at the final instant T for the random simulation.Finally, notice that the simulated system always had bounded outputs as a consequence of the model itself, which deals with limited densities on arcs (see Section 2).In particular, in our case study, the upper limits for arcs 13 and 23 (OUT 2 and OUT 1) were 0.1218 (36.54 kWh) and 0.180075 (54.0225 kWh), respectively.

Computational Cost
This subsection presents some details about the computational cost for the simulation of the presented energy system.We focused on a network represented by the couple (I, J ), with I ={I n } n=1,...,N and J ={j k } k=1,...,J .
In order to find a suitable numerical approximation in [0, T] for the density functions ρ n (t, x), n = 1, . . ., N, on the arcs and updates of the boundary data at the nodes, assume that each arc I n , n = 1, . . ., N has a length L n , and the space and time grid sizes are ∆x n and ∆t n , respectively.Using the Godunov method, the computational cost depends on for the densities and boundary data, respectively.
For simplicity, for the simulation of the described energy hub, we considered a constant space grid size ∆x, assuming that (∆x n , ∆t n ) = (∆x, 0.5∆x) ∀ n = 1, . . ., N, and computed the CPU times (measured in seconds and calculated by an Intel(R) Core (TM) i7-3630 QM CPU @2.40 GHz, 8 GB RAM) and convergence errors.The obtained results are in Table 2.  From the previous table, we simply find that the CPU time increased by about 0.30 s when ∆x decreased.As for the convergence error, it almost remained the same for different values of ∆x.
Notice that the computational effort increased linearly with T and was linearly dependent on the number N of total arcs and the number J of junctions.Hence, the computational time was not influenced by the geometries and topologies of possibly more complex energy networks in a meaningful way.Finally, another interesting topic deals with the scalability of the used approach.Indeed, grid computing and parallelization are a possible alternative for the simulation algorithm.At each iteration of the algorithm, the various solutions at nodes were computed by solving independently linear programming problems at the nodes that implied a separate resolution for the RPs on the arcs.Hence, in order to redistribute the computational load, further instruments dealing with multiprocessing programming are a possibility for simulations.Some further studies, as well as different numerical approaches for conservation laws on networks, were carefully analyzed in [18,25,33].

Conclusions
In this paper, the theoretical foundations of a novel computing paradigm based on the fluid dynamic theory for modeling, analysis, and optimization of complex and networked energy systems were presented.The proposed paradigm is based on the challenging idea of integrating in a unique framework both network modeling and the optimization features, which are traditionally treated as two separate problems and solved by using distinct solution techniques.In order to address this issue, the application of conservation laws and the definition of a cost functional, which represents a term proportional to the kinetic energy of the system, were proposed for network modeling and optimization, respectively.Thanks to these features, the functional maximization is directly obtained as a result of the network model process, which optimally tunes the network parameters ruling the distribution of the energy flows among the network arcs.
The benefits due to the application of the proposed approach were assessed in a realistic case study dealing with the solution of the optimal energy flow management problem for a complex energy hub designed in Waterloo, Ontario, Canada.The obtained results demonstrate the effectiveness of the fluid dynamic-based approach in the task of optimizing the energy flows of the hub components in order to drastically reduce the conversion losses.
Finally, on the basis of repetitive simulations, it could be argued that the obtained solution is globally optimal and robust against distributed parameter variations, and this can be considered relevant in terms of benefits.A rigorous theoretical justification of these intuitions is currently under investigation by the authors, as well as further empirical validation for other real-world scenarios.

Figure 1 .
Figure 1.Simplified structure of the energy hub shown in [11].

1 and η HPP 2 ,
respectively, transforms a part of the electricity β 1 P E η HPP 1 into hydrogen that feeds the fuel cell (FC) and another part β 1 P E η HPP 2 to generate heat.•

2 ,
respectively, transforms a part of the natural gas β 3 P G η CHP 1 into electricity and another part β 3 P G η CHP 2 into heat.•

Figure 2 .
Figure 2. Topology of the energy hub.

Figure 3 .
Figure 3. (left)Evolution in [0, T] of E (t) for optimal distribution coefficients (dashed line) and the first five different random choices (continuous lines).(right) Zoomed-in section near the asymptotic values.

Figure 4 .
Figure 4. (left):Behavior in [0, T] of E (t) in the case of optimal distribution coefficients (dashed line) and the remaining different random choices (continuous lines).(right) Zoomed-in section near the asymptotic values.

Figure 5 .
Figure 5. Histogram of random values of simulations at T = 150 for [0, T] of E (t).The black point represents OPTconf T = 150, and the dashed line represents RAND T = 150.

Table 1 .
Compared values of E (t).