「カオス制御」ってかっこいいよね
はじめに
決定論的なダイナミクスに支配されながらも複雑な挙動を見せるシステムをカオスシステムといいます。カオスシステムの特徴の一つに初期値の僅かな違いがその後の時間発展に大きな影響を与える「初期値鋭敏性」があります。
例として、下図のSprott回路[1]を示します。図中の抵抗器とキャパシタはすべて同じ回路定数$R, C$を持ちます。
図中の3点の電位$x, y, z$を状態変数としたダイナミクスは以下の常微分方程式で支配されます。ここで、指数を含む項は$z$へのフィードバック経路にあるダイオードの特性です。
\begin{aligned}
\eta &:= \begin{bmatrix}
x \\ y \\ z
\end{bmatrix}\\
\dot{\eta} &= f(\eta) \\
& = \begin{bmatrix}
\frac{1}{CR}y \\ \frac{1}{CR}z \\ -\frac{1}{CR}x-10^{-12}\frac
{1}{C}\left(\mathrm{e}^{y/0.026}-1\right)-\frac{1}{CR}z
\end{bmatrix}
\end{aligned}
パラメータとなる回路定数を$R=1\ \mathrm{k\Omega},\ C=1\ \mathrm{\mu F}$としたとき、初期値$x=0.0, y=0.4, z=0.0$および$x=0.0, y=0.400000001, z=0.0$に対する時間応答と解軌道を求めると下図のようになり、時間経過とともに挙動が乖離することがわかります。
本稿では、この回路に制御入力$u$を追加した以下のシステムを対象とします。
\begin{aligned}
\dot{\eta} &= f(\eta, u) \\
& = \begin{bmatrix}
\frac{1}{CR}y \\ \frac{1}{CR}z \\ -\frac{1}{CR}x-10^{-12}\frac
{1}{C}\left(\mathrm{e}^{y/0.026}-1\right)-\frac{1}{CR}z-\frac{1}{CR}u
\end{bmatrix}
\end{aligned}
周期軌道の安定化制御
カオスシステムは初期値鋭敏性を持つため、初期値の微小変動が時間とともに拡大します。しかし前出の図のとおり、数周期のごく短期間であれば軌道の乖離は小さく留まる可能性があり、かつ決定論的なシステムであれば同一初期値に対する応答は確定的です。
従って、もしこの軌道上に点$\eta_f$が存在し、適当な周期$T$で
\eta_f(t) = \eta_f(t+T)
となるのであれば、カオスシステムの不規則な軌道は周期$T$の安定軌道に拘束されます。
つまり、軌道全体を制御せずとも、軌道が毎回$\eta_f$を通るように制御するだけで軌道を安定化できます。
Poincaré写像
$\eta_f$の解析のためのツールとして、Poincaré写像を導入します。Poincaré写像は周期軌道が適当な断面(Poincaré断面)を横切るとき、断面上の点の遷移を表す写像です[2]。
下図は前出のSprott回路において、$x=0$($y-z$平面)をPoincaré断面としたときの図示で、一度Poincaré断面を横切った軌道が次に再度横切る際、交点の座標の遷移を表す写像$F$がPoincaré写像です。
Poincaré断面は元の状態空間の部分空間であり、上記の場合断面上の点は
\xi:=\begin{bmatrix}
y \\ z
\end{bmatrix}
で表すことができます。またもとの連続時間システムはPoincaré写像により断面を横切る点のダイナミクスとして離散化されます。これを以下の離散時間システムとして表します。
\xi[k+1]=F(\xi[k], u[k])
ここで、Poincaré断面を前項の$\eta_f$を通るように設定すれば、定義よりこの点はPoincaré写像の平衡点
\xi_f(u)=F(\xi_f(u), u)
に一致します。
従って、平衡点$\xi_f$の安定化制御を行うことで軌道の安定化が達成されます。たとえば、$\xi_f$周りでの$F$を線形化し、$\xi_f$近傍で状態フィードバックを行うなどです。
OGY法[3]
OGY制御はOtt, Grebori, Yorkeにより提案されたカオス制御の一手法で、通常の線形化に基づく制御と同様に平衡点$\xi_f$周りでの線形化モデルを用いるものの、平衡点が鞍点(双極型不安定)であるときに、状態を安定部分空間に拘束することで安定化を達成します。これにより線形制御と比べ小さな制御入力で安定化を行うことができます。
以下でOGY法を導出します。まず、状態変数が安定部分空間に属するための条件を導き、次にその条件を満足するための制御入力の決定方法を示します。
状態変数が安定部分空間に属するための条件
2次元のPoincaré写像が$u=u_f$のとき$\xi=\xi_f(u_f)$に鞍点を持つと仮定します。
\xi_f(u_f)=F(\xi_f(u_f), u_f)
$F$に対し$\xi_f$まわりの線形近似モデルを
\begin{aligned}
\tilde{\xi}[k+1]&=A\tilde{\xi}[k]+B\tilde{u}[k]\\
\tilde{\xi}&:=\xi-\xi_f\\
\tilde{u}&:=u-u_f
\end{aligned}
とおきます。$A,B$は$\xi_f$周辺の入出力データを用いることもできますし、$F$が明らかであれば
A=\frac{\partial F}{\partial \xi}(\xi_f, u_f),\ B=\frac{\partial F}{\partial u}(\xi_f, u_f)
と求めることもできます。鞍点を仮定したため、$A$の固有値はSchurの意味での安定固有値$\lambda_s\in\mathbb{R}, \lvert \lambda_s\rvert<1$と不安定固有値$\lambda_u\in\mathbb{R}, \lvert \lambda_u\rvert>1$の組み合わせとなります。対応する正規化固有ベクトルを$e_{s}, e_{u}\in\mathbb{R}^{2}$とおき、$e_s, e_u$をまとめた行列
E:=\begin{bmatrix}
e_{s}& e_{u}
\end{bmatrix}\in\mathbb{R}^{2\times 2}
と$E$の逆行列となる$2\times 2$行列を
V:=E^{-1}=\begin{bmatrix}
v_s^\mathsf{T}\\v_u^\mathsf{T}
\end{bmatrix}\\
v_s, v_u\in\mathbb{R}^{2}
と定義すると、
VE=\begin{bmatrix}
v_s^{\mathsf{T}}e_s & v_s^{\mathsf{T}}e_u\\
v_u^{\mathsf{T}}e_s & v_u^{\mathsf{T}}e_u
\end{bmatrix}
=\begin{bmatrix}
1&0\\0&1
\end{bmatrix}
が成り立ちます。従って、$v_u$と$e_s$は直交しており、$\tilde{\xi}[k+1]$が$e_s$上にあることを
v_{u}^{\mathsf{T}}\tilde{\xi}[k+1]=0
と表現できます。
制御入力の決定
$A$の固有値分解
A=E\begin{bmatrix}
\lambda_s&0\\
0&\lambda_u
\end{bmatrix}V=\begin{bmatrix}
e_s & e_u
\end{bmatrix}\begin{bmatrix}
\lambda_s&0\\
0&\lambda_u
\end{bmatrix}\begin{bmatrix}
v_s^\mathsf{T}\\v_u^\mathsf{T}
\end{bmatrix}=e_s\lambda_sv_s^{\mathsf{T}}+e_u\lambda_uv_u^{\mathsf{T}}
および状態方程式をもとに計算を続けると、$\tilde{\xi}[k+1]$を$e_s$上に配置するための制御入力が
\begin{aligned}
v_{u}^{\mathsf{T}}\tilde{\xi}[k+1]&=v_{u}^{\mathsf{T}}\left(A\tilde{\xi}[k+1]+B\tilde{u}[k]\right)\\
&=v_{u}^{\mathsf{T}}\left(\left(e_s\lambda_sv_s^{\mathsf{T}}+e_u\lambda_uv_u^{\mathsf{T}}\right)\tilde{\xi}[k+1]+B\tilde{u}[k]\right)\\
&=\left(v_{u}^{\mathsf{T}}e_s\lambda_sv_s^{\mathsf{T}}+v_{u}^{\mathsf{T}}e_u\lambda_uv_u^{\mathsf{T}}\right)\tilde{\xi}[k+1]+v_{u}^{\mathsf{T}}B\tilde{u}[k]\\
&=\lambda_uv_u^{\mathsf{T}}\tilde{\xi}[k+1]+v_{u}^{\mathsf{T}}B\tilde{u}[k]\\
&=0\\
\Rightarrow \tilde{u}[k]&=-\frac{\lambda_u v_u^\mathsf{T}}{v_u^\mathsf{T}B}\tilde{\xi}[k]
\end{aligned}
と求まります。線形近似が有効なのは$\xi_f$近傍においてのみであるため、ある閾値$\epsilon$をもとに
\tilde{u}[k]=\left\{
\begin{array}{cl}
-\frac{\lambda_u v_u^\mathsf{T}}{v_u^\mathsf{T}B}\tilde{\xi}[k]&\left(\left\lVert\tilde{\xi}\right\rVert<\epsilon\right)\\
0&(\mathrm{otherwise})
\end{array}
\right.
と制御します。ただし、$v_u^\mathsf{T}B\neq0$を仮定しました。
設計例
冒頭のSprott回路について、MATLABにてOGY法の設計を行いシミュレーションで検証します。
Poincaré写像の平衡点の発見
まずはPoincaré写像の平衡点$\xi_f$を発見する必要があります。このため、まずは$u\equiv0$、初期状態$x=0, y=0.4, z=0$のシミュレーション結果からPoincaré断面をまたぐデータのみを抽出し、次にPoincaré断面をまたいだ際の状態変化量$\xi[k+1]-\xi[k]$が最も小さいものを平衡点の候補とします。
下図は$\xi[k+1]-\xi[k]$をquiverでプロットしたもので、今回は赤枠の範囲で$\xi_f$を探索します。
赤枠内で初期値を変化させながら$F$の位相平面をプロットしたものが下図です$(u\equiv 0)$。きれいな鞍点が得られていますので、今回はこの平衡点を制御したいと思います。
この点の周辺で初期値・入力を変化させながら入出力データを取得し、以下の正規方程式をもとに最小二乗法で線形近似系のパラメータを推定しました。
\begin{bmatrix}
y[k+1] \\ z[k+1]
\end{bmatrix}
=
\begin{bmatrix}
a_{yy} & a_{yz}\\
a_{zy} & a_{zz}
\end{bmatrix}
\begin{bmatrix}
y[k] \\ z[k]
\end{bmatrix}
+
\begin{bmatrix}
b_y \\ b_z
\end{bmatrix}
u[k]
+
\begin{bmatrix}
c_y \\ c_z
\end{bmatrix}
\begin{bmatrix}
y[1] \\ y[2] \\ \vdots \\ y[K] \\ z[1] \\ z[2] \\ \vdots \\ z[K]
\end{bmatrix}
=
\begin{bmatrix}
H & O\\
O & H
\end{bmatrix}
\begin{bmatrix}
a_{yy} \\ a_{yz} \\ b_{y} \\ c_{y} \\ a_{zy} \\ a_{zz} \\ b_{z} \\ c_{z}
\end{bmatrix}\\
H:=\begin{bmatrix}y[0] & z[0] & u[0] & 1\\
y[1] & z[1] & u[1] & 1\\
\vdots & \vdots& \vdots& \vdots&\\
y[K-1] & z[K-1] & u[K-1] & 1
\end{bmatrix}
平衡点の座標は$c_y, c_z$から
\xi_f=\left(I-A\right)^{-1}
\begin{bmatrix}
c_y\\c_z
\end{bmatrix}
で求まります。ここでは
A=
\begin{bmatrix}
-2.9902 & -2.2062\\
-1.3872 & -1.0143
\end{bmatrix},
B = \begin{bmatrix}
0.8284\\-0.3016
\end{bmatrix}\\
\xi_f =
\begin{bmatrix}
-0.8023 \\ -0.3806
\end{bmatrix}\\
を得ました。本システムにOGY法を適用した制御結果を以下に示します。ここでは制御実行閾値$\epsilon = 0.08$、初期状態$x=0, y=0.4, z=0$を設定しています。
一定時間経過後、安定軌道への拘束が成功しています。一度制御が有効になったあとに再度無効になっていたり、定常状態でもやや変動が見られますが、これはPoincaré断面上の点ではなく、Poincaré断面をまたいだ次の点を用いて制御しているため誤差が生じているようです。
MATLABコード
システム同定および制御シミュレーションに用いたMATLABコードをいかに示します。
システム同定
clear;
close all;
% circuit constants
C = 1e-6;
R = 1e3;
% dynamics
f = @(t, x, u) [
1/C/R*x(2);
1/C/R*x(3);
-1/C/R*x(1)-1e-12/C*(exp(x(2)/0.026)-1)-1/C/R*x(3)-1/C/R*u
];
% initial value set
y0 = -0.8 + linspace(-0.025, 0.025, 20);
z0 = -0.3825 + linspace(-0.025, 0.025, 20);
% simulation length for one period
T = 0.02;
% simulation setting
odeopt = odeset("RelTol", 1e-9, "AbsTol", 1e-9);
% input value set
u = linspace(-3e-2, 3e-2, 10);
% identification data
idd_u = zeros(length(u)*length(y0)*length(z0), 1);
idd_y = zeros(length(u)*length(y0)*length(z0), 1);
idd_z = zeros(length(u)*length(y0)*length(z0), 1);
idd_y_next = zeros(length(u)*length(y0)*length(z0), 1);
idd_z_next = zeros(length(u)*length(y0)*length(z0), 1);
% simulation
idx_data = 1;
for i_u = 1:length(u)
for i_y = 1:length(y0)
for i_z =1:length(z0)
x0 = [0; y0(i_y); z0(i_z)];
[t, sol] = ode45(@(t, x) f(t, x, u(i_u)), [0, T], x0, odeopt);
% find index of Poincare section traverse
i_cross = find(sol(1:end-1,1) > 0 & sol(2:end,1) < 0, 1);
% store data
idd_u(idx_data) = u(i_u);
idd_y(idx_data) = y0(i_y);
idd_z(idx_data) = z0(i_z);
idd_y_next(idx_data) = sol(i_cross, 2);
idd_z_next(idx_data) = sol(i_cross, 3);
idx_data = idx_data + 1;
end
end
end
% identification by least squares
theta = [
idd_y idd_z idd_u ones(size(idd_u)) zeros(length(idd_u), 4);
zeros(length(idd_u), 4) idd_y idd_z idd_u ones(size(idd_u));
] \ [idd_y_next; idd_z_next];
A = [
theta(1) theta(2);
theta(5) theta(6)
];
B = [
theta(3);
theta(7);
];
xi_f = (eye(2)-A) \ [theta(4); theta(8)];
制御シミュレーション
clear;
close all;
% circuit constants
C = 1e-6;
R = 1e3;
% dynamics
f = @(t, x, u) [
1/C/R*x(2);
1/C/R*x(3);
-1/C/R*x(1)-1e-12/C*(exp(x(2)/0.026)-1)-1/C/R*x(3)-1/C/R*u
];
% load identification result
load("fp_local_model.mat");
% design feedback gain
[E, Lambda] = eig(A);
E = normalize(E, 1, "norm");
V = inv(E);
lambda_u = Lambda(1,1);
vu = V(1,:)';
g = -(lambda_u*vu') / (vu'*B);
% control threshold
epsilon = 8e-2;
% simulation length for one period
T = 0.02;
% simulation setting
odeopt = odeset("RelTol", 1e-9, "AbsTol", 1e-9);
% number of simulation period
K = 100;
% full state(3 dim)
eta = cell(1, K+1);
% state on Poincare section
xi = zeros(2, K);
eta{1} = [0; 0.4; 0];
u = zeros(1, K);
hp = plot3(nan, nan, nan);
xlabel('x');
ylabel('y');
zlabel('z');
% simulation
for k = 1:K-1
xi(:,k) = eta{k}(2:3, end);
% OGY method
if norm(xi(:,k) - xi_f) < epsilon
u(k) = g*(xi(:,k) - xi_f);
s_title = "CONTROL ON";
else
u(k) = 0;
s_title = "CONTROL OFF";
end
u(k) = 0;
[t, sol] = ode45(@(t, x) f(t, x, u(k)), [0, T], eta{k}(:,end), odeopt);
% find index of Poincare section traverse
i_cross = find(sol(1:end-1,1) > 0 & sol(2:end,1) < 0, 1);
eta{k+1} = sol(1:i_cross+1, :)';
hp.XData = eta{k}(1,:);
hp.YData = eta{k}(2,:);
hp.ZData = eta{k}(3,:);
title(s_title);
exportgraphics(gcf, "control_result_wocon.gif", "Append", true);
end
参考文献
[1] J. C. Sprott, "A New Chaotic Jerk Circuit"
[2] 機械工学辞典 ポアンカレ写像
[3] 山本 茂, 遅延フィードバックによるカオス制御