JP5761702B2 - Simulation apparatus, simulation method, and computer program - Google Patents

Description

  The present invention relates to a simulation apparatus, a simulation method, and a computer program for performing simulation.

  The Monte Carlo method (hereinafter referred to as “MC method”) is generally studied in a very wide range as a simulation method using random numbers. For example, in addition to cases where large-scale objects that are not easy to actually conduct experiments are targeted (for example, reliability prediction of structures), a large amount of experiment is required to find the target object (For example, computer simulation of chemical reaction), recently, the applicable fields such as computer graphics are very wide.

  The inventor of the present application broadly divides the MC method into a region integration type MC method and an event reproduction type MC method (for example, Non-Patent Document 1 and Non-Patent Document 2). The area integration type MC method is relatively easy to improve efficiency, and methods of improving efficiency include wide-ranging methods such as the Importance Sampling method and Latin Hypercube Sampling method.

  However, for problems where many parameters are involved in a complicated manner or problems where an analytical solution cannot be obtained, it may be necessary to fully study the description of the integration region.

  In addition, when evaluating the reliability of a structural system, the region integral MC method is often considered assuming a bilinear stress / strain relationship. There are cases where the relationship (or restoring force characteristics) must be assumed (buckling problem, building seismic response analysis problem, etc.), and it is difficult to apply the domain integral MC method.

  In such a case, it is convenient to accurately reproduce the attention event and obtain the required amount as an average value of the observed observation amount. Such a method is sometimes called a Hit-or-miss MC method in applied mathematics. However, since the Hit-or-miss MC method is also used for region integration type MC calculation, it is different from region integration type MC calculation. Separately, it is called the event reproduction type MC method. The event reproduction type MC method can obtain a stable solution in almost all cases, and is effective as a solution for a wide range of problems. However, the event reproduction type MC method is inefficient, and, except for some studies, the efficiency improvement method has not been studied much.

S.A.Matsuho, “Efficient hit-or-miss Monte Carlo methods simulating relevantevents for structural reliability evaluation”, Journal of the Society of Materials Science, Japan, Vol.54, No.3, pp.314-319 (2005). S. A. Matsuho, “Basic study on structuralreliability evaluation based on Hit-or-miss Monte-Carlo method”, Journal of Constructional Steel, Vol.16, No.1, pp.147-152 (2008).

  The present invention has been made in view of the above situation, and an object of the present invention is to provide a simulation apparatus, a simulation method, and a computer program capable of improving the efficiency of the event reproduction type MC method.

In order to solve the above problem, the simulation apparatus according to the present invention takes a value of 0 or 1 depending on a vector x representing n parameters depending on whether or not an event of interest occurs. A simulation apparatus for obtaining a probability P _f of occurrence of an event according to a probability density function f _X ( x ) indicating the likelihood of occurrence of an event of interest based on a state index function I ( x ) and the vector x , P _f is calculated by using the multiple definite integral form shown in the following (Equation 3).

  According to the simulation apparatus and the like according to the present invention, in particular, in the event reproduction type MC method, there is no need to formulate a limit state function, and there is an effect that it is suitable for simulation by a digital computer.

It is a block diagram for demonstrating an example of a structure of the simulation apparatus in this Embodiment. It is a schematic diagram for demonstrating the lower bridge type Warren truss bridge used for the numerical calculation. In the truss bridge of FIG. 2, it is a figure which shows the cross-sectional area of each member, and elongate ratio 1 / r. Table 2 shows the calculation result of failure probability P _f together with the number of trials based on the event reproduction type MC method, and table showing the result of reliability evaluation by applying the Crude MC method to the formula of the definite integral transformation method FIG. As a result of trying to apply simple stratified sampling without sampling of stratified cells (Table 4), for reference, we tried to apply the pseudo-random allocation (Event) method to the event reproduction MC method. It is a figure which shows the result (table 5). It is a figure which shows the result (table 6) examined about the case where the quasi-random number which is a kind of quasi-MC method is applied to a definite form.

Hereinafter, embodiments of the present invention will be described with reference to the drawings.
FIG. 1 is a block diagram for explaining an example of a configuration of a simulation apparatus according to the present embodiment.

  The simulation apparatus 100 includes a definite integral conversion unit 110, a probability estimation unit 120, and an estimation result output unit 130. The simulation apparatus of the present invention can also be realized by installing a computer program for performing processing as described below in a digital computer. The random number generation unit 150 generates a random number for MC method simulation and outputs the random number to the probability estimation unit 120. Of course, the random number generator 150 may be provided in the computer. The type of random number is not particularly limited. For example, in addition to a uniform random number such as Mersenne Twister, an embodiment using a quasi-random number that does not necessarily correspond to a “random number”, a low difference quantity sequence, or the like is also considered. It is done.

The function defining unit 140 defines a state index function I ( x ) and a probability density function f _X ( x ), which will be described in detail later, and sends them to the definite integral conversion unit 110. The state index function and the probability density function can be input to a computer in advance and stored in a storage unit such as a hard disk. However, a configuration in which the state index function and the probability density function are input every time a simulation is performed is also possible. The state metrics obtained as a result of a simulation function I (x), the configuration using the information of the probability density function f _{X (x)} is also possible. The definite integral conversion unit 110 converts the state index function and the probability density function into a definite integral form.

In the event reproduction type MC method, n random numbers η ⁽ⁱ⁾ = ( ⁾ generated according to a probability density function f _X ( x ) of a basic random variable X = (X ₁ , X ₂ ,..., X _n ) ^T. η ₁ ⁽ⁱ⁾ , η ₂ ⁽ⁱ⁾ ,..., η _n ⁽ⁱ⁾ ) Perform one trial using ^T , repeat this N times, and estimate the required amount as an average value. In the event reproduction type MC method classified by the present inventor, n random numbers η ⁽ⁱ⁾ are not uniform random numbers.

In the case of structural reliability evaluation, which will be described later, in each trial, it is observed whether the target system is safe. A state index that takes a value of 0 if it is safe and 1 if there is a possibility of damage. Formulate using function I ( x ). That is, the estimated value P ^ _f of the failure probability _Pf is given by the following (Equation 1).

Limit state function g (x) is the case of danger g (x) ≦ 0, since it is defined to be safe if g (x)> 0, the equation (1) state indicator function I (x) Can be determined from the value of the limit state function g ( x ), but in the present embodiment, the event of interest is accurately reproduced in the computer, and the event of interest occurs (breaks) in each trial. Whether or not to observe.

In general, the failure probability P _f is estimated by performing event reproduction type MC calculation according to the above (Equation 1), or by using a limit state function g ( x ) as shown in the following (Equation 2). Region integral MC calculation.

However, the method using (Equation 1) requires a large number of trials for the occurrence of a rare event, and the calculation efficiency is poor. In the method using (Equation 2), the function form of the limit state function g ( x ), Alternatively, the function value needs to be given, and considering the formulation burden, it is not suitable for simulation with a digital computer.

In the present embodiment, considering that the above (Equation 1) is a sample average of the state index function I ( x ), a so-called expected value definition formula is used to form a definite integral as shown in the following (Equation 3). Tried to convert. Specifically, the definite integral conversion unit 110 converts the state index function and the probability density function given from the function defining unit 140 into a definite integral form such as the following (Equation 3).

As described above, the state index function I ( x ) takes a value of 0 when the target system is safe and 1 when there is a possibility of damage in each trial. By converting to the form of definite integral in this way, it is not necessary to calculate and formulate the value of the limit state function g ( x ), and a simulation method suitable for simulation using a digital computer can be provided. . In addition, with respect to integration using the MC method, there is an advantage that it is possible to achieve further efficiency by applying effective methods that have been variously studied so far.

For example, if the above (Equation 3) is evaluated by the Crude (introductory) MC method, for example, uniform random numbers ξ ⁽ⁱ⁾ = (ξ ₁ ⁾ for n basic random variables required for one trial. ^(I) , ξ ₂ ⁽ⁱ⁾ ,..., Ξ _n ⁽ⁱ⁾ ) Reproduce the event of interest using ^T , observe the value of the state index function I ( x ), and try this N times From the following (Equation 4), the failure probability P ^ _f is estimated. The estimation result output unit 130 outputs an estimation result.

(Equation 4) is obtained by discretizing the definite integral form shown in (Equation 3) so as to be suitable for calculation by a digital computer, and {Π (b _j −a _j )} is the basic random variable n. It is the volume of the Hypercube (hypercube) determined by the domain in the dimensional space. As described above, when converted into a definite integral form as in (Equation 3), various efficiency improvement methods can be applied. However, the inventor of the present application has applied the limit state function g ( x The Latin Hypercube Sampling method, which does not depend much on the shape of the integrand, is considered preferable in order to make use of the feature of the method of the present embodiment in which it is not necessary to care about the function form of

(Example of numerical calculation)
Hereinafter, an example of actual numerical calculation by the simulation apparatus of the present embodiment will be described. FIG. 2 is a schematic diagram for explaining a lower bridge type Warren truss bridge used for numerical calculation. This truss bridge is bilaterally symmetric, and FIG. 3 is a diagram showing the cross-sectional area of each member and the slenderness ratio 1 / r (where l is the member length and r is the secondary radius of the cross section). The width is assumed to be 6.0 m. The steel material used is SM400.

In the reliability evaluation, a dead load and a live load are assumed as acting loads. The dead load is a fixed amount. In addition, it is assumed that one large trailer (two in total) acts simultaneously at the center of each vertical line as the most dangerous live load loading state of the road bridge. The weight per large trailer follows a normal distribution with an average of 154.35 kN and a variance of 10912.5 kN ² with reference to past traffic flow survey results, and the working stress is calculated in consideration of impact.

For the impact, the impact coefficient specified in the road bridge specifications was used as the definite value. Further, the member strength σ _R assumes a complete correlation between the members. The average value of the tensile strength is assumed to follow a normal distribution having σ _ta = 140 N / mm ² (allowable tensile stress degree) and a coefficient of variation of 0.02. For compression members,

σ _ca = 140 (when l / r ≦ 18)
σ _ca = 140−0.82 (l / r−18) (when 18 <l / r ≦ 92)
σ _ca = 1200000 / (6700+ (l / r) ² ) (when l / r ≧ 92)

Assuming a normal distribution having an allowable compressive stress degree σ _ca given by (1) as an average value and a coefficient of variation of 0.02. In order to simplify the analysis, if the acting stress σ _{s of} one member exceeds the allowable stress σ _R (assuming that it is at least disadvantageous even if it is not actually broken) Assume. At this time, the example can be considered as the weakest link system, the overall structural system if even a damaged state by one of the components is as to failure state, to assess the failure probability P _f of the structural system.

  Under the above conditions, the reliability of the truss structure was evaluated. The computing environment was Windows (registered trademark) XP on Pentium (registered trademark) IV CPU 2.80 GHz, the Fujitsu Fortran & C Academic Package V4.0 C language was used, and Mersenne Twister was used for uniform random numbers.

FIG. 4 shows the calculation result of the failure probability P _f along with the number of trials based on the event reproduction type MC method (see Equation 1), and the Crude MC method in the formula of the definite integral transformation method (see Equation 4). It is a figure which shows the table 3 which shows the result of having performed reliability evaluation by the method to apply.

In each table, the estimated value of P ^ _f the failure probability P _f, Time shows the required calculation time (in seconds). For accuracy of time measurement in the calculation processing system, less than 1 second is expressed as 0 second. For example, in Table 2, 0 second of 100,000 trials is about 0.58 seconds from the time required for 10 million trials. Can be estimated.

  In the table, “——” indicates that the number of trials is small and the calculation result is not obtained. Further, ρ represents a coefficient of variation of the estimated value P ^ f, and is expressed by the following (Equation 5) and (Equation 6).

From Table 2, it can be seen that as the number of trials increases, the coefficient of variation ρ of the estimated value P ^ _f decreases and the estimation accuracy improves.

Next, consider the results of Table 3. In the trial of Table 3, the upper and lower limits b _j and a _j of the integration were set as b _j = m _j + 5σ _j and a _j = m _j −5σ _j . However, the upper and lower limit values of the integration are not limited to this, and can be appropriately set depending on the application, purpose, and the like. Here, m _j and σ _j are the average value and standard deviation of each basic random variable X _j . From Table 3, a million times and later attempts, it can be seen that ^{_{P ^ f = 7.4 × 10 -4}} ( effective digit number 2 digits) is obtained. In other words, it has been found that by using the definite integral conversion method of the present embodiment, the solution converges with a smaller number of trials.

However, a clear tendency for the coefficient of variation ρ could not be confirmed. This may for example be a large error estimate, than if failure probability P if the estimated value of _f is larger ρ becomes smaller value (P _f is small, it is easy to estimate better when P _f is greater ) And the like. Therefore, when referring to ρ, it can be considered that it is not appropriate to determine whether the accuracy is good or not by simply comparing ρ.

  In the definite integral conversion method of the present embodiment, since the basic equation (Equation 3) is expressed by a definite integral, various effective integration methods can be applied. Although the above suggested the application of the Latin Hypercube Sampling method, this calculation example is a two-dimensional problem, and since the number of dimensions is small, a simple stratified sampling that does not sample stratified cells (Equation 3) Tried to apply to. For reference, an attempt was made to apply a pseudo-random assignment (Assigning-Pseudorandom-Numbers) method to the event reproduction MC method of (Equation 1). FIG. 5 shows the results. The results are shown in Table 4 and Table 5, respectively.

  For simple stratified sampling, see, for example, Takao Tsuda, “Monte Carlo Method and Simulation”, revised edition, Baifukan Co., Ltd., October 10, 1987, pp. 196 84-106. In application of simple stratified sampling to (Equation 3) (Table 4), the domain of each basic random variable (upper and lower limits of integration) was equally divided into 10 parts for simplicity.

As shown in Table 4, P ^ _f = 7.45 × 10 ⁻⁴ (3 significant digits) is obtained in 10 million or more trials. That is, it can be seen that the calculation with higher accuracy than the table 3 can be realized. However, the required calculation time is slightly longer than that of Table 3.

  In the case of Table 5, as compared with the result of Table 2, the required calculation time for the same number of trials was shorter, and a tendency to converge to the solution with a smaller number of trials was observed. However, when compared with the definite integral conversion method (Table 3 and Table 4) of the present embodiment, it was confirmed that the convergence to the solution was inferior.

Finally, the case where a quasi-random number which is a kind of the quasi-MC method is applied to the above-described definite integral form of (Equation 3) was examined. FIG. 6 is a diagram showing the calculation result (table 6). Since the quasi-random number is not a pseudo-random number, there may be a truncation error in formula application, but there is no statistical variation associated with the estimation of the required amount, so ρ is not described in Table 6. From the calculation result of Table 6, it can be seen that P ^ _f = 7.449 × 10 ⁻⁴ (the number of significant digits is 4 digits) has been obtained after one million trials. As described above, the definite integral conversion method of the present embodiment can use not only a uniform random number but also a quasi-random number, and it can be predicted that the present invention is not necessarily limited to a “random number”. is there.

  The simulation method according to the present invention is not limited to the reliability evaluation of structures, but can be applied to a wide variety of other simulations, computer graphics, and financial engineering, and its practical effects are extremely large.

  The present invention can be applied to, for example, a simulation by a digital computer.

DESCRIPTION OF SYMBOLS 100 Simulation apparatus 110 Definite integral conversion part 120 Probability estimation part 130 Estimation result output part 140 Function definition part 150 Random number generation part

Claims (6)

Based on a vector x representing n parameters, a state index function I (x) that takes a value of 0 or 1 depending on whether or not the event of interest occurs, and the event x of interest based on the vector x. A simulation apparatus for obtaining a probability Pf of occurrence of the event according to a probability density function fX (x) indicating the likelihood of occurrence,
When the state index function I (x) and the probability density function fX (x) are supplied, a definite integral conversion unit that forms the following (Equation 1) is provided.
A simulation apparatus having a function of calculating a probability Pf based on (Formula 1) formed by the definite integral conversion unit.

The definite integral converter is
It has a function of forming the following (Equation 2) obtained by discretizing (Equation 1),
The simulation apparatus according to claim 1, wherein the simulation apparatus has a function of calculating an estimated value of the probability P _f based on (Equation 2).

The definite integral converter is
The upper and lower limit values b _j and a _j of the integration of (Equation 1) and / or (Equation 2) are defined from the average value and standard deviation of each basic random variable X _j. 2. The simulation apparatus according to 2. 3. The simulation apparatus according to claim 2, wherein ξ shown in (Expression 2) is a uniform random number.
Based on a vector x representing n parameters, a state index function I ( x ) that takes a value of 0 or 1 depending on whether or not the event of interest occurs, and the event x of interest based on the vector x . According to a probability density function f _X ( x ) indicating the likelihood of occurrence, a simulation method by a simulation device for determining the probability P _f of occurrence of the event,
When the status indicator function I (x) and the probability density function f _{X (x)} is supplied, the constant integral transform unit forms the following equation (3),
A simulation method characterized in that the probability P _f is calculated based on (Equation 3) formed by the definite integral conversion unit.

Based on a vector x representing n parameters, a state index function I ( x ) that takes a value of 0 or 1 depending on whether or not the event of interest occurs, and the event x of interest based on the vector x . A computer program for determining a probability P _f of occurrence of the event according to a probability density function f _X ( x ) indicating the likelihood of occurrence,
When the status indicator function I (x) and the probability density function _f X (x) is supplied to form the following equation (4),
A computer program characterized in that the probability P _f is calculated based on the formed (Equation 4).

Images