Maxima で見る古典制御

□ ラプラス変換

   デルタ関数

    (%i ) laplace(delta(t), t, s);

   単位ステップ関数

    (%i ) laplace(1, t, s);

   三角関数

    (%i ) laplace(sin(t), t, s);

   指数関数とランプ関数

    (%i ) laplace(%e^(-a*t)*t^2, t, s);

□ ラプラス変換の部分分数展開

    (%i ) G: 1/((s+2)^2*(s+3));
    (%i ) partfrac(G, s);

    (%i ) G: 1/(s*((s+2)^2+1));
    (%i ) partfrac(G, s);

      複素極までは展開されないことに注意する必要がある．複素極に展開した時の係数
      はたとえば次のようにすれば求められる．

    (%i ) r: solve(denom(G) = 0, s);
    (%i ) c[1]: ratsimp(G*(s-ev(s,r[i])));
    (%i ) ratsimp(ev(c[1], r[1]));

      次のようにして全ての係数を求めることもできる．

    (%i ) for i:1 thru length(r) do 
            (c[i]: ratsimp(G*(s-ev(s,r[i]))),
             print(r[i], ":", ratsimp(ev(c[i],r[i]))));

      maxima 5.13.0 では上記ではなく，次のように計算する必要がある．

    (%i ) r: solve(denom(G) = 0, s);
    (%i ) for i:1 thru length(r) do 
            (c[i]: num(G)/(gfactor(denom(G))/(s-ev(s,r[i]))),
             print(r[i], ":", rectform(ev(c[i],r[i]))));
    

□ 逆ラプラス変換

    (%i ) ilt(1/(s+a), s, t);

    (%i ) assume(w>0);
    (%i ) ilt((s+a)/((s+a)^2+w^2), s, t);
    (%i ) forget(w>0);

    (%i ) ilt(1/((s+2)^2*(s+3)), s, t);

    (%i ) ilt(1/(s*((s+2)^2+1)), s, t);

□ 周波数伝達関数の計算

   伝達関数が
              K
     G(s) = ------
            (1+Ts)
   で与えられる１次遅れ要素の周波数伝達関数と，その極座標表示を求める．

    (%i ) G(s) := K/(1+T*s);
    (%i ) G(%i*w);
    (%i ) realpart(G(%i*w)));
    (%i ) imagpart(G(%i*w));
    (%i ) assume(K>0);
    (%i ) abs(G(%i*w));
    (%i ) carg(G(%i*w));
    (%i ) forget(K>0);

□ ２次系の過渡応答

   ２次系

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

   で ζ = 2/10, ω = 1, x(t) = 1 (単位ステップ入力)の場合の過渡応答を求める．
   このとき，伝達関数 G(s) は 
                                       1
                              -----------------
                                 2          
                               s   + 2 ζ s + 1
   で，ステップ入力のラプラス変換は 1/s で与えられる．

    (%i ) zeta: 2/10;
    (%i ) g: 1/(s^2+2*zeta*s+1);
    (%i ) define(f(t), ilt(g * (1/s), s, t));
    (%i ) plot2d(f, [t,0,30], [GNUPLOT_PREAMBLE, "set grid"]);

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

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

    (%i ) g(s):=1/((1+s)*(1+2*s)*(1+3*s));
    (%i ) define(x(w), ratsimp(realpart(g(%i*w))));
    (%i ) define(y(w), ratsimp(imagpart(g(%i*w))));
    (%i ) plot2d([parametric,x,y],[w,0,1000],[nticks,300]);

□ ボード線図

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

    (%i ) g(s):=(1+3*s)/((1+4*s)*(1+2*s));
    (%i ) define(gain(w), 20*log(ratsimp(abs(g(%i*w))))/log(10));
    (%i ) plot2d(gain, [w,0.01,10],[GNUPLOT_PREAMBLE, "set log x; set grid"]);
    (%i ) define(angle(w), (180/%pi)*atan(ratsimp(imagpart(g(%i*w))/realpart(g(%i*w)))));
    (%i ) plot2d(angle, [w,0.01,10],[GNUPLOT_PREAMBLE, "set log x;set grid"]);

□ 安定性の判別

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

    (%i ) expand(solve(4*s^3 + 6*s^2 + 6*s + 2 = 0, s));

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

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

    (%i )float(realpart(float(solve(s^4 + s^3 + 6*s^2 + 6*s + 2 = 0, s))));

   次の式でも求まるが，計算は遅い．
   
    (%i )float(realpart(solve(s^4 + s^3 + 6*s^2 + 6*s + 2 = 0, s)));

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