JP7496280B2

JP7496280B2 - 関数同定方法

Info

Publication number: JP7496280B2
Application number: JP2020170109A
Authority: JP
Inventors: ▲琢▼也鈴木
Original assignee: Takenaka Corp
Current assignee: Takenaka Corp
Priority date: 2020-10-07
Filing date: 2020-10-07
Publication date: 2024-06-06
Anticipated expiration: 2040-10-07
Also published as: JP2022061876A

Description

本発明は、関数同定方法に関する。

実験などから得られた観測データから、データ間の関係を関数の形で見つけることは、物理法則を発見する上で重要とされる。このように未知の関数を同定する場合、同定対象とする関数は、様々な関数が網羅的に表せるように表現力を高く設定する必要がある。

例えば、観測データとパラメータとの線形関係を表す計画行列の算出式が既知ではない場合でも最小二乗法をパラメータ推定に適用可能とする技術がある（特許文献１参照）。また、数式空間を探索する関数同定に関する技術として遺伝的プログラミングを用いた手法が知られている（非特許文献１及び非特許文献２参照）。

特開２００７－２１８８１１号公報

伊庭斉志,佐藤泰介, "システム同定アプローチに基づく遺伝的プログラミング"，人工知能学会誌，10,4 (1994) 又吉光邦, 名嘉村盛和, 宮城隼夫： "Tree 染色体構造を持ったGAでの関数同定のための局所探索手法"， IEEJ Trans. E IS, Vol.126,No.1, 2006

以上のように、観測データから得られる入出力関係を正解データとして、入出力の関係を数式で表現したい場合、関数の表現力をどのように定めるかが問題となる。同定対象とする関数の表現力が高い場合、同定された関数に不要な項が入ってしまうことが避けられず、重要度の高い単純な式表現が見えにくくなってしまう場合がある。

例えば同定したい式ｆ（）は、様々な式を表現できる必要があるが、一方で式ｆ（）の表現力を上げすぎると無駄な項が入りやすくなるという問題がある。以下の式ｆ（）はそのような無駄な項が入っている場合の一例である。

単に、入出力システムを再現したい、という要請に対してはニューラルネットワークで対応できるもののニューラルネットワークでは入出力システムの処理の理解が困難になるという問題がある。

本発明は上記事実を考慮して、不要な項を取り除き、データ間の関係を単純化した数式を得ることを目的とする。

上記目的を達成するために、本発明の関数同定方法は、数式中の項ごとの未知のパラメータを未知係数群として、前記未知係数群の候補係数ベクトルについて、計算対象の候補係数ベクトルを用いた前記数式と正解データの入力から計算された出力と前記正解データの出力との誤差に基づいて、所定の推定手法により、当該誤差を補正するように前記未知のパラメータを修正した第１候補係数ベクトルを計算し、前記第１候補係数ベクトルを用いた前記数式と前記正解データの入力から計算された出力と前記正解データの出力との誤差に基づいて、所定の推定手法により、当該誤差に影響を与えず、かつ、不要な前記未知のパラメータを除去することによりノルムを小さくするように修正した第２候補係数ベクトルを計算し、前記第２候補係数ベクトルを次の前記計算対象の候補係数ベクトルとして、前記第１候補係数ベクトル及び前記第２候補係数ベクトルの計算を所定の条件を満たすまで繰り返すことにより、前記数式の形式を導出する、処理をコンピュータに実行させる。これにより、不要な項を取り除き、データ間の関係を単純化した数式を得ることを可能とする。

また、関数同定方法において、前記未知係数群を、少なくとも係数及び指数部の未知数とし、前記繰り返しにより得られる前記未知係数群のパラメータが所定値以上である項を計算することにより、前記数式の形式を導出するようにしてもよい。これにより、係数及び指数部の未知数のうち、不要なパラメータを除外して、データ間の関係を単純化した数式を導出することができる。

また、関数同定方法において、前記推定手法は、特異値の大きさに応じて定められたモードを用いる手法として、第１選択手法により特異値が大きな値となる低次モードを用いて前記第１候補係数ベクトルを計算し、第２選択手法により特異値が小さな値となる高次モードを用いて前記第２候補係数ベクトルを計算するようにしてもよい。これにより、特異値を基準に用いて、誤差を小さくすると共に不要な未知数を除去していくことができる。

なお，前記正解データの入力データに、前記数式に不要な入力データが含まれる場合であっても、関数同定方法による導出により、当該入力データが寄与しないことを確認することができる。

本発明によれば、不要な項を取り除き、データ間の関係を単純化した数式を得ることができる、という効果が得られる。

本発明の実施形態の関数同定装置の構成を示すブロック図である。発明の実施形態に係る関数同定装置における関数同定処理を示すフローチャートである。第１候補計算処理の詳細を示すフローチャートである。第２候補計算処理の詳細を示すフローチャートである。例題の片持ち梁のたわみ式におけるパラメータの関係を示すグラフである。例題において本提案手法によって同定された候補係数の一覧を示す図である。例題において正解データの生成にあたり、片持ち梁のたわみ式におけるパラメータを変化させる範囲を示す図である。例題における順問題、逆問題の関係を示す図である。例題における未知係数群の初期値の一例を示す図である。例題における試行回数に対する誤差及び入力ベクトルノルムの検証結果を示すグラフである。例題において、不要なパラメータを含めた場合に、同定された候補係数の一覧を示す図である。

［本発明の実施形態］
以下、図面を参照して、本発明の実施形態の関数同定手法について説明する。本実施形態の手法の特徴として、数式の関数同定を行う際に「正解データとの誤差を小さくするために未知数を調整する手順」と、「正解データとの誤差を変えずに、不要な未知数を減らしていく手順」に分離し、繰り返し実施することが挙げられる。これにより、観測データに対して最も単純化された数式の同定が可能となる。

以下では、本実施形態の関数同定手法に係る構成及び作用を説明し、例題の検証について述べる。なお、誤差修正法、及びノイズ除去法として、ＭＩＥＣ（Modal Iterative Error Correction）法を用いる場合について説明するが、他の誤差修正法（非線形最小二乗法、Newton法、最急降下法、ＧＡ（Genetic Algorithm））、又はノイズ除去法（最小二乗法等）を用いてもよい。

（構成）
図１は、関数同定装置１００の構成を示すブロック図である。図１に示すように、関数同定装置１００は、設定部１１０と、第１候補計算部１２０と、第２候補計算部１３０と、判定部１４０と、同定部１５０とを含む。関数同定装置１００は、ＣＰＵと、ＲＡＭと、各処理部を実行するためのプログラム及び各種データを記憶したＲＯＭと、を含むコンピュータで構成することができる（図示省略）。なお、具体的な処理については作用の説明で述べる。

設定部１１０は、数式を定義する。数式は、項ごとの指数及び係数等を未知として、これらの未知のパラメータからなる未知係数群を含む形で定義する。未知係数群は候補係数として候補係数ベクトルとして表すことができる。設定部１１０は、未知数の初期値を設定する。また、設定部１１０は、正解データを設定する。正解データは観測等によって得られた入出力データの各々である。なお、未知係数群は指数及び係数だけでなく、任意のべき乗部（ｓｉｎ値等）を未知数としてもよい。

候補係数ベクトルは、以下、第１候補計算部１２０での計算結果を第１候補係数ベクトル、第２候補計算部１３０での計算結果を第２候補係数ベクトルと呼ぶ。

第１候補計算部１２０は、計算対象の候補係数ベクトルを用いた数式と正解データの入力から計算された出力と正解データの出力との誤差に基づいて、第１候補係数ベクトルを計算する。第１候補計算部１２０では、誤差を補正するように計算対象の候補係数ベクトルを修正する。

第２候補計算部１３０は、第１候補計算部１２０で計算された第１候補係数ベクトルを用いた前記数式と正解データの入力から計算された出力と正解データの出力との誤差の偏微分に基づいて、第１候補係数ベクトルを修正した第２候補係数ベクトルを計算する。なお、ＭＩＥＣ法以外の手法では誤差の偏微分に限定されるものではない。第２候補計算部１３０では、誤差に影響を与えず、かつ、不要な未知のパラメータを除去することによりノルムを小さくするように第１候補係数ベクトルを修正する。ここで除去するとは、繰り返しの処理によって該当するパラメータの成分を減らしていくということである。

判定部１４０は、第１候補計算部１２０及び第２候補計算部１３０の計算を所定の条件を満たすまで繰り返す。条件は指定回数とすればよい。条件は閾値等としてもよい。繰り返しでは、第２候補係数ベクトルを次の計算対象の候補係数ベクトルとする。

同定部１５０は、最終的に得られた候補係数ベクトルを用いて、未知数が所定値以上である項を抽出することにより未知係数群のうち微小なパラメータを除外し、数式の形式を導出する。

（作用）
次に、本発明の実施形態の関数同定装置１００の作用について説明する。図２は、本発明の実施形態に係る関数同定装置１００における関数同定処理を示すフローチャートである。なお、関数同定装置１００としてのコンピュータはＣＰＵがＲＯＭからプログラム及び各種データを読み出して実行することにより、各種処理を行い、ＣＰＵが、コンピュータの各部として機能する。

ステップＳ１００では、設定部１１０が、未知係数群を含む数式を定義する。

ステップＳ１０２では、設定部１１０が、未知係数群の初期値を設定する。

ステップＳ１０４では、設定部１１０が、定義した数式の正解データを設定する。

ステップＳ１０６では、第１候補計算部１２０が、計算対象の候補係数ベクトルを用いた数式と正解データの入力から計算された出力と正解データの出力との誤差に基づいて、第１候補係数ベクトルを計算する。具体的な処理内容については後述する。

ステップＳ１０８では、第２候補計算部１３０が、第１候補計算部１２０で計算された第１候補係数ベクトルを用いた前記数式と正解データの入力から計算された出力と正解データの出力との誤差の偏微分に基づいて、第１候補係数ベクトルを修正した第２候補係数ベクトルを計算する。具体的な処理内容については後述する。

ステップＳ１１０では、判定部１４０が、指定回数まで繰り返したか否かを判定する。指定回数まで繰り返したと判定した場合にはステップＳ１１４へ移行し、指定回数まで繰り返していないと判定した場合にはステップＳ１１２へ移行する。

ステップＳ１１２では、判定部１４０が、処理対象の候補係数ベクトルに第２候補係数ベクトルを設定し、ステップＳ１０６に戻って処理を繰り返す。

ステップＳ１１４では、同定部１５０が、最終的に得られた候補係数ベクトル（第２候補係数ベクトル）のうち、所定値以上である成分のみを用いて数式を設定することにより、未知係数群の微小なパラメータを除外し、数式の形式を導出する。

次に、ステップＳ１０６の第１候補計算処理の詳細について図３のフローチャートを用いて説明する。以下は第１候補計算部１２０が行う第１候補計算処理である。第１候補計算処理では、誤差を小さくするように候補係数を修正する。以下は誤差修正法による第１候補計算処理の一例である。

ステップＳ２００では、現在の処理対象の候補係数ベクトルを用いた数式と正解データの入力から計算される出力と、正解データの出力との出力データの誤差を全正解データ件数分計算し、誤差ベクトルを計算する。

ステップＳ２０２では、現在の処理対象の候補係数ベクトルの値に対する誤差ベクトルの偏微分マトリクス（ヤコビアンマトリクス）を作成する。この偏微分の計算方法は、解析微分でも、数値微分でもよい。

ステップＳ２０４では、偏微分マトリクスの特異値分解を実行する。

ステップＳ２０６では、特異値分解結果から低次モードを抽出する。低次モードとは、特異値が大きなモードである。特異値が大きいモードとは、未知係数群の少し変動によって、誤差が大きく変わるモードである。逆に特異値が小さいモードとは、未知係数群の変動が、誤差の変化への影響が小さいモードである。第１選択手法による低次モードの選択法は、（１）予め設定した特異値の閾値（つまり閾値以上の値）から選択する方法、（２）予め設定したモード数まで選択する方法、（３）現在の出力ベクトル（誤差ベクトル）のノルム比による方法、あるいはこれらを組み合わせた方法などユーザーが適用する問題に合わせて設定する。モードの選択手法に関しては、例えば参考文献１に記載の手法を参照し、実施することができる。参考文献１には、差分ベクトルのノルムを用いる手法、特異値がモードに関する閾値ＴＨ以上であるモードを選択する手法、及び繰り返し回数に応じて徐々に閾値を下げる手法等が記載されている。
［参考文献１］特開２０１９－１９４８２７号公報

ステップＳ２０８では、特異値分解結果から、抽出した低次モードに対応する部分行列を用いて逆行列を作成する。

ステップＳ２１０では、作成された逆行列に誤差ベクトルを乗じて、修正係数ベクトルを計算する。

ステップＳ２１２では、計算した修正係数ベクトルを、現在の処理対象の候補係数ベクトルに加算し、第１候補係数ベクトルを計算する。

次に、ステップＳ１０８の第２候補計算処理の詳細について図４のフローチャートを用いて説明する。以下は第２候補計算部１３０が行う第２候補計算処理である。第２候補計算処理では、誤差を変えずに、候補係数ベクトルを小さくする。以下はノイズ除去法による第２候補計算処理の一例である。

ステップＳ３００では、現在の処理対象の第１候補係数ベクトルを用いた数式と正解データの入力から計算される出力と、正解データの出力との出力データの誤差を全正解データ分計算し、誤差ベクトルとする。

ステップＳ３０２では、現在の処理対象の第１候補係数ベクトルの値に対する誤差ベクトルの偏微分マトリクス（ヤコビアンマトリクス）を作成する。この偏微分の計算方法は、解析微分でも、数値微分でもよい。

ステップＳ３０４では、偏微分マトリクスの特異値分解を実行する。

ステップＳ３０６では、特異値分解結果から高次モードを抽出する。高次モードとは、特異値が小さな値を選択するモードである。第２選択手法による高次モードの選択法は、（１）予め設定した特異値の閾値（つまり閾値以下の値）から選択する方法、（２）閾値を変えつつ予め設定したモード数まで選択する方法、（３）現在の入力ベクトルのノルム比による方法、あるいはこれらを組み合わせた方法などユーザーが適用する問題に合わせて設定する。第２選択手法も第１選択手法と同様に、参考文献１に記載の手法を参照する。

ステップＳ３０８では、特異値分解結果から、抽出した高次モードに対応する部分行列を作成する。

ステップＳ３１０では、部分行列を用いて、現在の処理対象の第１候補係数ベクトルに含まれる、これらの抽出された高次モードの各成分量を計算する。

ステップＳ３１２では、高次モードの各成分に除去率ηを乗じて、除去すべき候補係数ベクトルを計算し、現在の処理対象の第１候補係数ベクトルの入力パラメータから除去すべき修正係数ベクトルを計算する。

ステップＳ３１４では、計算した修正係数ベクトルを、現在の処理対象の第１候補係数ベクトルから減算し、第２候補係数ベクトルを計算する。

以上、説明したように、本発明の実施形態に係る関数同定装置１００によれば、不要な項を取り除き、データ間の関係を単純化した数式を得ることができる。

（例題の検証）
関数同定装置１００の関数同定処理の効果を、例題に適用して検証して確認を行った。例題の数式として、図５及び以下（１）式に示す片持ち梁のたわみδの算定式に適用した場合について説明する。図５は片持ち梁のたわみ式におけるパラメータの関係を示すグラフである。（１）式が未知数係数群から同定したいたわみ式である。

・・・（１）

上記たわみ式を同定する問題設定に対して、設定部１１０で定義する数式を以下（２）式で表す。（２）式は６つの級数で構成し、各項は６つの全変数のべき乗の積で構成する。なお級数の項数は一例である。変数は観測データのパラメータであるＥ、Ｇ、Ｌ、Ｑ、Ａｓ、及びＩである。Ｅはヤング率、Ｇはせん断弾性係数、Ｌはスパン長、Ｑはせん断力、Ａｓはせん断断面積、Ｉは断面２次モーメントである。また、各項の係数及び指数部を未知数（未知係数群）とする。これにより全部で未知係数群は４２個となる。

・・・（２）

（２）式に対して、別途作成した正解データ（３００の正解となる入出力データ）を教師データとして、本実施形態の提案手法を用いて関数同定を行った。適用の具体例については後述する。

図６に本提案手法によって同定された候補係数の一覧を示す。プラスのパラメータが分子、マイナスのパラメータが分母に対応する。べき乗数が０．１未満の項を無視すれば（つまり、べき乗数が０．１以上の項を計算対象とすれば）、残った第３項および第４項は正解となる（１）式のたわみ式と一致していることが確認できる。以上のように本実施形態の関数同定手法の効果が確認することができた。ここでべき乗数が０．１未満の項を考慮しないことが、所定値以上である項を計算することの一例である。なお、０．１以上とするのは一例であり、適宜、適切な値を設定すればよい。

適用の具体例について説明する。まず適用した正解データの生成手法について説明する。正解データは、図７に示した範囲で観測データのパラメータをランダムに変化させ、理論式を使って「Ｅ，Ｇ，Ｌ，Ｑ，Ａｓ，Ｉのときの正解δ」を３００個計算し、正解データとした。３００個の正解データについては、計算に用いたパラメータが本開示の「正解データの入力」の一例であり、正解δが本開示の「正解データの出力」の一例である。３００個の正解データの各々について、候補係数ベクトル（第１候補係数ベクトル）を用いた数式と当該正解データの入力から計算される出力と、当該正解データの出力との誤差を求める。そして、図８に示すような逆問題を解くことで、この正解データとの誤差を減らすような、関数の未知数（係数及び指数）を求めた。順解析を、４２次元の「式の未知係数群」から３００次元の「全テストデータの誤差」の計算とし、この逆解析により逆問題を解くことを行った。逆問題は上記図３～５に示したフローチャートの手法により解くことができる。

その他の解析条件について説明する。未知係数群の初期値は図９に示すように設定した。また、第１候補係数ベクトル計算時に用いるモード選択手法はノルム比による方法を用い、その閾値は０．２とし、繰返し回数を２０００、摂動量を０．０１とする。第２候補係数ベクトル計算時に用いるモード選択手法は特異値閾値による方法を用い、その閾値は１．０Ｅ－５、除去率ηのクリーニング比率は１０％とした。また試行回数に対する誤差及び入力ベクトルノルムの検証結果を示すグラフを図１０に示す。この検証結果ではＩｔｅｒ（試行回数）１８７４の結果では、Ｅｒｒｏｒ＿ｎｏｒｍ（誤差）＝０．００７２となったことが確認できた。

以上のようにたわみ式のような構造計算式を求めたい数式とし、未知係数群は、構造計算式に用いられ得る観測データのパラメータの未知係数群として、本実施形態の手法を適用できる。繰り返しにおいて、観測データのパラメータのうち構造計算式に寄与しない観測データのパラメータを未知係数群の計算結果をもとに除去することにより、数式の形式を同定する。

なお、不要な変数を追加した場合について説明する。例えば、観測データのパラメータに温度Ｔという変数があったとして、このパラメータが関数に影響を与えるかが不明な場合は、この変数も数式では定義され得る。以下（３）式に温度Ｔを定義する数式に含む場合を示す。

・・・（３）

温度Ｔのパラメータのとり得る範囲は０～６００として正解データを生成した。その他の解析条件は上記同様とした。この場合、図１１に同定された候補係数の一覧を示す。図１１に示すように温度Ｔの変数は０になっており、たわみ式において意味のない温度Ｔについては同定される数式には入らないように求められることが確認できた。

なお、本発明は、上述した実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。例えば、正解データに誤差がある場合において、正解データを、正解データを取得した際の環境等の要因を示す外部情報を用いて補正することも可能である。

１００関数同定装置
１１０設定部
１２０第１候補計算部
１３０第２候補計算部
１４０判定部
１５０同定部

Claims

数式中の項ごとの未知のパラメータを未知係数群として、前記未知係数群の候補係数ベクトルについて、計算対象の候補係数ベクトルを用いた前記数式と正解データの入力から計算される出力と前記正解データの出力との誤差に基づいて、所定の推定手法により、当該誤差を補正するように前記未知のパラメータを修正した第１候補係数ベクトルを計算し、
前記第１候補係数ベクトルを用いた前記数式と前記正解データの入力から計算される出力と前記正解データの出力との誤差に基づいて、所定の推定手法により、当該誤差に影響を与えず、かつ、不要な前記未知のパラメータを除去することによりノルムを小さくするように前記第１候補係数ベクトルを修正した第２候補係数ベクトルを計算し、
前記第２候補係数ベクトルを次の前記計算対象の候補係数ベクトルとして、前記第１候補係数ベクトル及び前記第２候補係数ベクトルの計算を所定の条件を満たすまで繰り返すことにより、前記数式の形式を導出する、
処理をコンピュータに実行させる関数同定方法。
前記未知係数群を、少なくとも係数及び指数部の未知数とし、
前記繰り返しにより得られる前記未知係数群のパラメータが所定値以上である項を計算することにより、前記数式の形式を導出する請求項１に記載の関数同定方法。
前記推定手法は、特異値の大きさに応じて定められたモードを用いる手法として、
第１選択手法においては特異値が大きな値を選択する低次モードを用いて前記第１候補係数ベクトルを計算し、第２選択手法においては特異値が小さな値を選択する高次モードを用いて前記第２候補係数ベクトルを計算する請求項１又は請求項２に記載の関数同定方法。