%図14.8, 図14.9，図14.10
%
clear all; close all;

%図14.8のプロット

%伝達関数の分子・分母多項式を与える
num = [ 0 0 1 ]; %分子多項式
den = [ 1 1 0 ]; %分母多項式

%コントローラC_{0}の分子・分母多項式を与える．
numc0 = [ 1 ]; %C_{0}の分子多項式
denc0 = [ 1 ]; %C_{0}の分母多項式

%コントローラC_{1}=K_pの分子・分母多項式を与える．
numc1 = [ 12 ]; %C_{1}の分子多項式
denc1 = [ 1 ]; %C_{1}の分母多項式

%位相進みコントローラのパラメータと分子・分母多項式を与える．
omega3 = 10; %位相進みコントローラのパラメータ
omega4 = 1; %位相進みコントローラのパラメータ
numPLE = [ omega3 omega3 * omega4 ]; %位相進みコントローラの分子多項式
denPLE = [ omega4 omega4 * omega3 ]; %位相進みコントローラの分母多項式

%制御対象とコントローラの伝達関数表現を与える
sys = tf( num, den ); %制御対象の伝達関数表現
c0 = tf( numc0, denc0 ); %C_{0}の伝達関数表現
c1 = tf( numc1, denc1 ); %C_{1}の伝達関数表現
cPLE = tf( numPLE, denPLE ); %位相進みコントローラの伝達関数表現

%開ループ伝達関数を求める
sysL0 = c0*sys; %L_0(s) = P(s)の場合の開ループ伝達関数
sysL1 = c1*sys; %L_1(s) = P(s)C_1(s)の場合の開ループ伝達関数
sysL2 = c1*cPLE*sys; %L_2(s) = P(s)C_2(s)の場合の開ループ伝達関数

%角周波数の範囲を指定
w = logspace(-1, 2, 1000); %対数的に等間隔なベクトルの生成(10^{-1}から10^{2}で1000点)

%図14.8のプロット

figure(1) %図のウィンドウを開く

%図14.8(a)のプロット
[ gainL0 phaseL0 wL0 ] = bode( sysL0, w ); %ゲインと位相
gainL0 = gainL0(:); %L_0(s) = P(s)の場合のゲインの配列を1次元ベクトルに変換
gainL1dB = 20*log10(gainL0); %L_0(s) = P(s)の場合のゲインをデシベル値に変換
phaseL0 = phaseL0(:); %L_0(s) = P(s)の場合の位相の配列を1次元ベクトルに変換

subplot(2,3,1) %複数の図を並べるためのコマンド．2行3列の1行1列目という意味
semilogx(wL0, gainL1dB); %ゲイン線図をプロット
xlim([10^(-1),10^2]) %横軸（角周波数）の範囲の指定
ylim([-30 45]) %縦軸の範囲の指定
xticks([10^(-1) 10^(0) 10^(1) 10^(2)]) %横軸の目盛りの値の設定
yticks([-20 0 20 40]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('Frequency (rad/s)'); %横軸のラベル表示
ylabel('Gain [dB]'); %縦軸のラベル表示
title('L_0(s) = P(s)') %タイトルの表示

subplot(2,3,4) %複数の図を並べるためのコマンド．2行3列の2行1列目という意味
semilogx(wL0, phaseL0); %位相線図をプロット
xlim([10^(-1),10^2]) %横軸（角周波数）の範囲の指定
ylim([-180 -90]) %縦軸の範囲の指定
xticks([10^(-1) 10^(0) 10^(1) 10^(2)]) %横軸の目盛りの値の設定
yticks([-180 -135 -90]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('Frequency (rad/s)'); %横軸のラベル表示
ylabel('Phase [deg]'); %縦軸のラベル表示

%図14.8(b)のプロット
[ gainL1 phaseL1 wL1 ] = bode( sysL1, w ); %ゲインと位相
gainL1 = gainL1(:); %L_1(s) = P(s)C_1(s)の場合のゲインの配列を1次元ベクトルに変換
gainL1dB = 20*log10(gainL1); %L_1(s) = P(s)C_1(s)の場合のゲインをデシベル値に変換
phaseL1 = phaseL1(:); %L_1(s) = P(s)C_1(s)の場合の位相の配列を1次元ベクトルに変換

subplot(2,3,2) %複数の図を並べるためのコマンド．2行3列の1行2列目という意味
semilogx(wL1, gainL1dB); %ゲイン線図をプロット
xlim([10^(-1),10^2]) %横軸（時間軸の範囲の指定）
ylim([-30 45]) %縦軸（縦軸の範囲の指定）
xticks([10^(-1) 10^(0) 10^(1) 10^(2)]) %横軸の目盛りの値の設定
yticks([-20 0 20 40]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('Frequency (rad/s)'); %横軸（時間軸）のラベル表示
%ylabel('Gain [dB]'); %縦軸のラベル表示
title('L_1(s) = P(s)C_1(s)')

subplot(2,3,5) %複数の図を並べるためのコマンド．2行3列の2行2列目という意味
semilogx(wL1, phaseL1); %位相線図をプロット
xlim([10^(-1),10^2]) %横軸（時間軸の範囲の指定）
ylim([-180 -90]) %縦軸（縦軸の範囲の指定）
xticks([10^(-1) 10^(0) 10^(1) 10^(2)]) %横軸の目盛りの値の設定
yticks([-180 -135 -90]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('Frequency (rad/s)'); %横軸（時間軸）のラベル表示
%ylabel('Phase [deg]'); %縦軸のラベル表示

%図14.8(c)のプロット
[ gainL2 phaseL2 wL2 ] = bode( sysL2, w ); %ゲインと位相
gainL2 = gainL2(:); %L_2(s) = P(s)C_2(s)の場合のゲインの配列を1次元ベクトルに変換
gainL2dB = 20*log10(gainL2); %L_2(s) = P(s)C_2(s)の場合のゲインをデシベル値に変換
phaseL2 = phaseL2(:); %L_2(s) = P(s)C_2(s)の場合の位相の配列を1次元ベクトルに変換

subplot(2,3,3) %複数の図を並べるためのコマンド．2行3列の1行3列目という意味
semilogx(wL2, gainL2dB); %ゲイン線図をプロット
xlim([10^(-1),10^2]) %横軸（角周波数）の範囲の指定
ylim([-30 45]) %縦軸の範囲の指定
xticks([10^(-1) 10^(0) 10^(1) 10^(2)]) %横軸の目盛りの値の設定
yticks([-20 0 20 40]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('Frequency (rad/s)'); %横軸のラベル表示
%ylabel('Gain [dB]'); %縦軸のラベル表示
title('L_2(s) = P(s)C_2(s)') %タイトルの表示

subplot(2,3,6) %複数の図を並べるためのコマンド．2行3列の2行3列目という意味
semilogx(wL2, phaseL2); %位相線図をプロット
xlim([10^(-1),10^2]) %横軸（角周波数）の範囲の指定
ylim([-180 -90]) %縦軸の範囲の指定
xticks([10^(-1) 10^(0) 10^(1) 10^(2)]) %横軸の目盛りの値の設定
yticks([-180 -135 -90]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('Frequency (rad/s)'); %横軸のラベル表示
%ylabel('Phase [deg]'); %縦軸のラベル表示

%図14.9のプロット

figure(2) %図のウィンドウを開く

%制御対象とコントローラのフィードバック結合を求める
sysc0 = feedback( sysL0, 1 ); %L_0(s) = P(s)の場合のフィードバック結合
sysc1 = feedback( sysL1, 1 ); %L_1(s) = P(s)C_1(s)の場合のフィードバック結合
sysc2 = feedback( sysL2, 1 ); %L_2(s) = P(s)C_2(s)の場合のフィードバック結合

%時間変数の定義
t = 0:0.01:6; %0から6まで0.01刻み

%ステップ応答の計算
[ y0, t0 ] = step( sysc0, t ); %L_0(s) = P(s)の場合のステップ応答
[ y1, t1 ] = step( sysc1, t ); %L_1(s) = P(s)C_1(s)場合のステップ応答
[ y2, t2 ] = step( sysc2, t ); %L_2(s) = P(s)C_2(s)の場合のステップ応答

%図14.9(a)のプロット
subplot(1,3,1) %複数の図を並べるためのコマンド．1行3列の1列目という意味
plot(t0,y0); %ステップ応答をプロット
xlim([0,6]) %横軸（時間軸）の範囲の指定
ylim([0.0 2.0]) %縦軸の範囲の指定
xticks([0 2.0 4.0]) %横軸の目盛りの値の設定
yticks([0.0 1.0 2.0]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('time  t [s]'); %横軸のラベル表示
ylabel('y(t)') %縦軸のラベル表示
title('C_0(s) = 1') %タイトルの表示

%図14.9(b)のプロット
subplot(1,3,2) %複数の図を並べるためのコマンド．1行3列の2列目という意味
plot(t1,y1); %ステップ応答をプロット
xlim([0,6]) %横軸（時間軸）の範囲の指定
ylim([0.0 2.0]) %縦軸の範囲の指定
xticks([0 2.0 4.0]) %横軸の目盛りの値の設定
yticks([0.0 1.0 2.0]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('time  t [s]'); %横軸のラベル表示
%ylabel('y(t)') %縦軸のラベル表示
title('C_1(s) = 12') %タイトルの表示

%図14.9(c)のプロット
subplot(1,3,3) %複数の図を並べるためのコマンド．1行3列の3列目という意味
plot(t2,y2); %ステップ応答をプロット
xlim([0,6]) %横軸（時間軸）の範囲の指定
ylim([0.0 2.0]) %縦軸の範囲の指定
xticks([0 2.0 4.0]) %横軸の目盛りの値の設定
yticks([0.0 1.0 2.0]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('time  t [s]'); %横軸のラベル表示
%ylabel('y(t)') %縦軸のラベル表示
title('C_2(s) = 12*(10/1)*(s+1)/(s+10)') %タイトルの表示

%図14.10のプロット
figure(3) %図のウィンドウを開く

%位相進みコントローラC_2(s)を与える
C2 = c1*cPLE; %C_2(s)

[ gainc2 phasec2 wc2 ] = bode( C2, w ); %ゲインと位相
gainc2 = gainc2(:); %C_2(s)のゲインの配列を1次元ベクトルに変換
gainc2dB = 20*log10(gainc2); %C_2(s)のゲインをデシベル値に変換
phasec2 = phasec2(:); %C_2(s)の位相の配列を1次元ベクトルに変換

subplot(2,1,1) %複数の図を並べるためのコマンド．2行1列の1行1列目という意味
semilogx(wc2, gainc2dB); %ゲイン線図をプロット
xlim([10^(-1),10^2]) %横軸（角周波数）の範囲の指定
ylim([-10 50]) %縦軸の範囲の指定
xticks([10^(-1) 10^(0) 10^(1) 10^(2)]) %横軸の目盛りの値の設定
yticks([0 20 40]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('Frequency (rad/s)'); %横軸のラベル表示
ylabel('Gain [dB]'); %縦軸のラベル表示

subplot(2,1,2) %複数の図を並べるためのコマンド．2行1列の2行1列目という意味
semilogx(wc2, phasec2); %位相線図をプロット
xlim([10^(-1),10^2]) %横軸（角周波数）の範囲の指定
ylim([-10 90]) %縦軸の範囲の指定
xticks([10^(-1) 10^(0) 10^(1) 10^(2)]) %横軸の目盛りの値の設定
yticks([0 45 90]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('Frequency (rad/s)'); %横軸のラベル表示
ylabel('Phase [deg]'); %縦軸のラベル表示

Published with MATLAB® R2016b