計測自動制御学会産業論文集
Vol. 10, No. 19, 160/166 (2011)

特集「移動体の誘導制御」

非線形計画法と Unscented 変換を用いた移動体のパラメータ推定†
横山

信 宏*

Parameter Estimation of a Vehicle Using Nonlinear Programming with Unscented Transform
Nobuhiro YOKOYAMA*This paper presents a new method for parameter estimation of the dynamics of a vehicle that is expressed by nonlinear
state and observation equations with unknown parameters. The developed method employs a smoothing technique based on
nonlinear programming with an unscented transform. It covers general classes of problems, including those which havenonlinear process noise, unknown covariances of the process and measurement noise. By applying the developed method tothe problems of estimating a characteristic parameter of a reentry body and the aerodynamic coefficients of a research
aircraft, good identification accuracies comparable to those of established methods as well as a stable convergence property
are demonstrated.
Key Words. Parameter estimation, nonlinear programming, smoothing, unscented transform

音とプロセス雑音の両方を考慮しなければならない場合も多1.はじめにく，その場合は Filter Error Method (FEM) ４) や拡大系に対す

移動体の運動シミュレーションや制御系設計においては，一般にその動特性を記述する数学モデルをあらかじめ構築しておく必要がある．数学モデルの構築においては，物理法則などに基づいてモデルの構造をあらかじめ決定できる場合とそうでない場合とがあり 1, 2)，それぞれにおいて異なるアプ
ローチが採られる．前者の場合は，あらかじめ構造が定められた数学モデル中の未知パラメータを何らかの手法を用いて決定するが，後者の場合は，数学モデルを表す式の構造と未知パラメータを，同時に決定する手法が用いられる．たとえば，航空機におけるシステム同定では，よく知られた運動方
程式を用いて精度の高い数学モデルを構築できるため，前者のアプローチが採用されることが多い 3)．本稿では，前者のアプローチに関し，数学モデルが非線形の場合にも適用可能なパラメータ推定手法について述べる．
非線形システムに対するパラメータ推定手法に関しては，いくつもの研究例がある．Equation Error Method や Output
Error Method は，実装が比較的容易で，よく普及している手
法である
†

3)

．しかしながら，実システムにおいては，観測雑るフィルタリング 5-7) およびスムージング 8, 9) が有効である
ことが知 られ ている． FEM で は ，Extended Kalman Filter(EKF)や Unscented Kalman Filter (UKF)6) を用いることで，プ
ロセス雑音の影響を除去した推定出力に対して最尤法を適用する．ただし，FEM においては，初期時刻から終端時刻までの状態方程式の推移を直接計算するため，扱う最適化問題の非線形性が強くなりやすく（たとえば，初期状態変数の推定値の僅かな更新が，終端時刻付近の状態変数において非線形性の強い変化を引き起こすなど），収束のロバスト性が比較的低いという問題点がある．一方，拡大系に対するフィル
タリングおよびスムージングにおいては，未知パラメータを
状態変数の一部に含めた拡大系に対して，EKF や UKF などの フ ィ ル タ ， あ る い は Unscented Raugh-Tang-Striebel
Smoother (URTSS)8) などのスムーザーを適用する．フィルタ
リングにおいては，時々刻々の観測データから逐次的に推定を行うのに対し，スムージングにおいては，全観測データが得られた後でそれらすべてに基づいて推定を行うことから，より高い推定精度を得られる傾向がある 8)．なお，フィルタリング・スムージングのいずれにおいても，初期状態変数や第 27 回誘導制御シンポジウム(2010・5)および 18th IFAC

雑音の共分散などを推定するために，別の手法が必要な点にSymposium on Automatic Control in Aerospace(2010・9)にて注意が必要である 9)．
本稿においても，観測雑音とプロセス雑音の両方に対応し一部発表*

防衛大学校航空宇宙工学科*

横須賀市走水 1-10-20

た，スムージングに基づく推定手法を提案する．ただし，提Department of Aerospace Engineering, National Defense案手法では，上述の手法のように拡大系に対してスムージン

Academyグを行うのではなく，非線形計画法を利用して，本来の系に(Received November 29, 2010)

対して直接的なスムージングを行う．文献 10 においても，TRIA 019/11/1019 © 2010SICE

非線形計画法に基づくスムージング手法が述べられているが，

N

N

k 0

k 1

log p( X N , θ Z N )   log p( zk | xk , θ )   log p ( xk | xk 1 , θ )提案手法は，非線形なプロセス雑音および雑音共分散の推定(6)

 log p ( x0 )  log p (θ )  log p ( Z N )を含むより一般的なクラスの問題に適用することが可能であ(2)式より，与えられた xk と θ のもとで zk は vk を線形変換

り，この点において，このクラスの問題に適用することがでだし，文献 9 の手法で非線形なプロセス雑音を扱う場合には，した Gaussian となる．すなわち， p( zk | xk , θ ) はexp[( zk  hk )T R 1 ( zk  hk ) / 2]p ( zk | xk , θ ) (2 ) p / 2 (det R)1/ 2

更にシステムを拡大して，推定すべき系の次元を増大させるとなり 11)，さらに log p ( zk | xk , θ ) はきない文献 10 の手法に対する優位性を有する．文献 9 の手法や FEM でも，同じクラスの問題を扱うことができる．た必要があるのに対し，提案手法ではその必要が無い．また，

log p ( zk | xk , θ )提案手法では，推定される状態変数の遷移を，FEM のよう1  [( zk  hk )T R 1 ( zk  hk )  log det R  p log(2 )]
2に状態方程式を直接数値積分して求めるのではなく，最適化

(7)

(8)

で表される．ただし， hk  h( xk , uk , θ ) とする．

計算における等式制約条件を通じて間接的に求めるため，
FEM と比較して収束のロバスト性が良いという利点もある．

つぎに， log p ( xk | xk 1 , θ ) に関しては， xk が Gaussian で近本稿の次節以降は以下のように構成される．第 2 節では提似可能であると仮定すれば，以下に述べる Unscented 変換 7)

案手法について述べる．第 3 節において，2 種類の移動体のにより計算することが可能である． (3) 式で定義されるプロ

パラメータ推定問題を扱い，提案手法の有効性を示す．第 4セス雑音 wk の分布を代表する点として，下記の (2q  1) 個の節で結論を述べる．

シグマ点を定義する．
: i  0,0
1/ 2
ωi   q   Qi
:1  i  q,
1/ 2
 q   Qi  n : q  1  i  2q.

2.推定手法
状態方程式および観測方程式は以下の離散系として与えられるとする．
xk 1  f ( xk , uk , θ , wk )(1)

zk  h( xk , uk , θ )  v k(2)

ただし， Qi1/ 2 は共分散行列 Q の平方根行列の i 列目を表す列
ベクトルである．また，    2 (q   )  q であり，  はシグ
マ点の広がりを調整するパラメータを，  は確率分布の形状

ただし，添え字 k ( 0,, N ) はステップ番号， xk   n は状態に応じて調整すべきパラメータを表す．さらに，確率分布の変数， uk   m は制御入力， zk   p は観測出力， θ   l は未平均および共分散を求める上で必要な重み係数を以下で定義

知パラメータ， wk   q はプロセス雑音, vk   r は観測雑音する．

をそれぞれ表す． w k , v k は，平均が零で共分散がそれぞれQ   q q , R   r  r で与えられる白色 Gaussian であると仮定する．すなわち，w k  N (0, Q )(3)v k  N (0, R)(4)とする．また，実験によって入出力の時歴 U N  {u0 , , uN } ,Z N  {z0 , , z N } が観測されているものとし， U N は不確かさ

2  / (q   )  1     : i  0
Wi ( c )  :1  i  2q
 1 / [2(q   )](11)

は，次式で定義される (2q  1) 個の点で代表する．

め除去しておくものとする．

χ k ,i  f ( xk 1 , uk 1 , θ , ωi ), (i  0,, 2q)

提案するパラメータ推定手法は，与えられた Z N のもとで，
未知パラメータ θ と状態変数の時歴 X N  { x0 , , x N } が実現(12)

これらの点を用いることで， p ( xk | xk 1 , θ ) に従う xk の平均される事後確率密度関数 p( X N , θ | Z N ) を最大化するアプローxˆ k および共分散 Pk は，次式から求められる．

チに基づく．Bayes の法則およびシステムの Markov 性によ2nり， p( X N , θ | Z N ) は以下のように変形できる．
p( X N , θ Z N )xˆ k  Wi ( m ) χ k ,i(13)

Pk  Wi ( c ) ( χ k ,i  xˆ k )( χ k ,i  xˆ k )T

(14)

i 0

2n

 p( Z N | X N , θ ) p( X N , θ ) / p( Z N )(5)式の両辺の対数をとると，次式を得る．

(10)る際に用いられる．与えられた xk 1 , θ のもとでの xk の分布雑音が存在する場合は，フィルタ処理などによってあらかじ

 N
 N

  p( zk | xk , θ )   p( xk | xk 1 , θ )  p( x0 ) p(θ ) / p ( Z N ) k 0  k 1
:i  0  / (q   )Wi ( m )  
 1 / [2(q   )] :1  i  2q

Wi ( m ) は平均値を求める際に用いられ， Wi ( c ) は共分散を求めのない確定的な変数の時系列であるとする．入力に強い観測 p ( Z N | X N , θ ) p ( X N | θ ) p (θ ) / p ( Z N )

(9)

i 0

(5)

f が非線形関数であっても，(13),(14)式で与えられる平均と共分散は真値を高い精度で近似できることが知られている 7)．
ただし，(14)式は半正定性を保証するが，正定性を保証しな
いため，微小な正の定数  を用いて次のように修正する．
161

2n

(15)

制約条件式 (19) で関係付けられているように， ek は状態変以上に基づき， xk 1 , θ が与えられたときの xk を Gaussian と(17) 式に現れる ( xk  xˆ k )T Pk1 ( xk  xˆ k ) の項は， Pk1 の条件数して近似すると， log p( xk | xk 1 , θ ) はが原因でしばしば最適化を行うことが困難となりやすいが，

Pk  Wi ( c ) ( χ k ,i  xˆ k )( χ k ,i  xˆ k )T   I n

数 xk の平均 xˆ k からの偏差を Pk1/ 2 で正規化した変数を表す．

i 0log p ( xk | xk 1 , θ )変数 e1 , , e N と制約条件式(19)を導入することで，(18)式の評1  ( xk  xˆ k )T Pk1 ( xk  xˆ k )  log det Pk  n log(2 ) 2(16)

価関数では，安定した最適化計算が容易な ek T ek の形式に置き換えられている．

で表すことができる．

本手法では，状態変数の推移を，初期時刻から終端時刻ま(8) 式および (16) 式を (6) 式に代入し， p( Z N ) は推定に無関

で状態方程式(1)の繰り返し計算で求めるのではなく，(19)式係であるため定数として扱うと，最大化すべき事後確率密度の等式制約条件によって間接的に求めている．これは軌道最p( X N , θ | Z N ) の対数関数として，適化の分野において良く用いられる Direct Collocation と呼ば

log p( X N , θ | Z N )れる手法と同じであり 12) ，最適化問題における状態方程式

1 N
    ( zk  hk )T R 1 ( zk  hk )  ( N  1)log det R
2 k 0の非線形性の影響を緩和することができる．また，最適化問 {( xk  xˆ k )T Pk1 ( xk  xˆ k )  log det Pk }
k 1 log p( x0 )  log p (θ )  const
N

題(18)-(21)における制約条件の Jacobian や Lagrange 関数の(17)

Hessian は，Direct Collocation に基づく軌道と静的パラメータの同時最適化問題 13)と同様のスパース構造を持っている．
そこで，数値最適化アルゴリズムとして，上述の同時最適化が得られる．ここで， p ( x0 ) および p(θ ) については，それ問題と同様に，Sparse SQP (Sequential Quadratic Programming)

らの確率密度関数が事前にわかっている場合は(17)式の中に12, 13)

取り込めば良いが，不明な場合は一様分布であると仮定し，のスパース構造を利用して，高速な最適化計算を行うことがを適用する．このアルゴリズムは，Jacobian や Hessian定数として扱う． x0 , θ に関する事前の確率密度関数の情報可能であり，収束のロバスト性も優れていることが知られては，一般に推定精度を高める上で有益であるが，一様分布といる. なお，アルゴリズム中における微係数の数値計算には，仮定する場合，推定精度を高めるのに寄与しないことに注意

有限差分法ではなく，アルゴリズムの自動微分 14) を用いる

されたい．

ことで，最適化計算を進める上で必要な計算精度を確保する．

実際には，共分散行列である R と Q が未知である場合がつぎに， R , Q の最適化ステップについて考える．(17)式

多いが，この場合， R と Q の各成分も (17) 式の最大化に基を以下のように書き直しておく．
log p( X N , θ | Z N )づいて推定を行う．ただし，数値実験の結果， R と Q を他

  J1 ( R) / 2  J 2 (Q ) / 2  log p ( x0 )  log p (θ )  constの変数と同時に最適化することは，数値計算における収束の観点から望ましくないことがわかったため， R , Q の最適化

(22)

ただし，
N

と X N , θ の最適化は同時に実行せず，交互に繰り返し実行すJ1 ( R)   ( zk  hk )T R 1 ( zk  hk )  ( N  1)log det R

るものとする．このアプローチは FEM においても同様の理

(23)

k 0

由から採用されている 3, 4)．

N

J 2 (Q )   { xk  xˆ k (Q )}T Pk (Q ) 1{ xk  xˆ k (Q )}  log det Pk (Q ) 

X N , θ の最適化ステップにおいては，与えられた R , Q のk 1

もとで以下の問題を数値的に解く．

(24)

・最適化変数: x0 , x1 , , x N , e1 , , e N , θである． J1 ( R) を R に関して偏微分すると，次式を得る．
J1 ( R)N  R 1   ( zk  hk )( zk  hk )T  R 1  ( N  1) R 1 (25)R k 0

・最小化すべき評価関数:
N

N

k 0

k 1

J   ( zk  hk )T R 1 ( zk  hk )   (ek T ek  2log det Pk1/ 2 )(18)(25)式および J1 ( R) が R に関して凸な関数であることから，

2log p( x0 )  2log p(θ )log p( X N , θ | Z N ) を最大化する R は次のようになる．

・制約条件:

R
xk  xˆ k  Pk1/ 2ek (k  1,, N ),(19)c ( x0 , θ )  0(20)s ( x0 , θ )  0

(21)

1 N
 ( zk  hk )( zk  hk )TN  1 k 0(26)

ただし，(26)式は共分散行列の正定性を保証しないため，こ
こでも正の微小な値  を用いて以下のように修正する．

R

こ こ で ， P (k  1,, N ) は Pk の 平 方 根 行 列 で あ り ，
1/ 2
k

1 N
 ( zk  hk )( zk  hk )T   I r
N  1 k 0(27)

Pk  Pk1/ 2 ( Pk1/ 2 )T を満たす．また，(20),(21)式は x0 , θ に関するまた，上式から求まる R の非対角成分は実用上推定精度に情報を取り込むための線形あるいは非線形の制約条件である．

ほとんど影響を与えないため，零に近似して計算を簡略化す
る．一方， Q に関しては解析的に最適解を求めることが一

162

般に困難である．ここでは，下記の最適化問題を数値的に解くことで， Q を求める．
Q ( diag[Q11 , Q22 , , Qqq ])  arg min J 2 (Q )Qii  (28)(28)式の数値最適化問題のアルゴリズムには，さまざまなものが利用可能であるが，ここでは方向集合法 15) を適用する．
方向集合法は，小規模最適化問題においては比較的ロバスト
な収束性を有する．
以上に基づき，提案手法のアルゴリズムは以下のようにまとめることができる．

Step 1: X N , θ , R , Q を初期化する．

Fig. 1. Absolute mean errors of states in the reentry body problem

Step 2: R , Q を固定し，(18)-(21)式で記述される最適化問題
を解いて X N , θ を更新する．

Table 1 Statistics of the parameters in the reentry body problem

Step 3: R , Q が未知の場合， X N , θ を固定し， (27),(28) 式か(true value of b : 103 )ら R , Q を求めて対角成分だけを残す．

Mean (Absolute error)Standard
deviationb

1.016  103 ( 1.615  105 )2.621  105Q

2.643  10 ( 1.435  10 )1.840  103

R

1.007  106 ( 6.880  103 )

1.857  105b

3

5.934  105

Step 4: R , Q の更新前からの変化が十分に小さければ終了．
そうでなければ Step 2 へ戻る．

Developed
algorithm

3. 計算例

URTSS

3.1 レーダートラッキング問題
最初に，非線形システムのパラメータ推定に関する例題として有名な，レーダートラッキングによる再突入物体のパラメータ推定問題 6) を考える．本問においては，高々度から鉛直に落下する物体までの距離をレーダーで測定し，その時歴

1.012  10 ( 1.180  10 )び観測雑音の分散の真値は，それぞれ次式で与える．

Q  2.5  103[ft 2 / s 2 ], R  106 [ft 2 ](32)以上に基づき，10 通りのランダムなシミュレーションを状態方程式:

hk 1  hk  

tk 1

実行して作成した観測データに対し，提案手法を適用した．

Vdt

(29)

[b exp( h)(V  w) 2 ] dt

(30)

tk

tk

5

V0  2  104 [ft/ s] ， b  103[1 / ft] で与える．プロセス雑音およ

は以下のように記述される．

tk 1

2

る．初期状態および未知パラメータの真値は， h0  3  105 [ft] ,から物体を特徴付けるパラメータの値を推定する．システム

Vk 1  Vk  

3Unscented 変換におけるパラメータは   102 ,   2 ,   2とした．最適化計算における初期解として， x0 を真値に，b を 3  105 [1 / ft] に設定し， w  0 とおいて式(29), (30)を順次計算することで状態変数時歴 X N を与えた．分散は未知とし，観測方程式:

その初期解については Q  1[ft 2 / s 2 ], R  1[ft 2 ] とした．比較のzmk  d  (hk  href )  vk
2

2

(31)

ただし，記号の定義は以下の通りである． t : 時刻[s], h : 高度[ft]， V : 鉛直方向速度[ft/s], zm : レンジ(観測出力) [ft]， b :再突入物体を特徴付ける未知パラメータ[1/ft] ， w : 突風の速度[ft/s] (プロセス雑音)， v : レンジの雑音（観測雑音）[ft]，ため，上述の 10 通りの観測データ列に対して，拡大系に対する URTSS によるスムージングを用いて，別途パラメータ推定を行った．以降，便宜上この手法を単に URTSS と呼ぶ．

URTSS における拡大状態変数は， h,V , b で構成されるが，その初期値に関し，平均をそれぞれ 3  105 [ft] , 2  104 [ft/ s] ,3  105 [1 / ft] に ， 分 散 を そ れ ぞ れ 106 [ft 2 ], 4  104 [ft 2 / s 2 ] , : 大気密度の特性を表す定数（= 5  105 [1/ft]）， href : レー
ダーの設置高度（  105 [ft] ）， d : 物体とレーダーの水平方向の距離の差（  105 [ft] ）．文献 6 と同様に，制御入力は考えず，重力加速度の影響も無視している．なお，(30)式において，非線形なプロセス雑音が印加されると仮定しており，
この点は文献 6 の定式化と異なる．観測は 1Hz で 60 秒間にわたって行われるものとし，状態方程式(29), (30)における積
分項は，刻み幅 1[s]の 4 次の Runge-Kutta 公式で数値積分す106 [1 / ft 2 ] に設定した．また，観測雑音・プロセス雑音の分
散は真値に固定した．
Fig. 1 に高度と速度に関する絶対誤差の平均値(10 ケース分)の時歴を示す．高々度では空気抵抗が小さく，重力を無視しているために物体はほぼ等速運動をするが，時刻 10[s]付近で空気抵抗が急増し，状態方程式の非線形性が強くなる．
これに伴い，URTSS においては時刻 10[s]付近における速度の誤差が大きくなっており，非線形性の影響が強く表れてい163

Table 4 Known parameters in the research aircraft problem

Table 2 Statistics of the parameters in the reentry body problem
4(true value of b : 2.0  10 )Developed
algorithm
URTSS

Mean (Absolute error)Standard
deviationb

1.996  104 ( 4.092  107 )7.206  106Q

2.567  103 ( 6.721  101 )2.264  103

R

1.062  106 ( 6.238  104 )1.612  105b

2.015  104 ( 1.479  106 )1.003  105

Parameter(true value of b : 5.0  103 [1 / ft] )Developed
algorithm
URTSS

Standard
deviationb

4.953  103 ( 4.671  103 )2.652  104Q

2.249  10 ( 2.512  10 )1.508  103

R

9.955  105 ( 4.471  103 )

1.916  105b

31.240  103

3

2

4

4.547  10 ( 4.527  10 )0.7920 [kg/m3]

Gravity acceleration g

9.80665 [m/s2]

Mass m

7472.0 [kg]2.2847 104 [kg m2]

Moment of inertia I yy

Table 3 Statistics of the parameters in the reentry body problem

Mean (Absolute error)

Value

Air density 

Reference area S

30.097 [m2]Mean aerodynamic chord c

0.6075 [m]Reference velocity Vref

104.67 [m/s]Thrust inclination angle  T

0.0525 [rad]Moment arm of the thrust0.1603 [m]ltx sin  T  ltz cos  T

ることが推察されるが，提案手法では URTSS と比較してこの影響が弱い．
Table 1 に，10 ケースの観測データから得られた b, Q, R の
平均値と標準偏差をまとめる．推定すべき未知パラメータでFig. 2 Measured inputs of the research aircraft problem

ある b の平均値は，URTSS・提案手法共にほぼ真値の付近に収束しており，その精度は同程度であるが，標準偏差が小さい点で提案手法の方が優れている．また， R と Q の平均値と同様に，提案手法では真値付近に収束していることが確認に関しても，提案手法では真値付近に収束していることが確

できる．

認できる．

3.2 実験用航空機の空力微係数推定問題

つぎに，収束のロバスト性を確認するため， b の初期解を102 , 104 , 105 , 106 , 107 [1 / ft] の各値に変えて，それぞれのもつぎに，実験用航空機 HFB-320 の縦運動に関する空力微とで状態変数の初期解も作り直して最適化計算を行った．計

係数の推定問題を考える．飛行試験データは，文献 3 で公開算の結果，いずれのケースにおいても，上記のケースと同様されているものを使用する．試験時間は 60 秒間，観測レーに b は真値の近傍に収束し，提案手法の優れたロバスト性を

トは 10Hz である．飛行試験においては，Fig. 2 に示すとおり，制御入力の 1 つである推力を一定に保持しつつ，もう 1

確認することができた．
さらに， b の真値を 2.0  104 [1 / ft] および 5.0  103[1 / ft] につの制御入力である昇降舵角を，航空機の縦の短周期モード

設定した場合（便宜上，これらをケース B および C と呼び，およびフゴイドモードを誘起するような時歴で動かしている．

3

真値が 10 [1 / ft] の場合をケース A と呼ぶ）における推定結果 を Table 2 お よ び Table 3 に 示 す ． b の 初 期 推 定 値

システムは以下のように表わす．
状態方程式:（URTSS においては b の初期値の平均）については，ケー6Vk 1  Vk  4

ス B で 6  10 [1 / ft] ， ケ ー ス C で 1.5  10 [1 / ft] と し ，URTSS における b の初期値の分散については，ケース B で
7

2

5

tk 1

tk

[

tk 1

T cos(   T )  D
 g sin(   )  w1 ]dt (33)m k 1   k   [ q 2tk

10 [1 / ft ] ，ケース C で 10 [1 / ft ] としている．その他の

L  T sin(   T ) g
 cos(   )  w2 ]dt
mV
V

(34)

問題設定は，ケース A と同じであり，観測データもそれぞれの真値のもとで，10 通りのランダムなシミュレーション k 1   k  

tk 1

tk

を実行して作成している．Table 2 および Table 3 より，推定された b に関する平均値の誤差および標準偏差の大きさを比qk 1  qk  

tk 1

tk

べると，提案手法の方が URTSS よりも優れている．また，
ケース B, C における R と Q の平均値に関しても，ケース A

164(q  w3 ) dt

(35)

[ M  T (ltx sin  T  ltz cos  T )  w4 ] / I yy dt (36)

ここで，記号の定義は以下の通りである． V : 速度[m/s],  :迎角[rad],  : ピッチ角[rad], q : ピッチ角速度[rad/s],  e : 昇降舵角[rad], T : 推力[N], D : 抗力[N], L : 揚力[N], M : ピッチン
グモーメント [Nm], Vm : 観測速度 [m/s],  m : 観測迎角 [rad], m : 観測ピッチ角[rad], qm : 観測ピッチ角速度[rad/s], qm : 観
測ピッチ角加速度 [rad/s2], axm : 機体 X 軸方向の観測加速度[m/s2], azm : 機体 Z 軸方向の観測加速度 [m/s2]， w1 , , w4 : プ
ロセス雑音， vk (1) ,, vk (7) : 観測雑音．状態変数は V ,  ,  ,q , 制御入力は  e , T , 観測出力は添え字 m のついた 7 つの変数である．Table 4 に定数パラメータの値をまとめる．また，
D , L , M は次式で与える．
D  qd S [CD 0  CDV (Vk / Vref )  CD  ]L  qd S [CL 0  CLV (Vk / Vref )  CL ]M  qd Sc [Cm 0  CmV (Vk / Vref )  Cm(44)

 Cmq (qc ) / (2Vref )  Cm e e ]ただし， q  V / 2 であり， CD 0 , CDV , CD , CL 0 , CLV , CL ,2

Cm 0 , CmV , Cm , Cmq , Cm e は空力微係数（推定すべき未知パラメータ）である．状態方程式(33)-(36)における積分項は，刻み幅 0.1[s]の 4 次の Runge-Kutta 公式で数値積分する．

Fig. 3 Outputs in the research aircraft problem

最適化計算における初期解として，未知パラメータの値はTable 5 Estimated parameters in the research aircraft problem

すべて 1 で与え，状態変数の時歴 X N は，観測データの端点である (Vm 0 , m 0 , m 0 , qm 0 ) と (VmN , mN , mN , qmN ) を線分で内挿す

Parameter

FEM

CD 0

0.123

Developed
algorithm
0.126

CDV

0.0645

0.0667

Unscented 変換におけるパラメータは   102 ,   2 ,   1

CD

0.320

0.303

とした．

CL 0

0.0929

0.0946

観測出力の時歴と推定された状態変数に基づく出力の時歴

CLV

0.149

0.153

を Fig. 3 に示す．全体的に出力はよく一致していることがわ

CL

4.328

4.285

かる．Table 5 に未知パラメータの推定結果をまとめる．

Cm 0

0.112

0.119

Table 5 における”FEM”の項目は，FEM による推定結果 3) を

CmV

0.0039

引用したものである．ただし，FEM の最適化計算においてCm0.968

1.024は，初期解が表中に示した最終的な解の近傍に設定されておCmq

34.710

35.188

り，FEM において提案手法と同様に未知パラメータの初期Cm e1.529

1.581

解をすべて 1 として最適化計算を行うと，計算を収束させる

ることによって与えた．一方，プロセス雑音・観測雑音の共分 散 の 初 期 解 は Q  I 4 , R  104  I 7 に 設 定 し た ． ま た ，ことができなかった．一方，提案手法では，未知パラメータの初期解をすべて 1 とした場合，FEM で収束した初期解と同じ値にした場合のいずれにおいても，収束させることがで観測方程式:

きた．未知パラメータの真値は不明であるが，Table 5 より，

Vmk  Vk  vk (1)(37)

 mk   k  vk (2)(38)

FEM と同じ程度の推定精度が提案手法によって達成できて mk   k  vk (3)(39)

いるものと解釈できる．

qmk  qk  vk (4)(40)

qmk  [ M k  Tk (ltx sin  T  ltz cos  T )] / I yy  vk (5)(41)

axmk  ( Lk sin  k  Dk cos  k  Tk cos  T ) / m  vk (6)(42)

azmk  ( Lk cos  k  Dk sin  k  Tk sin  T ) / m  vk (7)(43)

いずれも FEM と提案手法の推定結果が大体一致しており，

4. おわりに
本稿では，観測雑音とプロセス雑音の両方に対応した，移動体のパラメータ推定手法について提案した．提案手法は，
非線形計画法と Unscented 変換に基づいており，非線形なプ
ロセス雑音および雑音共分散の推定を含む一般的なクラスの165

問題に適用することが可能である．例題への適用を通じ，提[著 者 紹 介]

案手法が，安定した収束性と，従来手法と同等あるいはそれ以上の推定精度を有することを示した．

横

山

信

宏 （正会員）提案手法では， 収束の安定性を確保するために，状態変2005 年東京大学大学院工学系研究科博

数・未知パラメータに関する最適化と，雑音の共分散に関す

士課程航空宇宙工学専攻修了．同年宇宙航る最適化を交互に行っており，アルゴリズムがやや複雑であ空研究開発機構宇宙航空プロジェクト研究

る．推定精度や収束の安定性を確保しつつ，これらの最適化

員． 2007 年防衛大学校システム工学群航を同時に行えるような，より簡素化された手法を検討するこ

空宇宙工学科助教， 2011 年同講師，現在とが，今後の課題の 1 つである．

に至る．博士(工学)．航空機の誘導制御・システム最適化等の研究に従事．日本航空
宇宙学会，AIAA, IEEE の会員．

参考文献
1） L. Ljung: System identification – theory for the user, second
edition, Prentice Hall (1999)
2） 足立: MATLAB による制御のためのシステム同定，東京
電機大学出版局 (1996)
3） R. V. Jategaonkar: Flight vehicle system identification - a time
domain methodology, AIAA (2006)
4） R. V. Jategaonkar and E. Plaetschke: Algorithms for aircraft
parameter estimation accounting for process and measurement
noise, J. of Aircraft, 26-4, 360/372 (1989)
5） G. Chowdhary and R. Jategaonkar: Aerodynamic parameter
estimation from flight data applying extended and unscented
Kalman filter, AIAA Paper 2006-6146 (2006)
6） S., Julier, J. Uhlmann, and H. F. Durrant-Whyte: A new
method for the nonlinear transformation of means andcovariances in filters and estimators, IEEE Trans. on
Automatic Control, AC-45-3, 477/482 (2000)
7） S. Julier and J. K. Uhlmann,: Unscented filtering for nonlinear
estimation, Proceedings of the IEEE, 92-3, 401/422 (2004)
8） S. Särkkä: Unscented Rauch-Tung-Striebel smoother, IEEE
Trans. on Automatic Control, AC-53-3, 845/849 (2008)9） N. Yokoyama: Parameter estimation of aircraft dynamics via
unscented smoother with expectation-maximization algorithm,
J. of Guidance, Control, and Dynamics, 34-2, 426/436 (2011)10) B. M. Bell, J. V. Burke, and G. Pillonetto: An inequality
constrained nonlinear Kalman-Bucy smoother by interior
point likelihood maximization, Automatica, 45-1, 25/33(2009)
11) 加藤: 最適制御入門，126/129，東京大学出版会(1987)
12) J. T. Betts: Practical methods for optimal control using
nonlinear programming, SIAM (2001)
13) N. Yokoyama, S. Suzuki, and T. Tsuchiya: Convergenceacceleration of direct trajectory optimization using novel
Hessian calculation methods, J. of optimization theory andapplications, 136-3, 297/319 (2008)14) 久保田，伊理: アルゴリズムの自動微分と応用，コロナ
社 (1998)
15) R. P. Brent: Algorithms for minimization without derivatives.
Chapter 7, Prentice-Hall (1973)166