(57)【要約】【課題】非線形ダイナミクスを有する制御対象についてモデル化を容易にし、求められたモデルを用いて非線形ダイナミクスを有する制御対象の制御を可能とすることである。【解決手段】車両制御システム３０は、大別して２つの機能及び構成に分けて考えることができる。その１つは、モデル化コンピュータ４０によって実現されるもので、車両１０の制御特性を模擬した模擬車両２０から、制御に関するパラメータを多数収集し、それを用いて、非線形状態変数ｘと線形入力変数ｕについて2007141121.tif 000059の形式で車両１０のモデル化を行う。その２は、車両制御コンピュータ７０によって実現されるもので、モデル化コンピュータ４０によって得られた車両１０のモデルを用い、実際の車両１０について、姿勢制御が行われる。【選択図】図１
2007141121.tif 000002
【特許請求の範囲】
【請求項１】
非線形ダイナミクスを有する制御対象について、
制御対象の観察可能な非線形状態変数をｘとし、
制御対象に対する線形の入力変数をｕとし、
ｄｘ／ｄｔを算出対象とするときに、制御対象の状態方程式を
【数１】
2007141121.tif 000046
の形式にモデル化するモデル化装置であって、
制御対象について、ｄｘ／ｄｔ、ｘ、ｕを一組としたデータセットを複数収集し、その一部を複数のトレーニングセットとし、残りを複数のチェッキングセットとするデータセット収集手段と、
各ｕｉ
【数２】
2007141121.tif 000047
についてはたかだか１次形式とし、各ｘｉ
【数３】
2007141121.tif 000048
については高次形式を含み得るものとして、ｘ、ｕについての複数の部分多項式を生成する部分多項式生成手段と、
複数のトレーニングセットのデータに基き、複数の部分多項式のそれぞれの多項式係数を算出する多項式係数算出手段と、
算出された多項式係数を有する複数の部分多項式のそれぞれについてｄｘ／ｄｔの代表性を評価する手段であって、チェッキングセットの１のデータセットのｘ、ｕを部分多項式に適用して求められる値とそれに対応するｄｘ／ｄｔとの値との間の誤差を求め、これを複数のチェッキングセット全体について統計的に行うことで代表性を評価する評価手段と、
各部分多項式のそれぞれについて算出された代表性評価を比較し、制御対象の状態方程式として最も代表性の高い部分多項式を決定し、これを
【数４】
2007141121.tif 000049
の形式で表現するモデル化手段と、
を含むことを特徴とする非線形ダイナミクスを有する制御対象のためのモデル化装置。
【請求項２】
制御対象の観察可能な非線形状態変数をｘとし、制御対象に対する線形の入力変数をｕとし、ｄｘ／ｄｔを算出対象とする非線形ダイナミクスを有する制御対象を、任意に設定される目標制御状態に制御する制御装置であって、
請求項１に記載の非線形ダイナミクスを有する制御対象のためのモデル化装置によって求められる
【数５】
2007141121.tif 000050
のモデルのｆ（ｘ）とｇ（ｘ）の関数形と、データセット収集工程で収集されたデータセットの範囲に基いて規定されるｘ、ｕの制御適用範囲を表す集合であるＰｘ、Ｐｕと、を記憶する記憶部と、
記憶部から読み出されたｆ（ｘ）とｇ（ｘ）の関数形を用いて
【数６】
2007141121.tif 000051
の関係式に基いて制御対象のモデルを生成するモデル生成手段と、
制御対象から任意の時間ｋにおけるｘ（ｋ）を取得する観察手段と、
記憶部からＰｕを読み出してその要素である各ｕのそれぞれについて、
【数７】
2007141121.tif 000052
の関係式に基いて、任意の時間ｋの後の時間（ｋ＋１）におけるｘ（ｋ＋１）を算出し、算出されたｘ（ｋ＋１）が記憶部から読み出されたＰｘの要素であるか否かを判断する判断手段と、
を含み、
ｘ（ｋ＋１）がＰｘの要素であり、ｕがＰｕの要素であるときに、そのｕを制御対象に対する入力変数ｕ（ｋ）として用いることを特徴とする非線形ダイナミクスを有する制御対象のための制御装置。
【請求項３】
請求項２に記載の非線形ダイナミクスを有する制御対象のための制御装置において、さらに、
算出されたｘ（ｋ＋１）を目標制御状態に対して評価し、最も評価の高いｘ（ｋ＋１）に対応するｕをｕ（ｋ）として制御対象に出力する評価手段を含むことを特徴とする非線形ダイナミクスを有する制御対象のための制御装置。
【請求項４】
請求項３に記載の非線形ダイナミクスを有する制御対象のための制御装置において、
評価手段は、算出されたｘ（ｋ＋１）を目標制御状態に対し、二次計画法に基いて評価することを特徴とする非線形ダイナミクスを有する制御対象のための制御装置。
【請求項５】
車両の長手方向に沿った速度成分である縦方向速度ｘ１を検出する手段と、
車両の幅方向に沿った速度成分である横方向速度ｘ２を検出する手段と、
車両の長手方向を基準として車両の進行方向の角度ｕを変更するステアリング手段と、
を有する車両に対し、任意に設定された目標制御状態になるように車両の進行方向姿勢を制御する車両のための制御装置であって、
ｘ１、ｘ２を車両の観察可能な非線形状態変数とし、角度ｕを車両に対する線形の入力変数とし、ｄｘ１／ｄｔ、ｄｘ１／ｄｔを算出対象として、請求項１に記載の非線形ダイナミクスを有する制御対象のモデル化装置によって求められる
【数８】
2007141121.tif 000053
のモデルを、ｄｘ１／ｄｔ＝ｆ１（ｘ１，ｘ２）＋ｇ１（ｘ１，ｘ２）ｕ及びｄｘ２／ｄｔ＝ｆ２（ｘ１，ｘ２）＋ｇ２（ｘ１，ｘ２）ｕとして、ｆ１（ｘ１，ｘ２）、ｇ１（ｘ１，ｘ２）、ｆ２（ｘ１，ｘ２）、ｇ２（ｘ１，ｘ２）の各関数形と、データセット収集工程で収集されたデータセットの範囲に基いて規定されるｘ１、ｘ２、ｕの制御適用範囲を表す集合であるＰｘ１、Ｐｘ２、Ｐｕと、を記憶する記憶部と、
記憶部から読み出されたｆ１（ｘ１，ｘ２）、ｇ１（ｘ１，ｘ２）、ｆ２（ｘ１，ｘ２）、ｇ２（ｘ１，ｘ２）の各関数形を用いてｄｘ１／ｄｔ＝ｆ１（ｘ１，ｘ２）＋ｇ１（ｘ１，ｘ２）ｕ及びｄｘ２／ｄｔ＝ｆ２（ｘ１，ｘ２）＋ｇ２（ｘ１，ｘ２）ｕの関係式に基いて車両のモデルを生成するモデル生成手段と、
車両から任意の時間ｋにおけるｘ１（ｋ）、ｘ２（ｋ）を取得する観察手段と、
記憶部からＰｕを読み出してその要素である各ｕのそれぞれについて、ｄｘ１／ｄｔ＝ｆ１（ｘ１，ｘ２）＋ｇ１（ｘ１，ｘ２）ｕ及びｄｘ２／ｄｔ＝ｆ２（ｘ１，ｘ２）＋ｇ２（ｘ１，ｘ２）ｕの関係式に基いて、任意の時間ｋの後の時間（ｋ＋１）におけるｘ１（ｋ＋１）、ｘ２（ｋ＋１）を算出し、算出されたｘ１（ｋ＋１）、ｘ２（ｋ＋１）が記憶部から読み出されたＰｘ１、Ｐｘ２の要素であるか否かを判断する判断手段と、
ｘ１（ｋ＋１）がＰｘ１の要素であり、ｘ２（ｋ＋１）がＰｘ２の要素であり、ｕがＰｕの要素であるときに、そのｕを制御対象である車両に対する入力変数ｕ（ｋ）として用いることを特徴とする非線形ダイナミクスを有する車両のための制御装置。
【請求項６】
請求項５に記載の非線形ダイナミクスを有する車両のための制御装置において、さらに、
算出されたｘ１（ｋ＋１）、ｘ２（ｋ＋１）を目標制御状態に対して評価し、最も評価の高いｘ１（ｋ＋１）、ｘ２（ｋ＋１）に対応するｕをｕ（ｋ）として制御対象である車両に出力する評価手段を含むことを特徴とする非線形ダイナミクスを有する車両のための制御装置。
【請求項７】
請求項６に記載の非線形ダイナミクスを有する車両のための制御装置において、
目標制御状態は、駐車操作における目標縦方向速度ｘ１ｄと目標横方向速度ｘ２ｄであることを特徴とする非線形ダイナミクスを有する車両のための制御装置。
【請求項８】
請求項６に記載の非線形ダイナミクスを有する車両のための制御装置において、
モデル生成手段は、車両のモデルを、連続的値を有する状態変数及び入力変数と、離散的値を有する状態変数及び遷移変数とで表現されるハイブリッドモデルとして生成することを特徴とする非線形ダイナミクスを有する車両のための制御装置。
【請求項９】
非線形ダイナミクスを有する制御対象について、
制御対象の観察可能な非線形状態変数をｘとし、
制御対象に対する線形の入力変数をｕとし、
ｄｘ／ｄｔを算出対象とするときに、制御対象の状態方程式を
【数９】
2007141121.tif 000054
の形式にモデル化するモデル化装置上で実行されるモデル化プログラムであって、
制御対象について、ｄｘ／ｄｔ、ｘ、ｕを一組としたデータセットを複数収集し、その一部を複数のトレーニングセットとし、残りを複数のチェッキングセットとするデータセット収集処理手順と、
ｕについてはたかだか１次形式とし、ｘについては高次形式を含み得るものとして、ｘ、ｕについての複数の部分多項式を生成する部分多項式生成処理手順と、
複数のトレーニングセットのデータに基き、複数の部分多項式のそれぞれの多項式係数を算出する多項式係数算出処理手順と、
算出された多項式係数を有する複数の部分多項式のそれぞれについてｄｘ／ｄｔの代表性を評価する手段であって、チェッキングセットの１のデータセットのｘ、ｕを部分多項式に適用して求められる値とそれに対応するｄｘ／ｄｔとの値との間の誤差を求め、これを複数のチェッキングセット全体について統計的に行うことで代表性を評価する評価処理手順と、
各部分多項式のそれぞれについて算出された代表性評価を比較し、制御対象の状態方程式として最も代表性の高い部分多項式を決定し、これを
【数１０】
2007141121.tif 000055
の形式で表現するモデル化処理手順と、
を実行することを特徴とする非線形ダイナミクスを有する制御対象のためのモデル化プログラム。
【請求項１０】
制御対象の観察可能な非線形状態変数をｘとし、制御対象に対する線形の入力変数をｕとし、ｄｘ／ｄｔを算出対象とする非線形ダイナミクスを有する制御対象を、任意に設定される目標制御状態に制御する制御装置上で実行される制御プログラムであって、
請求項１に記載の非線形ダイナミクスを有する制御対象のためのモデル化装置によって求められる
【数１１】
2007141121.tif 000056
のモデルのｆ（ｘ）とｇ（ｘ）の関数形と、データセット収集工程で収集されたデータセットの範囲に基いて規定されるｘ、ｕの制御適用範囲を表す集合であるＰｘ、Ｐｕと、を記憶する記憶部から、ｆ（ｘ）とｇ（ｘ）の関数形とＰｘ、Ｐｕとを読み出す読出処理手順と、
記憶部から読み出されたｆ（ｘ）とｇ（ｘ）の関数形を用いて
【数１２】
2007141121.tif 000057
の関係式に基いて制御対象のモデルを生成するモデル生成処理手順と、
制御対象から任意の時間ｋにおけるｘ（ｋ）を取得する観察処理手順と、
記憶部からＰｕを読み出してその要素である各ｕのそれぞれについて、
【数１３】
2007141121.tif 000058
の関係式に基いて、任意の時間ｋの後の時間（ｋ＋１）におけるｘ（ｋ＋１）を算出し、算出されたｘ（ｋ＋１）が記憶部から読み出されたＰｘの要素であるか否かを判断する判断処理手順と、
を実行し、
ｘ（ｋ＋１）がＰｘの要素であり、ｕがＰｕの要素であるときに、そのｕを制御対象に対する入力変数ｕ（ｋ）として用いることを特徴とする非線形ダイナミクスを有する制御対象のための制御プログラム。
【請求項１１】
請求項１０に記載の非線形ダイナミクスを有する制御対象のための制御プログラムにおいて、さらに、
算出されたｘ（ｋ＋１）を目標制御状態に対して評価し、最も評価の高いｘ（ｋ＋１）に対応するｕをｕ（ｋ）として制御対象に出力する評価処理手順を含むことを特徴とする非線形ダイナミクスを有する制御対象のための制御プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、非線形ダイナミクスを有する制御対象のモデル化装置及びその上で実行されるプログラムと、モデル化装置等によって求められるモデルを用いた非線形ダイナミクスを有する制御対象のための制御装置、非線形ダイナミクスを有する制御対象のモデル化プログラム、及び非線形ダイナミクスを有する車両のための制御装置に係り、特に、制御対象の観察可能な非線形状態変数をｘとし、制御対象に対する線形の入力変数をｕとし、ｄｘ／ｄｔを算出対象とする非線形ダイナミクスを有する制御対象について、そのモデル化装置及びその上で実行されるプログラムと、モデル化装置等によって求められるモデルを用いた非線形ダイナミクスを有する制御対象のための制御装置、非線形ダイナミクスを有する制御対象のモデル化プログラム、及び非線形ダイナミクスを有する車両のための制御装置に関する。
【背景技術】
【０００２】
車両の操作と走行状態のように、非線形のダイナミクスを有するシステムが数多くあるが、これらのシステムを制御するコントローラ、すなわち制御装置を設計することは容易ではない。すなわち、従来の制御システムの構築を行う際のモデリング手法は、力学的な方程式表現を用いる解析的な手法であるが、この解析的手法を用いて非線形系に対するコントローラを設計しようとすると、次のような問題に直面する。第１に、非線形ダイナミクスをモデリングすることが困難で、常に精度のよい解析的モデルを得られる保証がないこと、第２に、仮に精度のよいモデリングが求められたとしても、そのモデルはコントローラの構築に適した形ではなく、コントローラの作成に容易な形への変換が困難である。
【０００３】
これらの問題に対し、１９６０年代に、アレクセイ・イヴァクネンコによりＧＭＤＨ（Ｇｒｏｕｐ  Ｍｅｔｈｏｄ  ｏｆ  Ｄａｔａ  Ｈａｎｄｌｉｎｇ）が提案された。このアルゴリズムは、未知のダイナミクスを自己組織化手法により多項式モデルに近似するものである。非特許文献１には、このＧＭＤＨの１つの適用例として、人間の運転技能が多項式表現によってモデル化されることが述べられており、また、非特許文献２には、他の適用例としては、気象のモデル化が報告されている。
【０００４】
【非特許文献１】Jong-Hae Kim, Yoshimichi Matsui, Shoichiro Hayakawa, Tatsuya Suzuki, Shigeru Okuma and Nuio Tsuchida, 「Acquisition and Modeling of Driving Skills by using Three Dimensional Driving Simulator」, IEICE Trans. Fundamentals, Vol.E88-A, No.3, Mar. 2005, p 770-778
【非特許文献２】William M. Lebow, Ramam K. Mehra and Paul M. Toldalagi, 「Forecasting Application of GMDH in Agriculture and Meteorological Time Series」,SELF-ORGANIZING METHODS IN MODELING, Dekker, vol. 54, p121-147
【発明の開示】
【発明が解決しようとする課題】
【０００５】
このように、ＧＭＤＨの手法は、非線形ダイナミクスを有するシステムにとって注目に値するものである。しかし、ＧＭＤＨは、例えば状態方程式の入力変数のように線形的に展開する変数と、例えば状態方程式の状態変数のように非線形的に展開する変数の指定がなく、すべての観測可能な変数に対して同じ部分多項式をもって展開するので、ＧＭＤＨによって得られる多項式モデルは非常に高次のものになりやすく、伝達関数や状態方程式に代表される現代制御理論による制御系設計の枠組の中での応用が困難である。したがって、これらは、特定の変数に対する解析、推定、認識、クラスタリング等、信号処理的な応用のみについて報告されており、制御系に直接応用された例はない。
【０００６】
本発明は、非線形ダイナミクスを有する制御対象についてモデル化を容易にする非線形ダイナミクスを有する制御対象のモデル化装置及びその上で実行されるプログラムを提供することである。また、求められたモデルを用いて非線形ダイナミクスを有する制御対象の制御を可能とする非線形ダイナミクスを有する制御対象のための制御装置、非線形ダイナミクスを有する制御対象のモデル化プログラム、及び非線形ダイナミクスを有する車両の制御装置を提供することである。
【課題を解決するための手段】
【０００７】
請求項１における本発明に係る非線形ダイナミクスを有する制御対象のモデル化装置は、非線形ダイナミクスを有する制御対象について、制御対象の観察可能な非線形状態変数をｘとし、制御対象に対する線形の入力変数をｕとし、ｄｘ／ｄｔを算出対象とするときに、制御対象の状態方程式を式（１）
【数１】
2007141121.tif 000003
の形式にモデル化するモデル化装置であって、制御対象について、ｄｘ／ｄｔ、ｘ、ｕを一組としたデータセットを複数収集し、その一部を複数のトレーニングセットとし、残りを複数のチェッキングセットとするデータセット収集手段と、式（２）に示されるｕについてはたかだか１次形式とし、
【数２】
2007141121.tif 000004
式（３）に示されるｘについては高次形式を含み得るものとして、
【数３】
2007141121.tif 000005
ｘ、ｕについての複数の部分多項式を生成する部分多項式生成手段と、複数のトレーニングセットのデータに基き、複数の部分多項式のそれぞれの多項式係数を算出する多項式係数算出手段と、算出された多項式係数を有する複数の部分多項式のそれぞれについてｄｘ／ｄｔの代表性を評価する手段であって、チェッキングセットの１のデータセットのｘ、ｕを部分多項式に適用して求められる値とそれに対応するｄｘ／ｄｔとの値との間の誤差を求め、これを複数のチェッキングセット全体について統計的に行うことで代表性を評価する評価手段と、各部分多項式のそれぞれについて算出された代表性評価を比較し、制御対象の状態方程式として最も代表性の高い部分多項式を決定し、これを式（４）
【数４】
2007141121.tif 000006
の形式で表現するモデル化手段と、を含むことを特徴とする。
【０００８】
また、本発明に係る非線形ダイナミクスを有する制御対象のための制御装置は、制御対象の観察可能な非線形状態変数をｘとし、制御対象に対する線形の入力変数をｕとし、ｄｘ／ｄｔを算出対象とする非線形ダイナミクスを有する制御対象を、任意に設定される目標制御状態に制御する制御装置であって、請求項１に記載の非線形ダイナミクスを有する制御対象のモデル化装置によって求められ、式（５）で示される
【数５】
2007141121.tif 000007
のモデルのｆ（ｘ）とｇ（ｘ）の関数形と、データセット収集工程で収集されたデータセットの範囲に基いて規定されるｘ、ｕの制御適用範囲を表す集合であるＰｘ、Ｐｕと、を記憶する記憶部と、記憶部から読み出されたｆ（ｘ）とｇ（ｘ）の関数形を用いて式（６）に示される
【数６】
2007141121.tif 000008
の関係式に基いて制御対象のモデルを生成するモデル生成手段と、制御対象から任意の時間ｋにおけるｘ（ｋ）を取得する観察手段と、記憶部からＰｕを読み出してその要素である各ｕのそれぞれについて、式（７）に示される
【数７】
2007141121.tif 000009
の関係式に基いて、任意の時間ｋの後の時間（ｋ＋１）におけるｘ（ｋ＋１）を算出し、算出されたｘ（ｋ＋１）が記憶部から読み出されたＰｘの要素であるか否かを判断する判断手段と、を含み、ｘ（ｋ＋１）がＰｘの要素であり、ｕがＰｕの要素であるときに、そのｕを制御対象に対する入力変数ｕ（ｋ）として用いることを特徴とする。
【０００９】
また、本発明に係る非線形ダイナミクスを有する制御対象のための制御装置において、さらに、算出されたｘ（ｋ＋１）を目標制御状態に対して評価し、最も評価の高いｘ（ｋ＋１）に対応するｕをｕ（ｋ）として制御対象に出力する評価手段を含むことが好ましい。
【００１０】
また、評価手段は、算出されたｘ（ｋ＋１）を目標制御状態に対し、二次計画法に基いて評価することが好ましい。
【００１１】
また、本発明に係る非線形ダイナミクスを有する車両のための制御装置は、車両の長手方向に沿った速度成分である縦方向速度ｘ１を検出する手段と、車両の幅方向に沿った速度成分である横方向速度ｘ２を検出する手段と、車両の長手方向を基準として車両の進行方向の角度ｕを変更するステアリング手段と、を有する車両に対し、任意に設定された目標制御状態になるように車両の進行方向姿勢を制御する車両の制御装置であって、ｘ１、ｘ２を車両の観察可能な非線形状態変数とし、角度ｕを車両に対する線形の入力変数とし、ｄｘ１／ｄｔ、ｄｘ１／ｄｔを算出対象として、請求項１に記載の非線形ダイナミクスを有する制御対象のモデル化装置によって求められ、式（８）で示される
【数８】
2007141121.tif 000010
のモデルを、ｄｘ１／ｄｔ＝ｆ１（ｘ１，ｘ２）＋ｇ１（ｘ１，ｘ２）ｕ及びｄｘ２／ｄｔ＝ｆ２（ｘ１，ｘ２）＋ｇ２（ｘ１，ｘ２）ｕとして、ｆ１（ｘ１，ｘ２）、ｇ１（ｘ１，ｘ２）、ｆ２（ｘ１，ｘ２）、ｇ２（ｘ１，ｘ２）の各関数形と、データセット収集工程で収集されたデータセットの範囲に基いて規定されるｘ１、ｘ２、ｕの制御適用範囲を表す集合であるＰｘ１、Ｐｘ２、Ｐｕと、を記憶する記憶部と、記憶部から読み出されたｆ１（ｘ１，ｘ２）、ｇ１（ｘ１，ｘ２）、ｆ２（ｘ１，ｘ２）、ｇ２（ｘ１，ｘ２）の各関数形を用いてｄｘ１／ｄｔ＝ｆ１（ｘ１，ｘ２）＋ｇ１（ｘ１，ｘ２）ｕ及びｄｘ２／ｄｔ＝ｆ２（ｘ１，ｘ２）＋ｇ２（ｘ１，ｘ２）ｕの関係式に基いて車両のモデルを生成するモデル生成手段と、車両から任意の時間ｋにおけるｘ１（ｋ）、ｘ２（ｋ）を取得する観察手段と、記憶部からＰｕを読み出してその要素である各ｕのそれぞれについて、ｄｘ１／ｄｔ＝ｆ１（ｘ１，ｘ２）＋ｇ１（ｘ１，ｘ２）ｕ及びｄｘ２／ｄｔ＝ｆ２（ｘ１，ｘ２）＋ｇ２（ｘ１，ｘ２）ｕの関係式に基いて、任意の時間ｋの後の時間（ｋ＋１）におけるｘ１（ｋ＋１）、ｘ２（ｋ＋１）を算出し、算出されたｘ１（ｋ＋１）、ｘ２（ｋ＋１）が記憶部から読み出されたＰｘ１、Ｐｘ２の要素であるか否かを判断する判断手段と、ｘ１（ｋ＋１）がＰｘ１の要素であり、ｘ２（ｋ＋１）がＰｘ２の要素であり、ｕがＰｕの要素であるときに、そのｕを制御対象である車両に対する入力変数ｕ（ｋ）として用いることを特徴とする。
【００１２】
また、本発明に係る非線形ダイナミクスを有する車両のための制御装置において、さらに、算出されたｘ１（ｋ＋１）、ｘ２（ｋ＋１）を目標制御状態に対して評価し、最も評価の高いｘ１（ｋ＋１）、ｘ２（ｋ＋１）に対応するｕをｕ（ｋ）として制御対象である車両に出力する評価手段を含むことが好ましい。
【００１３】
また、本発明に係る非線形ダイナミクスを有する車両のための制御装置において、目標制御状態は、駐車操作における目標縦方向速度ｘ１ｄと目標横方向速度ｘ２ｄであることが好ましい。
【００１４】
また、モデル生成手段は、車両のモデルを、連続的値を有する状態変数及び入力変数と、離散的値を有する状態変数及び遷移変数とで表現されるハイブリッドモデルとして生成することが好ましい。
【００１５】
また、本発明に係る非線形ダイナミクスを有する制御対象のモデル化プログラムは、非線形ダイナミクスを有する制御対象について、制御対象の観察可能な非線形状態変数をｘとし、制御対象に対する線形の入力変数をｕとし、ｄｘ／ｄｔを算出対象とするときに、制御対象の状態方程式を式（９）で示される
【数９】
2007141121.tif 000011
の形式にモデル化するモデル化装置上で実行されるモデル化プログラムであって、制御対象について、ｄｘ／ｄｔ、ｘ、ｕを一組としたデータセットを複数収集し、その一部を複数のトレーニングセットとし、残りを複数のチェッキングセットとするデータセット収集処理手順と、ｕについてはたかだか１次形式とし、ｘについては高次形式を含み得るものとして、ｘ、ｕについての複数の部分多項式を生成する部分多項式生成処理手順と、複数のトレーニングセットのデータに基き、複数の部分多項式のそれぞれの多項式係数を算出する多項式係数算出処理手順と、算出された多項式係数を有する複数の部分多項式のそれぞれについてｄｘ／ｄｔの代表性を評価する手段であって、チェッキングセットの１のデータセットのｘ、ｕを部分多項式に適用して求められる値とそれに対応するｄｘ／ｄｔとの値との間の誤差を求め、これを複数のチェッキングセット全体について統計的に行うことで代表性を評価する評価処理手順と、各部分多項式のそれぞれについて算出された代表性評価を比較し、制御対象の状態方程式として最も代表性の高い部分多項式を決定し、これを式（１０）で示される
【数１０】
2007141121.tif 000012
の形式で表現するモデル化処理手順と、を実行することを特徴とする。
【００１６】
また、本発明に係る非線形ダイナミクスを有する制御対象のための制御プログラムは、制御対象の観察可能な非線形状態変数をｘとし、制御対象に対する線形の入力変数をｕとし、ｄｘ／ｄｔを算出対象とする非線形ダイナミクスを有する制御対象を、任意に設定される目標制御状態に制御する制御装置上で実行される制御プログラムであって、請求項１に記載の非線形ダイナミクスを有する制御対象のモデル化装置によって求められ、式（１１）で示される
【数１１】
2007141121.tif 000013
のモデルのｆ（ｘ）とｇ（ｘ）の関数形と、データセット収集工程で収集されたデータセットの範囲に基いて規定されるｘ、ｕの制御適用範囲を表す集合であるＰｘ、Ｐｕと、を記憶する記憶部から、ｆ（ｘ）とｇ（ｘ）の関数形とＰｘ、Ｐｕとを読み出す読出処理手順と、記憶部から読み出されたｆ（ｘ）とｇ（ｘ）の関数形を用いて式（１２）で示される
【数１２】
2007141121.tif 000014
の関係式に基いて制御対象のモデルを生成するモデル生成処理手順と、制御対象から任意の時間ｋにおけるｘ（ｋ）を取得する観察処理手順と、記憶部からＰｕを読み出してその要素である各ｕのそれぞれについて、式（１３）で示される
【数１３】
2007141121.tif 000015
の関係式に基いて、任意の時間ｋの後の時間（ｋ＋１）におけるｘ（ｋ＋１）を算出し、算出されたｘ（ｋ＋１）が記憶部から読み出されたＰｘの要素であるか否かを判断する判断処理手順と、を実行し、ｘ（ｋ＋１）がＰｘの要素であり、ｕがＰｕの要素であるときに、そのｕを制御対象に対する入力変数ｕ（ｋ）として用いることを特徴とする。
【００１７】
また、本発明に係る非線形ダイナミクスを有する制御対象のための制御プログラムにおいて、さらに、算出されたｘ（ｋ＋１）を目標制御状態に対して評価し、最も評価の高いｘ（ｋ＋１）に対応するｕをｕ（ｋ）として制御対象に出力する評価処理手順を含むことが好ましい。
【発明の効果】
【００１８】
上記構成の少なくとも１つにより、制御対象の状態方程式を式（１４）で示される
【数１４】
2007141121.tif 000016
の形式にモデル化する。この形式は、従来の制御理論の枠組みの中で容易に解析でき、又制御装置として容易に実現できるものである。
【００１９】
この式（１５）で示される
【数１５】
2007141121.tif 000017
の形式にモデル化するために、ｘ、ｕについての複数の部分多項式を生成する。これらの部分多項式はいずれも、ｕについてはたかだか１次形式とし、ｘについては高次形式を含み得るものとする。ここで、ｕについてはたかだか１次形式とするとは、ｕの多項式展開が線形的であることを意味する。これに対しｘについては高次形式を含み得るものとするとは、ｘの多項式展開が非線形的であってもよいことを意味する。そして、制御対象からｄｘ／ｄｔ、ｘ、ｕを一組としたデータセットを複数収集する。これらのデータセットは、２つの目的に用いられる。１つ目の目的は、トレーニングセットとして各部分多項式の多項式係数を定めるために用いられ、２つ目の目的は、チェッキングセットとして多項式係数が定められた各部分多項式のそれぞれについてｄｘ／ｄｔの代表性を評価するために用いられる。そして、最もｄｘ／ｄｔの代表性の高い部分多項式が決定される。この決定された部分多項式は、上記のようにｕについて線形的に展開され、ｘについて非線形的に展開されているので、これは式（１６）で示される
【数１６】
2007141121.tif 000018
の形式に整理でき、かくてこの形式でのモデル化が実現される。
【００２０】
また、上記構成の少なくとも１つにより、上記のようにしてモデル化された式（１７）で示される
【数１７】
2007141121.tif 000019
のモデルのｆ（ｘ）とｇ（ｘ）の関数形を用いて、制御対象のシステムのモデルが生成される。そして、そのモデルの下で、データセット収集工程で収集されたデータセットの範囲に基いて規定されるｘ、ｕの制御適用範囲を表す集合であるＰｘ、Ｐｕを用いて、制御対象の制御が行われる。ここでＰｘ、Ｐｕを用いることとしたのは、Ｐｘの要素でないｘ、Ｐｕの要素でないｕを以後の解析に用いると、解析結果と実際の結果との乖離が生じうるので、以後の解析にはＰｘの要素であるｘ、Ｐｕの要素であるｕを用いることとするためである。
【００２１】
制御対象の制御を行う制御装置において、まず、制御対象から任意の時間ｋにおけるｘ（ｋ）が取得される。そして、Ｐｕの要素であるｕを変数として、式（１８）で示される
【数１８】
2007141121.tif 000020
の関係式に基いて、任意の時間ｋの後の時間（ｋ＋１）におけるｘ（ｋ＋１）を算出する。ここで算出されたｘ（ｋ＋１）がＰｘの要素であれば、そのｕは、制御対象に対する入力変数ｕ（ｋ）として用いることができることになる。そこでそのｕ（ｋ）を制御対象に入力し、制御対象が応答した結果を実際のｘ（ｋ＋１）として取得する。これを繰り返せば、制御対象に次々と入力変数ｕを与えることができる。入力変数ｕに対する制御対象の応答は、トレーニングセット及びチェッキングセットによってモデル化されているものであるから、Ｐｘの要素であるｘ、Ｐｕの要素であるｕを用いることで、制御対象の応答がモデル化された結果に近くなることが期待できる。このようにして、制御対象の制御が自動的に行われる。
【００２２】
ここで、Ｐｕの要素であるｕが複数あることがあり、したがって、ｘ（ｋ＋１）もまた複数ありえる。そこで、算出されたｘ（ｋ＋１）を目標制御状態に対して評価し、最も評価の高いｘ（ｋ＋１）に対応するｕをｕ（ｋ）として制御対象に出力することとする。これにより、かかるｕ（ｋ）が入力された制御対象は目標制御状態に対し近い結果を応答することが期待される。
【００２３】
また、評価手段は、算出されたｘ（ｋ＋１）を目標制御状態に対し、二次計画法に基いて評価することとするので、よく知られた方法で評価することができる。
【００２４】
また、上記構成の少なくとも１つにより、上記のようにしてモデル化された式（１９）で示される
【数１９】
2007141121.tif 000021
のモデルを車両に適用し、車両の進行方向姿勢を制御することができる。すなわち、観察可能な状態変数として車両の縦方向速度ｘ１と車両の横方向速度ｘ２を取り、入力変数として車両の進行方向の角度ｕを取る。そして、式（２０）で示される
【数２０】
2007141121.tif 000022
のモデルを、ｄｘ１／ｄｔ＝ｆ１（ｘ１，ｘ２）＋ｇ１（ｘ１，ｘ２）ｕ及びｄｘ２／ｄｔ＝ｆ２（ｘ１，ｘ２）＋ｇ２（ｘ１，ｘ２）ｕとする。したがって、ｆ１（ｘ１，ｘ２）、ｇ１（ｘ１，ｘ２）、ｆ２（ｘ１，ｘ２）、ｇ２（ｘ１，ｘ２）の各関数形を用いて車両のモデルを生成し、データセット収集工程で収集されたデータセットの範囲に基いて規定されるｘ１、ｘ２、ｕの制御適用範囲を表す集合であるＰｘ１、Ｐｘ２、Ｐｕを用いて車両の進行方向姿勢を自動的に制御することができる。
【００２５】
また、目標制御状態を、駐車操作における目標縦方向速度ｘ１ｄと目標横方向速度ｘ２ｄとするので、車両の駐車、並列駐車における車両の進行方向姿勢制御を自動的に行うことができる。
【００２６】
また、モデル生成にあたっては、車両のモデルを、連続的値を有する状態変数及び入力変数と、離散的値を有する状態変数及び遷移変数とで表現されるハイブリッドモデルとして生成する。例えば駐車制御における車両については、前進、停止、後退の３つの状態と、状態間の遷移とでモデル化ができるが、３つの状態の区別、各状態間の遷移は、連続的変数よりも離散的変数で表現する方が現実に近い。したがって、連続的値を有する変数と、離散的値を有する変数との表現を用いることで、現実の状態に適したシステムのモデル化を行うことができる。
【００２７】
上記のように、本発明に係る非線形ダイナミクスを有する制御対象のモデル化装置及びその上で実行されるプログラムによれば、非線形ダイナミクスを有する制御対象についてモデル化が容易となる。また、本発明に係る非線形ダイナミクスを有する制御対象のための制御装置、非線形ダイナミクスを有する制御対象のモデル化プログラムによれば、求められたモデルを用いて非線形ダイナミクスを有する制御対象の制御が可能となる。また、本発明に係る非線形ダイナミクスを有する車両のための制御装置によれば、求められたモデルを用いて非線形ダイナミクスを有する車両の制御が可能となる。
【発明を実施するための最良の形態】
【００２８】
以下に図面を用いて本発明に係る実施の形態につき詳細に説明する。以下において、制御対象は車両とし、より具体的には、並列駐車において前進、後退、停止を繰り返して車両の進行方向の姿勢を自動的に制御することを説明するが、駐車制御以外の車両制御であってもよく、また、車両以外の制御対象であってもよい。例えば、宇宙ロボットや医療用ロボットの高性能制御に用いられる制御装置及び制御プログラムのように、高精度の制御系の設計が必要なものについて適用されてもよい。
【００２９】
図１は、車両１０の進行方向における姿勢を制御する車両制御システム３０の全体構成を示す図である。具体的には、車両１０を駐車場において他の車両駐車スペースに対し並列駐車させるために、所定の低速度で、前進、停止、後退を繰り返させる無人操作自動制御を行う機能を車両制御システム３０は有する。ここで無人操作自動制御とは、運転者のステアリング操作なしで、車両１０を姿勢制御することを意味している。そして、低速下における駐車操作状態は、非線形のダイナミクスであるので、車両制御システム３０は、非線形ダイナミクスを制御対象とする制御システムである。
【００３０】
なお図１では、車両制御システム３０自体の構成要素ではないが、制御対象である車両１０と、車両制御システム３０を構築するために予め車両１０の制御に関するパラメータを収集するための模擬車両２０が図示されている。
【００３１】
制御対象である車両１０は、前輪１２及び後輪１４を備え、車両１０の進行方向の姿勢の変更を行うためのステアリング機構１６と、車両１０の長手方向に沿った速度成分である縦方向速度ｘ１と、車両１０の幅方向に沿った速度成分である横方向速度ｘ２とを検出する状態検出センサ１８とを有している。状態検出センサ１８はその他に車両１０の前進、停止、後退等を区別して検出する機能を有するものとすることができる。ステアリング機構１６は、ステアリング指令を受けて、車両１０の進行方向の姿勢を変更するものであるが、具体的には、車両の長手方向を基準として前輪１２の角度ｕを変更する機能を有する。
【００３２】
模擬車両２０は、もちろん制御対象である車両１０そのものを用いることができる。また、すでにシミュレーション実績が確立されている車両シミュレータを用いることができる。そのような車両シミュレータとして、米国ミシガン交通研究院大学（Ｕｎｉｖｅｒｓｉｔｙ  ｏｆ  Ｍｉｃｈｉｇａｎ  Ｔｒａｎｓｐｏｒｔａｔｉｏｎ  Ｒｅｓｅａｒｃｈ  Ｉｎｓｔｉｔｕｔｅ）によって開発された車両運転シミュレータであるＣａｒＳｉｍを用いることができる。
【００３３】
車両制御システム３０は、大別して２つの機能及び構成に分けて考えることができる。その１つは、車両１０の制御特性を模擬した模擬車両２０から、制御に関するパラメータを多数収集し、それを用いて車両１０の無人操作並列駐車の姿勢制御に関するモデル化を行う機能及び構成である。このモデル化機能は、モデル化コンピュータ４０によって実現される。その２は、モデル化コンピュータ４０によって得られた車両１０のモデルを用い、実際の車両１０について、無人操作並列駐車の姿勢制御を行う機能及び構成である。この車両制御機能は、車両制御コンピュータ７０によって実現される。以下では、最初にモデル化コンピュータ４０により実行されるモデル化機能について説明し、次に車両制御コンピュータ７０により実行される車両制御機能について説明する。なお、モデル化コンピュータ４０と車両制御コンピュータ７０とは別個独立のコンピュータとして説明するが、もちろん、モデル化機能と車両制御機能とを併せもつ１つのコンピュータで構成してもよい。
【００３４】
１．モデル化機能
図１にはモデル化機能を実行するためのモデル化コンピュータ４０の構成が示されている。モデル化コンピュータ４０は、ＣＰＵ４２と、模擬車両２０から収集されるデータセットを取得する等の機能を有する入力部４４と、モデル化の結果を車両制御コンピュータ７０に出力する等の機能を有する出力部４６と、プログラムやデータ等を記憶する記憶部４８を含んで構成され、これらの要素は内部バスによって相互に接続される。なお、取得されたデータセットはファイル５０として記憶部４８に記憶される。また、入力部４４は、他の入力データを入力するためのキーボード等を含むことができ、出力部４６はディスプレイやプリンタを含むものとすることができる。
【００３５】
ＣＰＵ４２は、模擬車両２０からモデル化に必要なデータセットを収集し取得するデータセット収集モジュール５２と、縦方向速度ｘ１、横方向速度ｘ２、車両姿勢角度ｕについて複数の部分多項式を設定し、収集されたデータセットに基いて各部分多項式についてその多項式係数を算出する多項式係数算出モジュール５４と、算出された多項式係数を有する複数の部分多項式について、算出すべきｄｘ１／ｄｔ、ｄｘ２／ｄｔに対する代表性を評価する部分多項式評価モジュール５６と、評価に基き、式（２１）で示される
【数２１】
2007141121.tif 000023
の形式でモデル化を行うモデル化モジュール５８と、モデル化された内容を出力する出力モジュール６０とを含んで構成される。これらの機能はソフトウェアで実現でき、具体的には、対応するモデル化プログラムを実行することで実現できる。これらの機能の一部をハードウェアで実現することもできる。
【００３６】
かかる構成の作用、特にＣＰＵ４２の各機能について図２のフローチャート及び図３のデータファイルを用いて詳細に説明する。図２は、非線形ダイナミクスを有する車両のモデル化の手順を示すフローチャートで、各手順は、モデル化プログラムの各処理手順に相当する。
【００３７】
モデル化コンピュータ４０においてモデル化プログラムを立ち上げると、最初に、模擬車両２０から、できるだけ多くのデータセットを収集し、これらを取得する（Ｓ１０）。具体的には、ＣＰＵ４２のデータセット収集モジュール５２の機能により、模擬車両２０であるシミュレータＣａｒＳｉｍに指令が与えられ、様々な条件の組み合わせのデータセットの収集を実行させ、その結果をモデル化コンピュータ４０に伝送させて、これを取得する。取得された各データセットは、一旦記憶部４８に記憶される。
【００３８】
シミュレータＣａｒＳｉｍには、車両１０を模擬するための車両データが予め与えられる。一例を上げると、車両１０の長手方向寸法として３８００ｍｍ、車両１０の質量として１５２７ｋｇ、車両１０の重心から前輪１２及び後輪１４までの距離としてそれぞれ１０１４ｍｍ及び１６７６ｍｍ、前輪１２及び後輪１４の半径としてそれぞれ２７７．５６６ｍｍ及び２８６．３８７ｍｍ、車両１０のロール慣性、ピッチ慣性、ヨー慣性としてそれぞれ６０６．１ｋｇ・ｍ２、２７４１．９ｋｇ・ｍ２、２７４１．９ｋｇ・ｍ２等がＣａｒＳｉｍに与えられる。ＣａｒＳｉｍはこれらの車両条件が与えられると、車両１０についてシミュレーションを実行し、車両１０の縦方向速度ｘ１、横方向速度ｘ２、車両１０の進行方向の角度ｕ、ｄｘ１／ｄｔ、ｄｘ２／ｄｔを１組のデータセットとして、上記車両条件を満たす多くのデータセットを算出することができる。
【００３９】
例えば、車両１０を所定の低速度下で走行状態とし、その条件の下で任意のステアリング操作を表す角度ｕをサンプリングサイクルに合わせて変化させ、その応答としてｘ１、ｘ２、ｄｘ１／ｄｔ、ｄｘ２／ｄｔをサンプリングサイクルに合わせて出力させ、同じサンプリングサイクルの角度ｕと対応付けてデータセットとすることができる。角度ｕの変化は、ステアリングの急旋回の限度範囲等を想定して設定することができる。
【００４０】
ここで、図１の車両１０で説明したように、車両１０のモデル化においては、ｘ１、ｘ２が観察可能な状態変数に対応し、ｕが車両１０に入力される入力変数に対応し、ｄｘ１／ｄｔ、ｄｘ２／ｄｔがモデル化において算出すべき状態変数に対応する。
【００４１】
図３のデータファイル１００は、収集されたデータセットを用いて車両１０のモデル化を行うときの各データの様子を模擬的に示す図である。ここでは、収集されたデータセットの用途を示すセット１０２の欄と、収集された各データセットを識別するための番号１０４の欄と、収集されたデータセットの内容が書き込まれるデータ収集１０６の欄等が設けられている。
【００４２】
なお、上記のように、モデル化において算出すべき状態変数は、ｄｘ１／ｄｔ、ｄｘ２／ｄｔの２つがあり、図３に示される例は、このうちのｄｘ１／ｄｔに関するデータファイルである。したがって、モデル化のためのデータファイルとしては、図３と同様の内容のｄｘ１／ｄｔに関するデータファイルがもう１つあることになる。また、車両１０の動作モデルが複数あるときは、各モデルについてそれぞれデータ収集を行うので、その場合には、動作モデルごとに図３のようなデータファイルがあることになる。例えば、無人並列駐車のモデルを、前進モデル、停止モデル、後退モデルのように分けて構築する場合には、それぞれのモデルについてデータセットの収集を行い、したがってそれぞれに図３のようなデータファイルがあることになる。
【００４３】
図３の例では、データ収集１０６の欄に、収集されたデータセットとして合計５０，０００セットが書き込まれている。収集されるデータセットの数は、できるだけ多い方が後の車両制御のときの適用範囲が広くなる。しかしデータ収集の時間がその分長くなるので、適当な数を予め定めておくことができる。所定の数のデータセットの収集がなされると、これをトレーニングセットとチェッキングセットに２分する（Ｓ１２）。トレーニングセットは、部分多項式の多項式係数を算出するために用いるものであり、チェッキングセットは算出された多項式係数を有する部分多項式の評価を行うものであるので、後者に比して前者のデータセットの数が多い方が好ましい。例えば、トレーニングセット：チェッキングセット＝７：３程度に二分することができる。図３の例では、トレーニングセットを３５，０００セット、チェッキングセットを１５，０００セットとして二分し、それぞれの区分がセット欄１０２に示されている。
【００４４】
次に、状態変数ｘ１、ｘ２、入力変数ｕ について部分多項式ｚを設定する（Ｓ１４）。この部分多項式ｚは、算出すべき状態変数ｄｘ／ｄｔを、観察可能な状態変数ｘと、入力変数ｕとの多項式ｚで表現しようとするものである。なお、一般的には、ｘと表記すればこれは位置であり、ｄｘ／ｄｔと表記すればこれは速度を示すことが多いが、ここでは、状態変数をｘとし、その時間微分をｄｘ／ｄｔとしているので、状態変数ｘを速度にとれば、ｄｘ／ｄｔは、加速度を示していることになるので、注意が必要である。
【００４５】
ここで、制御対象は、観察可能な状態変数は非線形であり、入力変数は線形である非線形ダイナミクスを有するものとする。この場合、多項式展開は、非線形状態変数については非線形展開となり、線形入力変数については線形展開にすることがよい。具体的には、部分多項式ｚは、２つの変数の多項式展開の形式とし、各変数が非線形の属性を有するか、線形の属性を有するかで、次のような４通りの形式の部分多項式を用いる。
【００４６】
以下の式では、ζｐ、ζｑを変数とし、ζＮとξＮを非線形の属性を有する変数により構成される集合、ζＬとξＬを線形の属性を有する変数により構成される集合とし、ａｋｉを多項式係数とする。ｋは、状態変数を区別する添字符号で、上記の例では、状態変数ｘ１に対応するものがｋ＝１、状態変数ｘ２に対応するものがｋ＝２と考えてよい。ｉは、複数ある多項式係数をそれぞれ区別する添字符号である。
【数２２】
2007141121.tif 000024
【数２３】
2007141121.tif 000025
【数２４】
2007141121.tif 000026
【数２５】
2007141121.tif 000027
【００４７】
式（２２）は、２つの変数とも非線形の属性を有する場合の部分多項式の形式を示す。式（２２）の次の記述は、この場合の部分多項式の結果物であるｚｋが、非線形の属性を有する集合に加えられることを示す。つまり、２つの非線形変数の部分多項式の結果物は、非線形である。
【００４８】
式（２３）、式（２４）は、一方の変数が非線形の属性を有し、他方の変数が線形の属性を有する場合の部分多項式の形式を示す。この場合に、線形の属性を有する変数は、たかだか１次形式で展開されるのに対し、非線形の属性を有する変数は２次形式まで展開される。このときの部分多項式の結果物は、線形の属性を有する集合に加えられる。
【００４９】
式（２５）は、２つの変数とも線形の属性を有する場合の部分多項式の形式を示す。この場合も、線形の属性を有する変数は、たかだか１次形式で展開される。このときの部分多項式の結果物は、線形の属性を有する集合に加えられる。
【００５０】
再び図３に戻って説明すると、第１世代部分多項式１０８の欄に、式（２２）から式（２５）の結果が書き込まれる。ここで、第１世代部分多項式では、変数として、データ収集１０６の欄の３つの変数ｘ１、ｘ２、ｕを用い、その３つの変数のうち２つを用いて部分多項式を構成する。図３において、線形・非線形１２０の欄は、部分多項式ｚの属性を示し、部分多項式１２２の欄は、複数の部分多項式を符号で示し、構成要素１２４の欄は、部分多項式の展開要素である２つの変数を示す。
【００５１】
例えば、部分多項式ｚ１は、変数ｘ１、ｘ２で展開され、ここで変数ｘ１、ｘ２はいずれも非線形の属性を有するので、式（１）の展開形式となり、その結果物の属性は非線形である。線形・非線形１２０の欄でＮとは非線形の属性を有することを示す。また、部分多項式ｚ２は、変数ｘ１、ｕで展開され、ここで変数ｘ１は非線形、ｕは線形の属性を有するので、式（２３）又は式（２４）の展開形式となり、その結果物の属性は線形である。線形・非線形１２０の欄でＬとは線形の属性を有することを示す。同様に、部分多項式ｚ３は、変数ｘ２、ｕで展開され、ここで変数ｘ２は非線形、ｕは線形の属性を有するので、式（２３）又は式（２４）の展開形式となり、その結果物の属性は線形である。
【００５２】
ここで、後述する第２世代、第３世代等の部分多項式について先に説明する。第２世代以降の部分多項式は、チェッキングセットでの評価において、第１世代部分多項式の形式では算出すべき状態変数ｄｘ／ｄｔを代表するのにまだ十分でないとされるときに用いられる、さらに高次の部分多項式である。例えば第２世代部分多項式は、データ収集された変数と、第１世代部分多項式の結果物との中の２つの変数を展開要素として構成される。すなわち、図３の例では、ｘ１、ｘ２、ｕに加え、ｚ１、ｚ２、ｚ３も部分多項式の展開要素となる。したがって、図３の第２世代部分多項式１１０の欄に示されるように、第２世代部分多項式はｚ４からｚ９までの６種類あることになる。
【００５３】
例えば、部分多項式ｚ４は、変数ｘ１、ｚ１で展開され、ここで変数ｘ１、ｚ１はいずれも非線形の属性を有するので、式（２２）の展開形式となり、その結果物の属性は非線形である。また、部分多項式ｚ５は、変数ｘ１、ｚ２で展開され、ここで変数ｘ１は非線形、ｚ２は上記のように線形の属性を有するので、式（２３）又は式（２４）の展開形式となり、その結果物の属性は線形である。この場合でも、ｚ５についての式（２３）又は式（２４）を注意深く観察してわかるように、線形の属性を有する変数ｕは、たかだか１次形式でしか展開されていない。他の部分多項式ｚ６以下の内容も同様に説明できる。図３の線形・非線形１２０の欄、部分多項式１２２の欄、構成要素１２４の欄には、これらの内容が書き込まれている。
【００５４】
第３世代多項式は、データ収集されたｘ１、ｘ２、ｕ、第１世代部分多項式の結果物ｚ１、ｚ２、ｚ３に加え、第２部分多項式の結果物ｚ４、ｚ５、ｚ６、ｚ７、ｚ８、ｚ９も部分多項式の展開要素となる。なお、ｚ１１のように、部分多項式の結果物が線形の属性を有するものは、その部分多項式に対応する式（２３）又は式（２４）、あるいは式（２５）を注意深く観察してわかるように、線形の属性を有する変数ｕは、たかだか１次形式でしか展開されていない。
【００５５】
再び図２に戻り、部分多項式の設定が終わると、トレーニングセットのデータを用いて、各部分多項式の多項式係数ａｋｉの算出が行われる（Ｓ１６）。具体的には、ＣＰＵ４２の多項式係数算出モジュール５４の機能により、各部分多項式の基本展開要素であるｘ１、ｘ２、ｕにトレーニングセットのデータセットの値を代入する。図３の例では３５，０００セットのデータセットの値を１つの部分多項式に対し適用し、未知の多項式係数を統計的に算出する。統計的算出には、算出すべき状態変数と、部分多項式の結果物との差の二乗の累積が最小になるような最小二乗法を用いることができる。すなわち、式（１）から式（４）と同様の表記を用いて、算出すべき状態変数をｄζ／ｄｔとして、トレーニングセットの全データセットについて式（２６）を最小にするような多項式係数を算出する。
【数２６】
2007141121.tif 000028
式（２６）において累積範囲の１からｊは、トレーニングセットの複数のデータセットの全部を示すものであり、図３の例ではｊ＝３５，０００である。
【００５６】
図３において、処理１３０はそのことを示している。このようにして、トレーニングセットのデータセットに基き、複数の部分多項式のすべてについて、その多項式係数を算出することができ、ここにおいて、すべての部分多項式の内容を統計的に確定できる。
【００５７】
次に、チェッキングセットによる各部分多項式の評価が行われる（Ｓ１８）。評価には、次の式（２７）で算出される評価指標ｒが用いられる。
【数２７】
2007141121.tif 000029
式（２７）において累積範囲のｊ＋１からｎは、チェッキングセットの複数のデータセットの全部を示すものであり、図３の例ではｊ＋１＝３５，００１であり、ｎ＝５０，０００である。式（２７）の評価指標ｒは、標準偏差に類似の無次元数であり、値が小さいほど算出すべき状態変数ｄζ／ｄｔからの誤差が少ない。
【００５８】
図３において処理１４０はそのことを示している。そして、評価１４２の欄には、算出されたｒを用いて、算出すべき状態変数ｄζ／ｄｔを最も代表する部分多項式を絞る手順が示されている。すなわち、各部分多項式ｚ１、ｚ２、ｚ３、ｚ４等について算出された評価指数ｒが欄１４４に書き込まれ、つぎに、適当な閾値Ｒを用いて、異常な値のｒを除外した結果が欄１４６に示される。欄１４６において×マークは、閾値Ｒ以上のｒであり、これらは以後の評価において考慮されないことを示す。そして、第１世代多項式の中で最も評価指標ｒの値がよいもの、第２世代多項式の中で最も評価指標ｒの値がよいものを選び、これを欄１４８に書き込む。図３の例では、第１世代多項式の中で最も評価指標ｒの値がよいものがｒ２で、第２世代多項式の中で、最も評価指標ｒの値がよいものがｒ６であることが示されている。
【００５９】
そして、評価が十分か否かを判断する（Ｓ２０）。具体的には、第１世代多項式の中で最も評価がよい部分多項式の評価指標の値と、第２世代多項式の中で最も評価がよい部分多項式の評価指標の値とを比較する。図３の例ではｒ２の値とｒ６の値を比較する。ここでｒ２＞ｒ６であると、評価が十分であるとは判断されず、さらに高次の部分多項式を設定する（Ｓ２２）。上記の例では、さらに第３世代部分多項式について評価指標ｒを求める。すなわち、その世代多項式の評価より次世代多項式の評価の方が良い場合は、さらに高次の世代の部分多項式を評価する。これを繰り返し、次世代多項式の評価が改善されなくなるまで上記の処理を行う。最も、適当な繰り返し数で評価を打ち切ることもできる。このようにして、最も評価指標が良い部分多項式を、算出すべき状態変数の代表性の高いものとして決定することができる。
【００６０】
評価が十分であるとして最も代表性の高い部分多項式が決定されると、これを整理して、式（２８）の形式に整え、制御対象を式（２８）でモデル化する（Ｓ２４）。具体的には、ＣＰＵ４２のモデル化モジュール５８の機能により、式（２８）の形式に整理されてモデル化が行われる。ここで、観察可能な状態変数はｘ１、ｘ２の２つであるが、式（２８）では、代表的に状態変数を単にｘで示してある。
【数２８】
2007141121.tif 000030
【００６１】
上記のように、各部分多項式の展開において、線形の属性を有する入力変数ｕはたかだか１次形式でしか展開されないので、どの部分多項式の形式が最も代表性の高い部分多項式として決定されても、式（２８）の形式に整理することができる。例えば図３の例で、評価指標ｒ６が最も良いとすると、最も代表性の高い部分多項式はｚ６であり、その部分多項式の形式は式（２９）のようになる。
【数２９】
2007141121.tif 000031
【００６２】
このようにしてモデル化が行われると、その成果として、式（７）における関数形ｆ（ｘ）、ｇ（ｘ）と、データ収集されたデータセットの集合を示すＰｘ、Ｐｕとが出力される（Ｓ２６）。具体的には、ＣＰＵ４２の出力モジュール６０の機能により、ｆ（ｘ）、ｇ（ｘ）と、Ｐｘ、Ｐｕとが、出力部４６から出力され、適当なデータ伝送媒体を介して車両制御コンピュータ７０に入力される。車両制御コンピュータ７０は、入力されたｆ（ｘ）、ｇ（ｘ）と、Ｐｘ、Ｐｕとを取得し、その記憶部７８に記憶する。
【００６３】
ここでＰｘ、Ｐｕを説明する。車両１０のモデル化として必要なものは、関数形ｆ（ｘ）、ｇ（ｘ）であるが、この関数形ｆ（ｘ）、ｇ（ｘ）は、収集された複数のデータセットに基いて統計的処理によって得られたものである。したがって、収集されたデータセットがカバーする範囲を超える状態変数、入力変数の場合には、この関数形ｆ（ｘ）、ｇ（ｘ）を用いたモデル化が十分でないことが生じえる。特に非線形ダイナミクスを取り扱う場合には、そのような注意が必要である。そこで、関数形ｆ（ｘ）、ｇ（ｘ）を用いたモデル化は、収集されたデータセットの範囲の状態変数ｘ、入力変数ｕの中で十分有効であることを示すため、状態変数ｘ、入力変数ｕの適用可能範囲として、Ｐｘ、Ｐｕを示すことが望ましい。
【００６４】
したがって、Ｐｘ、Ｐｕとは、制御モデルに適用可能な、いわゆるパーミッシブルな変数の集合を示すものである。具体的には、収集されたデータセットの中の状態変数ｘの最大値と最小値等でＰｘを表すことができ、収集されたデータセットの中の入力変数ｕの最大値と最小値等でＰｕを表すことができる。
【００６５】
さらに数式化してＰｘを示すには式（３０）を用いることができる。
【数３０】
2007141121.tif 000032
ここで、ｉは、上記の例でｘ１、ｘ２を区別する符号であり、ｉ＝１，２である。ｋはサンプリングの順番を示す符号であり、現在がｋ、次のサンプリングサイクルのときがｋ＋１、現在の１つ過去のサンプリングサイクルのときがｋ−１である。Ｔｓはサンプリング間隔である。添字のＭとｍとは、それぞれ、収集されたデータセットの範囲における最大値と最小値とを意味している。
【００６６】
同様に、数式化してＰｕを示すには式（３１）を用いることができる。
【数３１】
2007141121.tif 000033
ここで、ｋ−２は、ｋ−１よりも１つ過去のサンプリングタイムのときである。
【００６７】
以上がモデル化機能についての説明であり、次に、モデル化機能により得られた関数形ｆ（ｘ）、ｇ（ｘ）と、パーミッシブルなの集合を示すＰｘ、Ｐｕを用いた車両１０の制御について説明する。
【００６８】
２．車両制御機能
再び図１に戻ると、図１には、モデル化機能により得られたｆ（ｘ）、ｇ（ｘ）と、Ｐｘ、Ｐｕとを用いて車両１０の姿勢制御を実行するための車両制御コンピュータ７０の構成が示されている。車両制御コンピュータ７０は、ＣＰＵ７２と、モデル化コンピュータ４０から出力されるｆ（ｘ）、ｇ（ｘ）と、Ｐｘ、Ｐｕとを取得する等の機能を有する入力部７４と、ディスプレイやプリンタ等の出力部７６と、プログラムやデータ等を記憶する記憶部７８と、車両１０の状態検出センサ１８との間のインタフェース回路であるセンサＩ／Ｆ８２と、ステアリング機構１６との間のインタフェース回路であるステアリングＩ／Ｆ８４とを含んで構成され、これらの要素は内部バスによって相互に接続される。なお、取得されたｆ（ｘ）、ｇ（ｘ）、Ｐｘ、Ｐｕは、ファイル８０として、記憶部７８に記憶される。また、入力部４４は、他の入力データを入力するためのキーボード等を含むことができる。
【００６９】
ＣＰＵ７２は、記憶部７８からｆ（ｘ）、ｇ（ｘ）、Ｐｘ、Ｐｕを読み出す読出モジュール８６と、ｆ（ｘ）、ｇ（ｘ）に基いて車両１０のモデルを生成するモデル生成モジュール８８と、車両１０から状態変数ｘ（ｋ）を取得する取得モジュール９０と、生成された車両１０のモデルに従ってｘ（ｋ＋１）を算出するｘ（ｋ＋１）算出モジュール９２と、算出されたｘ（ｋ＋１）を目標制御状態に対し評価する評価モジュール９４と、評価に基き入力変数ｕ（ｋ）を車両１０に出力する入力変数ｕ（ｋ）出力モジュール９６とを含んで構成される。これらの機能はソフトウェアで実現でき、具体的には、対応する車両制御プログラムを実行することで実現できる。これらの機能の一部をハードウェアで実現することもできる。
【００７０】
かかる構成の作用、特にＣＰＵ７２の各機能について図４のフローチャートを用いて詳細に説明する。図４は、非線形ダイナミクスを有する車両の制御の手順を示すフローチャートで、各手順は、車両制御プログラムの各処理手順に相当する。
【００７１】
車両制御コンピュータ７０の記憶部７８には、予めモデル化コンピュータ４０から入力されたｆ（ｘ）、ｇ（ｘ）、Ｐｘ、Ｐｕがファイル８０に記憶されている（Ｓ３０）。Ｐｘ、Ｐｕは、上記のように、モデル化コンピュータ４０によって収集された多数のデータセットの範囲を表すものであるが、記憶部７８は、この収集された多数のデータセットの内容全部を記憶するのではなく、式（９）、式（１０）で表される内容のみを記憶する。すなわち、上記の例で言えば、５０，０００セットのデータセットの内容の記憶の代わりに、状態変数ｘ、入力変数ｕの最大値、最小値を表すデータのみを記憶する。このことにより、記憶部７８の記憶負荷を軽くすることができる。
【００７２】
そして、車両制御コンピュータ７０において車両制御プログラムが立ち上げられると、このｆ（ｘ）、ｇ（ｘ）、Ｐｘ、Ｐｕが読み出される（Ｓ３２）。この機能はＣＰＵ７２の読出モジュール８６によって実行される。ここで、車両１０の無人操作並列駐車の姿勢制御を行う場合には、並列駐車の制御モデルに対応するｆ（ｘ）、ｇ（ｘ）、Ｐｘ、Ｐｕが読み出される。すなわち、駐車の制御は、車両１０を前進、停止、後退させることを含むので、車両１０の前進の場合に用いられるｆ（ｘ）、ｇ（ｘ）、Ｐｘ、Ｐｕ、停止の場合に用いられるｆ（ｘ）、ｇ（ｘ）、Ｐｘ、Ｐｕ、後退の場合に用いられるｆ（ｘ）、ｇ（ｘ）、Ｐｘ、Ｐｕがそれぞれ読み出される。換言すれば、モデル化コンピュータ４０におけるモデル化は、車両１０の駐車の際の前進、停止、後退についてデータを収集し、モデル化を行っている。なお、状態変数はｘ１、ｘ２の２つがあるので、ｆ（ｘ）、ｇ（ｘ）、Ｐｘもそれぞれ２つあることになるが、ここでは説明の便宜上、ｘで代表させている。
【００７３】
次に、読み出されたｆ（ｘ）、ｇ（ｘ）を用いて、車両１０のモデルが生成され、いわゆるシステムのモデル化が行われる（Ｓ３４）。この機能は、ＣＰＵ７２のモデル生成モジュール８８によって実行される。並列駐車の制御モデルは、車両１０の前進、停止、後退の３つの状態と、その間の遷移とを含むように生成される。
【００７４】
図５は、並列駐車の制御のための車両１０のモデルを示す図である。ここで、前進、停止、後退の各状態における制御モデルには、式（２８）に従って、これをサンプリングサイクルで表現した式（３２）の形式を用いている。
【数３２】
2007141121.tif 000034
そして、前進、停止、後退の各状態を示すものとして、２値の状態変数ｘｄを用い、各状態の間の遷移を示すものとして、２値の遷移変数δを用いている。つまり、式（７）で用いられる連続的値を有する状態変数ｘ及び入力変数ｕとともに、離散的な値を有する状態変数ｘｄと遷移変数δとを合わせて用いて、車両１０のモデルが生成されている。連続的値の変数と離散的値の変数とを合わせ用いていることから、このモデルをハイブリッドモデルと呼ぶことができる。
【００７５】
２値の状態変数ｘｄは、目標制御状態を示すもので、ｘｄ＝［０  １］Ｔが前進を、ｘｄ＝［１  １］Ｔが停止を、ｘｄ＝［１  ０］Ｔが後退をそれぞれ表している。これを数式化して示すには、式（３３）を用いることができる。ここで、ｆ（ｘ（ｋ））、ｇ（ｘ（ｋ））に付けられた添字は、ｆが前進、ｓが停止、ｂが後退を示している。すなわち、ｆ（ｘ（ｋ））、ｇ（ｘ（ｋ））は、前進、停止、後退においてそれぞれ異なるものが用いられている。
【数３３】
2007141121.tif 000035
【００７６】
遷移変数δは、前進から停止へ遷移するとき、停止から後退へ遷移するときには、δ＝［０  １］Ｔを用い、逆に後退から停止へ遷移するとき、停止から前進へ遷移するときには、δ＝［１  ０］Ｔ  を用い、各状態で留まっているときには、δ＝［０  ０］Ｔを用いて、状態間の遷移を表している。これを数式化して示すには、式（３４）を用いることができる。ここで、ｓ（ｋ）は、車両１０の状態検出センサ１８において検出される方向ベクトルであり、ＹＳＦ、ＹＦＳ、ＹＳＢ、ＹＢＳ、ＹＳＳ、ＹＦＦ、ＹＢＢは、それぞれ車両１０が停止から前進への切り替え、前進から停止への切り替え、停止から後退への切り替え、後退から停止への切り替え、停止、前進、後退を続けるための必要な状態検出センサ１８の状態を示す。
【数３４】
2007141121.tif 000036
【００７７】
また、各状態の遷移を数式で示すには式（３５）を用いることができる。
【数３５】
2007141121.tif 000037
【００７８】
再び図４に戻り、車両１０のシステムモデル化がなされると、次に車両１０から状態変数ｘ（ｋ）が取得される（Ｓ３６）。具体的には、ＣＰＵ７２の状態変数ｘ（ｋ）取得モジュール９０の機能により、センサＩ／Ｆ８２を介して車両１０の状態検出センサ１８に指令が出され、状態変数ｘ（ｋ）として、車両１０の縦方向速度ｘ１（ｋ）、横方向速度ｘ２（ｋ）が各サンプルサイクルごとに検出され、再びセンサＩ／Ｆ８２を介し伝送されて、ＣＰＵ７２がこれを取得する。
【００７９】
次に、取得された状態変数ｘ（ｋ）のサンプルサイクルｋの次のサンプルサイクル（ｋ＋１）における状態変数ｘ（ｋ＋１）を算出する（Ｓ３８）。具体的には、ＣＰＵ７２のｘ（ｋ＋１）算出モジュール９２の機能により次のようにして実行される。すなわち、車両１０のモデルを示す式（７）をサンプルサイクルにおける形式に書き換え、式（３６）で示される形式を用い、Ｐｕの要素である入力変数ｕについてｘ（ｋ＋１）を算出する。
【数３６】
2007141121.tif 000038
そして、算出されたｘ（ｋ＋１）が、Ｐｘの要素であるか否か判断される（Ｓ４０）。ｘ（ｋ＋１）がＰｘの要素であると判断されると、次の評価工程（Ｓ４２）に進み、そうでなければ、入力変数ｕを変更してＳ３８、Ｓ４０の工程が繰り返される。このときも、入力変数ｕはＰｕの要素でなければならないことはもちろんである。
【００８０】
このようにして、入力変数ｕがＰｕの要素であり、状態変数ｘ（ｋ＋１）がＰｘの要素である条件を満たすｘ（ｋ＋１）が得られると、目標制御状態に対するｘ（ｋ＋１）の評価が行われる（Ｓ４２）。評価は、ＣＰＵ７２の評価モジュール９４の機能によって実行される。具体的には、ｘ１（ｋ＋１）、ｘ２（ｋ＋１）の総合評価が二次計画法基準によって行われる。すなわち、ｘ１の目標制御状態をｘ１ｄ、ｘ２の目標制御状態をｘ２ｄとして、式（３７）で示される二次計画法評価基準が計算され、その値を最小とするｘ１（ｋ＋１）、ｘ２（ｋ＋１）に対応する入力変数ｕ（ｋ）を求めることが行われる。
【数３７】
2007141121.tif 000039
【００８１】
ここで目標制御状態ｘ１ｄ、ｘ２ｄは、それぞれ縦方向速度及び横方向速度の目標値である。すなわち、並列駐車の場合には、ごく低速で制御が行われるので、その速度を目標制御状態とすることが好ましい。目標制御状態としてはこのほかに、駐車目標位置等を設定することができる。その場合には、駐車目標位置を速度等に変換して二次計画法基準を計算することができる。
【００８２】
求められたｕ（ｋ）は、次のサンプルサイクルにおける車両１０に対する入力変数として、ステアリングＩ／Ｆ８４を介し、車両１０のステアリング機構１６に出力される（Ｓ４４）。この機能は、ＣＰＵ７２の入力変数ｕ（ｋ）出力モジュール９６によって実行される。
【００８３】
このようにして、車両制御コンピュータ７０は、車両１０からの状態変数ｘ（ｋ）に基いて、次々と入力変数ｕを車両１０に出力する。これによって、車両１０の並列駐車の無人操作が行われる。
【実施例１】
【００８４】
車両として、その長手方向寸法が３８００ｍｍ、質量が１５２７ｋｇ、重心から前輪及び後輪までの距離がそれぞれ１０１４ｍｍ及び１６７６ｍｍ、前輪及び後輪の半径がそれぞれ２７７．５６６ｍｍ及び２８６．３８７ｍｍ、ロール慣性、ピッチ慣性、ヨー慣性がそれぞれ６０６．１ｋｇ・ｍ２、２７４１．９ｋｇ・ｍ２、２７４１．９ｋｇ・ｍ２を有するものについて、式（３８）に示される形式となるモデル化を行った。
【数３８】
2007141121.tif 000040
この諸元は、模擬車両としてのＣａｒＳｉｍに与えられた。データ収集に際し用いた入力変数ｕの様子を図６に示す。図６は、縦軸が車両の進行方向に沿った車両の角度ｕで、横軸は時間である。ここでは車両の並列駐車のときのステアリング機構の操作のように、時間と共に、右旋回、左旋回を交互に繰り返し、その角度ｕの大きさも時間と共に次第に大きくしその後次第に小さくして変化させている。このような入力変数ｕの下で、サンプリングサイクルを細かくして、５０，０００セットのｘ１、ｘ２、ｄｘ１／ｄｔ、ｄｘ２／ｄｔのデータセットを収集し、これに基いて、モデル化コンピュータ４０によってｆ（ｘ）、ｇ（ｘ）、Ｐｘ、Ｐｕを求めた。求められたｆ（ｘ）、ｇ（ｘ）、Ｐｘ、Ｐｕは車両制御コンピュータ７０に与えられた。具体的なｆ１（ｘ１、ｘ２）、ｇ１（ｘ１、ｘ２）、ｆ２（ｘ１、ｘ２）、ｇ２（ｘ１、ｘ２）は、式（３９）から（４２）に示す通りである。
【数３９】
2007141121.tif 000041
【数４０】
2007141121.tif 000042
【数４１】
2007141121.tif 000043
【数４２】
2007141121.tif 000044
【００８５】
車両制御コンピュータ７０では、求められたｆ（ｘ）、ｇ（ｘ）に基いて車両のモデルが生成され、そして、Ｐｘ、Ｐｕの下で、二次計画法評価基準が最小となる入力変数ｕ（ｋ）が求められて、各サンプルサイクルで車両に入力された。車両の状態検出センサによって、車両の実際の状態が実測され、その値と、式（４３）に示されるモデルで計算される値とが比較された。
【数４３】
2007141121.tif 000045
【００８６】
図７、図８は、それぞれ車両の縦方向加速度、横方向加速度の実測と計算との比較を示す図である。いずれの図も横軸は時間で、縦軸は図７が縦方向加速度、図８が横方向加速度である。実際には、実測と計算とはきわめて近く、その相関係数は図７において０．９９４６８、図８において０．９９９５６８であった。したがって、図７、図８上で実測と計算との差異を見出すことは実際上困難であった。
【図面の簡単な説明】
【００８７】
【図１】本発明に係る実施の形態において、車両の進行方向における姿勢を制御する車両制御システムの全体構成を示す図である。
【図２】本発明に係る実施の形態において、非線形ダイナミクスを有する車両のモデル化の手順を示すフローチャートである。
【図３】本発明に係る実施の形態において、収集されたデータセットを用いて車両のモデル化を行うときの各データの様子を示す図である。
【図４】本発明に係る実施の形態において、非線形ダイナミクスを有する車両の制御の手順を示すフローチャートである。
【図５】本発明に係る実施の形態において、並列駐車の制御のための車両のモデルを示す図である。
【図６】本発明に係る実施の形態において、データ収集に際し用いる入力変数の様子を示す図である。
【図７】本発明に係る実施の形態において、車両の縦方向加速度の実測と計算との比較を示す図である。
【図８】本発明に係る実施の形態において、車両の横方向加速度の実測と計算との比較を示す図である。
【符号の説明】
【００８８】
１０  車両、１２  前輪、１４  後輪、１６  ステアリング機構、１８  状態検出センサ、２０  模擬車両、３０  車両制御システム、４０  モデル化コンピュータ、４２，７２  ＣＰＵ、４４，７４  入力部、４６，７６  出力部、４８，７８  記憶部、５０，８０  ファイル、７０  車両制御コンピュータ、１００  データファイル、８２  センサＩ／Ｆ、８４  ステアリングＩ／Ｆ。
【図１】
2007141121.tif 000060
【図２】
2007141121.tif 000061
【図３】
2007141121.tif 000062
【図４】
2007141121.tif 000063
【図５】
2007141121.tif 000064
【図６】
2007141121.tif 000065
【図７】
2007141121.tif 000066
【図８】
2007141121.tif 000067