次のような微分方程式
で表される系の伝達関数を考えます。微分方程式をラプラス変換すると伝達関数は
のようになります。
##ボード線図 制御工学にはシステムを解析するためのボード線図があり、ゲイン線図と位相線図の組み合わせで使われることが多いです。 ゲイン線図を書くにはまず$s=j\omega$を代入すればよく
の後ゲインを
によって求めます。また、位相線図は
によって求めます。
今回は一般で解くわけにも行きませんので、バネマスダンパ系
を例に使用します。バネマスダンパ系については過去記事[1]を参照してください。
コマンドを叩くだけ。 bodeが関数、sysには伝達関数を引数で与えます。
bode(sys)
grid on
伝達関数は
num=[1];
den=[4 2 5];
sys=tf(num,den);
のように定義します。
結果は
%% 関数を使わずに作成
% 対数的に等間隔な行列を用意
omega = logspace(-1,2,1000);
figure;
plot(omega)
grid on
logspace(A,B,C)で対数的に等間隔な行列を用意します。 この例では-1から2までの間を1000個対数的に等間隔に分割しています。
Gjw = polyval(num,sqrt(-1)*omega)./polyval(den,sqrt(-1)*omega);
Mag = 20*log10(abs(Gjw));
Phase = rad2deg(unwrap(angle(Gjw)));
polyval(A,B)は多項式AにBを代入します。 あとは定義通りです。 結果は
clear
close all
%% 伝達関数を定義
num=[1];
den=[4 2 5];
sys=tf(num,den);
%% コマンドで描画
bode(sys)
grid on
%% 関数を使わずに作成
% 対数的に等間隔な行列を用意
omega = logspace(-1,2,1000);
figure;
plot(omega)
grid on
Gjw = polyval(num,sqrt(-1)*omega)./polyval(den,sqrt(-1)*omega);
Mag = 20*log10(abs(Gjw));
Phase = rad2deg(unwrap(angle(Gjw)));
figure;
subplot(211)
semilogx(omega,Mag)
ylabel('Magnitude (dB)');
xlim([-10^-1 10^1]);
ylim([-60 0]);
yticks([-60 -50 -40 -30 -20 -10 0]);
grid on
subplot(212)
semilogx(omega,Phase)
xlabel('Frequency (rad/s)');
ylabel('Phase (deg)');
xlim([-10^-1 10^1])
ylim([-182 0]);
yticks([-180 -135 -90 -45 0]);
grid on