JP2021043573A - 出力推定方法 - Google Patents

Abstract

【課題】一回の実験データのみに基づいて、ある制御器を導入した時の閉ループ系の時間応答を推定する。【解決手段】制御器Ｋを介して制御対象への入力に対する出力を推定する方法であって、１回の試行に基づいて制御対象への周期入出力応答（周期入力信号ｕ０と周期出力信号ｙ０）を求める工程と、仮想入力信号ｒに対する前記制御対象の応答（仮想出力信号ｙ）を（１）式で求める工程を有する。ここで、ｒ１（ｊω）は、（２）式に従う。ｒ１（ｊω）、ｒ（ｊω）、ｕ０（ｊω）、ｙ０（ｊω）は、それぞれ信号ｒ１、ｒ、ｕ０、ｙ０をフーリエ変換したもの、またＫ−１（ｊω）は１／Ｋ（ｊω）、Ｋ（ｊω）は制御器Ｋの伝達関数を表す。【選択図】図１

Description

本発明は、制御器Ｋを介して制御対象への入力に対する出力を推定する方法に関するものであり、特に制御対象への入出力を１度行った結果（以後「一回入出力応答」と呼ぶ。）だけを利用し、制御器をつないだ制御対象への任意の入力に対する出力を推定する方法に関するものである。

通常制御工学では、制御対象の正確なモデルシステムを同定し、そのモデルに基づいて良好な制御性能が期待できる制御器を設計し、その制御性能を実験で検証するアプローチがなされる。

しかし、実際の制御対象の多くは、激しく変化させる応答や、開ループ応答を取得できなかったり、長時間に渡る実験ができなかったりする。この場合、システム同定によって正確なモデルを得ることは難しい。

このようなケースに対して、モデルを求めることを経ず、実験データから直接制御器を調整するデータ駆動制御が近年活発に検討されている（非特許文献１〜４）。そのなかでも、Ｖｉｒｔｕａｌ Ｒｅｆｅｒｅｎｃｅ Ｆｅｅｄｂａｃｋ Ｔｕｎｉｎｇ（ＶＲＦＴ）と、Ｆｉｃｔｉｔｉｏｕｓ Ｒｅｆｅｒｅｎｃｅ Ｉｔｅｒａｔｉｖｅ Ｔｕｎｉｎｇ（ＦＲＩＴ）は、一回の実験データのみに基づいて、閉ループ系が参照モデルに近づくように制御ゲインを最適化できる。そのため、多くの応用例が報告されている。

Ｍ．Ｃ．Ｃａｍｐｉ，Ａ．Ｌｅｃｃｈｉｎｉ，Ｓ．Ｍ．Ｓａｖａｒｅｓｉ，Ｖｉｒｔｕａｌ ｒｅｆｅｒｅｎｃｅ ｆｅｅｄｂａｃｋ ｔｕｎｉｎｇ：ａ ｄｉｒｅｃｔ ｍｅｔｈｏｄ ｆｏｒ ｔｈｅ ｄｅｓｉｇｎ ｏｆ ｆｅｅｄｂａｃｋ ｃｏｎｔｒｏｌｌｅｒｓ，Ａｕｔｏｍａｔｉｃａ，Ｖｏｌ．３８，ｐｐ．１３３７−１３４６，２００２． 相馬将太郎，金子修，藤井隆雄，一回の実験データを用いた制御器パラメータチューニングの新しいアプローチ− Ｆｉｃｔｉｔｉｏｕｓ Ｒｅｆｅｒｅｎｃｅ Ｉｔｅｒａｔｉｖｅ Ｔｕｎｉｎｇ の提案，システム制御情報学会論文誌，Ｖｏｌ．１７，Ｎｏ．１２，ｐｐ．５２８−５３６，２００４ 金子修，山本透他，特集データ駆動制御：新機軸と新地平，計測と制御，Ｖｏｌ．５２，Ｎｏ．１０，ｐｐ．８９２−８９７，２０１３ 田坂謙一，加納学，小河守正，増田士朗，山本透，閉ループデータに基づく直接的ＰＩＤ調整とその不安定プロセスへの適用，システム制御情報学会論文誌，Ｖｏｌ．２２，Ｎｏ．４，ｐｐ．１３７−１４４，２００９

ＶＲＦＴやＦＲＩＴでは、参照モデルが適切でない場合は、制御性能が悪化するばかりか、制御系が不安定になるおそれすらあるという問題があった。

本発明は上記課題に鑑みて想到されたものであり、一回の実験データのみに基づいて、ある制御器を導入したときの閉ループ系の時間応答を推定する。

より具体的に、本発明に係る出力推定方法は、
制御器Ｋを介して制御対象への入力に対する出力を推定する方法であって、
１回の試行に基づいて制御対象への周期入出力応答（周期入力信号ｕ_０と周期出力信号ｙ_０）を求める工程と、
仮想入力信号ｒに対する前記制御対象の応答（仮想出力信号ｙ）を（１）式で求める工程を有することを特徴とする。

ここで、ｒ_１（ｊω）は、（２）式に従う。ｒ_１（ｊω）、ｒ（ｊω）、ｕ_０（ｊω）、ｙ_０（ｊω）は、それぞれ信号ｒ_１、ｒ、ｕ_０、ｙ_０をフーリエ変換したものを表す。またＫ^−１（ｊω）は１／Ｋ（ｊω）を表し、Ｋ（ｊω）は制御器Ｋの伝達関数を表す。Ｆ^−１は逆フーリエ変換を表す。

本発明に係る出力推定方法では、実験データのフーリエ変換と制御器の周波数特性からフィードバック系の入出力応答の周波数成分を求め、それを逆フーリエ変換して仮想時間応答を求める点に特徴がある。つまり、仮想出力信号（推定された出力信号）を求めるのに、伝達関数モデルを一切用いていないので、モデルの適不適や、予めモデルについての構造などの不知などに基づく課題を回避することができる。

本発明に係る出力推定方法の処理の流れを示す図である。 非周期入力信号ｕ_００と非周期出力信号ｙ_００を周期入力信号ｕ_０と周期出力信号ｙ_０へ変換する手順を示す図である。 図３（ａ）は、非周期入力信号ｕ_００と非周期出力信号ｙ_００の関係を示すブロック図であり、図３（ｂ）は、仮想入力信号ｒと仮想出力信号ｙの関係を示すブロック図であり、図３（ｃ）は制御系にｒ１を入力した場合の信号の流れを示すブロック図である。 Ｇ１からＧ９までの制御対象に対し、非周期入力信号としてステップ信号を入力したときの非周期出力信号の観測結果を示す図である。 制御対象への真の時間応答と仮想時間応答を比較したグラフである。 各制御器のパラメータを最適化した場合の真の時間応答と仮想時間応答を比較したグラフである。 一回入出力応答ｕ_００、ｙ_００が閉ループデータで、Ｇが不安定な場合のフィードバック系の真の出力応答と仮想時間応答を示す図である。 パラメータを最適化した場合の真の時間応答と仮想時間応答を比較したグラフである。 実施例に用いたモータの写真である。 モータに対する一回入出力応答と、真の時間応答（真の出力信号）と、仮想時間応答を比較したグラフである。 モータに対する一回入出力応答の真の時間応答（真の入力信号）を示した図である。

以下に本発明に係る出力推定方法について図面および実施例を示し説明を行う。なお、以下の説明は、本発明の一実施形態および一実施例を例示するものであり、本発明が以下の説明に限定されるものではない。以下の説明は本発明の趣旨を逸脱しない範囲で改変することができる。

また、非周期入力信号ｕ_００、非周期出力信号ｙ_００、周期入力信号ｕ_０、周期出力信号ｙ_０、制御系への仮想入力信号ｒ、制御系への仮想出力信号ｙ、制御対象への仮想入力信号ｕは、時間関数で連続信号である。しかし、公知の方法によっていつでも離散的処理として扱うことができる。すなわち、本明細書において、連続関数としての表現は離散的な表現を含むものである。以後離散的に信号を扱う際には、ｒ（ｋ）等のようにサンプル数を表す番号ｋが括弧付で用いられる場合もある。

図１に本発明に係る出力推定方法の処理の流れを示し、その後各工程について詳細を示す。本発明に係る出力推定方法では、制御対象へ周期入力信号ｕ_０を入力した時の周期出力信号ｙ_０を一回入出力応答（一回の試行）で求め、その周期入出力応答を使って仮想入力信号ｒに対する仮想出力信号ｙを推定する。ここで周期入出力応答を求める方法が２通りある。

図１の処理Ｐ１００を参照する。Ｇ（ｚ）は、加熱炉やモータといった制御対象を示す。ｚは進み演算子である。第１の方法は、まず周期入力信号ｕ_０を作製し、それを制御対象Ｇへ入力することで周期出力信号ｙ_０を求める。これは直接的な方法である。実際に入出力を測定するのはこの１回だけである。

つぎに処理Ｐ１０１および処理Ｐ１０２を参照する。第２の方法は、非周期入力信号ｕ_００を制御対象Ｇへ入力し非周期出力信号ｙ_００を得る（処理Ｐ１０１）。そして次に処理Ｐ１０１で得た非周期入力信号ｕ_００および非周期出力信号ｙ_００を周期応答信号へと変換する（処理Ｐ１０２）。周期応答信号への変換は、それぞれの信号を反転させた信号を元の信号に続けて加えることで得ることができる。

変換された周期入力信号ｕ_０および周期出力信号ｙ_０は、周期応答信号であり、一定の条件の下で、非周期入力信号ｕ_００および非周期出力信号ｙ_００とみなしてよい。この方法では、実際に制御対象Ｇで信号の入出力を確認するのは、非周期入力信号ｕ_００を制御対象Ｇへ入力し非周期出力信号ｙ_００を得る処理Ｐ１０１の１回だけである。なお、本明細書では、処理Ｐ１０１および処理Ｐ１０２の方法で周期入出力応答（周期入力信号ｕ_０、周期出力信号ｙ_０）を求める方法を主に説明を続ける。

次に処理Ｐ１０３を参照する。ここで制御対象Ｇ（ｚ）を制御する制御器Ｋを導入する。そして、これらの離散時間周波数伝達関数をそれぞれＧ（ｊω）およびＫ（ｊω）とする。また、制御器Ｋを含む制御系全体に対する仮想入力信号ｒをフーリエ変換したものをｒ（ｊω）とする。仮想入力信号ｒ（ｊω）による制御器Ｋ（ｊω）の出力を、制御対象Ｇ（ｊω）に対する仮想入力信号ｕ（ｊω）とし、制御対象Ｇ（ｊω）に仮想入力信号ｕ（ｊω）が入力されたときの仮想出力信号をｙ（ｊω）とする。すると制御器Ｋ（ｚ）に制御された制御対象Ｇ（ｚ）の仮想出力信号ｙは、（１）式で表される。

ここで、ｒ_１（ｊω）、ｒ（ｊω）、ｕ_０（ｊω）、ｙ_０（ｊω）は、それぞれ信号ｒ_１、ｒ、ｕ_０、ｙ_０をフーリエ変換したものを表す。またＫ^−１（ｊω）は１／Ｋ（ｊω）を表し、Ｋ（ｊω）は制御器Ｋの伝達関数を表す。Ｆ^−１は逆フーリエ変換を表す。

また、ｒ_１（ｊω）は、（２）式のように表される。ｕ_０（ｊω）、ｙ_０（ｊω）は、周期入力信号ｕ_０および周期出力信号ｙ_０をフーリエ変換したものである。仮想入力信号ｒは、もちろん既知の信号である。ｒ_１（ｊω）は、自ら設定した制御器Ｋの伝達関数と、周期入力信号ｕ_０（ｊω）および周期出力信号ｙ_０（ｊω）という既知の情報で表される。

したがって、（１）式は、制御系へ仮想入力信号ｒが入力された時の仮想出力信号ｙを既知の情報だけから推定できることを示している。このようにして求めた制御対象への仮想入力信号ｕと仮想出力信号ｙを仮想時間応答と呼ぶ。

なお、仮想入力信号ｒ（ｋ）は、評価したい閉ループ系（制御系）への目標入力である。周期入力信号ｕ_０、周期出力信号ｙ_０、仮想入力信号ｒ（ｋ）の離散フーリエ変換をそれぞれ周期入力信号ｕ_０（ｊω）、周期出力信号ｙ_０（ｊω）、仮想入力信号ｒ（ｊω）とする。ここでｊは虚数単位、ω［ｒａｄ／ｓ］はｔ_ｓ［ｓ］をサンプル時間として（３）式で与えられる角周波数である。

また、本発明の出力推定方法では周期入力信号ｕ_０、仮想入力信号ｒおよび離散時間制御器Ｋについて、次の仮定を満たすものとする。
Ａ１）ｕ_０、ｒは、周期的信号の整数周期分の時系列である。このとき、ｙ_０も周期的信号となる。
Ａ２）ｕ_０は、ｒに含まれる周波数成分と同じ周波数の成分を含む。つまり、あるω_ａでｒ（ｊω_ａ）≠０ならば、ｕ_０（ｊω_ａ）≠０である。その対偶をとると、あるω_ａでｕ_０（ｊω_ａ）＝０ならばｒ（ｊω_ａ）＝０である。
Ａ３）Ｋはフィードバック系を安定化する。
Ａ４）Ｋ^−１（ｊω）ｕ_０（ｊω）の逆離散フーリエ変換の最終値を定常状態とみなすことができる。これはＫ（ｚ）を逆Ｚ変換または逆双一次変換したＫ（ｓ）の零点が虚軸上になく、データ数Ｎが十分大きければ満たされる。
Ａ５）Ｇ（ｚ）を逆Ｚ変換または逆双一次変換したＧ（ｓ）の零点が虚軸上にあるとき、Ｋ（ｓ）はその零点を極としてもたない。
Ａ６）非周期入力信号ｕ_００、非周期出力信号ｙ_００をステップ応答データとした場合、非周期入力信号ｕ_００、非周期出力信号ｙ_００の初期値と最終値は定常状態とみなせる。

＜処理Ｐ１０１＞
次に処理Ｐ１０１から詳細に説明を行う。処理Ｐ１０１では、制御対象Ｇ（ｚ）に実際の非周期入力信号ｕ_００を入力し、その出力である非周期出力信号ｙ_００を得る。

制御対象Ｇ（ｚ）は、線形時不変な離散時間一入出力系で、一回の実験データはＧ（ｚ）の入出力の時系列である非周期入力信号ｕ_００（ｋ）、非周期出力信号ｙ_００（ｋ）であるとする。ｚはすでに述べたように進み演算子であり、ｋ（＝１、２、・・・・、Ｎ）はサンプル数、Ｎはデータの総数である。非周期入力信号ｕ_００、非周期出力信号ｙ_００は閉ループ応答でも開ループ応答でもどちらでもよい。

＜処理Ｐ１０２＞
仮定Ａ６のもと、非周期入力信号ｕ_００、非周期出力信号ｙ_００を用いて仮定Ａ１を満足する周期的信号ｕ_０、ｙ_０を以下のように作成する。

開ループや閉ループの非周期入力信号ｕ_００、非周期出力信号ｙ_００は、ステップ応答とした場合、仮定Ａ１を満たさないが、非周期入力信号ｕ_００、非周期出力信号ｙ_００の波形をコピーし、その上下をひっくり返してつなぎ合わせた波形は仮定Ａ１を満たす。具体的に、図２（ａ）は、非周期入出力応答ｕ_００、ｙ_００をステップ応答波形とした場合を示す。横軸は時間（離散化されているので、ステップ数）であり、縦軸は信号の強さである。また、図２（ｂ）は図２（ａ）の波形をコピーして反転し、つなぎ合わせた波形を示す。図２（ａ）同様に横軸はステップ数であり、縦軸は信号の強さである。

図２（ｂ）の信号は、図２（ａ）のステップ応答に（４）式および（５）式の処理を行うことで、得ることができる。図２（ｂ）の信号は１周期分の周期入出力応答（周期入力信号ｕ_０と周期出力信号ｙ_０）であり、これらは仮定Ａ１を満たす。そして、仮定Ａ６のもとでこれらの周期入力信号ｕ_０、周期出力信号ｙ_０は制御対象Ｇの真の応答と一致する。

以上のような手順に従うことで、一回入出力応答（非周期入力信号ｕ_００、非周期出力信号ｙ_００）を周期入出力応答（周期入力信号ｕ_０、周期出力信号ｙ_０）へ変換することができる。また、所定の仮定の下で、周期信号は、制御対象Ｇの真の応答と一致する。

＜処理Ｐ１０３＞
次に制御対象Ｇを制御器Ｋで制御する際に、制御器Ｋへ仮想入力信号ｒ（ｋ）を入力した時の仮想出力信号ｙ（ｋ）が（１）式で示せることを説明する。

周期信号は正弦波の和で表される。したがって、仮定Ａ１より周期入力信号ｕ_０、周期出力信号ｙ_０、仮想入力信号ｒは正弦波の和で表わすことができる。この時周波数応答は定常特性であり、過渡特性を含まない。制御対象Ｇの離散時間周波数伝達関数をＧ（ｊω）とすると、周期入力信号ｕ_０に対する周期出力信号ｙ_０は（６）式のように表される。

（６）式をブロック線図で表すと図３（ａ）となる。次に制御対象Ｇを制御する制御器Ｋを導入する。制御器Ｋを導入した時のブロック線図を図３（ｂ）に示す。図３（ｂ）より（７）式を得る。

（７）式に（６）式を代入し、（８）式を得る。

ここで（８）式の分母の部分をｒ_１（ｊω）と置けば（９）式を得る。なお、ｒ_１（ｊω）を顕に書けば（１０）式である。また、（１０）式中Ｋ^−１（ｊω）は１／Ｋ（ｊω）の意味である。

（９）式の右辺は、制御器Ｋへの仮想入力信号ｒ（ｊω）と周期入力信号ｕ_０および周期出力信号ｙ_０と、制御器Ｋの伝達関数であって、全て既知の情報である。すなわち、これら既知の情報に基づいて、制御器Ｋへの仮想入力信号ｒ（ｊω）に対する仮想出力信号ｙ（ｊω）が求まる。なお、（２）式は（１０）式と同じである。

（６）式および（１０）式をブロック線図で表すと、図３（ｃ）のように表される。同図より、ｒ_１は制御対象Ｇの入出力応答が周期入力信号ｕ_０、周期出力信号ｙ_０になるときのフィードバック系の目標入力である。仮定Ａ４の下で、ｒ_１は周期的信号とみなすことができる。

また、この時、制御対象Ｇへの仮想入力信号ｕ（ｊω）は、図３（ｂ）と（６）式および（９）式より（１１）式のように得られる。

仮想出力信号ｙ（ｋ）は、仮想入力信号ｒ（ｋ）を制御器Ｋに入力した時の出力として求められたが、これは実際に観測されたものではない。しかし、一回入出力応答である非周期入力信号ｕ_００、非周期出力信号ｙ_００から上記のように推定することができる。仮定Ａ１、Ａ３、Ａ４が満たされる条件下で、仮想入出力応答（制御対象Ｇへの仮想入力信号ｕ（ｊω）、仮想出力信号ｙ（ｊω））を逆離散フーリエ変換することで、仮想時間応答（制御対象Ｇへの仮想入力信号ｕ、仮想出力信号ｙ）を求めることができる。また、ｕは（１１）式を用いる代わりに、ｙを用いてｕ＝（Ｋｒ−ｙ）をシミュレーションしても求めることができる。

一方、この仮想時間応答では、制御器Ｋの下、仮想入力信号ｒ（ｋ）を固定すれば、仮想出力信号ｙ（ｋ）を推定することができるともいえる。したがって、制御器Ｋの構成を決めれば、仮想出力信号ｙ（ｋ）に一定の拘束条件を与えることで、制御器Ｋのパラメータを最適化することもできる。

例えば、制御器ＫをＰＩＤ制御器とする。この制御器Ｋの伝達関数は（１２）式のように表されれる。

この制御器Ｋのパラメータをθ＝（Ｋ_ｐ、Ｋ_ｉ、Ｋ_ｄ）と設定し、制御器Ｋ（θ）と表す。以下の手順で制御器Ｋ（θ）を最適化する。
１）系を安定化する制御器Ｋ（θ）の初期値、オーバーシュートの許容値を設定する。
２）オーバーシュートを拘束条件、整定時間を評価関数として、オーバーシュートが許容値以下で、評価関数を最小にするＫ（θ）を探索する最適化問題を解く。このとき、オーバーシュートと整定時間はＫ（θ）を離散化した制御器を用いた時のフィードバック系の仮想時間応答から求める。また、仮定Ａ４が満たされるように、Ｋ^−１（θ）が安定で遅い極を持たないことも拘束条件に入れる。

なお、本方法はＫ（θ）が系を不安定化する場合には適用できず、仮想時間応答は真の応答との間に誤差が発生してしまうという問題がある。しかし、フィードバック機構が効いている場合、安定限界に近づくと発振条件に近づいて時間応答が振動的になる。

本手順では、制御器の初期値が系を安定化していて、オーバーシュートを拘束条件に入れているため、オーバーシュートが小さくなる方向に最適化すれば振動が抑制されて安定限界から遠ざかると考えられる。

そのため、最適化手法は大域的最適解を求めるタイプではなく、初期値周りの局所的最適解を探索するタイプを用いる必要がある。

このような最適化問題の解法としては、例えば、ＭＡＴＬＡＢのｆｍｉｎｃｏｎを用いることができる。これはｉｎｔｅｒｉｏｒ−ｐｏｉｎｔアルゴリズムという制約付き非線形多変数関数の最小値を探索する手法であり、初期値周りの局所的最適解を求めることができる。もちろん、他の手法を用いてもよい。

＜一回入出力応答ｕ_００、ｙ_００が開ループデータの場合＞
（真の時間応答と仮想時間応答の比較）
本発明による出力推定方法でＰＩＤ制御時の仮想時間応答を計算し、その精度を確認した。以下の９つの連続時間伝達関数（（１３）式〜（２０）式）を制御対象とした。

これらの制御対象は限界感度法によって、同一のＰＩＤゲインを得ることができる。具体的な制御器Ｋの伝達関数は（２１）式のようになる。

Ｇ_１（ｓ）からＧ_９（ｓ）をサンプル時間ｔ_ｓ＝０．０１ｓｅｃでＺ変換し、それぞれに非周期入力信号ｕ_００をステップ信号として１００秒間印加したときの非周期出力信号ｙ_００を得た。

ただし、Ｇ_４については、無定位系でステップ応答の最終値が定常状態にならない。そこでＧ_４のみ非周期入力信号ｕ_００を１０秒までは１、それ以降は０のパルス入力とした。得られた非周期出力信号ｙ_００を図４に示す。図４を参照して、横軸は経過時間であり、縦軸は非周期出力信号ｙ_００の信号強度である。

Ｇ５は、他の制御対象と比較して収束が遅く、１００秒寸前でまだ変化しているが、１００秒時点を定常状態とみなした。

（２１）式に以下の（２２）式を代入して離散化した制御器Ｋを用いたフィードバック系を得た。制御器Ｋの目標値（仮想入力信号ｒ）としてステップ関数を入力した時の真の時間応答と、本発明に係る出力推定方法で得られる仮想時間応答（仮想出力信号ｙ、制御対象への仮想入力信号ｕ）および、その誤差をそれぞれ図５に示す。

図５を参照して、図５（ａ）は、制御対象の出力ｙであり、図５（ｂ）は、制御対象への入力信号（制御器Ｋの出力信号）ｕである。両図とも横軸は時間（ｓｅｃ）であり縦軸は信号強度（無単位）である。図５（ａ）および図５（ｂ）より、真の時間応答と仮想時間応答の値はほぼ重なっていて、真の時間応答と、仮想時間応答の違いが判らない。両者の誤差は、０．０１未満と小さかった。

（制御器Ｋの最適化）
次にオーバーシュート２０％以下を拘束条件として、定常値の２％以内への整定時間を評価関数とし、ＭＡＴＬＡＢのｆｍｉｎｃｏｎによって評価関数を最小にするパラメータの最適化問題を解いた。パラメータは（１２）式のθ＝（Ｋ_ｐ、Ｋ_ｉ、Ｋ_ｄ）である。

初期値は（２１）式の各項の係数とした。また、仮定Ａ４が満たされるように、Ｋ（θ）の零点が安定でその実部の大きさが１［ｒａｄ／ｓ］以上になることを拘束条件に加えた。

最適化して得られたＫ（θ）を用いた時のフィードバック系の真の時間応答と、仮想時間応答を図６に示す。図６（ａ）は、制御対象の出力ｙであり、図６（ｂ）は、制御対象への入力信号（制御器Ｋの出力信号）ｕである。両図とも横軸は時間（ｓｅｃ）であり縦軸は信号強度（無単位）である。図６（ａ）および図６（ｂ）より、真の時間応答と仮想時間応答の値はほぼ重なっており、両者の違いは判らないことがわかる。誤差は、出力ｙは０．０１未満であり、制御対象への入力ｕは０．０３未満と小さかった。

最適化前の図５と比較すると、Ｇ_９だけは初期値に近いところで極値に陥ってあまり変わらなかったが、Ｇ_９以外の全てにおいてオーバーシュートが抑制され、整定時間が小さくなった。

＜一回入出力応答ｕ_００、ｙ_００が閉ループデータでＧが不安定な非最小位相系の場合＞
（真の時間応答と仮想時間応答の比較）
不安定な非最小位相系としては、（２３）式の伝達関数（Ｇ_１０）とした。

状態フィードバックと同一次元オブザーバーを併合した連続時間制御器Ｋ（ｓ）を極配置法で設定し、Ｇ（ｓ）とＫ（ｓ）をサンプル時間０．０１ｓｅｃでＺ変換した。一回入出力応答は、閉ループ系の配置極を、−１０、−１０、−２０、−２０にして目標値をステップ関数にしたときの１００秒間の応答とした。閉ループ系の配置極を、ｐ、ｐ、−２０、−２０にして、ｐを−１、−５、−１０、−２５、−５０とした場合のそれぞれのフィードバック系の真の出力応答と仮想時間応答を図７に示す。

図７は横軸が時間（ｓｅｃ）であり、縦軸は信号強度（無単位）である。図７より真の出力応答と仮想時間応答はほぼ重なっており、その誤差は１０^−９未満と小さかった。

（制御器Ｋの最適化）
次にオーバーシュート２０％以下を拘束条件として、定常値の１０％以内への整定時間を評価関数とし、ＭＡＴＬＡＢのｆｍｉｎｃｏｎによって評価関数を最小にする制御器Ｋ（θ）を探す最適化問題を解いた。制御器Ｋ（θ）は、（２４）式とした。パラメータはθ＝（θ_１、θ_２、θ_３、θ_４）で与えた。

仮定Ａ４が満たされるように、Ｋ（θ）の零点が安定でその実部の大きさが１［ｒａｄ／ｓ］以上になることを拘束条件に加えた。

一回入出力応答（非周期入力信号ｕ_００、非周期出力信号ｙ_００）が得られた時のＫ（θ）を初期値とし、最適化して得られたＫ（θ）を用いた時のフィードバック系の真の時間応答と、仮想時間応答を図８に示す。図８は横軸が時間（ｓｅｃ）であり、縦軸は信号強度（無単位）である。図８よりｙ_００と比べてｙはオーバーシュートが抑制され、整定時間が小さくなった。また、真の出力応答と仮想時間応答はほぼ重なっており、その誤差は１０^−１０未満と小さかった。

＜制御対象をモータにした場合＞
制御対象を図９に示すＤＣモータに直径６８ｍｍ、重量１００ｇの車輪を付けたシステムとした。このＤＣモータ（Ｂｒｉｎｇｓｍａｒｔ ＪＧＡ２５−３７１−２０１）には減速比１：２１．３のギアヘッドとモータ一回転あたり１２カウントのエンコーダが取り付けられている。

エンコーダの出力パルス信号は、ＦＶコンバータ（ＪＲＣＮＪＭ４１５１）によって周波数−電圧変換され、マイコンＡｒｄｕｉｎｏ Ｍｅｇａ２５６０の１０ｂｉｔのＡＤ変換器で量子化される。サンプル時間は０．０１秒であり、制御器で計算した制御入力ｕは、８ｂｉｔで量子化される。

量子化された制御入力ｕは、チョッピング周波数９８０ＨｚＰＷＭでモータドライバ（ＤＦＲＯＢＯＴ Ｌ２９８Ｐ）に出力される。このモータには摩擦による不感帯があり、小さなモータ入力電圧ｖ［Ｖ］では回転しない。これを補償するために、プラス側とマイナス側の不感帯幅ｈ_１、ｈ_２を測定し、（２５）式のように制御入力ｕにｈ_１＝−ｈ_２＝２．０Ｖを加えた。

この式は、不感帯の逆関数の近似であるため、この補償によって制御入力ｕから制御出力ｙまでの特性はほぼ線形になると考えられる。ここでは、目標速度ｒを１１．６７ｒｐｍから２０ｒｐｍにした時のステップ応答を評価した。

ステップ応答法で設計したＰＩＤ制御器により、一回入出力応答データである非周期入力信号ｕ_００、非周期出力信号ｙ_００を１０秒間取得した。その初めの２秒間分の波形を図１０に示す。図１０（ａ）は出力信号を表し、図１０（ｂ）は、モータ（制御対象）への入力を表す。両図とも横軸は時間（ｓｅｃ）であり、縦軸は出力である。図１０（ａ）では回転数（ｒｐｍ）であり、図１０（ｂ）は電圧（Ｖ）である。非周期出力信号ｙ_００は、図１０（ａ）中で「ｙ_００」で表し、非周期入力信号ｕ_００は図１０（ｂ）中の「ｕ_００」で表した。

オーバーシュート１０％以下を拘束条件として、定常値の２％以内への整定時間を評価関数とし、ＭＡＴＬＡＢのｆｍｉｎｃｏｎによって評価関数を最小化した。具体的には、（２６）式のＫ（θ）を探す最適化問題を解いた。

調整パラメータθ＝（Ｋ_ｐ、Ｋ_ｉ、Ｋ_ｄ）とし、非周期出力信号ｙ_００を取得した時のＰＩＤパラメータを初期値とした。非周期出力信号ｙ_００に定常値の２％を超えるノイズが重畳していてこのままでは整定しないため、仮想出力信号ｙを（２７）式の移動平均フィルタに通してから整定時間とオーバーシュートを測定した。

図１０、図１１に最適化して得られたＫ（θ）を用いた時のフィードバック系の真の時間応答（図１０（ａ）中で「Ｔｒｕｅ ｙ」及び図１１中で「Ｔｒｕｅ ｕ」）と仮想時間応答（図１０（ａ）中で「Ｖｉｒｔｕａｌ ｙ」及び図１０（ｂ）中で「Ｖｉｒｔｕａｌ ｕ」）を示す。なお、真の時間応答とは、仮想時間応答を用いて最適化したＫ（θ）の制御器でモータを駆動させた場合のモータの入力電圧ｕと回転数ｙである。

同図より非周期出力信号ｙ_００に比べ仮想時間応答と真の応答はオーバーシュートが抑制され、整定時間が小さくなった。しかし、仮想時間応答にはおもに約５Ｈｚの振動が生じたが、これは真の出力応答には見られず、真の入力応答に高周波振動が発生した。このような誤差があるものの、周波数成分が豊富ではないステップ応答を基にしたものにも関わらず、本発明に係る出力推定方法を用いて、制御器のパラメータを最適化することで制御性能が向上されることを確認できた。

＜捕捉＞
（９）式と（１１）式において、ある角周波数ωａにおいて、ｒ_１（ｊω）＝０のとき、ゼロ割が生じる場合がある。その時はｙ（ｊω）＝ｕ（ｊω）＝０と設定すればよい。以下詳説する。
図３（ｃ）のフィードバック系の目標入力がｒ_１（ｊω）＝０であるので、ｙ_０（ｊω）、ｕ_０（ｊω）について（２８）式、（２９）式を得る。

（２８）式より［Ｇ（ｊω_ａ）Ｋ（ｊω_ａ）］^−１＋１≠０ならばｙ_０（ｊω_ａ）＝０である。ナイキストの安定判別より、安定余裕がゼロとなる安定限界のときに限りＧ（ｊω_ａ）Ｋ（ｊω_ａ）＝−１となる。しかし、仮定Ａ３よりフィードバック系が安定なのでそうはならない。したがってｙ_０（ｊω_ａ）＝０を得る。

（２９）式よりＫ^−１（ｊω_ａ）＋Ｇ（ｊω_ａ）≠０ならばｕ_０（ｊω_ａ）＝０である。Ｋ^−１（ｊω_ａ）＋Ｇ（ｊω_ａ）＝０となるのは、Ｇ（ｊω_ａ）Ｋ（ｊω_ａ）＝−１またはＫ^−１（ｊω_ａ）＝Ｇ（ｊω_ａ）＝０のときだけである。

仮定Ａ３よりＧ（ｊω_ａ）Ｋ（ｊω_ａ）≠−１であり、仮定Ａ５よりＧ（ｊω_ａ）＝０のときＫ^−１（ｊω_ａ）≠０であるので、ｕ_０（ｊω_ａ）＝０を得る。よって仮定Ａ２よりｒ（ｊω_ａ）＝０である。

ｒ（ｊω_ａ）＝０のとき、図３（ｂ）に対して同様の議論を行うと、ｙ（ｊω_ａ）＝ｕ（ｊω_ａ）＝０を得る。

本発明に係る出力推定方法は、制御対象への試行のための信号入力が制限されているような場合に有用である。

Claims (3)

1. 制御器Ｋを介して制御対象への入力に対する出力を推定する方法であって、
   １回の試行に基づいて制御対象への周期入出力応答（周期入力信号ｕ_０と周期出力信号ｙ_０）を求める工程と、
   仮想入力信号ｒに対する前記制御対象の応答（仮想出力信号ｙ）を（１）式で求める工程を有することを特徴とする出力推定方法。
   

   

   ここで、ｒ_１（ｊω）は、（２）式に従う。ｒ_１（ｊω）、ｒ（ｊω）、ｕ_０（ｊω）、ｙ_０（ｊω）は、それぞれ信号ｒ_１、ｒ、ｕ_０、ｙ_０をフーリエ変換したものを表す。またＫ^−１（ｊω）は１／Ｋ（ｊω）を表し、Ｋ（ｊω）は制御器Ｋの伝達関数を表す。
2. 前記周期入出力応答を求める工程は、前記１回の試行で、前記制御対象へ非周期入力信号ｕ_００を入力し、非周期出力信号ｙ_００を取得し、前記非周期信号（非周期入力信号ｕ_００、非周期出力信号ｙ_００）を周期信号（周期入力信号ｕ_０、周期出力信号ｙ_０）に変換する工程であることを特徴とする請求項１に記載された出力推定方法。
3. 制御対象を制御する制御器のパラメータをＫ（θ）として、請求項１または２の何れかの請求項に記載された出力推定方法による前記仮想出力信号および仮想入力信号を用いて、Ｋ（θ）を最適化する制御器の最適化方法。

Images