Health-Aware Economic MPC for Operational Management of Flow-Based Networks Using Bayesian Networks

: This paper presents a health-aware economic Model Predictive Control (EMPC) approach for the Prognostics and Health Management (PHM) of generalized ﬂow-based networks. The proposed approach consists of the integration of the network reliability model obtained from a Bayesian network in the control model. The controller is then able to optimally manage the supply taking into consideration the distribution of the control effort, to extend the life of the actuators by delaying the network reliability decay as much as possible. It also considers an optimal inventory replenishment policy based on a desired risk acceptability level, leading to the availability of safety stocks for unexpected excess demand in networks. The proposed implementation is illustrated with a real case study corresponding to an aggregate model of the Drinking Water transport Network of Barcelona.


Introduction
The management of flow-based networks is an interesting and increasing research subject owing to the fact that there are uncertain and complex systems with an important economic, environmental, and social impact. The use of control strategies that take into account the system health by monitoring their components reliability is necessary to guarantee a quality service, while minimizing the fault appearance and reducing the operational costs. However, the random behaviour of the network demands and the variability on the prices of the electricity (which directly influences the actuators operation costs) are critical characteristics when trying to control this type of system in real-time.
To deal with this problem, several approaches have been proposed from the research community. Most of them employ a Model Predictive Control (MPC) approach, due to the fact that it suitably fits with the need to take into account the demand forecasts to sufficiently fill the different reservoirs on time. This control strategy is based on solving a finite time-horizon optimization problem, given some future predictions, in order to minimize some considered control objectives and satisfy the set of constraints, including the system model and physical/operational limitations.
However, for managing water systems, MPC is not implemented in a classical way due to there being no reference volume to be tracked [1]. Instead of this, the standard MPC forces the system to follow a setpoint, but it does not guarantee the economic efficiency of the system evolution toward this setpoint, and this is just the general aim in the operation of several industrial processes: the reduction of costs associated to the energy consumption.
For this purpose, there exists a systematic method to optimize the cost of a system management: the Economic MPC (EMPC) [2]. However, there is a well-known problem related with this strategy, which is the aim of obtaining a family of optimal setpoints considering economic efficiency rather than aiming the controlled system to reach a certain setpoint [2]. Nevertheless, different approaches were suggested to enhance this concept, such as in Karimi Pour et al. [3], where they implemented an LPV-MPC controller based on a single-layer economic optimization problem with dynamic constraints, including predetermined safety levels to deal with the demand uncertainty. In Grosso et al. [4], where they also deal with an MPC based on a single-layer economic optimization problem with dynamic constraints to cope with the components degradation awareness and safety stock availability to satisfy non-stationary flow demands. In Karimi Pour et al. [5], where they proposed a health-aware LPV-MPC by using a chance-constraints approach of the reliability model. In Gokdere et al. [6], where the actuator lifespan is included as an additional goal parameterized in the linear quadratic optimal controller to reduce the cost of maintenance. Or in Pereira et al. [7], where it introduces an MPC strategy based on distributing the loads among redundant actuators, also with the use of constraints to ensure that the accumulated degradation of the actuators will not achieve an unsafe level within a prediction horizon.
Despite this, the reliability is the system's (or component) ability to perform its intended functions. In a Drinking Water transport Network (DWN) case, some conditions that would mainly influence it are: the capacity and the quality of the water accessible at the sources, as well as the pump/pipe failure rates [8,9]. Moreover, it is characterised according to the interdependence topology based on the network layout, given by the specific combination of the involved components. Subsequently, the system reliability leads to a nonlinear mathematical model when considering the control input.
In the literature, how the control affects the system health is considered by adding a damage index in the goals of the optimization-based control and by adjusting the trade-off between optimal operation and damage mitigation by weight tuning [10], or by defining constraints involving the actuators' reliabilities [8]. However, considering the reliability at the components level and not at the system level is the main drawback of these previous methods, as it implies a high computational cost. In addition, the existing gradient-based numerical algorithms do not certify that the obtained solution corresponds to the global one, because of the non-convexity of the associated optimization problem. This problem could be overcome by transforming the nonlinear optimization problem into a quadratic problem through a linearisation method.
The main contribution of this paper is to integrate a Bayesian network reliability model of a generalized flow-based network into the common economic MPC used in these cases and evaluate it together with the optimal inventory replenishment conditions. Thus, the system reliability obtained from the Bayesian inference is an online event-oriented performance criterion that measures the probability that all demands will be fully met within a given horizon from the available stock in time, under normal and adverse conditions. In addition, the method, instead of considering the reliability of all components, basically relies on the degradation of the actuators, which are the network components most exposed to stress when applying the control actions.
From this point on, the design of the proposed approach will be contextualized on a DWN, but it is important to note that the approach is generalized, and can be implemented in an analogous manner for other flow-based network cases as, e.g., power networks. The paper is structured as follows: Section 2 presents the EMPC of DWN. Section 3 describes the inclusion of the health-aware objective in the EMPC. Section 4 illustrates the proposed approach in an aggregate model of the Barcelona DWN. Finally, Section 5 draws the main conclusions and suggests future research topics.
Notation: R, R + , R n , R m×n represent the set of real numbers, of non-negative real numbers, of column real vectors of length n and of m by n real matrices, respectively. Similarly, I + presents the set of non-negative integer numbers including zero. The operator ⊕ denotes the direct matrix sum. . is the matrix spectral norm and ||.|| 2 is the squared 2-norm. The superscript represents the transpose of a vector/matrix and operators <, ≤, =, >, ≥ represents element-wise vector relations.

EMPC for DWN
As mentioned above, the proposed approach is based on an Economic MPC. The following subsections introduce the model of a DWN as the control-oriented one, and the designed optimization problem for the EMPC.

Control-Oriented Model
The control-oriented model to implement the MPC is expressed in discrete-time statespace system considering a flow-based modelling approach. The dynamics of the storage devices for all time instant k ∈ Z ≥0 lead to the following state equation where x ∈ R n x are the volumes and n x the number of storage devices on the network. The vector u ∈ R n u represent the control inputs associated to the flow rates through the actuators of the network, being n u the total number of them. The vector d ∈ R n d represents the disturbances corresponding to the consumer water demands, being n d the total number of them. A ∈ R n x ×n x , B ∈ R n x ×n u and B d ∈ R n x ×n d are the system time-invariant matrices that depend on the network topology. Furthermore, the system is subject to some constraints. First of all, it is subject to the flow-mass balance relations in the nodes (being n n the total number of them), leading to static equation for each node formulated in a matrix form as follows where E u ∈ R n n ×n u and E d ∈ R n n ×n d are the time-invariant matrices that depend on the network junctions connections. Besides, a DWN is subject to the physical inputs and states constraints, provided by convex and closed polytopic sets defined as: where G ∈ R n x ×n x , g ∈ R n x , H ∈ R n u ×n u , and h ∈ R n u are matrices collecting the system constraints.
Concerning the operation of the considered flow-based networks, it is assumed that the demands in d(k) and the states in x(k) are measurable at each time instant k ∈ Z ≥0 ; while the pair (A, B) is stabilizable.

Optimization Problem Formulation
As the MPC requires some criteria to obtain the control actions, an optimization problem must be defined. Then, the control aim is to minimize a convex stage cost function J : Z ≥0 × X × U −→ R ≥0 , which could include any functional relation with the system operation. Therefore, the control goal can be expressed as a convex multi-objective cost function to minimize. In this case, the objective terms considered to manage the flow-based network are: • Economic objective: the economical costs that involve the flow transport while providing the demanded volume should be minimized. This cost function corresponds to: where α(k) (α 1 + α 2 (k)) ∈ R n u is the price per volume unit, including the fixed costs α 1 ∈ R n u related to the supply flow and the variable costs α 2 (k) ∈ R n u related to the time-varying electricity associated cost. W e is a diagonal positive definite matrix that is used as a weight to prioritize the terms in the complete objective function.

•
Safety objective: the storage devices should guarantee some safety volume to deal with unexpected variations in the demand. This goal can be formulated as follows: where x s indicates the storage safety volume. However, this piecewise linear formulation can be avoided by considering that the safety cost function can be expressed through a soft constraint by using a slack variable ξ like: and also being introduced as an objective term to retain feasibility of the optimization problem: where W s is a diagonal positive definite matrix that is used as a weight to prioritize the terms in the complete objective function. • Smoothness objective: to avoid overloads in the pipes, and preserve the network components lifetime, the actuators are managed based on smooth control actions. To achieve this smoothing effect, the variation of the control actions among two consecutive time instants is penalized as follows: where ∆u(k) u(k) − u(k − 1), and W ∆u is a diagonal positive definite matrix that is used as a weight to prioritize the terms in the complete objective function.
Thus, the following multi-objective cost function to be minimized in the prediction horizon N p can be formulated as: Finally, to assure the feasibility of the obtained control actions, the cost function must be subject to the system constraints introduced in Section 2.1, as well as to the one related with the safety objective. So, at each time instant, the following optimization problem is solved online: min subject to: The sequences of optimal control actions u [1,Np ] , and safety levels ξ * (k) = {ξ(i|k)} i∈Z [1,Np ] are computed on-line, solving this optimization problem using the receding horizon philosophy [11]. That is, this problem is solved at the current time instant k by using x(0|k) as the initial condition that is obtained from measurements/state estimation at a time instant k. Then, the first value u * (0|k) of the optimal input sequence u * (k) is applied to the system. This procedure is repeated successively for the next time instants.

Reliability-Aware Economic MPC Using Bayesian Network Approach
As discussed in the introduction, the main contribution of this work is to integrate the information about system health in the MPC controller by using a Bayesian model and considering that the actuators are the only components that are affected by the degradation. A continuous evaluation of their states can then be performed to check the overall system reliability.

Bayesian Model
Before introducing the Bayesian Model (BM), it would be convenient to review what Bayesian probability and Bayes' theorem stands for. The first basically consists of an interpretation of the probability, in such a way that instead of frequency or propensity of some phenomenon, probability is interpreted as a degree of belief in an event, like quantifying a reasonable expectation. The second is a probability theory used to deal with Bayesian statistics, which describes a conditional probability for an event based on some data, such as prior information or beliefs about the event. Thus, a Bayesian Model is a very useful statistical model to manage the inference of complex systems with conditional dependencies, and to compute the reliability of any event on them. In addition, they are usually represented by an acyclic directed graph, which easily shows the parental relationships between the relevant system components.
In this work, a BM is aimed to be used for the reliability evaluation of the flow-based network, but only considering the actuators as components with considerable reliability influence. In order to lighten the computation of the proposed approach, and consequently, make it able to compute larger systems, a BM is evaluated for each demand independently. Once the reliabilities of all the demands are obtained, the global system reliability is calculated. Figure 1 shows a graph example of the BM for a single demand, from a simple network example illustrated in Figure 2.  Then, to integrate the BM in the MPC, an equivalent interpretation of the implied inference is required. This reinterpretation starts by finding the minimal paths from all the available sources, and, as indicated, by only taking into account the actuators involved in them. The pure parents are considered to be these actuators, corresponding to the first light blue nodes column of Figure 1. These are observable evidences, required to perform the conditional probabilities for the subsequent children (the next graph columns of Figure 1), which are, respectively:

1.
All paths leading to the relevant demand, whose reliabilities lie on the series multiplication of the actuators' reliabilties involved in each one, which means the direct product of these probabilities: P = ∏ n i=0 P i , implying that all of them must be operative to make the path feasible. In the network example of Figure 1 this is represented with the nodes corresponding to the available paths in the second purple nodes column. The parents for each one are the actuators involved in each path, and they are the same ones whose reliabilities will be multiplied to obtain each path's one; 2.
All sources with access to provide the relevant demand. In this case, the reliability for each one is the parallel multiplication of the reliabilities of those paths that provide supply from the relevant source, which means the complementary product of all the complementaries of these probabilities: P = 1 − ∏ n i=0 1 − P i , implying that one feasible path is enough to make the source available. In the network example of Figure 1, this is represented by the nodes corresponding to the available sources in the third light green nodes column. The parents for each one are the paths supplying from each source, and they are the same ones whose reliabilities will be computed to obtain each source's one; 3.
Finally, there is the relevant demand, whose reliability is also a parallel multiplication of the available sources, implying that one available source is enough to provide the supply. In this case, it corresponds to the last column with a single red node of the graph example of Figure 1, and all the sources are the parents to perform the relevant reliability calculation.
Afterwards, once the reliabilities of all the demands are computed separately, the reliability for the global system is evaluated as well as a series multiplication of all the demands, implying that all of them must be provided to satisfy the system.

Individual Reliability
To characterize the evolution of the reliability with time, many types of distributions can be found in the literature. The most commonly used ones are normal, log-normal, exponential, and Weibull distributions [12]. In this work, the exponential function is considered.
The failure rate λ, an important factor in reliability definition, stands for the fraction of the density of the stochastic lifetime to the remainder function (i.e., conditional probability). Particularly, systems are designed to operate under different load values. According to Jiang and Jardine [12], the load directly impacts the component failure rate. Therefore, in order to present the assessment of system reliability, the relationship between load and failure rate should be taken into account.
In this paper, actuator failure rates are considered to be in function of the load of the applied control input. The following exponential law is the most widely used relationship to characterize the variation of the actuator fault rates with respect to the load where λ 0 i represents the baseline failure rate (nominal failure rate), β i is a constant parameter that depends on the actuator characteristics, and u i (k) is the control action at a time instant k for the i-th actuator, being n the total number of actuators.
Accordingly to the failure rate definition, the reliability of a system or component can be described as: "Reliability is determined as the probability that a system (or a component) will perform their functioning satisfactorily for a certain period of time subject to operating conditions." [13] From a mathematical point of view, reliability R(t) is defined as the probability of the successful operation of a system in the intervening period from time 0 to time t where T is a stochastic variable that describes the time until failure. Furthermore, the unreliability of a system (or a component) F(t) is determined as the probability that the component or system encounters the first failure or has failed one or more times among the time interval 0 to time t.
Considering the system (or component) is always in one of the two possible states (operational or failed), the following relation is satisfied Then, the reliability of a component R i (t), in the useful life period, can be specified at a certain time t by exploiting the exponential function as follows In discrete-time, Equation (15) can be rewritten as where λ i (k) is the failure rate at a time instant k that is acquired from the i-th component under varying load levels u i ; and T s is the sampling time.

Overall System Reliability Modelling
The system state of health can be determined by means of the overall system reliability, R g (k) that can be determined from the elementary reliability components. Thus, R g (k) depends of the system component interconnections that can generally be modelled as a combination of parallel and/or series configurations [14].
To add the Bayesian inference in the MPC, the developed reinterpretation of the Bayesian model, introduced in Section 3.1, must be performed. We start by computing the reliabilities for the multiple paths P j subsystems as R p j (k) = ∏ i∈P j R i (k), P j ⊂ I j = 1, 2, . . . , n p being n p the total number of paths; and P j the subset of the i-indices corresponding to the actuators involved in each path.
Consecutively, the paths are subdivided in new subsystems. First, according to the demand which they provide and then according to the source from where the supply is provided. So, to obtain the reliabilities for each source, the following computation is performed being n s the total number of sources; and S l the subset of the j-indices corresponding to the paths providing from each source. To obtain the reliabilities for each demand being n d the total number of sources; and D h the subset of the l-indices corresponding to the available sources for each demand. Finally, to infer the overall system reliability, the reliabilities of all the demands are evaluated in series, since, as mentioned above, they must all be supplied for the system feasibility

Augmented System
The economic MPC formulation presented in Section 2 must be modified to include the preservation of the actuators' lifetime. This is achieved by adding a new term in the MPC objective function that aims to achieve reliability maximization, and by augmenting the system model according to the reliability model obtained using the Bayesian modelling approach presented above. For this purpose, the reliabilities evolution over time is included in the model in such a way that a conversion that allows computing them in a linear-like form is needed. This conversion is based on applying logarithms, starting from rewriting the individual actuators reliability evolution from Equation (16) as follows This expression allows the system to actualize the parent reliabilites, but, for computation simplicity, the extra objective term should only include the whole system reliability, computed from (20). So, to obtain the demands reliabilities using logarithms to preserve the linear computation convenience, the Equation (19) is expanded while taking into account the unreliability, according to (14) and leading to Then, to relate (21) with (22), the following variable change is introduced which, by applying logarithms, drives to and according to (22), which is equivalent to the following expression with β j = log z j (k) All in all, the structure of the augmented system is the following where the state vector will also be augmented by including the logarithms of the demands' unreliabilities and the logarithms of the actuators' reliability, in order to actualize them properly in each iteration, as follows and the system matrices corresponding to However, to be able to implement the Equation (21) in the augmented system, we should remove the exponential of the failure rate function (12); since it has the control input in its exponent, and it does not comply with the linear computation of (27). To achieve the required linearity, a good approach would be to approximate the exponential law to its corresponding first order Taylor series, as follows As can be appreciated in Figure 3, for small values, the approximation fits quite well. Taking into account that the real values of λ 0 i , β i and u i (k) will lead to small values satisfying the approximation.
Then, in terms of (12), the expression would end up being Therefore, the influence of the failure decay can be interpreted with the resulting terms of (31), which corresponds to a typical linear function: with a constant independent term, which would represent the unavoidable constant degradation over time; and a proportional term as a function of the control input, which would represent the relative degradation to the applied control. Assuming that the constant degradation computation does not provide us with any further advantage, owing to it only would imply the addition of a constant term in (27), which is removed as well to simplify the computation, (31) resulting in

Simple Network Example
This subsection includes an example to validate the method with a previously illustrated network (Figure 2). For this purpose, the network reliability (reliability of the demand to be supplied) is calculated in the three different ways that this work involves. The first one corresponds to obtain the inference by defining the Bayesian Model in a specific syntax for a package called py − bbn, which is a Python implementation of probabilistic and causal inference in Bayesian Belief Networks using exact inference algorithms (See documentation for further information: https://py-bbn.readthedocs.io/ (Accessed on: 28 March 2022)). The second way corresponds to the classical statistics reinterpretation of the BM, introduced in Section 3.1. Finally, the network reliability can also be computed by the linear formulation integrated in the MPC, derived from the classical statistics and discussed in the previous section. Moreover, it includes the demonstration of the expression concerning the β j summaries in the states matrix of the augmented system, since it is not obvious from the original Equation (26).
The example network is very simple, with only one demand, two sources, two tanks, and five actuators, leading up to three feasible paths: Then, for the linear formulation, the mathematical procedure would start by computing each of the β j for each path. For this, the actuators' reliabilities are involved (correspondence: A i =⇒ R i ). Once they are computed, they should be multiplied by the sum of the logarithms of the reliabilities of the same actuators involved added together, obtaining the following expanded expression (where F d1 stands for the unreliability of the unique demand), It is easy to realise the possibility of the expression rewriting by taking advantage of the commutative properties of addition and multiplication to get the logarithms of the actuators' reliabilities as the common factors, instead of the β j . So, the following equation is equivalent to (33), and it is the one implemented in the augmented system.
Finally, considering R 1 = R 2 = R 3 = R and R 4 = R 5 = 1, the results for different cases are included in the next table.
As can be seen from Table 1, the results from the classical statistics and the linear formulation are exactly the same, which makes sense since the latter is derived from the former. Furthermore, the results obtained by the Bayesian Networks Python package are almost the same, so the method is hereby validated. However, it differs a bit, and even increases the difference when the reliability takes lower values. This may happen due to the algorithms implemented in the library.

Optimization Problem
Regarding the objective function, the economic MPC formulation presented in (10) is modified by adding the next term: where W r is a positive scalar used as a weight to prioritize the terms in the complete objective function. Leading the following new multi-objective cost function to be minimized in the prediction horizon N p : To assure the feasibility of the obtained control actions, the cost function must be subject to the system constraints introduced in Section 2.1, as well as to the one related with the safety objective. So, at each time instant, the following optimization problem is solved online: subject to:

Control Scheme
The scheme of Figure 4 uses blocks to represents the different parts involved in the control system, and where the signals come from. As can be seen, the states come from two different blocks: one represents the plant of the real network, from where the tanks levels would be obtained by sensors;, while the other represents the computed dynamic Bayesian network, which recalculates the reliabilities given the actual control action.

Case Study
To evaluate the approach proposed in Section 3, a part of the Barcelona DWN, presented in Ocampo-Martínez et al. [15], is used as the case study. This network, corresponding to the one in Figure 5, includes nine sources, consisting of five underground and four surface sources, which currently provide an inflow of about 2 m 3 /s. It is composed of 17 tanks, 12 nodes, 25 demands, and 61 actuators (valves and pumps). Following the steps of the mentioned section, the network can be simplified to a graph where the nodes are the tanks/nodes and the links are the actuators (Figure 6) in order to find the minimal paths easily. Then, the corresponding reliability computation is realised.

Results and Discussion
To show and assess the effect of including the reliabilities awareness to the MPC optimization, its related weight of the relevant term in the objective function (W r ) can be set greater than 0; otherwise, the controller would not take them into account. Besides, the total economic cost and the final reliability of the whole system could be good indicators to evaluate their contributions. Then, a 4 days horizon simulation of the introduced case study was performed with some demand and cost values taken from [15]. The following figure includes the results of the mentioned simulation considering different objective function terms' weights, aiming to highlight their influence. In the combinations, only two of the weights were changing (W e and W r ); the other two were fixed: W ∆u = 50; W s = 10. Figure 7 is composed of two axes: the left one corresponds to the accumulated operational cost over time, and the functions related to this axis are those starting at 0 of the same axis; the right one corresponds to the whole system reliability evolution over time, and the functions related to this second axis are those starting at 1. Then, there is one pair of functions for each combination with the same style, detailed in the figure's legend, with one for each axis. The first observable evidence from these results, and an expected outcome, is that taking the overall reliability into account makes the whole system reliability remain closer to 1-otherwise it gets worse.
However, it was expected that the reliabilities integration would negatively affect the total operational costs; but surprisingly, if the economic term weight of the objective function is not 0, the cost remains practically the same and keeps the same trend.
Regarding the two weights modified to see their influence on the system, the economic term weight does not seem to offer a variability on the results but rather the choice of whether or not to optimize the production. This owes to the fact that, until it was disabled by a null weight, the trend was the same, and, once disabled, the trend was held. On the other hand, the reliability term weight seems to yield a minimal variance, but it looks like it is restricted by the economic term, since until it was not disabled by a null weight, and the reliability is not optimized in a considerably better way.

Conclusions
This paper proposes an innovative approach of a health-aware economic model predictive control (EMPC) for the management of generalized flow-based networks. The main enhancement with respect to some existing approaches relies on the dynamic integration of the Bayesian model of the whole system. It is included in the controller to manage the supply, by taking into consideration the distribution of the control effort, to extend the life of the network by delaying the reliability decay as much as possible. This is considered in the MPC optimization problem as a new objective term, and it is implemented with an equivalent interpretation of the Bayesian inference in a linear manner to be able to solve the MPC.
In addition, previous related works have been considered to keep some interesting features, such as the inclusion of an optimal inventory replenishment policy based on a desired risk acceptability level, leading to the availability of safety stocks for unexpected excess demand in networks.
The proposed implementation has been illustrated with a real case study corresponding to a sector of the water transport network of Barcelona. Although the results have a coherent behaviour, most likely, with other network examples, the results could become more evident and noticeable in terms of the influence of the weights; due to it depending intrinsically on the configuration of the network. For example, if there are multiple actuators that work in parallel, and which can then be alternated to distribute the effort, the reliability could be better preserved; otherwise, no matter how much you tune the weights, if the network does not offer versatility, the results are rigid.
In future work, the study of predicting the system reliability of water distribution networks by adding the consideration of the pressure model will be addressed. It will also be interesting to implement the approach to other types of flow-based network systems, such as power networks.
Author Contributions: Conceptualisation, J.P., V.P. and F.N.; methodology, J.P., V.P. and F.N.; software, J.P.; writing-original draft preparation, J.P.; writing-review and editing, J.P., V.P. and F.N.; supervision, V.P. and F.N. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: All the required data are included in the paper.

Conflicts of Interest:
The authors declare no conflict of interest.