跳到主要内容

逐次最小二乗法によるシステム同定シミュレーション

この文書では、逐次最小二乗法（Recursive Least Squares, RLS）を用いたシステム同定のシミュレーションについて詳しく説明します。2次線形システムのパラメータ推定、時変システム、および補助変数法を含む包括的な分析を行います。

3.5.1 同定対象システム​

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

xt=−a1xt−1−a2xt−2+b1ut−1+b2ut−2yt=xt+vt\begin{aligned} x_{t} &=-a_1x_{t-1}-a_2x_{t-2}+b_1u_{t-1}+b_2u_{t-2} \\ y_t &=x_t+v_t \end{aligned}

ここで、真のパラメータ値は：

a1=−1.5,a2=0.7,b1=1.0,b2=0.5a_1=-1.5,\quad a_2=0.7,\quad b_1=1.0, \quad b_2=0.5

システム特性​

  • 入力信号: [0,1][0,1] の範囲で一様分布する乱数
  • 観測雑音: ゼロ平均のガウス分布乱数
  • 無相関性: 雑音は入力信号と無相関
雑音信号比（NSR）

雑音のレベルは雑音信号比（Noise Signal Ratio, NSRNSR）で評価されます。

NSR=σvσxNSR = \frac{\sigma_v}{\sigma_x}

ここで、σv\sigma_vは雑音vtv_tの標準偏差、σx\sigma_xは真の出力信号xtx_tの標準偏差です。

3.5.2 逐次最小二乗法によるパラメータ推定​

数学的定式化​

観測方程式：

yt=−a1yt−1−a2yt−2+b1ut−1+b2ut−2+rty_t=-a_1y_{t-1}-a_2y_{t-2}+b_1u_{t-1}+b_2u_{t-2}+r_t

回帰ベクトル：

ztT=[−yt−1,−yt−2,ut−1,ut−2]z_t^T=[-y_{t-1},-y_{t-2},u_{t-1},u_{t-2}]

パラメータベクトル：

θT=[a1,a2,b1,b2]\theta^T=[a_1,a_2,b_1,b_2]

シミュレーション設定​

  • データ長: N=3000N = 3000（十分大きな整数）
  • NSR レベル: 0%, 10%, 20%
  • 目的: 各NSRレベルでのパラメータ推定性能の比較

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

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

% パラメータ設定
N = 3000;
rand('seed',1);

% 真のパラメータ
a1 = -1.5; a2 = 0.7; b1 = 1.0; b2 = 0.5;

% 初期化
t = zeros(1, N+1);
u = rand(1, N+1);
x = zeros(1, N+1);
theta1_hat = zeros(4,N+1); % NSR=0
theta2_hat = zeros(4,N+1); % NSR=0.1
theta3_hat = zeros(4,N+1); % NSR=0.2
alpha = 1000;
rho = 0.95*ones(1, N+1);

% システム応答生成
for k=3:N+1
    t(k) = k-1;
    x(k) = -a1*x(k-1)-a2*x(k-2)+b1*u(k-1)+b2*u(k-2);
end

% 各NSRレベルでの推定
for NSR = 0:0.1:0.2
    % ガウス雑音生成
    mu = 0;
    sigma = NSR*std(x);
    v = sigma*randn(1,N+1) + mu;
    
    % 観測信号
    y = x + v;
    
    % 逐次最小二乗法
    for k=3:N+1
        z(:,k) = [-y(k-1),-y(k-2),u(k-1),u(k-2)]';
        
        % 忘却係数更新
        rho(k) = (1-0.01)*rho(k-1)+0.01;
        
        % 共分散行列更新
        P(:,:,k) = (P(:,:,k-1) - (P(:,:,k-1)*psi(:,k)*phi(:,k)'*P(:,:,k-1))/(rho(k) + phi(:,k)'*P(:,:,k-1)*psi(:,k)))/rho(k);
        
        % ゲインベクトル
        L(:,k) = (P(:,:,k-1)*psi(:,k))/(rho(k) + phi(:,k)'*P(:,:,k-1)*psi(:,k));
        
        % 予測誤差
        E(k) = y(k) - phi(:,k)'*theta_hat(:,k-1);
        
        % パラメータ更新
        theta_hat(:,k) = theta_hat(:,k-1) + L(:,k)*E(k);
    end
end

3.5.3 時変システムのパラメータ推定​

時変システムモデル​

2次の時変システム：

xt=−a1,txt−1−a2,txt−2+b1,tut−1+b2,tut−2yt=xt+vt\begin{aligned} &x_t = -a_{1,t}x_{t-1}-a_{2,t}x_{t-2}+b_{1,t}u_{t-1}+b_{2,t}u_{t-2} \\ &y_t = x_t+v_t \end{aligned}

時変パラメータ：

a1,t=−1.0+0.5sin⁡(2πt/2000)a2,t=0.7b1,t=1.0b2,t=0.5+0.4cos⁡(2πt/2000)\begin{aligned} a_{1,t} &= -1.0+0.5\sin(2\pi t/2000) \\ a_{2,t} &= 0.7 \\ b_{1,t} &= 1.0 \\ b_{2,t} &= 0.5+0.4\cos(2\pi t/2000) \end{aligned}

時変パラメータ推定の特徴​

  1. 追従性能: 忘却係数により過去のデータの影響を減衰
  2. 適応能力: パラメータ変化に対する追従速度
  3. 雑音感度: 時変システムでは雑音の影響がより顕著

MATLABシミュレーション（時変システム）​

% 時変パラメータ生成
for k=1:N+1
    a1(k) = -1.0 + 0.5*sin(2*pi*k/2000);
    a2(k) = 0.7;
    b1(k) = 1.0;
    b2(k) = 0.5 + 0.4*cos(2*pi*k/2000);
end

% 時変システム応答
for k=3:N+1
    x(k) = -a1(k-1)*x(k-1)-a2(k-2)*x(k-2)+b1(k-1)*u(k-1)+b2(k-2)*u(k-2);
end

% 逐次推定（忘却係数 = 0.98）
for k=3:N+1
    rho(k) = 0.98; % 固定忘却係数
    % [推定アルゴリズムは上記と同様]
end

3.5.4 逐次補助変数法によるパラメータ推定​

補助変数法は、雑音と回帰ベクトルの相関を除去することでより正確な推定を実現します。

方法1: 遅延された入力を補助変数とする​

補助変数ベクトル：

mt=[ut−1−d,ut−2−d,ut−1,ut−2]Tm_t = [u_{t-1-d}, u_{t-2-d}, u_{t-1}, u_{t-2}]^T

ここで、d≥nd \geq n（nnはシステム次数）の遅延を設けます。

MATLABコード（遅延入力）​

d = 2; % 遅延パラメータ

% 遅延された入力の設定
for k=3+d:N+1
    m(:,k) = [u(k-1-d), u(k-2-d), u(k-1), u(k-2)]';
end

% 補助変数法推定
for k=3:N+1
    phi(:,k) = z(:,k);     % 回帰ベクトル
    psi(:,k) = m(:,k);     % 補助変数ベクトル
    
    % 共分散行列更新
    P(:,:,k) = (P(:,:,k-1) - (P(:,:,k-1)*psi(:,k)*phi(:,k)'*P(:,:,k-1))/(rho(k) + phi(:,k)'*P(:,:,k-1)*psi(:,k)))/rho(k);
    
    % パラメータ更新
    L(:,k) = (P(:,:,k-1)*psi(:,k))/(rho(k) + phi(:,k)'*P(:,:,k-1)*psi(:,k));
    E(k) = y(k) - phi(:,k)'*theta_hat(:,k-1);
    theta_hat(:,k) = theta_hat(:,k-1) + L(:,k)*E(k);
end

方法2: 遅延された出力を補助変数とする​

補助変数ベクトル：

mt=[yt−1−d,yt−2−d,ut−1,ut−2]Tm_t = [y_{t-1-d}, y_{t-2-d}, u_{t-1}, u_{t-2}]^T

方法3: 推定された出力を補助変数とする​

この方法では、推定されたパラメータから出力を予測し、それを補助変数として使用します。

% 安定性チェック
p = [1 theta_hat(1,k-1) theta_hat(2,k-1)];
s = roots(p);

if abs(s(1)) < 1 && abs(s(2)) < 1
    % 安定な場合
    bar_a1(k-1) = theta_hat(1,k-1);
    bar_a2(k-1) = theta_hat(2,k-1);
else
    % 不安定な場合は調整
    bar_a1(k-1) = bar_a1(k-2) + 0.5*(theta_hat(1,k-1) - bar_a1(k-2));
    bar_a2(k-1) = bar_a2(k-2) + 0.5*(theta_hat(2,k-1) - bar_a2(k-2));
end

% 推定出力計算
hat_x(k) = -bar_a1(k-1)*hat_x(k-1) - bar_a2(k-1)*hat_x(k-2) + theta_hat(3,k-1)*u(k-1) + theta_hat(4,k-1)*u(k-2);

% 補助変数
m(:,k) = [-hat_x(k-1), -hat_x(k-2), u(k-1), u(k-2)]';

性能比較と考察​

手法比較表​

手法利点欠点適用場面
逐次最小二乗法計算が簡単、リアルタイム実装容易雑音に対する感度が高い低雑音環境
補助変数法（遅延入力）雑音の影響を軽減遅延による情報損失入力信号が豊富
補助変数法（遅延出力）バランスの取れた性能実装がやや複雑一般的な応用
補助変数法（推定出力）最も高い推定精度安定性チェックが必要高精度が要求される場面
推奨事項
  • 低雑音環境: 逐次最小二乗法で十分
  • 中程度雑音: 補助変数法（遅延出力）を推奨
  • 高雑音環境: 補助変数法（推定出力）を使用
  • 時変システム: 忘却係数の適切な調整が重要

忘却係数の影響​

忘却係数ρ\rhoの設定は推定性能に大きく影響します：

  • ρ=1.0\rho = 1.0: 理想的だが数値的に不安定
  • ρ=0.98\rho = 0.98: 時変システムに適している
  • ρ=0.95\rho = 0.95: より高い追従性能だが雑音に敏感

実用的考慮事項​

  1. 初期値設定: P(0)=αIP(0) = \alpha I（α\alphaは大きな正数）
  2. 数値安定性: 定期的な行列の条件数チェック
  3. 計算量: リアルタイム実装時のCPU負荷考慮
  4. メモリ使用量: 長時間運転時のメモリ管理

まとめ​

逐次最小二乗法とその拡張である補助変数法は、システム同定において強力な手法です。雑音レベル、システムの時変性、計算資源などを総合的に考慮して最適な手法を選択することが重要です。

実際の応用では、複数の手法を組み合わせたハイブリッドアプローチや、適応的に手法を切り替える戦略も有効です。