はじめに
著者が制御器の設計および実装を行う際、基本的には連続時間で設計⇒離散化して実装を行ってきました。なので離散時間状態方程式はその存在すらほぼ知らずにいたのですが、本業の横分野について勉強する中でどうもそういうのがあるらしい、というのに気付きました。
その際に、「ん、これって連続→離散化した場合とで何が違うの？」という点が引っ掛かりました。これら２つのアプローチについて、その差異を明確に書いたドキュメントが意外と無かったので勉強がてらまとめるものとします。
結論だけ言えば、サンプリング周期間において入力uが一定条件のもと、連続時間状態方程式の前進差分による離散化と、離散時間状態方程式の指数関数の１次近似は完全に一致するという結果が得られました。離散時間状態方程式、これで怖くない！
また本稿では不随して２つの考察を行っていますが、これもかなり興味深い内容になりました。うーん、制御工学は奥が深い。
なお下記、少し手探りしつつ書いている面があるので間違いがあれば指摘頂ければ幸いです。
連続時間状態方程式とその離散化
下記で与えられます。
\dot{x}(t)=Ax(t)+Bu(t)
上記を離散化する場合は、ラプラス変換→z変換(離散化)を行うことで行います。
ラプラス変換は下記。なおx(0)=0であるものとします。
sX(s)=AX(s)+BU(s)\\
X(s)=(sI-A)^{-1}BU(s)
次にz変換ですが、具体的な議論を行う目的で以降では行列A,Bの内容を明示するものとします。
議論の単純化のためモデル化の対象をRL回路とし、このときのA,Bは下記で与えられます。
電圧方程式より\, v=Ri+L \dot{i} \, \\
i=x,v=uとして式を整理することで\\A = -R/L\\B=1/L
もはやA,Bは行列ですらありませんが(笑)特に問題はありません。
z変換は前進差分にて行うものとします。Tはサンプリング周期です。
s = \frac{z-1}{T}
ラプラス変換後の式に代入することで下記式が得られます。
(z-1+\frac{R}{L}T)i[n]= \frac{1}{L}Tv[n] \\
i[n+1] = (1-\frac{R}{L}T) i[n]+\frac{1}{L}Tv[n]
#離散時間状態方程式
下記で与えられます。
x((n+1)T)=e^{AT}x(nT)+\int_{0}^{T}e^{A\tau}Bu(T-\tau)d\tau
入力uが、サンプリング周期Tの間に一定であるものとすると下記変形できます。この前提は非常に大事です。
when \,u(\tau)=u_{0} \\ \,
subject \,to \,\, nT<\tau<(n+1)T\\
x((n+1)T)=e^{AT}x(nT)+ \int_{0}^{T}e^{A\tau}Bu_0d\tau \\
=e^{AT}x(nT)+A^{-1}[e^{A \tau}Bu_{0}]^{T}_0\\
=e^{AT}x(nT)+A^{-1}(e^{AT}-I)Bu_{0}
「離散」状態方程式なのに連続→離散で得られた式とは到底かけ離れていますね。近づけるには、
指数関数の項を１次まで近似します。このことで下記式が得られます。
when \,e^{AT}\simeq I+AT\\
x((n+1)T)=(1+AT)x(nT)+TBu_{0}
A,BをRL回路のものに置き換え、またu0=u(nT)=v(nt)と書き換えることで下記与えられます。
i((n+1)T)=(1-\frac{R}{L} T)i(nT)+T\frac{1}{L}v(nT)
上式は(nT)を[n]に置き換えることで連続時間状態方程式の前進差分による離散化結果と一致します。やったね。
ところで① 後退差分だとどうなる？
離散化といっても色々ありますよね、後退、双一次etc...
双一次は明らかに違う結果になることが予測できますが、後退だとどうでしょう。試してみます。
後退差分におけるz変換の式は下記です。
s = \frac{1-z^{-1}}{T}
ラプラス変換後の連続時間状態方程式に上式を代入することで下記式が得られます。
(1-z^{-1}+ \frac{R}{L}T)i[n]= \frac{1}{L}Tv[n] \\
i[n] = \frac{1}{1+\frac{R}{L}T} i[n-1]+\frac{T}{L}\frac{1}{1+\frac{R}{L}T}v[n]
なんと、全然違う式が出てしまいました。このへんが離散化の嫌らしいところと言えます。
なお前進差分と後退差分とで、解析結果は大差ありません。グラフは上が２条件の解析結果を重ねたもので、下が２条件の差です。このため前進差分と後退差分はどっちがどっちだかの理解がいい加減であっても実用上は問題にならなかったりします。
参考までにMATLABのコードを示します。
close all;
clear all;
%% 入出力関係の定義
Len = 1000; %データ長の定義
v = ones(1,Len); %ステップ応答を定義
i1 = zeros(1,Len);%電流iの初期化 前進差分用
i2 = zeros(1,Len);%後退差分用
%% パラメータ設定
L= 0.01;
R = 0.2;
T = 1e-3;
%% 前進差分の場合の結果を計算
a1 = (1-R/L*T);
b1 = 1/L*T;
for n = 1:Len-1
i1(n+1) = a1*i1(n)+b1*v(n);
end
plot(i1(1:500));
hold on;
%% 前進差分の場合の結果を計算
a2 = 1/(1+R/L*T);
b2 = T/L*a2;
for n = 2:Len
i2(n) = a2*i2(n-1)+b2*v(n);
end
figure(1);
subplot(2,1,1);hold on;grid on;
plot(i1);plot(i2);
subplot(2,1,2);
plot(i1-i2);
hold on;grid on;
#ところで② 逆ラプラス変換との関係は？
上記ではz変換によって周波数領域→時間領域への変換を行っていますが、逆ラプラス変換によっても行うことが出来ます。この場合と離散時間状態方程式との差異についても、気になったので明らかにしてみます。
ラプラス変換後のRL回路の方程式は下記。
X(s)=(sI-A)^{-1}BU(s) \,より\\
I(s)=\frac{1}{Ls+R}V(s)
逆ラプラス変換を行う場合、あらかじめV(s)の内容を与えてあげる必要があります。ここではステップ応答を想定し、電圧は一定値Eであるものとします。
ラプラス変換一覧表より
V(s)=\frac{E}{s}
代入して式変形します。
I(s)=\frac{1}{Ls+R} \frac{E}{s}=\frac{E}{L}(\frac{1}{s+\frac{R}{L}} \frac{1}{s}) \\
=\frac{E}{L}\frac{L}{R}(\frac{1}{s}-\frac{1}{s+\frac{R}{L}})=\frac{E}{R}(\frac{1}{s}-\frac{1}{s+\frac{R}{L}})
ラプラス変換一覧表を右から左に読むことで逆ラプラス変換します。結果は以下。
i(t)=\frac{E}{R}(1-e^{-\frac{R}{L}t})
もう一方の離散時間状態方程式ですが、まず下記に式を再掲します。
x((n+1)T)=e^{AT}x(nT)+A^{-1}(e^{AT}-I)Bu_{0}
以下ではnTを始点0と定義、(n+1)Tは始点から任意時間t経過後と非常にトリッキーな置き換えを行います。またここでx(0)=0と定義します。さらにはステップ応答なのでu0は時間によらずEと置き換えることが出来ます。このとき下記式が得られます。
i(t)=-\frac{L}{R}(e^{-\frac{R}{L}t}-1)\frac{E}{L}=\frac{E}{R}(1-e^{-\frac{R}{L}t})
逆ラプラスの結果と一致しました。
上記から分かることは、
・連続時間状態方程式と離散時間状態方程式は本質的には等価で、式の形が違う
・連続状態方程式→離散状態方程式への置き換え時には、何の情報の欠落も発生しない
ということが言えます。
情報の欠落が発生するのは、連続時間状態方程式をz変換により離散化するとき、および離散時間状態方程式の指数関数項に対し近似を適用するときになります。
情報の欠落を抑制する手法としてz変換における双一次近似、指数関数の２次以上の近似が挙げられ、それぞれ前進差分や１次近似と比べると如実な性能改善が望めます。ただし当然ながら計算量とのトレードオフになるので、どちらを使うかは設計者の意図が問われます。さすがにイマドキは離散化による計算量のひっ迫は「基本的に」起きにくいので、とりあえず双一次変換や２次近似を使ってしまってもよい気もしますが・・・