DigitalServo

Service

Blog

About

Library

    線形制御系設計

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

    モーションコントロール

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

    モータドライブ

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

    システム同定

    ホワイトボックスモデル
    パラメトリックモデル

一括最小二乗法によるパラメータ同定

パラメトリックモデルはシステムに関して事前情報が存在する場合に採用される。 ここでは線形回帰モデルに対して一括最小二乗法を用いてパラメータ同定を行う方法について説明する。

問題定義

目的変数，説明変数，パラメータをそれぞれ y∈R, ϕ∈Rn, θ∈Rny\in\mathbb{R},\ \bm{\phi}\in\mathbb{R}^{\rm n},\ \bm{\theta}\in\mathbb{R}^{\rm n}とし，目的変数の決定に関して以下の線形回帰モデルを仮定する。
y^=θTϕ\begin{align} \hat{y} &= \bm{\theta}^{\mathrm T}\bm{\phi} \end{align}
ただし，上添字 ^\hat{} は推定値であることを表す。説明変数の選択およびその選択数などはシステムの事前情報を用いて決定する。例えば，時系列モデルを作成する場合には，システムの次数(設定したい極の個数など)から説明変数の数を決定することもできる。

一括最小二乗法

本手法では与えられたデータセットの全てを使用して評価関数を最小化するパラメータを求める。目的変数と説明変数が N\rm N セット取得できたとして，最小二乗法の評価関数を θ\bm{\theta} の関数として次のように定める。
J(θ)=1N∑ k=1N(y[k]−y^[k])2=1N∑ k=1N(y[k]−θTϕ[k])2\begin{align} J(\bm{\theta}) &= \frac{1}{\rm N}\sum^{\rm N}_{\ k= 1} \left(y[{\rm k}] - \hat{y}[{\rm k}] \right)^2\\ &= \frac{1}{\rm N}\sum^{\rm N}_{\ k= 1} \left(y[{\rm k}] - \bm{\theta}^{\mathrm T}\bm{\phi}[{\rm k}]\right)^2 \end{align}
すなわち，求めるパラメータ θ^\hat{\bm{\theta}} は次のように記述される。
θ^=arg min⁡θJ(θ)\begin{align} \hat{\bm{\theta}} = \argmin_{\bm{\theta}} J(\bm{\theta}) \end{align}
これを解くために，評価関数を整理して以下を得る。
J(θ)=1N∑ k=1Ny2[k]−2θTf+θTRθf=1N∑ k=1Nϕ[k]y[k]R=1N∑ k=1Nϕ[k]ϕT[k]\begin{align} J(\bm{\theta}) &= \frac{1}{\rm N}\sum^{\rm N}_{\ k= 1}y^2[{\rm k}] - 2\bm{\theta}^{\mathrm T}\bm{f} + \bm{\theta}^{\mathrm T}\bm{R} \bm{\theta} \\ \bm{f} &= \frac{1}{\rm N}\sum^{\rm N}_{\ k= 1}\bm{\phi}[{\rm k}]y[{\rm k}]\\ \bm{R} &= \frac{1}{\rm N}\sum^{\rm N}_{\ k= 1}\bm{\phi}[{\rm k}]\bm{\phi}^{\mathrm T}[{\rm k}] \end{align}
この関数に停留点が存在し，かつその停留点におけるヘッセ行列が正定値であれば，その停留点が θ^\hat{\bm{\theta}} となる。評価関数の勾配は次のように与えられる。
∇J(θ)=∂∂θJ(θ)=−2f+2Rθ\begin{align} \nabla J(\theta) = \frac{\partial}{\partial\bm{\theta}}J(\bm{\theta}) &= - 2\bm{f} + 2\bm{R} \bm{\theta} \end{align}
以上より，θ^\hat{\bm{\theta}} が満たすべき正規方程式が求まる。
Rθ^=f\begin{align} \bm{R} \hat{\bm{\theta}} &= \bm{f} \end{align}
また，停留点におけるヘッセ行列は以下のように求められる。
∂∂θT∇J(θ^)=R\begin{align} \frac{\partial}{\partial\bm{\theta}^{\mathrm T}}\nabla J(\hat{\bm{\theta}}) &= \bm{R} \end{align}
すなわち，R\bm{R} が正定値行列であれば評価関数は停留点で最小値を取り，また正則であるため正規方程式を解くことができる。
θ^=R−1f\begin{align} \hat{\bm{\theta}} = \bm{R}^{-1}\bm{f} \end{align}

ARXモデルへの適用

システムをARXモデルで表現することについて考える。
A(z)y[k]=B(z)u[k]+w[k]A(z)≡1+∑k=1naa[k]z−kB(z)≡∑k=1nbb[k]z−k\begin{align} A(z)y[{\rm k}] &= B(z)u[{\rm k}] + w[{\rm k}] \\ A(z) &\equiv 1 + \sum_{{\rm k}=1}^{\rm n_{\rm a}} a[{\rm k}]z^{-{\rm k}}\\ B(z) &\equiv \sum_{{\rm k}=1}^{\rm n_{\rm b}} b[{\rm k}]z^{-{\rm k}}\\ \end{align}
ただし，u,y,wu, y, w はそれぞれ入力，出力，雑音を表し，ai,bj (i∈[1,na], j∈[1,nb])a_{\rm i}, b_{\rm j}\ ({\rm i}\in[1,{\rm n_{\rm a}}],\ {\rm j}\in[1,{\rm n_{\rm b}}]) は回帰モデルのパラメータである。これを差分方程式の形式とするために整理すれば，以下の表現を得る。
y[k]=[1−A(z)]y[k]+B(z)u[k]+w[k]=−∑k=1na[k]z−ky[k]+∑k=1mb[k]z−ku[k]+w[k]=θaTϕy[k]+θbTϕu[k]+w[k]θa≡[a1a2⋯an]Tθb≡[b1b2⋯bm]Tϕy[k]≡[−y[k−1]−y[k−2]⋯−y[k−n]]Tϕu[k]≡[u[k−1]u[k−2]⋯u[k−m]]T\begin{align} y[{\rm k}] &= \left[1-A(z)\right]y[{\rm k}] + B(z)u[k] + w[k]\\ &= -\sum_{{\rm k}=1}^{\rm n} a[{\rm k}]z^{-{\rm k}}y[{\rm k}] + \sum_{{\rm k}=1}^{\rm m} b[{\rm k}]z^{-{\rm k}}u[{\rm k}] + w[{\rm k}] \\ &= \bm{\theta}_{a}^{\mathrm T}\bm{\phi}_{y}[{\rm k}] + \bm{\theta}_{b}^{\mathrm T}\bm{\phi}_{u}[{\rm k}] + w[{\rm k}]\\ \bm{\theta}_{a} &\equiv \begin{bmatrix} a_{1} & a_{2} & \cdots & a_{\rm n} \end{bmatrix}^{\mathrm T} \\ \bm{\theta}_{b} &\equiv \begin{bmatrix} b_{1} & b_{2} & \cdots & b_{\rm m} \end{bmatrix}^{\mathrm T} \\ \bm{\phi}_{\rm y}[{\rm k}] &\equiv \begin{bmatrix} -y[{\rm k}-1] & -y[{\rm k}-2] & \cdots & -y[{\rm k}-{\rm n}] \end{bmatrix}^{\mathrm T} \\ \bm{\phi}_{\rm u}[{\rm k}] &\equiv \begin{bmatrix} u[{\rm k}-1] & u[{\rm k}-2] & \cdots & u[{\rm k}-{\rm m}] \end{bmatrix}^{\mathrm T} \\ \end{align}
説明変数 ϕy,ϕu\bm{\phi}_{y}, \bm{\phi}_{u} は前時刻までの利用可能な状態量から構成されるベクトルとなる。これをさらに整理すると，以下を得る。
y[k]=θTϕ+w[k]θ≡[θaθb]∈Rna+nbϕ[k]≡[ϕy[k]ϕu[k]]∈R(na+nb)×(na+nb)\begin{align} y[{\rm k}] &= \bm{\theta}^{\mathrm T}\bm{\phi} + w[{\rm k}]\\ \bm{\theta} &\equiv \begin{bmatrix} \bm{\theta}_{a} \\ \bm{\theta}_{b} \end{bmatrix} \in \mathbb{R}^{{\rm n}_{\rm a}+{\rm n}_{\rm b}}\\ \bm{\phi}[{\rm k}] &\equiv \begin{bmatrix} \bm{\phi}_{\rm y}[{\rm k}] \\ \bm{\phi}_{\rm u}[{\rm k}] \end{bmatrix} \in \mathbb{R}^{({\rm n}_{\rm a}+{\rm n}_{\rm b})\times({\rm n}_{\rm a}+{\rm n}_{\rm b})} \end{align}
これは定式化された形式の表現であるため，一括最小二乗法によるパラメータ推定を行うことができる。

可同定性

パラメータ θ\bm{\theta} が推定可能な条件は，R\bm{R} が正定値となることである。ARXモデルの同定において，R\bm{R} は次の値を取る。
R=1N∑k=1N[ϕy[k]ϕyT[k]ϕy[k]ϕuT[k]ϕu[k]ϕyT[k]ϕu[k]ϕuT[k]]=[1N∑k=1Nϕy[k]ϕyT[k]1N∑k=1Nϕy[k]ϕuT[k]1N∑k=1Nϕu[k]ϕyT[k]1N∑k=1Nϕu[k]ϕuT[k]]≡[RyRyuRuyRu]\begin{align} \bm{R} &= \frac{1}{\rm N}\sum^{\rm N}_{{\rm k}=1} \begin{bmatrix} \bm{\phi}_{\rm y}[{\rm k}]\bm{\phi}^{\mathrm T}_{\rm y}[{\rm k}] & \bm{\phi}_{\rm y}[{\rm k}]\bm{\phi}^{\mathrm T}_{\rm u}[{\rm k}] \\ \bm{\phi}_{\rm u}[{\rm k}]\bm{\phi}^{\mathrm T}_{\rm y}[{\rm k}] & \bm{\phi}_{\rm u}[{\rm k}]\bm{\phi}^{\mathrm T}_{\rm u}[{\rm k}] \\ \end{bmatrix} \\ &= \begin{bmatrix} \frac{1}{\rm N}\sum^{\rm N}_{{\rm k}=1}\bm{\phi}_{\rm y}[{\rm k}]\bm{\phi}^{\mathrm T}_{\rm y}[{\rm k}] & \frac{1}{\rm N}\sum^{\rm N}_{{\rm k}=1}\bm{\phi}_{\rm y}[{\rm k}]\bm{\phi}^{\mathrm T}_{\rm u}[{\rm k}] \\ \frac{1}{\rm N}\sum^{\rm N}_{{\rm k}=1}\bm{\phi}_{\rm u}[{\rm k}]\bm{\phi}^{\mathrm T}_{\rm y}[{\rm k}] & \frac{1}{\rm N}\sum^{\rm N}_{{\rm k}=1}\bm{\phi}_{\rm u}[{\rm k}]\bm{\phi}^{\mathrm T}_{\rm u}[{\rm k}] \\ \end{bmatrix}\\ &\equiv \begin{bmatrix} \bm{R}_{y} & \bm{R}_{yu} \\ \bm{R}_{uy} & \bm{R}_{u} \\ \end{bmatrix} \end{align}
ただし，Ru,Ry\bm{R}_{u}, \bm{R}_{y} は自己相関行列，Ruy,Ryu\bm{R}_{uy}, \bm{R}_{yu} は相互相関行列であり，Ruy=RyuT\bm{R}_{uy} = \bm{R}_{yu}^{\mathrm T} を満たす。この行列が正定値となる条件は，シューア補行列 R/Ru\bm{R}/\bm{R}_{u} を用いて次のように記述される。
R≻0⇔Ru≻0 ∧ R/Ru≻0R/Ru=Ry−RyuRu−1Ruy=Cov(y∣u)\begin{align} \bm{R} \succ 0 &\Leftrightarrow \bm{R}_{u} \succ 0\ \land\ \bm{R}/\bm{R}_{u} \succ 0 \\ \bm{R}/\bm{R}_{u} &= \bm{R}_{y} - \bm{R}_{yu}\bm{R}_{u}^{-1}\bm{R}_{uy} = {\rm Cov}(y|u) \end{align}
ここで，Ry\bm{R}_{y} が半正定値となると補行列も半正定値となるため，Ry\bm{R}_{y} も正定値となる必要がある。

PE性

この R\bm{R} の正定値性については持続的励振性 (Persistence of excition: PE) を基に議論がされている[1-4]。PE性は信号の情報の豊かさを示す指標であり，入力が豊かさを持つことでシステムの応答にも豊かさが与えられる。システムへの入力がPE性信号である場合，システムの持つ全てのモードが励起されることで R\bm{R} の列ベクトルが線型独立となり，正定値となる。

文献[1]では時間領域においてPE性信号の条件とその次数を定義しており，信号 uu についてその自己相関関数 ru(τ)r_{u}(\tau) を用いた以下の行列 Rm\bm{\mathcal{R}}_{\rm m} が正定値であり，かつ Rm+1\bm{\mathcal{R}}_{\rm m+1} が特異であるとき，信号は次数 m\rm m のPE性信号と定義している。
ru(T)=lim⁡N→∞1N∑t=1Nu(t)u(t+T)Rm=[ru(0)ru(1)⋯ru(m−1)ru(1)ru(0)⋯ru(m−2)⋮⋮⋱⋮ru(m−1)ru(m−2)⋯ru(0)]\begin{align} r_{u}(T) &= \lim_{{\rm N}\rightarrow\infty} \frac{1}{\rm N}\sum^{\rm N}_{t=1} u(t)u(t+T) \\ \bm{\mathcal{R}}_{\rm m} &= \begin{bmatrix} r_{u}(0) & r_{u}(1) & \cdots & r_{u}({\rm m}-1) \\ r_{u}(1) & r_{u}(0) & \cdots & r_{u}({\rm m}-2) \\ \vdots & \vdots & \ddots & \vdots\\ r_{u}({\rm m}-1) & r_{u}({\rm m}-2) & \cdots & r_{u}(0) \\ \end{bmatrix} \end{align}
Ljung[2, 4]は周波数領域におけるPE性信号の条件とその次数について与えており，期待値と自己相関が時間不変である信号 (準定常信号) uu についてスペクトル密度 Φu(ω)\Phi_{u}(\omega) が m\rm m 点で非零となるとき，その信号を次数 m\rm m のPE性信号と定義している。スペクトル密度は自己相関関数のフーリエ変換であるため，本質的に双方は等価である。

正弦波のPE性

次の単一周波数を持つ信号について考える。
u(t)=sin⁡(ω0t)\begin{align} u(t)=\sin(\omega_{0}t) \end{align}
この信号の自己相関関数は次のようになる。
ru(τ)=E[sin⁡(ω0t)sin⁡(ω0t−τ)]=E[12(cos⁡(ω0τ)−cos⁡(2ω0t−ω0τ))]=12E[cos⁡(ω0τ)]−12E[cos⁡(2ω0t−ω0τ))]\begin{align} r_{u}(\tau) &= E[\sin(\omega_{0}t)\sin(\omega_{0}t - \tau)] \\ &= E\left[\frac{1}{2} \left(\cos(\omega_{0}\tau) - \cos(2\omega_{0}t - \omega_{0}\tau)\right) \right] \\ &= \frac{1}{2} E\left[\cos(\omega_{0}\tau)\right] - \frac{1}{2}E\left[\cos(2\omega_{0}t - \omega_{0}\tau)) \right] \end{align}
右辺第一項の cos⁡(ω0τ)\cos(\omega_{0}\tau) は定数であること，第二項の周期関数のため期待値が 00 となることを考慮すると，以下を得る。
ru(τ)=12cos⁡(ω0τ)\begin{align} r_{u}(\tau) &= \frac{1}{2} \cos(\omega_{0}\tau) \end{align}
したがって，スペクトル密度は次のようになる。
Φu(ω)=∑τ=−∞∞ru(τ)e−jωτ=∑τ=−∞∞12cos⁡(ω0τ)e−jωτ=∑τ=−∞∞12⋅12(ejω0τ+e−jω0τ)e−jωτ=14∑τ=−∞∞(e−j(ω−ω0)τ+e−j(ω+ω0)τ)=14{D∞(−(ω−ω0))+D∞(−(ω+ω0))}\begin{align} \Phi_{u}(\omega) &= \sum_{\tau=-\infty}^{\infty}r_{u}(\tau) e^{-j\omega\tau} = \sum_{\tau=-\infty}^{\infty} \frac{1}{2} \cos(\omega_{0}\tau) e^{-j\omega\tau}\\ &= \sum_{\tau=-\infty}^{\infty} \frac{1}{2}\cdot \frac{1}{2}(e^{j\omega_{0}\tau} + e^{-j\omega_{0}\tau}) e^{-j\omega\tau} \\ &= \frac{1}{4}\sum_{\tau=-\infty}^{\infty} \left(e^{-j(\omega - \omega_{0})\tau} + e^{-j(\omega + \omega_{0})\tau}\right) \\ &= \frac{1}{4} \left\{D_{\infty}(-(\omega - \omega_{0})) + D_{\infty}(-(\omega + \omega_{0})) \right\} \end{align}
ただし，DD はディリクレ核とする。ディリクレ核はデルタ収束関数であり，スペクトル密度を次のように近似表現できる。
Φu(ω)∼14{2πδ(−(ω−ω0))+2πδ(−(ω+ω0))}=π2{δ(ω−ω0)+δ(ω+ω0)}\begin{align} \Phi_{u}(\omega) &\sim \frac{1}{4} \left\{2\pi\delta(-(\omega - \omega_{0})) + 2\pi\delta(-(\omega + \omega_{0})) \right\} \\ &= \frac{\pi}{2} \left\{\delta(\omega - \omega_{0}) + \delta(\omega + \omega_{0}) \right\} \end{align}
このスペクトル密度は ω=±ω0\omega=\pm\omega_{0} の2点において非零となるため，単一周波数を持つ正弦波のPE性次数は2となる。また，以下に示す n\rm n 個の周波数成分を持つマルチサイン信号は次数 2n2{\rm n} のPE性信号となる。
u(t)=∑k=1nμksin⁡(ωkt+ψk), ωi≠ωj,i≠j,ωi≠0,ωi≠π\begin{align} u(t) &= \sum^{\rm n}_{\rm k=1}\mu_{\rm k}\sin(\omega_{\rm k}t + \psi_{\rm k}),\ \omega_{i}\neq\omega_{\rm j}, {\rm i}\neq{\rm j}, \omega_{\rm i}\neq 0, \omega_{\rm i}\neq \pi \end{align}
ただし，μ,ψ\mu, \psi は振幅と位相を表す。
システム同定を行うために必要なPE性の次数については，R\bm{R} の構造について確認する必要がある。R\bm{R} を区分行列として表現した場合，対角成分には Ry\bm{R}_{y} と Ru\bm{R}_{u} が現れるが，先に示した通り R\bm{R} を正定値とするためにはこれらの行列を正定値都する必要がある。文献[3]では，入力信号 uu が十分なPE性を持つ場合，システムを通して出力 yy にもPE性が伝播され，R\bm{R} も正定値となることが示されている。このとき，PE性の次数は R\bm{R}のランクを保つため，na+nb\rm n_{\rm a} + n_{\rm b} 以上にすべきとの指標がある。ただし，情報豊富な信号は最小二乗法の安定性を高めるため，多くの信号を入力することが望ましい。

参考文献

[1]K.J. Åström, B. Torsten, "Numerical Identification of Linear Dynamic Systems from Normal Operating Records," in IFAC Proceedings Volumes, vol. 2, no. 2, pp. 96-111, Sept. 1965.
[2]L. Ljung, "Characterization of the Concept of `Persistently Exciting` in the Frequency Domain," in Department of Automatic Control, Lund Institute of Technology, no. Research Reports TFRT-3038, 1971.
[3]M. Green, J. B. Moore, "Persistence of excitation in linear systems," in Systems & Control Letters, vol. 7, no. 5, pp. 351-360, Mar. 1986.
[4]L. Ljung, "System Identification: Theory for the User (2nd ed.)," in Prentice Hall, Information and System Sciences Series, pp. 408-457 (Chapter 13), 1999.

CONTENTS