利得を指定して IIR フィルタ ( 感度関数も含む ) を設計する方式の SCILAB 用プログラムの紹介 木坂正志 MK サーボ開発 2010/3/16 NSS2

利得を指定して

IIRフィルタ（感度関数も含む）を設計する方

式の

SCILAB用プログラムの紹介

木坂正志 MKサーボ開発

最近の発表との関係 FIRによる位相設計 （_{Filtered-x LMS)} IIRによる位相設計 最小位相関数による_{利得の設計} 最小位相関数を使ったサーボ設計（SISO) 安定多項式の実数表現 適応制御 木坂 正志_{:”必ず安定となるIIRフィルタとそのハードディスク装置の外乱に対する適応抑制への応用”、第} ６回制御部門大会資料、_{Vol.2、469/474(2006)} 木坂 正志_{:”最小位相関数による適応制御とその位置決め制御への応用”、平成２０年電気学会産業応} 用部門大会 木坂 正志_{:”GKYP 補題を使った感度関数設計”、第8 回計測自動制御学会制御部門大会、2008.3}

SCILAB 5.2.0 アプリケーション モジュール管理（ATOMS) Minphase インストールすると Scilab ヘルプの中に minphase のヘルプが入る。 Scilabデモのなかにminphaseが入る。 log_fir(): 最小位相関数の周波数特性を求める。 piir(): ナイキスト周波数で進まない位相特性から安定なIIRフィルタを得る。 ソースコードは_{ATOMSにloginするととれる。}

アプリケーションを選ぶ

P 最小位相関数とサーボの関係 C D -+ W D N L i i q L i i P P p z z p P = -=

Õ

Õ

= -= 1 1 a b D N d L L i i d q L q L i i D D d z z d D D N = -=

Õ

Õ

+ + = + -+ -= 1 1 a b Dは積分器のように決まった関数。 Cを求める。 PDの分子の次数は、分母より低いとする。

感度関数Sを求める D N d L D D d L D D N N D D D D X X z D P X z D P C D pdP D P D P PDC S D D + + º º + = + = 1 1 Xを求めそこからCを解く。 分母のzはXで消去できるので安定な関数ならなんでも良い。 N N D D D d L N N D D d L N N X X D P X z D pdP X X D P z D pdP C = D - = D -+ + ₁ 1 S、Cは安定なのでXは分母、分子安定でなければならない。 （任意の対象で安定な_{Cが得られるわけではない 参考文献）} Xは最小位相関数。 ただし、以下の二つの制約条件が付く。

( ) ( ) ( )

=0 -+ k k D k D d L k P D X D b b b b _³1 k b

( ) ( )

_{k N k} = 0 N D P b b 等式条件 ¥ ® ® + = + z C D pdP D P z X N N D D d L _D 1 _{+ ® ® z ®¥} D P D P z D P D D N N d L D D D 1 0 X f z z D P X f D P D D D D L d d L d L D D d L D D + + + + =

Xは最小位相関数なので、単位円外にゼロ、極はない。 .. log 3 3 2 2 1 1 0 + + + + º _{a a z}- _{a z}- _{a z} -X .. 3 cos 2 cos cos log ) real(log log log = ₊ + = ₊ +a₀ +a₁ w +a₂ w +a₃ w + z D P X z D P S D D L d D D d L D D .. 3 sin 2 sin sin ) imag(log = Ð + ₁ + ₂ + ₃ + + Ð = Ð _{+ +} a w a w a w z D P X z D P S D D L d D D d L D D

( )

log cos cos2 cos3 ..

log w - ₊ ³ a₀ +a₁ w +a₂ w +a₃ w + z D P S D d L D D U S_U(ω)：感度関数の上限とする

感度関数の利得は logS = log

(

S ejÐS

)

より_{logをとった実部に対応するので}

.. 3 cos 2 cos cos ) real(log log X = X = a₀ +a₁ w +a₂ w +a₃ w + Log(利得) 周波数 π 周波数軸を時間軸とみるとフーリエ級数表現といえる。 利得を表す任意の連続関数を近似できるといえる。 （位相だけ与えた場合もそれを近似可能といえる。） 感度関数の制約である面積公式は、直流成分がゼロとみなせる。 _a0=0 （不安定極がある場合も） ¥ ® ® + + + º _{a z}- _{a z}- _{a z}- _z X exp( 3 ..) 1 3 2 2 1 1 このとき

( )

log cos cos2 cos3 .. log w - ₊ ³ a₁ w +a₂ w +a₃ w + z D P S D d L D D U

( ) ( ) ( )

= 0 -+ k k D k D d L k P D X D b b b b PN

( ) ( )

bk DN bk =0 1 ³ k b d £ - S 1 ロバスト条件がある場合 _S real imaginary 1 サーボ設計のまとめ（最小位相関数の特別な場合といえる） N N D D D d L N N D D d L N N X X D P X z D pdP X X D P z D pdP C = D - = D -+ + ₁ 1 ...) exp( 3 3 2 2 1 1 + + + = _{a z}- _{a z}- _{a z} X S Ð S S_U(ω)：感度関数の上限の仕様 コントローラ Xを有理多項式であらわせるなら。

最小位相関数_{Fの利得の仕様から、以下の式で係数を得る。}

L

U F F a a a a F

F log real(log ) cos cos2 cos3 .. log

log ³ = = ₀ + ₁ w + ₂ w + ₃ w + ³ 位相を含めた周波数特性を .. log 3 3 2 2 1 1 0 + + + + º _{a a z}- _{a z}- _{a z} -F で得る。 係数を得る方法は、Remez法が有名。 サーボを考えると、等式条件があるので、他の方法が望ましい。 本題：最小位相関数の伝達関数を求める。（一般の話）

.. ) (z º a₀ +a₁z-1+a₂z-2 +a₃z-3 + G )) ( ) ( ( 2 1 log _{F = G z +G}* _z GKYP補題で扱える。 サイズの大きなマトリクスをつかうので計算たいへん。 木坂 正志_{:”GKYP 補題を使った感度関数設計”、第8 回計測自動制御学会制御部門大会、2008.3}

) 1 log( log .. 3 cos 2 cos cos ) 1 log(

logg_i + +err ³ a₀ +a₁ w_i +a₂ w_i +a₃ w_i + ³ g_i - +err

log_fir()でやってること

( )

/(1 ) ) 1 ( err F g err g_i + ³ w_i ³ _i + いくつかの周波数で以下を指定 g_i目標 err許容誤差 LMIで以下の線形不等式を解く この解は無限にあるので、微分と差分の差を表す以下の評価関数Eを最小 にするようにしている。これで、点間をなだらかに結ぶようにしている。 .. ) ( ) ( logF z º G z = a₀ +a₁z-1 +a₂z-2 +a₃z-3 + T e z G z T e G e G D T j i T j T j i i i w w w w w ¶ ¶ ¶ ¶ -D -º ( +D ) ( ) ( ) ) 1 (

å

= = 100 0 2 i i D E

aa = log_fir(t0,forder,file,err) 関数の引数 t0: sample time forder: cosine級数の次数 file: 仕様を書いた外部ファイル名 err: 許容誤差 file f1 g1 p1 f2 g2 p2 .. 周波数f_i(Hz)で、利得がg_i(dBではない）と指定している。 p_iは使っていない。 exp(G)を301個の周波数で計算したapgain.txtを出力する。 LMIがとまる場合は次数を上げるか、errを大きくする。 .. ) ( 3 3 2 2 1 1 0 + + + + º _{a a z}- _{a z}- _{a z} -z G

piir()でやってること 木坂 正志_{:”IIRフィルタによる同定問題とそのサーボ設計への応用”、第36回制御理論シンポジウム、} 2007.9

å

å

= = º = _n l l l n i i i z b z a z g z f z p 0 0 ) ( ) ( ) (

å

å

å

= + -= -= º Ð º Ð = Ð = Ð = Ð n t t n t n l l l n i i iz b z x z XZ a z g z f z g z f z p 2 0 0 0 ) / 1 ( ) ( ) ( ) ( ) (

[

x x x n x n

]

X º _{0 1} .. _{2 -1 2}

[

_{n n n n}

]

t z z z z Z _º - - +1 _.. -1 ナイキスト周波数で進まない位相が与えられたとき、それを満たす安定な IIRフィルタを得る。（PRGはナイキスト周波数の条件をCheckしてはいない。） 等価な位相を持つ関数を定義。

( ) ( )

(

)

(

)

(

)

(

)

C S n n n n jL L n n n n j n n n n z z z z Z º + ú ú ú ú ú ú û ù ê ê ê ê ê ê ë é -+ ú ú ú ú ú ú û ù ê ê ê ê ê ê ë é -= ú ú ú ú ú ú û ù ê ê ê ê ê ê ë é º -w w w w w w w w sin 1 sin .. 1 sin sin cos 1 cos .. 1 cos cos .. 1 1

( )

z f z g z XZ h º ( ) (1/ ) = 仕様。各周波数で定義されたθに誤差α>0で近似する。

( )

w -q £a <p /2 Ð_{h e}j ベクトル L º

[

LC LS

]

との内積で仕様を表す。 a q q cos sin cos _{t t} X XLL XL _ú > û ù ê ë é

(

A B

)

B A AB = cos Ð -Ð

0 cos 0 0 sin 2 2 > ú û ù ê ë é -t t t_{L X} U XLU a a ú û ù ê ë é -º q q q q cos sin sin cos U ú û ù ê ë é -ú û ù ê ë é -ú û ù ê ë é -= ú û ù ê ë é - q q q q a a q q q q a a cos sin sin cos cos 0 0 sin cos sin sin cos cos 0 0 sin 2 2 2 2 t U U ú û ù ê ë é -ú û ù ê ë é -= q q q q a q a q a q a q cos sin sin cos cos cos sin sin cos sin sin cos 2 2 2 2 ú û ù ê ë é -+ + -= a q a q a q q a q q a q q a q q a q a q 2 2 2 2 2 2 2 2 2 2 2 2 cos cos sin sin cos sin cos sin sin cos cos sin cos sin sin cos cos sin sin cos

(

)

(

)

ú û ù ê ë é -= a q a q q q q q a q a q 2 2 2 2 2 2 2 2 cos cos cos 1 sin sin cos sin cos cos sin cos 1 cos

[

]

- ×I ú û ù ê ë é = ú û ù ê ë é -= q q a q q a q q q q q a q ₂ 2 2 2 2 cos sin cos sin cos cos sin sin cos sin cos cos cos

(

cos

)

0 sin cos 2 2 > -÷÷ ø ö çç è æ ú û ù ê ë é a q q _{t t} X XLL XL

[

u v

]

XLU º

(

sin cos

)(

sin cos

)

0 cos sin2 2 2 2 _{a -v a = u a -v a u a +v a >} u 0 cos sin cos > > ú û ù ê ë é = a q q _{t t} X XLL XL u 0 cos sina +v a > u 0 cos sina -v a > u Xに対する線形不等式となる。 LMIで解ける。 u v よりuは正なので

[

]

0 cos 0 0 sin 2 2 > ú û ù ê ë é ú û ù ê ë é - v u v u a a

( )

å

_Õ

(

) (

_Õ

)

_Õ

(

)

_Õ

= = = = = _÷ ø ö ç è æ -= -= = º n k n k k n i i k n i i n n t t t z z z z z z z x z g z f z h 1 1 1 1 2 0 ₁ ) / 1 ( ) ( l r l r gが安定のためにはh(z)の分子多項式で単位円外の根ρがn個以上必要 h(z)の分子は2n次、分母はn次の安定根なので、ωが0から2πまで増加したとき、 h(z)の軌跡が原点の周りを時計回りに0回以上回らなければならない。 実関数なので、これは、ナイキスト周波数で位相が進まないことを意味する。

Õ

Õ

= = -º _n i k n i i z z z p 0 0 1 ) ( r l 係数xを得た後、h(z)よりf(z)、g(z)を得る。 実際の計算では、分母多項式の 最高次数の係数が１となるようにしている。

安定根候補がn個以上ある場合、n個をとりあえず選び、残りを分子に使うとする。 ) / 1 ( ) ( ) / 1 ( ) ( ) ( ) ( ) ( ) ( ) ( z U z z U z U z z G z F z G z U z F z P m m m m m n m n n m m n- ₌ -= ) / 1 ( ) ( ) ( ) ( 0 2 0 z U z u z z u u z z G z F m m m i i m l k m l k k k l l m n m n ÷ ø ö ç è æ _{+ +} =

åå

= =

å

= -+ -) / 1 ( cos 2 ) ( ) ( 0 2 0 z U z u k u u z z G z F m m m i i m l k m l k l l m n m n ÷ ø ö ç è æ ₊ = -

åå

= = +

å

= w ) (cos ) ( ) ( 0 w m m i m i i m n m n _O z u z z G z F

å

= + -Þ

å

= º m i i i m z u z U 0 ) ( 絶対値が１以上の 根を集めた関数 分母は安定 O_mは正の実関数なら位相条件は満たされる。

å

å

= = + -= = m i i m m i m i i m n m n i a O z u z z G z F z P 0 0 cos ) (cos ) ( ) ( ) ( w w

å

å

= + -= -- ÷÷ø ö çç è æ + = _m i m i i m i i i i m n m n z u z z a z z G z F z P 0 0 2 ) ( ) ( ) ( 利得の仕様から、係数_{aを得る。} 得られる関数 O_mは負にならない関数から選ぶべきだが、PRGでは考慮していない。

fp=piir(t0,nf,file1,tol,[qx,[nfx]]) t0: sample time nf: 最小位相関数のときのIIRフィルタの次数 file1: 仕様を書いた外部ファイル名(“apgain.txt”) tol: 許容誤差 (degree) qx: 結果の関数の分母と分子の次数差 nfx: fp=x*nfx nfxをz^-nfとするとFIRとなる。 fp: 結果のIIRフィルタ 関数の引数 file1 f1 g1 p1 f2 g2 p2 .. 利得_giは、_{Kを決めるために必要} LMIが止まる場合は、次数を上げる。その後で止まる場合は次数を偶数とする。 最小位相関数でない周波数特性を与えた場合は、_{nfにマイナス値を} 使うと、O_m(cosω)の部分を使い利得を合わせる。 n+m次のIIRフィルタとなる。

Õ

Õ

= = -º _n i k n i i z z K z p 0 0 / 1 ) ( r l

Demo

WindowsでDemoを動かすためには、Current directoryを変える必要がある。 （メニューバー）より変えられる。

Known function case

2次関数の近似 Approximation of exp(1/z)

exp(1/z)を2次関数で近似

Discretization from a continuous transfer function F(s) = H(s)(1-1/z)/s

D(z)=(F(ω-Nω_s)+ F(ω-(N-1)ω_s)+.. F(ω+(N-1)ω_s)+ F(ω+Nω_s))/t_s Function with arbitrary gain

///// approximation of exp(1/z) z=%z; ts = 1; // sample time fz=1/z; //target function fs0=0.01*1/ts/2; // start freq. fsn=0.49*1/ts/2; // stop freq. ff = fs0+(fsn-fs0)*(0:20)/20; zf = exp(%i*ff'*2*%pi*ts); qf = exp(horner(fz,zf));

at =[ff',abs(qf),imag(log(qf))/%pi*180]; // prepare data in the file if MSDOS then unix('del foo');

else unix('rm -f foo'); end write('foo',at)

px = piir(ts,2,'foo',0.01); // The order of px is 2. Phase error is 0.01 // compare the result

clf() qa=horner(px,zf); bode(ff,[qf.';qa.']) // px is almost same as fz ２次IIRで近似 複素関数は正則領域内の一部の線で 一致すると全正則領域で一致する。

Function with arbitrary gain z=%z; ts = 1; // sample time at = [ // gain spec 0.01,0.01,0; // gain is 0.01 at 0.01Hz 0.1,1 ,0; // 0.1 at 1 0.2,2 ,0; 0.3,1.0 ,0; 0.4,0.01 ,0; 0.49,0.01,0];

if MSDOS then unix('del foo'); else unix('rm -f foo'); end write('foo',at)

q1=log_fir(ts,80,'foo',.03);

分母にz-1を持っていると分かっている場合

// the target function is stable except a predetermined factor // but not the minimum phase function

z=%z; ts = 1; // sample time fz=(z+4)*(z-0.2)/(z-1)/(z+0.1); //target function fs0=0.01*1/ts/2; // start freq. fsn=0.49*1/ts/2; // stop freq. ff = fs0+(fsn-fs0)*(0:20)/20; zf = exp(%i*ff'*2*%pi*ts); qf = horner(fz,zf);

at =[ff',abs(qf),imag(log(qf))/%pi*180]; // prepare data in the file if MSDOS then unix('del foo');

else unix('rm -f foo'); end write('foo',at)

px = piir(ts,-2,'foo',0.001,0,1/(z-1)); // the function has a predetermined factor z-1. 指定する次数を上げると、安定根のペアが分子分母に加わる。

まとめ SCILAB 5.2.0 アプリケーション モジュール管理（ATOMS) Minphase 利得の仕様から周波数特性を計算する。 aa = log_fir(t0,forder,file,err) ナイキスト周波数で進まない位相を与えたとき それを実現する安定な_{IIRを得る。} fp=piir(t0,nf,file1,tol,[qx,[nfx]]) サーボをサポートする_{Versionに更新予定}