図14.8のプロット
In [1]:
#計算とグラフプロットに必要なモジュールの読み込み
import numpy as np
from control import matlab
from matplotlib import pyplot as plt
from scipy import arange
from scipy import signal
In [2]:
#制御対象の分子・分母多項式を与える
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 = matlab.tf( num, den ) #制御対象の伝達関数表現
c0 = matlab.tf( numc0, denc0 ) #C_{0}の伝達関数表現
c1 = matlab.tf( numc1, denc1 ) #C_{1}の伝達関数表現
cPLE = matlab.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)の場合の開ループ伝達関数
In [3]:
#開ループ伝達関数L_0(s)の分子多項式と伝達関数表現を与える
[numL00, denL00] = matlab.tfdata(sysL0) #開ループ伝達関数の分子・分母多項式を抜き出す
numL0 = numL00[0][0] #分子多項式におけるPythonの配列の処理
denL0 = denL00[0][0] #分母多項式におけるPythonの配列の処理
#開ループ伝達関数L_1(s)の分子多項式と伝達関数表現を与える
[numL11, denL11] = matlab.tfdata(sysL1) #開ループ伝達関数の分子・分母多項式を抜き出す
numL1 = numL11[0][0] #分子多項式におけるPythonの配列の処理
denL1 = denL11[0][0] #分母多項式におけるPythonの配列の処理
#開ループ伝達関数L_2(s)の分子多項式と伝達関数表現を与える
[numL22, denL22] = matlab.tfdata(sysL2) #開ループ伝達関数の分子・分母多項式を抜き出す
numL2 = numL22[0][0] #分子多項式におけるPythonの配列の処理
denL2 = denL22[0][0] #分母多項式におけるPythonの配列の処理
#各コントローラに対応した開ループ伝達関数を与える
sysLG0 = signal.lti(numL0, denL0) #L_0の場合の伝達関数表現（signal.ltiの場合）
sysLG1 = signal.lti(numL1, denL1) #L_1の場合の伝達関数表現（signal.ltiの場合）
sysLG2 = signal.lti(numL2, denL2) #L_2の場合の伝達関数表現（signal.ltiの場合）
In [4]:
#角周波数の範囲を指定
w = np.logspace(-1, 2, 1000) #対数的に等間隔なベクトルの生成(10^{-1}から10^{2}で1000点)
#図14.8(a)
#ゲインと位相の計算
wL0, gainL0, phaseL0 = signal.bode(sysLG0, w) #ゲインと位相
#図14.8(a)のプロット
plt.subplot(2, 3, 1) #複数の図を並べるためのコマンド．2行3列の1行1列目という意味
plt.semilogx(wL0, gainL0) #ゲイン線図をプロット
plt.xlim([0.1,100]) #横軸（角周波数）の範囲の指定
plt.ylim([-20,60]) #縦軸の範囲の指定
plt.yticks([-20,0,20,40]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
#plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Gain[dB]") #縦軸のラベル表示
plt.title("L_0(s) = P(s)") #タイトルの表示
plt.subplot(2, 3, 4) #複数の図を並べるためのコマンド．2行3列の2行1列目という意味
plt.semilogx(wL0, phaseL0) #位相線図をプロット
plt.xlim([0.1,100]) #横軸（角周波数）の範囲の指定
plt.ylim([-180,-90]) #縦軸の範囲の指定
plt.yticks([-180,-135,-90]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Phase[deg]") #縦軸のラベル表示
#図14.8(b)
#ゲインと位相の計算
wL1, gainL1, phaseL1 = signal.bode(sysLG1, w) #ゲインと位相
#図14.8(b)のプロット
plt.subplot(2, 3, 2) #複数の図を並べるためのコマンド．2行3列の1行2列目という意味
plt.semilogx(wL1, gainL1) #ゲイン線図をプロット
plt.xlim([0.1,100]) #横軸（角周波数）の範囲の指定
plt.ylim([-20,60]) #縦軸の範囲の指定
plt.yticks([-20,0,20,40]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
#plt.xlabel("w[rad/s]") #横軸のラベル表示
#plt.ylabel("Gain[dB]") #縦軸のラベル表示
plt.title("L_1(s) = P(s)C_1(s)") #タイトルの表示
plt.subplot(2, 3, 5) #複数の図を並べるためのコマンド．2行3列の2行2列目という意味
plt.semilogx(wL1, phaseL1) #位相線図をプロット
plt.xlim([0.1,100]) #横軸（角周波数）の範囲の指定
plt.ylim([-180,-90]) #縦軸の範囲の指定
plt.yticks([-180,-135,-90]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Phase[deg]") #縦軸のラベル表示
#図14.8(c)
#ゲインと位相の計算
wL2, gainL2, phaseL2 = signal.bode(sysLG2, w) #ゲインと位相
#図14.8(c)のプロット
plt.subplot(2, 3, 3) #複数の図を並べるためのコマンド．2行3列の1行3列目という意味
plt.semilogx(wL2, gainL2) #ゲイン線図をプロット
plt.xlim([0.1,100]) #横軸（角周波数）の範囲の指定
plt.ylim([-20,60]) #縦軸の範囲の指定
plt.yticks([-20,0,20,40]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
#plt.xlabel("w[rad/s]") #横軸のラベル表示
#plt.ylabel("Gain[dB]") #縦軸のラベル表示
plt.title("L_2(s) = P(s)C_2(s)") #タイトルの表示
plt.subplot(2, 3, 6) #複数の図を並べるためのコマンド．2行3列の2行3列目という意味
plt.semilogx(wL2, phaseL2) #位相線図をプロット
plt.xlim([0.1,100]) #横軸（角周波数）の範囲の指定
plt.ylim([-180,-90]) #縦軸の範囲の指定
plt.yticks([-180,-135,-90]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("w[rad/s]") #横軸のラベル表示
#plt.ylabel("Phase[deg]") #縦軸のラベル表示
plt.show() #グラフの表示
In [5]:
#図14.8のプロット
# ゲイン線図のプロット
plt.subplot(2, 1, 1) #2つの図を縦に並べるためのコマンド．2行1列の1行目という意味
plt.semilogx(wL0, gainL0, label="L_0(s) = P(s)") #ゲイン線図をプロット
plt.semilogx(wL1, gainL1, label="L_1(s) = P(s)C_1(s)")#ゲイン線図をプロット
plt.semilogx(wL2, gainL2, label="L_2(s) = P(s)C_2(s)")#ゲイン線図をプロット
plt.xlim([0.1,100]) #横軸（角周波数）の範囲の指定
plt.ylim([-20,60]) #縦軸の範囲の指定
plt.grid(color='gray') #罫線を灰色で表示
#plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Gain[dB]") #縦軸のラベル表示
plt.legend() #凡例の表示
# 位相線図のプロット
plt.subplot(2, 1, 2) #2つの図を縦に並べるためのコマンド．2行1列の2行目という意味
plt.semilogx(wL0, phaseL0, label="L_0(s) = P(s)") #位相線図をプロット
plt.semilogx(wL1, phaseL1, label="L_1(s) = P(s)C_1(s)") #位相線図をプロット
plt.semilogx(wL2, phaseL2, label="L_2(s) = P(s)C_2(s)") #位相線図をプロット
plt.xlim([0.1,100]) #横軸（角周波数）の範囲の指定
plt.ylim([-180,-90]) #縦軸の範囲の指定
plt.yticks([-180,-135,-90]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Phase[deg]") #縦軸のラベル表示
plt.legend() #凡例の表示
plt.show() #グラフの表示
図14.9のプロット
In [6]:
#制御対象とコントローラのフィードバック結合を求める
sysc0 = matlab.feedback(sysL0, 1, -1) #L_0(s) = P(s)の場合のフィードバック結合
sysc1 = matlab.feedback(sysL1, 1, -1) #L_1(s) = P(s)C_1(s)の場合のフィードバック結合
sysc2 = matlab.feedback(sysL2, 1, -1) #L_2(s) = P(s)C_2(s)の場合のフィードバック結合
#時間変数の定義
t = arange(0, 6, 0.01) #0から6まで0.01刻み
#ステップ応答の計算
y0, t0 = matlab.step(sysc0, t) #L_0(s) = P(s)の場合のステップ応答
y1, t1 = matlab.step(sysc1, t) #L_1(s) = P(s)C_1(s)場合のステップ応答
y2, t2 = matlab.step(sysc2, t) #L_2(s) = P(s)C_2(s)の場合のステップ応答
#図14.9(a)のプロット
plt.subplot(1,3,1) #複数の図を並べるためのコマンド．1行3列の1列目という意味
plt.plot(t0, y0) #ステップ応答をプロット
plt.xlim([0,6]) #横軸(時間軸）の範囲の指定
plt.yticks([0,1.0,2.0]) #縦軸の範囲の指定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("time t[s]") #横軸のラベル表示
plt.ylabel("y(t)") #縦軸のラベル表示
plt.title("C_0(s) = 1") #タイトルの表示
#図14.9(b)のプロット
plt.subplot(1,3,2) #複数の図を並べるためのコマンド．1行3列の2列目という意味
plt.plot(t1, y1) #ステップ応答をプロット
plt.xlim([0,6]) #横軸(時間軸）の範囲の指定
plt.yticks([0,1.0,2.0]) #縦軸の範囲の指定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("time t[s]") #横軸のラベル表示
#plt.ylabel("y(t)") #縦軸のラベル表示
plt.title("C_1(s) = 12") #タイトルの表示
#図14.9(c)のプロット
plt.subplot(1,3,3) #複数の図を並べるためのコマンド．1行3列の3列目という意味
plt.plot(t2, y2) #ステップ応答をプロット
plt.xlim([0,6]) #横軸(時間軸）の範囲の指定
plt.yticks([0,1.0,2.0]) #縦軸の範囲の指定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("time t[s]") #横軸のラベル表示
#plt.ylabel("y(t)") #縦軸のラベル表示
plt.title("C_2(s) = 12*(10/1)*(s+1)/(s+10)") #タイトルの表示
plt.show() #グラフの表示
In [7]:
plt.plot(t0, y0, label = "C_0(s) = 1") #ステップ応答をプロット
plt.plot(t1, y1, label = "C_1(s) = 12") #ステップ応答をプロット
plt.plot(t2, y2, label = "C_2(s) = 12*(10/1)*(s+1)/(s+10)") #ステップ応答をプロット
plt.xlim([0,6]) #横軸(時間軸）の範囲の指定
plt.yticks([0,1.0,2.0]) #縦軸の範囲の指定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("time t[s]") #横軸のラベル表示
plt.ylabel("y(t)") #縦軸のラベル表示
plt.legend() #凡例の表示
plt.show() #グラフの表示
図14.10のプロット
In [8]:
#位相進みコントローラC_2(s)を与える
C2 = c1*cPLE #C_2(s)
#C_2(s)の分子多項式と伝達関数表現を与える
[numC22, denC22] = matlab.tfdata(C2) #C_2(s)の分子・分母多項式を抜き出す
numC2= numC22[0][0] #分子多項式におけるPythonの配列の処理
denC2 = denC22[0][0] #分母多項式におけるPythonの配列の処理
#C_2(s)の伝達関数を与える
sysC2 = signal.lti(numC2, denC2) #C_2の場合の伝達関数（signal.ltiの場合）
#ゲインと位相の計算
wC2, magC2, phaseC2 = signal.bode(sysC2, w) #ゲインと位相
#図14.10のプロット
plt.subplot(2, 1, 1) #複数の図を並べるためのコマンド．2行1列の1行目という意味
plt.semilogx(wC2, magC2) #ゲイン線図をプロット
plt.xlim([0.1,100]) #横軸（角周波数）の範囲の指定
plt.ylim([-10,50]) #縦軸の範囲の指定
plt.yticks([0,20,40]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
#plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Gain[dB]") #縦軸のラベル表示
plt.title("C_2(s)") #タイトルの表示
plt.subplot(2, 1, 2) #複数の図を並べるためのコマンド．2行1列の2行目という意味
plt.semilogx(wC2, phaseC2) #位相線図をプロット
plt.xlim([0.1,100]) #横軸（角周波数）の範囲の指定
plt.ylim([0,90]) #縦軸の範囲の指定
plt.yticks([0, 45, 90]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Phase[deg]") #縦軸のラベル表示
plt.show()
