ドライビングスティフネス同定を用いた 電気自動車のスリップ率制御に基づくヨーレート制御 狩野岳史 藤本博志 （横浜国立大学） Yaw-rate Control Based on Slip-ratio Control with Driving Stiffness Identification for Electric Vehicle T. Kanou, H. Fujimoto (Yokohama National University) Abstract– On low mu terrain, cooperative control of DYC with traction control can enhances performance than simple DYC. Traction control based on slip-ratio control needs driving stiffness of the terrain. In this paper, we propose the yaw-rate control based on slip-ratio control with driving stiffness identification in real time. Real time identification of driving stiffness achieves accurate driving force control for yaw-moment generation which is demanded by DYC controller. The proposed method is verified by the simulation and the experiment. Key Words: Electric Vehicle, Yaw-rate Control, Traction Control, Slip-ratio Control, Driving Stiffness 1 序論 近年環境問題への関心の高まりから，走行時に直接 CO2 を排出しない電気自動車が注目を集めている。し かし電気自動車の利点はこれだけに留まらない。モータ を動力とする電気自動車には以下のような特徴がある。 • 内燃機関に比べてトルク応答が速い • 発生するトルクを正確に把握できる • インホイールモータにより分散配置が可能 モータによってもたらされるこれらの利点によって，電 気自動車は高い制御性を持つ [1]。 これまでの著者らの研究グループでは，インホイー ルモータを持つ電気自動車に対して，低 µ 路において 車輪の空転を防ぐトラクション制御 (アンチスリップコ ントロール: ASC) や，走行中の車両のヨー運動をノミ ナル化するヨーレート制御 (ヨーモーメントオブザー バ: YMO) を提案し，シミュレーションや実験により その有効性を確認した [2]。また YMO に基づいたコー ナリングスティフネス推定法を提案し，さらに ASC と 組み合わせることで低 µ 路においても外乱の影響を抑 えられることも確認した [3]。 低 µ 路においては，車輪が空転しやすいため，ヨー レート制御器の要求するヨーモーメントが左右のイン ホイールモータの駆動力差により十分に発生できない 状況が起こりうる。このような時，トラクション制御 により車輪の空転を抑制することでヨーレート制御の 性能の維持が期待できる。 文献 [4] では車輪の空転によるヨーモーメントの飽和 をプラントの操作量飽和と捉え，アンチワインドアッ プ制御による改善が検討されている。 また文献 [5] では四輪の横力，制駆動力，接地荷重を 自由に制御可能な車両において，各輪の最適な µ 利用 率を実現する制御系の研究がなされている。 本稿ではドライビングスティフネス同定値を用いた スリップ率制御に基づくヨーレート制御法を提案する。 上位のヨーレート制御器より与えられた各輪の駆動力 指令を下位の車輪速度制御に基づくスリップ率制御器 [6] によって実現する。この方法の特徴は，路面のドラ 第 8 回制御部門大会（2008 年 3 月 5 日∼7 日・京都） Fig. 1: Experiment vehicle. Fig. 2: In-Wheel Motor. イビングスティフネスのリアルタイム同定を組み合わ せることによって，車輪がスリップしやすい低 µ 路に おいても，より正確な駆動力制御が可能となる点であ る。このため，正確なヨーモーメント生成と車輪空転の 抑制から，従来のヨーレート制御に比べてより安定し た姿勢制御を実現する。提案する制御系をシミュレー ションと小型の電気自動車を用いた走行実験により検 証する。 2 実験車両 実験機は市販の小型電気自動車 (トヨタ車体製 COMS ロング) をベースに改造している。この車両の特徴は後 輪二輪にそれぞれインホイールモータを搭載しており， 高い車両制御性を備えている点である。さらに Myway 技研株式会社と協力して制作したインバータシステム SY0004/08/0000-06414 c 2008 SICE Table 1: Specification of Experiment Vehicle. Drive Train 2 PM Motors Max. Power 2kW×2 Max. Torque 100Nm×2 Gear Ratio 6.267 Battery Lead Acid Total Voltage 72V (12V×6) Controller DSP TMS320VC33-150MHz Body Size Weight 480kg (with one person) Wheel Inertia 1.0Nms2 Wheel Radius 0.22m Ǭ ǩH ;H 3.1 ǩT NT ;T Fig. 3: Equivalent two-wheel vehicle model. δ γ∗ Kf τf.s+1 In.s + - + + Nin + P Nz + - N^ t ヨーモーメントオブザーバを用いた走行 安定化制御 [2] Vehicle Plant γ In.s ωc s+ωc 車両運動方程式 車両モデルとして Fig. 3 に示すように四輪を車軸に 関して対象と仮定した二輪車モデルを想定し，車体の ローリングを無視し，一定速度で走行していると仮定 する。この車両モデルの水平面内における運動方程式 は次式で表される。 M V (β˙ + γ) = 2Yf + 2Yr (1) I γ˙ = 2Yf lf − 2Yr lr (2) ただし，M は車体質量，V は車体速度，β は車体横滑 り角，γ はヨーレート，I はヨー軸周りの慣性モーメ ント，lf ，lr はそれぞれ車体重心から前輪軸，後輪軸ま での距離である。Yf ，Yr は前輪，後輪のコーナリング フォースであり，タイヤの特性に依存する非線形変数 である。 3.2 Ǫ ǫ を用いてモータを制御している。サンプリング周波数 は 10kHz である。Fig. 1 に実験車両の外観を，Fig. 2 に搭載されているインホイールモータを，Table. 1 に 車両の主要諸元を示す。 3 8 NH Fig. 4: Block diagram of vehicle stability control based on yaw-moment observer. 構成すれば外乱モーメントは抑圧され， γ= I dγ = Nz − Ndt dt (3) と簡単化される。式 (3) において，ヨーレート γ を検 出し，Fig. 4 に示されるヨーモーメントオブザーバを (4) とノミナル化される。本システムでは生成した駆動力 モーメント Nz を以下の方法で左右の車輪に駆動力分 配する。 Nz ∗ Fac = lp (FR∗ − FL∗ ) (5) = (6) FR∗ + FL∗ ∗ ただし，Fac はアクセル指令，FR∗ ，FL∗ は右輪，左輪の 駆動力指令である。式 (5)，(6) を解くことで，左右そ れぞれの駆動力を決定することができる。さらに， ヨーレート制御法 上記の β や V を測定するには車体の実速度ベクトル が必要となる。また，コーナリングフォースについて は非線形変数であるので測定や推定は困難である。文 献 [2] ではこれらのパラメータによる影響を外乱モーメ ントと捉え外乱オブザーバによる一括補償を提案して いる。 式 (2) の右辺で表されるタイヤに発生するコーナリ ングフォースによるモーメントを Nt ，式 (2) で考慮さ れていない路面状況の変化，横風などによる外乱モー メントを Nd ，実験機の後輪に装着されている二つのイ ンホイールモータの駆動力差により生成できるヨー軸 周りのモーメントを Nz とし，Ndt = Nd − Nt とおく と，式 (2) は， 1 Nin In s TR∗ = rFR∗ (7) = (8) TL∗ rFL∗ として，各車輪駆動モータに対するトルク指令を決定 する。ただし，r は車輪半径，lp は各モータから車軸 中心までの距離である。 実験機のモータには 100Nm の出力トルクリミッタが ∗ が共に大きいと，式 (5)，(6) の あるため，Nz と Fac 条件を満たせなくなる。その際に本稿ではアクセル指 ∗ を引き下げることで駆動力モーメント Nz の生 令 Fac 成を優先して確保している。 4 車輪速度制御に基づくスリップ率制御 [6] 4.1 並進方向運動方程式 モータの時定数が非常に小さいと仮定する時，車両 の運動方程式は以下の式で表現できる。 Jω ω˙ M V˙ Vω = T − rFd (9) = Fd − Fdr (10) = rω (11) ω GS + 2+ EQPVTQNNGT − 6 λ 1500 8ω 8 8GJKENG 2NCPV dry driving force[N] λ GS Fig. 5: Block diagram of slip-ratio control based on wheel rotation speed control. 4 (F 6CE 6&. 0\ . (F &U &U λ4 #FCRVKXG 2+EQPVTQNNGT λ. 4 6 1000 wet 500 ice 0 0 0.2 0.4 λ 4 ECNEWNCVQT 8ω 4 6 λ . ECNEWNCVQT 8ω . Fig. 6: Block diagram of yaw-rate control based on slip-ratio control. 各変数は，モータの回転角速度 ω ，車体速度 V ，車輪 速度 Vω ，モータトルク T ，駆動力 Fd ，走行抵抗 Fdr とする。各定数は，車体質量 M ，タイヤ半径 r，車輪 回転慣性モーメント Jω とする。またスリップ率は式 (12) で表されるが，本稿では駆動時のみ考えるものと し max(Vω , V ) = Vω とする。 但し，J(λ) は車輪と車体を合わせた回転部分の実際の 慣性モーメントである。 文献 [6] ではこの式に基づいた以下のスリップ率推定 器を提案している。 4.2 T rFˆdr ω˙ − 2 + 2 (14) ω r Mω r Mω 車輪速度制御 センサから得た車輪速度と推定して得られたスリッ プ率から車体速度を得ることが出来る。そして車体速 度と目標スリップ率から目標車輪速度を計算すること が出来る。 ω∗ = ˆ Vˆ rω(1 − λ) = ∗ ∗ r(1 − λ ) r(1 − λ ) T (15) これより，一般的に使用されているモータの電流制御 ループの外側に速度ループを組んだ回転速度制御を用 いることで，スリップ率制御を実現する。 @ (F Fig. 8: Block diagram of driving force observer. 制御器は PI 制御器を用いる。一般的に速度ループは 車輪の慣性モーメントを用いて以下の式に基づいて極 配置法により制御ゲインを決定する。 (12) T − rFdr + r2 M ω λ˙ T − rFdr + r2 M ω λ˙ := (13) Jω + r2 M (1 − λ) J(λ) ,ωU .2( ω= この式 (12) を両辺時間微分し，式 (9) から (11) を代 入し，Vω , V, Fd を消去すると以下の式を得ることが出 来る。 ω 8GJKENG 2NCPV ˆ + 1 + Jω ˆ˙ = − ω˙ λ λ ω r2 M 1 Fig. 7: Driving force curve and driving stiffness. 8 ω˙ = 0.8 slip ratio[-] . 6 #FCRVKXG 2+EQPVTQNNGT Vω − V λ= max(Vω , V ) 0.6 1 T Jω s (16) しかし式 (13) より車体全体の慣性モーメントはスリッ プ率に応じて変化するため車輪の慣性モーメントのみ で制御器を設計すると，スリップ率の変化に応じて相 対的に制御器の極が変化する。そこで式 (13) の J(λ) = Jω + r2 M (1 − λ) に対して極が一定になるように制御 ゲインをオンラインチューニングする。Fig. 5 に車輪 速度制御に基づくスリップ率制御ブロック図を示す。 ω= 5 1 T (Jω + r2 M (1 − λ))s (17) スリップ率制御に基づくヨーレート制御 5.1 スリップ率制御の指令値生成 ヨーレート制御では左右の駆動力差によりヨーモー メントを発生し車体の姿勢を制御するが，低 µ 路などに おいて駆動輪が空転するとコントローラの求めるヨー モーメントが得られなかったり，駆動輪のコーナリン グフォースが著しく減少し後輪駆動車の場合，スピン する事などが考えられる。 乾燥路において車体速度がほぼ一定のヨーレート制 御では，車輪角速度の変動 ω˙ が十分に小さいため式 (9) より式 (7),(8) も成り立つ。しかし低 µ 路において加減 速を伴う旋回を行うと車輪は容易に空転し，ω˙ が大き な値を持つためこの関係は成立しない。このときヨー レート制御のために正確に駆動力を制御するにはスリッ プ率制御が有効である。 そこで Fig. 6 に示すように，上位のヨーレート制御 系から与えられた駆動力指令を，路面のドライビング 100 0.1 7000 0.08 60 0.06 40 20 0.5 1 1.5 0.04 0.02 conventional method proposed method reference 0 0 driving stiffness 6500 driving stiffness[-] 80 slip-ratio[-] torque[Nm] slip-ratio reference 2 0 0 6000 5500 5000 4500 0.5 time[sec] 1 1.5 2 4000 0 0.5 1 1.5 2 time[sec] time[sec] Fig. 9: Driving force control based Fig. 10: Slip-ratio response in simu- Fig. 11: Identification of driving stiffon slip-ratio control. ness in simulation. lation. スティフネスを用いてスリップ率の指令値に換算し，下 位のスリップ率制御系に与えることでより正確な駆動 力配分を実現する。 λ∗ = Fd∗ Ds (18) 逐次型最小二乗法は観測信号 ξ が小さくなるとゼロ 割に近い演算が発生し推定精度は大きく低下する。そ のため本稿では条件付更新 [9] を用いて，ξ が ±0.005 より大きいときのみ同定値を更新している。 5.3 提案法による駆動力制御 ただし，λ∗ はスリップ率指令値，Fd∗ は駆動力指令 値，Ds はドライビングスティフネスである。この Ds 値としてリアルタイムに同定した値を制御に用いるこ とで，路面変化に対応したより正確な駆動力制御が可 能となる。 Fig. 7 はスリップ率に対する駆動力の発生分布であ る。駆動力はスリップ率が小さい領域ではほぼ線形に 増加する。この領域の線形増加の傾きをドライビング スティフネスといい，この比例係数を用いることで駆 動力指令をスリップ率指令に変換する。 一般的に駆動力はスリップ率 0.2 辺りのピークを越 えると緩やかに減衰する特性を持つ。この領域ではタ イヤは空転しているためモータトルクに対して十分な 駆動力が得られない事に加え，コーナリングフォース も著しく低下するため後輪駆動車の場合スピンを誘発 する。そのため駆動力はこのピークより手前の領域で 制御することが望ましい。 上記のスリップ率制御とドライビングスティフネス 同定を組み合わせた駆動力制御のシミュレーションを 行う。指令値を 80Nm から 50Nm にステップ状に変動 させたときに，実際に実現できる駆動力を比較する。駆 動力は駆動力オブザーバによって観測する。 Fig. 9 にシミュレーションによるドライビングスティ フネス同定値を用いた駆動力制御の結果を示す。但し 図中ではトルク次元で表示している。式 (7),(8) に基づ いたオープンループ制御では ω˙ = 0 のため厳密には指 令値どおりの駆動トルクを実現できないが，ドライビ ングスティフネスをリアルタイムに同定しスリップ率 制御を行うことで指令値どおりの駆動力を実現するこ とが出来る。その時のスリップ率応答を Fig. 10 に，ド ライビングスティフネス同定値を Fig. 11 に示す。車輪 速度制御の極は-50rad/s とした。スリップ率が小さく なるにつれドライビングスティフネスは大きな値が同 定されている。 5.2 ドライビングスティフネス同定 駆動輪一輪あたりの発生する駆動力は以下の式によっ て記述できる。 6 Fd = Ds λ (19) := θ · ξ (20) ただし，Fd は駆動力である。本稿では Fd は駆動力オブ ザーバ [7] を用いて推定されている。駆動力オブザーバ のブロック図を Fig. 8 に示す。Fd = M V˙ として加速度 センサを用いる方法は車体のピッチングやローリング の影響，走行抵抗が無視できないため不向きである。ま た，λ は推定もしくは測定が可能であるため，ξ := λ と して上記のように書き換え未知変数を Ds のみとする。 本稿ではドライビングスティフネスの同定に忘却係 数 λ∗ を用いた逐次型最小二乗法を用い，次式に基づき 推定する [8]。 Γ(k) = Γ(k − 1)ξ 2 (k)Γ(k − 1) 1 Γ(k − 1) − (21) λ∗ λ∗ + ξ(k)Γ(k − 1)ξ ˆ θ(k) = θ(k ˆ− 1) − Γ(k − 1)ξ(k) λ∗ + ξ(k)Γ(k − 1)ξ(k) ˆ − 1) − y(k) (22) · ξ(k)θ(k ヨーレート制御シミュレーション 乾燥路から 0.5sec に低 µ 路に進入して加速旋回する シミュレーションを行う。車両モデルは実験車両に基 づき，車体質量 480kg，ヨー慣性モーメント 351Nms2 ， ホイールベース 1.71m，トレッド 0.815m とした。時速 16km/h で走行中にハンドルを 180deg(実舵角 9.29deg に相当) 切って加速旋回する。駆動トルクは 140Nm(片 輪 70Nm) とし，ヨーレート指令値は乾燥路を同程度の 車体速度，操舵角で走行した時のヨーレート応答を基 に生成した。ヨーレート制御器はノミナルプラント In1 s に対して閉ループ極を −5rad/s に配置した比例制御と し，YMO のカットオフ周波数は 10rad/s とした。ス リップ率制御の極は −50rad/s とし，指令値には ±0.05 のリミッタを設定した。逐次型最小二乗法の忘却係数 は 0.99 とした。なお本稿で検証するにあたって作成し た車両シミュレータには坂井氏 [10] の 2 次元タイヤモ デル/発生力計算プログラム 2D-TireModel ver 1.0 を 使用している。 Fig. 12 にヨーレート応答を示す。非制御時とヨー レート制御時では後輪が滑りスピン傾向へ移行してい るが，提案法では乾燥路と同等の旋回性能を得られて いる。 7 1 without control yaw control yaw + slip control on dry road 6 lateral acceleration[m/s/s] yaw-rate[rad/s] 0.8 0.6 0.4 without control yaw control yaw + slip control reference 0.2 0 0 0.5 1 1.5 5 4 3 2 1 0 0 2 0.5 time[sec] Fig. 12: Yaw-rate response in simulation. 5 14 yaw-rate + slip-ratio control driving stiffness[-] y[m] yaw-rate control low mu terrain 3 2 1 5 10 15 20 x[m] Fig. 14: Trajectory of the vehicle in simulation. Fig. 13 に横加速度応答を示す。乾燥路では速度の増 加に伴って，横加速度も一定増加となっている。非制 御時とヨーレート制御のみの時では低 µ 路に進入する と車輪が空転し横滑りするため横加速度は急激に増加 するが，提案法では乾燥路ほどの旋回は実現できない ものの，車輪の空転が抑え横滑りも抑えることが出来 るため，横加速度の急激な増加も抑えられている。 車両の走行軌跡を Fig. 14 に示す。非制御時はヨー モーメントによる姿勢制御が行われないため車体は外 側に流れ，特に後輪駆動であるため後輪のコーナリン グフォースが著しく減少しスピンし始めている。ヨー レート制御のみ行ったときでは車体のスピンを抑える 方向にヨーモーメントを発生するために左右のモータ にそれぞれの方向に過大なトルク指令が与えられるが， 大空転した両輪が横滑りしスピンに至る。スリップ率 制御を加えたヨーレート制御では，DYC の駆動力配分 を活かしつつ車輪の空転を抑えることでスピンやドリ フトを抑えた旋回を実現している。 Fig. 15 にドライビングスティフネスの同定値を示す。 乾燥路では 45000 程度，低 µ 路では 5000 以下程度と， 路面の滑りやすさに応じてドライビングスティフネス が同定されている。 7 4 4 2 0 0 2 driving stiffness identification 6 4 x 10 without control 10 8 1.5 Fig. 13: Lateral acceleration response in simulation. 16 12 1 time[sec] 実験結果 実験環境はアクリル製の板の上に散水することでつ くった低 µ 路面上で行う。0.7sec より低 µ 路に進入し 1.5sec で乾燥路に復帰している。実験条件はシミュレー ションと同様とする。トルク指令値もシミュレーショ ンと同様としているが，実験では走行抵抗の影響によ り結果的に等速旋回となっている。なお実験では簡単 0 0 0.5 1 1.5 2 time[sec] Fig. 15: Identification of driving stiffness in simulation. のため車体速度は推定値ではなく，非駆動輪である前 輪の車輪速度をエンコーダで取得した値を使用して制 御を行っている。ヨーレート指令値は乾燥路を同程度 の車体速度，操舵角で走行した時のヨーレート応答を 基に生成した。 Fig. 16 にヨーレート応答を示す。非制御時とヨー レート制御時ではスピンしている。ヨーレート制御で は式 (5) より左右のトルクの和は非制御時と同様とな るが，ヨーモーメント生成の指令値によっては片輪あ たりのトルク指令は非制御時よりも大きくなる事があ るため，結果的に車輪の空転は大きくなり車両はスピ ンする。提案法では若干ヨーレートが増大しているが 非制御時に比べてその増加量は小さく抑えられている。 Fig. 17 に横加速度応答を示す。スピンする非制御時 とヨーレート制御時では横滑りが発生し，横加速度が 増加している。提案法では制御によるモータの振動に 若干センサが影響されているが，横滑りを最も小さく 抑えることが出来ている。 Fig. 18 に旋回内輪のスリップ率応答を示す。非制御 時とヨーレート制御時では同様にスピンしている。ヨー レート制御では上述の理由により非制御時よりもスリッ プ率が増大している。提案法ではトルク指令と同定し たドライビングスティフネスより生成したスリップ率 指令値どおりに制御され，発散を抑えている。 Fig. 17 にドライビングスティフネスを示す。乾燥路 でおよそ 18000 程度，低 µ 路で 5000 程度と同定され ている。それぞれの路面への進入に伴って速やかに同 定されている。 1 6 5 lateral acceleration[m/s/s] 0.8 yaw-rate[rad/s] without control yaw control yaw + slip control 0.6 0.4 without control yaw control yaw + slip control reference 0.2 0 0 0.5 1 1.5 4 3 2 1 0 0 2 0.5 time[sec] Fig. 16: Yaw-rate response in experiment. 3 without control yaw control yaw + slip control reference driving stiffness[-] slip-ratio[-] 0.1 0 4 2 High mu terrain 1.5 1 0.5 0.5 1 1.5 2 time[sec] Fig. 18: Slip-ratio response in experiment. 8 x 10 2.5 0.2 0 2 driving stiffness identification 0.3 -0.1 1.5 Fig. 17: Lateral acceleration response in experiment. 0.5 0.4 1 time[sec] 結論 本稿ではドライビングスティフネス同定値を用いた スリップ率制御に基づくヨーレート制御法を提案し，そ の有効性をシミュレーションと実験によって示した。 従来の DYC は低 µ 路では車輪が空転し，車両の操 縦安定性能を低下させる可能性があった。しかし提案 法では，DYC 制御器の各輪の駆動力指令を路面の情報 に基づいてスリップ率指令に変換し，スリップ率制御 で実現することにより，車輪の空転によるスピンを抑 制するとともにヨーレート制御性能の向上を実現した。 今後は現在非駆動輪の車輪速度より測定している車 体速度を推定によって取得し，四輪駆動車にも適用で きる制御系の研究を進める予定である。 謝辞 最後に本研究の一部は NEDO 産業技術研究助成 (プ ロジェクト ID:05A48701d) によって行われたことを付 記する。 参考文献 [1] Y. Hori: “Future Vehicle Driven by Electricity and Control-Research on Four-Wheel-Motored: UOT Electric March II”, IEEE Trans. IE, Vol. 51, No. 5, 2004. [2] 齋藤, 藤本, 野口, “スリップ及びヨーモーメントオブザー バを用いた電気自動車の操縦安定化制御”, 電気学会産 業計測制御研究会, IIC-03-52, pp. 41–46, 2003. [3] 裙坂, 藤本, 野口, “コーナリングスティフネス推定値を用 いた電気自動車の走行安定化制御”, 電気学会論文誌 D, IEEJ Trans. IA, Vol. 126, No. 7, pp. 996–1002, 2006. [4] 原, 堀, 坂井, “車輪の駆動力飽和にもとづく電気自動 車の 2 次元運動制御”, 電気学会産業計測制御研究会, IIC-01-20, pp. 115–120, 2001. 0 0 Low mu terrain 0.5 1 1.5 2 time[sec] Fig. 19: Identification of driving stiffness in experiment. [5] 小野, 服部, 門崎, “4 輪アクティブ操舵による安全性・快 適性向上”, 自動車技術会学術講演会前刷集, No. 11-07, pp. 1–6, 2007. [6] 藤井, 藤本, “車体速検出不要な電気自動車の車輪速制御 に基づくスリップ率制御”, 電気学会自動車研究会, VT07-05, pp. 27–32, 2007. [7] H. Sado, S. Sakai, Y. Hori: “Road Condition Estimation for Traction Control in Electic Vehicle”, In Proc. 1999 IEEE International Symposium on Industrial Electronics,pp.973-978, 1999. [8] 堀, 大西, 応用制御工学, 丸善株式会社, 1998. [9] H. Fujimoto, B. Yao, “Multirate Adaptive Robust Control for Discrete-Time Non-Minimum Phase System and Application to Linear Motors”, IEEE/ASME Transactions on Mechatronics, Vol. 10, No. 4, pp. 371–377, 2005. [10] Shin-ichiro Sakai, “2D-TireModel ver 1.0”, http://sakai.nnl.isas.ac.jp/index_j.html , 2000.
© Copyright 2024 ExpyDoc
