Skip to main content

線形離散時間システムモデルのシミュレーション

システムモデル​

2次の線形離散時間システムを考える：

xt=−a1xt−1−a2xt−2+b1ut−1+b2ut−2x_t=-a_1x_{t-1}-a_2x_{t-2}+b_1u_{t-1}+b_2u_{t-2}

Z変換による表現​

(1+a1z−1+a2z−2)X(z)=(b1z−1+b2z−2)U(z)(1+a_1z^{-1}+a_2z^{-2})X(z)=(b_1z^{-1}+b_2z^{-2})U(z)

伝達関数​

H(z)=X(z)U(z)=b1z−1+b2z−21+a1z−1+a2z−2=z−1+0.5z−21−1.5z−1+0.7z−2=z+0.5z2−1.5z+0.7H(z)=\frac{X(z)}{U(z)}=\frac{b_1z^{-1}+b_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}}=\frac{z^{-1}+0.5z^{-2}}{1-1.5z^{-1}+0.7z^{-2}}=\frac{z+0.5}{z^2-1.5z+0.7}

部分分数展開​

H(z)=z−1(1+0.5z−1)(1−1.5+j0.552z−1)(1−1.5−j0.552z−1)H(z)=\frac{ z^{-1}(1+0.5z^{-1}) }{ (1-\frac{1.5+j\sqrt{0.55}}{2}z^{-1})(1-\frac{1.5-j\sqrt{0.55}}{2}z^{-1})}
MATLABでの数値計算

MATLABでディジタル信号を扱う時に、そのまま差分方程式を使える。Z変換などが要らない。

MATLABシミュレーションコード​

初期設定​

clear
close all
set(0,'DefaultAxesFontName','Times New Roman')
set(0,'DefaultAxesFontSize',12)

%--------------------------------------------------------------------
%  同定ゼミ宿題pdf p.13 線形離散時間システムモデルのシミュレーション
%                                                       4/13 by YANG
%---------------------------------------------------------------------

パラメータ設定​

% 初期パラメータ
Tmax = 10;           % シミュレーション時間
samp = 0.01;         % サンプリング周期
n = Tmax/samp;       % サンプル数
u = rand(1,n+1);     % 入力（一様乱数）

NSR = 0.0;           % 雑音信号比（Noise Signal Ratio）
% NSR 0, 0.1, 0.2で比較可能

% システムパラメータ
a1 = -1.5;
a2 = 0.7;
b1 = 1.0;
b2 = 0.5;

配列初期化​

t = zeros(1, n+1);   % 時間配列
x = zeros(1, n+1);   % 真の出力
y = zeros(1, n+1);   % 観測値（雑音付き）

システム応答計算​

% 真の出力の計算
for k = 3:n+1
    t(k) = (k-1) * samp;
    x(k) = -a1*x(k-1) - a2*x(k-2) + b1*u(k-1) + b2*u(k-2);
end

観測雑音の追加​

% ゼロ平均のガウス分布乱数（正規分布）
mu = 0;                    % 平均値
sigma = NSR * std(x);      % 標準偏差 NSR = sigma/std(x)
v = sigma * randn(1,n+1) + mu;

% 観測値の計算
for k = 3:n+1
    y(k) = x(k) + v(k);
end

結果のプロット​

subplot(2,2,1); plot(t,u); title('入力信号 u(t)');
subplot(2,2,2); plot(t,x); title('真の出力 x(t)');
subplot(2,2,3); plot(t,y); title('観測出力 y(t)');
subplot(2,2,4); plot(t,v); title('観測雑音 v(t)');

シミュレーション結果の解析​

雑音信号比（NSR）の影響​

  • NSR = 0.0: 雑音なし、理想的な条件
  • NSR = 0.1: 10%の雑音レベル
  • NSR = 0.2: 20%の雑音レベル
雑音信号比とは

NSR（Noise Signal Ratio）は観測雑音の標準偏差と真の出力信号の標準偏差との比である： NSR=σvσxNSR = \frac{\sigma_v}{\sigma_x}

システムの安定性​

システムの極は次の方程式の解： z2−1.5z+0.7=0z^2 - 1.5z + 0.7 = 0

複素共役極： z=1.5±j0.552z = \frac{1.5 \pm j\sqrt{0.55}}{2}

極の大きさ： ∣z∣=0.7≈0.837<1|z| = \sqrt{0.7} \approx 0.837 < 1

安定性の確認

すべての極が単位円内にあるため、システムは安定である。

応用例​

このシミュレーションは以下の用途に応用できる：

  1. システム同定: パラメータa1,a2,b1,b2a_1, a_2, b_1, b_2の推定
  2. 雑音の影響評価: 異なるNSRレベルでの性能比較
  3. 制御器設計: システム特性の理解
  4. フィルタ設計: 雑音除去手法の検討

実装のポイント​

初期値の設定​

  • 配列の初期化でzerosを使用
  • ループはk=3から開始（過去2サンプルが必要）

雑音生成​

  • randnでガウス分布乱数生成
  • std(x)で真の出力の標準偏差を計算
  • NSRによって雑音レベルを調整

プロット設定​

  • subplotで複数のグラフを同時表示
  • フォント設定で可読性を向上
  • 日本語タイトルで理解しやすさを向上

まとめ​

本シミュレーションでは：

  1. 差分方程式: 離散時間システムの直接実装
  2. 雑音モデル: 実際的な観測条件のシミュレーション
  3. 可視化: 入力、出力、雑音の関係の理解
  4. パラメータ影響: NSRによる性能評価

これらの要素は、システム同定や制御系設計における基本的な数値実験の基盤となる。