Modeling and simulation in engineering

This review article will explore the innovative and popular theme of engineering modeling and simulation, predominantly in the manufacturing industry and cybersecurity world, citing severe challenges, advantages and time‐ and budget saving solutions and its future. The power of simulation is not an exaggeration but an understatement. The favorable outcomes since the advent of digital computers and software revolution could not have been achieved, especially without the multiple benefits of statistical simulation, which underlies the widespread use of modeling and simulation in engineering and sciences, stretching from A (Astronomy) to Z (Zoology). This refers not only to research findings in verifying a certain piece of theory, such as that of the recently discovered Higgs Boson, but in testing new products to innovate new discoveries so as to make our universe a more peaceful place by modeling and simulating the future projects and taking precautions before disasters occur. The review explores a cross section of engineering modeling and simulation practices illustrating a window of numerical examples. WIREs Comput Stat 2013, 5:239–266. doi: 10.1002/wics.1254


INTRODUCTION AND BRIEF HISTORY TO SIMULATION AND MOTIVATION
C omputer modeling and simulation (M&S), as programs or networks of computers mimicking the execution of an abstract model of many natural systems from physical and life sciences to social and managerial sciences, and primarily engineering, have become an integral part of digital experimentation. M&S proves useful to estimate the performance of complex engineering systems when too prohibitive for analytical solutions. A simulation is defined as the reproduction of an event with the use of scientific models. A model is a physical, mathematical, or other logical representation of a system, process, or phenomenon. Time-independent static Monte Carlo (MC) or conversely dynamic Discrete Event Simulation (DES) to manage events in real time for engineering applications will be reviewed. Taxonomywise, simulated computer models may be stochastic or deterministic, and dynamic or static, and discrete or continuous.
Modern computer simulation developed in parallel with the rapid-growth of computer use during the development of the Manhattan Project in WWII to nondestructively model and simulate the nuclear detonation before it was destructively dropped on Hiroshima and Nagasaki in Japan in 1945. Therefore, the history of simulation is interesting and intriguing. Some earliest pioneers can be observed in Ref. 1 Lord Rayleigh in 1899 showed that a onedimensional random walk without absorbing barriers could provide an approximate solution to a parabolic differential equation. In 1908 W.S. Gosset (with a nickname, Student) used experimental sampling to help him towards his discovery of the distribution of the correlation coefficient and to bolster his faith in his so-called t-distribution. 2 A.N. Kolmogorov in 1931 showed the relationship between Markov stochastic processes and certain integro-differential equations. Stanislaw Ulam at Los Alamos labs performed simulation in 1945 during the WWII in the bomb-building Manhattan Project before proposing the Teller-Ulam thermonuclear weapon design. Ulam suggested first the 'Russian Roulette' and 'splitting' methods, for evaluating complicated mathematical integrals for nuclear chain reactions that later led to the systematic Monte Carlo methods by von Neumann, Metropolis and others. John von Neumann explored the behavior of neutron chain reactions in fission devices using statistical sampling methods in 1948 (such as the acceptance-rejection method) employing the newly developed electronic computing techniques. Neumann proposed the agent-based Von Neumann Machine, 3 a theoretical machine capable of reproduction following detailed instructions to copy itself. Ulam suggested a machine as a collection of cells on a grid. The idea intrigued von Neumann, who created the first of the devices later termed, cellular automaton. 4 John Conway constructed the well-known Game of Life, 5 operated in a virtual world in the form of a two-dimensional checkerboard. A team headed by N. Metropolis using the ENIAC Computer in 1948 carried out what's contemporarily known as modern Monte Carlo calculations.
Computer simulation has been widely used in engineering systems to validate the effectiveness of tentative decisions regarding a new plan or schedule, or its outcomes, without actually experiencing the actual conditions, which could in actuality cost more resources or partial to full destruction such as in the simulation of the nuclear bomb. In a book titled, Simulation Engineering, by Jim Ledin in 2001, 6 the author outlines his twofold purpose as follows: (1) simulation engineering (SE) is the application of engineering discipline to the development of good simulations.
(2) Similarly, SE occurs when simulations become part of an engineering process when applied as tools to develop better products and test processes with a greater efficiency for different types of complex embedded systems. The latter purpose (2) is the subject matter of this review article. The IEEE June 2012 Spectrum issue, emphasized that the Modeling and Simulation effect is a creative and time-saving topic of interest ranging from automotive engineering of hybrid vehicles to finding solutions to treating nuclear waste, and upgrading the nuts and bolts of the Electrical Power (Smart) Grid and moreover, supercomputing research. 7

GENERIC THEORY-CASE STUDIES ON GOODNESS OF FIT FOR UNIFORM NUMBERS
A formal scientific theory of simulation, to verify a validated model so as to mimic a physical or a social system, does not exist in terms of conventional math-statistical theorems and their subsequent proofs. However, heuristic modeling formalisms at an advanced level for engineers through cellular automaton for Monte Carlo and Discrete Event Simulations are studied by Zeigler et al. 8 (Ch. 4), although these formalisms do not lend themselves to easy algorithmic implementations for practicing engineers or scientists as this review article purports to. Moreover, the fundamental process of verifying sequences of uniform deviates from an associated generator where Ho: Uniform Random Sequence is (quasi) random versus Ha: Sequence is not random, is an accepted technique. For instance, χ 2 tests, such as those by Leven and Wolfowitz 9 and Knuth, 10 are popularly well-accepted math-statistical scientific practices to theorize the verifiability of uniform random numbers essential to the realm of statistical simulation. In order to clarify the validation of the above stated Ho: Random Sequence versus Ha: Not Random Sequence, the commercial JAVA code's uniform number generator will be tested for randomness, as illustrated in a series of screenshots from Figures A1-A9 in Appendix A by using Stewart's JAVA program to implement Knuth's technique. 11 The results show that by law of large numbers for only n ≥ N ≈ 50,000; E(θ) → 0.5 with probability 1, for θ ∼ Uniform(0, 1) from the uniform number generator imbedded in the Java code, 'Ho: Random Sequence' is not rejected. Therefore n = 50,000 runs is a new standard for attaining quasirandomness; not 5000 anymore as practiced in 1980s.

WHY CRUCIAL TO ENGINEERING-MANUFACTURING AND CYBER DEFENSE ISSUES
The power of simulation is prevalent as the audiovisual Ref 12 favorably explains certain topics related to production and manufacturing engineering. In 'Modeling and Simulation in Manufacturing and Defense Systems Acquisition', the Board on Manufacturing and Engineering Design (BMED) emphasized the importance of modeling and simulation in not only making the right decisions but also incurring fewer expenses. 13 Similarly, the Wychavon (UK) council has adopted manufacturing industry's simulation model to reduce waste and improve performance. 14 Since the US manufacturing industry is challenged by increased global competition and price erosion, one can benefit from manufacturing simulation to eliminate bottlenecks, enhance lean manufacturing, optimize capacity planning and optimize production output. In an Annotated Discrete Event Simulation Bibliography, there exist 325 articles on manufacturing simulation as cited in Ref 15. A certain bibliography displays 112 publications on 'Load Models for Power Flow and Dynamic Performance Simulation' by the IEEE Transactions on Electric Power Systems. 16 On the contrary, there are fewer simulation studies in cybersecurity-and-defense-related theoretical and applied research. In their 2007 article as titled, 'Cyber Attack Modeling and Simulation for Network Security Analysis', the authors Kuhl et al. discuss a simulation modeling approach to represent computer networks and intrusion detection systems (IDS) to efficiently simulate cyber-attack scenarios in order to test and evaluate cyber security systems. 17 You-Tube-based audio-visual Cybersecurity Simulation roundtable 18 underlines the power of simulation in cybersecurity scenarios. Under 'War game reveals US lacks cyber-crisis skills' in a war game, 19 sponsored by a nonprofit group and attended by former top-ranking national security officials, laid bare that the US government lacks answers to such key questions. Former Clinton press secretary Joe Lockhart said that people would be scared by the simulation but he added, '. . .that's a good thing. ' 20 considers modeling and simulation of individual components and systems toward assessment of security risk, in addition to his publications where theoretical models are confirmed using Monte Carlo and Discrete event simulation runs. [21][22][23][24] Further, certain manufacturing-and cybersecurity-themed examples will be reviewed through working details of how the modeling should be validating the physical model and the subsequent simulation computationally verifying the solutions accurately and cost effectively. These reviews are the tips of the iceberg, as industries will continue to design and discover new products and services by M&S. Figure 1 displays the interaction between the process of building a model by focusing on the interplay between (1) experimental results, (2) simulation results, and (3) theoretical predictions as displayed in Ref 25. A favorable example of this interplay is presented in a recent WIREs article titled Cloud Computing 26 (Figure 8, p. 55), which displays an experimental scenario for a trivial Cyber Cloud.
On the other hand Ref 26 (Figure 9, p. 57) outlines the Markovian theoretical predictions followed by the simulation results for the same scenario of two 1-GB units serving a constant load of 1.5 GB for 13 cycles. The resulting availability of this small Cloud: (1) 0.307 for Experimental, (2) 0.305 for Simulation after 1 million runs, or trials and (3) 0.331 for Markov Theoretical, allowing a negligible error content, which diminishes to less than 3% as the size of the experiment increases from a few to many hundreds of units. In the event of large cyber CLOUDS such as those with 398 units, the authors showed that the experimental approach was infeasible, and theoretical result was not mathematically tractable. However, supercomputerdriven programming worked for days regarding the basic two-state assumptions, crunching 2 398 ( 10 100 ) Markov states to 93.8% reliability. DES result was a satisfactorily comparable 90.5%.

A CROSS SECTION OF MODELING AND SIMULATION ISSUES IN MANUFACTURING
Simulation use in production is not new. For the sake of a few examples, various authors from 25 years ago published articles on simulating flexible manufacturing systems (FMS), machine utilizations and production rates, and modeling of Automated Manufacturing Systems (AMS). [27][28][29] Given the advances in pervasive computing regarding communication networks, as well as recently popularized large scale cloud computing in cyber networks; quantitative risk assessment of a manufacturing unit and their network availability have become challenging tasks. An often overlooked fact is that many real-life grid units such as routers or servers in cyber physical systems to the manufacturing assemblies in automotive or avionics, etc. and the intricate telephony networks (wired or wireless), and water-supply networks or hydroelectric dams, do not operate in an idealized simple setting of either full or zero capacity. This fact therefore necessitates the inclusion of degrees of derated (inbetween UP and DOWN states) capacity. Because of lack of closed-form solutions in the three-state model including DERATED as opposed to that of the conventional UP-and-DOWN dichotomous two-state model, a summary of three or multistate system inferential analysis will be reviewed by using Monte Carlo simulations. This process will employ the empirical Bayesian principles to estimate the full and derated availability probability distributions. The historical failure and repair data, or operating (full or derated) and nonoperating hours, as the input data, will  32 The primary difference between the above listed three references and this review article is the empirical Bayesian treatment of the three states to estimate their probability distributions by Monte Carlo simulations based on the Sahinoglu-Libby probability density function, originally derived independently by both Sahinoglu and Libby in 1981. [33][34][35][36][37] The closest among these three articles, that is, by Rao et al. 30 uses only four transition rates in a three-state Markov model whereas Sahinoglu's model uses all six transitions. 38,39 However, this review article's simulation approach is even more powerful and flexible as the application can be extrapolated based on identical principles to four or more states, whereas reference by Rao et al. 30 deals solely with differential equations limited in scope. Others, Lins et al. 31 and Shah et al., 32 are on slightly different but not identical topics; all of which do not employ modeling and simulation techniques or generate closed form statistical probability density function (p.d.f.) expressions and derivations.

Modeling and Simulation of Multistate Production Units and Systems in Manufacturing
Most research articles or books on reliability theory are devoted to traditional binary reliability models allowing for only two possible states for a system and its components: perfect functionality or complete failure. However, many real-world systems are composed of multistate components which have different performance levels and several failure modes with varying effects on the entire system performance. Such systems are called multistate systems (MSS). 40 Examples of MSS are cyber systems where the unit performance is characterized by the data processing speed or server gigabyte capacity and similar to electric power systems, where the generating unit performance is depicted by its generating capacity. In the electric power supply system of generating facilities, each generator can function at different levels of capacity with a given probability. This may result from the outages of several auxiliaries such as pulverizers, water pumps, fans, boilers, etc. Billinton  The power of simulation once again flexes its muscle as a favorable exit out of this theoretical impasse. The Monte Carlo technique remains the only available feasible way to solve the proposed three-state problem, whose math-statistical closed-form solution does not actually exist. This is mainly because the three-state Markov model's random variables' (UP, DOWN, DER) probability distributions cannot be derived through math-statistical transformations due to mathematical intractability and lack of sufficient statistical theory. The probability density function of the Forced Outage Rate (FOR) was earlier analyzed in a textbook by the primary author, who designated that the Sahinoglu-Libby (SL) probability model can be used if certain underlying assumptions hold. 20,[33][34][35] Libby and Novick independently have studied multivariate generalized beta distributions for utility assessment; however their analysis was for only two-states similarly, not for multistate hence the term, G3B: Generalized 3-Parameter Beta. 36,37 The failure and repair rates were taken to be the generalized gamma random deviates where the corresponding gamma shape and scale parameters, respectively, were not equal. The two-state SL density was shown to default to that of a standard two-parameter beta density function when the shape parameters are identical. The stochastic method proposed was superior to estimating availability by dividing total uptime by exposure time. Examples had shown the validity of this method to avoid over-or underestimation of availability when only small samples or insufficient data exist for the historical life cycles of units. In this article, however, additionally we shall also review, similar to the two-state SL, a computational three-state simulation model. Because of the infeasibility of closed-form solutions, the analysis will be carried out using Monte Carlo simulations, obeying the Bayesian principles similar to Chapter 5 of the author's textbook. 20 In studying large capacity production units, it is necessary to consider the probabilities associated with one or more forced derated states rather than accepting the unit being either available or unavailable, according to Billinton. 44 Following the Monte Carlo simulations, analytical p.d.f.s of the multiple states will be approximated using their associated moments.

Two-State Sahinoglu-Libby Probability Model of Production Units (Closed-Form Solution)
In using the distribution function technique, the p.d.f. of FOR = q = λ/(λ + μ) is obtained first by deriving its cumulative density function Then, taking its derivative to obtain g Q (q) as per Eqs (5A.1)-(5A.18) in Appendix 5A, on pp. 26 (1) the usual two-parameter beta p.d.f. is obtained. An alternative original derivation of the same p.d.f. termed under generalized multivariate beta distribution is given by Libby in 1981 and1982. 36,37 The expression in Eq. (1) can also be reformulated in terms of SL(α = a + c, β = b + d, L = (β 1 /β 2 )), as follows: (3) distribution is not necessarily so, and can be skewed positively or negatively, depending on L > 1 and L < 1 respectively, because the mode, skewness, and kurtosis of SL random variable now also depend on L. For 0 < L < 1, the SL p.d.f. stays below the plot of the related standard Beta near zero but crosses the latter to become the greater of the two p.d.f.s at a point 36,37 : The reverse action holds true for L > 1 with the same crossing point, y 0 . The major drawback to the distribution is that there is no closed form for finite estimates of the moments. The moment generating function for the univariate SL distribution is an infinite series. 36,37

Three-State Sahinoglu Probability Model of Production Units (Monte Carlo Simulation)
In studying large capacity generation (power) or production (or cyber-physical) units, it may be necessary to consider the probabilities associated with one or more forced derated-outage states as in multistate, rather than considering the unit as being either available or unavailable. [40][41][42]44,45 In summary, there are gray areas or in-between capacities which are called derated or degraded states. However in this review article, we will only consider a single derated state rather than multiple ones, which may well exist in practice such as in 50%, 60%, or 75% derated capacity. But now, we have not only full-FOR but also derated-FOR (or DFOR), that will be equal to the total derated operating time over the total exposure time. That is, DFOR = DER time/(UP time + DER time + DOWN time). It is also well documented that any calculated FOR or DFOR is not only a constant but also a specific single realization of its random variable. 20 The probability density function of the FOR by empirical Bayesian analysis was identified in Section Two-State Sahinoglu-Libby Probability Model of Production Units (Closed-Form Solution) to be the Sahinoglu-Libby (SL) probability density, where certain underlying assumptions hold. However, we shall review above and beyond a traditional closed-form two-state SL; namely, a three-state SL where the transition rates are gamma distributed (see derivations in subsections of Three-State Sahinoglu Probability Model of Production Units (Monte Carlo Simulation)). Let us examine and review the following state space diagram in Figure 3 by Billinton from his textbook. 44 Let λ = transition rate from UP (fully operational) to DOWN (forced outage) state. Let μ = transition rate from DOWN to UP state; δ = transition rate from UP to DER (partially forced outage) state; β = transition rate from DER (partially forced outage) to UP state; α = transition rate from DER to DOWN state; γ = transition rate from DOWN to DER state. Using Figure 3 given, by changing to the Greek variables from the Latin originals (a-f) cited in the same reference 44 (p.156, Fig. 4.2); the time-dependent but steady state probabilities of occupying one of the three states are given as follow from (5)-(7), assuming negative exponential densities for each state's sojourn time, which will converge to: where DENOMINATOR = μβ + μγ + αβ + λβ + λγ A closed form solution of the three-state SL is intractable and analytically impossible in this setting with six random variables, as compared solely to the two variables in Section Two-State Sahinoglu-Libby Probability Model of Production Units (Closed-Form Solution). We will therefore have to simulate the P(UP), P(DER) and P(DOWN) from Eqs (5)-(7) by generating the recursive Monte Carlo simulated deviates of the state transition rates. Empirical Bayesian analysis will be pursued through deriving first the conditional posterior densities of the six transition rates from subsections of Three-State Sahinoglu Probability Model of Production Units (Monte Carlo Simulation), and using random uniforms for generating the transitions that constitute the probabilities in Eqs (5)- (7). See Figure 4 for a sample draft scenarios to illustrate transitions of Figure 3.

UP-to-DOWN Failure Transition Rate (λ)
, for example, from x1 to w1, or x2 to w2 in Figure 4 Let a = number of occurrences of UP (operating) times before DOWN (recovery) going DOWN (recovery) for 'a' such occurrences. λ = full UP-to-DOWN rate. c = shape parameter of gamma prior for the full failure rate λ. ξ = inverse scale parameter of gamma prior for the full failure rate λ. Now let the failure rate, λ have a gamma prior distribution: The joint likelihood of the UP-time random variables is the joint distribution of data and prior becomes: Thus, the posterior distribution for the random variable λ is Note that x is a vector.

DOWN-to-UP Recovery Transition Rate (μ), for example, from y1 to x1, or y2 to z1 in Figure 4
Let b = number of occurrences of DOWN (recovery) times before UP (operating) before going UP for 'b' many such occurrences μ = full recovery (DOWN-to-UP) rate d = shape parameter of gamma prior for the full recovery rate μ η = inverse scale parameter of gamma prior for the full recovery rate μ Now let the full recovery rate, μ, have a gamma prior distribution: The joint likelihood of the DOWN-time random variables is The joint distribution of data and prior becomes: Thus, similarly skipping two intermediate steps, the posterior distribution for μ is which is also distributed as Note that y is a vector.

UP-to-DER Failure Transition Rate (δ), e.g. from z1 to u1, or z2 to s2 in Figure 4
Let o = number of occurrences of UP times before e = shape parameter of gamma prior for the UP-to-DER failure rate δ.
= inverse scale parameter of gamma prior for the UP-to-DER failure rate δ.
Now let the UP-to-DER failure rate δ have a gamma prior distribution: Thus, similarly skipping two intermediate steps, the conditional posterior density of δ becomes: Note that z is a vector.

DER-to-UP Recovery Transition Rate (β), for example, from u1 to x2, or u2 to z2 in Figure 4
Let k = number of occurrences of DER times before UP Thus, similarly skipping two intermediate steps, the conditional posterior density of β: which is also distributed as Note that u is a vector.

DER-to-DOWN Failure Transition Rate (α); for example, from s1 to y2, or s2 to y3 in Figure 4
Let j = number of occurrences of DER failure times before DOWN s T = j 1 S i = total DER failure times before going DOWN for 'j' many such occurrences. S i ∼ αe −αS . α = DER-to-DOWN failure rate. g = shape parameter of gamma prior for DERto-DOWN failure rate α. ψ = inverse scale parameter of gamma prior for DER-to-DOWN failure rate α. Now let the DER-to-DOWN failure rate α have a gamma prior distribution: Thus, similarly skipping two intermediate steps, the conditional posterior density of α: which is also a Gamma [j + g, (s T + ψ) −1 ]. Note that s is a vector. (γ ), e.g. from w1 to s1, or w2 to u2 in Figure 4 Let p = number of occurrences of DOWN times before DER w T = p 1 W i = total DOWN times before going DER for 'p' many such occurrences. W i ∼ γ e −γ W . γ = DOWN-to-DER recovery rate. h = shape parameter of gamma prior for the DOWN-to-DER recovery rate γ . π = inverse scale parameter of gamma prior for the DOWN-to-DER recovery rate γ . Now let the DOWN-to-DER recovery rate γ have a gamma prior distribution:

DOWN-to-DER Recovery Transition Rate
Thus, similarly skipping two intermediate steps, the conditional posterior density of γ : which is also distributed as Gamma [p + h, (w T + π ) −1 ]. Note that w is a vector. Table 1 displays the input data as tabulated for the following example covering the first five episodes of six different sojourn times (see Figures 3 and 4). The cumulative probabilities of states are calculated by Monte Carlo Simulation method using input from Table 1 as follows in Tables 2-4 for UP, DER and DOWN states in 100, 1000, and 10,000 simulation runs, respectively. Figures 5-7 using Eqs (5)-(7) will convert these tabulations into cumulative frequency plots utilizing the six transitions of Figure 3 in subsections of Three-State      (9)- (24). Plots shown in Figure 8 are the extrapolated JAVA versions of the EXCEL applications in   That is, Mean (E) ≈ Median (M) where a spike or two for the Mode (m) will not violate the symmetric appearance, as evident in Figures 9 and 10. Note, Mean = E(q) and Median = q 0.5 if loss functions are assumed to be squared error and absolute error respectively, where mode is the maximum likelihood estimator. This follows from the fact that E(q −q) 2 , if it exists, is a minimum whenq = E(q), that is, the mean of the conditional (posterior) distribution of q. Then E(q) is the Bayes solution:

Statistical Simulation of Three-State Units to Estimate the Density of UP, DOWN, and DER
Similarly according to Hogg and Craig 46 (p. 262), the median of the random variable Q is the Bayes estimator using an informative prior when the loss function is given as L(q,q) = |q −q|. If E(|q −q|) exists, thenq = q 0.5 minimizes the loss function, i.e., the median of the conditional posterior distribution of q. The median is resistant to changes. Then, q 0.5 or median of q, that is, q M is the Bayes solution as the 50th percentile or 0.5 quantile, or second quartile for q, as follows:

Example of SL Simulation for Modeling Network of 2 Two-State (UP-DN) Units
Given the following simplest series system of two identical components in Figure 12, whose default operational probability for each is P(UP) = 0.9 and hence P(System) = 0.9 2 = 0.81. We now force these units have their unavailability r.v. distributed with SL displayed as in Figure 2's l.h.s. column, where g Q (q) is formulated as follows:  (Table 5). Historical failure and repair data are given in Figure 2. The flat deterministic outcome is 0.9 2 = 0.81 whereas SLdistributional input-output relationship is unknown due to the closed-form derivation of the product of random variables being not available. Since Eqs (25) and (26) are not closed form solutions and tedious numerical integration is needed, Monte Carlo simulation can be the only solution for much larger networks if analytical tools are not available 45 (pp. 196-197, Figures 4 and 5) where analytical integration becomes an impossible task.

Example of Another SL Simulation for Modeling Network of 7 Two-State (UP-DN) Units
For the sake of a convenient example, a feasible and probable seven-node complex architectural style is taken 20 (p. 254) with failure and repair history including the prior parameters displayed on the l.h.s. with 10 each ups and downs lasting 1000 and 111.11 h, respectively, in Figure 13. The author assumes for the hypothetical control architecture an identical SLdistribution for its unavailability as displayed in Table 6 employing historical data for its components simulated 1000 times in 100-tuples of networks. This means 100,000 simulation runs overall. The analytical result being unknown for a complex system as the 7-node network depicted in Figure 13, the resulting simulation is 0.785 as in Table 6.
In the June 2012 issue of the IEEE Computer with 'Computing in Asia' as the cover feature, an article titled 'Computing for the Next-Generation Automobile' displays three hybrid vehicle architectural styles: series, parallel and series-parallel, and then the (Toyota) Prius integrated THS II control architecture. It is mentioned that most vehicles today come with more than 50 embedded computer components, called electronic control units (ECUs) 7 (pp. [34][35].

A REVIEW OF MODELING AND SIMULATION IN CYBERSECURITY
Modeling and simulation (M&S) is a vital tool that can be leveraged for process improvement, and technology/capability development and evaluation. It is the process of designing a model of a system and conducting simulated experiments to preview and predict system behavior and evaluate optimal strategies for system operation. A short review of approaches will be covered in the world of cyber security on MC or DES. With the cyber security breaches rampant in the world, some of the most creative solutions to counteract these problems can be obtained by digital simulation faster, safer and cheaper than they can be resolved in the physical labs. In his related article, Rinaldi highlights M&S as a crosscutting initiative to increase the security of critical infrastructures. 48 Their Strategy states that modeling, simulation, and analysis must be employed to 'develop creative approaches and enable complex decision support, risk management, and resource investment activities to combat terrorism at home'. Rinaldi concludes that the multidisciplinary science of interdependent infrastructures is immature, and requires M&S to mature it, and adds that they are developing, among others, at Sandia Labs the following techniques. Aggregate Supply and Demand (What-if Analyses), Dynamic Simulations, and ABM, which at a macro level similar to cellular automata, is out of scope for this review. Some examples of M&S and DES, in the cybersecurity field will follow.

Monte-Carlo Value-at-Risk Approach by Kim et al. in Cloud Computing
Based on today's volatile market conditions, the ability to generate accurate and timely risk measures has become critical to operating successfully, and necessary for survival. Value-at-Risk (VaR) is a market standard risk measure used by senior management and regulators to quantify the risk level of a firm's holdings. However, the time-critical nature and dynamic computational workloads of VaR applications make it essential for computing infrastructures to handle bursts in computing and storage resources needs. This requires on-demand scalability, dynamic provisioning, and the integration of distributed resources. A VaR calculation will typically start after the end of the trading day, when market data and final positions have been verified. It must be complete, and updated risk numbers must be available, before the start of the next trading day. As the number and complexity of positions change, the computational requirements for the calculation can change significantly, however the completion deadline of the beginning of the next trading day remains fixed. Furthermore, as market conditions change, a firm may want to vary the number of Monte Carlo scenarios run (and thus the resolution of the calculation), which will add additional variability to the computation time. Specifically, the authors demonstrate how the Comet Cloud autonomic computing engine can support online multiresolution VaR analytics, a candidate for Cloud architecture by integrating of private and Internet cloud resources. 49

Monte Carlo and Discrete Event Simulations in Sahinoglu's Security-Meter (SM) Risk Model
Four examples will be studied regarding Monte Carlo (MC) and Digital Event Simulation in the field of Cybersecurity.

Example for Security Meter Risk Modeling and Simulation
Assume two vulnerabilities and two threats in a 2 × 2 × 2 set up as in Figure 14. 20,22 Let X (total number of cyber-attacks detected) = 360/year and let X 11 = 98, X 12 = 82, X 21 = 82, Let Y (total number of attacks undetected) = 10/year and let Y 11 = 2, Y 12 = 3, Y 21 = 3, When we keep Figure 14 in sight, we obtain the risk ratios and ECL (Expected Cost of Loss) as follow.
We place the estimated input values for the security meter in Figure 14 to calculate total residual risk. Therefore, once you build the probabilistic model from the empirical data, as above, which should verify the final results, you can forecast or predict any 'taxonomic' activity whether it is the number of vulnerabilities or threats or crashes as in Table 7. For the study above, the total number of crashes is 10 out of 370 total events, which gives a ratio of 10/370 = 0.0270 to verify the final results in Figure 14.
Using this probabilistically accurate model, we can predict what will happen in a different setting or year for a newly given explanatory set of data as in Table 7. If a clue suggests to us a future 1000 total episodes and 500 episodes of vulnerabilities of V 1 , then by the avalanche effect, we can fill in all the other blanks, such as for V 2 =500. Then  Table 7. If the asset is $2500 and the criticality constant is 0.4¸then the ECL (expected cost of loss) is demonstrated in Figure 14 Figure 15 and Figure 16 Given the Total Number of Attacks The analyst is expected to simulate a cybercomponent's (such as a server) tree-diagram 10 consecutive times from the beginning of the year (e.g., 1/1/2013) until the end of 1000 years (i.e., 12/31/3012) in an 8,760,000 h period, with a life cycle of crashes or saves for a total of 10 × 1000 = 10,000 simulation runs. The input data is tabulated in Table 7 to conduct the generation of random deviates. At the end of this planned time period, the analyst will fill in the elements of the tree diagram for a 2 × 2 × 2 security meter's tree diagram model as in Figure 14. Recall that the rates are the reciprocals of the means for the assumption of a negative exponential probability density function to represent the distribution of time to crash. For example, if λ = 98 per 8760 h, the mean time to crash is 8760/98 = 89.38 h. Use the input as in Table 7. 20,22 We observe a result of TRR = 0.0269 ≈ 0.027 in Figure 15.

Monte Carlo (Static Time-Independent) Simulation using Poisson p.d.f.
Using the identical information in Section Example for Security Meter Risk Modeling and Simulation, the analyst is expected to use the principles of Monte Carlo simulation to simulate the 2 × 2 × 2 security meter as in Table 7 and Figure 14 for 10 repeated trials. One employs the Poisson distribution for generating failure and repair rates for each leg in the tree diagram of the 2 × 2 × 2 model shown in Figure 14. The rates are given as the count of saves (repairs) or crashes (failures) annually. The necessary rates of occurrence for the Poisson distribution's random value generation were given in the empirical data in Table 7 above. For each security meter realization, get a risk value and average it over n = 10,000 in 1000 increments. When you average over n = 1000 runs, you should get the same value as in Figure 15. Using the same data, as projected, we get the same results in Figure 16 as in Figure 15. That is, TRR = 0.0269 ≈ 0.027. 20,24 Therefore, DES and MC results were identical to four decimals as expected.

Monte Carlo (Static Time-Independent) Simulation using a Continuous Uniform p.d.f.
In another cyber server setting, with three vulnerabilities (V 1 , V 2 , V 3 ) with four threat levels forV 1 , three threat levels each for V 2 and for V 3 ; the following upper and lower risk values for U(a,b) are assumed to be available for each vulnerability, threat and corresponding LCM variables, as follow in Table 8. Selected upper and lower uniformly distributed example values are very close to reduce variation. Theoretical derivations for TRR's mean and variance are by MAPLE software using are as follow. 24 M := 0.260432113, V := 0.0000144453852, as copied from MAPLE outcomes, are tabulated in Table 8:

DISCUSSION AND CONCLUSION
The power of simulation is evident from countless number of contemporary research works in addition to industrial and military undertakings to save time  Table 8. On the manufacturing or production front, the author in response to then-in-2007-unsolved homework question 5.5 on p. 256 from his Wiley textbook, 20 reviews the set of techniques to generate the multistate probability distribution model of an important pillar of trustworthiness, that is, availability. Namely, when the availability (or reliability) of a unit is at stake, and while the unit possesses three operational states with a derated state added beyond the usual two-state binary or dichotomous assumption, conventional applications do not suffice. Therefore, it is worth to review the fact that the primary difference between other related works 30-32 and author's empirical Bayesian treatment of the three states of a repairable hardware unit is to estimate the p.d.f.s of these three states by using Monte Carlo simulations. 38,39 The closest article to this one uses only four transition rates in a three-state Markov model whereas the reviewed Monte Carlo model uses all six transitions. 30 This reviewed statistical simulation approach is powerful and flexible, whereas Ref 30 deals with differential equations limited in scope. Other close references deal with different topics; however none use any simulation techniques. 31,32 It is currently infeasible to find closed form solutions for the random variables of UP, DER, and DOWN expressed by Eqs (5) 53 Currently, only deterministic probabilities can be calculated through Markov algebra, but not their probability densities. For example, a credit card is either closed (if less than a critical credit score), open (more than) or only conditionally usable for urgent cases (between lower and upper). The Bank actuarians may want to estimate the p.d.f.s of these three states to conduct statistical inference using customer-based empirical data by employing empirical Bayesian analysis. Multistate systems such as in the case of four multiple derated states representing electric power turbines, as cited in Refs 20 (p. 280, Figure 6.26) and 45 (p. 201, Figure 10) can be derived. These estimators for unit availability can further be propagated to simulate the source-target availability for troublesome complex networks.
Regarding the cyber security science and engineering issues however, implementation of modeling and simulation compared to manufacturing industry is fairly new progressing at an experimental stage. This fact is not only because of involvement of human life and death situations in adversity, as compared to accidental casualties in the production world, but also because of lack of theoretical and experiential data base dating back to only 1990s since the launch of public internet. The author, by following examples in this area proceeds with currently popular VaR technique by Kim et al. 49 and Security Meter and CLOUD simulation tools (CLOURA) by Sahinoglu et al. 26 Monte Carlo VaR is costly to execute; it does not incorporate cost comparisons when taking measures. Consider a medium size firm holding positions in 20,000 different financial instruments. Running a 100,000 simulation Monte Carlo VaR calculation requires generating 2 billion simulated instrument prices. With a conservative estimate of 10 milliseconds per pricing, this calculation requires more than 5500 h of processor time over an 8 h window. The capital cost of hardware plus the operational cost for data center space, power, cooling and maintenance makes this cost prohibitive to all but the largest firms. However, scalable CLOURA is a very fast algorithm, that is, it can simulate a CLOUD system with 430 servers for 1000 years in less than 4 min. 26 SM simulations as in Tables 7-9 and Figures 15-17 and are relatively fast and accurate, comparable to their analytical counterparts. [20][21][22] Overall M&S techniques abound, particularly face-saving in the case of theoretical impasses, and sometimes the only viable solutions in engineering and scientific applications. The multiples of positive results render M&S methods among the most useful and practical, as well as affordable algorithms of our time. If one day, humankind can make it to the surface of the red planet Mars, it will be possible because humans will have nondestructively travelled to Mars some tri-zillion times by riding on the cyber space Volume 5, May/June 2013 through digital simulation rather than on the outer space. The author contends that positive solutions will realize for cancer and currently incurable diseases by crunching computationally intensive and nonlethal M&S techniques. The application of M&S to engineering, cyberspace and health informatics, however, is not an easy task with much progress remains to be done. This overview also aims to provoke thoughts and stimulate ideas for such goals by exploring interdisciplinary avenues through M&S using supercomputing.
Finally, one exam question in a Cybersecurity M.S. program's midterm exam at AUM 54 asked, 'What would separate you in your future job if you took an M&S course, and others did not have any clue?' The following four responses in Appendix B were gathered from candidates invariably all with a military background, either on active duty or retired USAF. The responses, as quoted, did demonstrate an awakening of mind on the timely significance of Modeling and Simulation in cyberspace and reliability and security engineering.
Overview wires.wiley.com/compstats FIGURE A4 | Uniform numbers testing; Ho: random versus Ha: not random for 5000 runs. Ho is rejected. On the average, one out of 10 cycles of 5000 = 50,000 simulations will end up rejecting Ho: random. FIGURE A5 | Uniform numbers testing; Ho: random versus Ha: not random for 10000 runs. Ho is NOT rejected. FIGURE A6 | Uniform numbers testing; Ho: random versus Ha: not random for 10,000 runs. Ho is rejected. On the average, one out of 25 cycles of 10,000 = 250,000 simulations will end up rejecting Ho: random.