l型のナイキスト軌跡をscipyで書いてみた

scipyによるナイキスト軌跡の書き方

パッケージのインポート

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

開ループ関数
$L(s) = \frac{30}{(s + 1)(s + 2)(s + 3)}$
のナイキスト軌跡

s = signal.ZerosPolesGain([], [-1, -2, -3], 30)
w, H = signal.freqresp(s)
print(type(H))
plt.plot(H.real, H.imag, "b")
plt.plot(H.real, -H.imag, "r")
plt.grid()

ナイキスト軌跡から閉ループ系は安定であることがわかる．

1型の開ループ関数

$L(s) = \frac{3}{s(s + 1)(s + 2)}$
上と同様の方法で1型の開ループ関数を描いてみる．

s = signal.ZerosPolesGain([], [0, -1, -2], 3)
w = 10**np.linspace(-1, 4, 1024)
_, H = signal.freqresp(s, w)
plt.plot(H.real, H.imag)
plt.plot(H.real, -H.imag, "--")
plt.grid()

開ループ伝達関数に不安定極がないので，この場合閉ループ系は安定だと考えられるが，
l型の開ループ系でもナイキスト軌跡を閉じた軌跡にする場合，下図のような経路に沿って積分を行う．

この時，$s = 0$の極は積分経路の外なので，開ループ伝達関数の不安定極とは数えない．同様の伝達関数
$L(s) = \frac{3}{s(s + 1)(s + 2)}$
に対してナイキスト軌跡を描くと

s = signal.ZerosPolesGain([], [0, -1, -2], 3)
w = 10**np.linspace(-0.5, 4, 1024)
_, H = signal.freqresp(s, w)
plt.plot(H.real, H.imag, "b")
plt.plot(H.real, -H.imag, "b--")
z = 10**(-0.5)*np.exp(1j*np.linspace(0, np.pi/2, 1024))
Hn = K/((z*(z + 1)*(z + 2)))
plt.plot(Hn.real, Hn.imag, "r")
z = 10**(-0.5)*np.exp(1j*np.linspace(-np.pi/2, 0, 1024))
Hn = K/((z*(z + 1)*(z + 2)))
plt.plot(Hn.real, Hn.imag, "r--")
plt.plot(-1, 0, "kx")
plt.grid()

出力された図は閉ループ系のベクトル軌跡と(-1, 0)の位置関係がわかりやすくなっている．しかし，実際には積分経路に閉ループ系の極が全て入っている必要があることから，大きい半円の半径を十分大きく，小さい半円の半径を十分小さくしなければならない．出力された図に関しては，見やすいが，その意味では不正確なものである．

小さい円の半径を小さくすると青の軌跡は無限遠に飛んでいくので，図が見づらくなる．
matlabだと(-1, 0)付近を拡大できるので，正確だと思う．

2型の場合

$L(s) = \frac{1}{s^2(s + 1)}$
積分経路は1型の時と同様で，経路はうまく取れているとする．

s = signal.ZerosPolesGain([], [0, 0, -1], 1)
w = 10**np.linspace(-0.5, 4, 1024)
_, H = signal.freqresp(s, w)
plt.plot(H.real, H.imag, "b")
plt.plot(H.real, -H.imag, "b--")
z = 10**(-0.5)*np.exp(1j*np.linspace(0, np.pi/2, 1024))
Hn = 1/(z**2*(z + 1))
plt.plot(Hn.real, Hn.imag, "r")
z = 10**-0.5*np.exp(1j*np.linspace(-np.pi/2, 0, 1024))
Hn = 1/(z**2*(z + 1))
plt.plot(Hn.real, Hn.imag, "r--")
plt.plot(-1, 0, "kx")
plt.grid()

閉ループ経路は(-1, 0)の周りを時計回りに2回りしているので，不安定である．

3型の場合

$L(s) = \frac{K(s + 10)^2}{s^3}$
ゲインKはK > 5の時安定であることが，ラウス=フルビッツの方法などを使えばわかる．図はK = 10の場合である．

K = 10
m = 0.9
w = 10**np.linspace(m, 2, 1024)
_, H = signal.freqresp(([-10, -10], [0, 0, 0], K), w)
print(w)
plt.plot(H.real, H.imag, "b")
plt.plot(H.real, -H.imag, "b--")
plt.plot(-1, 0, "kx")
plt.grid()

z = 10**m*np.exp(1j*np.linspace(0, np.pi/2, 1024))
H = K*(z + 10)**2/z**3
plt.plot(H.real, H.imag, "r")
z = 10**m*np.exp(1j*np.linspace(-np.pi/2, 0, 1024))
H = K*(z + 10)**2/z**3
plt.plot(H.real, H.imag, "r--")

見やすくするために，小さい方の円を大きくとってある．
この図を見ると，(-1, 0)は閉ループの経路の進む向きに対して外側にあって，
その周りを回っていないことがわかる．

おまけ

積分経路の図の書き方

z = np.exp(1j*np.linspace(-np.pi/2, np.pi/2, 1024))
plt.plot(z.real, z.imag, "b", 5*z.real, 5*z.imag, "b")
y = np.linspace(1, 5, 1024)
x = np.zeros(len(y))
plt.plot(x, y, "b", x, -y, "b")
plt.grid()
plt.axis("equal")
plt.quiver(1, -0.3, 0, 1, scale=25, color="b")
plt.quiver(5, 0.2, 0, -1, scale=25, color="b")
plt.quiver(0, 3, 0, 1, scale=25, color="b")
plt.quiver(0, -3, 0, 1, scale=25, color="b")
plt.xticks([0])
plt.yticks([0])

参考文献

杉江俊治・藤田政之.フィードバック制御入門.コロナ社,1999.

感謝

読んでくれてありがとうございます．