:Amazonリンク:[システム同定の基礎|https://amzn.to/2OSoV79] !!! 2018年5月1日 Scilabでシステム同定をやりたいのだが、、どうすれば良い？ 検索するも、MATLABだったらすぐできそう、という情報しか得られない。 MATLAB欲しいなあ、と思いつつも、そんなお金はないので、Scilabで頑張る。 !!! 2018年5月2日 ところで、突然「システム同定」とか言い出しだが、定義は、、 :システム同定:計測データから動的モデルを構築するための数学的ツールやアルゴリズムを指す用語。 (Wikipedia「[システム同定|https://ja.wikipedia.org/wiki/%E3%82%B7%E3%82%B9%E3%83%86%E3%83%A0%E5%90%8C%E5%AE%9A]」より) システム同定で制御対象のモデルとして構築できれば、実物を用いる代わりにモデルを用いて制御の検証ができる。 Scilabでシステム同定をやるには、「CACSD」を用いる。「CACSD」とは「Computer-Aided Control System Design」の略らしい。 CACSDライブラリの中で、システム同定に関係ありそうなのは、、 [Scilab Online Help|https://help.scilab.org/]で調べた。 「CACSD」の中のさらに「identification」に分類されている関数が、システム同定に関するもの。 その中でさらに「同定」と説明があるのは、、 *armax: armax 同定 *armax1: armax 同定 *findABCD: 離散時間システムの部分空間の同定 *findAC: 離散時間システム部分空間の同定 *time_id: SISO 最小二乗同定 どれが良いだろうか、、？ !!! 2018年5月3日 どの「同定」を試してみるか、、ヘルプを見て「time_id」の記載が簡単に見えたので、試してみる。 「time_id」の説明のトップは「SISO最小二乗同定」とある。 「SISO」は「Single In Single Out」の略だろう。いわゆる「単入力単出力」なので、簡単そうだ。 ヘルプには簡単な例もある。これを動かしてみよう。 z = poly(0, 'z'); h = (1 - 2 * z)/(z^2 - 0.5 * z + 5) rep = [0;ldiv(h('num'), h('den'), 20)]; //インパルス応答 H = time_id(2, 'impuls', rep) // flts および uを用いる同じ例 u = zeros(1, 20);u(1) = 1; rep = flts(u, tf2ss(h)); //インパルス応答 H = time_id(2, u, rep) // step response u = ones(1, 20); rep = flts(u, tf2ss(h)); //ステップ応答. H = time_id(2, 'step', rep) H = time_id(3, u, rep) // uを入力として指定, 高次過ぎる系を指定 まず、 z = poly(0, 'z'); から。Polyは多項式を定義する関数で、zを用いた多項式を定義するための準備作業のようだ。 実行しても、反応はない。 h = (1 - 2 * z)/(z^2 - 0.5 * z + 5) h = 1 - 2z ------------- 2 5 - 0.5z + z 実行すると、Scilabのコンソールに出力された。 どうやら、zを用いた式を作り、それをhとしたようだ。 rep = [0;ldiv(h('num'), h('den'), 20)]; //インパルス応答 この中の、ldivは、「多項式行列の長除算」を行う、らしい。 ""x=ldiv(n,d,k) は, nのdによる長除算の 最初のk個の係数を出力します. すなわち,有理行列[nij(z)/dij(z)]の無限大近傍での テイラー展開を出力します. さっぱり分からない、、なお、「テイラー展開」とは、 ""数学において、テイラー級数 (Taylor series) は関数のある一点での導関数たちの値から計算される項の無限和として関数を表したものである。そのような級数を得ることをテイラー展開という。 コンソールにrepと入れてみると、テイラー展開の係数が保存されている。 rep = 0. -2. 0. 10. 5. -47.5 -48.75 213.125 350.3125 -890.46875 -2196.7969 3353.9453 12660.957 -10439.248 -68524.409 17934.036 351589.06 86124.354 -1714883.1 -1288063.3 7930384. まとめると、「rep = [0;ldiv(h('num'), h('den'), 20)]; //インパルス応答」は、、 * "[", "]" は、行列・ベクトルの定義 (ここではベクトル)。";" で区切ると列ベクトルになる * h('num') は、hの分子 (numerator)。h(2), h.num, numer(h) とも書ける * h('den) は、hの分母 (denominator)。h(3), h.den, denom(h) とも書ける * 20 とあるので、最初の20個の係数を出力。"0" と合わせてベクトル内の値は21個 コメントに「インパルス応答」とあるので、repはインパルス応答のようだ。 H = time_id(2, 'impuls', rep) H = 1 - 2z ------------- 2 5 - 0.5z + z time_id登場。引数のうち、 * 2 は出力する伝達関数の次元 * 'impuls' は指定したベクトルがインパルス応答であることを指定 * rep は先ほどのベクトル。インパルス応答 time_id() で計算された H は、、最初の h を同じ。 つまり、このサンプルは、 * 伝達関数 h を作成 * h からインパルス応答 rep を生成 * rep から伝達関数 H を同定 ということ？、、、すごくない？ // flts および uを用いる同じ例 次は、同じことを flts と u でやる、と、、 u = zeros(1, 20);u(1) = 1; zeros(1, 20) で、値が0の行ベクトル (20列) を作り、 u(1) = 1 で、最初だけ1、残りは0の行ベクトルになった。 u = column 1 to 11 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. column 12 to 20 0. 0. 0. 0. 0. 0. 0. 0. 0. 次の、 rep = flts(u, tf2ss(h)); //インパルス応答 * tf2ss() は、伝達関数を状態空間表現に変換する * flts() は、時間応答を生成する つまり、flts() は、状態空間表現にした h に対し、インパルス入力信号 u を適用し、時間応答を得る、、 ところで、これで得られた rep の値は、 rep = column 1 to 8 0. -2. 0. 10. 5. -47.5 -48.75 213.125 column 9 to 12 350.3125 -890.46875 -2196.7969 3353.9453 column 13 to 16 12660.957 -10439.248 -68524.409 17934.036 column 17 to 20 351589.06 86124.354 -1714883.1 -1288063.3 と、先ほどの ldiv() と異なり、最初から 0 が入っている。 H = time_id(2, u, rep) 先ほどと同じ、time_id() が登場。 今度は2番目の引数が u となっている。システムへの入力ベクトルをそのまま指定した模様。 H = 1 - 2z ------------- 2 5 - 0.5z + z なんと、先ほどと同じ結果が、、 まとめると、 * 伝達関数 h を状態空間表現に変換 (tf2ss() ) * インパルス入力 u を h に適用し、インパルス応答を生成 (flts() ) * u および rep から伝達関数を同定 (H) ただしよく考えると、1つ目と同じくインパルス応答から伝達関数を同定している。 ためしに2番目の引数 u を、1つ目と同じ 'impuls' にしてみても結果は同じだった。 // step response 今度はステップ応答らしい。 // u = ones(1, 20); ones() で値が1の行ベクトル (20列) を作る。コンソールで確認すると、確かに 1 だ。 u = column 1 to 11 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. column 12 to 20 1. 1. 1. 1. 1. 1. 1. 1. 1. rep = flts(u, tf2ss(h)); //ステップ応答. 先ほども出た。tf2ss() で状態空間表現にし、flts() で時間応答を生成する。 H = time_id(2, 'step', rep) H = 1 - 2z ------------- 2 5 - 0.5z + z ステップ応答でも、元の伝達関数が同定できた。すごい。 なお、次の例は実際より高次元 (3) を指定したとき。 この場合、正しい結果を返さない。でも分子分母を z で割ると、ほぼ同じになる。 H = time_id(3, u, rep) // uを入力として指定, 高次過ぎる系を指定 H = 2 -4.226D-10 + z - 2z -------------------- 2 3 5z - 0.5z + z :Amazonリンク:[システム同定 (計測・制御テクノロジーシリーズ)|https://amzn.to/2KAfypi] !!! 2018年5月4日 さて、昨日まで time_id() を例を使って体験してみた。 ただ、ちょっと分からないことがある。 time_id() は引数としてn, u, y をとる。 その中の応答ベクトル y だが、最初の例は ldiv() から、それ以降は flts() を用いて応答ベクトルを生成していた。 しかし、それら2つは、 * ldiv: 多項式行列の長除算。長除算の係数を出力 * flts: 時間応答 とある。 !!! 2018年5月5日 昨日のつづき。 time_id() の引数 y の値は、Scilabのヘルプでは「応答ベクトル」とされている。 ところでいま気付いたが、Scilabのローカルのヘルプ「ヘルプブラウザ」はオンラインヘルプと同じ内容だった。 さらに、ヘルプ上のサンプルコードをコンソールで実行したり、SciNotes に転記したり、色々便利。 で、応答ベクトルだが、私は「システムへの信号入力に対する応答が時系列で並んでいるもの」という風に考えていた。 flts() の説明は「時間応答」で「応答ベクトル」っぽいが、ldiv() の説明は「多項式行列の長除算」で、テイラー展開の係数出力とある。ちょっと応答ベクトルっぽくない。グラフにプロットするとこんな感じ。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image ScilabLdivPlot180505.png}} しかし、flts() の出力もこんな感じ。きっと、私の「応答ベクトル」の理解がおかしいのだろう、、、 ところで、「長除算」って何だ？と思って調べたが、よく分からない。 しかし、オンラインヘルプの言語を英語にしたところ、「long division」とあり、要は「筆算による除算」だった。 !!! 2018年5月6日 昨日のつづき。 多項式を筆算で除算、、うーんよく分からない。 ldiv(1, 7, 5) とやってみる。1 ÷ 7 だ。 ans = 0.1428571 0. 0. 0. 0. うん、ただの1 ÷ 7 だ。これが 1 ÷ 7 のテイラー展開の結果なのか、、？ この問題はこれ以上引っ張っても仕方ない。 「応答ベクトル」とは、「テイラー展開の係数のベクトル」だとしておこう。 さて、私が結局やりたいことは、 * ステップ入力やインパルス入力に対する時系列の応答データ列から、そのシステムの同定を行う だ。 time_id() では何となく出来なさそうな気がしてきた。 でもとりあえず、やってみよう。 :Amazonリンク:[システム同定―部分空間法からのアプローチ|https://amzn.to/2vtIhYe] !!! 2018年5月7日 昨日のつづき。時系列の入出力データを用意しよう。 何にしようか？ RC回路の電圧、これは1次遅れとなるため、これを使ってみる。 !!! 2018年5月8日 RC回路の電圧だが、[LTspice|http://www.analog.com/jp/design-center/design-tools-and-calculators/ltspice-simulator.html] で作り、グラフにプロットした。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image RcCircuit180508.png}} (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image RcGraph180508.png}} グラフのデータをテキストファイルにエクスポートした。 [RcData180508.txt|http://moi2.sakura.ne.jp/fswiki/wiki.cgi?page=Scilab%A4%F2%BB%C8%A4%C3%A4%C6%A4%DF%A4%EB&file=RcData180508%2Etxt&action=ATTACH] これを、Scilabに読み込ませる。 rc = fscanfMat("RcData180508.txt") !!! 2018年5月9日 [「Maxima」と「Scilab」で学ぶ古典制御 (I・O BOOKS)|https://amzn.to/2I06eyX]という本の、時定数と定常ゲインを算出するコードを参考に書いてみる。 !!! 2018年5月13日 続き。 function y = resp1(t, param) y = param(1) * (1 - exp(-t / param(2))); endfunction function er = error_resp1(x, t, data) er = sum((resp1(t, x) - data).^2); endfunction rc = fscanfMat("RcData180508.txt"); t = rc(:, 1); v = rc(:, 2); param0 = [2, 1]; [fer, popt, gropt] = leastsq(list(error_resp1, t, v), param0); 計算結果は popt に格納されているが、 popt = 5.0000609 0.1000043 定常ゲインが5、時定数は約0.1と出た。 入力したファイルは、定常値5V、時定数0.1 (100kΩ × 1uF = 0.1) なので合っている。 なお、電源電圧が5Vなので、本当は定常ゲインは1だ。 この方法は、最小二乗法関数 leastsq() を用いて定常ゲイン、時定数を同定している。 この方法はモデルが既に分かっているときは良いが、分からないときはどうしようもなさそうだ。 では、armax(), armax1() はどうだろうか？ armax() はn次元、armax1() は1次元のARXプロセスの係数を同定する、とある。 なら、armax1() の方が簡単だろう。サンプルコードもこちらのほうが短い。 まず、サンプルコードを動かしてみる。 a = [1, -2.851, 2.717, -0.865]; b = [0, 1, 1, 1]; d = [1, 0.7, 0.2]; ar = armac(a, b, d, 1, 1, 1); disp(_("Simulation of an ARMAX process:")); disp(ar); n = 300; u = -prbs_a(n, 1, int([2.5,5,10,17.5,20,22,27,35]*100/12)); zd = narsimul(ar, u); // armax1を使用: 有色ノイズ同定 // 同じ例を以下のように試すことができます [arc1, resid] = armax1(3,3,2,zd(1:n),u,1); disp(arc1); * 1行から3行は、ベクトル a, b, d を作っている。 * 4行は、、armac() 新しい関数だ。armax処理のためのデータ構造を作るらしい * 6行は、armax処理のデータ構造 ar を表示している * 8行。prbs_a() 。これも新しい関数。疑似バイナリ乱数列を作る。u には 1 と -1 がランダムに並んでいる * 9行。narsimul() 。モデル ar と 入力信号 u で、armaxシミュレーションを行う * 12行。armax1() 関数でシステム同定を行う。引数は先頭の 3, 3, 2 は自己回帰の次数、9行で生成した出力信号 zd、8行で生成した入力信号 u、最後は、1 のときに b0 を 0 と仮定する、らしい さて、元のモデル ar は、、 --> disp(ar); A(z^-1)y=B(z^-1)u + D(z^-1) e(t) A(x) = 2 3 1 -2.851x +2.717x -0.865x B(x) = 2 3 x +x +x D(x) 2 1 +0.7x +0.2x e(t)=Sig*w(t); w(t) 1-dim white noise Sig= | 1 | そして、同定したモデル arc1 は、、 --> disp(arc1); A(z^-1)y=B(z^-1)u + D(z^-1) e(t) A(x) = 2 3 1 -2.8567167x +2.7284292x -0.8707818x B(x) = 2 3 1.1990058x +1.2509158x +0.5501596x D(x) 2 1 +0.3460801x +0.0466038x e(t)=Sig*w(t); w(t) 1-dim white noise Sig= | 1.1431121 | A(x) は割りと近い。B(x)、D(x) は割りと違うとこもある。 とは言え、元々考えていた、時系列の入出力データからシステムの同定を行う方法だ。 :Amazonリンク:[MATLABによる制御のためのシステム同定|https://amzn.to/2Oh9gNW] !!! 2018年5月14日 5月2日に調べた、Scilab で利用できるシステム同定の関数を改めて列挙する。 *armax: armax 同定 *armax1: armax 同定 *findABCD: 離散時間システムの部分空間の同定 *findAC: 離散時間システム部分空間の同定 *time_id: SISO 最小二乗同定 time_id() は時系列データをそのまま使えなさそうだった。 armax1() は時系列データを元に同定できた。 armax() は、armax1() をn次元に拡張した方法らしいので、使い方はだいたい同じっぽい。 あとは、findAC() と findABCD() だ。 !!! 2018年5月15日 続き。 findAC() と findABCD() は、どちらも「離散時間システムの部分空間の同定」となっている。 ほぼ同じっぽい。 早速、サンプルコードを動かしてみる。findAC() から。 //指定した線形システムからデータを生成 A = [ 0.5, 0.1, -0.1, 0.2; 0.1, 0, -0.1, -0.1; -0.4, -0.6, -0.7, -0.1; 0.8, 0, -0.6, -0.6]; B = [0.8;0.1;1;-1]; C = [1 2 -1 0]; SYS = syslin(0.1, A, B, C); nsmp = 100; U = prbs_a(nsmp, nsmp / 5); Y = (flts(U, SYS) + 0.3 * rand(1, nsmp, 'normal')); // Rを計算 S = 15;L = 1; [R, N, SVAL] = findR(S, Y', U'); N = 3; METH = 3;TOL = -1; [A, C] = findAC(S, N, L, R, METH, TOL); * 2〜7行: A、B、C は、syslin() への引数 * 8行: syslin() は、線形システムの定義。サンプル時間0.1秒の状態空間を作る * 10行: prbs_a() は前にも出てきた、疑似バイナリ乱数列生成関数 * 11行: flts() 、これも前に出た。時間応答を生成する関数。状態空間 SYS に対し、入力信号 U を適用。rand() は乱数を作る * 14行: findR()、これは、「線形時不変システムの行列を指定する際のプリプロセッサ」とある。よく分からない。Y と U についてるシングルクォートは、行列の転置指定 * 17行: findAC() うん、分からない。 !!! 2018年5月16日 findAC(), findABCD() の中にある、A,B,C,D は、[状態空間表現|https://ja.wikipedia.org/wiki/%E7%8A%B6%E6%85%8B%E7%A9%BA%E9%96%93_(%E5%88%B6%E5%BE%A1%E7%90%86%E8%AB%96)] のために用いられる。 dx(t)/dt = Ax(t) + Bu(t) y(t) = Cx(t) + Du(t) x が状態ベクトル、y が出力ベクトル、u が入力ベクトル、A が状態行列、B が入力行列、C が出力行列、D が直達行列。 :Amazonリンク:[カルマンフィルタとシステムの同定 動的逆問題へのアプローチ|https://amzn.to/2OPZYJy] !!! 2018年5月17日 続き。 findAC() のサンプルコードをもう一度見てみる。 //指定した線形システムからデータを生成 A = [ 0.5, 0.1, -0.1, 0.2; 0.1, 0, -0.1, -0.1; -0.4, -0.6, -0.7, -0.1; 0.8, 0, -0.6, -0.6]; B = [0.8;0.1;1;-1]; C = [1 2 -1 0]; SYS = syslin(0.1, A, B, C); nsmp = 100; U = prbs_a(nsmp, nsmp / 5); Y = (flts(U, SYS) + 0.3 * rand(1, nsmp, 'normal')); // Rを計算 S = 15;L = 1; [R, N, SVAL] = findR(S, Y', U'); N = 3; METH = 3;TOL = -1; [A, C] = findAC(S, N, L, R, METH, TOL); A, B, C は、行列およびベクトル。それ以上でもそれ以下でもない。 これを、syslin() の引数にしている。 SYS = syslin(0.1, A, B, C); syslin() の結果は SYS に保存されたが、その内容を表示してみる。 SYS = SYS(1) (state-space system:) !lss A B C D X0 dt ! SYS(2)= A matrix = 0.5 0.1 -0.1 0.2 0.1 0. -0.1 -0.1 -0.4 -0.6 -0.7 -0.1 0.8 0. -0.6 -0.6 SYS(3)= B matrix = 0.8 0.1 1. -1. SYS(4)= C matrix = 1. 2. -1. 0. SYS(5)= D matrix = 0. SYS(6)= X0 (initial state) = 0. 0. 0. 0. SYS(7)= Time domain = 0.1 ヘルプには、syslin() は、A, B, C を元にリスト構造を作り、整合性を確認する、と書いてある。 そして、多くの関数がこの関数が作ったリスト構造を用いるらしい。 要は、データ構造を作っているだけのようだ。 nsmp = 100; U = prbs_a(nsmp, nsmp / 5); は、100個の、-1 または 1 を値としてとる、バイナリ乱数列を作る。 そして、最大20回、符号を反転させる。 U = column 1 to 11 1. 1. 1. 1. 1. 1. -1. -1. -1. -1. -1. column 12 to 22 -1. -1. -1. -1. -1. -1. -1. -1. 1. 1. -1. column 23 to 33 -1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. column 34 to 44 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. column 45 to 55 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. 1. column 56 to 66 1. -1. -1. -1. -1. -1. -1. 1. 1. 1. -1. column 67 to 77 1. 1. -1. -1. -1. -1. 1. 1. 1. -1. -1. column 78 to 88 -1. -1. -1. -1. -1. -1. -1. 1. 1. 1. -1. column 89 to 99 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. column 100 1. 符号の反転回数は、数えてない。 Y = (flts(U, SYS) + 0.3 * rand(1, nsmp, 'normal')); flts() は、前にも出た関数で、「時間応答」なるものを返すのだが、未だによく理解できていない。 とりあえず、2つに分けてみる。rand() を見てみよう。 r = 0.3 * rand(1, nsmp, 'normal')); まず、右半分を r とした。 r = column 1 to 4 0.089193 0.1592555 -0.4621402 -0.1189909 column 5 to 8 0.1548976 0.0022698 0.3126737 0.8011533 〜省略〜 column 94 to 97 0.224157 -0.017764 0.1135313 -0.4775086 column 98 to 100 -0.1171852 -0.0765337 -0.0440469 nsmp が 100 なので、100個のベクトルだ。'normal' は、平均0、分散1 の、ガウス乱数生成器を指定している。 y = flts(U, SYS) 続いて、左半分を y とした。 y = column 1 to 7 0. 0. 1.25 1.099 1.6896 1.57881 1.938271 column 8 to 11 1.8002288 -0.4304262 -0.2923465 -1.2425835 〜省略〜 column 93 to 96 2.9541454 0.3048054 3.1354941 0.5874792 column 97 to 100 3.2493289 0.7162335 3.3101171 0.7789876 あれ、なんだか以前の flts() と違い、時系列の応答みたいになっている？ --> plot(r, 'r') --> plot(y, 'b') --> plot(yr, 'g') (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image ScilabFltsPlot180517.png}} 赤が r、青が y、緑が y + r でプロットした。 つまりこれは、flts() の結果にノイズを追加したものだ。 :Amazonリンク:[新版 システム同定: ―部分空間法からのアプローチ―|https://amzn.to/2OhaTew] !!! 2018年5月18日 続き。 flts() の出力が、時系列の出力っぽい。 以前、time_id() のサンプルコードで使われていた flts() のときには、テイラー展開の係数、ではないかと考えていたが、、 さて、先を続ける。 findR() を見てみよう。 // Rを計算 S = 15;L = 1; [R, N, SVAL] = findR(S, Y', U'); findR() の説明は「線形時不変システムの行列を指定する際のプリプロセッサ」、、なんだろう？ さらに説明を見る。 findR は,コレスキーまたは(高速)QR分解および 部分空間同定法(MOESP または N4SID) 線形時不変システムの行列を推定する際の 入出力データの前処理を行い,システム次数を指定します. さらに分からない。 とりあえず、findAC()、findABCD() を呼び出すときに必ず必要になるようだ。 引数の Y', U' は、時系列の出力と入力だ。 S は、「ブロックハンケル行列のブロック行の数」とある。 「行列の行の数」ということは分かるが、「ブロックハンケル行列」とは、、？ !!! 2018年5月19日 「ブロックハンケル行列」とは何か？調べてみた。 まず「ハンケル行列」から調べてみよう。 「ハンケル行列」とは「左下から右上へ向かう対角線と平行となる行列成分が全て等しくなっている正方行列」である。 ちなみに「正方行列」は行と列の数が同じ行列。 「正方形」の「正方」と同じ。 こんな感じ。 a0 a1 a2 a3 a4 a5 a1 a2 a3 a4 a5 a6 a2 a3 a4 a5 a6 a7 a3 a4 a5 a6 a7 a8 a4 a5 a6 a7 a8 a9 a5 a6 a7 a8 a9 a10 で、「ブロックハンケル行列」とは、、結局よくわからず。 ハンケル行列が要素であるハンケル行列なのか、行列の一部分だけハンケル行列なのか、、 まあ、良い。 続き。findAC() のサンプルコードをもう一度見てみる。 //指定した線形システムからデータを生成 A = [ 0.5, 0.1, -0.1, 0.2; 0.1, 0, -0.1, -0.1; -0.4, -0.6, -0.7, -0.1; 0.8, 0, -0.6, -0.6]; B = [0.8;0.1;1;-1]; C = [1 2 -1 0]; SYS = syslin(0.1, A, B, C); nsmp = 100; U = prbs_a(nsmp, nsmp / 5); Y = (flts(U, SYS) + 0.3 * rand(1, nsmp, 'normal')); // Rを計算 S = 15;L = 1; [R, N, SVAL] = findR(S, Y', U'); N = 3; METH = 3;TOL = -1; [A, C] = findAC(S, N, L, R, METH, TOL); で、findR() は R と N と SVAL を返した。 // Rを計算 S = 15;L = 1; [R, N, SVAL] = findR(S, Y', U'); R は、60×60の巨大な行列だった。 N は、離散時間実現の次数。でも次の行で 3 を設定している。findR() が返した N の値は 2 だった。どういうこと？ SVAL は、次数の推定に使われる特異値。 !!! 2018年5月20日 続き。 N = 3; METH = 3;TOL = -1; [A, C] = findAC(S, N, L, R, METH, TOL); findAC() は、状態空間表現の 状態行列 A と 出力行列 C を見つける。 dx(t)/dt = Ax(t) + Bu(t) y(t) = Cx(t) + Du(t) その結果が、これだ。 A = 0.7192229 0.2058825 -0.0107432 -0.1481165 -0.899457 0.0510542 0.0344916 -0.0248388 0.1169426 C = 0.6758155 -0.4034059 0.0228385 それでもって、元々の A と C がこれ。 A = [ 0.5, 0.1, -0.1, 0.2; 0.1, 0, -0.1, -0.1; -0.4, -0.6, -0.7, -0.1; 0.8, 0, -0.6, -0.6]; B = [0.8;0.1;1;-1]; C = [1 2 -1 0]; 行列・ベクトルの大きさが違う、、 :Amazonリンク:[MATLAB/Simulinkとモデルベース設計による2足歩行ロボット・シミュレーション (プレミアムブックス版)|https://amzn.to/2Og0SOD] !!! 2018年5月21日 続き。 N = 3; METH = 3;TOL = -1; [A, C] = findAC(S, N, L, R, METH, TOL); 引数についてもう一度振り返る。 N = 3 は、findR() で計算された N の値が 2 だったので、再設定している。 METH は使用する手法で、3 はデフォルト設定を使うらしいが、デフォルト設定が書いてない。 が、findABCD() の説明を見ると、MOESP法らしい。 TOL は「行列のランクを推定するために使用される許容誤差」。「行列のランク」ってなんだ？と思ったら、「行列の階数」だった。 N、L は「整数」としか説明が書いてないが、これも findABCD() の説明には書いてあった。 N がシステムの次数で、L が出力の数、だ。 Nが、状態空間表現の x() のベクトル長さ、L が、y() のベクトル長さ。 S はブロックハンケル行列のブロック行の数、まだ良くわかってない。 R が findR() の計算結果。 !!! 2018年5月22日 findAC() の続き。 findAC() は、findR() の戻り値 R を引数にとる。 findAC() のサンプルコードの R は、60×60の巨大な行列だ。 さらに、左下隅から対角線直前まで値が全てゼロ。 ""[R,N] = findR(S,Y,U,METH,ALG,JOBD,TOL,PRINTW) は, 入出力データおよび離散時間実現の次数Nから構築された 結合ブロックハンケル行列の 上三角分解 R を返します. 「上三角形分解」とは「LU分解」のことのようだ。 LU分解を行うと、計算が簡単になるらしい。 R が何を表しているのかは見てもよく分からないが、、 * 状態空間表現 A, B, Cから線形システム SYS を syslin() で構成し、 * それと入力 U を元に flts() 時間応答を生成し、 * 前処理を findR() で行い、 * findAC() でシステム同定を行った ということか。 !!! 2018年5月23日 findAC() の続き。 とりあえず、状態行列 A と 出力行列 C は同定されたから、これを用い、flts() で時間応答を作ってみる。 //指定した線形システムからデータを生成 A = [ 0.5, 0.1, -0.1, 0.2; 0.1, 0, -0.1, -0.1; -0.4, -0.6, -0.7, -0.1; 0.8, 0, -0.6, -0.6]; B = [0.8;0.1;1;-1]; C = [1 2 -1 0]; SYS = syslin(0.1, A, B, C); nsmp = 100; U = prbs_a(nsmp, nsmp / 5); Y = (flts(U, SYS) + 0.3 * rand(1, nsmp, 'normal')); // Rを計算 S = 15;L = 1; [R, N, SVAL] = findR(S, Y', U'); N = 3; METH = 3;TOL = -1; // A2, C2 に入れるため、コメントアウト //[A, C] = findAC(S, N, L, R, METH, TOL); // 追加のコード // A2, C2 に入れる [A2,C2] = findAC(S,N,L,R,METH,TOL); // 次数を合わせる B2 = [1; 1; 1]; // 線形システムを構成する SYS2 = syslin(0.1, A2, B2, C2); // 時間応答を生成する Y2=flts(U, SYS2); // コンソールに A2、C2 を表示する disp(A2) disp(C2) // 時間応答をプロットする plot(Y, "r"); plot(Y2, "b"); 結果。 A2 = 0.7147033 0.05207 -0.4512613 -0.2256827 -0.9500092 -0.227936 0.1269412 -0.0745013 -0.6476534 C2 = -0.6484759 0.3669036 -0.5365736 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image ScilabFindACPlot180523.png}} 反転している、、 でも何回か繰り返していると、信号が反転しないときもある。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image ScilabFindACPlot180523_2.png}} B2 を適当に作ったからか？ :Amazonリンク:[カルマンフィルタの基礎|https://amzn.to/2M8Gguk] !!! 2018年5月24日 findAC() の続き。 昨日は同定した状態行列 A、出力行列 C を用いてプロットしたら、もとの信号と反転してしまった。 しかし、Scilab には findBD() という関数もある。 ""findBD 関数は, SLICOTルーチン IB01CD により, 状態量初期値および離散時間システムにおけるシステム行列 B と D を推定します. と書かれている。推定って、要は「同定」ってこと？ とりあえず、試してみる。 //指定した線形システムからデータを生成 A = [ 0.5, 0.1,-0.1, 0.2; 0.1, 0, -0.1,-0.1; -0.4,-0.6,-0.7,-0.1; 0.8, 0, -0.6,-0.6]; B = [0.8;0.1;1;-1]; C = [1 2 -1 0]; SYS=syslin(0.1,A,B,C); nsmp=100; U=prbs_a(nsmp,nsmp/5); Y=(flts(U,SYS)+0.3*rand(1,nsmp,'normal')); // Rを計算 S=15;L=1; [R,N,SVAL] = findR(S,Y',U'); N=3; METH=3;TOL=-1; [A2,C2] = findAC(S,N,L,R,METH,TOL); [B2, V, rcnd] = findBD(2, 1, 1, A, C, Y', U'); //B2 = [1; 1; 1]; SYS2 = syslin(0.1, A2, B2, C2); Y2=flts(U, SYS2); disp(A2) disp(C2) clf(); plot(Y, "r"); plot(Y2, "b"); 関数の 12 行目 flts ( C:\Program Files\scilab-6.0.0\modules\cacsd\macros\flts.sci 行 25 ) 実行されたファイルの 21 行目 C:\Users\takat7a\Desktop\aaa.sce error: 入力引数 #1 の型が間違っています: 文字列を指定してください. 動かない、、 !!! 2018年5月25日 findAC() の続き。 昨日は findBD() で B を計算しようと思ったら、出来なかった。 では、findBD() のサンプルコードを実行してみよう。 //　指定した線形システムからデータを生成 A = [ 0.5, 0.1,-0.1, 0.2; 0.1, 0, -0.1,-0.1; -0.4,-0.6,-0.7,-0.1; 0.8, 0, -0.6,-0.6]; B = [0.8;0.1;1;-1]; C = [1 2 -1 0]; SYS=syslin(0.1,A,B,C); nsmp=100; U=prbs_a(nsmp,nsmp/5); Y=(flts(U,SYS)+0.3*rand(1,nsmp,'normal')); // Rを計算 S=15;L=1; [R,N,SVAL] = findR(S,Y',U'); N=3; METH=3;TOL=-1; [A,C] = findAC(S,N,L,R,METH,TOL); [X0,B,D] = findBD(1,1,2,A,C,Y',U') SYS1=syslin(1,A,B,C,D,X0); Y1=flts(U,SYS1); clf();plot2d((1:nsmp)',[Y',Y1']) しかし、なぜかエラーが発生する。 APIエラー getScalarDouble: 入力引数 #7 の型が間違っています: スカラーを指定してください. の中に D = [] B = [] X0 = -1.0494476 -0.5147216 0.856072 しかし、コードを見ると、 * findAC() で 状態行列 A、出力行列 C を計算 * findBD() で入力行列 B、直達行列 D を計算 * 元の行列で作った時系列データ Y と、同定によって得られた行列で作った Y1 をプロット そのため、使い方は間違っていないようだ。 さて、以前、「システム同定」と思しき関数を列挙した。 *armax: armax 同定 *armax1: armax 同定 *findABCD: 離散時間システムの部分空間の同定 *findAC: 離散時間システム部分空間の同定 *time_id: SISO 最小二乗同定 あと残っているのは、findABCD() だけだ。 findAC() を使ってみて、結局 findBD() を使うのだったら、findABCD() でいいか、と思っている。 :Amazonリンク:[カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― (統計学One Point 2)|https://amzn.to/2M3MPio] !!! 2018年5月26日 早速、findABCD() のサンプルコードを実行する。 // 指定した線形システムからデータを生成 A = [ 0.5, 0.1,-0.1, 0.2; 0.1, 0, -0.1,-0.1; -0.4,-0.6,-0.7,-0.1; 0.8, 0, -0.6,-0.6]; B = [0.8;0.1;1;-1]; C = [1 2 -1 0]; SYS=syslin(0.1,A,B,C); nsmp=100; U=prbs_a(nsmp,nsmp/5); Y=(flts(U,SYS)+0.3*rand(1,nsmp,'normal')); // Rを計算 S=15; [R,N1,SVAL] = findR(S,Y',U'); N=3; SYS1 = findABCD(S,N,1,R) ;SYS1.dt=0.1; SYS1.X0 = inistate(SYS1,Y',U'); Y1=flts(U,SYS1); clf();plot2d((1:nsmp)',[Y',Y1']) コードを見ていく。 // 指定した線形システムからデータを生成 A = [ 0.5, 0.1,-0.1, 0.2; 0.1, 0, -0.1,-0.1; -0.4,-0.6,-0.7,-0.1; 0.8, 0, -0.6,-0.6]; B = [0.8;0.1;1;-1]; C = [1 2 -1 0]; SYS=syslin(0.1,A,B,C); syslin() で状態空間表現 SYS を構成している。 nsmp=100; U=prbs_a(nsmp,nsmp/5); 入力信号 U を生成。prbs_a() は、疑似バイナリ乱数列を生成する。 Y=(flts(U,SYS)+0.3*rand(1,nsmp,'normal')); flts() を用いて、SYS に 入力信号 U を適用し、時系列の出力信号 Y を生成。 rand() でノイズを追加している。 // Rを計算 S=15; [R,N1,SVAL] = findR(S,Y',U'); findR() は、findABCD() を呼び出す前の前処理。 N=3; SYS1 = findABCD(S,N,1,R) ;SYS1.dt=0.1; findABCD() によるシステム同定。 SYS1.X0 = inistate(SYS1,Y',U'); inistate() は、今まで出てこなかった。 「離散時間システムの状態量の初期値を推定」する。 Y1=flts(U,SYS1); 同定した状態空間表現 SYS1 を用いて、時系列の出力信号 Y を生成している。 clf();plot2d((1:nsmp)',[Y',Y1']) 最後に、元の信号 Y と 推定した信号 Y1 をグラフにプロットし、比較している。 いい感じに再現しているように思う。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image ScilabFindAbcdPlot180526.png}} !!! 2018年5月27日 さて、これまでシステム同定用関数と思われる関数、 *armax: armax 同定 *armax1: armax 同定 *findABCD: 離散時間システムの部分空間の同定 *findAC: 離散時間システム部分空間の同定 *time_id: SISO 最小二乗同定 の動作を一通り調べた。 次は、以前作ったCR回路の時系列データを armax1(), armax(), findABCD() で同定してみよう。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image RcCircuit180508.png}} (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image RcGraph180508.png}} まず、簡単そうな armax1() から。 文法は次の通り。 [arc, resid] = armax1(r, s, q, y, u [, b0f]) まず、r, s, q は、自己回帰の係数、、、 初っ端から分からない。 とりあえず、サンプルコードと同じ「3, 3, 2」にしてみよう。 ｙは出力信号。これは、以前作ったCR回路の時系列データの出力電圧とする。 rc = fscanfMat("RcData180508.txt"); v = rc(:, 2); u は入力信号。ステップ信号だから、最初から同じ値にする。 u = ones(size(v, "r"), 5); そして、armax1() を実行してみた。 // armax1() による同定 [arc1, resid] = armax1(3, 3, 2, v, u); 結果を表示してみる。 arc1 = A(z^-1)y=B(z^-1)u + D(z^-1) e(t) A(x) = 2 3 1 -0.832217x -0.2756286x +0.1641095x B(x) = 2 3 0.0142063 +0.0142063x +0.0142063x +0.0142063x D(x) 2 1 -0.1996464x -0.1755013x e(t)=Sig*w(t); w(t) 1-dim white noise Sig= | 0.0514156 | 先程の r, s, q は、それぞれ A, B, D の次数を示すらしい。 さて、narsimul() でシミュレーションをしてみる。 ヘルプを見ていて、arsimul() という関数もあることに気付いた。 ともあれ、narsimul() でやってみる。 rc = fscanfMat("RcData180508.txt"); v = rc(:, 2); u = ones(size(v, "r"), 1) * 5; // armax1() による同定 [arc1, resid] = armax1(3, 3, 2, v, u); y = narsimul(arc1, u'); clf(); plot(v, "b") plot(y, "r"); その結果が下のグラフ。 青が元データ、赤が ARMAX同定結果。 何かがおかしい、、 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image ScilabArmaxRcPlot180527.png}} :Amazonリンク:[非線形カルマンフィルタ|https://amzn.to/2Ol0dLB] !!! 2018年5月28日 続き。 青がRC回路の電圧なのだが、LTspice 上で描画されたグラフと立ち上がりが違う。 LTspice でエクスポートしたデータには時刻も記録されているのだが、実は時間間隔が等間隔ではないことが分かった。 一応、トリッキーな方法で出来なくはないようだが、普通の方法でやりたい。 では、と、[MPLAB Mindi|https://www.microchip.com/mplab/mplab-mindi] を試してみたが、等間隔で出力されない。 結局、[OpenModelica|https://openmodelica.org/] を試したところ、上手くいった。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image RcCircuitOm180528.png}} データ: [RcData180527.txt|http://moi2.sakura.ne.jp/fswiki/wiki.cgi?page=Scilab%A4%F2%BB%C8%A4%C3%A4%C6%A4%DF%A4%EB&file=RcData180508%2Etxt&action=ATTACH] このデータでやってみたところ、、 rc = fscanfMat("RcData180527.txt"); v = rc(:, 2); u = ones(size(v, "r"), 1) * 5; // armax1() による同定 [arc1, resid] = armax1(3, 3, 2, v, u); y = narsimul(arc1, u'); clf(); plot(v, "b") plot(y, "r"); (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image ScilabArmaxRcPlot180528.png}} 青が元データ、赤が ARMAX同定結果。 RCのグラフも出たし、armax1() の同定結果でシミュレーションした結果もピッタリあった。 すごい。 !!! 2018年5月29日 昨日のつづき。 ところで、armax1() の引数 r, s, q は、前回はサンプルコードと同じ 3, 3, 2 としたが、本来はどのように決めればよいのか？ 試しに、r, s, q = 0, 0, 0 でやってみる。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image ScilabArmaxRcPlot180529.png}} まったくだめだ。 次に、r, s, q = 1, 1, 1 で。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image ScilabArmaxRcPlot180529_2.png}} あれ、r, s, q = 3, 3, 2 とあまり変わらない。 !!! 2018年5月30日 続き。 armax() も試してみる。 // armax() による同定 rc = fscanfMat("RcData180527.txt"); v = rc(:, 2); u = ones(size(v, "r"), 1) * 5; [arc2, la2, lb2, sig2, resid2] = armax(1, 1, v', u'); y2 = narsimul(arc2, u'); clf(); plot(v, "b") plot(y2, "r"); (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180530ScilabArmaxRcPlot.png}} ほとんど同じみたいだ。 ただし、v, u は v', u' と指定しないといけない。 また、次のような警告が出た。 意味は分からない。 警告: armax: z*z' は数値的に特異です。 armax: 推定器の標準偏差 a: 0.00000 0.00123 armax: 推定器の標準偏差 b: 0.00057 0.00057 :Amazonリンク:[基礎からわかる時系列分析 ―Rで実践するカルマンフィルタ・MCMC・粒子フィルター (Data Science Library)|https://amzn.to/2Ma0HqA] !!! 2018年5月31日 続き。 armax1() と armax() で同定した値を、同じグラフに書いてみる。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180531ScilabArmaxRcPloty1y2.png}} ほぼ同じだ、、 同定したARMAX表現は、、まず armax1()。 arc1 = A(z^-1)y=B(z^-1)u + D(z^-1) e(t) A(x) = 1 -0.9899243x B(x) = 0.0050324 +0.0050324x D(x) 1 -0.0043229x e(t)=Sig*w(t); w(t) 1-dim white noise Sig= | 0.0001269 | 次に armax()。 arc2 = A(z^-1)y=B(z^-1)u + D(z^-1) e(t) A(x) = 1 -0.9902668x B(x) = 0.0048726 +0.0048726x D(x) 1 e(t)=Sig*w(t); w(t) 1-dim white noise Sig= | 0.0015570 | 若干異なるが、だいたい同じかな、、 状態空間表現は、、まず armax1()。 ss1 = ss1(1) (state-space system:) !lss A B C D X0 dt ! ss1(2)= A matrix = 0.9899243 0.0050324 0. 0. ss1(3)= B matrix = 0.0050324 1. ss1(4)= C matrix = 1. 0. ss1(5)= D matrix = 0. ss1(6)= X0 (initial state) = 0. 0. ss1(7)= Time domain = d armax() は、、 ss2 = ss2(1) (state-space system:) !lss A B C D X0 dt ! ss2(2)= A matrix = 0.9902668 0.0048726 0. 0. ss2(3)= B matrix = 0.0048726 1. ss2(4)= C matrix = 1. 0. ss2(5)= D matrix = 0. ss2(6)= X0 (initial state) = 0. 0. ss2(7)= Time domain = d :Amazonリンク:[信号・システム理論の基礎―フーリエ解析、ラプラス変換、z変換を系統的に学ぶ|https://amzn.to/2wf6cdM] !!! 2018年6月1日 続き。 armax()、armax1() で得られた状態空間表現を、伝達関数で表現してみる。 これは単純に tf = ss2tf(ss) で得られる。 armax1() の状態空間表現を伝達関数表現にすると、、 tf1 = 0.0050324 + 0.0050324z ----------------------- 2 -0.9899243z + z armax() の状態空間表現を伝達関数表現にすると、、 tf2 = 0.0048726 + 0.0048726z ----------------------- 2 -0.9902668z + z まあ、ほぼ同じだ。 ただ、RC回路って、１次遅れ、 K G(s) = -------- 1 + Ts になるはずなんだが、これは2次遅れになっている？ !!! 2018年6月2日 続き。 昨日までarmax()、armax1() を試していた。 次は、findABCD() を使ってみる。 こんなコード。 // findABCD() による同定 rc = fscanfMat("RcData180527.txt"); v = rc(:, 2)'; u = ones(size(v, "c"), 1)' * 5; S = 1; [R, N, SVAL] = findR(S, v, u); sys = findABCD(S, N, 1, R); sys.dt = 1; sys.X0 = inistate(sys, v, u); y3 = flts(u, sys); scf(); plot(v, "b"); plot(y2, "r"); 結果。 sorder: サンプル数を 4005 以上としてください サンプル数が少ない？ 1000データあるのに、、 仕方ない、データを増やそう、、 シミュレーションは 1 秒間だけど、5 秒間にしたら 1 msec単位だから5,000だ。 sorder: サンプル数を 20005 以上としてください だめか、、 あのRCの電圧は 1 秒で飽和してるからな、、 では、サンプリング周期を 0.2 msec 単位にしてみよう。 sorder: サンプル数を 20005 以上としてください だめだ、、 !!! 2018年6月3日 なぜか findR() で、4倍のサンプル数を要求される。 データ数を増やすと、さらにその4倍だ。 findR() 自体は、sorder() という関数を呼び出しているだけのラッパー関数である。 そこで、sorder() を直接呼び出してあれこれやってみたが、状況は変わらず、、 !!! 2018年6月4日 RC回路はステップ応答なので、情報量が少ないのだろうか？ もっと信号を増やしてみる。 !!! 2018年6月5日 続き。 信号を増やす方法として、サンプルコードでよく出てきた prbs_a() を用いる。 サンプルコードだと、信号反転回数をデータ数の1/5にしていたので、 u = prbs_a(1000, 200); これをCSVファイルに保存する。 csvWrite(u, 'rand.csv'); これをExcelで開き、OpenModelica の Modelica.Electrical.Analog.Sources.TableVoltage に設定するテーブルを作る。 !!! 2018年6月6日 続き。 prbs_a() を csvWrite() で出力すると、-1 と 1 を横一列に出力する。 そこで、Excelで経過秒と電圧を並べた形式、 0,0 0.001,5 のような 0.000 秒から 1.000 秒までで、0 or 5 の値をとるCSVファイルになるように整形する。 さらにテキストエディタ gvim で、 * 全ての末尾に ';' を追加する: %s/\(.*\)/\1; * 先頭行の頭に '[' を追加する * 最終行の末尾の ';' を削除し、 ']' を追加する * 全部1行につなげる: 先頭で qaJq 999@a これで、OpenModelica の Modelica.Electrical.Analog.Sources.TableVoltage に設定する値が完成、フィールドに設定する。 これで、シミュレーションを行い、結果をCSVで保存する。 !!! 2018年6月10日 続き。 時間がかかってしまった、、 OpenModelica の TableVoltage でのシミュレーションは上手くいったのだが、出力されたCSVファイルはなぜか同一行が3つあったため。 よくよく考えてみると、CR回路の伝達関数は K G(s) = -------- 1 + Ts と分かっているのだから、Scilab でシミュレーションすれば良い。 // 5Vのステップ入力 u2 = ones(1, 1000) * 5; // CR回路の伝達関数を作る。定常ゲイン 1, 時定数 0.1 秒 s = poly(0, 's'); G = 1 / (1 + 0.1 * s); // 線形システムを作る sys2 = syslin('c', G); // 線形システムを離散化する sl = dscr(tf2ss(sys2), 0.001); // ステップ入力を適用し、時間応答を生成する v2 = flts(u2, sl); (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180610RcCircuit.png}} できた。 元のCRグラフと重ねてプロットしてもほとんど重なっている。 では、ランダム信号を生成し、CR回路に適用してみよう。 変更するのは u2 のみ。 u2 = prbs_a(1000, 200); u2 = u2 + 1; u2 = u2 * 2.5; prbs_a() は -1 と 1 を出力するので * -1 を 0 にするために 1 を足す * 1 を足したので、1 は 2 になっている。だから、2.5 を掛ければ 5 になる (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180610RcCircuit_2.png}} それっぽい信号が出来た。 では、findABCD() に適用してみよう。 結果、 sorder: サンプル数を 4001 以上としてください だめだ、、と思っていたら、出来た。 今まで「サンプル数をxxxx以上としてください」と言われ、S を小さくすると数が減っていったので S = 1 にしていた。 しかし、ある程度大きくないといけないらしい。15 に戻した。なお、3より小さくすると Scilab が落ちた。 S = 15; サンプルコードにあった、乱数を追加した。これがないと、Scilab が落ちる。 ただし、乱数の大きさによっては同定レベルが下るので、あまり大きくないほうが良さそう。 また入力信号がステップ入力だと、やはり Scilab が落ちるので、ランダム信号が必要。 v2 = v2 + 0.01 * rand(1, 1000, 'normal'); 入力は列ベクトルでないとだめらしい。行ベクトルだと「サンプル数をxxxx以上としてください」と言われ続ける。 [R, N, SVAL] = findR(S, v2', u2'); sys = findABCD(S, N, 1, R); sys.dt = 0.001; ここも列ベクトルに sys.X0 = inistate(sys, v2', u2'); (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180610ScilabFindAbcdPlot.png}} !!! 2018年6月11日 続き。findABCD() で得られた sys は、次のようになった。 状態空間表現: sys = sys(1) (state-space system:) !lss A B C D X0 dt ! sys(2)= A matrix = 0.9900688 sys(3)= B matrix = -0.0359226 sys(4)= C matrix = -0.2763687 sys(5)= D matrix = -0.0001233 sys(6)= X0 (initial state) = -0.018516 sys(7)= Time domain = 0.001 伝達関数表現: --> ss2tf(sys) ans = 0.0100500 - 0.0001233z ----------------------- -0.9900688 + z 元の状態空間表現は、 伝達関数 (連続時間): 1 --------- 1 + 0.1s 状態空間表現: sl(1) (state-space system:) !lss A B C D X0 dt ! sl(2)= A matrix = 0.9900498 sl(3)= B matrix = 0.0031465 sl(4)= C matrix = 3.1622777 sl(5)= D matrix = 0. sl(6)= X0 (initial state) = 0. sl(7)= Time domain = 0.001 伝達関数 (離散時間): 0.0099502 --------------- -0.9900498 + z 状態空間表現は、状態行列 A しか合っていないように気もするが、、 元のグラフとぴったり重なっているということは、肝心なところは合っている、ということか？ 離散時間の伝達関数は、得られた伝達関数には分子に z が入っているのが気になるが、それ以外は係数などよく合っている。 ただ、不満なのが、ステップ入力とその応答データだと Scilab が落ちてしまう点。 なんとかならないものか、、 !!! 2018年6月12日 続き。 findABCD() で、ステップ入力でも同定ができるようにするには、、 現在のCR回路のデータは1ミリ秒単位で1秒のデータだが、多いようだ。 立ち上がり時間 (ステップ応答の10％から90％までの時間) に 5 〜 8 サンプリングで良いようだ。 調べたところ、立ち上がり時間が約220ミリ秒。 5 〜 8 サンプリングだと、27.5 〜 44 ミリ秒。 なので、30ミリ秒のサンプリングにしてみる。 !!! 2018年6月13日 続き。 CR回路を30ミリ秒サンプリングにしてみた。 smpl = 33; dt = 0.03; // CR回路, ステップ入力 u2 = ones(1, smpl) * 5; s = poly(0, 's'); G = 1 / (1 + 0.1 * s); sys2 = syslin('c', G); sl = dscr(tf2ss(sys2), dt); v2 = flts(u2, sl); 結果は、、「サンプル数をxxxx以上としてください」が出る、、 !!! 2018年6月16日 続き。 ステップ入力データだと findABCD() がどうしても実行できない。 そこで色々やってみたところ、、できた。 ステップ入力信号にもノイズを足してやると良い。 // findABCD (ステップ入力) // 入力・出力ともにノイズを加算する uStepNoize = u_step + 0.01 * rand(1, smpl, 'normal'); vStepNoize = v_step + 0.01 * rand(1, smpl, 'normal'); [R, N, SVAL] = findR(S, vStepNoize', uStepNoize'); sys2 = findABCD(S, N, 1, R); sys2.dt = dt; sys2.X0 = inistate(sys2, vStepNoize', uStepNoize'); vAbcdStep2 = flts(u_step, sys2); scf(); plot(x, vStepNoize, "b"); plot(x, vAbcdStep2, "r"); plot(x, uStepNoize, "c"); title("Subspace identification (Step)"); legend(["v_step", "v_findABCD", "u_step"], -4); xgrid(); 結果、同定できた。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180616ScilabFindAbcdStepPlot.png}} findABCD() とランダム信号を使って同定した状態空間表現は、 sys = sys(1) (state-space system:) !lss A B C D X0 dt ! sys(2)= A matrix = 0.9753136 sys(3)= B matrix = -0.0811942 sys(4)= C matrix = -0.3040028 sys(5)= D matrix = 0.0000259 sys(6)= X0 (initial state) = -0.0040826 sys(7)= Time domain = 0.0025 findABCD() と、ノイズを加えたステップ入力・応答信号を使って同定した状態空間表現は、 sys2 = sys2(1) (state-space system:) !lss A B C D X0 dt ! sys2(2)= A matrix = 0.9754135 sys2(3)= B matrix = -0.0820994 sys2(4)= C matrix = -0.3038203 sys2(5)= D matrix = -0.0145246 sys2(6)= X0 (initial state) = -0.2721375 sys2(7)= Time domain = 0.0025 大体合っているっぽい。 !!! 2018年6月17日 続き。 findABCD() へのステップ入力・応答信号にノイズを加えることで、システム同定ができるようになった。 ちなみにサンプリング数を極端に少なくすると「サンプル数を89以上としてください」が出るが、89以上に設定すれば問題ない。 こちらがサンプル数89のとき。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180617ScilabCrIdentificationPlot89.png}} サンプル数1000のときはこちら。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180617ScilabCrIdentificationPlot1000.png}} サンプル数の少ないほうが若干各グラフがずれているが、ほとんど変わらない。 !!! 2018年6月18日 続き。 ところで、今までCR回路のデータは1000サンプリング/秒、1秒間データだった。 これを、10サンプリング/秒、10秒間にするとどうなるだろう？ この場合、ステップ入力だと最初の1秒間のみ信号が変化し、あとは平坦になる。 結果はこのようになった。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180618ScilabCrIdentificationPlot10Smpl10Sweep.png}} ステップ入力の findABCD() は発振しているが、ランダム信号の findABCD() は大丈夫だった。 意外にも、ステップ入力の armax1()、armax() はそれなりに同定している。 !!! 2018年6月19日 続き。 ところで、findABCD() によるステップ入力による同定では、信号にノイズを加えないと Scilab が落ちてしまった。 では、armax1()、armax() の同定信号にノイズを加えたらどうなるか？ やってみた。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180619ScilabCrIdentificationPlotWithNoize.png}} 拡大したところ。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180619ScilabCrIdentificationPlotWithNoizeExpand.png}} armax1()、armax() はノイズに弱く、少サンプル数に強い。 findABCD() は、ノイズに強く、少サンプル数に弱い。 !!! 2018年6月21日 続き。 ところで、物の本 (システム同定の基礎, 足立修一著) によると、「直流成分を取り除く」ほうが良いらしい (Step 3 入出力データの前処理 P. 23)。 除去したらどの様になるか？ (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180621ScilabFindAbcdPlotAc.png}} あまり良くない、、 !!! 2018年6月22日 続き。 昨日のは、findABCD() だからダメなのかも、と思い、armax() で試してみる。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180622ScilabArmaxPlotAc.png}} 、、、なんか変だ。 !!! 2018年6月23日 続き。 直流成分の除去をしなくても問題なさそうだし、除去すると同定がうまく出来なかったので、放って置く。 ところで、今までCR回路にステップ入力を与えたときの応答を同定していた。 逆に、応答から入力を同定できるだろうか？ 試してみる。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180623ScilabArmaxPlotRev.png}} !!! 2018年6月26日 続き。 応答信号から入力信号へのシステム同定は難しそうなので、次に進む。 CR回路の応答に、むだ時間を追加した信号はどうだろう？ 実は、むだ時間を記述するには工夫が必要だと判明。 むだ時間は「パデ近似」で表現するのだとか。 MATLABでは pade() という関数がある。 Scilabでは、、ない。 というわけで、パデ近似を定義する。 !!! 2018年7月2日 続き。 パデ近似を実装してみよう。 MATLABのpade関数の説明サイト、「[むだ時間のパデ近似|https://jp.mathworks.com/help/control/ref/pade.html] 」によると、0.1秒遅延の3次パデ近似は、次のようになる。 sysx = -s^3 + 120 s^2 - 6000 s + 1.2e05 -------------------------------- s^3 + 120 s^2 + 6000 s + 1.2e05 [つまらない技術者の日記|https://blogs.yahoo.co.jp/der_winter] の、 [scilab vol.2|https://blogs.yahoo.co.jp/der_winter/10103902.html] に書かれたコードだと、 ans = 2 3 216 - 10.8s + 0.18s - 0.001s ----------------------------- 2 3 216 + 10.8s + 0.18s + 0.001s なんだか係数が異なる、、 [Scilabで学ぶフィードバック制御|http://programmer-life.net/index.html] の、 [周波数応答：pade関数の作成|http://programmer-life.net/tfunc_090.html] をそのまま使ってみると、 --> [norm, denom] = pade(0.1, 3) denom = 120000. 6000. 120. 1. norm = -120000. -6000. 120. -1. 係数はMATLABのと同じだ。ただし、伝達関数形式ではないが、、 というわけで、後者の関数を伝達関数表現にして用いる。 !!! 2018年7月3日 続き。 係数の配列を次のように変換したい。 denom: 12000 + 6000s + 120s^2 + s^3 norm: -12000 - 6000s + 120s^2 - s^3 どうやら、inv_coeff() がそれに相当するようだ。 A = 4. 6. 7. 0. 0. 8. 0. 6. 7. 1. 0. 0. --> inv_coeff(A, 5) ans = 2 5 4 +6x +7x +8x 2 3 6x +7x +x !!! 2018年7月4日 続き。 inv_coeff() を使ってみた。 --> disp(inv_coeff(pNorm',3, 's')) 2 3 -120000 -6000s +120s -s できた。 !!! 2018年7月6日 続き。 pade() によるむだ時間を伝達関数に追加してみた。 G = 1 / (1 + 0.1 * s); を G = padeSys / (1 + 0.1 * s); のようにするだけ。 padeSys は、pade() と inv_coeff() で組み立てた、むだ時間の伝達関数。 プロットしてみた。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180706ScilabPadeDelayPlot.png}} 、、グラフが逆さまになっているのと、振幅が異なっている。 あの pade() の係数はよく見ると、MATLABのそれとは異なっている。 ダメだ、、でも意外と、findABCD() の同定は頑張っている。 !!! 2018年7月9日 続き。 もう一つの pade()、 [つまらない技術者の日記|https://blogs.yahoo.co.jp/der_winter] の、 [scilab vol.2|https://blogs.yahoo.co.jp/der_winter/10103902.html] に書かれたコードを使ってみる。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180709ScilabPade2DelayPlot.png}} 遅れの最後の方で発振する感じになっているが、振幅は元と同じだし、まあ良いかもしれない。 とりあえこれを使おう。 青いグラフは findABCD() で同定したグラフだが、いまいち。 これでも、何回か繰り返した中で一番マシな方。 さて、遅れ時間を含む伝達関数を同定するにはどうすればいいだろうか？ !!! 2018年7月10日 続き。 Scilabでむだ時間を含んだシステムの同定について調べているが、ほぼ情報が出てこない。 !!! 2018年7月11日 続き。 現在のコード。 clc(); clear(); xdel(winsid()); //https://blogs.yahoo.co.jp/der_winter/10103902.html から拝借したコード function delay_sys = pade( L, order ) num = ( - L*%s + 2*order ) ^order; den = ( L*%s + 2*order ) ^order; delay_sys = num/den; endfunction // CR回路 (一次遅れ) smplPerSec = 1000; sweepTime = 1; smpl = sweepTime * smplPerSec; dt = 1 / smplPerSec; // 疑似乱数入力 u_rand = prbs_a(smpl, smpl / 20); u_rand = u_rand + 1; u_rand = u_rand * 2.5; // ステップ入力 u_step = ones(1, smpl) * 5; // 時間軸 x = linspace(0, sweepTime, smpl); s = poly(0, 's'); // むだ時間 deadTime = 0.1 order = 2; padeSys = pade(deadTime, order); disp(padeSys); if deadTime == 0 then padeSys = 1; end // CR: 定常ゲイン1, 時定数0.1 G = padeSys / (1 + 0.1 * s); //G = 1 / (1 + 0.1 * s); sysCr = syslin('c', G); sl = dscr(tf2ss(sysCr), dt); v_rand = flts(u_rand, sl); v_step = flts(u_step, sl); // ノイズ加算 uStepNoize = u_step + 0.001 * rand(1, smpl, 'normal'); vStepNoize = v_step + 0.001 * rand(1, smpl, 'normal'); uRandNoize = u_rand + 0.001 * rand(1, smpl, 'normal'); vRandNoize = v_rand + 0.001 * rand(1, smpl, 'normal'); // ステップ入力・ランダム入力の応答 scf(); plot(x, uStepNoize, "g"); plot(x, vStepNoize, "r"); plot(x, vRandNoize, "b"); plot(x, uRandNoize ,"c"); title('CR circuit'); legend(["u_step", "v_step", "v_rand", "u_rand"], -4); xgrid(); // findABCD() による同定 S = 15; [R, N, SVAL] = findR(S, vRandNoize', uRandNoize'); ssAbcd = findABCD(S, N, 1, R); ssAbcd.dt = dt; ssAbcd.X0 = inistate(ssAbcd, vRandNoize', uRandNoize'); vAbcdStep = flts(u_step, ssAbcd); scf(); plot(x, uRandNoize, "c"); plot(x, vRandNoize, "m"); plot(x, v_step, "r"); plot(x, vAbcdStep, "b"); title("Subspace identification"); legend(["u_rand", "v_random", "v_step", "v_findABCD"], -4); xgrid(); // findABCD (ステップ入力) S = 15; [R, N, SVAL] = findR(S, vStepNoize', uStepNoize'); ssAbcdStep = findABCD(S, N, 1, R); ssAbcdStep.dt = dt; ssAbcdStep.X0 = inistate(ssAbcdStep, vStepNoize', uStepNoize'); vAbcdStep2 = flts(u_step, ssAbcdStep); scf(); plot(x, uStepNoize, "c"); plot(x, v_step, "r"); plot(x, vAbcdStep2, "b"); title("Subspace identification (Step)"); legend(["u_step", "v_step", "v_findABCD"], -4); xgrid(); 現在までに分かっていることは、 *パデ関数 (pade() ) の次数を上げると、findABCD() がうまく同定できない *findR() の引数 S (ブロックハンケル行列のブロック行の数) の値を増やしていく (15→100) と、同定能力が上がる *パデ関数 (pade() ) の次数を上げると、ステップ信号を与えた方の findABCD() がうまく同定できない *ステップ信号による同定の場合、findR() が返す次数 N を findABCD() に与えるよりも、自分で次数を指定したほうが上手くいく !!! 2018年7月13日 続き。 むだ時間0.1秒、時定数0.1秒、定常値5のシステムのステップ信号を findABCD() で同定すると、次のようになる。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180713ScilabPade2FindAbcdStepN1.png}} findR() が返した次数 1 をそのまま使うと、こうなる。 そこで、自分で次数を上げてあげる。 次数を 2 にしたとき。実行するたびにグラフの形は変わるが、これは割と良いほう。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180713ScilabPade2FindAbcdStepN2.png}} さらに次数を 3 にしたとき。だいたい合ってきた。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180713ScilabPade2FindAbcdStepN3.png}} !!! 2018年7月17日 続き。 むだ時間0.1秒、時定数0.1秒、定常値5のシステムのステップ信号を findABCD() で同定する、の続き。 findR() の引数 S (ブロックハンケル行列のブロック行の数) の値は、現在は 15。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180713ScilabPade2FindAbcdStepN3.png}} この S を、15 から 30 に増やしてみる。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180717ScilabPade2FindAbcdStepN3S30.png}} S=15 のときは、元のグラフと合わないときが結構あるが、S=30 だと、ほぼ同じになる。 次は、2次遅れの同定をしてみよう。 !!! 2018年7月18日 2次遅れの伝達関数は、次の形になる。 K ωn^2 G = -------------------------- s^2 + 2 ζ ωn s + ωn^2 K は定常ゲイン、ωn は固有周波数、ζ は減衰係数。 とりあえず、グラフを描いてみよう。 !!! 2018年7月19日 続き。 2次遅れのグラフを描いてみる。 wn = 10; zeta = 1; G = wn^2 / (s^2 + 2 * zeta * wn * s + wn^2); sysCr = syslin('c', G); sl = dscr(tf2ss(sysCr), dt); v_rand = flts(u_rand, sl); v_step = flts(u_step, sl); ωn = 10, ζ = 1 で描いてみた。それっぽいのができた。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180719Scilab2ndOrderDelay.png}} 紫のグラフは、今までの1次遅れ。 !!! 2018年7月20日 続き。 2次遅れでも、ステップ信号を入力するよりも、擬似乱数信号を入力するほうが、同定能力が上がる。 同定信号として擬似乱数信号を入力したときのグラフ。 赤 (元の伝達関数のステップ応答) と青 (同定した状態空間表現のステップ応答) はほぼ重なっている。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180720Scilab2ndOrderDelayFindAbcd.png}} 同定信号としてステップ信号を入力したときのグラフ。 赤と青はずれている。推定された次数が、2ではなく1になっている。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180720Scilab2ndOrderDelayFindAbcdStep.png}} ただし、次数を正しく2と推定するときは、ステップ信号による同定でもそれなりのグラフとなる。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180720Scilab2ndOrderDelayFindAbcdStepN2.png}} !!! 2018年7月21日 続き。 2次遅れは、むだ時間に比べ同定が楽かもしれない。 同定された状態空間表現を伝達関数に変換する。 元の離散時間の伝達関数は、、 --> ss2tf(sl) ans = 0.0000493 + 0.0000497z --------------------------- 2 0.9801987 - 1.9800997z + z ステップ信号による同定・擬似乱数信号による同定結果を離散時間の伝達関数にすると --> ss2tf(ssAbcdStep) ans = 2 -0.0128824 + 0.0323398z - 0.0193583z ------------------------------------- 2 0.9801994 - 1.9801004z + z --> ss2tf(ssAbcd) ans = 2 0.0000055 + 0.0001289z - 0.0000359z ------------------------------------ 2 0.9803086 - 1.9802100z + z 分子の部分にz^2が現れているのが異なるところ。 !!! 2018年7月23日 離散時間の伝達関数から連続時間の伝達関数に変換しようと思ったが、Scilabには存在しないようだ。 MATLABだと、c2d()、d2c()という関数があるらしい。 !!! 2018年7月26日 続き。 離散時間伝達関数を差分方程式に変換するにはどうすれば良いだろうか？ ネットで調べても、引っかからない。 そもそもそういうことは必要性がない？ でも、連続時間伝達関数を微分方程式に戻し、そこから差分方程式にする、というのはよくあるのではないか？と思ったのだが、、 !!! 2018年7月30日 ちょっとz変換について調べる。 次の1次遅れの連続時間伝達関数、 s = poly(0, 's); Gc = 1 / (1 + 0.1 * s); --> Gc Gc = 1 --------- 1 + 0.1s これを、状態空間表現に変換してみる。 --> SS = tf2ss(Gc) SS = SS(1) (state-space system:) !lss A B C D X0 dt ! SS(2)= A matrix = -10. SS(3)= B matrix = 3.1622777 SS(4)= C matrix = 3.1622777 SS(5)= D matrix = 0. SS(6)= X0 (initial state) = 0. SS(7)= Time domain = [] さらに、離散時間の状態空間表現に変換してみる。 --> SSd = dscr(SS) 警告: dscr: 入力引数 #1 は連続時間としてください. SSd = SSd(1) (state-space system:) !lss A B C D X0 dt ! SSd(2)= A matrix = 0.9900498 SSd(3)= B matrix = 0.0031465 SSd(4)= C matrix = 3.1622777 SSd(5)= D matrix = 0. SSd(6)= X0 (initial state) = 0. SSd(7)= Time domain = 0.001 なぜか警告が出た。 さらに、離散時間伝達関数に変換してみる。 --> Gd = ss2tf(SSd) Gd = 0.0099502 --------------- -0.9900498 + z !!! 2018年8月2日 続き。 それぞれ、ボード線図を描いてみる。 Gc = 1 / (1 + 0.1 * s); SSc = tf2ss(Gc); SSd = dscr(SSc); Gd = ss2tf(SSd); scf(); bode(SSc); scf(); bode(Gd); すると、 警告: dscr: 入力引数 #1 は連続時間としてください. 関数の 17 行目 repfreq ( C:\Program Files\scilab-6.0.0\modules\cacsd\macros\repfreq.sci 行 29 ) 関数の 41 行目 bode ( C:\Program Files\scilab-6.0.0\modules\cacsd\macros\bode.sci 行 52 ) 実行されたファイルの 392 行目 C:\RcAnalyze180527.sce error: 入力引数 #2 の型が間違っています. どうやら、連続時間の状態空間表現の Time Domain が [] なのがダメらしい。 Gc = 1 / (1 + 0.1 * s); SSc = tf2ss(Gc); // 追加 SSc(7) = 'c'; SSd = dscr(SSc); Gd = ss2tf(SSd); scf(); bode(SSc); scf(); bode(Gd); すると、 警告: calfrq: ナイキスト周波数を越える周波数は、無視されます。 という警告が出て、グラフが表示された。 こちらが連続時間。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180803Scilab2ndOrderDelayBodeCont.png}} こちらが離散時間。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180803Scilab2ndOrderDelayBodeDesc.png}} グラフの形が違う？ と思ったら、グラフの軸の範囲が異なる。 「ナイキスト周波数〜」の警告は、離散時間の場合は意味がないからそれ以降のグラフは出していないようだ。 !!! 2018年8月6日 続き。 まとめて描いた。 (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180806Scilab2ndOrderDelayBodeBoth.png}} 重なっているっぽいが、、若干変な形。 と思ったら、dscr() で離散化するとき、サンプル周期を指定できる。 SSdの値を確認すると、サンプル周期は 0.001 だった。 SSd(7)= Time domain = 0.001 これを、もっと細かくしてみる。 0.0001だと、、 Gc = 1 / (1 + 0.1 * s); SSc = tf2ss(Gc); SSc(7) = 'c'; SSd = dscr(SSc, 0.0001); // ここ Gd = ss2tf(SSd); scf(); bode([SSc;Gd], 10^-3, 10^3, 0.0001, ['SSc';'Gd']); // ここ (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180806Scilab2ndOrderDelayBodeBoth2.png}} 変わった。 !!! 2018年8月7日 続き。 さらにサンプル周期を細かくしてみる。 0.00001 だと、、 Gc = 1 / (1 + 0.1 * s); SSc = tf2ss(Gc); SSc(7) = 'c'; SSd = dscr(SSc, 0.00001); // ここ Gd = ss2tf(SSd); scf(); bode([SSc;Gd], 10^-3, 10^3, 0.00001, ['SSc';'Gd']); // ここ (画像が表示されない場合があります。そのときは、右クリックで「画像を再読込み or 表示」を試してください) {{ref_image 180807Scilab2ndOrderDelayBodeBoth.png}} これで、連続時間と離散時間のボード線図が重なった。 逆に言えば、サンプル周期を細かくしないと、連続時間の周波数応答とは異なるということか。 !!! 2018年8月14日 続き。 離散時間の伝達関数と連続時間の伝達関数をボード線図に描き、ほぼ同じ軌道だった。 さて、連続時間伝達関数 Gc = 1 --------- 1 + 0.1s を、離散時間伝達関数に変換すると、 Gd = 0.0099502 --------------- -0.9900498 + z これを時間領域の関数に戻す。まず、連続時間伝達関数。Maximaで逆ラプラス変換する。 (%i8) exp1: 1 / (1 + 0.1 * s); (exp1) 1 / (0.1 * s + 1) (%i10) exp2: ilt(exp1, s, t); rat: replaced 0.1 by 1/10 = 0.1 (exp2) 10 * %e^(-10 * t) 10 * exp(-10 * t) となった。 !!! 2018年8月16日 続き。 離散時間伝達関数を時間領域の関数に戻す。離散時間領域の関数か。 、、、と思ったが、以前できなかったのを忘れていた。 MATLABなら、d2c() という関数があるのに、、 !!! 2018年8月21日 続き。 7月26日に ""離散時間伝達関数を差分方程式に変換するにはどうすれば良いだろうか？ と書いた。しかし、、 離散時間伝達関数の逆変換とは、逆z変換そのものだ。馬鹿だ、、 今まで、z変換は全くやっていない。 これを機会にやってみる。 変換する式は、 Gd = 0.0099502 --------------- -0.9900498 + z !!! 2018年8月23日 続き。 z変換の勉強から。 [z変換 (やる夫で学ぶディジタル信号処理)|http://www.ic.is.tohoku.ac.jp/~swk/lecture/yaruodsp/zt.html] を見ながら、、 Z変換の定義は、 ∞ X(z) = Σx[n] z^-n n=0 と定義し、次のように書く。 Z[x[]] = X(z) Z x[n] ---> X(z) そして、次の数列 x[-1] = 0, x[0] = 1, x[1] = 3, x[2] = -2, x[3] = -1, x[4] = 0 のz変換は、、 ∞ -n X(z) = Σx[n] z n=0 -1 -2 -3 = 1 +3z -2z -z となるようだ。コピペした。 * 最初の x[-1] は、定義が n=0 からだから、無視される。 * 次の x[0] = 1 は、そのまま 1 だ。定義からは 1z^-0 だから、z^0 = 1 * x[1] = 3 は、x[n] の部分が 3、z^-n の部分が -1 * x[2] = -2 は、x[n] の部分が 2、z^-n の部分が -2 * x[3] = -1 は、x[n] の部分が -1、z^-n の部分が -3 * x[4] = 0 は、x[n] の部分が 0、z^-n の部分が -4 だから、0z^-4 = 0 なので、式には現れない そうすると、一番単純なz変換は、数列の表現をただ変えただけ？ !!! 2018年9月11日 続き。間が空いてしまった。 z変換の勉強の続き。 [z変換 (やる夫で学ぶディジタル信号処理)|http://www.ic.is.tohoku.ac.jp/~swk/lecture/yaruodsp/zt.html] を見ながら、、 n x[n] = α のz変換は、、 ∞ n -n X(z) = Σα z n=0 ここで、a^r b^r = (a b)^r だから、a = α、b = z^-1 とすれば、、 ∞ -1 n X(z) = Σ(αz ) n=0 ここで、a = 1, r = αz^-1 とすると、等比級数だ。 Σ(ar)^n = a + ar + ar^2 + ... + ar^n-1 両辺をr倍すると、 rΣ(ar)^n = ar + ar^2 + ar^3 + ... + ar^n 最初の式をr倍した式で引くと、、 Σ(ar)^n - rΣ(ar)^n = (a + ar + ar^2 + ... + ar^n-1) - (ar + ar^2 + ar^3 + ... + ar^n) = a - ar^n さらに、Σ(ar)^n について解くと、、 Σ(ar)^n - rΣ(ar)^n = (1 - r)Σ(ar)^n = a - ar^n Σ(ar)^n = (a - ar^n) / (1 - r) さらに、n は ∞ になるので、 lim { (a - ar^n) / (1 - r) } n→∞ = lim { a / (1 - r) } - lim { ar^n / (1 - r) } n→∞ n→∞ ここで、r < 1 ならば、 = a / (1 - r) となる。さらに a = 1, r = αz^-1 だったので、 = 1 / (1 - αz^-1) となった。 !!! 2018年9月12日 z変換の勉強の続き。 x[n] = δ[n] のz変換は、、 ∞ X(x) = Σδ[n]^-n n=0 離散値のデルタ関数は、 x[0] = 1, x[1] = 0, x[2] = 0, x[3] = 0, x[n] = 0 ということは、、 ∞ -n X(z) = Σδ[n] z n=0 = 1z^-0 +0z^-1 +0z^-2 +0z^-3 = 1 となるようだ。 * δ[0]は、1 。定義からは 1z^-0 だから、1 * δ[1]以降は、0。 なお、デルタ関数は、連続時間でラプラス変換・フーリエ変換しても1、離散時間でz変換しても1、だ。 !!! 2018年9月13日 [z変換 (やる夫で学ぶディジタル信号処理)|http://www.ic.is.tohoku.ac.jp/~swk/lecture/yaruodsp/zt.html] での勉強の続き。 ところで、 1 / (1 - αz^-1) は、|αz^-1| < 1 でないと収束しない。 z変換の収束範囲は、|z| > |α| になる。 前に、z変換の定義は、 ∞ X(z) = Σx[n] z^-n n=0 と書いたが、このzは、e^(c+jω) だ。変形すると、e^c e^jω。 オイラーの公式が、 e^(iθ) = cos(θ) + i sin(θ) なので、この式 e^c e^jω は、半径 e^c の円を描く。 なお、ラプラス変換式の s は、s = c + jω である。 !!! 2018年9月15日 [z変換 (やる夫で学ぶディジタル信号処理)|http://www.ic.is.tohoku.ac.jp/~swk/lecture/yaruodsp/zt.html] での勉強の続き。 一気に「線形差分方程式とz変換」に行く。 差分方程式 N M Σ (a_k y[n - k]) = Σ (b_k x[n - k]) k=0 k=0 が題材。 微分方程式ではラプラス変換、差分方程式にはz変換を使う。 *微分方程式: d/dt → s倍 *差分方程式: 1時刻の遅延 演算子「D」を導入する。 D y[n] = y[n - 1] さきほどの差分方程式を書き直す。 N M Σ (a_k D^k y[n]) = Σ (b_k D^k x[n]) k=0 k=0 突然だが、離散時間フーリエ逆変換の式。 1 π jωn f[n] = --- ∫ F(ω) e dω 2π -π 離散時間フーリエ逆変換で置き換え、総和と積分を入れ替える。 1 N π 1 M π --- Σ (a_k D^k ∫ Y(ω) e^jωn dω) = --- Σ (b_k D^k ∫ X(ω) e^jωn dω) 2π k=0 -π 2π k=0 -π 1 π N 1 π M --- ∫ Σ (a_k D^k Y(ω) e^jωn dω) = --- ∫ Σ (b_k D^k X(ω) e^jωn dω) 2π -π k=0 2π -π k=0 遅延演算子 D を e^jωn に作用させると、 jωn jω(n - 1) D e = e -jω jωn = e e なので、D = e^-jω となる。 すると、 1 π N 1 π M --- ∫ [Σ (a_k (e^-jω)^k Y(ω) ] e^jωn dω) = --- ∫ [Σ (b_k (e^-jω)^k X(ω)] e^jωn dω) 2π -π k=0 2π -π k=0 総和の外側は同じなので、取り除く。 N M Σ (a_k (e^-jω)^k Y(ω)) = Σ (b_k (e^-jω)^k X(ω)) k=0 k=0 逆z変換で考えると、e^-jω が e^-(c+jω) になる。 (c+jω)n (c+jω)(n - 1) D e = e -(c+jω) (c+jω)n = e e なので、D = e^-(c+jω) となる。 すると、 N M Σ (a_k (e^-(c+jω))^k Y(ω)) = Σ (b_k (e^-(c+jω))^k X(ω)) k=0 k=0 z = e^(c+jω) なので、 N M Σ (a_k z^-k Y(z)) = Σ (b_k z^-k X(z)) k=0 k=0 Y(z) について解くと、 M Σ (b_k z^-k) k=0 Y(z) = --------------- X(z) N Σ (a_k z^-k) k=0 伝達関数 H(z) = Y(z) / X(z) とすれば、 M Σ (b_k z^-k) k=0 H(z) = --------------- N Σ (a_k z^-k) k=0 疲れた。 !!! 2018年10月16日 [z変換 (やる夫で学ぶディジタル信号処理)|http://www.ic.is.tohoku.ac.jp/~swk/lecture/yaruodsp/zt.html] をもとに変換をしてみる。 Gd = 0.0099502 --------------- -0.9900498 + z 先日、逆z変換について調べたが、実際はz変換表を用いる。 使えそうな変換は、 , ,時間領域,z領域 ,指数関数,a^n u0(t),1/(1 - az^-1) だが、上の式のzには冪乗の係数がない。 Gd = 0.0099502 --------------- -0.9900498 + z どうしようかと考えたが、Gd の分子分母に z^-1 をかけてみよう。 0.0099502 z^-1 --------------- * ------ -0.9900498 + z z^-1 0.0099502 * z^-1 -------------------------------- -0.9900498 * z^-1 + z^1 * z^-1 0.0099502 * z^-1 ------------------------------- -0.9900498 * z^-1 + z^(1-1) 0.0099502 * z^-1 ----------------------- -0.9900498 * z^-1 + 1 これで z^-1 ができた。 {{include_html amazon_match}} {{include_html adsense}}