Lyapunovの安定判別法 - 制御工学ブログ

状態方程式表現されたシステム

ここでは、制御系の安定性について説明したいと思います。まずは、以下のシステムが与えられているものとします。

\begin{equation} \dot{x}(t) = Ax(t) + Bu(t) \\ y(t) = Cx(t) \end{equation}

リアプノフの安定判別法

リアプノフの安定判別法について説明していきます。リアプノフの安定判別法は自律システムの安定性を解を求めることなく判別する有力な安定性判別手法になります。

非線形システムに対応できるなど、かなり有力なツールとなっています。基本的なアイディアは以下の通りです。

まずスカラ値関数 $V(x(t))$ を考えます。これは任意の $x(t)$ に対して常に正値をとる正定値関数です。例えば、

\begin{equation} V(x)=x^T P x\end{equation}

として $P$ を正定行列と設定してやるとこれは正定値関数となります。

このとき、 $V(x(t))$ の時間微分が常に負であれば $V(x(t))$ の値は時間の経過とともに常に減少していくことになります。これがリアプノフ関数を利用した安定性判別の基本的なアイディアになります。

線形システムに対するリアプノフ方程式

線形システム（自律系）に対してリアプノフ方程式は自律系の $A$ 行列および行列 $P, Q$ から構成される以下の式です。

\begin{equation} A^T P+ P A + Q = 0\end{equation}

この式を満たすような正定行列 $P$ および半正定行列 $Q$ の組が存在すれば、システムは漸近安定となります。もし、 $A$ の固有値の実部が全て負の場合は、適当に与えた $Q$ に対して、リアプノフ方程式を解き、得られた $P$ が正定行列となることが知られています。このとき、リアプノフ関数は

\begin{equation} V = x(t)^T P x(t)\end{equation}

と与えられ、その微分値は

$\dot{V} = x(t)^T (A^T P + PA) x(t)$

であり、さらに、リアプノフ方程式を満足する $Q$ は半正定であり

\begin{equation}\dot{V} = - x(t)^T Q x(t)\end{equation}

として $\dot{V}$ が常に非正となり単調減少します。

MATLABによる安定性解析の基本例

MATLABで安定性解析をするプログラム例を示します。

A = [-1 0 0; 0 -2 -1; 0 1 -2];       
Q = [1 0 0; 0 1 0; 0 0 1]; %単位行列
P = lyap(A', Q);

このとき、出力結果は次のようになります。

P = 
    0.5000         0         0
         0    0.2500         0
         0         0    0.2500

実際、得られた $P$ を用いて $A^T P + PA + Q$ を計算すると0になります。また、当然のように $P$ は正定行列となっています。

なお、不安定システムに対して上記の関数lyapを適用しても正定となるような $P$ は求まりません。

ここでは、 $n = 2$ のシステムにおいて行列 $A$ が次式で与えられ、 $u(t) = 0$ とした場合について図を示します。

\begin{equation} A = \begin{bmatrix}0&1\\-7.1 &-3.4 \end{bmatrix}\end{equation}

$Q = I$ と単位行列を与えてリアプノフ方程式を解いたとき、リアプノフ行列 $P$ は次式のように求まります。

\begin{equation}P = \begin{bmatrix}1.4306&0.0704\\0.0704 &0.1678 \end{bmatrix}\end{equation}

このとき、 $P$ は正定行列として求まっています。

\begin{equation} \dot{x} = A x\end{equation}

の解軌道を以下の図の破線で表します。また、実線でリアプノフ関数の等高線を表します。破線は常にリアプノフ関数値が小さくなる方向に進んでいることが確認できます。 $A$ 行列の与えられ方によってどの関数がリアプノフ関数になるかが変わってきます。

2次元システムにおける具体例と解軌道の可視化

ここでは、 $n = 2$ のシステムにおいて異なる特性を持つ4つのケースについて詳しく説明します。これらの例は、リアプノフ関数と解軌道の関係、固有値が与える動的特性の違い、そして座標変換の影響を理解するために設計されています。

ケース01: 標準的な減衰振動系

$A$ 行列が次式で与えられ、 $u(t) = 0$ とした場合について説明します。前節の $A$ を考えます。このとき、 $P$ は正定行列として求まります。

% ケース01: 標準的な減衰振動系
A = [0 1; -7.1 -3.4];
Q = eye(2);
P = lyap(A', Q);

% 固有値の計算
lambda = eig(A);
fprintf('固有値: %.4f ± %.4fi\n', real(lambda(1)), abs(imag(lambda(1))));

% P行列の表示
disp('P行列:');
disp(P);

% 検証: A'*P + P*A + Q = 0
residual = A'*P + P*A + Q;
fprintf('残差ノルム: %.2e\n', norm(residual));

固有値: $\lambda_{1,2} = -1.7 \pm 2.4495i$

特徴:

複素共役固有値を持つため、スパイラル状の軌道となる
実部が負( $-1.7$ )のため、減衰振動を示す
虚部( $\pm 2.4495$ )が振動周波数を決定
適度な減衰率と振動周波数を持つ標準的なケース

解軌道（破線）は常にリアプノフ関数の等高線（実線）に対して内側に向かって進み、リアプノフ関数値が小さくなる方向に進んでいることが確認できます。

ケース02: 緩やかな減衰振動系

$A$ 行列が次式で与えられる場合について説明します。

\begin{equation} A = \begin{bmatrix}-0.2&-0.01\\0.2 &-0.2 \end{bmatrix}\end{equation}

% ケース02: 緩やかな減衰振動系
A = [-0.2 -0.01; 0.2 -0.2];
Q = eye(2);
P = lyap(A', Q);

lambda = eig(A);
fprintf('固有値: %.4f ± %.4fi\n', real(lambda(1)), abs(imag(lambda(1))));
disp('P行列:');
disp(P);

固有値: $\lambda_{1,2} = -0.2 \pm 0.1414i$

特徴:

実部が小さい( $-0.2$ )ため、非常にゆっくりとした減衰
虚部も小さい( $\pm 0.1414$ )ため、振動周波数が低い
収束までに長時間を要する
リアプノフ関数の等高線は比較的円に近い形状

ケース01との比較: 減衰率が約1/8.5、振動周波数が約1/17と、どちらも大幅に小さくなっています。これにより、応答が非常に緩やかになることがわかります。

ケース03: 高周波振動系

$A$ 行列が次式で与えられる場合について説明します。

\begin{equation} A = \begin{bmatrix}-0.2&-5\\1 &-0.2 \end{bmatrix}\end{equation}

% ケース03: 高周波振動系
A = [-0.2 -5; 1 -0.2];
Q = eye(2);
P = lyap(A', Q);

lambda = eig(A);
fprintf('固有値: %.4f ± %.4fi\n', real(lambda(1)), abs(imag(lambda(1))));
disp('P行列:');
disp(P);

固有値: $\lambda_{1,2} = -0.2 \pm 2.2361i$

特徴:

実部はケース02と同じ( $-0.2$ )だが、虚部が大きい( $\pm 2.2361$ )
減衰率は同じでも、振動周波数が高い
細かい振動をしながらゆっくり収束する
リアプノフ関数の等高線は楕円形状が顕著

重要な示唆: このケースは、固有値の実部（減衰率）と虚部（振動周波数）が独立に変化できることを示しています。ケース02と03を比較することで、減衰率を変えずに振動周波数だけを変更した効果を理解できます。

ケース04: 相似変換されたシステム

ケース03のシステムに対して座標変換を施した場合について説明します。変換行列 $T_v$ を次のように設定します。

\begin{equation} T_v = \begin{bmatrix}1&3\\0 &0.4 \end{bmatrix}\end{equation}

変換後の $A$ 行列は次式で求められます。

\begin{equation} A = T_v^{-1} \begin{bmatrix}-0.2&-5\\1 &-0.2 \end{bmatrix} T_v \end{equation}

% ケース04: 相似変換されたシステム
Tv = [1 3; 0 0.4];
A_original = [-0.2 -5; 1 -0.2];
A = inv(Tv) * A_original * Tv;

Q = eye(2);
P = lyap(A', Q);

lambda = eig(A);
fprintf('固有値: %.4f ± %.4fi\n', real(lambda(1)), abs(imag(lambda(1))));
disp('変換後のA行列:');
disp(A);
disp('P行列:');
disp(P);

% ケース03のP行列と比較
P_case03 = lyap(A_original', Q);
fprintf('P行列の関係検証: norm(P - Tv^T * P_case03 * Tv) = %.2e\n', ...
        norm(P - Tv' * P_case03 * Tv));

固有値: $\lambda_{1,2} = -0.2 \pm 2.2361i$ （ケース03と同一）

特徴:

相似変換により $A$ 行列の要素は変化するが、固有値は不変
したがって、システムの安定性と動的特性は変わらない
リアプノフ関数の形状は座標系に依存して変化する
$P$ 行列は $\tilde{P} = T_v^T P_{03} T_v$ の関係を満たす

教育的意義: このケースは、座標変換がシステムの本質的な性質（固有値、安定性）を変えないことを示しています。異なる座標系で同じシステムを表現しても、リアプノフ関数による安定性判別は有効であることが確認できます。

4つのケースの比較分析

ケース	固有値	減衰率 (実部)	振動周波数 (虚部)	主な特徴
01	$-1.7 \pm 2.4495i$	-1.7	2.4495	標準的な減衰振動、適度な収束速度
02	$-0.2 \pm 0.1414i$	-0.2	0.1414	緩やかな減衰、低周波振動
03	$-0.2 \pm 2.2361i$	-0.2	2.2361	緩やかな減衰、高周波振動
04	$-0.2 \pm 2.2361i$	-0.2	2.2361	03の座標変換版、固有値は同一

これらの例から学べること

減衰率の影響（ケース01 vs 02, 03）: 固有値の実部が負で大きいほど、速く収束します。
振動周波数の影響（ケース02 vs 03）: 固有値の虚部が大きいほど、振動が速くなります。減衰率が同じでも振動周波数は独立に変化できます。
座標変換の不変性（ケース03 vs 04）: 相似変換によってシステム行列の表現は変わりますが、固有値と安定性は変わりません。
リアプノフ関数の汎用性: どのような安定なシステムに対しても、適切な正定行列 $P$ が存在し、リアプノフ関数を構成できます。
等高線と軌道の関係: 解軌道は常にリアプノフ関数の等高線に対して内側（関数値が小さくなる方向）に進みます。

統合解析プログラム

4つのケースをまとめて解析するMATLABプログラムを以下に示します。

% 4つのケースの統合解析
clear; clc;

% ケースの定義
cases = {
    [0 1; -7.1 -3.4], '01: 標準的な減衰振動系';
    [-0.2 -0.01; 0.2 -0.2], '02: 緩やかな減衰振動系';
    [-0.2 -5; 1 -0.2], '03: 高周波振動系';
};

% ケース04の追加
Tv = [1 3; 0 0.4];
A04 = inv(Tv) * [-0.2 -5; 1 -0.2] * Tv;
cases{4,1} = A04;
cases{4,2} = '04: 座標変換後（ケース03）';

% 各ケースの解析
fprintf('=================================================\n');
fprintf('リアプノフ安定性解析 - 4つのケース比較\n');
fprintf('=================================================\n\n');

for i = 1:4
    A = cases{i,1};
    fprintf('=== %s ===\n', cases{i,2});
    
    % 固有値の計算
    lambda = eig(A);
    fprintf('固有値: %.4f ± %.4fi\n', real(lambda(1)), abs(imag(lambda(1))));
    fprintf('減衰率（実部）: %.4f\n', real(lambda(1)));
    fprintf('振動周波数（虚部）: %.4f\n', abs(imag(lambda(1))));
    
    % リアプノフ方程式の解
    Q = eye(2);
    P = lyap(A', Q);
    fprintf('P行列:\n');
    disp(P);
    fprintf('P行列の最小固有値: %.6f （正定性確認）\n', min(eig(P)));
    
    % リアプノフ方程式の検証
    residual = A'*P + P*A + Q;
    fprintf('リアプノフ方程式の残差ノルム: %.2e\n', norm(residual));
    fprintf('\n');
end

% ケース03と04の関係確認
fprintf('=== ケース03とケース04の相似変換関係の検証 ===\n');
A03 = cases{3,1};
P03 = lyap(A03', eye(2));
P04 = lyap(A04', eye(2));
P04_transformed = Tv' * P03 * Tv;
fprintf('変換による誤差: %.2e\n', norm(P04 - P04_transformed));
fprintf('=================================================\n');