【０００１】
【発明の属する技術分野】
本発明は、特性不明なシステムに対する入力信号と応答信号とから、システムの伝達関数を推定する特性不明システムの伝達関数推定装置に関する。
【０００２】
【従来の技術】
従来の特性不明システムの伝達関数推定装置としては、ステップ状の入力信号を印加したときのステップ応答波形を用いてインパルス応答列を求め、インパルス応答列に基づいて算出した離散時間状態空間モデルを連続系状態空間モデルに変換し、この連続系状態空間モデルの次数をインパルス応答列で構成されるハンケル行列のランクから推定して伝達関数を算出するというものが知られている（例えば、特許文献１参照）。
【０００３】
また、システムの同定手法としては、既存の周波数応答データに基づいて、これに最小二乗法を適用することによりラプラス演算子ｓの有理関数形式の低次線形モデルを同定し、伝達関数推定の際には、算出した最小二乗解に対して出力誤差の標準偏差（相対誤差評価）で特性近似精度評価を行い、最も評価の良い伝達関数を選択するというものが知られている（例えば、非特許文献１参照）。
【０００４】
【特許文献１】
特公１９９４−７７２１１号（第２頁、図１）
【非特許文献１】
「周波数領域最小二乗近似による系統の低次線形モデル同定法」、電気学会論文誌Ｂ、平成１３年１月、１２１巻１号、ｐ５２−５９
【０００５】
【発明が解決しようとする課題】
しかしながら、上記特許文献１に記載の従来例にあっては、入力信号にステップ入力を適用しているので、システムに直流量を入力することになり、適用範囲が限定されるという未解決の課題がある。
また、伝達関数推定の際には、インパルス応答列から連続系状態空間モデルを求め、この連続系状態空間モデルの次数をインパルス応答列で構成されるハンケル行列のランクから推定して伝達関数を算出するので、システムの入力信号と出力応答信号から直接伝達関数を推定することができないという未解決の課題がある。
【０００６】
さらに、上記非特許文献１に記載の従来例にあっては、伝達関数を推定するためには、推定に必要な周波数領域に対する周波数応答データが予め用意されていることが前提であり、特性が不明なシステムに対しては適用できないという未解決の課題がある。
また、伝達関数推定の際に、出力誤差の標準偏差（相対誤差評価）で特性近似精度評価を行い、最も評価の良い伝達関数を選択するので、想定する複数個の伝達関数に対する相互評価を必要とするという未解決の課題がある。
そこで、本発明は、上記従来例の未解決の課題に着目してなされたものであり、システムの入力信号と過渡応答信号から直接伝達関数を推定することができる特性不明システムの伝達関数推定装置を提供することを目的としている。
【０００７】
【課題を解決するための手段】
上記目的を達成するために、請求項１に係る特性不明システムの伝達関数推定装置は、複数の異なる周波数の正弦波を合成した多重正弦波で構成される入力信号とそれに対する応答信号との関係を表す有理関数と、前記有理関数の周波数応答データとに基づいて、有理関数の未知係数を変数とした行列式を生成する行列式生成手段と、該行列式生成手段で生成した行列式から伝達関数を推定する伝達関数推定手段とを備えていることを特徴としている。
【０００８】
この請求項１に係る発明では、複数の異なる周波数の正弦波を合成した多重正弦波信号を入力信号として適用するので、システムに必要以上の直流量を入力することなく、短期間の測定でシステムの特性を近似する伝達関数を推定することができる。また、入力信号の周波数応答データに基づいて生成した行列式から伝達関数を推定するので、複数個の伝達関数に対する相互評価を必要とせず、システムの入力信号と過渡応答信号から直接伝達関数を推定することができる。
【０００９】
また、請求項２に係る特性不明システムの伝達関数推定装置は、請求項１に係る発明において、前記行列式生成手段は、入力信号とそれに対する応答信号との関係を所定の次数の差分方程式で表し、その差分方程式をｚ演算子の多項式からなる有理関数で表し、その有理関数とそれに対する入力信号の周波数帯域内の周波数応答データとに基づいて有理関数の未知係数を変数とした行列式を生成するように構成されていることを特徴としている。
【００１０】
この請求項２に係る発明では、入力信号と応答信号との関係を表す有理関数とそれに対する入力信号の周波数帯域内の周波数応答データに基づいて行列式を生成し、その行列式から伝達関数を推定するので、複数個の伝達関数に対する相互評価を必要とせず、システムの入力信号と過渡応答信号から直接伝達関数を推定することができる。
【００１１】
また、請求項３に係る特性不明システムの伝達関数推定装置は、請求項１又は２に係る発明において、前記行列式生成手段は、前記周波数応答データに基づいて、有理関数の未知係数を変数として生成した行列式に重み係数を乗算して、重み付きの行列式を生成する重み付け手段を備えていることを特徴としている。
この請求項３に係る発明では、周波数応答データに重み係数を乗算して得られる行列式から伝達関数を推定するので、伝達関数の推定誤差を低減することができる。
【００１２】
また、請求項４に係る特性不明システムの伝達関数推定装置は、請求項１乃至３の何れかに係る発明において、前記伝達関数推定手段は、伝達関数の次数を決定する次数決定手段と、伝達関数の係数を決定する係数決定手段とを有し、前記次数決定手段は、前記行列式生成手段で生成した行列の線形独立な列数が、その行列の列数と等しくなるように伝達関数の次数を決定し、前記係数決定手段は、前記次数決定手段で決定した次数に応じて、前記行列式生成手段で生成した行列式から伝達関数の係数を決定するように構成されていることを特徴としている。
【００１３】
この請求項４に係る発明では、データ行列の線形独立な列数がデータ行列の列数と等しくなるように伝達関数の次数を決定し、決定した次数に応じて伝達関数の係数を決定するので、システムの特性を推定するために十分な最小次数を有する伝達関数を推定することができる。
また、請求項５に係る特性不明システムの伝達関数推定装置は、請求項１乃至４の何れかに係る発明において、前記伝達関数推定手段で推定した伝達関数の安定性及び最小位相特性の何れかから、前記伝達関数推定手段で推定した伝達関数の補正の必要性を判断する伝達関数補正判断手段と、該伝達関数補正判断手段により伝達関数が不安定及び非最小位相特性の何れかであると判断されたときに、前記伝達関数推定手段で推定した伝達関数の次数を少なくとも１減らして伝達関数の係数を補正する伝達関数補正手段とを備えていることを特徴としている。
【００１４】
この請求項５に係る発明では、推定した伝達関数の安定性又は最小位相特性から補正の必要性を判断し、この判断結果に応じて、推定した伝達関数の次数を少なくとも１減らして伝達関数の係数を補正するので、システムの特性を推定するために十分な最小次数を有する伝達関数をより精度良く推定することができる。
【００１５】
【発明の実施の形態】
以下、本発明の実施の形態を図面に基づいて説明する。
図１は本発明の一実施形態を示す概略構成図であって、図中１は特性が不明なシステムであり、この特性不明システム１に入力信号ｕが入力されると、それに対して応答信号ｙを出力するように構成されている。
また、特性不明システム１の伝達関数を推定するための伝達関数推定装置２が備えられ、入力信号ｕ及び応答信号ｙのある一定期間のサンプリングデータｕ［ｎ］、ｙ［ｎ］を測定し、これらのデータの関係を表す差分方程式から導出した有理関数と、この有理関数をもとに得られる周波数応答データとから行列式を作成する。さらに、この行列の線形独立な列数を求め、線形独立な列数と同じ列数を持つ行列式が得られるように有理関数の次数を決定し、決定した次数とそのときの行列式の解とに基づいて伝達関数を推定する。そして、この伝達関数の安定性又は最小位相特性の判定を行い、判定結果に応じて、推定した伝達関数を補正し、システムの特性を推定するために十分な最小次数を有する伝達関数を推定する。
入力信号ｕは、複数の異なる周波数の正弦波を合成した多重正弦波信号であり、下記（１）式で表される。
【００１６】
【数１】
【００１７】
ここで、Ｍは合成する正弦波の数、ｆｋはｋ番目の正弦波の周波数成分、Ａｋはｋ番目の正弦波の周波数成分の振幅、Φｋはｋ番目の正弦波の周波数成分の位相である。ｋ番目の正弦波の周波数成分ｆｋは、伝達関数推定に必要な下限周波数成分の１周期をＴとすると、ｆｋ＝ｋ／Ｔである。また、合成正弦波の数Ｍは、伝達関数推定に必要な周波数領域の全域にわたる周波数成分を持つ正弦波を合成した多重正弦波信号とするために必要十分な数として設定し、ｋ番目の正弦波の周波数成分の位相Φｋは、合成正弦波の波高値が均一となるように設定する。
このように、複数の異なる周波数成分を含む多重正弦波信号を入力信号として適用することにより、特定の周波数成分を持つ正弦波信号を適用する場合と比較して、短時間で伝達関数の推定を行うことができる。
また、この入力信号ｕと応答信号ｙとの関係は、下記（２）式に示すような差分方程式で表すことができる。
【００１８】
【数２】
【００１９】
ここで、ｑは伝達関数の次数、ａｉ、ｂｉは差分方程式の未知係数である。
図２は、本発明の第１の実施形態における、伝達関数推定装置２で実行される伝達関数推定処理を示すフローチャートであって、先ず、ステップＳ１で、前記（１）式で表される入力信号と、それに対する応答信号のサンプリングデータｕ［ｎ］及びｙ［ｎ］を測定する。
次いでステップＳ２に移行して、入力信号ｕ［ｎ］と応答信号ｙ［ｎ］との関係を下記（３）式に示すように差分方程式で表す。ここで、ｐは初期設定次数である。
【００２０】
【数３】
【００２１】
そして、前記（３）式の関係を、次式に示すｚ演算子の多項式からなる有理関数に変換する。
【００２２】
【数４】
次にステップＳ３で、前記（４）式をもとに、多重正弦波信号ｕ（ｔ）の周波数帯域内の周波数応答Ｈ（ｚｋ）を算出する。ここで、ｚｋ＝ｅｘｐ（ｊ２πｆｋ／ｆｓ）、ｆｓ＝１／Δｔ、Δｔはサンプリング時間である。この周波数応答Ｈ（ｚｋ）をｍ点算出し、下記（５）式に示す周波数応答列ｄを得る。
ｄ＝［Ｈ（ｚ１），Ｈ（ｚ２），…，Ｈ（ｚｍ）］Ｔ ………（５）
次にステップＳ４で、前記（４）式の有理関数をもとに、周波数応答Ｈ（ｚｋ）に基づいて、下記（６）式に示すように、係数ａ’ｉ（ｉ＝１，…，ｐ）、ｂ’ｉ（ｉ＝０，…，ｐ）を変数ｘとした行列式Ａｘ＝ｄを作成する。
【００２３】
【数５】
【００２４】
次にステップＳ５で、行列Ａの線形独立な列数ＲＡ、即ちｒａｎｋＡを求め、ステップＳ６に移行し、行列Ａの線形独立な列数ＲＡが行列Ａの列数Ｒと等しいか否かを判定する。ＲＡ＜Ｒであるときには、ステップＳ７に移行して有理関数の次数ｐをデクリメントして前記ステップＳ４に移行する。
一方、ステップＳ６の判定結果が、ＲＡ＝Ｒであるときには、ｒａｎｋＡはフルランクであると判断してステップＳ８に移行する。
【００２５】
ステップＳ８では、ＲＡ＝Ｒであるときの次数ｐと、このときの行列式Ａｘ＝ｄとに基づいて解ｘ即ち係数ａ’ｉ、ｂ’ｉを算出してステップＳ９に移行し、算出した解ｘ即ち係数ａ’ｉ、ｂ’ｉから下記（７）式に示す伝達関数候補Ｈ’（ｚ）を求める。
Ｈ’（ｚ）＝（ｂ’０＋ｂ’１ｚ＋…＋ｂ’ｐｚｐ）／（１＋ａ’１ｚ＋…＋ａ’ｐｚｐ） ………（７）
次にステップＳ１０で、前記ステップＳ９で推定した伝達関数候補Ｈ’（ｚ）をもとにシステムが安定であるか、又は最小位相特性を満足するか否かを判定し、伝達関数候補Ｈ’（ｚ）で表されるシステムが安定でない又は最小位相特性を満足しないときには、前記ステップＳ７に移行する。
【００２６】
一方、ステップＳ１０の判定結果が、上記条件を満足しているときには、伝達関数候補Ｈ’（ｚ）を補正する必要はないと判断してステップＳ１１に移行し、前記ステップＳ９で求めた伝達関数候補Ｈ’（ｚ）を、システムの特性を推定するために十分な最小次数を有する伝達関数Ｈ（ｚ）として設定してから伝達関数推定処理を終了する。
この図２の処理において、ステップＳ２〜Ｓ４の処理が行列式生成手段に対応し、ステップＳ５〜Ｓ７の処理が次数決定手段に対応し、ステップＳ８の処理が係数決定手段に対応し、ステップＳ１０の処理が伝達関数補正判断手段に対応している。
【００２７】
したがって、今、特性が不明なシステムに入力信号としてＭ個の異なる周波数の正弦波を合成した多重正弦波信号ｕ（ｔ）を入力したとき、それに対する応答信号ｙ（ｔ）を出力するものとする。この場合には、先ず入力信号と応答信号のサンプリングデータｕ［ｎ］、ｙ［ｎ］を測定し、図２の伝達関数推定処理において、ステップＳ２で入力信号と応答信号との関係を表す差分方程式の次数を予め設定した初期設定次数ｐとして、前記（４）式をもとに有理関数Ｈ（ｚ）を導出し、ステップＳ３で周波数応答データ列ｄを求める。
【００２８】
次にステップＳ４で周波数応答データ列ｄに基づいて、次数ｐに応じて前記（６）式に示す行列式Ａｘ＝ｄを作成し、ステップＳ６で行列Ａの線形独立な列数ＲＡが行列Ａの列数Ｒと等しいか否かを判定する。線形独立な列数ＲＡが行列Ａの列数より少ないときには、ステップＳ７に移行して次数ｐをデクリメントしてからステップＳ４に移行し、初期設定次数ｐより１小さい次数で前記（６）式に示す行列式Ａｘ＝ｄを作成する。そして、ステップＳ６でこの行列Ａの線形独立な列数ＲＡが行列Ａの列数と等しいか否かを判定する。
【００２９】
このように、線形独立な列数ＲＡと行列の列数Ｒが等しくなるまで、即ちｒａｎｋＡがフルランクとなるまでステップＳ４からステップＳ７までの処理を繰り返し、ＲＡ＝ＲとなったときにステップＳ６からステップＳ８に移行して、そのときの次数ｐでの行列式Ａｘ＝ｄの解ｘ、即ち未知係数ａ’ｉ、ｂ’ｉを算出し、伝達関数候補Ｈ’（ｚ）を推定する。
【００３０】
次に、ステップＳ１０で、推定した伝達関数候補Ｈ’（ｚ）をもとにシステムが安定であるか、又は最小位相特性を満足するか否かを判定し、システムが不安定又は最小位相特性を満足しないときには、ステップＳ７に移行して次数ｐをデクリメントしてからステップＳ４に移行し、設定次数ｐより１小さい次数で前記（６）式に示す行列式Ａｘ＝ｄを作成する。
【００３１】
このように、伝達関数候補Ｈ’（ｚ）が安定であるか、又は最小位相特性を満足するまでステップＳ４からステップＳ１０までの処理を繰り返し、上記条件が成立したときに、ステップＳ１０からステップＳ１１に移行して、伝達関数候補Ｈ’（ｚ）を伝達関数Ｈ（ｚ）として設定することにより、特性不明システム１の特性を推定するために十分な最小次数を有する伝達関数が求まる。
【００３２】
したがって、上記第１の実施形態では、複数の異なる周波数成分を含む多重正弦波信号を入力信号として適用することにより、システムに必要以上の直流量を入力することなく、特定の周波数成分を持つ正弦波信号を適用する場合と比較して、短時間で伝達関数の推定を行うことができる。
また、周波数応答データに基づいて、伝達関数の未知係数を変数とした行列式を作成し、この行列式から伝達関数の次数及び係数を推定するので、複数個の伝達関数に対する相互評価を必要とせず、システムの入力信号と過渡応答信号から直接伝達関数を推定することができる。
【００３３】
さらに、周波数応答データに基づいて作成したデータ行列の線形独立な列数が、データ行列の列数と等しくなるように、即ちデータ行列のランクがフルランクとなるように伝達関数の次数を決定し、そのときの次数に応じて伝達関数の係数を決定するので、システムの特性を推定するために十分な最小次数を有する伝達関数を推定することができる。
【００３４】
また、推定した伝達関数候補に対して安定であるか又は最小位相特性を満足するか否かの判定を行い、伝達関数候補が安定でない又は最小位相特性を満足しないと判断したときには、推定した伝達関数候補の次数を少なくとも１減らして伝達関数候補の係数を補正し、最終的な伝達関数を設定するので、システムの特性を推定するために十分な最小次数を有する伝達関数をより精度良く推定することができる。
次に、本発明の第２の実施形態について説明する。
この第２の実施形態は、前述した第１の実施形態において、周波数応答データ列に基づいて生成した行列式に重み係数行列を乗算し、これにより得られる重み付き行列式に基づいて、伝達関数の次数及び係数を決定するようにしたものである。
【００３５】
図３は、第２の実施形態において、伝達関数推定装置２で実行される伝達関数推定処理の処理手順を示すフローチャートであって、図２に示す第１の実施形態における伝達関数推定処理において、ステップＳ４の後に、ステップＳ４で生成した行列式に重み係数行列Ｗを乗算して重み付き行列式Ｃｘ＝ｇを生成するステップＳ２１が追加され、ステップＳ５がデータ行列ＣのランクＲＣを算出するステップＳ２２に置換され、ステップＳ６がランクＲＣと行列Ｃの列数Ｒとが等しいか否かを判定するステップＳ２３に置換され、ステップＳ８が、前記ステップＳ２３でＲＣ＝Ｒであると判断したときに、重み付き行列式Ｃｘ＝ｇに基づいて伝達関数の次数及び係数を決定するステップＳ２４に置換されていることを除いては図２と同様の処理を行い、図２と同一部には同一符号を付与しその詳細な説明は省略する。
【００３６】
この図３の処理において、ステップＳ２１の処理が重み付け手段に対応している。
この第２の実施形態によると、ステップＳ２１で、前記ステップＳ４で生成した行列式Ａｘ＝ｄに重み係数行列Ｗを乗算し、重み付き行列式Ｃｘ＝ｇを生成する。
Ｃ＝ＷＡ ………（８）
ｇ＝Ｗｄ ………（９）
Ｗ＝ｄｉａｇ（ｗ１，ｗ２，…，ｗｍ） ………（１０）
ここで、ｗｉ（ｉ＝１，…，ｍ）は重み係数、ｄｉａｇ（ ）は括弧内の値を対角要素とする対角行列を表す関数である。
【００３７】
次にステップＳ２２で、行列Ｃの線形独立な列数ＲＣ、即ちｒａｎｋＣを求め、ステップＳ２３に移行し、行列Ｃの線形独立な列数ＲＣが行列Ｃの列数Ｒと等しいか否かを判定する。ＲＣ＜Ｒであるときには、前記ステップＳ７に移行する。
一方、ステップＳ２３の判定結果が、ＲＣ＝Ｒであるときには、ｒａｎｋＣはフルランクであると判断してステップＳ２４に移行し、ＲＣ＝Ｒであるときの次数ｐと、このときの重み付き行列式Ｃｘ＝ｇとに基づいて解ｘを算出し、算出した解ｘ即ち未知係数ａ’ｉ、ｂ’ｉから前記（７）式をもとに伝達関数候補Ｈ’（ｚ）を推定する。
【００３８】
このように、上記第２の実施形態では、周波数応答データに重み係数を乗算し、これにより得られる重み付き行列式に基づいて伝達関数の次数及び係数を決定するので、伝達関数の推定誤差を低減することができる。
なお、上記各実施形態においては、ステップＳ１０で、推定した伝達関数候補に対して伝達関数補正判断を行う場合について説明したが、これに限定されるものではなく、伝達関数補正判断を行わずに伝達関数候補を直接システムの特性を推定するために十分な最小次数を有する伝達関数として設定するようにしてもよい。
また、上記各実施形態においては、ステップＳ７で次数を１ずつ減らして最小次数を推定する場合について説明したが、これに限定されるものではなく、１以上初期設定次数未満の値ずつ減らして最小次数を推定するようにしてもよい。
【００３９】
【発明の効果】
以上説明したように、本発明の請求項１に係る発明によれば、複数の異なる周波数の正弦波を合成した多重正弦波信号を入力信号として適用するので、システムに必要以上の直流量を入力することなく、短期間の測定でシステムの特性を近似する伝達関数を推定することができると共に、有理関数の周波数応答データに基づいて生成した行列式から伝達関数を推定するので、複数個の伝達関数に対する相互評価を必要とせず、システムの入力信号と過渡応答信号から直接伝達関数を推定することができるという効果が得られる。
【００４０】
また、請求項２に係る発明によれば、入力信号と応答信号との関係を表す有理関数とそれに対する入力信号の周波数帯域内の周波数応答データに基づいて行列式を生成し、その行列式から伝達関数を推定するので、複数個の伝達関数に対する相互評価を必要とせず、システムの入力信号と過渡応答信号から直接伝達関数を推定することができるという効果が得られる。
【００４１】
さらに、請求項３に係る発明によれば、周波数応答データに重み係数を乗算して得られる行列式から伝達関数を推定するので、伝達関数の推定誤差を低減することができるという効果が得られる。
また、請求項４に係る発明によれば、データ行列の線形独立な列数がデータ行列の列数と等しくなるように伝達関数の次数を決定し、決定した次数に応じて伝達関数の係数を決定するので、システムの特性を推定するために十分な最小次数を有する伝達関数を推定することができるという効果が得られる。
【００４２】
さらにまた、請求項５に係る発明によれば、推定した伝達関数の安定性又は最小位相特性から補正の必要性を判断し、この判断結果に応じて、推定した伝達関数の次数を少なくとも１減らして伝達関数の係数を補正するので、システムの特性を推定するために十分な最小次数を有する伝達関数をより精度良く推定することができるという効果が得られる。
【図面の簡単な説明】
【図１】本発明の実施形態を示す概略構成図である。
【図２】本発明の第１の実施形態における図１の伝達関数推定装置２の伝達関数推定処理を示すフローチャートである。
【図３】本発明の第２の実施形態における図１の伝達関数推定装置２の伝達関数推定処理を示すフローチャートである。
【符号の説明】
１ 特定不明システム
２ 伝達関数推定装置[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a transfer function estimation device for a system with unknown characteristics, which estimates a transfer function of the system from an input signal and a response signal for the system with unknown characteristics.
[0002]
[Prior art]
As a transfer function estimation device for a conventional system with unknown characteristics, an impulse response sequence is obtained using a step response waveform when a stepped input signal is applied, and a discrete-time state space model calculated based on the impulse response sequence is continuously used. It is known that a transfer function is calculated by converting to a system state space model and estimating the order of the continuous system state space model from the rank of a Hankel matrix composed of impulse response sequences (for example, Patent Document 1). reference).
[0003]
As a system identification method, a low-order linear model in the rational function form of the Laplace operator s is identified by applying the least square method to existing frequency response data, and the transfer function is estimated. Is known to perform a characteristic approximation accuracy evaluation with a standard deviation of output error (relative error evaluation) for the calculated least squares solution and select a transfer function with the best evaluation (for example, non-patent) Reference 1).
[0004]
[Patent Document 1]
Japanese Patent Publication No. 1994-77211 (2nd page, Fig. 1)
[Non-Patent Document 1]
"Method for identifying low-order linear model of system by frequency domain least square approximation", IEEJ Transactions B, 2001, Vol. 121, No. 1, p52-59
[0005]
[Problems to be solved by the invention]
However, in the conventional example described in Patent Document 1, since the step input is applied to the input signal, the DC amount is input to the system, and the application range is limited. There is.
Also, when estimating the transfer function, a continuous state space model is obtained from the impulse response sequence, and the transfer function is calculated by estimating the order of the continuous state space model from the rank of the Hankel matrix composed of the impulse response sequence. Therefore, there is an unsolved problem that the transfer function cannot be estimated directly from the input signal and output response signal of the system.
[0006]
Furthermore, in the conventional example described in Non-Patent Document 1, in order to estimate the transfer function, it is assumed that frequency response data for the frequency region necessary for estimation is prepared in advance, and the characteristics are There is an unresolved issue that cannot be applied to unknown systems.
In addition, when estimating the transfer function, the characteristic approximation accuracy is evaluated using the standard deviation of output error (relative error evaluation), and the transfer function with the best evaluation is selected, so mutual evaluation of multiple assumed transfer functions is required. There is an unresolved issue.
Therefore, the present invention has been made paying attention to the above-mentioned unsolved problems of the conventional example, and a transfer function estimation device for a system with unknown characteristics that can directly estimate the transfer function from the input signal and transient response signal of the system. The purpose is to provide.
[0007]
[Means for Solving the Problems]
In order to achieve the above object, a transfer function estimation device for a system with unknown characteristics according to claim 1 is a relationship between an input signal composed of multiple sine waves obtained by synthesizing a plurality of sine waves of different frequencies and a response signal thereto. A determinant generating means for generating a determinant with a variable of an unknown coefficient of the rational function based on a rational function representing the rational function and the frequency response data of the rational function, and transmitting from the determinant generated by the determinant generating means And a transfer function estimating means for estimating the function.
[0008]
In the invention according to claim 1, since a multiple sine wave signal obtained by synthesizing a plurality of sine waves of different frequencies is applied as an input signal, the system can be measured in a short period of time without inputting an excessive amount of direct current to the system. A transfer function approximating the characteristics of can be estimated. Since the transfer function is estimated from the determinant generated based on the frequency response data of the input signal, the transfer function is estimated directly from the input signal and transient response signal of the system without the need for mutual evaluation of multiple transfer functions. can do.
[0009]
According to a second aspect of the present invention, the determinant generating means is configured to determine the relationship between the input signal and the response signal in accordance with a difference equation of a predetermined order. Determining the difference equation as a rational function consisting of a polynomial of the z operator, and using the rational function and the frequency response data in the frequency band of the input signal to the determinant with the unknown coefficient of the rational function as a variable It is characterized by being configured to generate.
[0010]
In the invention according to claim 2, a determinant is generated based on a rational function representing the relationship between the input signal and the response signal and frequency response data in the frequency band of the input signal corresponding to the rational function, and the transfer function is calculated from the determinant. Since the estimation is performed, the transfer function can be estimated directly from the input signal of the system and the transient response signal without requiring mutual evaluation for a plurality of transfer functions.
[0011]
Further, in the invention according to claim 1 or 2, the determinant generation means uses the unknown coefficient of the rational function as a variable based on the frequency response data. A weighting means for generating a weighted determinant by multiplying the generated determinant by a weighting coefficient is provided.
In the invention according to the third aspect, since the transfer function is estimated from the determinant obtained by multiplying the frequency response data by the weight coefficient, the estimation error of the transfer function can be reduced.
[0012]
According to a fourth aspect of the present invention, there is provided a transfer function estimation device for an unknown characteristic system according to any one of the first to third aspects, wherein the transfer function estimation means includes an order determination means for determining the order of the transfer function, and a transfer function. Coefficient determining means for determining the coefficient of the function, wherein the order determining means includes a transfer function so that the number of linearly independent columns of the matrix generated by the determinant generating means is equal to the number of columns of the matrix. The order is determined, and the coefficient determination unit is configured to determine a coefficient of a transfer function from a determinant generated by the determinant generation unit according to the order determined by the order determination unit. It is said.
[0013]
In the invention according to claim 4, the order of the transfer function is determined so that the number of linearly independent columns of the data matrix is equal to the number of columns of the data matrix, and the coefficient of the transfer function is determined according to the determined order. A transfer function with a minimum order sufficient to estimate the characteristics of the system can be estimated.
According to a fifth aspect of the present invention, there is provided a transfer function estimation device for a system with unknown characteristics according to any one of the first to fourth aspects, wherein either the stability of the transfer function or the minimum phase characteristic estimated by the transfer function estimation means is used. A transfer function correction determining means for determining the necessity of correcting the transfer function estimated by the transfer function estimating means, and the transfer function is determined to be either unstable or non-minimum phase characteristic by the transfer function correction determining means. The transfer function correcting unit corrects the coefficient of the transfer function by reducing the order of the transfer function estimated by the transfer function estimating unit by at least 1 when judged.
[0014]
In the invention according to claim 5, the necessity of correction is determined from the stability or minimum phase characteristic of the estimated transfer function, and the order of the estimated transfer function is reduced by at least 1 according to the determination result. Since the coefficients are corrected, a transfer function having a minimum order sufficient to estimate the characteristics of the system can be estimated with higher accuracy.
[0015]
DETAILED DESCRIPTION OF THE INVENTION
Hereinafter, embodiments of the present invention will be described with reference to the drawings.
FIG. 1 is a schematic configuration diagram showing an embodiment of the present invention. In FIG. 1, reference numeral 1 denotes a system with unknown characteristics. When an input signal u is input to the unknown characteristics system 1, a response signal is sent to the system. It is configured to output y.
In addition, a transfer function estimation device 2 for estimating the transfer function of the characteristic unknown system 1 is provided, and sampling data u [n], y [n] of a certain period of the input signal u and the response signal y are measured, A determinant is created from a rational function derived from a differential equation representing the relationship between these data and frequency response data obtained based on this rational function. Furthermore, the number of linearly independent columns of this matrix is obtained, the order of the rational function is determined so that a determinant having the same number of columns as the number of linearly independent columns is obtained, and the determined order and the solution of the determinant at that time are determined. Based on the above, the transfer function is estimated. Then, the stability of the transfer function or the minimum phase characteristic is determined, and the estimated transfer function is corrected according to the determination result, and a transfer function having a minimum order sufficient to estimate the system characteristics is estimated. .
The input signal u is a multiple sine wave signal obtained by synthesizing a plurality of sine waves having different frequencies, and is represented by the following equation (1).
[0016]
[Expression 1]
[0017]
Here, M is the number of sine waves to be combined, f k is the frequency component of the kth sine wave, A k is the amplitude of the frequency component of the kth sine wave, and Φ k is the frequency component of the kth sine wave. It is a phase. The frequency component f k of the kth sine wave is f k = k / T, where T is one period of the lower limit frequency component necessary for the transfer function estimation. Further, the number M of combined sine waves is set as a necessary and sufficient number to obtain a multiple sine wave signal obtained by combining sine waves having frequency components over the entire frequency region necessary for the transfer function estimation. The phase Φ k of the frequency component of the wave is set so that the peak value of the synthesized sine wave is uniform.
In this way, by applying a multiple sine wave signal including a plurality of different frequency components as an input signal, the transfer function can be estimated in a shorter time than when a sine wave signal having a specific frequency component is applied. It can be carried out.
The relationship between the input signal u and the response signal y can be expressed by a difference equation as shown in the following equation (2).
[0018]
[Expression 2]
[0019]
Here, q is the order of the transfer function, and a i and b i are unknown coefficients of the difference equation.
FIG. 2 is a flowchart showing a transfer function estimation process executed by the transfer function estimation device 2 in the first embodiment of the present invention. First, in step S1, the input represented by the above equation (1) is shown. The sampling data u [n] and y [n] of the signal and the corresponding response signal are measured.
Next, the process proceeds to step S2, and the relationship between the input signal u [n] and the response signal y [n] is expressed by a difference equation as shown in the following equation (3). Here, p is an initial set order.
[0020]
[Equation 3]
[0021]
Then, the relationship of the expression (3) is converted into a rational function composed of a polynomial of the z operator shown in the following expression.
[0022]
[Expression 4]
In step S3, a frequency response H (z k ) in the frequency band of the multiple sine wave signal u (t) is calculated based on the equation (4). Here, z k = exp (j2πf k / f s ), f s = 1 / Δt, and Δt are sampling times. This frequency response H (z k ) is calculated at m points to obtain a frequency response sequence d shown in the following equation (5).
d = [H (z 1 ), H (z 2 ),..., H (z m )] T (5)
Next, in step S4, based on the rational function of the equation (4), based on the frequency response H (z k ), as shown in the following equation (6), coefficients a ′ i (i = 1,... , P), b ′ i (i = 0,..., P) is created as a determinant Ax = d.
[0023]
[Equation 5]
[0024]
In step S5, linearly independent columns R A of the matrix A, that seek Ranka, the process proceeds to step S6, whether linearly independent columns R A of the matrix A is equal to the number of columns R of the matrix A Determine. When R A <R, the process proceeds to step S7, the order p of the rational function is decremented, and the process proceeds to step S4.
On the other hand, when the determination result of step S6 is R A = R, it is determined that rank A is full rank, and the process proceeds to step S8.
[0025]
In step S8, a solution x, that is, coefficients a ′ i and b ′ i are calculated based on the order p when R A = R and the determinant Ax = d at this time, and the process proceeds to step S9. A transfer function candidate H ′ (z) shown in the following equation (7) is obtained from the solution x, that is, the coefficients a ′ i and b ′ i .
H ′ (z) = (b ′ 0 + b ′ 1 z +... + B ′ p z p ) / (1 + a ′ 1 z +... + A ′ p z p ) (7)
Next, in step S10, it is determined whether the system is stable or satisfies the minimum phase characteristic based on the transfer function candidate H ′ (z) estimated in step S9, and the transfer function candidate H ′. When the system represented by (z) is not stable or does not satisfy the minimum phase characteristic, the process proceeds to step S7.
[0026]
On the other hand, if the determination result in step S10 satisfies the above condition, it is determined that there is no need to correct the transfer function candidate H ′ (z), the process proceeds to step S11, and the transfer function obtained in step S9 is determined. The candidate H ′ (z) is set as a transfer function H (z) having a minimum order sufficient to estimate the characteristics of the system, and then the transfer function estimation process is terminated.
In the processing of FIG. 2, the processing of steps S2 to S4 corresponds to the determinant generation means, the processing of steps S5 to S7 corresponds to the order determination means, the processing of step S8 corresponds to the coefficient determination means, and step S10 This process corresponds to the transfer function correction determination means.
[0027]
Therefore, when a multiple sine wave signal u (t) obtained by synthesizing M sine waves of different frequencies is input as an input signal to a system whose characteristics are unknown, a response signal y (t) is output as a response signal y (t). To do. In this case, first, sampling data u [n], y [n] of the input signal and the response signal are measured, and in the transfer function estimation process of FIG. 2, a difference representing the relationship between the input signal and the response signal in step S2. A rational function H (z) is derived on the basis of the above equation (4) with the order of the equation set as an initial setting order p, and a frequency response data string d is obtained in step S3.
[0028]
Next, in step S4, based on the frequency response data string d, a determinant Ax = d shown in the above equation (6) is created according to the order p, and in step S6, the linearly independent column number R A of the matrix A is a matrix. It is determined whether or not the number of columns R of A is equal. When the number of linearly independent columns RA is smaller than the number of columns of the matrix A, the process proceeds to step S7, the order p is decremented, and then the process proceeds to step S4. The determinant Ax = d shown in FIG. In step S6, it is determined whether or not the number of linearly independent columns RA of the matrix A is equal to the number of columns of the matrix A.
[0029]
In this way, the process from step S4 to step S7 is repeated until R A = R until the linearly independent column number R A and the column number R of the matrix become equal, that is, rank A reaches the full rank. The process proceeds from step S6 to step S8, and the solution x of the determinant Ax = d at the order p at that time, that is, the unknown coefficients a ′ i and b ′ i are calculated, and the transfer function candidate H ′ (z) is estimated. To do.
[0030]
Next, in step S10, it is determined whether the system is stable or satisfies the minimum phase characteristic based on the estimated transfer function candidate H ′ (z), and the system is unstable or has the minimum phase characteristic. If not satisfied, the process proceeds to step S7, the order p is decremented, then the process proceeds to step S4, and the determinant Ax = d shown in the above equation (6) is created with an order smaller than the set order p.
[0031]
In this way, the processes from step S4 to step S10 are repeated until the transfer function candidate H ′ (z) is stable or satisfies the minimum phase characteristic, and when the above condition is satisfied, step S10 to step S11 are performed. Then, the transfer function candidate H ′ (z) is set as the transfer function H (z), whereby a transfer function having a minimum order sufficient for estimating the characteristics of the characteristic unknown system 1 is obtained.
[0032]
Therefore, in the first embodiment, by applying a multiple sine wave signal including a plurality of different frequency components as an input signal, a sine having a specific frequency component can be obtained without inputting an excessive DC amount to the system. Compared with the case where a wave signal is applied, the transfer function can be estimated in a short time.
Also, based on the frequency response data, a determinant is created with the unknown coefficient of the transfer function as a variable, and the order and coefficient of the transfer function are estimated from this determinant. Therefore, mutual evaluation of multiple transfer functions is required. Instead, the transfer function can be estimated directly from the system input signal and the transient response signal.
[0033]
Further, the order of the transfer function is determined so that the number of linearly independent columns of the data matrix created based on the frequency response data is equal to the number of columns of the data matrix, that is, the rank of the data matrix is full rank. Since the coefficient of the transfer function is determined according to the order at that time, it is possible to estimate the transfer function having the minimum order sufficient for estimating the characteristics of the system.
[0034]
Further, it is determined whether the estimated transfer function candidate is stable or satisfies the minimum phase characteristic, and when it is determined that the transfer function candidate is not stable or does not satisfy the minimum phase characteristic, the estimated transfer function Since the function candidate order is reduced by at least 1 and the coefficients of the transfer function candidate are corrected and the final transfer function is set, the transfer function having the minimum order sufficient to estimate the characteristics of the system is estimated more accurately. be able to.
Next, a second embodiment of the present invention will be described.
In the second embodiment, the determinant generated based on the frequency response data sequence is multiplied by the weighting coefficient matrix in the first embodiment, and the transfer function is calculated based on the weighted determinant obtained thereby. The order and coefficient are determined.
[0035]
FIG. 3 is a flowchart showing a processing procedure of a transfer function estimation process executed by the transfer function estimation device 2 in the second embodiment. In the transfer function estimation process in the first embodiment shown in FIG. after step S4, step S21 is added for multiplying the weighting factor matrix W to a determinant generated in step S4 generates a weighting matrix equation Cx = g, step S5 calculates the rank R C of the data matrix C It is replaced with step S22, and step S6 is replaced with step S23 for determining whether rank RC is equal to the number of columns R of matrix C, and step S8 determines that RC = R in step S23. 2 is replaced with step S24 for determining the order and coefficient of the transfer function based on the weighted determinant Cx = g. The same reference numerals are given to the same parts as those in FIG. 2, and the detailed description thereof is omitted.
[0036]
In the process of FIG. 3, the process of step S21 corresponds to the weighting means.
According to the second embodiment, in step S21, the determinant Ax = d generated in step S4 is multiplied by the weight coefficient matrix W to generate a weighted determinant Cx = g.
C = WA (8)
g = Wd ......... (9)
W = diag (w 1 , w 2 ,..., W m ) (10)
Here, w i (i = 1,..., M) is a weighting coefficient, and diag () is a function representing a diagonal matrix having the values in parentheses as diagonal elements.
[0037]
In step S22, linearly independent columns R C matrix C, that seek RankC, the process proceeds to step S23, whether linearly independent columns R C matrix C is equal to the number of columns R of the matrix C Determine. When R C <R, the process proceeds to step S7.
On the other hand, when the determination result of step S23 is R C = R, it is determined that rankC is full rank, and the process proceeds to step S24. The order p when R C = R and the weighting at this time are determined. A solution x is calculated based on the determinant Cx = g, and a transfer function candidate H ′ (z) is estimated from the calculated solution x, that is, unknown coefficients a ′ i and b ′ i based on the equation (7). .
[0038]
Thus, in the second embodiment, the frequency response data is multiplied by the weighting coefficient, and the order and coefficient of the transfer function are determined based on the weighted determinant thus obtained. Can be reduced.
In each of the above embodiments, the case where the transfer function correction determination is performed on the estimated transfer function candidate in step S10 has been described. However, the present invention is not limited to this, and the transfer function correction determination is not performed. The transfer function candidate may be directly set as a transfer function having a minimum order sufficient to estimate the characteristics of the system.
In each of the above embodiments, the case has been described in which the minimum order is estimated by reducing the order by 1 in step S7. However, the present invention is not limited to this, and the minimum is reduced by one or more values less than the initial setting order. The order may be estimated.
[0039]
【The invention's effect】
As described above, according to the first aspect of the present invention, a multiple sine wave signal obtained by synthesizing a plurality of sine waves having different frequencies is applied as an input signal. Transfer function that approximates the characteristics of the system by short-term measurement, and the transfer function is estimated from the determinant generated based on the frequency response data of the rational function. There is an effect that the transfer function can be estimated directly from the input signal of the system and the transient response signal without requiring mutual evaluation of the function.
[0040]
According to the invention of claim 2, the determinant is generated based on the rational function representing the relationship between the input signal and the response signal and the frequency response data in the frequency band of the input signal corresponding to the rational function. Since the transfer function is estimated, there is an effect that the transfer function can be estimated directly from the input signal of the system and the transient response signal without requiring mutual evaluation for a plurality of transfer functions.
[0041]
Furthermore, according to the third aspect of the invention, since the transfer function is estimated from the determinant obtained by multiplying the frequency response data by the weighting coefficient, the effect that the estimation error of the transfer function can be reduced is obtained. .
According to the invention of claim 4, the order of the transfer function is determined so that the number of linearly independent columns of the data matrix is equal to the number of columns of the data matrix, and the coefficient of the transfer function is set according to the determined order. As a result, it is possible to estimate a transfer function having a minimum order sufficient to estimate the characteristics of the system.
[0042]
Furthermore, according to the invention of claim 5, the necessity of correction is determined from the stability or minimum phase characteristic of the estimated transfer function, and the estimated transfer function order is reduced by at least 1 according to the determination result. Thus, the transfer function coefficient is corrected, so that the transfer function having the minimum order sufficient for estimating the system characteristics can be estimated with higher accuracy.
[Brief description of the drawings]
FIG. 1 is a schematic configuration diagram showing an embodiment of the present invention.
FIG. 2 is a flowchart showing a transfer function estimation process of the transfer function estimation device 2 of FIG. 1 in the first embodiment of the present invention.
FIG. 3 is a flowchart showing a transfer function estimation process of the transfer function estimation device 2 of FIG. 1 in the second embodiment of the present invention.
[Explanation of symbols]
1 Unspecified system 2 Transfer function estimation device
