Modelling and Resolution of Dynamic Reliability Problems by the Coupling of Simulink and the Stochastic Hybrid Fault Tree Object Oriented (SHyFTOO) Library

: Dependability assessment is one of the most important activities for the analysis of complex systems. Classical analysis techniques of safety, risk, and dependability, like Fault Tree Analysis or Reliability Block Diagrams, are easy to implement, but they estimate inaccurate dependability results due to their simpliﬁed hypotheses that assume the components’ malfunctions to be independent from each other and from the system working conditions. Recent contributions within the umbrella of Dynamic Probabilistic Risk Assessment have shown the potential to improve the accuracy of classical dependability analysis methods. Among them, Stochastic Hybrid Fault Tree Automaton (SHyFTA) is a promising methodology because it can combine a Dynamic Fault Tree model with the physics-based deterministic model of a system process, and it can generate dependability metrics along with performance indicators of the physical variables. This paper presents the Stochastic Hybrid Fault Tree Object Oriented (SHyFTOO), a Matlab ® software library for the modelling and the resolution of a SHyFTA model. One of the novel features discussed in this contribution is the ease of coupling with a Matlab ® Simulink model that facilitates the design of complex system dynamics. To demonstrate the utilization of this software library and the augmented capability of generating further dependability indicators, three di ﬀ erent case studies are discussed and solved with a thorough description for the implementation of the corresponding SHyFTA models.


Introduction
Dependability is the ability to avoid service failures that are more frequent and severe than is acceptable, and it is comprised of the following attributes [1]: reliability: continuity of correct service; -availability: readiness for correct service; -safety: absence of catastrophic consequences on the user(s) and environment; -maintainability: ability to undergo modifications and repairs.
The dependability assessment of a system is a crucial activity for engineers and industrial stakeholders so as to ensure the dependable operation of the system. Several different stochastic formalisms have been developed to evaluate the different dependability attributes of complex ). Each of these formalisms has different advantages and limitations [2], and accordingly, the choice of the appropriate dependability assessment methodology depends on the complexity of the system, amount of available information, and required key performance indicators.
One of the main shortcomings in the above-mentioned formalisms is that they cannot model the dependency between the failure behavior of the components and the system working conditions. In fact, the probability distribution of the time to failure of the system components is fixed and must be provided as a static input of the model. Driven by the need for an accurate dependability analysis process, in order to tackle the previous limitation, during the last few decades, researchers have developed hybrid models that are able to account for the stochastic and deterministic failure and operation processes including random phenomena and physics-based evolutions. The branch of dependability modelling that considers stochastic and deterministic processes is known as Dynamic Probabilistic Risk Assessment (DPRA) or Dynamic Reliability [3]. In a DPRA model, the physics of a process can influence the stochastic behavior of the system components. One of the main difficulties of such a modelling process is the level of detail that influences the model size and its complexity [4].
The system dependability assessment is carried out by a team of experts that may eventually include risk practitioners, mathematicians, and engineers [5]. The common goal is the conception of a model that defines the system failure and operation processes to quantify the occurrence probability of undesired events. This process requires a detailed knowledge of the system processes and the ability to describe them with a stochastic model. Computer-aided tools can ease this activity, and the literature presents several contributions that demonstrate the potential of academic tools and their main scope of application with respect to the process under analysis and the stochastic modelling approach [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22].
Following this research stream, in this paper, the authors focus on the modelling of Stochastic Hybrid Fault Tree Automaton (SHyFTA) and illustrate how to use the SHyFTOO Matlab ® toolbox to implement and solve general DPRA problems. SHyFTA is a stochastic formalism of DPRA based on Repairable Dynamic Fault Tree (RDFT) [21], and SHyFTOO is a software library presented in a previous research paper by the authors [22]. The breakthrough of the framework and of the proposed tool, compared to the competitive tools that are now available in the industrial scenario, stands in its native hybrid feature, which is able to combine Time-Driven and Discrete-Event simulation paradigms to solve the continuous dynamic of a physical process and to evaluate at any relevant time event the reliability status of the system.
The aim of this contribution is twofold: (i) illustrate the integration of a SHyFTA model coded with the SHyFTOO library with Matlab ® Simulink toolbox; (ii) present new emerging capabilities of a simulated SHyFTA model that, thanks to the coupling of these two computer-aided components, allows the analysis of more complex dependencies and the measurement of several key performance indexes of a system process.
The rest of the paper is organized as follows. Section 2 resumes the evolution of the stochastic Fault Tree Analysis (FTA) toward the hybrid version that integrates stochastic and deterministic processes. Section 3 illustrates how to model a SHyFTA model with the SHyFTOO library and Simulink with the aid of a case study. Section 4 presents the resolution of three different models to illustrate a portfolio of applications that can be analysed, exploiting the SHyFTOO library to compute the main dependability attributes and other key performance indicators. In Section 5, conclusions and future works are summarized. Moreover, the Appendices A and B contain several information about the SHyFTOO library that can aid interested readers and practitioners in using this tool.  Figure 1 shows a breakdown of the three main categories of the stochastic methods used in the dependability assessment, and Table 1 displays their characteristics and modelling hypotheses. As shown in Figure 1, Fault Tree Analysis (FTA) is a common pattern of each category because it has been object of continuous modification over the last few years, improving its expressiveness and accuracy. A state-of-the-art survey demonstrates that several modifications of the original FTA formalism have been attempted to overcome the limitations of the SFT. These techniques are grouped into the second category of models, here named as "Advanced & General-Purpose Methods". The main improvements over the "Traditional Static Methods" focus on the capability to model complex failure dependencies such as time-event sequences, multi-state degradation, and standby policies [33]. However, various methodologies do not have the same properties, and the use of one or the other approach depends on the dependability problem under study.

The Evolution of Fault Tree Analysis
Among the various techniques, Dynamic Fault Tree (DFT) [34] is probably the most known method which addresses temporal events through dynamic gates. DFTs kick-started the proliferation of other formalisms that has continued for two decades now. For instance, a recent research paper [35] proposed a new fault tree modeling paradigm named Temporal Dynamic Fault Trees (TDFT) that is able to calculate the reliability and availability of a system exploiting the notion of soft-faults, i.e., temporary events that disappear after the source of the interference is no longer present. Other contributions [36,37] integrated the primitives of parametric, Repairable Dynamic Fault Trees (RDFT), and introduced a new formalism called Generalized Fault Trees (GFT) with the aim of reducing the cost of analysis and enabling the evaluation of different repair policies. Moreover, they generalized this methodology to handle parametric fault trees and dynamic fault trees. Authors in [38] suggested Boolean Logic Driven Markov Process (BDMP) as an equivalent formalism to evaluate DFT characterized by repair transitions. In [39], an ad-hoc approach for dealing with RDFTs is proposed. The main novelty is the introduction of a simplifying assumption of independency among gate inputs. In RDFTs, the occurrence of a basic event needs to be considered along with its repair process. There are system state-transitions, and these are handled using the concept of renewal  and future works are summarized. Moreover, the Appendices A and B contain several information about the SHyFTOO library that can aid interested readers and practitioners in using this tool. Figure 1 shows a breakdown of the three main categories of the stochastic methods used in the dependability assessment, and Table 1 displays their characteristics and modelling hypotheses. As shown in Figure 1, Fault Tree Analysis (FTA) is a common pattern of each category because it has been object of continuous modification over the last few years, improving its expressiveness and accuracy.

The Evolution of Fault Tree Analysis
Static Fault Tree (SFT) methodology was conceived in 1962 at Bell Laboratories by H.A. Watson, under a US Air Force Ballistics Systems Division contract to evaluate the Minuteman I Intercontinental Ballistic Missile (ICBM) Launch Control System [23]. The SFT formalism provides an easy way to model and understand the failure behaviour of a system. However, its underlying assumptions limit the application to systems with working and failed states. The construction of a fault tree follows a top-down approach that starts with the identification of a major undesired event (top event) as consequence of a cascade of intermediate events which can be modelled using the Boolean logic gates AND and OR. The leaves of the fault tree are called "basic events" and cannot be further decomposed. They are characterized by the probability distribution of the time to occur (or to fail) of the system components that are generally provided by the components manufacturer. Static Fault Tree belongs to the class of Boolean methods. Therefore, they can be solved with analytic algorithms that provide an exact result [24][25][26][27][28][29] also in the case of non-coherent systems, described with negative logic gates [30]. Static Fault Tree is currently one of the preferred methodologies in the industrial sector [31,32] and, in authors' opinion, this prevalence will be maintained as long as the "Advanced & General-Purpose" and the "Hybrid" methodologies do not become of easy access and understanding for the industrial risk practitioners. and future works are summarized. Moreover, the Appendices A and B contain several information about the SHyFTOO library that can aid interested readers and practitioners in using this tool. Figure 1 shows a breakdown of the three main categories of the stochastic methods used in the dependability assessment, and Table 1 displays their characteristics and modelling hypotheses. As shown in Figure 1, Fault Tree Analysis (FTA) is a common pattern of each category because it has been object of continuous modification over the last few years, improving its expressiveness and accuracy.

The Evolution of Fault Tree Analysis
Static Fault Tree (SFT) methodology was conceived in 1962 at Bell Laboratories by H.A. Watson, under a US Air Force Ballistics Systems Division contract to evaluate the Minuteman I Intercontinental Ballistic Missile (ICBM) Launch Control System [23]. The SFT formalism provides an easy way to model and understand the failure behaviour of a system. However, its underlying assumptions limit the application to systems with working and failed states. The construction of a fault tree follows a top-down approach that starts with the identification of a major undesired event (top event) as consequence of a cascade of intermediate events which can be modelled using the Boolean logic gates AND and OR. The leaves of the fault tree are called "basic events" and cannot be further decomposed. They are characterized by the probability distribution of the time to occur (or to fail) of the system components that are generally provided by the components manufacturer. Static Fault Tree belongs to the class of Boolean methods. Therefore, they can be solved with analytic algorithms that provide an exact result [24][25][26][27][28][29] also in the case of non-coherent systems, described with negative logic gates [30]. Static Fault Tree is currently one of the preferred methodologies in the industrial sector [31,32] and, in authors' opinion, this prevalence will be maintained as long as the "Advanced & General-Purpose" and the "Hybrid" methodologies do not become of easy access and understanding for the industrial risk practitioners. and future works are summarized. Moreover, the Appendices A and B contain several information about the SHyFTOO library that can aid interested readers and practitioners in using this tool. Figure 1 shows a breakdown of the three main categories of the stochastic methods used in the dependability assessment, and Table 1 displays their characteristics and modelling hypotheses. As shown in Figure 1, Fault Tree Analysis (FTA) is a common pattern of each category because it has been object of continuous modification over the last few years, improving its expressiveness and accuracy.

The Evolution of Fault Tree Analysis
Static Fault Tree (SFT) methodology was conceived in 1962 at Bell Laboratories by H.A. Watson, under a US Air Force Ballistics Systems Division contract to evaluate the Minuteman I Intercontinental Ballistic Missile (ICBM) Launch Control System [23]. The SFT formalism provides an easy way to model and understand the failure behaviour of a system. However, its underlying assumptions limit the application to systems with working and failed states. The construction of a fault tree follows a top-down approach that starts with the identification of a major undesired event (top event) as consequence of a cascade of intermediate events which can be modelled using the Boolean logic gates AND and OR. The leaves of the fault tree are called "basic events" and cannot be further decomposed. They are characterized by the probability distribution of the time to occur (or to fail) of the system components that are generally provided by the components manufacturer. Static Fault Tree belongs to the class of Boolean methods. Therefore, they can be solved with analytic algorithms that provide an exact result [24][25][26][27][28][29] also in the case of non-coherent systems, described with negative logic gates [30]. Static Fault Tree is currently one of the preferred methodologies in the industrial sector [31,32] and, in authors' opinion, this prevalence will be maintained as long as the "Advanced & General-Purpose" and the "Hybrid" methodologies do not become of easy access and understanding for the industrial risk practitioners.  [23]. The SFT formalism provides an easy way to model and understand the failure behaviour of a system. However, its underlying assumptions limit the application to systems with working and failed states. The construction of a fault tree follows a top-down approach that starts with the identification of a major undesired event (top event) as consequence of a cascade of intermediate events which can be modelled using the Boolean logic gates AND and OR. The leaves of the fault tree are called "basic events" and cannot be further decomposed. They are characterized by the probability distribution of the time to occur (or to fail) of the system components that are generally provided by the components manufacturer. Static Fault Tree belongs to the class of Boolean methods. Therefore, they can be solved with analytic algorithms that provide an exact result [24][25][26][27][28][29] also in the case of non-coherent systems, described with negative logic gates [30]. Static Fault Tree is currently one of the preferred methodologies in the industrial sector [31,32] and, in authors' opinion, this prevalence will be maintained as long as the "Advanced & General-Purpose" and the "Hybrid" methodologies do not become of easy access and understanding for the industrial risk practitioners.
A state-of-the-art survey demonstrates that several modifications of the original FTA formalism have been attempted to overcome the limitations of the SFT. These techniques are grouped into the second category of models, here named as "Advanced & General-Purpose Methods". The main improvements over the "Traditional Static Methods" focus on the capability to model complex failure dependencies such as time-event sequences, multi-state degradation, and standby policies [33]. However, various methodologies do not have the same properties, and the use of one or the other approach depends on the dependability problem under study.
Among the various techniques, Dynamic Fault Tree (DFT) [34] is probably the most known method which addresses temporal events through dynamic gates. DFTs kick-started the proliferation of other formalisms that has continued for two decades now. For instance, a recent research paper [35] proposed a new fault tree modeling paradigm named Temporal Dynamic Fault Trees (TDFT) that is able to calculate the reliability and availability of a system exploiting the notion of soft-faults, i.e., temporary events that disappear after the source of the interference is no longer present. Other contributions [36,37] integrated the primitives of parametric, Repairable Dynamic Fault Trees (RDFT), and introduced a new formalism called Generalized Fault Trees (GFT) with the aim of reducing the cost of analysis and enabling the evaluation of different repair policies. Moreover, they generalized this methodology to handle parametric fault trees and dynamic fault trees. Authors in [38] suggested Boolean Logic Driven Markov Process (BDMP) as an equivalent formalism to evaluate DFT characterized by repair transitions. In [39], an ad-hoc approach for dealing with RDFTs is proposed. The main novelty is the introduction of a simplifying assumption of independency among gate inputs. In RDFTs, the occurrence of a basic event needs to be considered along with its repair process. There are system state-transitions, and these are handled using the concept of renewal processes instead of the commonly used Markov chains.
Other suggested improvements can be found in the literature, such as [40], which integrated prognostic estimations into dependability assessment to avoid using average failure rate data, which leads to outdated system-level dependability estimation, and proposed a new framework, which results in more realistic and accurate dependability predictions. A methodology for combining the Fuzzy logic and FTA are discussed in [41] in order to cope with the information uncertainty due to limited data of the failure components probability, whereas [42] presents an improved methodology for analysing the qualitative importance of components in the functional and reliability structures of a system.
The resolution of the above-mentioned models is not as simple as SFT, and accordingly, the following three techniques are generally adopted: (i) algebraic methods, (ii) conversion into an equivalent model, and (iii) simulation. As far as it concerns the algebraic methods, [43,44] presented the algebraic relations describing the temporal priorities and the quantitative analysis calculation of the gates. Similar works are addressed in [45], which presents a special version of sequence algebras, representing failures of non-repairable components, and proposes a "sequence decision diagram" to encode sets of data.
Under some hypotheses, a fault tree can be converted into its equivalent state-space model like Markov chains or Stochastic Petri Nets. However, this approach has certain limitations. One is that the transformation can be made only if the model meets the needed requirements in the secondary modelling method, e.g., only components with exponential failures can be modeled using Markov chains. The other is the exponential state-space growth of the secondary models with a large number of components. Nevertheless, in a recent contribution [46], a simplified hierarchical approach shows that it is possible to convert a DFT by the mean of a limited state space constituted by five statuses and solve it through the Semi-Markov Process theory. In this way, the exclusive use of the exponential probability distributions is overcome. Unfortunately, the same approach is not much effective for RDFT.
Finally, the simulation approach based on Monte Carlo simulation (MCS) offers a modeling flexibility that overcomes the limitations of the previous methods. Many repeated iterations of the system evolutions are generated to achieve an acceptable accuracy in estimations. There are practical case studies that have used MCS to solve complex systems modeled by a fault tree. In [47], dynamic gates with repairs have been used to model an electrical power supply system of a nuclear power plant (NPP), and due to its complexity, MCS is used to solve this model. The main issues of the simulation are the computational speed and approximation of the results that can be unacceptable for real-time applications that require a fast response. Fortunately, this is not the case of offline dependability assessments whose results, generally, are included in the safety and risk report documentation of the enterprise. Recent contributions [48] propose relieving solutions for alleviating this issue and tackle rare events.
Although the category of "Advanced & General-Purpose Methods" is powerful, the main limitation of all these formalisms is that they are purely stochastic models. With these models, the boundary conditions must be fixed as input of the model, whereby the description of the stochastic events of the model gets limited to the usage of static probability density functions. Hybrid methods [49][50][51][52], falling under the umbrella of Dynamic Reliability, alleviate this issue, and, as shown in Figure 1, they embrace the scope covered by the other categories. Dynamic Reliability aims to relax the rigid hypotheses of traditional reliability, enabling the possibility to model multi-state systems and consider changes of the nominal design condition of a system. It can consider numerous characteristics of complex systems, such as inclusion of environmental dependencies, interactions between continuous process variables and system components, and stochastic and deterministic behaviours evolving in time. As matter of fact, a component does not always operate around the nominal design operative conditions, resulting in deviations of performance and wearing-out. The formulation and resolution of a Dynamic Reliability problem is not an easy task because, as shown in Equation (1), the mathematical relationship to write and solve it has to consider the stochastic events affecting the dynamic evolution of the working and operational conditions W (Equation (2)) and the variation of the physical dynamic behavior (Equation (3)), by means of a set of Partial Differential Algebraic Equations (PDAE).
where P is the probability that the time to failure of the system (S) is greater than the mission time T M and L is the aging variable. The working and operating conditions W can be described by the mean of a non-linear functional Γ that depends on time t, on the system dynamics X(t) and on a vector Λ of stochastic variables.
Finally, the general formulation of the PDAE problem, in terms of thermo-fluid, chemical, and rotational dynamic laws (e.g., mass energy balance, heat transfer equations) has the following form: where t denotes the time, y and u are the vectors of differential variables, algebraic variables and inputs, respectively. The vectors X and y are, respectively, the system and the process variables (state variables). F and B are vector functions. The partial derivative base vector, ζ, is an appropriate distribution domain, usually expressing geometry dimensions (e.g., length, width, radius) [53]. In order to simplify the modelling of such problems, recent studies have proposed methodologies based on the separation of concerns that allows the independent modelling of the deterministic and stochastic processes [54]. In this category of methods, the Stochastic Hybrid Fault Tree Automaton (SHyFTA) is the one that mostly resembles to the FTA [55]. A SHyFTA model is comprised of two interdependent models. The deterministic model and the stochastic model with the latter taking the form of a Repairable Dynamic Fault Tree. The formal definition states that a SHyFTA is a 13-uplet.
where t denotes the time, y and u are the vectors of differential variables, algebraic variables and inputs, respectively. The vectors X and y are, respectively, the system and the process variables (state variables). F and B are vector functions. The partial derivative base vector, ζ, is an appropriate distribution domain, usually expressing geometry dimensions (e.g., length, width, radius) [53]. In order to simplify the modelling of such problems, recent studies have proposed methodologies based on the separation of concerns that allows the independent modelling of the deterministic and stochastic processes [54]. In this category of methods, the Stochastic Hybrid Fault Tree Automaton (SHyFTA) is the one that mostly resembles to the FTA [55]. A SHyFTA model is comprised of two interdependent models. The deterministic model and the stochastic model with the latter taking the form of a Repairable Dynamic Fault Tree. The formal definition states that a SHyFTA is a 13-uplet. • is a finite set of arcs of the form ( , ɛj, Gk, ') where and ' are, respectively, the origin and the goal discrete states of the arc k, ɛj is the event associated with this arc, Gk is the guard condition on X in state ; • is a finite set of guard condition functions on each real variable Xi on the state j. • δ: × X → (ℝ → ℝ) is a function of "activities", describing the evolution of real variables in each discrete state; • is a finite set of clocks on ℝ that identify the firing of a deterministic or a stochastic event; is an application that associates a distribution function to the stochastic events ƐS, according to the clock H, the system evolution Χ and the discrete state ; • P is the instantaneous probability to be in i ∈ S; • GA is the finite set of gates of the fault tree model; • BE is the finite set of basic events of the fault tree model. The set BE contains a subset of novel conception, the Hybrid Basic Events (HBE), HBE ⊆ BE: ∀ ∈ HBE ∃ f : × × X → (ℝ → [0, 1]). This last characterization states that HBEs are those set of basic events whose failure distribution depends on the evolution of the system and vary with continuity in time; • T is the top event of the fault tree corresponding with the output of the main gate; • C is the set of connections between gates and basic events.
For a thorough description of the SHyFTA formalism, readers can refer to [55]. The Stochastic Hybrid Automaton is not easy to use and solve because it models, without constraints, any type of dependency. For this reason, the most appropriate resolution method is based on Monte Carlo simulation. In order to favour the design of a model, an ad-hoc software library, , X, Υ, δ, H, G, F, P, GA, BE, T, C

•
S is a finite set of discrete states {S D , S S }, where S D is the subset of deterministic states and S S of the stochastic one; • laws (e.g., mass energy balance, heat transfer equations) has the following form: e time, y and u are the vectors of differential variables, algebraic variables and . The vectors X and y are, respectively, the system and the process variables (state are vector functions. The partial derivative base vector, ζ, is an appropriate , usually expressing geometry dimensions (e.g., length, width, radius) [53]. simplify the modelling of such problems, recent studies have proposed ed on the separation of concerns that allows the independent modelling of the tochastic processes [54]. In this category of methods, the Stochastic Hybrid Fault HyFTA) is the one that mostly resembles to the FTA [55]. st characterization states that HBEs are those set of basic events whose failure pends on the evolution of the system and vary with continuity in time; ent of the fault tree corresponding with the output of the main gate; onnections between gates and basic events.
description of the SHyFTA formalism, readers can refer to [55]. Hybrid Automaton is not easy to use and solve because it models, without e of dependency. For this reason, the most appropriate resolution method is based ulation. In order to favour the design of a model, an ad-hoc software library, is a finite set of events { rotational dynamic laws (e.g., mass energy balance, heat transfer equations) has the following form: where t denotes the time, y and u are the vectors of differential variables, algebraic variables and inputs, respectively. The vectors X and y are, respectively, the system and the process variables (state variables). F and B are vector functions. The partial derivative base vector, ζ, is an appropriate distribution domain, usually expressing geometry dimensions (e.g., length, width, radius) [53]. In order to simplify the modelling of such problems, recent studies have proposed methodologies based on the separation of concerns that allows the independent modelling of the deterministic and stochastic processes [54]. In this category of methods, the Stochastic Hybrid Fault Tree Automaton (SHyFTA) is the one that mostly resembles to the FTA [55]. A SHyFTA model is comprised of two interdependent models. The deterministic model and the stochastic model with the latter taking the form of a Repairable Dynamic Fault Tree. The formal definition states that a SHyFTA is a 13-uplet. is a finite set of arcs of the form ( , ɛj, Gk, ') where and ' are, respectively, the origin and the goal discrete states of the arc k, ɛj is the event associated with this arc, Gk is the guard condition on X in state ; • is a finite set of guard condition functions on each real variable Xi on the state j. • δ: × X → (ℝ → ℝ) is a function of "activities", describing the evolution of real variables in each discrete state; • is a finite set of clocks on ℝ that identify the firing of a deterministic or a stochastic event; • F: × × X → (ℝ → [0,1]) is an application that associates a distribution function to the stochastic events ƐS, according to the clock H, the system evolution Χ and the discrete state ; • P is the instantaneous probability to be in i ∈ S; • GA is the finite set of gates of the fault tree model; • BE is the finite set of basic events of the fault tree model. The set BE contains a subset of novel conception, the Hybrid Basic Events (HBE), HBE ⊆ BE: ∀ ∈ HBE ∃ f : × × X → (ℝ → [0, 1]). This last characterization states that HBEs are those set of basic events whose failure distribution depends on the evolution of the system and vary with continuity in time; • T is the top event of the fault tree corresponding with the output of the main gate; • C is the set of connections between gates and basic events.
For a thorough description of the SHyFTA formalism, readers can refer to [55]. The Stochastic Hybrid Automaton is not easy to use and solve because it models, without constraints, any type of dependency. For this reason, the most appropriate resolution method is based on Monte Carlo simulation. In order to favour the design of a model, an ad-hoc software library, D , rotational dynamic laws (e.g., mass energy balance, heat transfer equations) has the following form: where t denotes the time, y and u are the vectors of differential variables, algebraic variables and inputs, respectively. The vectors X and y are, respectively, the system and the process variables (state variables). F and B are vector functions. The partial derivative base vector, ζ, is an appropriate distribution domain, usually expressing geometry dimensions (e.g., length, width, radius) [53]. In order to simplify the modelling of such problems, recent studies have proposed methodologies based on the separation of concerns that allows the independent modelling of the deterministic and stochastic processes [54]. In this category of methods, the Stochastic Hybrid Fault Tree Automaton (SHyFTA) is the one that mostly resembles to the FTA [55]. A SHyFTA model is comprised of two interdependent models. The deterministic model and the stochastic model with the latter taking the form of a Repairable Dynamic Fault Tree. The formal definition states that a SHyFTA is a 13-uplet. • is a finite set of arcs of the form ( , ɛj, Gk, ') where and ' are, respectively, the origin and the goal discrete states of the arc k, ɛj is the event associated with this arc, Gk is the guard condition on X in state ; • is a finite set of guard condition functions on each real variable Xi on the state j. • δ: × X → (ℝ → ℝ) is a function of "activities", describing the evolution of real variables in each discrete state; • is a finite set of clocks on ℝ that identify the firing of a deterministic or a stochastic event; • F: × × X → (ℝ → [0,1]) is an application that associates a distribution function to the stochastic events ƐS, according to the clock H, the system evolution Χ and the discrete state ; • P is the instantaneous probability to be in i ∈ S; • GA is the finite set of gates of the fault tree model; • BE is the finite set of basic events of the fault tree model. The set BE contains a subset of novel conception, the Hybrid Basic Events (HBE), HBE ⊆ BE: ∀ ∈ HBE ∃ f : × × X → (ℝ → [0, 1]). This last characterization states that HBEs are those set of basic events whose failure distribution depends on the evolution of the system and vary with continuity in time; • T is the top event of the fault tree corresponding with the output of the main gate; • C is the set of connections between gates and basic events.
For a thorough description of the SHyFTA formalism, readers can refer to [55]. The Stochastic Hybrid Automaton is not easy to use and solve because it models, without constraints, any type of dependency. For this reason, the most appropriate resolution method is based on Monte Carlo simulation. In order to favour the design of a model, an ad-hoc software library, S }, where rotational dynamic laws (e.g., mass energy balance, heat transfer equations) has the following form: where t denotes the time, y and u are the vectors of differential variables, algebraic variables and inputs, respectively. The vectors X and y are, respectively, the system and the process variables (state variables). F and B are vector functions. The partial derivative base vector, ζ, is an appropriate distribution domain, usually expressing geometry dimensions (e.g., length, width, radius) [53]. In order to simplify the modelling of such problems, recent studies have proposed methodologies based on the separation of concerns that allows the independent modelling of the deterministic and stochastic processes [54]. In this category of methods, the Stochastic Hybrid Fault Tree Automaton (SHyFTA) is the one that mostly resembles to the FTA [55]. A SHyFTA model is comprised of two interdependent models. The deterministic model and the stochastic model with the latter taking the form of a Repairable Dynamic Fault Tree. The formal definition states that a SHyFTA is a 13-uplet. • is a finite set of arcs of the form ( , ɛj, Gk, ') where and ' are, respectively, the origin and the goal discrete states of the arc k, ɛj is the event associated with this arc, Gk is the guard condition on X in state ; • is a finite set of guard condition functions on each real variable Xi on the state j. • δ: × X → (ℝ → ℝ) is a function of "activities", describing the evolution of real variables in each discrete state; • is a finite set of clocks on ℝ that identify the firing of a deterministic or a stochastic event; • F: × × X → (ℝ → [0,1]) is an application that associates a distribution function to the stochastic events ƐS, according to the clock H, the system evolution Χ and the discrete state ; • P is the instantaneous probability to be in i ∈ S; • GA is the finite set of gates of the fault tree model; • BE is the finite set of basic events of the fault tree model. The set BE contains a subset of novel conception, the Hybrid Basic Events (HBE), HBE ⊆ BE: ∀ ∈ HBE ∃ f : × × X → (ℝ → [0, 1]). This last characterization states that HBEs are those set of basic events whose failure distribution depends on the evolution of the system and vary with continuity in time; • T is the top event of the fault tree corresponding with the output of the main gate; • C is the set of connections between gates and basic events.
For a thorough description of the SHyFTA formalism, readers can refer to [55]. The Stochastic Hybrid Automaton is not easy to use and solve because it models, without constraints, any type of dependency. For this reason, the most appropriate resolution method is based on Monte Carlo simulation. In order to favour the design of a model, an ad-hoc software library, D is the subset of deterministic events and rotational dynamic laws (e.g., mass energy balance, heat transfer eq where t denotes the time, y and u are the vectors of differential v inputs, respectively. The vectors X and y are, respectively, the system variables). F and B are vector functions. The partial derivative b distribution domain, usually expressing geometry dimensions (e.g., In order to simplify the modelling of such problems, methodologies based on the separation of concerns that allows th deterministic and stochastic processes [54]. In this category of meth Tree Automaton (SHyFTA) is the one that mostly resembles to the comprised of two interdependent models. The deterministic model a latter taking the form of a Repairable Dynamic Fault Tree. The forma is a 13-uplet. Υ is a finite set of arcs of the form (s, ε j , G k , s') where s and s' are, respectively, the origin and the goal discrete states of the arc k, ε j is the event associated with this arc, G k is the guard condition on X in state s; • G is a finite set of guard condition functions on each real variable X i on the state s j.
is a function of "activities", describing the evolution of real variables in each discrete state; • H is a finite set of clocks on R that identify the firing of a deterministic or a stochastic event; is an application that associates a distribution function to the stochastic events ormation 2019, 10, x FOR PEER REVIEW 6 of 39 here ℙ is the probability that the time to failure of the system (S) is greater than the mission time TM d L is the aging variable. The working and operating conditions W can be described by the mean a non-linear functional  that depends on time t, on the system dynamics X(t) and on a vector  of ochastic variables.
Finally, the general formulation of the PDAE problem, in terms of thermo-fluid, chemical, and tational dynamic laws (e.g., mass energy balance, heat transfer equations) has the following form: here t denotes the time, y and u are the vectors of differential variables, algebraic variables and puts, respectively. The vectors X and y are, respectively, the system and the process variables (state riables). F and B are vector functions. The partial derivative base vector, ζ, is an appropriate stribution domain, usually expressing geometry dimensions (e.g., length, width, radius) [53]. In order to simplify the modelling of such problems, recent studies have proposed ethodologies based on the separation of concerns that allows the independent modelling of the terministic and stochastic processes [54]. In this category of methods, the Stochastic Hybrid Fault ee Automaton (SHyFTA) is the one that mostly resembles to the FTA [55]. Gk, ') where and ' are, respectively, the origin and the goal discrete states of the arc k, ɛj is the event associated with this arc, Gk is the guard condition on X in state ; is a finite set of guard condition functions on each real variable Xi on the state j. δ: × X → (ℝ → ℝ) is a function of "activities", describing the evolution of real variables in each discrete state; is a finite set of clocks on ℝ that identify the firing of a deterministic or a stochastic event; F: × × X → (ℝ → [0,1]) is an application that associates a distribution function to the stochastic events ƐS, according to the clock H, the system evolution Χ and the discrete state ; P is the instantaneous probability to be in i ∈ S; GA is the finite set of gates of the fault tree model; BE is the finite set of basic events of the fault tree model. The set BE contains a subset of novel conception, the Hybrid Basic Events (HBE), HBE ⊆ BE: ∀ ∈ HBE ∃ f : × × X → (ℝ → [0, 1]). This last characterization states that HBEs are those set of basic events whose failure distribution depends on the evolution of the system and vary with continuity in time; T is the top event of the fault tree corresponding with the output of the main gate; C is the set of connections between gates and basic events.
For a thorough description of the SHyFTA formalism, readers can refer to [55]. The Stochastic Hybrid Automaton is not easy to use and solve because it models, without nstraints, any type of dependency. For this reason, the most appropriate resolution method is based Monte Carlo simulation. In order to favour the design of a model, an ad-hoc software library, S, according to the clock H, the system evolution X and the discrete state S; • P is the instantaneous probability to be in s i ∈ S S ; • GA is the finite set of gates of the fault tree model; • BE is the finite set of basic events of the fault tree model. The set BE contains a subset of novel conception, the Hybrid Basic Events (HBE), HBE ⊆ BE : . This last characterization states that HBEs are those set of basic events whose failure distribution depends on the evolution of the system and vary with continuity in time; • T is the top event of the fault tree corresponding with the output of the main gate; • C is the set of connections between gates and basic events.
For a thorough description of the SHyFTA formalism, readers can refer to [55]. The Stochastic Hybrid Automaton is not easy to use and solve because it models, without constraints, any type of dependency. For this reason, the most appropriate resolution method is based on Monte Carlo simulation. In order to favour the design of a model, an ad-hoc software library, named SHyFTOO, has been specifically developed under the Matlab ® framework [22]. SHyFTOO is characterized by a computer algorithm that combines Time-Driven and Discrete-Event simulation [22,55,56]. The main disadvantage of this library is that it can require numerous iterations and long runs to retrieve accurate results. However, the authors are confident that the speed of penetration of new paradigms based on parallel and cloud computing [57] can alleviate this issue. Moreover, the deepness of a model [4], meant as the capability to describe the accuracy of a phenomena within a correct sizing of the boundary conditions (and to avoid too complex models), must be considered in order to design a model that is effective enough for the type of analysis to achieve. In fact, while traditional models are based on fixed probability density functions, in the Dynamic Reliability, the specific experience and knowledge of the system equipment and of the processes by the modelers is an important factor that affects the accuracy of the final model (e.g., the relationships between components-failure dynamics and the system working conditions). The proposed version of SHyFTOO has been coded in Matlab ® , which offers a high-level user-friendly interface and facilitates the modelling of dynamic systems.
In the Appendix A, all the main components of the SHyFTOO library are described. Moreover, the novel standard template Simulink blocks that belong to the SHyFTOO library package and that allow a DFT to couple with a Simulink dynamic system of Matlab ® are presented.

Case Study
This section describes the main case study of this manuscript and the SHyFTA model implemented with the SHyFTOO library. Its resolution is then presented in Section 4 together with other case study results. To better understand the process under analysis, Figure 2 shows a representation of the SHyFTA model hereby described. As it will be shown, it is implemented by using several built-in Simulink blocks that are combined with the standard template blocks of the SHyFTOO library discussed in the Appendix A.
correct sizing of the boundary conditions (and to avoid too complex models), must be considered in order to design a model that is effective enough for the type of analysis to achieve. In fact, while traditional models are based on fixed probability density functions, in the Dynamic Reliability, the specific experience and knowledge of the system equipment and of the processes by the modelers is an important factor that affects the accuracy of the final model (e.g., the relationships between components-failure dynamics and the system working conditions). The proposed version of SHyFTOO has been coded in Matlab ® , which offers a high-level user-friendly interface and facilitates the modelling of dynamic systems.
In the Appendix A, all the main components of the SHyFTOO library are described. Moreover, the novel standard template Simulink blocks that belong to the SHyFTOO library package and that allow a DFT to couple with a Simulink dynamic system of Matlab ® are presented.

Case Study
This section describes the main case study of this manuscript and the SHyFTA model implemented with the SHyFTOO library. Its resolution is then presented in Section 4 together with other case study results. To better understand the process under analysis, Figure 2 shows a representation of the SHyFTA model hereby described. As it will be shown, it is implemented by using several built-in Simulink blocks that are combined with the standard template blocks of the SHyFTOO library discussed in the Appendix A. This system is constituted by a distillation column, working continuously 24 h per 365 days a year, that has to purify the output mixture of an industrial process, separating the solid part from the liquid. The undesired top event corresponds with the unavailability of the system, which involves the stopping of the purifying process. The input section of the distillation column (identified with the subsystem IN) is constituted by two electric pumps (HBE1 and HBE2) and by a valve (BE3) in charge of pumping the mixture to purify to the distillation column. The electric pumps are in a standby configuration; therefore, they are modelled with a SPARE gate (SP1). The suction system (OUT) of the distillation column is made up of a valve (BE4) and by an electric pump (HBE5). Moreover, it is assumed that the tank containing the solid part needs to be cleaned when it is full. This event is This system is constituted by a distillation column, working continuously 24 h per 365 days a year, that has to purify the output mixture of an industrial process, separating the solid part from the liquid. The undesired top event corresponds with the unavailability of the system, which involves the stopping of the purifying process. The input section of the distillation column (identified with the subsystem IN) is constituted by two electric pumps (HBE1 and HBE2) and by a valve (BE3) in charge of pumping the mixture to purify to the distillation column. The electric pumps are in a standby configuration; therefore, they are modelled with a SPARE gate (SP1). The suction system (OUT) of the distillation column is made up of a valve (BE4) and by an electric pump (HBE5). Moreover, it is assumed that the tank containing the solid part needs to be cleaned when it is full. This event is represented by the hybrid basic event HBE6. As shown in the fault tree, the top event is verified if any among the IN, OUT or HBE6 occurs. For this reason, an OR gate (TOP) is used.
The electric pumps (HBE1, HBE2, and HBE5), modelled as hybrid basic events, are characterized by a dynamic failure rate according to the formula of the MIL-HBDK-217 standard [58]. This standard has the great merit to identify a mathematical relationship between the failure rate and the aging L. Therefore, even the correctness of this relationship cannot be guaranteed by the authors. It represents a valuable example for illustrating how to use the SHyFTA methodology.
where α B is the characteristic life of the bearings and α W is the characteristic life of the coils. Instead, the hybrid event HBE6 depends on the physical characteristics of the solid tank of the distillation column and on the physical characteristics of the input mixture to separate. In fact, the volume of the solid part can be computed by the following formula: where S(τ) is the fraction of solid part contained in the flow mixture Q(τ), input of the distillation column. When the volume of the solid part deposited in the tank reaches the Vs threshold of the tank (V ST ), Vs = V ST , the distillation process must be stopped to perform the operation of cleaning of the solid tank of the distillation column.
Moreover, a quality requirement of the cleaning process specifies that the liquid part achieves the standard purity if, after a stop of the system, the process restarts within 16 hours. In other words, the mixture inside the distillation column cannot lie idle longer than 16 hours, otherwise it must be dumped. Finally, at the end of the year, the company has to determine the amount of mixture that could not be purified. If it exceeds a certain threshold, the safety of the environment is compromised, and an economic penalty has to be paid.
The process described so far contains many elements that cannot be modelled with any of the modelling "Traditional Static" or the "Advanced & General-Purpose" methods. In fact, the dynamic failure rate specified in Equation (4), the condition for triggering the HBE6 in Equation (5) and the constraints concerning the quality requirement and the safety of the environment have to be modelled considering the dynamics of the physical process. A SHyFTA model of the presented case study is developed using the SHyFTOO, and the next subsections contains the main information for implementing it. Several information that will be presented refer to the main components of the SHyFTOO library that are discussed, thoroughly, in the Appendix A. To improve the reading of the manuscript, throughout the manuscript, the elements of the SHyFTOO library that are part of the Appendix A are referred with the following notation: << component >>.

Definition of the Fault Tree
The definition of a DFT requires the customization of the <<initFaultTree.m>> script. As is shown in Algorithm 1, the code for the specification of the fault tree model is intuitive. Line 2 defines the mission time Tm, whereas from line 5 to 18 the fault tree components are set. Initially, the basic events to be used as parameters of the gate definition are defined. Line 18 is very important because the SHyFTOO library expects to recognize a Top Event by the use of a variable TOP that is equals to the top event gate defined by the user.

Definition of the Physical Process (with Simulink)
As already anticipated, the Simulink toolbox simplifies the design of complex dynamic systems. Thanks to the SHyFTOO library, it can be easily coupled with the stochastic process by accessing the parameters of the DFT components that must be incorporated within the physical process (e.g., the status of a BE). The Simulink implementation of the physical process for the case study proposed is shown in Figure 3. Line 2 defines the mission time Tm, whereas from line 5 to 18 the fault tree components are set. Initially, the basic events to be used as parameters of the gate definition are defined. Line 18 is very important because the SHyFTOO library expects to recognize a Top Event by the use of a variable TOP that is equals to the top event gate defined by the user.

Definition of the Physical Process (with Simulink)
As already anticipated, the Simulink toolbox simplifies the design of complex dynamic systems. Thanks to the SHyFTOO library, it can be easily coupled with the stochastic process by accessing the parameters of the DFT components that must be incorporated within the physical process (e.g., the status of a BE). The Simulink implementation of the physical process for the case study proposed is shown in Figure 3. Among the displayed blocks, blocks 1 and block 2 belong to the standard built-in blocks of the SHyFTOO Simulink (presented in the Appendix A). They are respectively the <<ITER EVOLUTION >> and <<RACE CONDITION >> and they must be included without any customization in the Simulink implementation of the physical process. Blocks 3, 4, and 5 inherit from the <<GENERIC HYBRID BASIC EVENT>> block, and only a few modifications are needed to adapt the internal <<HBE_SAMPLING>> block to the characteristics of the probability density function of the corresponding events. Block 6 models the physical process that determines the triggering of the HBE6. In fact, it depends on the condition Vs = VST, knowing that Vs follows Equation (5). Finally, block 7 models the quality requirement of the cleaning process, which is not a dependability attribute but a physical property that depends on the availability of the distillation column. Among the displayed blocks, blocks 1 and block 2 belong to the standard built-in blocks of the SHyFTOO Simulink (presented in the Appendix A). They are respectively the <<ITER EVOLUTION >> and <<RACE CONDITION >> and they must be included without any customization in the Simulink implementation of the physical process. Blocks 3, 4, and 5 inherit from the <<GENERIC HYBRID BASIC EVENT>> block, and only a few modifications are needed to adapt the internal <<HBE_SAMPLING>> block to the characteristics of the probability density function of the corresponding events. Block 6 models the physical process that determines the triggering of the HBE6. In fact, it depends on the condition Vs = V ST , knowing that Vs follows Equation (5). Finally, block 7 models the quality requirement of the cleaning process, which is not a dependability attribute but a physical property that depends on the availability of the distillation column.
To better understand the custom modification of all these blocks, in the next subsections, the Simulink models are discussed into more detail.

Customization of the HBE1 (and HBE5) Block
The Simulink ensemble shown in Figure 4 is devoted to the failure sampling of the hybrid basic events HBE1, identified with the label '4 , and with the same logic to the HBE5, identified with the label '5 .

Customization of the HBE1 (and HBE5) Block
The Simulink ensemble shown in Figure 4 is devoted to the failure sampling of the hybrid basic events HBE1, identified with the label '4′, and with the same logic to the HBE5, identified with the label '5′. The block HBE1_SAMPLING contains the logic to simulate, at any timestep, the failure of the electric pump according to the dynamic failure rate specified in Equation (4). This equation depends on the aging of the electric pump, and one of the inputs to provide to the "Interpreted Matlab function" is the aging variable. As shown in Figure 4, the aging is computed using an integrator block that takes as input a Boolean signal (converted into a double value), the output of a complex logic that considers the status of the electric pump itself (<<HBE1.Status>>), and the status of the Top Event. If both are good, it means that the pump is operating correctly, so its aging grows linearly. Clearly, if the pump has failed, the aging can no longer increase. With the same reasoning, the aging does not increase even if the Top Event has occurred because it is assumed that, if the cleaning process has paused all, the components stop operating until the system is restored. The integrator block takes another input as a reset condition so as to reset the value of the aging to 0 when the pump is restored as new.

Customization of the HBE2 Block
The Simulink ensemble shown in Figure 5 models the failure sampling of the hybrid basic events HBE2 identified with the label '3′. The block HBE1_SAMPLING contains the logic to simulate, at any timestep, the failure of the electric pump according to the dynamic failure rate specified in Equation (4). This equation depends on the aging of the electric pump, and one of the inputs to provide to the "Interpreted Matlab function" is the aging variable. As shown in Figure 4, the aging is computed using an integrator block that takes as input a Boolean signal (converted into a double value), the output of a complex logic that considers the status of the electric pump itself (<<HBE1.Status>>), and the status of the Top Event. If both are good, it means that the pump is operating correctly, so its aging grows linearly. Clearly, if the pump has failed, the aging can no longer increase. With the same reasoning, the aging does not increase even if the Top Event has occurred because it is assumed that, if the cleaning process has paused all, the components stop operating until the system is restored. The integrator block takes another input as a reset condition so as to reset the value of the aging to 0 when the pump is restored as new.

Customization of the HBE2 Block
The Simulink ensemble shown in Figure 5 models the failure sampling of the hybrid basic events HBE2 identified with the label '3 .

Customization of the HBE6 Block
The Simulink ensemble shown in Figure 6 is devoted to the failure sampling of the hybrid basic events HBE6, identified with the label '6 .

Customization of the HBE6 Block
The Simulink ensemble shown in Figure 6 is devoted to the failure sampling of the hybrid basic events HBE6, identified with the label '6′. The blocks "% particle" and "QI" represent respectively the two physical variables constituting the input mixture to the process, where it is assumed that the former is the percentage of the solid part inside the infinitesimal volume QI, constituting the instantaneous flowrate of mass processed in the distillation column. Those are random variables that vary at any timestep of the simulation with a uniform probability density function distribution. It was assumed that the percentage of solid part can take at most 2% of the instantaneous QI mass. The QI vary in the range [400, 600] m 3 /h. The blocks "% particle" and "QI" represent respectively the two physical variables constituting the input mixture to the process, where it is assumed that the former is the percentage of the solid part inside the infinitesimal volume QI, constituting the instantaneous flowrate of mass processed in the distillation column. Those are random variables that vary at any timestep of the simulation with a uniform probability density function distribution. It was assumed that the percentage of solid part can take at most 2% of the instantaneous QI mass. The QI vary in the range [400, 600] m 3 /h. Back to Figure 6, it is possible to see that these two variables are multiplied and integrated so as to implement the formula of Equation (2). The product block takes as input also the <<TOP.Status>> because, if the system is unavailable, it is assumed to stop the cleaning process with a loss in terms of mixture mass, which never stops flowing (in that case it is lost and must be recorded to evaluate the environmental damage). Finally, the HBE6_SAMPLING block, shown in Figure 7, modifies a bit the template model of Figure A2 (see Appendix A) as the event modelling the filling of the solid tank is not governed by a stochastic probability density function, but by a physical condition. Therefore, the "Interpreted Matlab function" to use in this block has been codified ad-hoc, as shown in Algorithm 2.
of mixture mass, which never stops flowing (in that case it is lost and must be recorded to evaluate the environmental damage). Finally, the HBE6_SAMPLING block, shown in Figure 7, modifies a bit the template model of Figure A2 (see Appendix A) as the event modelling the filling of the solid tank is not governed by a stochastic probability density function, but by a physical condition. Therefore, the "Interpreted Matlab function" to use in this block has been codified ad-hoc, as shown in Algorithm 2. Line 2 is a Matlab ® API that enables the assignment and the evaluation of the workspace variables (e.g., the global variables of the main script <<SHyFTAmain.m>>) within a Matlab ® function. For instance, in line 5, the variable FTA is initialized with the value of the workspace variable FTA that contains all the element of the stochastic model, including basic events and gates. In line 6, the variable HBEid is initialized with the identifier of the hybrid basic events whose name is contained in the input u (2), passed to the function. Line 7 initializes the local variable decanterCapacity with the value set in the script <<SHyFTAmain.m>> and it represents the capacity VST of the solid tank. From line 9 to 13, the function performs the comparison of the current volume at time t passed as input u (1) and the value of the solid tank (decanterCapacity). If the inequality is verified, the function sets the property <<FailureTime>> of the hybrid basic event with the value of the current simulation time, currentTime, that was retrieved in line 4 using the Matlab ® API get_param (). Afterward, in line 12, the global variable FTA is set with the local variable FTA that contains the change on the hybrid basic event declared as failed. In line 13 the output of the function is set equal to 0. This value is returned and used in the block assertion <<HBEevaluateFT>> that pausing the iteration performs the evaluation Line 2 is a Matlab ® API that enables the assignment and the evaluation of the workspace variables (e.g., the global variables of the main script <<SHyFTAmain.m>>) within a Matlab ® function. For instance, in line 5, the variable FTA is initialized with the value of the workspace variable FTA that contains all the element of the stochastic model, including basic events and gates. In line 6, the variable HBEid is initialized with the identifier of the hybrid basic events whose name is contained in the input u (2), passed to the function. Line 7 initializes the local variable decanterCapacity with the value set in the script <<SHyFTAmain.m>> and it represents the capacity V ST of the solid tank. From line 9 to 13, the function performs the comparison of the current volume at time t passed as input u (1) and the value of the solid tank (decanterCapacity). If the inequality is verified, the function sets the property <<FailureTime>> of the hybrid basic event with the value of the current simulation time, currentTime, that was retrieved in line 4 using the Matlab ® API get_param (). Afterward, in line 12, the global variable FTA is set with the local variable FTA that contains the change on the hybrid basic event declared as failed. In line 13 the output of the function is set equal to 0. This value is returned and used in the block assertion <<HBEevaluateFT>> that pausing the iteration performs the evaluation of the stochastic model. If the returning value y is equal to 1, the assertion block does not raise, and the iteration can continue without interruptions.

Customization of the Physical Process
The Simulink ensemble shown in Figure 8 models the part of the physical process devoted to the measure of the quality requirement of the cleaned liquid. the iteration can continue without interruptions.

Customization of the Physical Process
The Simulink ensemble shown in Figure 8 models the part of the physical process devoted to the measure of the quality requirement of the cleaned liquid. The output variable QI_Mixture records the flow of mixture which is correctly processed. It is obtained multiplying the input QI with the <<Top.Status>>. In other words, the variable QI_Mixture will record the value 0 for those instants of time in which the system is not working (i.e., namely it is unavailable). On the other hand, with the same logic, the output variable QI_Lost records the flow of mixture which cannot be processed during the system unavailability. This value is obtained multiplying the input variable QI with the complement to one of the <<Top.Status>> (= "⌐TOP EVENT STATUS").
The other part of the ensemble implements the logic to verify the quality requirement of the liquid part, separated in the distillation column. It makes use of the <<Top.Status>>. Any change from bad to good and vice versa triggers an assertion (SHA_TOP) that verifies if the amount of time from the last stop (due to the unavailability of the system) exceeds the threshold dictated by the quality process variable time4GoodProcess. The following code, shown in Algorithm 3, is used. At first, it can be noticed that, in Matlab ® , the script has complete visibility of the workspace global variables and therefore it is not required any coder.extrinsic API to access them. Line 2 checks if the Top Event has occurred and assigns to the global variable lastStop the value of the current time, which has been retrieved using the get_param API of Matlab ® .

Algorithm 3
Code for the definition of the quality requirement of the case study.  The output variable QI_Mixture records the flow of mixture which is correctly processed. It is obtained multiplying the input QI with the <<Top.Status>>. In other words, the variable QI_Mixture will record the value 0 for those instants of time in which the system is not working (i.e., namely it is unavailable). On the other hand, with the same logic, the output variable QI_Lost records the flow of mixture which cannot be processed during the system unavailability. This value is obtained multiplying the input variable QI with the complement to one of the <<Top.Status>> (= "

TOP EVENT STATUS").
The other part of the ensemble implements the logic to verify the quality requirement of the liquid part, separated in the distillation column. It makes use of the <<Top.Status>>. Any change from bad to good and vice versa triggers an assertion (SHA_TOP) that verifies if the amount of time from the last stop (due to the unavailability of the system) exceeds the threshold dictated by the quality process variable time4GoodProcess. The following code, shown in Algorithm 3, is used. At first, it can be noticed that, in Matlab ® , the script has complete visibility of the workspace global variables and therefore it is not required any coder.extrinsic API to access them. Line 2 checks if the Top Event has occurred and assigns to the global variable lastStop the value of the current time, which has been retrieved using the get_param API of Matlab ® .  Otherwise, when the system has recovered, the algorithm checks the time needed to restart the process, obtained as the difference between the current time and the last stop variable. If this value is longer than the constant process variable time4GoodProcess (that specifies the max amount of time that the mixture can be idle to remain good) then the global variable volumeNok-keeping track of the volume of mixture so far dumped-has to be updated.
In order to record the other relevant physical, the script <<shyftaMetrics.m>> is used. Since, the <<shyftaMetrics.m>> is launched at the end of each iteration, it is possible to compute the volume of mixture cleaned and the volume of mixture that could not be processed because of the unavailability of the system using the Matlab ® function sum, passing as parameter the arrays of the QI_Mixture and QI_Lost set in the Simulink block shown in Figure 8.

Case Studies Resolution
This section contains the results of several case studies, including the one described in Section 3, selected to demonstrate the applicability of the SHyFTA methodology to a wide range of industrial applications. The resolution of these case studies has been carried out with the SHyFTOO library, but the corresponding models are not included for space limits. Figure 9 shows the simplified SHyFTA model of an electric motor. This system consists of several components where, among all, it is possible to identify the BE1 and BE3 (bearings), and BE2 (winding or coil), characterized by the dynamic failure rates of Equations (6) and (7).

Electric Motor
In Equations (6) and (7), L(t) represents the aging of the component, and T A is the ambient temperature. The parameters of these equations are taken from the standard MIL-HDBK-217 [58]. This standard has the great merit to identify a mathematical relationship between the failure rates and the real-time variables of the aging and of the ambient temperature; even the correctness of this relationship cannot be guaranteed by the authors, it represents a valuable example for illustrating how to use the SHyFTA methodology. Unfortunately, we could not find a different data source that could have allowed us to perform a double check of the final results.
The failure of one of these components causes the failure of the entire system. To account for further types of failure, a generic BE4 event has been included in the fault tree model. As it is possible to see in the Figure 9, BE4 is the only basic event modelled as a traditional basic event. The main goal is to assess the reliability of the electric motor and compare it with the one obtained using a traditional FTA model. Due to its limitations, the static FTA cannot take as input a dynamic failure rate and therefore the parameters shown in Table 2 have been fixed with constant values using the method of computation described in the [58]. Moreover, for λBE1 and λBE2, the worstcase failure rates (with L = Tm = 8760 h) have been used. The main goal is to assess the reliability of the electric motor and compare it with the one obtained using a traditional FTA model. Due to its limitations, the static FTA cannot take as input a dynamic failure rate and therefore the parameters shown in Table 2 have been fixed with constant values using the method of computation described in the [58]. Moreover, for λ BE1 and λ BE2 , the worst-case failure rates (with L = T m = 8760 h) have been used. Table 2. Failure and repair rates of the Fault Tree model of the electric motor case study.

Event
Failure Rate (MIL-HDBK-217) 1.00 × 10 −6 (h −1 ) Figure 10 shows the unreliability results for the two models-the Fault Tree and the SHyFTA, respectively. It is possible to notice that the unreliability computed with the Fault Tree is higher than the one simulated with the SHyFTA model. The former increases linearly whereas the latter presents a pronounced non-linear behaviour after the 7000th hour from the beginning of the mission. The main goal is to assess the reliability of the electric motor and compare it with the one obtained using a traditional FTA model. Due to its limitations, the static FTA cannot take as input a dynamic failure rate and therefore the parameters shown in Table 2 have been fixed with constant values using the method of computation described in the [58]. Moreover, for λBE1 and λBE2, the worstcase failure rates (with L = Tm = 8760 h) have been used.

Event
Failure Rate (MIL-HDBK-217) BE1/BE2 5.98 × 10 −6 (h −1 ) BE3 5.89 × 10 −7 (h −1 ) BE4 1.00 × 10 −6 (h −1 ) Figure 10 shows the unreliability results for the two models-the Fault Tree and the SHyFTA, respectively. It is possible to notice that the unreliability computed with the Fault Tree is higher than the one simulated with the SHyFTA model. The former increases linearly whereas the latter presents a pronounced non-linear behaviour after the 7000th hour from the beginning of the mission. To better understand the effect of the bearing over the unreliability of motor, it is important to recall that its failure rate increases with the square of the aging L and decreases with the cube of αB To better understand the effect of the bearing over the unreliability of motor, it is important to recall that its failure rate increases with the square of the aging L and decreases with the cube of α B as shown in Equation (6). On the other hand, α B grows with the increasing of the ambient temperature whose trend is shown in Figure 11. The combination of these effects on the dynamic failure rate of the bearing can be analysed by observing the trend in Figure 12. At the beginning of the mission, the failure rate is low, although around the 1000th hour, it presents an interesting peak of about 3.5 × 10 −5 caused by a sudden decrease of the ambient temperature in wintertime.
as shown in Equation (6). On the other hand, αB grows with the increasing of the ambient temperature whose trend is shown in Figure 11. The combination of these effects on the dynamic failure rate of the bearing can be analysed by observing the trend in Figure 12. At the beginning of the mission, the failure rate is low, although around the 1000th hour, it presents an interesting peak of about 3.5 × 10 −5 caused by a sudden decrease of the ambient temperature in wintertime.  During the hot season (indicatively from the 1 June-1 October, in the interval comprised between the 3624th and the 6552th h) the failure rate is constant because of the prevalence of the parameter αB over the aging L. However, when the ambient temperature starts decreasing again (after as shown in Equation (6). On the other hand, αB grows with the increasing of the ambient temperature whose trend is shown in Figure 11. The combination of these effects on the dynamic failure rate of the bearing can be analysed by observing the trend in Figure 12. At the beginning of the mission, the failure rate is low, although around the 1000th hour, it presents an interesting peak of about 3.5 × 10 −5 caused by a sudden decrease of the ambient temperature in wintertime.  During the hot season (indicatively from the 1 June-1 October, in the interval comprised between the 3624th and the 6552th h) the failure rate is constant because of the prevalence of the parameter αB over the aging L. However, when the ambient temperature starts decreasing again (after During the hot season (indicatively from the 1 June-1 October, in the interval comprised between the 3624th and the 6552th h) the failure rate is constant because of the prevalence of the parameter α B over the aging L. However, when the ambient temperature starts decreasing again (after 1 October), the aging effect begins to dominate and the failure rate presents several peaks that are revealed during the night-time hours, when the ambient temperature decreases. This explains the important spike of the unreliability of the motor, after the 7000th h. These results demonstrate the huge difference between a static model and a hybrid dynamic reliability model, with this latter able to model and capture non-trivial physics-based dynamic operation features of a real electro-mechanical process.

Domestic Photovoltaic Power Plant with Storage System
In this case study [59], a SHyFTA model of a domestic grid-connected photovoltaic power plant has been simulated in order to compare the performance and the economic benefits of the same configuration, with and without a storage system. The main technical information of the system is shown in Table 3. The single-line diagram of the grid-connected photovoltaic power plant is shown in Figure 13. The SHyFTA model has to simulate the process of energy supply for a generic household equipped with a grid-connected photovoltaic power plant and a storage system. The single-line diagram, in Figure 13, allows the identification of the main sub-systems: the photovoltaic power plant (PV Generator), the storage system (BAT), and the equipment of the grid connection coupling (GCC) that allow the coupling with the electrical grid. 1 October), the aging effect begins to dominate and the failure rate presents several peaks that are revealed during the night-time hours, when the ambient temperature decreases. This explains the important spike of the unreliability of the motor, after the 7000th h. These results demonstrate the huge difference between a static model and a hybrid dynamic reliability model, with this latter able to model and capture non-trivial physics-based dynamic operation features of a real electromechanical process.

Domestic Photovoltaic Power Plant with Storage System
In this case study [59], a SHyFTA model of a domestic grid-connected photovoltaic power plant has been simulated in order to compare the performance and the economic benefits of the same configuration, with and without a storage system. The main technical information of the system is shown in Table 3. The single-line diagram of the grid-connected photovoltaic power plant is shown in Figure 13. The SHyFTA model has to simulate the process of energy supply for a generic household equipped with a grid-connected photovoltaic power plant and a storage system. The single-line diagram, in Figure 13, allows the identification of the main sub-systems: the photovoltaic power plant (PV Generator), the storage system (BAT), and the equipment of the grid connection coupling (GCC) that allow the coupling with the electrical grid.  It is possible to identify the following main sub-sections: 1.
Direct Current Section (DCS), made up of string protection diodes (SPR), a DC disconnector (DCD), and a surge protection device (SPD); 3.
Alternating Current Section (ACS), made up of an inverter (INV) and an AC circuit breaker (ACB);

4.
Grid Connector Coupling (GCC), made up of an AC disconnector (ACD), a differential circuit breaker (DCB), and a generic sub-system representing the electrical grid (GRD).

5.
Battery (BAT) that is connected in parallel in the AC section.
The technical regulation of the grid-connected photovoltaic [60] power plant states that if the electrical grid fails, the power plant must be disconnected, stopping the production and the energy supply of the household. In this case, and any time the photovoltaic power plant is disconnected, the battery is also forbidden to supply energy to the household. This scenario corresponds with the top event of the DFT of Figure 14. The DFT is constituted by the Top Event AND gate (TE) that takes as input the output of the OR gate GCC (OR (ACD, DCB, GRD)) and the OR gate "PV Down" (OR (ACS, DCS)). The former models any type of disconnection of the electrical grid, whereas the latter the unavailability of the photovoltaic power plant that occurs if the electrical circuit of the PV Generator gets open (any failure of the ACS or DCS components). The AND Gate PVM models the failure of the photovoltaic strings; although the modules are connected in series, the by-pass diodes guarantee the electrical isolation of those modules that are not working properly. As far as it concerns the battery, it must be pointed out that its unavailability does not cause a stop of the household energy supply because the grid can fulfill the energy request. On the other hand, the battery is disconnected if the PV power plant is down. These behaviors can be modelled with two FDEP Gates. The first takes as primary input the output of the GCC and as secondary input the OR Gate "PV down" [59]. The second FDEP takes as primary input the "PV Down" output and as second input the gate "BAT down". The failure and repair rates of the DFT model are shown in Table 4.
It is possible to identify the following main sub-sections: 1. PV Module (PVM), made up by ten photovoltaic modules (M1-M10); 2. Direct Current Section (DCS), made up of string protection diodes (SPR), a DC disconnector (DCD), and a surge protection device (SPD); 3. Alternating Current Section (ACS), made up of an inverter (INV) and an AC circuit breaker (ACB); 4. Grid Connector Coupling (GCC), made up of an AC disconnector (ACD), a differential circuit breaker (DCB), and a generic sub-system representing the electrical grid (GRD).

Battery (BAT) that is connected in parallel in the AC section.
The technical regulation of the grid-connected photovoltaic [60] power plant states that if the electrical grid fails, the power plant must be disconnected, stopping the production and the energy supply of the household. In this case, and any time the photovoltaic power plant is disconnected, the battery is also forbidden to supply energy to the household. This scenario corresponds with the top event of the DFT of Figure 14. The DFT is constituted by the Top Event AND gate (TE) that takes as input the output of the OR gate GCC (OR (ACD, DCB, GRD)) and the OR gate "PV Down" (OR (ACS, DCS)). The former models any type of disconnection of the electrical grid, whereas the latter the unavailability of the photovoltaic power plant that occurs if the electrical circuit of the PV Generator gets open (any failure of the ACS or DCS components). The AND Gate PVM models the failure of the photovoltaic strings; although the modules are connected in series, the by-pass diodes guarantee the electrical isolation of those modules that are not working properly. As far as it concerns the battery, it must be pointed out that its unavailability does not cause a stop of the household energy supply because the grid can fulfill the energy request. On the other hand, the battery is disconnected if the PV power plant is down. These behaviors can be modelled with two FDEP Gates. The first takes as primary input the output of the GCC and as secondary input the OR Gate "PV down" [59]. The second FDEP takes as primary input the "PV Down" output and as second input the gate "BAT down". The failure and repair rates of the DFT model are shown in Table 4. The inverter and the battery are the most sensitive components to the aging. Therefore, they have been modelled with a dynamic failure rate that allows to consider the wearing-out. To this end, a Weibull pdf with shape factor β > 1 (i.e., the failure rate is increasing with respect to time) [54] and The inverter and the battery are the most sensitive components to the aging. Therefore, they have been modelled with a dynamic failure rate that allows to consider the wearing-out. To this end, a Weibull pdf with shape factor β > 1 (i.e., the failure rate is increasing with respect to time) [54] and a scale parameter γ (that is generally set with the corresponding constant failure rate of an exponential distribution) has been used. Equation (8) defines the failure rate as   (5), β = 1.5, γ = 1.14 × 10 5 2.1 × 10 −3

BAT Energy Shortage
This event is modelled in the physical process. See [59] for more info.
The economic assessment has been carried out using the cash-flow method along the 20 years of operations, considering the main characteristics of the Italian Market. In order to compute the net present value (NPV) of the investments, it was required that we evaluate the service availability. This attribute of the dependability assessment depends on the failure/repair behaviour described by the DFT of Figure 14 and also on the capability of the system to provide the expected service when it is needed. This latter corresponds with the energy provisioning. Therefore, in this case study, the hourly household energy demand (Table 5) must be provided as input to the SHyFTA model. Comparing the two design configurations, with and without storage system, the results of the SHyFTA model allowed to estimate the (1) expected energy production that is injected to the grid ( Figure 15) and the (2) energy required (from the grid) to satisfy the energy consumptions of the household (Figure 16). Comparing the two design configurations, with and without storage system, the results of the SHyFTA model allowed to estimate the (1) expected energy production that is injected to the grid ( Figure 15) and the (2) energy required (from the grid) to satisfy the energy consumptions of the household (Figure 16).    Comparing the two design configurations, with and without storage system, the results of the SHyFTA model allowed to estimate the (1) expected energy production that is injected to the grid ( Figure 15) and the (2) energy required (from the grid) to satisfy the energy consumptions of the household (Figure 16).   These two variables depend on the profile of the hourly energy demand of Table 5. Clearly, the investment of the PV system (in both the configurations) reduces the energy demand required from the grid of the household and has an impact into the electrical bill. This latter information, together with the cost of the investments, is used to perform a cash-flow analysis that, in this example, has been carried out considering the Italian market scenario. Figures 17 and 18 show the results of the economic analysis. In particular, Figure 17 presents the cumulated discounted cash-flow along the investment horizon time and Figure 18 the payback time (in years) as a function of the battery cost (expressed in €/kWh). These two variables depend on the profile of the hourly energy demand of Table 5. Clearly, the investment of the PV system (in both the configurations) reduces the energy demand required from the grid of the household and has an impact into the electrical bill. This latter information, together with the cost of the investments, is used to perform a cash-flow analysis that, in this example, has been carried out considering the Italian market scenario. Figures 17 to 18 show the results of the economic analysis. In particular, Figure 17 presents the cumulated discounted cash-flow along the investment horizon time and Figure 18 the payback time (in years) as a function of the battery cost (expressed in €/kWh).   These two variables depend on the profile of the hourly energy demand of Table 5. Clearly, the investment of the PV system (in both the configurations) reduces the energy demand required from the grid of the household and has an impact into the electrical bill. This latter information, together with the cost of the investments, is used to perform a cash-flow analysis that, in this example, has been carried out considering the Italian market scenario. Figures 17 to 18 show the results of the economic analysis. In particular, Figure 17 presents the cumulated discounted cash-flow along the investment horizon time and Figure 18 the payback time (in years) as a function of the battery cost (expressed in €/kWh).   In Figure 17, it is clearly demonstrated that the investment of a PV system without battery is more convenient than an investment of a PV system with battery. In fact, the final NPV (20y, 3%) of the investment without battery amounts to 4496 €, versus 4410 € of the investment with battery. On the other hand, in Figure 18 it is shown that the payback time of the investment without battery is nine years, whereas with the battery amounts to 11 years. Only a drastic cost decrease of the battery system can modify the results in favor of the battery installation. In fact, as shown in Figure 18, the payback time of the PV system with battery becomes equal only if the cost of the battery is lowered to 150 €/kWh.

Distillation Column
The results presented in this section refer to the case study discussed in Section 3. Whereas in the previous case studies the SHyFTA was used to evaluate the system reliability, the service availability and the NPV, in this final example the SHyFTA model is coded to dimension the capacity of the decanter that maximizes the yearly profit (Incomes-Costs) along the useful life of the system (T UL = 10 years). First of all, knowing that the cost of installation for a solid tank with a capacity of 0.1 m 3 amounts to $15,000 USD, it is possible to compute the cost of installation for different sizes of the solid tank. In fact, the well-known "0.6 rule" states that the relationship between the increase in equipment cost (C) and the increase in capacity (V) is given by the following scaling law [61]: where m denotes the scale coefficient (= 0.6). Therefore, Table 6 resumes the installation costs and the corresponding yearly expense (amortization) spread during the useful life, assuming a discount rate Therefore, it is possible to extract from the SHyFTA simulation the yearly volume of the mixture processed in a good way (OK), dumped (NOK), or not processed at all (Lost). The simulation results are shown in Figure 19 (the left axis is assigned to the OK, whereas the right to the NOK and Lost volumes).
Knowing that the purified liquid is sold at the price of $ 0.01 /m 3 , it is possible to compute the yearly profit and the yearly missed income. Moreover, among the yearly costs, it is necessary to consider a penalty (of $ 50,000 /year) that must be sustained if the safety environment requirements are not fulfilled. In this case, the safety is the probability of the system to avoid, during a year of activity, that the yearly loss of mixture in the environment exceeds a certain threshold, fixed by the authority. Therefore, the SHyFTA model can be coded so as to compute this probability on the basis of the volume of mixture not processed by the distillation column (Volume Lost) during a year of operations and of a maximum threshold (V THmax = 4.85 × 10 5 m 3 ). Figure 20 shows the unsafety simulated for a year of activity. . Figure 19. Volume of the mixture processed in a good way (OK), dumped (NOK), or not processed at all (Lost).
Knowing that the purified liquid is sold at the price of $ 0.01 /m 3 , it is possible to compute the yearly profit and the yearly missed income. Moreover, among the yearly costs, it is necessary to consider a penalty (of $ 50,000 /year) that must be sustained if the safety environment requirements are not fulfilled. In this case, the safety is the probability of the system to avoid, during a year of activity, that the yearly loss of mixture in the environment exceeds a certain threshold, fixed by the authority. Therefore, the SHyFTA model can be coded so as to compute this probability on the basis of the volume of mixture not processed by the distillation column (Volume Lost) during a year of operations and of a maximum threshold (VTHmax = 4.85 × 10 5 m 3 ). Figure 20 shows the unsafety simulated for a year of activity.
. Figure 20. Unsafety of the process of distillation. . Figure 19. Volume of the mixture processed in a good way (OK), dumped (NOK), or not processed at all (Lost).
Knowing that the purified liquid is sold at the price of $ 0.01 /m 3 , it is possible to compute the yearly profit and the yearly missed income. Moreover, among the yearly costs, it is necessary to consider a penalty (of $ 50,000 /year) that must be sustained if the safety environment requirements are not fulfilled. In this case, the safety is the probability of the system to avoid, during a year of activity, that the yearly loss of mixture in the environment exceeds a certain threshold, fixed by the authority. Therefore, the SHyFTA model can be coded so as to compute this probability on the basis of the volume of mixture not processed by the distillation column (Volume Lost) during a year of operations and of a maximum threshold (VTHmax = 4.85 × 10 5 m 3 ). Figure 20 shows the unsafety  It is possible to see that the system with a solid tank of 0.1 m 3 has a probability equals to 1 to exceed the threshold V THmax . Therefore, this solution must be certainly discarded. In order to evaluate the other solutions and account also for the expected penalty, the corresponding unsafety of each solution must be multiplied with the penalty. Table 7 resumes the values of the yearly profits, expenses and expected penalty for the possible solutions of the solid tank. Therefore, on the basis of the dependability analysis performed, it is possible to compare the yearly net profit of the various solutions.
As shown in Table 8 the most convenient technical solution, with an annual net profit of $ 19,302/year, is the one with a solid tank of 0.9 m 3 of capacity. According to the results the installation of a larger solid tanks can decrease the unsafety of the system and increase the mean availability but, from the economic viewpoint, is not justified.

Conclusions
With the advent of hybrid methodologies, the capability and expressiveness of the stochastic models have increased to such an extent that they require new tools able to assist engineers and risk practitioners in the practice of the dependability assessment. Therefore, the need for a ready-to-use computer-aided tool for the dependability assessment of industrial systems has motivated the research proposed in this paper.
To demonstrate the importance of the topic, a survey of the state-of-the-art has been presented showing the evolution of the Fault Tree Analysis (FTA) related methods. Nowadays, the conception of the Stochastic Hybrid Fault Tree Automaton (SHyFTA) represents a mature alternative to FTA to analyse complex dependable systems under the viewpoint of the Dynamic Reliability Probabilistic Assessment.
The SHyFTOO Matlab ® library is a software toolbox that allows for the modelling and the simulation of such hybrid models. The main novelty presented in this paper, over the previous implementation of the SHyFTOO [22], is the capability to couple SHyFTOO with Simulink. Simulink is a toolbox integrated in Matlab ® , specifically conceived for the modelling of dynamic systems. This simplifies tremendously the coding of the physical process of a SHyFTA model and provides the user with a toolbox able to compute not only the typical means of the dependability assessment but also other custom performance indexes of the process under analysis.
In this paper, the core concepts of the SHyFTA methodology have been summarized, and the main features and components of the SHyFTOO Simulink library have been presented, in order to offer to the reader a practical step-by-step guide to implement and solve, in an easy way, a SHyFTA model. Afterward, several case studies have been discussed and solved to demonstrate the flexibility of the tool in modelling and solving non-trivial industrial applications. Moreover, the tool allowed to evaluate many properties of an industrial process that cannot be described with non-hybrid methodologies, thus obtaining results and process information that would not be achieved with a different approach.
In future studies, the authors aim to improve the performance of the tool investigating the possibility to integrate it with an algorithm for the variance reduction and for the implementation of rare-event approximations.
The software library developed in Matlab, including some SHyFTA case studies can be freely downloaded from the github webpage of the corresponding author at the following link: https: //github.com/chiacchiof/SHyFTOO-Matlab.  -SEQ ('SEQ'): the SEQ gate is a dynamic gate that forces the input at failing from the right to the left. It triggers as soon the last input has failed; -FDEP ('FDEP'): the FDEP gate is a dynamic gate that imposes the failure of the inputs connected to the gate on the basis of the state of its primary input.
The parameter IsFailureGate is a boolean value that allows to specify if the gate can be restored to the good status if the inputs are repairable [21]. As for the inputs of the gate, the SHyFTOO allows the codification of extended fault tree. The parameter IsFailureGate is a boolean value that allows to specify if the gate can be restored to the good status if the inputs are repairable [21]. As for the inputs of the gate, the SHyFTOO allows the codification of extended fault tree. These fault trees support the cascade of dynamic gates, (i.e., the inputs of the dynamic gates can be of basic events or gates) [36]. These relevant blocks allow the coupling with the stochastic model (the fault tree components defined in the script initFaultTree) and they have to be used as explained in the following: - The "ITER EVOLUTION" must copied as it is and must not be modified. It is an ensemble of blocks that control, for each iteration, the simulation time. When the mission time Tm has reached, a new iteration is automatically restarted. The marking of the "ITER EVOLUTION"components are summarized in Table A1.

Mission Time Tm
It represents the mission time of the system. This parameter is set in the script SHyFTAmain.m The Simulink block must not be deleted or modified.

Clock
It is the Simulink block representing the simulation clock. The Simulink block must not be deleted.

TimeSimulationClock
This block evaluates the time evolution until the end of the mission time. The Simulink block must not be deleted or modified. - The "RACE CONDITION" must be copied as it is and must not be modified. The marking of the "RACE CONDITION" components are summarized in Table A2. This ensemble controls the occurrence of the next discrete event time-point of status change for the basic events of the fault tree. If the nextEvent time-point is higher than the current clock, the iteration is paused, and the fault tree status is evaluated. This block is in race condition against all the "GENERIC HYBRID BASIC EVENT" ensembles, as shown in Figure A5.  These relevant blocks allow the coupling with the stochastic model (the fault tree components defined in the script initFaultTree) and they have to be used as explained in the following: - The "ITER EVOLUTION" must copied as it is and must not be modified. It is an ensemble of blocks that control, for each iteration, the simulation time. When the mission time Tm has reached, a new iteration is automatically restarted. The marking of the "ITER EVOLUTION" components are summarized in Table A1.

Element Description
Mission Time T m It represents the mission time of the system. This parameter is set in the script SHyFTAmain.m The Simulink block must not be deleted or modified. These relevant blocks allow the coupling with the stochastic model (the fault tree components defined in the script initFaultTree) and they have to be used as explained in the following: - The "ITER EVOLUTION" must copied as it is and must not be modified. It is an ensemble of blocks that control, for each iteration, the simulation time. When the mission time Tm has reached, a new iteration is automatically restarted. The marking of the "ITER EVOLUTION"components are summarized in Table A1.

Mission Time Tm
It represents the mission time of the system. This parameter is set in the script SHyFTAmain.m The Simulink block must not be deleted or modified.

Clock
It is the Simulink block representing the simulation clock. The Simulink block must not be deleted.

TimeSimulationClock
This block evaluates the time evolution until the end of the mission time. The Simulink block must not be deleted or modified. - The "RACE CONDITION" must be copied as it is and must not be modified. The marking of the "RACE CONDITION" components are summarized in Table A2. This ensemble controls the occurrence of the next discrete event time-point of status change for the basic events of the fault tree. If the nextEvent time-point is higher than the current clock, the iteration is paused, and the fault tree status is evaluated. This block is in race condition against all the "GENERIC HYBRID BASIC EVENT" ensembles, as shown in Figure A5.

Element Description
nextEvent It contains a global variable array that keeps track of the events (of the basic events) that must occur. The Simulink block must not be deleted or modified.
It is the Simulink block representing the simulation clock.
The Simulink block must not be deleted.

TimeSimulationClock
This block evaluates the time evolution until the end of the mission time.
The Simulink block must not be deleted or modified.
-The "RACE CONDITION" must be copied as it is and must not be modified. The marking of the "RACE CONDITION" components are summarized in Table A2. This ensemble controls the occurrence of the next discrete event time-point of status change for the basic events of the fault tree. If the nextEvent time-point is higher than the current clock, the iteration is paused, and the fault tree status is evaluated. This block is in race condition against all the "GENERIC HYBRID BASIC EVENT" ensembles, as shown in Figure A5.

Relation Operator
It is the Simulink block that compares two inputs. The output is a Boolean value (true if the input 1 is greater than input 2, false vice versa). The Simulink block must not be deleted.

BEevaluteFT
This block is a Simulink Assertion that evaluates the occurrence of the next basic event in race condition with all the existing hybrid basic events handled in the GENERIC HYBRID BASIC EVENT blocks. The Simulink block must not be deleted or modified.
-For each hybrid basic event of the fault tree, the Simulink model SHyFTA(.slx) has to contain a block ensemble of type "GENERIC HYBRID BASIC EVENT". The task of this block is to verify the status of the generic hybrid basic event that is in race condition against all the other hybrid basic events and the regular basic events. The race condition is handled in the "RACE CONDITION" ensemble.
Conversely from the previous blocks, for each "GENERIC HYBRID BASIC EVENT", the setting of the HBE_Status and HBE blocks must be modified. For instance, as shown in Figure A4, these settings have to specify the name of the hybrid basic event, as defined in the initFaultTree script. In Figure A4 it is possible to see that the Constant value takes the value BasicEventName.Status (e.g., in that case the name of the hybrid basic event specified in the initFaultTree is HBE3). The same reasoning applies for the block HBE.
The marking of the "GENERIC HYBRID BASIC EVENT" components are summarized in Table  A3.

HBE_Status
It contains the SHyFTOO variable of the status of the hybrid basic event of the corresponding "GENERIC HYBRID BASIC EVENT" block. The value of the HBE_Status of each GENERIC HYBRID BASIC EVENT" block must be modified.

Logical Port
It is the Simulink block that invert the logical input. The output is a Boolean value (true if the input is false and false if the input is true). The Simulink block must not be deleted.

HBE_Name
It contains the SHyFTOO variable of the name of the hybrid basic event of the corresponding "GENERIC HYBRID BASIC EVENT" block.
It is the Simulink block that compares two inputs. The output is a Boolean value (true if the input 1 is greater than input 2, false vice versa). The Simulink block must not be deleted.

BEevaluteFT
This block is a Simulink Assertion that evaluates the occurrence of the next basic event in race condition with all the existing hybrid basic events handled in the GENERIC HYBRID BASIC EVENT blocks. The Simulink block must not be deleted or modified.
-For each hybrid basic event of the fault tree, the Simulink model SHyFTA(.slx) has to contain a block ensemble of type "GENERIC HYBRID BASIC EVENT". The task of this block is to verify the status of the generic hybrid basic event that is in race condition against all the other hybrid basic events and the regular basic events. The race condition is handled in the "RACE CONDITION" ensemble.
Conversely from the previous blocks, for each "GENERIC HYBRID BASIC EVENT", the setting of the HBE_Status and HBE blocks must be modified. For instance, as shown in Figure A4, these settings have to specify the name of the hybrid basic event, as defined in the initFaultTree script. In Figure A4 it is possible to see that the Constant value takes the value BasicEventName.Status (e.g., in that case the name of the hybrid basic event specified in the initFaultTree is HBE3). The same reasoning applies for the block HBE.
The marking of the "GENERIC HYBRID BASIC EVENT" components are summarized in Table A3.
Information 2019, 10, x FOR PEER REVIEW 30 of 39

Relation Operator
It is the Simulink block that compares two inputs. The output is a Boolean value (true if the input 1 is greater than input 2, false vice versa). The Simulink block must not be deleted.

BEevaluteFT
This block is a Simulink Assertion that evaluates the occurrence of the next basic event in race condition with all the existing hybrid basic events handled in the GENERIC HYBRID BASIC EVENT blocks. The Simulink block must not be deleted or modified.
-For each hybrid basic event of the fault tree, the Simulink model SHyFTA(.slx) has to contain a block ensemble of type "GENERIC HYBRID BASIC EVENT". The task of this block is to verify the status of the generic hybrid basic event that is in race condition against all the other hybrid basic events and the regular basic events. The race condition is handled in the "RACE CONDITION" ensemble.
Conversely from the previous blocks, for each "GENERIC HYBRID BASIC EVENT", the setting of the HBE_Status and HBE blocks must be modified. For instance, as shown in Figure A4, these settings have to specify the name of the hybrid basic event, as defined in the initFaultTree script. In Figure A4 it is possible to see that the Constant value takes the value BasicEventName.Status (e.g., in that case the name of the hybrid basic event specified in the initFaultTree is HBE3). The same reasoning applies for the block HBE.
The marking of the "GENERIC HYBRID BASIC EVENT" components are summarized in Table  A3.

Element Description
HBE_Status It contains the SHyFTOO variable of the status of the hybrid basic event of the corresponding "GENERIC HYBRID BASIC EVENT" block. The value of the HBE_Status of each GENERIC HYBRID BASIC EVENT" block must be modified. It is the Simulink block that invert the logical input. The output is a

HBE_Status
It contains the SHyFTOO variable of the status of the hybrid basic event of the corresponding "GENERIC HYBRID BASIC EVENT" block. The value of the HBE_Status of each GENERIC HYBRID BASIC EVENT" block must be modified.
Logical Port Figure A4. Customization of a Simulink block of type Constant.

HBE_Status
It contains the SHyFTOO variable of the status of the hybrid basic event of the corresponding "GENERIC HYBRID BASIC EVENT" block. The value of the HBE_Status of each GENERIC HYBRID BASIC EVENT" block must be modified.

Logical Port
It is the Simulink block that invert the logical input. The output is a Boolean value (true if the input is false and false if the input is true). The Simulink block must not be deleted.

HBE_Name
It contains the SHyFTOO variable of the name of the hybrid basic event of the corresponding "GENERIC HYBRID BASIC EVENT" block.
It is the Simulink block that invert the logical input. The output is a Boolean value (true if the input is false and false if the input is true). The Simulink block must not be deleted.

HBE_Name
It contains the SHyFTOO variable of the name of the hybrid basic event of the corresponding "GENERIC HYBRID BASIC EVENT" block. The value of the HBE_Name of each GENERIC HYBRID BASIC EVENT" block must be modified.
The HBE_B block of the "GENERIC HYBRID BASIC EVENT" shown in Figure A5 controls the next time-event of the corresponding hybrid basic event. The marking of the components of this block are described in Table A4. The value of the HBE_Name of each GENERIC HYBRID BASIC EVENT" block must be modified.
The HBE_B block of the "GENERIC HYBRID BASIC EVENT" shown in Figure A5 controls the next time-event of the corresponding hybrid basic event. The marking of the components of this block are described in Table A4. The physical process is the most suitable simulation scope for the hybrid events that are characterized by dynamic parameters that change throughout the evolution of the physical process. This block contains a Matlab ® function (Dynamic parameters) that updates the failure parameters of the component at each time-step of the iteration and pass it, together with other inputs (e.g., clock, mission time Tm, and id of the HBE) to an "Interpreted Matlab Function", named "Sample HBE", that samples whether the hybrid basic event has occurred or not. In fact, at any timestep, the "Interpreted Matlab Function" inverts the probability density function characterizing the hybrid basic event and compares it with a real random value in (0, 1]. If the comparison is positive, an exception is raised and handled in the HBEevaluateFT assertion block. In this case, the Simulink iteration is paused, and the control is passed to the SHyFTOO library to verify the fault tree status. Otherwise, the Simulink iteration continues normally. The mechanism hereby described is the one that allows the race condition together with the block "RACE CONDITION".

Element Description
Clock It takes the input of the Simulink clock (refer to Table A1). This input must not be deleted or modified.

Tm
It takes the input of the mission time (refer to Table A1). This input must not be deleted or modified. The physical process is the most suitable simulation scope for the hybrid events that are characterized by dynamic parameters that change throughout the evolution of the physical process. This block contains a Matlab ® function (Dynamic parameters) that updates the failure parameters of the component at each time-step of the iteration and pass it, together with other inputs (e.g., clock, mission time Tm, and id of the HBE) to an "Interpreted Matlab Function", named "Sample HBE", that samples whether the hybrid basic event has occurred or not. In fact, at any timestep, the "Interpreted Matlab Function" inverts the probability density function characterizing the hybrid basic event and compares it with a real random value in (0, 1]. If the comparison is positive, an exception is raised and handled in the HBEevaluateFT assertion block. In this case, the Simulink iteration is paused, and the control is passed to the SHyFTOO library to verify the fault tree status. Otherwise, the Simulink iteration continues normally. The mechanism hereby described is the one that allows the race condition together with the block "RACE CONDITION".

Element Description
Clock It takes the input of the Simulink clock (refer to Table A1). This input must not be deleted or modified.
T m It takes the input of the mission time (refer to Table A1). This input must not be deleted or modified.

ID_HBE
It takes the input of the HBE_Name (refer to Table A3). This input must not be deleted or modified.

Dynamic parameters
Additional input block that can vary (depends on the physical process). The input must be modified according to the dependencies between dynamic failure parameters and the physical process.

SAMPLE_HBE
It is a function that samples the hybrid basic event occurrence. This input must not be deleted or modified.

Appendix A.4. Script shyftaMetrics
The script shyftaMetrics must be modified on the basis of the custom user variables that have been defined in the main script SHyFTAmain and that need to be updated at the end of each iteration. This can be, for example, the case of the continuous variables of the physical process that are averaged at the end of the simulation. The script in the Algorithm 4 applies for the generic variable name variableName defined in the SHyFTAmain script. For each custom variable, the lines from 2 to 8 must be copied and pasted, modifying the variable of interest. As already said, in order to keep track of the evolution of these variables, a dedicated Simulink block ToWorkspace has to be used in the Simulink model SHyFTA.slx, setting, as shown in Figure A5, the parameter VariableName = variableNameCum and the Save format = Timeseries.

Appendix A.5. Properties of the Basic Events and Gates
Tables A5 and A6 contain the list of attributes (or properties) of the SHyFTOO library defined for the basic events and gates. They can be easily accessed in the Matlab ® workspace using the syntax Component.Property (e.g., BE3.Name).

Appendix A.4. Script shyftaMetrics
The script shyftaMetrics must be modified on the basis of the custom user variables that have been defined in the main script SHyFTAmain and that need to be updated at the end of each iteration. This can be, for example, the case of the continuous variables of the physical process that are averaged at the end of the simulation. The script in the Algorithm 4 applies for the generic variable name variableName defined in the SHyFTAmain script. For each custom variable, the lines from 2 to 8 must be copied and pasted, modifying the variable of interest. As already said, in order to keep track of the evolution of these variables, a dedicated Simulink block ToWorkspace has to be used in the Simulink model SHyFTA.slx, setting, as shown in Figure A5, the parameter VariableName = variableNameCum and the Save format = Timeseries.

Appendix A.5. Properties of the Basic Events and Gates
Tables A5 and A6 contain the list of attributes (or properties) of the SHyFTOO library defined for the basic events and gates. They can be easily accessed in the Matlab ® workspace using the syntax Component.Property (e.g., BE3.Name).

Appendix A.4. Script shyftaMetrics
The script shyftaMetrics must be modified on the basis of the custom user variables that have been defined in the main script SHyFTAmain and that need to be updated at the end of each iteration. This can be, for example, the case of the continuous variables of the physical process that are averaged at the end of the simulation. The script in the Algorithm A1 applies for the generic variable name variableName defined in the SHyFTAmain script. For each custom variable, the lines from 2 to 8 must be copied and pasted, modifying the variable of interest. As already said, in order to keep track of the evolution of these variables, a dedicated Simulink block ToWorkspace has to be used in the Simulink model SHyFTA.slx, setting, as shown in Figure A5, the parameter VariableName = variableNameCum and the Save format = Timeseries.

Appendix A.5. Properties of the Basic Events and Gates
Tables A5 and A6 contain the list of attributes (or properties) of the SHyFTOO library defined for the basic events and gates. They can be easily accessed in the Matlab ® workspace using the syntax Component.Property (e.g., BE3.Name). This section contains information about the Simulink models of the case study discussed in Section 3. They have been organized in different Tables A7-A10 that refer respectively to the manuscript Figures 4-6 and Figure 8. Table A7. Marking of the components of the hybrid basic events HBE1 and HBE5 (refer to Figure 4).

Clock
It takes the input of the Simulink clock of Table A1.

Mission Time
It takes the input of the mission time of Table A1. TOP_EVENT_STATUS It is the signal of Top Event status (Top.Status).

HBE1_Status
Refer to the HBE_Status (refer to Table A3). The block must be set = HBE1.Status (or HBE5.Status)

HBE1
Refer to the HBE_Name (refer to Table A3). The block must be set = 'HBE1 (or 'HBE5 ) NOT Refer to the Logical Port "NOT" (refer to Table A3).

Logical Port
Information 2019, 10, x FOR PEER REVIEW 33 of 39

Appendix B
This section contains information about the Simulink models of the case study discussed in Section 3. They have been organized in different Tables A7-A10 that refer respectively to the manuscript Figures 4-6 and 8. Table A7. Marking of the components of the hybrid basic events HBE1 and HBE5 (refer to Figure 4).

Element Description
Clock It takes the input of the Simulink clock of Table  A1.

Mission Time
It takes the input of the mission time of Table A1.

TOP_EVENT_STATUS
It is the signal of Top Event status (Top.Status).

HBE1_Status
Refer to the HBE_Status (refer to Table A3). The block must be set = HBE1.Status (or HBE5.Status)

NOT
Refer to the Logical Port "NOT" (refer to

Appendix B
This section contains information about the Simulink models of the case study discussed in Section 3. They have been organized in different Tables A7-A10 that refer respectively to the manuscript Figures 4-6 and 8. Table A7. Marking of the components of the hybrid basic events HBE1 and HBE5 (refer to Figure 4).

Element Description
Clock It takes the input of the Simulink clock of Table  A1.

Mission Time
It takes the input of the mission time of Table A1. TOP_EVENT_STATUS It is the signal of Top Event status (Top.Status).

HBE1_Status
Refer to the HBE_Status (refer to Table A3). The block must be set = HBE1.Status (or HBE5.Status)

NOT
Refer to the Logical Port "NOT" (refer to Table  A3).

Logical Port
It is the Simulink block "AND" gate.   Figure 5).

Element Description
It is the Simulink block that converts its input in a double type variable.

Integrator
Information 2019, 10, x FOR PEER REVIEW 33 of 39

Appendix B
This section contains information about the Simulink models of the case study discussed in Section 3. They have been organized in different Tables A7-A10 that refer respectively to the manuscript Figures 4-6 and 8. Table A7. Marking of the components of the hybrid basic events HBE1 and HBE5 (refer to Figure 4).

Clock
It takes the input of the Simulink clock of Table  A1.

Mission Time
It takes the input of the mission time of Table A1.

TOP_EVENT_STATUS
It is the signal of Top Event status (Top.Status).

HBE1_Status
Refer to the HBE_Status (refer to Table A3). The block must be set = HBE1.Status (or HBE5.Status)

NOT
Refer to the Logical Port "NOT" (refer to Table  A3).

Logical Port
It is the Simulink block "AND" gate.   Figure 5).

Element Description Clock
It takes the input of the Simulink clock (refer to Table A1).
It is the Simulink block integrator. It returns the aging of the component. Input 1 is the on/off of the component Input 2 is the reset condition corresponding with the status of the HBE1_Status of the component. In fact, it is assumed that the components aging is reset to zero when the component is restored after a fault. Table A8. Marking of the components of the hybrid basic events HBE2 (referring to Figure 5).

Element Description
Clock It takes the input of the Simulink clock (refer to Table A1).

Mission Time
It takes the input of the mission time (refer to Table A1).

TOP_EVENT_STATUS
It is the signal of Top Event status (Top.Status).

HBE2_Status
Refer to the HBE_Status (refer to Table A3). The block must be set = HBE2.Status

HBE2_InUseBy
It contains the SHyFTOO variable of the property "InUseBy" of the hybrid basic event HBE2. The block must be set = HBE2.InUseBy

HBE2
Refer to the HBE_Name (refer to Table A3). The block must be set = 'HBE2 NOT Refer to the Logical Port "NOT" (refer to Table A3).

Logical Port (left AND)
Section 3. They have been organized in different Tables A7-A10 that refer respectively to the manuscript Figures 4-6 and 8. Table A7. Marking of the components of the hybrid basic events HBE1 and HBE5 (refer to Figure 4).

Element Description
Clock It takes the input of the Simulink clock of Table  A1.  Mission Time  It takes the input of the mission time of Table A1. TOP_EVENT_STATUS It is the signal of Top Event status (Top.Status).

HBE1_Status
Refer to the HBE_Status (refer to Table A3). The block must be set = HBE1.Status (or HBE5.Status)

NOT
Refer to the Logical Port "NOT" (refer to Table  A3).

Logical Port
It is the Simulink block "AND" gate.  .   Table A8. Marking of the components of the hybrid basic events HBE2 (referring to Figure 5).

Element Description Clock
It takes the input of the Simulink clock (refer to

Appendix B
This section contains information about the Simulink models of the case study discussed in Section 3. They have been organized in different Tables A7-A10 that refer respectively to the manuscript Figures 4-6 and 8. Table A7. Marking of the components of the hybrid basic events HBE1 and HBE5 (refer to Figure 4).

Element Description
Clock It takes the input of the Simulink clock of Table  A1.

Mission Time
It takes the input of the mission time of Table A1.

TOP_EVENT_STATUS
It is the signal of Top Event status (Top.Status).

HBE1_Status
Refer to the HBE_Status (refer to Table A3). The block must be set = HBE1.Status (or HBE5.Status)

HBE1
Refer to the HBE_Name (refer to Table A3).
The block must be set = 'HBE1′ (or 'HBE5′) NOT Refer to the Logical Port "NOT" (refer to Table  A3).

Logical Port
It is the Simulink block "AND" gate.  .   Table A8. Marking of the components of the hybrid basic events HBE2 (referring to Figure 5).

Element Description Clock
It takes the input of the Simulink clock (refer to Table A1).
It is the Simulink block "AND" gate. It returns a Boolean value (true if both the inputs are true, false in all the other cases). The output of this block allows to integrate the aging of the component represented by the hybrid basic event. The aging increases only if the component is good (and not in standby) and the system is available (if the TOP EVENT STATUS is false the component is supposed not to work).

Appendix B
This section contains information about the Simulink models of the case study discussed in Section 3. They have been organized in different Tables A7-A10 that refer respectively to the manuscript Figures 4-6 and 8. Table A7. Marking of the components of the hybrid basic events HBE1 and HBE5 (refer to Figure 4).

Element Description
Clock It takes the input of the Simulink clock of Table  A1.  Mission Time  It takes the input of the mission time of Table A1. TOP_EVENT_STATUS It is the signal of Top Event status (Top.Status).

HBE1_Status
Refer to the HBE_Status (refer to Table A3). The block must be set = HBE1.Status (or HBE5.Status)

NOT
Refer to the Logical Port "NOT" (refer to Table  A3).

Logical Port
It is the Simulink block "AND" gate.  .   Table A8. Marking of the components of the hybrid basic events HBE2 (referring to Figure 5).

Element Description Clock
It takes the input of the Simulink clock (refer to Table A1).
It is the Simulink block that converts its input in a double type variable.

Integrator
This section contains information about the Simulink models of the case study discussed in Section 3. They have been organized in different Tables A7-A10 that refer respectively to the manuscript Figures 4-6 and 8. Table A7. Marking of the components of the hybrid basic events HBE1 and HBE5 (refer to Figure 4).

Element Description
Clock It takes the input of the Simulink clock of Table  A1.  Mission Time  It takes the input of the mission time of Table A1. TOP_EVENT_STATUS It is the signal of Top Event status (Top.Status).

HBE1_Status
Refer to the HBE_Status (refer to Table A3). The block must be set = HBE1.Status (or HBE5.Status)

NOT
Refer to the Logical Port "NOT" (refer to Table  A3).

Logical Port
It is the Simulink block "AND" gate.  .   Table A8. Marking of the components of the hybrid basic events HBE2 (referring to Figure 5).

Element Description Clock
It takes the input of the Simulink clock (refer to Table A1).
It is the Simulink block integrator. It returns the aging of the component. Input 1 is the on/off of the component. Input 2 is the reset condition corresponding with the status of the HBE2_Status of the component. It is assumed that the components aging is reset to zero when the component is restored after a fault. Table A9. Marking of the components of the hybrid basic events HBE6 (referring to Figure 6).

% Particle
It is a Simulink "uniform random number" generator that models the percentage of solid particles of the mixture. The block is set as shown in Figure A6b QI It is a Simulink "uniform random number" generator that models the instantaneous volume (flow rate) of mixture. The block is set as shown in Figure A6a  Integrator It is the Simulink block integrator. It returns the volume of the solid particle deposited in the solid tank. This output is the input of the HBE6_SAMPLING block that verifies when the solid tank threshold is reached (event occurrence of the hybrid basic event HBE6). Input 1 is the output of the product block above described. Input 2 is the reset condition, represented by the status of the top event. When the top event occurs, the process is stopped (and the solid tank is emptied). Figure A6 shows the setting of the block "QI" and "% Particle" according to the characteristics of the physical process described in the case study. In order to avoid the repetition of iterations with the same evolution (for the correct implementation of a Monte Carlo simulation), the Seed parameter must be set with a random number. This feature can be obtained by using a global array variable randomSeed that is initialized at the beginning of each new iteration (as previously described in the script <<SHyFTAmain.m>>). Table A10. Marking of the components of the physical process (referring to Figure 8).

Element
Description It is the signal of QI (refer to It is the Simulink block of the product operation. It returns the instantaneous flow QI processed by the distillation column when the system is working (if the Top Event status is bad, the process is stopped and the product is null). The output is the input of the "To Workspace" block QI_Mixture.

TOP_EVENT_STATUS
It is the neglection of Top Event status (Top.Status).

HBE2_Status
Refer to the HBE_Status (refer to Table A3). The block must be set = HBE6.Status

HBE6
Refer to the HBE_Name (refer to Table A3). The block must be set = 'HBE6 Product

HBE2_Status
Refer to the HBE_Status (refer to Table A3). The block must be set = HBE2.Status

HBE2_InUseBy
It contains the SHyFTOO variable of the property "InUseBy" of the hybrid basic event HBE2. The block must be set = HBE2.InUseBy

HBE2
Refer to the HBE_Name (refer to Table A3). The block must be set = 'HBE2′ NOT Refer to the Logical Port "NOT" (refer to Table A3).

Logical Port (left AND)
It is the Simulink block "AND" gate. It returns a Boolean value (true if both the inputs are true, false in all the other cases). The output of this block returns the working (on/off) condition of the standby component. It is true if the component status is OK and if the standby condition is not false.

Logical Port (right AND)
It is the Simulink block "AND" gate. It returns a Boolean value (true if both the inputs are true, false in all the other cases). The output of this block allows to integrate the aging of the component represented by the hybrid basic event. The aging increases only if the component is good (and not in standby) and the system is available (if the TOP EVENT STATUS is false the component is supposed not to work).

Data type conversion
It is the Simulink block that converts its input in a double type variable.

Integrator
It is the Simulink block integrator. It returns the aging of the component. Input 1 is the on/off of the component. Input 2 is the reset condition corresponding with the status of the HBE2_Status of the component. It is assumed that the components aging is reset to zero when the component is restored after a fault. Table A9. Marking of the components of the hybrid basic events HBE6 (referring to Figure 6).

% Particle
It is a Simulink "uniform random number" generator that models the percentage of solid particles of the mixture. The block is set as shown in Figure A6b QI It is a Simulink "uniform random number" generator that models the instantaneous volume (flow rate) of mixture. The block is set as shown in Figure

HBE2_Status
Refer to the HBE_Status (refer to Table A3). The block must be set = HBE6.Status

HBE6
Refer to the HBE_Name (refer to Table A3). The block must be set = 'HBE6′

Product
It is the Simulink block of the product operation. It returns the quantity of solid particle in the instantaneous infinitesimal volume of mixture (if the Top Event status is bad, the process is stopped and the product is null). The output is the input of the integrator block.
It is the Simulink block of the product operation. It returns the quantity of solid particle in the instantaneous infinitesimal volume of mixture (if the Top Event status is bad, the process is stopped and the product is null). The output is the input of the integrator block.

Integrator
Information 2019, 10, x FOR PEER REVIEW 35 of 39

Integrator
It is the Simulink block integrator. It returns the volume of the solid particle deposited in the solid tank. This output is the input of the HBE6_SAMPLING block that verifies when the solid tank threshold is reached (event occurrence of the hybrid basic event HBE6). Input 1 is the output of the product block above described. Input 2 is the reset condition, represented by the status of the top event. When the top event occurs, the process is stopped (and the solid tank is emptied). Figure A6 shows the setting of the block "QI" and "% Particle" according to the characteristics of the physical process described in the case study.
It is the Simulink block integrator. It returns the volume of the solid particle deposited in the solid tank. This output is the input of the HBE6_SAMPLING block that verifies when the solid tank threshold is reached (event occurrence of the hybrid basic event HBE6). Input 1 is the output of the product block above described. Input 2 is the reset condition, represented by the status of the top event. When the top event occurs, the process is stopped (and the solid tank is emptied). Figure A6 shows the setting of the block "QI" and "% Particle" according to the characteristics of the physical process described in the case study.
Integrator of the HBE6_SAMPLING block that verifies when the solid tank threshold is reached (event occurrence of the hybrid basic event HBE6). Input 1 is the output of the product block above described. Input 2 is the reset condition, represented by the status of the top event. When the top event occurs, the process is stopped (and the solid tank is emptied). Figure A6 shows the setting of the block "QI" and "% Particle" according to the characteristics of the physical process described in the case study. In order to avoid the repetition of iterations with the same evolution (for the correct implementation of a Monte Carlo simulation), the Seed parameter must be set with a random number. This feature can be obtained by using a global array variable randomSeed that is initialized at the beginning of each new iteration (as previously described in the script <<SHyFTAmain.m>>). Table A10. Marking of the components of the physical process (referring to Figure 8).

Element
Description QI It is the signal of QI (refer to It is the Simulink block of the product operation. It returns the instantaneous flow QI processed by the distillation column when the system is working (if the Top Event status is bad, the process is stopped and the product is null). The output is the input of the "To Workspace" block QI_Mixture. In order to avoid the repetition of iterations with the same evolution (for the correct implementation of a Monte Carlo simulation), the Seed parameter must be set with a random number. This feature can be obtained by using a global array variable randomSeed that is initialized at the beginning of each new iteration (as previously described in the script <<SHyFTAmain.m>>). Table A10. Marking of the components of the physical process (referring to Figure 8).

QI
It is the signal of QI (refer to Table A9) TOP_EVENT_STATUS It is the signal of Top Event status (Top.Status) in the same block Integrator solid particle deposited in the solid tank. This output is the input of the HBE6_SAMPLING block that verifies when the solid tank threshold is reached (event occurrence of the hybrid basic event HBE6). Input 1 is the output of the product block above described. Input 2 is the reset condition, represented by the status of the top event. When the top event occurs, the process is stopped (and the solid tank is emptied). Figure A6 shows the setting of the block "QI" and "% Particle" according to the characteristics of the physical process described in the case study. In order to avoid the repetition of iterations with the same evolution (for the correct implementation of a Monte Carlo simulation), the Seed parameter must be set with a random number. This feature can be obtained by using a global array variable randomSeed that is initialized at the beginning of each new iteration (as previously described in the script <<SHyFTAmain.m>>). Table A10. Marking of the components of the physical process (referring to Figure 8).

Element
Description It is the signal of QI (refer to It is the Simulink block of the product operation. It returns the instantaneous flow QI processed by the distillation column when the system is working (if the Top Event status is bad, the process is stopped and the product is null). The output is the input of the "To Workspace" block QI_Mixture. Integrator of the HBE6_SAMPLING block that verifies when the solid tank threshold is reached (event occurrence of the hybrid basic event HBE6). Input 1 is the output of the product block above described. Input 2 is the reset condition, represented by the status of the top event. When the top event occurs, the process is stopped (and the solid tank is emptied). Figure A6 shows the setting of the block "QI" and "% Particle" according to the characteristics of the physical process described in the case study. In order to avoid the repetition of iterations with the same evolution (for the correct implementation of a Monte Carlo simulation), the Seed parameter must be set with a random number. This feature can be obtained by using a global array variable randomSeed that is initialized at the beginning of each new iteration (as previously described in the script <<SHyFTAmain.m>>). Table A10. Marking of the components of the physical process (referring to Figure 8).

Element
Description QI It is the signal of QI (refer to It is the Simulink block of the product operation. It returns the instantaneous flow QI processed by the distillation column when the system is working (if the Top Event status is bad, the process is stopped and the product is null). The output is the input of the "To Workspace" block QI_Mixture.
It is the Simulink block of the product operation. It returns the instantaneous flow QI processed by the distillation column when the system is working (if the Top Event status is bad, the process is stopped and the product is null). The output is the input of the "To Workspace" block QI_Mixture.

QI_Mixture
It is the Matlab "To Workspace" block to store the simulation variable. It is used to store the instantaneous flow QI processed when the system is working.

Product (at the bottom)
Information 2019, 10, x FOR PEER REVIEW 36 of 39

QI_Mixture
It is the Matlab "To Workspace" block to store the simulation variable. It is used to store the instantaneous flow QI processed when the system is working.

Product (at the bottom)
It is the Simulink block of the product operation. It returns the instantaneous flow QI lost when the system is not working (if the Top Event status is bad, the variable ⌐TOP_EVENT_STATUS is true and the product is non null). The output is the input of the "To Workspace" block QI_Lost.

QI_Lost
It is the Matlab "To Workspace" block to store the simulation variable. It is used to store the instantaneous flow QI lost when the system is stopped (or not working).

Detect Change
It is the Simulink block "Detect Change". It triggers a Boolean true value as soon its input changes. It is used to activate the SHA_TOP assertion The input is the top event status.

Assertion Block
It is the Simulink block of the product operation. It returns the instantaneous flow QI lost when the system is not working (if the Top Event status is bad, the variable Information 2019, 10, x FOR PEER REVIEW Integrator It is the Simulink block integrator. It returns the v solid particle deposited in the solid tank. This outp of the HBE6_SAMPLING block that verifies when threshold is reached (event occurrence of the hybr HBE6). Input 1 is the output of the product block above Input 2 is the reset condition, represented by the st event. When the top event occurs, the process is the solid tank is emptied). Figure A6 shows the setting of the block "QI" and "% Particle" according to the ch of the physical process described in the case study. In order to avoid the repetition of iterations with the same evolution (for implementation of a Monte Carlo simulation), the Seed parameter must be set with a rand This feature can be obtained by using a global array variable randomSeed that is initia beginning of each new iteration (as previously described in the script <<SHyFTAmain.m> Table A10. Marking of the components of the physical process (referring to Figure 8)

Element
Description It is the signal of QI (refer to It is the Simulink block of the product operation. It re instantaneous flow QI processed by the distillation colu the system is working (if the Top Event status is bad, th is stopped and the product is null). The output is the input of the "To Workspace" block Q TOP_EVENT_STATUS is true and the product is non null). The output is the input of the "To Workspace" block QI_Lost.

QI_Lost
It is the Matlab "To Workspace" block to store the simulation variable. It is used to store the instantaneous flow QI lost when the system is stopped (or not working).

Detect Change
Information 2019, 10, x FOR PEER REVIEW 36 of 39

QI_Mixture
It is the Matlab "To Workspace" block to store the simulation variable. It is used to store the instantaneous flow QI processed when the system is working.

Product (at the bottom)
It is the Simulink block of the product operation. It returns the instantaneous flow QI lost when the system is not working (if the Top Event status is bad, the variable ⌐TOP_EVENT_STATUS is true and the product is non null). The output is the input of the "To Workspace" block QI_Lost.

QI_Lost
It is the Matlab "To Workspace" block to store the simulation variable. It is used to store the instantaneous flow QI lost when the system is stopped (or not working).

Detect Change
It is the Simulink block "Detect Change". It triggers a Boolean true value as soon its input changes. It is used to activate the SHA_TOP assertion The input is the top event status.

Assertion
Block