目次
第1 章 序論 3 1.1 研究背景・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・3 1.2 研究目的・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・5 第2 章 非線形システム同定法 6 2.1 従来のシステム同定理論・ ・・・・・・・・・・・・・・・・・・・・・6 2.1.1 システム同定とは・・・・・・・・・・・・・・・・・・・・・6 2.1.2 ARX モデル・・・ ・・・・・・・・・・・・・・・・・・・・ 7 2.1.3 最小二乗法に基づくシステム同定法 ・・・・・・・・・・・・ 10 2.2 問題の定式化 ・・・・・・・・・・・・・・・・・・・・・・・・・・・13 2.3 準ニュートン法 ・・・・・・・・・・・・・・・・・・・・・・・・・・15 2.3.1 直線探索法の概要 ・・・・・・・・・・・・・・・・・・・・16 2.3.2 Armijo の基準 ・・・・・・・・・・・・・・・・・・・・・17 2.3.3 降下方向 ・・・・・・・・・・・・・・・・・・・・・・・・17 2.3.4 セカント条件 ・・・・・・・・・・・・・・・・・・・・・・19 2.3.5 BFGS 公式 ・・・・・・・・・・・・・・・・・・・・・・21 第3 章 非線形システムの同定 22 3.1 飽和特性を有する非線形センサの場合 ・・・・・・・・・・・・・・・・22 3.1.1 問題の定式化・・・・・・・・・・・・・・・・・・・・・・・23 3.1.2 シミュレーション結果・・・・・・・・・・・・・・・・・・・24 3.2 分解能を有する非線形センサの場合 ・・・・・・・・・・・・・・・・・28 3.2.1 問題の定式化・・・・・・・・・・・・・・・・・・・・・・・28 3.2.2 分解能を考慮した場合・・・・・・・・・・・・・・・・・・・30 3.2.2 シミュレーション結果・・・・・・・・・・・・・・・・・・・31 第4 章 制御対象の概要と同定実験 36 4.1 SPIDER 駆動精密ステージの概要・・・・・・・・・・・・・・・・・・36 4.2 制御対象の事前同定実験 ・・・・・・・・・・・・・・・・・・・・・・39 4.3 従来法との比較・・・・・・・・・・・・・・・・・・・・・・・・・・42 第5 章 精密ステージ制御への応用 45 5.1 精密ステージに及ぼす摩擦力 ・・・・・・・・・・・・・・・・・・・・455.2 シミュレーション結果・・・・・・・・・・・・・・・・・・・・・・・48 5.3 実験結果・・・・・・・・・・・・・・・・・・・・・・・・・・・・・49 第6 章 まとめ 51 参考文献 52 発表論文 53 謝辞 53
第１章 序論
1.1 研究背景
線形システム同定に関する古典的な理論研究は、Ljung の著書によってほぼ体系づけられ、 また、システム同定理論を実行するソフトウェアも用意されている。 しかしながら、世の中に実在する自動車エンジン、超精密ステージ、ロボットなど、全 てのシステムはなんらかの非線形特性を有している。したがって、このようなシステムに 対しては、制御ばかりでなくシステムの予測、解析、故障診断、さらに高性能化などを行 う際に、非線形特性を考慮したモデル化手法及び制御法が必要となる。そのため、システ ムのモデリングにあたって非線形システム同定は近年注目を集めており、1990 年代以降に シグモイド型ニューラルネットワーク、RBF（Radial Basis Function）ネットワーク、フ ァジィモデルなどによるアプローチが精力的に研究されているが、非線形と呼ぶ対象があ まりに広大なため、まだ体系的な非線形同定理論は完備されていない。 また、従来のシステム同定理論においては入出力信号を測定するセンサは線形センサが 用いられ,観測雑音としてたとえば白色雑音が混入する場合を仮定することがほとんどで あった。しかしながら,産業界ではコスト面の理由から,つねに（高価な）線形センサを利 用できず,（安価な）静的非線形性をもつセンサを利用しなければならない状況も存在する。 非線形センサとして例えば、自動車エンジンの空燃比を検出する際に用いられる安価な ジルコニア O２センサやバイナリセンサ等が挙げられる。そのような測定値を用いて従来の 線形システム同定法を行った場合、同定精度は劣化してしまう。 また、センサの静的非線形性として量子化誤差がある。近年、ナノテクノロジーの開発 が活発に行われており、その中で機械システムをナノメートルオーダーで精密に制御する ことが必要とされている。例えば、半導体製造装置に使用される精密ステージには、ナノ レベルの精度が要求される。 機械システムをナノメートルオーダーで制御するためには、ナノレベルでの制御対象の 挙動を表すモデルが必要である。なぜなら、微小動作領域での運動のモデルと動作範囲の 広い運動のモデルは異なるからである。そのため、ナノレベルでの制御にはナノレベルで の運動に特化したモデルを同定する必要がある システム同定においては、一般に、アナログであるプラント入出力の信号を A/D 変換器 でディジタル化したデータを用いる。このとき、得られるディジタルデータは有限語調な ので量子化誤差が生じる。この量子化誤差は、特にナノレベルの運動モデルの運動モデル の同定精度に悪影響を及ぼす。なぜなら、A/D 変換器の分解能はそのフルレンジに対応して 決まるため、広範囲での動作が必要なプラントでは量子化単位がナノレベルより大きくな1.2 研究目的
同定対象は線形システムであると仮定するが、飽和、分解能などの静的非線形特性を有 するセンサにより入出力データを測定する場合、測定値は非線形誤差を含むことになる。 ここで従来のシステム同定理論では係数パラメータのみを推定するだけで非線形誤差は推 定されない。したがって同定の精度が劣化してしまう。そこで、本研究の目的として、シ ステム同定のための評価関数の中に、予測誤差の二乗和だけでなく、この非線形誤差に関 する項も組み込むことによって、係数パラメータだけでなく、非線形誤差も同時に推定す る手法を数値シミュレーションによりその有効性を検証する。また、精密ステージの位置 決め制御に対してこの手法を応用し,その有効性をシミュレーションと実験によって検証 する。 本論文の構成は以下となっている。 第２章では非線形システム同定法について述べる。ここでは、まず従来のシステム同定 理論について述べ、次に非線形誤差を考慮した非線形システム同定法を述べた後、最適化 問題を解くために用いた準ニュートン法のアルゴリズムについて述べる。第３章では飽和 特性を有する場合と分解能を有する場合に対して問題の定式化を行い、それぞれの場合に 対してシミュレーションを行い有効性を検証する。第 4 章では、精密ステージの静止摩擦、 分解能誤差を考慮した事前同定実験を行い、これを真値の同定モデルとした。この真値の 同定モデルを用いて非線形システム同定法により導出した同定モデル、従来法により導出 した同定モデルの比較を行う。第 5 章では精密ステージの位置決め制御への応用として本 研究に用いた実験装置の概要の説明をし、ステージの真値、従来のシステム同定法による 同定モデルとのシミュレーション結果の比較、そして実験結果を示す。第 6 章で本論文の まとめを述べる。第２章 非線形システム同定法
2.1 従来のシステム同定理論
2.1.1 システム同定とは
対象 とするシ ステムのふる まいを特 徴づけるモデ ルを構築 することをモ デリン グ (modeling)という。対象とするシステムや目的に応じて様々なモデリングの手法が存在する が、モデルとして数学モデルを用いるのか、図形モデルを用いるのかにより大分類できる。 数学モデルとは、システムのふるまいを代数方程式、微分方程式、あるいは論理式のよう な数学的な表現を用いたものである。制御系設計を行う場合には、一般に制御対象である システムの数学モデルが必要になる。このとき、数学モデルを構成する方法には次のよう なものがある。 • ホワイトボックスモデリング これは、モデルの構造が完全に既知である場合のモデリング方法である。ここでホ ワイトボックスとは、中が全て見える箱という意味で、システムを白い箱とみなし ている。この方法は物理モデリングとも呼ばれ、対象を支配する物理法則 (例えば 運動方程式、回路方程式、電磁界方程式、物質収支、熱収支、あるいは化学反応方 程式など) に基づいてモデリングを行う方法である。 • ブラックボックスモデリング これは、対象に関する物理的な情報や事前情報をいっさい利用しない方法である。 ここでは対象を黒い箱とみなしている。この方法が、システム同定と呼ばれるもの である。 ここで、ホワイトボックスモデリングは、対象が比較的小規模で、物理法則のみによっ てその動特性が記述できる場合には非常に有効であるが、対象が複雑になるにつれてその 適用は難しくなる。そのような状況で有効な方法が、ブラックボックスモデリング、即ち システム同定である。2.1.2 ARX モデル
ここでは、ARX モデルについて説明する。ARX モデルは式誤差モデルの一つである。式 誤差モデルとは、差分方程式)
(
)
1
(
)
(
k
a
1y
k
a
nay
k
n
ay
)
(
)
(
)
1
(
1u
k
b
u
k
n
e
k
b
nb
b
…(2.1.2.1) で表されるものである。ただし、u
は入力、y
は出力。a
、b
はパラメータである。また、)
(k
e
は外乱項で、差分方程式に直接、誤差として入っている。この式誤差モデルにおいて、 ARX モデルは外乱項e
(k
)
を白色雑音w
(k
)
と仮定したものである。すなわち、)
(
)
1
(
)
(
k
a
1y
k
a
nay
k
n
ay
)
(
)
(
)
1
(
1u
k
b
u
k
n
w
k
b
nb
b
…(2.1.2.2) と表される。 ここで、ARX モデルの離散時間 LTI システムの一般的な表現とその 1 段予測誤差につい て述べる。 まず、離散時間LTI システムとしては一般的に)
(
)
(
)
(
)
(
)
(
k
G
q
u
k
H
q
w
k
y
…(2.1.2.3) と表される。ここで、u
(k
)
は入力、y
(k
)
は出力、G
(q
)
は伝達関数、H
(q
)
は雑音モデル、)
(k
w
は白色雑音である。また、離散時間LTI システムの 1 段予測誤差は、上式で定義した 離散時間LTI モデルにおいて、時刻(
k
1
)
までに測定された入出力データに基づいた出力)
(k
y
の1 段予測誤差y
ˆ
(
k
|
)
として)
(
)
,
(
)
,
(
)
(
)]
,
(
1
[
)
|
(
ˆ
k
H
q
y
k
H
q
G
q
u
k
y
q
q
…(2.1.2.4) と与えられる。ただし、
はモデルを記述するパラメータから構成されるベクトルである。 以上のことから、ARX モデルの離散時間 LTI システムは、パラメータベクトルを]
,
,
,
,
,
[
a
1
a
nab
1
b
nb
…(2.1.2.5) とし、データベクトル(回帰ベクトル)をT b a
u
k
u
k
n
n
k
y
k
y
k
)
[
(
1
),
,
(
),
(
1
),
,
(
)]
(
…(2.1.2.6) と定義すれば、ARX モデルの出力y
(k
)
は)
(
)
(
)
(
k
k
w
k
y
T
…(2.1.2.7) となり、次式で定義される既約なシフトオペレータq
)
1
(
)
(
1
k
u
k
u
q
…(2.1.2.8) を用いた二つの多項式 na naq
a
q
a
q
A
1
11
)
(
…(2.1.2.9) nb nbq
b
q
b
q
B
1
1)
(
…(2.1.2.10) を用いて)
(
)
(
)
(
)
(
)
(
q
y
k
B
q
u
k
w
k
A
…(2.1.2.11) と書き直された後、伝達関数G
(q
)
と雑音モデルH
(q
)
を)
(
)
(
)
(
q
A
q
B
q
G
,)
(
1
)
(
q
A
q
H
…(2.1.2.12) とおくことで、ARX モデルの離散時間 LTI システムは)
(
)
(
1
)
(
)
(
)
(
)
(
w
k
q
A
k
u
q
A
q
B
k
y
…(2.1.2.13) とかけ、(
na
nb
)
個のパラメータを表現したパラメトリックモデルとなる。 さらに、ARX モデルの 1 段予測誤差は)
(
)
(
)
(
)
(
)]
(
1
[
)
|
(
k
A
q
y
k
B
q
u
k
k
y
T
…(2.1.2.14) と表され、この式からわかるようにARX モデルでは 1 段予測誤差は
に関して線形な関係 式となる。このことからARX モデルは線形回帰モデルといわれ、ARX モデルは線形なモ デルといえる。 また、図2.2.1 に ARX モデルのブロック線図を示しておく。 図2.1.1 ARX モデルのブロック線図2.1.3 最小二乗法
LS[Least-Squares method]法は、評価規範 (2.1.3.1) を最小化するθを求めることに対応しており、同定手順は以下のようになっている。 Step1: 推定誤差の計算 (2.1.3.2) Step2: θのパラメータの推定値
を 求める。 (2.1.3.3) なぜこうなるのかについて、詳しく説明する。 出力の一段階予測値y
(
k
)
がθに関して線形である線形回帰モデル (2.1.3.4) について考える。このときの予測誤差は、 (2.1.3.5) となり、このモデルに最小二乗法を適用する。すなわち、パラメータ推定のための評価規 範は次式のようになる。 (2.1.3.6) さらに、計算すると、 (2.1.3.7) また、
N k Nk
N
J
1 2)
(
1
)
(
)
1
(
)
(
)
(
)
(
k
y
k
Tk
k
N k N k Tk
y
k
k
k
k
1 1 1)
(
)
(
)
(
)
(
)
(
ˆ
)
(
)
(
k
k
y
T
)
(
)
(
)
(
k
y
k
T
k
N k T N k Ny
k
k
N
k
N
J
1 2 1 2)
(
)
(
1
)
(
1
)
(
N k T T T Ny
k
y
k
k
N
J
1 2)
(
)
(
2
)
(
1
)
(
とおくと、
J
N(
)
は以下のようになる。 (2.1.3.8) （1） 未知パラメータが１つの場合 このとき、J
N(
)
の変数はすべてスカラになり、θをx
とおくと、 (2.1.3.9) となり、x
に関する二次方程式となる。ここで、R
(N
)
＞0 であれば、J
N(
)
は下に凸にな るので、J
N(
)
の最小値は、以下のように求めることが可能である。 (2.1.3.10) よって、 (2.1.3.11) のとき、J
N(
)
は最小値をとる。 （2） 未知パラメータが２つ以上の場合 (1)と同様であるが、下に凸になるための条件が、R
(N
)
＞0 から、行列R
(N
)
が正定値行 列（固有値が全て正である行列）であるとういう条件に変わる。 これは、正定値行列であれば、逆行列が存在するので、逆行列を用いてパラメータの推定 値を求めることが出来るからである。)
(
NJ
をθで微分したものを０とすると、下のようになり、 (2.1.3.12)
N k N k N kk
k
N
N
R
k
y
k
N
N
f
k
y
N
N
C
1 1 1 2)
(
)
(
1
)
(
)
(
)
(
1
)
(
)
(
1
)
(
(
)
2
)
(
N
f
N
R
N
C
J
N
T
T
2)
(
2
)
(
N
f
N
x
R
N
x
C
J
N
0
)
(
2
)
(
2
f
N
R
N
J
dx
d
N)
(
)
(
N
R
N
f
x
(
)
0
2
J
f
N
R
N
d
d
N)
(
NJ
を最小とするパラメータの推定値は、 (2.1.3.13) となる。 LS 法は、一括処理最小二乗法、あるいは、オフライン最小二乗法と呼ばれる。 オフラインの意味は、同定と実験を別々に行うということである。)
(
)
(
)
(
N
R
N
1f
N
2.2 問題の定式化
本論文では、図 2.1.1 にブロック線図を示したような非線形特性をもつセンサを用いたシ ステム同定問題について考える。ここで考える非線形特性は、飽和、不感帯などである。 線形システムの入力側に（静的）非線形性を有する場合はHammerstein モデルと呼ばれ、 非線形特性をもつアクチュエーターによって線形システムを制御する際に有用であること が知られている。それに対して、ここで考えている同定問題では、対象に飽和が存在しな いそのままの入力信号が印加されているがその入力信号は正確には測定できず、飽和特性 をもつセンサによって測定される状況を想定している。出力信号も非線形特性をもつセン サを用いて測定されるが、これは通常用いられるWeiner モデルと同じ形式をとっている。Input signal u(t) Output signal y(t)
u(k) y(k)
Nonlinear sensor
Measured input signal Measured output signal
)
(k
u
)
(k
y
Linear system
Sampling
Sampling
Static nonlinearity Static nonlinearity 図 2.2.1 非線形センサによる線形システムの同定問題例えば、コントローラの特性が既知でないときの閉ループ同定実験の場合や、入力信号 が他のシステムの出力として与えられている場合などは、開ループ同定実験の場合と異な り、入力信号自身をセンサによって測定しなければならず、このような問題設定が必要と なる。 図 2.1.1 における同定対象は線形時不変システムであると仮定する。すると、同定対象の 離散時間における入出力関係は、差分方程式
)
(
)
1
(
)
(
k
a
1y
k
a
y
k
n
y
n
b
1u
(
k
1
)
b
nu
(
k
n
)
(2.2.1) で記述される。ただし、u(k)は入力信号、y(k)は出力信号である。本論文では、システムの 次数ｎは既知であると仮定する。 また、本論文では、センサの非線形特性として飽和、不感帯などを考えるが、非線形特 性の種類とその特性を与えるパラメータは既知であると仮定する。従って、本論文では、 何らかの方法によって、センサの非線形特性があらかじめ同定できている状況を想定する。 もちろん、センサの非線形特性を入出力データから同時に同定する方法を開発することは、 今後の重要な課題である。2.3 準ニュートン法
この節では、直線探索法の基本的な概念を述べ、それを基に準ニュートン法の概要を述 べる。準ニュートン法の全体的手順は図 2.3.1 のように構成される。また、イメージ図を図 2.3.2 に示す。 ②探索方向の計算 ① 初期化 ④直線探索 ⑤Hkの更新 ③終了条件 停止 図 2.3.1 準ニュートン法による探索 図 2.3.2 イメージ図① 初期化：初期点 x0、初期行列 H0(通常単位行列) を与え、k = 0 とする ② 探索方向の計算：探索方向 dk = −Hk∇f(xk) を求める ③ 終了条件：停止条件が満たされていれば、xk を解として出力し停止。そうでなけれ ば④へ ④ 直線探索： Armijo の基準を用いた直線探索によりステップ幅 αk を求める xk+1 = xk + αkdk とおく ⑤ Hｋの更新：sk = xk+1 − xk, yk = ∇f(xk+1) − ∇f(xk) をそれぞれ求める BFGS 公式 (9) を用いて、行列 Hk を更新する k = k + 1 とおいて ②へ
2.3.1 直線探索法の概要
制約無し最適化問題を解くためのアルゴリズムを構築する際には、初期点 x0 を与え、 系列 {xk} を生成する。そしてこれ以上の更新が起こらない、または十分な精度の近似が出 来ている際に、終了する。基本的に xk→ xk+1 と更新する際に、f(xk), ∇f(xk) 等の情報を 用いる。現在の点 xk から次の点 xk+1 へ移動する際に、直線探索法と信頼領域法の２つ の基本的なアプローチがある。以下では、本論文で用いる直線探索法について記述する。 直線探索法では、探索方向 pk を決定し、この向きに現在いる点 xk から進む。xk から 方向ベクトルの向きへ移動する距離 α を、次の一次の最小化問題を解くことで求めるこ とになる。 min ( ) 0 f xk
kpk （2.3.1.1） 式 （2.3.1.1） を正確に解くことで、最も有用なステップ幅 αk を見つけることができ るが、この大域的最小解を特定するのは困難である。そのため、典型的な直線探索アルゴ リズムでは、候補となる α の系列を生成し、ある条件が満たされた時に、値の１つを受 理して終了する。その条件としては、Armijo の基準、Wolfe の条件、Goldstein の条件等 があるが、以下では最も基本的でかつ本論文で使用した Armijo の基準について紹介する。2.3.2 Armijo の基準
k
の満たすべき条件として、十分に目的関数の値を減少させるような値をとることが望 まれる。Armijo の基準は、
kが次の式を満たす時に受理する手法である。 T k k k k k k kp
f
x
c
f
x
p
x
f
(
)
(
)
1
(
)
（2.3.2.1） ここでc1 はc1(0,1)となる定数である。式 （2.3.2.1）の左辺・右辺をそれぞれ
kの関数 と見て
(
),
l
(
)
とおく。この時、関数l
(
)
は、傾きがc
1
f
(
x
k)
Tp
k となり、負の値を とる。十分に降下する条件として、
(
)
l
(
)
を満たす時に、
を受理するようなアルゴ リズムとなる。 ただし、0 に近い全ての α を式（2.3.2.1）は満たすことになるので、このアルゴリズ ムは合理的に進行しない可能性がある。α が小さくなり過ぎるのを防ぐために、Wolfe の 条件や Goldstein の条件などがある。本論文におけるシミュレーションでは、Armijo の 基準を使用して実装している。2.3.3 降下方向
探索方向の選択には、様々な手法があるが、探索方向を決定する際に、目的関数を何ら かの意味で近似したモデル関数を局所的に最小化することが考えられる。モデル関数とし て、1 次のモデル (linear model)p
x
f
x
f
l
p
x
f
(
k
k k)
(
)
(
k)
(
k)
T と 2 次モデル (quadratic model)
f
x
k kp
kq
f
x
kf
x
k Tp
p
Tf
(
x
k)
p
2
1
)
(
)
(
)
(
)
(
2 を用いるのが一般的である。1 次モデルに基づいた方法として最急降下法が、また 2 次モ デルに基づいた方法としてニュートン法が挙げられる。探索方向 を選択する際には、目的 関数の値が単調に減少する、すなわちf
(
x
k
1
)
f
(
x
k)
となるような系列が生成されるこ とが望ましい。このような系列の大域的収束性を保証するためには、k 回目の近似解 での 探索方向（p
k方向）での方向微係数が負になることが望ましい。f(x) が微分可能の場合に は、このことは
f
(
x
k)
Tp
0
（2.3.3.1） と書ける 従って、式（2.3.3.1）が成り立つような探索方向p
kを見つけることが重要になる。 このような探索方向を、f
(x
)
の降下方向 (decent direction) と呼ぶ。最急降下法では、探 索方向を、p
k
f
(
x
k)
とする。この時、 f(xk)T pf(xk)2 0 となるので、これは式（2.3.3.1）を満たす。また、ニュートン法では、探索方向を、)
(
k kf
x
p
とする。ヘッセ行列
2f
(
x
k)
が正定値行列となる時、0
)
(
)
(
)
(
)
(
2 1
k k k k T kp
f
x
f
x
f
x
x
f
となり、式（2.3.3.1） を満たす。ニュートン法は、局所的に二次収束することが知られ ているが、
2f
(
x
k)
が正定値行列でない時、式（2.3.3.1）、（2.3.2.1）を満たさないので、 ニュートン方向は定義されない。この点を補うために、修正ニュートン法, 信頼領域法, 準ニュートン法などが提案されている。 ニュートン法では前述したように、ヘッセ行列が正定値行列である時にしか探索方向が 定義されない。その短所を補うべく代替案として、準ニュートン法が提案された。準ニュ ートン法は、最急降下法の「大域的収束性」とニュートン法の「局所的に速い収束性」と いう長所をそれぞれ合わせ持っており、制約無し問題に対する数値解法の中では、現在最 も有力な手法の 1 つとなっている。 準ニュートン法では、ニュートン方向におけるヘッセ行列の代わりに、その近似行列と して、正定値対称行列
B
kを用いる。探索方向は、以下のようになる。p
k
B
k 1
f
(
x
k)
（2.3.3.2） が正定値行列である、という条件の元で、0
)
(
)
(
)
(
1
k k k T kp
f
x
B
f
x
x
f
となるので、探索方向p
kはf
(x
)
の降下方向となる。2.3.4 セカント条件
行列B
kよって、
2f
(
x
k)
を近似したいので、ヘッセ行列の情報を取り込んだなんらかの 基準をB
kが満足することが必要となる。そこで、目的関数f
(
x
k)
の一次導関数のテイラー 展開より、
f
x
tp
f
x
p
p
dt
p
p
x
f
x
f
dt
p
tp
x
f
x
f
p
x
f
1 0 2 2 2 1 0 2)
(
)
(
)
(
)
(
)
(
)
(
)
(
x
x
k 、p
x
k1
x
k とすると、f
(x
)
が連続であることより、)
(
)
)(
1
(
)
(
)
1
(
x
kf
x
k 2f
x
kx
k 1x
ko
x
k 1x
kf
と書け、o
(
x
k1
x
k)
0
とすると、
f
(
x
k
1
)
f
(
x
k)
2f
(
x
k
1
)(
x
k1
x
k)
となり、このときs
k
x
k1
x
k、y
k
f
(
x
k1)
f
(
x
k)
とおいた時に、近似行列B
k1 がy
k
B
k1s
k （2.3.4.1） を満たすことを条件とすることで、行列 は、ヘッセ行列
2f
(
x
k)
を近似できる。式 （4.1.4） はセカント条件と呼ばれている。従って、セカント条件（4.1.4）を満たすようにB
kを更新 してB
k1を作ることが考えられる。これが準ニュートン法の基本的な考え方となる。行列 1 kB
の成分が未知数であることから、セカント条件は不定方程式になる。そのためセカント 条件を満たすB
k1は無数に存在し、様々な更新公式が考えられている。2.3.5 BFGS 公式
セカント条件（4.1.4）を満たす更新公式の代表的なものとして、Davidon(1959) と
Fletcher and Powell(1963) に よ っ て 提 案 さ れ た DFP 公 式 と Broyden(1970),Fletcher(1970) Goldfarb(1970),Shanno(1970) に よ っ て 提 案 さ れ た BFGS 公式の２つがある。 DFP 公式 k T k T k k k T k k k T k k T k k k k k k
y
s
y
y
y
s
s
B
s
y
s
y
s
B
B
B
1
(
1
)
BFGS 公式 k T k T k k k k T k T k k k k k ky
s
y
y
s
B
s
s
B
s
B
B
B
1
(
)
実践的な準ニュートン法のアルゴリズムでは、ヘッセ行列
2f
(
x
k)
を行列B
kで近似するよ りも、ヘッセ行列の逆行列
2f
(
x
k)
1を行列H
kで近似する方が一般的である。理由は更 新公式を利用して得られたB
kの更に逆行列B
k1を求める、という手間を省くためである。 このヘッセ行列の逆行列の近似行列 の BFGS 公式は以下のようになる。 k T k T k k k T k k k T k k T k T k k k T k k k k ky
s
s
s
y
s
y
H
y
y
s
y
H
s
s
y
H
H
H
1
(
)
(
1
)
（2.3.5.1） 今日、理論的にも実用的にも最も有効な公式は BFGS 公式であることが、広く知られてい る。第
3 章 非線形システムの同定
3.1 飽和特性を有する非線形センサへの応用
3.1.1 問題の定式化
図 2.1.1 における非線形特性として、図 3.1.1 に示した飽和特性の場合を考える。図 3.1.1 において、非線形センサによって測定された出力データをy
(k
)
、出力における飽和量を)
(k
とすると、真の出力信号の推定値y
ˆ k
(
)
は次式のように表わされる。
L
k
y
if
k
k
y
L
k
y
if
k
y
L
k
y
if
k
k
y
k
y
)
(
),
(
)
(
)
(
),
(
)
(
),
(
)
(
)
(
ˆ
ここで出力の線形領域は
L
とし、L の値は既知であると仮定する。これより。出力信号 の飽和量に対する制約条件
(
k
)
0
（3.1.1.1） が得られる。この制約条件より、飽和量を推定する方向（+あるいはー）を限定することが できる。 いま、時刻ｋにおける式誤差e(k)は、次式のように記述できる。)
(
)
1
(
)
(
ˆ
)
1
(
ˆ
)
(
ˆ
)
(
1 1n
k
u
b
k
u
b
n
k
y
a
k
y
a
k
y
k
e
n n
（3.1.1.2） ここで、式誤差を計算するとき、真の出力データが利用できないことに注意する。
Output
δ
(k)L
－
L
L input
―L
図 3.1.1 飽和特性3.1.2 飽和特性に対する解法
飽和特性を考慮したシステム同定問題は次の制約付き最適化問題に帰着できる。min
(
)
,
.
(
)
0
1 2
k
to
subj
k
e
N k
（3.1.2.1） ただし、N はデータ数である。本論文では、飽和誤差 に関する項を導入した評価関数
N k N kk
h
k
e
J
1 1 2)
(
)
(
（3.1.2.2） を設定する。ただし、h
(k
)
は次式で定義される
(k
)
に対するペナルティ関数である。
0
)
(
,
0
0
)
(
,
)
(
2k
if
k
if
k
h
)
(k
y
)
(
ˆ k
y
ただし、
は定数重みである。ペナルティ関数を導入したことによって、制約つき最適化問 題を制約なし最適化問題に変換した点が重要である。 ここで、（3.2.2.1）式の目的関数は未知パラメータに関して２次形式とはならず、（3.2.2.1） 式は２次計画問題とはならない。そこで本論文では、ペナルティ関数を導入し、非線形計 画法の一つである準ニュートン法を用いて（3.2.2.2）式を最小化することにより係数パラ メータと飽和量を同時に推定することにする。本論文で準ニュートン法を用いた理由は以 下の通りである。まず、アルゴリズムが簡単で広く用いられている。つぎに、ニュートン 法では目的関数のヘッセ行列が悪条件となるような数値的問題が発生する可能性があるが、 準 ニ ュ ー ト ン 法 に お い て 、 ヘ ッ セ 行 列 の 逆 行 列 を 近 似 す る 更 新 公 式 と し て BFGS(Broyden-Fletcher-Goldfarb-Shanno)公式を利用することにより、数値的な悪条件の問題 を回避することができる。 なお、飽和を考慮しない通常の最小２乗法による推定値をパラメータ推定値の初期値と し、
(k
)
の初期値はともに０とする。 いま、目的関数の次数が２であり、ペナルティ関数の次数も２であるため、ペナルティ関 数の効果が現れない可能性もある。そのような場合には、ペナルティ関数を対数関数や指 数関数などに変更すればよい。3.1.3 シミュレーション結果
対 象 シ ス テ ム と し て は １ 次 遅 れ 系 を 仮 定 し 、 そ の パ ラ メ ー タ 値 は9708738
.
0
1
a
,b
1
0
.
0291262
とし、出力端飽和を考慮した。入力信号として正規性白 色雑音を用いた。また、センサの線形領域をL
0
.
07
と設定した。すなわち、これらを超 える範囲ではデータは飽和するものとした。また、同定に用いる式誤差の数はN
100
と した。ただし、定数重みは
10
とした。ここで終了条件としては、目的関数値の変化が 510
1
以下になったとき、更新を終了した。図 3.1.2 に入出力信号を示す。図中、real が真の出力信号、measured は飽和を持つ非線形センサにより計測された出力信号、newton が推定された飽和量を用いて再構成された出力信号、LS は飽和された出力信号を従来の最 小二乗法に適用し同定したモデルの出力信号を示している。図より、LS では誤差が見られ るが、newton は real とほぼ重なっており、真の出力信号を推定できていることがわかる。図 3.1.2 入出力信号 また、図 3.1.3 に周波数伝達関数の同定結果を示す。図より、飽和特性を考慮しない従来の 最小二乗法（LS）では、周波数特性を同定できていないが、非線形システム同定法(newton) を用いると、真値（real）とほぼ等しい周波数特性が同定出来ていることがわかる。表 3.1.1 に定量的評価とパラメータ推定誤差率の比較を示す。また、図 3.1.4 に出力との誤差を示す。 図より LS では誤差が見られるが、newton は誤差が 0 となっている。これらからも、非線 形システム同定法(newton)は対象のパラメータを精度よく同定できていることが確認でき る。
図 3.1.3 周波数伝達関数の同定結果 表 3.1.1 定量的評価 未知パラメータ 真値 最小２乗法 （推定誤差率） 非線形同定法 （推定誤差率） a1 -0.9708738 -0.9456740(
2.60 %
) -0.9709301(0.01 %
) b1 0.0291262 0.0239106(17.91 %
) 0.0291270(0.00 %
)3.2 分解能を有する非線形センサの場合
3.2.1 問題の定式化
まず、図 3.2.1 に示した不感帯特性の場合を考える。ただし、不感帯領域は
D
とし、0
D
の値は既知であるとする。図 3.2.1 から明らかなように、非線形センサによって測 定された出力データy
(k
)
は不感帯による値のずれを有しているので、それらを次式のよう に補正する必要がある。
0
)
(
,
)
(
0
)
(
),
(
0
)
(
,
)
(
)
(
'
k
y
if
D
k
y
k
y
if
k
y
k
y
if
D
k
y
k
y
（3.2.1.1） これらの処理により、不感帯の範囲外で測定されたデータから、不感帯によるずれを除去 することができる。 すると、真の出力信号の推定値y
ˆ k
(
)
は次式のように表わすことができる。
D
k
y
if
k
y
D
k
y
if
k
k
y
)
(
'
),
(
)
(
'
),
(
)
(
ˆ
（3.2.1.2） ただし、
(k
)
は、出力端において不感帯により測定できない微小データを推定するために、 新たに導入された変数であり、次の制約条件を満たさなければならない。
(
k
)
D
（3.2.1.3） 不感帯特性を考慮したシステム同定問題は次の制約付き最適化問題に帰着できる。e
k
subj
to
k
D
N k
)
(
.
,
)
(
min
1 2
（3.2.2.1） ただし、e(k)は（3.1.1.2）式により計算される式誤差である。この問題を解くために、次のペナルティ関数を導入する。
D
k
if
D
k
if
D
k
k
h
)
(
,
0
)
(
,
)
(
)
(
2
（3.2.2.2） ただし、
は定数重みである。ペナルティ関数を導入したことにより、制約付き最適化問題 を、次の制約なし最適化問題に変形することができる。
N k N kk
h
k
e
J
1 1 ( 2)
(
)
(
min
min
output
－Ｄ
Ｄ input
図 3.2.1 不感帯特性)
(k
y
)
(k
y
3.2.2 分解能を考慮した場合
図 3.2.2 にセンサの分解能を 0.05mm とした時の出力信号を示す。図中、青色（実線で 表示）が真値、黒色（点線で表示）が分解能を有した非線形センサで測定した信号である。 図より測定値はこの分解能より微小なデータは測定できず丸められていることが確認でき る。本論文では 3.2.1 で示した不感帯領域をセンサの分解能により測定できない範囲と考 え,非線形システム同定法を応用する。ただし、係数パラメータの初期値としては、分解能 で丸められた測定値に基づいて最小二乗法から推定された値を用い、
(k
)
の初期値は 0 と する。 図 3.2.2 分解能を有する場合3.2.3 シミュレーション結果
3.1.3 項と同じ１次遅れ系を仮定し、入力信号としてN
(
0
,
0
.
8
2)
の正規性白色雑音を用い た。また、センサの分解能を 0.05mm と設定し、シミュレーションを行った。さらに同定に 用いる式誤差の数はN
1000
とした。図3.2.3 にシミュレーション結果を示す。図中,real は分解能を持たない真値,measured は分解能で丸められた測定値,newton は準ニュートン法 により推定された微小データ β(k)を用いて再構成された出力信号を示している。図より newton と real とがほぼ一致しており,真の出力信号を推定できていることがわかる。図3.2.4 に周波数伝達関数の同定結果を示す。図中,LS は丸められた出力に従来の最小二乗法を適用 し同定したモデルの周波数特性である。この結果からも周波数特性が改善されていること がわかる。また,表 3.2.1 に各同定パラメータに対する定量的比較結果を示す。 図3.2.3 出力信号図3.2.4 周波数伝達関数の同定結果 表 3.2.1 定量的評価 未知パラメータ 真値 最小２乗法 （推定誤差率） 非線形同定法 （推定誤差率） a1 -0.9708738 -0.9696640 (
0.12 %
) -0.9697482 (0.12 %
) b1 0.0291262 0.0252060 (13.46 %
0.0292092 (0.29 %
)また、対象システムとして次式のような 2 次遅れ系を仮定し,制御入力信号として
)
8
.
0
,
0
(
2N
の正規性白色信号を用いてシミュレーションを行った。 (3.2.3.1) また,シミュレーション条件としてサンプリング時間を 0.5ms,センサの分解能を 0.05mm とした。つまりこの分解能より微小なデータは測定できず丸められる。図3.2.5 にシミュレ ーション結果を示す。図中,real は分解能を持たない真値,measured は分解能で丸められた測 定値,newton は準ニュートン法により推定された微小データ β(k)を用いて再構成された出力 信号を示している。図より newton と real とがほぼ一致しており,真の出力信号を推定できて いることがわかる。図3.2.6 に周波数伝達関数の同定結果を示す。図中,LS は丸められた出 力に従来の最小二乗法を適用し同定したモデルの周波数特性である。この結果からも周波 数特性が改善されていることがわかる。また,表 3.2.2 に各同定パラメータに対する定量的 比較結果を示す。また, 図 3.2.7 に真値との出力誤差を示す。LS では誤差が見られるが、 newton はほぼ 0 に近い結果となっている。これらからも対象の特性をよく近似しているこ とが確認できる。 図3.2.5 出力信号
1
.
0
150
2
)
(
2 2 2
n n n ns
s
s
P
図3.2.6 周波数伝達関数の同定結果 表 3.2.2 同定パラメータの比較 真値 最小２乗法(推定誤差) 非線形同定法 （推定誤差） a1 1.948 -1.919(1.46 %) -1.931(0.86 %) a2 0.970 0.923(2.89 %) 0.953(1.72 %) b1 0.011 0.009(10.68 %) 0.010(1.42 %) b2 0.011 0.010(1.10 %) 0.011(0.41 %)
第
4 章 制御対象の概要と同定実験
本章では、実験装置として用いる精密ステージについて簡単に説明し、分解能を考慮し た非線形システム同定法で導出した同定モデルと従来法で導出した同定モデルの比較検証 のために同定実験を行い制御対象の真値の同定モデルを導出した。4.1 SPIDER 駆動精密ステージの概要
本研究に用いた実験装置の構成図を図4.1.1 に示す。ホスト PC から送られた入力指令は、 PCI スロットに装着したパラレル IO カードを利用して、サーボインターフェスユニット（モ ーションコントローラ）、アンプを通してSPIDER に送信される。ここのアンプにおいて入 力指令電圧は、130/10 倍（以下 13 倍）されて送られる。また、エンコーダ及びリミットセ ンサよりステージの位置情報ストロークリミット情報が読み込まれ、ホストPC に送られる。 ホストPC には OS として Windows98 を搭載した PC を用い、VisualC++により実行プログラムを作成している。I/F カードには Interface Corporation 製 16/16bitI/O PCI ボード PCI-2735 を使用している。圧電素子を使用しているアクチュエータを含むステージシステ ムは熊本テクノロジー、太平洋セメント社らの共同研究により開発された。
Linear encoder
Linear
guide
Limit
sensor
SPIDER
Stage
Scale
PC
with servo system
Motion controller
control input
position
signal
operatorlimit
signal
Guide
Plate
図4.1.1 ハードウェア構成図SPIDER
Slide
plate
Stage
SPIDER
Slide
plate
Stage
図 4.1.2 精密ステージpreload mechanism
20mm
piezoelectric actuator
preload mechanism
20mm
piezoelectric actuator
図4.1.3 アクチュエータ部拡大図 ステージシステムの写真を図4.1.2 に、アクチュエータ拡大図を図 4.1.3 に示す。位置を 測定するために、リニアエンコーダがステージ稼動部の下面に取り付けてある。リニアエ ンコーダはミツトヨ製で、計測分解能は電気分割ユニットのスイッチ切り替えにより、最 小10nm となっている。制御入力となる圧電素子への最大印加電圧は±130V である。駆動 周波数は1kHz～60kHz まで設定変更が可能であり、ステージストロークは 4 インチウェ ハ対応の約100mm となっている。ステージの仕様を表 2.1 に示し、また圧電素子の性能表 を表2.2 に示す。また、圧電素子の静的な発生力は最大伸縮素子発生力 660N（４脚同時） となる。予圧力が50N、摺動面の摩擦力が 15N であることから足の運動に十分な駆動力が 得られている。また、ガイドプレートの平均表面粗さが約0.2μm であることから、伸縮方 向に39V のオフセット電圧を印加することで表面粗さの影響を軽減している。表4.1.1 ステージの仕様 可動部質量 1kg 駆動周波数設定 1～60kHz 最大推力 13N 最大印加電圧 ±130V ストローク 100mm 位置分解能 100nm 表4.1.2 アクチュエータ（SPIDER）の仕様 材質 PB(Zr,Ti)O3 密度 7.8×103kg/m3 伸縮率 660×10-12m/V 剪断率 1010×10-12m/V 積層枚数 4（伸縮）×4（剪断）
4.2 制御対象の事前同定実験
非線形システム同定法により導出した同定モデルと従来法により導出した同定モデルの 比較のために、精密ステージの真値の同定モデルを導出する必要がある。そのため、この モデルを得るためにシステム同定実験を行う。ここではМ系列信号を入力し,その出力応答 から同定を行うが、静止摩擦を考慮し、分解能誤差の影響を受けないよう十分なオフセッ トを入力電圧に加えた状態での実験を行った。 実験条件は、同定に用いたМ系列信号のオフセットを 6.5V、サンプリング時間を 0.5ms、 センサ分解能を 100nm と設定した。実験より得られた入出力信号を図 4.2.1 に示す。図中、 上部はМ系列信号、下部は位置出力である。分解能に対し位置出力が十分に大きいため分 解能の影響を受けていないと考えられる。また、図 4.2.2 に出力の速度の一部を示す。図 より、速度が 0 より大きいため、静止摩擦の影響を受けていないことが確認できる。また、 図 4.2.3 に周波数伝達関数の同定結果、図 4.2.4 に実験出力とモデル出力の比較を示す。 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 6.4 6.45 6.5 6.55 6.6 6.65 Time [s] Input 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 20 40 60 80 100 Time [s] P o si ti o n [m m ] 図 4.2.1 入出力信号1000 1500 2000 2500 3000 3500 30 32 34 36 38 40 42 44 46 48 50 Number of samples V el oci ty [ m m /s ] 図 4.2.2 速度データ -150 -100 -50 0 50 100 振幅 (d B) 10-4 10-2 100 102 104 -180 0 180 360 540 720 位相 (d e g) ボード線図 周波数 (rad/sec) 図 4.2.3 周波数伝達関数の同定結果
0 200 400 600 800 1000 1200 1400 1600 1800 2000 -10 0 10 20 30 40 50 60 70 80 90 Number of samples O ut put v ol ta ge [V ] experimental simulation 図 4.2.4 実験出力（実線）とモデル出力（点線）
4.3 従来法との比較
4.2 節で求めた真値モデルを用いて、分解能を考慮した非線形システム同定法により導出 した同定モデルと、従来法により導出した同定モデルとの比較を行う。以下にその実験結 果を示す。 実験条件は、同定に用いたМ系列信号のオフセットを 1.7V、サンプリング時間を 0.5ms、 センサ分解能を 100nm と設定した。実験より得られた入出力信号を図 4.3.1 に示す。図中、 上部はМ系列信号、下部は位置出力である。位置出力が小さいため、分解能の影響を受け ていると考えられる。図 4.3.2 に出力の速度を示す。図より、静止摩擦の影響を受けてい ないことが確認できる。また、図 4.3.3 に実験出力と各同定モデルとの出力誤差の同定結 果を示す。すなわち,分解能の影響を受けない状態での同定モデル出力（real）,同一入力に 対する非線形システム同定法により導出したモデルの出力（newton）,最小二乗法により導 出したモデル出力（LS）を比較して示す。また, 図 4.3.4 には対応する各同定モデルの周 波数特性を示す。これらから newton は時間応答,周波数特性ともに real とほぼ重なっている ことがわかる。0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 1.6 1.65 1.7 1.75 1.8 Time [s] Input 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.5 1 1.5 2 Time [s] P os it io n [m m ] 図 4.3.1 入出力信号 0 500 1000 1500 2000 2500 3000 3500 4000 0 0.5 1 1.5 2 2.5 Number of samples V el oci ty [ m m /s ] 図 4.3.2 速度データ
図 4.3.3 出力誤差
![図 5.1.1 摩擦特性 図 5.1.2 スティックスリップ現象 -40-30-20-1001020 30 40- 80- 6002040Friction u [V]- 40- 206080161301拡大図Steady-state velocity v[mm/s]-40-30-20-10010203040u [V]161301161301Steady-state velocity v[mm/s]拡大図](https://thumb-ap.123doks.com/thumbv2/123deta/6324805.1097798/46.892.208.645.152.488/図511摩擦特性図512スティックスリップ現象Friction大図SteadystateSteadystate拡大.webp)
