Python で見る古典制御

    numpy, scipy.linalg, scipy.integrate.odeint, matplotlib.pyplot を使
    うため、まず import する。
  
    >>> import numpy
    >>> import scipy.linalg
    >>> import matplotlib.pyplot as plt
    >>> from scipy.integrate import odeint
   
□ ２次系の過渡応答

   ２次系

   y''(t) + 2ζω y'(t) + ω^2 y(t) = x(t)

   で ζ = 0.2, ω = 1, x(t) = 1 (単位ステップ入力)の場合の過渡応答を求める．

   y1' = y2
   y2' = -2*0.2*y2 - y1 + 1

    >>> def ydot(y, t):
            return (y[1], -2*0.2*y[1] - y[0] + 1)
    >>> y0 = (0, 0)
    >>> t = numpy.linspace(0, 14, 1000)
    >>> y = odeint(ydot, y0, t)
    >>> plt.plot(t, y[:,0])
    >>> plt.grid('on')
    >>> plt.show()

□ ベクトル軌跡（ナイキスト線図）

   伝達関数が
                   1
     G(s) = -----------------
            (1+s)(1+2s)(1+3s)
   で与えられる制御系のベクトル軌跡を描く．

    >>> def G(s):
            return 1 / ((1 + s)*(1 + 2 * s)*(1 + 3 * s))
    >>> w = numpy.logspace(-2, 4, 300)
    >>> s = map(G, 1j * w)
    >>> plt.plot(numpy.real(s), numpy.imag(s))
    >>> plt.grid('on')
    >>> plt.show()

□ ボード線図

   伝達関数が
               (1+3s)
     G(s) = ------------
            (1+s)(1+2s)
   で与えられる制御系のボード線図を描く．

    >>> def G(s):
            return (1 + 3 * s) / ((1 + s)*(1 + 2 * s))
    >>> w = numpy.logspace(-2, 1, 100)
    >>> s = map(G, 1j*w)

    ゲイン線図
    >>> plt.semilogx(w, 20*numpy.log10(numpy.abs(s)))
    >>> plt.grid('on')
    >>> plt.show()

    位相線図
    >>> plt.semilogx(w, 180.*numpy.angle(s)/numpy.pi)
    >>> plt.grid('on')
    >>> plt.show()

     位相が 180(-180)度を越えるとグラフが不連続になることがあるので注意が
     必要である．

□ 安定性の判別

   特性方程式が
     4s^3 + 6s^2 + 6s + 2 = 0
   で与えられる制御系の安定判別を行う．

    >>> p = numpy.array([4, 6, 6, 2])
    >>> numpy.roots(p)

    r の実部が全て負なので，制御系は安定である．

   特性方程式が
     s^4 + s^3 + 6s^2 + 6s + 2 = 0
   で与えられる制御系の安定判別を行う．

    >>> p = numpy.array([1, 1, 6, 6, 2])
    >>> numpy.roots(p)

    r の実部に正の値を取るものがあるので，制御系は不安定である．

□ 根軌跡

   一巡伝達関数が
                     K
     G(s)H(s) = ------------
                 s(s+1)(s+2)
   で与えられる制御系の根軌跡を描く．

   特性方程式は
     s^3 + 3s^2 + 2s + K = 0
   で与えられる．

    >>> def f(k):
            return numpy.roots(numpy.array([1, 3, 2, k]))
    >>> k = numpy.linspace(0, 30, 61)
    >>> nc = numpy.size(k.shape)
    >>> ran = range(nc)
    >>> r = f(k[0])
    >>> for i in ran[1:]:
            r = numpy.concatenate((r, f(k[i])))
    >>> plt.plot(numpy.real(r), numpy.imag(r), 'x')
    >>> plt.grid('on')
    >>> plt.show()