Monte Carlo simulation of a 3 binary-state component system in Python, part 1 of 3.
A quantum mechanics approximation.
Abstract
The content of this article focuses on an area of industrial engineering called RAMS (Reliability, Availability, Maintainability & Supportability) and on which is predictive maintenance the main target to achieve. So the methodology described here is an excellent resource to make an analysis of reliability in critical complex systems such as thermo-nuclear plants, oil&gas industry, and aerospace projects, among others.
Introduction
When it comes to simulating stochastic processes, Monte Carlo methods play a fundamental role. We define as stochastic processes all those whose realizations (execution of an event directly from the process) depend on time and, therefore, do not repeat themselves over time, in other words, a single realization of the process corresponds exclusively to each moment, it never will be repeated again. Stochastic processes are known as probabilistic processes also.
Monte Carlo methods have demonstrated to be useful for estimating solutions of integral equations such as those presented in solving real-world problems regarding reliability diagnostics.
With the application of Monte Carlo methods in solving the calculation of the reliability of a system, what we are doing is essentially the randomization of deterministic optimization methods, which gives rise to Randomized Computing. In my opinion, this is the key step, because randomization infers variance in some natural way and it makes up for math models, being more realistic.
Randomized Computing, a brief introduction
The algorithms that describe simulation models can be viewed as part of Randomized Computing. The introduction of these probabilistic algorithms has enabled numerous advances in engineering areas. a probabilistic algorithm is one that receives, among other input data, a source of random bits.
The three main basic advantages of using probabilistic algorithms are the following:
1.- The execution time or memory requirements of a probabilistic algorithm are typically less than those of the best known deterministic algorithm for a problem.
2.- They are usually easy to understand and implement.
3.- Randomizing a badly behaving deterministic algorithm in the worst case usually leads to a well-behaving probabilistic algorithm with a high probability of success.
At this point, it can already be intuited that one of the keys in probabilistic algorithms is to have an excellent engine for generating (pseudo)random numbers. The fact that countless commercial environments with pseudo-random number generators are available, makes us think that we should not worry about this matter, but nothing is further from the truth. Experience shows us that this matter must be treated with care, especially due to the complexity of the simulation models that are carried out today, hence this area is an active field of research.
One way to generate pseudo-random numbers is based on the use of physical mechanisms. For instance, we can highlight the white noise produced by electronic circuits, the count of radioactive particles emitted, the atmospheric noise, etc. However, the use of physical mechanisms for simulation is discouraged, since reproducibility is a desirable characteristic, so that we can replicate the experiments again and again under the same conditions. Therefore, we must look for algorithmic procedures for generating pseudo-random numbers. The idea, introduced by John Von Neumann, is to produce numbers that appear random, using the arithmetic operations of the computer and starting from an initial seed.
The following is the algorithm that the main programming languages have adopted when computing pseudo-random numbers. However, it should be noted that not all commercial generators produce good quality results and this can have negative consequences on the simulation process.
The Mersenne Twister algorithm, specifically the MT19937 for 32-bit and the MT19937–64 for 64-bit respectively, is usually used by languages such as Dyalog APL, Apache Commons, C++, CUDA Library, Microsoft Excel, GAUSS, GLib, GNU Multiple Precision Arithmetic Library, GNU Octave, GNU Scientific Library, Gretl, IDL, Julia, CMU Common Lisp, Embeddable Common Lisp, Steel Bank Common Lisp, Maple, MATLAB, Mathematica, Free Pascal, PHP, Python, R, Ruby, SAS, SageMath, Scilab, Stata y SPSS. It is based on the property of Mersenne primes that are of the type:
Mersenne prime numbers are a fundamental part of the perfect even numbers found out by Euler that are of the type:
How to transform pseudo-random numbers into stochastic variables
At this point, the aim is to transform the random numbers obtained into random variables from probabilistic distributions of interest for the problem under study. The method is applied to both discrete and continuous, univariate, and multivariate distributions (exponential, Weibull, multivariate normal, Wishart, Dirichlet, and so on). To do so we use the Inverse transformation method (also known as inversion sampling, the inverse probability integral transform, the inverse transform sampling, Smirnov transform, or the golden rule). See:
Inverse transformation method
The inverse transformation method states that if a random variable, x, has a cumulative probability function that admits inversion
then it is verified that the transformed variable
always follows a continuous uniform distribution. Hence, if we solve x as a function of U we have:
what if it is generated a pseudo-random variable, u, that follows:
it leads to:
with what we have pulled off random variables of the distribution of interest, as shown in the following figure:
The actual problem that the simulation engineers usually face and moreover a minimum condition to apply the inverse transformation method is to know the explicit form of the distribution function of the stochastic process. Sometimes the density function is known but the distribution function is not. This can be a tough problem, both analytically and computationally, hence for these cases, it is necessary to use Von Neumann’s Rejection method.
Von Neumann's Rejection method
This method is based on taking samples from simple and known distributions, but not from the distribution that governs the phenomenon that we want to simulate and which is much more complex to analyze, so random sampling will be used under certain conditions to calculate integrals defined with Monte Carlo simulation.
The main condition is that the probability density function of the known function is an envelope of the density function to be calculated, as shown in the following figure:
We must try to make the enveloping function fit the objective function as best as possible to minimize the rejected candidate variables and maximize the candidate variables to be accepted.
Simulation of a system by using the equations of quantum mechanics.
What it means “Transport” in a simulation?
Monte Carlo Simulation (MCS) offers a powerful tool for evaluating the reliability of a complex system, due to the modeling flexibility that it offers regardless of the type and size of the problem. The method is based on the repeated sampling of the realizations of the system configurations, which, however, rarely fail, so a large number of realizations must be simulated to achieve an acceptable precision in the estimated probability of failure, with large computation times. For this reason, efficient sampling techniques for system failure realizations are of interest, in order to reduce computational effort. For this purpose, techniques such as Subset Simulation and Line Sampling have been developed, which make it possible to improve the efficiency of the Monte Carlo simulation in estimating the probability of system failure.
Here are some definitions to better understand the mathematical theory:
Plant: system of N components properly connected.
Component: a subsystem of the plant that can remain on different (multi) exclusive states (nominal, fault, standby mode, …). Stochastic transitions from state to state occur for stochastic times.
Plant State for t: represents the set of states in which the N components remain on t. The plant states are marked by a scale that lists all the possible combinations of all the component states.
Plant Transition: When any of the plant components make a state transition, the plant is said to have made a transition. The moment in which the plant carries out the n-th transition is called t_n and the state in which the plant will enter is called k_n.
Mission Time: it is the lifetime for which the system is expected to operate correctly and it is the one that will provide us with a metric on the reliability parameter (R) of the system. Once this time has been passed, if the system fails, it is mainly due to a low level of attention to the maintainability parameter (M) of the complete set of components.
As can be seen in the following figure, we can represent the Monte Carlo tests by means of a two-dimensional trajectory between the states reached by the system and the entry phases (time) in the transitions. Specifically, the image represents three simulations, two of which have entered a system failure state and only one has exceeded the mission time. Technically each of these representations is a simulation of random walks with the Markovian property.
Just as an example, if we have a system composed of 10 components with binary states (On / Off), there are 2 ^10 = 1024 possible states of the system.
To carry out the simulation, we must first know the equations that will govern the stochastic behavior of the system. Such equations are presented below and simulate the ‘transports’ of the states of a plant.
Quantum equations to simulate a plant states transport
The equations that govern stochastic changes come from quantum mechanics. Specifically, the following are used:
Free flight kernel: the conditional probability of a transition at t in dt given that the preceding transition occurred at t’ and that the state before the transition occurred was k’.
Collision kernel: the conditional probability that the plant enters state k, given a transition that occurred at time t when the system was in state k’.
Transport kernel: the combination of both previous probabilities form the conditional probability that the plant enters a state k at time t in dt given that the preceding transition occurred at t’ and that the state before the transition occurred was k’.
The following figure summarizes the transport kernel schematic:
It can be observed that the dependent variable of the transport kernel is a probability density function of a transition of the system at t giving rise to the entry into state k, and it is called the Boltzmann integral equation, which is solved by Neumann’s serial expansion:
For now, that’s enough…
In subsequent parts of this article, I’ll talk about how to compute all the critical states for a system using a technique named Minimal Cut Sets, and in the last part, I’ll talk about how I’ve programmed all this stuff in Python.
Coming soon: Part 2: System failure states
If you liked this, please support clicking the 💚 below so other people will see this here on Medium.
