Reduction of Petri net maintenance modeling complexity via Approximate Bayesian Computation

Reliability Engineering and System Safety 222 (2022) 108365 Available online 19 February 2022 0951-8320/© 2022 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Contents lists available at ScienceDirect Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress Reduction of Petri net maintenance modeling complexity via Approximate Bayesian Computation Manuel Chiachío a,b,∗, Ali Saleh a,b, Susannah Naybour c, Juan Chiachío a,b, John Andrews c aAndalusian Research Institute in Data Science and Computational Intelligence, University of Granada, 18071 Granada, Spain bDept. Structural Mechanics & Hydraulics Engineering, University of Granada, 18071 Granada, Spain cResilience Engineering Group, University of Nottingham, University Park, Nottingham, UK ARTICLE INFO Keywords: Petri nets Model similarity Bayesian inference Approximate Bayesian Computation Maintenance models ABSTRACT The accurate modeling of engineering systems and processes using Petri nets often results in complex graph representations that are computationally intensive, limiting the potential of this modeling tool in real life applications. This paper presents a methodology to properly define the optimal structure and properties of a reduced Petri net that mimic the output of a reference Petri net model. The methodology is based on Approximate Bayesian Computation to infer the plausible values of the model parameters of the reduced model in a rigorous probabilistic way. Also, the method provides a numerical measure of the level of approximation of the reduced model structure, thus allowing the selection of the optimal reduced structure among a set of potential candidates. The suitability of the proposed methodology is illustrated using a simple illustrative example and a system reliability engineering case study, showing satisfactory results. The results also show that the method allows flexible reduction of the structure of the complex Petri net model taken as reference, and provides numerical justification for the choice of the reduced model structure. 1. Introduction In most countries of the developed world, the lifespan of a large part of critical infrastructures (such as transportation and energy assets) is approaching their end. Methods and models for optimal asset management are thus a necessity to avoid an uncontrolled increase of system failures and unexpected downtimes, whilst keeping maintenance and inspection costs under reasonably lower levels [1]. In the literature, a number of mathematical and computational techniques for optimal asset management have matured over the last decades aiming at optimizing the maintenance and operational costs [2]. Some of these methods include soft-computing techniques such as genetic algorithms [3], neural networks [4], and support vector machines [5]; along with event-based models like Markov chains [6], Bayesian networks [7], and Petri nets [8–11], to name but the most relevant ones. Among the aforementioned approaches, Petri nets (PNs) [12] are typically regarded as powerful modeling tools for maintenance and inspection modeling due to their ability to account for resource availability, concurrency, and synchronization, which are common aspects that underline the majority of the asset management models [13,14]. The basic concepts relative to the theory of PNs are summarized in [15], whereas practical reliability engineering applications of PNs ∗Corresponding author at: Andalusian Research Institute in Data Science and Computational Intelligence, University of Granada, 18071 Granada, Spain. E-mail address: [email protected] (M. Chiachío). can be found in [16,17]. Moreover, PNs can be combined with other computational techniques such as object-oriented programming [9], fuzzy sets [18], expert system modeling [13,19], etc., which greatly increase their suitability for maintenance and inspection modeling in engineering problems. Irrespectively, one of the main challenges in adopting PNs to model maintenance of engineering systems is the handling of complex processes, including their inter-dependencies and workflows, such as those appearing in real-world reliability engineering and maintenance domain. The latter requires large PN representations which are difficult (if not impossible) to interpret graphically, and also heavy computational resources for their numerical evaluation [31,32]. To efficiently approach this practical challenge, an increasing number of researchers have reported progress towards the improvement of the efficiency of PN simulation of large maintenance models [33]. These methods require specialist computational hardware and still can be expensive in terms of cost and computational time. Other existing alternatives consist in reducing the PN complexity to partially alleviate their computational cost [21]. These techniques share the common principle of using a set of rules to shorten specific graphical sub-structures commonly found in PNs while preserving the overall PN behavior unviolated. Table 1 lists the main recent contributions to the literature on this topic and https://doi.org/10.1016/j.ress.2022.108365 Received 2 September 2021; Received in revised form 27 January 2022; Accepted 30 January 2022 Reliability Engineering and System Safety 222 (2022) 108365 2 M. Chiachío et al. Table 1 Summary of recent relevant works about PN reduction methodologies. Method Description Benefits Limitations Refs. Reduction of PN structure following defined rules. PN reduction in a rule-based way, considering model sub-structures. Reduction rules are well understood preventing illogical reductions by the user. Reductions limited to the pre-defined rules. Reduction error not always quantified and limited transitions types considered. [20–24] Symmetry reduction for time Petri net state classes Explore the symmetries in the system to enumerate equivalence states instead of all the reachable states It can be used with state class abstraction. Can be applied to construct State Class Graph of TPNs Beneficial only for the cases where symmetries exists [25] Counting PN markings from reduction equations. Reduction by removing redundant places and transitions and representing the effects of reductions by linear equations Can be combined with other techniques like symmetry reductions. Keep track between the markings of initial and reduced PN. The method can only reduce PNs with specific properties (redundant transitions, redundant places, or place agglomerations). [26,27] Discovering PNs from event logs The firing sequence of transitions in PN structures is used to recreate event logs. PN structures and firing times can be discovered. Not limited to PNs with exponential probability models. Methods can yield different equivalent PN structures. Limited probability models governing the PN delay times are considered. [28] Identification of PN time parameters through comparison to event sequences. PN time parameters are discovered for two normal transitions, in conflict, to match a pre-known event sequence. The algorithm shows good convergence for simple event sequences. Normal distribution for the transitions of the PN is assumed. Large event sequences lead to a low computational efficiency. [29] Reconstructing PN model parameters through comparison to data. PN model parameters are discovered by matching the model output to an expected output. Allows constant probability values within the PN model to be discovered. The approach is not used to discover probability distributions governing PN transition firing times. [30] briefly indicates their key benefits and limitations. From a generic perspective, the existing methods are capable of effectively reducing the model complexity, however this is at the cost of adopting a number of assumptions about the types of nodes that can be used and their parameters, which limit their applicability to real-world cases, precisely where the benefits of these methods can be fully exploited. Hence, there is a need for more versatile and efficient PN reduction techniques to deal with large and complex system reliability and maintenance models. In this context, this paper proposes a novel methodology for PN complexity reduction based on model similarity with respect to a reference PN model through Bayesian learning. More specifically, the proposed approach finds the most plausible values of a set of model parameters governing the reduced PN by comparing both model outputs, namely the reduced and the reference one, within an Approximate Bayesian Computation framework [34,35]. The resulting inference equations are analytically intractable, hence the ABC-SubSim algorithm [36] is adopted whereby the probability density functions (PDFs) of the model parameters are approximated by samples which represent random realizations of the parameter values from the reduced network. The ABC-SubSim method, which is used in this work due to its ability to perform the parameter inference within multi-dimensional parameter spaces and under very low approximation error [36], produces parameter samples through Markov Chain Monte Carlo (MCMC) sampling over a nested sequence of subsets until a predefined approximation error is reached. To overcome the well-known manual scaling required in MCMC sampling [37], this paper also provides an improved version of the ABC-SubSim algorithm with respect to the original one in [36] by adaptive scaling of the MCMC sampling across the subsets. Finally, and as by product, the proposed inference method also produces a quantitative measure of the reduced model suitability through Bayesian evidence computation [38], which aids the task of model selection in the frequent case of multiple reduced models available. The resulting methodology is versatile enough to be applied to any generic PN model with firing times governed by any probability distribution, or even to high-level PN variants like the referred highlevel Petri Nets (HLPN) [39], which allow consideration of nodes with complex behavior like reset and probability transitions, inhibitor arcs, numerical and predicate places, etc. Here, the proposed approach is tested on a real-world case study about a two-pump lube oil hydraulic system whose inspection, maintenance, and renewal process is modeled through a reference, yet complex, HLPN. Three candidate reduced PN models are defined following a criteria of increasing complexity allowing their suitability to be evaluated. The results show that the proposed reduction methodology is able to infer the parameters governing the inspection, maintenance and renewal processes of the equivalent reduced PN models in liaison with the reference HLPN, demonstrating a high level of accuracy in approximation to the reference model alongside a significant reduction in computational cost. The results also show that the reduced models can provide the reliability curves of the two-pump engineering system with high accuracy with respect to the one obtained using the reference HLPN. The remainder of the paper is organized as follows. Section 2 presents the fundamentals behind the proposed computational method from basic principles of PNs to the ABC-SubSim method. Section 3provides description of the adaptive scaling method for the ABC-SubSim algorithm. The mathematical basis and computational aspects of the PN reduction by ABC-SubSim method is given in Section 4. Section 5 illustrates the approach over a simple PN architecture before presenting the approach in application to an engineering case study in Section 6. Section 7discusses the main findings and provides a comparative analysis, and finally Section 8gives concluding remarks. 2. Fundamentals 2.1. Petri net modeling A PN is a directed bipartite graph consisting of two types of nodes, known as places and transitions. The transitions, which are symbolized using boxes, are responsible of the dynamic behavior of the PN, and enable the system to move from one state to another. The places, which are symbolized using circles, represent particular states of the system (e.g. a component in a failure state or an inspection activity which is currently in progress) [15]. The distribution of tokens over the places at a specific execution time determines the state of the PN, and is expressed as a vector of integer values referred to as marking. Mathematically, a PN can be defined as a tuple N=⟨𝐏,𝐓,𝐅,𝐖,𝐌0⟩ [15], where 𝐏and 𝐓denote the set of 𝑛𝑝places and 𝑛𝑡transitions Reliability Engineering and System Safety 222 (2022) 108365 3 M. Chiachío et al. of the PN, whose connections are expressed through the set of edges 𝐅⊆(𝐏×𝐓) ∪ (𝐓×𝐏).𝐖∶𝐅→N>0is a weight function which assigns a natural number (1by default) to each arc within 𝐅, and finally 𝐌0∈N𝑛𝑝is the initial marking. Additionally, in this work the term ∙𝐏𝑡will be used to denote the pre-set of 𝑡, namely the set of input places of transition 𝑡, and 𝐏∙ 𝑡for the set of output places of transition 𝑡, also referred to as the post-set of 𝑡. The dynamics of a PN can be mathematically described through a state equation defined as follows: 𝐌𝑘+1 =𝐌𝑘+𝐀𝑇𝐮𝑘(1) where 𝐌𝑘is the marking at time step 𝑘, being 𝑘∈Na time index. 𝐀is the incidence matrix, which can be obtained as the result of subtracting the backward incidence matrix 𝐀−=[𝑎− 𝑖𝑗]from the forward incidence matrix 𝐀+=[𝑎+ 𝑖𝑗]. The terms 𝑎− 𝑖𝑗 and 𝑎+ 𝑖𝑗 coincide with the weights 𝑤− 𝑖𝑗 and 𝑤+ 𝑖𝑗 of the arcs of transition 𝑡𝑖from its pre-set and post-set places 𝑝𝑗, respectively, where 𝑖= 1,…, 𝑛𝑡, 𝑗 = 1,…, 𝑛𝑝. Finally, the term 𝐮𝑘= (𝑢1,𝑘, 𝑢2,𝑘,…, 𝑢𝑛𝑡,𝑘)𝑇in Eq. (1) is the firing vector whose elements are binary values which adopt the unity if the transition is fired and 0 otherwise, according to the firing rule [15]. In practical applications dealing with the modeling of dynamical systems, transitions are typically assigned time delays, referred to as 𝜏. The resulting PNs are called Timed Petri Nets (TPN) if the delays are deterministic, and Stochastic Petri Nets (SPN) if the delays are specified by a probability model [15]. In this work, PN variants referred to as high-level Petri nets (HLPN) [39] are used to consider higher complexities in the systems. One of the variants is the probability transitions whose production of tokens to their output arcs is carried out under a specified probability. In mathematical terms, if 𝑡𝑖is a so-called probability transition, then a probability value 𝜋𝑖𝑧 is assigned to every 𝑧th output arc from its post-set 𝐏∙ 𝑡𝑖. Thus, the associated elements to 𝑡𝑖from the forward matrix can be expressed as: 𝑎+ 𝑖𝑧 =𝑤+ 𝑖𝑧I𝑧(2) where 𝑤+ 𝑖𝑧 is the output arc weight from transition 𝑡𝑖to place 𝑝𝑧∈𝐏∙ 𝑡𝑖, and I𝑧is a Boolean indicator given as follows: I𝑧={1,if 𝑧=𝑧′ 0,otherwise (3) In the last equation, 𝑧′is the index of the output arc chosen according to the associated probability of each output arc. Finally note that HLPN models used in this work make also use of inhibitor arcs [15] (represented with a small circle as ending, like those depicted in Fig. 4) which produce the opposite effect of the firing rule, i.e. they prevent a transition from being enabled once its pre-set places are marked. Fig. 1 is provided to illustrate a sample HLPN of five places and four transitions to model the maintenance of an engineering component under abnormal condition. Places 𝑝1,𝑝2,𝑝3and 𝑝4describe the normal, degraded, critical, and failed states of a component, respectively, whereas place 𝑝5accounts for the number of times the component changes to any malfunction state. To represent the probability of moving from the normal to any abnormal state, the probability transition 𝑡1 is used, whose delay time 𝜏1accounts for the time required before the component malfunction. Transitions 𝑡2,𝑡3and 𝑡4describe maintenance actions to return the component to its normal condition, and their delay times 𝜏2,𝜏3and 𝜏4, describe the times of the maintenance operations. 2.2. Approximate Bayesian computation method The approximate Bayesian computation (ABC) method [34] is under the category of Bayesian statistics and related to the Bayesian model updating approach [38] in those cases where the likelihood function is difficult, or analytically intractable. In Bayesian model updating, the interest is about updating the prior knowledge about a set parameters 𝜃∈𝜣⊂R𝑑of a parameterized model based on the information gained from a set of data, 𝑦∈⊂R𝓁. The term is the observation space which contains all possible observational outcomes according to a specific model class . If multiple candidate models classes are available, this methodology also provides a framework to assess the plausibility of each model class. The ABC technique produces posterior samples of the pairs (𝜃, 𝑥) ∈ ⊆𝜣×which makes the model response 𝑥=𝑥(𝜃) ∼ 𝑝(𝑥|𝜃)𝑝(𝜃)lay within a defined region around the data 𝑦. This region can be formally defined as: 𝜖(𝑦) = {𝑥∈∶𝜌(𝜂(𝑥), 𝜂(𝑦))⩽𝜖}(4) where 𝜌(⋅) ∶ R𝑑×𝓁→Ris a metric function to evaluate the closeness of the simulated model output 𝑥to the data 𝑦,𝜖is a tolerance parameter, and 𝜂(⋅)is a summary statistic [40] which, if required, allows the comparison between 𝑥and 𝑦in a weak manner. The referred posterior samples produce an approximation of the posterior PDF 𝑝𝜖(𝜃, 𝑥|𝑦) = 𝑝(𝜃, 𝑥|𝑥∈𝜖(𝑦))which assigns higher probability density values to those pairs (𝜃, 𝑥) ∈ that satisfy the condition 𝜌(𝜂(𝑥), 𝜂(𝑦))⩽𝜖. 2.3. ABC-SubSim method As previously mentioned, the success of an ABC method is dependent on the choice of certain hyper-parameters like the summary statistic 𝜂(⋅), the metric 𝜌, and a tolerance value 𝜖. Among them, the tolerance hyper-parameter is of special importance since it entails a critical trade-off between accuracy of the posterior approximation and computational cost. In the literature, a number of techniques can be found to address this trade-off by combining the ABC principles with efficient sampling algorithms. The reader is referred to [35] for a comprehensive overview of the ABC methods. Irrespective of the efficiency demonstrated by some of the referred methods, using 𝜖→0still translates into heavy computation except for those which sequentially adapt 𝜖until reaching a desirable small value [41]. Among those, the so called approximate Bayesian computation by subset simulation, also referred to as ABC-SubSim algorithm [36], has been proved to be one of the most efficient ABC algorithms in the literature having been included in several well-known ABC user-platforms such as ABCpy [42] and Pi4U [43] (https://github.com/cselab/pi4u), and thus is the algorithm adopted for the PN reduction methodology presented in this paper. The ABC-SubSim depends on the SubSet Simulation, originally presented in [44] as a method to avoid costly or inaccurate rare event simulation that arises due to the existence of one or several small failure regions. By Subset Simulation, the rare event simulation is avoided by the introduction of several intermediate regions 𝑗⊂with larger probabilities by generating conditional samples (𝜃, 𝑥) ∈ 𝑗, 𝑗 = 1,…, 𝑚, where: 𝑗={(𝜃, 𝑥) ∈ ∶𝜌(𝜂(𝑥), 𝜂(𝑦)) ≤𝜖𝑗}(5) In the last equation, the intermediate tolerances 𝜖𝑗follow a decreasing sequence 𝜖1⩾𝜖2…⩾𝜖𝑚≡𝜖. A practical way to implement the ABC-SubSim is by fixing the size of the samples generated in 𝑗−1 that lie in 𝑗as 𝑁𝑃0, where 𝑁 is the amount of samples that populate the 𝑗th simulation level, and 𝑃0is a conditional probability acting as algorithm hyper-parameter defined by the modeler [45]. The outcome of this algorithm is a set of (𝜃, 𝑥) ∈ samples that populate the posterior PDF 𝑝𝜖(𝜃, 𝑥|𝑦). Also, and as by-product, the algorithm straightforwardly produces a measure of the evidence (also referred to as marginal likelihood) of a candidate reduced PN model via the conditional probabilities involved in the simulation [36], as follows: 𝑃𝜖(𝑦|𝑗) = 𝑃(1) 𝑚 ∏ 𝑗=2 𝑃(𝑗|𝑗−1) ≈ 𝑃𝑚 0(6) where 𝑗indicates the 𝑗th model class within a set of available models, which in our case, corresponds to any available reduced PN model. In Section 4, a pseudo-code implementation of the ABC-SubSim algorithm is given within the context of the PN reduction methodology proposed in this work. Reliability Engineering and System Safety 222 (2022) 108365 4 M. Chiachío et al. 3. ABC-SubSim with adaptive sampling As stated in the last section, in ABC-SubSim algorithm 𝑁𝑃0samples are selected at every 𝑗th simulation level acting as seeds of the 𝑗subset, which generate new samples via MCMC sampling. In this work, the Modified Metropolis Hastings algorithm [44] (MMA) is used as MCMC method to generate successive chain values from each seed. Irrespectively, this method requires a tricky selection of a proposal PDF used to evolve values in each Markov chain due its strong impact on the quality of the ABC posterior approximation along with on the speed of convergence of the algorithm [46]. In the literature, the optimal scaling of the standard deviation of the proposal PDF has been considered within the context of Subset Simulation by manual calibration and preliminary trials, thus increasing the computational cost [36,37,45]. To avoid the aforementioned drawback, Vakilzadeh et al. [47] proposed an adaptive method to tune the proposal standard deviation per simulation level for the ABC-SubSim algorithm with soundly results. However the performance of the method in [47] requires the introduction of three new sensible hyper-parameter to the algorithm (see inputs 𝑁𝑎,𝜉1and 𝜆𝑗in [47, Algorithm 2]), which again needs modeler expertise to tune them. More recently, [48] proposed an adaptive method to tune the conditional probability hyper-parameter 𝑃0which also helps to scale the proposal standard deviation, as by product. As in [47], the resulting sampling scheme in [48] is automatic, however it is not efficient in avoiding high sampling rejections in subsets of small size. In this section, a method based on samples’ values distributions is presented to obtain appropriate standard deviation automatically, thus avoiding the modeler’s intervention. The method is based on the idea that each dimension of the seeds has its own standard deviation. Also the method assumes that the best region to search the value of the proposals is located around their seed, but should not include another seeds. To obtain the seed 𝑘′closest to the 𝑘th seed, the distances from each seed to all the remaining seeds are calculated. Note that the closest seed to another is the seed with the minimum-strictly positive distance to it, because otherwise a zero distance would imply repetition of the seed, which cannot be properly considered as a seed. In this work, the L2-distance is used for calculating distance 𝚥(⋅)between two seeds as follows: 𝚥𝑘−𝑘′=[𝑑 ∑ 𝑞=1 ⎛⎜⎜⎝ 𝜃(𝑘) 𝑗(𝑞)−𝜃(𝑘′) 𝑗(𝑞)  𝜃𝑗(𝑞)⎞⎟⎟⎠ 2]1∕2 (7) where 𝜃(𝑘) 𝑗(𝑞)and 𝜃(𝑘′) 𝑗(𝑞)are the 𝑞th-components of the 𝜃-values of seeds 𝑘 and 𝑘′, respectively, where 𝑞= 1,…, 𝑑, and 𝑘= 1,…, 𝑘′,…, 𝑁𝑃0. The term  𝜃𝑞is the range of the seed parameter values in the 𝑞th-dimension, which is calculated as follows:  𝜃(𝑞)=‖max{𝜃(𝑘) (𝑞)}𝑁𝑃0 𝑘=1 − min{𝜃(𝑘) (𝑞)}𝑁𝑃0 𝑘=1 ‖(8) where the subscript 𝑗has been avoided for simpler notation. Then, for each of the seeds of the 𝑗th simulation level, the standard deviation in each dimension is obtained as the component of the distance from the seed to its closest seed in that dimension, as follows: 𝜎(𝑘) (𝑞)=𝜃(𝑘) (𝑞)−𝜃(𝑘′) (𝑞)(9) Remark 1. Since the 𝜃-values are not equally distributed within the parameter space 𝜣, hence the normalization in Eqs. (7) through  𝜃𝑞 enables the parameter components 𝑞= 1,…, 𝑑 to equally contribute to the distance 𝚥. Otherwise, if  𝜃𝑞were discarded, then the effect of each dimension 𝑞on 𝚥would be larger as the 𝜃-values in that dimension increase. Remark 2. A limitation of Eq. (9) is that in case that the closest seed 𝑘′has the same value as seed 𝑘within any dimension 𝑞, then 𝜎(𝑘) 𝑞will result in zero leading to a lack of exploration for this seed in this dimension. To avoid this critical situation, and in order for the proposals to be fairly distributed in all dimensions around any 𝑘th seed, then the proposal standard deviations for each seed in all dimension are normalized, as follows: 𝜎(𝑘) 𝑞=𝜎(𝑘) 𝑞  𝜃𝑞 (10) Finally, the standard deviations are recalculated based on the maximum normalized standard deviation among all dimensions. Mathematically: 𝜎(𝑘) 𝑞= 𝜃𝑞⋅max{𝜎(𝑘) 𝑞}𝑑 𝑞=1 (11) The application of this automatic scaling method for the proposal standard deviation results in an algorithm that gives a step-wise improvement to the intermediate posterior regions which are broadly sampled at increasingly higher resolutions. This differs from existing versions of the ABC-SubSim algorithm and MCMC algorithms for Subset Simulation [37], where the whole possible region is sampled uniformly resulting in repeated sampling of unnecessary areas, thus increasing the amount of dependent samples and decreasing the acceptance rate. The advantages of using the proposed method for calculating standard deviation are shown in the illustrative example given in Section 5. 4. Petri Net model similarity by ABC-SubSim This section introduces the method that allows a reduced PN mimics a bigger reference one with quantified uncertainty. The proposed PN reduction method implements the ABC-Subsim inference method explained above to allow the reduced PN, referred to as (𝑟), to closely approximate the behavior of a reference one, denoted here as (𝑠), so that its signal outputs from a selection of places can be closely replicated over time by those from (𝑟). The method assumes that a reduced PN (𝑟), whose behavior wants to be mimicked with respect to (𝑠), is available. At this standpoint is important to note that (𝑟)must be chosen in a manner that retains the capability to sufficiently reproduce the reliability estimators along with the maintenance activities to be modeled whilst reducing the complexity of the reference one (𝑠). Indeed, the reduced PN must have the capacity to reproduce the key outputs of (𝑠)and incorporate any behaviors for optimization of the system under consideration, thus it is recommended that the net structure of (𝑟)is decided based on the knowledge of the modeler. Irrespective to the previous comments, this work is not focused on the methods for net reduction of PNs and the reader is referred to [49] for rules for optimal simplification of PN models. Thus, let us assume that (𝑟)is available and that its behavior is intended to reproduce the reference one based on a selection of comparison places 𝐏(𝑟)= {𝑝(𝑟) 𝑢}𝑛𝑟 𝑢=1 for the reduced PN, and 𝐏(𝑠)= {𝑝(𝑠) 𝑣}𝑛𝑠 𝑣=1 for the reference one, where 𝐏(𝑟)and 𝐏(𝑠)are not necessarily of same size. The basis of comparison for both PN models is their marking at a specific prediction time 𝑛restricted to the referred comparison places, namely 𝐌(𝑟) 𝑛=(𝑀(𝑟) 𝑛(1),…, 𝑀(𝑟) 𝑛(𝑢),…, 𝑀(𝑟) 𝑛(𝑛𝑟))and 𝐌(𝑠) 𝑛= (𝑀(𝑠) 𝑛(1),…, 𝑀(𝑠) 𝑛(𝑣),…, 𝑀(𝑠) 𝑛(𝑛𝑠))for the reduced and reference PNs, respectively. Next, let us assume that a subset of the 𝑛(𝑟) 𝑡transitions of (𝑟)are selected such that their delay times are parameterized through a set of uncertain parameter vector 𝜃= {𝜃1, 𝜃2,…, 𝜃𝑑} ∈ 𝜣⊂R𝑑whose prior information is known through their component-wise prior PDF 𝑝(𝜃1),…, 𝑝(𝜃𝑑), so that 𝑝(𝜃) = 𝑝(𝜃1)⋅𝑝(𝜃2)⋅…⋅𝑝(𝜃𝑑). Therefore, the marking of the reduced PN depends upon the uncertain parameters 𝜃, i.e., 𝐌(𝑟) 𝑛=𝐌(𝑟) 𝑛(𝜃)henceforth the challenge relies on obtaining the set of 𝜃-values which makes 𝐌(𝑟) 𝑛(𝜃) ≈ 𝐌(𝑠) 𝑛,∀𝑛. Reliability Engineering and System Safety 222 (2022) 108365 5 M. Chiachío et al. Algorithm 1 Petri net similarity algorithm by ABC-SubSim method Inputs: 𝑝0∈ [0,1] {𝑃0= 0.2is recommended}. 𝑁, {number of samples per intermediate level}; 𝜖, {tolerance}; 𝑛𝑟,{number of output places to compare}; 𝑑, {dimension of the 𝜃-parameter space}. Algorithm: 1: Sample [(𝜃(1) 0, 𝑥(1) 0),…,(𝜃(𝑛) 0, 𝑥(𝑛) 0),…,(𝜃(𝑁) 0, 𝑥(𝑁) 0)], where (𝜃, 𝑥) ∼ 𝑝(𝜃)𝑝(𝑥|𝜃) 2: Set 𝑗= 0,𝜖𝑗= ∞ 3: while 𝜖𝑗> 𝜖 do 4: 𝑗←𝑗+ 1 5: for 𝑛∶ 1,…, 𝑁 do 6: for 𝑢∶ 1,…, 𝑛𝑟do 7: Evaluate the centroid linkage distance 𝑑(𝑢)according to Eq. (12) 8: end for 9: Obtain 𝑚𝑎𝑥{𝑑(𝑢)} and 𝑚𝑖𝑛{𝑑(𝑢)} and normalize the distances according to Eq. (13) to obtain 𝜌𝑗(𝜃(𝑛)) 10: end for 11: Renumber [(𝜃(𝑛) 𝑗−1, 𝑥(𝑛) 𝑗−1), 𝑛 ∶ 1,…, 𝑁]so that 𝜌(1) 𝑗⩽𝜌(2) 𝑗⩽…𝜌(𝑁) 𝑗 12: Fix 𝜖𝑗=1 2(𝜌(𝑁𝑃0) 𝑗+𝜌(𝑁𝑃0+1) 𝑗) 13: if 𝜖𝑗⩽𝜖then 14: End algorithm 15: end if 16: for 𝑘= 1,…, 𝑁𝑃0do 17: Select as a seed (𝜃(𝑘),1 𝑗, 𝑥(𝑘),1 𝑗)=(𝜃(𝑘) 𝑗−1, 𝑥(𝑘) 𝑗−1)∼𝑝(𝜃, 𝑥|(𝜃, 𝑥) ∈ 𝑗) 18: (Requires 𝑘 > 1) Calculate the distance 𝚥𝑘−𝑘′, according to Eq. (7), where 𝑘′= 1,…, 𝑁𝑃0, 𝑘′≠𝑘 19: Set 𝜃𝑘′ 𝑗as the closest seed by min(𝚥𝑘−𝑘′), excluding repeated seeds 20: for q= 1,..., 𝑑do 21: Obtain 𝜎(𝑘) (𝑞)=𝜃(𝑘) (𝑞)−𝜃(𝑘′) (𝑞) 22: Obtain  𝜃𝑞according to Eq. (8) and normalize 𝜎(𝑘) 𝑞=𝜎(𝑘) 𝑞∕ 𝜃𝑞 23: end for 24: Set component-wise proposal standard deviation 𝜎(𝑘) (𝑞)= 𝜃𝑞⋅ max {𝜎(𝑘) 𝑞}, 𝑞 = 1,…, 𝑑 25: Run MMA [44] with proposal standard deviation 𝜎(𝑘) (𝑞)to generate 1∕𝑃0 states of a Markov chain lying in 𝑗:[(𝜃(𝑘),1 𝑗, 𝑥(𝑘),1 𝑗),…,(𝜃(𝑘),1∕𝑃0 𝑗, 𝑥(𝑘),1∕𝑃0 𝑗)] 26: end for 27: Renumber [(𝜃(𝑘),𝑖 𝑗, 𝑥(𝑘),𝑖 𝑗) ∶ 𝑘= 1,…, 𝑁𝑃0;𝑖= 1,…,1∕𝑃0]as 28: [(𝜃(1) 𝑗, 𝑥(1) 𝑗),…,(𝜃(𝑁) 𝑗, 𝑥(𝑁) 𝑗)] 29: end while Under the Approximate Bayesian computation perspective, the latter implies obtaining the posterior PDF of the model parameters and their corresponding marking response (𝜃, 𝑥) ∈ as 𝑝𝜖(𝜃, 𝑥|𝑦), where 𝑥=𝑥(𝜃)≐𝐌(𝑟) 0−𝑛, the simulated marking of the reduced PN (𝑟)from initial time to 𝑛, and 𝑦refers to the data which, in this case, correspond to the marking of the reference one (𝑠)from 0to 𝑛, i.e, 𝑦≐𝐌(𝑠) 0−𝑛. As with any ABC method, the marking similarity 𝑥(𝜃) ≈ 𝑦is established under a tolerance value 𝜖, along with through a metric 𝜌(⋅)and a summary statistic 𝜂(⋅). Thus, the next step to apply the ABC inference method is to define 𝜌(⋅)and 𝜂(⋅)to properly compare the response markings from (𝑟)and (𝑠). Several distance measures are available in the literature, and among them, the centroid linkage distance [50] was revealed as the most suitable one since it involves quantifying the distance between the mean values of two comparable clusters. To apply this specific metric, the amount of comparison places from 𝐏(𝑟)and 𝐏(𝑠)should be equal, thus the comparison is carried out through the marking of the respective places in a one-by-one basis. Since 𝑛𝑟=𝑛𝑠in this particular case, the indexes 𝑢and 𝑣can be fused into a single one (e.g., 𝑢), thus 𝑥(𝑢, 𝜃)and 𝑦(𝑢)denote the markings up to time 𝑛for the reduced and reference PN model, respectively. Let us now split the markings 𝑥(𝑢, 𝜃) and 𝑦(𝑢)into a number of overlapping sub-vectors 𝑥ℎ(𝑢, 𝜃)and 𝑦ℎ(𝑢),ℎ= 1,…, 𝐻, each of same length and with an overlapping percentage of 40%. Then, the centroid linkage distance applied to places 𝑝(𝑟) 𝑢and 𝑝(𝑠) 𝑢, Fig. 1. High Level Petri net used as illustrative example. The figure provides indication of the arc weights and production probabilities of probability transition 𝑡1. referred to here as 𝑑(𝑢)≡𝑑(𝑢, 𝑢) = 𝑑(𝑣, 𝑣),𝑢= 1,…, 𝑛𝑟, compares the mean values of the sliding markings under a L2-norm, mathematically: 𝑑(𝑢) = (𝐻 ∑ ℎ=1 (𝑥ℎ(𝑢, 𝜃) − 𝑦ℎ(𝑢))2)1∕2 (12) where 𝑥(⋅)and 𝑦(⋅)refer to the mean of the ℎth sliding marking of the reduced and reference PN models, respectively. Note that calculating the metric 𝑑(𝑢)from each comparison places 𝑝(𝑟) 𝑢and 𝑝(𝑠) 𝑢, where 𝑢= 1,…, 𝑛𝑟, will result in a vector of metrics whereby an overall metric 𝜌=𝜌(𝜃)can be obtained as follows: 𝜌(𝜃) = 𝑛𝑟 ∑ 𝑢=1 𝑑(𝑢) − min{𝑑(𝑢)}𝑛𝑟 𝑢=1 ∥ max{𝑑(𝑢)}𝑛𝑟 𝑢=1 − min{𝑑(𝑢)}𝑛𝑟 𝑢=1 ∥(13) where 𝜌only depends upon the model parameters 𝜃through 𝑑(𝑢), as can be seen from Eq. (12). Note also that in Eq. (13), the individual distances 𝑑(𝑢)are normalized such that their influences for the overall metric 𝜌vary within [0,1]. At this standpoint, the ABC-SubSim method can be applied whereby posterior samples of the model parameters 𝜃and their corresponding marking response 𝑥=𝑥(𝜃), are obtained as output. A pseudo-code implementation of the PN model similarity method by ABC-SubSim has been provided as Algorithm 1. To produce the initial subset (see Step 1 of the Algorithm 1), samples of the model parameters are drawn from the prior PDF 𝑝(𝜃)by which (𝑟)can be simulated via Monte Carlo simulation to obtain 𝑥=𝑥(𝜃)over a time duration 𝑛. The process continues with the production of subsets until the metric falls below a specified threshold 𝜖, meaning that ABC-SubSim has reached the posterior region of the model parameters. Once a posterior region in the parameter space has been found, the error associated with the reduced PN can be straightforwardly identified from the ABC-SubSim tolerance value 𝜖. In summary, the proposed method updates a number of uncertain parameters from a proposed reduced PN such that the signal outputs over time of a reference PN can be closely replicated. In this method, care must be taken on which parameters to update within the reduced PN, since it is possible to fit multiple parameters involving high model complexity and over-fitting to the reference PN. On the contrary, using few of them can imply very low accuracy of the marking from the reduced model with respect to the reference one, thus obtaining meaningless results from a reliability engineering point of view. Through the evidence calculation presented in Eq. (6), a trade-off between complexity and accuracy can be established, not only at the level of possible parameterizations of a reduced PN model, but also at the level of different potential reduced PN schemes with respect to a reference one. This will be explored in the example presented in Section 6. Reliability Engineering and System Safety 222 (2022) 108365 6 M. Chiachío et al. Fig. 2. Parameter updating of the 𝜃1and 𝜃2parameters from the illustrative Petri net example of Section 5. Top panels show the results when using the classical proposal standard deviation method, whereas the bottom ones present the results using the proposed adaptive method for scaling the standard deviation. The shaded areas drawn using yellow color indicate the proposal’s areas. 5. Illustrative example This section gives an example to illustrate the use of the proposed method on a simple PN architecture. The example also illustrates the advantages of using the proposed adaptive sampling method for the ABC-SubSim algorithm. To this end, the PN structure shown in Fig. 1 is considered again. The reader is referred to Section 2.1 for a description of the PN structure. Here, all arcs’ weights are set to 1, and the production probabilities from 𝑡1to {𝑝2, 𝑝3, 𝑝4}are chosen so that {𝜋12 = 0.9999, 𝜋13 = 0.00009, 𝜋14 = 0.00001}, respectively. The delay times of transitions 𝑡1, 𝑡2,𝑡3and 𝑡4are chosen as {𝜏1= 40, 𝜏2= 0.05, 𝜏3= 0.1, 𝜏4= 10}, respectively, which are used to obtain synthetic reference data. Then, {𝜏1, 𝜏2}are assumed to be unknown and updatable parameters using the proposed method. Their prior distributions are assumed as uniforms, as follows: 𝜃1∼[10,100],𝜃2∼[0.01,1]. In this example, the delay time of 𝑡1is chosen to be much greater than that of 𝑡2in the reference case, to show the ability of the model in solving critical cases where the algorithm has to infer parameters whose components belong to very different ranges of values. In this example, the places 𝑝1,𝑝2and 𝑝5were considered as comparison places thus their markings 𝐌(𝑠) 𝑛and 𝐌(𝑟) 𝑛(𝜃)for the reference and reduced PN model, respectively, were used to compute the ABCSubSim metric. The centroid linkage metric with one sliding region for each signal was used for places 𝑝1and 𝑝2. However, since that 𝑝5is a counter place, only its last marking was considered when calculating its metric. The original ABC-SubSim algorithm and the one with adaptive sampling presented in Section 2.3, were used to comparatively examine their performance. Each ABC-SubSim run employed 𝑁= 1000 samples per simulation level and a maximum of eight simulation levels. Fig. 2 shows plots of posterior samples of the 𝜃1, 𝜃2parameters produced at the second, fifth, and eighth simulation level of the ABCSubSim algorithm. The plots also depict the seeds with their accepted and rejected proposals, along with their proposals’ areas. The proposal area for each seed is represented by a yellow color shade ellipse whose dimension-wise semi-axis is twice the standard deviation in that dimension, in order for the ellipse to contain 95th percentile of the proposals. The results show that the posterior samples lie around 𝜃1= 40 and 𝜃2= 0.05 for both, the classical proposal standard deviation method and the proposed method to adaptive select the aforementioned algorithm hyper-parameter. Top panels (obtained using the classical standard deviation method) reveal that the amount of rejected samples highly increases as the algorithm progresses towards deeper simulation levels, which is a well-known drawback of the original ABC-SubSim algorithm. These rejected samples, which are drawn as repeating samples, conform the majority of samples within the final posterior PDF, thus significantly decreasing the acceptance rate and the statistical information of the posterior. However, posterior samples obtained using the proposed adaptive method for the ABC-SubSim are drawn around the seeds’ area leading to an increased acceptance rate of new samples even when the range of values is very small, as shown in the 6th simulation level from Fig. 2 (bottom-right panel). As overall consequence, higher resolution for the ABC-SubSim inference is achieved by adopting lower tolerance values 𝜖, since the statistical noise of the sampling method of the algorithm is reduced at a minimum. 6. Case study The previous section illustrated a case where a SPN model with two uncertain parameters was compared by itself using synthetic data obtained by taking the parameters as known. In this section, the proposed Bayesian methodology is applied to infer the model parameters of a set of reduced PN models based on the response of a larger reference PN model of an engineering system. 6.1. Description of the engineering system and its maintenance model The engineering system to be modeled corresponds to a two-pumps system to support the pressure of a lube oil supply facility which operates in a continuous regimen, as depicted in Fig. 3. The system contains a main pump, that is usually in operation except when a failure occurs whereupon a drop in pressure would be observed since the supply flow would be less than the demand flow. In this case, an auxiliary pump will then come into play to recover the pressure of the system. Further details on this system along with their physical data can be found in [51, Chapter 7]. In this case study, an insufficient system pressure is classified as a failed system state, which can occur only if both pumps fail. The Reliability Engineering and System Safety 222 (2022) 108365 7 M. Chiachío et al. Fig. 3. Engineering diagram of the two pumps system considered as case study. The diagram shows a duty/standby arrangement of the pumps with indication of the auxiliary components of the system, like valves, starters, pressure sensor controllers, and connections. failure of a pump occurs when there is insufficient outflow from the pump to maintain the oil pressure, and it is revealed through a pressure sensor over the main pump. In this example, the back-up pump is assumed to have unrevealed failures, since it is usually not in operation. Both pumps are assumed to be inspected on a periodic basis whose results can reveal either a failure or a degraded condition in any or both pumps. Here, degradation can happen at a different rate in each pump and is defined as a state where the pump has sufficient capacity to support the system, but maintenance is required. Mathematically, the pump degradation is assumed to be modeled as a two parameter Weibull PDF representing the time to failure, namely (12,1) for the main pump, and (30,2) for the auxiliary one. After inspection, a maintenance action is activated if both pumps are in the failed state or, one of them is in the failed state whilst the other is in a degraded state. Another option is that any of the pumps was in a degraded state, then an opportunistic maintenance action is activated under a probability, which in this case study equals 0.2. After maintenance, each pump is returned to the working state, thus recovering the overall system to normal operation. The following assumptions are also adopted for modeling the maintenance of the system under study: (a) The rate of degradation of the back-up pump does not change when the main pump fails; (b) The pressure sensor associated with the main pump will always immediately reveal a failed state of the main pump; (c) The components in the remaining parts of the system are always in the working state; (d) Inspection is always successful in discovering a degraded or failed state; (e) If a degraded state is identified but not rectified in time, then a failure can occur; (f) If there is a failed state in either of the pumps, then maintenance occurs as soon as a maintenance team is available. 6.2. Petri net maintenance model acting as reference In this section, a HLPN is proposed to model the system behavior along with its inspection and maintenance actions, as previously described in Section 6.1. The HLPN model is depicted in Fig. 4, and it will be used as the reference net for the reduction by ABC-SubSim numerical experiments that will be shown further below. In this HLPN, places 𝑝1and 𝑝2model the overall working and failed states of the system, respectively. Since the pumps are assumed to operate in parallel, thus both must be in the failed state for a system failure to occur. In such a case, transition 𝑡1immediately marks place 𝑝2. Then, if any of the pumps returns to the working state, the marking of place 𝑝2is removed and place 𝑝1is marked again by either transition 𝑡2or 𝑡3. For the back-up pump, places 𝑝3, 𝑝4, and 𝑝5represent the working, degraded, and failed states, respectively. When place 𝑝9is marked then an inspection action is enabled which lasts until place 𝑝8is marked. If degradation is discovered, then place 𝑝7is marked and maintenance is requested. Analogously, when place 𝑝6is marked then there is a discovered failure of the pump and, as well, maintenance is required. The main pump states are represented on the right hand side of the graph. In this part of the model, places 𝑝10, 𝑝11 and 𝑝12 corresponds to the good,degraded, and fail states, respectively. Similarly to the architecture described for the back-up pump, a token in place 𝑝15 indicates that inspection is enabled for the main pump, which is underway until place 𝑝14 is marked. If a degradation is discovered, then place 𝑝13 is marked whereupon maintenance is requested. The remaining nodes of the model govern the maintenance scheduling for the pumps. Specifically, place 𝑝18 and transition 𝑡25 model the availability of maintenance resources which produce a token in 𝑝19 indicating that maintenance is possible. When place 𝑝17 is marked then maintenance is requested for a pump in the system due to an identified failure in one of the pumps. A token in 𝑝22 indicates that maintenance is available to any of the pumps to improve a degraded state, whereas, in an analogous way, a token in 𝑝23 then indicates maintenance to repair a failure state of any pump. If neither of the pumps are in the failed state, then the available maintenance resources can be assigned to repair a degraded state, which is modeled by the marking of place 𝑝21 after enabling of transitions 𝑡26 and 𝑡27. In this model, inspections are assumed to be carried out every 6 weeks with a maintenance interval of 3 weeks, taken both terms in average. The rest of parameters are described in Table 5. Initially, the system starts in the working state represented by one token in places 𝑝1,𝑝3,𝑝10,𝑝9,𝑝15 and 𝑝18. The Monte Carlo simulation has been used to simulate the response of the reference PN model described above using 𝑁= 4000 samples for every stochastic transitions of the model. The results, which have been produced in week units, are shown using dark line in Fig. 6 for places 𝑝1and 𝑝2, and also in Fig. 9 for places 𝑝5and 𝑝12. 6.3. Reduced Petri net maintenance models This section presents three different reduced SPN models, namely (𝑟) 1,(𝑟) 2, and (𝑟) 3, to approximate the reference one, (𝑠), defined in the previous section. The reduced models are parameterized with a set of uncertain parameters that will be described in the following subsections, which together with their prior information, conform the model classes 1,2, and 3. The uncertain parameters are updated using the proposed Bayesian methodology using the marking of the places 𝑝1,𝑝2,𝑝5, and 𝑝12, from (𝑠), as data. Algorithm 1is applied using 𝑃0= 0.2,𝑁= 4000 samples per simulation level, and a tolerance 𝜖= 0.035. Gaussian densities are used as proposal PDFs for the MCMC sampling of Algorithm 1using adaptive scaling of their standard deviations, as proposed in Section 3. 6.3.1. Reduced model 1 Fig. 1 shows the first reduced model (𝑟) 1which is based on a simple net structure of two transitions, 𝑡1and 𝑡2, and two places, 𝑝1 and 𝑝2. This reduced SPN extracts the overall behavior of both pumps in the system along with their maintenance and inspection strategies to solely provide an estimate of the system states, namely working or failed state, represented by places 𝑝1and 𝑝2, respectively. The transition 𝑡1governs the time to total system failure whereas 𝑡2stands for the Reliability Engineering and System Safety 222 (2022) 108365 8 M. Chiachío et al. Fig. 4. Reference HLPN model of the case study presented in Section 6. The dashed rectangles indicate nodes representing the state of the main and back-up pumps. The remaining places and transitions model maintenance scheduling on the system level. The dark small rectangles indicate instant transitions whereas the blue text labels provide some explanatory information about key places. Table 2 Description of the transitions from the PN model shown in Fig. 1. Transition Type Distribution Parameters Description of action 𝑡1Stochastic Normal 𝜃1,20 Initiates system failure 𝑡2Stochastic Normal 𝜃2,1Initiates system repair time to repair of the system. Both transition times, namely 𝜏1and 𝜏2(expressed in week units), are modeled as stochastic transitions assumed as Gaussians, i.e., 𝜏1∼(𝜃1,20) and 𝜏2∼(𝜃2,20), where 𝜃1, 𝜃2are considered as the uncertain parameters of this reduced model (Table 2 provides an overview about the transitions). The prior PDFs of the model parameters are modeled as uniform PDFs within the [0,200] and [0,18] intervals, namely 𝑝(𝜃1|1) = (0,200) and 𝑝(𝜃2|1) = (0,18), respectively. The metric values are obtained by comparing the marking of places 𝑝1and 𝑝2of (𝑠)with respect to places 𝑝1and 𝑝2in the reduced one. Fig. 5 shows the resulting posterior PDF of the uncertain parameters 𝜃1, 𝜃2after the application of algorithm 1. The values within the posterior region that represent the Maximum a posteriori (MAP) were selected as point-valued estimates of 𝜃that make the reduced model better match the behavior of the reference one. These values are 𝜃1,𝑀𝐴𝑃 = 36.9, and 𝜃2,𝑀𝐴𝑃 = 1.0in weeks units, for transitions 𝑡1 and 𝑡2, respectively. The MAP values were then used to comparatively reproduce the marking sequence of the output places of (𝑟) 1during a500 weeks time. Additionally, the 10th and 90th percentiles of the results were calculated after simulating the reduced model using all the 𝜃-samples from the posterior region given in Fig. 5. The results are shown in Fig. 6 for reference places 𝑝1(left panel) and 𝑝2(right one). These results demonstrate that the inferred reduced model (𝑟) 1 evaluated using 𝜃𝑀𝐴𝑃 is able to capture the throughput of tokens passing through the comparison places 𝑝1and 𝑝2with good accuracy. Note also from Fig. 6 that the overall response of (𝑟) 1when using all possible model parameters values taken from the inferred posterior PDF, under-predicts the reference model response since there are some Reliability Engineering and System Safety 222 (2022) 108365 9 M. Chiachío et al. Fig. 5. Scatter plot representation of posterior samples in the 𝜃-space as ABC-SubSim output for the first reduced SPN model. On the diagonal, kernel density estimates are shown for the marginal posterior PDF of 𝜃1and 𝜃2at different simulation levels, differentiated using a color scale indicated at the right side of the plot. Table 3 Description of the transitions from the SPN model shown in Fig. 7. Transition Type Distribution Parameters Description of action 𝑡1Symbolic None – Fires instantaneously 𝑡2Symbolic None – Fires instantaneously 𝑡3Symbolic None – Fires instantaneously 𝑡4Stochastic Gaussian 𝜃1,2Initiates failure (back-up pump) 𝑡5Stochastic Gaussian 𝜃2,1Initiates repair (back-up pump) 𝑡6Stochastic Gaussian 𝜃3,1Initiates failure (main pump) 𝑡7Stochastic Gaussian 𝜃4,1Initiates repair (main pump) outliers which are not well covered within the 10th and 90th probability band. The latter can be viewed as a consequence of the extremely reduced net structure of (𝑟) 1with respect to the reference one, since this model summarizes the behavior of the entire system into a single cycle (failure–repair) of transitions for both pumps. 6.3.2. Reduced model 2 The second reduced model slightly increases the complexity of the previous one by considering simple sub-nets for each of the pumps in the system, specifically {𝑝3, 𝑝4, 𝑡4, 𝑡5}for the back-up pump, and {𝑝5, 𝑝6, 𝑡6, 𝑡7}for the main one, as depicted in Fig. 7. These sub-nets are complemented with three transitions {𝑡1, 𝑡2, 𝑡3}that govern the overall system states, namely working and failed states, represented by places 𝑝1and 𝑝2, respectively. Similarly to the reference model, the pumps are operating in parallel such that both must be in the failed state for a failure to occur. In this reduced structure, each of the pump’s sub-net only consists of the working and failed states, corresponding to 𝑝3and 𝑝4for the backup pump, and to 𝑝5and 𝑝6for the main one, respectively. Transitions 𝑡4, 𝑡5and 𝑡6, 𝑡7govern the times to failure and times to repair of the respective pumps, which are assumed to be modeled as Gaussian PDFs as follows: 𝜏4∼(𝜃1,2); 𝜏5∼(𝜃2,1); 𝜏6∼(𝜃3,1); 𝜏7∼ (𝜃4,1), where 𝜃= {𝜃1,…, 𝜃4}are considered as the uncertain model parameters of this model. The component-wise prior PDFs of these uncertain parameters are known and considered as uniforms, given by: 𝜃1∼[0,100],𝜃2∼[0,20],𝜃3∼[0,80] and 𝜃4∼[0,20], respectively, where these values have been considered sufficiently large to allow the method to obtain the posterior values within them. Table 3 gives further information about the referred transitions. In this model the marking sequences of places 𝐏(𝑠)= {𝑝1, 𝑝2, 𝑝5, 𝑝12}in the reference model were compared with the sequences of places 𝐏(𝑟) 2= {𝑝1, 𝑝2, 𝑝4, 𝑝6}in the reduced model respectively, since they hold the same interpretation. The data for the Bayesian inference are gathered here as the vector 𝑦=[𝑀(1)(𝑠) 0−500, 𝑀(2)(𝑠) 0−500, 𝑀(5)(𝑠) 0−500, 𝑀(12)(𝑠) 0−500], where 𝑀(⋅)(𝑠) 0−500 refers to the marking of a specific place of the reference PN for a 500 weeks-period. Algorithm 1was used to update the uncertain parameters 𝜃so as to make (𝑟) 2(𝜃)working similar to (𝑠).Fig. 8 shows the output of the Bayesian inference of 𝜃with respect to the data 𝑦taken from (𝑠). These results disclose MAP values for 𝜃given by 𝜃1,𝑀𝐴𝑃 = 11.96, 𝜃2,𝑀𝐴𝑃 = 6.07, 𝜃3,𝑀𝐴𝑃 = 4.38 and 𝜃4,𝑀𝐴𝑃 = 0.002 weeks, respectively. As with the previous section, the reduced model was simulated using the MAP values to comparatively obtain the marking sequences of the output places from (𝑟) 2with respect to those from (𝑠), and the results are shown in Fig. 9. This figure also represent the 90th and 10th percentile for the average markings sequences obtained through the use of all posterior samples from Fig. 8. The results show that the simulations from the reduced model can capture the response from the reference one when performing a posterior statistical averaging from the posterior samples. However, when observing the model response from the 𝜃𝑀𝐴𝑃 (in green), there are periodic peaks in the marking sequence of the reduced model that cannot be seen in that of the reference one (in gray). This reflects that this reduced model has limited capacity to reliably reproduce the behavior of the reference model, perhaps due to lack of extra nodes to better capture the system behavior including not only failure and recovering, but also maintenance actions. The latter will be tested in the following section. 6.3.3. Reduced model 3 The third reduced model (𝑟) 3considered in this study is based on an extension of the reduced model (𝑟) 2presented before by adding three extra nodes to the sub-nets on the lower left hand side and the lower right hand side of the graph, as depicted in Fig. 10. The new nodes correspond to {𝑡5, 𝑝5, 𝑡7}for the sub-net of the back-up pump and to {𝑡9, 𝑝8, 𝑡11}for the one of main pump, and allow each pump operation being modeled through three states: working, degraded, and the failed state, respectively. As with the previous reduced models presented above, the pumps operate in parallel such that both must be in the failed state for a system failure to occur. The system failure is manifested through a token in 𝑝2, otherwise the system is in the working state, represented by one token in 𝑝1. In this reduced model, the times to degradation and failure are governed by the transitions 𝑡4and 𝑡5for the back-up pump, and by transitions 𝑡8and 𝑡9for the main one. Transitions 𝑡6, 𝑡7and 𝑡10, 𝑡11, model the respective maintenance after a degraded and failed state are revealed in the back-up and main pump, respectively. The time delays associated to the aforementioned maintenance transitions are assumed to be modeled as Gaussian PDFs as follows: 𝜏6∼(𝜃1,1); 𝜏7∼ (𝜃2,1); 𝜏10 ∼(𝜃3,1); 𝜏11 ∼(𝜃4,1), being 𝜃= {𝜃1,…, 𝜃4}the uncertain model parameters. Similarly to the previous reduced model, the component-wise prior PDFs of these uncertain parameters are considered as uniforms, which in this case are given by: 𝜃1∼[0,20], 𝜃2∼[0,40],𝜃3∼[0,40] and 𝜃4∼[0,4], respectively. Note that the prior interval for 𝜃4is significantly smaller than the rest, and the reason is due to shorter repair times expected for the main pump under failure. Under these settings, Algorithm 1has been applied to obtain the posterior PDF of 𝜃that makes the response of model (𝑟) 3similar to the that from the reference one, (𝑠). The model response is evaluated through the marking sequences of places 𝐏(𝑟) 3= {𝑝1, 𝑝2, 𝑝5, 𝑝8}in comparison to the marking from places 𝐏(𝑠)= {𝑝1, 𝑝2, 𝑝5, 𝑝12}, in the reference model. Likewise the previous cases, the data for the Bayesian inference are given by the marking of the comparison places 𝐏(𝑠)over 500 weeks-period simulation, i.e, 𝑦=[𝑀(1)(𝑠) 0−500, 𝑀(2)(𝑠) 0−500, 𝑀(5)(𝑠) 0−500, 𝑀(12)(𝑠) 0−500]. The ABC-SubSim algorithm led to the posterior inference Reliability Engineering and System Safety 222 (2022) 108365 16 M. Chiachío et al. [13] Taleb-Berrouane M, Khan F, Amyotte P. BayesIan stochastic Petri nets (BSPN)-a new modelling tool for dynamic safety and reliability analysis. Reliab Eng Syst Saf 2020;193:106587. [14] Andrews J, Fecarotti C. System design and maintenance modelling for safety in extended life operation. Reliab Eng Syst Saf 2017;163:95–108. [15] Murata T. Petri Nets: Properties, analysis and applications. Proc IEEE 1989;77:541–80. [16] Murad CA, de Melani AHA, de Michalski MAC, Netto AC, de Souza GFM. Estimation of operational and maintenance tasks influence on equipment availability through Petri net modeling. In: 2020 annual reliability and maintainability symposium (rams). 2020, p. 1–6. http://dx.doi.org/10.1109/RAMS48030.2020. 9153582. [17] Latsou C, Dunnett SJ, Jackson LM, et al. A new methodology for automated Petri net generation: Method application. Reliab Eng Syst Saf 2019;185:113–23. [18] Knezevic J, Odoom E. Reliability modelling of repairable systems using Petri nets and fuzzy lambda–tau methodology. Reliab Eng Syst Saf 2001;73(1):1–17. [19] Chiachío M, Chiachío J, Prescott D, Andrews J. A new paradigm for uncertain knowledge representation by plausible Petri nets. Inform Sci 2018;453:323–45. [20] Senderovich A, Rogge-Solti A, Gal A, Mendling J, Mandelbaum A, Kadish S, Bunnell CA. Data-driven performance analysis of scheduled processes. In: International conference on business process management. Springer; 2016, p. 35–52. [21] Senderovich A, Shleyfman A, M. Weidlich M, A. Gal A, Mandelbaum A. To aggregate or to eliminate? Optimal model simplification for improved process performance prediction. Inf Syst 2018;78:96–111. [22] Berthomieu B, Le Botlan D, Dal Zilio S. Counting Petri net markings from reduction equations. Int J Softw Tools Technol Transf 2020;22:163–81. [23] Bønneland FM, Dyhr J, Jensen PG, Johannsen M, Srba J. Stubborn versus structural reductions for Petri nets. J Log Algebr Methods Program 2019;102:46–63. [24] Padma V. A study on flexible manufacturing systems using petr net reduction. IOP Conf Ser: Mater Sci Eng 2021;1130. [25] Bourdil P-A, Berthomieu B, Dal Zilio S, Vernadat F. Symmetry reduction for time Petri net state classes. Sci Comput Program 2016;132:209–25. [26] Berthomieu B, Le Botlan D, Dal Zilio S. Counting Petri net markings from reduction equations. Int J Softw Tools Technol Transf 2020;22(2):163–81. [27] Berthomieu B, Le Botlan D, Dal Zilio S. Petri Net reductions for counting markings. In: International symposium on model checking software. Springer; 2018, p. 65–84. [28] Rogge-Solti A, Mans RS, van der Aalst WM, Weske M. Repairing event logs using timed process models. In: OTM confederated international conferences‘‘ on the move to meaningful internet systems’’. Springer; 2013, p. 705–8. [29] Leclercq E, Lefebvre D, Medhi S. Identification of timed stochastic Petri net models with normal distributions of firing periods. 13, 2009, p. 948–53, [30] Buchholz R, Krull C, Horton G. Reconstructing model parameters in partially-observable discrete stochastic systems. 2011, p. 159–74. [31] Li X-Y, Liu Y, Lin Y-H, Xiao L-H, Zio E, Kang R. A generalized Petri net-based modeling framework for service reliability evaluation and management of cloud data centers. Reliab Eng Syst Saf 2021;207:107381. [32] Yan R, Dunnett S. Improving the strategy of maintaining offshore wind turbines through Petri net modelling. Appl Sci 2021;11(2):574. [33] Yianni PC, Neves L, Andrews J. Accelerating Petri-net simulations using NVIDIA graphics processing units. European J Oper Res 2018;265:361–71. [34] Marjoram P, Molitor J, Plagnol V, Tavare S. Markov chain Monte Carlo without likelihoods. Proc Natl Acad Sci 2003;100(26):15324–8. [35] Sisson SA, Fan Y, Beaumont M. Handbook of Approximate Bayesian Computation. CRC Press; 2018. [36] Chiachio M, Beck J, Chiachio J, Rus G. Approximate Bayesian computation by subset simulation. SIAM J Sci Comput 2014;36:A1339–58. [37] Papaioannou I, Betz W, Zwirglmaier K, Straub D. MCMC algorithms for subset simulation. Probab Eng Mech 2015;41:89–103. [38] Beck JL. BayesIan system identification based on probability logic. Struct Control Health Monit 2010;17(7):825–47. [39] Jensen K, Rozenberg G. High-level petri nets: theory and application. Springer Science & Business Media; 2012. [40] Fearnhead P, Prangle D. Constructing summary statistics for approximate Bayesian computation: Semi-automatic approximate Bayesian computation. J R Stat Soc Ser B Stat Methodol 2012;74(3):419–74. [41] Beaumont M, Cornuet J-M, Marin J-M, Robert C. Adaptive approximate Bayesian computation. Biometrika 2009;96:983–90. [42] Dutta R, Schoengens M, Pacchiardi L, Ummadisingu A, Widmer N, Onnela JP, Mira A. ABCpy: A high-performance computing perspective to approximate Bayesian computation. 2020, arXiv:1711.04694. [43] Hadjidoukas P. High performance framework for Bayesian uncertainty quantification and optimization. 2021, URL: https://github.com/cselab/pi4u/tree/master/ inference (Accessed Feb 25, 2021). [44] Au S-K, Beck J. Estimation of small failure probabilities in high dimensions by subset simulation. Probab Eng Mech 2001;16:263–77. [45] Zuev KM, Beck JL, Au S-K, Katafygiotis LS. BayesIan post-processor and other enhancements of subset simulation for estimating failure probabilities in high dimensions. Comput Struct 2012;92:283–96. [46] Gelman A, Roberts G, Gilks W. Efficient metropolis jumping rules. Bayesian statistics, Vol. 5, Oxford University Press Oxford; 1996. [47] Vakilzadeh MK, Huang Y, Beck JL, Abrahamsson T. Approximate Bayesian computation by subset simulation using hierarchical state-space models. Mech Syst Signal Process 2017;84:2–20. [48] Barros J, Chiachio M, Chiachio J, Cabanillas F. Adaptive approximate Bayesian computation by subset simulation for structural model calibration. Comput Aided Civ Infrastruct Eng 2021. http://dx.doi.org/10.1111/mice.12762, (in press). [49] Pham D, Karaboga D. Intelligent optimisation techniques. 2nd ed.. SpringerVerlag Berlin Heidelberg; 2000. [50] Desa MM, Desa E. Encyclopedia of distances. Springer-Verlag Berlin Heidelberg; 2009. [51] Forsthoffer MS. More best practices for rotating equipment. ButterworthHeinemann; 2017. [52] Ciardo G, Trivedi KS. A decomposition approach for stochastic Petri net models. In: Proceedings of the fourth international workshop on petri nets and performance models. IEEE Computer Society; 1991, p. 74–5. [53] Woodside CM, Li Y. Performance Petri net analysis of communications protocol software by delay-equivalent aggregation. In: Proceedings of the fourth international workshop on petri nets and performance models. IEEE Computer Society; 1991, p. 64–5.