(57)【要約】【課題】微分方程式で表されるシステムの同定に要する計算量を減らす。【解決手段】本係数決定方法は、システムを表す微分方程式の次数を第１の次数に設定した場合の第１の微分係数を、システムへのディジタル入力及びシステムからのディジタル出力に対する離散フーリエ変換の結果に基づき算出し、微分方程式の次数を第２の次数に設定した場合の第２の微分係数を記憶装置から読み出し、第１の微分係数と読み出した第２の微分係数との比較に基づき、微分方程式の微分係数の収束について判定し、判定の結果に基づき、微分方程式の微分係数を決定する処理を含む。【選択図】図４
2018092218.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータに、
システムを表す微分方程式の次数を第１の次数に設定した場合の第１の微分係数を、前記システムへのディジタル入力及び前記システムからのディジタル出力に対する離散フーリエ変換の結果に基づき算出し、
前記微分方程式の次数を第２の次数に設定した場合の第２の微分係数を記憶装置から読み出し、前記第１の微分係数と読み出した前記第２の微分係数との比較に基づき、前記微分方程式の微分係数の収束について判定し、
判定の結果に基づき、前記微分方程式の微分係数を決定する、
処理を実行させる係数決定プログラム。
【請求項２】
前記第１の次数が前記第２の次数より大きい、
請求項１記載の係数決定プログラム。
【請求項３】
前記コンピュータに、
前記微分方程式に複数の微分演算子が含まれる場合、前記複数の微分演算子の次数の組合せを、前記複数の微分演算子の次数の総和が前記第２の次数以上になり、且つ、前記微分方程式の次数を前記第２の次数に設定した場合における前記複数の微分演算子の次数の組合せと異なるように設定する、
処理をさらに実行させる請求項１記載の係数決定プログラム。
【請求項４】
前記微分方程式の微分係数の収束について判定を行う処理において、
前記第１の微分係数及び前記第２の微分係数のいずれにも含まれる微分係数が、前記微分方程式の次数を前記第２の実数に設定した場合と前記微分方程式の次数を前記第２の次数に設定した場合とで同じ値であり、且つ、前記第１の微分係数に含まれ且つ前記第２の微分係数に含まれない微分係数の値がゼロであるか判定する、
請求項１記載の係数決定プログラム。
【請求項５】
前記微分方程式の微分係数を決定する処理において、
前記判定の結果が前記微分方程式の微分係数が収束したことを示している場合、前記微分方程式の微分係数を決定し、前記判定の結果が前記部分方程式の微分係数が収束していないことを示している場合、算出した前記微分係数を前記記憶装置に格納し且つ前記微分方程式の次数を第３の次数に設定する、
請求項１記載の係数決定プログラム。
【請求項６】
コンピュータが、
システムを表す微分方程式の次数を第１の次数に設定した場合の第１の微分係数を、前記システムへのディジタル入力及び前記システムからのディジタル出力に対する離散フーリエ変換の結果に基づき算出し、
前記微分方程式の次数を第２の次数に設定した場合の第２の微分係数を記憶装置から読み出し、前記第１の微分係数と読み出した前記第２の微分係数との比較に基づき、前記微分方程式の微分係数の収束について判定し、
判定の結果に基づき、前記微分方程式の微分係数を決定する、
処理を実行する係数決定方法。
【請求項７】
システムを表す微分方程式の次数を第１の次数に設定した場合の第１の微分係数を、前記システムへのディジタル入力及び前記システムからのディジタル出力に対する離散フーリエ変換の結果に基づき算出する算出部と、
前記微分方程式の次数を第２の次数に設定した場合の第２の微分係数を記憶装置から読み出し、前記第１の微分係数と読み出した前記第２の微分係数との比較に基づき、前記微分方程式の微分係数の収束について判定し、判定の結果に基づき、前記微分方程式の微分係数を決定する決定部と、
を有する係数決定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、システム同定技術に関する。
【背景技術】
【０００２】
システムへの入力とシステムからの出力とに基づいてシステムを同定する技術が存在する。例えば図１に示すように、常微分方程式で表される物理的なシステムについてアナログの入力ｕ（ｔ）及び出力ｙ（ｔ）が与えられたとき、以下に示すような関係が成り立つように微分演算子Ｐ（ｄ／ｄｔ）の微分係数Ｐ0乃至Ｐmと微分演算子Ｑ（ｄ／ｄｔ）の微分係数Ｑ0乃至Ｑl-1とを決定する問題を解くことで同定が行われる。なお、本願では、Ｐ（ｄ／ｄｔ）ｙ（ｔ）とＱ（ｄ／ｄｔ）ｕ（ｔ）とがおおよそ一致する場合も「＝」で表すものとする。
【０００３】
【数１】
2018092218.tif 000003
【０００４】
システムへの入力がアナログデータではなくディジタルデータ（すなわち離散データ）であることもある。この場合、例えば図２に示すように、常微分方程式で表される物理的なシステムについてディジタルの入力ｕ0、ｕ1...ｕ2nと出力ｙ0、ｙ1...ｙ2nとが与えられたとき、（１）式のような関係が成り立つように微分演算子Ｐ（ｄ／ｄｔ）の微分係数Ｐ0乃至Ｐmと微分演算子Ｑ（ｄ／ｄｔ）の微分係数Ｑ0乃至Ｑl-1とを決定する問題を解くことで同定が行われる。
【０００５】
離散データを用いたシステム同定に関して、或る文献は、以下のような技術を開示する。具体的には、情報処理装置が、システムのディジタル入力及びディジタル出力のそれぞれについて、指数関数の線型和で表されるフィッティング曲線を生成する。そして、情報処理装置は、微分演算子とディジタル入力についてのフィッティング曲線との積の結果に含まれる各指数関数の係数と、ディジタル出力についてのフィッティング曲線に含まれる同一の指数関数の係数とが一致するように、微分演算子の係数を算出する。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２０１５−１３５６６２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００７】
ところで、人の運動や車両の移動などについての力学的な数理モデルを同定する場合、計測される離散データ（例えば、加速度の離散データ）に含まれる高周波成分はノイズに関連するものであるため、低周波成分を用いて同定が行われる。すなわち、数理モデルは低階微分の項のみを含む微分方程式で表される。
【０００８】
上記技術は汎用的な技術であり、微分方程式の次数の設定について決まりが無いので、離散データの全周波数成分を使用して解析が行われる。従って、上記のような数理モデルの同定には必ずしも適していない。また、微分方程式の次数を変えて同定を行うたびに同定結果について検証が行われるので、検証が完了するまでの計算量が多くなることがある。
【０００９】
本発明の目的は、一側面によれば、微分方程式で表されるシステムの同定に要する計算量を減らす技術を提供することである。
【課題を解決するための手段】
【００１０】
一態様に係る係数決定方法は、システムを表す微分方程式の次数を第１の次数に設定した場合の第１の微分係数を、システムへのディジタル入力及びシステムからのディジタル出力に対する離散フーリエ変換の結果に基づき算出し、微分方程式の次数を第２の次数に設定した場合の第２の微分係数を記憶装置から読み出し、第１の微分係数と読み出した第２の微分係数との比較に基づき、微分方程式の微分係数の収束について判定し、判定の結果に基づき、微分方程式の微分係数を決定する処理を含む。
【発明の効果】
【００１１】
微分方程式で表されるシステムの同定に要する計算量を減らせるようになる。
【図面の簡単な説明】
【００１２】
【図１】図１は、アナログ入出力のシステムを模式的に示す図である。
【図２】図２は、ディジタル入出力のシステムを模式的に示す図である。
【図３】図３は、本実施の形態に係るシステムの概要を示す図である。
【図４】図４は、本実施の形態に係る情報処理装置の機能ブロック図である。
【図５】図５は、メインの処理フローを示す図である。
【図６】図６は、ディジタル入出力のシステムを模式的に示す図である。
【図７】図７は、第１データ格納部に格納されるデータの一例を示す図である。
【図８】図８は、第１データ格納部に格納されるデータの一例を示す図である。
【図９】図９は、第２データ格納部に格納されるデータの一例を示す図である。
【図１０】図１０は、第２データ格納部に格納されるデータの一例を示す図である。
【図１１】図１１は、線形回帰を説明するための図である。
【図１２】図１２は、ｍの値及びｌの値の更新について説明するための図である。
【図１３】図１３は、微分係数の算出結果を示す図である。
【図１４】図１４は、予測処理の処理フローを示す図である。
【図１５】図１５は、コンピュータの機能ブロック図である。
【発明を実施するための形態】
【００１３】
図３に、本実施の形態に係るシステムの一例を示す。車両１０００には、検出装置１００が設けられており、当該検出装置１００は、車両１０００の加速度を測定する加速度センサ１０１と、車両１０００の走行距離を測定する距離センサ１０２と、通信部１０３とを有する。なお、車両１０００は、１台だけではなく、複数台の場合もある。
【００１４】
検出装置１００の通信部１０３は、インターネットなどのネットワーク２００を介して、例えばデータセンタやクラウド内の情報処理装置３００へ測定データを送信する。車両１０００は、ｍ次（ｍは０以上の整数）の微分演算子とｌ次（ｌは０以上の整数）の微分演算子とを含む以下の微分方程式で表される物理的なシステムである。
【００１５】
【数２】
2018092218.tif 000004
【００１６】
図４に、情報処理装置３００の構成例を示す。
【００１７】
情報処理装置３００は、通信部３１０と、入出力データ格納部３２０と、学習処理部３３０と、予測処理部３４０とを有する。
【００１８】
通信部３１０は、車両１０００における検出装置１００からの測定データを受信し、入出力データ格納部３２０に格納する。学習処理部３３０は、ディジタル入力及びディジタル出力から、車両１０００を表す微分方程式を同定するための処理を実行する。学習処理部３３０は、ＤＦＴ（Discrete Fourier Transformation）計算部３３１と、第１データ格納部３３２と、第１演算部３３３と、第２データ格納部３３４と、更新部３３５と、比較部３３６とを有する。第１演算部３３３は、第１係数算出部３３３１と、線型回帰処理部３３３２とを有する。
【００１９】
ＤＦＴ計算部３３１は、入出力データ格納部３２０に格納されているデータに基づき処理を実行し、処理結果を第１データ格納部３３２に格納する。第１係数算出部３３３１は、第１データ格納部３３２に格納されているデータに基づき、線型回帰処理部３３３２で用いられるデータを生成し、第１データ格納部３３２に格納する。線型回帰処理部３３３２は、第１データ格納部３３２に格納されているデータを用いて線型回帰処理を実行し、処理結果を第２データ格納部３３４に格納する。比較部３３６は、第２データ格納部３３４に格納されているデータに基づき比較処理を行い、処理結果に応じて、第２演算部３４３又は更新部３３５に通知を行う。更新部３３５は、ｍの値およびｌの値を更新し、更新後のｍの値およびｌの値をＤＦＴ計算部３３１に通知する。
【００２０】
線型回帰処理の処理結果は、ｍ次の微分演算子の微分係数およびｌ次の微分演算子の微分係数を含む。ｍ次の微分演算子の微分係数およびｌ次の微分演算子の微分係数が得られれば、車両１０００を表す微分方程式が得られることになる。
【００２１】
予測処理部３４０は、同定された微分方程式及び新たなディジタル入力からディジタル出力を予測する処理を実行する。予測処理部３４０は、フィッティング曲線生成部３４１と、第３データ格納部３４２と、第２演算部３４３と、第４データ格納部３４４とを有する。第２演算部３４３は、第２係数算出部３４３１と、出力データ算出部３４３２とを有する。
【００２２】
フィッティング曲線生成部３４１は、離散フーリエ変換（ＤＦＴ）を行って、入出力データ格納部３２０に格納されたディジタル入力についてのフィッティング曲線を生成し、フィッティング曲線のデータを第３データ格納部３４２に格納する。第２係数算出部３４３１は、比較部３３６から通知された微分係数と、第３データ格納部３４２に格納されているフィッティング曲線のデータとから、出力関数の係数を算出する。出力データ算出部３４３２は、第２係数算出部３４３１により算出された係数で特定される出力関数により、出力データを算出し、第４データ格納部３４４に格納する。
【００２３】
次に、図５乃至図９を用いて、情報処理装置３００が実行する処理について説明する。なお、車両１０００の検出装置１００の通信部１０３は、任意のタイミングで、ディジタル入力（例えば加速度センサ１０１の測定データ）及びディジタル出力（例えば距離センサ１０２の測定データ）のデータを、情報処理装置３００に送信するものとする。情報処理装置３００の通信部３１０は、車両１０００からの受信データを、入出力データ格納部３２０に格納するものとする。
【００２４】
より具体的には、図６に示すように、ディジタル入力Ｕj（ｊ＝１乃至ｒ）に対して、システムの出力Ｙjが得られるものとする。なお、Ｕj及びＹjは、以下に示すように、２ｎ＋１個の離散点を含む。
【００２５】
【数３】
2018092218.tif 000005
【００２６】
このシステムは、上でも述べたように、ｍ次の微分演算子Ｐ（ｄ／ｄｔ）とｌ次の微分演算子Ｑ（ｄ／ｄｔ）とを含む微分方程式Ｐ（ｄ／ｄｔ）ｙj（ｔ）＝Ｑ（ｄ／ｄｔ）ｕj（ｔ）で表される。
【００２７】
そうすると、学習処理部３３０のＤＦＴ計算部３３１は、ｍ＝０と設定し且つｌ＝０と設定し（図５：ステップＳ１）、さらにＤＦＴ計算部３３１は、ｐ＝ｍ＋ｌと設定する（ステップＳ３）。
【００２８】
ＤＦＴ計算部３３１は、係数Ｙj,p、Ｕj,p、Ｙj,-p及びＵj,-pが既に算出されたか判定する（ステップＳ５）。
【００２９】
係数Ｙj,p、Ｕj,p、Ｙj,-p及びＵj,-pが既に算出された場合（ステップＳ５：Ｙｅｓルート）、処理はステップＳ９に移行する。一方、Ｙj,p、Ｕj,p、Ｙj,-p及びＵj,-pが未だ算出されていない場合（ステップＳ５：Ｎｏルート）、ＤＦＴ計算部３３１は、入出力データ格納部３２０に格納されているデータに基づき、係数Ｙj,p、Ｕj,p、Ｙj,-p及びＵj,-pを算出する（ステップＳ７）。ＤＦＴ計算部３３１は、算出結果を第１データ格納部３３２に格納する。
【００３０】
図７に、ステップＳ７の処理によって第１データ格納部３３２に格納されるデータの一例を示す。図７には、ｐ＝０の場合にステップＳ７の処理によって第１データ格納部３３２に格納される算出結果が示されている。なお、例えば図８に示すように、ステップＳ７の処理が行われる度に第１データ格納部３３２に処理結果が追加されていくようになる。
【００３１】
ステップＳ７においては、後述する（１４）式によって係数Ｙj,p、Ｕj,p、Ｙj,-p及びＵj,-pが算出されるが、この（１４）式の導出過程等について以下で説明する。なお、ここでは説明を分かりやすくするためｊの記載が一部で省略されている。
【００３２】
離散点の実数値ｙ0，ｙ1，．．．，ｙ2nに対してＤＦＴを行って、ｎ＋１個の周波数を算出すると、以下のようになる。
【００３３】
【数４】
2018092218.tif 000006
【００３４】
各周波数Ｙ0，Ｙ1，．．．，Ｙnは、複素数になるので、実部ａq及び虚部ｂqにより以下のように表される。
【００３５】
【数５】
2018092218.tif 000007
【００３６】
なお、離散点は実数値であるのでＤＦＴの性質から、ｂ0＝０となることに注意する。
【００３７】
このように得られたａq及びｂqを用いて、以下に示す有限の（すなわちｎ次の）フーリエ級数を生成することができる。
【００３８】
【数６】
2018092218.tif 000008
【００３９】
ここで得られたフーリエ級数が離散点ｙ0，ｙ1，．．．，ｙnのフィッティング曲線となる。このフィッティング曲線は、以下で説明するように、これらの離散点を通過する。
【００４０】
なお、２ｎ＋１個の離散点ｙ0，ｙ1，．．．，ｙ2nが得られると、ＤＦＴによって２ｎ＋１個の周波数を得ることができる。
【００４１】
【数７】
2018092218.tif 000009
【００４２】
この離散点が実数値である場合には、以下のような関係式が得られる。
【００４３】
【数８】
2018092218.tif 000010
【００４４】
バー付きのＹは、Ｙの共役数を表す。このことから、以下の関係式も得られる。
【００４５】
【数９】
2018092218.tif 000011
【００４６】
ここで周波数Ｙ0，Ｙ1，．．．，Ｙ2nから元の離散点のデータを復元するために逆離散フーリエ変換（ＩＤＦＴ：Inverse DFT）を行う場合には、以下の数式の１行目のように表される。
【００４７】
【数１０】
2018092218.tif 000012
【００４８】
この２行目は、上で示したＹとＹの共役数との関係から、１行目の右辺を分解したものに相当する。３行目以降は、パラメータｑ等を調整して、周波数の実部ａq及び虚部ｂqで表現して整理を進めると、最終行が得られる。
【００４９】
一方、（６）式で示したｎ次のフーリエ級数ｙ（ｔ）をｔ＝０，１，．．．，２ｎでサンプリングすると、以下の式が得られる。
【００５０】
【数１１】
2018092218.tif 000013
【００５１】
この（１１）式の右辺と逆離散フーリエ変換の式の最終行は等しいので、これによってｙi＝ｙ（ｉ）（ｉ＝０，１，．．．，２ｎ）となることが分かる。すなわち、（６）式のｎ次のフーリエ級数ｙ（ｔ）は、元の離散点の各々を通過することが分かる。
【００５２】
また、任意のｎ次のフーリエ級数ｇ（ｔ）が与えられたとき、ｇ（ｔ）をｔ＝０，１，．．．，２ｎでサンプリングして離散点ｇi＝ｇ（ｉ）（ｉ＝０，１，．．．，２ｎ）を求めてから、この離散点ｇi（ｉ＝０，１，．．．，２ｎ）のフィッティング曲線を求めると、元のｇ（ｔ）に戻る。これは上で述べたのと同様の手法を用いれば証明できるが、ここではその詳細については省略する。以上のことから、このフィッティング法は、離散点と有限のフーリエ級数との間での一対一且つ上の対応になっている。
【００５３】
従って、ディジタル出力についてのフィッティング曲線ｙj（ｔ）及びディジタル入力についてのフィッティング曲線ｕj（ｔ）は、以下のように得られる。
【００５４】
【数１２】
2018092218.tif 000014
【００５５】
但し、このフィッティング曲線ｙj（ｔ）及びｕj（ｔ）は、ネイピア数ｅを底とする指数関数の線型和で表すことができる。このように指数関数の線型和で表される場合には、ｍ次の微分を簡単に行うことができる。
【００５６】
【数１３】
2018092218.tif 000015
【００５７】
このように、係数Ｙj,q及びＵj,qは、ネイピア数ｅを底とする指数関数の線型和でフィッティング曲線を表す場合に使用される係数に相当する。係数Ｙj,q及びＵj,qは以下のように表される。
【００５８】
【数１４】
2018092218.tif 000016
【００５９】
図５の説明に戻り、第１演算部３３３は、ｑ＝−ｐ，−ｐ＋１，．．．，ｐについて、第１データ格納部３３２に格納されている係数Ｙj,q、Ｕj,q、Ｙj,-q及びＵj,-qを読み出し、読み出した係数を用いてｍ次の微分演算子Ｐ（ｄ／ｄｔ）の微分係数Ｐ0乃至Ｐm及びｌ次の微分演算子Ｑ（ｄ／ｄｔ）の微分係数Ｑ1乃至Ｑlを算出する（ステップＳ９）。第１演算部３３３は、算出した微分係数を第２データ格納部３３４に格納する。
【００６０】
図９に、第２データ格納部３３４に格納されるデータの一例を示す。図９の例では、ｍ＝０且つｌ＝０の場合にステップＳ９の処理によって第２データ格納部３３４に格納される算出結果が示されている。なお、例えば図１０に示すように、ステップＳ９の処理が行われる度に第２データ格納部３３４に処理結果が追加されていくようになる。
【００６１】
ステップＳ９の処理について説明を追加する。ステップＳ９においては、まず第１係数算出部３３３１が、係数Ｙj,q及びＵj,qから、ｍ次の微分演算子の各微分係数およびｌ次の微分演算子の各微分係数を算出するためのαβ平面の点群（Ｃ，Ｄ）を複数生成する。第１係数算出部３３３１は、処理結果を第１データ格納部３３２に格納する。
【００６２】
上でも述べたように、Ｐ（ｄ／ｄｔ）ｙj（ｔ）＝Ｑ（ｄ／ｄｔ）ｕj（ｔ）であるから、左辺のＰ（ｄ／ｄｔ）ｙj（ｔ）は、以下のように表される。
【００６３】
【数１５】
2018092218.tif 000017
【００６４】
ｙjは、上でも述べたように、指数関数の線型和となるので、１階微分であれば係数が、指数関数の指数部の係数となり、２階微分であれば係数が（指数関数の指数部の係数）2となり、３階微分であれば係数が（指数関数の指数部の係数）3となり、ｍ階微分であれば係数が（指数関数の指数部の係数）mとなる。このように微分演算が非常に簡単になる。
【００６５】
すなわち、以下のように表される。
【００６６】
【数１６】
2018092218.tif 000018
【００６７】
従って、上で述べた（１５）式が得られる。
【００６８】
同様に、Ｑ（ｄ／ｄｔ）も以下のように表される。
【００６９】
【数１７】
2018092218.tif 000019
【００７０】
そうすると、Ｐ（ｄ／ｄｔ）ｙj（ｔ）＝Ｑ（ｄ／ｄｔ）ｕj（ｔ）は、以下のように表される。
【００７１】
【数１８】
2018092218.tif 000020
【００７２】
なお、（１８）式の各指数関数は、直交しているので、各指数関数の係数が左辺と右辺とで一致していれば、左辺全体と右辺全体とが一致することになる。すなわち、以下の関係が成り立つ。
【００７３】
【数１９】
2018092218.tif 000021
【００７４】
ここでＵj,q（２√（−１）ｑπ／（２ｎ＋１））lは定数項であるから右辺に置き、左辺を未知の微分係数Ｐ0乃至Ｐm及びＱ1乃至Ｑlにまとめると、以下のように変形される。
【００７５】
【数２０】
2018092218.tif 000022
【００７６】
ステップＳ９では、この（２０）式における未知の微分係数に乗じられる係数群がｍ次元ベクトルＣj,qとして算出され、（２０）式における定数項の値がＤj,qとして以下のように算出される。
【００７７】
【数２１】
2018092218.tif 000023
【００７８】
そうすると、線型回帰処理部３３３２は、このような点群（Ｃj,q，Ｄj,q）に対して最小二乗法による線型回帰を行って微分係数Ｐ0乃至ＰmとＱ1乃至Ｑlとを算出し、第２データ格納部３３４に格納する。ここでは、第１データ格納部３３２に格納されたデータを用いて処理が行われる。
【００７９】
この算出においては、β＝Ｐ（α）＝α0Ｐ0＋α1Ｐ1＋α2Ｐ2＋・・・・＋αm-1Ｐm-1＋αmＰm＋αm+1Ｑ1＋αm+2Ｑ2＋・・・＋αm+lＱlという回帰直線の未知の係数Ｐ0乃至ＰmとＱ1乃至Ｑlとが算出される。すなわち、模式的に示せば図１１に示すように、αβ平面において、ステップＳ９で示される点群が配置された場合に、誤差が最小となるようなβ＝Ｐ（α）に含まれる係数Ｐ0乃至ＰmとＱ1乃至Ｑlとを算出する。最小二乗法については周知であるから、これ以上述べない。
【００８０】
以上のような処理を実行することによって、最も誤差が小さくＰ（ｄ／ｄｔ）ｙj（ｔ）≒Ｑ（ｄ／ｄｔ）ｕj（ｔ）となるような微分演算子の微分係数が算出される。
【００８１】
上でも述べたように、フィッティング曲線が指数関数の線型和で表されるので、高次の微分も簡単に行うことができ、上で述べた処理を高速に実行することができる。
【００８２】
図５の説明に戻り、比較部３３６は、第２データ格納部３３４に格納されたデータに基づき、今回算出され且つ前回も算出された微分係数が、前回と今回とで同じ値であるか判定する（ステップＳ１１）。例えば、前回算出された微分係数が図９に示すようにＰ0であり、今回算出された微分係数が図１０に示すようにＰ0及びＰ1である場合、Ｐ0が前回と今回とで同じ値であるか判定される。
【００８３】
今回算出され且つ前回も算出された微分係数が、前回と今回とで同じ値ではない場合（ステップＳ１１：Ｎｏルート）、比較部３３６は、更新部３３５に対して通知を行う。これに応じ、更新部３３５は、ｍの値及びｌの値を更新する（ステップＳ１５）。処理はステップＳ３に戻る。なお、ステップＳ１１の処理を初めて実行する場合も処理はステップＳ１１のＮｏルートに進む。
【００８４】
本実施の形態においては、更新部３３５は、例えば図１２に示した情報に従ってｍの値及びｌの値を更新する。図１２の例では、まずｍ＋ｌ＝０が成立するように更新が行われ、次にｍ＋ｌ＝１が成立するように更新が行われ、・・・というように、ｍ＋ｌ＝Ｘを満たすｍとｌとの組合せがなくなった場合に、ｍ＋ｌ＝Ｘ＋１を満たすｍとｌとの組合せが採用されることが繰り返し行われる。このように、ｍ＋ｌの値がなるべく変わらないように更新を行うことで、ステップＳ７の処理を行う回数を減らせるようになる。
【００８５】
一方、今回算出され且つ前回も算出された微分係数が、前回と今回とで同じ値である場合（ステップＳ１１：Ｙｅｓルート）、比較部３３６は、以下の処理を実行する。具体的には、比較部３３６は、第２データ格納部３３４に格納されたデータに基づき、今回算出され且つ前回算出されていない微分係数の値が０であるか判定する（ステップＳ１３）。例えば、前回算出された微分係数が図９に示すようにＰ0であり、今回算出された微分係数が図１０に示すようにＰ0及びＰ1である場合、Ｐ1が０であるか判定される。
【００８６】
今回算出され且つ前回算出されていない微分係数の値が０ではない場合（ステップＳ１３：Ｎｏルート）、処理はステップＳ１５に移行する。なお、ステップＳ１３の処理を初めて実行する場合も処理はステップＳ１３のＮｏルートに進む。
【００８７】
一方、今回算出され且つ前回算出されていない微分係数の値が０である場合（ステップＳ１３：Ｙｅｓルート）、比較部３３６は、今回又は前回算出された微分係数の値を含む通知を学習結果として第２演算部３４３に送信する。これに応じ、第２演算部３４３は、通知された微分係数の値を用いた予測処理の実行を開始する。
【００８８】
ここまでの処理について、具体例を用いて説明を追加する。例えばｍ＝０且つｌ＝０の場合にはＹj,0及びＵj,0が算出され、Ｙj,0及びＵj,0に基づきＰ0が算出される。そして初回の処理なのでｍの値およびｌの値が更新される。ここでは、ｍ＝１且つｌ＝０に更新されたとする。この場合にはＹj,1、Ｕj,1、Ｙj,-1及びＵj,-1が算出され、算出されたＹj,1、Ｕj,1、Ｙj,-1及びＵj,-1とＹj,0及びＵj,0とに基づき、Ｐ0及びＰ1が算出される。
【００８９】
ここで、ｍ＝０且つｌ＝０のときに算出されたＰ0と、ｍ＝１且つｌ＝０のときに算出されたＰ0とが同じであるか判定され、同じである場合にはｍ＝１且つｌ＝０のときに算出されたＰ1＝０であるか判定される。Ｐ0の値が同じであると判定され且つＰ1＝０である場合には、予測処理部３４０の処理が開始される。一方、Ｐ0の値が同じではない又はＰ1＝０ではない場合、ｍの値及びｌの値が更新される。ここでは、ｍ＝０且つｌ＝１に更新される。そして、ｍ＝０且つｌ＝１のときに算出された微分係数について、上記と同様の処理が行われる。
【００９０】
別の例として、以下の（２２）式に示す関数ｕ（ｔ）を（２ｎ＋１）点でサンプルした離散データ（ｕ（０），ｕ（１），．．．，ｕ（２ｎ））を入力とし、以下の（２３）式に示す関数ｙ（ｔ）を（２ｎ＋１）点でサンプルした離散データ（ｙ（０），ｙ（１），．．．，ｙ（２ｎ））を出力とするシステムの同定を行うことを考える。
【００９１】
【数２２】
2018092218.tif 000024
【００９２】
【数２３】
2018092218.tif 000025
【００９３】
この場合、ｙ（ｔ）＝（（ｄ／ｄｔ）＋１）ｕ（ｔ）が成立する。
【００９４】
図１３に示すように、ｍ＝０且つｌ＝１のときにＰ0＝１及びＱ1＝１が算出され、ｍの値及びｌの値が更新され、ｍ＝０且つｌ＝２のときにＰ0＝１、Ｑ1＝１およびＱ2＝０が算出される。すると、正しく同定が行われたと判定され、予測処理部３４０の処理が開始される。この判定には、或るｍの値およびｌの値で同定が正しく同定が行われた場合、更新後のｍの値及びｌの値によって新たに算出された微分係数の値は０になり且つ元々算出された微分係数の値については更新前後で差が無いという性質が利用されている。本実施の形態においては、（２）式に示したようにＱ（ｄ／ｄｔ）を定義したことによってこの性質を利用することが可能になっている。
【００９５】
背景技術の欄で述べたように、２０１５−１３５６６２号公報の技術を利用する場合には、ｍとｌとの組合せの各々について、予測処理部３４０は同定が適切か否かの検証を行うことになる。一方、本実施の形態の方法によれば、予測処理部３４０は同定が適切か否かの検証を毎回行わなくてもよいので、同定にかかる時間を短縮できるようになる。
【００９６】
また、高周波成分に対応する高階微分の項が不要であり且つ低階微分の項のみを含む微分方程式で表されるシステムに対しても、同定を適切に行うことができるようになる。
【００９７】
図５の説明に戻り、予測処理部３４０は、学習処理部３３０の学習結果に基づき、予測処理を実行する（ステップＳ１７）。予測処理については、図１４を用いて詳細に説明する。
【００９８】
なお、学習処理を行った後は、車両１０００の検出装置１００は、加速度センサ１０１の測定データのみ、すなわちディジタル入力のデータを通信部１０３から情報処理装置３００に送信することになる。これによって情報処理装置３００への送信データ量が削減される。ここでは、学習処理で用いたディジタル入力ではなく、別のディジタル入力として新たに加速度センサ１０１の測定データを取得する。
【００９９】
従って、情報処理装置３００の通信部３１０は、ディジタル入力のデータを、車両１０００から受信すると、入出力データ格納部３２０に格納する。
【０１００】
そして、予測処理部３４０のフィッティング曲線生成部３４１は、入出力データ格納部３２０からディジタル入力のデータＵ（＝（ｕ0，ｕ1，．．，ｕ2n））を読み出す（図１４：ステップＳ２１）。そして、フィッティング曲線生成部３４１は、ＤＦＴにより、ディジタル入力Ｕについてのフィッティング曲線ｕ（ｔ）を生成し、当該曲線のデータを、第３データ格納部３４２に格納する（ステップＳ２３）。（１３）式におけるＵq（ｑは−ｎ乃至＋ｎ）を算出する。なお、ｊ＝１のみなので、サフィックスｊは省略されている。
【０１０１】
そして、第２演算部３４３の第２係数算出部３４３１は、フィッティング処理の結果Ｕq及び学習処理部３３０から通知された微分演算子の微分係数Ｐ0乃至ＰmとＱ1乃至Ｑlとを用いて、出力関数ｙ（ｔ）の係数Ｙq（ｑは−ｎ乃至＋ｎ）を算出する（ステップＳ２５）。
【０１０２】
具体的には、以下の式にてＹqを算出する。
【０１０３】
【数２４】
2018092218.tif 000026
【０１０４】
そうすると、係数Ｙqから出力関数ｙ（ｔ）は、以下のように表される。
【０１０５】
【数２５】
2018092218.tif 000027
【０１０６】
よって、出力データ算出部３４３２は、（２５）式に従って、予め定められたｔ（＝０，１，．．．，２ｎ）について、出力データＹ＝（ｙ0，ｙ1，．．．，ｙ2n）を算出し、第４データ格納部３４４に格納する（ステップＳ２７）。
【０１０７】
このようにすれば学習処理によって得られた微分演算子の係数を用いて、走行距離を推定し、出力することが可能となる。
【０１０８】
本実施の形態の方法によれば、例えば人の動作の力学的な物理モデルのような、低階微分を含む微分方程式と離散データの低周波部分とが重要なモデルに対して、同定にかかる時間を短縮することができ、また、適切な微分方程式を特定できる。例えば、２階微分を含む（２（ｄ／ｄｔ）2＋（ｄ／ｄｔ）＋１）ｙ（ｔ）＝ｕ（ｔ）という微分方程式のシステムについて、サンプル数２ｎ＋１＝２００１の入出力離散データを使用して同定を行った場合、既存の方法と比べて計算量を１／２０にすることができるようになる。
【０１０９】
以上本発明の実施の形態を説明したが、本発明はこれに限定されるものではない。例えば、車両１０００に設けた加速度センサ１０１及び距離センサ１０２を前提にした説明を行ったが、同定を行うべきシステムの入力及び出力に応じた他のセンサによって測定されるデータであっても良い。
【０１１０】
また、図３及び図４に示した機能ブロック図は一例であって、プログラムモジュール構成とは一致しない場合もある。さらに処理フローについても、処理結果が変わらない限り、処理順番を入れ替えたり、複数ステップを並列実行するようにしても良い。
【０１１１】
さらに、図３及び図４では、検出装置１００と情報処理装置３００とが分けて実装されていたが、用途によっては一体化された装置を実装する場合もある。さらに、検出装置１００は、例えば携帯電話機やカーナビゲーション装置などの装置に組み込まれる場合もある。
【０１１２】
さらに、システム同定だけではなく、人工知能の学習問題にも応用が可能である。
【０１１３】
なお、上で述べた情報処理装置３００は、コンピュータ装置であって、図１５に示すように、メモリ２５０１とＣＰＵ（Central Processing Unit）２５０３とハードディスク・ドライブ（ＨＤＤ：Hard Disk Drive）２５０５と表示装置２５０９に接続される表示制御部２５０７とリムーバブル・ディスク２５１１用のドライブ装置２５１３と入力装置２５１５とネットワークに接続するための通信制御部２５１７とがバス２５１９で接続されている。オペレーティング・システム（ＯＳ：Operating System）及び本実施例における処理を実施するためのアプリケーション・プログラムは、ＨＤＤ２５０５に格納されており、ＣＰＵ２５０３により実行される際にはＨＤＤ２５０５からメモリ２５０１に読み出される。ＣＰＵ２５０３は、アプリケーション・プログラムの処理内容に応じて表示制御部２５０７、通信制御部２５１７、ドライブ装置２５１３を制御して、所定の動作を行わせる。また、処理途中のデータについては、主としてメモリ２５０１に格納されるが、ＨＤＤ２５０５に格納されるようにしてもよい。本技術の実施例では、上で述べた処理を実施するためのアプリケーション・プログラムはコンピュータ読み取り可能なリムーバブル・ディスク２５１１に格納されて頒布され、ドライブ装置２５１３からＨＤＤ２５０５にインストールされる。インターネットなどのネットワーク及び通信制御部２５１７を経由して、ＨＤＤ２５０５にインストールされる場合もある。このようなコンピュータ装置は、上で述べたＣＰＵ２５０３、メモリ２５０１などのハードウエアとＯＳ及びアプリケーション・プログラムなどのプログラムとが有機的に協働することにより、上で述べたような各種機能を実現する。
【０１１４】
以上述べた本実施の形態をまとめると、以下のようになる。
【０１１５】
本実施の形態の第１の態様に係る係数決定方法は、システムを表す微分方程式の次数を第１の次数に設定した場合の第１の微分係数を、システムへのディジタル入力及びシステムからのディジタル出力に対する離散フーリエ変換の結果に基づき算出し、微分方程式の次数を第２の次数に設定した場合の第２の微分係数を記憶装置から読み出し、第１の微分係数と読み出した第２の微分係数との比較に基づき、微分方程式の微分係数の収束について判定し、判定の結果に基づき、微分方程式の微分係数を決定する処理を含む。
【０１１６】
このようにすれば、算出した微分係数が適切であるか否かの検証を毎回行わなくてもよいので、微分方程式で表されるシステムの同定に要する計算量を減らすことができるようになる。
【０１１７】
また、第１の次数が第２の次数より大きくてもよい。
【０１１８】
例えば人の運動や車両の移動など、力学的な要因で発生する事象は、低次の微分方程式で表されることが多い。上で述べたようにすれば、低次の微分方程式から先に処理が行われるようになるので、特に上記のような事象を対象とする場合に同定に要する時間を短くすることができる。すなわち、同定に要する時間を減らすことができるようになる。
【０１１９】
また、本係数決定方法は、（Ｄ）微分方程式に複数の微分演算子が含まれる場合、複数の微分演算子の次数の組合せを、複数の微分演算子の次数の総和が第２の次数以上になり、且つ、微分方程式の次数を第２の次数に設定した場合における複数の微分演算子の次数の組合せと異なるように設定する処理をさらに含んでもよい。
【０１２０】
離散フーリエ変換が無駄に行われて計算量が増えることを抑制できるようになる。
【０１２１】
また、微分方程式の微分係数の収束について判定を行う処理において、（ｂ１）第１の微分係数及び第２の微分係数のいずれにも含まれる微分係数が、微分方程式の次数を第２の実数に設定した場合と微分方程式の次数を第２の次数に設定した場合とで同じ値であり、且つ、第１の微分係数に含まれ且つ第２の微分係数に含まれない微分係数の値がゼロであるか判定してもよい。
【０１２２】
適切な微分係数が算出されたことを確認できるようになる。
【０１２３】
また、微分方程式の微分係数を決定する処理において、（ｃ１）判定の結果が微分方程式の微分係数が収束したことを示している場合、微分方程式の微分係数を決定し、判定の結果が部分方程式の微分係数が収束していないことを示している場合、算出した微分係数を記憶装置に格納し且つ微分方程式の次数を第３の次数に設定してもよい。
【０１２４】
算出結果を使用して次の処理を実行できるので、計算量の増加を抑制できるようになる。
【０１２５】
本実施の形態の第２の態様に係る係数決定装置（例えば、情報処理装置３００）は、（Ｅ）システムを表す微分方程式の次数を第１の次数に設定した場合の微分係数を、システムへのディジタル入力及びシステムからのディジタル出力に対する離散フーリエ変換の結果に基づき算出する算出部（例えば、第１演算部３３３）と、（Ｆ）微分方程式の次数を第２の次数に設定した場合の微分係数を記憶装置から読み出し、算出した微分係数と読み出した微分係数との比較に基づき、微分方程式の微分係数の収束について判定を行い、判定の結果に基づき、微分方程式の微分係数を決定する決定部（例えば、比較部３３６）とを有する。
【０１２６】
なお、上で述べたような処理をコンピュータに実行させるためのプログラムを作成することができ、当該プログラムは、例えばフレキシブル・ディスク、ＣＤ−ＲＯＭなどの光ディスク、光磁気ディスク、半導体メモリ（例えばＲＯＭ）、ハードディスク等のコンピュータ読み取り可能な記憶媒体又は記憶装置に格納される。なお、処理途中のデータについては、ＲＡＭ等の記憶装置に一時保管される。
【０１２７】
以上の実施例を含む実施形態に関し、さらに以下の付記を開示する。
【０１２８】
（付記１）
コンピュータに、
システムを表す微分方程式の次数を第１の次数に設定した場合の第１の微分係数を、前記システムへのディジタル入力及び前記システムからのディジタル出力に対する離散フーリエ変換の結果に基づき算出し、
前記微分方程式の次数を第２の次数に設定した場合の第２の微分係数を記憶装置から読み出し、前記第１の微分係数と読み出した前記第２の微分係数との比較に基づき、前記微分方程式の微分係数の収束について判定し、
判定の結果に基づき、前記微分方程式の微分係数を決定する、
処理を実行させる係数決定プログラム。
【０１２９】
（付記２）
前記第１の次数が前記第２の次数より大きい、
付記１記載の係数決定プログラム。
【０１３０】
（付記３）
前記コンピュータに、
前記微分方程式に複数の微分演算子が含まれる場合、前記複数の微分演算子の次数の組合せを、前記複数の微分演算子の次数の総和が前記第２の次数以上になり、且つ、前記微分方程式の次数を前記第２の次数に設定した場合における前記複数の微分演算子の次数の組合せと異なるように設定する、
処理をさらに実行させる付記１記載の係数決定プログラム。
【０１３１】
（付記４）
前記微分方程式の微分係数の収束について判定を行う処理において、
前記第１の微分係数及び前記第２の微分係数のいずれにも含まれる微分係数が、前記微分方程式の次数を前記第２の実数に設定した場合と前記微分方程式の次数を前記第２の次数に設定した場合とで同じ値であり、且つ、前記第１の微分係数に含まれ且つ前記第２の微分係数に含まれない微分係数の値がゼロであるか判定する、
付記１記載の係数決定プログラム。
【０１３２】
（付記５）
前記微分方程式の微分係数を決定する処理において、
前記判定の結果が前記微分方程式の微分係数が収束したことを示している場合、前記微分方程式の微分係数を決定し、前記判定の結果が前記部分方程式の微分係数が収束していないことを示している場合、算出した前記微分係数を前記記憶装置に格納し且つ前記微分方程式の次数を第３の次数に設定する、
付記１記載の係数決定プログラム。
【０１３３】
（付記６）
コンピュータが、
システムを表す微分方程式の次数を第１の次数に設定した場合の第１の微分係数を、前記システムへのディジタル入力及び前記システムからのディジタル出力に対する離散フーリエ変換の結果に基づき算出し、
前記微分方程式の次数を第２の次数に設定した場合の第２の微分係数を記憶装置から読み出し、前記第１の微分係数と読み出した前記第２の微分係数との比較に基づき、前記微分方程式の微分係数の収束について判定し、
判定の結果に基づき、前記微分方程式の微分係数を決定する、
処理を実行する係数決定方法。
【０１３４】
（付記７）
システムを表す微分方程式の次数を第１の次数に設定した場合の第１の微分係数を、前記システムへのディジタル入力及び前記システムからのディジタル出力に対する離散フーリエ変換の結果に基づき算出する算出部と、
前記微分方程式の次数を第２の次数に設定した場合の第２の微分係数を記憶装置から読み出し、前記第１の微分係数と読み出した前記第２の微分係数との比較に基づき、前記微分方程式の微分係数の収束について判定し、判定の結果に基づき、前記微分方程式の微分係数を決定する決定部と、
を有する係数決定装置。
【符号の説明】
【０１３５】
１００  検出装置
１０１  加速度センサ
１０２  距離センサ
１０３  通信部
２００  ネットワーク
３００  情報処理装置
３１０  通信部
３２０  入出力データ格納部
３３０  学習処理部
３４０  予測処理部
３３１  ＤＦＴ計算部
３３２  第１データ格納部
３３３  第１演算部
３３４  第２データ格納部
３３５  更新部
３３６  比較部
３４１  フィッティング曲線生成部
３４２  第３データ格納部
３４３  第２演算部
３４４  第４データ格納部
【図１】
2018092218.tif 000028
【図２】
2018092218.tif 000029
【図３】
2018092218.tif 000030
【図４】
2018092218.tif 000031
【図５】
2018092218.tif 000032
【図６】
2018092218.tif 000033
【図７】
2018092218.tif 000034
【図８】
2018092218.tif 000035
【図９】
2018092218.tif 000036
【図１０】
2018092218.tif 000037
【図１１】
2018092218.tif 000038
【図１２】
2018092218.tif 000039
【図１３】
2018092218.tif 000040
【図１４】
2018092218.tif 000041
【図１５】
2018092218.tif 000042