特許 6989099 周波数不感帯を抑制する周波数応答解析システム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6989099

(24)【登録日】2021年12月6日

(45)【発行日】2022年1月5日

(54)【発明の名称】周波数不感帯を抑制する周波数応答解析システム

(51)【国際特許分類】

   G05B 11/36 20060101AFI20211220BHJP

   G05B 23/02 20060101ALI20211220BHJP

【ＦＩ】

   G05B11/36 501C

   G05B23/02 P

【請求項の数】1

【全頁数】11

(21)【出願番号】特願2017-40246(P2017-40246)

(22)【出願日】2017年3月3日

(65)【公開番号】特開2018-147149(P2018-147149A)

(43)【公開日】2018年9月20日

【審査請求日】2020年2月28日

(73)【特許権者】

【識別番号】304021277

【氏名又は名称】国立大学法人 名古屋工業大学

(72)【発明者】

【氏名】前田 佳弘

(72)【発明者】

【氏名】立花 弘貴

【審査官】 影山 直洋

(56)【参考文献】

【文献】 特開２００５−２１４７１１（ＪＰ，Ａ）

【文献】 特開２０１６−０９０６７０（ＪＰ，Ａ）

【文献】 特開２０００−２７８９９０（ＪＰ，Ａ）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０５Ｂ １１／３６

Ｇ０５Ｂ ２３／０２

(57)【特許請求の範囲】

【請求項1】

プラント（３）の周波数特性を算出するＦＲＡシステム（２）において、
前記プラント（３）の力次元の入力（Ｕ）に対して、
第１差分器（１４）を介して現在のサンプル時刻と１サンプル前における入力の差分値Ｕ’を算出し、
続いて第１ＤＦＴ演算部（１６）を介して周波数応答Ｕ’（ω）を算出し、
前記プラント（３）の位置次元の出力（Ｙ）に対して、
第２差分器（１５）を介して現在のサンプル時刻と１サンプル前における出力の差分値Ｙ’を算出し、
続いて第２ＤＦＴ演算部（１７）を介して周波数応答Ｙ’（ω）を算出し、
乗除算器（１８）において、前記周波数応答Ｕ’（ω）及び前記周波数応答Ｙ’（ω）より前記プラント（３）の周波数特性を算出することを特徴とするＦＲＡシステム（２）。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、プラントの周波数特性を同定するＦＲＡ（周波数応答解析）システムに関する。

【背景技術】

【0002】

通常サーボシステムは共振モードを有しており、「稼働中」のプラント入出力信号を用いてプラントの周波数応答解析の結果を得て作動状況を確認し制御することは、サーボシステムを位置決め性能等を確保して運用するために重要である。

【0003】

サーボシステムに対する、線形ダイナミクスの周波数特性を同定する周波数応答解析のアルゴリズムに関する従来技術は主に２つある。

【0004】

従来技術１は、例えば特許文献１のようにサーボアナライザなどの周波数応答解析（以下、ＦＲＡと記す）を行うアルゴリズムを備えた機器を用いて、プラントを含む制御系に正弦波信号を掃引し、その際のプラント入出力信号を離散フーリエ変換（以下、ＤＦＴ）することで周波数応答解析結果を得るものである。本技術は、掃引周波数に対して高いＳ／Ｎ比の入出力信号を得るために、サーボシステムが稼働していない状態で正弦波掃引試験を実施することを前提とした技術である。従って、仮に稼働中に正弦波信号を掃引した場合、それに起因する位置決め応答劣化などの問題が生じる。

【0005】

一方、例えば非特許文献２のように離散フーリエ変換（以下、ＤＦＴと記す）によって、稼働中の入出力信号からＦＲＡを行う従来技術２においては、単純にＤＦＴを実施するだけでは入出力信号のＤＦＴ結果に表れる周期的な周波数不感帯（周期的不感帯）によって、ＦＲＡ結果において実際には存在しない周期的なスパイク状の誤差（ＦＲＡ誤差）が生じる。これにより、プラントモデルパラメータの同定精度が劣化する、モデルベース制御による位置決め性能が劣化する、という問題点を有する。

【先行技術文献】

【特許文献】

【0006】

【特許文献1】特開平８−０９４６９０

【非特許文献】

【0007】

【非特許文献1】電気学会論文誌Ｃ, 「閉ループステップ応答データを用いた機械系の周波数応答推定」, Vol.131, No.4, pp.751-757, 2011

【発明の概要】

【発明が解決しようとする課題】

【0008】

本発明の課題は、「稼働中」の入出力信号を用いたＦＲＡ技術を前提として、ＤＦＴにより生じるスパイク状のＦＲＡ誤差を発生させない、プラントの周波数特性を同定するＦＲＡシステムを提供することである。

【課題を解決するための手段】

【0009】

課題を解決するために本発明は、プラント（３）の周波数特性を同定するＦＲＡシステム（２）において、プラント（３）の力次元の入力（Ｕ）に対して、第１差分器（１４）を介して差分値Ｕ’を算出し、続いて第１ＤＦＴ演算部（１６）を介して周波数応答Ｕ’（ω）を算出し、プラント（３）の位置次元の出力（Ｙ）に対して、第２差分器（１５）を介して差分値Ｙ’を算出し、続いて第２ＤＦＴ演算部（１７）を介して周波数応答Ｙ’（ω）を算出し、乗除算器（１８）において、周波数応答Ｕ’（ω）及び周波数応答Ｙ’（ω）よりプラント（３）の周波数特性を算出することを特徴とするＦＲＡシステム（２）である。
なお、上記及び特許請求の範囲における括弧内の符号は、特許請求の範囲に記載された用語と後述の実施形態に記載されて当該用語の例となる具体物等との対応関係を示すものである。

【発明の効果】

【0010】

発明によると、サーボシステムが「稼働中」の状態でもＦＲＡを行えるので、ＦＲＡのためにサーボシステムを停止する必要がなく、サーボシステムの稼働率を低下させることはない。また、稼働中のプラント周波数特性をリアルタイムに観測し、また、ＦＲＡ精度を飛躍的に向上させてＦＲＡを同定できるため、プラント特性変動に起因する位置決め性能劣化を抑制できる。更に、故障が生じる前に異常診断でき制御器更新ができる等の付加価値機能の追加が可能となる。また、従来技術２に対して差分器２個の追加で実現でき、低実装・低演算コストでＦＲＡ精度を飛躍的に向上できる。

【図面の簡単な説明】

【0011】

【図1】本実施形態のＦＲＡシステム構成を示すブロック線図。

【図2】本実施形態の入力の時間応答シミュレーション結果。

【図3】本実施形態の入力差分値の時間応答シミュレーション結果。

【図4】本実施形態の出力の時間応答シミュレーション結果。

【図5】本実施形態の出力差分値の時間応答シミュレーション結果。

【図6】本実施形態の入力差分値の周波数応答シミュレーション結果。

【図7】本実施形態の出力差分値の周波数応答シミュレーション結果。

【図8】本実施形態の推定プラント周波数応答シミュレーション結果。

【図9】従来技術２の推定プラント周波数応答シミュレーション結果。

【図10】本実施形態のプラントモデル同定シミュレーション結果。

【図11】従来技術２のプラントモデル同定シミュレーション結果。

【図12】本実施形態のＦＲＡ結果から同定したプラントモデルを制御器設計モデルに用いた位置決め応答シミュレーション結果（比較：従来技術２）。

【図13】従来技術１のＦＲＡシステム構成を示すブロック線図。

【図14】従来技術２のＦＲＡシステム構成を示すブロック線図。

【図15】従来技術２の入力の周波数応答シミュレーション結果。

【図16】従来技術２の出力の周波数応答シミュレーション結果。

【発明を実施するための形態】

【0012】

以下、図面を参照しつつ本発明の実施の形態について説明する。本発明は、以下の実施形態に限定されるものではなく、発明の範囲を逸脱しない限りにおいて、変更、修正、改良を加え得るものである。
また、図８から図１２において、データ線に符号（ア）、（イ）、（ウ）、（エ）、（オ）、（カ）、（キ）を表示しているが、同一符号は同一データ線を示す。

【0013】

（従来技術の課題１）
図１３は、従来技術１のＦＲＡシステム（２００）の構成を示すブロック線図である。位置指令を零（Ｒ＝０）としたサーボシステムに対し、ＦＲＡシステムから加算器２０１を介して加振信号として正弦波Ｕｓｉｎ（２０２）を制御系に入力し、プラントの入力Ｕ、出力ＹからＤＦＴによって入力の周波数応答Ｕ（ω）（２０３）、出力の周波数応答Ｙ（ω）（２０４）を算出する。そして、Ｕ（ω）とＹ（ω）を用いて推定プラント周波数特性Ｐｅ（ω）を求める。ＤＦＴ演算部の前段における差分器は有しない。原理的に、ＦＲＡのために加振信号を掃引するため、稼働中に実行すると位置決め性能劣化を招く。

【0014】

加振信号ｕｓｉｎは、単一周波数の正弦波を用いたが、これに限らずＭ系列や合成正弦波等の複数の周波数成分から構成される信号でもよい。

【0015】

（従来技術の課題２）
図１４は、従来技術２のＦＲＡシステム３００の構成を示すブロック線図である。本発明と同様に、稼働中のプラントの入力Ｕ、出力ＹからＤＦＴによって入力の周波数応答Ｕ（ω）（３０１）、出力の周波数応答Ｙ（ω）（３０２）を算出する。そして、Ｕ（ω）とＹ（ω）を用いて推定プラント周波数特性Ｐｅ（ω）を求める。ＤＦＴ演算部の前段における差分器は有しない。

【0016】

従来技術２において、ＤＦＴ演算に用いる入力Ｕと出力Ｙは、繰返し信号（初期値と最終値が一致した時系列信号）とするため、実施形態に示す数式（２）で定義される過渡応答時系列信号を用いてＮ＋１サンプル以降で折り返した信号（データ数２Ｎ）をＵ、Ｙとして用いる。ＵとＹは、それぞれ数式（９）、（１０）として定義される。

【数9】

【数10】

【0017】

図１５に従来技術２における入力の周波数応答Ｕ（ω）、図１６に従来技術２における出力の周波数応答Ｙ（ω）を示す。両図において、周期的不感帯（一定の周波数間隔でゲインが急峻に減少する）が見られる。これが入出力信号のＤＦＴ結果に表れる周期的な周波数不感帯（周期的不感帯）であり、これによりＦＲＡ結果において実際には存在しないスパイク状の誤差（ＦＲＡ誤差）として生じる。

【0018】

ここで、数式（９）に対するＤＦＴ変換は数式（１１）で表せる。

【数11】

数式（３）と同様に、入力Ｕの差分値Ｕ’をデータ長２Ｎに対して定義するとき、数式（１１）は数式（９）の関係を考慮して、数式（１２）のように変形できる。

【数12】

ここで、数式（１２）の１つめの指数関数の和の部分を、正弦関数を用いて表すと数式（１３）となる。

【数13】

分子の項に着目すると、Ｕ（ω）はω＝ｐ／ＮＴ（ｐ＝０，１，・・・）の周期で零となり、これが図１５における周期的不感帯の要因である。なお、出力の周波数応答Ｙ（ω）についても同様である。

【0019】

Ｕ（ω）とＹ（ω）を用いて、推定プラント周波数特性は数式（１４）で求められる。

【数14】

分子と分母の１つめの指数関数の和の部分は、理論的には相殺されるが、実際のＤＦＴを行う計算機は有限の演算精度となる。これによって、理想的に相殺されなかった影響として、図９の太破線（ウ）に見られるスパイク状のＦＲＡ誤差が生じる。

【0020】

数式（１４）に着目すると、分子と分母の１つめの指数関数の和の部分が理想的に相殺される場合、分子分母共に２つめの指数関数の和の部分のみが残る。それらは、入力、出力のそれぞれの差分値のＤＦＴ結果を示しおり、数式（６）、（７）、（８）で表される本発明と等価の式となる。すなわち、本発明における差分値に対するＤＦＴは、スパイク状のＦＲＡ誤差を招く数値演算を根本的に排除している。

【0021】

（実施形態）
サーボシステムは、プラント部（可動部を備える揺動あるいは直動機械、可動部に対する推力を発生するアクチュエータ、アクチュエータを駆動するためのサーボアンプ及び可動部位置を検出するセンサ等）と、コントローラ部（制御演算を行う補償器）と、周波数応答解析（以下、ＦＲＡ）部と、によって構成される。
本実施形態は、このようなメカトロニクス機器のプラント線形ダイナミクスに対する、稼働中の入出力信号を用いた周波数応答解析アルゴリズムに関するものである。

【0022】

図１は、ＦＲＡシステムの実施形態のブロック線図である。サーボシステム１のＦＲＡシステム２は、ＤＦＴ演算部１６、１７と、差分器１４、１５と、乗除算器１８を有する。プラント部とコントローラ部は、プラント３と、制御演算を行うフィードバック（ＦＢ）補償器Ｃ（ｚ）４と、フィードフォワード（ＦＦ）補償器ＦＦｃｏｍｐ５と、を有する。プラント３には、機器部のサーボアンプやセンサに起因する入出力の遅れ要素、機械系における共振モードなどの線形ダイナミクスが含まれる。
サーボシステム１に対して、位置指令Ｒ（８）が入力されると、ＦＦ補償器５は制御アルゴリズムに基づいてＦＦ位置指令軌道Ｒｆｆ（９）とＦＦ入力Ｕｆｆ（１０）を算出する。ＦＦ位置指令軌道Ｒｆｆは、第一加算器１２に入力され、センサによって検出される位置出力Ｙ（７）との差分が計算されてＦＢ補償器Ｃ（ｚ）に入力される。ＦＢ補償器Ｃ（ｚ）は、制御アルゴリズムに基づき位置の次元から力の次元に変換し、算出されるＦＢ入力Ｕｆｂ（１１）は、第二加算器１３に入力され、ＦＦ入力Ｕｆｆとの加算によってプラント３への入力Ｕ（６）が計算される。プラント３は、入力Ｕによって位置出力Ｙを生じる。入力Ｕは、第一差分器１４へ入力され、入力差分値Ｕ’ （１９）を出力し、第一ＤＦＴ演算部（１６）を介して入力差分値Ｕ’の周波数応答Ｕ’（ω）を算出する。同様に、出力Ｙは、第二差分器（１５）へ入力され、出力差分値Ｙ’ （２１）を出力し、第二ＤＦＴ演算部１７を介して出力差分値Ｙ’の周波数応答Ｙ’（ω）を算出する。乗除算器１８を用いて、Ｙ’（ω）に対するＵ’（ω）の除算値を計算し、推定プラント特性Ｐｅ（ω）（２３）となる。
本実施形態は、ＦＲＡシステム２において、入力Ｕと位置Ｙから差分値Ｕ’とＹ’が出力され、ＤＦＴ演算部によってそれぞれのＤＦＴであるＵ’（ω）とＹ’（ω）を出力する。Ｕ’（ω）とＹ’（ω）は、後述の原理によってＦＲＡ誤差を招くスパイク状の周期的誤差は除去されており、Ｕ’（ω）とＹ’（ω）を用いて推定プラント周波数特性Ｐｅ（ω）を計算することで、線形プラント特性に対する高精度なＦＲＡアルゴリズムとなる。

【0023】

プラント特性Ｐ（ｚ）は、共振モードを有する離散時間伝達関数であり、例えば連続時間伝達関数表現で数式（１）のように示される。

【数1】

【0024】

ＦＲＡシステム１２０にて、力次元の入力Ｕの各サンプル時刻（サンプル周期Ｔ）における値を数式（２）のように表すとき（データ長Ｎ）、Ｕと第一差分器１１１との積として、差分値Ｕ’は数式（３）のように算出される。
図２に入力Ｕの、図３に入力差分値Ｕ’の時間応答シミュレーション結果を示す。

【数2】

【数3】

【0025】

同様に、位置次元の出力Ｙの各サンプル時刻における値を数式（４）のように表すとき、Ｙと第二差分器１１１との積として、差分値Ｙ’は数式（５）のように算出される。
図４に出力Ｙの、図５に出力差分値Ｙ’の時間応答シミュレーション結果を示す。

【数4】

【数5】

【0026】

その後、第一ＤＦＴ演算部１６と第二ＤＦＴ演算部１７において、Ｕ’の周波数応答Ｕ’（ω）とＹ’の周波数応答Ｙ’（ω）を、数式（６）と（７）に従ってそれぞれ算出する。

【数6】

【数7】

図６にＵ’（ω）の、図７にＹ’（ω）の周波数応答シミュレーション結果を示す。
従来技術２である図１５に示すＵ（ω）、図１６に示すＹ（ω）と比較すると、周期的不感帯（一定の周波数間隔でゲインが急峻に減少する）が発生していないことが分かる。

【0027】

乗除算器１１３によりＹ’（ω）に対するＵ’（ω）の除算として推定プラント周波数応答Ｐｅ（ω）を求める。この関係を数式（８）に示す。

【数8】

【0028】

図８は、本実施形態の推定プラント周波数特性のシミュレーション結果（Ｐｅ（ω）のゲイン）を示す。（ア）の太破線は、同定すべきプラント特性Ｐ（ω）であり、目標値に相当する。本実施形態の結果を、（イ）の細実線で示す。（ア）及び（イ）は、ゲイン、位相共に良く一致しており、一本の線のように見える。すわなち、稼働中の入出力信号を用いて高精度なＦＲＡ結果を得ている。

【0029】

一方、稼働時におけるＦＲＡ法の図８の比較として、図９に従来技術２による推定プラント周波数特性のシミュレーション結果（Ｐｅ（ω）のゲイン）を示す。（ア）の太破線の同定すべきプラント特性Ｐ（ω）に対して、（ウ）の細実線で示すＦＲＡ結果は周期的なスパイク状のＦＲＡ誤差が生じ、（ア）と（ウ）は広い周波数範囲に渡って一致していない。

【0030】

以上、本実施形態によれば、稼働中の入出力信号を用いたＦＲＡにおいて、周期的なスパイク状誤差が発生しない。よって、プラント周波数特性を高精度に同定するＦＲＡシステムを提供することができる。また、ＦＲＡのためにサーボシステムを停止する必要がなく、サーボシステムの稼働率を低下させることはない。更に、故障が生じる前に異常診断でき制御器更新ができる等の付加価値機能の追加が可能となる。また、従来技術２に対して差分器２個の追加で実現でき、低実装・低演算コストでＦＲＡ精度を飛躍的に向上できる。

【0031】

図１０は、本実施形態のＦＲＡ結果に対して、偏分反復法などの最適化アルゴリズムを用いて同定した、数式（１）のプラントモデル伝達関数の周波数特性シミュレーション結果を示している。（イ）の太破線は本実施形態によるプラントのＦＲＡ結果Ｐｅ（ω）、（エ）の細実線は同定したプラントモデルである。（イ）及び（エ）はよく一致し、一本の線のように見える。これは、図で示されるように、Ｐｅ（ω）において周期的なＦＲＡ誤差が抑制されているからである。

【0032】

一方、稼働時におけるＦＲＡ法の図１０と比較として、図１１に従来技術２のＦＲＡ結果に対して同定した数式（１）のプラントモデル伝達関数の周波数特性シミュレーション結果を示す。（ウ）の太破線は本実施形態によるプラントのＦＲＡ結果Ｐｅ（ω）、（オ）の細実線は同定したプラントモデルである。（ウ）における周期的なＦＲＡ誤差により（オ）は共振モードに対して明確な誤差を生じている。これは、（ウ）における共振モード付近のスパイク状のＦＲＡ結果に対して、モデルパラメータ同定アルゴリズムが局所解に陥りやすくなるためである。

【0033】

図１２は、本実施形態のＦＲＡ結果に対して同定したプラントモデルを、図１のＦＦ制御器の設計モデルに用いた、位置決め応答シミュレーション結果を示す。比較として、従来技術２の場合を示す。（カ）の本実施形態では、共振モードに起因する残留振動を抑制し、水平破線で示される目標整定精度内に収束している。これは、図１０で示されるように、同定すべきプラント周波数特性Ｐ（ω）を、本実施形態が良く同定できるためである。
一方、（キ）の従来技術２では、残留振動を伴うオ−バーシュート及びアンダーシュートが発生している。これは、同定すべきプラント周波数特性Ｐ（ω）を、従来技術２は精度良く同定できないからである。このような同定誤差が生じた場合、サーボシステムの位置決め精度は劣化する。よって、プラント周波数特性Ｐ（ω）に対して高精度なモデルを獲得するために、数値シミュレーションや実験による試行錯誤的なモデルパラメータの調整・更新に多大な労力を必要とし、その人的労力や時間も多大となる。

【産業上の利用可能性】

【0034】

周波数応答解析機器を扱うメーカにとって、精密サーボシステム等で問題視されている正弦波掃引時と稼働中のプラントダイナミクス変化に対し、稼働中の入出力信号を用いたオンラインＦＲＡ技術として有望である。

【0035】

精密サーボシステムを扱うメーカにとって、サーボシミュレーションや有効な制御器設計に用いるプラントモデルの高精度同定手法として、合わせ込みのための試行錯誤的なシミュレーション及び実験に係る工数削減と、サーボ性能の向上技術として有望である。

【符号の説明】

【0036】

１ サーボシステム
２ ＦＲＡシステム
３ プラントＰ（ｚ）
４ ＦＢ補償器Ｃ（ｚ）
５ ＦＦ補償器ＦＦｃｏｍｐ
１１ 第１加算器
１３ 第２加算器
１４ 第１差分器：（ｚ−１）／ｚ
１５ 第２差分器：（ｚ−１）／ｚ
１６ 第１ＤＦＴ演算部
１７ 第２ＤＦＴ演算部
１８ 乗除算器

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】