JP3881510B2 - 曲面モデルの同定方法及びプログラム - Google Patents
曲面モデルの同定方法及びプログラム Download PDFInfo
- Publication number
- JP3881510B2 JP3881510B2 JP2000378625A JP2000378625A JP3881510B2 JP 3881510 B2 JP3881510 B2 JP 3881510B2 JP 2000378625 A JP2000378625 A JP 2000378625A JP 2000378625 A JP2000378625 A JP 2000378625A JP 3881510 B2 JP3881510 B2 JP 3881510B2
- Authority
- JP
- Japan
- Prior art keywords
- data
- identification
- calculated
- distance
- calculation step
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Landscapes
- Complex Calculations (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、製剤設計、材料設計、システムの運転条件・製品製造条件の設計／調整等に使用される応答曲面法などに利用可能な、曲面モデルの同定方法及びプログラムに関するものである。
【０００２】
【従来の技術】
従来、実験計画法に基づいて収集されたデータを用いて、効率よく製品の製造条件・システムの運転条件等を設計・調整するために、応答曲面法（Response Surface Methodology；RSM ）が用いられている。以下、従来の応答曲面法（文献１「轟 章、応答曲面法、平11-12月日本機械学会講習会「応答曲面法による非線形問題の最適設計入門」資料、URL：http://ueno.mes.titech.ac.jp/rec-res.html にて公開中」、文献２「長谷川ほか、実数型交差モデルと応答曲面モデルを用いた逐次的近似最適化に関する一試行、日本計算工学会論文集、２０００年５月２４日発行」、文献３「Takayama、Artificial Neural Network as a Novel Method to Optimize Pharmaceutical Formulations、Pharmaceutical Research, Vol.16, No.1, 1999」など）について説明する。
【０００３】
図１１は応答曲面法における処理の流れを示すフローチャート図である。応答曲面法には、実験計画法に基づき必要最低限のデータを収集する処理（ステップ４０１）と、収集されたデータをもとに各要因（設計条件）とその応答（最適にしたい特性値）を結びつける近似関数（曲面モデル）を同定する処理（ステップ４０２）と、得られた近似関数（曲面モデル）を利用して、所望の特性値となるような最適設計条件を各種最適化アルゴリズムを利用して求める処理（ステップ４０３）の３つの処理ステップがある。
【０００４】
次に、上記ステップ４０１〜４０３において、従来利用されている技術を以下に説明する。ステップ４０１では、実験計画法に基づき、実験計画データを収集する。これにより必要最低限なデータを効率よく収集する。利用される実験計画法の要因配置（実験点配置）としては、完備要因計画（全因子計画；fully factorial design）法、直交計画（orthogonal design ）法、球形中心２次複合計画（中央複合計画；central composite design）法などが代表的である（文献１、文献４「橋田 充 編、経口投与製剤の処方設計（第7章 処方・製法の最適化）、薬業時報社、１９９８」）。
【０００５】
ステップ４０２では、ステップ４０１で収集したデータを使って、各要因（設計条件）ｘi （ｉ＝１，２，・・・，ｋ）とその応答（最適にしたい特性値）ｙとを結びつける近似関数（曲面モデル）ｙ＝ｆ（ｘ1，ｘ2，・・・，ｘk ）を以下のような手法を利用して求める。
【０００６】
（Ｉ）多項式を用いた手法（文献１、文献２）は、各要因（入力変数）と特性値（出力変数）とを基本的には線形の関数で近似する手法である。多くの場合は、以下の２次多項式が利用される。
【０００７】
【数１】
【０００８】
式（１８）において、ａ0，ｂi，ｃi，ｃij は定数である。非線形な対象に対応する場合には、相互作用項（式（１８）の右辺第４項）を導入したり、べき乗関数、指数関数あるいは対数関数などを利用して、入力変数の変換を行ったりする。前記２次多項式が、応答曲面法の近似関数としては最も頻繁に利用される手法である。
【０００９】
（ＩＩ）ニューラルネットワークを用いた手法（文献３）は、非線形な入出力関係を近似する手法である。入力を各要因（設計条件）、出力を特性値として、階層型のニューラルネットワークにバックプロパゲーション法を利用して学習させる。入力層から中間層、出力層に至る信号は、前層からの信号を次の式（１９）により統合し、式（２０）のようなシグモイド関数に代表される伝達関数により、次の層へと出力される。
【００１０】
【数２】
【００１１】
【数３】
【００１２】
なお、ｗpqは次層ユニットｑと前層ユニットｐ間の荷重、ｘp は前層からの出力であり、ｆ（ｙq ）は次層への出力として伝達される。また、ａはシグモイド関数の勾配係数である。
【００１３】
（ＩＩＩ）スプライン関数を用いた手法（文献５「桜井 明著、スプライン関数入門、東京電機大学出版局、１９８１」、文献６「桜井 明監修、パソコンによるスプライン関数、東京電機大学出版局、１９８８」）は、非線形な特性であっても関数化することができ、滑らかな曲線/曲面が得られる手法である。通常広く認知されているスプラインは、文献「マグローヒル科学技術用語大辞典」に記載されている通り、導関数に基づくスプライン関数のことである。この導関数に基づくスプライン関数の数学的定義は、文献４に記載されている通り、「ｍ次のスプライン関数は、そのｍ階微分が階段関数で、ｍ−１階以下の微分が連続であるような関数である」というものである。また、ｎ個の接点をもつｍ次のスプライン関数は以下のようになる。
【００１４】
【数４】
【００１５】
なお、（ｘ−ｑi）m +は切断べき関数と呼ばれる関数であり、以下のように表される。
【００１６】
【数５】
【００１７】
（ＩＩＩＩ）多変数適応回帰スプライン（MARS）を用いた手法（文献７「T.J.Hastie and R.J.Tibshirani、Nonparametric Regression and Classification Part I - Nonparametric Regression, NATO ASI Series. F. Computer and System Sciences, Vol.136, pp.62-69, 1994」）は、多項式の回帰近似を区分に分割して行う手法である。分割する領域の大きさは計算により自動的に決定される。
【００１８】
最後に、ステップ４０３では、得られた曲面モデルを利用して、公知の各種最適化アルゴリズムにより、所望の最適値探索を行い、最適設計条件を見出す。連続量に対する公知の最適化アルゴリズムとしては、シンプレックス法、勾配法、共役勾配法、ニュートン法、準ニュートン法、シミュレーテッド・アニーリング法、トンネリング・アルゴリズム、遺伝的アルゴリズム（GA）などが代表的である（文献８「茨木俊秀ほか、最適化の手法、共立出版、１９９３」、文献９「相吉 英太郎、「最適化」の展望、計測と制御、Vol.29,No.12、１９９０」）。
【００１９】
また、制約条件下での最適化手法として、ペナルティ関数を利用した方法（文献３、文献８）などもある。さらに、複数の曲面モデルを利用することで、複数特性の最適値（妥協解）を見出す多目的最適化（文献３、文献１０「西川ほか、岩波講座 情報科学１９「最適化」、岩波書店、１９８２」）も実施することができる。多目的最適解を求める方法としては、重みパラメータ法、ε制約法、目標計画法、対話型の方法などがある。
【００２０】
【発明が解決しようとする課題】
応答曲面法は、数少ないデータ（必要最低限のデータ）から、効率よく最適条件を見出すために利用される。しかしながら、従来の応答曲面法では、前述のステップ４０２における近似関数（曲面モデル）の同定に労力を要することがあり、実際には効率が悪い可能性があるという問題点があった。例えば、前述の従来の近似関数を利用した場合には、以下のような問題がある。
【００２１】
（Ｉ）多項式を用いた手法は、基本的には、入出力関係が線形な対象をモデル化する手法であるため、各要因と特性値が非線形な関係の場合（現実にはこのケースが多い）は、効率よくモデル化することが困難となる場合がある。非線形な対象をモデル化する手段として、相互作用項を追加したり、べき乗関数、指数関数あるいは対数関数などを導入して近似することもあるが、必要な精度のモデルを得るためには多大な労力を要することがある。例えば、相互作用項を導入する場合には、すべての要因の相互作用項を追加すべきか、あるいは一部の相互作用項とするか、一部の相互作用項とする場合には、どの要因の組み合わせを追加すればよいかということは試してみないと分からない。また、どのような関数を導入すればよいかについても試してみないと分からないことがある。つまり、近似関数を得るためには、試行錯誤的となり効率が悪い。
【００２２】
さらに、非線形性の強い対象の場合は、全領域を１つの関数で近似することが困難な場合がある。そのため、小さな領域に分けてそれぞれに異なる多項式近似関数をあてはめることもある。しかし、対象の特性が未知なものについて、どの領域がどのような関数に相当するかを決めることは困難である。そのため，結果的に所望の精度のモデルを得ることができない。また、２次多項式による近似では、現実的にはあり得ない値（特性値の限界値以上あるいは限界値以下の値）の応答が現れる場合があり、所望の応答曲面にならない場合がある。この場合、最適な領域の分割や関数の選択にさらに十分な時間をかければ、所望のモデルが得られる可能性もあるが、非常に多大な労力が必要となり、モデル同定の効率が悪化する。
【００２３】
（ＩＩ）ニューラルネットワークを用いた手法では、ユーザが指定すべきパラメータが多すぎ、また、パラメータを決定するための指標がないので、ネットワークモデルを試行錯誤的に求めざるを得なくなり、実際には効率よく所望の近似関数を求めることは困難である。例えば、ユーザが予め設定する必要があるパラメータとして代表的なものには、ネットワークの階層数、中間層のユニット数、伝達関数の関数式（伝達関数をシグモイド関数とする場合には、その関数の傾き）、適用する学習法則、学習停止条件、学習スタート時の荷重係数の初期値などがある。
【００２４】
このようなパラメータを変化させると、応答曲面が大きく変化するが、パラメータをどう変化させればよりよいモデルが得られるかを知るための明確な指標がないため、結果的には試行錯誤的にモデルを作らざるを得ない。さらに、学習計算は繰り返し計算（収束計算）となるため、計算時間が必要となることがある。ニューラルネットワークで近似した応答曲面では、応答曲面が伝達関数であるシグモイド関数の形状となってしまい、所望の応答曲面にならない場合がある。この場合、よりよい応答曲面にするにはさらに時間がかかり、モデル同定の効率が悪化する。
【００２５】
（ＩＩＩ）通常、工学分野で広く認知・利用されている導関数に基づくスプライン関数では、１要因（１入力変数）に対する出力値の補間曲線（近似関数）を描くことは容易にできるが、２要因（２入力変数）以上の場合には、以下のような制限が発生する。
（α）補間計算（近似関数の同定）をするために与えられるデータが、入力空間内で格子状となるデータであり、各座標点に順序をつける必要がある（文献１１「William H. Pressほか、Numerical Recipes in C、技術評論社、１９９３」）。
（β）与えられたデータに囲まれる区分領域（２入力変数の場合は、３角領域）を切り出す必要がある（文献１２「曲渕 英邦、スプライン関数を用いたランダム布置標本データからの２次元補間、日本建築学会大会学術講演梗概集（九州）１９９８年９月」）。
【００２６】
実験データの近似を目的とするならば、（α）の格子状データが得られることは期待できないため、不規則分布のデータを前提とした（β）を適用せざるを得ない。しかし、２入力変数以上の場合、補間計算以外に領域切り出しのアルゴリズム（例えば、Delaunay三角分割；文献１３「MATLAB ver5.3 オンラインマニュアル、Function Reference "delaunay" 」）が必要であり、また各次元ごとに（入力変数の数に応じて）領域切り出しのアルゴリズムが必要である（２入力変数の場合には例えば３角形分割、３入力変数の場合には例えば４面体領域切り出しなど）。つまり、アルゴリズムが高次元になればなるほど、アルゴリズムが複雑となり、プログラムも大きくなり、計算時間も必要となるため、現実的には３入力変数以上の補間にはほとんど利用されていない。
【００２７】
（ＩＩＩＩ）多変数適応回帰スプライン（MARS）を用いた手法は、基本的には多項式近似と同様なアルゴリズムであるため、多項式と同様な問題がある。ただし、区分領域での近似を採用していることから、非線形近似にもある程度対応可能である。また、近似の基本となる関数もあらかじめ候補をいくつか用意しておけば、そのなかで最適化なものを自動的に選択することができる。しかし、近似の基本となる関数の候補として何を用意すればよいかわからない場合は、多数の関数候補を用意する必要があり、これによって計算負荷（計算時間、要求メモリサイズ）が増大する。また、非線形性が強い場合などには、より多数の区分領域に分けて回帰近似する必要があり、これによっても計算負荷が大きくなるため、実用上問題である。また、領域を区分して多項式近似をするため、実験データを細分して利用することとなる。このため、全領域において、実験計画法に基づき収集したデータでは、多項式近似に必要なデータ量が得られないことがある。また、細分した領域ごとに実験計画法を適用しデータを収集することは、実験コストの増大や効率の悪化となるため、実用上問題である。
【００２８】
なお、他の補間手法（近似関数）として、ラグランジュ補間多項式もあるが、この手法を利用した場合には、データ点数によっては、ルンゲの現象という振動現象が現れる。これにより、なめらかではない不自然な、応答の激しい曲面となってしまうので（文献１４「桜井 明監修、パソコンによるスプライン関数、東京電機大学出版局、１９８８、pp.4-6」）、実用上問題である。
【００２９】
以上のように、従来の応答曲面法では、モデル同定が試行錯誤的となり、モデル同定に多大な時間がかかるという問題点があった。また、従来の応答曲面法では、モデル同定を誤る可能性があった。
本発明は、上記課題を解決するためになされたもので、高精度な曲面モデルの同定を効率良く行うことができる曲面モデルの同定方法及びプログラムを提供することを目的とする。
【００３０】
【課題を解決するための手段】
本発明は、ｎ（ｎは２以上の整数）種類の変量間の関係を近似したｎ次元の曲面モデルをコンピュータを使用して同定する曲面モデルの同定方法であって、予めデータが収集された前記ｎ種類の変量のうちｎ−１種類の変量で定義されるｎ−１次元の空間において各データ間の距離を計算して記憶装置に記憶させる第１の同定時距離計算ステップ（３０２）と、この第１の同定時距離計算ステップで計算された各距離に対応するグリーン関数値を計算して記憶装置に記憶させる第１の同定時グリーン関数計算ステップ（３０３）と、この第１の同定時グリーン関数計算ステップで計算されたグリーン関数値から前記収集されたデータが前記曲面モデルの同定にとって有効なデータか否かを示す評価指標値を計算して記憶装置に記憶させる評価指標値計算ステップ（３０４）と、この評価指標値が最大となるようなデータを前記収集されたデータ中から選択するデータ選択ステップ（３０５，３０６）と、この選択されたデータについて、前記ｎ種類の変量のうちｎ−１種類の変量で定義されるｎ−１次元の空間において各データ間の距離を計算して記憶装置に記憶させる第２の同定時距離計算ステップ（１０２）と、この第２の同定時距離計算ステップで計算された各距離に対応するグリーン関数値を計算して記憶装置に記憶させる第２の同定時グリーン関数計算ステップ（１０３）と、前記選択されたデータのうち前記ｎ−１種類を除く１種類のデータと前記第２の同定時グリーン関数計算ステップで計算されたグリーン関数値とから前記ｎ次元の曲面モデルを計算して記憶装置に記憶させるモデル同定ステップ（１０４）とからなるものである。このように、本発明は、従来の応答曲面法（実験計画法、近似関数、既存最適化アルゴリズムを組み合わせて、少ない実験データから効率よく最適解を求める手法）の近似関数として、グリーン関数に基づく重調和スプライン補間を適用している。重調和スプライン補間は、与えられたデータのみを使って非線形な特性であっても応答曲面を求めることができ、ユーザは近似関数の次数、構造などの各種パラメータの指定をする必要がない。また、重調和スプライン補間は、多次元への拡張が容易であり、多要因のモデルにも柔軟に対応することができる。
【００３１】
また、本発明の曲面モデルの同定方法の１構成例において、前記ｎ種類の変量のデータは、実験計画法によって収集されるものである。
【００３２】
また、本発明のプログラムは、予めデータが収集された前記ｎ種類の変量のうちｎ−１種類の変量で定義されるｎ−１次元の空間において各データ間の距離を計算して記憶装置に記憶させる第１の同定時距離計算ステップと、この第１の同定時距離計算ステップで計算された各距離に対応するグリーン関数値を計算して記憶装置に記憶させる第１の同定時グリーン関数計算ステップと、この第１の同定時グリーン関数計算ステップで計算されたグリーン関数値から前記収集されたデータが前記曲面モデルの同定にとって有効なデータか否かを示す評価指標値を計算して記憶装置に記憶させる評価指標値計算ステップと、この評価指標値が最大となるようなデータを前記収集されたデータ中から選択するデータ選択ステップと、この選択されたデータについて、前記ｎ種類の変量のうちｎ−１種類の変量で定義されるｎ−１次元の空間において各データ間の距離を計算して記憶装置に記憶させる第２の同定時距離計算ステップと、この第２の同定時距離計算ステップで計算された各距離に対応するグリーン関数値を計算して記憶装置に記憶させる第２の同定時グリーン関数計算ステップと、前記選択されたデータのうち前記ｎ−１種類を除く１種類のデータと前記第２の同定時グリーン関数計算ステップで計算されたグリーン関数値とから前記ｎ次元の曲面モデルを計算して記憶装置に記憶させるモデル同定ステップとをコンピュータに実行させるものである。
【００３３】
【発明の実施の形態】
［実施の形態の１］
本発明では、応答曲面法において効率よく最適条件を探索するため、効率よく近似関数（曲面モデル）を得ることを目的としており、衛星データから等高線を得るために利用されている、グリーン関数に基づく重調和スプラインを近似関数として適用している。
【００３４】
グリーン関数に基づく重調和スプラインについては、例えば文献１５「David T. Sandwell、Biharmonic Spline Interpolation of GEOS-3 and SEASAT Altimeter Data、Geophysical Research Letters, Vol.14, No.2, pp.139-142, 1987」、また文献１６「Paul Wessel, David Bercovici, Interpolation with Splines in Tension: A Green's Function Approach, Mathmatical Geology, Vol.30, No.1, 1998 」に記載されている。
【００３５】
重調和スプラインは、工学分野で広く認知・利用されている、導関数に基づくスプライン関数とは、理論的にまったく異なる補間方法である。しかし、補間曲面は、導関数に基づくスプライン関数と同様に滑らかな曲面であり、非線形な特性でもモデル化することができる。
【００３６】
従来工学の分野で利用されているスプライン関数は、導関数に基づく補間値計算アルゴリズムを有するスプライン関数である。これに対し、重調和スプラインは、積分定理（グリーン定理）に基づく補間値計算アルゴリズムとなっている。そのため、従来のスプラインとは異なり、不規則分布データであっても特別なアルゴリズムを用意することなく、補間計算可能である。さらに、この手法は、データを与えるだけで、データの分布のみから関数を計算するため、ユーザが予め設定する必要があるパラメータは存在しない。また、多次元の入力変数（多要因）でも２入力変数（２要因）以下の場合と同様なアルゴリズムで動作するため、多次元（多要因）な対象にも、柔軟に適用可能である。
【００３７】
以下、この重調和スプライン補間のアルゴリズムを記述する。文献１５によると、任意の入力変数ベクトルｘに対する重調和スプライン関数ｆ（ｘ）は次式のようになる。
【００３８】
【数６】
【００３９】
ただし、Ｎは予め収集されたデータｘj（ｊ＝１〜Ｎ） の数、ｚは任意のデータｘと予め収集されたデータｘj との距離である。また、φm（ｚ） はグリーン関数であり、距離ｚとデータｘ，ｘj の次元数（データｘ，ｘj を特徴付ける変量の数）ｍとに応じて、次式のように定められる。
【００４０】
【数７】
【００４１】
式（１）の係数ｃj は、データが与えられると線形マトリクスを解くことにより自動的に算出できる。これにより求めたい入力変数ベクトルｘi に応じた出力値（補間値）ｙi は、次式により算出できる。
【００４２】
【数８】
【００４３】
式（３）において、ｚijはデータｘi と予め収集されたデータｘj との距離である。なお、式（１）は、次式より導出される。
【００４４】
【数９】
【００４５】
ただし、∇4 は重調和演算子、δ（ｚ）はデルタ関数である。
以上のように、重調和スプラインでは、データが与えられると、入力空間におけるデータ間の距離と、各次元（入力変数の数）に応じたグリーン関数とから、自動的に関数値が計算できる。そのため、不規則分布データであっても、さらに多次元（多要因）データであっても、特別なアルゴリズムを用意する必要はなく、柔軟に多次元拡張可能である。また線形マトリクス計算であるため、計算負荷も大きくない。
【００４６】
次に、以上のような重調和スプライン補間のアルゴリズムを用いた曲面モデルの同定方法について詳細に説明する。図１は本発明の第１の実施の形態となる曲面モデルの同定装置の構成を示すブロック図、図２は図１の同定装置の動作を示すフローチャート図である。図２（ａ）は同定装置のモデル同定時の動作を示し、図２（ｂ）は出力データ計算時の動作を示している。
【００４７】
図１の同定装置は、例えば実験計画法によって予めデータが収集されたｎ（ｎは２以上の整数）種類の変量のうちｎ−１種類の変量で定義されるｎ−１次元の空間において各データ間の距離を計算する距離計算部１と、この距離計算部１で計算された各距離に対応するグリーン関数値を計算するグリーン関数計算部２と、前記収集されたデータのうち前記ｎ−１種類を除く１種類のデータと前記グリーン関数値とから前記ｎ種類の変量間の関係を近似したｎ次元の曲面モデルを求めるモデル係数計算部３と、グリーン関数値と同定されたｎ次元の曲面モデルとを用いて、特定の入力データに対応する１次元の出力データを補間計算する出力データ計算部４とを有している。
【００４８】
図１の同定装置は、コンピュータ上で実現することができる。コンピュータは、ＣＰＵ、ＲＯＭ、ＲＡＭ、ディスプレイ装置やキーボードあるいは外部記憶装置とのインタフェースをとるための回路などを備えた周知の構成のものでよい。ＣＰＵは、ＲＯＭ若しくはＲＡＭに記憶されたプログラム、又はキーボードから入力されたコマンドに従って処理を実行する。また、ＣＰＵは、外部記憶装置にデータを書き込んだり、外部記憶装置からデータを読み出したりすることができる。
【００４９】
本実施の形態の曲面モデルの同定方法を実現させるためのプログラムは、フレキシブルディスク、ＣＤ−ＲＯＭ、メモリカード等の記録媒体に記録された状態で提供される。この記録媒体を外部記憶装置に挿入すると、媒体に記録されたプログラムが読み取られ、コンピュータに転送される。そして、ＣＰＵは、読み込んだプログラムをＲＡＭ等に書き込む。こうして、同定装置は、以下で説明するような処理を実行する。
【００５０】
以下、このような同定装置のモデル同定時の動作を図２（ａ）を参照して説明する。まず、同定装置のユーザは、ｎ種類の変量によって特徴付けられるｎ次元データを例えば実験計画法によって複数個収集し、前記ｎ種類の変量をｎ−１種類と１種類の変量に分け、これに対応して前記ｎ次元データをｎ−１次元のデータと１次元のデータとに分けて同定装置に入力する。
【００５１】
一般的には、前記ｎ−１種類の変量はユーザが求める解の要因となる変量（例えば、曲面モデルの生成理由が製剤処方の最適化にあるとすれば、要因は各原材料の量）であり、前記１種類の変量は解となる変量（処方された製剤の特性値）である。本発明では、ｎ−１種類の変量が要因（入力）で、残りの１種類の変量が解（出力）である必要は必ずしもないが、ここでは便宜的に、ｎ−１種類の変量によって特徴付けられるｎ−１次元データを入力データ、残りの１種類の変量によって特徴付けられる１次元データを出力データと呼ぶ。
【００５２】
同定装置の距離計算部１は、ユーザによってデータが入力されると（ステップ１０１においてＹＥＳ）、ｎ−１次元の空間においてｎ−１次元の各入力データ間の距離ｚを入力データの全ての組み合わせについて計算して、図示しない記憶装置に格納する（ステップ１０２）。ここでの距離ｚは、ユークリッド距離である。
【００５３】
例えば、ｎ＝３として、入力の２変量をＸ1，Ｘ2、出力の１変量をＹとし、図３に示すように、２次元の入力データｘ1，ｘ2，ｘ3 の値をそれぞれ（Ｘ11，Ｘ21），（Ｘ12，Ｘ22），（Ｘ13，Ｘ23）とし、入力データｘ1，ｘ2，ｘ3 に対応する出力データをそれぞれｙ1，ｙ2，ｙ3 とすれば、入力データｘ1，ｘ2間の距離ｚ12（又はＺ21）、入力データｘ2，ｘ3間の距離ｚ23（又はＺ32）、入力データｘ3，ｘ1間の距離ｚ31（又はＺ13）はそれぞれ次式のように得られる。
【００５４】
ｚ12＝Ｚ21＝（（Ｘ11−Ｘ12）2＋（Ｘ21−Ｘ22）2）1/2 ・・・（５）
ｚ23＝Ｚ32＝（（Ｘ12−Ｘ13）2＋（Ｘ22−Ｘ23）2）1/2 ・・・（６）
ｚ31＝Ｚ13＝（（Ｘ13−Ｘ11）2＋（Ｘ23−Ｘ21）2）1/2 ・・・（７）
【００５５】
続いて、グリーン関数計算部２は、距離計算部１で計算された距離ｚ毎に対応するグリーン関数値φm（ｚ） を計算して、図示しない記憶装置に格納する（ステップ１０３）。前述のｎ＝３の例では、距離ｚ12，ｚ23，ｚ31に応じてグリーン関数値φm（ｚ12），φm（ｚ23），φm（ｚ31） を計算することになる。前述のとおり、グリーン関数値φm（ｚ） は式（２）によって計算することができる。なお、本実施の形態では、ｍ＝ｎ−１である。
【００５６】
モデル係数計算部３は、出力データｙとグリーン関数値φm（ｚ） とからｎ次元の曲面モデルのモデル係数ｃj （ｊ＝１〜Ｎ）を計算して、図示しない記憶装置に格納する（ステップ１０４）。モデル係数ｃj を求めるには、出力データｙとグリーン関数値φm（ｚ） とを式（３）に代入すればよい。前述のｎ＝３の例ではＮ＝３なので、以下のような連立方程式が得られる。
【００５７】
ｙ1 ＝ｃ2φm（ｚ12）＋ｃ3φm（ｚ13） ・・・（８）
ｙ2 ＝ｃ1φm（ｚ21）＋ｃ3φm（ｚ23） ・・・（９）
ｙ3 ＝ｃ1φm（ｚ31）＋ｃ2φm（ｚ32） ・・・（１０）
【００５８】
モデル係数計算部３は、式（８）〜式（１０）の連立方程式から最小二乗法により、モデル係数ｃ1，ｃ2，ｃ3 を計算する。
以上により、モデル同定が終了する。同定装置によって同定されたｎ＝３次元の曲面モデルは、式（１）とモデル係数ｃ1，ｃ2，ｃ3 とにより、次式のように表すことができる。
ｙ＝ｃ1φm（ｚi1）＋ｃ2φm（ｚi2）＋ｃ3φm（ｚi3） ・・・（１１）
【００５９】
次に、ｎ−１（＝ｍ）次元の空間内の特定の入力データｘi に対応する出力データ（補関値）ｙi を同定装置に計算させる場合の動作を図２（ｂ）を参照して説明する。
【００６０】
同定装置の距離計算部１は、ユーザによって入力データｘi が入力されると（図２ステップ２０１においてＹＥＳ）、ｎ−１次元の空間において入力データｘi とモデル同定時に入力された前記収集された各入力データｘ1，ｘ2，ｘ3 との間の距離ｚを図４に示すように入力データｘ1，ｘ2，ｘ3 毎に計算して、図示しない記憶装置に格納する（ステップ２０２）。
【００６１】
入力データｘiの値を（Ｘ1i，Ｘ2i）とすると、入力データｘi，ｘ１間の距離ｚi1、入力データｘi，ｘ2間の距離ｚi2、入力データｘi，ｘ3間の距離ｚi3はそれぞれ次式のように得られる。
ｚi1＝（（Ｘ1i−Ｘ11）2＋（Ｘ2i−Ｘ21）2）1/2 ・・・（１２）
ｚi2＝（（Ｘ1i−Ｘ12）2＋（Ｘ2i−Ｘ22）2）1/2 ・・・（１３）
ｚi3＝（（Ｘ1i−Ｘ13）2＋（Ｘ2i−Ｘ23）2）1/2 ・・・（１４）
【００６２】
続いて、グリーン関数計算部２は、距離計算部１で計算された距離ｚi1，ｚi2，ｚi3毎に対応するグリーン関数値φm（ｚi1），φm（ｚi2），φm（ｚi3） を計算して、図示しない記憶装置に格納する（ステップ２０３）。
【００６３】
出力データ計算部４は、特定の入力データｘi に対応する出力データ（補関値）ｙi を計算して、図示しない記憶装置に格納する（ステップ２０４）。出力データｙi を計算するには、グリーン関数計算部２で計算されたグリーン関数値φm（ｚi1），φm（ｚi2），φm（ｚi3） を以下に示すように式（３）に代入すればよい。
ｙi ＝ｃ1φm（ｚi1）＋ｃ2φm（ｚi2）＋ｃ3φm（ｚi3） ・・・（１５）
こうして、特定の入力データｘi に対応する出力データｙi 、すなわち任意の条件に対する特性値の計算が可能となり、最適解の探索が可能となる。
【００６４】
次に、本実施の形態の曲面モデルの同定方法を適用した結果を従来手法と比較して示す。図５〜図７は、実験計画法によって収集した実験データ（製剤特性データ；３要因）に従来手法と本実施の形態の手法を適用した結果得られた曲面モデルを応答曲面として図示した図である。図５（ａ）、図６（ａ）、図７（ａ）は従来の２次多項式を用いて同定した曲面を示し、図５（ｂ）、図６（ｂ）、図７（ｂ）は従来のニューラルネットワークを用いて同定した曲面を示し、図５（ｃ）、図６（ｃ）、図７（ｃ）は本実施の形態の同定方法を用いて同定した曲面を示している。
【００６５】
図５〜図７の例では、喘息の治療用の気管支拡張剤の最適製剤設計に応答曲面法を適用している。気管支拡張剤は、患者が錠剤を服用したときに錠剤の崩壊を抑えるゲル化剤と崩壊を促進する崩壊剤とを錠剤の主成分に混合したものである。図５〜図７において、入力変量である要因パラメータ１はゲル化剤の量であり、同じく入力変量である要因パラメータ２は崩壊剤の量である。また、図５における出力変量である特性値１は、速放率であり、図６における出力変量である特性値２は速放性（錠剤の崩壊が促進される特性）の速度であり、図６における出力変量である特性値３は徐放性（錠剤の崩壊が抑制される特性）の速度である。
【００６６】
図５（ａ）、図６（ａ）、図７（ａ）に示す２次多項式の応答曲面では、製剤設計の専門家の経験に反する、現実ではありえないような特性値が出力される。また、図５（ｂ）、図６（ｂ）、図７（ｂ）に示すニューラルネットワークの応答曲面では、伝達関数となっているシグモイド関数の形状がそのまま現れており、所望の応答曲面とはなっていない。これに対して、図５（ｃ）、図６（ｃ）、図７（ｃ）に示す本実施の形態の曲面では、専門家の経験と合致する、所望の滑らかな応答曲面が得られている。
【００６７】
以上のように、重調和スプラインを利用すれば、非線形な最適化対象に対し、ユーザがパラメータの指定をすることなく、自動的に最適設計に必要な近似関数（曲面モデル）を得ることができる。このため、従来手法と比べて格段に効率のよい最適設計が可能となる。また、従来手法と異なり、現実にありえないような応答となることが少ないため、最適化対象の知識があまりない人であっても、効率よく精度の高い曲面モデルを得ることができ、条件設計の誤りを回避することができる。
なお、本実施の形態では、ｎ＝３の場合を例にとって説明しているが、これに限るものではなく、ｎが２以上であれば本発明を適用することができる。
【００６８】
［実施の形態の２］
実施の形態の１の曲面モデルの同定方法では、データのスムージング（収集したすべてのデータを必ずしも通らない関数近似）が可能である。データのスムージングのアルゴリズムとしては、以下のようなものがある。
【００６９】
（Ｉ）収集したデータを間引いて（サンプリングして）補間関数を求める方法（文献１５）。
（ＩＩ）特異値分解を利用する方法（文献１６）。
（ＩＩＩ）前述の式（４）の代わりに、以下の式（５）を利用する方法（文献１６）。
【００７０】
【数１０】
【００７１】
式（１６）において、∇2 はラプラシアン演算子である。また、Ｔ／Ｄの比を変えてスムージングの程度を制御することができる。
図８に本実施の形態のデータスムージングを適用した結果を示す。図８（ａ）は入力変量が１種類（ｎ＝２）の場合、図８（ｂ）、図８（ｃ）は入力変数が２種類（ｎ＝３）の場合の結果である。なお、図８（ａ）において、破線は全収集データを補間した結果を示し、実線はスムージングをしながら補間した結果を示す。
【００７２】
このように、収集したデータにノイズが含まれる場合には、全ての収集データを必ず通るような曲面モデルを求めるよりも、本実施の形態のスムージングを適用して曲面モデルを求める方が有効であり、所望の応答曲面を得ることができる。なお、収集したデータの信頼性が高い場合はこの限りでない。
【００７３】
［実施の形態の３］
図９は本発明の第３の実施の形態となる曲面モデルの同定装置の構成を示すブロック図であり、図１と同一の構成には同一の符号を付してある。
本実施の形態は、図１の同定装置に、距離計算部５、グリーン関数計算部６、評価指標計算部７及びデータ選択部８を加えたものである。実施の形態の１と同様に、本実施の形態の同定装置もコンピュータで実現することができる。
【００７４】
以下、本実施の形態の同定装置の動作を図１０を参照して説明する。まず、同定装置の距離計算部５は、ユーザによって収集データが入力されると（ステップ３０１においてＹＥＳ）、ｎ−１次元の空間においてｎ−１次元の各入力データ間の距離ｚを入力データの全ての組み合わせについて計算して、図示しない記憶装置に格納する（ステップ３０２）。このときの計算は、実施の形態の１の距離計算部１における処理と全く同じである。
【００７５】
続いて、グリーン関数計算部６は、距離計算部５で計算された距離ｚ毎に対応するグリーン関数値φm（ｚ） を計算して、図示しない記憶装置に格納する（ステップ３０３）。このときの計算も、実施の形態の１のグリーン関数計算部２における処理と全く同じである。
【００７６】
評価指標計算部７は、グリーン関数計算部６で計算されたグリーン関数値φm（ｚ） を基に、ユーザによって収集され同定装置に入力されたデータが曲面モデルの同定にとって有効なデータか否かを示す評価指標値を計算して、図示しない記憶装置に格納する（ステップ３０４）。
【００７７】
評価指標値の計算には、多項式を利用する場合に適用されている各種最適基準（Ａ，Ｄ，Ｇ，Ｑなどの最適基準；文献３「轟 章、応答曲面法、平11-12月日本機械学会講習会「応答曲面法による非線形問題の最適設計入門」資料、URL：http://ueno.mes.titech.ac.jp/rec-res.html にて公開中」）を利用する。ただし、文献３で利用されている入力変数値の行列ではなく、グリーン関数値の行列を使用する。例えば、Ａ最適基準の場合、評価指標値は、（ＧT Ｇ）／ｎの逆行列の対角成分の和となる。なお、ここでのｎは実験データ数（有効な実験データ数）である。また、行列Ｇは以下のような式を利用する。
【００７８】
【数１１】
【００７９】
次に、データ選択部８は、評価指標計算部７で計算された評価指標値が最適（最大もしくは最小）かどうかを判定し（ステップ３０５）、評価指標値が最適でない場合（すなわち、距離計算部５に入力されたデータが曲面モデルの同定にとって最適でない場合）、ユーザによって収集され同定装置に入力されたｎ次元データ中から一部のデータを選択的に取り出して、距離計算部５に与える（ステップ３０６）。これにより、ステップ３０２〜３０５の処理が再び実行される。
【００８０】
このように、ユーザによって収集され同定装置に入力されたｎ次元データ中から一部のｎ次元データのみを選択して評価指標値を計算することを繰り返すことにより、評価指標値が最適となるようなデータを選択する。そして、データ選択部８は、評価指標値が最適となったｎ次元データ（ｎ−１次元の入力データと１次元の出力データ）を距離計算部１に与える。以降の動作は実施の形態の１と全く同じであるので、説明は省略する。なお、Ａ最適基準を用いる場合は、対角成分の和が最小のものが最適であるが、別の基準を用いる場合は、最大のものが最適である場合もある。
【００８１】
以上により、本実施の形態では、過去に収集したデータの中からモデル同定にとって有効なデータのみを適切に選択でき、曲面モデルを誤って同定することを回避することができる。
また、評価指標計算部７で計算された評価指標値を用いて、ユーザがデータの収集方法を見直すことも可能である。
【００８２】
なお、実施の形態の１〜３の曲面モデルの同定方法は、少ないデータから効率よく最適設計を行うための応答曲面法として利用することができ、以下のような様々な分野に適用することができる。
（Ａ）医薬品の最適製剤設計（文献３）。
（Ｂ）複合材（文献１）。
（Ｃ）構造設計（文献２）。
（Ｄ）塗料（文献１７「長倉 稔、塗料・塗装技術者のための実験計画法（第２０回）、塗装技術、1998年２月号」）、化粧品などの配合率最適設計。
（Ｅ）食品加工（文献１８「相島 鐵郎、最適化手法のフレーバー開発への応用、月刊フードケミカル、1996-3」）、機械加工分野などの製造設備の最適設計、運転条件設計／調整。
（Ｆ）その他一般のシステム（プロセス分野、ビル空調、水処理、ごみ処理施設などの設備）の最適設計・運転条件設計／調整。
（Ｇ）経済分野などにおける意思決定問題（状況に応じて採るべき行動の支援など）。
（Ｈ）その他一般の最適設計問題、意思決定問題。
【００８３】
【発明の効果】
本発明によれば、従来の応答曲面法の近似関数として、グリーン関数に基づく重調和スプライン補間を適用しているので、ユーザは近似関数の次数や構造などの各種パラメータの指定をすることなく、収集したデータのみを使って非線形な最適化対象であっても所望の曲面モデルを求めることができる。これにより、より効率（労力削減、時間短縮）のよい応答曲面法による最適設計が可能となる。また、少ない時間でも精度よく所望の曲面モデルが得られるので、多項式近似のような現実的にはありえない近似を防ぐことができ、最適化対象の知識が少ない人であっても、条件設計のミスを回避することができる。さらに、重調和スプラインを用いたことにより、多次元への拡張が容易となり、多入力（多要因）の対象にも柔軟に対応することができる。また、アルゴリズムが単純であり、同一アルゴリズムで多次元化が可能であるため、プログラムがよりコンパクトにでき、さらに最適設計のための計算のパフォーマンスを向上させることができる。
【００８４】
また、出力計算時距離計算ステップ、出力計算時グリーン関数計算ステップ及び出力データ計算ステップを実行することにより、特定の入力データに対応する出力データ（補間値）を計算することができる。
【００８５】
また、第２の同定時距離計算ステップ、第２の同定時グリーン関数計算ステップ、評価指標値計算ステップ及びデータ選択ステップを実行した後に、この選択されたデータを前記収集されたデータとして、第１の同定時距離計算ステップ、第１の同定時グリーン関数計算ステップ及びモデル同定ステップを実行することにより、過去に収集したデータの中からモデル同定にとって有効なデータのみを適切に選択でき、曲面モデルを誤って同定することを回避することができる。
【図面の簡単な説明】
【図１】 本発明の第１の実施の形態となる曲面モデルの同定装置の構成を示すブロック図である。
【図２】 図１の曲面モデルの同定装置の動作を示すフローチャート図である。
【図３】 本発明の第１の実施の形態における距離計算部のモデル同定時の動作を説明するための説明図である。
【図４】 本発明の第１の実施の形態における距離計算部の出力データ計算時の動作を説明するための説明図である。
【図５】 実験データに従来の手法と本発明の第１の実施の形態の手法を適用した結果得られた曲面モデルを応答曲面として図示した図である。
【図６】 実験データに従来の手法と本発明の第１の実施の形態の手法を適用した結果得られた曲面モデルを応答曲面として図示した図である。
【図７】 実験データに従来の手法と本発明の第１の実施の形態の手法を適用した結果得られた曲面モデルを応答曲面として図示した図である。
【図８】 本発明の第２の実施の形態におけるデータスムージングの効果を示す図である。
【図９】 本発明の第３の実施の形態となる曲面モデルの同定装置の構成を示すブロック図である。
【図１０】 図９の曲面モデルの同定装置の動作を示すフローチャート図である。
【図１１】 応答曲面法における処理の流れを示すフローチャート図である。
【符号の説明】
１…距離計算部、２…グリーン関数計算部、３…モデル係数計算部、４…出力データ計算部、５…距離計算部、６…グリーン関数計算部、７…評価指標計算部、８…データ選択部。
【発明の属する技術分野】
本発明は、製剤設計、材料設計、システムの運転条件・製品製造条件の設計／調整等に使用される応答曲面法などに利用可能な、曲面モデルの同定方法及びプログラムに関するものである。
【０００２】
【従来の技術】
従来、実験計画法に基づいて収集されたデータを用いて、効率よく製品の製造条件・システムの運転条件等を設計・調整するために、応答曲面法（Response Surface Methodology；RSM ）が用いられている。以下、従来の応答曲面法（文献１「轟 章、応答曲面法、平11-12月日本機械学会講習会「応答曲面法による非線形問題の最適設計入門」資料、URL：http://ueno.mes.titech.ac.jp/rec-res.html にて公開中」、文献２「長谷川ほか、実数型交差モデルと応答曲面モデルを用いた逐次的近似最適化に関する一試行、日本計算工学会論文集、２０００年５月２４日発行」、文献３「Takayama、Artificial Neural Network as a Novel Method to Optimize Pharmaceutical Formulations、Pharmaceutical Research, Vol.16, No.1, 1999」など）について説明する。
【０００３】
図１１は応答曲面法における処理の流れを示すフローチャート図である。応答曲面法には、実験計画法に基づき必要最低限のデータを収集する処理（ステップ４０１）と、収集されたデータをもとに各要因（設計条件）とその応答（最適にしたい特性値）を結びつける近似関数（曲面モデル）を同定する処理（ステップ４０２）と、得られた近似関数（曲面モデル）を利用して、所望の特性値となるような最適設計条件を各種最適化アルゴリズムを利用して求める処理（ステップ４０３）の３つの処理ステップがある。
【０００４】
次に、上記ステップ４０１〜４０３において、従来利用されている技術を以下に説明する。ステップ４０１では、実験計画法に基づき、実験計画データを収集する。これにより必要最低限なデータを効率よく収集する。利用される実験計画法の要因配置（実験点配置）としては、完備要因計画（全因子計画；fully factorial design）法、直交計画（orthogonal design ）法、球形中心２次複合計画（中央複合計画；central composite design）法などが代表的である（文献１、文献４「橋田 充 編、経口投与製剤の処方設計（第7章 処方・製法の最適化）、薬業時報社、１９９８」）。
【０００５】
ステップ４０２では、ステップ４０１で収集したデータを使って、各要因（設計条件）ｘi （ｉ＝１，２，・・・，ｋ）とその応答（最適にしたい特性値）ｙとを結びつける近似関数（曲面モデル）ｙ＝ｆ（ｘ1，ｘ2，・・・，ｘk ）を以下のような手法を利用して求める。
【０００６】
（Ｉ）多項式を用いた手法（文献１、文献２）は、各要因（入力変数）と特性値（出力変数）とを基本的には線形の関数で近似する手法である。多くの場合は、以下の２次多項式が利用される。
【０００７】
【数１】
式（１８）において、ａ0，ｂi，ｃi，ｃij は定数である。非線形な対象に対応する場合には、相互作用項（式（１８）の右辺第４項）を導入したり、べき乗関数、指数関数あるいは対数関数などを利用して、入力変数の変換を行ったりする。前記２次多項式が、応答曲面法の近似関数としては最も頻繁に利用される手法である。
【０００９】
（ＩＩ）ニューラルネットワークを用いた手法（文献３）は、非線形な入出力関係を近似する手法である。入力を各要因（設計条件）、出力を特性値として、階層型のニューラルネットワークにバックプロパゲーション法を利用して学習させる。入力層から中間層、出力層に至る信号は、前層からの信号を次の式（１９）により統合し、式（２０）のようなシグモイド関数に代表される伝達関数により、次の層へと出力される。
【００１０】
【数２】
【数３】
なお、ｗpqは次層ユニットｑと前層ユニットｐ間の荷重、ｘp は前層からの出力であり、ｆ（ｙq ）は次層への出力として伝達される。また、ａはシグモイド関数の勾配係数である。
【００１３】
（ＩＩＩ）スプライン関数を用いた手法（文献５「桜井 明著、スプライン関数入門、東京電機大学出版局、１９８１」、文献６「桜井 明監修、パソコンによるスプライン関数、東京電機大学出版局、１９８８」）は、非線形な特性であっても関数化することができ、滑らかな曲線/曲面が得られる手法である。通常広く認知されているスプラインは、文献「マグローヒル科学技術用語大辞典」に記載されている通り、導関数に基づくスプライン関数のことである。この導関数に基づくスプライン関数の数学的定義は、文献４に記載されている通り、「ｍ次のスプライン関数は、そのｍ階微分が階段関数で、ｍ−１階以下の微分が連続であるような関数である」というものである。また、ｎ個の接点をもつｍ次のスプライン関数は以下のようになる。
【００１４】
【数４】
なお、（ｘ−ｑi）m +は切断べき関数と呼ばれる関数であり、以下のように表される。
【００１６】
【数５】
（ＩＩＩＩ）多変数適応回帰スプライン（MARS）を用いた手法（文献７「T.J.Hastie and R.J.Tibshirani、Nonparametric Regression and Classification Part I - Nonparametric Regression, NATO ASI Series. F. Computer and System Sciences, Vol.136, pp.62-69, 1994」）は、多項式の回帰近似を区分に分割して行う手法である。分割する領域の大きさは計算により自動的に決定される。
【００１８】
最後に、ステップ４０３では、得られた曲面モデルを利用して、公知の各種最適化アルゴリズムにより、所望の最適値探索を行い、最適設計条件を見出す。連続量に対する公知の最適化アルゴリズムとしては、シンプレックス法、勾配法、共役勾配法、ニュートン法、準ニュートン法、シミュレーテッド・アニーリング法、トンネリング・アルゴリズム、遺伝的アルゴリズム（GA）などが代表的である（文献８「茨木俊秀ほか、最適化の手法、共立出版、１９９３」、文献９「相吉 英太郎、「最適化」の展望、計測と制御、Vol.29,No.12、１９９０」）。
【００１９】
また、制約条件下での最適化手法として、ペナルティ関数を利用した方法（文献３、文献８）などもある。さらに、複数の曲面モデルを利用することで、複数特性の最適値（妥協解）を見出す多目的最適化（文献３、文献１０「西川ほか、岩波講座 情報科学１９「最適化」、岩波書店、１９８２」）も実施することができる。多目的最適解を求める方法としては、重みパラメータ法、ε制約法、目標計画法、対話型の方法などがある。
【００２０】
【発明が解決しようとする課題】
応答曲面法は、数少ないデータ（必要最低限のデータ）から、効率よく最適条件を見出すために利用される。しかしながら、従来の応答曲面法では、前述のステップ４０２における近似関数（曲面モデル）の同定に労力を要することがあり、実際には効率が悪い可能性があるという問題点があった。例えば、前述の従来の近似関数を利用した場合には、以下のような問題がある。
【００２１】
（Ｉ）多項式を用いた手法は、基本的には、入出力関係が線形な対象をモデル化する手法であるため、各要因と特性値が非線形な関係の場合（現実にはこのケースが多い）は、効率よくモデル化することが困難となる場合がある。非線形な対象をモデル化する手段として、相互作用項を追加したり、べき乗関数、指数関数あるいは対数関数などを導入して近似することもあるが、必要な精度のモデルを得るためには多大な労力を要することがある。例えば、相互作用項を導入する場合には、すべての要因の相互作用項を追加すべきか、あるいは一部の相互作用項とするか、一部の相互作用項とする場合には、どの要因の組み合わせを追加すればよいかということは試してみないと分からない。また、どのような関数を導入すればよいかについても試してみないと分からないことがある。つまり、近似関数を得るためには、試行錯誤的となり効率が悪い。
【００２２】
さらに、非線形性の強い対象の場合は、全領域を１つの関数で近似することが困難な場合がある。そのため、小さな領域に分けてそれぞれに異なる多項式近似関数をあてはめることもある。しかし、対象の特性が未知なものについて、どの領域がどのような関数に相当するかを決めることは困難である。そのため，結果的に所望の精度のモデルを得ることができない。また、２次多項式による近似では、現実的にはあり得ない値（特性値の限界値以上あるいは限界値以下の値）の応答が現れる場合があり、所望の応答曲面にならない場合がある。この場合、最適な領域の分割や関数の選択にさらに十分な時間をかければ、所望のモデルが得られる可能性もあるが、非常に多大な労力が必要となり、モデル同定の効率が悪化する。
【００２３】
（ＩＩ）ニューラルネットワークを用いた手法では、ユーザが指定すべきパラメータが多すぎ、また、パラメータを決定するための指標がないので、ネットワークモデルを試行錯誤的に求めざるを得なくなり、実際には効率よく所望の近似関数を求めることは困難である。例えば、ユーザが予め設定する必要があるパラメータとして代表的なものには、ネットワークの階層数、中間層のユニット数、伝達関数の関数式（伝達関数をシグモイド関数とする場合には、その関数の傾き）、適用する学習法則、学習停止条件、学習スタート時の荷重係数の初期値などがある。
【００２４】
このようなパラメータを変化させると、応答曲面が大きく変化するが、パラメータをどう変化させればよりよいモデルが得られるかを知るための明確な指標がないため、結果的には試行錯誤的にモデルを作らざるを得ない。さらに、学習計算は繰り返し計算（収束計算）となるため、計算時間が必要となることがある。ニューラルネットワークで近似した応答曲面では、応答曲面が伝達関数であるシグモイド関数の形状となってしまい、所望の応答曲面にならない場合がある。この場合、よりよい応答曲面にするにはさらに時間がかかり、モデル同定の効率が悪化する。
【００２５】
（ＩＩＩ）通常、工学分野で広く認知・利用されている導関数に基づくスプライン関数では、１要因（１入力変数）に対する出力値の補間曲線（近似関数）を描くことは容易にできるが、２要因（２入力変数）以上の場合には、以下のような制限が発生する。
（α）補間計算（近似関数の同定）をするために与えられるデータが、入力空間内で格子状となるデータであり、各座標点に順序をつける必要がある（文献１１「William H. Pressほか、Numerical Recipes in C、技術評論社、１９９３」）。
（β）与えられたデータに囲まれる区分領域（２入力変数の場合は、３角領域）を切り出す必要がある（文献１２「曲渕 英邦、スプライン関数を用いたランダム布置標本データからの２次元補間、日本建築学会大会学術講演梗概集（九州）１９９８年９月」）。
【００２６】
実験データの近似を目的とするならば、（α）の格子状データが得られることは期待できないため、不規則分布のデータを前提とした（β）を適用せざるを得ない。しかし、２入力変数以上の場合、補間計算以外に領域切り出しのアルゴリズム（例えば、Delaunay三角分割；文献１３「MATLAB ver5.3 オンラインマニュアル、Function Reference "delaunay" 」）が必要であり、また各次元ごとに（入力変数の数に応じて）領域切り出しのアルゴリズムが必要である（２入力変数の場合には例えば３角形分割、３入力変数の場合には例えば４面体領域切り出しなど）。つまり、アルゴリズムが高次元になればなるほど、アルゴリズムが複雑となり、プログラムも大きくなり、計算時間も必要となるため、現実的には３入力変数以上の補間にはほとんど利用されていない。
【００２７】
（ＩＩＩＩ）多変数適応回帰スプライン（MARS）を用いた手法は、基本的には多項式近似と同様なアルゴリズムであるため、多項式と同様な問題がある。ただし、区分領域での近似を採用していることから、非線形近似にもある程度対応可能である。また、近似の基本となる関数もあらかじめ候補をいくつか用意しておけば、そのなかで最適化なものを自動的に選択することができる。しかし、近似の基本となる関数の候補として何を用意すればよいかわからない場合は、多数の関数候補を用意する必要があり、これによって計算負荷（計算時間、要求メモリサイズ）が増大する。また、非線形性が強い場合などには、より多数の区分領域に分けて回帰近似する必要があり、これによっても計算負荷が大きくなるため、実用上問題である。また、領域を区分して多項式近似をするため、実験データを細分して利用することとなる。このため、全領域において、実験計画法に基づき収集したデータでは、多項式近似に必要なデータ量が得られないことがある。また、細分した領域ごとに実験計画法を適用しデータを収集することは、実験コストの増大や効率の悪化となるため、実用上問題である。
【００２８】
なお、他の補間手法（近似関数）として、ラグランジュ補間多項式もあるが、この手法を利用した場合には、データ点数によっては、ルンゲの現象という振動現象が現れる。これにより、なめらかではない不自然な、応答の激しい曲面となってしまうので（文献１４「桜井 明監修、パソコンによるスプライン関数、東京電機大学出版局、１９８８、pp.4-6」）、実用上問題である。
【００２９】
以上のように、従来の応答曲面法では、モデル同定が試行錯誤的となり、モデル同定に多大な時間がかかるという問題点があった。また、従来の応答曲面法では、モデル同定を誤る可能性があった。
本発明は、上記課題を解決するためになされたもので、高精度な曲面モデルの同定を効率良く行うことができる曲面モデルの同定方法及びプログラムを提供することを目的とする。
【００３０】
【課題を解決するための手段】
本発明は、ｎ（ｎは２以上の整数）種類の変量間の関係を近似したｎ次元の曲面モデルをコンピュータを使用して同定する曲面モデルの同定方法であって、予めデータが収集された前記ｎ種類の変量のうちｎ−１種類の変量で定義されるｎ−１次元の空間において各データ間の距離を計算して記憶装置に記憶させる第１の同定時距離計算ステップ（３０２）と、この第１の同定時距離計算ステップで計算された各距離に対応するグリーン関数値を計算して記憶装置に記憶させる第１の同定時グリーン関数計算ステップ（３０３）と、この第１の同定時グリーン関数計算ステップで計算されたグリーン関数値から前記収集されたデータが前記曲面モデルの同定にとって有効なデータか否かを示す評価指標値を計算して記憶装置に記憶させる評価指標値計算ステップ（３０４）と、この評価指標値が最大となるようなデータを前記収集されたデータ中から選択するデータ選択ステップ（３０５，３０６）と、この選択されたデータについて、前記ｎ種類の変量のうちｎ−１種類の変量で定義されるｎ−１次元の空間において各データ間の距離を計算して記憶装置に記憶させる第２の同定時距離計算ステップ（１０２）と、この第２の同定時距離計算ステップで計算された各距離に対応するグリーン関数値を計算して記憶装置に記憶させる第２の同定時グリーン関数計算ステップ（１０３）と、前記選択されたデータのうち前記ｎ−１種類を除く１種類のデータと前記第２の同定時グリーン関数計算ステップで計算されたグリーン関数値とから前記ｎ次元の曲面モデルを計算して記憶装置に記憶させるモデル同定ステップ（１０４）とからなるものである。このように、本発明は、従来の応答曲面法（実験計画法、近似関数、既存最適化アルゴリズムを組み合わせて、少ない実験データから効率よく最適解を求める手法）の近似関数として、グリーン関数に基づく重調和スプライン補間を適用している。重調和スプライン補間は、与えられたデータのみを使って非線形な特性であっても応答曲面を求めることができ、ユーザは近似関数の次数、構造などの各種パラメータの指定をする必要がない。また、重調和スプライン補間は、多次元への拡張が容易であり、多要因のモデルにも柔軟に対応することができる。
【００３１】
また、本発明の曲面モデルの同定方法の１構成例において、前記ｎ種類の変量のデータは、実験計画法によって収集されるものである。
【００３２】
また、本発明のプログラムは、予めデータが収集された前記ｎ種類の変量のうちｎ−１種類の変量で定義されるｎ−１次元の空間において各データ間の距離を計算して記憶装置に記憶させる第１の同定時距離計算ステップと、この第１の同定時距離計算ステップで計算された各距離に対応するグリーン関数値を計算して記憶装置に記憶させる第１の同定時グリーン関数計算ステップと、この第１の同定時グリーン関数計算ステップで計算されたグリーン関数値から前記収集されたデータが前記曲面モデルの同定にとって有効なデータか否かを示す評価指標値を計算して記憶装置に記憶させる評価指標値計算ステップと、この評価指標値が最大となるようなデータを前記収集されたデータ中から選択するデータ選択ステップと、この選択されたデータについて、前記ｎ種類の変量のうちｎ−１種類の変量で定義されるｎ−１次元の空間において各データ間の距離を計算して記憶装置に記憶させる第２の同定時距離計算ステップと、この第２の同定時距離計算ステップで計算された各距離に対応するグリーン関数値を計算して記憶装置に記憶させる第２の同定時グリーン関数計算ステップと、前記選択されたデータのうち前記ｎ−１種類を除く１種類のデータと前記第２の同定時グリーン関数計算ステップで計算されたグリーン関数値とから前記ｎ次元の曲面モデルを計算して記憶装置に記憶させるモデル同定ステップとをコンピュータに実行させるものである。
【００３３】
【発明の実施の形態】
［実施の形態の１］
本発明では、応答曲面法において効率よく最適条件を探索するため、効率よく近似関数（曲面モデル）を得ることを目的としており、衛星データから等高線を得るために利用されている、グリーン関数に基づく重調和スプラインを近似関数として適用している。
【００３４】
グリーン関数に基づく重調和スプラインについては、例えば文献１５「David T. Sandwell、Biharmonic Spline Interpolation of GEOS-3 and SEASAT Altimeter Data、Geophysical Research Letters, Vol.14, No.2, pp.139-142, 1987」、また文献１６「Paul Wessel, David Bercovici, Interpolation with Splines in Tension: A Green's Function Approach, Mathmatical Geology, Vol.30, No.1, 1998 」に記載されている。
【００３５】
重調和スプラインは、工学分野で広く認知・利用されている、導関数に基づくスプライン関数とは、理論的にまったく異なる補間方法である。しかし、補間曲面は、導関数に基づくスプライン関数と同様に滑らかな曲面であり、非線形な特性でもモデル化することができる。
【００３６】
従来工学の分野で利用されているスプライン関数は、導関数に基づく補間値計算アルゴリズムを有するスプライン関数である。これに対し、重調和スプラインは、積分定理（グリーン定理）に基づく補間値計算アルゴリズムとなっている。そのため、従来のスプラインとは異なり、不規則分布データであっても特別なアルゴリズムを用意することなく、補間計算可能である。さらに、この手法は、データを与えるだけで、データの分布のみから関数を計算するため、ユーザが予め設定する必要があるパラメータは存在しない。また、多次元の入力変数（多要因）でも２入力変数（２要因）以下の場合と同様なアルゴリズムで動作するため、多次元（多要因）な対象にも、柔軟に適用可能である。
【００３７】
以下、この重調和スプライン補間のアルゴリズムを記述する。文献１５によると、任意の入力変数ベクトルｘに対する重調和スプライン関数ｆ（ｘ）は次式のようになる。
【００３８】
【数６】
ただし、Ｎは予め収集されたデータｘj（ｊ＝１〜Ｎ） の数、ｚは任意のデータｘと予め収集されたデータｘj との距離である。また、φm（ｚ） はグリーン関数であり、距離ｚとデータｘ，ｘj の次元数（データｘ，ｘj を特徴付ける変量の数）ｍとに応じて、次式のように定められる。
【００４０】
【数７】
式（１）の係数ｃj は、データが与えられると線形マトリクスを解くことにより自動的に算出できる。これにより求めたい入力変数ベクトルｘi に応じた出力値（補間値）ｙi は、次式により算出できる。
【００４２】
【数８】
式（３）において、ｚijはデータｘi と予め収集されたデータｘj との距離である。なお、式（１）は、次式より導出される。
【００４４】
【数９】
ただし、∇4 は重調和演算子、δ（ｚ）はデルタ関数である。
以上のように、重調和スプラインでは、データが与えられると、入力空間におけるデータ間の距離と、各次元（入力変数の数）に応じたグリーン関数とから、自動的に関数値が計算できる。そのため、不規則分布データであっても、さらに多次元（多要因）データであっても、特別なアルゴリズムを用意する必要はなく、柔軟に多次元拡張可能である。また線形マトリクス計算であるため、計算負荷も大きくない。
【００４６】
次に、以上のような重調和スプライン補間のアルゴリズムを用いた曲面モデルの同定方法について詳細に説明する。図１は本発明の第１の実施の形態となる曲面モデルの同定装置の構成を示すブロック図、図２は図１の同定装置の動作を示すフローチャート図である。図２（ａ）は同定装置のモデル同定時の動作を示し、図２（ｂ）は出力データ計算時の動作を示している。
【００４７】
図１の同定装置は、例えば実験計画法によって予めデータが収集されたｎ（ｎは２以上の整数）種類の変量のうちｎ−１種類の変量で定義されるｎ−１次元の空間において各データ間の距離を計算する距離計算部１と、この距離計算部１で計算された各距離に対応するグリーン関数値を計算するグリーン関数計算部２と、前記収集されたデータのうち前記ｎ−１種類を除く１種類のデータと前記グリーン関数値とから前記ｎ種類の変量間の関係を近似したｎ次元の曲面モデルを求めるモデル係数計算部３と、グリーン関数値と同定されたｎ次元の曲面モデルとを用いて、特定の入力データに対応する１次元の出力データを補間計算する出力データ計算部４とを有している。
【００４８】
図１の同定装置は、コンピュータ上で実現することができる。コンピュータは、ＣＰＵ、ＲＯＭ、ＲＡＭ、ディスプレイ装置やキーボードあるいは外部記憶装置とのインタフェースをとるための回路などを備えた周知の構成のものでよい。ＣＰＵは、ＲＯＭ若しくはＲＡＭに記憶されたプログラム、又はキーボードから入力されたコマンドに従って処理を実行する。また、ＣＰＵは、外部記憶装置にデータを書き込んだり、外部記憶装置からデータを読み出したりすることができる。
【００４９】
本実施の形態の曲面モデルの同定方法を実現させるためのプログラムは、フレキシブルディスク、ＣＤ−ＲＯＭ、メモリカード等の記録媒体に記録された状態で提供される。この記録媒体を外部記憶装置に挿入すると、媒体に記録されたプログラムが読み取られ、コンピュータに転送される。そして、ＣＰＵは、読み込んだプログラムをＲＡＭ等に書き込む。こうして、同定装置は、以下で説明するような処理を実行する。
【００５０】
以下、このような同定装置のモデル同定時の動作を図２（ａ）を参照して説明する。まず、同定装置のユーザは、ｎ種類の変量によって特徴付けられるｎ次元データを例えば実験計画法によって複数個収集し、前記ｎ種類の変量をｎ−１種類と１種類の変量に分け、これに対応して前記ｎ次元データをｎ−１次元のデータと１次元のデータとに分けて同定装置に入力する。
【００５１】
一般的には、前記ｎ−１種類の変量はユーザが求める解の要因となる変量（例えば、曲面モデルの生成理由が製剤処方の最適化にあるとすれば、要因は各原材料の量）であり、前記１種類の変量は解となる変量（処方された製剤の特性値）である。本発明では、ｎ−１種類の変量が要因（入力）で、残りの１種類の変量が解（出力）である必要は必ずしもないが、ここでは便宜的に、ｎ−１種類の変量によって特徴付けられるｎ−１次元データを入力データ、残りの１種類の変量によって特徴付けられる１次元データを出力データと呼ぶ。
【００５２】
同定装置の距離計算部１は、ユーザによってデータが入力されると（ステップ１０１においてＹＥＳ）、ｎ−１次元の空間においてｎ−１次元の各入力データ間の距離ｚを入力データの全ての組み合わせについて計算して、図示しない記憶装置に格納する（ステップ１０２）。ここでの距離ｚは、ユークリッド距離である。
【００５３】
例えば、ｎ＝３として、入力の２変量をＸ1，Ｘ2、出力の１変量をＹとし、図３に示すように、２次元の入力データｘ1，ｘ2，ｘ3 の値をそれぞれ（Ｘ11，Ｘ21），（Ｘ12，Ｘ22），（Ｘ13，Ｘ23）とし、入力データｘ1，ｘ2，ｘ3 に対応する出力データをそれぞれｙ1，ｙ2，ｙ3 とすれば、入力データｘ1，ｘ2間の距離ｚ12（又はＺ21）、入力データｘ2，ｘ3間の距離ｚ23（又はＺ32）、入力データｘ3，ｘ1間の距離ｚ31（又はＺ13）はそれぞれ次式のように得られる。
【００５４】
ｚ12＝Ｚ21＝（（Ｘ11−Ｘ12）2＋（Ｘ21−Ｘ22）2）1/2 ・・・（５）
ｚ23＝Ｚ32＝（（Ｘ12−Ｘ13）2＋（Ｘ22−Ｘ23）2）1/2 ・・・（６）
ｚ31＝Ｚ13＝（（Ｘ13−Ｘ11）2＋（Ｘ23−Ｘ21）2）1/2 ・・・（７）
【００５５】
続いて、グリーン関数計算部２は、距離計算部１で計算された距離ｚ毎に対応するグリーン関数値φm（ｚ） を計算して、図示しない記憶装置に格納する（ステップ１０３）。前述のｎ＝３の例では、距離ｚ12，ｚ23，ｚ31に応じてグリーン関数値φm（ｚ12），φm（ｚ23），φm（ｚ31） を計算することになる。前述のとおり、グリーン関数値φm（ｚ） は式（２）によって計算することができる。なお、本実施の形態では、ｍ＝ｎ−１である。
【００５６】
モデル係数計算部３は、出力データｙとグリーン関数値φm（ｚ） とからｎ次元の曲面モデルのモデル係数ｃj （ｊ＝１〜Ｎ）を計算して、図示しない記憶装置に格納する（ステップ１０４）。モデル係数ｃj を求めるには、出力データｙとグリーン関数値φm（ｚ） とを式（３）に代入すればよい。前述のｎ＝３の例ではＮ＝３なので、以下のような連立方程式が得られる。
【００５７】
ｙ1 ＝ｃ2φm（ｚ12）＋ｃ3φm（ｚ13） ・・・（８）
ｙ2 ＝ｃ1φm（ｚ21）＋ｃ3φm（ｚ23） ・・・（９）
ｙ3 ＝ｃ1φm（ｚ31）＋ｃ2φm（ｚ32） ・・・（１０）
【００５８】
モデル係数計算部３は、式（８）〜式（１０）の連立方程式から最小二乗法により、モデル係数ｃ1，ｃ2，ｃ3 を計算する。
以上により、モデル同定が終了する。同定装置によって同定されたｎ＝３次元の曲面モデルは、式（１）とモデル係数ｃ1，ｃ2，ｃ3 とにより、次式のように表すことができる。
ｙ＝ｃ1φm（ｚi1）＋ｃ2φm（ｚi2）＋ｃ3φm（ｚi3） ・・・（１１）
【００５９】
次に、ｎ−１（＝ｍ）次元の空間内の特定の入力データｘi に対応する出力データ（補関値）ｙi を同定装置に計算させる場合の動作を図２（ｂ）を参照して説明する。
【００６０】
同定装置の距離計算部１は、ユーザによって入力データｘi が入力されると（図２ステップ２０１においてＹＥＳ）、ｎ−１次元の空間において入力データｘi とモデル同定時に入力された前記収集された各入力データｘ1，ｘ2，ｘ3 との間の距離ｚを図４に示すように入力データｘ1，ｘ2，ｘ3 毎に計算して、図示しない記憶装置に格納する（ステップ２０２）。
【００６１】
入力データｘiの値を（Ｘ1i，Ｘ2i）とすると、入力データｘi，ｘ１間の距離ｚi1、入力データｘi，ｘ2間の距離ｚi2、入力データｘi，ｘ3間の距離ｚi3はそれぞれ次式のように得られる。
ｚi1＝（（Ｘ1i−Ｘ11）2＋（Ｘ2i−Ｘ21）2）1/2 ・・・（１２）
ｚi2＝（（Ｘ1i−Ｘ12）2＋（Ｘ2i−Ｘ22）2）1/2 ・・・（１３）
ｚi3＝（（Ｘ1i−Ｘ13）2＋（Ｘ2i−Ｘ23）2）1/2 ・・・（１４）
【００６２】
続いて、グリーン関数計算部２は、距離計算部１で計算された距離ｚi1，ｚi2，ｚi3毎に対応するグリーン関数値φm（ｚi1），φm（ｚi2），φm（ｚi3） を計算して、図示しない記憶装置に格納する（ステップ２０３）。
【００６３】
出力データ計算部４は、特定の入力データｘi に対応する出力データ（補関値）ｙi を計算して、図示しない記憶装置に格納する（ステップ２０４）。出力データｙi を計算するには、グリーン関数計算部２で計算されたグリーン関数値φm（ｚi1），φm（ｚi2），φm（ｚi3） を以下に示すように式（３）に代入すればよい。
ｙi ＝ｃ1φm（ｚi1）＋ｃ2φm（ｚi2）＋ｃ3φm（ｚi3） ・・・（１５）
こうして、特定の入力データｘi に対応する出力データｙi 、すなわち任意の条件に対する特性値の計算が可能となり、最適解の探索が可能となる。
【００６４】
次に、本実施の形態の曲面モデルの同定方法を適用した結果を従来手法と比較して示す。図５〜図７は、実験計画法によって収集した実験データ（製剤特性データ；３要因）に従来手法と本実施の形態の手法を適用した結果得られた曲面モデルを応答曲面として図示した図である。図５（ａ）、図６（ａ）、図７（ａ）は従来の２次多項式を用いて同定した曲面を示し、図５（ｂ）、図６（ｂ）、図７（ｂ）は従来のニューラルネットワークを用いて同定した曲面を示し、図５（ｃ）、図６（ｃ）、図７（ｃ）は本実施の形態の同定方法を用いて同定した曲面を示している。
【００６５】
図５〜図７の例では、喘息の治療用の気管支拡張剤の最適製剤設計に応答曲面法を適用している。気管支拡張剤は、患者が錠剤を服用したときに錠剤の崩壊を抑えるゲル化剤と崩壊を促進する崩壊剤とを錠剤の主成分に混合したものである。図５〜図７において、入力変量である要因パラメータ１はゲル化剤の量であり、同じく入力変量である要因パラメータ２は崩壊剤の量である。また、図５における出力変量である特性値１は、速放率であり、図６における出力変量である特性値２は速放性（錠剤の崩壊が促進される特性）の速度であり、図６における出力変量である特性値３は徐放性（錠剤の崩壊が抑制される特性）の速度である。
【００６６】
図５（ａ）、図６（ａ）、図７（ａ）に示す２次多項式の応答曲面では、製剤設計の専門家の経験に反する、現実ではありえないような特性値が出力される。また、図５（ｂ）、図６（ｂ）、図７（ｂ）に示すニューラルネットワークの応答曲面では、伝達関数となっているシグモイド関数の形状がそのまま現れており、所望の応答曲面とはなっていない。これに対して、図５（ｃ）、図６（ｃ）、図７（ｃ）に示す本実施の形態の曲面では、専門家の経験と合致する、所望の滑らかな応答曲面が得られている。
【００６７】
以上のように、重調和スプラインを利用すれば、非線形な最適化対象に対し、ユーザがパラメータの指定をすることなく、自動的に最適設計に必要な近似関数（曲面モデル）を得ることができる。このため、従来手法と比べて格段に効率のよい最適設計が可能となる。また、従来手法と異なり、現実にありえないような応答となることが少ないため、最適化対象の知識があまりない人であっても、効率よく精度の高い曲面モデルを得ることができ、条件設計の誤りを回避することができる。
なお、本実施の形態では、ｎ＝３の場合を例にとって説明しているが、これに限るものではなく、ｎが２以上であれば本発明を適用することができる。
【００６８】
［実施の形態の２］
実施の形態の１の曲面モデルの同定方法では、データのスムージング（収集したすべてのデータを必ずしも通らない関数近似）が可能である。データのスムージングのアルゴリズムとしては、以下のようなものがある。
【００６９】
（Ｉ）収集したデータを間引いて（サンプリングして）補間関数を求める方法（文献１５）。
（ＩＩ）特異値分解を利用する方法（文献１６）。
（ＩＩＩ）前述の式（４）の代わりに、以下の式（５）を利用する方法（文献１６）。
【００７０】
【数１０】
式（１６）において、∇2 はラプラシアン演算子である。また、Ｔ／Ｄの比を変えてスムージングの程度を制御することができる。
図８に本実施の形態のデータスムージングを適用した結果を示す。図８（ａ）は入力変量が１種類（ｎ＝２）の場合、図８（ｂ）、図８（ｃ）は入力変数が２種類（ｎ＝３）の場合の結果である。なお、図８（ａ）において、破線は全収集データを補間した結果を示し、実線はスムージングをしながら補間した結果を示す。
【００７２】
このように、収集したデータにノイズが含まれる場合には、全ての収集データを必ず通るような曲面モデルを求めるよりも、本実施の形態のスムージングを適用して曲面モデルを求める方が有効であり、所望の応答曲面を得ることができる。なお、収集したデータの信頼性が高い場合はこの限りでない。
【００７３】
［実施の形態の３］
図９は本発明の第３の実施の形態となる曲面モデルの同定装置の構成を示すブロック図であり、図１と同一の構成には同一の符号を付してある。
本実施の形態は、図１の同定装置に、距離計算部５、グリーン関数計算部６、評価指標計算部７及びデータ選択部８を加えたものである。実施の形態の１と同様に、本実施の形態の同定装置もコンピュータで実現することができる。
【００７４】
以下、本実施の形態の同定装置の動作を図１０を参照して説明する。まず、同定装置の距離計算部５は、ユーザによって収集データが入力されると（ステップ３０１においてＹＥＳ）、ｎ−１次元の空間においてｎ−１次元の各入力データ間の距離ｚを入力データの全ての組み合わせについて計算して、図示しない記憶装置に格納する（ステップ３０２）。このときの計算は、実施の形態の１の距離計算部１における処理と全く同じである。
【００７５】
続いて、グリーン関数計算部６は、距離計算部５で計算された距離ｚ毎に対応するグリーン関数値φm（ｚ） を計算して、図示しない記憶装置に格納する（ステップ３０３）。このときの計算も、実施の形態の１のグリーン関数計算部２における処理と全く同じである。
【００７６】
評価指標計算部７は、グリーン関数計算部６で計算されたグリーン関数値φm（ｚ） を基に、ユーザによって収集され同定装置に入力されたデータが曲面モデルの同定にとって有効なデータか否かを示す評価指標値を計算して、図示しない記憶装置に格納する（ステップ３０４）。
【００７７】
評価指標値の計算には、多項式を利用する場合に適用されている各種最適基準（Ａ，Ｄ，Ｇ，Ｑなどの最適基準；文献３「轟 章、応答曲面法、平11-12月日本機械学会講習会「応答曲面法による非線形問題の最適設計入門」資料、URL：http://ueno.mes.titech.ac.jp/rec-res.html にて公開中」）を利用する。ただし、文献３で利用されている入力変数値の行列ではなく、グリーン関数値の行列を使用する。例えば、Ａ最適基準の場合、評価指標値は、（ＧT Ｇ）／ｎの逆行列の対角成分の和となる。なお、ここでのｎは実験データ数（有効な実験データ数）である。また、行列Ｇは以下のような式を利用する。
【００７８】
【数１１】
次に、データ選択部８は、評価指標計算部７で計算された評価指標値が最適（最大もしくは最小）かどうかを判定し（ステップ３０５）、評価指標値が最適でない場合（すなわち、距離計算部５に入力されたデータが曲面モデルの同定にとって最適でない場合）、ユーザによって収集され同定装置に入力されたｎ次元データ中から一部のデータを選択的に取り出して、距離計算部５に与える（ステップ３０６）。これにより、ステップ３０２〜３０５の処理が再び実行される。
【００８０】
このように、ユーザによって収集され同定装置に入力されたｎ次元データ中から一部のｎ次元データのみを選択して評価指標値を計算することを繰り返すことにより、評価指標値が最適となるようなデータを選択する。そして、データ選択部８は、評価指標値が最適となったｎ次元データ（ｎ−１次元の入力データと１次元の出力データ）を距離計算部１に与える。以降の動作は実施の形態の１と全く同じであるので、説明は省略する。なお、Ａ最適基準を用いる場合は、対角成分の和が最小のものが最適であるが、別の基準を用いる場合は、最大のものが最適である場合もある。
【００８１】
以上により、本実施の形態では、過去に収集したデータの中からモデル同定にとって有効なデータのみを適切に選択でき、曲面モデルを誤って同定することを回避することができる。
また、評価指標計算部７で計算された評価指標値を用いて、ユーザがデータの収集方法を見直すことも可能である。
【００８２】
なお、実施の形態の１〜３の曲面モデルの同定方法は、少ないデータから効率よく最適設計を行うための応答曲面法として利用することができ、以下のような様々な分野に適用することができる。
（Ａ）医薬品の最適製剤設計（文献３）。
（Ｂ）複合材（文献１）。
（Ｃ）構造設計（文献２）。
（Ｄ）塗料（文献１７「長倉 稔、塗料・塗装技術者のための実験計画法（第２０回）、塗装技術、1998年２月号」）、化粧品などの配合率最適設計。
（Ｅ）食品加工（文献１８「相島 鐵郎、最適化手法のフレーバー開発への応用、月刊フードケミカル、1996-3」）、機械加工分野などの製造設備の最適設計、運転条件設計／調整。
（Ｆ）その他一般のシステム（プロセス分野、ビル空調、水処理、ごみ処理施設などの設備）の最適設計・運転条件設計／調整。
（Ｇ）経済分野などにおける意思決定問題（状況に応じて採るべき行動の支援など）。
（Ｈ）その他一般の最適設計問題、意思決定問題。
【００８３】
【発明の効果】
本発明によれば、従来の応答曲面法の近似関数として、グリーン関数に基づく重調和スプライン補間を適用しているので、ユーザは近似関数の次数や構造などの各種パラメータの指定をすることなく、収集したデータのみを使って非線形な最適化対象であっても所望の曲面モデルを求めることができる。これにより、より効率（労力削減、時間短縮）のよい応答曲面法による最適設計が可能となる。また、少ない時間でも精度よく所望の曲面モデルが得られるので、多項式近似のような現実的にはありえない近似を防ぐことができ、最適化対象の知識が少ない人であっても、条件設計のミスを回避することができる。さらに、重調和スプラインを用いたことにより、多次元への拡張が容易となり、多入力（多要因）の対象にも柔軟に対応することができる。また、アルゴリズムが単純であり、同一アルゴリズムで多次元化が可能であるため、プログラムがよりコンパクトにでき、さらに最適設計のための計算のパフォーマンスを向上させることができる。
【００８４】
また、出力計算時距離計算ステップ、出力計算時グリーン関数計算ステップ及び出力データ計算ステップを実行することにより、特定の入力データに対応する出力データ（補間値）を計算することができる。
【００８５】
また、第２の同定時距離計算ステップ、第２の同定時グリーン関数計算ステップ、評価指標値計算ステップ及びデータ選択ステップを実行した後に、この選択されたデータを前記収集されたデータとして、第１の同定時距離計算ステップ、第１の同定時グリーン関数計算ステップ及びモデル同定ステップを実行することにより、過去に収集したデータの中からモデル同定にとって有効なデータのみを適切に選択でき、曲面モデルを誤って同定することを回避することができる。
【図面の簡単な説明】
【図１】 本発明の第１の実施の形態となる曲面モデルの同定装置の構成を示すブロック図である。
【図２】 図１の曲面モデルの同定装置の動作を示すフローチャート図である。
【図３】 本発明の第１の実施の形態における距離計算部のモデル同定時の動作を説明するための説明図である。
【図４】 本発明の第１の実施の形態における距離計算部の出力データ計算時の動作を説明するための説明図である。
【図５】 実験データに従来の手法と本発明の第１の実施の形態の手法を適用した結果得られた曲面モデルを応答曲面として図示した図である。
【図６】 実験データに従来の手法と本発明の第１の実施の形態の手法を適用した結果得られた曲面モデルを応答曲面として図示した図である。
【図７】 実験データに従来の手法と本発明の第１の実施の形態の手法を適用した結果得られた曲面モデルを応答曲面として図示した図である。
【図８】 本発明の第２の実施の形態におけるデータスムージングの効果を示す図である。
【図９】 本発明の第３の実施の形態となる曲面モデルの同定装置の構成を示すブロック図である。
【図１０】 図９の曲面モデルの同定装置の動作を示すフローチャート図である。
【図１１】 応答曲面法における処理の流れを示すフローチャート図である。
【符号の説明】
１…距離計算部、２…グリーン関数計算部、３…モデル係数計算部、４…出力データ計算部、５…距離計算部、６…グリーン関数計算部、７…評価指標計算部、８…データ選択部。
Claims (3)
- ｎ（ｎは２以上の整数）種類の変量間の関係を近似したｎ次元の曲面モデルをコンピュータを使用して同定する曲面モデルの同定方法であって、
予めデータが収集された前記ｎ種類の変量のうちｎ−１種類の変量で定義されるｎ−１次元の空間において各データ間の距離を計算して記憶装置に記憶させる第１の同定時距離計算ステップと、
この第１の同定時距離計算ステップで計算された各距離に対応するグリーン関数値を計算して記憶装置に記憶させる第１の同定時グリーン関数計算ステップと、
この第１の同定時グリーン関数計算ステップで計算されたグリーン関数値から前記収集されたデータが前記曲面モデルの同定にとって有効なデータか否かを示す評価指標値を計算して記憶装置に記憶させる評価指標値計算ステップと、
この評価指標値が最大となるようなデータを前記収集されたデータ中から選択するデータ選択ステップと、
この選択されたデータについて、前記ｎ種類の変量のうちｎ−１種類の変量で定義されるｎ−１次元の空間において各データ間の距離を計算して記憶装置に記憶させる第２の同定時距離計算ステップと、
この第２の同定時距離計算ステップで計算された各距離に対応するグリーン関数値を計算して記憶装置に記憶させる第２の同定時グリーン関数計算ステップと、
前記選択されたデータのうち前記ｎ−１種類を除く１種類のデータと前記第２の同定時グリーン関数計算ステップで計算されたグリーン関数値とから前記ｎ次元の曲面モデルを計算して記憶装置に記憶させるモデル同定ステップとからなることを特徴とする曲面モデルの同定方法。 - 請求項１記載の曲面モデルの同定方法において、
前記ｎ種類の変量のデータは、実験計画法によって収集されることを特徴とする曲面モデルの同定方法。 - ｎ（ｎは２以上の整数）種類の変量間の関係を近似したｎ次元の曲面モデルの同定をコンピュータに実行させるプログラムであって、
予めデータが収集された前記ｎ種類の変量のうちｎ−１種類の変量で定義されるｎ−１次元の空間において各データ間の距離を計算して記憶装置に記憶させる第１の同定時距離計算ステップと、
この第１の同定時距離計算ステップで計算された各距離に対応するグリーン関数値を計算して記憶装置に記憶させる第１の同定時グリーン関数計算ステップと、
この第１の同定時グリーン関数計算ステップで計算されたグリーン関数値から前記収集されたデータが前記曲面モデルの同定にとって有効なデータか否かを示す評価指標値を計算して記憶装置に記憶させる評価指標値計算ステップと、
この評価指標値が最大となるようなデータを前記収集されたデータ中から選択するデータ選択ステップと、
この選択されたデータについて、前記ｎ種類の変量のうちｎ−１種類の変量で定義されるｎ−１次元の空間において各データ間の距離を計算して記憶装置に記憶させる第２の同定時距離計算ステップと、
この第２の同定時距離計算ステップで計算された各距離に対応するグリーン関数値を計算して記憶装置に記憶させる第２の同定時グリーン関数計算ステップと、
前記選択されたデータのうち前記ｎ−１種類を除く１種類のデータと前記第２の同定時グリーン関数計算ステップで計算されたグリーン関数値とから前記ｎ次元の曲面モデルを計算して記憶装置に記憶させるモデル同定ステップとをコンピュータに実行させることを特徴とするプログラム。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2000378625A JP3881510B2 (ja) | 2000-12-13 | 2000-12-13 | 曲面モデルの同定方法及びプログラム |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2000378625A JP3881510B2 (ja) | 2000-12-13 | 2000-12-13 | 曲面モデルの同定方法及びプログラム |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2002183111A JP2002183111A (ja) | 2002-06-28 |
| JP3881510B2 true JP3881510B2 (ja) | 2007-02-14 |
Family
ID=18847160
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2000378625A Expired - Fee Related JP3881510B2 (ja) | 2000-12-13 | 2000-12-13 | 曲面モデルの同定方法及びプログラム |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP3881510B2 (ja) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2010135370A (ja) * | 2008-12-02 | 2010-06-17 | Yamatake Corp | 検査位置決定方法、検査情報管理システム及び検査方法 |
Families Citing this family (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4535811B2 (ja) * | 2004-09-06 | 2010-09-01 | 株式会社山武 | データ処理方法及びプログラム |
| JP2006084732A (ja) * | 2004-09-15 | 2006-03-30 | Univ Of Tokyo | 多項式近似に基づく雑音下音声認識のためのモデル適応法 |
| JP4627674B2 (ja) * | 2005-03-24 | 2011-02-09 | 株式会社山武 | データ処理方法及びプログラム |
| JP4918569B2 (ja) | 2009-06-02 | 2012-04-18 | 株式会社小野測器 | 応答曲面モデル作成装置及び応答曲面モデル方法及び応答曲面モデルプログラム |
| JP6292566B2 (ja) * | 2013-08-20 | 2018-03-14 | 国立大学法人東京工業大学 | 推定方法およびそれを利用した推定装置 |
| KR101664795B1 (ko) * | 2015-10-13 | 2016-10-13 | 한국생산기술연구원 | 중간 곡면 패널 최적 설계방법, 이를 이용한 중간 곡면 패널 제조 방법 및 이를 이용한 곡면 패널 제조 방법 |
| JP6973370B2 (ja) * | 2016-02-26 | 2021-11-24 | 日本電気株式会社 | 最適化システム、最適化方法および最適化プログラム |
| US11586907B2 (en) | 2018-02-27 | 2023-02-21 | Stmicroelectronics S.R.L. | Arithmetic unit for deep learning acceleration |
| US11687762B2 (en) * | 2018-02-27 | 2023-06-27 | Stmicroelectronics S.R.L. | Acceleration unit for a deep learning engine |
| US10977854B2 (en) | 2018-02-27 | 2021-04-13 | Stmicroelectronics International N.V. | Data volume sculptor for deep learning acceleration |
| CN109902428B (zh) * | 2019-03-12 | 2023-05-02 | 武汉轻工大学 | 椭圆锥面方程的模型识别方法、装置、终端设备及可读存储介质 |
| CN110009747B (zh) * | 2019-04-11 | 2023-03-31 | 武汉轻工大学 | 单叶双曲面方程识别方法、设备、存储介质及装置 |
| JPWO2022244142A1 (ja) * | 2021-05-19 | 2022-11-24 |
-
2000
- 2000-12-13 JP JP2000378625A patent/JP3881510B2/ja not_active Expired - Fee Related
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2010135370A (ja) * | 2008-12-02 | 2010-06-17 | Yamatake Corp | 検査位置決定方法、検査情報管理システム及び検査方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| JP2002183111A (ja) | 2002-06-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP3881510B2 (ja) | 曲面モデルの同定方法及びプログラム | |
| Gálvez et al. | A new iterative mutually coupled hybrid GA–PSO approach for curve fitting in manufacturing | |
| Yu | A GP-AHP method for solving group decision-making fuzzy AHP problems | |
| Wang et al. | Adaptive MLS-HDMR metamodeling techniques for high dimensional problems | |
| EP1840841B1 (en) | Evolutionary direct manipulation of free form deformation representations for design optimisation | |
| Jemai et al. | Pyramidal hybrid approach: Wavelet network with OLS algorithm-based image classification | |
| Ebrahimi et al. | A composite iterative procedure with fast convergence rate for the progressive-iteration approximation of curves | |
| Liu et al. | Non-iterative approach for global mesh optimization | |
| Furukawa et al. | Accurate cyclic plastic analysis using a neural network material model | |
| Yang et al. | Designing fuzzy supply chain network problem by mean-risk optimization method | |
| Wan et al. | How to determine the minimum number of fuzzy rules to achieve given accuracy: a computational geometric approach to SISO case | |
| Khamdamov et al. | Algorithms of multidimensional signals processing based on cubic basis splines for information systems and processes | |
| US20080147352A1 (en) | Jacobian Scaling Parameter System and Method for Automatic Grid Generation | |
| Balzer et al. | Isogeometric finite-elements methods and variational reconstruction tasks in vision—A perfect match | |
| US6304888B1 (en) | High speed numerical integration method and system | |
| JP2001325294A (ja) | 類似画像検索方法および類似画像検索装置 | |
| CN112669127B (zh) | 用于商品推荐的方法、装置及设备 | |
| Yang et al. | Finite element mesh deformation with the skeleton-section template | |
| Kronfeld et al. | Representation of motion spaces using spline functions and fourier series | |
| Marco et al. | Structural shape optimization using Cartesian grids and automatic h-adaptive mesh projection | |
| Viúdez-Moreiras et al. | Performance comparison of Kriging and SVR surrogate models applied to the objective function prediction within aerodynamic shape optimization | |
| Du | Interactive shape design using volumetric implicit PDEs | |
| Gerzen | Analysis and applications of variational sensitivity information in structural optimisation | |
| Winter et al. | NURBS-based shape and parameter optimization of structural components with an adaptive amount of control points | |
| Beliakov | On a Method of Multivariate Density Estimate Basedon Nearest Neighbours Graphs |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20040601 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20060718 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20060725 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20060824 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20060926 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20061006 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20061107 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20061110 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 Ref document number: 3881510 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20091117 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20101117 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20101117 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111117 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121117 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20131117 Year of fee payment: 7 |
|
| LAPS | Cancellation because of no payment of annual fees |