実際に使える制御工学：寄り道編（４）「システム・パラメータ同定。１回につき100万円の技術。」

＜この記事で伝えたい事＞ 

• システム・パラメータ同定は、制御対象をハッキリさせる技術です
  制御設計において、最初にやるべきことは制御対象をハッキリさせることです。
  まずは、数学・物理学知識を用いて制御対象を導きだすのがセオリーですが、実験データに基づいて詳細に制御対象特性を把握したいときは、システム/パラメータ同定をします。 

1. はじめに

自己紹介はこちらです。
 
この記事は、制御工学の周辺技術としてのシステム・パラメータ同定技術の紹介・解説のために書いています。

2. システム・パラメータ同定とは？

システム同定は、実際に計測した実験データに基づき、システムの数理モデル（伝達関数や状態関数）を導き出す技術です。

パラメータ同定は、数理モデル内の未知パラメータを導き出す技術です。
 

運動方程式や電圧方程式の形は分かるけど、摩擦などの数値が分からない（または数値の精度に不安がある）場合、パラメータ同定が使われます。
 
一方で、運動方程式や電圧方程式の形からして
「本当にこれで良いのか？」
と不安な場合は、システム同定が使われます。

ただ、同定したシステムから物理的な意味を読み取るのは難度が高いです。
それには、ある程度の知識と経験が必要です。

素人には、そうそうできることではありません。
（おそらく、単価が高いのはこのためです）

3. フーリエ変換によるシステム同定

実験データを測定し、それに対して
「フーリエ変換」（正確には「離散フーリエ変換」）
と呼ばれる信号処理をかけると、そのデータがどの周波数の成分をどれだけ持つかが分かります。
（フーリエ変換の詳細については、付録1,2で説明します）
 
フーリエ変換の結果は、下式のように複素数（＝実数＋虚数）の形で得られます。
（$${j}$$は虚数、$${n}$$は整数、$${\omega}$$は角周波数） 

$$
F_{(\omega_n)}= R_{(\omega_n)}+j\cdot I_{(\omega_n)}
$$

よって、角周波数$${\omega_n}$$における大きさ$${|F_(\omega_n)|}$$・位相$${\angle F_(\omega_n)}$$は、それぞれ下式で求められます。

$$
\begin{array}{ll}
|F_{(\omega_n)}|&=\sqrt{{R_{(\omega_n)}}^2+{I_{(\omega_n)}}^2}\\
\angle F_{(\omega_n)}&=tan^{-1}\left(\cfrac{I_{(\omega_n)}}{R_{(\omega_n)}}\right)
\end{array}
$$

そのデータがどの周波数成分を多く含むかさえ分かれば良いときは、位相の計算は省略できます。
（これは周波数スペクトル解析と呼ばれます）
 

システムの伝達関数を明らかにしたいときは、大きさ・位相を両方計算した方が良いです。

システム入力データを$${u}$$、出力データを$${y}$$とし、それぞれの大きさ・位相を$${|F_{u(\omega_n)}|}$$，$${\angle F_{u(\omega_n)}}$$，$${|F_{y(\omega_n)}|}$$，$${\angle F_{y(\omega_n)}}$$とします。
このとき、システム伝達関数$${G}$$のゲイン$${|G|}$$と、位相$${\angle G}$$は、以下のように求められます。
（ゲインの単位を[dB]にするには、$${20\cdot log10(|G|)}$$を計算する必要あり）

$$
\begin{array}{ll}
|G_{(\omega_n)}|&=\cfrac{|F_{y(\omega_n)}|}{|F_{u(\omega_n)}|}\\
\angle G_{(\omega_n)}&=\angle F_{y(\omega_n)}-\angle F_{u(\omega_n)}
\end{array}
$$

それぞれの角周波数$${\omega_n}$$について計算すれば、システム伝達関数の特性が数値的に分かります。

ただ、実験データといえど、以下の要因がノイズになります。

• センサ絡みの特性（絶対・相対精度やAD変換分解能）

• データの桁落ち（特に、固定小数点演算のとき）

 これらの影響で、多少は真の特性に対してズレが出たり、キレイな曲線が得られなかったりします。

キレイな曲線が欲しい場合、得られた計算結果に対しカーブフィッティングするのをオススメします。

4. 最小二乗法によるパラメータ同定

例えば、モータの運動方程式が下式なのは分かっているとします。

$$
J_m\cdot \ddot{\theta}+D_m\cdot \dot{\theta}=\tau
$$

ここで、
トルク$${\tau}$$、角加速度$${\ddot{\theta}}$$、角速度$${\dot{\theta}}$$は測定可能（または測定値から計算可能）なものの、
モータイナーシャ$${J_m}$$や粘性摩擦係数$${D_m}$$の正確な数値は分からない
とします。

そんなときは、パラメータ同定の出番です。

まずは、モータ動作中に$${\tau,\ddot{\theta},\dot{\theta}}$$のデータを$${n}$$回測定して、それをただ縦に並べます。
すると、以下のようになります。
（$${n}$$は整数）

$$
J_m\cdot \ddot{\theta}_{(1)}+D_m\cdot \dot{\theta}_{(1)}=\tau_{(1)}\\
J_m\cdot \ddot{\theta}_{(2)}+D_m\cdot \dot{\theta}_{(2)}=\tau_{(2)}\\
\vdots\\
J_m\cdot \ddot{\theta}_{(n)}+D_m\cdot \dot{\theta}_{(n)}=\tau_{(n)}\\
$$

動作中に$${J_m,D_m}$$は変わらないと仮定します。
そして、行列を使ってまとめると下式になります。

$$
\left[
\begin{array}{cc}
\ddot{\theta}_{(1)}&\dot{\theta}_{(1)}\\
\ddot{\theta}_{(2)}&\dot{\theta}_{(2)}\\
\vdots&\vdots\\
\ddot{\theta}_{(n)}&\dot{\theta}_{(n)}\\
\end{array}
\right]
\left[
\begin{array}{c}
J_m\\
D_m
\end{array}
\right]
=
\left[
\begin{array}{c}
\tau_{(1)}\\
\tau_{(2)}\\
\vdots\\
\tau_{(n)}\\
\end{array}
\right]
$$

ここで、$${\ddot{\theta}}$$と$${\dot{\theta}}$$で構成される行列の転置行列を、両辺の左側から掛けます。
（転置行列：その行列の行と列を入れ替えた行列）
その結果、下式が得られます。

$$
\left[
\begin{array}{cccc}
\ddot{\theta}_{(1)}&\ddot{\theta}_{(2)}&\cdots&\ddot{\theta}_{(n)}\\
\dot{\theta}_{(1)}&\dot{\theta}_{(2)}&\cdots&\dot{\theta}_{(n)}
\end{array}
\right]
\left[
\begin{array}{cc}
\ddot{\theta}_{(1)}&\dot{\theta}_{(1)}\\
\ddot{\theta}_{(2)}&\dot{\theta}_{(2)}\\
\vdots&\vdots\\
\ddot{\theta}_{(n)}&\dot{\theta}_{(n)}\\
\end{array}
\right]
\left[
\begin{array}{c}
J_m\\
D_m
\end{array}
\right]
=
\left[
\begin{array}{cccc}
\ddot{\theta}_{(1)}&\ddot{\theta}_{(2)}&\cdots&\ddot{\theta}_{(n)}\\
\dot{\theta}_{(1)}&\dot{\theta}_{(2)}&\cdots&\dot{\theta}_{(n)}
\end{array}
\right]
\left[
\begin{array}{c}
\tau_{(1)}\\
\tau_{(2)}\\
\vdots\\
\tau_{(n)}\\
\end{array}
\right]
$$

$$
\begin{array}{l}
\left[
\begin{array}{cc}
({\ddot{\theta}_{(1)}}^2+{\ddot{\theta}_{(2)}}^2+\cdots+{\ddot{\theta}_{(n)}}^2)
&(\ddot{\theta}_{(1)}\cdot\dot{\theta}_{(1)}+\ddot{\theta}_{(2)}\cdot\dot{\theta}_{(2)}+\cdots+\ddot{\theta}_{(n)}\cdot\dot{\theta}_{(n)})\\
(\ddot{\theta}_{(1)}\cdot\dot{\theta}_{(1)}+\ddot{\theta}_{(2)}\cdot\dot{\theta}_{(2)}+\cdots+\ddot{\theta}_{(n)}\cdot\dot{\theta}_{(n)})
&({\dot{\theta}_{(1)}}^2+{\dot{\theta}_{(2)}}^2+\cdots+{\dot{\theta}_{(n)}}^2)
\end{array}
\right]
\left[
\begin{array}{c}
J_m\\
D_m
\end{array}
\right]\\
=
\left[
\begin{array}{c}
(\ddot{\theta}_{(1)}\cdot\tau_{(1)}+\ddot{\theta}_{(2)}\cdot\tau_{(2)}+\cdots+\ddot{\theta}_{(n)}\cdot\tau_{(n)})\\
(\dot{\theta}_{(1)}\cdot\tau_{(1)}+\dot{\theta}_{(2)}\cdot\tau_{(2)}+\cdots+\dot{\theta}_{(n)}\cdot\tau_{(n)})
\end{array}
\right]
\end{array}
$$

$$
\left[
\begin{array}{cc}
\displaystyle\sum_{k=1}^n\left({\ddot{\theta}_{(k)}}^2\right)
&
\displaystyle\sum_{k=1}^n\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)
\\
\displaystyle\sum_{k=1}^n\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)
&
\displaystyle\sum_{k=1}^n\left({\dot{\theta}_{(k)}}^2\right)
\end{array}
\right]
\left[
\begin{array}{c}
J_m\\
D_m
\end{array}
\right]
=
\left[
\begin{array}{c}
\displaystyle\sum_{k=1}^n\left(\ddot{\theta}_{(k)}\cdot\tau_{(k)}\right)
\\
\displaystyle\sum_{k=1}^n\left(\dot{\theta}_{(k)}\cdot\tau_{(k)}\right)
\end{array}
\right]
$$

ある行列に対して、その転置行列をかけると、必ず正方行列になります。
（正方行列：縦と横のサイズが同じ行列）

しかし、正方行列が、必ず正則行列だとは断言できません。
（正則行列：逆行列のある正方行列）

しかし、それなりに多くのデータを使えば、上式中の正方行列はまず正則になります。
よって、以下のように未知パラメータ$${J_m,D_m}$$が得られます。

$$
\left[
\begin{array}{c}
J_m\\
D_m
\end{array}
\right]
=
\left[
\begin{array}{cc}
\displaystyle\sum_{k=1}^n\left({\ddot{\theta}_{(k)}}^2\right)
&
\displaystyle\sum_{k=1}^n\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)
\\
\displaystyle\sum_{k=1}^n\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)
&
\displaystyle\sum_{k=1}^n\left({\dot{\theta}_{(k)}}^2\right)
\end{array}
\right]^{-1}
\left[
\begin{array}{c}
\displaystyle\sum_{k=1}^n\left(\ddot{\theta}_{(k)}\cdot\tau_{(k)}\right)
\\
\displaystyle\sum_{k=1}^n\left(\dot{\theta}_{(k)}\cdot\tau_{(k)}\right)
\end{array}
\right]
$$

運動方程式や電圧方程式中の未知パラメータの数が増えても、同じ要領でパラメータを同定できます。

なお、ここでは$${J_m,D_m}$$などの未知パラメータは動作中に変わらないとしましたが、実際には変わるケースもあります。

そうした場合にも対応できる方法として、逐次最小二乗法と呼ばれる技術があります。
（詳細については、付録3で説明します）

5. おわりに

この記事では、システム/パラメータ同定を紹介・解説しました。
 
世の中には、システムの数理モデルやパラメータが不確かであっても、高精度な制御を実現する方法を研究する方もおられます。
（モデルフリー制御など）

きっと、そこに学術的な価値はあるのだろうと思います。
 
しかし、今のところは、制御対象がハッキリと分かっていた方が、圧倒的に高い性能を実現できます。
遠い未来はともかく、今後100年くらいは制御技術にとって、システム/パラメータ同定は重要な周辺技術であり続けると思います。

なお、記事タイトルの「１回につき100万円」は誇張ではありません。

私が学生のころ、先生から
「今度のシステム同定には100万円くらいかかる（から成果出してね）」
と言われたことがあります。

2026年現在は、当時より物価上昇しており、技術サービスの価格も上がっているはずです。
これを考慮すると、現在は100万円以上かかってもおかしくありません。

前回：寄り道編（３）「状態方程式」はこちら
次回：寄り道編（５）「システム同定時の入力信号」はこちら

付録1. フーリエ変換

フーリエ変換とは、あらゆる周期波形をSIN波の足し合わせに分解し、どの周波数のSIN波の大きさがどれだけなのかを明らかにする方法です。
 
「あらゆる周期波形」とは、
「変換対象の波形は自由に決められるけど、それが周期的に繰り返されるものとして扱われますよ」
という意味の言葉と思ってもらえば、それでだいたい合っています。

ちなみに、変換対象波形の最初と最後の数値が違うと、急激な波形の立ち上がり（または立ち下がり）があるとして、フーリエ変換することになります。

すると、本来より高周波数の成分を多く含む結果が出やすいです。

窓関数というものを使えば、強制的に最初と最後の数値を０にそろえて、この問題を回避できます。

詳細説明は省略しますが、知っておけば役立つことがあると思います。

話を本題に戻します。

フーリエ変換を理解するには、
「どんな波形でも、SIN波の足し合わせで表現できる」
という結論ありきで考えると分かりやすいです。
 
このとき、我々が解析対象とする波形$${f_{(t)}}$$は、下式のように表せます。
（$${T}$$は周期、$${K}$$は整数、オイラーの公式より$${e^{\pm j\cdot\theta}=cos\theta\pm j\cdot sin\theta}$$）

$$
\begin{array}{ll}
f_{(t)}&=\displaystyle\sum_{K=-\infty}^{\infty}
\left(
c_{(K)}\cdot e^{j\cdot\frac{2\pi\cdot K}{T}\cdot t}
\right)\\
&=\displaystyle\sum_{K=-\infty}^{\infty}
\left(
c_{(K)}\cdot
\left(
cos\frac{2\pi\cdot K}{T}\cdot t + j\cdot sin\frac{2\pi\cdot K}{T}\cdot t
\right)
\right)
\end{array}
$$

我々が求めたいのは、係数$${c_{(K)}}$$です。
これは、下式によって求められます。
（$${k}$$は整数、$${c_{(k)}}$$は未知だが定数）

$$
\begin{array}{ll}
c_{(k)}&=\cfrac{1}{T}\cdot\displaystyle\int_0^T f_{(t)}\cdot e^{-j\cdot\frac{2\pi\cdot k}{T}\cdot  t}dt\\
&=
\cfrac{1}{T}\cdot
\displaystyle\int_0^T f_{(t)}\cdot
\left(
cos\frac{2\pi\cdot k}{T}\cdot t + j\cdot sin\frac{2\pi\cdot k}{T}\cdot t
\right)dt
\end{array}
$$

$${k}$$は無数にある$${K}$$のうちの１つです。
$${k=1,2,3,\cdots}$$と変えていくことで、それに応じた角周波数$${\omega=\frac{2\cdot \pi\cdot k}{T}}$$のSIN・COS波にかかる係数$${c_{(k)}}$$を求められます。

ここからは、なぜそうなるかを解説します。
これまでの数式から、下式が成り立ちます。

$$
\begin{array}{ll}
c_{(k)}&=\cfrac{1}{T}\cdot\displaystyle\int_0^T f_{(t)}\cdot e^{-j\cdot\frac{2\pi\cdot k}{T}\cdot  t}dt\\
&=
\cfrac{1}{T}\cdot\displaystyle\int_0^T
\displaystyle\sum_{K=-\infty}^{\infty}
\left(
c_{(K)}\cdot e^{j\cdot\frac{2\pi\cdot K}{T}\cdot t}
\right)
\cdot e^{-j\cdot\frac{2\pi\cdot k}{T}\cdot  t}dt\\
&=
\cfrac{1}{T}\cdot\displaystyle\int_0^T
\displaystyle\sum_{K=-\infty}^{\infty}
\left(
c_{(K)}\cdot e^{j\cdot\frac{2\pi\cdot (K-k)}{T}\cdot t}
\right)dt\\
&=
\cfrac{1}{T}\cdot\displaystyle\int_0^T
\left(
c_{(k)}\cdot e^{j\cdot\frac{2\pi\cdot (k-k)}{T}\cdot t}
\right)dt\\
&     +\cfrac{1}{T}\cdot\displaystyle\int_0^T
\displaystyle\sum_{K=-\infty}^{k-1}
\left(
c_{(K)}\cdot e^{j\cdot\frac{2\pi\cdot (K-k)}{T}\cdot t}
\right)dt\\
&          +\cfrac{1}{T}\cdot\displaystyle\int_0^T
\displaystyle\sum_{K=k+1}^{\infty}
\left(
c_{(K)}\cdot e^{j\cdot\frac{2\pi\cdot (K-k)}{T}\cdot t}
\right)dt
\end{array}
$$

上式右辺は、$${K=k}$$の項と、$${K\neq k}$$の項に大別できます。
まずは、$${K=k}$$の項に着目します。

$$
\begin{array}{ll}
\cfrac{1}{T}\cdot\displaystyle\int_0^T
\left(
c_{(k)}\cdot e^{j\cdot\frac{2\pi\cdot (k-k)}{T}\cdot t}
\right)dt
&=
\cfrac{1}{T}\cdot\displaystyle\int_0^T
\left(
c_{(k)}\cdot e^{j\cdot\frac{2\pi\cdot 0}{T}\cdot t}
\right)dt\\
&=\cfrac{1}{T}\cdot\displaystyle\int_0^T
\left(
c_{(k)}\cdot e^{0}
\right)dt\\
&=\cfrac{1}{T}\cdot\displaystyle\int_0^T
c_{(k)}dt\\
&=\cfrac{c_{(k)}}{T}\cdot[t]_0^T\\
&=\cfrac{c_{(k)}}{T}\cdot(T-0)\\
&=c_{(k)}
\end{array}
$$

最終的に、$${K=k}$$の項はそのまま$${c_{(k)}}$$です。

となれば、$${K\neq k}$$の項が０になってくれれば、正しく$${c_{(k)}}$$が求められることになります。

さしあたり、$${K=(k+1)～\infty}$$の方に着目してみましょう。

$$
\begin{array}{ll}
&\cfrac{1}{T}\cdot\displaystyle\int_0^T
\displaystyle\sum_{K=k+1}^{\infty}
\left(
c_{(K)}\cdot e^{j\cdot\frac{2\pi\cdot (K-k)}{T}\cdot t}
\right)dt\\
&     =\cfrac{1}{T}\cdot\displaystyle\int_0^T
\left(c_{(k+1)}\cdot(cos\frac{2\pi}{T}\cdot t+j\cdot sin\frac{2\pi}{T}\cdot t)
\right)dt\\
&               +\cfrac{1}{T}\cdot\displaystyle\int_0^T
\left(c_{(k+2)}\cdot(cos\frac{4\pi}{T}\cdot t+j\cdot sin\frac{4\pi}{T}\cdot t)
\right)dt\\
&                    +\cdots\\
&     =\cfrac{c_{(k+1)}}{T}\cdot\displaystyle\int_0^T
\left(cos\frac{2\pi}{T}\cdot t+j\cdot sin\frac{2\pi}{T}\cdot t
\right)dt\\
&               +\cfrac{c_{(k+2)}}{T}\cdot\displaystyle\int_0^T
\left(cos\frac{4\pi}{T}\cdot t+j\cdot sin\frac{4\pi}{T}\cdot t
\right)dt\\
&                    +\cdots\\
&     =0
\end{array}
$$

ここで重要なのは、$${K,k}$$がいずれも整数なことです。
$${K,k}$$がいずれも整数だと、上式中のSIN,COSの定積分の区間は、必ずキッチリ整数周期分になります。

下図は、２周期分のSIN波を表したものです。

積分とは、波形の面積を求める計算です。
このため、SIN波をキッチリ２周期分の区間で定積分すると、その結果は０になります。
（青部分の面積はプラス・赤部分の面積はマイナス扱いなので）

積分区間を２周期から３周期に変えても、定積分の結果は０です。
（積分区間がキッチリ整数周期分なら、SIN,COSの定積分値は０）

このため、$${K=(k+1)～\infty}$$の項＝０です。
同様に、$${K=-\infty～(k-1)}$$の項＝０です。

以上のことから、

$$
\begin{array}{ll}
c_{(k)}&=\cfrac{1}{T}\cdot\displaystyle\int_0^T f_{(t)}\cdot e^{-j\cdot\frac{2\pi\cdot k}{T}\cdot  t}dt\\
&=
\cfrac{1}{T}\cdot
\displaystyle\int_0^T f_{(t)}\cdot
\left(
cos\frac{2\pi\cdot k}{T}\cdot t + j\cdot sin\frac{2\pi\cdot k}{T}\cdot t
\right)dt
\end{array}
$$

が成り立つと分かってもらえると思います。

付録2. 離散フーリエ変換

離散フーリエ変換（DFT：Discrete Fourier Transform）は、変換対象波形がディジタル（離散時間）データなことを前提としたフーリエ変換です。

基本的に、DFTは通常のフーリエ変換と同じです。

DFTが通常のフーリエ変換と異なる点は、

• 変換対象波形$${\bm{f_{(n)}}}$$はサンプル時間$${\bm{t_s}}$$毎に変化する

• Nサンプルで１周期として扱う
  （周期を$${T}$$とすると$${N=T/t_s}$$、$${T}$$は$${t_s}$$の整数倍にする必要あり）

• フーリエ変換の周波数範囲が有限
  （最大周波数：$${0.5/{t_s}_{[Hz]}}$$）

• 積分計算の代わりに積算計算を使う

です。

解析対象とするディジタルデータ列$${f_{(n)}}$$は、下式のように表せます。
（$${k,n}$$は整数、$${n}$$はサンプリング番号でもある）

$$
\begin{array}{ll}
f_{(n)}&=\displaystyle\sum_{k=0}^{N-1}
\left(
F_{(k)}\cdot e^{j\cdot\frac{2\pi\cdot k}{T}\cdot (t_s\cdot n)}
\right)\\
&=\displaystyle\sum_{k=0}^{N-1}
\left(
F_{(k)}\cdot e^{j\cdot\frac{2\pi\cdot k}{N}\cdot n}
\right)\\
&=\displaystyle\sum_{k=0}^{N-1}
\left(
F_{(k)}\cdot
\left(
cos\frac{2\pi\cdot k}{N}\cdot n + j\cdot sin\frac{2\pi\cdot k}{N}\cdot n
\right)
\right)\\
&=\displaystyle\sum_{k=0}^{N-1}
\left(
F_{(k)}\cdot
cos\frac{2\pi\cdot k}{N}\cdot n
\right)
+j\cdot\displaystyle\sum_{k=0}^{N-1}
\left(
F_{(k)}\cdot
sin\frac{2\pi\cdot k}{N}\cdot n
\right)
\end{array}
$$

DFTでの$${F_{(k)}}$$は、通常のフーリエ変換における$${c_{(k)}}$$に相当します。
よって、$${c_{(k)}}$$を求める計算式中の積分を積算に置き換えたものが、$${F_{(k)}}$$を求める計算式になります。

$$
\begin{array}{ll}
F_{(k)}&=\cfrac{1}{T}\cdot\displaystyle\sum_{n=0}^{N-1}
\left(
f_{(n)}\cdot e^{-j\cdot\frac{2\pi\cdot k}{T}\cdot (t_s\cdot n)}
\right)\cdot t_s\\
&=\cfrac{1}{N}\cdot\displaystyle\sum_{n=0}^{N-1}
\left(
f_{(n)}\cdot e^{-j\cdot\frac{2\pi\cdot k}{N}\cdot n}
\right)\\
&=\cfrac{1}{N}\cdot\displaystyle\sum_{n=0}^{N-1}
\left(
f_{(n)}\cdot \left(cos\frac{2\pi\cdot k}{N}\cdot n-j\cdot sin\frac{2\pi\cdot k}{N}\cdot n\right)
\right)\\
&=\cfrac{1}{N}\cdot\displaystyle\sum_{n=0}^{N-1}
\left(
f_{(n)}\cdot cos\frac{2\pi\cdot k}{N}\cdot n\right)
-j\cdot\cfrac{1}{N}\cdot\displaystyle\sum_{n=0}^{N-1}
\left(
f_{(n)}\cdot sin\frac{2\pi\cdot k}{N}\cdot n
\right)
\end{array}
$$

後は、$${k=0,1,\cdots,N-1}$$と$${k}$$を変えながら、$${F_{(k)}}$$を求めていけばDFTは完了です。

スペクトルグラフを描く場合、横軸は$${\omega=\frac{2\pi\cdot k}{N\cdot t_s}_{[rad/s]}}$$または$${f=\frac{k}{N\cdot t_s}_{[Hz]}}$$とすれば良いです。
（$${k=0,1,\cdots,N-1}$$）

なお、積分を積算に置き換えた結果、DFTは行列形式で表現することが可能になっています。

具体的な数式は以下です。

$$
\left[
\begin{array}{c}
F_{(0)}\\
F_{(1)}\\
F_{(2)}\\
\vdots\\
F_{(N-1)}
\end{array}
\right]
=
\cfrac{1}{N}
\left[
\begin{array}{ccccc}
1&1&1&\cdots &1\\
1&e^{-j\cdot \frac{2\pi}{N}}&e^{-j\cdot \frac{4\pi}{N}}&\cdots &e^{-j\cdot \frac{2\pi\cdot(N-1)}{N}}\\
1&e^{-j\cdot \frac{4\pi}{N}}&e^{-j\cdot \frac{8\pi}{N}}&\cdots &e^{-j\cdot \frac{4\pi\cdot(N-1)}{N}}\\
\vdots & \vdots & \vdots & &\vdots\\
1&e^{-j\cdot \frac{2\pi\cdot(N-1)}{N}}&e^{-j\cdot \frac{4\pi\cdot(N-1)}{N}}&\cdots&e^{-j\cdot \frac{2\pi\cdot(N-1)\cdot(N-1)}{N}}
\end{array}
\right]
\left[
\begin{array}{c}
f_{(0)}\\
f_{(1)}\\
f_{(2)}\\
\vdots\\
f_{(N-1)}
\end{array}
\right]
$$

右辺の左側行列の中身は、全て$${e^{-j\cdot\frac{2\pi}{N}}}$$のべき乗になります。
（$${1}$$は$${e^{-j\cdot\frac{2\pi}{N}}}$$の$${0}$$乗）

また、$${e}$$のべき乗には、以下のような数学的性質があります。
（オイラーの公式から導き出せます）

• $${\bm{e^{j\cdot(\theta+\alpha\cdot\pi)}=e^{j\cdot \theta}}}$$（$${\alpha}$$は偶数）

• $${\bm{e^{j\cdot(\theta+\beta\pi)}=e^{-j\cdot \theta}=-e^{j\cdot\theta}}}$$（$${\beta}$$は奇数）

これらを利用して、少ない計算量でDFTを実行する方法を
「高速フーリエ変換」（FFT：Fast Fourier Transform）
と呼びます。

DFTとFFTとで得られる結果は同じなので、FFTの詳細説明は省略します。

ただ、データ数が膨大だったり、限られた計算リソースでフーリエ変換を実行する必要がある場合は、FFTを使うと良いです。

なお、複素数計算に対応しないプログラミング言語でDFTを実装する場合、実数部分と虚数部分をぞれぞれ個別に計算させれば問題なく実装できます。
（変数管理が面倒なのは仕方ないです）

付録3. 逐次最小二乗法

モータイナーシャ$${J_m}$$や粘性摩擦係数$${D_m}$$などの同定対象パラメータは、時々刻々変化する場合があります。

そうしたとき、最小二乗法に用いるデータを逐次追加し、最小二乗法の計算を逐次実行していくのが効果的です。

ただ、毎度毎度、4章の計算を丸ごとやるのは、コンピュータにとってさえ大変なことです。
（時間が経つほど計算量が増えてしまうこともありますし）

前回の最小二乗法の計算結果に対し、新たに追加されたデータによる計算結果を付け足せば、少ない計算量で逐次的な最小二乗法の計算ができます。

これが、逐次最小二乗法です。 

実際には、ある程度データが溜まったら通常の最小二乗法の計算をし、その次の回から逐次計算部分を付け足していくという使い方をします。

それでは、具体的なやり方を説明します。
まず、4章の最小二乗法の計算式を再掲します。

$$
\left[
\begin{array}{c}
J_m\\
D_m
\end{array}
\right]
=
\left[
\begin{array}{cc}
\displaystyle\sum_{k=1}^n\left({\ddot{\theta}_{(k)}}^2\right)
&
\displaystyle\sum_{k=1}^n\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)
\\
\displaystyle\sum_{k=1}^n\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)
&
\displaystyle\sum_{k=1}^n\left({\dot{\theta}_{(k)}}^2\right)
\end{array}
\right]^{-1}
\left[
\begin{array}{c}
\displaystyle\sum_{k=1}^n\left(\ddot{\theta}_{(k)}\cdot\tau_{(k)}\right)
\\
\displaystyle\sum_{k=1}^n\left(\dot{\theta}_{(k)}\cdot\tau_{(k)}\right)
\end{array}
\right]
$$

これを最新のn回めのデータに関する部分と、そうでない部分に分けます。

$$
\begin{array}{ll}
\left[
\begin{array}{c}
J_m\\
D_m
\end{array}
\right]
&=
\left(
\left[
\begin{array}{cc}
\displaystyle\sum_{k=1}^{n-1}\left({\ddot{\theta}_{(k)}}^2\right)+{\ddot{\theta}_{(k)}}^2
&
\displaystyle\sum_{k=1}^{n-1}\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)+\ddot{\theta}_{(n)}\cdot\dot{\theta}_{(n)}
\\
\displaystyle\sum_{k=1}^{n-1}\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)+\ddot{\theta}_{(n)}\cdot\dot{\theta}_{(n)}
&
\displaystyle\sum_{k=1}^{n-1}\left({\dot{\theta}_{(k)}}^2\right)+{\dot{\theta}_{(n)}}^2
\end{array}
\right]
\right)^{-1}
\left[
\begin{array}{c}
\displaystyle\sum_{k=1}^{n-1}\left(\ddot{\theta}_{(k)}\cdot\tau_{(k)}\right)+\ddot{\theta}_{(n)}\cdot\tau_{(n)}
\\
\displaystyle\sum_{k=1}^{n-1}\left(\dot{\theta}_{(k)}\cdot\tau_{(k)}\right)+\dot{\theta}_{(n)}\cdot\tau_{(n)}
\end{array}
\right]\\
&=
\left(
\left[
\begin{array}{cc}
\displaystyle\sum_{k=1}^{n-1}\left({\ddot{\theta}_{(k)}}^2\right)
&
\displaystyle\sum_{k=1}^{n-1}\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)
\\
\displaystyle\sum_{k=1}^{n-1}\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)
&
\displaystyle\sum_{k=1}^{n-1}\left({\dot{\theta}_{(k)}}^2\right)
\end{array}
\right]
+
\left[
\begin{array}{cc}
{\ddot{\theta}_{(n)}}^2&\ddot{\theta}_{(n)}\cdot\dot{\theta}_{(n)}\\
\ddot{\theta}_{(n)}\cdot\dot{\theta}_{(n)}&{\dot{\theta}_{(n)}}^2
\end{array}
\right]
\right)^{-1}
\left[
\begin{array}{c}
\displaystyle\sum_{k=1}^{n-1}\left(\ddot{\theta}_{(k)}\cdot\tau_{(k)}\right)+\ddot{\theta}_{(n)}\cdot\tau_{(n)}
\\
\displaystyle\sum_{k=1}^{n-1}\left(\dot{\theta}_{(k)}\cdot\tau_{(k)}\right)+\dot{\theta}_{(n)}\cdot\tau_{(n)}
\end{array}
\right]\\
&=
\left(
\left[
\begin{array}{cc}
\displaystyle\sum_{k=1}^{n-1}\left({\ddot{\theta}_{(k)}}^2\right)
&
\displaystyle\sum_{k=1}^{n-1}\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)
\\
\displaystyle\sum_{k=1}^{n-1}\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)
&
\displaystyle\sum_{k=1}^{n-1}\left({\dot{\theta}_{(k)}}^2\right)
\end{array}
\right]
+
\left[
\begin{array}{c}
\ddot{\theta}_{(n)}\\
\dot{\theta}_{(n)}
\end{array}
\right]
\cdot 1 \cdot
\left[
\begin{array}{cc}
\ddot{\theta}_{(n)}&\dot{\theta}_{(n)}
\end{array}
\right]
\right)^{-1}
\left[
\begin{array}{c}
\displaystyle\sum_{k=1}^{n-1}\left(\ddot{\theta}_{(k)}\cdot\tau_{(k)}\right)+\ddot{\theta}_{(n)}\cdot\tau_{(n)}
\\
\displaystyle\sum_{k=1}^{n-1}\left(\dot{\theta}_{(k)}\cdot\tau_{(k)}\right)+\dot{\theta}_{(n)}\cdot\tau_{(n)}
\end{array}
\right]
\end{array}
$$

ここで、逆行列の補助定理$${(A+BDC)^{-1}=A^{-1}B(D^{-1}+CA^{-1}B)^{-1}CA^{-1}}$$を適用することを考えます。
（中学・高校では教わらないようなので、知らないのが普通な定理だと思います）

$$
A=
\left[
\begin{array}{cc}
\displaystyle\sum_{k=1}^{n-1}\left({\ddot{\theta}_{(k)}}^2\right)
&
\displaystyle\sum_{k=1}^{n-1}\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)
\\
\displaystyle\sum_{k=1}^{n-1}\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)
&
\displaystyle\sum_{k=1}^{n-1}\left({\dot{\theta}_{(k)}}^2\right)
\end{array}
\right]
$$

$$
B=
\left[
\begin{array}{c}
\ddot{\theta}_{(n)}\\
\dot{\theta}_{(n)}
\end{array}
\right]
$$

$$
C=
\left[
\begin{array}{cc}
\ddot{\theta}_{(n)}&
\dot{\theta}_{(n)}
\end{array}
\right]
$$

$$
D=1
$$

$$
\begin{array}{l}
\left(
\left[
\begin{array}{cc}
\displaystyle\sum_{k=1}^{n-1}\left({\ddot{\theta}_{(k)}}^2\right)
&
\displaystyle\sum_{k=1}^{n-1}\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)
\\
\displaystyle\sum_{k=1}^{n-1}\left(\ddot{\theta}_{(k)}\cdot\dot{\theta}_{(k)}\right)
&
\displaystyle\sum_{k=1}^{n-1}\left({\dot{\theta}_{(k)}}^2\right)
\end{array}
\right]
+
\left[
\begin{array}{c}
\ddot{\theta}_{(n)}\\
\dot{\theta}_{(n)}
\end{array}
\right]
\cdot 1 \cdot
\left[
\begin{array}{cc}
\ddot{\theta}_{(n)}&\dot{\theta}_{(n)}
\end{array}
\right]
\right)^{-1}\\
=A^{-1}+A^{-1}
\left[
\begin{array}{c}
\ddot{\theta}_{(n)}\\
\dot{\theta}_{(n)}
\end{array}
\right]
\left(
1+
\left[
\begin{array}{cc}
\ddot{\theta}_{(n)}&\dot{\theta}_{(n)}
\end{array}
\right]A^{-1}
\left[
\begin{array}{c}
\ddot{\theta}_{(n)}\\ \dot{\theta}_{(n)}
\end{array}
\right]
\right)^{-1}
\left[
\begin{array}{cc}
\ddot{\theta}_{(n)}&\dot{\theta}_{(n)}
\end{array}
\right] A^{-1}
\end{array}
$$

逆行列の補助定理の良い点は、逆行列計算をすることなく、結果が得られるところです。

ここでは、未知パラメータが２つしかないので、あまり有難みを感じないかもしれません。

しかし、未知パラメータが増えて行列サイズが大きくなると、有難みが分かると思います。

なお、$${A^{-1}}$$の初期値には、逐次最小二乗法に入る前の計算結果を使えます。
 
$${n+1}$$回めでは、$${n}$$回めで逆行列の補助定理で計算した$${(A+BDC)^{-1}}$$を$${A^{-1}}$$として使えます。

こうして、新たなデータを付け足していく形でのパラメータ同定が可能になります。

しかし、このままでは逐次計算しても狙った効果は得にくいです。

なぜなら、パラメータ同定結果が、開始時から現時刻までの全データに対する平均値的な数値となるためです。

パラメータが時々刻々変化するなら、古いデータがパラメータ同定結果に与える影響を小さくする工夫が必要です。

具体的な方法としては、

• 忘却係数λを導入する（0<λ<1）

• 新しいデータを1つ付け足したら、玉突きのように古いデータを差し引く

の２とおりが挙げられます。
（他にもあるかもしれませんが）

今回は、忘却係数$${\lambda}$$の導入のやり方のみを紹介しようと思います。

これについては、数式よりコードで説明した方が分かりやすいと思います。
行列計算なので、MATLABスクリプト形式で説明します。

Lambda=0.99;   // 忘却係数λ（0<λ<1）

// 以下、パラメータ同定のメイン処理部分（サンプル毎に繰り返す）
// Trqはトルク、Accは角加速度、Velは角速度
Ainv=Ainv/Lambda;   // 過去データにλ掛け（逆行列なのでλの逆数を掛ける）
B1=Lambda*B1+Acc*Trq;   // 過去データにλ掛けして、現在データ分を加算
B2=Lambda*B2+Vel*Trq;   // 過去データにλ掛けして、現在データ分を加算

B=[B1; B2];
BT=B’;   // 行列Bの転置行列

Ainv=Ainv-Ainv＊B/(1+BT*Ainv*B)*BT*Ainv;   // 現在データを使ったAinv計算（逆行列の補助定理）

Result=Ainv*B;   // 最小二乗法でのパラメータ同定結果（Resultの１行目がJm、２行目がDm）

こうすると、繰り返し計算１回毎に過去データに$${\lambda}$$を１回掛けられます。
 
仮に、$${\lambda=0.99}$$に設定したとして、１サンプル前のデータには$${0.99}$$がそのまま掛かるだけです。

しかし、２サンプル前のデータには$${0.99^2\fallingdotseq0.98}$$が、１０0サンプル前のデータなら$${0.99^{100}\fallingdotseq0.37}$$が掛かるのと同じになります。

このように、忘却係数$${\lambda}$$を導入すれば、過去データほど数値が小さくなって、パラメータ同定結果への影響度合いが低くなっていきます。 

下図は、忘却係数なし/ありのパラメータ同定シミュレーション結果です。
（動作中にモータに連結されたボールねじテーブル上におもりが載せられるなどして、モータからみたイナーシャが増加する想定のシミュレーション）

忘却係数なしだと、これまでの全データの平均値的な結果が出ます。
上図の場合、全データの半分はイナーシャが増える前のデータなため、同定値が真値に行き着いていません。

一方、忘却係数$${\lambda}$$を導入すると、過去データのパラメータ同定結果への影響度合いが小さくなって、短い時間（＝少ないサンプル数）で同定値が真値に近づきます。

この結果から、パラメータが時々刻々変化する場合は、忘却係数$${\lambda}$$を導入した方が良いと分かります。 

ただ、短時間でパラメータ同定するために$${\lambda}$$を小さくしすぎると、ほぼ現在データのみからパラメータ同定をするのに近い状態になります。

そうすると、各種ノイズの影響を受けやすく、同定精度が低下する可能性が高くなります。
 
このトレードオフを考慮して、$${\lambda}$$は$${0.9～0.999}$$の範囲に設定するケースが多いようです。
（データのサンプルレートにもよりますが）

前回：寄り道編（３）「状態方程式」はこちら
次回：寄り道編（５）「システム同定時の入力信号」はこちら