図12.6〜図12.9のプロット

#計算とグラフプロットに必要なモジュールの読み込み
import numpy as np
from control import matlab
from matplotlib import pyplot as plt
from scipy import arange 
from scipy import signal

#比例要素Kの分子・分母多項式を与える
numk = [10] #比例要素の分子多項式
denk = [1] #比例要素の分母多項式

#G_1(s)の分子・分母多項式を与える
num1 = [0, 1] #G_1(s)の分子多項式
den1 = [10, 1] #G_1(s)の分母多項式

#G_2(s)の分子・分母多項式を与える
num2 = [1, 1] #G_2(s)の分子多項式
den2 = [0, 1] #G_2(s)の分母多項式

#G_3(s)の分子・分母多項式を与える
num3 = [0, 1] #G_3(s)の分子多項式
den3 = [0.1, 1] #G_3(s)の分母多項式

#伝達関数表現を与える
sysK = signal.lti(numk, denk) #比例要素の伝達関数表現（signal.ltiの場合）
sys1 = signal.lti(num1, den1) #G_1(s)の伝達関数表現（signal.ltiの場合）
sys2 = signal.lti(num2, den2) #G_2(s)の伝達関数表現（signal.ltiの場合）
sys3 = signal.lti(num3, den3) #G_3(s)の伝達関数表現（signal.ltiの場合）

/Users/kaz/.pyenv/versions/anaconda3-4.3.1/lib/python3.6/site-packages/scipy/signal/filter_design.py:1549: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless
  "results may be meaningless", BadCoefficients)

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

# ゲインと位相の計算
wK, gainK, phaseK = signal.bode(sysK, w) #比例要素のゲインと位相
w1, gain1, phase1 = signal.bode(sys1, w) #G_1(s)ゲインと位相
w2, gain2, phase2 = signal.bode(sys2, w) #G_2(s)ゲインと位相
w3, gain3, phase3 = signal.bode(sys3, w) #G_3(s)ゲインと位相

# 図12.6のプロット

plt.subplot(2, 1, 1) #複数の図を並べるためのコマンド．2行1列の1行目という意味
plt.semilogx(wK, gainK) #ゲイン線図をプロット
plt.xlim([0.001,1000])#横軸（角周波数）の範囲の指定
plt.ylim([-40,40]) #縦軸の範囲の指定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Gain[dB]") #縦軸のラベル表示
plt.title("Bode diagram of K") #タイトルの表示
plt.show() #グラフの表示

# 位相線図のプロット
plt.subplot(2, 1, 2) #複数の図を縦に並べるためのコマンド．2行1列の2行目という意味
plt.semilogx(wK, phaseK) #位相線図をプロット
plt.xlim([0.001,1000])#横軸（角周波数）の範囲の指定
plt.ylim([-90,0]) #縦軸の範囲の指定
plt.yticks([-90,-45,0]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Phase[deg]") #縦軸のラベル表示
plt.show() #グラフの表示

# 図12.7のプロット

plt.subplot(2, 1, 1) #複数の図を並べるためのコマンド．2行1列の1行目という意味
plt.semilogx(w1, gain1) #ゲイン線図をプロット
plt.xlim([0.001,1000])#横軸（角周波数）の範囲の指定
plt.ylim([-40,40]) #縦軸の範囲の指定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Gain[dB]") #縦軸のラベル表示
plt.title("Bode diagram of G_1(s)") #タイトルの表示
plt.show() #グラフの表示

# 位相線図のプロット
plt.subplot(2, 1, 2) #複数の図を縦に並べるためのコマンド．2行1列の2行目という意味
plt.semilogx(w1, phase1) #位相線図をプロット
plt.xlim([0.001,1000])#横軸（角周波数）の範囲の指定
plt.ylim([-90,0]) #縦軸の範囲の指定
plt.yticks([-90,-45,0]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Phase[deg]") #縦軸のラベル表示
plt.show() #グラフの表示

# 図12.8のプロット

plt.subplot(2, 1, 1) #複数の図を並べるためのコマンド．2行1列の1行目という意味
plt.semilogx(w2, gain2) #ゲイン線図をプロット
plt.xlim([0.001,1000])#横軸（角周波数）の範囲の指定
plt.ylim([-40,40]) #縦軸の範囲の指定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Gain[dB]") #縦軸のラベル表示
plt.title("Bode diagram of G_2(s)") #タイトルの表示
plt.show() #グラフの表示

# 位相線図のプロット
plt.subplot(2, 1, 2) #複数の図を縦に並べるためのコマンド．2行1列の2行目という意味
plt.semilogx(w2, phase2) #位相線図をプロット
plt.xlim([0.001,1000])#横軸（角周波数）の範囲の指定
plt.ylim([-90,0]) #縦軸の範囲の指定
plt.yticks([-90,-45,0]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Phase[deg]") #縦軸のラベル表示
plt.show() #グラフの表示

# 図12.9のプロット

plt.subplot(2, 1, 1) #複数の図を並べるためのコマンド．2行1列の1行目という意味
plt.semilogx(w3, gain3) #ゲイン線図をプロット
plt.xlim([0.001,1000])#横軸（角周波数）の範囲の指定
plt.ylim([-40,40]) #縦軸の範囲の指定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Gain[dB]") #縦軸のラベル表示
plt.title("Bode diagram of G_3(s)") #タイトルの表示
plt.show() #グラフの表示

# 位相線図のプロット
plt.subplot(2, 1, 2) #複数の図を縦に並べるためのコマンド．2行1列の2行目という意味
plt.semilogx(w3, phase3) #位相線図をプロット
plt.xlim([0.001,1000])#横軸（角周波数）の範囲の指定
plt.ylim([-90,0]) #縦軸の範囲の指定
plt.yticks([-90,-45,0]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Phase[deg]") #縦軸のラベル表示
plt.show() #グラフの表示

図12.10のプロット

#図12.10のプロット

# ゲイン線図のプロット
plt.subplot(2, 1, 1) #2つの図を縦に並べるためのコマンド．2行1列の1行目という意味
plt.semilogx(w3, gainK + gain1 + gain2 + gain3) #ゲイン線図をプロット
plt.xlim([0.001,1000])#横軸（角周波数）の範囲の指定
plt.ylim([-40,40]) #縦軸の範囲の指定
plt.grid(color='gray') #罫線を灰色で表示
#plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("gain[dB]") #縦軸のラベル表示
plt.title("Bode diagram of G(s)")

# 位相線図のプロット
plt.subplot(2, 1, 2) #2つの図を縦に並べるためのコマンド．2行1列の2行目という意味
plt.semilogx(w3, phaseK + phase1 + phase2 + phase3) #位相線図をプロット
plt.xlim([0.001,1000])#横軸（角周波数）の範囲の指定
plt.ylim([-90,0]) #縦軸の範囲の指定
plt.yticks([-90,-45,0]) #縦軸の目盛りの値の設定
plt.grid(color='gray') #罫線を灰色で表示
plt.xlabel("w[rad/s]") #横軸のラベル表示
plt.ylabel("Phase[deg]") #縦軸のラベル表示
plt.show() #グラフの表示