(57)【特許請求の範囲】
【請求項１】
動的システムの周波数特性を示す周波数応答関数からｓ演算子もしくはｚ演算子の有理関数形式の伝達関数を同定するモデル同定装置において、
前記伝達関数の周波数特性と前記周波数応答関数が示す周波数特性との各周波数点での差の二乗を最小にするために、前記伝達関数を構成する分母多項式及び分子多項式の未知係数を未知変数とし、前記周波数応答関数が示す周波数特性を定数とした第１の連立方程式を作成し、この第１の連立方程式に対して総合最小二乗法を用いて得た第１の解ｘを導出し、この第１の解ｘを前記第１の連立方程式に代入して得られる残余を定数項として作成した第２の連立方程式に対して総合最小二乗法を用いて得た第２の解ｅを導出し、これら２つの解ｘ，ｅの和ｙを伝達関数の係数とする伝達関数同定手段
を備えたことを特徴とするモデル同定装置。
【請求項２】
前記伝達関数同定手段は、前記第１の連立方程式に対して重み付き最小二乗法の反復計算を適用し、
その重みを、前記第１の連立方程式に対し総合最小二乗法を用いて得た第１の解ｘと、この第１の解ｘを前記第２の連立方程式に代入して得られる残余を定数項として作成した第２の連立方程式に対して総合最小二乗法を用いて得た解ｅとの和ｙを基に設定し、
前記反復計算を、前記二つの解ｘ，ｅの和ｙを代入して得られる伝達関数の周波数応答特性と前記周波数応答関数の周波数特性との相対偏差が設定した閾値を下回っていた場合、もしくは反復計算の前回得られた解と今回得られた解の相対誤差が設定した閾値を下回っていた場合に終了する
ことを特徴とする請求項１に記載のモデル同定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、動的システムの周波数特性を示す周波数応答関数から数学的モデルを同定するモデル同定装置に関する。
【背景技術】
【０００２】
測定により事前情報として得られた周波数応答関数Ｈ（ｚ）からその特性を示す伝達関数Ｇ（ｚ）を導出する問題は、数式１に示すように、上記の伝達関数を構成する分母多項式及び分子多項式の未知係数を未知変数とし、上記の周波数応答関数が示す周波数特性を定数とした連立方程式を定義し、線形連立方程式を解く問題に帰着することができる。
【数１】
0004470855.tif 000002
【０００３】
そして、未知変数よりも式の数が多い場合、その解は一意に求まらず、上記の伝達関数の周波数特性と上記の周波数応答関数が示す周波数特性との各周波数点での差の二乗が最小となる線形連立方程式の解を求める最小二乗問題（数式２）となる。
【数２】
0004470855.tif 000003
【０００４】
また、数式２で問題となる線形化による解の偏り（線形化誤差）を解消するために、例えば数式３に示すように解を用いた重みを付加し、重みを更新しながら反復計算を行って収束解を求める。この収束解を伝達関数の係数とする。
【数３】
0004470855.tif 000004
【０００５】
以上説明した特許文献１における従来の最小二乗問題では、数式１で示される連立方程式を構成する係数項行列Ａにはエラーがなく、測定により得られた各周波数点での周波数特性データで構成される列ベクトルｄにのみエラーが存在することを前提として最小二乗解を求めるようになっている。
しかしながら、数式１、数式２に示したように、行列Ａには、線形化を行ったことにより、測定で得られた周波数特性が含まれており、行列Ａにエラーがないと仮定して得られる最小二乗解には測定誤差に起因する偏り誤差が含まれる欠点がある。
【０００６】
また、数式３に基づいた反復計算においては、偏り誤差を含んだ解を基に重みを設定することから、測定エラー以外のエラーが行列Ａに含まれることになり、反復計算が収束しないか、もしくは収束解が得られても線形化誤差を解消することができないという欠点がある。
これらの欠点を補うために、近年、列ベクトルｄだけではなく係数項行列Ａにもエラーが含まれることを考慮して連立方程式の最小二乗解を得ることができる手法として総合最小二乗法（Ｔｏｔａｌ  Ｌｅａｓｔ  Ｓｑｕａｒｅ）という手法が提案されている（非特許文献２参照）。
【０００７】
【非特許文献１】野田  琢、長岡  直人、「線形化最小二乗方によるＡＲＭＡ過渡現象計算モデルの開発」、電気学会論文誌Ｂ部門、１１４巻４号、１９９４年、Ｐ−３９６−４０１
【非特許文献２】Ｓａｂｉｎｅ  Ｖａｎ  Ｈｕｆｆｅｌ，  Ｊｏｏｓ  Ｖａｎｄｅｗａｌｌｅ，  “Ｔｈｅ  Ｔｏｔａｌ  Ｌｅａｓｔ  Ｓｑｕａｒｅｓ  Ｐｒｏｂｌｅｍ  Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ａｓｐｅｃｔｓ  ａｎｄ  Ａｎａｌｙｓｉｓ”，  Ｓｏｃｉｅｔｙ  ｆｏｒ  Ｉｎｄｕｓｔｒｉａｌ  ａｎｄ  Ａｐｐｌｉｅｄ  Ｍａｔｈｅｍａｔｉｃｓ，  １９９１
【発明の開示】
【発明が解決しようとする課題】
【０００８】
上記特許文献における総合最小二乗法を数式２、数式３に適用することは、係数項行列Ａと列ベクトルｂで構成される行列に対して特異値分解を行うことで容易に実行できる。しかし、総合最小二乗法はデータに含まれる過大なエラーに対して従来の最小二乗法に比べて影響を受けやすくなるため、解の精度の劣化や、反復計算を行う際の重みの変化に過敏になり収束解が得られないといった問題がある。
【０００９】
本発明は、このような課題に鑑みてなされたものであり、エラーを含む周波数応答関数を基に伝達関数を同定する場合でも、エラーを考慮した精度の良い最小二乗解を得ることができ、また反復計算の過程での線形化誤差に影響されず、最小二乗解の精度の劣化を回避することができ、更に安定して収束解を得ることができるモデル同定装置を提供することを目的としている。
【課題を解決するための手段】
【００１０】
上記目的を達成するために、本発明の請求項１によるモデル同定装置は、動的システムの周波数特性を示す周波数応答関数からｓ演算子もしくはｚ演算子の有理関数形式の伝達関数を同定するモデル同定装置において、前記伝達関数の周波数特性と前記周波数応答関数が示す周波数特性との各周波数点での差の二乗を最小にするために、前記伝達関数を構成する分母多項式及び分子多項式の未知係数を未知変数とし、前記周波数応答関数が示す周波数特性を定数とした第１の連立方程式を作成し、この第１の連立方程式に対して総合最小二乗法を用いて得た第１の解ｘを導出し、この第１の解ｘを前記第１の連立方程式に代入して得られる残余を定数項として作成した第２の連立方程式に対して総合最小二乗法を用いて得た第２の解ｅを導出し、これら２つの解ｘ，ｅの和ｙを伝達関数の係数とする伝達関数同定手段を備えたことを特徴とする。
この構成によれば、２つの解ｘ，ｅの和ｙを伝達関数の係数とする同定処理、即ち、２つの解ｘ，ｅの和を真の総合最小二乗解ｙとする同定処理を行うことで、エラーを含む周波数応答関数Ｈ（ｚ）を基に伝達関数Ｇ（ｚ）を同定する際においても、エラーを考慮した精度のよい最小二乗解を得ることができる。
【００１１】
また、本発明の請求項２によるモデル同定装置は、前記伝達関数同定手段は、前記第１の連立方程式に対して重み付き最小二乗法の反復計算を適用し、その重みを、前記第１の連立方程式に対し総合最小二乗法を用いて得た第１の解ｘと、この第１の解ｘを前記第２の連立方程式に代入して得られる残余を定数項として作成した第２の連立方程式に対して総合最小二乗法を用いて得た解ｅとの和ｙを基に設定し、前記反復計算を、前記二つの解ｘ，ｅの和ｙを代入して得られる伝達関数の周波数応答特性と前記周波数応答関数の周波数特性との相対偏差が設定した閾値を下回っていた場合、もしくは反復計算の前回得られた解と今回得られた解の相対誤差が設定した閾値を下回っていた場合に終了することを特徴とする。
この構成によれば、重みｗを更新して再度、伝達関数の係数（真の総合最小二乗解）ｙを求める反復計算の過程での線形化誤差に影響されず、解の精度の劣化を回避することができ、更に、真の総合最小二乗解ｙのみで重みｗの更新を行うことで、安定に収束解を得ることができる。
【発明の効果】
【００１２】
以上説明したように本発明によれば、エラーを含む周波数応答関数を基に伝達関数を同定する場合でも、エラーを考慮した精度の良い最小二乗解を得ることができ、また反復計算の過程での線形化誤差に影響されず、最小二乗解の精度の劣化を回避することができ、更に安定して収束解を得ることができるという効果がある。
【発明を実施するための最良の形態】
【００１３】
以下、本発明の実施の形態を、図面を参照して説明する。但し、本明細書中の全図において相互に対応する部分には同一符号を付し、重複部分においては後述での説明を適時省略する。
（第１の実施の形態）
図１は、本発明の第１の実施の形態に係るモデル同定装置の構成を示す図である。
図１に示すモデル同定装置１０は、関数読込部１２及び伝達関数同定部１４を備えて構成され、伝達関数同定部１４は、第１の連立方程式作成部１６と、第１の行列作成部１８と、第１の特異値分解部２０と、第１の総合最小二乗解計算部２２と、残余計算部２４と、第２の連立方程式作成部２６と、第２の行列作成部２８と、第２の特異値分解部３０と、第２の総合最小二乗解計算部３２と、総合最小二乗解加算部３４とを備えて構成されている。
【００１４】
関数読込部１２は、対象となる動的システムの周波数特性を示す周波数応答関数Ｈ（ｚ）を読み込むものである。
第１の連立方程式作成部１６は、その読み込まれた周波数応答関数Ｈ（ｚ）を基に、次数をｐと仮定した伝達関数Ｇ（ｚ）を定義し、この伝達関数Ｇ（ｚ）を構成する分母多項式及び分子多項式の未知係数を未知変数とし、周波数応答関数Ｈ（ｚ）が示す周波数特性を定数とした連立方程式（数式１）を作成するものである。
第１の行列作成部１８は、数式１に示す係数項行列Ａと定数項列ベクトルｄを用いてｍ行２ｐ＋２列の行列Ｂ（数式４）を作成するものである。
【数４】
0004470855.tif 000005
【００１５】
第１の特異値分解部２０は、数式４に示す行列Ｂの特異値分解を行うものである（数式５）。
【数５】
0004470855.tif 000006
【００１６】
ここで、行列Ｕは左特異列ベクトルで構成されるｍ行ｍ列の直交行列、行列Ｖは右特異ベクトルで構成される２ｐ＋２行２ｐ＋２列の直交行列、行列Ｓは対角要素が特異値で構成されるｍ行２ｐ＋２列の対角行列である。
第１の総合最小二乗解計算部２２は、行列Ｖの２ｐ＋２列目の列ベクトルを用いて、数式１の第１の総合最小二乗解ｘを計算するものである（数式６）。
【数６】
0004470855.tif 000007
【００１７】
ここで、ｖi,jは行列Ｖのｉ行ｊ列目の要素を示す。
残余計算部２４は、第１の総合最小二乗解計算部２２で得られた解ｘを数式１に代入し、数式１の残余ｒを計算するものである（数式７）。
【数７】
0004470855.tif 000008
【００１８】
第２の連立方程式作成部２６は、数式７で求まった残余ｒと数式１の係数項行列Ａを用いて、変数ｅとした新しい連立方程式を作成するものである（数式８）。
【数８】
0004470855.tif 000009
【００１９】
第２の行列作成部２８は、第１の行列作成部１８のように、係数項行列Ａと定数項列ベクトルｒを用いてｍ行２ｐ＋２列の行列Ｃ（数式９）を作成するものである。
【数９】
0004470855.tif 000010
【００２０】
第２の特異値分解部３０は、第１の特異値分解部２０のように、数式９に示す行列Ｃの特異値分解を行うものである（数式１０）。
【数１０】
0004470855.tif 000011
【００２１】
ここで、行列Ｕｃは左特異列ベクトルで構成されるｍ行ｍ列の直交行列、行列Ｖｃは右特異ベクトルで構成される２ｐ＋２行２ｐ＋２列の直交行列、行列Ｓｃは対角要素が特異値で構成されるｍ行２ｐ＋２列の対角行列である。
第２の総合最小二乗解計算部３２は、行列Ｖｃの２ｐ＋２列目の列ベクトルを用いて、数式８の第２の総合最小二乗解ｅを計算するものである（数式１１）。
【数１１】
0004470855.tif 000012
【００２２】
ここで、ｖci,jは行列Ｖｃのｉ行ｊ列目の要素を示す。
総合最小二乗解加算部３４は、第１と第２の２つの総合最小二乗解ｘ，ｅの和ｙを計算し、このｙを数式１の真の総合最小二乗解とするものである。
このような構成のモデル同定装置１０による伝達関数Ｇ（ｚ）の同定処理の動作を、図２に示すフローチャートを参照して説明する。
【００２３】
ステップＳ１において、関数読込部１２で対象となる動的システムの周波数特性を示す周波数応答関数Ｈ（ｚ）が読み込まれる。ステップＳ２において、上記ステップＳ１で読み込まれた周波数応答関数Ｈ（ｚ）を基に、次数をｐと仮定した伝達関数Ｇ（ｚ）が定義され、この伝達関数を構成する分母多項式及び分子多項式の未知係数が未知変数とされ、先の周波数応答関数Ｈ（ｚ）が示す周波数特性を定数とした数式１の連立方程式Ａｘ＝ｄが作成される。
【００２４】
続いて、ステップＳ３において、第１の行列作成部１８で数式１に示す係数項行列Ａと定数項列ベクトルｄが用いられて、数式４のｍ行２ｐ＋２列の行列Ｂ＝[Ａ  ｄ]が作成される。ステップＳ４において、数式４に示す行列Ｂの特異値分解が数式５のＢ＝ＵＳＶTように行われる。
そして、ステップＳ５において、第１の総合最小二乗解計算部２２にて、数式５における前述した行列Ｖの２ｐ＋２列目の列ベクトルが用いられて、数式１の第１の総合最小二乗解ｘが数式６のように計算される。
【００２５】
次に、ステップＳ６において、上記ステップＳ５で得られた解ｘが数式１に代入され、数式１の残余ｒが数式７のｒ＝ｄ−Ａｘのように計算される。
ステップＳ７において、数式７で求まった残余ｒと数式１の係数項行列Ａとが用いられて、変数ｅとされた新しい連立方程式が数式８のＡｅ＝ｒのように作成される。
次に、ステップＳ８において、上記ステップＳ３のように、第２の行列作成部２８にて係数項行列Ａと定数項列ベクトルｒが用いられて、数式９のＣ＝[Ａ  ｒ]の行列が作成される。
【００２６】
ステップＳ９において、上記ステップＳ４のように、数式９に示す行列Ｃの特異値分解が数式１０のＣ＝ＵｃＳｃＶｃTように行われる。
ステップＳ１０において、数式９における前述した行列Ｖｃの２ｐ＋２列目の列ベクトルが用いられて、数式８の第２の総合最小二乗解ｅが数式１１のように計算される。
そして、ステップＳ１１において、総合最小二乗解加算部３４にて上記ステップＳ５とＳ１０で得られた２つの解ｘ、ｅの和である真の総合最小二乗解ｙが計算されて処理が終了する。
【００２７】
このように第１の実施の形態のモデル同定装置１０によれば、周波数応答関数Ｈ（ｚ）が示す周波数特性を定数とし、同定する伝達関数Ｇ（ｚ）を構成する分母多項式及び分子多項式の未知係数を未知変数とした連立方程式に対して第１の総合最小二乗解ｘを求め、更に、その解ｘを基に残余を計算し、その残余から残余に対する第２の総合最小二乗解ｅを計算し、２つの解ｘ，ｅの和を真の総合最小二乗解ｙとする処理を行うことで、エラーを含む周波数応答関数Ｈ（ｚ）を基に伝達関数Ｇ（ｚ）を同定する際においても、エラーを考慮した精度のよい最小二乗解を得ることができる。
【００２８】
（第２の実施の形態）
図５は、本発明の第２の実施の形態に係るモデル同定装置の構成を示す図である。
図５に示すモデル同定装置４０は、関数読込部１２と、重み初期設定部４２と、伝達関数同定部４４を備えて構成され、伝達関数同定部４４は、第１の連立方程式作成部１６と、第１の行列作成部１８と、第１の特異値分解部２０と、第１の総合最小二乗解計算部２２と、残余計算部２４と、第２の連立方程式作成部２６と、第２の行列作成部２８と、第２の特異値分解部３０と、第２の総合最小二乗解計算部３２と、総合最小二乗解加算部３４と、計算収束判定部４６と、重み更新部４８とを備えて構成されている。
【００２９】
重み初期設定部４２は、重み（重み係数）の初期設定を行うものである。但し、重みの初期設定は単位行列を使う。
計算収束判定部４６は、真の総合最小二乗解ｙを求める計算の収束判定を後述で説明するように行い、収束していないと判定された場合に、その未収束の判定結果を重み更新部４８へ出力する。
反復計算の収束判定は、同定の収束性もしくは解の収束性を判別して打ち切り判定を行う。例えば、真の総合最小二乗解ｙを代入して得られる伝達関数の周波数応答特性と周波数応答関数の周波数特性との相対偏差が設定された閾値を下回っていた場合（数式１２）、又は反復過程の前回の計算で得られた解と今回得られた解の相対誤差が設定した閾値を下回っていた場合（数式１３）に反復計算を打ち切る。
【数１２】
0004470855.tif 000013
【００３０】
【数１３】
0004470855.tif 000014
【００３１】
重み更新部４８は、計算収束判定部４６から未収束の判定結果が入力された場合に、重み（重み係数）の更新を行って第１の連立方程式作成部１６へ出力する。その重みは数式３に示したように、伝達関数同定部４４の構成要素１６〜３４，４６，４８の順に行われる反復計算の過程での前回の計算で得られた真の総合最小二乗解ｙを基に設定する。この重みを更新した後、再度、第１の連立方程式作成部１６に戻って順次同様の処理にて真の総合最小二乗解ｙの計算を行うようになっている。
【００３２】
このような構成のモデル同定装置４０による伝達関数Ｇ（ｚ）の同定処理の動作を、図４に示すフローチャートを参照して説明する。
ステップＳ２１にて関数読込部１２で周波数応答関数Ｈ（ｚ）が読み込まれ、ステップＳ２２にて、重み初期設定部４２にて重みｗ＝１とする初期設定が行われる。この後、ステップＳ２３〜Ｓ３２においては、上記第１の実施の形態のステップＳ２〜Ｓ１１と同様に処理が行われる。この処理にて真の総合最小二乗解ｙが求められると、ステップＳ３３において、計算収束判定部４６にて収束判定が行われる。この結果、収束していないと判断された場合は、ステップＳ３４において、重み更新部４８にて重みｗが数式３に示すように、反復過程の前回の計算で得られた真の総合最小二乗解ｙを基に更新される。この更新後、ステップＳ２３に戻り、その更新された重みｗが第１の連立方程式作成部１６に設定された後、上記同様に再度、真の総合最小二乗解ｙが求められる。そして、ステップＳ３３にて計算の収束判定が行われ、この結果、何れかの上記閾値を下回っていた場合、反復計算が打ち切られ、所定の伝達関数Ｇ（ｚ）が得られる。
【００３３】
このように第２の実施の形態のモデル同定装置４０によれば、上記第１の実施の形態と同様にエラーを含む周波数応答関数を基に伝達関数を同定する際においても、エラーを考慮した精度のよい最小二乗解を得ることができる。また、重みｗを更新して再度、真の総合最小二乗解ｙを求める反復計算の過程での線形化誤差に影響されず、解の精度の劣化を回避することができ、更に、真の総合最小二乗解ｙのみで重みｗの更新を行うことで、安定に収束解を得ることができる。
【図面の簡単な説明】
【００３４】
【図１】本発明の第１の実施の形態に係るモデル同定装置の構成を示す図である。
【図２】第１の実施の形態に係るモデル同定装置による伝達関数の同定処理の動作を説明するためのフローチャートである。
【図３】本発明の第２の実施の形態に係るモデル同定装置の構成を示す図である。
【図４】第２の実施の形態に係るモデル同定装置による伝達関数の同定処理の動作を説明するためのフローチャートである。
【符号の説明】
【００３５】
１０，４０  モデル同定装置
１２  関数読込部
１４，４４  伝達関数同定部
１６  第１の連立方程式作成部
１８  第１の行列作成部
２０  第１の特異値分解部
２２  第１の総合最小二乗解計算部
２４  残余計算部
２６  第２の連立方程式作成部
２８  第２の行列作成部
３０  第２の特異値分解部
３２  第２の総合最小二乗解計算部
３４  総合最小二乗解加算部
４２  重み初期設定部
４６  計算収束判定部
４８  重み更新部
【図１】
0004470855.tif 000015
【図２】
0004470855.tif 000016
【図３】
0004470855.tif 000017
【図４】
0004470855.tif 000018