判断・操作統合型運転行動モデルに基づく革新的運転行動支援の創出 タカタ財団助成研究論文集 ISSN 2185

判断・操作統合型運転行動モデルに基づく

革新的運転行動支援の創出

― 平成 22 年度（中間報告） タカタ財団助成研究論文 ―

研究代表者

鈴木 達也

ISSN 2185-8950

研究代表者

名古屋大学大学院工学研究科教授

鈴木

達也

第1章 序論 1 1.1 研究背景 . . . 1 1.1.1 運転支援システムの現状 . . . 1 1.1.2 運転行動の数理モデル化 . . . 3 1.2 研究目的と概要. . . 4 1.3 本論文の構成 . . . 4 第2章 実験装置・条件 6 2.1 実車およびドライビングシミュレータ（DS）の特長 . . . 6 2.2 シミュレータ実験の内容 . . . 7 2.2.1 概要 . . . 7 2.2.2 制御系(PC1) . . . 9 2.2.3 スクリーン表示演算用PC(PC2) . . . 10 2.2.4 最適アシスト量計算用PC(PC3) . . . 10 2.2.5 走行環境(DS) . . . 10 2.2.6 実験手順(DS) . . . 11 第3章 ハイブリッドシステムモデルとその同定手法 14 3.1 ARXモデル . . . 14 3.2 PWARXモデル . . . 15 3.2.1 PWARXモデルの定義 . . . 15 3.2.2 PWARXモデルによるパラメータ同定の流れ . . . 15 3.2.3 K-means法を用いたPWARXモデル同定 . . . 16 3.2.4 サポートベクターマシン(SVM) . . . 20 3.3 PrARXモデル . . . 23

3.3.1 PrARXモデルの定義 . . . 23 3.3.2 モデル例 . . . 23 3.3.3 PrARXモデルのパラメータ同定法 . . . 24 3.3.4 パラメータの逐次更新手法 . . . 25 3.4 ハイブリッドシステムによる前方車追従行動のモデル化 . . . 27 第4章 モデル予測型アシストシステム 31 4.1 モデル予測型運転アシストシステム . . . 31 4.2 アシスト量最適化問題の定式化 . . . 33 4.3 MLDSを用いたドライバモデルの表現 . . . 34 4.3.1 論理式から線形不等式への置き換え . . . 34 4.3.2 MLDSシステムの一般形 . . . 36 4.3.3 運転行動モデルのMLDS表現 . . . 37 4.4 KdBの線形化 . . . 38 4.5 前方車モデルおよび自車モデルの設計 . . . 39 4.6 アシスト量に関する制約 . . . 39 4.6.1 アシスト量の最大値の制約 . . . 39 4.6.2 アシスト変化量の最大値の制約 . . . 40 4.6.3 ペダル操作とアシストの干渉抑制 . . . 40 4.6.4 相対速度によるアシスト量の制約 . . . 40 第5章 実験的検証 41 5.1 PWARXモデルを用いたオフラインシステム同定 . . . 41 5.2 PrARXモデルを用いたオンラインシステム同定 . . . 45 5.3 モデル予測型アシストシステムの評価 . . . 46 5.3.1 実装上の詳細 . . . 46 5.3.2 アシストシステムを実装した実験データ . . . 48 第6章 結論 56 6.1 まとめ . . . 56 謝辞 参考文献

1.1 交通事故発生件数・死者数・負傷者数の推移(昭和25年∼平成20年)(文 献[1]より転載) . . . 2 1.2 ITSの開発分野(文献[2]より転載) . . . 2 2.1 DSの構成図 . . . 8 2.2 DSの運転席 . . . 9 2.3 統合制御プログラムのコントロール画面. . . 10 2.4 走行環境の鳥瞰図 . . . 11 2.5 前方車の走行速度パターン . . . 12 3.1 クラスタリングの流れ . . . 16 3.2 SVMの概略 . . . 20 3.3 ソフトマージン. . . 22 3.4 PrARXモデルの例 . . . 24 3.5 学習用モデルパラメータの時間変化 . . . 26 3.6 逐次更新により推定されたパラメータの時間変化と真値との比較 . . . . 27 3.7 出力信号の比較. . . 28 3.8 KdBのイメージ図 . . . 29

4.1 Model-Predictive Driving Assist System . . . 32

4.2 モデル予測制御の流れ . . . 32

4.3 MLDSの置換例１ . . . 35

4.4 MLDSの置換例２ . . . 35

5.1 クラスタリング結果(E-1：時系列) . . . 43

5.3 クラスタリング結果  (E-1：車間距離−相対速度) . . . 43 5.4 クラスタリング結果  (E-1：相対速度−ペダル操作量) . . . 44 5.5 SVMによるモード分離平面(E-1) . . . 45 5.6 オンライン推定によるモデルパラメータの時間変化. . . 47 5.7 観測データ例 . . . 48 5.8 クラスタリング結果例 . . . 48 5.9 ２つのモードに分離した例 . . . 48 5.10 パラメータ同定用データの分布図(E-1)(車間距離-相対速度) . . . 49 5.11 パラメータ同定用データの分布図(E-1)(相対速度-ペダル操作量) . . . 49 5.12 ω2= 0.05のアシスト車データ分布図(E-1)(車間距離-相対速度) . . . 50 5.13 ω2= 0.05のアシスト車データ分布図(E-1)(相対速度-ペダル操作量) . . 50 5.14 ω2= 0.5のアシスト車データ分布図(E-1)(車間距離-相対速度) . . . 50 5.15 ω2= 0.5のアシスト車データ分布図(E-1)(相対速度-ペダル操作量) . . . 50 5.16 パラメータ同定用データの時系列図(E-1) . . . 51 5.17 ω2= 0.05のアシスト車データの時系列図(E-1) . . . 51 5.18 ω2= 0.5のアシスト車データの時系列図(E-1) . . . 52 5.19 パラメータ同定用データの分布図(E-2)(車間距離-相対速度) . . . 52 5.20 パラメータ同定用データの分布図(E-2)(相対速度-ペダル操作量) . . . 52 5.21 ω2= 0.05のアシスト車データ分布図(E-2)(車間距離-相対速度) . . . 53 5.22 ω2= 0.05のアシスト車データ分布図(E-2)(相対速度-ペダル操作量) . . 53 5.23 ω2= 0.5のアシスト車データ分布図(E-2)(車間距離-相対速度) . . . 53 5.24 ω2= 0.5のアシスト車データ分布図(E-2)(相対速度-ペダル操作量) . . . 53 5.25 パラメータ同定用データの時系列図(E-2) . . . 54 5.26 ω2= 0.05のアシスト車データの時系列図(E-2) . . . 54 5.27 ω2= 0.5のアシスト車データの時系列図(E-2) . . . 55

2.1 DSの構成 . . . 8 3.1 KdB導出に用いる変数一覧 . . . 28 4.1 置換例の変数一覧 . . . 35 5.1 被験者の個人情報(E-1∼E-2) . . . 42 5.2 PWARXモデルパラメータ1 . . . 42 5.3 SVMパラメータ推定結果 . . . 44 5.4 Initial Parametor θ . . . . 46 5.5 Initial Parametor η . . . . 46 5.6 出力誤差の比較. . . 46 5.7 MILPの計算に要する時間 . . . 48 5.8 危険モードのデータ点数(全3042点) . . . 50 5.9 KdB値の二乗平均平方根 . . . 50 5.10 相対速度の二乗平均平方根 . . . 51

序論

1.1

研究背景

1.1.1

運転支援システムの現状

近年，わが国の交通事故死者数は減少傾向にあり，平成21年の死者数は4914人である ことが発表された．これは，第2次ピークである平成4年の半数を下回る数値であり，昭 和27年以来の4千人台となった．しかし，一日平均死者数は13.46人と決して少ない数 ではなく，事故を1件でも減らすことが強く望まれている(図1.1)． このため，世間では交通違反に対する罰則の強化や飲酒運転根絶の呼びかけなどの動き が高まり，それと同時に，事故を未然に防止するための様々な技術の研究もなされてい る．道路や自車両に取り付けられたセンサーやカメラによって周辺を監視し，危険を感知 したらドライバに注意を促すなどの安全技術は，すでに実用化されている．今後，道路交 通の安全化，効率化を目的としたITS（Intelligent Transport Systems: 高度交通システ

ム，図1.2）計画の一環であり，先進技術を活用して安全性を格段に高める先進安全自動

車（ASV: Advanced Safety Vehicle）の開発研究[3]がさらに推し進められ，搭載される

運転支援システムの種類が増えるとともに，さらに普及していくことが予想される． ここで，これまでの運転支援システム開発では，その機能をまず実現することが先決で あり，運転支援システムの機能・性能は平均的な運転者に対応するものであった．車間距 離維持システムや車線維持支援システムなどもその例である． しかし今後，より自動車交通の安全性を向上させると共に自動車交通利用者や社会が受 容できる運転支援システムの開発が重要となっており，個々の運転者に適合したシステム 開発が求められる．この実現のためには，運転者の行動の種々の特性を解析し，その特性 をベースに運転支援システムを開発・設計すること，つまり車を「人間−機械系」として 総合的に捉えることが一つの方策になる[4][5]．

図1.1 交通事故発生件数・死者数・負傷者数の推移(昭和25年∼平成20年)(文献[1]より転載)

1.1.2

運転行動の数理モデル化

個々のドライバに適合したアシストシステムを開発するためには，ドライバの運転行動 を解析し，モデル化することが不可欠である．しかし，現状で実現されている運転支援シ ステムの多くはドライバの個人特性を明示的に考慮しているとは言い難い．この最も大き な要因は，ドライバの個人特性を定量的に反映したモデルが十分に整備されていない点に 起因していると考えられる． 運転行動のモデル化の研究として，従来からファジー理論やニューラルネットワークを

用いたモデルの研究[6][7][8]や，隠れマルコフモデル(Hidden Markov Model: HMM)

を用いた研究[9][10]がある．しかしこれら非線形モデルによる手法は(1) モデル自体が 複雑すぎる，(2) 運転行動を物理的に解釈できない，(3) アシストシステムに応用するに はその有用性に疑問が残る，等の問題がある． そこで，これらの問題を解決するために他の人間行動を取り扱った従来研究に目を移す と，神経科学，システム工学，制御工学，情報処理，認知科学等幅広い分野で様々な視点に 立ったモデルが提案されている．その中でも，神経科学分野におけるE.Bizziら[11][12] による研究や，パターン認識分野におけるC.Breglerら[13]やL.Goncalvesら[14]の研 究の中で実に興味深い報告がなされている．彼らの報告の中で注目すべきは，「人間行動 をプリミティブな動作の合成として捉える」，すなわち「人間行動は単純な動作の組み合 わせで表現できる」という点にある．E.Bizziらによる研究では，プリミティブな動作と いう概念を裏付けるいくつかの動物実験がなされている．またC.Breglerらの研究では行 動のモード分割の考え方に基づき，行動を記録した動画像データに対してパターン認識手 法を適用することで行動の分割を試みている． この考え方を運転行動に応用し，連続モデルと離散モデルの組み合わせによるハイブ リッドシステムとして運転行動を表現する提案がこれまでに筆者らによりなされており [15]，有用性の検証が行われてきた．文献[16]では，人間の前方車追従行動をハイブリッ ドシステムとしてモデル化し，クラスタリング手法を用いてパラメータの同定を行う手法 が考案された．また，運電データをもとに運転者の制御モードの推定を行い，現在の運転 状況が危険なモードにあると判別された場合にブレーキを加えるアシストが提案された． しかしながら，ブレーキを加えるタイミングは運転行動モデルに基づくものの，ブレーキ 量自体の与え方は設計者の経験にもとづき決められている． 運転行動モデルを基にした支援システムにおいては，モデルによる人間行動の再現性が 非常に重要な要素となるが，人間の自動車運転行動における「操作」「判断」の特性は，支 援システムへの慣れ，運転への習熟，疲労などの様々な要因によって常に変化すると考え られる. このため，より積極的な個人適応型運転支援を実現するためには，運転行動モデ ルのパラメータを常に再現性の高い値に更新することが必要となる.

1.2

研究目的と概要

本研究では，個人特性を考慮したアシストシステムの統一的な設計法の確立を目指す． はじめに，人間の運転行動はいくつかの簡単な操作とそれらの切り換えよって実現され ている仮説にもとづき，運転行動をハイブリッド動的システムとして表現することを提案 する．運転行動モデルのパラメータは，実際の運転データをもとにシステム同定を通じて 決定される．これにより，個々の運転者の特性を捉えた運転行動のモデリングおよび解析 が可能となる．具体的には，ダイナミクスの切り換わりを確定的な論理条件で表現する区 分的ARX（PWARX）モデルと，切り換わりを確率的に表現することで判断のあいまい

さが反映される確率重みARX（PrARX）モデルの2種類を考える．まずPWARXモデ

ルに対して統計的手法を用いてパラメータ推定を行い，得られたモデルパラメータを通じ て運転行動を分析する．また，PrARXに対して最急降下法にもとづきパラメータを推定 する方法を提案し，これを用いて時間とともに変化する人間の運転特性を捉えることが可 能なオンラインシステム同定を行う． 次に，同定した運転行動モデルを用いて人間の判断特性を考慮したアシスト系の設計を 行う．その手法として，人間の運転行動モデルと車両モデルをMLDS（混合論理動的シ ステム)を用いて線形表現し，混合整数計画問題として定式化することで最適なアシスト 量を計算する方法を提案する．これにより，“ドライバが危険だと感じる瞬間”や”加速を 行いたいと感じる瞬間”を車両モデルが素早く感知し，評価関数に基づいてアシスト量を 加えることで，ドライバにとって適切な支援を行うことが可能となる．提案するアシスト システムを運転シミュレータ上に実装し，その有用性を検証する．

1.3

本論文の構成

本論文は，以下のように構成される． 第2章では，本研究で用いる実験装置であるドライビングシミュレータ(DS)および実 車の特長を比較し，それぞれの構成について説明する． 第3章では，ハイブリッドシステムモデルのPWARXモデルおよびその確率的近似モ デルであるPrARXモデルを導入し，そのシステム同定手法を述べる．次に，これらを用 いてドライバの前方車追従行動表現する方法を示す． 第4章では，最適化計算によるアシスト量の決定手法を説明する．また，ドライバモデ ルや自車モデルを最適化計算の制約条件として用いる手法を説明する．

上で得られた運転データをもとにシステム同定を行った結果を示し，前方車追従行動の定

量的説明を行う．また，PWARXとして表現された運転行動モデルをアシストシステム

に組み入れて最適化計算を行った結果を示し，アシストの有無による結果の違いを検証す る．

実験装置・条件

本章では，本研究で用いる実験装置であるドライビングシミュレータ(以下DS)および 実車の特長について述べ，その後，それぞれの実験装置の構成について述べる．また，DS を用いた実験の内容について説明する． 2.1節では，データを収集する際のDSを用いる利点と，実車を用いる利点について データの信憑性等の観点から比較する． 2.2節では，DSの構成や，実験の条件，手順について説明する． 2.3節では，実車実験に用いた実験装置の構成や，実験の条件，手順などについて説明 する．

2.1

実車およびドライビングシミュレータ（DS）の特長

一般に，DSとは実際に道路を走行するのではなく，視覚情報，加速度情報および音情 報を疑似的に創り出し，その仮想世界で走行する装置のことをいう．そのため実車では行 えないような危険，もしくは特殊な条件での走行データの収集を行なうことが可能となる が，一方で，実際にアシスト系を運用する場合は，シミュレータだけではなく，実車を用 いて安全性の検証などを十分にすべきである．以下に，実車を用いた場合のデータ収集と DSを用いた場合のデータ収集の特長を述べる. • 実車でのデータ収集 – 得られるデータの信憑性が高い． – 自車の挙動を直に体感することが出来る． – 多様性が低い． – データを習得するために各種センサを車両に取り付ける必要がある． – 汎用性が低い． – 得られる運転情報はプロのドライバが主流である．本来欲しい一般ドライバも

– 再現性が低い． – 同じシチュエーションの実験を行うことが困難である． • DSでのデータ収集 – 実車に比べて信憑性が低い． – 自車の挙動を直に体感することは出来ない． – シミュレーションにおける環境と実際の環境，自動車モデルと実車を比較する 必要がある． – 多様性が高い． – 車のモデルを変えることで，あらゆる車種にも対応が可能である． – 汎用性が高い． – 被験者の安全が保証されており，一般のドライバの運転情報も収集可能である． – 再現性が高い． – コンピュータで処理しているので，同じ走行環境，条件で何度でも同じ実験が 可能である．

2.2

シミュレータ実験の内容

2.2.1

概要

本研究で使用するDriving Simulator (DS) のシステム構成を図2.1に示す．本システ ムは，3面スクリーンによる前方180◦の視覚呈示が可能な運転環境シミュレータであり， 以下のような特徴を持つ． • 3面，180◦スクリーンによる広い画像呈示 • 3面ミラーと後面による広い画像呈示 • 実車インパネを使用することによる実車の再現 • 反力呈示可能なステアリング • 油圧システムのブレーキを持つペダル操作系 • 5.1chサラウンドシステムを使用した音響 • 物理シミュレータを使用することによる実車の挙動の再現 そして，シミュレータに使用しているPCの役割はそれぞれ • PC1：制御系(車両ダイナミクスの計算部) • PC2：スクリーン表示演算用 • PC3：最適アシスト量演算用

MILP SOLVER Memolink 図2.1 DSの構成図 となっている．また，それぞれのPCはInterface社製Memolinkによってデータの受け 渡しが可能となっている．制御系と表示演算系の構成の詳細を表2.1に示す．制御系に用 いている接続デバイスは，ステアリング，アクセル，ブレーキ，スピーカである．シート は三菱自動車エンジニアリング社から購入した実車のものを運転席内に設置した．ステア リングは市販のものを用いており，左右にそれぞれ450◦ 回転する．アクセル,ブレーキ はステアリングとセットになっているものを用いている．インパネ内には速度メータが設 置されており，PC1により制御することでドライバは実際の走行速度を知ることができ る．ただし，本実験では速度メーターを使用していない． 表2.1 DSの構成 制御系(PC1) 表示演算用(PC2) 最適化演算用(PC3)

OS Windows XP Windows XP Windows XP

professional professional professional

CPU Core2 DUO E8500 Core i7 950 Pentium 4

3.16GHz 3.07GHz 3.00GHz メモリ 3.00GB RAM 3.00GB RAM グラフィック NVIDIA GeForce ボード GTX 295 接続デバイス ステアリング プロジェクタ ×3 アクセル （前方3画面） ブレーキ スピーカ 速度メータ

図2.2 DSの運転席 （注：画面ははめ込み合成） 以下に，それぞれのPCの詳細な構成とその働きを述べる．

2.2.2

制御系

(PC1)

制御系のPCには，統合制御プログラムが実装されており，車両のダイナミクス計算や スピーカからエンジン音の出力，USBを介してのステアリング，アクセル，ブレーキ操作 情報の収集を行う．自動車への入力装置としてステアリング，アクセル，ブレーキをシー トとともに持ち込みコクピットを形成している(図2.2参照)． CarSim 車両ダイナミクスはCarSim(ver 5.15)と呼ばれる車両運動ダイナミクス計算ソフト ウェアを用いている．CarSimには様々な入力があり，これを自由に変えることでその条 件を満たす車両の挙動をシミュレーションすることができる．またタイヤを1本ずつ独立 に計算しているため，歩道への乗り上げ等もCarSimを用いることによって忠実に再現し ている． 統合制御プログラム CarSim とその他デバイスのデータの橋渡しのためにC++Builder にて作成した． CarSimとは共有メモリでつながっており，CarSimは統合制御プログラムが共有メモリ

図2.3 統合制御プログラムのコントロール画面 に書き込んだ情報を基に車両ダイナミクスを計算する．その計算結果は共有メモリに書き 込まれ，統合制御プログラムがメモリンクを用いて表示系のPC，および速度メータ制御 用のPCに送信する．またこのプログラムは対向車などの他車の制御もしており，図2.3 に示すコントロール画面により他車の位置・速度等の制御を行なっている．エンジン音は CarSimの計算結果であるエンジン回転数に応じて統合制御プログラムがDirectXを用い てエンジン音を生成し，スピーカより出力する．

2.2.3

スクリーン表示演算用

PC(PC2)

前方スクリーン映像の計算にPCを１台(PC2)割り当て，表示系プログラムを用意し た．仮想街環境の3次元モデル，車両の3次元モデル，各3次元モデルに描くテクスチャ 画像をあらかじめ用意しておく．表示系プログラムは，Memolinkを通してPC1から受 け取ったシミュレーション結果を基に，現時刻における3次元モデルを組み立てる．組み 立てた3次元モデルを基に，前方スクリーンに投影する映像を作成し，プロジェクターを 通して前方スクリーンに表示する．

2.2.4

最適アシスト量計算用

PC(PC3)

後述する最適化計算を行うために，PCを1台(PC3)割り当てた．RS-232Cを用いて PC1から環境情報を取得し，後の章で示す最適化問題をNU-OPTを用いて解き，PC1 へとその解を転送する．

2.2.5

走行環境

(DS)

本研究で用いた走行環境の全体像を図2.4に示す．

図2.4 走行環境の鳥瞰図 想定環境として，片側二車線で，無限遠に伸びる直進*1_{のみの高速道路環境を用意し} た．また追越車線には車は走行せず，対向車線には85km/h(走行車線)と105km/h(追越 車線)で約200m間隔*2_{で常に走行している．} 前方車は，全試行について図2.5に示した速度パターンに従い走行し，自車の10m前 方に存在した状態でスタートする．

2.2.6

実験手順

(DS)

DSを用いた実験では，まずは平常車における前方車追従の行動モデルを構築し，その 後にアシストシステムを搭載した車を運転して，アシストの有無による違いを考察する． 実験の手順を以下に示す． 1. 個人情報の口頭質問 2. DS説明・実験走行内容の説明 3. 平常車練習走行 4. 平常車走行 5. 平常車運転モデル構築 6. アシスト車実験 7. 操作感に関する口頭質問 *1_{全長6kmの直進高速道路環境を繰り返し表示している} *2_{走行車線：90m,240m,350m 追越車線：120m,150m,160m,260mの間隔で6台が繰り返し走行}

0 100 200 300 400 500 600 0 5 10 15 20 25 30 35 40

Velocity of forward vehicle(m/s)]

time(sec) 図2.5 前方車の走行速度パターン 上記の実験手順において実際に行うことを順に記載する． 個人情報の口頭質問 実験前に，性別・年齢・視力・免許取得からの年月・年間走行距離・運転頻度等につい て質問する． DS説明・実験走行内容の説明 DSの操作の説明と，今回用いる実験環境と走行手順について説明し，実験内容を理解 してもらう． 平常車走行 平常車を用いて10分間の前方車追従を行ってもらう．前方車は前節でも記述した通り， 図2.5の速度パターンに従って走行する． 平常車運転モデル構築 あらかじめ用意した解析プログラムにより，平常車運転モデルを構築する．ここで用い る解析理論は第3章に記述する．

アシスト車実験 前手順において構築された運転モデルを用いて被験者の運転行動に応じたアシストシス テムを構築し，アシスト車を用いて10分間の前方車追従を行ってもらう．前方車は図2.5 の速度パターンに従って走行する．実際に実験で用いたアシスト機能については第4章に 記す． アシスト車アンケート記入 実験終了後，運転の感覚について以下の質問に答えてもらった． 1. 運転中に違和感があったか 2. 運転しやすかったか 3. アシストのタイミングは適当であったか

ハイブリッドシステムモデルとその

同定手法

本章では，ハイブリッドシステムモデルの一種であるPWARXモデルおよびその確率 的近似モデルであるPrARXモデルを導入し，そのシステム同定手法を述べる． 3.1節では，古典的なシステム同定モデルである線形自己回帰(ARX)モデルについて 述べる．ARXモデルはハイブリッドシステムモデル中の1つの動作モードに対応する． 次に3.2節では，複数のARXモデルから成る区分的ARX(PWARX)モデルについて 述べ，そのシステム同定に必要な K-meansクラスタリングとサポートベクターマシン (SVM)のアルゴリズムについて詳しく説明する． 3.3節では，モードの切り換わりを確率的に表現する確率重みARX(PrARX)モデルを 導入し，その特性を利用したオンラインパラメータ推定手法を述べる． 最後に3.4では前方車追従タスクにおける人間の運転行動をハイブリッドシステムとし て表現する方法について述べる．

3.1

ARX

モデル

ARX(Auto Rgressive and eXogenous)モデルとは，過去の入出力信号の線形和として

次の出力値が与えられる線形モデルである．まず以下のような差分方程式を考える． yk+ a1yk−1+· · · + anyk−n= b1uk−1+· · · + bmuk−m+ ek (3.1) 式 (3.1) にお いて，yk はサン プリン グ時刻 k における システム の出力 ，uk はサン プリング時刻 k における入力，n, m はそれぞれ過去の出力，入力の数である．また a1, a2, . . . , an, b1, b2, . . . , bmはそれぞれ出力，入力のパラメータを表しており，ekは式誤 差で白色雑音とする． ここで a1, a2, . . . , an, b1, b2, . . . , bm をパラメータベクトルとして以下のように定義

θ = [a1, . . . , an, b1, . . . , bm]T (3.2) また，回帰ベクトルを xk = [−yk−1, . . . ,−yk−n, uk−1, . . . , uk−m]T (3.3) と定義すると，出力ykは次式のように表現できる． yk = θTxk+ ek (3.4) ARXモデルは最も基本的なシステム同定モデルの1つであり，そのパラメータ同定も 最小二乗法にもとづいて容易に行える．しかしながら，人間の自動車運転行動などの非線 形性を示す挙動のモデリングには適さない．次節以降では，複数のARXを組み合わせる ことでより複雑な挙動を表現することが可能な2つのモデルについて議論する．

3.2

PWARX

モデル

3.2.1

PWARX

モデルの定義

PWARXモデルは離散的な状態遷移と連続的なモデル(ARXモデル)を組み合わせた モデルで，ハイブリッドシステムの中でも代表的な表現手法の一つである[20]．具体的に は，前節で導入した式(3.3)の回帰ベクトルがs個の多面体領域Ui={ui, θi∈ R(n+m)} のどの領域に属するかによってダイナミクスが切り替わるもので，ARXの元の式である 式(3.4)は yk =            θT1xk+ ek if xk ∈ U1 θT 2xk+ ek if xk ∈ U2    ... θT sxk+ ek if xk ∈ Us (3.5) と表現される．

3.2.2

PWARX

モデルによるパラメータ同定の流れ

PWARXモデルのパラメータ推定手法は様々である[21]．例えばFerrariら[20] によ る確定的クラスタリングによる手法，Nakadaら[22]による統計的クラスタリングに基づ く方法，A.Lj.Juloskiら[23]によるベイズ推定に基づく方法そしてBemporadら[24]に よる確定的誤差バウンドに基づく方法などがある．これらの手法は一長一短であり，同定 対象によって使い分ける必要がある．本研究では前方車追従行動の同定問題を扱うため，

事前知識のない物理システムに対するパラメータ同定に適した確定的クラスタリングよ る手法を採用した．クラスタリングとはデータ解析の諸技法の中で，外的基準なしに自 動的に分類する方法，いいかえれば，データ以外にあらかじめ基準を設定することなく， データの集まりをいくつかのグループに分ける方法の事を言う．以降，本研究で用いた K-meansクラスタリングの具体的な理論について説明しながら，PWARXモデルのパラ メータ同定の大まかな流れを述べる．図3.1にその概略図を示す．

Clustering the feature space (K-means Method)

(c)

x y

Transform to sample space

(d)

x y Pi mi :Mean Value feature vector LDs i i :Estimated Parameter

Transform to feature space

(a)

(b)

図3.1 クラスタリングの流れ

3.2.3

K-means

法を用いた

PWARX

モデル同定

K-means法ではクラスターの数をあらかじめ指定し，特徴量をs個のクラスターに分 割する．本研究では運転行動のダイナミクスを考慮するため，以下に示す手順に従いクラ スタリングを行った． 1. 各データ点の特徴量を，その周辺の局所データ集合のダイナミクスを用いて計算 する 2. 各データ点の持つ特徴量の信頼度指標を計算する 3. 特徴量とその信頼度指標を用いてクラスタリングを行う

特徴量の生成 入出力から特徴量を抽出し，特徴量に基づきデータのクラスタリングを行う．しかし特 徴量の選定方法は様々であり，解析対象によって異なる．本研究では前方車追従行動のダ イナミクスを考慮したモデル化を目指すため，文献[20]で紹介された局所集合(ローカル データセット)による特徴量の抽出を行う．以下にその手順を示す． 1. 実験により得られた入出力データは，単位がそれぞれ異なり，そのまま用いること は望ましくない．このため式(3.6)のように正規化を行う． ¯ di= di max|di| (3.6) 2. 正規化された入出力信号から回帰ベクトルxkを得る． 3. xj(j = 1, 2, . . . , n)に対し最も距離の近いc個のデータ(x1LDs, xLDs2 ,· · · , xLDsc ) を集め，局所データ集合LDsj = (uLDsj, yLDsj)を作る． 4. 局所データ集合LDsj より、θ LDsj j を以下のように計算する． θLDsi j = (Φ T jΦj)−1ΦTjyLDsj Φj = [ xLDsj 1 x LDsj 2 · · · x LDsj c ]T (3.7) これは，yLDsj を出力系列，uLDsj を入力系列として，最小二乗推定を行ったとき に得られるパラメータである． 5. 局所データ集合LDsj のxLDsj の平均mj を計算する． mj = 1 c ∑ xLDsj, j = 1, . . . , N (3.8) 以上の操作によって導出されたθLDsj _とm j から特徴量 ξj = [(θLDsj ) T , mTj] T , ∀j = 1, . . . , N (3.9) を定義する． 特徴量の信頼度指標の計算 上記の手順に従い特徴量を用いてクラスタリングを行う．クラスタリングのアルゴリズ ムを以下で説明する． 1. クラスタリングを行う際，一般的にはユークリッド距離を用いるが，アウトライア の影響を受けやすいという欠点を持つ．そこで距離の定義には特徴量の共分散を用 いたマハラノビス距離を取り入れる．式(3.7)のΦj を用いてθLDsj の共分散Vj を求める．

Vj = SSRj c− (n + 1)(Φ T jΦj)−1 SSRj = yTLDsj(I− Φj(Φ T j Φj)−1ΦTj )yLDsj (3.10) SSRj はローカルデータセットにおけるパラメータベクトルの信頼度を表して いる． 2. 同様にして式(3.8)のmj の共分散Qj を求める． Q_j = ∑ x)∈LDsj (xLDs− mj)(xLDs− mj)T, mj = 1 c ∑ (x)∈LDs xLDs, j = 1, . . . , N (3.11) 3. Vj とQj を分散Rj としてまとめる． Rj = [ Vj 0 0 Qj ] (3.12) 4. マハラノビス距離を用いた式(3.15)の評価関数を用いて，特徴ベクトルをs個の 部分集合Diにクラスタリングをする． 但し，µj は後述する部分集合Diの重み 付き重心である． J ({Di}si=1,{µi}si=1) = s ∑ i=1 ∑ ξj∈Di ||ξj − µi||2_R−1 j = s ∑ i=1 ∑ ξj∈Di (ξj− µi)TR−1j (ξj − µi) (3.13) 重み付きk-means法のアルゴリズム クラスタリング手法は，データマイニングとして多く用いられる一般的な手法である． データ点間の距離を用いることでデータ点間の類似度を評価し，類似度の高いデータ点群 によるクラスタを構築する．k-means法はクラスタリング手法の中でも，分割的なクラス タリングと呼ばれ，全データ点を既知のクラスタ数N に分割するための手法である． 本論文では，通常のk-means法の距離計算において，信頼度指標の重み付き距離を導 入した重み付きk-means手法に基づいて各データ点の属するモードを判別する．本手法 は以下に示すようなアルゴリズムに基づいて行われる． Step.1 部分集合の中心（以下，代表点）µ(0)_i (i = 1, 2, . . . , N )を乱数により初期化す る．

← 0 Step.3 以下，繰り返す． Step.4 評価関数JKM({D (r) i }si=1,{µ (r) i }si=1)を最小にするような集合D (r) i を計算す る． ただし， JKM({Di}si=1,{µi}si=1) = s ∑ i=1 ∑ ξj∈Di ||ξj− µi||2_R−1 j (3.14) = s ∑ i=1 ∑ ξj∈Di (ξj − µi)TR−1j (ξj − µi) (3.15) Step.5 次式を満たすように代表点を更新する．    ∑ j:ξj∈Di(r) R−1_j    µ(r+1) i = ∑ j:ξ_∈Di(r) R−1_j ξj (3.16) Step.6 r← r + 1とする． Step.7 書き条件を満たすまで上記を繰り返す 代表点を更新する前後の評価関数が一致する JKM({D (r−1) i } s i=1,{µ (r) i } s i=1) = JKM({D (r−1) i } s i=1,{µ (r−1) i } s i=1) (3.17) Step.8 得られたD(r)_i ，µ(r)_i を最適解とし，終了する． Di∗= D (r) i , µ∗i = µ (r) i (3.18) Step.9 評価関数JKM が最小となるような部分集合Di∗が得られたので，以下の式に 従って各データ点が所属するモードstを決定する st = i s. t. ξt∈ Di∗ (3.19) Step.10 終わり 本アルゴリズムにおける評価関数JKM({Di}si=1,{µi}si=1)は，各特徴量ベクトルξj か らそれに最も近い代表点µiまでの，R−1_j によって重み付けられた特徴量空間上の二乗距 離の和となっており，代表点を用いて各データ点を表現する際の歪み量に対応する．すな わち，各データ点の特徴量が精度よく求まってさえいれば，この評価関数が小さいほど代 表点による各データ点の近似精度が高く，精度を落とすことなくPWARX表現が得られ ることとなる．また，Step 5では，R_j−1によって重みづけられた重み付き平均を用いる ことにより，信頼度の低い点であるアウトライアの影響を低減することができる．

以上の手法により，データ中の各モードに対して最小二乗推定を施すことで式(3.5)に おける未知パラメータ，θi(i = 1, 2, 3, 4)を求めることが出来る．

3.2.4

サポートベクターマシン

(SVM)

本研究では，上述のクラスタリング手法によって分けられた各基本タスクに対し，タス クの切り替え条件を求めるために，SVMを用い推定を行った．以下にSVMについて説 明を加える． SVMはパターン認識における学習モデルの1つである.基本的に2つのクラスを識別 する識別器を構成するための学習法であり，文字認識など多クラスの識別器を構成するた めには，複数のSVMを組み合わせるなどの必要がある. SVMの特徴は，学習データか ら「マージン最大化」という基準で識別関数のパラメータを決定する，という点である （図3.2の 1 ||w||がマージンである）. 入力ベクトルをx，識別関数を決定するパラメータ をw,hとし，識別関数を以下のように定義する. f (x) = sign(g(x)) (3.20) g(x) = wTx− h (3.21) ここで関数sign(u)は，u > 0のとき1，u≤ 0のとき−1をとる. 学習データはn個与えられているとしxi(i = 1, 2, . . . , n)と表す. これらのデータを2 つのクラスX1, X2に分類するとし，学習データ集合に対しg(x)がつぎの条件を満たす ようなパラメータをもとめる. g(xi) = { wTxi− h ≥ 1 if xi∈ X1 wT_x i − h ≤ −1 if xi∈ X2 (3.22) wTx-h=1 wT_{x-h = -1} 1/||w|| 1/||w|| 図3.2 SVMの概略

ti= { 1 if xi∈ X1 −1 if xi∈ X2 (3.23) これを用いて式(3.22)を書き直すと次式のようになる. ti(wTx− h) − 1 ≥ 0 (3.24) ここで（3.24）式を満たすようなパラメータw, hが存在すれば，図3.2のようにg(x) = 1 とg(x) =−1の2枚の超平面で学習データが完全に分類され，2枚の超平面の間には学 習データが1つも存在しないことを示している. このとき，識別平面とこれらの超平面との距離（マージンの大きさ）は 1 ||w||となる.し たがって，マージンを最大化するパラメータw, hを求める問題は，制約条件 ti(wTx− h) − 1 ≥ 0 (3.25) のもとで，目的関数 L(w) = 1 2||w|| 2 _(3.26) を最小化するパラメータを求める問題として置き換えられる. ソフトマージン 上述のSVMは，学習データが線形分離可能にのみ利用できる．しかし，実際の問題で 線形分離可能である場合は稀である．したがって，線形分離不可能な場合に対して，図 3.3のように多少の識別誤りを許すように制約を緩める方法がとられる．これを「ソフト マージン」と呼ぶ．ソフトマージンでは，マージン _||w||1 を最大としながら，いくつかの サンプルが識別面を超え，別のクラスに入ってしまうことを許す．反対側にどの程度入り 込んだかという距離を，パラメータζi(≥ 0)を用いて，_||w||ζi と表すと，その和である N ∑ i=1 ζi ||w|| (3.27) はなるべく小さいことが望ましい．これらの条件から最適な識別面を求める問題は，制約 条件 ζi≥ 0, ti(wTxi− h) ≥ 1 − ζi, (i = 1, . . . , N ) (3.28) のもとで，目的関数 L(w, ζ) = 1 2||w|| 2_{+ γ} N ∑ i=1 ζi (3.29)

を最小とするパラメータを求める問題と置き換えられる．ここでパラメータγは，マージ ンの大きさとはみ出しの程度とのバランスを決める定数である． 本研究では，ソフトマージン法を用いて線形分離不可能な場合に対処したうえで分離面 のパラメータを求めた． wT_x-h=1 wT_{x-h = -1} 1/||w|| 1/||w|| 図3.3 ソフトマージン

3.3

PrARX

モデル

3.2節で紹介したPWARXモデルは人間の行動を確定的な境界面で分割しているが, 本節ではPWARXと多項ロジスティック回帰の考え方を融合した PrARX(Probability-weighted ARX)モデルについて説明する.

3.3.1

PrARX

モデルの定義

PrARXモデルの出力は次式で与えられる． yk = s ∑ i=1 Pi(xk)θiTxk (3.30) PrARXモデルはs個のARXモデルが確率Piによって重みづけされる形となっている． Piはモードiが選択される確率を表し，以下のようなソフトマックス関数によって与えら れる. Pi(x) = exp(ηT i x) ∑s j=1exp(ηTj x) (3.31) ηs= 0 (3.32) 確率Piはパラメータηによって決定され,多項ロジスティック回帰モデルにおける出力 と同じ形式となっている. このように，PrARXモデルは多項ロジスティック回帰モデル とPWARXモデルを融合した形となっている． また，PrARXモデルの確率重み関数Piをもとに，次式によって確定的なモード条件 を得ることができる． x∈ Ui i = argmaxPi(x) (3.33) これにより，PrARXモデルを機械的にPWARXモデルへ変換することができる.

3.3.2

モデル例

PrARXモデルの例として,3モード1変数のPrARXモデルのサンプルを以下に示す. このモデルのパラメータは, θ1= [ 0.5 5 ] , θ2= [ −0.1 3 ] , θ3= [ −0.4 15 ] η1= [ −3 45 ] , η2= [ −1.5 30 ]

0 5 10 15 20 25 30 0 5 10 y 0 5 10 15 20 25 30 0 0.5 1 Probability (P i ) u 図3.4 PrARXモデルの例 のように設定した. 図のように,赤の破線で示された3つのARXモデルが確率Pに従って,青の実線のよ うに滑らかに切り替わっていることが確認できる. このことから,PrARXモデルは,ダイ ナミクスを表すパラメータθ,判断を表すパラメータηで決定され,モードの切り替わりを 滑らかにしたハイブリッドシステムとみなすことができる.

3.3.3

PrARX

モデルのパラメータ同定法

PrARXモデルのパラメータは，以下の誤差関数を最小にするθ，ηを求めることに よって得られる．  = 1 N N ∑ k=1 kyk− ˆykk2, (3.34) ˆ yk= s ∑ i=1 Pi(xk)θTixk (3.35) この目的関数はモデルによって推定された出力yˆkと実際に観測された出力yk との誤差 のノルムの二乗和を表している. 本研究ではこの最適化問題を最急降下法を用いて解く. 最急降下法は以下の更新式に従い,得られた最急降下方向にパラメータを更新することに よって最適解を求める手法である.

θ_i(t+1)= θ(t)_i − α ∂ ∂θ_i(t) (3.36) η(t+1)_i = η(t)_i − β ∂ ∂η_i(t) (3.37) 目的関数のパラメータによる偏微分は, ∂ ∂θi =−1 N N ∑ k=1 ekPi(xk)xk (3.38) ∂ ∂ηi =−1 N N ∑ k=1 2ekPi(xk)xk(θTi x− ˆyk) (3.39) ここで,θ_i(t),η(t)_i はt回更新後のパラメータである. α,βは微小な正の数を任意で使用する. また,初期値依存性の問題から局所最適解を回避するため多くの初期値の組を用意し,検 証を行う必要がある. PrARXモデルの利点として,この最急降下法というアルゴリズム一つで,ダイナミクス に関するパラメータθiとモードの境界を決定するパラメータηiが同時に同定できるとい う点がある.

3.3.4

パラメータの逐次更新手法

人間の操作・判断特性は，支援システムの機能やタスクへの慣れなどによって時々刻々 変化していくと考えられる．したがって，PrARXモデルのパラメータθ，η も常にそれ に合わせて最適な値に調整することが望ましい. 前述の最急降下法を用いたパラメータ同 定方法は局所最適解を回避するために大量の初期値が必要となり，結果パラメータ更新に 時間がかかり，実時間でのモデルの最適化が難しいという問題がある．支援システムの機 能やタスクへの慣れに起因する操作・判断特性の変化は,時刻とともに徐々に変化してい くと仮定すると，現時刻のパラメータを何らかの規則に従って修正することで，時刻に対 応した最適パラメータを求めることが可能になると考えられる．そこで，PrARXモデル のパラメータが最急降下法に基づくパラメータ更新によって決定されていることを利用 し，現時刻のパラメータを最急降下法に基づくパラメータ更新の初期値として用いること で計算時間の短縮と，局所最適解の回避を両立させる手法で解析を行う． 逐次更新の手順 1. 取得したデータの内，現時刻から一定のステップ数τ 前までのデータのみを用い， 最急降下法を用いてそのデータ点内で最適なパラメータを計算する．このとき，最 急降下法に用いる初期値は前時刻でのθ,ηを用いる. 2. 前回の開始点から1ステップだけずらし,そこからτ ステップ前の解析区間につい て,現時刻でのモデルパラメータを初期値として,最急降下法を行う.

0 500 1000 1500 2000 0 0.05 0.1 0.15 0.2 θ 11 0 500 1000 1500 2000 −0.04 −0.02 0 k θ 12 0 500 1000 1500 2000 −0.1 −0.05 0 0.05 0.1 θ 21 0 500 1000 1500 2000 −0.05 −0.03 −0.01 0.01 k θ 22 (a) θ1 (b) θ2 0 500 1000 1500 2000 0 5 10 η 11 0 500 1000 1500 2000 −5 0 5 k η 12 (c) η1 図3.5 学習用モデルパラメータの時間変化 3. 2.の手順を繰り返す. 逐次更新の適用例 以下では，パラメータの時間変化が既知であるPrARXモデルに対して，その入出力信 号を学習データとして逐次更新手法により求めたモデルパラメータの時間変化を真値と比 較する．学習データとして,式(3.30),(3.32)の形で表現されるモード数2のデータ2000 点を用意する. 入力信号xは式(3.40),モデルパラメータθ,ηは図3.5のように与える. ま た，この例においてはデータ長τ は200とした． φk= [uk1]T, uk = sin ( k 50π ) (3.40) 以上の学習データを逐次更新の手順に従い,解析を行った結果,推定されたパラメータ

0 500 1000 1500 2000 −0.04 −0.03 −0.02 −0.01 0 θ 12 k 0 500 1000 1500 2000 0 0.05 0.1 0.15 0.2 θ 11 0 500 1000 1500 2000 −0.1 −0.05 0 0.05 0.1 θ 21 0 500 1000 1500 2000 −0.1 −0.05 0 0.05 0.1 k θ 22 (a) θ1 (b) θ2 0 500 1000 1500 2000 0 5 10 η 11 0 500 1000 1500 2000 −5 0 5 k η 12 (c) η1 図3.6 逐次更新により推定されたパラメータの時間変化と真値との比較 図3.6に示された推定パラメータの様子から,正解パラメータの変化に対応したパラ メータ推定が行われていることが確認でき,図3.7で示されるように,実際の出力とほぼ相 違ない形で出力の推定が行われている.

3.4

ハイブリッドシステムによる前方車追従行動のモデル化

前方車追従行動をモデル化するための入出力変数を以下に示す．入力信号を以下のよう に定義する． • u1：接近離間状態評価指標(KdB)[dB][19] • u2：前方車と自車の車間距離[m] • u3：前方車と自車の相対速度[m/s]

0 500 1000 1500 2000 −0.02 −0.01 0 0.01 0.02 0.03 0.04 Output k True model Identified model 図3.7 出力信号の比較 表3.1 KdB導出に用いる変数一覧 Hp=前方車の車高 H=網膜に映る前方車の車高 Wp=前方車の車幅 W =網膜に映る前方車の車幅 Sp=前方車の背面積 S=網膜に映る前方車の背面積 D=車間距離 f =焦点距離 次に出力信号は以下で与える． • y：自車加速度[m/s2] ここで接近離間状態指標(以下KdB)とは，伊佐治ら[19]によって提案されたドライバ の状態推定を行うための指標であり，「ドライバは前方車との相対関係の変化を前方車の 視覚的面積変化として認識し加減速操作を行っている」という仮定で定義されている．具 体的にはドライバの網膜上に投影される前方車の見かけ上の面積の時間変化率を用いて， 以下のように導出される．まず変数の定義を表3.1に示す． 前方車に追従走行する場合を考えると，ドライバの網膜上に投影される前方車の見かけ 上の面積Sは式(3.41)で表現される． W = Wp× f D, H = Hp× f D, S = W × H = Wp× Hp× f2× 1 D2. (3.41) この時のドライバの網膜上に投影される前方車の見かけ上の面積の時間変化率は式(3.42) のようにあらわすことができ，これを前方車の面積の時間変化率Kとする． dS dt = d(W × H) dt ∝ d dt( 1 D2) =− 2 D3 × V r =: K (3.42)

図3.8 KdBのイメージ図 この前方車の面積の時間変化率Kは，たとえば車間距離D = 1∼ 100[m]の間で106の オーダで大きく変化する指標となるため，扱いやすいデシベル表示とする．自車の100[m] 前方に存在する相対速度V r = −0.1[km/h]で接近してくる前方車の面積の時間変化率 K0をドライバが面積変化に気づく事ができる最小検出限界と仮定し，このときの値を 0[dB]と定義する(式(3.43))． V r =−0.1 3.6[m/s] D = 100[m] K0=− 2 D3 × V r = − 2 1003 − ×(− 0.1 3.6)' 5 × 10 −8 _(3.43) つまり，前方車の面積の時間変化率K0= 5× 10−8の時のデシベル値を0[dB]とし，式 (3.44)のようにあらわす指標をあらためて接近離間状態評価指標KdBとして定義する． なお，接近離間状態評価指標は前方車が近づいてくるときを正の値，離れていくときを負 の値とした． u1=      10× log10(−κ) if κ < −1 −10 × log10(κ) if κ > 1 0 if − 1 ≤ κ ≤ 1 (3.44) ただし， κ = 2× Vr D3 × 1 5× 10−8 = 4× 10 7_× Vr D3 (3.45) である． また，人間が入力情報を認知してペダルの操作するまでの遅れを1ステップ，ペダルの 操作の後に車両の加速度に反映されるまでの遅れを2ステップと仮定する．すると，前方

車追従行動はPWARXモデルとして以下のように表せる． yk =                        θ11u1,k−3+ θ12u2,k−3+ θ13u3,k−3+ θ14yk−1 if u∈ U1 θ21u1,k−3+ θ22u2,k−3+ θ23u3,k−3+ θ24yk−1 if u∈ U2    ... θi1u1,k−3+ θi2u2,k−3+ θi3u3,k−3+ θi4yk−1 if u∈ Ui    ... θs1u1,k−3+ θs2u2,k−3+ θs3u3,k−3+ θs4yk−1 if u∈ Us (3.46) また，PrARXモデルに関しても同様に表現可能である． モード数sに関しては，一般にsが大きいほど精緻なモデルであるといえるが，同時に パラメータ推定のコストの増大，解釈の難しさ，および後述するアシストシステムにおけ る計算コストの増大を招く．この点を考慮し，本研究では2モードないし4モードのモデ ルを採用することとする．

モデル予測型アシストシステム

本章では，モデル予測型運転アシストシステムについて述べる．モデル予測型アシスト システムとは，モデル予測制御のコンセプトにもとづき，運転行動モデルから予測される 運転者の将来の運転行動に対して最適なアシストをリアルタイムで計算するシステムであ る．個々の運転者のデータをもとに運転行動モデルを同定することで，個人に適合したア シストが実現される． はじめに5.3.2においてモデル予測型アシストシステムの全体像を示す． 次に4.2節において，モデル予測制御の各制御周期で計算するアシスト量最適化問題を 定式化する． 4.3節では，モード分割された人間の運転行動モデルをMLDSで表現する手法につい て説明する． 4.4節では非線形関数であるKdBの線形化について述べる． 4.5節では，MILPの制約条件として用いる，自車モデルおよび前方車モデルについて 説明する． 4.6節では，アシストを加える条件およびアシストの線形不等式表現を説明する．

4.1

モデル予測型運転アシストシステム

図4.1にモデル予測型運転アシストシステムの全体図を示す．アシストシステムは，各 制御周期において後述するアシスト量最適化問題を計算し，得られたアシスト量を運転者 のペダル操作に足し合わせることで運転支援を実現する． 一般にモデル予測制御とは，プラントモデルをもとに最適な制御入力の時系列を計算 し，その先頭要素のみを実際にプラントに印加するという手続きを所定の制御周期でく り返す制御手法である[27]．モデル予測制御の処理の流れを図4.2に示す．環境情報を t, t + 1,· · · と1ステップごとに観測できるとし，現時点の時刻をtとする．まず，現時刻 tにおいて前方車の速度，自車の位置，速度，加速度といった環境情報を計測し，後述す

＋ ＋ ＋ ＋ ＋ ＋ ＋ ＋

Assist optimizer based on MILP

Operation

Assistance

Human

Driver

S

e

ns

or

y

inf

or

m

a

ti

on

Constraints

Environment

and vehicles

Driver model

Physical and

Logical constraints Vehicle models

J

Objective function

Sensory Information

and parameters

Given

図4.1 Model-Predictive Driving Assist System

(a) 時刻tにおける信号 (b) 1ステップ後の時刻t0における信号 図4.2 モデル予測制御の流れ る最適化問題をを解くことによりN ステップ先までの最適アシスト解を求める．得られ たN ステップのアシスト時系列のうち，現時刻に相当するyas[t]だけを次のサンプル時 刻t + 1まで実際に系に加え続ける．1ステップ後の時刻t + 1において環境情報を再び 観測する．このとき，観測値はモデルの不完全さや外乱などの影響のため，時刻tで予測 した値とは必ずしも一致しない．そこで，時刻t + 1を新たに現時刻t0とし，観測値y(t0)

プ最適化問題を解き決め直す．このように，モデル予測制御はステップごとに予測区間， 制御入力区間をスライドし，開ループの最適化問題を解くことにより，制御入力を定める 手法である．

4.2

アシスト量最適化問題の定式化

本節では，アシスト量を決定する最適化問題を定式化する．アシスト量最適化問題はお おまかに以下のような形をとる． Given u1[1], u2[1], u3[1], y[1], Find u1[2],· · · , u1[K], u2[2],· · · , u2[K], u3[2],· · · , u3[K], y[2],· · · , y[K], y+_as[1],· · · , y_as+[K], y−_as[1],· · · , y_as−[K], and (M LDSに用いる論理変数：4.3節参照) which minimize J = K ∑ i=1 [

ω1ζ2[i] + ω2(u3[i]− ˆu3)2+ y−

2 as[i] + y +2 as[i] ] (4.1) subject to (ドライバモデル：4.3節参照)， (アシストに関する制約: 4.6節参照)， (自車モデル，前方車モデル：4.5節参照) ここでu1, u2, u3はそれぞれKdB，車間距離および相対速度であり，yは自車加速度で ある．また，y+_{, y}−_{はそれぞれ加速アシストおよび減速アシストを表す．} この最適化問題は，自車の位置，速度，加速度の現在値（初期値）を入力値とし，ドラ イバモデルや自車モデルなどを用いることで数ステップ後までの出力(自車加速度)を予 測し，ステップごとの最適アシスト量(m/s2)を求める． ここで，評価関数について説明する．ζ は危険度指標を表しており，TTC（Time To Collision）と呼ばれる指標の逆数を用いている（ζ = u3/u2）．Kは最適化の対象となる

時間区間の長さを表す．右辺の第1項は減速アシストに関するコストであり，危険度を低 減する効果を持つ．第2項は加速アシストのコストであり，相対速度を極力目標値に近づ ける効果を持つ．最後に第3項と第4項はアシスト量それ自体を節約するためのコスト である．ω1およびω2は重みであり，ω1に比してω2が大きい場合は，より相対速度を目 標値に近付けようとする効果が強まり加速アシストが働きやすくなり，逆にω2が小さい 場合は危険度を重視し減速アシストが優位となる． アシスト信号は加速度信号として運転者の操作により生じる自車加速度に加算される． したがってアシストが加えられた車両の加速度は，式(4.2)のように表される． ytotal[i] = yi+ yas−[i− 1] + y + as[i− 1] (4.2) 式(3.46)で表されたドライバモデルには論理条件が含まれており，また，式(3.44)で 表わされるKdBは非線形である．このため，上述の最適化問題はそのままの形では解 くことが非常に難しい．そこで，4.3節で述べる方法でこの問題を混合整数線形計画問題 （MILP）と呼ばれる最適化問題に近似的に変形する．MILPに対しては優れたソルバが 開発されているので，これらを用いて効率的に計算が可能である．また，4.4節以降では， KdBの線形化手法や，自車モデル，前方車モデルの設計手法，その他付加的な拘束条件 について説明する．

4.3

MLDS

を用いたドライバモデルの表現

4.3.1

論理式から線形不等式への置き換え

混合論理動的システム(Mixed Logical Dynamical System, MLDS)とは，PWARXモ

デルのような論理条件を含むシステムを，0-1変数と不等式制約条件によって等価的に表 すシステム表現の一種である．本節で使用する変数を表4.1にまとめる． 置き換えの例として，次の論理命題Xを考える． X ≡ {f(x) ≤ 0} (4.3) この命題の真偽を0-1変数δで関連づけることを考える．つまり， { {f(x) ≤ 0} ↔ {δ = 1} {f(x) > 0} ↔ {δ = 0} (4.4) が成り立つとする．式(4.4)を幾何学的に考えると，図4.3のようになる． 図4.3のf (x)について，δを用いた不等式で表わすと， { f (x)≤ M(1 − δ) f (x)≥  + (m − )δ (4.5)

表4.1 置換例の変数一覧 M : f (x)の最大値 m : f (x)の最小値 :微小正数 δ : 0-1変数 z:連続値補助変数 図4.3 MLDSの置換例１ 図4.4 MLDSの置換例２ となる．したがって，式(4.4)で表わされる関連付けは，式(4.5)の線形不等式群と等価 であるといえる． もう一つの例として，0-1変数δと，線形関数f (x)の積 z≡ δf(x) (4.6) により新たな補助変数zを導入する場合に適用することについて考える．式(4.6)は， { {δ = 0} → {z = 0} {δ = 1} → {z = f(x)} (4.7) と等価である．式(4.7)は，幾何学的に考えると図4.4のようになる． 図4.4のzについて，f (x)およびδを用いた不等式で表わすと，          z≤ Mδ z≥ mδ z≤ f(x) − m(1 − δ) z≥ f(x) − M(1 − δ) (4.8) が成り立つ．したがって，非線形な表現である式(4.6)が，線形不等式群である式(4.8) と等価であるといえる． ここで，MLDSの簡単な例を紹介する (文献[25] から抜粋)．次のように，連続変数

x(t)のとる値によりダイナミクスが異なるシステムを取り上げる． x(t + 1) = { 0.8x(t) + u(t) if x(t)≥ 0, −0.8x(t) + u(t) if x(t) < 0, (4.9) ここで，m ≤ x(t) ≤ M とする．つぎに，x(t)≥ 0 の条件は，表現変換として論理変数 δ(t)を用いることで _{ [δ(t) = 1]↔ [x(t) ≥ 0] [δ(t) = 0]↔ [x(t) < 0] (4.10) と関連付ける．これにより，場合分けを含む式(4.9)を，単一式(4.11)のように書き換え ることが可能となる． x(t + 1) = 1.6δ(t)x(t)− 0.8x(t) + u(t) (4.11) これで式(4.11)は単一式になったが，その前提となる式(4.10)が新たに導入されてし まった．そこで，式(4.10)を式(4.5)と同様に { −mδ(t) ≤ x(t) − m −(M + )δ ≤ −x(t) −  (4.12) と不等式表現に書き換える．ここで，は微小正数である． また，式(4.11)にはδ(t)x(t)という非線形項を含むので，これにも対処する必要があ る．そこで，z(t) = δ(t)x(t) と定義される連続値補助変数を用いると，式(4.8)と同様に          z(t)≤ Mδ z(t)≥ mδ z(t)≤ x(t) − m(1 − δ(t)) z(t)≥ x(t) − M(1 − δ(t)) (4.13) となり，これらの制約を組み合わせることで，式(4.9)は次のように表現可能となる． x(t + 1) = 1.6z(t)− 0.8x(t) + u(t) (4.14) 以上をまとめると，式(4.9)のシステム方程式は，線形等式表現である式(4.14)と線形 不等式表現である式(4.12), 式(4.13)の組み合わせで記述できることになる．

4.3.2

MLDS

システムの一般形

MLDSシステムは次のような状態方程式と線形不等式によって表される．      x(t + 1) = Atx(t) + B1tu(t) + B2tδ(t) + B3tz(t) y(t) = Ctx(t) + D1tu(t) + D2tδ(t) + D3tz(t) E2tδ(t) + E3tz(t)≤ E1tu(t) + E4tx(t) + E5t (4.15)

し，δ(t)は，離散値入力あるいは離散値補助変数を表している．z(t)は連続値補助変数で ある．A，B，C，D，Eは適当な次数の行列であり，添え字tは時変を意味している． 結局,式(4.15)は通常の状態方程式表現にδ，zといった補助変数が加わり，さらにδ，z を一意に定めるために第3式の線形不等式が加わった表現方法となっている．後の計算負 荷を小さくするためには，第3式をできるだけ簡潔に表現する必要がある．

4.3.3

運転行動モデルの

MLDS

表現

この節では，式(4.16)に示すシステムをMLDSで表現する手法を示す．なお，5.3.1節 で記したように，計算負荷を小さくするために運転行動モデルのMode数を2として表現 する．以下，iは0.2秒ごとのステップ数を表す． {

y[i] = θ11· u1[i− 3] + θ12· u2[i− 3] + θ13· u3[i− 3] + θ14· y[i − 1] if u ∈ [ModeA] y[i] = θ21· u1[i− 3] + θ22· u2[i− 3] + θ23· u3[i− 3] + θ24· y[i − 1] if u ∈ [ModeB]

(4.16)

モードの切り替え条件を表す分離面の線形式を式(4.17)に示す．

S1[i] = η11· u1[i− 3] + η12· u2[i− 3] + η13· u3[i− 3] + η14· y[i − 1] + η15 (4.17)

S1[i] > 0とS1[i]≤ 0でModeが分離されるので，表現変換として論理変数δ[i]を用いる

ことで _{ [δ1[i] = 1]↔ [S1[i] > 0] [δ1[i] = 0]↔ [S1[i]≤ 0] (4.18) と関連付ける．これを式(4.5)と同様にして不等式表現に書き換えると， { S1[i]≤ S1,M· δ1[i] S1[i]≥  + (s1,m− ) · (1 − δ1[i]) (4.19) となる．したがって，式(4.18)で表わされる関連付けを，式(4.19)の線形不等式群で表 現することが可能となる．式(4.16)を，δ1を用いた表現に書き直してみる． {

y[i] = θ11· u1[i− 3] + θ12· u2[i− 3] + θ13· u3[i− 3] + θ14· y[i − 1] if δ1[i] = 1 y[i] = θ21· u1[i− 3] + θ22· u2[i− 3] + θ23· u3[i− 3] + θ24· y[i − 1] if δ1[i] = 0

(4.20) ここで，0− 1変数δと変数u, yの積zを次のように定義する．         

z11[i] = δ1· u1[i] z12[i] = δ1· u2[i] z13[i] = δ1· u3[i] z14[i] = δ1· y[i] z21[i] = (1− δ1)· u1[i] z22[i] = (1− δ1)· u2[i] z23[i] = (1− δ1)· u3[i] z24[i] = (1− δ1)· y[i]