2 posts tagged with "後退差分"

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

May 26, 2024 · 9 min read

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. 双一次変換の特性

特徴：

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

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

周波数特性の保持

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

精度と安定性の評価

サンプリング周期の影響

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

適用場面

前進差分

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

後退差分

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

双一次変換

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

まとめ

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

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

手法選択の指針

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

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

制御対象​

状態空間表現への変換​

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. 双一次変換の特性​

精度と安定性の評価​

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

適用場面​

前進差分​

後退差分​

双一次変換​

まとめ​

制御対象

状態空間表現への変換

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. 双一次変換の特性

精度と安定性の評価

サンプリング周期の影響

適用場面

前進差分

後退差分

双一次変換

まとめ