5 posts tagged with "MATLAB"

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

June 1, 2024 · 8 min read

R&D Researcher, Control & Robotics

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

3.5.1 同定対象システム

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

\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}

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

a_1=-1.5,\quad a_2=0.7,\quad b_1=1.0, \quad b_2=0.5

システム特性

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

雑音信号比（NSR）

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

$NSR = \frac{\sigma_v}{\sigma_x}$

ここで、 $\sigma_v$ は雑音 $v_t$ の標準偏差、 $\sigma_x$ は真の出力信号 $x_t$ の標準偏差です。

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

数学的定式化

観測方程式：

y_t=-a_1y_{t-1}-a_2y_{t-2}+b_1u_{t-1}+b_2u_{t-2}+r_t

回帰ベクトル：

z_t^T=[-y_{t-1},-y_{t-2},u_{t-1},u_{t-2}]

パラメータベクトル：

\theta^T=[a_1,a_2,b_1,b_2]

シミュレーション設定

データ長: $N = 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次の時変システム：

\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}

時変パラメータ：

\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}

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

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

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: 遅延された入力を補助変数とする

補助変数ベクトル：

m_t = [u_{t-1-d}, u_{t-2-d}, u_{t-1}, u_{t-2}]^T

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

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: 遅延された出力を補助変数とする

補助変数ベクトル：

m_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$ の設定は推定性能に大きく影響します：

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

実用的考慮事項

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

まとめ

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

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

PID制御器のシミュレーション

May 26, 2024 · 9 min read

Guangze Yang

R&D Researcher, Control & Robotics

制御対象

4次の伝達関数で表される制御対象を考える：

G(s)=\frac{12s+8}{20s^4+112s^3+147s^2+62s+8}

状態空間表現への変換

num = [0 0 12 8];
den = [20 113 147 62 8];
sysg = tf(num, den);                    % 伝達関数
[Ao,Bo,Co,Do] = tf2ss(num,den);        % 状態空間表現への変換

PID制御器

PID制御ブロック図

離散時間での表現（ $k \to n$ ）：

e(n)=r(n)-y(n-1)

偏差の定義

偏差 = 目標値 - 一個前の制御量

U(z)=E(z)C(z)

PID制御器伝達関数

C(s)=c_p+c_i\frac{1}{s}+c_d\frac{s}{\gamma s+1}, \quad \gamma=4\times T

gamma = 4*T;  % γ（微分項のフィルタ時定数）
cp = 6;       % 比例ゲイン
ci = 1;       % 積分ゲイン  
cd = 7;       % 微分ゲイン

数値積分手法による離散化

1. 前進差分

PID制御器のZ変換表現

C(z)=c_p+c_i\frac{Tz^{-1}}{1-z^{-1}}+c_d\frac{1-z^{-1}}{\gamma+(T-\gamma)z^{-1}}

P-I-D分解実装（推奨方法）

比例制御（P）：

u_p(n)=c_pe(n)

積分制御（I）：

u_i(n)-u_i(n-1)=c_iTe(n-1) \Rightarrow u_i(n)=u_i(n-1)+c_iTe(n-1)

積分制御の意味

偏差の累積で制御する

微分制御（D）：

\gamma u_d(n)+(T-\gamma)u_d(n-1)=c_d(e(n)-e(n-1))

\Rightarrow u_d(n)=\frac{c_d}{\gamma}\left(e(n)-e(n-1)-\frac{T-\gamma}{\gamma}u_d(n-1)\right)

微分制御の意味

偏差の変化で制御するつもりだが、初期状態には変化がない。微分の $s$ × ローパスフィルタ $\frac{1}{\gamma s+1}$

状態空間表現（前進差分）

x(n+1)=x(n)+T(Ax(n)+bu(n))

y(n)=Cx(n)+Du(n)

2. 後退差分

PID制御器のZ変換表現

C(z)=c_p+c_i\frac{T}{1-z^{-1}}+c_d\frac{1-z^{-1}}{\gamma(1-z^{-1})+T}

P-I-D分解実装

比例制御（P）：

u_p(n)=c_pe(n)

積分制御（I）：

u_i(n)-u_i(n-1)=c_iTe(n) \Rightarrow u_i(n)=c_iTe(n)+u_i(n-1)

微分制御（D）：

(\gamma+T)u_d(n)-\gamma u_d(n-1)=c_d(e(n)-e(n-1))

\Rightarrow u_d(n)=\frac{c_d(e(n)-e(n-1)) +\gamma u_d(n-1)}{(\gamma+T)}

状態空間表現（後退差分）

x(n)-x(n-1)=T(Ax(n)+bu(n)) \Rightarrow x(n)=(I-TA)^{-1}(x(n-1)+Tbu(n))

y(n)=Cx(n)+Du(n)

3. 双一次変換法

PID制御器のZ変換表現

\begin{aligned} C(z)&=c_p+c_i\frac{T}{2}\frac{1+z^{-1}}{1-z^{-1}}+c_d\frac{1}{\gamma+\frac{T(1+z^{-1})}{2(1-z^{-1})}}\\ &=c_p+\frac{Tc_i}{2}\frac{1+z^{-1}}{1-z^{-1}}+2c_d\frac{1-z^{-1}}{2\gamma+T+(T-2\gamma)z^{-1}} \end{aligned}

P-I-D分解実装

比例制御（P）：

u_p(n)=c_pe(n)

積分制御（I）：

u_i(n)-u_i(n-1)=\frac{c_iT}{2}(e(n)+e(n-1)) \Rightarrow u_i(n)=\frac{c_iT}{2}(e(n)+e(n-1))+u_i(n-1)

微分制御（D）：

(2\gamma+T)u_d(n)+(T-2\gamma)u_d(n-1)=2c_d(e(n)-e(n-1))

\Rightarrow u_d(n)= \frac{2c_d(e(n)-e(n-1))+(2\gamma-T)u_d(n-1)}{2\gamma+T}

状態空間表現（双一次変換）

\frac{2}{T}(x(n)-x(n-1))=A(x(n)+x(n-1))+b(u(n)+u(n-1))

\Rightarrow x(n)=(2I-TA)^{-1}[(2I+AT)x(n-1)+bT(u(n)+u(n-1))]

y(n)=Cx(n)+Du(n)

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

プログラミング上の注意点

行列計算の注意事項

行列計算する時、行列の次元に注意

inv()  % 逆行列
eye()  % 単位行列I
./     % 行列除算
.*     % 行列要素ごとの乗算

最も注意すべきなのは、括弧と正負の符号である。理由は、コードの数式が読みづらい。

解決策：

コードをLaTeX公式に入れて、確認する
目と手でもう一回計算して確認する

制御ループの実装

for k=1:N+1
    y(k) = C*x(k) + D*u(k)     % 出力計算
    e(k) = r(k) - y(k)         % 偏差計算
    
    % 離散方法
    up(k) = ...                % 比例制御
    ui(k) = ...                % 積分制御  
    ud(k) = ...                % 微分制御
    u(k) = up(k) + ui(k) + ud(k)
    
    % 前進差分で状態更新
    x(k+1) = x(k) + T*(A*x(k) + B*u(k))
end

制御順序の重要性

実際の場合、センサからの観測値をもらうので、 $y(t)$ が先の方である
しかし、 $y(t)$ を先にすると、 $y=Cx+Du$ における $u$ がずっと0で計算している。普通に $D=0$ なので影響がない
また、現在の入力 $u(t)$ を求めたい。 $u(t)$ を求めるために、現在の出力（観察値） $y(t)$ が必要である。出力を得るために、現在まだ求めていない入力を使うべきではない

完全なシミュレーションコード

PID制御器比較

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

%--------------------------------------------------------------------
%  同定ゼミ宿題 PID制御器のシミュレーション
%                                                5/26 by YANG
%---------------------------------------------------------------------

Tmax = 30;    % シミュレーション時間
T = 0.05;     % サンプリング周期
n = Tmax/T;   % サンプル数
gamma = 4*T;  % フィルタ時定数

% PIDパラメータ
cp = 6;       % 比例ゲイン
ci = 1;       % 積分ゲイン
cd = 7;       % 微分ゲイン

% 制御対象
num = [0 0 12 8];
den = [20 113 147 62 8];
sysg = tf(num, den);
[Ao,Bo,Co,Do] = tf2ss(num,den);

% 配列初期化
t = zeros(1, n+1);
r = ones(1, n+1);      % 目標値（ステップ入力）

% 前進差分用配列
upF = zeros(1, n+1);   % P制御入力
uiF = zeros(1, n+1);   % I制御入力
udF = zeros(1, n+1);   % D制御入力
uF = zeros(1, n+1);    % 総制御入力
xF = zeros(4, n+1);    % 状態変数
yF = zeros(1, n+1);    % 出力
eF = zeros(1, n+1);    % 偏差

% 前進差分によるPID制御
for k = 2:n+1
    t(k) = (k-1)*T;
    yF(k) = Co*xF(:,k) + Do*uF(k);
    eF(k) = r(k) - yF(k);
    
    % PID制御計算
    upF(k) = cp*eF(k);
    uiF(k) = uiF(k-1) + ci*T*eF(k-1);
    udF(k) = (cd*(eF(k)-eF(k-1)) - (T-gamma)*udF(k-1))/gamma;
    uF(k) = upF(k) + uiF(k) + udF(k);
    
    % 状態更新
    xF(:,k+1) = xF(:,k) + T*(Ao*xF(:,k) + Bo*uF(k));
end

% 他の手法（後退差分、双一次変換）も同様に実装...

P-I-D個別効果比較

% P, I, D制御の個別効果を比較するシミュレーション
for k = 2:n+1
    t(k) = (k-1)*T;
    
    % P制御のみ
    ypF(k) = Co*xpF(:,k) + Do*up(k);
    epF(k) = r(k) - ypF(k);
    up(k) = cp*epF(k);
    xpF(:,k+1) = xpF(:,k) + T*(Ao*xpF(:,k) + Bo*up(k));
    
    % I制御のみ
    yiF(k) = Co*xiF(:,k) + Do*ui(k);
    eiF(k) = r(k) - yiF(k);
    ui(k) = ui(k-1) + ci*T*eiF(k-1);
    xiF(:,k+1) = xiF(:,k) + T*(Ao*xiF(:,k) + Bo*ui(k));
    
    % D制御のみ
    ydF(k) = Co*xdF(:,k) + Do*ud(k);
    edF(k) = r(k) - ydF(k);
    ud(k) = (cd*(edF(k)-edF(k-1)) - (T-gamma)*ud(k-1))/gamma;
    xdF(:,k+1) = xdF(:,k) + T*(Ao*xdF(:,k) + Bo*ud(k));
end

結果可視化

% 制御入力の比較
subplot(2,1,1); 
plot(t,uF,t,uB,t,uD); 
title('制御入力 u(t)');
legend('前進差分','後退差分','双一次変換'); 
ylabel('u(t)'); 
xlabel('Time[sec]'); 
grid on

% 出力応答の比較
subplot(2,1,2);
plot(t,yF,t,yB,t,yD,t,r,'--'); 
title('目標値r(t)と出力y(t)');
legend('前進差分','後退差分','双一次変換','目標値'); 
ylabel('y(t)'); 
xlabel('Time[sec]'); 
grid on

性能比較と考察

各手法の特性

手法	安定性	精度	計算負荷	実装容易さ
前進差分	条件付き	低	低	高
後退差分	良好	中	中	中
双一次変換	優秀	高	中	中

PID制御の各要素の役割

P制御（比例制御）

効果: 偏差に比例した制御
特徴: 応答速度の向上、定常偏差の残存
調整: ゲインを大きくすると応答が速くなるが、振動しやすくなる

I制御（積分制御）

効果: 偏差の累積を解消
特徴: 定常偏差の除去、応答の遅れ
調整: ゲインを大きくすると定常偏差は早く除去されるが、オーバーシュートが増加

D制御（微分制御）

効果: 偏差の変化率に基づく制御
特徴: 応答の改善、雑音の増幅
調整: ゲインを大きくすると安定性は向上するが、雑音に敏感になる

実装上のポイント

ゼロオーダーホールドの適用: 因果関係の保持
フォント設定: 可読性の向上
軸ラベルと単位: グラフの明確化
初期値設定: 適切な初期条件の設定

まとめ

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

離散化手法: 3つの数値積分手法による比較
PID制御: 各制御要素の個別効果と統合効果
実装技法: MATLABによる効率的なプログラミング手法
性能評価: 安定性、精度、実用性の観点からの評価

これらの知識は、実際の制御系設計における重要な基盤となる。特に、離散化手法の選択は制御性能に大きく影響するため、システムの特性と要求性能に応じた適切な選択が重要である。

数値積分手法の比較シミュレーション

April 21, 2024 · 5 min read

Guangze Yang

R&D Researcher, Control & Robotics

前進差分、後退差分、双一次変換をそれぞれ用いて、例題1.2.1の電気回路のステップ応答を数値計算で求め、数値計算の結果と真の応答をプロットし、比較する。

電気回路モデル

RL回路の微分方程式：

L\frac{di(t)}{dt} + Ri(t) = u(t)

真の解（ステップ応答）：

i(t) = \frac{1}{R}(1 - e^{-\frac{R}{L}t})

数値積分手法

1. 前進差分（Forward Difference）

微分の近似：

sf(n) = \frac{f(n+1)-f(n)}{T}

よって：

L\frac{i(n+1)-i(n)}{T} + Ri(n) = u(n)

整理すると：

i(n+1) = \frac{1}{L}(Tu(n) + (L-RT)i(n))

時間をシフトして：

i(n) = \frac{1}{L}(Tu(n-1) + (L-RT)i(n-1))

2. 後退差分（Backward Difference）

微分の近似：

sf(n) = \frac{f(n)-f(n-1)}{T}

よって：

L\frac{i(n)-i(n-1)}{T} + Ri(n) = u(n)

整理すると：

i(n) = \frac{1}{RT+L}(Tu(n) + Li(n-1))

3. 双一次変換（Bilinear Transform）

台形積分による近似：

sf(n) = \frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}f(n)

微分方程式に適用：

L\frac{2}{T}(i(n)-i(n-1)) + R(i(n)+i(n-1)) = u(n)+u(n-1)

整理すると：

i(n) = \frac{1}{2L+RT} \left( Tu(n)+Tu(n-1) + (2L-RT)i(n-1) \right)

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

初期設定

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

%--------------------------------------------------------------------
%  同定ゼミ宿題 電気回路応答のシミュレーション
%                                                4/21 by YANG
%---------------------------------------------------------------------

パラメータ設定

Tmax = 10;           % 最大時間
samp = 0.5;          % サンプリング周期
T = samp;
n = Tmax/samp;       % サンプル数
u = ones(1,n+1);     % ステップ入力

% 回路パラメータ
R = 1;               % 抵抗値
L = 1;               % インダクタンス値

配列初期化

% 数列
t  = zeros(1, n+1);  % 時間配列
i  = zeros(1, n+1);  % 真の信号
i1 = zeros(1, n+1);  % 前進差分
i2 = zeros(1, n+1);  % 後退差分
i3 = zeros(1, n+1);  % 双一次変換

各手法による数値計算

真の解析解

% 真の信号
for k = 2:n+1
    t(k) = (k-1) * samp;
    i(k) = (1 - exp((-t(k)*R)/L)) / R;
end

前進差分

% 前進差分
for k = 2:n+1
    t(k) = (k-1) * samp;
    i1(k) = (samp*u(k-1) + (L-R*samp)*i1(k-1)) / L;
end

後退差分

% 後退差分
for k = 2:n+1
    t(k) = (k-1) * samp;
    i2(k) = (T*u(k) + L*i2(k-1)) / (R*samp+L);
end

双一次変換

% 双一次変換
for k = 2:n+1
    t(k) = (k-1) * samp;
    i3(k) = (T*u(k) + T*u(k-1) + (2*L-R*T)*i3(k-1)) / (2*L+R*samp);
end

結果のプロット

% 比較プロット
subplot(3,1,1); 
plot(t,i,t,i1); 
title('前進差分');
legend('真の解','前進差分');  
xlabel('t'); 
ylabel('i(t)');
grid on

subplot(3,1,2); 
plot(t,i,t,i2); 
title('後退差分');
legend('真の解','後退差分');      
xlabel('t'); 
ylabel('i(t)');
grid on

subplot(3,1,3); 
plot(t,i,t,i3); 
title('双一次変換');
legend('真の解','双一次変換');  
xlabel('t'); 
ylabel('i(t)');
grid on

数値積分手法の特性比較

1. 前進差分の特性

特徴：

明示的（explicit）手法
計算が簡単
安定性に制限がある

安定性条件： $|1 - \frac{RT}{L}| \leq 1$

これより： $T \leq \frac{2L}{R}$

安定性の制限

前進差分は大きなサンプリング周期で不安定になる可能性がある。

2. 後退差分の特性

特徴：

暗示的（implicit）手法
無条件安定
計算がやや複雑

安定性： 分母が $RT + L > 0$ で常に正のため、無条件安定。

安定性の優位性

後退差分は大きなサンプリング周期でも安定である。

3. 双一次変換の特性

特徴：

台形積分による近似
優れた周波数特性
バランスの取れた精度

安定性： 連続時間システムが安定なら、離散化後も安定性が保持される。

周波数特性の保持

双一次変換は周波数応答特性を最もよく保持する。

精度と安定性の評価

サンプリング周期の影響

手法	精度	安定性	計算コスト
前進差分	低〜中	条件付き	低
後退差分	中	無条件	中
双一次変換	高	無条件	中

適用場面

前進差分

高速な計算が必要
小さなサンプリング周期
リアルタイム処理

後退差分

安定性が重要
大きなサンプリング周期
ロバスト性が要求される場合

双一次変換

高精度が必要
周波数特性の保持が重要
制御系設計

まとめ

本シミュレーションにより以下のことが確認できる：

精度比較: 双一次変換が最も高精度
安定性: 後退差分と双一次変換が優秀
計算効率: 前進差分が最も簡単
実用性: 双一次変換が総合的に優れている

手法選択の指針

高精度が必要 → 双一次変換
安定性重視 → 後退差分
高速計算 → 前進差分（安定性条件下）

これらの数値積分手法の理解は、システム同定や制御系設計において重要な基礎知識となる。

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

April 13, 2024 · 4 min read

Guangze Yang

R&D Researcher, Control & Robotics

システムモデル

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

x_t=-a_1x_{t-1}-a_2x_{t-2}+b_1u_{t-1}+b_2u_{t-2}

Z変換による表現

(1+a_1z^{-1}+a_2z^{-2})X(z)=(b_1z^{-1}+b_2z^{-2})U(z)

伝達関数

H(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)=\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 = \frac{\sigma_v}{\sigma_x}$

システムの安定性

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

複素共役極： $z = \frac{1.5 \pm j\sqrt{0.55}}{2}$

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

安定性の確認

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

応用例

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

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

実装のポイント

初期値の設定

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

雑音生成

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

プロット設定

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

まとめ

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

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

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

VSCodeによるMATLABの使用

January 1, 2024 · 4 min read

Guangze Yang

R&D Researcher, Control & Robotics

MATLABは、科学技術計算やデータ分析、アルゴリズム開発などに利用される強力なツールですが、専用のMATLABエディタ以外で作業したい場合もあります。VSCodeは、とても人気なエディタであり、たくさんの開発環境に適用できます。例え、MATLAB、LaTex, Python, Cなどがあります。VSCode (Visual Studio Code) を使用してMATLABコードを編集、実行する方法について解説します。

1. 必要な準備

VSCodeでMATLABを使用するためには、以下の手順を行う必要があります。

必要なソフトウェア

MATLAB: MATLABがインストールされている必要があります。
VSCode: 最新バージョンのVisual Studio Codeをインストールしてください。
MATLAB Extension: VSCodeでMATLABコードをサポートするための拡張機能。
推奨拡張機能:
- Code Runner: ワンクリックでコードを実行。
- Bracket Pair Colorizer: 括弧の対応を色分け。

2. VSCodeのセットアップ

2.1 MATLAB拡張機能のインストール

VSCodeを起動します。
サイドバーの「Extensions(拡張機能)」アイコンをクリックします（ショートカットキー: Ctrl+Shift+X）。
検索バーに「MATLAB」と入力し、MathWorksが提供する公式のMATLAB拡張機能をインストールします。

MATLAB Extension

2.2 MATLAB CLIパスの設定

MATLABコードを実行するには、MATLABのコマンドラインインターフェイス (CLI) へのパスを指定する必要があります。

MATLABがインストールされているディレクトリを確認します:
- Windowsの場合: C:\Program Files\MATLAB\R202Xx\bin\matlab.exe
- macOSの場合: /Applications/MATLAB_R202Xx.app/bin/matlab
- Linuxの場合: /usr/local/MATLAB/R202Xx/bin/matlab
VSCodeの設定を開きます（Ctrl+,）。
検索バーにMATLAB Pathと入力し、MATLAB CLIのフルパスを設定します。

MATLAB Path Setting

3. MATLABコードの編集と実行

3.1 新しいMATLABスクリプトの作成

VSCodeで新しいファイルを作成します（Ctrl+N）。
ファイルをexample.mとして保存します。
MATLABコードを書きます。以下は簡単な例です：

x = 0:0.1:10;
y = sin(x);
plot(x, y);
title('Sine Wave');
xlabel('x');
ylabel('sin(x)');

3.2 MATLABコードの実行

ターミナルを開きます（`Ctrl+``）。
MATLAB CLIを使用してスクリプトを実行します：

matlab -r "example"

MATLABのデスクトップ環境が起動し、プロットが表示されます。

3.3 拡張機能 Code Runner による実行

VSCodeで「Code Runner」拡張機能をインストールします。
- 拡張機能の検索バーに「Code Runner」と入力し、インストールをクリック。
MATLABスクリプトを開き、右上の「Run Code」ボタンをクリックするだけでコードが実行されます。
- または、Ctrl+Alt+Nで実行可能。
実行結果はVSCodeの出力ターミナルに表示されます。

Code Runner

4. 補足情報

MATLAB Live Scriptsの注意点

VSCodeではMATLAB Live Scripts (.mlxファイル) の編集は直接サポートされていません。これらを編集する場合は、MATLAB専用のエディタを使用してください。

まとめ

VSCodeでMATLABを使用することで、柔軟な作業環境を構築できます。特に、他のプログラミング言語と組み合わせたプロジェクトを扱う際には便利です。ぜひ試してみてください！

3.5.1 同定対象システム​

システム特性​

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

数学的定式化​

シミュレーション設定​

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

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

時変システムモデル​

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

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

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

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

MATLABコード（遅延入力）​

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

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

性能比較と考察​

手法比較表​

忘却係数の影響​

実用的考慮事項​

まとめ​

制御対象​

状態空間表現への変換​

PID制御器​

PID制御器伝達関数​

数値積分手法による離散化​

1. 前進差分​

PID制御器のZ変換表現​

P-I-D分解実装（推奨方法）​

状態空間表現（前進差分）​

2. 後退差分​

PID制御器のZ変換表現​

P-I-D分解実装​

状態空間表現（後退差分）​

3. 双一次変換法​

PID制御器のZ変換表現​

P-I-D分解実装​

状態空間表現（双一次変換）​

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

プログラミング上の注意点​

制御ループの実装​

完全なシミュレーションコード​

PID制御器比較​

P-I-D個別効果比較​

結果可視化​

性能比較と考察​

各手法の特性​

PID制御の各要素の役割​

P制御（比例制御）​

I制御（積分制御）​

D制御（微分制御）​

実装上のポイント​

まとめ​

電気回路モデル​

数値積分手法​

1. 前進差分（Forward Difference）​

2. 後退差分（Backward Difference）​

3. 双一次変換（Bilinear Transform）​

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

初期設定​

パラメータ設定​

配列初期化​

各手法による数値計算​

真の解析解​

前進差分​

後退差分​

双一次変換​

結果のプロット​

数値積分手法の特性比較​

1. 前進差分の特性​

2. 後退差分の特性​

3. 双一次変換の特性​

精度と安定性の評価​

サンプリング周期の影響​

適用場面​

前進差分​

後退差分​

双一次変換​

まとめ​

システムモデル​

Z変換による表現​

3.5.1 同定対象システム

システム特性

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

数学的定式化

シミュレーション設定

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

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

時変システムモデル

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

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

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

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

MATLABコード（遅延入力）

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

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

性能比較と考察

手法比較表

忘却係数の影響

実用的考慮事項

まとめ

制御対象

状態空間表現への変換

PID制御器

PID制御器伝達関数

数値積分手法による離散化

1. 前進差分

PID制御器のZ変換表現

P-I-D分解実装（推奨方法）

状態空間表現（前進差分）

2. 後退差分

PID制御器のZ変換表現

P-I-D分解実装

状態空間表現（後退差分）

3. 双一次変換法

PID制御器のZ変換表現

P-I-D分解実装

状態空間表現（双一次変換）

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

プログラミング上の注意点

制御ループの実装

完全なシミュレーションコード

PID制御器比較

P-I-D個別効果比較

結果可視化

性能比較と考察

各手法の特性

PID制御の各要素の役割

P制御（比例制御）

I制御（積分制御）

D制御（微分制御）

実装上のポイント

まとめ

電気回路モデル

数値積分手法

1. 前進差分（Forward Difference）

2. 後退差分（Backward Difference）

3. 双一次変換（Bilinear Transform）

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

初期設定

パラメータ設定

配列初期化

各手法による数値計算

真の解析解

前進差分

後退差分

双一次変換

結果のプロット

数値積分手法の特性比較

1. 前進差分の特性

2. 後退差分の特性

3. 双一次変換の特性

精度と安定性の評価

サンプリング周期の影響

適用場面

前進差分

後退差分

双一次変換

まとめ

システムモデル

Z変換による表現