A Fluid Dynamic Approach To Model And Optimize Energy Flows In Networked Systems

1. Introduction

In the context of energy sectors, a primary importance is given to problems of energy conversion/management ([1,2]), as well as possible interactions between renewable sources and environment (see [3] for power planning, [4] for control issues and [5] for a survey). This occurs especially in cases of multi-generation systems, useful to produce electricity as well as hydrogen, heat and cooling power, with consequent advantages for high efficiencies and reduced CO₂ 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] and [7] for possible applications as well as [8] for a complete vision. In particular, energy hubs are a class of multi-generation systems where multiple energy carriers are converted, stored and dissipated, [9,10].

In this paper, following the ideas proposed in [11] and [12], an energy hub, designed in Waterloo, Canada, is modeled and energy flows at its nodes are optimized. Indeed, this problem is already considered in [11], where the authors provide a control-oriented methodology, based on a mixed integer dynamic model and an optimal scheduling (see also [13]) which is robust to uncertainties in specific scenarios. In this case, the main discussion is different, as we want to guarantee a robust optimization focusing on this fundamental issue: a description of the space-time behavior of the energy flows in order to explicitly mitigate the power fluctuations that could affect the 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 Lightwill-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 the car traffic (queues formation and backward propagation), and that, in cases of energy systems, offers the possibility of analyzing energy flows in each part of the network and for every 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 to focus on scenarios that are not always in the steady states. Indeed, the superiority of the LWR model is also due to results of existence and uniqueness of solutions for large networks, guaranteeing a solid analytical theory for numerical approximations and optimization problems (similar drawbacks are also described for different types of PDEs in [20] and [21] and for energy issues in [22]). For instance, in [23], efficient numerical algorithms are described to treat complex networks in acceptable computational times. In particular, this is the main issue to address optimization problems, which otherwise would be computationally too expensive. Similar results are given in [24] and [25]; [26] and [27] deal with possible optimization techniques via genetic algorithms, while [28] considers optimization and control problems through machine learning techniques.

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 space-time behaviors of energy flows. The same does not happen for richer fluid dynamic models (see, for instance, [29]) and/or static models, which consider only steady states.
• Analysis. As for the theory on networks, the LWR model has fundamental and detailed results. To the best knowledge of the authors, there are not similar and complete theoretical developments for other models, especially of fluid dynamic type.
• Numerics and Optimization. The robust theory for LWR also gives rise to fast numerical algorithms (an example is in [30]), which allow to consider complicated optimization strategies.

Following the just described approach, an energy hub is modeled with a finite set of arcs, that meet at some nodes, i.e. junctions. Space-time evolutions of energy flows within the arcs are found using conservation laws. The dynamics at junctions is solved via linear programming problems that consider the maximization of the through flows under constraints that foresee: bounded incoming and outgoing flows; distribution coefficients that determine how flows on incoming arcs distribute to outgoing ones.

Once assumed the model for energy flows, we consider a cost functional to estimate a term proportional to the kinetic energy on the energy hub. In particular, we want to maximize the functional with respect to (w.r.t.) distribution coefficients at nodes. Unfortunately, as it is difficult to foresee a-priori an exact evolution of energy flows on the hub, the 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 is used in [31]):

• For network topologies with a unique node and every initial data, compute the optimal distribution coefficients. Then, assuming infinite length arcs so as to avoid boundary data effects, consider the asymptotic solution.
• For generic networks with complex topologies, use the (locally) optimal distribution coefficients at each node, updating the values of the parameters at every time instant through the actual flows on the arcs near the junction.
• Using simulations, verify the performances obtained by (locally) optimal distribution coefficients through comparisons with random choices of 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 junctions. For this reason, and also considering the topology of the energy hub under discussion, we focus on the particular case of nodes of $1\times 2$ type, i.e. one incoming arc and two outgoing ones. The second step is done 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 simulations are obtained via some methods described in [32,33] and [34].

Considering that the optimization approach considers a local optimization in the asymptotic state, the obtained results are very good: using 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 a 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 the energy hub control, as well as for a possible scheduling of resources over a long time interval.

The paper has the following structure. Section 2 presents a model for energy flows on networks. The subsequent Section 3 is devoted to the optimization results on energy networks via the cost functional described above. Section 4 provides an example of a real energy hub designed in Waterloo, Canada. Subsection 4.2 presents some simulations of energy transitions over the energy hub under discussion for two possible choices of distribution coefficients: optimal and random. Conclusions (Section 5) end the paper.

2. Theoretical Foundations

An energy network is described by a couple $\left(\mathcal{I},\mathcal{J}\right)$, where: $\mathcal{I}={\left\{I_n\right\}}_{n=1,\dots ,N}$ is the set of arcs $I_n$, $n=1,\dots ,N$, each one represented by $\left[a_n,b_n\right]\subset \mathbb{R}$; $\mathcal{J}={\left\{j_k\right\}}_{k=1,\dots ,J}$ is the collection of nodes $j_k$, $k=1,\dots ,$ $J.$

Each arc $I_n\in \mathcal{I},n=1,\dots ,N$, is described by:

• a density function ${\rho }_n:={\rho }_n\left(t,x\right)\in \left[0,{\rho }_{max}^n\right]$, $\left(t,x\right)\in \left[0,+\infty \right[\times \left[a_n,b_n\right]$, where ${\rho }_{max}^n$ is the maximal allowed density for arc $I_n$;
• a velocity function $v_n:=v_n\left({\rho }_n\right)\in \left[0,v_{max}^n\right]$, where $v_{max}^n$ indicates the maximal velocity for particles travelling on arc $I_n$;
• a flux function defined as $f_n\left({\rho }_n\right):={\rho }_n{v}_n$.

The three above quantities have the following interpretation: ${\rho }_n$ is the energy density at time t in the point x of arc $I_i$, $v_n$ is the average velocity of each energy particle, while $f_n\left({\rho }_n\right)$ is the flux associated to ${\rho }_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, that deal with gas, heat, electricity, and so on.

For $I_n\in \mathcal{I},n=1,\dots ,N$, the evolution of ${\rho }_n\left(t,x\right)$ follows the Lighthill-Whitham-Richards (LWR) model ([15,16]), expressed by the conservation law:

$$\frac{\partial }{\partial t}{\rho }_n+\frac{\partial }{\partial x}{f}_n\left({\rho }_n\right)=0.$$

(1)

Without loss of generality, choosing, ∀$n=1,\dots ,N$, ${\rho }_{max}^n={\rho }_{max}$, $v_{max}^n=v_{max}$ and a decreasing velocity function $v_n\left({\rho }_n\right)=v_{max}\left(1-\frac{{\rho }_n}{{\rho }_{max}}\right),$ ${\rho }_n\in \left[0{\rho }_{max}\right]$, the flux function $f\left({\rho }_n\right):=f_n\left({\rho }_n\right)$ simply reads as:

$$f\left({\rho }_n\right)=v_{max}{\rho }_n\left(1-{\displaystyle \frac{{\rho }_n}{{\rho }_{max}}}\right),{\rho }_n\in \left[0,{\rho }_{max}\right].$$

(2)

The evolution at a node $j_k\in \mathcal{J},k=1,\dots ,J$, obeys Riemann Problems (RPs), i.e. Cauchy Problems that have constant initial data for incoming and outgoing arcs.

Fix a node $j_k$ of $r\times s$ type, namely r incoming arcs $I_{\phi }^k$, $\phi =1,\dots ,r$, and s outgoing ones $I_{\psi }^k$, $\psi =r+1,\dots ,r+s$, and indicate by ${\rho }_0^k=\left({\rho }_{1,0}^k,\dots ,{\rho }_{r,0}^k,{\rho }_{r+1,0}^k,\dots ,{\rho }_{r+s,0}^k\right)\in {\left[0,{\rho }_{max}\right]}^{r+s}$ the initial datum at $j_k$.

Definition 1.

For the node $j_k$ a Riemann Solver (RS) is a function

$$RS:{\left[0,{\rho }_{max}\right]}^r\times {\left[0,{\rho }_{max}\right]}^s\to {\left[0,{\rho }_{max}\right]}^r\times {\left[0,{\rho }_{max}\right]}^s$$

such that

$${\rho }_i^k(t,x):=\left\{\begin{array}{cc}({\rho }_{i,0}^k,{\widehat{\rho }}_i^k),\hfill & I_i^k,i=1,\dots ,r,\hfill \\ ({\widehat{\rho }}_i^k,{\rho }_{i,0}^k),\hfill & I_i^k,i=r+1,\dots ,r+s,\hfill \end{array}\right.$$

is a weak solution [14] of (1) with initial datum ${\rho }_{i,0}^k$ and boundary condition ${\widehat{\rho }}_i^k$, and such that furthermore: (C1) $RS\left(RS\left({\rho }_0^k\right)\right)=RS\left({\rho }_0^k\right)$; (C2) on arc $I_{\phi }^k$, $\phi =1,\dots ,r$, the wave $\left({\rho }_{\phi ,0}^k,{\widehat{\rho }}_{\phi }^k\right)$ has negative speed, and on arc $I_{\psi }^k$, $\psi =r+1,\dots ,r+s$, the wave $\left({\widehat{\rho }}_{\psi }^k,{\rho }_{\psi ,0}^k\right)$ has positive speed.

Intuitively, for the assigned initial datum ${\rho }_0^k$ at $j_k$, a RS associates a vector ${\widehat{\rho }}^k=\left({\widehat{\rho }}_1^k,\dots ,{\widehat{\rho }}_n^k,{\widehat{\rho }}_{n+1}^k,\dots ,{\widehat{\rho }}_{n+m}^k\right)\in {\left[0,{\rho }_{max}\right]}^{r+s}$ such that on $I_i^k$, $i=1,\dots ,r+s$, ${\rho }_i^k$ is a solution of (1) with initial datum ${\rho }_{i,0}^k$ and boundary condition ${\widehat{\rho }}_i^k$. More precisely, for $\phi \in \{1,\dots ,r\}$ and $\psi \in \{r+1,\dots ,r+s\}$, the solution consists of the wave $\left({\rho }_{\phi ,0}^k,{\widehat{\rho }}_{\phi }^k\right)$ on $I_{\phi }^k$ and the wave $\left({\widehat{\rho }}_{\psi }^k,{\rho }_{\psi ,0}^k\right)$ on $I_{\psi }^k$.

If $r\le s$, a possible RS at the node $j_k$ is constructed using the rules (see [18]):

(A) 

The traffic of particles distributes at $j_k$ according to some coefficients, collected in a matrix $A_{j_k}=\left({\alpha }_{\psi ,\phi }^{j_k}\right),$ $\phi =1,\dots ,r,$ $\psi =r+1,\dots ,r+s$, $0<{\alpha }_{\psi ,\phi }^{j_k}<1,$ ${\displaystyle \sum _{\psi =r+1}^{r+s}}{\alpha }_{\psi ,\phi }^{j_k}=1$. The $\phi -$th column of $A_{j_k}$ represents the percentages of particles that, from $I_{\phi }^k$, distribute to the outgoing arcs;

(B) 

The flux through $j_k$ is maximized respecting rule (A).

If $r>s$, beside rules (A) and (B), a further criterion is needed. For instance, if $j_k$ is of $r\times 1$ type, a possible rule is the following:

(C_r×1) 

Not all particles enter the outgoing arc, and assume that Q is the quantity that can do it. Then, $p_{\phi }^{k}Q$ particles come from $I_{\phi }^k$ to cross the arc junction, where $0<p_{\phi }^k<1$, ${\displaystyle \sum _{\phi =1}^r}{p}_{\phi }^k=1$, is the priority parameter of $I_{\phi }^k$, $\phi =1,\dots ,r$.

Remark 1.

Notice that, if $r=1$ and $s=2$, i. e. the junction $j_k$ is of $1\times 2$ type, the matrix $A_{j_k}$ has only the parameters ${\alpha }_k:={\alpha }_{2,1}^{j_k}$ and $1-{\alpha }_k:={\alpha }_{3,1}^{j_k}$. If $r=2$ and $s=1$, i. e. $j_k$ of $2\times 1$ type, the priority parameters are $p_k:=p_1^k$ and $1-p_k:=p_2^k$, while $A_{j_k}=\left(1,1\right)$.

From rules (A), (B) and (C) (this last one if necessary), for a junction $j_k$ of $r\times s$ type, with initial datum ${\rho }_0^k$ and the flux function (2), the solution ${\widehat{\rho }}^k$ to the RP at $j_k$ is as follows. Consider the function $\omega :\left[0,{\rho }_{max}\right]\to \left[0,{\rho }_{max}\right]$ that, for each arc $I_i^k$, $i=1,\dots ,r+s,$ satisfies the properties:

• $f\left(\omega \left({\rho }_i^k\right)\right)=f\left({\rho }_i^k\right)$∀${\rho }_i^k\in \left[0,{\rho }_{max}\right]$;
• $\omega \left({\rho }_i^k\right)\ne {\rho }_i^k$∀${\rho }_i^k\in \left[0,{\rho }_{max}\right]\setminus \left\{\frac{{\rho }_{max}}{2}\right\}$.

Then, for $I_{\phi }^k$, $\phi =1,\dots ,r$:

• ${\widehat{\rho }}_{\phi }^k\in \left\{{\rho }_{\phi ,0}^k\right\}\cup \left]\omega \left({\rho }_{\phi ,0}^k\right),{\rho }_{max}\right]$ if $0\le {\rho }_{\phi ,0}^k\le \frac{{\rho }_{max}}{2}$;
• ${\widehat{\rho }}_{\phi }^k\in \left[\frac{{\rho }_{max}}{2},{\rho }_{max}\right]$ if $\frac{{\rho }_{max}}{2}\le {\rho }_{\phi ,0}^k\le {\rho }_{max}$.

For $I_{\psi }^k$, $\psi =r+1,\dots ,r+s$:

• ${\widehat{\rho }}_{\psi }^k\in \left[0,\frac{{\rho }_{max}}{2}\right]$ if $0\le {\rho }_{\psi ,0}^k\le \frac{{\rho }_{max}}{2}$;
• ${\widehat{\rho }}_{\psi }^k\in \left\{{\rho }_{\psi ,0}^k\right\}\cup \left[0,\omega \left({\rho }_{\psi ,0}^k\right)\right[$ if $\frac{{\rho }_{max}}{2}\le {\rho }_{\psi ,0}^k\le {\rho }_{max}$.

From ${\widehat{\rho }}_{\phi }^k$ and ${\widehat{\rho }}_{\psi }^k$, we get the maximal flux values on $I_{\phi }^k$, $\phi =1,\dots ,r$ and $I_{\psi }^k$, $\psi =r+1,\dots ,r+s$, i. e.:

$${\gamma }_{\phi }^{k,max}=\left\{\begin{array}{cc}f\left({\rho }_{\phi ,0}^k\right),\hfill & \mathrm{if}0\le {\rho }_{\phi ,0}^k\le \frac{{\rho }_{max}}{2},\\ f\left(\frac{{\rho }_{max}}{2}\right),\hfill & \mathrm{if}\frac{{\rho }_{max}}{2}\le {\rho }_{\phi ,0}^k\le {\rho }_{max},\end{array}\right.$$

(3)

$${\gamma }_{\psi }^{k,max}=\left\{\begin{array}{cc}f\left(\frac{{\rho }_{max}}{2}\right),\hfill & \mathrm{if}0\le {\rho }_{\psi ,0}^k\le \frac{{\rho }_{max}}{2},\\ f\left({\rho }_{\psi ,0}^k\right),\hfill & \mathrm{if}\frac{{\rho }_{max}}{2}\le {\rho }_{\psi ,0}^k\le {\rho }_{max}.\end{array}\right.$$

(4)

Remark 2.

Notice that such an approach allows to define solutions to the Cauchy problem on the network $\left(\mathcal{I},\mathcal{J}\right)$ via a wave-front tracking algorithm (see [17,18]). Refer to Appendix A for details about the construction of solutions.

3. Energy Optimization

Consider an energy network $\left(\mathcal{I},\mathcal{J}\right)$ as described in Section 2. For the optimization of the network performances, we define the cost functional:

$$\overline{\mathcal{E}}\left(t\right):={\left(\sum _{n=1}^N\underset{I_n}{\int }{v}_n\left({\rho }_n\left(t,x\right)\right)dx\right)}^2,$$

that represents a term proportional to the kinetic energy on the whole network.

Considering bounded ${\rho }_n\left(t,x\right)$, $n=1,\dots ,N$, the aim is to maximize $\overline{\mathcal{E}}\left(t\right)$ w.r.t. the distribution coefficients of matrices $A_{j_k}$∀$j_k\in$ $\mathcal{J}$.

As the solution of such optimization control problem involves space-time variables and the optimization itself refers only to $1\times 2$ nodes for the energy hub described in Section 4, we consider an approach defined by the steps:

• Consider a node $j_k\in$ $\mathcal{J}$ of $1\times 2$ type (one incoming arc, $I_1^k$, and two outgoing arcs, $I_2^k$ and $I_3^k$) for which only one distribution coefficient ${\alpha }_k$ is considered, see Remark 1. Assuming an initial datum $\left({\rho }_{1,0}^k,{\rho }_{2,0}^k,{\rho }_{3,0}^k\right)$ at $j_k$, fix the local cost functional:

  $${\mathcal{E}}_{j_k}\left(t\right):={\left(\sum _{m=1}^3\underset{I_m^k}{\int }{v}_m\left({\rho }_m^k\left(t,x\right)\right)dx\right)}^2.$$
• For a time horizon $\left[0,T\right]$, with T quite big, assume the traffic distribution coefficient ${\alpha }_k$ as control, and maximize ${\mathcal{E}}_{j_k}\left(T\right)$ w.r.t. ${\alpha }_k$.
• Construct the optimal solution of the overall network by localization, i.e by using the single optimization solutions at each node $j_k\in$ $\mathcal{J}$ of $1\times 2$ type.

For step 2, assume the conditions:

• ${\mathbf{T}}_1:{\gamma }_3^{k,max}\le \frac{{\gamma }_1^{k,max}}{2}<{\gamma }_1^{k,max}\le {\gamma }_2^{k,max}$;
• ${\mathbf{T}}_2:{\gamma }_2^{k,max}<\frac{{\gamma }_1^{k,max}}{2}<{\gamma }_1^{k,max}\le {\gamma }_3^{k,max}$;
• ${\mathbf{T}}_3:{\gamma }_2^{k,max}<{\gamma }_3^{k,max}<{\gamma }_1^{k,max}$;
• ${\mathbf{T}}_{3A}:{\gamma }_1^{k,max}-{\gamma }_3^{k,max}\ge {\gamma }_2^{k,max}$;
• ${\mathbf{T}}_{3B}:{\gamma }_1^{k,max}-{\gamma }_3^{k,max}<{\gamma }_2^{k,max}\le \frac{{\gamma }_1^{k,max}}{2}$;
• ${\mathbf{T}}_4:{\gamma }_3^{k,max}<{\gamma }_2^{k,max}<{\gamma }_1^{k,max}$;
• ${\mathbf{T}}_5:\frac{{\gamma }_1^{k,max}}{2}\le {\gamma }_1^{k,max}-{\gamma }_3^{k,max}<{\gamma }_2^{k,max}$,

and define:

$$g_{uv}^{k,max}:=\frac{{\gamma }_u^{k,max}}{{\gamma }_v^{k,max}}.$$

We get the following:

Theorem 1.

Consider a node $j_k\in$ $\mathcal{J}$ of $1\times 2$ type, and assume T sufficiently big. Then, ${\mathcal{E}}_{j_k}\left(T\right)$ is maximized for the value ${\alpha }_k^{opt}$ (for some cases, the optimal control does not exist and it is approximated through a positive and small constant ε):

$${\alpha }_k^{opt}=\left\{\begin{array}{cc}1-g_{31}^{k,max}+\epsilon ,\hfill & \mathrm{if}{\mathbf{T}}_{1}holds;\\ g_{21}^{k,max},\hfill & \mathrm{if}{\mathbf{T}}_{2}issatisfied;\\ {\left(1+g_{32}^{k,max}\right)}^{-1},\hfill & \mathrm{if}{\mathbf{T}}_3\mathrm{and}{\mathbf{T}}_{3A}bothhold;\\ g_{21}^{k,max}-\epsilon ,\hfill & \mathrm{if}{\mathbf{T}}_3\mathrm{and}{\mathbf{T}}_{3B}arebothtrue;\\ {\left(1+g_{32}^{k,max}\right)}^{-1}+\epsilon ,\hfill & \mathrm{if}{\mathbf{T}}_{3A}\mathrm{and}{\mathbf{T}}_{4}bothhold;\\ 1-g_{31}^{k,max}-\epsilon ,\hfill & \mathrm{if}{\mathbf{T}}_{5}issatisfied;\\ \frac{1}{2},\hfill & otherwise.\end{array}\right.$$

Proof. 

Assume that $v_{max}={\rho }_{max}=1$. Other cases for which $v_{max}\ne 1$ and/or ${\rho }_{max}\ne 1$ lead to the same results. Fix a node $j_k\in$ $\mathcal{J}$ of $1\times 2$ type and T quite big. Considering the solution to the RP at $j_k$, ${\mathcal{E}}_{j_k}\left(T\right)$ writes as:

$${\mathcal{E}}_{j_k}\left(T\right)=\left(3-s_1^k\sqrt{1-4{\widehat{\gamma }}_1^k}-s_2^k\sqrt{1-4{\alpha }_k{\widehat{\gamma }}_1^k}-s_3^k\sqrt{1-4\left(1-{\alpha }_k\right){\widehat{\gamma }}_1^k}\right),$$

(5)

where coefficients $s_1^k$ and $s_{\psi }^k,$ $\psi =2,3,$ are:

$$s_1^k=\left\{\begin{array}{cc}+1,\hfill & \mathrm{if}{\rho }_{1,0}^k\ge \frac{1}{2},\\ \hfill & \mathrm{or}{\rho }_{1,0}^k<\frac{1}{2}\mathrm{and}{\gamma }_1^{k,max}>min\left\{\frac{{\gamma }_2^{k,max}}{{\alpha }_k},\frac{{\gamma }_3^{k,max}}{1-{\alpha }_k}\right\},\\ -1,\hfill & \mathrm{if}{\rho }_{1,0}^k<\frac{1}{2}\mathrm{and}{\gamma }_1^{k,max}\le min\left\{\frac{{\gamma }_2^{k,max}}{{\alpha }_k},\frac{{\gamma }_3^{k,max}}{1-{\alpha }_k}\right\},\end{array}\right.$$

$$s_{\psi }^k=\left\{\begin{array}{cc}+1,\hfill & \mathrm{if}{\rho }_{\psi ,0}^k>\frac{1}{2}\mathrm{and}\frac{{\gamma }_{\psi }^{k,max}}{{\alpha }_{\psi }}\le min\left\{{\gamma }_1^{k,max},\frac{{\gamma }_{{\psi }^{\prime }}^{k,max}}{{\alpha }_{{\psi }^{\prime }}}\right\},{\psi }^{\prime }\ne \psi ,\\ -1,\hfill & \mathrm{if}{\rho }_{\psi ,0}^k\le \frac{1}{2},\\ \hfill & \mathrm{or}{\rho }_{\psi ,0}^k>\frac{1}{2}\mathrm{and}\frac{{\gamma }_{\psi }^{k,max}}{{\alpha }_{\psi }}>min\left\{{\gamma }_1^{k,max},\frac{{\gamma }_{{\psi }^{\prime }}^{k,max}}{{\alpha }_{{\psi }^{\prime }}}\right\},{\psi }^{\prime }\ne \psi ,\end{array}\right.$$

with

$${\alpha }_{\psi }=\left\{\begin{array}{cc}{\alpha }_k,\hfill & \mathrm{if}\psi =2,\\ 1-{\alpha }_k,\hfill & \mathrm{if}\psi =3.\end{array}\right.$$

For simplicity, from now on we drop the dependence on $j_k$ and T from $\mathcal{E}$. As the solution to the RP at $j_k$ depends on ${\alpha }_k$, we have various cases. Here, for seek of brevity, we consider only two of them, as the proof for other cases is similar.

Assume ${\gamma }_1^{k,max}<{\gamma }_3^{k,max}<{\gamma }_2^{k,max}$. In this case, ${\widehat{\gamma }}_1^k={\gamma }_1^{k,max}$, $s_1^k=s_2^k=s_3^k=-1$ and we have to maximize:

$$\mathcal{E}={\left(3+\sqrt{1-4{\gamma }_1^{k,max}}+\sqrt{1-4{\alpha }_k{\gamma }_1^{k,max}}+\sqrt{1-4\left(1-{\alpha }_k\right){\gamma }_1^{k,max}}\right)}^2,$$

defined for ${\alpha }_k\in \left]0,1\right[$. We get that:

$$\frac{\partial \mathcal{E}}{\partial {\alpha }_k}=4\mathcal{E}{\gamma }_1^{k,max}\left[\Phi \left(1-{\alpha }_k\right)-\Phi \left({\alpha }_k\right)\right],$$

(6)

where $\Phi \left({\alpha }_k\right)$ is:

$$\Phi \left({\alpha }_k\right):=\frac{1}{\sqrt{1-4{\alpha }_k{\gamma }_1^{k,max}}}.$$

Since

$$\frac{\partial \mathcal{E}}{\partial {\alpha }_k}\ge 0\iff {\alpha }_k\le \frac{1}{2},$$

then ${\alpha }_k^{opt}={\displaystyle \frac{1}{2}}$.

Consider now the case ${\gamma }_3^{k,max}<{\gamma }_1^{k,max}<{\gamma }_2^{k,max}$. We have that:

• if $0<{\alpha }_k\le 1-g_{31}^{k,max}$, ${\widehat{\gamma }}_1^k={\displaystyle \frac{{\gamma }_3^{k,max}}{1-{\alpha }_k}}$ and $s_1^k=s_3^k=+1$, $s_2^k=-1$;
• if $1-g_{31}^{k,max}<{\alpha }_k<1$, ${\widehat{\gamma }}_1^k={\gamma }_1^{k,max}$ and $s_1^k=s_2^k=s_3^k=-1$.

Hence, (5) becomes:

$$\mathcal{E}=\left\{\begin{array}{cc}{\mathcal{E}}_1,\hfill & 0<{\alpha }_k\le 1-g_{31}^{k,max},\\ {\mathcal{E}}_2,\hfill & 1-g_{31}^{k,max}<{\alpha }_k<1,\end{array}\right.$$

where:

$${\mathcal{E}}_1={\left(3-\sqrt{1-4{\displaystyle \frac{{\gamma }_3^{k,max}}{1-{\alpha }_k}}}+\sqrt{1-4{\alpha }_k{\displaystyle \frac{{\gamma }_3^{k,max}}{1-{\alpha }_k}}}-\sqrt{1-4{\gamma }_3^{k,max}}\right)}^2,$$

$${\mathcal{E}}_2={\left(3+\sqrt{1-4{\gamma }_1^{k,max}}+\sqrt{1-4{\alpha }_k{\gamma }_1^{k,max}}+\sqrt{1-4\left(1-{\alpha }_k\right){\gamma }_1^{k,max}}\right)}^2.$$

If $0<{\alpha }_k\le 1-g_{31}^{k,max}$, then

$$\frac{\partial \mathcal{E}}{\partial {\alpha }_k}=\frac{4\mathcal{E}{\gamma }_3^{k,max}}{{\left(1-{\alpha }_k\right)}^2}\left[\Psi \left(\frac{1}{1-{\alpha }_k}\right)-\Psi \left(\frac{{\alpha }_k}{1-{\alpha }_k}\right)\right],$$

where:

$$\Psi \left({\alpha }_k\right):=\frac{1}{\sqrt{1-4{\alpha }_k{\gamma }_3^{k,max}}}.$$

It is possible to verify that $\mathcal{E}$ is an increasing function. If $1-g_{31}^{k,max}<{\alpha }_k<1$, ${\displaystyle \frac{\partial \mathcal{E}}{\partial {\alpha }_k}}$ is the expression (6). Hence, we conclude that:

• if $1-g_{31}^{k,max}<{\displaystyle \frac{1}{2}},$ $\mathcal{E}$ is optimized for ${\alpha }_k^{opt}={\displaystyle \frac{1}{2}};$
• if $1-g_{31}^{k,max}\ge {\displaystyle \frac{1}{2}},$ $\mathcal{E}$ has not an optimal value, hence ${\alpha }_k^{opt}$ is chosen as ${\alpha }_k^{opt}=1-g_{31}^{k,max}+\epsilon$, where $\epsilon$ is a positive and small constant.

□

4. 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 here considered. This is a relevant problem in modern energy systems, where the increasing inter-dependencies between heterogeneous energy infrastructures is introducing new and more complex vulnerabilities, requiring effective modeling and optimization tools aimed at improving the accuracy and the robustness of coordinated control actions. In this context, the large scale deployment for the energy hub could have a strategic role, since it allows improving the energy networks 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, Canada, for the supply of commercial load, has been analyzed. More details on this system, which is schematically depicted in Figure 1, can be found in [11].

From Figure 1 it is worth noting as the input power flows are the electricity, $P_E$, and the natural gas, $P_G$, while the output power flows are electricity, $L_E$, and heat, $L_H$. Notice that $P_E$ is splitted into ${\beta }_1{P}_E$ and ${\beta }_2{P}_E$, while $P_G$ is divided into ${\beta }_3{P}_G$ and ${\beta }_4{P}_G$. Precisely, ${\beta }_1{P}_E$ and ${\beta }_2{P}_E$ are, respectively, inputs of the Hydrogen Production Plant (HPP) and of the Transformer (T); ${\beta }_3{P}_G$ is input of the Combined Heat Power (CHP), while ${\beta }_4{P}_G$ is input of the Furnace (F). The system works as follows:

• The subsystem HPP, characterized by electric-hydrogen and heat-hydrogen efficiences ${\eta }_1^{HPP}$ and ${\eta }_2^{HPP}$, transforms a part of the electricity, ${\beta }_1{P}_E{\eta }_1^{HPP}$, into hydrogen that feeds the Fuel Cell (FC) and another part, ${\beta }_1{P}_E{\eta }_2^{HPP}$, into heat.
• Using the hydrogen-electricity and hydrogen-heat efficiencies ${\eta }_1^{FC}$ and ${\eta }_2^{FC}$, the subsystem FC transform a part of the hydrogen, ${\beta }_1{P}_E{\eta }_1^{HPP}{\eta }_1^{FC}$, into electricity and another part, ${\beta }_1{P}_E{\eta }_1^{HPP}{\eta }_2^{FC}$, into heat.
• The subsystem T, due to its efficiency ${\eta }^T$, has the electricity power flow ${\beta }_2{P}_E{\eta }^T$ as output.
• The subsystem CHP, considering the gas-electric and gas-heat efficiencies ${\eta }_1^{CHP}$ and ${\eta }_2^{CHP}$, transforms a part of natural gas, ${\beta }_3{P}_G{\eta }_1^{CHP}$, into electricity and another part, ${\beta }_3{P}_G{\eta }_2^{CHP}$, into heat.
• Finally, the subsystem F, characterized by its efficiency ${\eta }^F$, has output ${\beta }_4{P}_G{\eta }^F$, that is heat.

As for the outputs of the energy hub, we simply get that:

$$L_E={\beta }_1{P}_E{\eta }_1^{HPP}{\eta }_1^{FC}+{\beta }_2{P}_E{\eta }^T+{\beta }_3{P}_G{\eta }_1^{CHP},$$

(7)

$$L_H={\beta }_1{P}_E{\eta }_1^{HPP}{\eta }_2^{FC}+{\beta }_1{P}_E{\eta }_2^{HPP}+{\beta }_3{P}_G{\eta }_2^{CHP}+{\beta }_4{P}_G{\eta }^F.$$

(8)

Notice that (7) and (8) are obtained considering some losses that depend on coefficients ${\eta }_1^{HPP}$, ${\eta }_2^{HPP}$, ${\eta }_1^{FC}$, ${\eta }_2^{FC}$, ${\eta }^T$, ${\eta }_1^{CHP}$, ${\eta }_2^{CHP}$ and ${\eta }^F$. Such parameters are usually fixed and, for the topology described in Figure 1, they are the following: ${\eta }_1^{HPP}=0.7,$ ${\eta }_2^{HPP}=0.2,$ ${\eta }_1^{FC}=0.55,$ ${\eta }_2^{FC}=0.4,$ ${\eta }^T=0.98,$ ${\eta }_1^{CHP}=0.45,$ ${\eta }_2^{CHP}=0.35,$ ${\eta }^F=0.48$. Moreover, note that ${\beta }_2=1-{\beta }_1,$ ${\beta }_4=1-{\beta }_3,$ $0<{\beta }_1<1,$ $0<{\beta }_2<1$, 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 ${\beta }_1$ and ${\beta }_2$ (of 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 choose 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 “output” different from $L_E$ and $L_H$.

4.2. Numerical Results

This section presents the simulations for the energy hub, represented in Figure 2 as a couple $\left(\mathcal{I},\mathcal{J}\right)$, with $\mathcal{I}={\left\{I_n\right\}}_{n=1,..,23}$ and $\mathcal{J}={\left\{j_k\right\}}_{k=1,..,11}$, see Section 2. The features of the network in Figure 2 are as follows. External arcs: $I_1$, $I_5$, $I_7$, $I_{12}$, $I_{13}$, $I_{14}$, $I_{17}$, $I_{19}$, $I_{21}$, $I_{23}$. Inner arcs: $I_2$, $I_3$, $I_4$, $I_6$, $I_8$, $I_9$, $I_{10}$, $I_{11}$, $I_{15}$, $I_{16}$, $I_{18}$, $I_{20}$, $I_{22}$. Nodes of: $1\times 2$ type, $j_1$, $j_4$, $j_7$, $j_9$; $2\times 1$ type: $j_5$, $j_6$, $j_{10}$, $j_{11}$; $1\times 3$ type: $j_2$, $j_3$, $j_8$.

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 with arc $I_m$, named simply m.

Notice that arcs 1 and 14 are the inputs of the systems, while outputs are arcs 13 and 23, indicated, respectively, by $OUT$ 2 and $OUT$ 1.

The performances of the network are evaluated through the cost functional $\overline{\mathcal{E}}\left(t\right)$, 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, distribution matrices $A_J$, $J\in \left\{2,3,4,8,9\right\}$ assume the form:

$$A_2=\left(\begin{array}{c}c{\alpha }_{4,2}\hfill \\ {\alpha }_{5,2}\hfill \\ {\alpha }_{6,2}\hfill \end{array}\right)=\left(\begin{array}{c}c0.7\hfill \\ 0.1\hfill \\ 0.2\hfill \end{array}\right),$$

$$A_3=\left(\begin{array}{c}c{\alpha }_{7,4}\hfill \\ {\alpha }_{8,4}\hfill \\ {\alpha }_{9,4}\hfill \end{array}\right)=\left(\begin{array}{c}c0.05\hfill \\ 0.55\hfill \\ 0.4\hfill \end{array}\right),$$

$$A_4=\left(\begin{array}{c}c{\alpha }_{11,3}\hfill \\ {\alpha }_{12,3}\hfill \end{array}\right)=\left(\begin{array}{c}c0.98\hfill \\ 0.02\hfill \end{array}\right),$$

$$A_8=\left(\begin{array}{c}c{\alpha }_{17,15}\hfill \\ {\alpha }_{18,15}\hfill \\ {\alpha }_{19,15}\hfill \end{array}\right)=\left(\begin{array}{c}c0.2\hfill \\ 0.45\hfill \\ 0.35\hfill \end{array}\right),$$

$$A_9=\left(\begin{array}{c}c{\alpha }_{20,16}\hfill \\ {\alpha }_{21,16}\hfill \end{array}\right)=\left(\begin{array}{c}c0.68\hfill \\ 0.32\hfill \end{array}\right).$$

Moreover, always referring to measures on the real hub, priority parameters are all chosen $0.5$ for incoming arcs of nodes 5, 6, 10 and 11. Hence, no control is considered for junctions 2, 3, 4, 5, 6, 8, 9, 10 and 11, and the optimization of $\overline{\mathcal{E}}\left(t\right)$ deals only with the distribution coefficients at nodes 1 and 7 of $1\times 2$ type.

As for the numerical construction of $\overline{\mathcal{E}}\left(t\right)$, it is necessary a suitable approximation for densities ${\rho }_i\left(t,x\right)$, $i=1,\dots ,23$, whose evolution is ruled by equation (1).

In this paper, we apply the Godunov scheme (see [32,33,34]),using a numerical grid with constant space and time sizes, $\Delta x=0.0125$ and $\Delta t=0.5\Delta x$, respectively (see subSection 4.4 for details about the computational cost). The network of Figure 2 is simulated in such conditions: time interval of simulation: $\left[0,T\right]$, with $T=150min$; empty arcs when the simulation starts ($t=0$); boundary data of Dirichlet type, equal to $0.3$ for arcs 1 and 14 while, for arcs 5, 7, 12, 13, 17, 19, 21 and 23, we choose a Dirichlet boundary data equal to $0.9$.

Notice that typical maximal values for inputs of our hub are 15 MWh and 20 MWh for arcs 1 and 14, respectively. In our case, associating at ${\rho }_{max}=1$ the quantity 15 MWh, we simply get that boundary data $0.3$ and $0.9$ correspond to $4.5$ MWh and $13.5$ MWh.

Two different choices of the distribution parameters are assumed for nodes 1 and 7:

• optimal case: parameters that optimize locally the asymptotic behaviour of $\overline{\mathcal{E}}\left(t\right)$, i.e. distribution coefficients that refer to Theorem 1 for junctions 1 and 7. Such type of simulation is useful to test the global performance, starting from analytical results that consider only a part of nodes of the network.
• random case: parameters at nodes 1 and 7 are chosen in a random way at $t=0$ and then are kept constant in $\left[0,T\right]$. A random simulation allows comparisons with network performances obtained via local optimal distribution coefficients.

4.3. Results Discussion

In Figure 3, Figure 4, Figure 5 and Figure 6 we show some simulation results for the energy hub. More precisely, the values of $\overline{\mathcal{E}}\left(t\right)$, computed on the whole network, is represented as function of the time.

In particular, the behavior of $\overline{\mathcal{E}}\left(t\right)$ in the optimal case is compared with the ones obtained via ten different simulation random studies. As for these last cases, Figure 3 and Figure 4 show the first five ones, while Figure 5 and Figure 6 the remaining other ones.

The optimization algorithm for $1\times 2$ nodes, which is of local type, can be applied to the complex topology of the energy hub without compromising the possibility of a global optimization. Such a situation is evident in the optimal case for $\overline{\mathcal{E}}\left(t\right)$, that is compared to the behaviors in ten different random cases. In Figure 3, Figure 4, Figure 5 and Figure 6, it is shown that the optimal case is always higher than random cases.

Indeed, 100 random cases have been simulated and compared with the optimal behavior for $\overline{\mathcal{E}}\left(t\right)$. Table 1 reports the value of $\overline{\mathcal{E}}\left(t\right)$ in the optimal configuration at $T=150$, i.e. OPTconf T, with the average value (RAND T) of random simulations at $T=150$.

Notice that OPTconf T is, as expected, higher than RAND T, namely the global optimization of local type has a strong robustness. Such result is also represented in Figure 7, where a histogram reports the values of $\overline{\mathcal{E}}\left(t\right)$ at the final instant T for the random simulation.

Finally, notice that the simulated system has always 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) are, respectively, $0.1218$ ($36.54$ kWh) and $0.180075$ ($54.0225$ kWh).

4.4. Computational Cost

This subsection presents some details about the computational cost for the simulation of the presented energy system. Focus on a network represented by the couple $\left(\mathcal{I},\mathcal{J}\right)$, with $\mathcal{I}={\left\{I_n\right\}}_{n=1,..,N}$ and $\mathcal{J}={\left\{j_k\right\}}_{k=1,..,J}$.

In order to find a suitable numerical approximation in $\left[0,T\right]$ for the density functions ${\rho }_n\left(t,x\right)$, $n=1,..,N$, on arcs and update of boundary data at nodes, assume that each arc $I_n$, $n=1,..,N$, has length $L_n$; space and time grid sizes are, respectively, $\Delta x_n$ and $\Delta t_n$. Using the Godunov method, the computational cost depends on ${\displaystyle \sum _{n=1}^N}{\displaystyle \frac{L_n}{\Delta x_n}}$ and ${\displaystyle \sum _{k=1}^J}{\displaystyle \frac{T}{\Delta t_k}}$ for densities and boundary data, respectively.

For simplicity, for the simulation of the described energy hub we consider a constant space grid size $\Delta x$, assume $\left(\Delta x_n,\Delta t_n\right)=\left(\Delta x,0.5\Delta x\right)$∀$n=1,\dots ,N$, and compute the CPU times (measured in seconds and calculated by an Intel(R) Core (TM) i7-3630 QM CPU @2.40 GHz, RAM 8 GB) and convergence errors. The obtained results are in Table 2.

From the previous table, we simply get that the CPU time increases of about 0.30 seconds when $\Delta x$ decreases. As for the convergence error, it almost remains the same for different values of $\Delta x$. Some further studies, as well as different numerical approaches for conservation laws on networks, are carefully analyzed in [18,23] and [30].

5. 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 have been 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. To address this issue the application of conservation laws and the definition of cost functional, which represents a term proportional to the kinetic energy of the system, have been proposed for dealing with the 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 have been assessed on a realistic case study, dealing with the solution of the optimal energy flow management problem for a complex energy hub designed in Waterloo, Canada. The obtained results demonstrated 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 argue that the obtained solution is globally optimal, and robust against distributed parameters variations, which can be considered further, and certainly relevant, benefits. A rigorous theoretical justification of these intuitions are currently under investigation by the authors.