5章 安定性¶
In [1]:
from control.matlab import *
import matplotlib.pyplot as plt
import numpy as np
#plt.rcParams['font.family'] ='sans-serif' #使用するフォント
plt.rcParams['font.family'] = 'Times New Roman' # font familyの設定
plt.rcParams['mathtext.fontset'] = 'cm' # math fontの設定
plt.rcParams['xtick.direction'] = 'in' #x軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams['ytick.direction'] = 'in' #y軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams['xtick.major.width'] = 1.0 #x軸主目盛り線の線幅
plt.rcParams['ytick.major.width'] = 1.0 #y軸主目盛り線の線幅
plt.rcParams['font.size'] = 11 #フォントの大きさ
plt.rcParams['axes.linewidth'] = 0.5 # 軸の線幅edge linewidth。囲みの太さ
plt.rcParams['mathtext.default'] = 'it'#'regular'
plt.rcParams['axes.xmargin'] = '0'
plt.rcParams['axes.ymargin'] = '0.05'
plt.rcParams['savefig.facecolor'] = 'None'
plt.rcParams['savefig.edgecolor'] = 'None'
plt.rcParams["legend.fancybox"] = True # 丸角
# plt.rcParams["legend.framealpha"] = 1 # 透明度の指定、0で塗りつぶしなし
# plt.rcParams["legend.edgecolor"] = 'gray' # edgeの色を変更
plt.rcParams["legend.handlelength"] = 1.8 # 凡例の線の長さを調節
plt.rcParams["legend.labelspacing"] = 0.4 # 垂直方向（縦）の距離の各凡例の距離
plt.rcParams["legend.handletextpad"] = 0.7 # 凡例の線と文字の距離の長さ
plt.rcParams["legend.markerscale"] = 1.0 # 点がある場合のmarker scale
In [2]:
def linestyle_generator():
linestyle = ['-', '--', '-.', ':']
lineID = 0
while True:
yield linestyle[lineID]
lineID = (lineID + 1) % len(linestyle)
In [3]:
def plot_set(fig_ax, *args):
fig_ax.set_xlabel(args[0])
fig_ax.set_ylabel(args[1])
fig_ax.grid(ls=':', lw=0.5)
if len(args)==3:
fig_ax.legend(loc=args[2])
In [4]:
def bodeplot_set(fig_ax, *args):
fig_ax[0].grid(which="both", ls=':', lw=0.5)
fig_ax[0].set_ylabel('Gain [dB]')
fig_ax[1].grid(which="both", ls=':', lw=0.5)
fig_ax[1].set_xlabel('$\omega$ [rad/s]')
fig_ax[1].set_ylabel('Phase [deg]')
if len(args) > 0:
fig_ax[1].legend(loc=args[0])
if len(args) > 1:
fig_ax[0].legend(loc=args[1])
In [5]:
# 図を保存するかどうか
is_savefig = False
# 図の保存パス
figpath="./notebook_output/"
In [6]:
# 数式処理のためにsympyをインポート
import sympy as sp
from sympy.matrices import *
例5.1: 伝達関数と状態方程式の安定性の関係¶
状態方程式を定義
In [7]:
A = np.matrix([[4,3],[-6,-5]])
B = np.matrix([[0],[1]])
C1 = np.matrix([1,0])
C2 = np.matrix([1,1])
D = np.matrix([0])
ss1 = ss(A, B, C1, D)
ss2 = ss(A, B, C2, D)
In [8]:
ss1
Out[8]:
\[
\left(
\begin{array}{rllrll|rll}
4\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&3\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\
-6\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&-5\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\
\hline
1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\
\end{array}\right)
\]
In [9]:
ss2
Out[9]:
\[
\left(
\begin{array}{rllrll|rll}
4\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&3\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\
-6\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&-5\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\
\hline
1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\
\end{array}\right)
\]
Aの固有値を調べてみると1と-2で，正のものを一つ含むので不安定である．
In [10]:
l, _ = np.linalg.eig(A)
l
Out[10]:
array([ 1., -2.])
一方，伝達関数を調べてみると，G1の極はAの固有値を二つとも保存しているが，G2の極は安定な固有値の方しか含まない．Aの不安定極は，出力行列C2によって隠されてしまっている．
In [11]:
G1 = ss2tf(ss1)
G1
Out[11]:
$$\frac{3}{s^2 + s - 2}$$
In [12]:
pole(G1)
Out[12]:
array([-2., 1.])
In [13]:
G2 = ss2tf(ss2)
G2
Out[13]:
$$\frac{1}{s + 2}$$
In [14]:
pole(G2)
Out[14]:
array([-2.])
コーヒーブレイク：安定性の用語に注意¶
積分器を含むシステム＝原点に極を持つ場合の挙動
In [15]:
G1=tf([1],[1,1,0])
G1
Out[15]:
$$\frac{1}{s^2 + s}$$
In [16]:
pole(G1)
Out[16]:
array([-1., 0.])
In [17]:
ss1 = tf2ss(G1)
ss1
Out[17]:
\[
\left(
\begin{array}{rllrll|rll}
-1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&-1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\
1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\
\hline
0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&-1\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}&0\phantom{.}&\hspace{-1em}&\hspace{-1em}\phantom{\cdot}\\
\end{array}\right)
\]
In [18]:
s, v = np.linalg.eig(ss1.A)
s
Out[18]:
array([ 0., -1.])
In [19]:
fig, ax = plt.subplots(figsize=(6, 6))
w = 1.5
Y, X = np.mgrid[-w:w:5j, -w:w:5j]
U = ss1.A[0,0]*X + ss1.A[0,1]*Y
V = ss1.A[1,0]*X + ss1.A[1,1]*Y
# 固有空間のプロット
vv = v * w * np.sqrt(2)
if s.imag[0] == 0 and s.imag[1] == 0: #固有値が複素数の場合はプロットできない
ax.plot([-vv[0,0],vv[0,0]], [-vv[1,0],vv[1,0]], c='k', ls='-.', lw=2)
ax.plot([-vv[0,1],vv[0,1]], [-vv[1,1],vv[1,1]], c='k', ls='--', lw=2)
ax.streamplot(X, Y, U, V, density=0.7, color='k', linewidth=0.5)
ax.set_xlim(-1.5,1.5)
ax.set_ylim(-1.5,1.5)
ax.set_xticks(np.arange(-1.5, 1.51, step=0.5))
ax.set_yticks(np.arange(-1.5, 1.51, step=0.5))
ax.grid(ls = ':')
plot_set(ax, '$x_1$', '$x_2$')
_, _, x = initial(ss1, X0=[1,-0.5], return_x=True)
ax.plot(x[:,0], x[:,1], lw=2, color='r',)
Out[19]:
[<matplotlib.lines.Line2D at 0x7f9b204cd220>]
In [20]:
fig, ax = plt.subplots(figsize=(6, 4.6))
x0 = [1,-0.5]
_, trange, x = initial(ss1, X0=x0, return_x=True) # シミュレーション
ax.plot(trange, x[:,0], lw=2, color='r', label='$x_1$')
ax.plot(trange, x[:,1], lw=2, color='b', label='$x_2$')
ax.axhline(0, color="k", linewidth=0.5)
ax.grid(ls = ':')
plot_set(ax,'t','x','best')
例5.2 Lyapunovの安定判別法¶
In [21]:
A = np.matrix([[1, 2],[-3, -4]])
B = np.matrix([[0],[1]])
C = np.matrix([1, 0])
D = np.matrix([0])
sys = ss(A, B, C, D)
print(sys)
A = [[ 1. 2.]
[-3. -4.]]
B = [[0.]
[1.]]
C = [[1. 0.]]
D = [[0.]]
In [22]:
Q = np.eye(A.shape[0])*12
Q
Out[22]:
array([[12., 0.],
[ 0., 12.]])
In [23]:
#関数の仕様上，A^Tを代入していることに注意
P = lyap(np.transpose(A),Q)
print(P)
[[27. 11.] [11. 7.]]
In [24]:
P*A + np.transpose(A)*P
Out[24]:
matrix([[-1.20000000e+01, -3.55271368e-15],
[-3.55271368e-15, -1.20000000e+01]])
固有値の確認
In [25]:
l, v = np.linalg.eig(A)
l
Out[25]:
array([-1., -2.])
※数式処理による理論解
In [26]:
import sympy
sympy.init_printing()
sympy.var('p11, p12, p22')
Asym=sympy.Matrix(A)
Psym=sympy.Matrix([[p11, p12],[p12, p22]])
Qsym=sympy.Matrix(Q)
Lyap = Psym*Asym + Asym.transpose()*Psym + Qsym
Lyap
Out[26]:
$\displaystyle \left[\begin{matrix}2 p_{11} - 6 p_{12} + 12.0 & 2 p_{11} - 3 p_{12} - 3 p_{22}\\2 p_{11} - 3 p_{12} - 3 p_{22} & 4 p_{12} - 8 p_{22} + 12.0\end{matrix}\right]$
In [27]:
eq1 = sympy.Eq(Lyap[0,0], 0)
eq2 = sympy.Eq(Lyap[0,1], 0)
eq3 = sympy.Eq(Lyap[1,1], 0)
sympy.solve([eq1, eq2, eq3], [p11, p12, p22])
Out[27]:
$\displaystyle \left\{ p_{11} : 27.0, \ p_{12} : 11.0, \ p_{22} : 7.0\right\}$
章末問題【２】も参照のこと．
例5.3 フィードバック系の内部安定性と極零相殺¶
数式処理による理論解析
In [28]:
s=sp.symbols('s')
K=(s-1)/(s+1)
P=1/(s-1)
K,P
Out[28]:
$\displaystyle \left( \frac{s - 1}{s + 1}, \ \frac{1}{s - 1}\right)$
In [29]:
Gyr=sp.factor(P*K/(1+P*K))
Gyr
Out[29]:
$\displaystyle \frac{1}{s + 2}$
In [30]:
Gyd=sp.factor(P/(1+P*K))
Gyd #不安定！
Out[30]:
$\displaystyle \frac{s + 1}{\left(s - 1\right) \left(s + 2\right)}$
In [31]:
Gur=sp.factor(K/(1+P*K))
Gur
Out[31]:
$\displaystyle \frac{s - 1}{s + 2}$
章末問題【３】も参照のこと．
例5.4(a) L(s)が積分器を持つ場合¶
In [32]:
L=tf([1],[1,0]) * tf([1],[1,1])**3
In [33]:
# 全体図(左)，拡大図（右）
fig, ax = plt.subplots(1,2,figsize=(12, 6))
x, y, _ = nyquist(L, logspace(-3,3,1000), plot=False)
for i in [0,1]:
ax[i].plot(x, y, c='k', ls='-', lw=2)
ax[i].plot(x, -y, c='k', ls=':', lw=2)
ax[i].grid(ls=':', lw=0.5)
ax[i].set_xlabel('Re')
ax[i].set_ylabel('Im')
ax[0].set_xlim(-5,15)
ax[0].set_ylim(-10, 10)
ax[1].set_xlim(-3,1)
ax[1].set_ylim(-2,2)
ax[1].scatter(0, 0, color='k')
ax[1].annotate('$O$', xy=(0.2, 0), size=10)
Out[33]:
Text(0.2, 0, '$O$')
In [34]:
# Bode線図（参考）
fig, ax = plt.subplots(2, 1, figsize=(8,7))
gain, phase, w = bode(L, logspace(-3,3), plot=False)
ax[0].semilogx(w, 20*np.log10(gain), color='k', lw=1)
ax[1].semilogx(w, phase*180/np.pi, color='k', lw=1)
bodeplot_set(ax)
ax[1].set_yticks([-360,-270,-180,-90,0])
fig.tight_layout()
L(s)が虚軸上に極を持つ場合¶
In [35]:
omega = 10.0
L=tf([omega],[1,0,omega**2]) * tf([1],[1,1])
In [36]:
# 全体図(左)，拡大図（右）
fig, ax = plt.subplots(1,2,figsize=(12, 6))
log_eps = 0.1
x, y, _ = nyquist(L, logspace(-3, 1-log_eps,1000), plot=False)
x2, y2, _ = nyquist(L, logspace(1+log_eps,3,1000), plot=False)
for i in [0,1]:
ax[i].plot(x, y, c='k', ls='-', lw=2)
ax[i].plot(x, -y, c='k', ls=':', lw=2)
ax[i].plot(x2, y2, c='k', ls='-', lw=2)
ax[i].plot(x2, -y2, c='k', ls=':', lw=2)
ax[i].grid(ls=':', lw=0.5)
ax[i].set_xlabel('Re')
ax[i].set_ylabel('Im')
# ax[0].set_xlim(-5,15)
# ax[0].set_ylim(-10, 10)
# ax[1].set_xlim(-3,1)
# ax[1].set_ylim(-2,2)
In [37]:
# Bode線図（参考）
fig, ax = plt.subplots(2, 1, figsize=(8,7))
gain, phase, w = bode(L, logspace(-3,3,1000), plot=False)
ax[0].semilogx(w, 20*np.log10(gain), color='k', lw=1)
ax[1].semilogx(w, phase*180/np.pi, color='k', lw=1)
bodeplot_set(ax)
ax[1].set_yticks([-90,-45,0])
fig.tight_layout()
/Users/minami/opt/anaconda3/lib/python3.8/site-packages/control/xferfcn.py:317: RuntimeWarning: divide by zero encountered in true_divide out[i][j] = (polyval(self.num[i][j], x_arr) / /Users/minami/opt/anaconda3/lib/python3.8/site-packages/control/xferfcn.py:317: RuntimeWarning: invalid value encountered in true_divide out[i][j] = (polyval(self.num[i][j], x_arr) /
章末問題¶
【１】 Routh-Hurwitzの方法による安定判別¶
In [38]:
import sympy as sp
s = sp.symbols('s')
a3,a2,a1,a0 = sp.symbols('a_3, a_2, a_1, a_0')
# b3,b2,b1 = sp.symbols('b_3, b_2, b_1')
# c3,c2,c1 = sp.symbols('c_3, c_2, c_1')
# d3,d2,d1 = sp.symbols('d_3, d_2, d_1')
sp.init_printing()
(1) Routhの方法¶
In [39]:
b1 = (a3*a2 - 1*a1) / a3
b2 = (a3*a0 - 1*0) / a3 # 係数がないところは0を埋める
b3 = (a3* 0 - 1*0) / a3 # 係数がないところは0を埋める
b1,b2,b3
Out[39]:
$\displaystyle \left( \frac{- a_{1} + a_{2} a_{3}}{a_{3}}, \ a_{0}, \ 0\right)$
In [40]:
c1 = (b1*a1 - a3*b2) / b1
c2 = (b1*0 - a3*b3) / b1
c3 = (b1*0 - a3*0 ) / b1
c1, c2, c3
Out[40]:
$\displaystyle \left( \frac{a_{3} \left(- a_{0} a_{3} + \frac{a_{1} \left(- a_{1} + a_{2} a_{3}\right)}{a_{3}}\right)}{- a_{1} + a_{2} a_{3}}, \ 0, \ 0\right)$
In [41]:
c1 = sp.factor(sp.expand(c1))
c1
Out[41]:
$\displaystyle \frac{a_{0} a_{3}^{2} + a_{1}^{2} - a_{1} a_{2} a_{3}}{a_{1} - a_{2} a_{3}}$
In [42]:
d1 = (c1*b2 - b1*c2) / c1
d2 = (c1*b3 - b1*c3) / c1
d1,d2
Out[42]:
$\displaystyle \left( a_{0}, \ 0\right)$
再評価
In [43]:
Routh = sp.Matrix([
[1, a2, a0],
[a3, a1, 0],
[b1, b2, b3],
[c1, c2, 0],
[d1, d2, 0]
])
Routh
Out[43]:
$\displaystyle \left[\begin{matrix}1 & a_{2} & a_{0}\\a_{3} & a_{1} & 0\\\frac{- a_{1} + a_{2} a_{3}}{a_{3}} & a_{0} & 0\\\frac{a_{0} a_{3}^{2} + a_{1}^{2} - a_{1} a_{2} a_{3}}{a_{1} - a_{2} a_{3}} & 0 & 0\\a_{0} & 0 & 0\end{matrix}\right]$
係数を代入して確認
In [44]:
values=[(a3,2), (a2,3), (a1,4), (a0, 5)]
In [45]:
Routh.subs(values)
Out[45]:
$\displaystyle \left[\begin{matrix}1 & 3 & 5\\2 & 4 & 0\\1 & 5 & 0\\-6 & 0 & 0\\5 & 0 & 0\end{matrix}\right]$
最左列をみると，２回符号が変わっているので，不安定（不安定極を二つもつ）ことがわかる．
(2) Hurwitzの方法¶
In [46]:
H4 = sp.Matrix([[a1, a3, 0, 0],[a0, a2, 1, 0],[0, a1, a3, 0], [0, a0, a2, 1]])
H3 = H4[:3,:3]
H2 = H4[:2,:2]
H1 = H4[:1,:1]
In [47]:
H1, H2, H3, H4
Out[47]:
$\displaystyle \left( \left[\begin{matrix}a_{1}\end{matrix}\right], \ \left[\begin{matrix}a_{1} & a_{3}\\a_{0} & a_{2}\end{matrix}\right], \ \left[\begin{matrix}a_{1} & a_{3} & 0\\a_{0} & a_{2} & 1\\0 & a_{1} & a_{3}\end{matrix}\right], \ \left[\begin{matrix}a_{1} & a_{3} & 0 & 0\\a_{0} & a_{2} & 1 & 0\\0 & a_{1} & a_{3} & 0\\0 & a_{0} & a_{2} & 1\end{matrix}\right]\right)$
In [48]:
sp.det(H1), sp.det(H2), sp.det(H3), sp.det(H4)
Out[48]:
$\displaystyle \left( a_{1}, \ - a_{0} a_{3} + a_{1} a_{2}, \ - a_{0} a_{3}^{2} - a_{1}^{2} + a_{1} a_{2} a_{3}, \ - a_{0} a_{3}^{2} - a_{1}^{2} + a_{1} a_{2} a_{3}\right)$
In [49]:
sp.Matrix([sp.det(H1), sp.det(H2), sp.det(H3), sp.det(H4)]).subs(values)
Out[49]:
$\displaystyle \left[\begin{matrix}4\\2\\-12\\-12\end{matrix}\right]$
(3) 確認¶
In [50]:
P = tf([1],[1,2,3,4,5])
P
Out[50]:
$$\frac{1}{s^4 + 2 s^3 + 3 s^2 + 4 s + 5}$$
In [51]:
pole(P)
Out[51]:
array([-1.28781548+0.85789676j, -1.28781548-0.85789676j,
0.28781548+1.41609308j, 0.28781548-1.41609308j])
【２】 Lyapunovの方法による安定判別¶
In [52]:
A = np.matrix([[-2, 1],[2, -3]])
B = np.matrix([[0],[1]])
C = np.matrix([1, 0])
D = np.matrix([0])
sys = ss(A, B, C, D)
print(sys)
A = [[-2. 1.]
[ 2. -3.]]
B = [[0.]
[1.]]
C = [[1. 0.]]
D = [[0.]]
In [53]:
Q = np.eye(A.shape[0])
Q
Out[53]:
array([[1., 0.],
[0., 1.]])
In [54]:
#関数の仕様上，A^Tを代入していることに注意
P = lyap(np.transpose(A),Q)
P*40
Out[54]:
matrix([[17., 7.],
[ 7., 9.]])
In [55]:
P*A + np.transpose(A)*P
Out[55]:
matrix([[-1.00000000e+00, 5.55111512e-17],
[ 5.55111512e-17, -1.00000000e+00]])
In [56]:
l, v = np.linalg.eig(A)
l
Out[56]:
array([-1., -4.])
※数式処理による理論解
In [57]:
import sympy
sympy.init_printing()
sympy.var('p, q, r')
Asym=sympy.Matrix(A)
Psym=sympy.Matrix([[p, q],[q, r]])
Qsym=sympy.Matrix(Q)
Lyap = Psym*Asym + Asym.transpose()*Psym + Qsym
Lyap
Out[57]:
$\displaystyle \left[\begin{matrix}- 4 p + 4 q + 1.0 & p - 5 q + 2 r\\p - 5 q + 2 r & 2 q - 6 r + 1.0\end{matrix}\right]$
In [58]:
eq1 = sympy.Eq(Lyap[0,0], 0)
eq2 = sympy.Eq(Lyap[0,1], 0)
eq3 = sympy.Eq(Lyap[1,1], 0)
sympy.solve([eq1, eq2, eq3], [p, q, r])
Out[58]:
$\displaystyle \left\{ p : 0.425, \ q : 0.175, \ r : 0.225\right\}$
【３】 フィードバック系の内部安定性と極零相殺¶
数式処理による理論解析
In [59]:
s=sp.symbols('s')
K=(s**2-3*s+2)/(s**2+2*s+1)
P=(s+1)/(s**2-4)
K,P
Out[59]:
$\displaystyle \left( \frac{s^{2} - 3 s + 2}{s^{2} + 2 s + 1}, \ \frac{s + 1}{s^{2} - 4}\right)$
In [60]:
K=sp.factor(K)
P=sp.factor(P)
K,P
Out[60]:
$\displaystyle \left( \frac{\left(s - 2\right) \left(s - 1\right)}{\left(s + 1\right)^{2}}, \ \frac{s + 1}{\left(s - 2\right) \left(s + 2\right)}\right)$
In [61]:
Gyr=sp.factor(P*K/(1+P*K))
Gyr
Out[61]:
$\displaystyle \frac{s - 1}{s^{2} + 4 s + 1}$
In [62]:
Gyd=sp.factor(P/(1+P*K))
Gyd
Out[62]:
$\displaystyle \frac{\left(s + 1\right)^{2}}{\left(s - 2\right) \left(s^{2} + 4 s + 1\right)}$
In [63]:
Gur=sp.factor(K/(1+P*K))
Gur
Out[63]:
$\displaystyle \frac{\left(s - 2\right) \left(s - 1\right) \left(s + 2\right)}{\left(s + 1\right) \left(s^{2} + 4 s + 1\right)}$
In [64]:
phi_FB=sp.factor(sp.denom(P)*sp.denom(K) + sp.numer(P)*sp.numer(K))
phi_FB
Out[64]:
$\displaystyle \left(s - 2\right) \left(s + 1\right) \left(s^{2} + 4 s + 1\right)$
【４】 $P(s)$が安定な場合のフィードバック制御系解析¶
(1) $r$から$y$までの伝達関数$G(s)$¶
In [65]:
import sympy as sp
sp.init_printing()
s=sp.Symbol('s')
Psym= sp.expand( 1 / ((s+1)**2 * (0.1*s+1)) )
Psym
Out[65]:
$\displaystyle \frac{1}{0.1 s^{3} + 1.2 s^{2} + 2.1 s + 1}$
In [66]:
Ksym=sp.Symbol('k')
Ksym
Out[66]:
$\displaystyle k$
In [67]:
Gsym = sp.simplify( Psym * Ksym / (1 + Psym * Ksym) )
Gsym
Out[67]:
$\displaystyle \frac{k}{k + 0.1 s^{3} + 1.2 s^{2} + 2.1 s + 1}$
(2) Bode線図とゲイン余裕¶
In [68]:
P = tf([1],[1,1])**2 * tf(1,[0.1, 1])
In [69]:
fig, ax = plt.subplots(2, 1, figsize=(4, 3.5))
k=1
gain, phase, w = bode(P*k, logspace(-2,3), plot=False)
ax[0].semilogx(w, 20*np.log10(gain), color='k', lw=1)
ax[1].semilogx(w, phase*180/np.pi, color='k', lw=1)
bodeplot_set(ax)
ax[0].semilogx(w, 0*w, c='c', ls='-', lw=1)
ax[1].semilogx(w, -180*np.ones_like(w), c='c', ls='-', lw=1)
ax[1].set_yticks([-270,-180,-90,0])
# ゲイン余裕，位相余裕，位相交差周波数，ゲイン交差周波数
gm, pm, wpc, wgc = margin(P*k)
GM = 20*np.log10(gm) #デシベル値に変換
ax[0].plot([wpc,wpc],[-GM,0], c='r', ls='-', lw=2)
ax[0].scatter(wpc,-GM, alpha=0.5, c='k')
ax[1].scatter(wpc,-180, alpha=0.5, c='k')
fig.tight_layout()
if (is_savefig):
fig.savefig(figpath+"ans/ch5_4_bode.pdf", transparent=True, bbox_inches="tight", pad_inches=0.0)