JP6874239B2 - Rare event analysis device, rare event analysis method and rare event analysis program recording medium - Google Patents

Description

The present invention relates to a rare event analysis device, a rare event analysis method, and a rare event analysis program recording medium.

When designing an optical system such as a telescope, it is necessary to design so that stray light does not occur. Here, "stray light" is light that is potentially generated in optical design and is caused by a rare non-interfering light ray that does not follow the design optical path. That is, light rays that enter the light receiving surface through a non-regular path (optical path) are called stray light, which causes a decrease in image contrast and noise.

However, the conditions under which non-interfering rays occur in complex and large systems are very complex, and stray light is rare, even if the expert has extensive experience and good intuition. Experts sometimes overlook them. This fact not only reduces the performance of the designed system, but also reduces its credibility and reliability.

On the other hand, it is practically difficult to automatically search for all very rare conditions by applying all possible conditions by optical simulation.

Simulation methods for analyzing stray light are known. For example, Patent Document 1 is an "optical simulation device" based on a more accurate ray tracing method, in which one ray is divided into a transmitted ray and a reflected ray so that stray light can be handled by tracing both rays. Has been realized.

Japanese Unexamined Patent Publication No. 06-23769

However, Patent Document 1 does not suggest or disclose any device for determining the combination of parameters (events) to be input to the optical simulation device in order to efficiently discover the rare event of stray light.

An object of the present invention is to provide a rare event analysis device, a rare event analysis method, and a rare event analysis program recording medium that solves the above-mentioned problems.

The rare event analysis device according to the present invention is a rare event analysis device that discovers the possibility that a rare event may occur, and is N in each test k (1 ≦ k ≦ K) based on the proposed distribution Q (X). _{An evaluation function that determines the N parameters X k} , each of which defines an event, and sequentially sets the N parameters X _k one by one to define the possibility that the rare event will occur. A parameter determination unit supplied to a simulator having C (X); an evaluation value acquisition unit that acquires _{N evaluation values C (X k} ) output from the simulator corresponding to the _{N parameters X k. And;, the parameter determination unit includes N past parameters X k-1} used in N simulations of the past test (k-1) and the N past parameters X _k-1. Based on the N past evaluation values C (X _k-1 ) obtained in the N simulations for, the N parameters X _k to be used in the N simulations in this test k were determined. The parameter determination unit _{inputs the N past parameters X k-1} and the N past evaluation values C (X _k-1 ), sets the evaluation value as a class, and sets the evaluation value as a class, and sets the parameters of each class. A histogram calculation unit that calculates a histogram with the number as a frequency; based on the calculated histogram, the distribution of random numbers for generating _{N parameters X k used for N simulations in the test k this time.} , frequency parameters of the threshold value C _th or more of each class defined with respect to classes of the histogram, a random number to provide a histogram as larger than the frequency of the parameters in the threshold C _th following each class is obtained As such, it has a defined random number distribution definition unit; and a parameter generation unit that generates _{each parameter X k} to be used for each simulation in the present test k based on the defined random number distribution.

The rare event analysis method of the present invention is a rare event analysis method for discovering the possibility that a rare event may occur, and is N in each test k (1 ≦ k ≦ K) based on the proposed distribution Q (X). _{An evaluation function that determines the N parameters X k} , each of which defines an event, and sequentially sets the N parameters X _k one by one to define the possibility that the rare event will occur. Determining the parameters to be supplied to the simulator having C (X); and acquiring the _{N evaluation values C (X k} ) output from the simulator corresponding to the _{N parameters X k.} When; comprises, determining the parameters, the N number of past parameters X _k-1 used in the N times of the simulation of past test (k-1), the N pieces of past parameters X _{k Based} on the N past evaluation values C (X _{k-1 ) obtained in the N simulations for -1, the N parameters X k} to be used in the N simulations in this test k are determined. To determine the parameters, the N past parameters X _k-1 and the N past evaluation values C (X _k-1 ) are input, and the evaluation values are classified into classes. Then, to calculate a histogram in which the number of parameters of each class is a frequency; to generate _{N parameters X k to be used for N times of simulation in the present test k based on the calculated histogram.} The distribution of random numbers in the above is a histogram in which the frequency of the parameter in each class above the _{threshold C th} determined for the class of the histogram is greater than the frequency of the parameter in each class below the _{threshold C th.} To obtain a random number that gives the above; to generate _{each parameter X k} to be used for each simulation in the present test k based on the distribution of the defined random number.

The rare event analysis program recording medium of the present invention is a recording medium on which a rare event analysis program for discovering the possibility of a rare event occurring is recorded on a computer, and the rare event analysis program is proposed on the computer. _{Based on Q (X), N parameters X k} , each of which defines an event, to be used for N simulations in each test k (1 ≦ k ≦ K) are determined, and the N parameters are determined. A parameter determination procedure for sequentially supplying X _k one by one to a simulator having an evaluation function C (X) that defines the possibility of the rare event occurring; the simulator corresponding to the _{N parameters X k.} The evaluation value acquisition procedure for acquiring N evaluation values C (X _k ) output from is executed, and the parameter determination procedure causes the computer to perform N simulations of the past test (k-1). _{The N past parameters X k-1} used in the above and the N _{past evaluation values C (X k-1} ) obtained in N simulations for the N past parameters X _{k-1. Based on the above, N parameters X k} to be used for N simulations in this test k are determined, and the parameter determination procedure causes the computer to determine the N past parameters X _k-1. And the above N past evaluation values C (X _k-1 ), and a histogram calculation procedure for calculating a histogram in which the evaluation value is a class and the number of parameters of each class is a frequency; Based on the histogram, the distribution of random numbers for generating the _{N parameters X k} used for the N simulations in the current test k is determined for each of the _{thresholds C th or more determined for the class of the histogram.} A procedure for defining a random distribution distribution and defining a random number so as to obtain a random number that gives a histogram in which the frequency of the parameter in the class is higher than the frequency of the parameter in each class below the threshold C _th; Based on the distribution, the parameter generation procedure for generating _{each parameter X k} used for each simulation in the test k of this time and; are executed.

According to the present invention, extremely rare and fatal events in the target system can be automatically and reliably analyzed.

It is a figure which shows the schematic structure of the optical system designed by the telescope to which this invention is applied. It is explanatory drawing for demonstrating the Multicanonical MCMC method used in this invention. It is explanatory drawing which shows four variables which define the event X on the light receiving surface of the design optical system shown in FIG. It is a figure which shows the result of the search of a stray light in a four-dimensional parameter space. It is a figure which shows an example of the histogram which shows the sampled event on the fatality degree (evaluation function) C (X) in the Multicanonical MCMC method. It is a figure for demonstrating the Multicanonical MCMC method specialized for the rare event search. It is a figure which shows the comparison of the number of rays obtained which cause a stray light with respect to the number of histogram updates between the standard Multicanonical MCMC method (sm-MCMC) and the proposed Focused Multicanonical MCMC method (fm-MCMC). .. It is a block diagram which shows the optical system design support system by 1st Example of this invention. It is a block diagram which shows the structure of the parameter determination part used in the optical system design support system shown in FIG. 8 in detail. It is a block diagram which shows the optical system design support system by 2nd Example of this invention. It is a block diagram which shows the structure of the parameter determination part used in the optical system design support system shown in FIG. 10 in detail. It is a block which shows the rare event analysis apparatus according to 3rd Example of this invention. It is a block diagram which shows the structure of the parameter determination part used in the rare event analysis apparatus shown in FIG. 12 in detail.

First, in order to facilitate the understanding of the present invention, the meanings of the terms used in the present specification will be briefly described.

Markov chain Monte Carlo methods (MCMC method: Markov chain Monte Carlo methods) is a general term for algorithms that sample probability distributions based on creating Markov chains that have the desired probability distribution as a balanced distribution. In other words, the MCMC method is, as the name implies, a method that combines the "Markov chain" and the "Monte Carlo method". Specifically, the MCMC method is a method of performing various calculations by sampling from a probability distribution using the property of a stochastic process called a Markov chain.

The multicanonical algorithm is one of the sampling methods in the MCMC method, and is used when integrating an integrand having an arbitrary shape by the metropolis hasting method. In the multicanonical method, sampling is performed according to the reciprocal of the density of states. In the multicanonical method, the density of states is calculated in advance by another method such as the Wang and Landau method.

The Metropolis-Hastings algorithm is used to build Markov chains used to generate sequences of statistical samples from probability distributions that are difficult to sample directly. It is a method. This array is used as an approximation (histogram) of the target distribution in the MCMC method, or is used for those that require integral calculation such as the expected value.

The Wang and Landau algorithm is one of the Monte Carlo methods used to calculate the density of states of a system. In the Wang and Landau method, a non-Markov chain random walk is performed in order to quickly calculate the energy of all possible states of the system, which is necessary for calculating the density of states. The Wang and Landau algorithm is important for calculating the density of states required to carry out the multicanocal method.

One of the gist of the present invention is that the Multicanonical MCMC method is applied to the purpose of discovering stray light (rare event) in the design stage of an optical system such as a telescope.

However, it is not easy to combine the two for the following reasons.

This is because the Multicanonical MCMC method is used in some scientific and mathematical fields to generate random numbers that follow a long-tailed distribution, but is used in engineering for rare event analysis that generally does not have a long tail. Because it has not been done.

In addition, the Multicanonical MCMC method has been applied in various fields in statistical mathematics, physics, and chemistry. For example, the Multicanonical MCMC method was used to sample very rare examples of random matrices, random graphs, and chaotic dynamics. In a more applied aspect, the Multicanonical MCMC method was used to sample rare noise to theoretically evaluate the ideal performance of error-correcting codes in telecommunications. However, the applicability and scalability of the Multicanonical MCMC method for practical real-world problems has not been demonstrated.

In general, in the simulation method for stray light analysis, there was an optical path (ray) that started from the light receiving surface of the telescope designed by backpropagation and ended at a place other than the regular light entrance. In the case (that is, when there is an optical path other than the correct optical path), it is judged that the optical path (light ray) may be stray light.

The light receiving surface is represented by points (x, y) having two-dimensional coordinates in the horizontal and vertical directions, and the optical path (ray) is represented by angles α and β in the latitude and longitude directions with respect to each point. That is, one optical path (ray) is represented by a set of parameter values composed of four variables (x, y, α, β). In other words, each optical path (ray) is defined by four variables (a set of parameter values) (defined by a vector in four-dimensional space whose coordinates are the four parameters). Therefore, in order to search for stray light, it is necessary to perform a simulation for a huge number of optical paths (light rays) while sequentially changing the above four variables (a set of parameter values).

[Explanation of Outline of Invention]
We apply this statistical sampling technique to the search for extremely rare and deadly events for reliable and reliable engineering design. The proposal of the present invention is summarized as follows.

The present inventors propose a novel method for automatically and surely analyzing extremely rare and fatal events in a target system using the Multicanonical MCMC method.

We also propose an extension of the Multicanonical MCMC method, in particular, for the purpose of achieving a more efficient search for rare and fatal events.

The present inventors have found that in experiments using actual scale design of the telescope, check the prowess of the method of the present inventors, the inventors of the analysis, in a 10 ³ simulation tests of the order, I was able to discover all the stray lights with success. This is five orders of magnitude less than the number of tests required by simulations based on known grid searches. Here, the "grid search" refers to a method of searching for a parameter that causes an event (stray light in the example of this stray light analysis) to be searched by brute force searching for a combination of parameters.

[Analysis for search for stray light]
The purpose of stray light analysis is to determine if unwanted light rays have propagated to the light receiving surface of the image sensor and if their intensity is large enough to cause noise in the image. For example, the telescope optical system 10 shown in FIG. 1 has four mirrors.

In FIG. 1, the telescope optical system 10 includes a first mirror 11, a second mirror 12, a third mirror 13, a plane mirror 14, and a light receiving surface 15.

In FIG. 1, the solid line 16 represents an optical path of normal light designed to be focused on the light receiving surface 15 to form an image. On the other hand, the broken line 17 represents the optical path of the stray light that has reached the light receiving surface 15 by a non-regular optical path.

The purpose of stray light analysis is to search for an optical path that occurs in a designed telescope and reaches the light receiving surface with a non-negligible amount of light that traces other than the optical path of normal light.

Tools for stray light analysis have been developed primarily by using one of the two Monte Carlo principles. A simple technique is to use the forward propagation (forward propagation) principle, in which light rays are emitted from a random position at a random angle in a given light source, to evaluate the arrival of the telescope at the light receiving surface. This principle requires a wide range of positions and angles of emission to be simulated. Because the light source potentially surrounds the telescope. Therefore, analyzes that use forward propagation (forward propagation) generally tend to be inefficient. Due to this difficulty, most tools for analyzing stray light use the backward propagation (backpropagation) principle, in which light rays are virtually emitted at random angles from random locations on the telescope's light receiving surface. And check their presence outside the telescope. All optical elements of a telescope are symmetric with respect to the ingress and egress of light rays, except for optical converters such as spectrometers, which have asymmetric light input and output characteristics. Therefore, backward propagation (backpropagation) provides the same optical path as forward propagation (forward propagation). The latter principle is more efficient than the former. This is because the range of angles and positions at which light is emitted from the light receiving surface is limited in most cases.

When the probability of stray light occurrence is very low, the efficiency of tools that use the backward propagation (backpropagation) principle is insufficient to enable stray light brute force search. Therefore, these tools are used by experts to simulate some optical paths that heuristically reach the light receiving surface towards the generation of stray light. We are trying to increase the reliability and efficiency of stray light analysis by proposing a synthesis method that combines expert knowledge with this tool. However, the analysis performed by a human expert depends on the skill and knowledge of the expert and, in some cases, cannot search for all stray lights. This is because the conditions under which stray light occurs are extremely complex and rare. In addition, human analysis cannot handle the limited usefulness of skilled professionals.

[Rare and deadly event search for engineering design]
The following characteristics are required to discover the rare events that cause fatal consequences given in many engineering design tasks.

(1) Many types of rare and fatal events are efficiently searched for in real-world systems characterized by discontinuity, non-linearity and complexity.
(2) Events that are implicitly given and follow very low occurrences are efficiently searched.
(3) The probability of event occurrence is quantitatively evaluated.

The third feature cannot be ignored in optical design tasks. This is because a quantitative assessment of the probability of a fatal event occurring is often required in engineering design tasks. If the probability of an event occurring is very low, it can be ignored in optical design tasks.

With respect to (1), methods such as convex optimization, which assume certain features of the objective function in the search space, are not applicable to this problem. This is because the type of feature is unknown. In this regard, evolutionary optimization frameworks, stochastic optimization frameworks, and information science optimization frameworks are candidates for exploring stray light conditions. Because they can search for a large number of local maxima. For example, genetic algorithms and simulated annealing are typical algorithms and are often performed to find a large number of local maxima in various problems. This is because these algorithms ideally do not require any assumptions about the underlying conforming features. In recent years, Bayesian optimization has often been used in various engineering applications. This assumes that the unknown and underlying fitted features are well approximated by surrogate models generated by some principles such as the Gaussian process, and points to the highest uncertainty of feature estimation. By evaluation, the local maximum of the feature is searched. This method is effective and efficient when the characteristic function of interest is expensive to evaluate, non-differentiable, non-convex, and multimodal.

However, because the exploratory principles of these frameworks are opportunistic and probabilistic, they need to be found by contingencies due to their very low unknown probabilities, high isolation and high localization. It is easy to overlook the converted maximal. In addition, optimization techniques generally do not provide information about the probability of occurrence of rare and fatal events in the system. This is because they do not evaluate the probability distribution of event occurrence in their search mechanism. Therefore, these methods do not satisfy requirements (2) and (3).

In contrast, the stochastic sampling framework described below also requires no assumptions about the underlying conforming features, but by weighting the search space according to its critical degree, it is of interest and unknown rarity. Efficiently accomplish your search by focusing on fatal events. In addition, it efficiently obtains fair probability information for rare and fatal events by preserving their probability distributions throughout the search. Therefore, this framework is sufficient to meet all three requirements.

[Efficient stochastic sampling of rare and fatal events]
Histogram method, Multicanonical MCMC method, and replica exchange MCMC method are typical methods of stochastic sampling for efficiently generating events including rare and fatal events. These methods have been developed primarily in the field of statistical physics, where events characterized by various features with their weights that accurately preserve a given probability distribution need to be analyzed. It was.

The histogram method simply samples events from events that have a distribution independent of their probabilities according to the reciprocal of their probability density function. However, the method requires a probability density function explicitly known as prior knowledge and is supported in space with only a few dimensions for search efficiency.

The Multicanonical MCMC method uses a single Markov chain Monte Carlo (MCMC) method, and the single MCMC method sequentially seeks out events from an implicitly given probability distribution in the highest possible dimensional space. During the iterative sampling process, the Multicanonical MCMC method iteratively updates the weighting functions of the sampled events to ensure that they are evenly distributed with respect to the degree of lethality. Therefore, the Multicanonical MCMC method efficiently enumerates various events, including targeted rare and fatal events, with those weights that maintain their probability distributions accurately, without any bias. become.

The replica exchange MCMC method uses the MCMC method in parallel in layers, but is also known as parallel tempering or Metropolis-coupled MCMC. Each MCMC method follows an implicitly given probability distribution, the top MCMC method follows a target probability distribution, and the deep MCMC method follows a flatter distribution. This technique sequentially and randomly exchanges events generated in one layer with other events in adjacent layers while satisfying detailed balance of it. Since the upper layer introduces events that follow a flatter distribution from the lower layer without disturbing its distribution, it acquires various events, including rare and fatal events in the highest MCMC method that follows the probability distribution of interest. .. Although this exchange mechanism can generate rare and fatal events more efficiently than the Multicanonical MCMC method, each step of the individual MCMC method is used to check the degree of fatality of the generated events. Due to the parallel MCMC method including simulation, it takes a lot of calculation time.

Throughout these techniques, these techniques have been applied in various fields in statistical mathematics, physics, and chemistry. For example, they were used to sample very rare examples of random matrices, random graphs, and chaotic dynamics. In a more applied aspect, these techniques have been used to sample rare noise to theoretically evaluate the ideal performance of error-correcting codes in telecommunications. However, the applicability and scalability of these techniques for practical real-world problems has not been demonstrated.

The very rare and deadly event search method proposed by the present inventors for industrial design requiring high reliability uses the Multicanonical MCMC method, and this study is particularly complicated and unsuitable. It is applied to the stray light analysis of a telescope with continuous features. This method is suitable for problems in which the probability distribution of optical events generated by the Monte Carlo simulator is not clearly known and the generation of events by the simulator is expensive.

[Multicanonical MCMC method]
We first describe the use of importance sampling to locate events by applying weights based on the proposed distribution. Subsequently, MCMC sampling with multicanonical weights for efficiently sampling events, including rare and fatal events, will be described. Further, an algorithm for obtaining an appropriate proposed distribution in the Multicanonical MCMC method will be described.

およびＣ_ｔｈ∈［Ｃ_ｍｉｎ，Ｃ_ｍａｘ］は、事象の致命的な度合いおよび致命的な事象を規定するために与えられたその閾値レベル（以下、単に「閾値」とも呼ぶ。）である。ここで、本発明者らは、［Ｃ_ｍｉｎ，Ｃ_ｍａｘ］上での値域Ｃ_ｉ、ｉ＝１、・・・、Ｂから成るヒストグラムを使用する。ここで、値域Ｃ_ｉの高さは、次式で表される。Given the discrete event space S, X ∈ S is an event. Note that the event space S can be discrete or continuous. When the event space S is continuous, the sum Σ on the subset of the event space S is subsequently rewritten by ∫. Further, the probability distribution at which event X occurs is given by P (X). The probability distribution P (X) is explicitly given in many engineering design problems whose purpose is an artificial system, while the probability distribution P (X) is a record observed in the past in many problems involving natural objects. Implicitly given by. We define a set of fatal events R = { _{X ∈ S | C th} ≤ C (X)}. here,

And C _th ∈ [C _min , C _max ] is the lethal degree of the event and its threshold level given to define the fatal event (hereinafter, also simply referred to as “threshold”). Here, the present inventors use a histogram consisting of the range C _i , i = 1, ..., B _{on [C min} , C _max]. Here, the height of the range C _i is expressed by the following equation.

与えられた事象Ｘに対して、Ｐ（Ｃ（Ｘ））はＰ（Ｃ_ｉ）を表す。ここで、致命的な度合いＣ（Ｘ）∈Ｃ_ｉである。致命的な度合いは評価関数やエネルギ関数とも呼ばれる。また、値域Ｃ_ｉは階級とも呼ばれ、値域Ｃ_ｉの高さ（事象の個数）は度数とも呼ばれる。

For a given event X, P (C (X)) represents P (C _i ). Here, a fatal degree C (X) ∈C _i. The degree of fatality is also called the evaluation function or energy function. The range C _i is also called the class, and _{the height of the range C i} (the number of events) is also called the frequency.

By assuming that the event X is a fatal event having a probability distribution P (X) that occurs with a very low probability and searching directly from the probability distribution P (X), it is possible to obtain almost a rare event X. Can not. Therefore, the purpose is to efficiently derive rare and fatal events that satisfy the _{threshold C th} ≤ C (X) under a given probability distribution P (X), and to efficiently derive the overall rare and fatal events. Probability of R: To obtain a good estimate of P (R) = P ( _Cth ≤ C (C)).

[Importance sampling]
In order to efficiently generate a fatal event for a given probability distribution P (X), we artificially raise the event in R by a multicanonical weight G (C (X)). Importance sampling, sampled from the proposed distribution Q (X) ∝ G (C (X)) P (X), is applied. Importance sampling refers to a method of obtaining an expected value using a weighted sample. Here, the present inventors assume that the proposed distribution Q (X) ≠ 0 whenever the probability distribution P (X) ≠ 0.

[MCMC sampling with multicanonical weights]
Choosing an appropriate proposed distribution Q (X) is key to the efficient occurrence of events, including rare and fatal events, in importance sampling. The Multicanonical MCMC method selects the proposed distribution Q (X) approximately represented by the following equations 3 and 4.

ここで、ｃは定数であり、各Ｃ_ｉ、ｉ＝１、・・・、Ｂに対応するＧ（Ｃ（Ｘ））はマルチカノニカル重みとして知られている。

Here, c is a constant, each _{C i, i = 1, ···} , corresponding to B G (C (X)) is known as a multi-canonical weights.

By substituting the number 4 into the number 3, the number 5 is obtained on {X ∈ S | C (X) ∈ C _{i } for all ranges C i, i = 1, ···, B.} ..

ここで、Ｚ＝Σ_ＳＰ（Ｘ）／Ｐ（Ｃ（Ｘ））である。従って、Multicanonical ＭＣＭＣ法は、図２に示すように、領域［Ｃ_ｍｉｎ、Ｃ_ｍａｘ］上の一様な分布でサンプルすることによってＲ＝{Ｘ∈Ｓ｜Ｃ_ｔｈ≦Ｃ（Ｘ）}における希少で致命的な事象のサンプルを多く得る。事象Ｘが一般に高次元のベクトルであるので、ＭＣＭＣ法は、提案分布Ｑ（Ｘ）から事象Ｘを効率的にサンプルするために使用される。

Here, a _{Z = Σ S P (X)} / P (C (X)). Therefore, the Multicanonical MCMC method is _{rare in R = {X ∈ S | C th} ≤ C (X)} by sampling with a uniform distribution over the _{region [C min} , C _{max], as shown in FIG.} Get many samples of fatal events. Since event X is generally a high dimensional vector, the MCMC method is used to efficiently sample event X from the proposed distribution Q (X).

[Iterative estimation of proposed distribution]
As mentioned above, the probability distribution P (X) is given, but the probability distribution P (X) is not given directly in many practical problems, so the number 4 is the multicanonical weight G (C). (X)) cannot be applied directly to obtain. Therefore, we apply the Wang and Landau method to efficiently estimate the multicanonical weight G (C (X)) using iterative modification at every MCMC step.

The key idea is _{to gradually modify G (C (X)) in each range C i} , i = 1, ..., B, and the corresponding proposed distribution Q (C (X)) is the interval [C. It is to guarantee that it is uniform on [ _min , C _max]. The simplest algorithm, all multi-canonical weights G (C _i) is, after a sufficiently large number of sampled event is accumulated in range, is updated, an entropic sampling method. In contrast, the Wang and Landau method updates the range directly after the sample event is generated in the range by the MCMC method. This frequent update more efficiently estimates the desired multicanonical weight G (C (X)).

The Wang and Landau method works as follows. After each step in the MCMC process, if the generated event belongs to range C _i, the corresponding multi-canonical weights G (C (X)) is multiplied by a constant factor 0 <F <1. This update is repeated until the histogram is almost flat. At this point, all histogram ranges are reset to 0 and F is corrected to √F. Then, the update of the multi-canonical weights G (C _i), until the histogram is sufficiently flat, again repeated. Finally, we obtain the multicanonical weights G (C _i ), i = 1, ..., B, which guarantees a sufficiently uniform proposed distribution Q (C (X)).

[Probability of rare and fatal events]
_{Given the samples X i} , i = 1, ..., M found from the proposed distribution Q (X), the probability of the set R of fatal events is given as the next number 6.

ここで、Ｉ（Ｃ_ｔｈ≦Ｃ（Ｘ_ｉ））は指示関数であって、次の数７として規定される。

Here, I ( _Cth ≤ C (X _i )) is an indicator function and is defined as the following equation 7.

この数式は計算するのに不便である。何故なら、提案分布Ｑ（Ｘ）は、実際問題として、Ｇ（Ｃ（Ｘ））Ｐ（Ｘ）から容易に取得されないからである。従って、本発明者らは、それを数３および数４で置換することによってうまく処理し、マルチカノニカル重みＧ（Ｃ（Ｘ））のみから成る次の数８の数式を導出して、容易に計算する。

This formula is inconvenient to calculate. This is because the proposed distribution Q (X) is not easily obtained from G (C (X)) P (X) as a practical matter. Therefore, we have successfully processed it by substituting it with equations 3 and 4, and easily derived the following equation of equation 8 consisting only of the multicanonical weight G (C (X)). calculate.

ここで、次の数９の関連が使用される。

Here, the following association of equation 9 is used.

[Methodology for searching for rare stray lights]
We propose a new methodology for the search for rare stray light based on simulations for the engineering design of telescopes. This methodology consists of the Multicanonical MCMC method described above and the target telescope simulator.

We define event X by a vector that provides the operating conditions of the target system. It is assumed that event X follows a given probability distribution P (X). MCMC methods are used to generate the multi-canonical weights G (C _i) event X following the same proposal distribution Q (X) to an initial probability distribution P (X) under. Here, i = 1, ..., B, and all histogram classes are uniform (all bin widths are uniform). Then, the simulation of the objective system Sim (X) was derived using the event X generated as the input operating condition, and the critical degree in the situation of the simulation result from the input event X was given heuristically. It is evaluated by the measure (evaluation function) C (X). The fatal degree updates the count in the range of the histogram Q (C _i) corresponding the value of the initial zero and the multi-canonical weights G (C _i) of the range by intervening said one-Landau method Used to update. These updates provide updates to the proposed distribution as Q (X) ∝G (C (X)) P (X), and generate and fatally the next event X according to the new proposed distribution Q (X) by the MCMC method. The simulation of the objective system that provides the degree (evaluation function) C (X) is repeated. This process continues until the histograms Q (C _i ), i = 1, ..., B are almost flat, and the multicanonical weights G (C _i ), i = 1, ..., B are satisfactory. Converges to, and a sufficient number of events X are accumulated.

Algorithm 1 in Table 1 below represents the pseudo code for this method. It contains a double'while'loop, with the inner loop in _{the proposed distribution Q (X) until the output histogram Q (C i} ), i = 1, ..., B are sufficiently uniform. On the other hand, the outer loop intends to continue the MCMC method, _{and repeats the modification of the multicanonical weights G (C i} ), i = 1, ..., B until it finally reaches the optimum. Intended. The present inventors apply the transition function T _δ (X) to the metropolis method, which is a standard MCMC algorithm, to generate the next candidate event X'from the current event X. After convergence of the multi-canonical weights G (C _i) has become optimum, the data set D consisting of M events is generated with a weight, the probability P of catastrophic event (R) is estimated from the data set D.

Hereinafter, embodiments of the present invention will be described in detail with reference to the drawings. The configurations described in each of the following embodiments are merely examples, and the technical scope of the present invention is not limited thereto.
[Outline of the invention]

The present invention relates to an optical system simulator Sim (X), which is commercially available and standard in the field of optical system design, including telescopes, as described below. In this simulator, event X indicates a condition in which light rays are emitted from the light receiving surface 15 (see FIG. 1) of the image sensor by a method of backward propagation (back propagation). Therefore, event X is also called a condition. It consists of four variables X = [x, y, α, β] ^T as shown in FIG. In FIG. 3, x and y are the horizontal and vertical positions of the emission on the light receiving surface 15 defined as _{x ∈ [−x range} , x _range ], y ∈ [−y _{range , y range], respectively.} .. In addition, α and β are the latitude and longitude angles of the emitted light rays in orthogonal space, respectively. The normal vector of the light receiving surface 15 is 0 degrees. The range of angles is α ∈ [-90, 90] and β ∈ [-90, 90], respectively. The event (four variables) X is also called a parameter (input parameter).

The probability distribution P (X) is given by the product of four independent Gaussian distributions as the next number 10.

ここで、ｄ＝４であり、Σ_Ｐは、次の数１１で示される。

Here, d = 4, and Σ _P is represented by the following equation 11.

この分布は、座標空間の原点で中心に置かれる。σ_ｘおよびσ_ｙは、それぞれ、ｘ_{ｒａｎｇｅ}およびｙ_{ｒａｎｇｅ}として設定され、およびσ_αおよびσ_βは、それぞれ、９０度である。

This distribution is centered at the origin of the coordinate space. σ _x and σ _y are set as _{x range} and y _range , respectively, _{and σ α} and σ _β are 90 degrees, respectively.

This distribution is very wide, it is slightly higher in the center, but almost uniform in the direction perpendicular to the light receiving surface. The choice of this probability distribution P (X) reflects the telescope's symmetric geometry and mirror configuration.

The critical degree (evaluation function) C (X) is heuristic but needs to be carefully specified. Because it defines a fatal event. In the analysis of the embodiments of the present invention, the critical degree is due to the fact that the light beam emitted from the light receiving surface 15 under condition X reaches the outside of the telescope via backward propagation (back propagation). Mainly stipulated. cr (X) is 1 if this happens and 0 otherwise. However, cr (X) is not definitive for a given condition X due to random reflections and refractions on the optical path. Therefore, the intrinsic lethal degree of condition X is assessed by obtaining the average of cr (X) for n tests such as the following equation 12.

ここで、cr_i（Ｘ）は第ｉの試験における結果である。この致命的な度合い（評価関数）Ｃ（Ｘ）は、致命的な状況への確率的な近似に関して単調に増加し、閾値レベルＣ_ｔｈ以上の致命的な度合い（評価関数）Ｃ（Ｘ）を持つ条件Ｘは、致命的な事象を表す。閾値レベルＣ_ｔｈは、０．６に設定され、ｎ個の試験における条件Ｘの大部分は、結果として、このシミュレーションにおいて致命的な状況となる。ｎの値は、致命的な度合い（評価関数）Ｃ（Ｘ）のコンピュータ計算コストと統計上の安定性との間の関係性を考慮して決定される必要がある。致命的な度合い（評価関数）Ｃ（Ｘ）は、２項式の分布に従い、その相対的標準誤差は、近似的に√（１−Ｃ（Ｘ））／ｎＣ（Ｘ）であり、それは、ｎと致命的な度合い（評価関数）Ｃ（Ｘ）とで単調に増加する。ｎは、本発明者らの実験において１００に設定される。何故なら、誤差は、致命的な度合い（評価関数）Ｃ（Ｘ）が１に近接する致命的な条件Ｘに対して十分に抑制されるからである。

Here, cr _i (X) is the result in the i-th test. This fatal degree (evaluation function) C (X) increases monotonically with respect to a stochastic approximation to a fatal situation, and the fatal degree (evaluation function) C (X) above the _{threshold level C th is increased.} The condition X to have represents a fatal event. The threshold level C _th is set to 0.6, and most of the conditions X in the n tests result in a fatal situation in this simulation. The value of n needs to be determined in consideration of the relationship between the computer computational cost of the critical degree (evaluation function) C (X) and the statistical stability. The critical degree (evaluation function) C (X) follows the distribution of the binomial equation, and its relative standard error is approximately √ (1-C (X)) / nC (X). It increases monotonically with n and the fatal degree (evaluation function) C (X). n is set to 100 in our experiments. This is because the error is sufficiently suppressed for the fatal condition X whose fatal degree (evaluation function) C (X) is close to 1.

We use the following transition function Tδ (X), where the output X'is randomly given by a multivariable independently distributed Gaussian centered under condition X, as in equation 13.

ここで、ｄ＝４であり、Σ_Ｔは、次の数１４で示される。

Here, d = 4, and Σ _T is represented by the following equation 14.

ｄ_ｘ、ｄ_ｙ、ｄ_α、ｄ_βは、それぞれ、光挙動の空間又は角度的分解能であって、それらは、可視光線の波長と望遠鏡のスケールとに基づいて、光学において知られている。ハイパーパラメータδは、条件Ｘから新しい条件Ｘ’の偏差を支配する。もしハイパーパラメータδが大き過ぎれば、新しい条件Ｘ’が条件Ｘからそれる範囲が、新しい条件Ｘ’が条件Ｘに近接する致命的な事象を見落としてしまう。もしハイパーパラメータδが小さ過ぎれば、新しい条件Ｘ’はほとんど条件Ｘの近傍に残って、条件Ｘから離れた致命的な事象をほとんど見つけられない。従って、試行錯誤接近方法を使用して、ハイパーパラメータδの値を適切に調整する必要がある。

d _x , _dy , d _α , d _β are the spatial or angular resolutions of light behavior, respectively, and they are known in optics based on the wavelength of visible light and the scale of the telescope. The hyperparameter δ governs the deviation of the new condition X'from the condition X. If the hyperparameter δ is too large, a fatal event in which the new condition X'deviates from the condition X and the new condition X'is close to the condition X is overlooked. If the hyperparameter δ is too small, the new condition X'remains almost in the vicinity of the condition X, and few fatal events away from the condition X can be found. Therefore, it is necessary to adjust the value of hyperparameter δ appropriately by using the trial and error approach method.

The maximum number of histogram updates, K _max, and the maximum number of MCMC steps, M _max , were set to 10 and 5000, respectively, throughout all of our experiments. The former is sufficient to approximate the F to unity, to fine-tune Multicanonical weight G (C _i). The latter also allows for sufficient relaxation of the MCMC process. The last parameter to be adjusted is the number B of the histogram range on the critical degree (evaluation function) C (X). Events below the maximum number M _max are uniformly dispersed over these B ranges, so the probability of all events occurring in one range is the same, and the number of events in one range is relative to each other. It follows a Poisson distribution with a typical standard deviation ε ≤ √ (B / M _max). We set B = 20, which maintains a reasonably small ε = 0.05 to 0.10.

[Embodiment]
Next, an embodiment for carrying out the present invention will be described using the embodiment.

[Experimental setup]
We used a lighting design analysis software called Light Tools released by Synopsys, Inc. This software is widely used in optical design and analysis. This software is capable of Monte Carlo ray tracing. The target optical system is a full-scale model of the designed telescope on a real scale. We have empirically known that the two optical paths cause stray light of similar intensity in the two symmetric regions of the light receiving surface 15 due to the geometric symmetry of the telescope.

The assumptions and conditions of this simulation are as follows.
We assume that the model used in the simulation is a perfect replica of the system of a life-size model in the real world.
-The mechanism for producing stray light is well approximated only by geometric optics.
• All tracked rays are well approximated by a single wavelength.
• The "backpropagation method" is used for efficient analysis in all experiments.
-Ray tracing is stopped under either of the following conditions:

- Configuration normalized energy xi] of is smaller than the threshold value xi] _th. This threshold ξ _th is considered to be the lower limit of the non-negligible relative light intensity. In practice, ξ _th is set to ^10-4 times the original ray energy.

-The rays reach the infinite space outside the telescope without encountering any obstacles.
-Reflection and refraction are calculated as follows.

-The mode of reflection and refraction on an object is selected according to the parameters given to that object and its surface.

-The angle of reflection and refraction at the object is stochastically selected according to the given parameters.

-A ray does not branch at any reflection point and any refraction point, and the incident light is attenuated to become reflected light.

The main programs including the above Algorithm 1 were implemented by the programming language Python. Each step of the MCMC method was performed by the main program in the Light Tools simulator via the Excel VBA and Light Tools user interfaces. All calculations in our experiments were performed using a general purpose personal computer with a 3.4 GHz Intel Core i7 core processor and 16 MB of memory, without additional hardware.

[Experimental result]
The present inventors performed a stray light analysis using the above [methodology] and the above [setup]. FIG. 4 shows the result of searching for stray light in a four-dimensional parameter space.

FIG. 4A shows a result histogram of the light rays 17 passing through the non-normal path and the arrival positions (x, y) on the light receiving surface 15. Higher ray counts are represented by brighter colors. This histogram shows two important peaks that represent a non-negligible non-normal ray, a set of stray lights.

4 (b-1) and (b-2) show histograms of the rays belonging to the selected region around the two peaks, respectively. These are represented in the space of the incident direction (α, β), and clearly indicate that there are two optical paths symmetrical with respect to β, which are the paths that produced one of the two peaks. The results obtained for these positions and paths and the intensity of the stray light 17 correspond well with the stray light obtained in the experimental results of a full-scale model telescope.

FIG. 5 (a) is a histogram showing sampled events on the critical degree (evaluation function) C (X) in the Multicanonical MCMC method of this analysis. In FIG. 5A, the horizontal axis represents the evaluation function C (X) and the vertical axis represents the sampled event. The total number of simulation tests and their calculation time to achieve this convergence and obtain the results shown in FIG. 4 were 5 × 10 ⁴ and 17.3 hours, respectively.

FIG. 5A shows a histogram in which the evaluation value is a class and the number of sampled events (parameters) of each class is a frequency. The histogram of FIG. 5A is normalized so that the maximum value of the class is 1.0.

[Multicanonical MCMC method specialized for rare event search]
As mentioned above, this technique follows standard Multicanonical MCMC methods to sample events and includes rare events with weights that accurately maintain a given probability distribution. This requires accurate analysis of the rare event probabilities P (R) and their associated statistics in statistical physics and other fields of application. The exact weights are obtained by generating various events over a region of critical degree [C _min , C _max ] to avoid statistical deviations, and therefore the objective histogram Q (C _i ), i = 1. , ..., B are set uniformly.

An alternative goal in practical problems, including our current applications in engineering design, is to discover many rare and deadly events, rather than derivatives of events with precise weights. is there. This can be achieved efficiently by treating the region of fatal events [C _th , C _max ] more heavily than the other regions [C _min , C _th ], as shown in FIG. The following analysis shows that the range in which a fatal event is observed is intensively and most efficiently reduced without exploring the event space S evenly, thereby focusing on the fatal event. deal with.

[Optimal proposal distribution]
A fully relaxed MCMC process meets detailed balance and the generated event X is in steady state. Thus, the optimal proposed distribution Q (C (X)) condition makes it possible to infer fatal events that should occur efficiently from the analysis of a single transition step in the MCMC process.

で表され、全てのｉ≠ｊ、ｊ、ｊ＝１、・・・、ｒに対して、Ｒ_ｉ∩Ｒ_ｊ＝０であり、／Ｒ＝Ｓ＼Ｒである。Ｒ_ｉは、致命的な事象の集合であり、／Ｒは、非致命的な事象の集合である。ここで、簡単のためのメタ事象として各事象集合を参照する。各致命的なメタ事象Ｒｉは、十分に特定の場所に局在化されており、他の致命的なメタ事象から十分に分離されているとする。Since the fatal event of the goal of the present inventors is positioned in a large number of regions belonging to the event space S as demonstrated by the analysis of the present inventors in the above [experimental result], the event space S. Can be modeled as S = R∪ / R. Here, R is the next number 15

It is represented by, and for all i ≠ j, j, j = 1, ..., R, R _i ∩ _{R j} = 0 and / R = S \ R. _Ri is a set of fatal events and / R is a set of non-fatal events. Here, each event set is referred to as a meta event for simplicity. It is assumed that each fatal meta-event Ri is well localized and well separated from other fatal meta-events.

It is assumed that Q (Y'| Y) is the transition probability of event X from meta-event Y to another meta-event Y'in the MCMC process that produces an event that follows the proposed distribution Q (X). From the above assumptions, it is further assumed that discrete transitions between any fatal metaevents are ignored. Therefore, only non-zero transition probabilities between _{all Ri and / R are considered.} The following conditions are retained by detailed balance between the state of the overall meta-event probability and _{Ri and / R.}

Considering the event space S in the most efficient way to avoid overlooking any rare and fatal event requires the highest transition probability between _{Ri and / R in total.}

By substituting the _{equation 16 and differentiating it by Q (R i} ), the following conditions are obtained for the maximum value.

Further, by substituting this equation with the equation 17, taking the sum of all i's, and applying the constraint of the equation 16, the following simple condition is obtained.

This result derives Q (R) = Q (/ R) = 0.5, that is, the following number 21.

Since the uniform distribution of Q (C (X)) provides the minimum deviation of any statistical inference, it is _{constant at Ri} and / R, respectively, while satisfying the above conditions, Q ( The following target histogram of C (X)) is used.

Are provided to define a fatal _{events, C th is, [C min,} C _max] Since generally close to the upper limit of the value of Q for catastrophic events (C (X)), non Greater than that for fatal events.

[Focused Multicanonical MCMC method]
Based on the results given in the previous [Optimal Proposed Distribution], we have named the "Focused Multicanonical MCMC method" to efficiently discover rare and fatal events. We propose an extension of the MCMC method. In Algorithm 1 of Table 1 above, the shape of the target histogram required for Q (Ci), i = 1, ..., B to converge is not flat, but a step function given by Equation 22. .. Other than this change, no other changes are required for the above algorithm.

The probability of a rare and fatal event P (R) is also easily calculated. The stepped histogram of Q (C (X)) in Equation 22 is provided by Equation 23 below.

ここで、ｃ（Ｃ（Ｘ））は、次の数２４で表される。

Here, c (C (X)) is represented by the following number 24.

そして、Ｚ_Ｃ＝Σ_Ｓｃ（Ｃ（Ｘ））Ｐ（Ｘ）／Ｐ（Ｃ（Ｘ））である。これは、数２３を数３に代入し、それをすべてのＣ_ｉ、ｉ＝１、・・・、Ｂに対して{Ｘ∈Ｓ｜Ｃ（Ｘ）∈Ｃ_ｉ}において周辺化することによって、容易に確認される。

Then, a _{Z C = Σ S c (C} (X)) P (X) / P (C (X)). This is done by substituting the number 23 into the number 3 and _{marginalizing it for all C i} , i = 1, ···, B in {X ∈ S | C (X) ∈ C _i }. , Easily confirmed.

By substituting the numbers 3, 23 and 24 into the number 6, the next number 26 is obtained.

ここで、次の数２７の関係が使用される。

Here, the relationship of the following equation 27 is used.

Therefore, once G (C (Xi)) is obtained by the Wang and Landau method, the probability P (R) of a rare and fatal event is calculated in a manner similar to the standard method.

[Experimental results of Focused Multicanonical MCMC method]
FIG. 5B shows a histogram of sampled events on C (X) in our Focused Multicanonical MCMC method using the _{threshold C th = 0.6 as described above.} This histogram has a stepped shape and is balanced by the _{threshold Cth.} To achieve this convergence, the total number of simulation experiments and their calculation time, respectively, was 6.1 × 10 ⁴ and 20.7 hours.

FIG. 7 shows a comparison of the number of rays obtained that cause stray light on histogram updates between the standard Multicanonical MCMC method (sm-MCMC) and the proposed Focused Multicanonical MCMC method (fm-MCMC). As can be easily confirmed from this figure, the efficiency of obtaining rare stray light by the proposed method converges to almost twice as high efficiency as that of the standard Multicanonical MCMC method.

[Comparison with other methods]
Table 2 shows a performance comparison of grid search, Bayesian optimization, standard Multicanonical MCMC method, and Focused Multicanonical MCMC method. Grid search is a comprehensive search algorithm for the entire event space and provides baseline performance. We have not evaluated the genetic algorithm and simulated annealing. Because they are infeasible by requiring a huge number of simulations.

The first and second columns of this Table 2 show the probability of finding a rare and fatal event for each simulation test, that is, the efficiency of the search. The first column shows the efficiency when the position [x, y] is fixed at the position where the stray light is observed and the direction [α, β] of the incident light causing the stray light is searched. The efficiency of Bayesian optimization is an order of magnitude higher than that of grid search, while the efficiency of standard Multicanonical MCMC and Focused Multicanonical MCMC methods is two orders of magnitude higher. The second column of Table 2 shows the efficiency when searching the entire four-dimensional space. Due to the large search space, Bayesian optimization, in dealing with 10 ⁵ number of tests, I could not find any fatal events. In contrast, the two Multicanonical MCMC methods can successfully discover one rare and fatal event within 1000 trials, an efficiency five orders of magnitude higher than grid search. In both 2D and 4D cases, the efficiency of the Focused Multicanonical MCMC method is twice that of the standard Multicanonical MCMC method. These results follow the discussion in [Multicanonical MCMC method specialized for rare event search].

The last column of Table 2 shows an estimate of the probability P (R) of a rare and fatal event, i.e. the risk of stray light. Grid search estimates are expected to contain high uncertainty. This is because grid search can only find a very limited number of rays that cause stray light, thus leading to large statistical errors. Bayesian optimization cannot provide this probability. Because it does not have stochastic information in the search. The standard Multicanonical MCMC and Focused Multicanonical MCMC methods numbers 8 and 26 provide estimates as large as the probabilities obtained by grid search.

FIG. 8 is a block diagram showing an optical system design support system 100 according to the first embodiment of the present invention using the optical system simulator Sim (X). The optical system design support system 100 is a system for discovering the possibility of stray light occurring in the optical system 10 designed like a telescope as shown in FIG.

The optical system design support system 100 includes a processing unit 200 and a storage unit 300. The processing unit 200 includes, for example, a processor such as a CPU (Central Processing Unit), but is not limited to this, of course. The storage unit 300 includes, for example, a RAM (random access memory), a ROM (read only memory), and the like, but of course, the storage unit 300 is not limited thereto.

The processing unit 200 includes a parameter determination unit 210, a locus calculation unit 220, and an evaluation value calculation unit 230. The storage unit 300 includes a parameter storage unit 310, an optical model storage unit 320, and an evaluation value storage unit 330.

The parameter determination unit 210 determines the parameters as described later. As each parameter, the parameter determination unit 210 has a two-dimensional position (x, y) on the light receiving surface 15 of the designed optical system 10 and an angle (α) in two directions in which the light beam enters the light receiving surface 15. , Β), a total of four-dimensional parameter combinations are determined. The parameter storage unit 310 stores the determined parameter.

The locus calculation unit 220 calculates the locus of light rays in the designed optical system 10 by using the determined parameters, as will be described later. The locus calculation unit 220 calculates the locus in which the light beam propagates back from the light receiving surface 15 toward the light inlet. The optical model storage unit 320 stores an optical model in which a plurality of optical components are arranged on the optical path from the light inlet to the light receiving surface 15 in the designed optical system 10. This optical model corresponds to the designed optical system 10. The locus calculation unit 220 receives the input of the optical model stored in the optical model storage unit 320, and calculates a locus indicating the behavior of the light beam in the optical model.

The evaluation value calculation unit 230 calculates an evaluation value indicating the possibility of stray light being generated for the parameter based on the comparison between the calculated trajectory and the trajectory expected under the parameter in the ideal optical system. To do. In other words, the evaluation value calculation unit 230 has an evaluation function C (X) uniquely defined according to the above-designed optical system 10 (optical model). That is, the optical system simulator Sim (X) has an evaluation function C (X). The evaluation value calculation unit 230 receives the calculated locus and the determined parameter, and calculates the evaluation value according to the evaluation function C (X). The evaluation value storage unit 330 stores the calculated evaluation value.

Therefore, the parameter storage unit 310 stores the parameters used for calculating the past locus. The evaluation value storage unit 330 stores the evaluation value calculated for the parameter.

The parameter determination unit 210 is used for the trajectory calculation from the next time onward so that the higher the evaluation value is, the closer to the parameter is based on the parameter stored in the parameter storage unit 310 and the evaluation value stored in the evaluation value storage unit 330. Determine the parameters.

The parameter determination unit 210 defines a distribution of random numbers so that a random number closer to the parameter can be obtained as the evaluation value is higher, and generates a parameter to be used for trajectory calculation from the next time onward based on this distribution.

FIG. 9 is a block diagram showing in detail the configuration of the parameter determination unit 210 used in FIG. The parameter determination unit 210 includes a histogram calculation unit 212, a random number distribution definition unit 214, and a parameter generation unit 216.

Hereinafter, the operation of the optical system design support system 100 will be described in more detail with reference to FIGS. 8 and 9.

The parameter storage unit 310 stores N parameters used for calculating the past trajectory.

Based on the above distribution, the parameter determination unit 210 determines N parameters to be used for the trajectory calculation N times from the next time onward.

The locus calculation unit 220 calculates N loci for each of the determined N parameters.

The evaluation value calculation unit 230 calculates N evaluation values for each of the determined N parameters and the calculated N loci. The evaluation value storage unit 330 stores N evaluation values calculated for N parameters.

The histogram calculation unit 212 receives N parameters used for calculating the past trajectory from the parameter storage unit 310, and receives N evaluation values calculated for the N parameters from the evaluation value storage unit 330. The histogram calculation unit 212 calculates a histogram in which the evaluation value is a class and the number of parameters of each class is a frequency.

Based on the calculated histogram, the random number distribution definition unit 214 determines the distribution of random numbers for generating N parameters to be used for N times of trajectory calculation from the next time onward, so that the histogram becomes substantially flat. Defined to be obtained (see figure on the left side of FIG. 6).

The parameter generation unit 216 generates each parameter to be used for each trajectory calculation from the next time onward according to the defined distribution of random numbers. That is, as shown in FIG. 5A, the parameter generation unit 216 generates N parameters so that the histogram is substantially flat.

Therefore, the parameter determination unit 210 determines N parameters to be used for the trajectory calculation of N times from the next time onward by using the standard Multicanonical MCMC method.

More specifically, when random numbers (parameters) are randomly generated, as shown in the figure on the left side of FIG. 2, a large number of parameters with evaluation values smaller than the _{threshold C th = 0.6 are obtained.} Will be. Therefore, the parameter determination unit 210 adopts all random numbers (parameters) that can obtain an evaluation value _{of the threshold C th} = 0.6 or more, but a random number (parameter) that can obtain an evaluation value smaller than the _{threshold C th = 0.6.} As shown in the figure on the right side of FIG. 2, N parameters are determined so that the histogram becomes substantially flat by discarding most of the parameters and adopting them occasionally.

By repeating the above operation a plurality of times, the optical system design support system 100 can determine N parameters so that the histogram is finally substantially flat as shown in FIG. 5 (a). it can. Here, "substantially flat" means that the frequency in each class is within a range of a predetermined error with respect to the average frequency. This will be described by taking the case of FIG. 5A as an example. It is assumed that the average frequency is 75 and the predetermined error is 15. In this case, the frequency in each class is in the range of 60 to 90.

As is clear from the above description, according to the first embodiment, the stray light in the optical system can be automatically and surely analyzed.

FIG. 10 is a block diagram showing an optical system design support system 100A according to a second embodiment of the present invention using the optical system simulator Sim (X). The optical system design support system 100A is also a system for discovering the possibility of stray light occurring in the optical system 10 designed like a telescope as shown in FIG.

The illustrated optical system design support system 100A has the same configuration and operation as the optical system design support system 100 shown in FIG. 8, except that the configuration and operation of the processing unit are different as described later. Therefore, a reference code of 200A is attached to the processing unit. Those having the same functions as those shown in FIG. 8 are designated by the same reference numerals, and the description thereof will be omitted for the sake of brevity.

The processing unit 200A has the same configuration and operation as the processing unit 200 shown in FIG. 8, except that the configuration and operation of the parameter determination unit are different as described later. Therefore, a reference code of 210A is attached to the parameter determination unit.

FIG. 11 is a block diagram showing in detail the configuration of the parameter determination unit 210A used in the optical system design support system 100A of FIG. The parameter determination unit 210A has the same configuration and operation as the parameter determination unit 210 shown in FIG. 9, except that the configuration and operation of the random number distribution definition unit are different as described later. Therefore, the reference code of 214A is attached to the random number distribution definition unit.

Random distribution defining unit 214A based on the calculated histogram, the distribution of the random numbers for generating N parameters to be used for path calculation of the next time N times, the threshold C _th defined with respect to the class of the histogram It is defined so that a random number that gives a histogram that gives a histogram in which the frequency of the parameter in each of the above classes is higher than the frequency of the parameter in each class below the _{threshold Cth can be obtained.}

In the illustrated example, the maximum value of the histogram class is normalized to be 1.0, and the threshold value C _th is equal to 0.6.

Random distribution defining unit 214A is a random number that the total of the frequencies of the sum and the threshold value C _th following parameters of the frequency threshold C _th or more parameters determined for the class of the histogram and gives histogram as substantially equal Define the distribution so that Here, "substantially equal", within a predetermined error range, it means that the threshold value C _th or more in total with a threshold value C _th following parameters of the frequency parameter frequency total and are equal.

The parameter generation unit 216 generates each parameter to be used for each trajectory calculation from the next time onward, as shown in FIG. 5B, according to the defined distribution of random numbers. That is, the parameter generating unit 216, as shown in FIG. 5 (b), the histogram is, the threshold value C _th or more in total with a threshold value C _th following parameters of the frequency parameter frequency total and of substantially equal As described above, N parameters are generated.

Therefore, the parameter determination unit 210A determines N parameters to be used for the trajectory calculation of N times from the next time onward by using the Focused Multicanonical MCMC method.

More specifically, when random numbers (parameters) are randomly generated, the evaluation values are those with a high parameter frequency at an evaluation value smaller than the _{threshold value C th = 0.6, as shown in the figure on the left side of FIG.} Will be obtained in large numbers. Therefore, the parameter determination unit 210A adopts all random numbers (parameters) that can obtain an evaluation value _{of the threshold value C th} = 0.6 or more, but a random number (parameter) that can obtain an evaluation value smaller than the _{threshold value C th = 0.6.} Most of the items are discarded and occasionally adopted. Thus, as shown in FIG. 5 (b), the parameter determination unit 210A, the histogram, the sum of the frequencies of the sum and the threshold value C _th following parameters of the frequency threshold C _th or more parameters are substantially Determine N parameters so that they are equal.

Optical design support system 100A, by repeating several times the operation, as shown in FIG. 5 (b), histogram of the frequency threshold C _th or more parameters sum threshold value C _th following the frequency parameters N parameters can be generated such that the sum is substantially equal.

In this example, the random number distribution definition unit 214A sets the threshold value C _th to 0.6, but the random number distribution definition unit 214A may set the threshold value C _th to a value larger than 0.5. However, the random number distribution defining unit 214A, it is preferable to set the threshold value _{C th} between the range of 0.6 to 0.8. However, it is more preferable that the random number distribution definition unit 214A sets the threshold value _Cth to a value of 0.6. The reason is as follows.

In general, the number of stray lights is not limited to one, and in most cases, there are a plurality of stray lights due to the above-mentioned symmetry and the like. In such a situation, if the threshold C _th is set higher than 0.8, most of the random numbers (parameters) of the evaluation values smaller than the _{threshold C th are discarded, and it is possible to search for other stray lights.} This is because the sex (probability) becomes very low. Further, if the threshold value C _th is set lower than 0.6, it becomes impossible to efficiently search for a plurality of stray lights. Therefore, it is preferable to set the _{threshold value Cth to a value of 0.6.}

As is clear from the above description, according to the second embodiment, a plurality of stray lights in the optical system can be automatically and efficiently analyzed.

FIG. 12 is a block diagram showing a rare event analysis device 100B according to a third embodiment of the present invention. The rare event analysis device 100B is a device that discovers the possibility that a rare event may occur. Specifically, the rare event analysis device 100B discovers, for example, stray light generated in an optical system 10 designed like a telescope as shown in FIG. 1 as a rare event. However, the rare event analysis device 100B can be applied not only to the simulation of the optical system 10 designed as such, but also to the simulation of other target systems.

The rare event analysis device 100B includes a processing unit 200B and a storage unit 300B. The processing unit 200B includes, for example, a processor such as a CPU (Central Processing Unit), but is not limited to this, of course. The storage unit 300B includes, for example, a RAM (random access memory), a ROM (read only memory), and the like, but of course, the storage unit 300B is not limited thereto.

The processing unit 200B includes a parameter determination unit 210B, a simulator 220B, and an evaluation value acquisition unit 230B. The storage unit 300 includes a parameter storage unit 310B and an evaluation value storage unit 330B.

The simulator 220B may be, for example, the above-mentioned optical system simulator Sim (X). In this case, the simulator 220B can be realized by the combination of the optical model storage unit 320, the locus calculation unit 220, and the evaluation value calculation unit 230 shown in FIG. However, the simulator 220B is not limited to the optical system simulator Sim (X), and may be another simulator. The simulator 220B has an evaluation function C (X) that defines the possibility that a rare event will occur.

_{The parameter determination unit 210B determines N parameters X k} to be used for N simulations in each test k (1 ≦ k ≦ K) based on the proposed distribution Q (X). Each of the N parameters X _k defines an event. The parameter determination unit 210B _{sequentially supplies N parameters X k} to the simulator 220B one by one. The parameter storage unit 310B stores the determined parameter.

The evaluation value acquisition unit 230B acquires _{N evaluation values C (X k} ) output from the simulator 220B corresponding to _{N parameters X k.} The evaluation value storage unit 330B stores these N evaluation values.

_{The parameter storage unit 310B stores N past parameters X k-1} used in N simulations of the past test (k-1). The evaluation value storage unit 330B stores N past evaluation values C (X _k-1 ) _{acquired in N simulations for those N past parameters X k-1.}

The parameter determination unit 210B has N past parameter X _{k-1 stored in the parameter storage unit 310B and N past evaluation values C (X k-1} ) stored in the evaluation value storage unit 330B. ) And, as will be described later, N parameters X _k to be used for N simulations in this test k are determined.

FIG. 13 is a block diagram showing in detail the configuration of the parameter determination unit 210B used in FIG. The parameter determination unit 210B includes a histogram calculation unit 212B, a random number distribution definition unit 214B, and a parameter generation unit 216B.

The histogram calculation unit 212B includes N past _{evaluation values X k-1 stored in the parameter storage unit 310B and N past evaluation values C (X k-1} ) stored in the evaluation value storage unit 330B. Receive. The histogram calculation unit 212B calculates a histogram in which the evaluation value is a class and the number of parameters of each class is a frequency.

Based on the calculated histogram, the random number distribution definition unit 214B determines the distribution of random numbers for generating the _{N parameters X k used in the N simulations of the test k this time with respect to the histogram class.} It is defined so that a random number can be obtained that gives a histogram in which the frequency of the parameter in each class above the threshold C _{th is higher than the frequency of the parameter in each class below the threshold C th.}

In the illustrated example, the maximum value of the histogram class is normalized to be 1.0, and the threshold value C _th is equal to 0.6.

Random distribution definition unit 214B generates a random number that the sum of the frequencies of the sum and the threshold value C _th following parameters of the frequency threshold C _th or more parameters determined for the class of the histogram and gives histogram as substantially equal Define the distribution so that Here, "substantially equal", within a predetermined error range, it means that the threshold value C _th or more in total with a threshold value C _th following parameters of the frequency parameter frequency total and are equal.

_{The parameter generation unit 216B generates each parameter X k} to be used for each simulation in the present test k, as shown in FIG. 5B, based on the defined random number distribution. That is, the parameter generating unit 216B, as shown in FIG. 5 (b), the histogram is, the threshold value C _th or more in total with a threshold value C _th following parameters of the frequency parameter frequency total and of substantially equal As described above, N parameters X _k used for N simulations in this test k are generated.

Therefore, the parameter determination unit 210B uses the Focused Multicanonical MCMC method to determine _{N parameters X k to be used for N simulations in the test k this time.}

More specifically, when random numbers (parameters) are randomly generated, the evaluation values are those with a high parameter frequency at an evaluation value smaller than the _{threshold value C th = 0.6, as shown in the figure on the left side of FIG.} Will be obtained in large numbers. Therefore, the parameter determination unit 210B adopts all random numbers (parameters) that can obtain an evaluation value _{of the threshold value C th} = 0.6 or more, but a random number (parameter) that can obtain an evaluation value smaller than the _{threshold value C th = 0.6.} Most of the items are discarded and occasionally adopted. Thus, as shown in FIG. 5 (b), the parameter determination unit 210B, the histogram, the sum of the frequencies of the sum and the threshold value C _th following parameters of the frequency threshold C _th or more parameters are substantially _{N parameters X k} used for N simulations in this test k are generated so as to be equal.

Rare event analysis apparatus 100B, by repeating K times the test operation, as shown in FIG. 5 (b), histogram of the frequency threshold C _th or more parameters sum threshold value C _th following the frequency parameters In the final test K, N final parameters X _K can be generated so that the sums are substantially equal.

In this example, the random number distribution definition unit 214B sets the threshold value C _th to 0.6, but the random number distribution definition unit 214B may set the threshold value C _th to a value larger than 0.5. However, the random number distribution definition section 214B, it is preferable to set the threshold value _{C th} between the range of 0.6 to 0.8. However, it is more preferable that the random number distribution definition unit 214B sets the threshold value _Cth to a value of 0.6. The reason is as follows.

In general, the number of rare events is not limited to one, and in most cases, there are a plurality of rare events. In such a situation, if the threshold C _th is set higher than 0.8, most of the random numbers (parameters) of the evaluation value (class) smaller than the _{threshold C th are discarded, and other rare events.} This is because the possibility (probability) of being able to search for is very low. Further, if the threshold value C _th is set lower than 0.6, it becomes impossible to efficiently search for a plurality of rare events. Therefore, it is preferable to set the _{threshold value Cth to a value of 0.6.}

As is clear from the above description, according to the third embodiment, it is possible to automatically and efficiently analyze a plurality of rare events in the target system.

The present invention is not limited to the above-described embodiment as it is, and at the implementation stage, the components can be modified and embodied within a range that does not deviate from the gist thereof. In addition, various inventions can be formed by appropriately combining a plurality of components.

Each part of the optical system design support system (rare event analysis device) may be realized by using a combination of hardware and software. In the form of combining hardware and software, an optical system design support (rare event analysis) program is deployed in RAM (random access memory), and a CPU (central) is developed based on the optical system design support (rare event analysis) program. By operating the hardware of the processor such as processing unit), each part is realized as various means. Further, the optical system design support (rare event analysis) program may be recorded on a recording medium and distributed. The optical system design support (rare event analysis) program recorded on the recording medium is read into the memory via a wired, wireless, or recording medium itself, and operates a processor (CPU) or the like. Examples of recording media include optical disks, magnetic disks, semiconductor memory devices, hard disks, and the like.

To explain the first and second embodiments in another expression, the parameter determination units 210, 210A and the locus of a computer operating as an optical system design support system are based on an optical system design support program developed in RAM. This can be achieved by operating as a combination of the calculation unit 220 and the evaluation value calculation unit 230.

To explain the third embodiment in another expression, the parameter determination unit 210B, the simulator 220B, and the evaluation value acquisition unit are based on the rare event analysis program developed in the RAM of the computer operating as the rare event analysis device. It can be realized by operating as a combination of 230B.

Further, the specific configuration of the present invention is not limited to the above-described embodiment, and is included in the present invention even if there is a change within a range not deviating from the gist of the present invention.

Although the present invention has been described above with reference to the embodiments, the present invention is not limited to the above examples. Various changes that can be understood by those skilled in the art can be made within the scope of the present invention in terms of the structure and details of the present invention.

The present invention can be used in combination with an optical system simulator such as an optical communication or a camera. The present invention can also be used in combination with other simulators such as ultimate condition search in engineering design.

This application claims priority on the basis of Japanese application Japanese Patent Application No. 2017-129002 filed on June 30, 2017, and incorporates all of its disclosures herein.

10 Designed optical system (telescope)
11 1st mirror 12 2nd mirror 13 3rd mirror 14 Plane mirror 15 Light receiving surface 16 Normal light path 17 Stray light path 100, 100A Optical system design support system 100B Rare event analyzer 200, 200A, 200B Processing unit 210 , 210A Parameter determination unit 212, 212B histogram calculation unit 214, 214A, 214B Random distribution definition unit 216, 216B Parameter generation unit 220 Trajectory calculation unit 220B Simulator 230 Evaluation value calculation unit 230B Evaluation value acquisition unit 300, 300B Storage unit 310, 310B Parameter storage unit 320 Optical model storage unit 330, 330B Evaluation value storage unit

Claims (10)

A rare event analysis device that discovers the possibility of rare events occurring.
_{Based on the proposed distribution Q (X), N parameters X k} , each of which defines an event, to be used for N simulations in each test k (1 ≦ k ≦ K) are determined, and the above N parameters are determined. A parameter determination unit that sequentially supplies the parameters X _k of the above to a simulator having an evaluation function C (X) that defines the possibility of the rare event occurring.
Wherein in response to N parameters X _k, and an evaluation value acquiring unit that acquires an evaluation value of N that is output from the simulator C (X _{k),
The parameter determination unit has N simulations of the N past parameters X k-1} used in the N simulations of the past test (k-1) and the N past parameters X _k-1. Based on the N past evaluation values C (X _k-1 ) obtained in, the N parameters X _k to be used for N simulations in this test k were determined.
The parameter determination unit
A histogram in which the N past parameters X _k-1 and the N past evaluation values C (X _k-1 ) are input, the evaluation value is a class, and the number of parameters of each class is a frequency. Histogram calculation unit that calculates
Based on the calculated histogram, the distribution of random numbers for generating the _{N parameters X k} used in the N simulations in the current test k is _{set to the threshold Cth determined for the class of the histogram.} A random number distribution definition unit that defines a random number that gives a histogram that gives a histogram in which the frequency of the parameters in each of the above classes is greater than the frequency of the parameters in each class below the threshold _Cth.
Based on the defined random number distribution, a parameter generation unit that generates _{each parameter X k used for each simulation in this test k, and a parameter generation unit.}
That we have a scarce event analysis apparatus. The random number distribution definition unit, within a predetermined error range, such as the sum of the frequencies of the sum and the threshold value the following parameters of the frequency threshold or more parameters determined for the class of the histogram is equal properly as random numbers to provide a histogram is obtained, defining the distribution, rare event analysis apparatus according to claim 1. The rare event analysis device according to claim 2 , wherein when the maximum value of the class is normalized to be 1.0, the random number distribution definition unit sets the threshold value to a value larger than 0.5. .. The rare event analysis device according to claim 3 , wherein the random number distribution definition unit sets the threshold value in the range of 0.6 or more and 0.8 or less. The rare event analysis device according to claim 4 , wherein the random number distribution definition unit sets the threshold value to a value of 0.6. A rare event analysis method that discovers the possibility of rare events occurring by computer.
_{Based on the proposed distribution Q (X), N parameters X k} , each of which defines an event, to be used for N simulations in each test k (1 ≦ k ≦ K) are determined, and the above N parameters are determined. The parameters X _k of the above are sequentially supplied to the simulator having the evaluation function C (X) that defines the possibility of the rare event occurring, and the parameters are determined.
Including the acquisition of _{N evaluation values C (X k} ) output from the simulator corresponding to the N parameters X _{k.
Determining the parameters is N times for the N past parameters X k-1} used in the N simulations of the past test (k-1) and the N past parameters X _k-1 . based on the past evaluation value C (X _k-1) of the N acquired by the simulation state, and are able to determine the N parameters X _k to be used in the simulation of N times in this study k ,
Determining the parameters
A histogram in which the N past parameters X _k-1 and the N past evaluation values C (X _k-1 ) are input, the evaluation value is a class, and the number of parameters of each class is a frequency. To calculate and
Based on the calculated histogram, the distribution of random numbers for generating the _{N parameters X k} used in the N simulations in the current test k is _{set to the threshold value C th determined for the class of the histogram.} It is defined so that a random number can be obtained that gives a histogram in which the frequency of the parameter in each of the above classes is greater than the frequency of the parameter in each class below the threshold _Cth.
Based on the defined random number distribution, each parameter X _k used for each simulation in the present test k is generated, and
Including scarce event analysis methods. The definition is a histogram in which the total frequency of the parameters above the threshold value and the total frequency of the parameters below the threshold value are equal within a predetermined error range. The rare event analysis method according to claim 6, wherein the distribution is defined so that a given random number can be obtained. The rare event analysis method according to claim 7, wherein the definition sets the threshold value to a value greater than 0.5 when the maximum value of the class is normalized to be 1.0. The definition is the rare event analysis method according to claim 8, wherein the threshold value is set in the range of 0.6 or more and 0.8 or less. A rare event analysis program that discovers the possibility that a rare event may occur in a computer, and the rare event analysis program is used in the computer.
_{Based on the proposed distribution Q (X), N parameters X k} , each of which defines an event, to be used for N simulations in each test k (1 ≦ k ≦ K) are determined, and the above N parameters are determined. A parameter determination procedure for sequentially supplying the parameters X _k of the above to a simulator having an evaluation function C (X) that defines the possibility of the rare event occurring, and
Wherein in response to N parameters X _k, the evaluation value obtaining step of obtaining the evaluation value of N which is output from the simulator C (X _k), is executed,
The parameter determination procedure is performed on the computer with respect _{to the N past parameters X k-1} used in the N simulations of the past test (k-1) and the N past parameters X _k-1. Based on the N past evaluation values C (X _k-1 ) obtained in the N simulations, the N parameters X _k used for the N simulations in the current test k are determined. Yes,
The parameter determination procedure is performed on the computer.
A histogram in which the N past parameters X _k-1 and the N past evaluation values C (X _k-1 ) are input, the evaluation value is a class, and the number of parameters of each class is a frequency. Histogram calculation procedure to calculate
Based on the calculated histogram, the distribution of random numbers for generating the _{N parameters X k} used in the N simulations in the current test k is _{set to the threshold Cth determined for the class of the histogram.} A random number distribution definition procedure and a procedure for defining a random number distribution so as to obtain a random number that gives a histogram in which the frequency of the parameter in each of the above classes is greater than the frequency of the parameter in each class below the threshold _Cth.
Based on the defined random number distribution, a parameter generation procedure for generating _{each parameter X k used for each simulation in the present test k, and a parameter generation procedure.}
Nozomi small event analysis program that Ru allowed to run.

Images