System reliability aware model predictive control framework

System Reliability Aware Model Predictive Control Framework Jean C. Salazara,∗, Philippe Weberb, Fatiha Nejjaria, Ramon Sarratea, Didier Theilliolb aUniversitat Polit`ecnica de Catalunya, Research Center for Supervision, Safety and Automatic Control (CS2AC), 10 Rambla Sant Nebridi, Terrassa, Spain bUniversit´e de Lorraine, Centre de Recherche en Automatique de Nancy (CRAN), FST - B.P. 70239 54506 Vandoeuvre-l`es-Nancy, France Abstract This work presents a Model Predictive Control (MPC) framework taking into account the usage of the actuators to preserve system reliability while maximizing control performance. Two approaches are proposed to preserve system reliability: a global approach that integrates in the control algorithm a representation of system reliability, and a local approach that integrates a representation of component reliability. The trade-off between the system reliability and the control performance should be taken into account. A methodology for MPC tuning is proposed to handle this trade-off. System and component reliability are computed based on Dynamic Bayesian Network. The effectiveness and benefits of the proposed control framework are discussed through its application to an over-actuated system. Keywords: Reliability, Dynamic Bayesian Networks, Model Predictive Control, Reliability Importance Measures 1. Introduction The degradation of physical components in engineering systems is generally inevitable. In particular, the degradation of actuators in a closed-loop control system can lead to poor performance and sometimes in a loss of controllability. Therefore, actuator health is of great importance for the safety and reliability of the controlled system. Thus, to avoid failures it is important to enhance system safety by taking into consideration the reliability of components in the controller design [12]. If the design objective is still to maintain the original system performance, this may force the remaining actuators to work beyond their normal duty to compensate the handicaps caused by the fault. Therefore, the tradeoff between achievable performance and available actuator capability should be carefully considered in all control designs. Considerable research has been carried out in order to enhance the system reliability from the manufacturing and system structure point of view. Recently, system reliability has been taken into account in the system control process through a Prognosis and Health Management (PHM) framework. Mainly because reliability is as a measure of how long the system will perform its function correctly [7] and can be used to predict future failures in the system given the state of its components. In some cases the control effort can be redistributed among the available actuators to relieve the work load and the stress factors on assets with worst conditions avoiding in this manner their deterioration. For this purpose, ∗Corresponding author Email address: [email protected] (Jean C. Salazar) an appropriate policy should be developed to redistribute the control effort until maintenance actions can be taken. In this work, the actuator and system reliability are integrated in the control design and used as a policy which serves to redistribute the control effort among the actuators [23, 24]. This paper presents a method to manage the system reliability while performing the control of the system within a framework based on Model Predictive Control (MPC). For this purpose, the paper presents two approaches, the local approach, based on the components reliability and the global approach focused on the system reliability. As a result, a trade-off between system reliability and control performance arises. The growing importance of maintenance has generated an increasing interest in the development and implementation of optimal maintenance strategies for improving system reliability, preventing the occurrence of system failures, and reducing maintenance costs of deteriorating systems [9]. MPC is an efficient technique to manage this kind of objectives. It allows the incorporation of several criteria in the optimization problem. By instance, in [8] the authors present an application of MPC to a Drinking Water Network (DWN) which includes in the optimization problem several criteria such as a criteria for economic cost, level of service and level of components degradation. In [21] and [22], the MPC formulation includes the actuator usage as constraints, whose objective is to maintain the accumulated usage under a safety level at the end of the mission. Reliability can be modelled as an exponential function [6, 33], as a Weibull function [2, 11] or a Gamma function Postprint submitted to Reliability Engineering & System Safety December 11, 2017 List of acronyms and notation BNs Bayesian Networks CAU Cumulative Actuator Usage CPT Conditional Probability Table DAG Directed Acyclic Graph DBNs Dynamic Bayesian Networks DWN Drinking Water Network IBBirnbaum’s measure IMs Importance Measures ISE Integrated Square Error JCR Joint Component Reliability MCs Markov chains MPC Model Predictive Control MIF Marginal Importance Factor PHM Prognosis and Health Management RBD Reliability Block Diagram RIMs Reliability Importance Measures SIMs Structural Importance Measures IβIdentity matrix of size β×β [Λ]βBlock column matrix composed by β×1 blocks of matrices Λ TΛ βBlock lower triangular matrix composed by β×β blocks of matrices Λ uiLower bound of the control input for the iactuator uiUpper bound of the control input for the iactuator HcControl horizon HpPrediction horizon ρEnergy weighting parameter εWeighting parameter of cost function objectives λiFailure rate of the ith component λ0 iBaseline failure rate of the ith component Ckkth minimal cut set Pssth minimal path set P(A) Probability of event A P(A|B) Probability of event A given B RiReliability of the ith component RReliability of the system TsSampling time Φ(·) Structure function eiRandom variable representing the state of the ith component SRandom variable representing the state of the system XiRandom variable representing the state of the ith component of the system XSystem state vector composed of Xielements TMMission time [17, 18, 26], among others. It can also be expressed as a stochastic process [20], using a Markov chain (MC) to model the reliability of the components with the drawback that in practice and depending on the complexity of the system this can lead to a combinatorial explosion of states [32]. The use of Bayesian Networks (BN) to model and computing the reliability of the system, taking into account observations (evidences) about the state of the components has been considered recently in some works [1, 2, 27, 30]. In [25] the authors present the advantages of BNs in comparison with reliability block diagrams. In [4], a fault tree is modelled using a BN. A comparison between MC and Dynamic Bayesian Networks (DBN) applied to reliability is presented in [31]. DBNs are interesting because they allow to model the system reliability with a factorization of the MC states leading to a compact model representation. The purpose of this paper is to present the benefit of taking into account system and component reliability as criteria in the MPC algorithm through an illustration dedicated to DWN as a part of a PHM strategy. The paper presents a study on local and global approaches to manage the system reliability, using Reliability Importance Measures (RIMs). The global system reliability is modelled using a DBN. The rest of the work is organized as follows: Section 2 deals with the components and system reliability modelling. Section 3 presents the RIM used in the global approach. In Section 4 a description of the MPC framework is presented. Section 5 presents the application of the MPC framework on a DWN system and some results of the control application are discussed. Finally in Section 6 some conclusions are given. 2. Reliability Modelling 2.1. Reliability concept Before introducing the concept of reliability, let us first define the concept of failure rate. The general definition of failure rate, denoted as λ, is stated as the fraction of the density of the random lifetime to the survival function (i.e. conditional probability). From this definition, the failure rate can be interpreted as the probability that a component will fail at the next instant time (k+ ∆t) given that it has survive until the current instant time (k). In this work, the definition given by [5] is used, where the failure rate for the ith item is modelled as: λi(k) = λ0 i×g(`, ϑ) (1) where λ0 irepresents the baseline failure rate or nominal failure rate and g(`, ϑ) is a load function also known as 2 covariate, that represents the effect of stress on the component failure rate. `represents an image of the applied load and ϑis a component parameter. Different definitions of function g(`, ϑ) exists in the literature. However, the exponential form is the most commonly used. In [13, 14] the authors propose a load function based on the root-mean-square of the applied control input (ui) until the end of the mission (TM), and an actuator parameter defined from the upper and lower bounds of control input. Moreover, the operating conditions and environmental factors, such as stress, load, temperature, pressure, etc., under which device operates must be taken into account. In this work, it is assumed that the failure rate is provided by: λi(k) = λ0 igi(ui(k)) (2) In (2) gi(ui(k)) represents the amount of load on the actuator and corresponds to the following normalized control action: gi(ui(k)) = ui(k)−ui ui−ui (3) where ui(k) is the control effort at time k,uiand uiare the minimum and maximum control efforts allowed for the ith actuator. The major actuator load corresponds to ui(k) = ui, which leads to the worst failure rate λi=λ0 i. Accordingly to the failure rate definition, reliability is the probability that a system will perform its functioning satisfactorily for a given period of time and under stated operating conditions. In this paper an exponential model for the component reliability is used: Ri(t) = e−Rt 0λi(v)dv (4) and in discrete time it follows: Ri(k+ 1) = e−Ts k+1 P v=0 λi(v)(5) where λi(v) is the failure rate that is obtained from the ith component under different levels of load and Tsis the sampling time. For analysis simplification but without loss of generality, it is assumed that all the components are mutually independent. This means that the joint distribution of the system state vector (X= (X1, . . . , Xn)) is determined by the reliability of its components. The dependence of the state of the system regards to the state of the components can be determined by means of its structure function (Φ(X)). This function indicates the status of the system (success or failure) given the states of each component. For instance, the structure function can be computed using either the minimal path sets (Ps) which are the minimal sets of elements of the system whose functioning (i.e., being up -Up-) ensures that the system is up (6), or the minimal cut sets (Ck) which are the minimal sets of elements of the system whose failure (i.e. being down -Dn-) causes the failure of the system (7). Φp(X)=1− s Y j=1  1−Y i∈Pj Xi (6) Φc(X) = k Y j=1  1−Y i∈Cj (1 −Xi) (7) where Xiis a random binary variable that represents the state of the component: Xi= 1 the component is Up, Xi= 0 the component is Dn. As both of methods are equivalent in this work the minimal path set approach is used. Hence, the reliability of the system is defined as the expectation of the structure function, this is: R=E[Φ(X)] (8) which can be expressed in probability terms as: R= P(Φ(X) = 1) (9) Therefore, if Φ(X) = 1 the system is up (Up), and if Φ(X) = 0 the system is down (Dn). Generally, the reliability is determined by the way the components are interconnected in the system, i.e., serial or parallel or a combination of them. 2.2. Bayesian Networks & Dynamic Bayesian Networks Basically, Bayesian Networks compute the probability distribution in a set of variables according to the prior knowledge of some variables and the observation of others [10]. BNs can be considered as a subset of Directed Acyclic Graph (DAG). To illustrate this concept, let Aand Bbe two nodes with two possible states (S1and S2) (Figure 1). A probability is associated to each state of the node Aand this probability is defined a priori for a root node and computed by inference for the others. The a priori probabilities of node Aare P(A=SA1) and P(A=SA2). Figure 1: Basic Bayesian Network. A Conditional Probability Table (CPT) is associated with node Band defines the conditional probability of the state of Bgiven the state of A(P(B|A)). Thus, the BN inference computes the marginal distribution e.g. P(B=SB1) as: P(B=SB1) =P(B=SB1|A=SA1)P(A=SA1) + P(B=SB1|A=SA2)P(A=SA2) (10) 3 A special case of Bayesian Network called Dynamic Bayesian Network (DBN) is used to model the time dependency of a process. In a DBN the random variables are time indexed, so the process state SA:{sA 1, ..., sA M}is represented by nodes Akand Ak+∆tat time instants kand k+ ∆t, respectively. The time dependence is represented by an arc and the temporal evolution of the variables is represented by successive time slices as shows Figure 2 (for simplicity let ∆tbe equal to 1). The transition probability between the states of the variable at time instant kto k+ 1 is defined by a CPT intertime slices. Given the probability distribution P(Ak) at time instant k, the network leads by inference to a unique distribution P(Ak+1) at time instant k+ 1. The computation over time of P(AT) is performed by iterative inferences starting from k= 0 and P(A0) [1]. Figure 2: DBN model for the ith component. 2.3. Components reliability modelling The decay of the components reliability can be modelled using a Markov chain process. This type of process is very useful in system reliability modelling due to its memoryless property, which means, that the state transitions only depend on the current and next state, and it does not depend on the previous state or the amount of time that the process has stayed in the current state, making it an effective tool for system reliability analysis [28]. Let eibe a discrete random variable in the Markov process representing the state of the ith component with two possible mutually exclusive states {Up,Dn}. The probabilistic state transition between the states is defined by: PMC =1−p12 p12 0 1 (11) where p12 ∼ =λ∆t:λrepresents a constant failure rate and ∆tthe time interval, and p12 can be interpreted as the probability that the component goes from state Up to Dn after ∆t. In the case of components whose failure rate depends on time, their reliability decay can be modeled by a semi- Markov chain (Figure 3). The MC is homogeneous if the transition probabilities of the states are independent on time. In a semi-Markov chain, the memoryless property is relaxed, which means that the state transition depends on the current state and the transition time. Therefore, the probability matrix (11) becomes dependent on time: PMC(eik+1|eik) = 1−p12(k)p12(k) 0 1 (12) Figure 3: Semi-Markov chain for the component reliability. From (12) it is clear that the failure rate is the probability of component state to be Dn at instant time k+ 1 given that its state was Up at time k, this is: λi(k) = P(ei(k+ 1) = Dn|ei(k) = Up) (13) The semi-Markov chain is modelled by a Dynamic Bayesian Network with CPT evolving according to time. In this case the evolution is represented in two time slides (Figure 4) [1]. Figure 4: DBN model for the ith component. Therefore, DBN computes the component reliability using: Ri(k+ 1) = P(ei(k+ 1) = Up) (14) which follows expression (5). 2.4. System reliability modelling A system is defined as a set of icomponents whose states are described by the random variable ei, as: P(Xi= 1) = P(ei=Up); P(Xi= 0) = P(ei=Dn) (15) where Xi:Xi= 1 if the ith component is up, Xi= 0 if the ith component is down. Assume that the state of the components are known. Therefore, in this work the system reliability is computed from its components reliability using a DBN. For instance, once the minimal path sets of the system are identified, its states (i.e., if the minimal path set is up or down) are represented by a random variable as a node in the DBN. The state of the components of the system are also represented as nodes in the DBN and are connected to their respective minimal path set node. Note that the probabilities of the component states are computed by inferences in the DBN modelling of the semi Markov chain described in Section 2.3. Finally, the state of the system (i.e., system reliability) is given by the top node of the BN which is connected to the minimal path set nodes. 4 For example, consider the reliability block diagram (RBD) of Figure 5, which shows the contribution of component reliability to the success or failure of a three components system. It is clear that with a minimum of two components the system can performs its function satisfactorily, being {1, 3}and {2, 3}the minimal path sets of the system. Although, in this example a RBD is used to define the structure of the BN, the minimal path sets can also be obtained by other methods or directly from the system. This will be illustrated in the case study in Section 5. Figure 5: RBD of a three components system example. Then, it is possible to build a BN of the system reliability from its minimal path sets see Figure 6, being ei the probability of the ith component state, Pithe probability of the state of the ith minimal path set and Sthe probability of the system state [27]. Figure 6: Bayesian Network of system reliability. Tables 1 and 2 present the CPTs of nodes P1and S, respectively. P1depends on the states of the components e1 and e3and its behaviour corresponds to an AND gate, i.e., all the components in a success path should be available for the system to be available. Sdepends on the state of success path nodes P1and P2and has the behaviour of an OR gate, i.e., if there is at least one success path available, then the system will be available. Table 1: CPT for node P1. e1e3P1 Up Dn Up Up 1 0 Up Dn 0 1 Dn Up 0 1 Dn Dn 0 1 Table 2: CPT for node S. P1P2S Up Dn Up Up 1 0 Up Dn 1 0 Dn Up 1 0 Dn Dn 0 1 It is possible then to compute the probability distribution for each variable conditioned by the values of the other variables in the graph. This feature is particularly important in case a control system must work in real-time, because in that case evidences acquired about a state variable must be propagated to update the state of the rest of the domain. Therefore, the reliability of the system (R) is computed using the BN as: R(k+ 1) = P(S(k+ 1) = Up) (16) In the case of complex structure systems with high amount of components the compute of the structure function becomes non trivial. The use of BNs to reliability modelling allows several modelling structures based on minimal path sets, or logical combination of components states (using AND or OR gates) [15]. 3. Reliability Importance Measures One of the objectives of system reliability analysis is to identify the weakness in a system and to quantify the impact of components failures over system functioning. Several indicators of reliability importance exist, each expressing the importance from a slightly different point of view. Importance Measures (IMs) were first introduced by [3] and are classified in two groups: Reliability Importance Measures (RIMs) and Structural Importance Measures (SIMs). The RIMs evaluate the relative importance of a component taking into account its contribution to the system reliability and the SIMs provide the relative importance of a component taking into account its position into the system structure. These metrics can be defined according to their functional aspect, leading to their modelling by the minimal path sets, and according to their dysfunctional aspect, leading to their modelling by the minimal cut sets. As both are equivalent, in this study only the functional aspect is used. The aim is to quantify the importance of the ith component for the reliability of the total system and how the changes in the component reliability impact in the system reliability. One of the most used RIMs is the Birnbaum’s importance measure [3] also known as Marginal Importance Factor (MIF). MIF is related to the probability that a component is critical for the system. Definition 1. The Birnbaum’s importance measure of the ith component for the functioning of the system, denoted as IBi(i;p), for a coherent system with independent components it is defined as: IB(i, p) =∂R(p) ∂pi =P(S=Up|ei=Up)−P(S=Up|ei=Dn)(17) This measure is well known from classical sensitivity analysis. It represents the probability that the failure and functioning of the ith component coincide with the failure 5 and functioning of the system. This can be interpreted as the maximum decrease of system reliability when the ith component changes from the condition of perfect functioning to failed condition. Note that the Birnbaum’s importance measure of the ith component only depends on the structure of the system and the reliabilities of the other components, and that it is independent of the actual reliability of the ith component. The Birnbaum’s importance measure of component i can be computed by inference using the DBN system model by introducing evidences about the states of each component (i.e., being Up or Dn) and computing the difference between the system reliabilities in both cases. 4. Reliability aware MPC Framework 4.1. MPC formulation Consider the following linear discrete-time dynamic model of a system described in the state-space form as: x(k+ 1) = Ax(k) + Bu(k) y(k) = Cx(k)(20) where for each k∈Z+, x(k)∈Rnis the state, u(k)∈Rp is the control input, y(k)∈Rqis the measured output, A∈Rn×nis the state matrix, B∈Rq×pis the input matrix and C∈Rq×nis the output matrix. Model-based predictive control (MPC) is a discrete time technique where an explicit dynamic model of the plant is considered to predict the system outputs by applying (20) recursively along the prediction horizon, obtaining (18) (Hp) where ˆu(k+i|k) is the control input corresponding to k+icalculated at time instant k, and the hat notation represents a predicted value. The control actions are manipulated throughout a finite control horizon (Hc) in order to minimize a given cost function, with Hc< Hp. This methodology is represented in Figure 7, where the past outputs and control input sequences (showed in the left side) are used to compute the future output sequence over Hpand future control inputs      ˆy(k+ 1|k) ˆy(k+ 2|k) . . . ˆy(k+Hp|k)      | {z } ˆ Y =     CA CA2 . . . CAHp      | {z } Fu x(k) +        CB 0q×p· · · 0q×p CAB CB · · · 0q×p . . .. . ....0q×p CAHp−1B CAHp−2B· · · Hp−Hc P i=0 CAiB        | {z } H      ˆu(k|k) ˆu(k+ 1|k) . . . ˆu(k+Hc−1|k)      | {z } ˆ U (18)      ˆu(k|k) ˆu(k+ 1|k) . . . ˆu(k+Hc−1|k)      | {z } ˆ U =     Ip0p×p· · · 0p×p IpIp· · · 0p×p . . .. . ..... . . IpIp· · · Ip      | {z } TIp Hc      ∆ˆu(k|k) ∆ˆu(k+ 1|k) . . . ∆ˆu(k+Hc−1|k)      | {z } ∆ˆ U +     Ip×p Ip×p . . . Ip×p      | {z } [Ip]Hc ˆu(k−1) (19) Figure 7: Model Predictive Control basic concept sequences over Hc. Such that the cost function is minimized. The cost function generally includes a term for the tracking error, minimizing the square difference between the predicted output ( ˆ Y) and the set-point (Yref ), and a term for the minimization of the square of the energy ( ˆ U) and its variation (∆ ˆ U) according to (18) and (19). In this work, the cost function Jis defined as: min (ˆ U(k|k), ..., ˆ U(k+Hc −1|k), ∆ˆ U(k|k), ..., ∆ˆ U(k+Hc−1|k)) J=εhˆ Y−Yref i>hˆ Y−Yref i + (1 −ε)ˆ U>ρˆ U+ ∆ ˆ U>∆ˆ U (21) subject to U≤ˆ U≤U X≤ˆ X≤X where: Yref = [yref (k+ 1), ..., yref (k+Hp)]>is the output reference, ρ=diag(ρ, ρ, ..., ρ) of size Hc×Hcand ρ=diag(ρ1, ρ2, ..., ρp) is the weighting for the energy and ε is a weighting parameter between tracking and energy objectives, U,U,Xand X, are the lower and upper bounds of the control effort and system states, respectively. The resulting control action ˆu(k|k) is injected to the system while ˆu(k+i|k)∀i= 1, ..., Hc−1 are discarded. At 6 next instant time k+ 1, y(k+1) is measured and the optimization problem is solved again. Thus, ˆu(k+ 1|k+ 1) is calculated moving the prediction horizon using the concept of sliding window. 4.2. MPC tuning methodology MPC tuning involves finding the appropriate values for the weighting parameters in the cost function, as well as the prediction and control horizons that depend on the application and the amount of data that can be handled. This work focuses on proposing a methodology to find ρand εin (21) as follows: Step 1: Select ε= 0 and through simulations find the optimal value of ρin which the reliability of the system at the end of the mission time is the highest. The comparison and selection of the best approach is performed based on the following criteria. Let the Joint Component Reliability index (JCR) be a measure of the remained overall reliability of the system components at the end of the mission time: JCR = p Y i=1 Ri(22) and, let the Cumulative Actuator Usage (CAU) be a measure of the actuator energy consumption defined as: CAU = Ts TM/Ts X k=0 U(k)>U(k) (23) where TMis the mission time. Step 2: Having selected ρas mentioned in the previous step, find the optimal value of εthrough simulations, starting from ε= 0 to 1. The idea is to study the effect produced by the variation of parameter εin the system reliability or in other words, the impact of the tracking error in the reliability. To measure this impact, an index called the integral square error (ISE) is defined: ISE = Ts TM/Ts X k=0 hˆ Y(k)−Yref (k)i>hˆ Y(k)−Yref (k)i (24) Therefore, the εvalue which corresponds to the highest system reliability and the lowest ISE. 4.3. Control action redistribution The weights ρin the cost function (21) redistribute the control effort among the actuators [23]. In this work, a global and a local approach are proposed to handle the reliability in the optimization problem. The global approach consists in considering the Birnbaum’s importance measure viewed in Section 3 (computed using the DBN reliability model) by setting: ρ(k) = diag(IBi(k)) ∀i= 1,2, . . . , p (25) in this case, it is expected that components with a greater contribution to the system reliability are used less than the others. On the other hand, the local approach tries to preserve system reliability by preserving component reliability by setting the weight as: ρ(k) = diag(1 −Ri(k)) ∀i= 1,2, . . . , p (26) This criteria aims at finding the optimal control actions and distribute it among the available actuators in such a way that actuators with lower reliability level are relieved. Hence, the use of highly reliable components is prioritized. The local approach assumes an equivalent contribution of component reliability to system reliability. However, this is hardly ever true. In fact, the DBN reliability model can intrinsically explain this relation. The control strategy scheme is presented in Figure 8. On the one hand, the MPC computes the control inputs according to: the cost function, a set of bounds, the current system state and the weight ρ. Then, the control input is injected to the system and used to compute the component failure rates. Figure 8: Block diagram of the approach On the other hand, those failure rates are used in the DBN to compute: the components reliability, the overall system reliability and the Birnbaum’s importance measure. This data is used then to update the weight ρ, which is used in the MPC algorithm, closing in this way the loop. 5. Application: Drinking Water Network 5.1. System description A DWN is composed by sources (water supplies), sinks (water demand sectors) and pipelines that link sources to sinks. It also contains active elements like pumps, valves. 7 Figure 9: Drinking water network diagram. The considered DWN corresponds to a subsystem of the Barcelona water transport network [30] (Figure 9). The network consists of 5 sources and 1 sink. It is assumed that the demand forecast at the sink (dm(k)) is known (Figure 10), and that any single source can satisfy this required water demand. It is also assumed that the volume of the tanks should follow a given set point (Figure 11). 0 50 100 150 200 250 300 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Time [h] Flow [m3/s] Figure 10: Drinking water demand. The DWN is modelled by applying mass balance to each tank and the following linear discrete-time model is obtained: x(k+ 1) = Ax(k) + Bu(k) + Bddm(k) y(k) = Cx(k)(27) where x(k)∈Rnare tanks volume, u(k)∈Rpare the control inputs (pump commands) with u(k)≥0∀k,y(k) ∈Rqare the measured tanks volume, dm(k)∈Rmis the water demand, A∈Rn×nis the state matrix, B∈Rq×p input matrix, Bd∈Rn×mis the disturbance matrix, and C∈Rq×nis the output matrix. 0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 1 2 3 4 5 x 104 Volume [m3] Time [h] Setpoint tank 1 Setpoint tank 2 Setpoint tank 3 Setpoint tank 4 Figure 11: Volume reference for the 4 reservoirs 5.2. Dynamic Bayesian Network modelling The DWN reliability is modelled using a DBN following the methodology explained in Section 2.4. Firstly, system components must be identified. In this case there are 10 pumps, 5 sources, 4 tanks and several pipes. Secondly, the minimal path sets should be determined. A minimal path set is composed by those components which allow a flow path between sources and sinks, such as pipes, tanks and pumps. A list of the components that correspond to each minimal path set is presented in Table 3. There are nine minimal path sets in the system of Figure 9. Each minimal path set is available depending on the reliability of its components. Pipes and tanks are considered perfectly reliable so they do not provide significant information to the network. Nevertheless, sources are included in the minimal path sets merely for illustrating the procedure. Table 3: Components and minimal path sets relationship. A1A2A3A4A5e1e2e3e4e5e6e7e8e9e10 P1× × × P2× × × P3× × × P4× × × × P5× × × × P6× × × × × P7× × × × P8× × × × P9× × × × × Provided the information of Table 3, the DBN presented in Figure 12 is built as follows: nodes eiand Aiare drawn for each component. Note that nodes eihave two time slices in time kand k+ 1 following the approach of Section 2.3. Then, these nodes are interconnected to their corresponding minimal path set nodes Piusing arcs. Finally, each minimal path set node is interconnected to the system reliability node S[30]. Initially, at instant k= 0, the pumps and the system are assumed to be fully reliable, i.e. their reliability is 1. Then, the probability of each node is computed using its CPT. 8 Figure 12: Dynamic Bayesian network model of the DWN. At each sampling time, the reliability Riof each pump is computed according to its failure rate using a MC (Figure 12). Its behaviour follows an exponential distribution as stated in (4). Note that it is independent of the previous states of the component. It only depends on its present state. In the DBN, this corresponds to the CPT shown in Table 4. Table 4: Inter-time slices CPT for node eik+1 . ei(k+ 1) ei(k)Up Dn Up 1-(λ0 iTs×gi(ui(k))) λ0 iTs×gi(ui(k)) Dn 0 1 Remark 1. The failure rate is computed according to gi(ui(k)). The CPT shown in Table 4 defines the discretized stochastic process of using (2) and (3) in (4). The CPT of node P1is shown in Table 5. This CPT depends on the states of the source 1 (A1) and pumps 1 and 5 (e1, e5). Its behaviour corresponds to an AND gate. Table 5: CPTs for nodes P1. A1e1k+1 e6k+1 P1 Up Dn Up Up Up 1 0 Up Up Dn 0 1 Up Dn Up 0 1 Up Dn Dn 0 1 Dn Up Up 0 1 Dn Up Dn 0 1 Dn Dn Up 0 1 Dn Dn Dn 0 1 It is assumed that with one source it is possible to satisfy the water demand. Thus, the availability of the system can be assured as long as at least one of paths Piis available, which corresponds to the CPT of node Sshown in Table 6. It depends on the state of nodes P1to P9and has the behaviour of an OR gate. Table 6: CPT for node S. P1P2P3. . . P9S Up Dn Up Up Up . . . Up 1 0 Up Up Up . . . Dn 1 0 - - - . . . Up 1 0 . . .. . .. . ..... . .. . .. . . Dn Dn Dn . . . Dn 0 1 The system reliability of the DWN is computed by implementing the DBN presented in the Figure 12 using the BNT toolbox for Matlab [19]. 5.3. MPC with reliability optimization The control of the DWN system is performed applying the MPC formulation of Section 4. Table 7 provides the simulation parameters used. A hierarchical control structure is assumed, where the MPC layer produces a set of set-points for the lower level flow controllers. Hence, a 1 hour sampling time is assumed for the upper level MPC. Since the water demand (Figure 10) exhibits a daily profile, a 24 h prediction horizon has been chosen and the initial tank volumes have been set to X0. Table 7: Simulation parameters parameter value Hp|Hc24 |8 Ts|TM[h] 1 |2000 ρ|ε{1, 1 −Ri, IBi} | {10−18, ..., 1} U[m3/s]0.75 0.75 0.75 1.20 0.85 1.60 1.70 0.85 1.70 1.60 U[m3/s]0 0 0 0 0 0 0 0 0 0 λ0[h−1×10−4]9.85 10.70 10.50 1.40 0.85 0.80 11.70 0.60 0.74 0.78 X[m3] 65200 3100 14450 11745 X[m3] 25000 2200 5200 3500 X0[m3] 45100 2650 9825 7622 To find the appropriate values for the weighting parameters the procedure presented in Section 4 is followed. First, Step 1 is applied and the following simulation scenarios are considered: ρi= 1 −Ri,ρi= IBiand ρi= 1 where no dependency on system reliability is assumed. 9