外乱オブザーバと未知入力オブザーバの違いについて


# はじめに
制御対象に印可される外乱を推定する手法として、外乱オブザーバが広く知られています。外乱オブザーバの応用事例は多く、特にモータ、ロボットなどのメカトロ制御を中心として様々な場面で活躍しています。

一方、似たような名前の推定手法として、未知入力オブザーバがあります。こちらの手法は外乱オブザーバに比べて文献が少なく、両者を比較している内容をあまり見かけなかったためこの記事を作成しました。特に未知入力オブザーバに重きをおいて基本的な考え方を述べた後、MATLAB/Simulinkによるシミュレーションで動きを確認します。
外乱オブザーバは推定した外乱を抑制することに用いることも多いですが、抑制(相殺)させるにあたって両者でやることは変わらないので、本記事では推定の様子のみに着目します。

以下、未知入力オブザーバをUIO(Unknown Input Observer)と表記します。
UIOは、入力端に印可された未知信号の影響を受けずに、制御対象の状態変数を推定することを目的としたものです。状態変数の推定を行うことで、入力端への未知信号も間接的に推定することができます。外乱オブザーバは考え方によっては出力端外乱も推定できますが、あまり一般的な使われ方ではないように思えます。

この記事では入力端への未知信号をすべて外乱、と呼ぶことにして、入力端の外乱の推定の様子を比較します。
UIOの文献はこちら[^1]、外乱オブザーバの文献はこちら[^2]です。細かい議論は飛ばすので、正確な内容は元文献を参照ください。

# 制御対象
制御対象は下式の線形の連続時間状態方程式であらわされるものとします
(この記事では制御まで行わないので、観測対象といった方がいいかもしれませんが)。
```math
\begin{align}
 \dot{x}(t) &= A_c x(t)+ B_c u(t) + E_c d(t) \\
y(t) &= C x(t) 
\end{align}
```

両手法ともに連続時間系としての適用を考えることができますが、後述の適用条件を緩和するため、基本的にはUIOは離散時間系で構築することが多いようです。
上記連続時間系に対し、サンプリング周期$T$内で入力$u$、外乱$d$がともに一定と仮定すると(いわゆるゼロ次ホールド)、下記の離散時間系が得られます。
```math
\begin{align}
 x_{k+1} &= A x_k + B u_k + E d_k \\
y_k &= C x_k
\end{align}
```
ここで、状態$x_k \in \mathbb{R}^n$、入力$u_k \in \mathbb{R}^p$、外乱$d_k \in \mathbb{R}^r$、出力$y_k \in \mathbb{R}^m$です。各係数行列$A$、$B$、$E$の導出方法は省略しますが、各変数のサイズに合った行列とします。

# 未知入力オブザーバの構成
結論を先に言ってしまうと、UIOは観測値を時間軸で2点取ってくることで、外乱の影響を完全に無くすことを目的としています。外乱の影響を除去した結果、外乱の復元も行うことができます。以下でその詳細を見ていきます。

下記の同一次元オブザーバを構成することを考えます。
執筆者の造語ですが、後述の少し異なったUIOと区別するため、これを便宜上シングルレートUIOと呼ぶことにします。
```math
\begin{align}
 \hat{x}_{k+1} &= (A- L_1 CA - L_2 C) \hat{x}_k + (B-L_1 C B) u_k + L_1 y_{k+1} + L_2 y_k \\
\end{align}
```
ここで、$\hat{x}_k \in \mathbb{R}^n$は推定状態、$L_1$、$L_2$は適当なサイズのオブザーバゲイン行列です。

右辺の添え字が通常の状態オブザーバと異なり、$y_{k+1}$が出てきています。いきなり出てきた$L_1 CA$の項や$L_1 CB$の項の意義がこの段階では伝わりづらいですが、$y\_{k+1}$の存在に合わせていると考えると合点がいくはずです。
通常のオブザーバはステップ$k$にてステップ$k+1$の予測の式となりますが、この場合はステップ$k+1$にてステップ$k+1$の$x$の推定を行います。つまり通常のオブザーバは一手先の予測式だったのに対し、未知入力オブザーバは現在の状況を推定するかたちとなります(現在の状況、を演算しようとした段階ですでに微妙に過去になってるわけですが)。
実装上はすべての$k$を-1して$\hat{x}\_k = \hat{A} \hat{x}\_{k-1} + \hat{B} u\_{k-1}+ L_1 y\_{k} + L_2 y\_{k-1}$のように考えるのがいいと思います。

このオブザーバに対し、推定誤差$e_{k+1} = x\_{k+1} - \hat{x}\_{k+1}$と置くと、推定誤差ダイナミクスは下記のようにあらわされます。
```math
\begin{align}
 e_{k+1}  = (A- L_1 CA - L_2 C) e_k + (E - L_1 C E) d_k 
\end{align}
```
したがって、
1. $E - L_1 C E=0$となる$L_1$が選べる
2. $(A- L_1 CA - L_2 C)$が安定極に極配置できる

場合に、推定誤差系が安定、かつ外乱$d_k$に不感となります。つまり、外乱の影響を一切受けずに推定誤差が時間とともに0へ向かうようになります。

前者1.の条件は、${\rm rank}(CE)=r$と等価であることが示されています。記号の意味を再掲すると、出力の係数行列$C$と外乱の係数行列$E$をかけた階数は外乱の次数$r$と同じ、という条件です。
後者2.の条件は、$(C, \ A-E(CE)^{+}CA)$のペアが可検出であることと等価であることが示されています。なお、$(CE)^{+}=[(CE)^T(CE)]^{-1} (CE)^T$はムーア―ペンローズの疑似逆行列です。前者1.と被ってる気がしますが、$(CE)^T(CE)$が非正則な場合は当然この条件も満たされません。


<!-- 
### 未知入力オブザーバにおける外乱の推定式
もともとの離散時間系のダイナミクスの両辺に$C$をかけ、$C x\_{k+1}=y\_{k+1}$を用いると外乱$d_k$は下式であらわせます。
```math
\begin{align}
 d_{k} &= (CE)^{+} (y_{k+1} - CA x_k - CBu_k) 
\end{align}
```

これを踏まえ、前述で得られた状態の推定値$\hat{x}_k$を用いて、外乱の推定を下式で行うことを考えます。
```math
\begin{align}
 \hat{d}_{k+1} &= (CE)^{+} (y_{k+1} - CA\hat{x}_k - CBu_k) 
\end{align}
```

外乱の推定誤差を書き下すと$\hat{d}\_{k+1}-d_k＝(CE)^{+} CA (\hat{x}_k - x_k)$となります。
これでは、サンプリング周期が短くなったことで$(CE)^{+}$が大きくなる・・・などによって過渡的に大きな推定誤差が生まれることが想像されます。

そこで、実際には下式を使って外乱を推定します。

```math
\begin{align}
 \hat{d}_{k+1} &= (CE)^{+} (y_{k+1} - CA\hat{x}_k - CBu_k) + L_3 (C \hat{x}_k - y_k)
\end{align}
```
$L_3$はオブザーバゲインで、再び外乱の推定誤差を書き下すと$\hat{d}\_{k+1}-d_k＝(L_3 C - (CE)^{+} CA) (\hat{x}_k - x_k)$となります。
$L_3=(CE)^{+}(CAC^T)(CC^T)^{-1}$と選ぶことで、この誤差を最小化することが提案されています。
-->

### シングルレートUIOにおける外乱の推定式

外乱の推定式は、状態の推定値$\hat{x}_k$を用いて、下式で得られます。
(もともとの離散時間系のダイナミクスの両辺に$C$をかけ、$C x\_{k+1}=y\_{k+1}$を用いて導出できます)

```math
\begin{align}
 \hat{d}_{k} &= (CE)^{+} (y_{k+1} - CA\hat{x}_k - CBu_k) 
\end{align}
```

元文献では左辺を$d\_{k+1}$としていましたが、あえて上記の表現にしました。$y_{k+1}$よりも過去の外乱が推定できる、という様子が正しいと思われます。
$L_3 (C \hat{x}_k - y_k)$という項を追加して過渡特性を改善できるという情報もありましたが、いまいちシミュレーションで効果のほどがわからなかったので、この記事では使わないことにします。



まとめると、条件1.、条件2.を満たす
(もしくは等価な条件1. ${\rm rank}(CE)=r$となる、条件2.$(C, \ A-E(CE)^{+}CA)$が可検出となる)場合のシングルレートUIOの設計手順は下記のとおりです。
1. $L_1 = E(CE)^{+}$ で$L_1$を求める。
2. 先に求めた$L_1$から、$(A- L_1 CA - L_2 C)$が安定となるように$L_2$を極配置法などで求める。
3. 上述の$\hat{x}\_{k+1}$、$\hat{d}\_{k}$についての推定式を実装する。

### シングルレートUIOの構成できる条件
<!--
さらっと条件を書きましたが、意外とこの条件2つは絡み合っています。例えば後述のバネマスダンパ系のシミュレーションを例に考えてみましょう。
位置、速度の2次元系にて外乱として力が加わり、位置のみが観測できるケースを考えてみます。すると、$C = [1 \ 0]$、$E = [0 \ 1]^T$となり、$r=1$で ${\rm rank}([1 \ 0] [ 0 \ 1]^T)=0$なので未知入力オブザーバは構成できません。
感覚的には、外乱の一つ遅れた位相の物理量まで観測できないとだめって感じですかね。この例はわかりやすいです。

では、速度のみが観測できるケースは未知入力オブザーバが構成できるのかというと、それも違います。$C = [0 \ 1]$、$E = [0 \ 1]^T$となり、$r=1$で ${\rm rank}([0 \ 1] [ 0 \ 1]^T)=1$なので条件1.的にはOKなのですが、実は条件2.が満たされません。 
$m=r$、つまり観測できる次数と外乱の次数が等しい場合、かならず条件2.が満たされなくなります。
執筆者の力及ばずでこの点詳しく説明できませんが、後述のシミュレーションで簡易的に可観測性についてはチェックをしてみます。
-->

さらっと条件を書きましたが、意外とこの2条件はやっかいです。

条件1.は連続時間系ではなく離散時間系の$E$行列についての判定です。もともとの$E=(\int^T_0 e^{A\tau} d\tau) E_c$に立ち返って$E_c$と$C$のみで判定したいところですが、この点は深堀りせず先に進みたいと思います。

また、条件1.を満たしても条件2.は満たさない場合もあります。特に、観測の次数$m$と外乱の次数$r$が一致しているとき、必ず条件2.は満たされなくなります。この点も長くなりそうなので詳細は割愛して、結果だけ記載します。


観測の次数$m$と外乱の次数$r$の大小関係で分けてUIOが構成できる条件をまとめると
 $m<r$の場合 → 明らかにUIOは構成できない
 $m>r$の場合 → 上述の条件.1、条件.2を満たせばシングルレートUIOが構成できる
 $m=r$の場合 → 上述のシングルレートUIOは構成できない

となっています。これに対し、$m=r$へ適用範囲を広げたUIOが提案されています。次項で紹介します。

<!-- 
ふとここで脱線して${\rm rank}(CE)=r$について考えてみると、やや厳しい条件であることがわかります。
例えば後述するDCモータのシミュレーションで位置、速度の2次元系にて外乱として力が加わり、位置のみが観測できるケースを考えると、$C = [1 \ 0]$、$E = [0 \ 1]^T$となります。
このとき$r=1$で ${\rm rank}([1 \ 0] [ 0 \ 1]^T)=0$となるので、位置のみ観測できDCモータに外力がかかる、という条件では適用できません。
(速度のみ観測できる場合($C = [0 \ 1]$)、位置速度両方観測できる場合($C$が2次元単位行列)はどちらも適用できます。感覚的には、外乱の一つ遅れた位相の物理量まで観測できないとだめって感じですかね。）
-->

# より適用範囲を広げた未知入力オブザーバの構成

上述の$m=r$の条件下での適用を考えていきます。以下の議論は$m=r$を前提にします($m>r$なら前述のUIOを使えばいいだけなので)。

ここで、一度もともとの制御対象の連続時間系に立ち返り、サンプリング周期$T$内でもう一点観測データ$z_k$が得られる、という状況を考えてみます。
$0<i<1$で、$i=1/2$ならセンサーの出力だけ2倍の周波数でサンプリングしている、という状況です。

![未知入力オブザーバ.drawio.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/701220/7a56b3e1-b74b-93d3-2c7f-f3b3ced35b83.png)



このとき、もともとの離散時間方程式に補助出力$z_k$が追加された下記の離散時間系を作ることができます。
いわゆるマルチサンプリングレートの2-delay出力離散時間系です。
```math
\begin{align}
 x_{k+1} &= A x_k + B u_k + E d_k \\
\left[
\begin{array}{cc}
y_k \\
z_k
\end{array}
\right] &= \left[
\begin{array}{cc}
C \\
C \tilde{A}
\end{array}
\right] x_k +
\left[
\begin{array}{cc}
0 \\
C \tilde{B}
\end{array}
\right] u_k + 
\left[
\begin{array}{cc}
0 \\
C \tilde{E}
\end{array}
\right] d_k 
\end{align}
```

ここで、$\tilde{A}=e^{A_c i T}$、$\tilde{B}=(\int^{iT}_0 e^{A_c \tau} d \tau) B_c$、$\tilde{E}=(\int^{iT}_0 e^{A_c \tau} d \tau) E_c$です。

この系に対し、先ほどとよく似た構成のオブザーバの適用を考えます。
これも執筆者の造語ですが、ダブルレートUIOと呼ぶことにします。
```math
\begin{align}
 \hat{x}_{k+1} &= (A- L_1 C \tilde{A} - L_2 C) \hat{x}_k + (B-L_1 C \tilde{B}) u_k + L_1 z_k + L_2 y_k \\
\end{align}
```

このオブザーバに対し、再び推定誤差$e_{k+1} = x\_{k+1} - \hat{x}\_{k+1}$と置くと、推定誤差ダイナミクスは下記のようにあらわされます。
```math
\begin{align}
 e_{k+1}  &= x_{k+1} - \hat{x}_{k+1} \\
 &= (A x_k + B u_k + E d_k) - ((A- L_1 C \tilde{A} - L_2 C) \hat{x}_k + (B-L_1 C \tilde{B}) u_k + L_1 z_k + L_2 y_k)\\
 &= (A x_k + B u_k + E d_k) - ((A- L_1 C \tilde{A} - L_2 C) \hat{x}_k + (B-L_1 C \tilde{B}) u_k + L_1 (C \tilde{A} 
 x_k + C \tilde{B} u_k +  C \tilde{E} d_k) + L_2 C x_k )
\end{align}
```
これを整理すると、$u_k$に関する項が早々に消え、先ほどとよく似た下式が得られます。
```math
\begin{align}
 e_{k+1}  &= (A- L_1 C \tilde{A} - L_2 C) e_k  + (E-L_1 C \tilde{E})  d_k
\end{align}
```

したがって、
1. $E - L_1 C \tilde{E}=0$となる$L_1$が選べる
2. $(A- L_1 C\tilde{A} - L_2 C)$が安定極に極配置できる


場合、外乱$d_k$に不感となるオブザーバが構成できます。
先ほどと違うのは、$m=r$前提なので疑似逆行列を持ち出さなくても良い点です。
したがって、上の条件1.は$C \tilde{E}$が正則かどうかと等価で、正則な場合には$L_1 = E(C \tilde{E})^{-1}$と選べばよいです。
条件2.は$L_1 = E(C \tilde{E})^{-1}$を代入の上、$(C, A- E(C \tilde{E})^{-1}C\tilde{A})$のペアが可検出であればよいです。



ここでまだしっくりこないのは、各々の条件に記号~、つまり$i$が含まれているからです。$i$を含まない条件に変えられれば、もともとの周期$T$内のどの位置でサンプリングさせるか、を気にせずUIOが構成できるか判定できます。

この点についても心配なく、それぞれに対応した条件として
1. ほとんどすべての$i$について、$CE$が正則なら$C \tilde{E}$も正則となる
2. もともとの連続時間系の$(C,A_c,E_c)$の組み合わせに関して(つまり$d$から$y$の経路に)原点の零点が存在しないとき、ほとんどすべての$i$にて$(C, A- E(C \tilde{E})^{-1}C\tilde{A})$が可検出となる

ことが示されています。やはり条件2がちょっとわかりにくいですが、これも後述の数値例で確認することとします。

### ダブルレートUIOにおける外乱の推定式
$m=r$ではないときと同様、外乱の推定式は下式となります。やはり、$z_k$より過去に位置する$d_k$が推定される式です。
```math
\begin{align}
 \hat{d}_{k} &= (C\tilde{E})^{-1} (z_k - C \tilde{A} \hat{x}_k - C \tilde{B} u_k) 
\end{align}
```

まとめると、$m=r$のときも、条件1.、条件2.を満たす場合にはダブルレートUIOであれば構成でき、設計手順は下記のとおりです。
1. $i$、つまりサンプリング周期$T$間でもう一つ測定値をとるタイミングを設定する。
2. $L_1 = E(C\tilde{E})^{-1}$ で$L_1$を求める。
3. 先に求めた$L_1$より、$(A- L_1 C\tilde{A} - L_2 C)$が安定となるように$L_2$を極配置法などで求める。
4. 上述の$\hat{x}\_{k+1}$、$\hat{d}\_{k}$についての推定式を実装する。


# 外乱オブザーバの構成

外乱オブザーバについては解説も豊富なので、結果のみ示します。状態変数などの記号は未知入力オブザーバと被りますがご了承ください。

離散時間系の外乱オブザーバは$d\_{k+1}=d_k$のステップ状の外乱を想定する場合、拡大系を構成して下式で得られます。
```math
\begin{align}
 \left(
\begin{array}{c}
\hat{x}_{k+1} \\
\hat{d}_{k+1}
\end{array}
\right) &= \left(
\begin{array}{cc}
A & E \\
0 & 1
\end{array}
\right)  \left(
\begin{array}{c}
\hat{x}_k \\
\hat{d}_k
\end{array}
\right) + \left(
\begin{array}{cc}
B \\
0
\end{array}
\right) u_k + L\left(C \left(
\begin{array}{c}
\hat{x}_k \\
\hat{d}_k
\end{array}
\right)  - y_k \right)
\end{align}
```
この状態方程式はもともとより次数が1大きく、オブザーバゲイン$L$もそれに合ったサイズの行列であることにご注意ください。

# シミュレーションによる比較

MATLAB/Simulinkコードは下記に置きます。
Control System Toolboxを使用しています。

https://github.com/SKMAAX/Control/tree/main/MATLAB/UIO
```MATLAB:DOB_vs_UIO.m
%% 物理パラメータの設定
J = 1;
c_th = 0.1;
k = 1;

%% 連続時間＆離散時間の状態空間モデルの設定
Ac = [0 1;-c_th/J -k/J ];
Bc = [0;1/J];
Ec = [0;1/J]; % 外乱の係数行列は入力と同じ。マッチング条件を満たす。
C = [1 0];
% C = [0 1]; %% 速度観測ではこのUIOは構成できない
Dc = 0;
T = 0.001;
[A, B, ~, D] = c2dm(Ac,Bc, C, Dc, T);
[~, E, ~, ~] = c2dm(Ac,Ec, C, Dc, T);

%% シングルレートのUIO構成条件のチェック
rank([eye(2); eye(2)*(A-E*inv((eye(2)*E)'*(eye(2)*E))*(eye(2)*E)'*eye(2)*A)]) %ランク2なのでOK
rank([[0 1]; [0 1]*(A-E*inv(([0 1]*E)'*([0 1]*E))*([0 1]*E)'* [0 1]*A)]) %ランク1なのでNG
rank([[1 0]; [1 0]*(A-E*inv(([1 0]*E)'*([1 0]*E))*([1 0]*E)'* [1 0]*A)]) %ランク1なのでNG

%% 未知入力オブザーバの設計
i = 0.5; %(0<i<1)の範囲で設定する
[A_til, B_til, C_til, D_til] = c2dm(Ac,Bc, C, Dc, T*i);
[~, E_til, ~, ~] = c2dm(Ac,Ec, C, Dc, T*i);

invCE_til = inv(C*E_til);

L1 = E * invCE_til;
pole_uio = [0.9 0.8];
L2 = place((A - L1*C*A_til)', C', pole_uio)';
%L3 = invCE_til * (C*A_til*C') * inv(C*C');
A_hat = A - L1*C*A_til - L2*C;
B_hat = B - L1*C*B_til;

%% ダブルレートのUIO構成条件のチェック
rank([[0 1]; [0 1]*(A-E*inv([0 1]*E_til)*[0 1]*A_til)]) %ランク1なのでNG
rank([[1 0]; [1 0]*(A-E*inv([1 0]*E_til)*[1 0]*A_til)]) %ランク2なのでOK

%% 外乱オブザーバ 
A_dob = [A E;0 0 1]; B_dob = [B;0]; C_dob = [C 0]; D_dob=0; % 離散時間の拡大系
pole_dob = [0.9 0.8 0.7];
L_dob = place(A_dob', C_dob', pole_dob)';

%% シミュレーション実行
out = sim('sim_DOB_vs_UIO');
```


簡単なシミュレーションとして、下記の連続時間系への適用を考えます。$x$は位置・速度からなる二次元の状態ベクトルで、例えばバネマスダンパ系を想定したような状態方程式です。
入力$u$、外力$d$ですが、今回入力に値が入ることはありません(今後のことを考え、一応残してあるだけです)。
観測出力は後から考えます。①位置・速度両方観測できる、②速度のみ観測できる、③位置のみ観測できる、のそれぞれを考えてみます。①は$m>r$の例、②と③は$m=r=1$となり、ダブルレートUIOを適用することとなります。この記事ではシングルレートUIOについては割愛して、ダブルレートUIOの挙動のみ検証してみます。

```math
\begin{align}
\dot{x}(t) &= \left(
\begin{array}{cc}
0 & 1 \\
-0.1 & -1
\end{array}
\right) x(t) + \left(
\begin{array}{cc}
0 \\
1
\end{array}
\right) u(t) + \left(
\begin{array}{cc}
0 \\
1
\end{array}
\right) d(t) 
\end{align}
```

この連続時間系を、サンプリング周期$T=0.001$[s]でゼロ次ホールド離散化します。
実際にはMATLABのc2dm関数で行ったので、コード内を確認ください。

<!--
実際にはMATLABのc2dm関数で行ったので、数値結果のみ載せます。
```math
\begin{align}
x_{k+1} &= \left(
\begin{array}{cc}
   1.0000\rm{e}+00 & 9.9950\rm{e}-04 \\
  -9.9950\rm{e}-05 & 9.9900\rm{e}-01
\end{array}
\right) x_k + \left(
\begin{array}{cc}
   4.9983\rm{e}-07 \\
   9.9950\rm{e}-04
\end{array}
\right) u_k + \left(
\begin{array}{cc}  
   \rm{4.9983e-07} \\
   \rm{9.9950e-04}
\end{array}
\right) d_k
\end{align}
```
-->

同様に、補助出力$z_k$の観測タイミングを$i=0.5$として、補助出力に関わる各行列を求めます。
こちらもMATLABのc2dm関数を使用しました。
$E$と$\tilde{E}$のみ、条件チェックに関わるところなので数値を載せます。下記の結果となりました。
```math
\begin{align}
\tilde{E}= \left(
\begin{array}{cc}
   0.0005\rm{e}-03 \\
   0.9995\rm{e}-03
\end{array}
\right), \ \tilde{E}= \left(
\begin{array}{cc}
   1.2498\rm{e}-07 \\
   4.9988\rm{e}-04
\end{array}
\right)
\end{align}
```

### 位置・速度両方観測の場合
位置・速度両方観測できる場合は、$C$が2次元の単位行列となります。
シングルレートUIOの条件1.は
```math
\begin{align}
{\rm rank } (CE) = {\rm rank } \left(
\left(
\begin{array}{cc}
 1 & 0 \\
 0 & 1 
\end{array}
\right)
\left(
\begin{array}{cc}
   0.0005\rm{e}-03 \\
   0.9995\rm{e}-03
\end{array}
\right)
\right) = 1
\end{align}
```
より、${\rm rank} (CE) = r$なので満たされます。

シングルレートUIOの条件2. $(C, \ A-E(CE)^{+}CA)$が可検出となる・・・のチェックは行列がごちゃごちゃしてくるので過程は省略します。気になる方はコードを確認ください。結局、可観測性行列を作ってみるとランク2となるので、この条件も満たされます。

したがって、位置・速度両方観測できる場合にはシングルレートUIOが構成できます。


### 速度のみ観測の場合
速度観測の場合は、$C=[0 \ 1]$となります。

シングルレートUIOの条件1.は $CE$が非ゼロの定数になるので、${\rm rank} (CE) = 1$となって満たされます。
シングルレートUIOの条件2.は $(C, \ A-E(CE)^{+}CA)$の可観測性行列を作ってみるとランク1となってしまい、満たされません(本来は可検出性を見るべきですが省略)。 やはり、$m=r$の条件下ではシングルレートUIOは構成できないようです。

続けて、ダブルレートUIOの方の条件を見てみます。
条件1.は、$C \tilde{E}$は非ゼロの定数となり明らかに正則なので満たされます。

条件2.の方を確認していきます。
$d$から$y$までの伝達関数は
```math
\begin{align}
G(s) &= C(sI-A_c)^{-1} E_c \\
&= C \ \frac{{\rm adj}(sI-A_c)}{{\rm det}(sI-A_c)} \ E_c
\end{align}
```
で得られます。実際に数値を入れていくと
```math
\begin{align}
G(s) &= (0 \ 1) \ \frac{{\rm adj}
\left(
\begin{array}{cc}
s & -1\\
0.1 & s+1
\end{array}
\right)
}{{\rm det} 
\left(
\begin{array}{cc}
s & -1\\
0.1 & s+1
\end{array}
\right)
} \ \left(
\begin{array}{cc}
0 \\
1
\end{array}
\right) \\
&= (0 \ 1) \ \frac{
\left(
\begin{array}{cc}
s+1 & 1\\
-0.1 & s
\end{array}
\right)
}{
s^2 + s + 0.1
} \ \left(
\begin{array}{cc}
0 \\
1
\end{array}
\right) \\
&= \frac{s}{
s^2 + s + 0.1
} 
\end{align}
```
最後の行で、$C$と$E_c$によって余因子行列の右下が抽出されました。
したがって、$d$から$y$の特性に原点ゼロが含まれており、速度のみ観測した場合には条件2.を満たさないことがわかります。
試しに$(C, \ A- E(C\tilde{E})^{-1} C\tilde{A})$の可観測性行列を作ってみると、やはりランク1となっています。

結局、速度のみ観測できる場合にはUIOが構成できないことがわかりました。

### 位置のみ観測の場合

位置観測の場合は、$C=[1 \ 0]$となります。

シングルレートUIOの条件1.は $CE$が非ゼロの定数になるので、${\rm rank} (CE) = 1$となって満たされます。
シングルレートUIOの条件2.は $(C, \ A-E(CE)^{+}CA)$の可観測性行列を作ってみるとランク1となってしまい、満たされません(本来は可検出性(以下略))。 やはり、$m=r$の条件下ではシングルレートUIOは構成できないようです。

続けて、ダブルレートUIOの方の条件を見てみます。
条件1.は、$C \tilde{E}$が非ゼロの定数となり明らかに正則なので満たされます。

条件2.は、速度観測の場合での計算過程を流用して下記の結果となります。
```math
\begin{align}
G(s) 
&= (1 \ 0) \ \frac{
\left(
\begin{array}{cc}
s+1 & 1\\
-0.1 & s
\end{array}
\right)
}{
s^2 + s + 0.1
} \ \left(
\begin{array}{cc}
0 \\
1
\end{array}
\right) \\
&= \frac{1}{
s^2 + s + 0.1
} 
\end{align}
```
余因子行列の右上が抽出されました。
こちらは、零点が含まれておらず、条件2.が満たされます。

したがって、位置観測の場合なら、ダブルレートUIOが構成できることがわかりました。
なぜ速度観測ではだめなのか？をもっと掘り下げたいところですが、執筆者の力不足もあり、何か補足できそうなら後日追記します。

## 未知入力オブザーバと外乱オブザーバの構成比較
以上で、位置観測であればダブルレートUIOが構成できるのがわかったので、実際に外乱オブザーバと比較してみます。
(外乱オブザーバの方も、位置観測にて拡大系が可観測であることが言えますが省略します)

Simulinkでの接続構成は下図のとおりです。入力はなく、外乱のみを印可します。
UIO、外乱オブザーバともに$y(t)$のサンプル値のみから外乱を推定する構成です。

外乱オブザーバは特に補足する点はなく、$T=0.001$[s]の各ステップで離散時間状態方程式を更新していく形です。
UIOについても、おおよそ列挙してきた式の通りの接続です。並べてみると、UIOでは2出力をサンプリングするオブザーバである、という点が強調されていると思います。

処理周期は外乱オブザーバとUIOで変わらず$T=0.001$で動いていることに注意ください。
赤丸部分が若干ハマったポイントでした。
この部分で、$y_k$はサンプリング後に$T$秒遅らせて、$z_k$は$(1-i)T$秒遅らせた後にサンプリングしています。

![未知入力オブザーバのSimulink接続図.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/701220/7b5064a6-e9d0-2c3d-c438-2711ca6929d5.png)

若干わかりにくいので、下記でもう少し補足します。

### 未知入力オブザーバの演算のタイミング 
演算のタイミングは下記のようなイメージになっています。$z_k$が得られた瞬間から$d_k$の推定はできますが、このモデルでは$k+1$の時点で$d_k$と$\hat{x}_k$を推定しています。

![未知入力オブザーバ  演算タイミング.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/701220/31d9522e-90be-bf2c-1204-8e812c0c751b.png)

マイコン等で実装する場合には、「演算時間」と書いた矢印の区間内で計算をすることになると思われます。したがって、考え方的には$i=0.5$とするのが最もわかりやすいですが、演算時間の余裕的には$i$は小さい方が有利です。ただし、$i$が小さいほど$(C\tilde{E})^{-1}$が大きくなります。つまり全体的にハイゲインになり、ノイズ等々の影響を受けやすくなります。
この点は、$i$を小さくするほど微小区間内で$d$が$z$に現れてくる様子をシビアにみる必要がありそう、という様子からも直感的にわかるかと思います。


## シミュレーション結果

だいぶ長くなってしまいましたが、ようやく結果を見てみます。
UIOのゲイン$L_2$は$[0.9 \ 0.8]$へ、外乱オブザーバのオブザーバゲインは$[0.9 \ 0.8 \ 0.7]$へ極配置することで得ました。

### ステップ状外乱の推定の様子
ステップ状外乱を印可し、それを両手法で推定した様子が下記です。
左上が位置、右上が速度で、左下が外乱推定の様子を表しています、それぞれの下部のグラフは、上部の誤差をプロットしたものです。外乱オブザーバには推定遅れがあるのに対し、UIOは$T=0.001$後に完全に外乱と同一の値となっています。

位置、速度の推定の様子がUIOだと悪いように見えますが、これは推定の値が1ステップ遅れているからです。そのため、元信号の変化量が大きければ大きいほど誤差が大きくなってしまっています。

1ステップ遅れていることを前提に時間をずらすと、外乱オブザーバと比較しても圧倒的に小さな誤差に収まります。時間をずらして比較しても良かったのですが、1ステップ遅れた値しかわからない、というのもデメリットと思ったためあえてこのように比較しました。以降のシミュレーション結果でも同様です。

![未知入力オブザーバ.drawio (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/701220/28995c71-8c31-3807-ac49-f6a632f5847f.png)


### 乱数状外乱の推定の様子
－1から1の範囲で変動する一様乱数を外乱として与えた結果も比較しました。
外乱オブザーバには推定遅れであまり外乱の特徴を捉えられていませんが、UIOは1ステップ後に完全に外乱と同一の値となっています。

ここでもUIOの波形は1ステップ遅れていることの補正はしていないため、信号の変化量が大きいところでは見かけ上誤差が大きく出ています。
信号全体を1ステップずらすと、どの波形も外乱オブザーバと比較して圧倒的に小さな誤差に収まります。
![未知入力オブザーバ.drawio (2).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/701220/fd50efc7-71c6-327e-236c-3b3761e9b246.png)

### ステップ状外乱の推定の様子(センサノイズあり)
先に見せたステップ状外乱の推定にて、センサ出力に-0.00001～0.00001の一様乱数を加えました。
結果、外乱オブザーバはノイズがあってもある程度の推定誤差に収まっているのに対し、UIOはほぼまともに推定できていない様子がわかります。右下のグラフは、左下の拡大図です。

シミュレーション内での位置の変化はせいぜい0.1以内で、それの1/10000程度のノイズを入れただけでUIOは推定が行えない様子が確認されました。

![未知入力オブザーバ.drawio (3).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/701220/b7e8ccc5-1a4d-003e-d5e4-e2b3d24243f0.png)

## 結果の解釈
上述のシミュレーション結果から、UIOはたしかに次ステップ(最短で$iT$後も可能)で外乱の正確な値を推定し、結果は外乱の性質によらないことがわかりました。

外乱オブザーバの方は、初めにステップ状外乱であることを仮定して、外乱モデルを状態方程式の中に組み込んでいる特徴があります。外乱オブザーバはステップ状の他にもランプ状、正弦波状の外乱等も推定できますが、いずれにしてもあらかじめそれらを観測モデルの中に組み込む必要があり、そこで想定した信号以外は推定できなかったり、時間がかかったりします。
これと比較したUIOは、まったく同じオブザーバで様々な外乱を推定することができます。

しかし、UIOには実用上の弱点も多いと感じました。観測値にわずかでもノイズが乗る場合、シミュレーションの限りではまともに推定値が得られない様子でした。今回試してはいませんが、逆モデルを構築して$d$を逆算するような形になっている都合上、モデル化誤差の影響も多分に受けそうな印象です。

# おわりに

この記事では、未知入力オブザーバ(UIO)について紹介し、外乱オブザーバと外乱推定の様子を比較しました。

どんな外乱でも同じオブザーバでピッタリ推定できてしまう、というUIOの挙動は非常に魅力的でした。しかし、マルチレートでやや複雑な推定則を実装しないといけない点、ノイズ等々の影響を受けやすい点から、外乱を推定する用途ではまだまだ産業界では外乱オブザーバが活躍し続けるのかなと感じました。

またこれも執筆者の私見ですが、一般にAD変換でセンサ値を取得する方がDA変換でアナログ値を外部出力するより時間がかかる、という性質があります。制御系の周期はセンサ値がどの周期でとってこれるか、で決まると言っても過言ではないと思っています。これらを考慮して制御入力(アナログ出力)の方のレートを上げるマルチレート制御もあるぐらいなので、その点踏まえてUIOはやや実装ハードルが高い印象でした。

これらを克服している、実際にこんな分野で応用されている、というUIOの事例があればぜひご紹介いただきたいです。また、間違い等ありましたらコメントでお知らせください。

ご覧いただきありがとうございました。



[^1]: 「Design of general structured observers for linear systems with unknown inputs」, Shao-Kung Chang, Wen-Tong You, Pau-Lo Hsu, Journal of the Franklin Institute, 1997, https://www.sciencedirect.com/science/article/abs/pii/S0016003296000774

[^2]: 「外乱オブザーバ」、島田明 著、コロナ社、 2021、https://www.coronasha.co.jp/np/isbn/9784339033823/