線形制御系設計

    線形システムの表現
    線形システムの特性解析
    制御器設計
    観測器設計

    モーションコントロール

    加速度制御
    剛体の姿勢制御

    モータドライブ

    駆動制御法
    磁極位置推定法

    システム同定

    ホワイトボックスモデル

姿勢推定器の設計

システムモデル

姿勢を表すクオタニオン q\bm{q} に関する支配方程式は、剛体系の角速度を ω\bm{\omega} として次のように記述される。
q˙=12[r1r2r3r0r3r2r3r0r1r2r1r0]ωJq(q)ω, \begin{align} \dot{\bm{q}} &= \frac{1}{2} \begin{bmatrix} -r_1 & -r_2 & -r_3 \\ r_0 & -r_3 & r_2\\ r_3 & r_0 & -r_1\\ -r_2 & r_1 & r_0 \end{bmatrix} \bm{\omega} \\ &\equiv\bm{J}_{q}(\bm{q})\bm{\omega}, \end{align}
ただし、Jq(q)\bm{J}_{q}(\bm{q}) はヤコビ行列を表す。

離散モデル

計算機で状態の予測と修正を行うために、離散時間モデルを導出する。計算機のサンプリング時間を TsT_{\rm s} として、クオタニオンの 11 ステップ先予測モデルを次のように設計する。
qk+1=qk+TsJq(qk)ωk\begin{align} \bm{q}_{k+1} &= \bm{q}_{k} + T_{\rm s}\bm{J}_{q}(\bm{q}_{k})\bm{\omega}_{k} \end{align}
ただし、下添字 kk はサンプリング点を表すものとする。このモデルでは矩形積分による予測を行っており、積分誤差はプロセス雑音と見做すものとする。本誤差はサンプリング時間が小さくなるにつれて減少し、観測器の誤差修正によってこれによる影響を抑圧するものとする。

観測モデル

観測器の修正は、基準となる状態量と予測された状態量を比較し、その誤差が減少するように行われる。回転するシステムでは機体座標系が基準座標系に対して回転するため、状態量をいずれかの座標系に射影して観測量の比較を行う必要がある。機体の姿勢を表すクオタニオン q\bm{q} が与えられた場合、機体座標系で xbR3\bm{x}^{\rm b} \in \mathbb{R}^{3} と表されるベクトルは基準座標系で次のように表される。
xr=Rbr(q)xb\begin{align} \bm{x}^{\rm r} &= \bm{R}_{\rm br}(\bm{q}) \bm{x}^{\rm b} \end{align}
ただし、Rbr\bm{R}_{\rm br} は以下に示す回転行列とする。
Rbr(q)[12(r22+r32)2(r1r2r0r3)2(r1r3+r0r2)2(r1r2+r0r3)12(r22+r32)2(r2r3r0r1)2(r1r3r0r2)2(r2r3+r0r1)12(r12+r22)]\begin{align} \bm{R}_{\rm br}(\bm{q}) \equiv \begin{bmatrix} 1 - 2\left(r_2^2 + r_3^2\right) & 2\left(r_1r_2 - r_0r_3\right) & 2\left(r_1r_3 + r_0r_2\right)\\ 2\left(r_1r_2 + r_0r_3\right) & 1 - 2\left(r_2^2 + r_3^2\right) & 2\left(r_2r_3 - r_0r_1\right)\\ 2\left(r_1r_3 - r_0r_2\right) & 2\left(r_2r_3 + r_0r_1\right) & 1 - 2\left(r_1^2 + r_2^2\right) \end{bmatrix} \end{align}
また、回転行列の性質より、逆変換は次のように表現される。
xb=RbrT(q)xr\begin{align} \bm{x}^{\rm b} &= \bm{R}_{\rm br}^{\mathrm T}(\bm{q}) \bm{x}^{\rm r} \end{align}

推定量の修正

非線形システムの予測誤差の二乗和を最小化することを目的として、ガウス・ニュートン法を使用する。ガウス・ニュートン法は非線形写像のパラメータ推定に使用することができ、元および像が x,yRm\bm{x}, \bm{y} \in \mathbb{R}^{m} として与えられた時、以下で与えられる予測残差 ek\bm{e}_{k}L2L_2 ノルムを最小化する写像 ff を反復計算によって求める。
ek=yky^ky^k=f(x,p)\begin{align} \bm{e}_{k} &= \bm{y}_{k} - \hat{\bm{y}}_{k} \\ \hat{\bm{y}}_{k} &= f(\bm{x}, \bm{p}) \end{align}
ただし、kk は反復試行数、pRn (nm)\bm{p} \in \mathbb{R}^{n}\ (n\le m) は写像のパラメータであり、写像 ff を求めることはこのパラメータ p\bm{p} を求めることと等価となる。ここで、関数 ff のヤコビ行列 Jf\bm{J}_{f} を用いて次のようにパラメータ更新を行うことで、予測残差の L2L_2 ノルムが減少する方向に写像が修正される。
JffpRm×npk+1=pk+Jf+(pk)ek\begin{align} \bm{J}_{f} &\equiv \frac{\partial f}{\partial \bm{p}} \in \mathbb{R}^{m\times n}\\ \bm{p}_{k+1} &= \bm{p}_{k} + \bm{J}_{f}^{+}(\bm{p}_{k})\bm{e}_{k} \end{align}
ただし、Jf+\bm{J}_{f}^{+}Jf\bm{J}_{f} の一般化逆行列であり、nmn\le m の条件下で次のように定義される。
Jf+=(JfTJf)1JfT\begin{align} \bm{J}_{f}^{+} = \left(\bm{J}_{f}^{\mathrm T}\bm{J}_{f}\right)^{-1}\bm{J}_{f}^{\mathrm T} \end{align}
この一般化逆行列 Jf+\bm{J}_{f}^{+}Jf\bm{J}_{f} の左逆元となることに留意する。この方法は、初期推定値が真値から離れている場合や写像 ff の非線形性が大きい場合にはパラメータ推定値が発散する恐れがある。そこで、縮小因子 γ\bm{\gamma} を導入して更新方程式を次のように修正する。
pk+1=pk+γJf+ek\begin{align} \bm{p}_{k+1} &= \bm{p}_{k} + \bm{\gamma}\bm{J}_{f}^{+}\bm{e}_{k} \end{align}

退化と不可観測性

ガウス・ニュートン法を使用する条件として、参照データ y\bm{y} の次元 mm が パラメータ p\bm{p} の次元 nn 以上であることとしている。これは、写像によって情報が退化した場合にパラメータ修正則が定まらないためと考えられる。

便宜上、写像 ff のヤコビ行列 Jf\bm{J}_{f} を用いて、y^k\hat{\bm{y}}_{k} を近似的に次のように表現する。
y~k=Jfpk\begin{align} \tilde{\bm{y}}_{k} &= \bm{J}_{f}\bm{p}_{k} \end{align}
このとき、パラメータ修正項を Δ\bm{\Delta} とすると、これは近似的に次のように表される。
ΔkJf+(yky~k)=Jf+ykJf+Jfpk\begin{align} \bm{\Delta}_{k} &\approx \bm{J}_{f}^{+}(\bm{y}_{k} - \tilde{\bm{y}}_{k}) \\ &= \bm{J}_{f}^{+}\bm{y}_{k} - \bm{J}_{f}^{+}\bm{J}_{f}\bm{p}_{k} \end{align}
ここで、Jf+\bm{J}_{f}^{+}Jf\bm{J}_{f} の左逆元となる場合には、右辺第二項は pk\bm{p}_{k} に復元される。一方で、n>mn > m となる場合には、一般可逆行列は左逆元とならないために復元されない。これは Jf\bm{J}_{f} の零空間 (退化空間) に起因する退化が原因と考えられる。パラメータ pk\bm{p}_{k} が 零空間 Ker(Jf)\rm Ker({\it \bm{J}_{f}}) の成分 pnull\bm{p}_{\rm null} を持ち、次のように記述される場合について考える。
pk=pm+pnullpmKer(Jf)\begin{align} \bm{p}_{k} &= \bm{p}_{\rm m} + \bm{p}_{\rm null}\\ \bm{p}_{\rm m} &\notin {\rm Ker}(\bm{J}_{f}) \end{align}
ここで、零空間成分 pnull\bm{p}_{\rm null} は写像 ff によって退化し、像の観測情報から元を逆推定することができない。実際に、パラメータ修正項 Δ\bm{\Delta} の近似値は次のように記述される。
ΔkJf+ykpm\begin{align} \bm{\Delta}_{k} &\approx \bm{J}_{f}^{+}\bm{y}_{k} - \bm{p}_{\rm m} \end{align}
このように、零空間の成分が修正項に反映されず、零空間成分が不定となる。したがって、n>mn > m の条件下ではパラメータ決定に冗長自由度が発生し、真値を推定することが困難となる。

オブザーバ設計

上記で設計した状態更新モデル、観測モデル、および推定量修正則に従って姿勢オブザーバを設計する。ここではシステムが慣性計測ユニットおよび地磁気センサを搭載することを想定し、機体角速度 ω\bm{\omega}、機体加速度 a\bm{a}、および地磁気 m\bm{m} が測定されるものとする。また。地磁気センサの出力は適切な写像によって単位球上に移されているものとする。

予測ステップ

システムの状態予測に当たっては以下の離散モデルを使用する。
qk+1=qk+TsJq(qk)ωk\begin{align} \bm{q}_{k+1} &= \bm{q}_{k} + T_{\rm s}\bm{J}_{q}(\bm{q}_{k})\bm{\omega}_{k} \end{align}
ここで、ωk\bm{\omega}_{k} は角速度の真値であるが、予測ステップではセンサ計測値からこの値を推定する必要がある。角速度センサはオフセットを持つことがあり、その値が時間の経過に伴い変動することから、角速度の推定値 ω^\hat{\bm{\omega}} を測定値 ωs\bm{\omega}^{\rm s} および バイアスモデル ωb\bm{\omega}^{\rm b} を用いて次のように表す。
ω^=ωs+ωb\begin{align} \hat{\bm{\omega}} &= \bm{\omega}^{\rm s} + \bm{\omega}^{\rm b} \end{align}
また、バイアスの時間変動が十分に小さいものとし、動的モデルを次のように設定する。
ddtωb=0 \begin{align} \frac{d}{dt}\bm{\omega}^{\rm b} = \bm{0} \end{align}
以上より、システムの状態 x\bm{x} を次のように設定し、予測ステップを設計する。
x=[qωb]R7x~k+1=[I4TsJq(q^k)00]x^k+[TsJq(q^k)0]ωks\begin{align} \bm{x} &= \begin{bmatrix} \bm{q} \\ \bm{\omega}_{\rm b} \end{bmatrix} \in \mathbb{R}^7\\ \tilde{\bm{x}}_{k+1} &= \begin{bmatrix} \bm{I}_{4} & T_{\rm s}\bm{J}_{q}(\hat{\bm{q}}_{k}) \\ \bm{0} & \bm{0} \end{bmatrix} \hat{\bm{x}}_{k} + \begin{bmatrix} T_{\rm s}\bm{J}_{q}(\hat{\bm{q}}_{k}) \\ \bm{0} \end{bmatrix} \bm{\omega}^{\rm s}_{k} \end{align}
だたし、上添字付き変数 X^,X~\hat{X}, \tilde{X} はそれぞれ XX の推定値および予測値を表すものとする。

修正ステップ

観測情報から得られる修正勾配を用いて、状態を修正する。観測系における基準情報は、基準座標系において計測される重力加速度 agr\bm{a}^{\rm r}_{\rm g} と地磁気 mer\bm{m}^{\rm r}_{\rm e}とする。
agr=[00g]Tmer=[100]T\begin{align} \bm{a}^{\rm r}_{\rm g} &= \begin{bmatrix} 0 & 0 & g \end{bmatrix}^{\mathrm T} \\ \bm{m}^{\rm r}_{\rm e} &= \begin{bmatrix} 1 & 0 & 0 \end{bmatrix}^{\mathrm T} \end{align}
ただし、gg は重力加速度を表す。これらを纏めて基準観測量 y\bm{y} とする。
y=[agrmer]R6\begin{align} \bm{y} = \begin{bmatrix} \bm{a}^{\rm r}_{\rm g} \\ \bm{m}^{\rm r}_{\rm e} \end{bmatrix} \in \mathbb{R}^{6} \end{align}
機体の姿勢を表すクオタニオン が q\bm{q} と表される場合、機体座標系で計測された重力加速度 および地磁気を agb,meb\bm{a}^{\rm b}_{\rm g}, \bm{m}^{\rm b}_{\rm e} とすると、基準座標系の重力加速度 agr\bm{a}^{\rm r}_{\rm g} および地磁気 mer\bm{m}^{\rm r}_{\rm e} は次のように推定される。
a^gr=Rbr(x^k)agbm^er=Rbr(x^k)meb\begin{align} \hat{\bm{a}}^{\rm r}_{\rm g} &= \bm{R}_{\rm br}(\hat{\bm{x}}_{k}) \bm{a}^{\rm b}_{\rm g}\\ \hat{\bm{m}}^{\rm r}_{\rm e} &= \bm{R}_{\rm br}(\hat{\bm{x}}_{k}) \bm{m}^{\rm b}_{\rm e} \end{align}
ただし、Rrb\bm{R}_{\rm rb} は以下に示す回転行列とする。
Rrb(x)[12(r22+r32)2(r1r2+r0r3)2(r1r3r0r2)2(r1r2r0r3)12(r22+r32)2(r2r3+r0r1)2(r1r3+r0r2)2(r2r3r0r1)12(r12+r22)]\begin{align} \bm{R}_{\rm rb}(\bm{x}) \equiv \begin{bmatrix} 1 - 2\left(r_2^2 + r_3^2\right) & 2\left(r_1r_2 + r_0r_3\right) & 2\left(r_1r_3 - r_0r_2\right)\\ 2\left(r_1r_2 - r_0r_3\right) & 1 - 2\left(r_2^2 + r_3^2\right) & 2\left(r_2r_3 + r_0r_1\right)\\ 2\left(r_1r_3 + r_0r_2\right) & 2\left(r_2r_3 - r_0r_1\right) & 1 - 2\left(r_1^2 + r_2^2\right) \end{bmatrix} \end{align}
以上より、予測観測量 y^k\hat{\bm{y}}_{k} は次のように記述される。
y^k=[a^grm^er]=[2(r^1r^3r^0r^2)g2(r^2r^3+r^0r^1)g(12(r^12+r^22))g12(r^22+r^32)2(r^1r^2r^0r^3)2(r^1r^3+r^0r^2)]\begin{align} \hat{\bm{y}}_{k} = \begin{bmatrix} \hat{\bm{a}}^{\rm r}_{\rm g} \\ \hat{\bm{m}}^{\rm r}_{\rm e} \end{bmatrix} = \begin{bmatrix} 2\left(\hat{r}_1\hat{r}_3 - \hat{r}_0\hat{r}_2\right)g\\ 2\left(\hat{r}_2\hat{r}_3 + \hat{r}_0\hat{r}_1\right)g\\ (1 - 2\left(\hat{r}_1^2 + \hat{r}_2^2\right))g\\ 1 - 2\left(\hat{r}_2^2 + \hat{r}_3^2\right)\\ 2\left(\hat{r}_1\hat{r}_2 - \hat{r}_0\hat{r}_3\right)\\ 2\left(\hat{r}_1\hat{r}_3 + \hat{r}_0\hat{r}_2\right) \end{bmatrix} \end{align}
設計の便宜上、y^k\hat{\bm{y}}_{k} の算出に関して写像 hh を用いて次のように記述する。
y^k=h(x^k)\begin{align} \hat{\bm{y}}_{k} = h(\hat{\bm{x}}_{k}) \end{align}
以上の y,y^k\bm{y}, \hat{\bm{y}}_{k} を使用して予測残差 ek\bm{e}_{k} が算出される。この予測残差 ek\bm{e}_{k} を小さくするようにパラメータ x\bm{x} を修正するためガウス・ニュートン法を適用したいが、今回の設計条件では写像 hh のヤコビ行列には零空間が存在する。
Jhhx^R6×7\bm{J}_{h} \equiv \frac{\partial h}{\partial \hat{\bm{x}}} \in \mathbb{R}^{6\times 7}
そこで、ωb\bm{\omega}^{\rm b} について以下の仮定を置き、観測系のヤコビ行列のランクを下げる。
hωb0Jh[hq^0] \begin{align} \frac{\partial h}{\partial \bm{\omega}^{\rm b}} &\approx \bm{0} \\ \rightarrow \bm{J}_{h} &\approx \begin{bmatrix} \cfrac{\partial h}{\partial \hat{\bm{q}}} & \bm{0} \end{bmatrix} \end{align}
低次元化されたヤコビ行列 Jh\bm{J}_{h'} を用いて、予測誤差 ek\bm{e}_{k} から q\bm{q} に関する修正項 Δkq\bm{\Delta}^{q}_{k} を得る。
Jhhq^=[2gr^22gr^32gr^02gr^12gr^12gr^02gr^32gr^22gr^02gr^12gr^22gr^32r^02r^12r^22r^32r^32r^22r^12r^02r^22r^32r^02r^1]Jh+(JhTJh)1JhTΔkq=Jh+ek \begin{align} \bm{J}_{h'} &\equiv \frac{\partial h}{\partial \hat{\bm{q}}} = \begin{bmatrix} -2g \hat{r}_2 & 2g \hat{r}_3 & -2g \hat{r}_0 & 2g \hat{r}_1 \\ 2g \hat{r}_1 & 2g \hat{r}_0 & 2g \hat{r}_3 & 2g \hat{r}_2 \\ 2g \hat{r}_0 & -2g \hat{r}_1 & -2g \hat{r}_2 & 2g \hat{r}_3 \\ 2 \hat{r}_0 & 2 \hat{r}_1 & -2 \hat{r}_2 & -2 \hat{r}_3 \\ -2 \hat{r}_3 & 2 \hat{r}_2 & 2 \hat{r}_1 & -2 \hat{r}_0 \\ 2 \hat{r}_2 & 2 \hat{r}_3 & 2 \hat{r}_0 & 2 \hat{r}_1 \\ \end{bmatrix}\\ \bm{J}^{+}_{h'} &\equiv \left(\bm{J}_{h'}^{\mathrm T}\bm{J}_{h'}\right)^{-1}\bm{J}_{h'}^{\mathrm T} \\ \bm{\Delta}^{q}_{k} &= \bm{J}^{+}_{h'} \bm{e}_{k} \end{align}
続いて、q\bm{q}ωb\bm{\omega}^{\rm b} の関係から、ωb\bm{\omega}^{\rm b} の修正項 Δkω\bm{\Delta}^{\omega}_{k} を次のように算出する。
Δkω=Jω+ΔkqJω+(JωTJω)1JωTJωqω^kbTsJq(q^k)\begin{align} \bm{\Delta}^{\omega}_{k} &= \bm{J}^{+}_{\omega} \bm{\Delta}^{q}_{k} \\ \bm{J}^{+}_{\omega} &\equiv \left(\bm{J}_{\omega}^{\mathrm T}\bm{J}_{\omega}\right)^{-1}\bm{J}_{\omega}^{\mathrm T} \\ \bm{J}_{\omega} &\equiv \frac{\partial \bm{q}}{\partial \hat{\bm{\omega}}^{\rm b}_{k}} \fallingdotseq T_{\rm s}\bm{J}_{q}(\hat{\bm{q}}_{k}) \end{align}
以上を用いて、状態 x\bm{x} の修正項 Δ\bm{\Delta} を次のように設計する。
Δk=[ΔkqΔkω]R7\begin{align} \bm{\Delta}_{k} &= \begin{bmatrix} \bm{\Delta}^{q}_{k} \\ \bm{\Delta}^{\omega}_{k} \end{bmatrix} \in \mathbb{R}^{7} \end{align}
上記の修正項に縮小因子 γ\bm{\gamma} を組み合わせ、予測値 x~k+1\tilde{\bm{x}}_{k+1} を修正することで推定値 x^k+1\hat{\bm{x}}_{k+1} とする。
x^k+1=x~k+1+γΔk\begin{align} \hat{\bm{x}}_{k+1} = \tilde{\bm{x}}_{k+1} + \bm{\gamma}\bm{\Delta}_{k} \end{align}
縮小因子 γ\bm{\gamma} はオブザーバゲインに相当し,収束速度を決定するパラメータとなる。

© DigitalServo