(57)【要約】【課題】負荷外力の存在するシステムのパラメータを精確に同定することのできるパラメータ同定方法を提供する。【解決手段】本発明の同定方法は、入力が駆動モータに対する目標値と外乱であって、出力が駆動モータの回転に関する２つの情報量となるような同定対象閉ループモデルで表現しておき、同定対象閉ループモデルは、パラメータを備える動特性を有すると共に、パラメータを備える動特性の出力側と外乱が入力される点との間から情報量が出力され、外乱が入力される点とパラメータの動特性の入力側との間から、情報量が出力されるように構成しておいて、同定対象閉ループモデルを構成するパルス伝達関数が、応答を規定するパラメータで表現されるものとなっており、同定対象モデルへの入力値とパルス伝達関数とで構成される誤差関数が最小値を取るように、応答を規定するパラメータを同定する。【選択図】図２０
2015111384.tif 000002
【特許請求の範囲】
【請求項１】
駆動モータによって構成され且つ負荷外力が付与される慣性系のシステムにおける応答を規定するパラメータを同定するパラメータ同定方法において、
外乱と、パラメータの変化量とが分離可能である際には、パラメータが同定可能であることを特徴とするシステムのパラメータ同定方法。
【請求項２】
前記慣性系のシステムは、パラメータα1を有するものとされ、
前記慣性系のシステムを、入力が駆動モータに対する目標値uと外乱wであって、出力が前記駆動モータの回転に関する２つの情報量y1、y2となるような同定対象閉ループモデルで表現しておき、
前記同定対象閉ループモデルは、パラメータα1を備える動特性を有すると共に、パラメータα1を備える動特性の出力側と外乱wが入力される点との間から、情報量y2が出力され、外乱wが入力される点とパラメータα1の動特性の入力側との間から、情報量y1が出力されるように構成しておいて、
前記同定対象閉ループモデルを構成するパルス伝達関数が、前記応答を規定するパラメータα1で表現されるものとなっており、  
前記同定対象モデルへの入力値とパルス伝達関数とで構成される誤差関数が最小値を取るように、前記応答を規定するパラメータα1を同定することを特徴とする請求項１に記載のシステムのパラメータ同定方法。 
【請求項３】
前記応答を規定するパラメータα1は、駆動モータに取り付けられた電流計測センサのゲイン特性であることを特徴とする請求項２に記載のシステムのパラメータ同定方法。
【請求項４】
前記慣性系のシステムにコントローラの出力値が入力され、且つ前記コントローラの入力に目標値と慣性系のシステムの出力値がフィードバックされる制御系においては、次式を満たす際に、慣性系のシステムのパラメータが同定可能とされていることを特徴とする請求項１〜３のいずれかに記載のシステムのパラメータ同定方法。
【数４９】
2015111384.tif 000053
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、回転機器などから構成されるシステムの入出力データに基づくパラメータ同定方法に関する。
【背景技術】
【０００２】
近年、環境問題や節電への意識が高まり、今まで以上に省エネルギーへの要求が強くなっている。これにともない、既存の制御系のエネルギー効率の劣化を検出するために、エネルギー効率を評価することもより重要な課題となっている。そのためには、稼働中の制御系の入出力データを用いてシステムのパラメータを同定することが有効と考えられる。
例えば、工場等で使用されるモータや油圧ポンプなどの回転機器に目を向けると、少ないエネルギーで高い効率で稼働させることを目指し、例えば、回転機器の各種パラメータ（例えば、摩擦係数）などを正確に知る必要が生じている。
【０００３】
あるシステムに関して、当該システムのパラメータを同定する手法としては、特許文献１に開示されたものがある。
特許文献１は、油圧サーボ弁と、油圧サーボ弁により駆動される油圧シリンダと、油圧サーボ弁に制御信号を供給し油圧シリンダの変位を測定する制御部とを備えた油圧システムにおける油圧システムパラメータ同定方法において、油圧サーボ弁に制御信号を供給する入力チャンネルと、油圧サーボ弁のスプール変位又は油圧シリンダのピストン変位を検出する出力チャンネルと、油圧サーボ弁に印加される入力波形を設定する処理と、前記設定する処理により設定された入力波形を入力チャンネルにより油圧サーボ弁へ入力し、油圧サーボ弁のスプール変位又は油圧シリンダのピストン変位を出力チャンネルにより計測する処理と、前記計測する処理により計測された油圧サーボ弁のスプール変位又は油圧シリンダのピストン変位に基づいて、油圧サーボ弁のスプールの粘性減衰率又は油圧シリンダのピストンの摩擦特性又は油圧サーボ弁のゲイン等の油圧システムパラメータを同定する処理とを含む油圧システムパラメータ同定方法を開示する。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２００１−１１７６２７号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
特許文献１は、あるシステム（油圧システム）におけるパラメータ同定方法を開示するものの、同定対象となっているシステムは、負荷が外乱として加わるものとなっていない。現実に使用されるシステムには、負荷が加わることになり、この負荷を考慮しない状況下でのパラメータ同定は、正確なパラメータの値を求めることができない。
例えば、工場等で使用されるモータや油圧ポンプなどの回転機器を考えた場合、この回転機器の稼働中には、負荷側からの反力などに伴う負荷トルクが回転機器に作用し、外乱として働く。斯かる状況下において、回転機器の各種パラメータを同定することは、特許文献１の技術を用いても困難であると思われる。つまり、回転機器に対して負荷トルクが外乱として加わるようなシステムのパラメータを精度よく推定する技術は、未だ開発されるには至っていない。
【０００６】
そこで、本発明は、上記問題点を鑑み、負荷外力の存在するシステムのパラメータを精確に同定することのできるパラメータ同定方法を提供することを目的とする。
【課題を解決するための手段】
【０００７】
上述の目的を達成するため、本発明においては以下の技術的手段を講じた。
本発明に係るシステムのパラメータ同定方法は、駆動モータによって構成され且つ負荷外力が付与される慣性系のシステムにおける応答を規定するパラメータを同定するパラメータ同定方法において、外乱と、パラメータの変化量とが分離可能である際には、パラメータが同定可能であることを特徴とする。
【０００８】
好ましくは、前記慣性系のシステムは、パラメータα1を有するものとされ、前記慣性
系のシステムを、入力が駆動モータに対する目標値uと外乱wであって、出力が前記駆動モータの回転に関する２つの情報量y1、y2となるような同定対象閉ループモデルで表現しておき、前記同定対象閉ループモデルは、パラメータα1を備える動特性を有すると共に、パラメータα1を備える動特性の出力側と外乱wが入力される点との間から、情報量y2が出力され、外乱wが入力される点とパラメータα1の動特性の入力側との間から、情報量y1が出力されるように構成しておいて、前記同定対象閉ループモデルを構成するパルス伝達関数が、前記応答を規定するパラメータα1で表現されるものとなっており、前記同定対象モデルへの入力値とパルス伝達関数とで構成される誤差関数が最小値を取るように、前記応答を規定するパラメータα1を同定することを特徴とする。
【０００９】
好ましくは、 前記応答を規定するパラメータα1は、駆動モータに取り付けられた電流計測センサのゲイン特性であるとよい。
好ましくは、前記慣性系のシステムにコントローラの出力値が入力され、且つ前記コントローラの入力に目標値と慣性系のシステムの出力値がフィードバックされる制御系においては、次式を満たす際に、慣性系のシステムのパラメータが同定可能とされているとよい。
【００１０】
【数５０】
2015111384.tif 000003
【発明の効果】
【００１１】
本発明によれば、負荷外力の存在するシステムのパラメータを精確に同定することが可能となる。
【図面の簡単な説明】
【００１２】
【図１】同定対象モデルを示した図である。
【図２】一慣性系に対するフィードバック制御系を示した図である。
【図３】モータとオペアンプの等価回路を示した図である。
【図４】実機実験を行うためのＳｉｍｕｌｉｎｋモデルを示したものである。
【図５】αと評価関数との関係を示した図である。
【図６】シミュレーション応答を示した図である。
【図７】シミュレーションに基づく評価関数を示した図である。
【図８】外乱を０としたときの応答を示した図である。
【図９】外乱を０としたときの評価関数を示した図である。
【図１０】Ｋ２（ｓ）に対する応答を示した図である。
【図１１】シミュレーションに基づく評価関数を示した図である。
【図１２】個別に離散化した際の閉ループ系を示した図である。
【図１３】φに依存する出力を示した図である。
【図１４】モータ電機子電流をｙ２としたときの評価関数を示した図である。
【図１５】Ｆをインナー化したときの評価関数を示した図である。
【図１６】ｗとｑの変化を示した図である。
【図１７】ｗと１ステップシフトしたｑの変化を示した図である。
【図１８】ｗと２ステップシフトしたｑの変化を示した図である。
【図１９】ｗとシフトしたｑの誤差を示した図である。
【図２０】パラメータが同定可能となるための必要十分条件を検証するためのフィードバック制御系を示した図である。
【発明を実施するための形態】
【００１３】
以下、本発明の実施形態を、図を基に説明する。
本願発明は、駆動モータによって構成され且つ負荷外力（外乱）が付与される慣性系のシステムにおける応答を規定するパラメータα1を同定するパラメータ同定方法である。このパラメータ同定方法は、慣性系のシステムを、入力が駆動モータに対する目標値uと外乱wであって、出力が前記駆動モータの回転に関する２つの情報量y1、y2となるような
同定対象閉ループモデルで表現しておくものである。この同定対象閉ループモデルは、パラメータα1を備える動特性を有すると共に、パラメータα1を備える動特性の出力側と外乱wが入力される点との間から、情報量y2が出力され、外乱wが入力される点とパラメータα1の動特性の入力側との間から、情報量y1が出力されるように構成してある。この同定対象閉ループモデルを構成するパルス伝達関数は、応答を規定するパラメータα1で表現されるものとなっており、本発明のパラメータ同定方法は同定対象閉ループモデルへの入力値とパルス伝達関数とで構成される誤差関数が最小値を取るように、応答を規定するパラメータα1を同定することを特徴としている。
【００１４】
一般に同定はシステムの入出力関係を調べることで行われる。このとき、システムの入出力データをすべて把握できることが望ましい。しかしながら、通常は制御系には負荷外力が加わる。多くの場合、負荷外力を測定して入力データとして扱うことは困難である。負荷外力を避けるために対象となる制御系を停止し、同定のためのテストを行うことも可能であるが、その場合には設備の生産性が低下してしまう。
【００１５】
そこで、本実施形態では、複数の出力信号を考えることで負荷外力に依存しない信号成分を取り出し、その信号を用いてパラメータの同定を行う手法を開示する。この手法では、負荷外力の存在する状況下においても、システムが稼働している状態での入出力データを用いてパラメータを同定することが可能である。
以降では、「パラメータ同定手法の概略」、「実際のモータ制御系にパラメータ同定手法を適応した実施形態」の順で本実施形態の説明を行う。
［パラメータ同定手法の概要］
まず、初めに、本実施形態の核となるパラメータ同定手法の概要を説明しておく。
図１に示されるパラメータに依存する２入力２出力の閉ループの離散時間システム（(1)式）を考える。
【００１６】
【数１】
2015111384.tif 000004
【００１７】
ここで、φは同定するパラメータ、u は既知の入力、w は外乱、y1 、y2 は同定に用いる出力である。なお、システムの劣化によるパラメータ変化を想定するので、パラメータφの同定は真値φ*の近傍で考える。また、φの変化はシステムの動特性に比べて遅く、同定においては一定と仮定する。
図１、(1)式より、予測出力 y(φ) は(2)式で記述される。
【００１８】
【数２】
2015111384.tif 000005
【００１９】
ただし、G(z)◇uはシステムG(z)に入力uを印加したときの出力をあらわす。一方、同定実験によって得られた出力を、次式のように表記する。
【００２０】
【数３】
2015111384.tif 000006
【００２１】
ここで、次式のような定義を行い、フィルタ F(z;φ)を定める。
【００２２】
【数４】
2015111384.tif 000007
【００２３】
上記したG11、G21がスカラーの場合、（不安定な極零相殺を許せば）一般性を失うことなくFは、(3)式で与えられる。
【００２４】
【数５】
2015111384.tif 000008
【００２５】
ただし、W(z;φ) は周波数重み関数である。
F(z;φ) を用いて信号vを、(4)式のように定める。
【００２６】
【数６】
2015111384.tif 000009
【００２７】
このとき、(5)式のようになるので、v(φ) はwに依存しない。
【００２８】
【数７】
2015111384.tif 000010
【００２９】
ただし、Gvuには、次式の関係がある。
【００３０】
【数８】
2015111384.tif 000011
【００３１】
ここでさらに、フィルタリングされた出力に関する予測誤差 e(φ)を次式のように定める。
【００３２】
【数９】
2015111384.tif 000012
【００３３】
すると、φ＝φ*のとき、e(φ)=0 が成り立つ。そこで、数値最適化手法を用いて、(6)式で定義される評価関数を最小にするφを求めることでφの同定を行う。
【００３４】
【数１０】
2015111384.tif 000013
【００３５】
以上までが、パラメータ同定手法の概略についての説明である。
［モータ制御系におけるパラメータ同定］
次に、上述したパラメータ同定手法を、実際にモータ制御系に適応した実施形態、言い
換えればモータ制御系におけるパラメータ同定について説明する。
まず、モータ制御系におけるパラメータ同定の手法について述べることとする。
（同定実験の結果）
同定実験において、対象とする制御系はモータに対するＰＤ制御系である（図２を参照）。なお、使用したモータ実験装置は、図３のように電流フィードバックで駆動される装置である。各パラメータのカタログ値を表１に示す。
【００３６】
【表１】
2015111384.tif 000014
【００３７】
図２の制御系においてトルク外乱Tf および目標信号rから、モータ角度θmおよび制御入力である印加電圧Vi  までのシステムは、(7)式〜(9)式によって与えられる。
【００３８】
【数１１】
2015111384.tif 000015
【００３９】
パラメータ同定を行う際には、rを既知入力、θm、Viを同定出力とする。
同定実験は、Real-Time Workshop （Mathworks社のＣ言語のコード発生器、商標）を用いて、図４に示すSimulinkモデル（Mathworks社のモデリング、シミュレーション、及び解析のためのシステム、登録商標）を利用して行った。ＰＤ補償器を事前に離散化することは行わず、Simulinkの固定ステップシミュレーションで実験を行った。サンプル周期は4msecである。
【００４０】
シミュレーションでの検討により、Rsの同定精度が他のパラメータの同定精度に大きな影響を与えることが判明しているので、まずは、Rsの同定を試みた。同定問題は離散時間システムに対して定義されているが、(7)式は連続時間システムであるので、そのままでは同定手法を適用できない。そこで、(7)式をゼロ次ホールド等価な意味で離散化して P(z; Rs) を求めた。離散化は各Rsに対して数値的に行うので、Rsに対して解析的な表現が得られるわけではない。しかしながら、与えられたRsに対する離散時間システムを求めることができれば数値最適化は可能である。
【００４１】
(3)式において W(z;φ)=1とし、得られた同定実験データを用いてパラメータ同定を行
ったところ、収束値が負の値となった。そこで、Rsをカタログ値のα倍に設定し、各α に対して評価関数Γ(αRs)の値を計算した。
図５の「●」のプロットはその結果である。図の横軸は、αの対数軸である。(1)式に基づいて理論的に検討する限りにおいては、評価関数はα=1の近傍で下に凸になるはずである。しかしながら、図５より、評価関数はα=1の近傍で下に凸にはなっていないことがわかる。また、離散化の方法を双一次変換、一次ホールド等価、インパルス不変法などに変更しても結果は改善されなかった（詳しくは、図５を参照）。
（シミュレーションでの原因調査）
上述したように、同定実験の結果は芳しいものとはなっていない。
そこで、シミュレーションに基づいて、その原因の究明を行った。考えられる原因はいくつかあるが、ここでは以下の２点について調べた。
【００４２】
・離散化による影響
・外乱周波数とサンプリング周波数の関係
図６は、外乱に白色外乱を与えて行ったシミュレーションの応答である。応答は閉ループ系を
・すべて連続時間のまま
・閉ループ系をゼロ次ホールド等価に離散化する
・対象は連続時間のままとし、コントローラのみ Tustin変換を用いて離散化する
の３通りで計算した。
【００４３】
ただし、「連続時間のまま」は連続時間の応答を求めるという意味ではなく、Simulink
に実装されている微分方程式の求解アルゴリズムを用いてシミュレーションを行うという意味である。また、白色外乱のサンプリング周期を 4msec、400msec とした。図６より、コントローラのみを離散化した場合には他の応答との間にわずかに誤差が生じていることがわかる。
【００４４】
以上の応答に対して前節のフィルタを適用した際の評価関数値を図７に示す。
図７より、すべて連続時間あるいはすべて離散時間の場合には、外乱の周波数によらず想定される評価関数になっていることが確認できる。一方、コントローラのみを離散化した場合には、真値の近傍で下に凸にならないような評価関数になっている。後者は同定実験で得られた結果と同様である。また、この場合にも外乱の周波数はほとんど影響を与えていない。
【００４５】
さらに、図８はRsが真値の場合に外乱を０としたときの目標値応答を求めたものである。図８の緑の応答がコントローラのみ離散化したシステムの応答であり、図８の赤の応答は連続時間システムの応答である。図８のような応答を用いて求めた評価関数を図９に示す。図９より、この場合にも評価関数値が真値の近傍で下に凸になっていないことが確認できる。
【００４６】
以上のように、パラメータ同定に問題が生じるのはコントローラの離散化によることがわかる。しかしながら、連続時間系の一巡伝達関数のゲイン線図を見ると、サンプリング周波数500π≒1500rad/sec 近傍でのゲインは-70dB 程度であり、十分小さい。よって、一巡伝達関数のゲインが高すぎるわけでも、制御系としてサンプリング周波数が低いわけでもない。
【００４７】
コントローラは近似微分を用いたＰＤ制御を用いているので、高周波ゲインは０にはなっていない。これがエイリアシングと同じような状況を引き起こしている可能性がある。そこで、コントローラを次式のような厳密にプロパーなものに変更した。
【００４８】
【数１２】
2015111384.tif 000016
【００４９】
K2(s) に対応する制御系に対して、図６と同じ条件でシミュレーションを行った結果を図１０に示す。白色外乱のサンプリング周期はいずれも4msである。また、このときの評
価関数の計算結果を図１１に示す。
図１１と図７を比較すると、図１１では真値の近傍で評価関数値がやや減少していることが確認できる。また、評価関数の形状も純連続・純離散の場合のものに近づいている。ただし、評価関数の値が純連続・純離散の場合とは大きく異なっていることから、コントローラが厳密にプロパーか否かは問題の本質的原因ではないと考えられる。すなわち、離散化ではあっても、エイリアシングが根本的な原因とは考えられない。
（理論的解析）
以上のように、コントローラのみ離散化することで、パラメータ同定が失敗に終わることが再現された。しかしながら、同時に、その原因はエイリアシングなどの離散化特有の性質とは考えにくいことも示唆されている。そこで、コントローラのみ離散化した場合と全て連続時間・全て離散時間の場合との違いを改めて考えてみる。
【００５０】
全て離散化したシステムでは連続時間で閉ループ伝達関数を求めた後に離散化を行っている。全て連続のシステムにおいてもシミュレーションを行う段階では本質的に同じ操作を行っている。一方、コントローラのみ離散化したシステムでは制御対象やコントローラを個別に離散化した後に、それらを接続することで閉ループ系を構成している。一般に、連続時間の閉ループ系を離散化したシステムと個々のシステムを離散化した後にそれらを接続して構成された離散時間系は異なるシステムとなる。
【００５１】
そこで、個々のシステムが離散化された状況を想定して、パラメータ同定問題における評価関数を理論的に解析する。
図２のモータ制御系において、コントローラ、制御対象をそれぞれ個別に離散化すると、閉ループ系は図１２のように表現できる。ただし、K(z)はコントローラ、α(z;φ)は制御入力からトルク、P(z;φ) はトルクから出力までの伝達関数である。α、P はいずれもパラメータに依存するものとする。前述のシミュレーションではφ=Rsであり、Pはφに依存しないものであった。ここではより一般的に、Pもφに依存する場合を考える。
このとき、入力u1、u2と出力y1 、y2は(10)式、(11)式の関係を満たすものとする。
【００５２】
【数１３】
2015111384.tif 000017
【００５３】
(10)式、(11)式を伝達関数行列を用いて表現すると、(12)式が得られる。
【００５４】
【数１４】
2015111384.tif 000018
【００５５】
さらに、(12)式の両辺に、(12)式の左辺の一項目の行列の逆行列を左からかけると、(13)式が得られる。
【００５６】
【数１５】
2015111384.tif 000019
【００５７】
ここで、φが真値φ*をとるときの α(z;φ)、P(z;φ) をそれぞれα*(z)、P*(z) とする。すなわち、
【００５８】
【数１６】
2015111384.tif 000020
【００５９】
である。
(13)式より、与えられた外乱u1（図１におけるwに対応） 、外生入力u2（図１におけるuに対応）に対して、実際に測定される出力 y*は(14)式で記述される。
【００６０】
【数１７】
2015111384.tif 000021
【００６１】
(13)式より、フィルターＦは周波数重み関数Ｗを用いて
【００６２】
【数１８】
2015111384.tif 000022
【００６３】
と与えられる。このとき、(13) 式より v = F◇u は、
【００６４】
【数１９】
2015111384.tif 000023
【００６５】
である。一方、(14)式より F◇y*は、
【００６６】
【数２０】
2015111384.tif 000024
【００６７】
となり、真値φ*に依存しない。(16)式、(17)式より、誤差eはφによらず e≡0 となる。すなわち、評価関数値はφによらずゼロである。
以上より、図１２のシステムを考える限り、誤差 e、評価関数はφによらずゼロであるためにφの同定が不可能であることが明らかとなった。これは、K(z) が既知の伝達関数であるために、パラメータφに関して y1とy2は本質的に同じ情報しか持っていないことが原因であると考えられる。このことはさらに、２つ目の出力として、制御入力ではなくy1とは独立にφに依存する出力を得る必要があることも意味している。
【００６８】
図７の全て連続・離散のシミュレーションでは、図１２と同じシステムを用いているにもかかわらず、評価関数は真値の近傍で下に凸になっている。これは既に述べたように、閉ループシステムが離散化と接続の順序に依存するためである。離散化してから接続した場合、y1とy2はφに関して本質的に同じ情報しか持たない。一方、接続してから離散化した場合、対象とコントローラの動特性が渾然一体となったまま離散化されるため、y1とy 2にはφに関して異なる情報が得られる可能性がある。全て連続の場合も、シミュレーションの際に離散化を行うことになるので状況は同じである。
（同定可能な出力）
前述したように、図１２の構成ではφを同定することはできない。φを同定するには y1とは独立にφに依存する出力を新たにとる必要がある。そこで、ここでは図１２のαを α2・α1に分割し、α1の出力をy2にとることを考える。これは、図２のモータ制御系では電流フィードバックに用いる電流センス抵抗の電圧値を測定すること、すなわち、モータの電機子電流を測定することに対応する。
【００６９】
図１３のシステムにおいては、入力u1、u2と出力y1、y2は、(18)式、(19)式の関係を満足する。
【００７０】
【数２１】
2015111384.tif 000025
【００７１】
これより、入出力関係は、(20)式で与えられる。
【００７２】
【数２２】
2015111384.tif 000026
【００７３】
ただし、α=α2・α1である。このときフィルターFは(21)式となる。
【００７４】
【数２３】
2015111384.tif 000027
【００７５】
このとき、v、F◇y*は以下のようになる。
【００７６】
【数２４】
2015111384.tif 000028
【００７７】
ただし、α*1、α*は、それぞれφが真値φ*であるときのα1 、αである。さらに、e は次式で与えられる。
【００７８】
【数２５】
2015111384.tif 000029
【００７９】
(24) 式より、e は（外乱u1が特別な場合を除いて）ゼロにはならないことがわかる。さらに、α がφに関して連続であれば、φがφ*に近づくことで α1がα*1に近づき、e はゼロに近づくことがわかる。よって、図１３の構成をとることができれば φの同定は可能と考えられる。
実際、モータ制御系において電流センス抵抗電圧をy2にとり、シミュレーションに基づいて評価関数を求めたものが図１４である。図１４より、φが真値の近傍で評価関数値が下に凸な準凸関数になっており、評価関数の最小化によってφを同定可能であることが確認できる。
（Fのインナー化）
さて、図１４の評価関数はφが小さくなると急激に大きくなる。これは、φが小さくな
るとFのゲインが大きくなるためである。仮にφ≠φ*のφで Fのゲインが小さくなるようなフィルタになっていると、評価関数はφ=φ*以外のφで下に凸になる可能性もある。評価関数の最小化によって同定を行うことを考えると、φによって Fのゲインが変化することは望ましくない。とくに、F のゲインは周波数に依存するため、F そのものを用いると設計者の意図とは無関係に周波数重み関数を導入していることになる。
【００８０】
よって、F のゲインはすべてのφとすべての周波数で一定であることが望ましい。そこで、F をそのまま用いるのではなく、各φに対して F を F = Fo Fiとインナ・アウタ分解し、Fiを用いて評価関数の計算を行うほうがよい。ただし、Fiはインナ因子、すなわち
Fi(ejθ) Fi(ejθ)*= 1を満足する安定なシステムである。一方、Foはアウタ因子、すなわち |Fo(jθ)|=|F(ejθ)| を満足する安定かつ最小位相なシステムである。
【００８１】
Fのインナアウタ分解はリカッチ方程式を用いて求めることができる。具体的には Fの状態空間実現を
【００８２】
【数２６】
2015111384.tif 000030
【００８３】
とする。このとき、リカッチ方程式
【００８４】
【数２７】
2015111384.tif 000031
【００８５】
の安定化解Xを用いると、Fのインナ因子は、
【００８６】
【数２８】
2015111384.tif 000032
【００８７】
で与えられる。ただし、  
【００８８】
【数２９】
2015111384.tif 000033
【００８９】
であり、Π は ΠΠT = CXCT + DDT を満足する正則行列である。Π は唯一ではない。例えばスカラーの場合には ± 1 倍だけの自由度がある。しかし、これらの自由度はインナアウタ分解の自由度と同じである。よって、Π は ΠΠT = CXCT + DDT を満足する正則行列ならば、どれをとっても構わない。
上記のアルゴリズムでインナアウタ分解を行うためには、F = [F1 F2] のF1、F2の少なくとも一方は双プロパーである必要がある。F1 、F2が厳密にプロパー、すなわち、F自体が厳密にプロパーな場合には、例えば F(z) zｌを改めて F(z) とすればよい。ただし、ｌはF1 、F2それぞれの相対次数のうち、小さいほうの値である。F(z) zｌは原点に ｌ 個の零点を有するが、これらの零点はアウタ因子に含まれるため、Fiは zｌには依存しない。
【００９０】
Fの代わりに Fiを用いて評価関数を計算したものが図１５である。図１４よりも評価関数のカーブが緩やかなものになっている。図１４と図１５の違いはさほど明確ではないが、これはもともと図 １４のグラフはφの真値の近傍で下に凸になっているためである。
（作用効果）
以上述べたように、本発明のシステムのパラメータ同定方法を用いることで、負荷外力の存在するシステムのパラメータを精確に同定することが可能となる。
［外乱の推定］
前述した「パラメータ同定手法の概要」において、φが定まったとき、外乱 wを計算したが、その具体的な計算方法について、以下に述べる。
（推定手法）
まず、真値φ*に対して y*= y(φ*) となるはずなので、(2)式より、
【００９１】
【数３０】
2015111384.tif 000034
【００９２】
が成り立つ。
G11(z;φ*)、G21 (z;φ*) の相対次数の小さいほうを r、両者に共通する不安定零点を
ζ1、・・・、ζm とすると、G11(z;φ*)、G21 (z;φ*)を以下のように分解することができる。
【００９３】
【数３１】
2015111384.tif 000035
【００９４】
であり、c1、c2は安定有理関数環上で互いに素な伝達関数である。c1、c2が互いに素であることから、ベズー方程式
【００９５】
【数３２】
2015111384.tif 000036
【００９６】
を満足する安定な伝達関数 x(z)、y(z) が存在する。このとき、(25)式、(26)式より、
【００９７】
【数３３】
2015111384.tif 000037
【００９８】
が得られる。よって、真値 φ*が得られれば、q = n(z)◇wは計算可能である。
qからwを求める際、n(z) の相対次数と不安定零点が障害となる。このうち、相対次数は時間軸がシフトするだけなので大きな問題とはならない。一方、不安定零点については、単純にqに1/n(z)を作用させるとn(z)の不安定零点が不安定極となるために、得られる
信号は時間の増大に対して発散するものとなる。
【００９９】
q=n(z)◇w を満足するwは、一般に
【０１００】
【数３４】
2015111384.tif 000038
【０１０１】
と書かれる。ただし、w0 は q = n(z)◇w を満足する一つの信号、ψiは n(z) のゼロダイナクスの基底となる信号である。n(z) の零点はすべて不安定零点であることから、ψiはいずれも（また、一般には w0も）時間の増大に対して発散する信号となる。(28)式からわかるようにq = n(z)◇w を満足する w は唯一には定まらない。しかしながら、(28)式を踏まえれば、例えば、｜w｜ が最小となるような w を求めることは可能である。もちろん、G11 (z;φ*)、G21 (z;φ*) が不安定零点を共有しなければ、このような不都合は生じない。
（モータ制御系を用いた有効性検証）
前節の推定手法をモータ制御系に適用し、その有効性を検証する。出力はモータ角度と電機子電流である。また、φの真値が得られたものと仮定する。このとき、G11 、G21は、
【０１０２】
【数３５】
2015111384.tif 000039
【０１０３】
である。ただし、z 1 = -1 + 4.29 × 10-6である。これより、G11 、G21 に共通する不安定零点は存在しないことがわかる。G11 、G21 双方とも相対次数は 1 なので、c 1 、c
2を以下のように定める。
【０１０４】
【数３６】
2015111384.tif 000040
【０１０５】
このとき、c1 、c2 とも双プロパーである。上記の c1 、c2 に対して、(27) 式のベズー方程式を満足する x、y を求める。ベズー方程式の解はリカッチ方程式を用いて求めることができる。ただし、今回は、上記の式の右辺の項の左逆システムが安定であったので、それを用いることにした。具体的には G11 、G21 の状態空間実現
【０１０６】
【数３７】
2015111384.tif 000041
【０１０７】
を用いて、x、y を
【０１０８】
【数３８】
2015111384.tif 000042
【０１０９】
で定めた。Ｄ†は Ｄの疑似逆行列である。この x、y に対して次式が成り立つ。
【０１１０】
【数３９】
2015111384.tif 000043
【０１１１】
求めたx、y を用いて
【０１１２】
【数４０】
2015111384.tif 000044
【０１１３】
を求めれば q =(1/z)◇w  となるので、1 サンプルステップだけ時間を遅らせた外乱信号が求まるはずである。
図１６は外乱を与えてシミュレーションを行った際の応答から q を求めた結果である。図１７、図１８はそれぞれq を１ステップ、２ステップシフトした結果である。図１９は１ステップ、２ステップシフトした場合それぞれについて、印加されている外乱との誤差を求めたものである。q を２ステップシフトしたものは実際に印加されている外乱によく一致している。この結果から、提案手法の有効性が確認できる。
【０１１４】
一方、理論的には q は２ステップではなく１ステップシフトすることで w に一致するはずである。しかしながら、上記の結果はそれとは異なっている。現時点では、この原因は不明であり、今後、その原因をあきらかにする必要がある。
（可同定性の十分条件）
なお、本実施形態では、図１２のシステムでは可同定性を満足せず、図１３のシステムにすることで可同定性が満足することを示した。より理論的な式で十分条件を示すと、与えられたG(z;φ)に対し、次式が成り立てばφ*は可同定である。
【０１１５】
【数４１】
2015111384.tif 000045
【０１１６】
実際、図１３のシステムではα1(1)≠0であれば、Hr(z)≠0となりφ*は可同定となる。なお、＊(1)は＊のφによる微分を示しており、α1(1)≠0は一般的な図１３のシステムでは成り立っている。(29)式は閉ループ伝達関数Gの外乱の影響分（(29)式のdet内の行列の第１列）とパラメータφの変化の影響（(29)式のdet内の行列の第２列）が独立であることを示している。また、独立であるということは、外乱の影響分とパラメータφの変化の影響分とが分離可能であることを示している。
（パラメータが同定可能となるための必要十分条件）
次に、図１に示すシステムでの可同定性、言い換えれば、図１に示すシステムにおいて、パラメータが同定可能となる(29)式が成り立つための必要十分条件を考える。
【０１１７】
図１において、G(z;φ) は閉ループ系の伝達関数行列となる。図２０は、この伝達関数行列G(z;φ)をより詳細に記したものであり、P(z;φ)，K(z) はそれぞれ制御対象、補償器（コントローラ）である。また w、u、v はそれぞれ外乱、目標信号、制御入力であり、y1、y2は観測出力である。
図２０におけるy1、y2を
【０１１８】
【数４２】
2015111384.tif 000046
【０１１９】
とし、各信号に対応して P(z;φ)，K(z) を下記のように分割しておく。
【０１２０】
【数４３】
2015111384.tif 000047
【０１２１】
このとき、
【０１２２】
【数４４】
2015111384.tif 000048
【０１２３】
である。それ故、
【０１２４】
【数４５】
2015111384.tif 000049
【０１２５】
が成り立つので、
【０１２６】
【数４６】
2015111384.tif 000050
【０１２７】
となる。
よって、図２０に対応するG(z;φ) は次式で与えられる。
【０１２８】
【数４７】
2015111384.tif 000051
【０１２９】
この G(z;φ) が、(29)式を満たすための必要十分条件は、Ku(z)≠0 のとき、
【０１３０】
【数４８】
2015111384.tif 000052
【０１３１】
が成り立つことである。
Ku(z)≠0 は目標信号uが制御入力vに反映されることを示しており、通常のフィードバック制御系では成り立っている。したがって、補償器（コントローラ）K(z)には関係しない制御対象P(z;φ)が(30)式を満たせば、ほとんど全ての任意の外乱に対してパラメータφが可同定であるための十分条件となる。
【０１３２】
(30)式は制御対象の伝達関数Pの外乱の影響分（(30)式のdet内の行列の第１列Pw(z;φ)）とパラメータφの変化の影響（(30)式のdet内の行列の第２列Pv(z;φ) (1) ）が独立であることを示している。また独立であるということは、外乱の影響分と、パラメータφの変化の影響分とが分離可能であることを示している。
以上、今回開示された実施形態はすべての点で例示であって制限的なものではないと考えられるべきである。本発明の範囲は上記した説明ではなくて特許請求の範囲によって示され、特許請求の範囲と均等の意味及び範囲内でのすべての変更が含まれることが意図される。
【図１】
2015111384.tif 000054
【図２】
2015111384.tif 000055
【図３】
2015111384.tif 000056
【図４】
2015111384.tif 000057
【図１２】
2015111384.tif 000058
【図１３】
2015111384.tif 000059
【図５】
2015111384.tif 000060
【図６】
2015111384.tif 000061
【図７】
2015111384.tif 000062
【図８】
2015111384.tif 000063
【図９】
2015111384.tif 000064
【図１０】
2015111384.tif 000065
【図１１】
2015111384.tif 000066
【図１４】
2015111384.tif 000067
【図１５】
2015111384.tif 000068
【図１６】
2015111384.tif 000069
【図１７】
2015111384.tif 000070
【図１８】
2015111384.tif 000071
【図１９】
2015111384.tif 000072
【図２０】
2015111384.tif 000073