平成１４年度 国立天文台 ALMA 共同開発研究 研究成果報告書
研究課題名：外乱に強い高速高精度アンテナ制御技術の開発
研究代表者： 岐阜大学工学部 佐々木実
目 的 ： ALMA の ア ン テ ナ で の ト ッ プ ３ の 課 題 は 、「 鏡 面 精 度 ＜ ２ ０ ミ ク ロ ン 」「Pointing/Tracking 精度＜０．６ 」「Fast Switch < 3 deg within 1.5 second」という性能の達成である。後者２項目について、新 10m 鏡に載せた加速度計や角度エンコーダで、計測を中心にして基礎的な評価が行われた。
Preliminary な結果として、（１）風外乱等に対して主鏡部は十分に剛構造である。
（２）通常追尾時の Pointing/Tracking はほぼ仕様を満足するものの、Fast Switching 時に 発 生 す る 振 動 の 収 束 が 十 分 で は な い た め Tracking 精 度 が 仕 様 を 満 た さ ず 、 FastSwitching 観測での時間効率が高くない、などが指摘されている。
つまり、新 10ｍアンテナは、制御部の性能が不十分であり、これを改善すれば観測性能の向上が期待できる。
本研究の目的は新 10m アンテナの駆動制御性能の向上である。このために現代制御理論を取り入れた高性能な制御の開発を行う。
アンテナ特性診断、最適制御回路の設計、その実装と一連の技術開発が完成すれば、単に新１０ｍアンテナのみならず、他の高速スイッチング性能が必要なアンテナにも応用できる。Fast switching 時に発生する振動の収束が十分ではないため Tracking 精度が仕様を満たさず，Fast Switching 観測での時間効率が高くない．
よって本研究では 10m アンテナの駆動制御の向上である．
方法

上記目的を達成するために，新１０ｍアンテナの加速度計並びにエンコーダ出力データからアンテナシステムの同定・解析を行い、Fast Switching 時の振動追尾誤差の収束の早い入力軌道の制御系の設計を行う．
システム同定方法（予測誤差法）
あるシステムに対して線形離散時間確率モデル

y (t ) = G ( z,θ )u (t ) + H ( z,θ )e(t )(1.1)

を仮定する．ここで，e(t)は，平均値 0，分散σ，のホワイトノイズ，θは未知パラメータベクトルで，d=dimθとするとき，θは Rd の部分集合Θ⊂Rｄに属するものとする．時刻 tまでの入出力データの集合を Zt=u(t)，y(t)，k=t，t-1，･･･とおく．式(1.1)のモデルを用いて，出力の一段予測値を計算し，式(1.1)を変形して移項すると，

y (t ) = G ( z ,θ )u (t ) + [ H ( z ,θ ) − 1] + e(t ) + e(t )(1.2)となる．式(1.2)を e(t)について求めると，e(t ) = H −1 ( z ,θ ) y (t ) − H −1 ( z ,θ )G ( z ,θ )u (t )(1.3)となり，その式(1.3)の e(t)を式(1.2)中に代入すると，y (t ) = H −1 ( z ,θ )G ( z ,θ )u (t ) + [1 − H −1 ( z ,θ )] y (t ) + e(t )(1.4)となる．式(1.4)の右辺の第一項，第 2 項の伝達関数は，次のような関係で表される．

lim H −1 ( z ,θ )G ( z ,θ ) = 0(1.5)

lim[1 − H −1 ( z ,θ )] = 0

(1.6)

z →∞

z →∞

式(1.5)，式(1.6)より，H-1(z,θ)G(z,θ)，1-H-1(z,θ)はともに厳密にプロパーであることになる．従って，式(1.4)の右辺第一項，第 2 項は，時刻 t-1 までの入出力のデータ Zt-1 の関数である．さらに，e(t)は仮定より Zt-1 には無相関であるから Zt-1 に基づく y(t)の一段階予測値
は，∧y (t / θ ) = H −1 ( z ,θ )G ( z ,θ )u (t ) + [1 − H −1 ( z ,θ )] y (t )(1.7)

で与えられる．従って，一段階予測誤差は，

ε (t ,θ ) H −1 ( z ,θ )[ y (t ) − G ( z,θ )u (t )](1.8)となる．ここで，形式的には式(1.4)での e(t)をε(t,θ)で置き換えたものに等しい．式(1.8)
から予測誤差は H(q)と G(q)の関数になっていることがわかる．ここで，入出力データ

Z N = {u (t ), y (t ), t = 0,1,･･･, N }(1.9)

が与えられた時，式(1.7)から予測誤差 e(t)，t=0，1， ･･･，N が計算できる．これよりモデ
ルを観測データ ZN にあてはめた時，明らかに小さな予測誤差を与えたほうがよいモデルであることがわかる．このようなことから，予測誤差法は予測誤差の分散を小さくするようにパラメータを決定する方法であることがわかる．つまり，この H(q)と G(q)をパラメータ化したものを調整することにより予測誤差の分散を減少させる方法である．これは，次の評価関数の値を減少させることを意味している．

N

V (G, H ) = ∑ e(t ) 2(1.10)

t =1

ARMAX モデル構造では，パラメータ化すると予測誤差は，e(t ) =

A(q)
B(q)
y (t ) −u (t )
C (q)
C (q)

(1.11)となり，このときの評価関数は，
N

V = ∑ e(t ) 2(1.12)

t =1となる．
１入力１出力の ARMAX モデル

y (t ) + a1 y (t − 1) + ･･･+ a na y (t − na) = b1u (t − nk ) + ･･･+ bnb u (t − nk − nb + 1) + e(t ) + c1 (t − 1) + ･･･+ c nc (t − nc)
2 入力１出力の ARMAX モデル

y (t ) + a1 y (t − 1) + ･･･ + a na y (t − na) = b11u (t − nk1 ) + ･･･ + b1nb u (t − nk1 − 1nb + 1) + b21u (t − nk 2 ) + ･･･+ b2 nb (t − nk 2 − 2nb + 1) + e(t ) + c1 (t − 1) + ･･･ + c nc (t − nc)

アンテナエンコーダー・センサーの位置

図 1.az 駆動角 0.3 度 el10 度

図 3.1.1 入力 1 出力の同定結果

図 3.2..2 入力 1 出力の同定結果

同定結果１
１入力時は AZ 方向のエンコーダーの位置，２入力時は AZ 方向，EL 方向のエンコーダーの位置，出力は①のセンサー振動を用いた．青が実際のシステムの応答,緑が同定結果．評
価基準として,以下の式を用いる.このとき時間がたつにつれ振動の隆起がみられなくなるので,200 項までを対象とする.
200Re = ∑ ( y i (t ) − θ i (t )) 2i =1

1 入力 1 出力の同定結果
Re=6.3373e+007
2 入力 1 出力の同定結果
Re=6.5076e+007
求められた同定結果の伝達関数の応答を図 4，図 5 に示す az1.5 度 el10 度,az6 度 el10 度の新たな入力と出力を使い,同定精度を確かめる．図 5.1 は az1.5 度 el10 度の時の 2 入力 1 出力による結果，図 5.2 は az6 度 el10 度の時の 2 入力 1 出力の結果である．又，同様に図 5.3
は az1.5 度 el10 度の時の 1 入力 1 出力の結果, 図 5.4 は az6 度 el10 度の時の 1 入力 1 出力の結果である.

図 4. az 駆動角 1.5 度 el10 度

図 4. az 駆動角 6 度 el10 度

図 5.1 2 入力系 az1.5 度 el10 度同定結果

図 5.2 2 入力系 az6 度 el10 度同定結果

図 5.3 1 入力系 az1.5 度 el10 度同定結果

図 5.4 1 入力系 az6 度 el10 度同定結果

以上の結果からわかるように，系の挙動が安定していない結果がほとんどである．ここで，図 3 のシステム結果を見てみると，振動挙動の始まりの部分が同定結果に反映していない，それは，入力より出力が早く始まっている事からだと推測できる．そこで，入力結果に対し，時間軸に修正を加えて同定した結果を以下に示す．

図 6.1 2 入力 az0.3 度 el10 度同定結果

図 6.3 2 入力 az6 度 el10 度同定結果

図 6.2 2 入力 az1.5 度 el10 度同定結果

図 6.4 2 入力 az1.5 度 el30 度同定結果

図 6.5 2 入力 az6 度 el60 度同定結果

同定結果２
図 6.1 は２入力によるプラントの同定結果である．図 6.2 は図 6.1 で求められた同定シス
テムの精度を確かめるために新しく 2 入力 az1.5 度 el10 度の入力を使用した結果である．
同様に，図 6.3 は 2 入力 az6 度 el10 度を使用した時の結果である.入力はそれぞれ 0.08 秒速くなるように時間軸に修正がかけられている.
結果として，同定精度は以前のものより向上したと思われる．
続いて，この同定システムの適用範囲を調べる．図 6.4 は駆動角 az1.5 度 el30 度(固定)，図6.5 は駆動角 az6 度 el60 度(固定)の結果である．青い線が実際の出力，緑の線が同定結果であり，この結果より，同定器は el の角度が異なると適用は難しくなる．よって，角度位置による，別の同定システムが必要となるかもしれない．

ニューラルネットワーク
ニューラルネットワークは，脳・神経系による情報処理方式に学び，多数の神経細胞(ニューロン)が網目状に結びついたネットワーク・システムをモデル化したものである．
階層型ニューラルネットワーク：信号が入力層(Input layer)，隠れ層(Hidden layer)，出力層(Output layer)の順に一方向に流れる．そして，ニューロンのモデルをユニットと呼び，ユニットはニューラルネットワーク上の交点にあたり，結合の強さを表す結合係数をつけられてユニットへ入力される．結合係数をつけられたそれぞれの入力の総和は，応答関数によって変形を受けた後出力される．ここでは，4 層の階層型ニューラルネットワークを用いた．ここで，Ii(i=1･･･n)は各入力層ユニットへの入力， Hj(j=1･･･m1)は１番目の各隠れ層ユニットの出力，Pk(k=1･･･m2)は 2 番目の各隠れ層ユニット出力，Ol(l=1･･･r)はニューラルネットワークの出力である．(以下では，添え字 i は入力層のユニット番号，j は隠れ層のユニット番号を表す)．隠れ層の各ユニットへの入力(内部ポテンシャル)shj 及び出力層への入力 s0 は
ns hj = ∑ I i･Wij + θ j(2.1)

i =0m2s 0l = ∑ Pk･Vk(2.2)

k =0となる．ただし，Wij は入力層から隠れ層への結合係数，Vk は隠れ層から出力層への結合係数であり，θj は閾値である．また，隠れ層の応答関数 fh(x)及び出力層の関数 f0(x)には次式に示すものを用いた．

f h ( x) =1
1 + exp(− x)f 0 ( x) = x

(2.3)
(2.4)

隠れ層に用いた関数はシグモイド関数と呼ばれる．このようなシグモイド関数を用いることによりニューラルネットワークは非線形連続関数を実現できることが知られている．また、この関数の特徴として，式(2.5)に示すように，シグモイド関数の微分式がシグモイド関数そのもので表現できることが挙げられる．

∂f h ( x)exp(− x)== f h {1 − f h ( x)}∂x{1 + exp(− x)}2(2.5)

式(2.1)，(2.2)，(2.3)，(2.4)より，隠れ層，出力層の出力は

H j = f h ( s hj ), O = f 0 ( s 0 )(2.6)となる．
次に，ニューラルネットワークが学習を行うには，何らかの評価基準が必要になり，その結果をフィードバックしてユニット間の結合係数を調整する．ここでは，学習アルゴリズムとして BP(Backpropagation)法を用い，前述のニューラルネットワークに適用する．
教師信号を d とすると，このときの誤差信号は

e = d −O(2.7)となる．これより誤差関数 Ep を 2 乗誤差を用いてEp =1 2e2(2.8)と 定 義す る． BP 法 は最 急 降下 法 (Steepest descent method)を 利 用 し て ， Ep の 勾 配

∂E p
∂Wij，

∂E p
∂V jから，Ep が減少する修正量を求めるものである．結合係数 Vj の微小変化に対する誤差関数の影響は

∂E p
∂V j=

∂E p ∂o ∂s 0・ ･∂o ∂s 0 ∂V j= −e･1･H j(2.9)(2.10)となり，これより更新式は，学習係数をαとして

V jnew = V jold + ∆V j(2.11)∆V j = −α･e･H j(2.12)となる．これを用いて逐次的に修正することによって

−

∂E p
∂V j∝ ∆V j

(2.13)

を得る．つまり，式(2.8)の 2 乗誤差は式(2.11)，(2.12)に基づく学習によって極小値に向かって減少していくことがわかる．また，結合係数 Wij についても同様に，誤差関数を結合係数で偏微分すると

∂E p
∂Wij=

∂E p ∂o ∂s 0 ∂H j ∂s hj･ ･･･∂o ∂s 0 ∂H j ∂s hj ∂Wij(2.14)= −e･V j･f ' h ( s hj )･I i(2.15)となる．式(2.14)，(2.15)は，式(2.5)のシグモイド関数の微分を用いると

∂E p
∂Wij

= −e･V j･f h ( s hj )･{1 − f h ( s hj )}･I i(2.16)となる．従って，更新式は学習係数を同様にαとして次式で求められる．

Wijnew = Wijold + ∆Wij(2.17)∂E p∆Wij = −α・∂Wij(2.18)= α・e・V・) {1 − f h ( s h )}・I ij f h (sh ・(2.19)

また，式(2.1)，(2.2)の閾値θＪは上式の入力層からの入力 I i を 1 とおいてθ jnew = θ old+ ∆θ jj(2.20)∆θ j = α・e・V・) {1 − f ( s hj )}j f ( s hj ・(2.21)と修正され，結合係数と同様に扱うことができる．
以上から，BP 法における結合係数と閾値の修正量は出力層における誤差の寄与量から求められている．
ニューラルネットワークによるシステム同定と制御

図７に示すとおり u が入力，y が出力，gy が推定出力となり，その誤差をニューラルネットワークで学習させることで，Plant の同定システムをニューラルネットワーク内に構築できる．また，その同定システムを用いて図 7.1 のような形の系を用いて適応的に入力を整形させ，制御する．

uyPlantgy

NeuralNetwork１

E=(y-gy)^2
図 7．ニューラルネットワーク同定器
ugu

NeuralyPlant
e=（y-gy）^2Network２

NeuralNetwork１gy
E=gy^2
図 7.1．ニューラルネットワーク同定・制御器

図 8．ニューラルネットによる同定結果

図 9．ニューラルネットワークモータの位置と振動図 10．ニューラルネットワークモータの位置

図 11．ニューラルネットワーク角速度と振動図 12．ニューラルネットワーク角速度と入力
図 8 はニューラルネットワークを同定器として使用した結果である.
青が az モータの位置(入力)，緑が el モータの位置(入力)，赤が実際の振動，水色がニューラルネットの推定値である．よく追従していることが分かる．
次に，ニューラルネットワークの同定結果に対してのコントローラ取り付け時の結果である．
図 9 の青い線は az モータの位置入力，緑の線は el モータの位置入力，赤い線が実振動，水色の線は結果である．
図 10 の青い線は az モータの位置入力，緑の線は el モータの位置入力，赤い線が az の推定位置入力である．
図 11 の青い線は az モータの角速度入力，緑の線は el モータの角速度入力，赤い線が実振動，水色の線は結果である．
図 12 の青い線は az モータの角速度入力，緑の線は el モータの角速度入力，赤い線が azの推定角速度入力である．
コントローラはまだ上手く設計できているとはいえず，よい結果が得られているとはいえない．コントローラの入力となる部分がサーボ設計と同じになるため，それを考慮に入れながらのコントローラの設計にはまだ検討が十分ではなく，時間がかかると思われる．
図 13 にシステム同定結果を用いて逆システムをフィードフォワードコントローラとして用いた結果を示す．青い線は，目標加速度値に振動データ成分をいれた場合で，.緑の線は,目標加速度振動を０としたときの結果である．

図１３．逆システムを用いた振動制御結果
以上より，同定器，ニューラルネット共にプラントのモデル化は拡張性の問題は残っているものの，推定できるものと思われる．今後の問題点としてはコントローラの設計時に目標値入力を修正しながら，振動抑制のできるものにしなければならない．そのために，ニューラルネットワークで使用される評価関数の見直し，コントローラ再設計をしていきたいと考えている．