JP6945492B2 - Analysis parameter estimation method - Google Patents

Description

本発明は、数値シミュレーションにおける解析パラメータの推定方法に関する。 The present invention relates to a method for estimating analysis parameters in numerical simulation.

近年の計算機技術の進展により、物理現象を模擬するための数値解析あるいは数値シミュレーションが、様々な分野で広く用いられるようになっている。工業的には、ある装置や機械が、その目的に適った性能を発揮できるかを事前に評価することを目的として、装置や機械の設計段階において事前に数値シミュレーションを実施し、その結果に基づいて装置や機械の設計を最適化することが広く行われている。 Due to recent advances in computer technology, numerical analysis or numerical simulation for simulating physical phenomena has come to be widely used in various fields. Industrially, in order to evaluate in advance whether a certain device or machine can exhibit the performance suitable for the purpose, numerical simulation is performed in advance at the design stage of the device or machine, and based on the result. It is widely practiced to optimize the design of equipment and machines.

こうした工業的な装置や機械における物理現象を対象とした数値シミュレーションを精度よく実行するためには、解析プログラムに与える初期条件や境界条件、ないしは解析プログラムの中で物理現象を記述するために用いられる物理モデルが含有するモデルパラメータなどの解析パラメータを適切に与える必要がある。 In order to accurately execute numerical simulations targeting physical phenomena in such industrial equipment and machines, it is used to describe the initial conditions and boundary conditions given to the analysis program, or to describe the physical phenomena in the analysis program. It is necessary to appropriately give analysis parameters such as model parameters contained in the physical model.

例えば、ある燃焼装置内の温度分布を評価するために、ナビエ・ストークス方程式や熱伝導方程式、燃焼反応のモデル式などに基づいた燃焼・流体シミュレーションを行うことができるが、その場合、境界条件としての燃焼装置入口での燃料ガスの流速および化学組成の分布や、燃焼反応のモデル式が含有するモデルパラメータ、例えば反応速度定数などの解析パラメータを適切に与えない限り、実用に足る解析精度を得ることはできない。 For example, in order to evaluate the temperature distribution in a certain combustion device, combustion / fluid simulation based on the Navier-Stokes equation, the heat conduction equation, the model formula of the combustion reaction, etc. can be performed, but in that case, as a boundary condition. Unless the distribution of the flow velocity and chemical composition of the fuel gas at the inlet of the combustion device and the model parameters contained in the model formula of the combustion reaction, for example, the analysis parameters such as the reaction rate constant are appropriately given, the analysis accuracy sufficient for practical use can be obtained. It is not possible.

しかしながら、一般に実際の工業的な装置や機械を対象とする場合、これらの解析パラメータを事前に全て正確に知ることは困難であることが多い。その場合、事前に正確に知ることができない解析パラメータについては、何らかの仮定値を与えて数値シミュレーションを実行することになる。そのため、数値シミュレーション自体を実行することは可能であっても、得られた結果が目的に適う十分な精度を有しているとは限らない。 However, in general, when targeting an actual industrial device or machine, it is often difficult to know all of these analysis parameters accurately in advance. In that case, for the analysis parameters that cannot be accurately known in advance, some hypothetical values are given and the numerical simulation is executed. Therefore, even if it is possible to execute the numerical simulation itself, the obtained result does not always have sufficient accuracy suitable for the purpose.

すなわち、装置や機器の設計に足る十分な精度において物理現象を模擬するための数値シミュレーションを実行するためには、未知の解析パラメータの値を的確に推定する必要がある。 That is, in order to execute a numerical simulation for simulating a physical phenomenon with sufficient accuracy for designing a device or equipment, it is necessary to accurately estimate the value of an unknown analysis parameter.

上記の課題に対しては、逆問題（inverse problem）あるいは逆解析（inverse analysis）の有用性が近年注目されている。逆問題あるいは逆解析とは、入力から出力を、または原因から結果を直接的に求める順問題（direct problem）あるいは順解析（direct analysis）とは逆の概念であり、出力から入力を、または結果から原因を求めようとするものであり、その概要については例えば非特許文献１に記述されている。 In recent years, attention has been paid to the usefulness of an inverse problem or an inverse analysis for the above problems. Inverse problem or inverse analysis is the opposite concept of direct problem or direct analysis, which seeks the effect from the input or the result directly from the cause, and the input from the output or the result. The cause is sought from the above, and the outline thereof is described in, for example, Non-Patent Document 1.

前記の燃焼装置内の温度分布を評価する問題についても、こうした逆解析手法を適用することにより、出力である燃焼装置の出口境界での燃焼ガスの温度および化学組成などの測定結果から、入力である燃焼装置の入口境界での燃料ガスの流速および化学組成の分布や、燃焼反応のモデルパラメータなどを求めることができると期待される。 By applying such an inverse analysis method to the problem of evaluating the temperature distribution in the combustion device, the measurement results such as the temperature and chemical composition of the combustion gas at the outlet boundary of the combustion device, which is the output, can be input. It is expected that the flow velocity of fuel gas and the distribution of chemical composition at the inlet boundary of a certain combustion device, the model parameters of the combustion reaction, etc. can be obtained.

実際、逆解析手法は様々な工業分野への応用が行われている。例えば非特許文献２では、自動車用ディスクブレーキに用いられるブレーキパッド摩擦材に関して、直接的に計測することが困難である摩擦材の弾性定数を、逆解析手法を適用することにより同定する手法について記述されている。具体的には、まず打撃試験によりブレーキパッドの固有振動数とモード形状を測定した上で、弾性定数を適当に仮定した数値シミュレーション（有限
要素解析）を多数実行して、計算結果を入力（推定対象パラメータ）である弾性定数と出力（計測可能な観測量）であるブレーキパッドの固有振動数およびモード形状の関係性として定式化し、最後に最尤推定法を用いて測定データと計算結果の誤差が最小となるよう弾性定数を同定する手法について記述されている。 In fact, the inverse analysis method has been applied to various industrial fields. For example, Non-Patent Document 2 describes a method for identifying a brake pad friction material used for an automobile disc brake by applying an inverse analysis method to identify the elastic constant of the friction material, which is difficult to measure directly. Has been done. Specifically, after first measuring the natural frequency and mode shape of the brake pad by a striking test, a large number of numerical simulations (finite element analysis) assuming elastic constants are executed, and the calculation results are input (estimated). Formulated as the relationship between the elastic constant (target parameter) and the natural frequency and mode shape of the brake pad, which is the output (measurable observation amount), and finally the error between the measured data and the calculation result using the most probable estimation method. The method of identifying the elastic constant so that is minimized is described.

また非特許文献３では、逆解析手法によるパラメータ推定問題における実験方法あるいは計測条件、例えば計測点の位置や計測する物理量などを最適化する手法について記述されている。具体的には、事後推定誤差の共分散行列の最大固有値に着目し、その値を最小化する実験方法あるいは計測条件を採用することで、効果的なパラメータ推定が実現できることが記述されている。 Further, Non-Patent Document 3 describes an experimental method or a measurement condition in a parameter estimation problem by an inverse analysis method, for example, a method of optimizing the position of a measurement point or a physical quantity to be measured. Specifically, it is described that effective parameter estimation can be realized by focusing on the maximum eigenvalue of the covariance matrix of the ex-post estimation error and adopting an experimental method or measurement condition that minimizes the value.

久保司郎、「逆問題の解析手法」、材料、日本材料学会、１９９２年、第４１巻、第４７０号、ｐ．１５９５−１６０４Shiro Kubo, "Inverse Problem Analysis Method", Materials, Japan Society of Materials Science, 1992, Vol. 41, No. 470, p. 1595-1604 井上正則ほか、「逆解析を用いたディスクブレーキパッド摩擦材の異方性弾性定数の非侵襲同定」、日本機械学会論文集、２０１６年、第８２巻、第８４０号、ｐ．１−１５Masanori Inoue et al., "Non-invasive identification of anisotropic elastic constants of disc brake pad friction materials using inverse analysis", Proceedings of the Japan Society of Mechanical Engineers, 2016, Vol. 82, No. 840, p. 1-15 Ameyaほか、「Optimization of Measurements For Inverse Problem」、3rd International Conference on Inverse Problems in Engineering、The American Society of Mechanical Engineers、１９９９年Ameya et al., "Optimization of Measurements For Inverse Problem", 3rd International Conference on Inverse Problems in Engineering, The American Society of Mechanical Engineers, 1999

しかしながら、一般に実際の工業的な装置や機械を対象とする場合、これらの解析パラメータを事前に全て正確に知ることは困難であることが多い。例えば前記の例では、装置や機械の形状によって、境界条件としての流路入口の流速および圧力分布を正確に測定することが物理的に困難な場合もあるし、また物性の良くわかっていない新たな流体などを対象とする場合のように、密度や粘性係数などのパラメータを事前に把握することが困難な場合もある。 However, in general, when targeting an actual industrial device or machine, it is often difficult to know all of these analysis parameters accurately in advance. For example, in the above example, depending on the shape of the device or machine, it may be physically difficult to accurately measure the flow velocity and pressure distribution at the flow path inlet as boundary conditions, and the physical properties are not well understood. In some cases, it may be difficult to grasp parameters such as density and viscosity coefficient in advance, as in the case of targeting various fluids.

こうした場合、事前に正確に知ることができない解析パラメータについては、何らかの仮定値を与えて数値シミュレーションを実行することになる。そのため、数値シミュレーション自体を実行することは可能であっても、得られた結果が目的に適う十分な精度を有しているとは限らない、という問題があった。すなわち、装置や機器の設計に足る十分な精度において物理現象を模擬するための数値シミュレーションを実行するためには、未知の解析パラメータの値を的確に推定する必要がある、という課題があった。 In such a case, numerical simulations are performed by giving some hypothetical values for analysis parameters that cannot be accurately known in advance. Therefore, even if it is possible to execute the numerical simulation itself, there is a problem that the obtained result does not always have sufficient accuracy suitable for the purpose. That is, there is a problem that it is necessary to accurately estimate the value of an unknown analysis parameter in order to execute a numerical simulation for simulating a physical phenomenon with sufficient accuracy for designing a device or an apparatus.

これらの非特許文献に記述されている概念や手法は、様々な工業的な装置や機械を対象とした数値シミュレーションにおいて、目的に適う十分な精度が得られるように、未知の解析パラメータの値を的確に推定するために適用することが可能と考えられるが、実用上はいくつかの課題がある。 The concepts and methods described in these non-patent documents use the values of unknown analytical parameters to obtain sufficient accuracy for the purpose in numerical simulations of various industrial devices and machines. Although it is considered possible to apply it for accurate estimation, there are some problems in practical use.

まず、前記の非特許文献２においては、有限要素解析の計算コスト削減のために、入力（推定対象パラメータ）と出力（計測可能な観測量）の関係性を、予め実行しておいた数値シミュレーションの結果に基づいて、線形あるいは非線形の多項式を用いて応答曲面として定式化しておき、それを実際の観測結果と照らし合わせることで、未知の（推定対象の）パラメータの同定が行われている。 First, in Non-Patent Document 2 described above, in order to reduce the calculation cost of finite element analysis, a numerical simulation in which the relationship between an input (estimated parameter) and an output (measurable observation amount) is executed in advance is performed. Based on the result of, an unknown (estimated target) parameter is identified by formulating it as a response curved surface using a linear or non-linear polynomial and comparing it with the actual observation result.

ここで、非特許文献２にも記載されているように、応答曲面の精度はパラメータの同定
精度に大きく影響するため、問題の性質に適した定式化の手法を選択することが重要となる。具体的には、問題における入力と出力の関係性を適切に表現できるよう、例えば応答曲面の多項式の次数を適切に選択することが必要となる。 Here, as described in Non-Patent Document 2, the accuracy of the response surface greatly affects the identification accuracy of the parameters, so it is important to select a formulation method suitable for the nature of the problem. Specifically, it is necessary to appropriately select the degree of the polynomial of the response surface, for example, so that the relationship between the input and the output in the problem can be appropriately expressed.

しかしながら、推定対象であるパラメータの数（次元）が大きい場合、応答曲面の多項式の次数を増やしていくと、応答曲面を構成する多項式の係数の数が爆発的に増大してしまう。例えば、推定対象であるパラメータの数が２０個ある場合、それらのパラメータと出力（観測量）の関係性を多項式で表すとすると、１次多項式の場合は係数の数は２１個となるが、２次多項式の場合は２３１個、３次多項式の場合は１５８１個となる。 However, when the number (dimension) of the parameter to be estimated is large, increasing the degree of the polynomial of the response surface explosively increases the number of coefficients of the polynomial constituting the response surface. For example, if the number of parameters to be estimated is 20, and the relationship between those parameters and the output (observed amount) is expressed by a polynomial, the number of coefficients is 21 in the case of a first-order polynomial. In the case of a quadratic polynomial, the number is 231, and in the case of a cubic polynomial, the number is 1581.

応答曲面の定式化にあたって、単純にこれらの多項式の係数を数値シミュレーションの計算結果に基づいて、例えば最小二乗回帰などの手法により決定するためには、係数の数に応じた多数の数値シミュレーションを、入力条件を変化させた上で実行する必要がある。すなわち、１次多項式を採用した場合には最低２１回の数値シミュレーションを実行すれば係数を決定することができるが、２次多項式を採用した場合は最低でも２３１回、３次多項式を採用した場合には最低でも１５８１回以上の数値シミュレーションを実行しない限り、応答曲面の係数を決定することができない。 In formulating the response surface, in order to simply determine the coefficients of these polynomials based on the calculation results of the numerical simulation by a method such as least squares regression, a large number of numerical simulations according to the number of coefficients are performed. It is necessary to execute after changing the input conditions. That is, when the first-order polynomial is adopted, the coefficient can be determined by executing the numerical simulation at least 21 times, but when the second-order polynomial is adopted, at least 231 times when the third-order polynomial is adopted. The coefficient of the response surface cannot be determined unless a numerical simulation of at least 1581 times is performed.

また、非特許文献３においては、こうした逆解析手法によるパラメータ推定問題において、パラメータの推定精度を高めるために観測量（計測点の位置や計測する物理量など）を最適化する手法について記述されている。しかしながら、「物理現象を模擬するための数値シミュレーションを所望の解析精度において実行する」という本来の目的に対しては、パラメータの推定精度を過度に高めることを志向することは必ずしも現実的ではない。つまり、パラメータの推定精度を高めるためには、より高い精度で前述の応答曲面を作成する必要が生じるため、その作成のために必要な数値シミュレーションの実行回数も増やす必要がある。 Further, Non-Patent Document 3 describes a method of optimizing an observable quantity (position of a measurement point, a physical quantity to be measured, etc.) in order to improve the estimation accuracy of the parameter in the parameter estimation problem by such an inverse analysis method. .. However, for the original purpose of "performing a numerical simulation for simulating a physical phenomenon with a desired analysis accuracy", it is not always realistic to aim at excessively increasing the estimation accuracy of parameters. That is, in order to improve the estimation accuracy of the parameters, it is necessary to create the above-mentioned response surface with higher accuracy, and therefore it is necessary to increase the number of executions of the numerical simulation required for the creation.

したがって、現実的に実行可能な数値シミュレーションの実行回数のもとで、できるだけ高次の非線形性を表現できる応答曲面の定式化手法と、得られた応答曲面が、目的とする数値シミュレーションにおいて所望の解析精度を得る上で十分な精度を有しているか否かを適切に評価できる、応答曲面および観測量の最適化手法とが必要となるが、そのような手法はこれまでに実用化されていなかった。 Therefore, a method for formulating a response surface that can express as high-order non-linearity as possible based on the number of executions of a realistically feasible numerical simulation and the obtained response surface are desired in the target numerical simulation. A response surface and observation amount optimization method that can appropriately evaluate whether or not the analysis accuracy is sufficient is required, and such a method has been put into practical use so far. There wasn't.

そのため、数値シミュレーションの計算コストが大きく、かつ問題を支配する物理現象の非線形性が強い問題、例えば乱流解析や燃焼解析などの問題に対して、このような応答曲面法に基づく逆解析手法を適用して、目的とする数値シミュレーションを所望の解析精度において実行するために最低限必要な精度において、未知の解析パラメータの値を的確に推定することは、これまでは実用的に困難であった。 Therefore, for problems with high computational cost of numerical simulation and strong non-linearity of physical phenomena that govern the problem, such as turbulence analysis and combustion analysis, an inverse analysis method based on such a response curved surface method is used. It has been practically difficult to accurately estimate the values of unknown analysis parameters with the minimum accuracy required to apply and execute the desired numerical simulation with the desired analysis accuracy. ..

そこで本発明では、物理現象を模擬するための数値シミュレーションを所望の解析精度において実行するために必要な精度において、未知の解析パラメータの値を的確に推定できるとともに、数値シミュレーションの計算コストが大きく、かつ問題を支配する物理現象の非線形性が強い問題に対しても適用可能な、解析パラメータの推定方法を提供することを課題とする。 Therefore, in the present invention, the values of unknown analysis parameters can be accurately estimated with the accuracy required to execute the numerical simulation for simulating the physical phenomenon with the desired analysis accuracy, and the calculation cost of the numerical simulation is large. Moreover, it is an object of the present invention to provide a method for estimating analysis parameters, which can be applied to a problem in which the physical phenomenon that governs the problem is strongly non-linear.

上記課題を解決するため、本発明の第１の態様である解析パラメータの推定方法は、物理現象を模擬するための数値シミュレーションを所望の解析精度において実行するために必要な、１つ以上の未知の解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）の推定値Ｘ^＊ _ｉを推定する手法であって、対象とする物理現象における、１つ以上の観測点における１
つ以上の観測変数Ｙ_ｊ（ｊ＝１，２，・・・，Ｍ）について、前記１つ以上の解析パラメータＸ_ｉの値を変化させた多数の数値シミュレーションを実行して得られた結果に基づいて、前記Ｘ_ｉと前記Ｙ_ｊとの関係性を定式化し、さらに前記物理現象について実験的計測を行うことにより、前記Ｙ_ｊの観測値Ｙ^＊ _ｊ（ｊ＝１，２，・・・，Ｍ）を得て、前記Ｙ^＊ _ｊの統計情報と、前記定式化された前記Ｘ_ｉと前記Ｙ_ｊの関係性に基づいて、ベイズ推定を行うことにより、前記観測値Ｙ^＊ _ｊが得られたもとでの、前記パラメータＸ_ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）＝Ｐ（Ｘ_ｉ｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）として、前記パラメータＸ_ｉの推定値Ｘ^＊ _ｉの確率分布を得ることを特徴とする。 In order to solve the above problems, the method for estimating analysis parameters according to the first aspect of the present invention is one or more unknowns necessary for executing a numerical simulation for simulating a physical phenomenon with a desired analysis accuracy. ^{It is a method of estimating the estimated value X *} _i of the analysis parameter X _i (i = 1, 2, ..., N) of the above, and 1 at one or more observation points in the target physical phenomenon.
One or more of the observed variables _{Y j (j = 1,2, ···} , M) for, on the results obtained by executing a number of numerical simulation of changing the value of the one or more analysis parameters X _i based on the _{X i} and the relationship between the _{Y j} formulated by performing an experimental measurement for further said physical phenomenon, said _{Y j} of observations _{^{Y * j (j = 1,2,}} ··· , M) is obtained, ^{and the observed value Y *} _j is obtained by performing Bayesian estimation based on the statistical information of the ^{Y *} _{j and the relationship between the formulated X i} and the Y _j. in Moto that is, the parameter _{X i} of the posterior probability distribution ^{_{P (X * i) = P}} (X i | Y * 1, Y * 2, ···, Y * M) as the estimated value of the parameter _{X i} It is characterized by obtaining a probability distribution of X ^* _i.

このように構成すると、解析パラメータＸ_ｉの値を変化させた多数の数値シミュレーションを実行して得られた結果に基づいて定式化された、解析パラメータＸ_ｉと観測変数Ｙ_ｊの関係性に基づいて、実際に得られた観測値Ｙ^＊ _ｊの統計情報を適切に説明できる解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの確率分布Ｐ（Ｘ^＊ _ｉ）を得ることができる。 With this configuration, it is based on the relationship between the _{analysis parameter X i} and the observation variable Y _j , which is formulated based on the results obtained by executing a large number of numerical simulations in which the values of the _{analysis parameter X i are changed.} Therefore, it is possible to obtain the probability distribution P (X ^* _i ) of the estimated value X ^* _{i of the analysis parameter X i that} can appropriately explain the statistical information of the actually obtained observed value Y ^* _j.

また、本発明の第２の態様である解析パラメータの推定方法では、前記Ｘ_ｉと前記Ｙ_ｊとの関係性は、Ｙ_ｊ（ｊ＝１，２，・・・，Ｍ）を目的変数とし、Ｘ_ｉ（ｉ＝１，２，・・・，Ｎ）を説明変数とした多変数多項式回帰により定式化されることを特徴とする。 Further, in the analysis parameter estimation method according to the second aspect of the present invention, the relationship between the _{X i} and the Y _{j has Y j} (j = 1, 2, ..., M) as the objective variable. , X _i (i = 1, 2, ..., N) as explanatory variables, and is formulated by multivariate polynomial regression.

このように構成すると、解析パラメータＸ_ｉと観測変数Ｙ_ｊの関係性を、その関係性に応じて適切な次数および／または項数の多変数多項式を選択することにより、精度よく表現することができる。 With this configuration _{, the relationship between the analysis parameter X i} and the observed variable Y _j can be accurately expressed by selecting an appropriate multivariate polynomial of degree and / or number of terms according to the relationship. can.

また、本発明の第３の態様である解析パラメータの推定方法は、前記多変数多項式回帰において、リッジ回帰、ラッソ回帰またはエラスティックネットのいずれかからなるスパース推定手法を用いて変数選択を行うことを特徴とする。 Further, in the method for estimating analysis parameters according to the third aspect of the present invention, in the multivariate polynomial regression, variable selection is performed using a sparse estimation method consisting of ridge regression, lasso regression, or elastic net. It is characterized by.

このように構成すると、前記多変数多項式を構成する変数および／または項の候補となる数が多い場合においても、解析パラメータＸ_ｉと観測変数Ｙ_ｊの関係性を、変数の数および／または項の数を最小限に抑制しつつ、かつ精度よく表現することができるため、計算負荷を抑制しつつ所望の解析精度を得ることができる。 With this configuration, even when there are many candidates for the variables and / or terms that make up the multivariate polynomial, _{the relationship between the analysis parameter X i} and the observed variable Y _j can be determined by the number of variables and / or terms. Since it is possible to express accurately while suppressing the number of the above to the minimum, it is possible to obtain a desired analysis accuracy while suppressing the calculation load.

また、本発明の第４の態様である解析パラメータの推定方法は、前記ベイズ推定において、マルコフ連鎖モンテカルロ法またはハミルトニアンモンテカルロ法を用いることを特徴とする。 The analysis parameter estimation method according to the fourth aspect of the present invention is characterized in that the Markov chain Monte Carlo method or the Hamiltonian Monte Carlo method is used in the Bayesian estimation.

このように構成すると、前記多変数多項式における変数の数および／または項の数が多い場合や、あるいは解析パラメータＸ_ｉと観測変数Ｙ_ｊの関係性が強い非線形性を有する場合、さらには観測値Ｙ^＊ _ｊの統計分布が正規分布に従わないような場合など、観測値Ｙ^＊ _ｊの統計情報から直接的に前記パラメータＸ_ｉの推定値Ｘ^＊ _ｉを推定することが困難である場合であっても、観測値Ｙ^＊ _ｊの統計情報に基づいて、前記パラメータＸ_ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）＝Ｐ（Ｘ_ｉ｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）を数値的に求めることが可能となる。 With this configuration, the multi-variable or when the number of the number and / or section of the variable is large in the polynomial, or analysis when the parameter X _i and observation relation of the variable Y _j has strong nonlinearity, further observations When it is difficult to estimate the estimated value X ^* _{i of the parameter X i} directly from the statistical information of the ^{observed value Y *} _j , such as when the statistical distribution of ^{Y *} _{j does not follow the normal distribution. However} , based on the statistical information of the observed value Y ^* j, the posterior probability distribution P (X ^* _i ) = P (X _i | Y ^* _1, Y ^* _2, ... _, Y ^* _{M ) of the parameter X i.} ) Can be calculated numerically.

また、本発明の第５の態様である解析パラメータの推定方法は、前記１つ以上の観測点は、前記パラメータＸ_ｉと相関の高い観測変数Ｙ_ｊ１（ｊ_１＝１，２，・・・，Ｍ_１）が得られるように選ばれた観測点と、前記数値シミュレーションにおいて、計算結果の評価に用いられる評価変数Ｙ_ｊ２（ｊ_２＝Ｍ_１＋１，Ｍ_１＋２，・・・，Ｍ）が得られるように選ばれた観測点とからなることを特徴とする。 Further, a method of estimating a fifth aspect of which analysis parameters of the present invention, the one or more observation points, the parameter _{X i} having high correlation with the observed variables _{Y j1 (j 1 = 1,2,} ··· , M _{1 ) and the evaluation variable Y j2} (j ₂ = M ₁ + 1, M ₁ + 2, ..., M) used to evaluate the calculation results in the numerical simulation. It is characterized by consisting of observation points selected so that

このように構成すると、前記未知の解析パラメータＸ_ｉを精度よく推定するために選ば
れた観測変数Ｙ_ｊ１と、前記数値シミュレーションにおける計算結果の精度を評価するために選ばれた評価変数Ｙ_ｊ２に基づいて、前記解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの確率分布Ｐ（Ｘ^＊ _ｉ）を効率的に得ることができる。 With this configuration, the observation variables Y _j1 chosen in order to estimate accurately the unknown analysis parameters X _i, the evaluation variable Y _j2 chosen to assess the accuracy of the calculation result in the numerical simulation Based on this, the probability distribution P (X ^* _i ) of the estimated value X ^* _{i of the analysis parameter X i} can be efficiently obtained.

さらに、本発明の第６の態様である解析パラメータの推定方法は、前記事後確率分布Ｐ（Ｘ_ｉ｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）に基づく、前記パラメータＸ_ｉの推定値Ｘ^＊ _ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）を、前記定式化された前記Ｘ_ｉと前記Ｙ_ｊとの関係性に入力することにより得られる、前記評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）＝Ｐ（Ｙ_ｊ２｜Ｘ^＊ _ｉ）＝Ｐ（Ｙ_ｊ２｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）について、前記事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間が、いずれかのｊ_２について所定の閾値を上回る場合に、前記観測変数の数Ｍ_１、Ｍ_２および／ないしＭを増やすことにより、前記事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間が、いずれのｊ_２についても所定の閾値を下回るように構成したことを特徴とする。 Further, the method for estimating the analysis parameter according to the sixth aspect of the present invention is based on the posterior probability distribution P (X _i | Y ^* _1, Y ^* _2, ... _, Y ^* _M ). _i of the estimated value ^X _{* i} posterior probability distribution P ^(X _{* i),} obtained by inputting to the relationship between the said formalized the _{X i Y j,} post of the evaluation variables _{Y j2} Regarding the probability distribution P (Y ^** _j2 ) = P (Y _j2 | X ^* _i ) = P (Y _j2 | Y ^* _1, Y ^* _2, ... _, Y ^* _M ), the posterior probability distribution P ( Y ^** _j2) credit period of, if above a predetermined threshold for any _{j 2,} to the number _M 1, _{M 2} and / absence of the observed variables by increasing the M, the posterior probability distribution P ( Y ^{** _j2)} credit section, characterized by being configured to be below a predetermined threshold for any j _2.

このように構成すると、前記評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間が所定の閾値を下回るように、すなわち目的とする数値シミュレーションの実行結果において前記評価変数Ｙ_ｊ２を所望の解析精度において得るために必要な精度において前記の解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの確率分布Ｐ（Ｘ^＊ _ｉ）を得ることができる。 With this configuration, the so trust interval evaluation variables Y _j2 posterior probability distribution P (Y ** _^j2) is below a predetermined threshold, i.e., the evaluation variable Y _j2 in the execution results of the numerical simulations of interest it is possible to obtain a probability distribution P of the estimated value X ^* _i of the analysis parameters X _i in precision needed to obtain the desired accuracy of analysis (X * ^_i).

一実施形態による未知の解析パラメータの推定方法のステップを示す図である。It is a figure which shows the step of the method of estimating the unknown analysis parameter by one Embodiment. 境界条件としての燃焼装置入口境界での燃料ガスの化学組成の分布が未知である場合における燃焼シミュレーションの一実施例を示す。An example of a combustion simulation in the case where the distribution of the chemical composition of the fuel gas at the boundary of the combustion device inlet as a boundary condition is unknown is shown. 解析パラメータＸ_１からＸ_１０までの事後確率分布Ｐ（Ｘ^＊ _ｉ）の推定結果の一例を示す。An example of the estimation result of the posterior probability distribution P (X ^* _i ) from the analysis parameters X _{1 to X 10 is shown.} 各観測点および評価点におけるガス温度の事後確率分布Ｐ（Ｙ^＊＊ _１）、Ｐ（Ｙ^＊＊ _５）、Ｐ（Ｙ^＊＊ _９）、Ｐ（Ｙ^＊＊ _１３）、Ｐ（Ｙ^＊＊ _１７）の一例を示す。Posterior probability distributions of gas temperature at each observation point and evaluation point P (Y ^** ₁ ), P (Y ^** ₅ ), P (Y ^** ₉ ), P (Y ^** ₁₃ ), P (Y ^**) _{An example of 17} ) is shown.

以下、図面を参照して、本発明の好適な実施形態について詳細に説明する。 Hereinafter, preferred embodiments of the present invention will be described in detail with reference to the drawings.

図１は、本発明の一実施形態による未知の解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）の推定方法のステップを示している。 Figure 1 is unknown analysis parameters according to an embodiment of the present invention _{X i (i = 1,2, ···} , N) shows a step of estimating methods.

はじめに、ステップ１として、解析対象に観測点および評価点を設定し、観測点における観測変数Ｙ_ｊ１（ｊ_１＝１，２，・・・，Ｍ_１）と、評価点における評価変数Ｙ_ｊ２（ｊ_２＝Ｍ_１＋１，Ｍ_１＋２，・・・，Ｍ）を定義する。解析対象は、例えば、ごみ焼却プラントやバイオマス燃焼プラントなどを含むが、これらに限られず、装置や機器の設計の過程において数値シミュレーションが用いられるあらゆる装置や機器が、解析対象となりうる。観測点の数と観測変数の数との関係は任意である。例えば、Ｍ_１個の観測点に各１つずつ合計Ｍ_１個の観測変数を定義してもよいし、１個の観測点にＭ_１個の観測変数を定義してもよいし、５個の観測点に各Ｍ_１／５個ずつ合計Ｍ_１個の観測変数を定義してもよい。 First, as step 1, the observation point and the evaluation point are set as the analysis target, and the observation variable Y _j1 (j ₁ = 1, 2, ..., M _{1 ) at the observation point and the evaluation variable Y j2} (M 1) at the evaluation point are set. Define j ₂ = M ₁ + 1, M ₁ + 2, ..., M). The analysis target includes, for example, a waste incineration plant and a biomass combustion plant, but is not limited to these, and any device or device for which numerical simulation is used in the process of designing the device or device can be analyzed. The relationship between the number of observation points and the number of observation variables is arbitrary. For example, it is possible to define each one at a total of M ₁ or observed variables to a _single observation point M, may be defined to M ₁ or observed variables in one observation point, five the observation points may define a total M ₁ or observed variables by each M _1/5 pieces.

後述するように、本実施形態における未知の解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）の推定は、観測点において得られる観測変数Ｙ_ｊ（ｊ＝１，２，・・・，Ｍ）の統計情報に基づいてなされるため、原理的には観測変数Ｙ_ｊの数Ｍが大きいほど、精度よく解析パラメータＸ_ｉを推定することができるが、一方で解析パラメータＸ_ｉを推定するための計算量が増大するという問題が生じる。 As will be described later, the estimation of the unknown analysis parameter X _{i (i = 1, 2, ..., N) in this embodiment is performed by the observation variable Y j} (j = 1, 2, ..., N) obtained at the observation point. - because they are made based on statistical information of M), as in principle a large number M of observed variables Y _j, can be accurately estimated analysis parameters X _i, the while the analysis parameters X _i The problem arises that the amount of calculation for estimation increases.

そのため本実施形態においては、観測変数Ｙ_ｊの数Ｍをできるだけ少なくした上で精度よく解析パラメータＸ_ｉを推定するために、以下に記すように、観測変数Ｙ_ｊ（ｊ＝１，２，・・・，Ｍ）を、観測変数Ｙ_ｊ１（ｊ_１＝１，２，・・・，Ｍ_１）と、評価変数Ｙ_ｊ２（ｊ_２＝Ｍ_１＋１，Ｍ_１＋２，・・・，Ｍ）の２つに分けて定義する。 Therefore, in the present embodiment, in order to estimate the analysis parameter X _{i accurately after reducing} the number M of the observation variable Y j as much as possible, the observation variable Y _j (j = 1, 2, ...・ ・, M) is the observation variable Y _j1 (j ₁ = 1, 2, ···, M ₁ ) and the evaluation variable Y _j2 (j ₂ = M ₁ + 1, M ₁ + 2, ···, M). It is defined in two parts.

観測変数Ｙ_ｊ１（ｊ_１＝１，２，・・・，Ｍ_１）については、対象とする物理現象を考慮して、推定すべき未知の解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）との相関が高い観測量として定義する。観測変数Ｙ_ｊ１の数Ｍ_１を適切に決定する方法については、後述する。 For the observed variable Y _j1 (j ₁ = 1, 2, ..., M _{1 ), an unknown analysis parameter X i} (i = 1, 2, ..., M 1) to be estimated in consideration of the target physical phenomenon.・, Defined as an observed quantity with a high correlation with N). A method for appropriately determining the number M ₁ of the observation variable Y _{j1 will be described later.}

評価変数Ｙ_ｊ２（ｊ_２＝Ｍ_１＋１，Ｍ_１＋２，・・・，Ｍ）については、対象とする数値シミュレーションにおいて、計算結果の評価を行う上で着目する変数、すなわち計算結果の精度を評価する対象となる変数として定義する。評価変数Ｙ_ｊ２の数Ｍ_２は、対象とする数値シミュレーションにおいて計算結果の評価を行う上で必要となる数とすればよい。 Regarding the evaluation variable Y _j2 (j ₂ = M ₁ + 1, M ₁ + 2, ..., M), the variable to be focused on when evaluating the calculation result in the target numerical simulation, that is, the accuracy of the calculation result is determined. Define it as a variable to be evaluated. The number M ₂ of the evaluation variable Y _j2 may be a number required for evaluating the calculation result in the target numerical simulation.

例えば前記の燃焼装置内の温度分布を評価する問題において、燃焼装置の入口境界での燃料ガスの流速および化学組成の分布が未知の解析パラメータである場合、観測変数Ｙ_ｊ１としては燃焼装置の入口境界にできるだけ近く、物理的に観測が可能な位置における装置内のガス流速および／または化学組成の計測値を採用することができる。また、評価変数Ｙ_ｊ２としては、装置内の代表的な箇所であり、かつ物理的に観測が可能である箇所における温度の計測値や、装置出口境界でのガス流速および／または化学組成の計測値を採用することができる。 For example, in the problem of evaluating the temperature distribution in the combustion device, when the flow velocity of the fuel gas and the distribution of the chemical composition at the inlet boundary of the combustion device are unknown analysis parameters, the observation variable Y _j1 is the inlet of the combustion device. Measurements of gas flow velocity and / or chemical composition within the device at locations as close to the boundary as possible and physically observable can be employed. Further, as the evaluation variable Y _j2 , the measured value of the temperature at a typical location in the apparatus and physically observable, and the measurement of the gas flow velocity and / or the chemical composition at the outlet boundary of the apparatus. The value can be adopted.

次に、ステップ２として、解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）を変化させた多数の数値シミュレーションを実行する。以下では、ステップ２において、解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）をＫ通りに変化させた数値シミュレーションを実行することとし、ｋ番目の計算条件に対応する解析パラメータＸ_ｉの値をＸ_ｉ ^ｋ（ｋ＝１，２，・・・，Ｋ）と記述する。また、ｋ番目の計算条件での数値シミュレーション結果における観測変数Ｙ_ｊ１および評価変数Ｙ_ｊ２の値を、それぞれＹ_ｊ１ ^ｋおよびＹ_ｊ２ ^ｋと記述する。計算条件の数Ｋを適切に決定する方法については、後述する。 Next, as step 2, _{a large number of numerical simulations in which the analysis parameters X i} (i = 1, 2, ..., N) are changed are executed. In the following, in step 2, _{the numerical simulation in which the analysis parameter X i} (i = 1, 2, ..., N) is changed in K ways is executed, and the analysis parameter X corresponding to the kth calculation condition is executed. describing values ^{_{i X i k (k = 1,2}} , ···, K) and. _{Further, the values of the observation variable Y j1} and the evaluation variable Y _j2 in the numerical simulation result under the kth calculation condition are described as _{Y j1} ^k and Y _j2 ^k , respectively. The method for appropriately determining the number K of the calculation conditions will be described later.

次に、ステップ３として、ステップ２で実行した数値シミュレーションの実行結果に基づいて、解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性を定式化する。 Next, in step 3, based on the execution results of numerical simulation performed in step 2, to formulate the relationship between the analysis parameters X _i observed variables Y _j1, and analysis parameters X _i and evaluation variables Y _j2.

本実施形態においては、解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性は、解析パラメータＸ_ｉについての多変数多項式回帰により定式化される。ここで、多変数多項式回帰式の次数については、模擬対象となる物理現象の複雑さや、解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性の複雑さに応じて決めるのがよい。 In the present embodiment, the relationship between the analysis parameters X _i observed variables Y _j1, and analysis parameters X _i and evaluation variables Y _j2 is formulated by multivariate polynomial regression of the analysis parameters X _i. Here, the order of the multivariate polynomial regression equation, complexity and physical phenomena to be simulated target, analysis parameters X _i and observed variables Y _j1, and the relationship of the complexity of evaluation variables Y _j2 and analysis parameters X _i It is better to decide according to it.

例えば、多変数多項式回帰式の次数を２とした場合において、解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性は次式のように定式化することができ、数値シミュレーションの計算条件であるＸ_ｉ ^ｋと、その実行結果であるＹ_ｊ１ ^ｋおよびＹ_ｊ２ ^ｋに基づいて、解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性を最もうまく説明できる次式の係数ａ、ｂ、ｃおよびｄを、最小二乗回帰や最尤推定法などの手法によって求める回帰問題に帰着する。多変数多項式回帰式の次数を１とした場合や、３以上とした場合についても、同様に定
式化することができる。 For example, in the case where the degree of a polynomial regression equation and 2, is that the relationship of the analysis parameters X _i and observed variables Y _j1, and analysis parameters X _i and evaluation variables Y _j2 is formulated as follows can, and _X ^{i k} is the calculation conditions of numerical simulation, based on the _{Y j1} ^k and _{Y j2} ^k which is the execution result, the analysis parameter _{X i} and observed variables _{Y j1,} and analysis parameters _{X i} and evaluation variables _{Y j2} The variables a, b, c and d of the following equations that can best explain the relationship between The same formula can be used when the degree of the multivariate polynomial regression equation is 1 or 3 or more.

しかしながら、解析パラメータＸ_ｉの数が大きいと、回帰係数ａ、ｂ、ｃおよびｄの数が増大する。あるいは、多変数多項式の次数が大きくなった場合にも、求めるべき回帰係数の数が増大する。その場合、後述するステップ６において、観測変数Ｙ_ｊ１および評価変数Ｙ_ｊ２の観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の統計情報に基づいて、解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの確率分布を得る処理の計算負荷が増大してしまうという問題が生じる。 However, when _{the number of analysis parameters X i} is large, the number of regression coefficients a, b, c and d increases. Alternatively, when the degree of the multivariable polynomial increases, the number of regression coefficients to be obtained also increases. In that case, in step 6 described later, the probability distribution of the estimated value X ^* _{i of the analysis parameter X i} is obtained based on the statistical information of the observed values Y ^* _j1 and Y ^* _{j2 of the observed variable Y j1} and the evaluation variable Y _j2. There arises a problem that the calculation load of processing increases.

この問題に対処するため、本実施形態においては、ステップ３における解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性の定式化に際して、リッジ回帰、ラッソ回帰またはエラスティックネットのいずれかからなるスパース推定手法を用いて説明変数の選択が行われる。 In order to deal with this problem, in the present embodiment, _{in formulating the relationship between the analysis parameter X i} and the observation variable Y _j1 and the analysis parameter X _i and the evaluation variable Y _j2 in step 3, ridge regression, lasso regression or Explanatory variables are selected using a sparse estimation method consisting of one of the elastic nets.

スパース推定手法とは、回帰モデルの損失関数に回帰係数のノルムに基づく正則化項を加えた正則化損失関数を最小化することにより、変数の数が多い場合において、推定の安定化とともに、真に意味のある変数を選択することが可能となる手法であり、リッジ回帰、ラッソ回帰またはエラスティックネット等の手法が広く用いられている。これらのスパース推定手法については公知であるのでその詳細は省略する。 The sparse estimation method is a method of minimizing the regularization loss function by adding a regularization term based on the norm of the regression coefficient to the loss function of the regression model, thereby stabilizing the estimation and making it true when the number of variables is large. It is a method that enables selection of meaningful variables, and methods such as ridge regression, lasso regression, and elastic net are widely used. Since these sparse estimation methods are known, the details thereof will be omitted.

次に、ステップ４として、スパース推定手法による多変数多項式回帰により定式化された解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の回帰式の精度を評価する。具体的には、ステップ２において実行した数値シミュレーションの全ての計算条件Ｘ_ｉ ^ｋについて、ステップ３で得られた回帰式に各Ｘ_ｉ ^ｋを代入することにより得られる予測値Ｙ^〜 _ｊ１ ^ｋおよびＹ^〜 _ｊ２ ^ｋの値と、シミュレーション結果として得られたＹ_ｊ１ ^ｋおよびＹ_ｊ２ ^ｋの値とを比較する。 Next, in step 4, to evaluate the regression accuracy of sparse estimation analysis parameters were formulated by multivariate polynomial regression by Method X _i and observed variables Y _j1, and the analysis parameters X _i evaluation variables Y _{j2. Specifically, for all the calculation conditions X i} ^k of the numerical simulation executed in step 2, the predicted values Y ^to _j1 ^k and Y _{obtained by substituting each X i} ^k into the regression equation obtained in step 3. ^{The values of ~} _j2 ^k are compared with the values of _{Y j1} ^k and Y _j2 ^k obtained as simulation results.

回帰式の精度の評価指標としては、例えば上記のＹ_ｊ１ ^ｋおよびＹ_ｊ２ ^ｋの値と、対応するＹ^〜 _ｊ１ ^ｋおよびＹ^〜 _ｊ２ ^ｋの値との二乗平均誤差を用いることができる。二乗平均誤差の値が所定の閾値よりも大きい場合は、シミュレーションの実行ケース数Ｋを増やすか、または多変量多項式における次数を増やした上で、再度ステップ２からステップ４を実行することにより、回帰式の精度を高めるのがよい。 As an evaluation index of the accuracy of the regression equation, for example, the root mean square error between the _{above values of Y j1} ^k and Y _j2 ^k and the corresponding values of Y ^to _j1 ^k and Y ^to _j2 ^{k can be used.} If the value of the root mean square error is larger than the predetermined threshold, the number of simulation execution cases K is increased, or the degree in the multivariate polynomial is increased, and then steps 2 to 4 are executed again to perform regression. It is better to improve the accuracy of the formula.

次に、ステップ５として、観測点における観測変数Ｙ_ｊ１の観測値Ｙ^＊ _ｊ１および評価点における評価変数Ｙ_ｊ２の観測値Ｙ^＊ _ｊ２を実験的に得る。 Next, as step 5, the observed value Y ^* _{j1 of the observation variable Y j1 at} the observation point and the observation value Y ^* _j2 of the evaluation variable Y _{j2 at} the evaluation point are experimentally obtained.

次に、ステップ６として、ステップ５で得られた観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の統計情報と、ステップ３で定式化された解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性とに基づいて、ベイズ推定を行うことにより、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２が得られたもとでの、解析パラメータＸ_ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）＝Ｐ（Ｘ_ｉ｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）を求める。 Next, as step 6, the statistical information of ^{the observed values Y *} _j1 and Y ^* _{j2 obtained in step 5, the analysis parameter X i} and the observation variable Y _j1 formulated in step 3, and the analysis parameter X _i The posterior probability distribution P (X ^* _i ) of the analysis parameter X _i ^{under the observation values Y *} _j1 and Y ^* _j2 obtained by performing Bayesian estimation based on the relationship of the evaluation variable Y _{j2 =} Find P (X _i | Y ^* _1, Y ^* _2, ... _, Y ^* _M ).

本ステップ６の目的は、実験的に得られた観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の情報から、それらを与える解析パラメータＸ^＊ _ｉを逆推定することである。 The purpose of this step 6 is to back-estimate the analysis parameters X ^* _i ^{that give them from the experimentally obtained information of the observed values Y *} _j1 and Y ^* _j2.

ここで、ステップ５において得られる観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２は、実際の装置または機器において実験的に計測される量であるため、実験上あるいは計測上の様々な要因を反映した観測誤差を含んでいる。そのため、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２は確定値として得られるのではなく、バラつきをもった統計分布として得られることが通常である。 Here, since the observed values Y ^* _j1 and Y ^* _j2 obtained in step 5 are quantities that are experimentally measured in an actual device or instrument, observation errors that reflect various experimental or measurement factors are reflected. Includes. Therefore, the observed values Y ^* _j1 and Y ^* _j2 are not usually obtained as definite values, but are usually obtained as statistical distributions with variations.

観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の分布におけるバラつきが小さい場合や、その分布が正規分布に従う場合は、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の各々の分布の平均値ないし中央値をそれらの代表値として確定的に用いることができる。この場合、解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性が線形関係にある場合は、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２を与える解析パラメータＸ^＊ _ｉを解析的に算出することができる。 If the variations in the distributions of the observed values Y ^* _j1 and Y ^* _j2 are small, or if the distribution follows a normal distribution, the average or median of the distributions of ^{the observed values Y *} _j1 and Y ^* _{j2 are representative of them.} It can be used deterministically as a value. In this case, if _{the relationship between the analysis parameter X i} and the observation variable Y _j1 and the relationship between the analysis parameter X _i and the evaluation variable Y _j2 are linear, the analysis parameter X ^* _i ^{that gives the observation values Y *} _j1 and Y ^* _j2. Can be calculated analytically.

しかしながら、解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性がより複雑である場合、例えばこれらが２次以上の多項式で表されるような非線形関係にある場合は、前記のように観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の各々の分布の平均値ないし中央値をそれらの代表値として確定的に用いることができる場合であっても、それらの観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２を与える解析パラメータＸ^＊ _ｉを解析的に算出することは不可能であることが多い。 However, if the relationship between the analysis parameters X _i and observed variables Y _j1, and analysis parameters X _i evaluation variable Y _j2 is more complex, for example, they are in non-linear relationship as represented by a quadratic or more polynomial In the case, even if ^{the average value or the median value of each distribution of the observed values Y *} _j1 and Y ^* _j2 can be deterministically used as their representative values as described above, those observed values Y It is often impossible to analytically calculate the analytical parameters X ^* _i ^{that give *} _j1 and Y ^* _j2.

あるいは、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の分布におけるバラつきが大きい場合や、その分布が正規分布に従わない場合には、そもそも観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の代表値を確定的に定義することができないため、たとえ解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性が線形関係にあったとしても、そのような観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２に対応する解析パラメータＸ^＊ _ｉを解析的に算出することは不可能である。 Alternatively, if there is a large variation in the distribution of the observed values Y ^* _j1 and Y ^* _j2 , or if the distribution does not follow the normal distribution, ^{the representative values of the observed values Y *} _j1 and Y ^* _j2 are definitely defined in the first place. Therefore, even if _{the relationship between the analysis parameter X i} and the observed variable Y _j1 and the analysis parameter X _i and the evaluation variable Y _j2 are linear, such observed values Y ^* _j1 and Y ^* It is impossible to analytically calculate the analysis parameter X ^* _{i corresponding to j2.}

そこで本実施形態による解析パラメータの推定方法は、ベイズ推定を用いることにより、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の統計情報に基づいて、それらを与える解析パラメータＸ^＊ _ｉの確率分布を、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２が得られたもとでの、解析パラメータＸ_ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）＝Ｐ（Ｘ_ｉ｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）として求めるように構成されている。 Therefore, the method of estimating the analysis parameters according to the present embodiment uses Bayesian estimation to obtain the probability distribution of the ^{analysis parameters X *} _i that give them based on the statistical information ^{of the observed values Y *} _j1 and Y ^* _j2. The posterior probability distribution P (X ^* _i ) = P (X _i | Y ^* _1, Y ^* _2, ... _, Y ^* _M ) of the analysis parameter X _i based on the obtained Y ^* _j1 and Y ^* _j2. It is configured to ask for.

具体的には、本実施形態においては、ステップ６におけるベイズ推定手法として、マルコフ連鎖モンテカルロ法またはハミルトニアンモンテカルロ法が用いられている。 Specifically, in the present embodiment, the Markov chain Monte Carlo method or the Hamiltonian Monte Carlo method is used as the Bayesian estimation method in step 6.

マルコフ連鎖モンテカルロ法は、ベイズ推定の枠組みにおいて、事後確率分布の計算を数値的に行うことができる手法として広く用いられており、メトロポリス−ヘイスティングアルゴリズム、ギブスサンプラーなどの手法の総称である。またハミルトニアンモンテカルロ法は、マルコフ連鎖モンテカルロ法における目標分布への収束性の改善を目的として近年開発され広く利用されている。これらの手法については公知であるのでその詳細は省略する。 The Markov chain Monte Carlo method is widely used in the framework of Bayesian estimation as a method that can calculate posterior probability distributions numerically, and is a general term for methods such as the Metropolis-Hastings algorithm and Gibbs sampler. The Hamiltonian Monte Carlo method has been recently developed and widely used for the purpose of improving the convergence to the target distribution in the Markov chain Monte Carlo method. Since these methods are known, the details thereof will be omitted.

これらのベイズ推定手法を用いることにより、前記多変数多項式における変数の数および／または項の数が多い場合や、あるいは解析パラメータＸ_ｉと観測変数Ｙ_ｊ１および／または評価変数Ｙ_ｊ２との関係性が強い非線形性を有する場合や、さらには観測変数Ｙ_ｊ１および／または評価変数Ｙ_ｊ２の統計分布が正規分布に従わないような場合など、観測変数Ｙ_ｊ１および評価変数Ｙ_ｊ２の観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の統計情報から直接的に前記解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉを推定することが困難な場合であっても、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の統計情報に基づいて、解析パラメータＸ_ｉの事後確率分布Ｐ
（Ｘ^＊ _ｉ）＝Ｐ（Ｘ_ｉ｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）を数値的に求めることが可能となる。 By using these Bayesian estimation technique, the multi-case number of variables the number of variables in the polynomial and / or terms are often and, or analyze relationship parameter X _i and the observed variables Y _j1 and / or evaluation variable Y _j2 If you have a strong nonlinearity, etc. Furthermore if the statistical distribution of the observed variables Y _j1 and / or evaluation variable Y _j2 is that do not follow a normal distribution, the observed values of the observation variables Y _j1 and evaluation variables Y _j2 Y ^* even directly the analysis parameters if _{X i} estimate ^X _{* i} is difficult to estimate the from _j1 and ^Y _{* j2} statistics, based on the observed value ^Y _{* j1} and ^Y _{* j2} statistics Te, the posterior probability distribution P of analysis parameters _{X i}
(X ^* _i ) = P (X _i | Y ^* _1, Y ^* _2, ... _, Y ^* _M ) can be calculated numerically.

次に、ステップ７として、ステップ６で得られた解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）を、ステップ３で定式化された解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性に入力することにより、評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）＝Ｐ（Ｙ_ｊ２｜Ｘ^＊ _ｉ）＝Ｐ（Ｙ_ｊ２｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）を得る。 Next, as step 7, the posterior probability distribution P (X ^* _i ) of the estimated value X ^* _{i of the analysis parameter X i} obtained in step 6 is subjected to the analysis parameter X _i and the evaluation variable Y formulated in step 3. by entering into relationships of _j2, the posterior probability distribution ^P of the evaluation variable ^{_{Y j2 (Y ** j2) =}} P (Y j2 | X * i) = P (Y j2 | Y * 1, Y * 2, ·・ ・_, Y ^* _M ) is obtained.

次に、ステップ８として、ステップ７で得られた評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間の評価を行う。 Next, as step 8, the credit interval of the posterior probability distribution P (Y ^** _{j2 ) of the evaluation variable Y j2 obtained in step 7 is evaluated.}

前述のように、本実施形態における評価変数Ｙ_ｊ２は、計算結果の評価を行う上で着目する変数、すなわち計算結果の精度を評価する対象となる変数として定義されている。したがって、観測変数Ｙ_ｊ１および評価変数Ｙ_ｊ２の観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２が得られたもとでの、評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間が、所定の閾値よりも小さくなっていれば、対象とする数値シミュレーションにおいて、所望の解析精度を得るために必要な精度において未知の解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）が求まったといえることになる。 As described above, the evaluation variable Y _j2 in the present embodiment is defined as a variable of interest in evaluating the calculation result, that is, a variable to be evaluated for the accuracy of the calculation result. Therefore, the credible interval of the posterior probability distribution P (Y ^** _j2 ) of the evaluation variable Y _j2 under the observation values Y ^* _j1 and Y ^* _{j2 of} the observation variable Y _j1 and the evaluation variable Y _{j2 is predetermined.} if smaller than the threshold value, the numerical simulation of interest, the posterior probability distribution P of the estimated value X ^* _i of unknown analysis parameters X _i in precision needed to obtain the desired accuracy of analysis (X ^* _i) Can be said to have been obtained.

ここで、評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間が、所定の閾値よりも小さくならなかった場合は、対象とする数値シミュレーションにおいて、解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）を、所望の解析精度を得るために必要な精度で求めるためには、情報が不足していることを意味している。したがって、そのような場合は、観測変数Ｙ_ｊ１の数Ｍ_１を増やした上で、再度ステップ１からステップ８までの手順を実行する。これにより、観測変数Ｙ_ｊ１から得られる情報を増やすことができ、その結果として、数値シミュレーションにおいて所望の解析精度を得るために必要な精度で、未知の解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）を求めることができる。 Here, if the credible interval of the posterior probability distribution P (Y ^** _j2 ) of the evaluation variable Y _j2 is not smaller than a predetermined threshold value, in the target numerical simulation, the estimated value X _{of the analysis parameter X i} ^* _i of the posterior probability distribution P a (X * _^i), in order to obtain accuracy needed to obtain the desired accuracy of analysis means that the information is missing. Therefore, in such a case, the procedure from step 1 to step 8 is executed again after increasing the number M _{1 of the observation variable Y j1.} As a result, the _{information obtained from the observation variable Y j1} can be increased, and as a result, the estimated value X ^* _i of the _{unknown analysis parameter X i has} the accuracy required to obtain the desired analysis accuracy in the numerical simulation. The posterior probability distribution P (X ^* _i ) can be obtained.

以下、図面を参照して、本発明の実施例について説明する。 Hereinafter, examples of the present invention will be described with reference to the drawings.

図２は、ある燃焼装置内の温度分布やガス組成分布を評価するために行う燃焼シミュレーションにおいて、境界条件としての燃焼装置入口境界での燃料ガスの化学組成の分布が未知である場合の一実施例を示している。 FIG. 2 shows a case where the distribution of the chemical composition of the fuel gas at the boundary of the combustion device inlet as a boundary condition is unknown in the combustion simulation performed to evaluate the temperature distribution and the gas composition distribution in a certain combustion device. An example is shown.

ここでは、図２に示されているように４つに分割されている燃焼装置の入口境界の各々の領域Ｆ_１、Ｆ_２、Ｆ_３、Ｆ_４における燃料ガスＧ_１、Ｇ_２、Ｇ_３、Ｇ_４の組成を表現するためのパラメータとして、合計２０個の解析パラメータＸ_ｉ（ｉ＝１〜２０）が定義されている。すなわち、ｆをある関数として、Ｇ_ｎ＝ｆ（Ｘ_ｉ）（ｎ＝１〜４、ｉ＝１〜２０）であり、これらＸ_ｉの値を、上述した実施形態による解析パラメータの推定方法によって推定するものとする。また、本実施例における数値シミュレーションの要求精度としては、燃焼装置出口のガス温度を±１０℃の精度で求めることが要求されているものとする。 Here, the fuel gases G ₁ , G ₂ , and G _{3 in the respective regions F 1} , F ₂ , F ₃ , and F ₄ of the inlet boundary of the combustion apparatus divided into four as shown in FIG. , as a parameter for expressing the composition of _{G 4,} a total of 20 pieces of analysis parameters _X i (i = _1~20) are defined. That is, with f as a function, G _n = f (X _i ) (n = 1 to 4, i = 1 to 20), and _{the values of these X i} are determined by the method of estimating the analysis parameters according to the above-described embodiment. It shall be estimated. Further, as the required accuracy of the numerical simulation in this embodiment, it is assumed that the gas temperature at the outlet of the combustion device is required to be obtained with an accuracy of ± 10 ° C.

はじめに、ステップ１として、観測点および観測変数、ならびに評価点および評価変数を設定する。具体的には、燃焼装置の入口境界の各領域Ｆ_１、Ｆ_２、Ｆ_３、Ｆ_４の直上であって、かつ燃焼装置内のガス温度およびガス組成の計測が実際に可能な箇所に、４つの観測点Ｒ_１、Ｒ_２、Ｒ_３、Ｒ_４を定義し、それらの位置における燃焼装置内のガス温度およびガス組成（未燃ガス濃度、酸素濃度、水分濃度の３つ）を観測変数として定義する。
また、燃焼装置の出口境界に１つの評価点Ｒ_５を定義し、その位置におけるガス温度およびガス組成（同じく未燃ガス濃度、酸素濃度、水分濃度の３つ）を評価変数として定義する。 First, as step 1, observation points and observation variables, and evaluation points and evaluation variables are set. Specifically, in each region of the inlet boundary of the combustion device, _{directly above F 1} , F ₂ , F ₃ , and F ₄ , and in a place where the gas temperature and gas composition can be actually measured in the combustion device. Four observation points R ₁ , R ₂ , R ₃ , and R ₄ are defined, and the gas temperature and gas composition (unburned gas concentration, oxygen concentration, and water concentration) in the combustion device at those positions are observed variables. Defined as.
Also, define one evaluation point R ₅ at the exit boundary of the combustion device, defined gas temperature and gas composition (also unburned gas concentration, oxygen concentration, three water concentration) at its position as an evaluation variable.

つまり、１つの観測点ないし評価点において、４つずつの観測変数ないし評価変数が定義される。以下では本実施例における観測変数Ｙ_ｊ１（ｊ_１＝１〜１６）および評価変数Ｙ_ｊ２（ｊ_２＝１７〜２０）を、それぞれ次の表１のように定義する。 That is, four observation variables or evaluation variables are defined at one observation point or evaluation point. In the following, the observation variable Y _j1 (j _{1 = 1 to} 16) and the evaluation variable Y _j2 (j ₂ = 17 to 20) in this embodiment are defined as shown in Table 1 below.

次に、ステップ２として、解析パラメータＸ_ｉ（ｉ＝１〜２０）を変化させた多数の数値シミュレーションを実行する。具体的には、ステップ３で各観測変数Ｙ_ｊ１（ｊ_１＝１〜１６）および評価変数Ｙ_ｊ２（ｊ_２＝１７〜２０）を各解析パラメータＸ_ｉの２次多変数多項式として表現することで応答曲面の定式化を行うことを前提に、２次多変数多項式の係数の数（２３１個）を上回るように合計２５０点（ｋ＝２５０）の数値シミュレーションを実行した。 Next, as step 2, _{a large number of numerical simulations in which the analysis parameters X i} (i = 1 to 20) are changed are executed. Specifically, it each observation in Step 3 variables _Y j1 _{to (j} 1 = _1~16) and evaluated variables _Y j2 the _(j 2 = 17 to 20) expressed as a quadratic multivariate polynomial of each analysis parameter _{X i} On the premise that the response surface was formulated in, a total of 250 points (k = 250) were numerically simulated so as to exceed the number of coefficients (231) of the quadratic multivariate polynomial.

次に、ステップ３として、ステップ２で実行した数値シミュレーションの実行結果に基づいて、解析パラメータＸ_ｉ（ｉ＝１〜２０）と、観測変数Ｙ_ｊ１（ｊ_１＝１〜１６）および評価変数Ｙ_ｊ２（ｊ_２＝１７〜２０）との関係性を定式化する。ここでは、２次の多変数多項式を用い、さらにラッソ回帰によるスパース推定手法を適用することで、説明変数の選択を行った。変数選択の結果の詳細は省略するが、もとの多変量多項式の係数の数（２３１個）に対して、選択された変数の数は各観測変数および評価変数に対して、１０個から９０個程度となり、効果的に変数選択が行われたことが分かる。 Next, as step 3, based on the execution result of the numerical simulation executed in step 2, the analysis parameter X _i (i = 1 to 20), the observation variable Y _j1 (j _{1 = 1 to} 16), and the evaluation variable Y _{Formulate the relationship with j2} (j ₂ = 17 to 20). Here, the explanatory variables were selected by using a quadratic multivariate polynomial and further applying a sparse estimation method by lasso regression. Although the details of the result of variable selection are omitted, the number of selected variables is 10 to 90 for each observed variable and evaluation variable for the number of coefficients (231) of the original multivariate polynomial. It can be seen that the variables were effectively selected.

次に、ステップ４として、上記により定式化された解析パラメータＸ_ｉ（ｉ＝１〜２０）と観測変数Ｙ_ｊ１（ｊ_１＝１〜１６）および評価変数Ｙ_ｊ２（ｊ_２＝１７〜２０）の回帰式の精度を評価する。具体的には、ステップ２において実行した数値シミュレーションの全ての計算条件Ｘ_ｉ ^ｋ（ｉ＝１〜２０、ｋ＝１〜２５０）について、ステップ３で得られた回帰式に各Ｘ_ｉ ^ｋを代入することにより得られる予測値Ｙ^〜 _ｊ１ ^ｋおよびＹ^〜 _ｊ２ ^ｋの値と、シミュレーション結果として得られたＹ_ｊ１ ^ｋおよびＹ_ｊ２ ^ｋの値との二乗平均誤差を比較し、十分な精度が得られていることを確認した。より高い精度が必要である場合には、図１に示したフローチャートに従って、適宜数値シミュレーションの実行回数な
いし多変量多項式の次数を増やした上で、再度ステップ２からステップ４までの手順を実行すればよい。 Next, as step 4, the analysis parameter X _i (i = 1 to 20), the observation variable Y _j1 (j _{1 = 1 to} 16) and the evaluation variable Y _j2 (j ₂ = 17 to 20) formulated as described above. Evaluate the accuracy of the regression equation of. Specifically, all the calculation conditions of numerical simulation performed in step ^{_{2 X i k (i = 1~20}} , k = 1~250) for, substituting each _X ^{i k} the obtained regression equation in Step 3 Sufficient accuracy was obtained by comparing the squared average error between the predicted values Y ^to _j1 ^k and Y ^to _j2 ^k obtained by the above and the values of Y _j1 ^k and Y _j2 ^{k obtained as simulation results.} I confirmed that. If higher accuracy is required, the number of times the numerical simulation is executed or the degree of the multivariate polynomial is increased as appropriate according to the flowchart shown in FIG. 1, and then the steps 2 to 4 are executed again. good.

次に、ステップ５として、観測点における観測変数Ｙ_ｊ１の観測値Ｙ^＊ _ｊ１および評価点における評価変数Ｙ_ｊ２の観測値Ｙ^＊ _ｊ２を実験的に得た。 Next, as step 5, the observed value Y ^* _{j1 of the observation variable Y j1 at} the observation point and the observation value Y ^* _j2 of the evaluation variable Y _{j2 at} the evaluation point were experimentally obtained.

次に、ステップ６として、ステップ５で得られた観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２の統計情報と、ステップ３で定式化された解析パラメータＸ_ｉと観測変数Ｙ_ｊ１、および解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性とに基づいて、ベイズ推定を行うことにより、観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２が得られたもとでの、解析パラメータＸ_ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）＝Ｐ（Ｘ_ｉ｜Ｙ^＊ _ｊ）をベイズ推定手法（ハミルトニアンモンテカルロ法）により求めた。 Next, as step 6, the statistical information of ^{the observed values Y *} _j1 and Y ^* _{j2 obtained in step 5, the analysis parameter X i} and the observation variable Y _j1 formulated in step 3, and the analysis parameter X _i The posterior probability distribution P (X ^* _i ) of the analysis parameter X _i ^{under the observation values Y *} _j1 and Y ^* _j2 obtained by Bayesian estimation based on the relationship of the evaluation variable Y _{j2 =} P (X _i | Y ^* _j ) was obtained by the Bayesian estimation method (Hamiltonian Monte Carlo method).

その結果の一例として、解析パラメータＸ_１からＸ_１０までの事後確率分布Ｐ（Ｘ^＊ _ｉ）の推定結果を図３に示す。これらの図に示されるように、実験的に得られる観測値の統計情報に対応する解析パラメータの値の事後確率分布は、一般に幅広い分布を示すことが分かる。 As an example of the result, FIG. 3 shows the estimation result of the posterior probability distribution P (X ^* _{i ) from the analysis parameters X 1} to X _10. As shown in these figures, it can be seen that the posterior probability distribution of the values of the analysis parameters corresponding to the statistical information of the observed values obtained experimentally generally shows a wide distribution.

次に、ステップ７として、ステップ６で得られた解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）を、ステップ３で定式化された解析パラメータＸ_ｉと評価変数Ｙ_ｊ２の関係性に入力することにより、評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）＝Ｐ（Ｙ_ｊ２｜Ｘ^＊ _ｉ）＝Ｐ（Ｙ_ｊ２｜Ｙ^＊ _ｊ）を得た。その結果の一例として、各観測点および評価点におけるガス温度の事後確率分布Ｐ（Ｙ^＊＊ _１）、Ｐ（Ｙ^＊＊ _５）、Ｐ（Ｙ^＊＊ _９）、Ｐ（Ｙ^＊＊ _１３）、Ｐ（Ｙ^＊＊ _１７）を図４に示す。 Next, as step 7, the posterior probability distribution P (X ^* _i ) of the estimated value X ^* _{i of the analysis parameter X i} obtained in step 6 is subjected to the analysis parameter X _i and the evaluation variable Y formulated in step 3. by inputting the relationship _j2, the posterior probability distribution ^P of the evaluation variables ^{_{Y j2 (Y ** j2) =}} P was obtained _{^{| | (Y * j Y j2}} ) (Y j2 X * i) = P. As an example of the result, posterior probability distributions of gas temperature at each observation point and evaluation point P (Y ^** ₁ ), P (Y ^** ₅ ), P (Y ^** ₉ ), P (Y ^** ₁₃ ) , P (Y ^** ₁₇ ) are shown in FIG.

前述のように、本手法により同定された未知の解析パラメータＸ_ｉ（ｉ＝１〜２０）の事後確率分布Ｐ（Ｘ^＊ _ｉ）は、図３に示すように幅広い分布を示すものの、各観測点および評価点におけるガス温度の事後確率分布は、図４に示すように観測変数および評価変数の実際の観測値（図中の太線）を含むシャープな分布となっており、それらの実際の観測値Ｙ^＊ _ｊ１およびＹ^＊ _ｊ２をうまく説明できていることが分かる。 As described above, the posterior probability distribution P (X ^* _i ) of the unknown analysis parameter X _i (i = 1 to 20) identified by this method shows a wide distribution as shown in FIG. 3, but each observation. As shown in Fig. 4, the posterior probability distribution of the gas temperature at the points and evaluation points is a sharp distribution including the observation variables and the actual observation values (thick lines in the figure) of the observation variables, and those actual observations. It can be seen that the values Y ^* _j1 and Y ^* _{j2 are well explained.}

次に、ステップ８として、ステップ７で得られた評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間の評価を行う。前述したように本実施例では、評価変数である燃焼装置出口のガス温度を±１０℃の精度で求めることが要求されていたが、図４で示される燃焼装置出口のガス温度Ｙ_１７の事後確率分布は、真値（実際の測定値）８３６℃に対して、約８３４〜８４１℃の範囲にシャープに分布しており、十分な精度が得られていることが分かる。より高い精度が必要である場合には、図１に示したフローチャートに従って、適宜観測変数の数を増やした上で、再度ステップ１からステップ８までの手順を実行すればよい。 Next, as step 8, the credit interval of the posterior probability distribution P (Y ^** _{j2 ) of the evaluation variable Y j2 obtained in step 7 is evaluated.} As described above, in this embodiment, it was required to obtain the gas temperature at the outlet of the combustion device, which is an evaluation variable, with an accuracy of ± 10 ° C., but after _{the gas temperature Y 17 at the outlet of the combustion device shown in FIG.} The probability distribution is sharply distributed in the range of about 834 to 841 ° C. with respect to the true value (actual measured value) of 836 ° C., and it can be seen that sufficient accuracy is obtained. When higher accuracy is required, the steps from step 1 to step 8 may be executed again after increasing the number of observation variables as appropriate according to the flowchart shown in FIG.

以上の実施例では、本発明を燃焼装置における未知のガス組成の同定に適用しているが、本発明はこの例に限定されず、その他の様々な工業的装置あるいは機械における未知の解析パラメータの推定、例えば機器や装置における温度分布や濃度分布の同定、境界表面力の同定、電流密度の同定などにも用いることができる。 In the above examples, the present invention is applied to the identification of an unknown gas composition in a combustion device, but the present invention is not limited to this example, and the unknown analytical parameters in various other industrial devices or machines are used. It can also be used for estimation, for example, identification of temperature distribution and concentration distribution in equipment and devices, identification of boundary surface force, identification of current density, and the like.

以上、本発明の実施形態および実施例についてその詳細を説明したが、本発明は上記の実施形態および実施例に限定されるものではなく、特許請求の範囲、及び明細書と図面に記載された技術的思想を逸脱しない範囲において種々変更又は修正を行って実施することが可能である。 Although the details of the embodiments and examples of the present invention have been described above, the present invention is not limited to the above-described embodiments and examples, and is described in the claims, the specification and the drawings. It is possible to make various changes or modifications within the range that does not deviate from the technical idea.

Claims (6)

物理現象を模擬するための数値シミュレーションを所望の解析精度において実行するために必要な、１つ以上の未知の解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）の推定値Ｘ^＊ _ｉを推定する方法であって、
対象とする物理現象における、１つ以上の観測点における１つ以上の観測変数Ｙ_ｊ（ｊ＝１，２，・・・，Ｍ）について、前記１つ以上の解析パラメータＸ_ｉの値を変化させた複数の数値シミュレーションを実行するステップと、
前記複数の数値シミュレーションを実行して得られた結果に基づいて、前記解析パラメータＸ_ｉと前記観測変数Ｙ_ｊとの関係性を定式化するステップと、
前記物理現象について実験的計測を行うことにより、前記観測変数Ｙ_ｊの観測値Ｙ^＊ _ｊ（ｊ＝１，２，・・・，Ｍ）を得るステップと、
前記観測変数Ｙ_ｊの観測値Ｙ^＊ _ｊの統計情報と、前記定式化された前記解析パラメータＸ_ｉと前記観測変数Ｙ_ｊの関係性とに基づいて、ベイズ推定を行うことにより、前記観測値Ｙ^＊ _ｊが得られたもとでの、前記解析パラメータＸ_ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）＝Ｐ（Ｘ_ｉ｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）として、前記解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの確率分布を得るステップと、
を含むことを特徴とする、解析パラメータの推定方法。 ^{Estimated value X * of} _{one or more unknown analysis parameters X i} (i = 1, 2, ..., N) required to perform a numerical simulation for simulating a physical phenomenon with desired analysis accuracy. _It is a method of estimating i
_{For one or more observation variables Y j} (j = 1, 2, ..., M) at one or more observation points in the target physical phenomenon _{, the values of the one or more analysis parameters X i} are changed. Steps to execute multiple numerical simulations
A step of, based on the result obtained by executing a plurality of numerical simulation, to formulate the relationship between the analysis parameters X _i and the observed variable Y _j,
The step of obtaining the observed value Y ^* _j (j = 1, 2, ..., M) of the _{observed variable Y j} by performing experimental measurement on the physical phenomenon, and
Wherein the observed variables Y _j observations Y ^* _j statistics, based on the a formalized the analysis parameters X _i and the relationship of the observed variables Y _j, by performing Bayesian estimation, the observation value _Given that Y ^* j is obtained, the posterior probability distribution P (X ^* _i ) = P (X _i | Y ^* _1, Y ^* _2, ... _, Y ^* _{M ) of the analysis parameter X i} is set as described above. Steps to obtain the probability distribution of the estimated value X ^* _i of the analysis parameter X _i,
A method of estimating analysis parameters, which comprises. 前記解析パラメータＸ_ｉと前記観測変数Ｙ_ｊとの関係性は、前記観測変数Ｙ_ｊ（ｊ＝１，２，・・・，Ｍ）を目的変数とし、前記解析パラメータＸ_ｉ（ｉ＝１，２，・・・，Ｎ）を説明変数とした多変数多項式回帰により定式化されることを特徴とする、請求項１に記載の解析パラメータの推定方法。 The relationship between the analysis parameter X _i and the observation variable Y _j is determined by using the observation variable Y _j (j = 1, 2, ..., M) as the objective variable and the analysis parameter X _i (i = 1,). The method for estimating analysis parameters according to claim 1, wherein the method is formulated by a multivariate polynomial regression using 2, ..., N) as an explanatory variable. 前記多変数多項式回帰において、リッジ回帰、ラッソ回帰またはエラスティックネットのいずれかからなるスパース推定手法を用いて前記説明変数の選択を行うことを特徴とする、請求項２に記載の解析パラメータの推定方法。 The analysis parameter estimation according to claim 2, wherein in the multivariate polynomial regression, the explanatory variables are selected by using a sparse estimation method including ridge regression, lasso regression, or elastic net. Method. 前記ベイズ推定において、マルコフ連鎖モンテカルロ法またはハミルトニアンモンテカルロ法を用いることを特徴とする、請求項１から３のいずれか１項に記載の解析パラメータの推定方法。 The method for estimating analysis parameters according to any one of claims 1 to 3, wherein the Markov chain Monte Carlo method or the Hamiltonian Monte Carlo method is used in the Bayesian estimation. 前記１つ以上の観測点は、前記解析パラメータＸ_ｉとの相関性の高い観測変数Ｙ_ｊ１（ｊ_１＝１，２，・・・，Ｍ_１）が得られるように選ばれた観測点と、前記数値シミュレーションにおいて計算結果の評価に用いられる評価変数Ｙ_ｊ２（ｊ_２＝Ｍ_１＋１，Ｍ_１＋２，・・・，Ｍ）が得られるように選ばれた観測点とからなることを特徴とする、請求項１から４のいずれか１項に記載の解析パラメータの推定方法。 Wherein the one or more observation points, the analysis parameters _{X i} and higher observed variables correlated of _{Y j1 (j 1 = 1,2,} ···, M 1) and observation point chosen to obtain _{, The observation point is selected so that the evaluation variable Y j2} (j ₂ = M ₁ + 1, M ₁ + 2, ..., M) used for evaluating the calculation result in the numerical simulation can be obtained. The method for estimating analysis parameters according to any one of claims 1 to 4. 前記解析パラメータＸ_ｉの推定値Ｘ^＊ _ｉの事後確率分布Ｐ（Ｘ^＊ _ｉ）を、前記定式化された前記解析パラメータＸ_ｉと前記観測変数Ｙ_ｊとの関係性に入力することにより、前記評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）＝Ｐ（Ｙ_ｊ２｜Ｘ^＊ _ｉ）＝Ｐ（Ｙ_ｊ２｜Ｙ^＊ _１，Ｙ^＊ _２，・・・_，Ｙ^＊ _Ｍ）を得るステップと、
前記評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間が、いずれかのｊ_２について所定の閾値を上回る場合に、前記観測変数Ｙ_ｊ１の数Ｍ_１を増やすことにより、前記評価変数Ｙ_ｊ２の事後確率分布Ｐ（Ｙ^＊＊ _ｊ２）の信用区間が、いずれのｊ_２についても前記所定の閾値を下回るようにするステップと、
をさらに含むことを特徴とする、請求項５に記載の解析パラメータの推定方法。
By inputting the relationship between the analysis parameters X _i estimate X ^* _i the posterior probability distribution P (X * _ⁱ⁾ of the analysis the was formulated parameter X _i and the observed variable Y _j of the Step to obtain posterior probability distribution P (Y ^** _j2 ) = P (Y _j2 | X ^* _i ) = P (Y _j2 | Y ^* _1, Y ^* _2, ... _, Y ^* _M ) of evaluation variable Y _j2 When,
When the credit interval of the posterior probability distribution P (Y ^** _j2 ) of the evaluation variable Y _j2 exceeds a predetermined threshold for _{any j 2} , the number M _{1 of the observation variable Y j1} is increased to obtain the above. A step of making the credit interval of the posterior probability distribution P (Y ^** _j2 ) of the evaluation variable Y _j2 fall below the predetermined threshold for _{any j 2.}
5. The method for estimating analysis parameters according to claim 5, further comprising.

Images