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 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による性能評価

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

制御対象​

状態空間表現への変換​

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制御（微分制御）​

実装上のポイント​

まとめ​

システムモデル​

Z変換による表現​

伝達関数​

部分分数展開​

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

初期設定​

パラメータ設定​

配列初期化​

システム応答計算​

観測雑音の追加​

結果のプロット​

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

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

システムの安定性​

応用例​

実装のポイント​

初期値の設定​

雑音生成​

プロット設定​

まとめ​

制御対象

状態空間表現への変換

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制御（微分制御）

実装上のポイント

まとめ

システムモデル

Z変換による表現

伝達関数

部分分数展開

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

初期設定

パラメータ設定

配列初期化

システム応答計算

観測雑音の追加

結果のプロット

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

雑音信号比（NSR）の影響

システムの安定性

応用例

実装のポイント

初期値の設定

雑音生成

プロット設定

まとめ