4. システムの安定判別

4-1. 安定なシステム

特性方程式が式(1)のとき、このシステムの安定判別を行え。$$s^5 +8s^4 + 25s^3 + 40s^2 + 34 s + 12=0 \;\;\; \cdots (1)$$

解答例：　式(1)は、安定であるための必要条件が成り立っている。つまり特性方程式の係数が全て「正」である。次に、ラウス表を作成する。
※ラウス表の作成方法に関しては、17. 安定判別 を参照願います。

\(s^5\)行	1	25	34
\(s^4\)行	8	40	12
\(s^3\)行	20	65/2	0
\(s^2\)行	27	12	0
\(s^1\)行	425/18	0
\(s^0\)行	12

ラウス表の第１列の要素が全て「正」であるから、このシステムは安定であると判定できる。また、式(1)の特性方程式は、式(2)と変形できる。その結果、特性根は、\(\lambda_1 = -1,\; \lambda_2 = -2,\; \lambda_3 = -3,\; \lambda_{4,5} = -1 \pm j\)となり、特性根の実部が負なので、システムは安定であることがわかる。$$s^5 +8s^4 + 25s^3 + 40s^2 + 34 s + 12=(s+1)(s+2)(s+3)(s+1-j)(s+1+j)=0 \;\;\; \cdots (2)$$

※Scilabによりラウス表を作成するスクリプトと実行結果を示す。

// 多項式の微分を計算する関数
function dpoly=polyder(p)
    n = length(p);
    if n == 1 then
        dpoly = 0; // 定数の微分は 0
    else
        dpoly = zeros(1, n - 1);
        for i = 1:n-1
            dpoly(i) = p(i) * (n - i);
        end
    end
endfunction

// ラウス表を作成するScilabスクリプト
function RouthArray=routh_hurwitz(coeffs)
    n = length(coeffs);
    cols = ceil(n / 2); // 列数を計算
    RouthArray = zeros(n, cols);

    // 最初の2行を初期化
    RouthArray(1, 1:length(coeffs(1:2:$))) = coeffs(1:2:$);
    RouthArray(2, 1:length(coeffs(2:2:$))) = coeffs(2:2:$);

    // 残りの行を計算
    for i = 3:n
        for j = 1:min(n - i + 1, cols)
            a = RouthArray(i - 2, 1);
            b = 0;
            d = 0;

            if j + 1 <= size(RouthArray, 2) then
                b = RouthArray(i - 2, j + 1);
                d = RouthArray(i - 1, j + 1);
            end

            c = RouthArray(i - 1, 1);

            if c == 0 then
                c = 1e-6;  // ゼロ回避のための微小値
            end

            RouthArray(i, j) = (c * b - a * d) / c;
        end

        // すべての要素が0の場合の特別処理
        if norm(RouthArray(i, :)) == 0 then
            aux = polyder(RouthArray(i - 1, :)); // 多項式微分
            RouthArray(i, 1:length(aux)) = aux;
        end
    end
endfunction

// 与えられた特性方程式の係数
coeffs = [1, 8, 25, 40, 34, 12]; // 5次の例

// ラウス表を計算して表示
RouthArray = routh_hurwitz(coeffs);
disp(RouthArray);

4-2. 減衰特性を指定する

特性方程式を式(3)としたとき、特性根が全て\(Re[\lambda] \le -1\)となるように、\(a,\;b\)の値を求めよ。$$s^3 + 4s^2 + a s + b=0 \;\;\; \cdots (3)$$

解答例：　特性方程式の特性根\(\lambda\)の全ての実部が負（\(\Re[\lambda] \leq -1\)）となるような \(a,\; b\) の範囲を求める。元の特性方程式は、$$s^3 + 4s^2 + a s + b = 0$$なので、実部が\(-1\)以下になる条件を求めるのに変数変換 を使う。\(s = w - 1 \quad \text{(実部を -1 にシフト)}\)とすればよい。$$(w-1)^3 + 4(w-1)^2 + a(w-1)+b=0 \\ w^3 + w^2 +(a-5)w + (3-a+b)=0$$を得るので、この特性方程式でラウス表を作成する。

\(w^3\)行	1	\(a-5\)
\(w^2\)行	1	\(3-a+b\)
\(w^1\)行	\(2a-b-8\)	0
\(w^0\)行	\(3-a+b\)

安定であるための必要条件は$$a-5 \gt 0 ,\quad 3-a+b \gt 0$$である。また、必要十分条件は、$$2a-b-8 \gt 0,\quad 3-a+b \gt 0$$である。以上より、\(a,\; b\)を満たすべき条件は、$$b \lt 2a-8 , \quad b \gt a-3$$となる。

※Scilabにより、\(a,\;b\)の範囲を求めてグラフ化するスクリプトと実行結果を図３に示す。

// Scilabスクリプト
// K の安定範囲を求める関数
function isStable=check_stability(K)
    // ラウス表の行を作成
    R1 = [1, K];
    R2 = [9, 60];
    R3 = [(9*K - 60) / 9, 0];
    R4 = [60, 0];

    // 安定性判定（すべての行の先頭要素が正）
    isStable = (R1(1) > 0) & (R2(1) > 0) & (R3(1) > 0) & (R4(1) > 0);
endfunction

// K の範囲を探索
K_vals = linspace(-50, 50, 500);
stable_K = [];

for i = 1:length(K_vals)
    K = K_vals(i);
    
    if check_stability(K) then
        stable_K = [stable_K; K];
    end
end

// 安定領域をプロット
scf();
clf();
plot(stable_K, zeros(stable_K), 'b.');
xlabel("K");
ylabel("");
title("安定となる K の範囲");

4-3. 制御定数\(K\)と制御系の安定性

特性方程式を式(4)とするとき、システムが安定となる\(K\)の範囲を求めよ。$$s^3+9s^2+Ks+60=0 \;\;\; \cdots (4)$$

解答例：　特性方程式よりラウス表を作成する。

\(s^3\)行	1	\(K\)
\(s^2\)行	9	60
\(s^1\)行	\((9K-60)/9\)	0
\(s^0\)行	60

ラウス表より、システムが安定となるには、$$\frac{9K-60}{9} \gt 0$$が必要十分条件なので、\(K\)の範囲は、$$K \gt \frac{60}{9}$$となる。

4-4. 制御定数\(K\)と制御系の安定性

解答例：　図４より、$$L(s) = \frac{K}{s(s^2+2s+4)} \;\;\; \cdots (5)$$なので、特性方程式は、$$1+L(s) = 0 \text{より、}s^3+2s^2+4s+K=0$$となる。この特性方程式についてラウス表を作成する。

\(s^3\)行	1	4
\(s^2\)行	2	\(K\)
\(s^1\)行	\((8-K)/2\)	0
\(s^0\)行	\(K\)

ラウス表から、安定であるための必要十分条件は、$$8-K \gt 0 , \quad K \gt 0$$となる。よって、フィードバック制御系が安定になる\(K\)の範囲は、$$0 \lt K \lt 8$$となる。

※Scilabにより、\(K\)の範囲を求めて、ナイキスト軌跡を表示するスクリプトと実行結果を図５，６，７，８に示す。

// Scilab スクリプト
// ラウス・フルビッツ法で安定性判定
function isStable=check_stability(K)
    // 特性方程式の係数
    coeffs = [1, 2, 4, K];

    // ラウス表の計算
    RouthArray = zeros(4, 2);
    RouthArray(1, :) = coeffs(1:2); 
    RouthArray(2, :) = coeffs(3:4);

    // 3行目と4行目の計算
    RouthArray(3, 1) = (RouthArray(2, 1) * RouthArray(1, 2) - RouthArray(1, 1) * RouthArray(2, 2)) / RouthArray(2, 1);
    RouthArray(3, 2) = 0;

    // 4行目の計算
    RouthArray(4, 1) = coeffs(4);
    RouthArray(4, 2) = 0;

    // 安定性条件（すべての行の先頭要素が正）
    isStable = and(RouthArray(:, 1) > 0);
endfunction

// K の範囲を探索
K_vals = linspace(-10, 50, 1000);
stable_K = [];

for i = 1:length(K_vals)
    K = K_vals(i);
    
    if check_stability(K) then
        stable_K = [stable_K; K];
    end
end

// ラウス表の結果を表示
disp("ラウス判別法による安定な K の範囲:");
if size(stable_K, 1) > 0 then
    disp([min(stable_K), max(stable_K)]);
else
    disp("安定な K の範囲が見つかりませんでした。");
end

// ナイキスト軌跡の描画
function draw_nyquist(K)
    s = %s;
    L = K / (s * (s^2 + 2 * s + 4));
    Ls = syslin('c', L);

    // 周波数応答とナイキスト軌跡
    clf();
    nyquist(Ls, 0.01, 100);
    xgrid();

    // -1点のプロット
    plot(-1, 0, 'ro');
    xlabel("Re(L(jω))");
    ylabel("Im(L(jω))");
    title("ナイキスト軌跡 (K = " + string(K) + ")");
    
    // 表示範囲を設定 (x軸, y軸とも ±2)
    ax = gca();                       // 現在の軸を取得
    ax.data_bounds = [-2, -2; 2, 2];  // 表示範囲を指定
endfunction

// 安定領域とナイキスト軌跡を表示
if size(stable_K, 1) > 0 then
   draw_nyquist(max(stable_K));
end
scf(1);
   draw_nyquist(4);
scf(2);
   draw_nyquist(10);

4-5. 位相余裕、利得余裕

4-4.（図４）のフィードバック制御系において、\(K = 6,8,10\)についての位相余裕、利得余裕を計算せよ。

解答例：　式(5)の\(L(s)\)で、\(s \to j \omega\)として、開ループ周波数伝達関数\(L(j\omega)\)を得る。$$L(j \omega) = \frac{K}{j \omega(-\omega^2 + j 2 \omega +4)}=\frac{-K}{2\omega^2 + j(\omega^3 - 4 \omega)} \\= \frac{-2K\omega^2 +jK(\omega^3 - 4\omega)}{(2\omega^2)^2 + (\omega^3 - 4\omega)^2} \;\;\; \cdots (6)$$ 利得が１となる角周波数は、$$(2 \omega^2)^2 + (\omega^3 - 4\omega)^2 = (-K)^2$$を解けば良い。\(\Omega = \omega^2\)とおくと、$$\Omega^3 - 4\Omega^2 + 16\Omega - K^2=0$$となる。\(K=6\)のとき、\(Omega = 2.835, \; 0.5824 \pm j3.5154\)なので、ゲイン交差角周波数（ゲインが１（0 dB))となる角周波数\(\omega_{cg}\)は、\(\omega_{cg} = 1.684\)[rad/s]である。この時の位相\(\theta_{\omega_{cg}}\)は、$$\theta_{\omega_{cg}} = \tan^{-1}\frac{\omega_{cg}^2 - 4}{-2 \omega_{cg}} = \tan^{-1}\frac{-1.165}{-3.368} = 160.92^{\circ}$$位相余裕は、$$\phi_m = \theta_{\omega_{cg}} - (-180^{\circ})=19.08^{\circ}$$である。
開ループ伝達関数の位相は、$$\theta(\omega) = \tan^{-1} \frac{\omega^2 - 4}{-2 \omega}$$なので、位相が\(-\pi\)となる角周波数は、\(\omega_{cp} =2.0\)である。利得は式(6)より、$$|L(j\omega)| = \frac{K}{\sqrt{(2\omega_{cp}^2)^2+(\omega_{cp}^3 - 4 \omega_{cp})^2}}$$となる。よって、利得余裕は、$$g_m = -20\log_{10} |L(j \omega)| = -20\log_{10} \left[\frac{K}{\sqrt{(2\omega_{cp}^2)^2+(\omega_{cp}^3 - 4 \omega_{cp})^2}}\right] = -20\log_{10} \frac{K}{8}$$となる。\(K=6\)のとき、\(g_m = 2.50\)[dB]となる。

※Scilabにより、\(K = 6,8,10\)についての位相余裕、利得余裕を計算するスクリプトと実行結果を図９に示す。\(K=8\)において、安定限界となることが分かる。

// Scilabスクリプト
// 位相余裕・利得余裕を計算する関数
function compute_margins(K)
    s = %s;
    L = K / (s * (s^2 + 2 * s + 4));
    sys = syslin('c', L);

    // 位相余裕と位相交差周波数
    [pm, wcp] = p_margin(sys);
    
    // 利得余裕とゲイン交差周波数
    [gm, wcg] = g_margin(sys);

    // 数値誤差対策の閾値
    threshold = 1e-6;

    // 利得余裕の dB 変換と誤差処理
    if abs(gm) < threshold then
        gm_str = "安定限界 (Critical Stability)";
    elseif gm > 0 then
        gm_db = 20 * log10(gm);
        gm_str = string(gm_db) + " dB";
    else
        gm_db = -20 * log10(abs(gm));  
        gm_str = string(gm_db) + " dB (不安定)";
    end

    // 結果を表示
    disp("========================");
    disp("K = " + string(K));
    disp("========================");
    disp("利得余裕 (Gain Margin) : " + gm_str);
    disp("位相余裕 (Phase Margin): " + string(pm) + "°");
    disp("ゲイン交差周波数 (wg) : " + string(wcg) + " rad/s");
    disp("位相交差周波数 (wp)  : " + string(wcp) + " rad/s");
endfunction

// K = 6, 8, 10 の場合の余裕を計算
K_values = [6, 8, 10];

for i = 1:length(K_values)
    compute_margins(K_values(i));
end