非線形システムの平衡点を求めて線形化する方法と、非線形システムに単一周波数の正弦波を入力する方法を紹介します。
1. 平衡点を算出して、平衡点を使用して線形化
例えば、単振子の運動方程式は、下記の式で表されます。
..
mlθ+mgsin(θ)=T .....(1)
この運動方程式を表した Simulink モデルが、pend_sys.mdl です。このモデルは、下記の関連ドキュメントからダウンロードできます。
上記の運動方程式は、θ=0 近傍で
sin(θ)=θ
と近似できることから、運動方程式は下記のように書き換えられます。
..
mlθ+mgθ=T
ラプラス変換して伝達関数を求めると、下式になります。
T 1
G(s)= --- = ------------ .....(2)
θ mls^2+mg
次に、運動方程式 (1) について、モデル pend_sys.mdl を用いて線形化を行います。
まず、TRIM 関数で、システムの入力と状態の操作点を求めます。
[x,u]=trim('pend_sys')
x =
0
0
u =
0
得られた操作点を、線形化を行う LINMOD 関数の入力引数に指定します。
[a,b,c,d]=linmod('pend_sys',x,u);
確認のため、伝達関数形式に変換すると（Control System Toolboxが必要）、下記結果が得られます。
tf(ss(a,b,c,d))
伝達関数:
0.1
------------- .....(3)
s^2 + 9.8
モデルでは、m=10, l=1, g=9.8 と置いているので、(2)と(3)の伝達関数が同じであることが分かります。
2. 非線形システムに単一周波数の正弦波を入力
非線形システムにある周波数のサイン波を入力し、出力の応答から振幅・位相を求め、その周波数の振幅・位相のデータとします。これを、複数の周波数のサイン波を入力することにより、ボード線図を描画します。サンプルモデル (pend_sys_sin.mdl) とサンプルプログラム(sin_bode.m) は、下記の関連ドキュメントからダウンロードできます。
モデルでは、Sine Wave ブロックの周波数を変数 w と設定しており、sin_bode.m で変数 w の値を変更しながら、入出力データを取得しています。得られた入出力データを用いて、FFT 関数で複素周波数応答 Pxy を算出し、入力となる正弦波の周波数 w に関する Pxy の値を取得して、最後に振幅・位相を求めてます。下記に、sin_bode.m のコードを記述します。
% 物理パラメータ
m=10;l=1;g=9.8;
% ボード線図を求める周波数ベクトル[rad/s]
wv = [logspace(-1,0,10),logspace(0,1,30),logspace(1,2,10)];
% シミュレーション・FFT用のパラメータ
num_data=32768; % FFTの点数
ts=0.001; % データのサンプル時間
t_end=(num_data-1)*ts; % シミュレーションの最終時刻
% 各周波数を入力とした時のシミュレーションと複素周波数応答を算出
for k = 1:length(wv)
w = wv(k);
sim('pend_sys_sin') % モデルの実行
Pxx = fft(indata,num_data); % 入力データをFFT
Pyy = fft(outdata,num_data); % 出力データをFFT
Pxy = (Pyy.*conj(Pxx))./(Pxx.*conj(Pxx)); % 複素周波数応答
Pxy = Pxy(1:num_data/2);
F = (0:num_data/2-1)*(1/ts)/num_data; % 周波数ベクトル
Pxy_w(k) =interp1(2*pi*F,Pxy,w); % 入力周波数に対する複素周波数応答を計算
end
% 振幅・位相を計算
Mag_dB = 20*log10(abs(Pxy_w));
Ph = unwrap(angle(Pxy_w))/(pi/180);
上記結果を、線形化したシステムの結果と比較するには、続けて下記コマンドを実行します。
% 線形化したシステムのボード線図を計算
sys = tf(1,[m*l 1 m*g]);
wb = logspace(-1,3,100);
[magb,phaseb]=bode(sys,wb);
magb = squeeze(magb); magb_db = 20*log10(magb);
phaseb = squeeze(phaseb);
% 結果の表示
figure(1)
subplot(211)
semilogx(wv,Mag_dB, wb, magb_db)
xlabel('Frequency (rad/s)'); ylabel('Magnitude (dB)');
xlim([0.1 100]); ylim([-100 20]);
subplot(212)
semilogx(wv,Ph, wb, phaseb)
xlabel('Frequency (rad/s)'); ylabel('Phase (deg)');
xlim([0.1 100]); ylim([-180 0]);
legend('estimation','linearized system')
ただし、各周波数においてシミュレーションされた出力応答が、周期性を持つ信号になっている場合に限られますので、出力応答を確認しつつボード線図を求める必要があります。
