超低周波雑音を含む操業データからの動特性推定
（昭和51年5月31日 原亭己i怪付）
制御工学教室古谷忠義
Identificatioll of Linear Systems in the Presence of Very Low Frequency Noise
by Tadayoshi FURUYA
If・nly伽ct岨ti・g P・・t・f i叩・t・・d・・tp・t・ignal i・amphfi・d f・・d・t・mi。、ti。，。f dyn、mi。
characteristics． very low frequency components are relatively large， which caused by drift of instru−
ment・and・。・・t・U・n・・y・1迦g・・f・y・t・m va・1・bl・． Th・・ef…a1・ng P・・i・d d・t・，h。。ld b，，eq。i−
red to esUmate parameters of dynamic characteristics．
This paper is concerned with the estimation of parameteTs of dynamic characteristics ln the P「e・en…f…yl。wf・eque・cy・。i・e・N・isei・・ep・・at・di・t。tw・P・・t・，・n・。fwhi・hi、v。，yl。w f・equency p・・L・・d th・・th・・i・蜘d・m・。ise p・・t． V・・y l。w f・eq・en・y p・・t田。 b， apP，。。im。t，d by p・ly・・mi・1・・f、time…i・bl・a・d…d・m・・i・・p・・t can b・・ep・e・e・t・d by aut。，，g，e、、ive f。，m．
We proposed an iterative algorithm for e5timation of parameter呂of dynamic characteristics and noise
・ha・a・t・・i・ti…It i・p。ssibl・t・d・t・・mine dy・ami…f・ub・・y・t・m i・th・feed b、、k l。。P、nd i。 tl1。
presenceof very low frequency compone耐s，since both of dynamics and noise chal acteristics are estimated．
Experimental results are included to demonst口te the fea5ibility of the proposed method．
1．緒 言 増幅すれば非定常性の外乱や計測器のドリフト雑音が相
対的に大きくなり．限られた時聞内においても，これら 最近のようにシステム運用に対して厳しい仕様の最適 の超低周波成分による影冊が顕著である，この対策とし 化や最適制御を行なうことが課せられるようになってく て平均値変動とも考えられる超低周波成分の除去を時間 ると・一般の工業プロセスでも常時制御系を最適に保つ 平均をとることによって行なうと｝乏時間のデ＿タを必要 ために計算機を導入し・オン・ラインでデ・一タを計算機 とし，逆に短時間のデータから動特性推定を行なうと超 に読み込み・データ処理を行なうことが普通になってき 低周波以外の雑音によって推定値にぱらつきが生じる。た・さらに高度にシステムを運用するためには静的な特性 したがって超低周波成分が大きいデータからの動特性推 のみならず動的な特性ことにシステム内部のサブ・シス 定は困難と考えられてきた。
テムの動特性やシステムに入ってくる外乱の特性を従来 超低周波成分の除去を行なう動特性推定法として高域 よ1〕正確に推定する必要がある。 i戸波器を用いる方法とW．D、 T． Davisコ）やNikiforvk3）ら
ところが運用状態にあるシステムにおいては監枇や が提案しているような峙間の多項式を用いる方法等があ データ集録に設けられている入出力の計装設備から得ら る．これらの方法は超低周波雑i昔を除去することには優 れるシステム変数の時間的変劃⑦は僅少であるうえに動特 れているが，通常の不規則雑音の除去すなわち推定が不 性推定に関係のない低周波成分が大部分であってどうし 十分であり，推定値の分散が大きいという欠点がある。
ても平均値成分を取リ去った変動成分を増幅し取り出さ さらに内部雑枠が存在するシステムのサブ・システムの ねばならない。11これは計測器などの感度限界の増大と 動特性推定においては推定値が不偏推定植とならないこ 計算機入力のためのAD変換器の鼠子化誤差を減少さ とかある．
せるたあに不可欠な手段である。しかし変動成分のみを これに対して通常の不規則雑音だけが存在する埼合に
100
は動特性の最小分散不握推定値を与える方法が数多く提 となる。ここでTをサンプル周期として
蕊蕊鷲を當㌶㌶：熟 力・；ぱ・（・）輌一＝煎醐
y缶二y（力丁）． 恥＝m（占丁）， 虜h＝，鵡τ）
いる。しかし通常の不規則雑音のみを仮定しているので
超低周灘音の存在する場合にはやはり長時聞のテータ である・②式につr5も喘畷散化すれば
を躍とする欠点がある． 炉゜・＋尼α劇 〔4）
欄文で述べ肪蹴ド・1フトのような超低周波麺 となる・
、 H｛のの離散値π▲は の存在する場合にも適用可能なよっに拡張した ものであ
る。ここでは超低周波成分を時間の多項式で近似し，そ 識㈲＝Φ識ほ一D＋b〜白 ㈲ のパラメータと動特性を同時に推定することによって超 加＝C転H㈲
低周波成分の除去を行なっている。 ここでψは還移行列である。
、 ＿ と書きあらわすことができる。15拭をZ変換し，整理する 2． ソステム・モァル
と
ここで対象とするシステムは線形であると仮定する。 識（ZX1一砂ガ1〕＝b1U（Z） ｛6｝
システムに重畳している雑音はその周波数成分によって 識（z）＝〔∬一φガ1） 1占蝋z）
高周波成分，入出力変動と同帯域の雑音とドリフトのよ となる。（7）式を㈲式に代入すると
うな超低周波成分に分けて考えることができる。このう 川訂＝c「（∬＿φガリー膓6zσ（訂 ｛7｝
ち高周波成分はアナログ・フィルタを挿入して連断する ＝B（三りτo（z）
ことができるカC超低周波成分の混入してくることは避 双z1）
けられない。したがってここではシステム．モデルを となる・A輌〒）カζ示しているように・ωが麟であれ
瞳1のよう蹴る・
@ ぱ。（。）＿1犯（。） ｛副
1一Σα，z f ホ
H田5E と仮定できる。すなわち堀f｝の離散値恥は ガロ
ぬニΣ四〜白一∫＋lu白 ｛9）
DRIF了 ∫■1
占
なる自己回帰モデルで生起しているものと仮定する。こ ・ 耐PUT Y、 こで職は正規白色性の雑音である。動特性推定の問題 脈 Su8−SY5榔 は与えられた入出力データ加と臨からパラメータ∫〜 ， Fi昏1 Mathematical Model 凸，疏を推定することである。
システムの入力を姻出力を。ω超｛鍋麟音を 3・酬雌定法
mω不規則性雑音を力（ ）とあらわすと入出力の関係は 3−1 原理
インバル碑答、ωを用いて 超低周灘音・蠕齢よぴ雑齢性それぞれをあち
・吋㌔（・輌一）＋川 ）｛1｝㌶二1の推定方法を考える・③式蹴（91式
と書くことができる・超低周灘音mω1塒間に閥す ，，、一。一堂∫。。』一呈繊τP壌訓、．、llO｝
る多項式 司 」司 国
杣＝。．増耐 ② の関係力fあり・…を正規白雛音と仮定して・ るので雌
」叶 尤推定を行なうことができる。尤度関数から で近似的にあらわす。ここでロiが超低周波雑音の時間変
化をあらわ抹知パラメータである。 ノー菖（・菖声・一璃・・冊一覧一一・アω
①式を嵯散値系に書き換えると を最小とすれば良い。ところが川｝式にはα、と♪島．、のとも ぼ
治＝Σん臨． ＋伽＋取 口｝ に未知項の植が存在しているので通常の線形推定を行な r口t
うことができない。そこで・なんらかの方法で得られた く得られたβ1，の差のノルムを規準化したもの 岳，品を用いて｛3｝，14｝式を変形した
‖βκ一β汗一1｛1／llβκll ｛151
エぞ ロ バ
瓜＝撫一君・ノr匡f 白一一品』男島伎7『）」 ｛121 が十分に0に近ければ収束したものとしてStep
から，‥の推定値茄を求める。そして評価関係ノのなかの 5に進み・収束していなければStep 2に戻る・
，1、を厄でi置き換えた」綱職 St・p 5・収束したβ棚の動雛と雑音繊をあ肋
炉鷺（玖菖∫・臨一蓋醐丁巨貢缶五』」〔13） ㍍㍍：窺灘墓：㌶：鷲
求め1 古が最小となるよコに疏，あ，畠を計算する。
以上の計算を収束するまで繰り返すことによって動特
が最小となるように岳，あ，亘、を決定する。この∫証は し祖・が正規白色性でなければパラメータ個数 パラメータに閲して線形であるので通常の最小2乗法と 」Vρ・N九鵡を変えて計算を繰り返す・
榊・パラメータを勅ることができる． 以上の計算手順によって糠剛撲∫・躍舗幽
次に繰り返し計算を考える。す軸ち∫。を最小とす と超低周波齢のを求めることができる・ ・
るノ2、、元｝．岳のうち∫〜．，命を用いて｛12）式から再び取を
、 ＿ ＿ START
性鵡剖描超酬灘勧それぞオ、のパラメータを DAT融UT
決定する。
3−2澗方法 育1・・倉1…⊇・Y・
の計筋法から考える．通常の最小2乗法と同様に（1ヨ｝式 C°HPUTほX』DヱX・YK
まず㈹式を最小とするように抗．蚕Lδ，を求めるため ．Tを最小とする届，め，応は次の線形方程式によって与え ［H端E
られる。 SOLVE Eo・IF Np NR No
リロ
ふ ば ふ ば
碧1Xぷ三江1≧ 〔1▲｝ n・一日N−ll光・・ C旺CK
ここで 1「部1 0K
」rl『＝〔r ・．曲．，…f 、．。・，i1，々丁，…（片丁〕・ i ヨ、．1…η、．Nα〕
YES STOPβ丁＝〔」；1∬元… ；．，…緬1…亘．，1δ1…δ。。〕 ［＾LCULATE W・u・1・G E・．uq）
したがってこの線形代数方程式を解けばパラメータ が 決定される。この計算において命の値は繰り返し計算
ごと1・ 2｝式から勅ら．れるのでX．xrと㌔、の汀、1、 醜2趾㎞ati°n A』・・i輪
聞する部分は計算が繰り返されるごとに異っているが，汀・の酬しない部分は常に同＿であるので繰り返し酬
4・実酬から除くことができる。 以上述べた方法によって，動特性推定の摸擬実験を行 以上の計算手順をまとめてみるとFig．2に示すよう なった。実験データは、ディジタル計口：機かアナログ計 になる。 算機によってシステムを膜擬1し，作製している。精度や
Step l．まず｛12）式の 抗，雨をすぺて0と置いて 求められた数値が亜要であるものについてはディジタル 瓜を求め｛14｝式から1回目の推定値βを得る。 計算機により行ない，一方統計的な理論と現実のシステ
Step 2．βのパラメータを用いて汀，を求める。 ムに差異の生じることが考えられる場合にはアナログ計 Step 3．再び（10式のΣ》5X「とΣ］X．治における ・ 算機を用いた。
ぬの関係する部分を計算し，llO式の代数方程式を 実験は圭に次の3点について行なわれている。まず超
解き，新しく鼻を得る。 低周披雑音の影韓とその除去，つぎにシステムにワィー
Step 4．収束の判定を行う。すなわち前回のβヤー1と新し ド・パック径路があり．その内部に雑音が存在している
102
場合その内部雑音の影響とその除去である。最後に超低 力をデータ処理装置によって計算機に読み込んだ。雑i音 周波雑音と不規則雑音がともに存在している場合につい と入力の分散の比を0．47になるようにκ∫を決定し・
てである。 サンプル周期を2秒，データ個数を300個に選び・実験 4−1 超低周波成分の除去 を行なった。パラメータ個数は残差が白色に近いように
以上述べた方法によって，超低周披成分が除去される 決定し，Nβ＝10， N4＝4，」V・＝0とし・推定した結果ことを示すために動特性推定を行なった。Fig．3に示す がFig．5である。通常の最小2乗法すなわち1鴫＝oと システムをディジタル計算機によって模擬し，入出力 した場合の推定値を参考のために示してある。この図か データとした。ここでは超低周波雑音η のみが重畳し らあきらかなように通常の最小2乗法では推定値にパイ ているものとし．m身は アスが存在している。これに対して，ここで提案した方 ，〃白＝0．0助丁＿5．10−5伍丁）コ 法によると十分な精度で推定できる。
と仮定した。このとき入力の分散は2．5である。
1脚1 °二：7，＿
1 1 Y貝 田目川輪15【 冊 UR 示「
Fi9．3 Sim血ted Sy8tem 鞠 1苫 ． Table lはこれらのパラメータの真値と推定値を示し
たものである，以上の結果から超低周波成分の影響の大
なること，多項式近似によって超低周波成分が除去され Fig・4 Simurated Sub8y8tem in the Feed B血￠k LooP ることがわかる。
T且BLE l
E・TI酪…V・…（・K−0．・，一・。・0、02・・0．00005・2） 9 巨
TRuE vAu」E NK＝2 83 NK＝1 NK＝O o㎡（
Hエ Ol32968 0、32968 0，2ア1q o．2708 n．ら
H2 022099 022099 0・210ア n2106
Hユ 0．1捌q O・ユ正131q O・1203 0・U9ア ゜ 2
H4 0．09929 0．09929 〔｝．0812 010808 0
H5 ｛〕．06656 り．0〔…656 0、〔］1｛68 0．046叫
H6 010q562 0．04552 0．〔〕29q O、0293 ．う，コ
1；；：；鑑1：器器；；：器；；；：1；；； Fi島5−・dR・・u1・・
Hg O』ユ54q o・013qq−0・0058−9』052
il iiiil°11iiiii』iiii＝：元二内灘三：㌻鴛：時に
c25・10−5 @5・10−5 擁している場合について錨を行なうた，対象のシス
・・… （3x・r12） テムはFig．6 、示すものでディ〃ル計轍・よって
模擬を行なった，サンプル周期を2秒に選び，超低周波
4−2 内部雑音のあるシステム 雑音ぬをO．005力丁とし重畳させ，データ個数，パフ
これまでに述べた方法はシステムに重畳している不規 メータ個数をそれぞれ∫Vp＝ユ0，1Vd＝400，1V・＝LIV・則雑音が白色雑音になるように推定しているので，内部 ＝2と選案し，推定を行なった。Figコが対象のシズテムの
継の存鱈るフ仁ド…。ク制御系においてもその ステ。プ酪の齪結果である．この齢雑剖寺性の推
内部のサプ学ステムの酬酷不雛定できると考え 定繍をイン・・ルス応馴の形式で獅したのが・Fig・
られる．これを示すため1・Fig．41・示す江テムをアナ 8である．これらの結果から対象のシステムの動雛と
ログ計算機によって模擬し，そのサブ・システムの入出 雑音特性を合わせて同時に推定できることがわかる・恥…。1，、 雅 H、 参考文献
一コー一一
〔U凸脅■1） Sふ25S i
l｝鈴木 肝：，古谷忠邊，藤井克彦：動特性推定に川いるプロ 怜
セスの操業データの処理方法． 1 i 測自亜》制御学会描文集．4．
出IT［Ho：s［
@ 1 1 ． ◆ ° 2、 175／1B1 （1967）．
ヘヰ
マロけn・引 f2s 1「 u 5s卓1 2｝W． D、 T． Davi5，」． L．1〕ouce：On・line Sy5tem ldenti伍ca．
，o
1．5
口
1．o
o．〜
tbn irl the Presence of Drift， PreprinLs o臼h〔｝IFAC Symp．，
Fi酔6 S㎞皿bted Proce88 3−12， Prague Czecho（1967），
3｝P．N． Nikiforuk， M． M． Gupta， R、 Ho「fman：Ident甫cation of Linear Discrete Systcm5 in the Presence of DnfL u5ing ． 一 内 ＾ ㌔ 「 Pseudo Random Sequences， Prep加t50f the 2nd Prague ． ● IFAC Symp．、7−4， Prague Czecho U970｝．
4｝D．Q． Mayne：Parameter Estimation， AutomaUca，3，
245／255（1966）．
1 5）S、G．T囲fe5tas：SDme compu民F且lded Estimation in S1㏄ha§ti〔 Corltrol Systems Iderltification， Int、」． Con甘ol，
電 ［｝叫唱r コv肱頂
12，3，385／399｛］970）．
1o 力（、［d 6）H・Akaike：Annals of the lnstiしute of Stat已tical Mathe・
Fi島7−・・dR…1t・・fF貼6 ，鷲驚、2ぽ迦、．S、㏄h田、，C．＿＿
98／102，Acadernic Press（1970）．
19 。 白 苫《玉 ） O o
Fig．8 Noi賠Dynamic8
5．むすび
プロセスの動特性の際に，どうしても超低周波成分の 混入．が避けられないことから、その成分を多項式近似で あちわし，不規則雑音のパラメータとともに推定する方 法について報告した。このなかで超低周波成分の存在す る場合やシステムの内部に雑音が存在する場合の動特性 推定が可能なことを実験例によって示した。ここで述べ た方法は偲号に含まれている雑音成分の特性をそれぞれ 2つの成分に分けて除去するようにしている。すなわち 多項式近似と雑音過程の自己回帰モデルによる表現であ る。したがって計算のうえで隔通性に富む形式であり，
残差が白色となるようにパラメータ個数を決定してやる だけでよいという特徴がある。
最後に，九州工業大学 米澤洋教授ならびに大阪大 学 藤井克彦教授に謝意を表するとともに，本研究が昭 和46年度文部省科学研究費｛奨励A｝によって行なわれ たことを付記する。
