(57)【特許請求の範囲】
【請求項１】
むだ時間を含む動特性が所定条件下で変化する特性を有し、参照パラメータが所定範囲内で変化するのに応じて、当該むだ時間が最大値および最小値を含むＭ（Ｍは値２以上の整数）個の整数値の間で連続的に変化するようにモデリングされる制御対象において、当該制御対象の制御量を制御入力によって制御する制御装置であって、
当該制御量の目標となる目標制御量を設定する目標制御量設定手段と、
前記参照パラメータを検出する参照パラメータ検出手段と、
前記制御量と前記制御入力の関係を定義した制御対象モデルを用い、前記Ｍ個のむだ時間がそれぞれ経過したタイミングでの前記制御量として、Ｍ個の予測値を算出する予測値算出手段と、
前記検出された参照パラメータに基づき、当該参照パラメータに対応するＭ個の重み関数値をそれぞれ算出する重み関数値算出手段と、
当該算出されたＭ個の重み関数値を前記算出されたＭ個の予測値にそれぞれ乗算することにより、Ｍ個の第１乗算値を算出するとともに、当該Ｍ個の第１乗算値の総和を前記制御量の予測値である予測制御量として設定する予測制御量設定手段と、
当該予測制御量が前記目標制御量になるように、前記制御入力を算出する制御入力算出手段と、
を備え、
前記Ｍ個の重み関数値は、前記参照パラメータの前記所定範囲におけるＭ個の領域にそれぞれ対応し、対応する領域において値０以外の値に設定され、かつ当該対応する領域以外において値０に設定されており、
当該Ｍ個の領域は隣り合う領域が互いに重なり合っており、
前記Ｍ個の重み関数値は、当該重なり合う領域における前記参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が所定値になるように設定されていることを特徴とする制御装置。
【請求項２】
前記Ｍ個の重み関数値を、前記Ｍ個のむだ時間前のタイミングでのＭ個の前記制御入力にそれぞれ乗算することにより、Ｍ個の第２乗算値を算出するとともに、当該Ｍ個の第２乗算値の総和を修正後制御入力として設定する修正後制御入力設定手段と、
前記制御量と当該修正後制御入力との関係を定義した修正後モデルを用いて導出した所定の同定アルゴリズムにより、当該修正後モデルのモデルパラメータをオンボード同定する同定手段と、
をさらに備え、
前記予測値算出手段は、当該同定されたモデルパラメータを前記制御対象モデルのモデルパラメータとして用いることを特徴とする請求項１に記載の制御装置。
【請求項３】
前記制御入力算出手段は、所定の周波数特性が得られるように設定された感度関数、相補感度関数および伝達関数のいずれかに基づいて導出された制御アルゴリズムを用いて、前記制御入力を算出することを特徴とする請求項１または２に記載の制御装置。
【請求項４】
むだ時間を含む動特性が所定条件下で変化する特性を有し、参照パラメータが所定範囲内で変化するのに応じて、当該むだ時間が最大値および最小値を含むＭ（Ｍは値２以上の整数）個の整数値の間で連続的に変化するようにモデリングされる制御対象において、当該制御対象の制御量を制御入力によって制御する制御装置であって、
前記参照パラメータを検出する参照パラメータ検出手段と、
前記検出された参照パラメータに基づき、当該参照パラメータに対応するＭ個の重み関数値をそれぞれ算出する重み関数値算出手段と、
当該算出されたＭ個の重み関数値を、前記Ｍ個のむだ時間前のタイミングでのＭ個の前記制御入力にそれぞれ乗算することにより、Ｍ個の乗算値を算出するとともに、当該Ｍ個の乗算値の総和を修正後制御入力として設定する修正後制御入力設定手段と、
前記制御量と当該修正後制御入力との関係を定義した修正後モデルを用いて導出した所定の同定アルゴリズムにより、当該修正後モデルのモデルパラメータをオンボード同定する同定手段と、
所定の制御アルゴリズムおよび制御対象モデルを用いて、前記制御入力を算出するとともに、当該制御対象モデルのモデルパラメータとして前記同定されたモデルパラメータを用いる制御入力算出手段と、
を備え、
前記Ｍ個の重み関数値は、前記参照パラメータの前記所定範囲におけるＭ個の領域にそれぞれ対応し、対応する領域において値０以外の値に設定され、かつ当該対応する領域以外において値０に設定されており、
当該Ｍ個の領域は隣り合う領域が互いに重なり合っており、
前記Ｍ個の重み関数値は、当該重なり合う領域における前記参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が所定値になるように設定されていることを特徴とする制御装置。
【請求項５】
前記所定の制御アルゴリズムは、所定の周波数特性が得られるように設定された感度関数、相補感度関数および伝達関数のいずれかに基づいて導出されたアルゴリズムであることを特徴とする請求項４に記載の制御装置。
【請求項６】
むだ時間を含む動特性が所定条件下で変化する特性を有し、参照パラメータが所定範囲内で変化するのに応じて、当該むだ時間が最大値および最小値を含むＭ（Ｍは値２以上の整数）個の整数値の間で連続的に変化するようにモデリングされる制御対象において、当該制御対象の制御量を制御入力によって制御する制御装置であって、
当該制御量の目標となる目標制御量を設定する目標制御量設定手段と、
前記参照パラメータを検出する参照パラメータ検出手段と、
当該検出された参照パラメータに基づき、当該参照パラメータに対応するＭ個の重み関数値をそれぞれ算出する重み関数値算出手段と、
当該算出されたＭ個の重み関数値を、前記Ｍ個のむだ時間前のタイミングでのＭ個の前記制御入力にそれぞれ乗算することにより、Ｍ個の乗算値を算出するとともに、当該Ｍ個の乗算値の総和を修正後制御入力として設定する修正後制御入力設定手段と、
当該修正後制御入力および前記制御量を用いて、外乱推定値を算出する外乱推定値算出手段と、
当該算出された外乱推定値を用い、前記制御量が前記目標制御量になるように、前記制御入力を算出する制御入力算出手段と、
を備え、
前記Ｍ個の重み関数値は、前記参照パラメータの前記所定範囲におけるＭ個の領域にそれぞれ対応し、対応する領域において値０以外の値に設定され、かつ当該対応する領域以外において値０に設定されており、
当該Ｍ個の領域は隣り合う領域が互いに重なり合っており、
前記Ｍ個の重み関数値は、当該重なり合う領域における前記参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が所定値になるように設定されていることを特徴とする制御装置。
【請求項７】
前記外乱推定値算出手段は、前記制御量の推定値である推定制御量を、当該推定制御量と前記修正後制御入力と外乱推定値と前記制御量との関係を定義したモデルを用いて算出するとともに、当該推定制御量と前記制御量との偏差が最小となるように、前記外乱推定値を算出することを特徴とする請求項６に記載の制御装置。
【請求項８】
前記制御量は、内燃機関の混合気の空燃比を表す値であり、前記制御入力は、当該内燃機関の燃料量を補正する補正係数であることを特徴とする請求項１ないし７のいずれかに記載の制御装置。
【請求項９】
前記制御量は、自動変速機の伝達トルク調整機構におけ出力回転数を表す値であり、前記制御入力は当該伝達トルク調整機構のアクチュエータへの入力であることを特徴とする請求項１ないし７のいずれかに記載の制御装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、むだ時間が変化する特性を備えた制御対象の制御量を制御入力によって制御する制御装置に関する。
【背景技術】
【０００２】
従来、内燃機関の混合気の空燃比を制御する制御装置として、特許文献１に記載されたものを本出願人は既に提案している。制御装置は、ＬＡＦセンサ、酸素濃度センサ、状態予測器、オンボード同定器、スライディングモードコントローラおよび目標空燃比算出部などを備えている。これらのＬＡＦセンサおよび酸素濃度センサはいずれも、内燃機関の排気通路内の排気ガス中の酸素濃度を表す値すなわち空燃比を検出するためのものであり、排気通路の集合部よりも下流側に設けられている。また、内燃機関は、ガソリンを燃料とするガソリンエンジンタイプのものであり、排気通路の集合部よりも下流側設けられた第１触媒装置と、その下流側に設けられた第２触媒装置を備えている。ＬＡＦセンサは、第１触媒装置の上流側に配置され、酸素濃度センサは、第１触媒装置と第２触媒装置の間に配置されている。
【０００３】
この制御装置では、制御対象モデルとして、ＬＡＦセンサの検出空燃比ＫＡＣＴと空燃比基準値ＦＬＡＦＢＡＳＥとの偏差（以下「空燃比偏差」という）ｋａｃｔを入力とし、酸素濃度センサの出力ＶＯＵＴと所定の目標値ＶＯＵＴ＿ＴＡＲＧＥＴとの偏差（以下「出力偏差」という）ＶＯ２を出力とする離散時間系モデルと、ＬＡＦセンサで検出された排ガスの空燃比が酸素濃度センサに検出されるまでのむだ時間ｄ１と、目標空燃比ＫＣＭＤがＬＡＦセンサの検出結果に反映されるまでのむだ時間ｄ２とを用いて、所定の制御アルゴリズムにより、制御入力としての目標空燃比ＫＣＭＤが算出される。また、これらのむだ時間ｄ１，ｄ２はいずれも一定値に設定されている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２０００−２３４５５０号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
特許文献１のように構成された内燃機関の場合、内燃機関の運転状態の変化、経年変化および個体間のばらつきなどに起因して、上記２つのむだ時間ｄ１，ｄ２の実際の値が変化する。その場合、特許文献１の制御装置によれば、２つのむだ時間ｄ１，ｄ２として一定の設定値を用いている関係上、制御精度の低下を招いてしまう。このような問題は、特許文献１のように、空燃比を制御する場合に限らず、むだ時間や応答遅れが変化する特性を備えた制御対象を制御する際にも生じてしまう。例えば、自動変速機のクラッチの接続・遮断状態を制御する場合にも生じてしまう。
【０００６】
本発明は、上記課題を解決するためになされたもので、むだ時間や応答遅れが変化する特性を備えた制御対象を制御する場合において、制御精度を向上させることができる制御装置を提供することを目的とする。
【課題を解決するための手段】
【０００７】
上記目的を達成するために、請求項１に係る発明は、むだ時間ｄ，ｄ”を含む動特性が所定条件下で変化する特性を有し、参照パラメータ（排ガスボリュームＶｅｘ、油温Ｔｏｉｌ）が所定範囲（０≦Ｖｅｘ≦ＶｅｘＭＡＸ、Ｔｏｉｌ≦ＴｏｉｌＭＡＸ）内で変化するのに応じて、むだ時間ｄ，ｄ”が最大値（値３）および最小値（値０）を含むＭ（Ｍは値２以上の整数）個の整数値の間で連続的に変化するようにモデリングされる制御対象において、制御対象の制御量（検出当量比ＫＡＣＴ，主軸回転数ＮＭ）を制御入力（空燃比補正係数ＫＡＦ、制御入力Ｕａｃｔ）によって制御する制御装置１，１Ａ〜１Ｅであって、制御量の目標となる目標制御量（目標当量比ＫＣＭＤ、目標主軸回転数ＮＭ＿ｃｍｄ）を設定する目標制御量設定手段（ＥＣＵ２、目標当量比算出部３０，２３０、目標主軸回転数算出部５１０）と、参照パラメータを検出する参照パラメータ検出手段（ＥＣＵ２、クランク角センサ２０、吸気圧センサ２２、油温センサ２６、ステップ１０）と、制御量と制御入力の関係を定義した制御対象モデル（式（２），（６８），（１０５））を用い、Ｍ個のむだ時間ｄ，ｄ”がそれぞれ経過したタイミングでの制御量として、Ｍ個の予測値（第０〜第３予測値ＰＲＥ＿ＫＡＣＴ＿０〜３、第０〜第３予測値ＰＲＥ＿ＮＭ＿０〜３）を算出する予測値算出手段（ＥＣＵ２、可変むだ時間状態予測器４０，２４０，５２０）と、検出された参照パラメータに基づき、参照パラメータに対応するＭ個（４個）の重み関数値Ｗｄｉ，Ｗｄｉ”（ｉ＝１〜４）をそれぞれ算出する重み関数値算出手段（ＥＣＵ２、可変むだ時間状態予測器４０，２４０，５２０）と、算出されたＭ個の重み関数値を算出されたＭ個の予測値にそれぞれ乗算することにより、Ｍ個の第１乗算値（Ｗｄｉ・ＰＲＥ＿ＫＡＣＴ＿４−ｉ、Ｗｄｉ”・ＰＲＥ＿ＮＭ＿４−ｉ）を算出するとともに、Ｍ個の第１乗算値の総和を制御量の予測値である予測制御量（予測当量比ＰＲＥ＿ＫＡＣＴ、予測主軸回転数ＰＲＥ＿ＮＭ）として設定する予測制御量設定手段（ＥＣＵ２、可変むだ時間状態予測器４０，２４０，５２０）と、予測制御量が目標制御量になるように、制御入力を算出する制御入力算出手段（ＥＣＵ２、周波数整形コントローラ１３０，３３０，５４０、２自由度応答指定型コントローラ３５０，３８０，７５０）と、を備え、Ｍ個（４個）の重み関数値Ｗｄｉ，Ｗｄｉ”は、参照パラメータの所定範囲におけるＭ個の領域にそれぞれ対応し、対応する領域において値０以外の値に設定され、かつ対応する領域以外において値０に設定されており、Ｍ個の領域は隣り合う領域が互いに重なり合っており、Ｍ個の重み関数値は、重なり合う領域における参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が所定値（値１）になるように設定されていることを特徴とする。
【０００８】
この制御装置によれば、制御量と制御入力の関係を定義した制御対象モデルを用い、Ｍ個のむだ時間がそれぞれ経過したタイミングでの制御量として、Ｍ個の予測値が算出され、検出された参照パラメータに基づき、参照パラメータに対応するＭ個の重み関数値がそれぞれ算出されるとともに、算出されたＭ個の重み関数値を算出されたＭ個の予測値にそれぞれ乗算することにより、Ｍ個の第１乗算値が算出される。さらに、Ｍ個の第１乗算値の総和が制御量の予測値である予測制御量として設定され、予測制御量が目標制御量になるように、制御入力が算出される。この場合、Ｍ個の重み関数値は、参照パラメータの所定範囲におけるＭ個の領域にそれぞれ対応し、対応する領域において値０以外の値に設定され、かつ対応する領域以外において値０に設定されており、Ｍ個の領域は隣り合う領域が互いに重なり合っており、Ｍ個の重み関数値は、重なり合う領域における参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が所定値になるように設定されている。
【０００９】
したがって、そのようなＭ個の重み関数値をＭ個の予測値にそれぞれ乗算したＭ個の第１乗算値は、Ｍ個の予測値を互いに連続するように重み付けされた値として算出されるとともに、そのような第１乗算値の総和が予測制御量として設定されるので、予測制御量を、Ｍ個の予測値を連続的に結合させた値として算出することができる。それにより、参照パラメータの変化に応じてむだ時間が変化した場合でも、予測制御量を、そのようなむだ時間の変化を適切に補償しながら精度よく算出することができる。特に、参照パラメータの急変に伴ってむだ時間が急変したときでも、予測制御量を、むだ時間の急変を適切に補償しながら、段差なく円滑に変化するように算出することができる。以上のように、予測制御量を精度よく算出することができる。そして、そのように算出された予測制御量が目標制御量になるように、制御入力が算出されるので、この制御入力によって、制御量を目標制御量に精度よく制御することができる。特に、制御入力の算出アルゴリズムとして、フィードバック制御アルゴリズムを用いた場合には、フィードバックゲインを高い状態に維持することができ、高精度かつ高応答性を確保しながら、制御量を目標制御量に追従させることができる。
【００１０】
請求項２に係る発明は、請求項１に記載の制御装置１，１Ａ，１Ｄにおいて、Ｍ個の重み関数値Ｗｄｉ，Ｗｄｉ”を、Ｍ個のむだ時間ｄ前のタイミングでのＭ個の制御入力（空燃比補正係数ＫＡＦ（ｋ−４＋ｉ），制御入力Ｕａｃｔ（ｋ−４＋ｉ））にそれぞれ乗算することにより、Ｍ個の第２乗算値（Ｗｄｉ・ＫＡＦ（ｋ−４＋ｉ），Ｗｄｉ”・Ｕａｃｔ（ｋ−４＋ｉ））を算出するとともに、Ｍ個の第２乗算値の総和を修正後制御入力ＫＡＦ＿ｍｏｄ，Ｕａｃｔ＿ｍｏｄとして設定する修正後制御入力設定手段（ＥＣＵ２、オンボード同定器６０，２６０，５３０）と、制御量と修正後制御入力との関係を定義した修正後モデル［式（３０）］を用いて導出した所定の同定アルゴリズム［式（１７）〜（２９），（７４）〜（８６），（１１２）〜（１２４）］により、修正後モデルのモデルパラメータ（同定値αｉｄ，αｉｄ”、モデルパラメータδ，α）をオンボード同定する同定手段（ＥＣＵ２、オンボード同定器６０，２６０，５３０）と、をさらに備え、予測値算出手段は、同定されたモデルパラメータを制御対象モデルのモデルパラメータとして用いることを特徴とする。
【００１１】
この制御装置によれば、Ｍ個の重み関数値を、Ｍ個のむだ時間前のタイミングでのＭ個の制御入力にそれぞれ乗算することにより、Ｍ個の第２乗算値が算出され、Ｍ個の第２乗算値の総和が修正後制御入力として設定される。この場合、Ｍ個の重み関数値は参照パラメータに対して前述したように設定されているので、参照パラメータの変化に応じてむだ時間が連続的に変化する場合でも、修正後制御入力を、そのようなむだ時間の変化を適切に補償しながら精度よく算出することができる。特に、参照パラメータの急変に伴ってむだ時間が急変したときでも、修正後制御入力を、むだ時間の急変を適切に補償しながら、段差なく円滑に変化するように算出することができる。さらに、制御量と修正後制御入力との関係を定義した修正後モデルを用いて導出した所定の同定アルゴリズムにより、修正後モデルのモデルパラメータがオンボード同定されるので、参照パラメータの変化に応じてむだ時間が変化したときでも、その影響を抑制しながら、制御入力モデルのモデルパラメータを精度よく同定することができる。そして、そのようなモデルパラメータを制御対象モデルのモデルパラメータとして用いられるので、制御性や、個体間のばらつきおよび経年変化の影響に対する制御のロバスト性を飛躍的に向上させることができる。
【００１２】
請求項３に係る発明は、請求項１または２に記載の制御装置１，１Ａ〜１Ｅにおいて、制御入力算出手段は、所定の周波数特性が得られるように設定された感度関数Ｓｄ、相補感度関数Ｔｄおよび伝達関数Ｃのいずれかに基づいて導出された制御アルゴリズム［（３４），（３５），（８７），（８８），（１２５），（１２６）］を用いて、制御入力を算出することを特徴とする。
【００１３】
この制御装置によれば、所定の周波数特性が得られるように設定された感度関数、相補感度関数および伝達関数のいずれかに基づいて導出された制御アルゴリズムを用いて、制御入力が算出されるので、前述したようにむだ時間の変化を適切に補償しながら、制御装置における外乱抑制特性やロバスト性を周波数軸上でダイレクトに指定（設定）することができる。それにより、外乱による制御量の変動を抑制したい周波数域において、外乱抑制能力やロバスト性を飛躍的に向上させることができる。
【００１４】
請求項４に係る発明は、むだ時間ｄ，ｄ”を含む動特性が所定条件下で変化する特性を有し、参照パラメータ（排ガスボリュームＶｅｘ、油温Ｔｏｉｌ）が所定範囲（０≦Ｖｅｘ≦ＶｅｘＭＡＸ、Ｔｏｉｌ≦ＴｏｉｌＭＡＸ）内で変化するのに応じて、むだ時間ｄ，ｄ”が最大値（値３）および最小値（値０）を含むＭ（Ｍは値２以上の整数）個の整数値の間で連続的に変化するようにモデリングされる制御対象において、制御対象の制御量（検出当量比ＫＡＣＴ，主軸回転数ＮＭ）を制御入力（空燃比補正係数ＫＡＦ、制御入力Ｕａｃｔ）によって制御する制御装置１，１Ａ，１Ｄであって、参照パラメータを検出する参照パラメータ検出手段（ＥＣＵ２、クランク角センサ２０、吸気圧センサ２２、油温センサ２６、ステップ１０）と、検出された参照パラメータに基づき、参照パラメータに対応するＭ個（４個）の重み関数値Ｗｄｉ，Ｗｄｉ”（ｉ＝１〜４）をそれぞれ算出する重み関数値算出手段（ＥＣＵ２、オンボード同定器６０，２６０，５３０）と、算出されたＭ個の重み関数値Ｗｄｉ，Ｗｄｉ”を、Ｍ個のむだ時間ｄ前のタイミングでのＭ個の制御入力（空燃比補正係数ＫＡＦ（ｋ−４＋ｉ），制御入力Ｕａｃｔ（ｋ−４＋ｉ））にそれぞれ乗算することにより、Ｍ個の乗算値（Ｗｄｉ・ＫＡＦ（ｋ−４＋ｉ），Ｗｄｉ”・Ｕａｃｔ（ｋ−４＋ｉ））を算出するとともに、Ｍ個の乗算値の総和を修正後制御入力ＫＡＦ＿ｍｏｄ，Ｕａｃｔ＿ｍｏｄとして設定する修正後制御入力設定手段（ＥＣＵ２、オンボード同定器６０，２６０，５３０）と、制御量と修正後制御入力との関係を定義した修正後モデル［式（３０）］を用いて導出した所定の同定アルゴリズム［式（１７）〜（２９），（７４）〜（８６），（１１２）〜（１２４）］により、修正後モデルのモデルパラメータ（同定値αｉｄ，αｉｄ”、モデルパラメータδ，α）をオンボード同定する同定手段（ＥＣＵ２、オンボード同定器６０，２６０，５３０）と、所定の制御アルゴリズムおよび制御対象モデルを用いて、制御入力を算出するとともに、制御対象モデルのモデルパラメータとして同定されたモデルパラメータを用いる制御入力算出手段（ＥＣＵ２、周波数整形コントローラ３３０）と、を備え、Ｍ個（４個）の重み関数値Ｗｄｉ，Ｗｄｉ”は、参照パラメータの所定範囲におけるＭ個の領域にそれぞれ対応し、対応する領域において値０以外の値に設定され、かつ対応する領域以外において値０に設定されており、Ｍ個の領域は隣り合う領域が互いに重なり合っており、Ｍ個の重み関数値は、重なり合う領域における参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が所定値（値１）になるように設定されていることを特徴とする。
【００１５】
この制御装置によれば、検出された参照パラメータに基づき、参照パラメータに対応するＭ個の重み関数値がそれぞれ算出され、Ｍ個の重み関数値を、Ｍ個のむだ時間前のタイミングでのＭ個の制御入力にそれぞれ乗算することにより、Ｍ個の乗算値が算出され、Ｍ個の乗算値の総和が修正後制御入力として設定される。この場合、Ｍ個の重み関数値は、参照パラメータの所定範囲におけるＭ個の領域にそれぞれ対応し、対応する領域において値０以外の値に設定され、かつ対応する領域以外において値０に設定されており、Ｍ個の領域は隣り合う領域が互いに重なり合っており、Ｍ個の重み関数値は、重なり合う領域における参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が所定値（値１）になるように設定されている。したがって、そのようなＭ個の重み関数値をＭ個のむだ時間前のタイミングでのＭ個の制御入力にそれぞれ乗算したＭ個の乗算値の総和が修正後制御入力として設定されるので、参照パラメータの変化に応じてむだ時間が連続的に変化する場合でも、修正後制御入力を、そのようなむだ時間の変化を適切に補償しながら精度よく算出することができる。特に、参照パラメータの急変に伴ってむだ時間が急変したときでも、修正後制御入力を、むだ時間の急変を適切に補償しながら、段差なく円滑に変化するように算出することができる。
【００１６】
また、制御量と修正後制御入力との関係を定義した修正後モデルを用いて導出した所定の同定アルゴリズムにより、修正後モデルのモデルパラメータがオンボード同定されるので、参照パラメータの変化に応じてむだ時間が変化したときでも、その影響を抑制しながら、制御入力モデルのモデルパラメータを精度よく同定することができる。さらに、所定の制御アルゴリズムおよび制御対象モデルを用いて、制御入力が算出されるとともに、上記のように同定されたモデルパラメータが制御対象モデルのモデルパラメータとして用いられるので、制御性や、個体間のばらつきおよび経年変化の影響に対する制御のロバスト性も飛躍的に向上させることができる。
【００１７】
請求項５に係る発明は、請求項４に記載の制御装置１，１Ａ，１Ｄにおいて、所定の制御アルゴリズムは、所定の周波数特性が得られるように設定された感度関数Ｓｄ、相補感度関数Ｔｄおよび伝達関数Ｃのいずれかに基づいて導出されたアルゴリズムであることを特徴とする。
【００１８】
この制御装置によれば、所定の周波数特性が得られるように設定された感度関数、相補感度関数および伝達関数のいずれかに基づいて導出された制御アルゴリズムを用いて制御入力が算出されるので、制御装置における外乱抑制特性やロバスト性を周波数軸上でダイレクトに指定（設定）することができる。それにより、外乱による制御量の変動を抑制したい周波数域において、外乱抑制能力やロバスト性を飛躍的に向上させることができる。
【００１９】
請求項６に係る発明は、むだ時間ｄ，ｄ”を含む動特性が所定条件下で変化する特性を有し、参照パラメータ（排ガスボリュームＶｅｘ、油温Ｔｏｉｌ）が所定範囲（０≦Ｖｅｘ≦ＶｅｘＭＡＸ、Ｔｏｉｌ≦ＴｏｉｌＭＡＸ）内で変化するのに応じて、むだ時間ｄ，ｄ”が最大値（値３）および最小値（値０）を含むＭ（Ｍは値２以上の整数）個の整数値の間で連続的に変化するようにモデリングされる制御対象において、制御対象の制御量（検出当量比ＫＡＣＴ，主軸回転数ＮＭ）を制御入力（空燃比補正係数ＫＡＦ、制御入力Ｕａｃｔ）によって制御する制御装置１Ｃ，１Ｅであって、制御量の目標となる目標制御量（目標当量比ＫＣＭＤ、目標主軸回転数ＮＭ＿ｃｍｄ）を設定する目標制御量設定手段（ＥＣＵ２、目標当量比算出部３０，２３０、目標主軸回転数算出部５１０）と、参照パラメータを検出する参照パラメータ検出手段（ＥＣＵ２、クランク角センサ２０、吸気圧センサ２２、油温センサ２６、ステップ１０）と、検出された参照パラメータに基づき、参照パラメータに対応するＭ個（４個）の重み関数値Ｗｄｉ（ｉ＝１〜４）をそれぞれ算出する重み関数値算出手段（ＥＣＵ２、オンボード同定器６０，５３０）と、算出されたＭ個の重み関数値Ｗｄｉ，Ｗｄｉ”を、Ｍ個のむだ時間ｄ前のタイミングでのＭ個の制御入力（空燃比補正係数ＫＡＦ（ｋ−４＋ｉ），制御入力Ｕａｃｔ（ｋ−４＋ｉ））にそれぞれ乗算することにより、Ｍ個の乗算値（Ｗｄｉ・ＫＡＦ（ｋ−４＋ｉ），Ｗｄｉ”・Ｕａｃｔ（ｋ−４＋ｉ））を算出するとともに、Ｍ個の乗算値の総和を修正後制御入力ＫＡＦ＿ｍｏｄ，Ｕａｃｔ＿ｍｏｄとして設定する修正後制御入力設定手段（ＥＣＵ２、オンボード同定器６０，２６０，５３０）と、修正後制御入力および制御量を用いて、外乱推定値εを算出する外乱推定値算出手段（ＥＣＵ２、適応外乱オブザーバ３７０，７４０）と、算出された外乱推定値εを用い、制御量が目標制御量になるように、制御入力を算出する制御入力算出手段（ＥＣＵ２、２自由度応答指定型コントローラ３８０，７５０）と、を備え、Ｍ個（４個）の重み関数値Ｗｄｉ，Ｗｄｉ”は、参照パラメータの所定範囲におけるＭ個の領域にそれぞれ対応し、対応する領域において値０以外の値に設定され、かつ対応する領域以外において値０に設定されており、Ｍ個の領域は隣り合う領域が互いに重なり合っており、Ｍ個の重み関数値は、重なり合う領域における参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が所定値（値１）になるように設定されていることを特徴とする。
【００２０】
この制御装置によれば、検出された参照パラメータに基づき、参照パラメータに対応するＭ個の重み関数値がそれぞれ算出され、Ｍ個の重み関数値を、Ｍ個のむだ時間前のタイミングでのＭ個の制御入力にそれぞれ乗算することにより、Ｍ個の乗算値が算出され、Ｍ個の乗算値の総和が修正後制御入力として設定される。この場合、Ｍ個の重み関数値は、参照パラメータの所定範囲におけるＭ個の領域にそれぞれ対応し、対応する領域において値０以外の値に設定され、かつ対応する領域以外において値０に設定されており、Ｍ個の領域は隣り合う領域が互いに重なり合っており、Ｍ個の重み関数値は、重なり合う領域における参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が所定値になるように設定されている。したがって、そのようなＭ個の重み関数値をＭ個のむだ時間前のタイミングでのＭ個の制御入力にそれぞれ乗算したＭ個の乗算値の総和が修正後制御入力として設定されるので、参照パラメータの変化に応じてむだ時間が連続的に変化する場合でも、修正後制御入力を、そのようなむだ時間の変化を適切に補償しながら精度よく算出することができる。特に、参照パラメータの急変に伴ってむだ時間が急変したときでも、修正後制御入力を、むだ時間の急変を適切に補償しながら、段差なく円滑に変化するように算出することができる。
【００２１】
さらに、そのように算出された修正後制御入力および制御量を用いて、外乱推定値が算出されるので、参照パラメータの変化に応じてむだ時間が連続的に変化する場合でも、そのようなむだ時間の変化を適切に補償しながら、外乱推定値を外乱を精度よく表す値として算出することができる。これに加えて、そのように算出された外乱推定値を用い、制御量が目標制御量になるように、制御入力が算出されるので、参照パラメータの変化に応じてむだ時間が連続的に変化する場合でも、制御入力を、そのようなむだ時間の変化を適切に補償しながら精度よく算出することができるとともに、外乱抑制能力すなわちロバスト性を向上させることができる。以上により、制御量と目標制御量との偏差の積分値を使用する制御アルゴリズムによって、制御入力を算出したときでも、制御量の振動的な挙動やオーバーシュート挙動の発生を回避しながら、制御量を目標制御量に精度よく制御することができる。特に、制御入力の算出アルゴリズムとして、フィードバック制御アルゴリズムを用いた場合には、フィードバックゲインを高い状態に維持することができ、高精度かつ高応答性を確保しながら、制御量を目標制御量に追従させることができる。
【００２２】
請求項７に係る発明は、請求項６に記載の制御装置１Ｃ，１Ｅにおいて、外乱推定値算出手段は、制御量の推定値である推定制御量（外乱推定用の推定検出当量比ＫＡＣＴ＿ａｄｖ、外乱推定用の推定主軸回転数ＮＭ＿ａｄｖ）を、推定制御量と修正後制御入力と外乱推定値と制御量との関係を定義したモデルを用いて算出するとともに、推定制御量と制御量との偏差が最小となるように、外乱推定値εを算出することを特徴とする。
【００２３】
この制御装置によれば、制御量の推定値である推定制御量が、推定制御量と修正後制御入力と外乱推定値と制御量との関係を定義したモデルを用いて算出される。この場合、修正後制御入力および外乱推定値は、前述したように、むだ時間の変化を適切に補償しながら精度よく算出されるので、参照パラメータの変化に応じてむだ時間が変化する場合でも、そのようなむだ時間の変化を適切に補償しながら、推定制御量を精度よく算出することができる。これに加えて、そのように算出された推定制御量と制御量との偏差が最小となるように、外乱推定値が算出されるので、外乱推定値の算出精度をさらに向上させることができ、それにより、制御量の目標制御量への制御精度をさらに向上させることができる。
【００２４】
請求項８に係る発明は、請求項１ないし７のいずれかに記載の制御装置１，１Ａ〜１Ｃにおいて、制御量は、内燃機関３の混合気の空燃比を表す値（検出当量比ＫＡＣＴ）であり、制御入力は、内燃機関３の燃料量を補正する補正係数（空燃比補正係数ＫＡＦ）であることを特徴とする。
【００２５】
この制御装置によれば、内燃機関の燃料量を補正する補正係数を制御入力として用いながら、内燃機関の混合気の空燃比を表す値を制御量として制御する場合において、請求項１ないし６の制御装置と同様の作用効果を得ることができる。
【００２６】
請求項９に係る発明は、請求項１ないし７のいずれかに記載の制御装置１Ｄ，１Ｅにおいて、制御量は、自動変速機４００の伝達トルク調整機構（クラッチ４１０）における出力回転数を表す値（主軸回転数ＮＭ）であり、制御入力は伝達トルク調整機構のアクチュエータへの入力Ｕａｃｔであることを特徴とする。
【００２７】
この制御装置によれば、伝達トルク調整機構のアクチュエータへの入力を制御入力として用いながら、自動変速機の伝達トルク調整機構における出力回転数を表す値を制御量として制御する場合において、請求項１ないし６の制御装置と同様の作用効果を得ることができる。
【図面の簡単な説明】
【００２８】
【図１】本発明の第１実施形態に係る制御装置およびそれを適用した内燃機関の概略的な構成を示す図である。
【図２】むだ時間ｄと排ガスボリュームＶｅｘの関係をモデリングした図である。
【図３】制御装置の構成を示すブロック図である。
【図４】要求トルクＴＲＱＤＲＶの算出に用いるマップの一例を示す図である。
【図５】可変むだ時間状態予測器の構成を示すブロック図である。
【図６】重み関数値Ｗｄｉの算出に用いるマップの一例を示す図である。
【図７】オンボード・スケジュールド・モデルパラメータ同定器の構成を示すブロック図である。
【図８】修正後制御入力算出部の構成を示すブロック図である。
【図９】同定値算出部の構成を示すブロック図である。
【図１０】基準モデルパラメータαｂｓの算出に用いるマップの一例を示す図である。
【図１１】重み関数値Ｗａｉの算出に用いるマップの一例を示す図である。
【図１２】制御装置におけるフィードバック制御系の構成をｚ領域で表したブロック図である。
【図１３】最適感度関数Ｓｏｐｔのゲイン曲線を表す図である。
【図１４】スライディングモード制御アルゴリズムの感度関数Ｓｓｌｄのゲイン曲線を示す図である。
【図１５】式（４２）の感度関数Ｓｄのゲイン曲線を示す図である。
【図１６】相補感度関数Ｔｄのゲイン曲線を示す図である。
【図１７】一次遅れ系のモデル化誤差Δｌのゲイン曲線を示す図である。
【図１８】式（５０）の伝達関数Ｐのボード線図である。
【図１９】式（４１）の伝達関数Ｐのボード線図である。
【図２０】空燃比制御処理を示すフローチャートである。
【図２１】モデル化誤差がないシミュレーション条件下での第１実施形態の制御装置による空燃比制御のシミュレーション結果の一例を示すタイミングチャートである。
【図２２】比較のために、モデル化誤差がないシミュレーション条件下で、制御装置における同定値αｉｄおよび予測当量比ＰＲＥ＿ＫＡＣＴの演算を停止したときの制御シミュレーション結果を示すタイミングチャートである。
【図２３】比較のために、モデル化誤差がないシミュレーション条件下で、制御装置における同定値αｉｄおよび予測当量比ＰＲＥ＿ＫＡＣＴの演算を停止しかつ感度設定パラメータβの値を変更したときの制御シミュレーション結果を示すタイミングチャートである。
【図２４】モデル化誤差が存在するシミュレーション条件下での第１実施形態の制御装置による空燃比制御のシミュレーション結果の一例を示すタイミングチャートである。
【図２５】比較のために、モデル化誤差が存在するシミュレーション条件下で、制御装置における同定値αｉｄおよび予測当量比ＰＲＥ＿ＫＡＣＴの演算を停止したときの制御シミュレーション結果を示すタイミングチャートである。
【図２６】比較のために、モデル化誤差が存在するシミュレーション条件下で、制御装置における同定値αｉｄおよび予測当量比ＰＲＥ＿ＫＡＣＴの演算を停止しかつ感度設定パラメータβの値を変更したときの制御シミュレーション結果を示すタイミングチャートである。
【図２７】比較のために、モデル化誤差が存在するシミュレーション条件下で、制御装置における同定値αｉｄの演算のみを停止したときの制御シミュレーション結果を示すタイミングチャートである。
【図２８】補正係数Ｋαｂｓの算出に用いるマップの一例を示す図である。
【図２９】重み関数値Ｗａｎｊの算出に用いるマップの一例を示す図である。
【図３０】重み関数値Ｗａａｈの算出に用いるマップの一例を示す図である。
【図３１】第２実施形態に係る制御装置の構成を示すブロック図である。
【図３２】第２実施形態の可変むだ時間状態予測器の構成を示すブロック図である。
【図３３】第２実施形態のオンボード・スケジュールド・モデルパラメータ同定器の構成を示すブロック図である。
【図３４】モデルパラメータベクトル算出部の構成を示すブロック図である。
【図３５】第３実施形態に係る制御装置の構成を示すブロック図である。
【図３６】第４実施形態に係る制御装置の構成を示すブロック図である。
【図３７】本発明の第５実施形態に係る制御装置およびそれを適用した内燃機関の駆動系の概略的な構成を示す図である。
【図３８】むだ時間ｄ”と油温Ｔｏｉｌの関係をモデリングした図である。
【図３９】クラッチコントローラの構成を示すブロック図である。
【図４０】目標クラッチすべり率Ｒｓｌｉｐ＿ｃｍｄの算出に用いるマップの一例を示す図である。
【図４１】重み関数値Ｗｄｉ”の算出に用いるマップの一例を示す図である。
【図４２】重み関数値Ｗａｉ”の算出に用いるマップの一例を示す図である。
【図４３】基準モデルパラメータαｂｓ”の算出に用いるマップの一例を示す図である。
【図４４】スロットル弁コントローラの構成を示すブロック図である。
【図４５】目標エンジントルクＴＲＱ＿ＥＮＧ＿ｃｍｄの算出に用いるマップの一例を示す図である。
【図４６】目標ＴＨ開度ＴＨ＿ｃｍｄの算出に用いるマップの一例を示す図である。
【図４７】第５実施形態の制御装置によるクラッチ制御のシミュレーション結果の一例を示すタイミングチャートである。
【図４８】第６実施形態の制御装置の構成を示すブロック図である。
【発明を実施するための形態】
【００２９】
以下、図面を参照しながら、本発明の第１実施形態に係る制御装置について説明する。本実施形態の図１に示す制御装置１は、後述する制御アルゴリズムによって、内燃機関（以下「エンジン」という）３の混合気の空燃比を制御するものであり、ＥＣＵ２を備えている。
【００３０】
エンジン３は、図示しない車両に搭載された直噴式のガソリンエンジンであり、このエンジン３には、燃料噴射弁４（１つのみ図示）が気筒ごとに取り付けられている。この燃料噴射弁４は、ＥＣＵ２に電気的に接続されており、その開弁時間および開弁タイミングがＥＣＵ２によって制御され、それにより、燃料噴射制御が実行される。この場合、通常運転状態では、空燃比がリーン側の値になるように、燃料噴射制御が実行され、それによって、エンジン３はリーンバーン運転される。
【００３１】
ＥＣＵ２には、クランク角センサ２０およびアクセル開度センサ２１がそれぞれ接続されている。このクランク角センサ２０（参照パラメータ検出手段）は、マグネットロータおよびＭＲＥピックアップで構成されており、クランクシャフト（図示せず）の回転に伴い、いずれもパルス信号であるＣＲＫ信号およびＴＤＣ信号をＥＣＵ２に出力する。
【００３２】
このＣＲＫ信号は、所定クランク角（例えば１゜）毎に１パルスが出力され、ＥＣＵ２は、このＣＲＫ信号に基づき、エンジン３の機関回転数（以下「エンジン回転数」という）ＮＥを算出する。また、ＴＤＣ信号は、各気筒のピストン（いずれも図示せず）が吸気行程のＴＤＣ位置よりも若干、手前の所定のクランク角位置にあることを表す信号であり、所定クランク角毎に１パルスが出力される。
【００３３】
また、アクセル開度センサ２１は、車両の図示しないアクセルペダルの踏み込み量（以下「アクセル開度」という）ＡＰを検出して、それを表す検出信号をＥＣＵ２に出力する。
【００３４】
一方、エンジン３の吸気通路５には、上流側から順に、スロットル弁機構６および吸気圧センサ２２が設けられている。このスロットル弁機構６は、スロットル弁６ａおよびこれを開閉駆動するＴＨアクチュエータ６ｂなどを備えている。スロットル弁６ａは、吸気通路５の途中に回動自在に設けられており、当該回動に伴う開度の変化によりスロットル弁６ａを通過する空気量を変化させる。ＴＨアクチュエータ６ｂは、ＥＣＵ２に接続されたモータにギヤ機構（いずれも図示せず）を組み合わせたものであり、ＥＣＵ２からの制御入力信号によって制御されることにより、スロットル弁６ａの開度を変化させる。
【００３５】
さらに、吸気圧センサ２２（参照パラメータ検出手段）は、吸気通路５のスロットル弁６ａよりも下流側のサージタンクの部分に配置されており、吸気通路５内の圧力（以下「吸気圧」という）を検出して、それを表す検出信号をＥＣＵ２に出力する。ＥＣＵ２は、この吸気圧センサ２２の検出信号に基づいて、吸気圧ＰＢを算出する。なお、この吸気圧ＰＢは絶対圧として算出される。
【００３６】
一方、エンジン３の排気通路１０には、上流側から順に、ＬＡＦセンサ２３、上流側三元触媒１１、酸素濃度センサ２４、下流側三元触媒１２、ユリア噴射弁１３、上流側選択還元触媒１４、ＮＨ３濃度センサ２５および下流側選択還元触媒１５が設けられている。
【００３７】
このＬＡＦセンサ２３は、ジルコニアおよび白金電極などで構成され、理論空燃比よりもリッチなリッチ領域から極リーン領域までの広範囲な空燃比の領域において、排気通路１０内を流れる排ガス中の酸素濃度をリニアに検出し、それを表す検出信号をＥＣＵ２に出力する。ＥＣＵ２は、このＬＡＦセンサ２３の検出信号の値に基づき、排ガス中の当量比を表す検出当量比ＫＡＣＴを算出する。なお、本実施形態では、検出当量比ＫＡＣＴが制御量および空燃比を表す値に相当する。
【００３８】
また、上流側三元触媒１１は、その温度が所定の活性化温度よりも高い領域にあるときに活性化し、排ガス中の有害な未燃成分を浄化する。さらに、下流側三元触媒１２は、上流側三元触媒１１と同じタイプのものであり、上流側選択還元触媒１４に流入する排ガス成分をＮＯｘ浄化に最適な状態に調整し、上流側選択還元触媒１４において高いＮＯｘ浄化率を確保することを目的として、上流側選択還元触媒１４の上流側に設けられている。なお、下流側三元触媒１２として、上流側三元触媒１１と異なるタイプのもの、例えば、リーンバーン運転中のＨＣやＣＯの酸化能力をより高めたものや、ＮＯを酸化してＮＯ２化する能力をより高めたものを用いてもよい。
【００３９】
さらに、酸素濃度センサ２４は、ジルコニアおよび白金電極などで構成され、上流側三元触媒１１を通過した排ガス中の酸素濃度に基づいた出力をＥＣＵ２に送る。この酸素濃度センサ２４の出力は、理論空燃比よりもリッチな混合気が燃焼したときには、ハイレベルの電圧値（例えば０．８Ｖ）となり、混合気がリーンのときには、ローレベルの電圧値（例えば０．２Ｖ）となるとともに、混合気が理論空燃比付近のときには、ハイレベルとローレベルの間の所定の目標値（例えば０．６Ｖ）となる。
【００４０】
一方、ユリア噴射弁１３は、ＥＣＵ２に電気的に接続されており、ＥＣＵ２からの制御入力信号によって開弁状態に駆動されると、ユリアタンク（図示せず）からの尿素水を排気通路１０内に噴射する。この場合、ユリア噴射弁１３から噴射された尿素水の尿素は、その一部が、排ガスの熱および選択還元触媒１４との接触によって、アンモニアに変化する。
【００４１】
また、上流側選択還元触媒１４は、還元剤としての尿素（Ｕｒｅａ）が存在する雰囲気下で、排ガス中の窒素酸化物（ＮＯｘ）を選択的に還元する。上流側選択還元触媒１４では、そのＮＯｘの還元作用において、尿素水の噴射時に尿素から変化したアンモニアも、尿素と一緒に消費されるとともに、消費されなかった分のアンモニアは、上流側選択還元触媒１４内に貯蔵される。
【００４２】
さらに、下流側選択還元触媒１５は、上流側選択還元触媒１４と同じタイプのものであり、排ガス中のＮＯｘを浄化することに加えて、上流側選択還元触媒１５を通過したアンモニアを捕捉することを目的として、上流側選択還元触媒１４の下流側に設けられている。本実施形態の場合、以上のユリア噴射弁１３および選択還元触媒１４，１５によって、尿素ＳＣＲ装置が構成されている。なお、下流側選択還元触媒１５として、上流側選択還元触媒１４と比べて低温でのＮＯｘ浄化性能をより高めたもの、例えば、Ｃｕゼオライトタイプのものや、後側に酸化触媒をゾーンコート化したものを用いてもよい。
【００４３】
さらに、ＮＨ３濃度センサ２５は、上流側選択還元触媒１４を通過した排ガス中のアンモニア濃度を検出し、それを表す検出信号をＥＣＵ２に出力する。ＥＣＵ２は、このＮＨ３濃度センサ２５の検出信号に基づき、ユリア噴射弁１３を介してユリア噴射量を制御し、それにより、尿素ＳＣＲ装置によるＮＯｘ浄化率やＮＯｘ浄化量を制御する。
【００４４】
一方、ＥＣＵ２は、ＣＰＵ、ＲＡＭ、ＲＯＭ、Ｉ／Ｏインターフェースおよび駆動回路（いずれも図示せず）などからなるマイクロコンピュータで構成されており、前述した各種のセンサ２０〜２５の検出信号などに応じて、エンジン３の運転状態を判別するとともに、後述する空燃比制御処理などを実行する。
【００４５】
なお、本実施形態では、ＥＣＵ２が、目標制御量設定手段、参照パラメータ検出手段、予測値算出手段、重み関数値算出手段、予測制御量設定手段、制御入力算出手段、修正後制御入力設定手段、同定手段、および外乱推定値算出手段に相当する。
【００４６】
次に、本実施形態の制御装置１について説明する。まず、本実施形態の制御装置１で用いる制御対象モデルについて説明する。エンジン３の燃料噴射弁４からＬＡＦセンサ２３までの系を、空燃比補正係数ＫＡＦを制御入力とし、検出当量比ＫＡＣＴを制御量とする一次遅れ系の制御対象と見なしてモデリングすると、下式（１）が得られる。この場合、空燃比補正係数ＫＡＦは、後述する制御アルゴリズムを用いて、当量比と同じ次元の値として算出される。
【数１】
0005400743.tif 000002
【００４７】
この式（１）のαはモデルパラメータである。また、式（１）において、記号（ｋ）付きの各離散データは、所定の制御周期ΔＴ（本実施形態ではＴＤＣ信号の発生周期）でサンプリングまたは算出されたデータであることを示しており、記号ｋ（ｋは正の整数）は各離散データのサンプリングまたは算出サイクルの順番を表している。この点は、以下の離散データにおいても同様である。なお、以下の説明では、各離散データにおける記号（ｋ）を適宜省略する。
【００４８】
上式（１）の場合、空燃比補正係数ＫＡＦと検出当量比ＫＡＣＴとの間に存在するむだ時間ｄが考慮されていないので、このむだ時間ｄを上式（１）に反映させると、下式（２）が得られる。この式（２）を制御対象モデルとして用いた理由については後述する。
【数２】
0005400743.tif 000003
【００４９】
この場合、上式（２）のむだ時間ｄは、エンジン３の運転状態に応じて変化するものであり、このむだ時間ｄと排ガスボリュームＶｅｘとの関係をモデリング（マッピング）すると、図２に示すモデル（マップ）が得られる。この排ガスボリュームＶｅｘ（参照パラメータ）は、排ガスの空間速度に相当する値であり、具体的には、エンジン回転数ＮＥおよび吸気圧ＰＢに応じて、図示しないマップを検索することにより算出される。
【００５０】
図２において、Ｖｅｘ１〜Ｖｅｘ４，ＶｅｘＭＡＸは、排ガスボリュームＶｅｘの所定値であり、０＜Ｖｅｘ１＜Ｖｅｘ２＜Ｖｅｘ３＜Ｖｅｘ４＜ＶｅｘＭＡＸが成立するように設定されている。また、所定値ＶｅｘＭＡＸは、エンジン３の運転中に排ガスボリュームＶｅｘが変化し得る領域の最大値に設定されている。言い換えれば、排ガスボリュームＶｅｘは、エンジン３の運転中、０〜ＶｅｘＭＡＸの領域内で変化する特性を有している。
【００５１】
本実施形態の制御装置１では、以上のむだ時間ｄを含む上式（２）の制御対象モデルを用いて、以下に述べるように、空燃比補正係数ＫＡＦなどの各種の演算値が算出される。図３に示すように、制御装置１は、目標当量比算出部３０、可変むだ時間状態予測器（以下「状態予測器」という）４０、オンボード・スケジュールド・モデルパラメータ同定器（以下「オンボード同定器」という）６０、および周波数整形コントローラ１３０を備えており、これらはいずれもＥＣＵ２によって構成されている。
【００５２】
この目標当量比算出部３０では、前述した検出当量比ＫＡＣＴの目標となる値として、目標当量比ＫＣＭＤが算出される。具体的には、エンジン回転数ＮＥおよびアクセル開度ＡＰに応じて、図示しないマップを検索することにより、要求トルクＴＲＱＤＲＶを算出し、この要求トルクＴＲＱＤＲＶおよびエンジン回転数ＮＥに応じて、図４に示すマップを検索することにより、目標当量比ＫＣＭＤが算出される。同図において、ＫＣＭＤ１〜ＫＣＭＤ４は、目標当量比ＫＣＭＤの所定値であり、ＫＣＭＤ１＝１，ＫＣＭＤ１＞ＫＣＭＤ２＞ＫＣＭＤ３＞ＫＣＭＤ４が成立するように設定される。
【００５３】
また、状態予測器４０では、後述する予測アルゴリズムを用いて、検出当量比ＫＡＣＴの予測値として、予測当量比ＰＲＥ＿ＫＡＣＴが算出され、オンボード同定器６０では、後述する同定アルゴリズムを用いて、前述したモデルパラメータαをオンボード同定した値として、同定値αｉｄが算出される。さらに、周波数整形コントローラ１３０では、後述する制御アルゴリズムを用いて、制御入力としての空燃比補正係数ＫＡＦが算出される。
【００５４】
なお、本実施形態では、目標当量比算出部３０が目標制御量設定手段に相当し、目標当量比ＫＣＭＤが目標制御量に相当する。また、状態予測器４０が、予測値算出手段、重み関数値算出手段および予測制御量設定手段に相当し、予測当量比ＰＲＥ＿ＫＡＣＴが予測制御量に相当する。さらに、オンボード同定器６０が、修正後制御入力設定手段、同定手段および重み関数値算出手段に相当し、周波数整形コントローラ１３０が制御入力算出手段に相当する。
【００５５】
次に、前述した状態予測器４０について説明する。この状態予測器４０は、以下に述べる予測アルゴリズムにより、予測当量比ＰＲＥ＿ＫＡＣＴを算出するものであり、この予測当量比ＰＲＥ＿ＫＡＣＴは、現在の制御系におけるむだ時間ｄが経過した後の制御タイミングにおける検出当量比ＫＡＣＴを予測した値に相当する。
【００５６】
図５に示すように、状態予測器４０は、３つの遅延要素４１〜４３と、増幅器４４と、３つの予測値算出部４５〜４７と、４つの重み関数値算出部４８〜５１と、４つの乗算器５２〜５５と、加算器５６とを備えている。
【００５７】
まず、増幅器４４で、下式（３）により、第０予測値ＰＲＥ＿ＫＡＣＴ＿０が算出される。すなわち、第０予測値ＰＲＥ＿ＫＡＣＴ＿０は、むだ時間ｄ＝０のときの検出当量比ＫＡＣＴ（ｋ）として算出される。
【数３】
0005400743.tif 000004
【００５８】
また、第１予測値算出部４５では、遅延要素４１で１制御サイクル分遅延された空燃比補正係数の値ＫＡＦ（ｋ−１）を用い、下式（４）により、第１予測値ＰＲＥ＿ＫＡＣＴ＿１が算出される。
【数４】
0005400743.tif 000005
【００５９】
この第１予測値ＰＲＥ＿ＫＡＣＴ＿１は、ｄ＝１のときのむだ時間ｄが経過したときの検出当量比ＫＡＣＴを予測した値に相当する。なお、上式（４）の導出手法については後述する。
【００６０】
さらに、第２予測値算出部４６では、２つの遅延要素４１，４２でそれぞれ１，２制御サイクル分ずつ、遅延された空燃比補正係数の値ＫＡＦ（ｋ−１），ＫＡＦ（ｋ−２）を用い、下式（５）により、第２予測値ＰＲＥ＿ＫＡＣＴ＿２が算出される。
【数５】
0005400743.tif 000006
【００６１】
この第２予測値ＰＲＥ＿ＫＡＣＴ＿２は、ｄ＝２のときのむだ時間ｄが経過したときの検出当量比ＫＡＣＴを予測した値に相当する。なお、上式（５）の導出手法については後述する。
【００６２】
また、第３予測値算出部４７では、３つの遅延要素４１〜４３でそれぞれ１〜３制御サイクル分ずつ、遅延された空燃比補正係数の値ＫＡＦ（ｋ−１），ＫＡＦ（ｋ−２），ＫＡＦ（ｋ−３）を用い、下式（６）により、第３予測値ＰＲＥ＿ＫＡＣＴ＿３が算出される。
【数６】
0005400743.tif 000007
【００６３】
この第３予測値ＰＲＥ＿ＫＡＣＴ＿３は、ｄ＝３のときのむだ時間ｄが経過したときの検出当量比ＫＡＣＴを予測した値に相当する。なお、上式（６）の導出手法については後述する。
【００６４】
さらに、４つの重み関数値算出部４８〜５１では、排ガスボリュームＶｅｘに応じて、図６に示すマップを検索することにより、４つの重み関数値Ｗｄ１〜Ｗｄ４がそれぞれ算出される。同図に示すように、４つの重み関数値Ｗｄ１〜Ｗｄ４はそれぞれ、排ガスボリュームＶｅｘが変化し得る領域を０≦Ｖｅｘ≦Ｖｅｘ２，Ｖｅｘ１≦Ｖｅｘ≦Ｖｅｘ３，Ｖｅｘ２≦Ｖｅｘ≦Ｖｅｘ４，Ｖｅｘ３≦Ｖｅｘ≦ＶｅｘＭＡＸの４つの領域に区分した場合において、これらの４つの領域に対応するように設定されているとともに、対応する領域では、値１以下の正の値に設定され、対応する領域以外では、値０に設定されている。
【００６５】
具体的には、重み関数値Ｗｄ１は、これが対応する領域（０≦Ｖｅｘ≦Ｖｅｘ２）では、Ｖｅｘ≦Ｖｅｘ１のときの値１を最大値として、Ｖｅｘ１＜Ｖｅｘの領域で排ガスボリュームＶｅｘが大きいほど、より小さい正の値に設定されているとともに、それ以外の領域では、値０に設定されている。また、重み関数値Ｗｄ２は、これが対応する領域（Ｖｅｘ１≦Ｖｅｘ≦Ｖｅｘ３）では、Ｖｅｘ＝Ｖｅｘ２のときの値１を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値０に設定されている。
【００６６】
さらに、重み関数値Ｗｄ３は、これが対応する領域（Ｖｅｘ２≦Ｖｅｘ≦Ｖｅｘ４）では、Ｖｅｘ＝Ｖｅｘ３のときの値１を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値０に設定されている。一方、重み関数値Ｗｄ４は、これが対応する領域（Ｖｅｘ３≦Ｖｅｘ≦ＶｅｘＭＡＸ）では、Ｖｅｘ４≦Ｖｅｘのときの値１を最大値として、排ガスボリュームＶｅｘが大きいほど、より大きい正の値に設定されているとともに、それ以外の領域では、値０に設定されている。
【００６７】
以上に加えて、４つの重み関数値Ｗｄｉ（ｉ＝１〜４）がそれぞれ対応する４つの領域は、上述したような、隣り合う領域が互いに重なり合うように設定されており、これらの互いに重なり合う領域における排ガスボリュームＶｅｘのいずれかの値に対応する重み関数値Ｗｄｉの値の和は、各重み関数値Ｗｄｉにおける最大値１と等しくなるように設定されている。
【００６８】
さらに、図６と前述した図２を比較すると明らかなように、互いに重なり合う３つの領域は、むだ時間ｄの勾配が一定状態に保持される３つの領域に対応するように設定されている。これに加えて、重み関数値Ｗｄ１はむだ時間ｄ＝３に対して、重み関数値Ｗｄ２はむだ時間ｄ＝２に対して、重み関数値Ｗｄ３はむだ時間ｄ＝１に対して、重み関数値Ｗｄ４はむだ時間ｄ＝０に対して重みが最も大きくなるようにそれぞれ設定されている。
【００６９】
また、乗算器５２では、重み関数値Ｗｄ４に第０予測値ＰＲＥ＿ＫＡＣＴ＿０を乗算することにより、乗算値Ｗｄ４・ＰＲＥ＿ＫＡＣＴ＿０が算出され、乗算器５３では、重み関数値Ｗｄ３に第１予測値ＰＲＥ＿ＫＡＣＴ＿１を乗算することにより、乗算値Ｗｄ３・ＰＲＥ＿ＫＡＣＴ＿１が算出される。さらに、乗算器５３では、重み関数値Ｗｄ２に第２予測値ＰＲＥ＿ＫＡＣＴ＿２を乗算することにより、乗算値Ｗｄ２・ＰＲＥ＿ＫＡＣＴ＿２が算出され、乗算器５４では、重み関数値Ｗｄ１に第３予測値ＰＲＥ＿ＫＡＣＴ＿３を乗算することにより、乗算値Ｗｄ１・ＰＲＥ＿ＫＡＣＴ＿３が算出される。
【００７０】
そして、加算器５６で、以上の４つの乗算値を互いに加算することにより、予測当量比ＰＲＥ＿ＫＡＣＴが算出される。すなわち、下式（７）により、予測当量比ＰＲＥ＿ＫＡＣが算出される。
【数７】
0005400743.tif 000008
【００７１】
以上のように、前述した４つの重み関数値Ｗｄｉを、４つの予測値ＰＲＥ＿ＫＡＣＴ＿４−ｉにそれぞれ乗算した値の総和として、予測当量比ＰＲＥ＿ＫＡＣＴが算出されるので、排ガスボリュームＶｅｘの変化に応じて、むだ時間ｄが図２に示すように値０から値３の間で連続的に変化したときでも、そのようなむだ時間ｄの変化を適切に反映させながら、予測当量比ＰＲＥ＿ＫＡＣＴを円滑に段差なく変化するような値として算出することができる。
【００７２】
なお、前述した第１〜３予測値ＰＲＥ＿ＫＡＣＴ＿１〜３の算出式（４）〜（６）は、以下に述べるように導出される。まず、前述した式（２）において、ｄ＝１とすると、下式（８）が得られる。
【数８】
0005400743.tif 000009
【００７３】
上式（８）において、右辺のＫＡＣＴ（ｋ＋１）をＰＲＥ＿ＫＡＣＴ＿１（ｋ）に、左辺のαをαｉｄ（ｋ）にそれぞれ置き換えると、前述した式（４）が得られる。
【００７４】
また、前述した式（２）において、ｄ＝２とすると、下式（９）が得られる。
【数９】
0005400743.tif 000010
【００７５】
上式（９）の各変数を１制御サイクル分、未来側にシフトさせると、下式（１０）が得られる。
【数１０】
0005400743.tif 000011
【００７６】
上式（１０）に上式（９）を代入すると、下式（１１）が得られる。
【数１１】
0005400743.tif 000012
【００７７】
上式（１１）において、右辺のＫＡＣＴ（ｋ＋２）をＰＲＥ＿ＫＡＣＴ＿２（ｋ）に、左辺のαをαｉｄ（ｋ）にそれぞれ置き換えると、前述した式（５）が得られる。
【００７８】
また、前述した式（２）において、ｄ＝３とすると、下式（１２）が得られる。
【数１２】
0005400743.tif 000013
【００７９】
上式（１２）の各変数を１制御サイクル分、未来側にシフトさせると、下式（１３）が得られる。
【数１３】
0005400743.tif 000014
【００８０】
上式（１３）に上式（１２）を代入すると、下式（１４）が得られる。
【数１４】
0005400743.tif 000015
【００８１】
さらに、上式（１３）の各変数を１制御サイクル分、未来側にシフトさせると、下式（１５）が得られる。
【数１５】
0005400743.tif 000016
【００８２】
この式（１５）に上式（１４）を代入すると、下式（１６）が得られる。
【数１６】
0005400743.tif 000017
【００８３】
この式（１６）において、右辺のＫＡＣＴ（ｋ＋３）をＰＲＥ＿ＫＡＣＴ＿３（ｋ）に、左辺のαをαｉｄ（ｋ）にそれぞれ置き換えると、前述した式（６）が得られる。
【００８４】
次に、前述したオンボード同定器６０について説明する。このオンボード同定器６０の場合、本実施形態の制御対象のように、むだ時間ｄが排ガスボリュームＶｅｘに応じて連続的に変化する場合において、以下に述べる拘束条件付きスケジュールド補正型の同定アルゴリズムを用いて、そのようなむだ時間ｄの変化を反映させながら、同定値αｉｄが算出される。なお、このオンボード同定器６０の同定アルゴリズムは、後述するように、前述した式（２）の右辺の値ＫＡＦ（ｋ−ｄ）を後述する修正後制御入力ＫＡＦ＿ｍｏｄ（ｋ）に置き換えた修正後モデル（後述する式（３０））に基づいて導出される。
【００８５】
このオンボード同定器６０は、図７に示すように、修正後制御入力算出部７０、３つの遅延要素６１〜６３、合成信号値算出部６４、推定合成信号値算出部６５、同定ゲイン算出部６６、減算器６７、乗算器６８および同定値算出部９０を備えている。
【００８６】
まず、修正後制御入力算出部７０について説明する。この修正後制御入力算出部７０は、修正後制御入力ＫＡＦ＿ｍｏｄを算出するものであり、図８に示すように、修正後制御入力算出部７０は、３つの遅延要素７１〜７３、４つの重み関数値算出部７４〜７７、４つの乗算器７８〜８１および加算器８２を備えている。
【００８７】
まず、４つの重み関数値算出部７４〜７７で、前述した４つの重み関数値算出部４８〜５１と同様に、排ガスボリュームＶｅｘに応じて、前述した図６に示すマップを検索することにより、４つの重み関数値Ｗｄ１〜Ｗｄ４がそれぞれ算出される。
【００８８】
また、乗算器７８では、空燃比補正係数の今回値ＫＡＦ（ｋ）を重み関数値Ｗｄ４（ｋ）に乗算することにより、乗算値Ｗｄ４（ｋ）・ＫＡＦ（ｋ）が算出され、乗算器７９では、遅延要素７１で１制御サイクル分遅延された空燃比補正係数の値ＫＡＦ（ｋ−１）を重み関数値Ｗｄ３（ｋ）に乗算することにより、乗算値Ｗｄ３（ｋ）・ＫＡＦ（ｋ−１）が算出される。
【００８９】
さらに、乗算器８０では、２つの遅延要素７１，７２によって２制御サイクル分遅延された空燃比補正係数の値ＫＡＦ（ｋ−２）を重み関数値Ｗｄ２（ｋ）に乗算することにより、乗算値Ｗｄ２（ｋ）・ＫＡＦ（ｋ−２）が算出され、乗算器８１では、３つの遅延要素７１〜７３によって３制御サイクル分遅延された空燃比補正係数の値ＫＡＦ（ｋ−３）を重み関数値Ｗｄ１（ｋ）に乗算することにより、乗算値Ｗｄ１（ｋ）・ＫＡＦ（ｋ−３）が算出される。
【００９０】
そして、加算器８２で、以上の４つの乗算値を用い、下式（１７）により、修正後制御入力ＫＡＦ＿ｍｏｄが算出される。
【数１７】
0005400743.tif 000018
【００９１】
図７に戻り、合成信号値算出部６４では、検出当量比ＫＡＣＴと、遅延要素６１によって１制御サイクル分遅延された検出当量比の値ＫＡＣＴ（ｋ−１）を用い、下式（１８）により、合成信号値Ｗ＿ａｃｔが算出される。
【数１８】
0005400743.tif 000019
【００９２】
また、推定合成信号値算出部６５では、遅延要素６１によって１制御サイクル分遅延された検出当量比の値ＫＡＣＴ（ｋ−１）と、遅延要素６２によって１制御サイクル分遅延された修正後制御入力の値ＫＡＦ＿ｍｏｄ（ｋ−１）とを用いて、下式（１９）により、偏差ζ’（ｋ−１）を算出した後、これと遅延要素６３によって１制御サイクル分遅延された同定値の値αｉｄ（ｋ−１）を用い、下式（２０）により、推定合成信号値Ｗ＿ｈａｔが算出される。
【数１９】
0005400743.tif 000020
【数２０】
0005400743.tif 000021
【００９３】
さらに、減算器６７では、下式（２１）により、同定誤差ｅｉｄ’が算出される。
【数２１】
0005400743.tif 000022
【００９４】
一方、同定ゲイン算出部６８では、下式（２２），（２３）により、同定ゲインＫｐ’が算出される。この同定ゲインＫｐ’は、同定値αｉｄの修正方向（正負）および修正量を規定するものである。
【数２２】
0005400743.tif 000023
【数２３】
0005400743.tif 000024
【００９５】
上式（２２）のゲインＰ’（ｋ）の初期値Ｐ’（０）は、下式（２４）のように定義される。
【数２４】
0005400743.tif 000025
ここで、Ｐ０は所定値に設定されている。
【００９６】
また、上式（２２）において、λ１，λ２は重みパラメータであり、これらの値λ１、λ２を下記のように設定することにより、同定アルゴリズムとして、３つのアルゴリズムのいずれかを選択することができる。
λ１＝１，λ２＝０  ；固定ゲインアルゴリズム
λ１＝１，λ２＝１  ；最小２乗法アルゴリズム 
λ１＝λ，λ２＝１  ；重み付き最小２乗法アルゴリズム
ここで、λは０＜λ＜１に設定される所定値である。本実施形態の場合、同定精度および制御精度を適切に確保するために、重み付き最小２乗法アルゴリズムを用いている。
【００９７】
さらに、乗算器６８で、同定ゲインＫｐ’と同定誤差ｅｉｄ’の乗算値Ｋｐ’・ｅｉｄ’が算出される。
【００９８】
次いで、同定値算出部９０では、上記の乗算値Ｋｐ’・ｅｉｄ’と排ガスボリュームＶｅｘを用いて、以下に述べるように同定値αｉｄが算出される。図９に示すように、この同定値算出部９０は、基準モデルパラメータ算出部９１、４つの重み関数値算出部９２〜９５、８つの乗算器９６〜１０３、５つの加算器１０４〜１０８、４つの遅延要素１０９〜１１２および４つの増幅器１１３〜１１６を備えている。
【００９９】
まず、基準モデルパラメータ算出部９１では、排ガスボリュームＶｅｘに応じて、図１０に示すマップを検索することにより、基準モデルパラメータαｂｓが算出される。同図において、Ｖｅｘ５〜Ｖｅｘ８は、排ガスボリュームＶｅｘの所定値であり、０＜Ｖｅｘ５＜Ｖｅｘ６＜Ｖｅｘ７＜Ｖｅｘ８＜ＶｅｘＭＡＸが成立するように設定されている。このマップでは、基準モデルパラメータαｂｓは、排ガスボリュームＶｅｘが大きいほど、より大きい値に設定されている。これは、排ガスボリュームＶｅｘが大きいほど、ＬＡＦセンサ２３におけるセンサカバーの孔を介した排ガスの交換が促進され、ＬＡＦセンサ２３の遅れ特性が小さくなることで、空燃比補正係数ＫＡＦが検出当量比ＫＡＣＴに及ぼす影響の度合がより大きくなることによる。
【０１００】
また、４つの重み関数値算出部９２〜９５では、排ガスボリュームＶｅｘに応じて、図１１に示すマップを検索することにより、４つの重み関数値Ｗａ１〜Ｗａ４がそれぞれ算出される。同図に示すように、４つの重み関数値Ｗａ１〜Ｗａ４はそれぞれ、排ガスボリュームＶｅｘが変化し得る領域を０≦Ｖｅｘ≦Ｖｅｘ６，Ｖｅｘ５≦Ｖｅｘ≦Ｖｅｘ７，Ｖｅｘ６≦Ｖｅｘ≦Ｖｅｘ８，Ｖｅｘ７≦Ｖｅｘ≦ＶｅｘＭＡＸの４つの領域に区分した場合において、これらの４つの領域に対応するように設定されているとともに、対応する領域では、値１以下の正の値に設定され、対応する領域以外では、値０に設定されている。
【０１０１】
すなわち、重み関数値Ｗａ１は、これが対応する領域（０≦Ｖｅｘ≦Ｖｅｘ６）では、Ｖｅｘ≦Ｖｅｘ５のときの値１を最大値として、排ガスボリュームＶｅｘが大きいほど、より小さい正の値に設定されているとともに、それ以外の領域では、値０に設定されている。また、重み関数値Ｗａ２は、これが対応する領域（Ｖｅｘ５≦Ｖｅｘ≦Ｖｅｘ７）では、Ｖｅｘ＝Ｖｅｘ６のときの値１を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値０に設定されている。
【０１０２】
さらに、重み関数値Ｗａ３は、これが対応する領域（Ｖｅｘ６≦Ｖｅｘ≦Ｖｅｘ８）では、Ｖｅｘ＝Ｖｅｘ７のときの値１を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値０に設定されている。一方、重み関数値Ｗａ４は、これが対応する領域（Ｖｅｘ７≦Ｖｅｘ≦ＶｅｘＭＡＸ）では、Ｖｅｘ８≦Ｖｅｘのときの値１を最大値として、排ガスボリュームＶｅｘが大きいほど、より大きい正の値に設定されているとともに、それ以外の領域では、値０に設定されている。
【０１０３】
以上に加えて、４つの重み関数値Ｗａｉ（ｉ＝１〜４）がそれぞれ対応する４つの領域は、上述したような、隣り合う領域が互いに重なり合うように設定されており、これらの互いに重なり合う領域における排ガスボリュームＶｅｘのいずれかの値に対応する重み関数値Ｗａｉの値の和は、各重み関数値Ｗａｉにおける最大値１と等しくなるように設定されている。図１１と前述した図１０を比較すると明らかなように、これらの互いに重なり合う３つの領域は、基準モデルパラメータαｂｓの勾配が一定状態に保持される３つの領域に対応するように設定されている。
【０１０４】
また、乗算器９６で、重み関数値Ｗａ１に値Ｋｐ’・ｅｉｄ’を乗算することにより、乗算値Ｗａ１・Ｋｐ’・ｅｉｄ’が算出され、遅延要素１０９で１制御サイクル分遅延された修正項ｄα１（ｋ−１）に、増幅器１１３でゲイン係数Ｈ（ｋ）を乗算することにより、乗算値Ｈ（ｋ）・ｄα１（ｋ−１）が算出される。なお、このゲイン係数Ｈについては後述する。そして、加算器１０４で、値Ｗａ１・Ｋｐ’・ｅｉｄ’に値Ｈ（ｋ）・ｄα１（ｋ−１）を加算することにより、修正項ｄα１が算出される。
【０１０５】
さらに、乗算器９７で、重み関数値Ｗａ２に値Ｋｐ’・ｅｉｄ’を乗算することにより、乗算値Ｗａ２・Ｋｐ’・ｅｉｄ’が算出され、遅延要素１１０で１制御サイクル分遅延された修正項ｄα２（ｋ−１）に、増幅器１１４でゲイン係数Ｈ（ｋ）を乗算することにより、乗算値Ｈ（ｋ）・ｄα２（ｋ−１）が算出される。そして、加算器１０５で、値Ｗａ２・Ｋｐ’・ｅｉｄ’に値Ｈ（ｋ）・ｄα２（ｋ−１）を加算することにより、修正項ｄα２が算出される。
【０１０６】
また、乗算器９８で、重み関数値Ｗａ３に値Ｋｐ’・ｅｉｄ’を乗算することにより、乗算値Ｗａ３・Ｋｐ’・ｅｉｄ’が算出され、遅延要素１１１で１制御サイクル分遅延された修正項ｄα３（ｋ−１）に、増幅器１１５でゲイン係数Ｈ（ｋ）を乗算することにより、乗算値Ｈ（ｋ）・ｄα３（ｋ−１）が算出される。そして、加算器１０６で、値Ｗａ３・Ｋｐ’・ｅｉｄ’に値Ｈ（ｋ）・ｄα３（ｋ−１）を加算することにより、修正項ｄα３が算出される。
【０１０７】
さらに、乗算器９９で、重み関数値Ｗａ４に値Ｋｐ’・ｅｉｄ’を乗算することにより、乗算値Ｗａ４・Ｋｐ’・ｅｉｄ’が算出され、遅延要素１１２で１制御サイクル分遅延された修正項ｄα４（ｋ−１）に、増幅器１１６でゲイン係数Ｈ（ｋ）を乗算することにより、乗算値Ｈ（ｋ）・ｄα４（ｋ−１）が算出される。そして、加算器１０７で、値Ｗａ４・Ｋｐ’・ｅｉｄ’に値Ｈ（ｋ）・ｄα４（ｋ−１）を加算することにより、修正項ｄα４が算出される。
【０１０８】
以上の増幅器１１３〜１１６では、ゲイン係数Ｈが以下の式（２５）〜（２７）に示すように算出される。
【数２５】
0005400743.tif 000026
【数２６】
0005400743.tif 000027
【数２７】
0005400743.tif 000028
【０１０９】
上式（２５）〜（２７）において、α＿Ｌは所定の下限値であり、α＿Ｈは所定の上限値である。また、η’は忘却係数であり、０＜η’≦１が成立するように設定される。同定値αｉｄの算出において、この忘却係数η’を用いた理由は、エンジン３の定常運転状態が長時間継続した場合、同定値αｉｄが増大化し、不適切な値となるおそれがあるので、それを回避するためである。さらに、上式（２６）に示すように、同定値αｉｄが下限値α＿Ｌと上限値α＿Ｈとの間の領域にあるときに、忘却係数η’による忘却効果を停止している理由は、このオンボード同定器６０で用いている同定アルゴリズムの場合、後述する同定条件１（拘束条件）を満たすように、同定値αｉｄを常に同定することができるので、拘束条件を満たすことを目的として、同定値αｉｄを後述する基準モデルパラメータαｂｓの近傍に無理に拘束する必要がないことによる。
【０１１０】
また、以上の４つの加算器１０４〜１０７における演算は、下式（２８）で表される。
【数２８】
0005400743.tif 000029
【０１１１】
さらに、乗算器１００〜１０３ではそれぞれ、４つの重み関数値Ｗａｉを４つの修正項ｄαｉに乗算することにより、４つの乗算値Ｗａｉ・ｄαｉが算出される。
【０１１２】
そして、加算器１０８で、下式（２９）により、同定値αｉｄが最終的に算出される。
【数２９】
0005400743.tif 000030
【０１１３】
以上のように、オンボード同定器６０において、修正後制御入力ＫＡＦ＿ｍｏｄが、４つの制御タイミングでの検出当量比ＫＡＣＴに４つの重み関数値Ｗｄｉをそれぞれ乗算した値の総和として算出され、４つの修正項ｄαｉが、この修正後制御入力ＫＡＦ＿ｍｏｄを用いて算出した同定誤差ｅｉｄ’と同定ゲインＫｐ’との積Ｋｐ’・ｅｉｄ’に４つの重み関数値Ｗａｉをそれぞれ乗算した値の総和として算出される。そして、その総和を基準モデルパラメータαｂｓに加算することにより、同定値αｉｄが算出されるので、排ガスボリュームＶｅｘの変化に応じて遅れ特性やむだ時間ｄが連続的に変化したときでも、２種類の重み関数値Ｗｄｉ，Ｗａｉの効果によって、その影響を抑制しながら、同定値αｉｄを円滑に変化するような値として同定することができる。
【０１１４】
同定値αｉｄの算出において、以上の式（１７）〜（２９）の同定アルゴリズムを用いた理由は以下による。まず、本実施形態の制御装置１における制御系は、空燃比補正係数ＫＡＦを制御入力とし、検出当量比ＫＡＣＴを制御量とするものであり、外乱がない状態では、定常偏差を生じない系となる。それに起因して、前述した式（２）の制御対象モデルの場合、入力と出力との間に定常偏差を生じないようにするために、入力項および出力項の乗算係数すなわちモデルパラメータα，１−αは、両者の和が値１になるように設定されている。
【０１１５】
この場合、２つのモデルパラメータα，１−αは、互いに独立した値を取ることができず、一方の増大に従って他方が減少するという互いに拘束し合う関係にある。したがって、これらのモデルパラメータα，１−αを同定する場合、一方の増大に従って他方が減少するという互いに拘束し合う条件（以下「拘束条件」という）を満たすように、２つのモデルパラメータα，１−αを同定する必要があり、以下、これを「同定条件１」という。ここで、最小２乗法などの一般的な同定アルゴリズムをそのまま用いた場合、この同定条件１を満たすのは困難である。
【０１１６】
これに加えて、前述したように、遅れ特性やむだ時間ｄが排ガスボリュームＶｅｘに応じて変化する特性を有しているので、最小２乗法などの一般的な同定アルゴリズムをそのまま用いた場合、遅れ特性やむだ時間ｄの変化を反映させながら、２つのモデルパラメータα，１−αを同定することができず、結果的に、２つのモデルパラメータα，１−αの同定精度が低下してしまう。したがって、同定精度を高めるために、遅れ特性やむだ時間ｄが変化した場合でも、それを適切に反映させながら、モデルパラメータα，１−αを同定する必要があり、以下、これを「同定条件２」という。
【０１１７】
まず、上記の同定条件２を満たすために、前述した式（２）に代えて、下式（３０）を制御対象モデルとして用いる。
【数３０】
0005400743.tif 000031
【０１１８】
この式（３０）は、前述した式（２）の右辺の値ＫＡＦ（ｋ−ｄ）を値ＫＡＦ＿ｍｏｄ（ｋ）に置き換えたものに相当する。この修正後制御入力ＫＡＦ＿ｍｏｄは、前述した式（１７）に示すように、４つの重み関数値Ｗｄｉと４つの空燃比補正係数ＫＡＦとの乗算値の和として算出され、４つの重み関数値Ｗｄｉは、前述した手法により算出されるので、修正後制御入力ＫＡＦ＿ｍｏｄを、むだ時間ｄが変化したときでも、それを適切に反映させながら算出することができる。これに加えて、重み関数値Ｗａｉを用いることによって、遅れ特性の変化を反映させながら、４つの修正項ｄαｉを算出することができる。すなわち、上記の同定条件２を満たすことができる。
【０１１９】
上式（３０）を変形すると、下式（３１）が得られる。
【数３１】
0005400743.tif 000032
【０１２０】
上式（３１）の左辺および右辺をそれぞれ、下式（３２），（３３）に示すように、合成信号値Ｗ＿ａｃｔおよび推定合成信号値Ｗ＿ｈａｔとして定義する。
【数３２】
0005400743.tif 000033
【数３３】
0005400743.tif 000034
【０１２１】
このように定義した場合、上記の同定条件１を満たすには、合成信号値Ｗ＿ａｃｔおよび推定合成信号値Ｗ＿ｈａｔが互いに一致するように、制御対象モデルのモデルパラメータを同定すればよいことになる。すなわち、前述した同定誤差ｅｉｄ’が値０になるように、同定値αｉｄを同定（算出）すればよいことになる。以上の理由により、前述した式（１７）〜（２９）の同定アルゴリズムを用いて、同定値αｉｄが算出される。
【０１２２】
また、制御対象モデルのモデルパラメータαと排ガスボリュームＶｅｘとが前述した図１０の関係を有している場合、最小２乗法などの一般的な同定アルゴリズムでは、図１０の関係を基準としてモデルパラメータを同定することができない。これに対して、本実施形態のオンボード同定器６０の場合、排ガスボリュームＶｅｘに応じて図１０のマップを検索することによって、基準モデルパラメータαｂｓを算出し、この基準モデルパラメータαｂｓを、前述した重み関数値Ｗａｉと修正項ｄαｉとの積の総和で修正することにより、同定値αｉｄが算出されるので、高い同定精度を確保することができる。
【０１２３】
次に、前述した周波数整形コントローラ１３０について説明する。この周波数整形コントローラ１３０は、予測当量比ＰＲＥ＿ＫＡＣＴが目標当量比ＫＣＭＤに収束するように、空燃比補正係数ＫＡＦを算出するもの、言い換えれば、検出当量比ＫＡＣＴが目標当量比ＫＣＭＤに収束するように、空燃比補正係数ＫＡＦを算出するものである。この周波数整形コントローラ１３０では、まず、下式（３４）に示すように、前述した予測当量比ＰＲＥ＿ＫＡＣＴから目標当量比ＫＣＭＤを減算することにより、予測追従誤差ＰＲＥ＿ｅが算出される。
【数３４】
0005400743.tif 000035
【０１２４】
次いで、下式（３５）により、制御入力としての空燃比補正係数ＫＡＦが算出される。
【数３５】
0005400743.tif 000036
この式（３５）において、βは感度設定パラメータであり、以下に述べる手法により所定値（例えば値０．６）に設定されている。
【０１２５】
次に、上記周波数整形コントローラ１３０の制御アルゴリズムの導出原理について説明する。本実施形態の場合、良好な排ガス特性と良好な燃費とを両立させることを目的として、ガソリンエンジンタイプのエンジン３において、その空燃比をリーン側に制御し、リーンバーン運転するとともに、尿素ＳＣＲ装置によって排ガス中のＮＯｘを浄化するように構成されている。
【０１２６】
このように構成した場合、ガソリンエンジンでは、リーンバーン運転中の燃焼安定性が低く、燃焼可能な混合気の空燃比が所定範囲内に制限されることになるので、空燃比が一時的に過度にリーン化する現象を抑制する必要がある。この現象は、過渡運転状条件下で特に発生しやすい。これに加えて、リーンバーン運転中は、燃焼変動に起因して、サージング現象が発生しやすくなるので、それを抑制するために、燃料量が過度に変動しないように制御する必要がある。これらの要求を満たすには、低周波数の外乱に対する抑制能力が低く、かつ高周波数の外乱に対する抑制能力が高くなるように、空燃比を制御することが必要となる。以下、これを「制御条件φ」という。
【０１２７】
ここで、本実施形態の制御装置１のようなフィードバック制御系の構成、すなわち制御入力としての空燃比補正係数ＫＡＦを制御対象に入力することにより、検出当量比ＫＡＣＴを目標当量比ＫＣＭＤに収束させるようにフィードバック制御する系の構成を、ｚ領域のブロック図で表すと、図１２のようになる。同図において、Ｃ（ｚ）はコントローラの伝達関数を、Ｐ（ｚ）は制御対象の伝達関数を、Ｄ（ｚ）は外乱をそれぞれ表している。なお、以下の説明では、各データにおける記号（ｚ）を適宜省略する。
【０１２８】
この制御系の場合、外乱Ｄと検出当量比ＫＡＣＴとの間の伝達関数すなわち感度関数Ｓは、下式（３６）のようになる。
【数３６】
0005400743.tif 000037
【０１２９】
この場合、上記制御条件φを満たすには、感度関数Ｓにおいて、そのゲイン特性（すなわち周波数応答特性）を表すゲイン曲線が図１３に示すものになることが必要になる。以下、制御条件φを満たす図１３のゲイン曲線が得られる感度関数を「最適感度関数Ｓｏｐｔ」という。同図において、ＦＱ１は、所定周波数を表しており、この所定周波数ＦＱ１は予め実験により設定される。同図に示すように、この最適感度関数Ｓｏｐｔの場合、所定周波数ＦＱ１以上の、外乱を抑制する必要性が高い周波数領域（以下「外乱抑制領域」という）では、高ゲインに設定されている一方、所定周波数ＦＱ１未満の、外乱を抑制する必要性が低い周波数域（以下「外乱非抑制領域」という）では、外乱抑制領域よりも低いゲインに設定されている。より具体的には、最適感度関数Ｓｏｐｔは、外乱抑制領域では、高ゲインかつフラットな状態のゲイン特性を示すとともに、外乱非抑制領域では、周波数が低いほど、ゲインが連続的に急減するような特性を示すように設定されている。なお、本実施形態では、制御条件φを満たすように設定された最適感度関数Ｓｏｐｔが、所定の周波数特性が得られるように設定された感度関数に相当する。
【０１３０】
ここで、前述した特許文献１におけるスライディングモード制御アルゴリズムを図１２の制御系に適用した場合、以下のようになる。すなわち、スライディングモード制御アルゴリズムの場合、追従誤差ｅおよび切換関数σが以下の式（３７），（３８）のように定義される。
【数３７】
0005400743.tif 000038
【数３８】
0005400743.tif 000039
ここで、ＰＯＬＥ＿Ｅは、−１＜ＰＯＬＥ＿Ｅ＜０が成立するように設定される切換関数設定パラメータである。
【０１３１】
スライディングモード制御アルゴリズムの場合、制御対象の動特性をσ＝０になるように拘束する制御手法であり、σ＝０を上式（３８）に適用すると、下式（３９）が得られるとともに、この式（３９）を整理すると、下式（４０）が得られる。
【数３９】
0005400743.tif 000040
【数４０】
0005400743.tif 000041
【０１３２】
上式（４０）は、入力の存在しない一次遅れ系を表している。すなわち、スライディングモード制御アルゴリズムは、制御対象の動特性を入力の存在しない一次遅れ系に拘束する制御アルゴリズムであり、そのような一次遅れ系の感度関数Ｓｓｌｄのゲイン曲線は、図１４に実線で示すものとなる。同図を参照すると明らかなように、この感度関数Ｓｓｌｄの場合、そのゲイン曲線は同図に破線で示す最適感度関数Ｓｏｐｔのゲイン曲線とかなり近似しており、前述した制御条件φを満たすものとなることが判る。
【０１３３】
ここで、スライディングモード制御アルゴリズムの場合、追従誤差ｅが切換直線上の値となる（すなわちσ＝０となる）までの到達モードと、追従誤差ｅが切換直線上の値となった（すなわち、制御対象の動特性が入力の存在しない一次遅れ系に拘束された）以降のスライディングモードとが存在している。そのため、スライディングモードにあるときには、制御条件φを満たすことができるものの、到達モードにあるときには、制御条件φを満たすことができないことになる。すなわち、スライディングモード制御アルゴリズムの場合、制御条件φを常に満たすことができない。
【０１３４】
そこで、本実施形態では、制御条件φを常に満たす制御アルゴリズムとして、以下に述べるように、感度関数Ｓｄを制御条件φを満たすように予め設定する制御アルゴリズムを採用した。まず、当量比の次元を備える空燃比補正係数ＫＡＦを制御入力とし、検出当量比ＫＡＣＴを制御量とする系を一次遅れ系と想定すると、その制御対象モデルは前述した式（１）となるとともに、その制御対象モデルのｚ領域の伝達関数Ｐは、下式（４１）となる。
【数４１】
0005400743.tif 000042
【０１３５】
一方、制御条件φを満たす感度関数Ｓｄを下式（４２）に示すように定義する。
【数４２】
0005400743.tif 000043
【０１３６】
上式（４２）のβは、感度関数設定パラメータであり、０＜β＜１となる所定値に設定される。上式（４２）の感度関数Ｓｄにおいて、β＝０．６と設定したときのゲイン曲線は、図１５に示す実線で示すものとなる。同図を参照すると明らかなように、この感度関数Ｓｄの場合、そのゲイン曲線は、同図に破線で示す最適感度関数Ｓｏｐｔのゲイン曲線とかなり近似しており、前述した制御条件φを満たすものとなることが判る。
【０１３７】
この感度関数Ｓｄと、コントローラの伝達関数Ｃおよび制御対象の伝達関数Ｐの関係は、下式（４３）に示すものとなる。
【数４３】
0005400743.tif 000044
【０１３８】
上式（４３）を変形し、感度関数Ｓｄの定義式をコントローラＣについて解くと、下式（４４）が得られる。
【数４４】
0005400743.tif 000045
【０１３９】
この式（４４）に式（４２）を代入すると、下式（４５）が得られる。
【数４５】
0005400743.tif 000046
【０１４０】
この式（４５）を離散時間系の漸化式で表すと、下式（４６）が得られる。
【数４６】
0005400743.tif 000047
【０１４１】
この式（４６）を参照すると、制御対象モデルのモデルパラメータαと、感度関数Ｓｄの周波数応答特性（ゲイン特性）を決定する感度関数設定パラメータβとによって、コントローラのフィードバックゲインを指定（設定）できることが判る。
【０１４２】
一方、前述した図１２の制御系の場合、相補感度関数Ｔは下式（４７）で表される。
【数４７】
0005400743.tif 000048
【０１４３】
ここで、相補感度関数Ｔと感度関数Ｓの関係は下式（４８）で表されることが知られている。
【数４８】
0005400743.tif 000049
【０１４４】
以上の式（４７），（４８）から明らかなように、前述した式（４６）を導出した手法は、外乱Ｄと検出当量比ＫＡＣＴとの間の周波数応答特性（ゲイン特性）を決定すると同時に、目標当量比ＫＣＭＤと検出当量比ＫＡＣＴとの間の周波数応答特性（ゲイン特性）を決定していることになる。
【０１４５】
ここで、図１５の感度関数Ｓｄに対応する相補感度関数をＴｄとした場合、この相補感度関数Ｔｄのゲイン曲線は、図１６に示すものとなる。なお、同図において、破線で示す曲線は、後述するモデル化誤差Δｌを相補感度関数Ｔｄに反映させた場合のゲイン曲線である。
【０１４６】
以上のように、式（４６）の制御アルゴリズムは、制御条件φを満たす感度関数Ｓｄを用いて導出されたものであるが、この式（４６）をそのまま用いて制御を実行しようとした場合、以下に述べる問題点１，２が生じてしまう。
【０１４７】
＜問題点１＞：制御対象モデルのモデルパラメータαの変動やばらつきに対応することができず、高いロバスト性を確保することができない。例えば、従来のＰＩＤ制御アルゴリズムや最適制御アルゴリズムと同様のロバスト性しか確保することができない。
＜問題点２＞：制御対象がむだ時間特性を備えている場合、それに対応することができず、制御精度の低下を招いてしまう。
【０１４８】
まず、＜問題点１＞について具体的に説明する。式（１）の制御対象モデルと、実際の制御対象とのモデル式誤差をΔｌ（ｚ）とした場合、制御系の安定条件として、下式（４９）の成立が必要であることが知られている。
【数４９】
0005400743.tif 000050
【０１４９】
ここで、式（１）の一次遅れ系モデルのような遅れ系モデルの場合、そのモデル化誤差Δｌは、図１７に示すように、周波数域が高くなるほど、より増大するという特性を備えているので、このモデル化誤差Δｌを前述した相補感度関数Ｔｄに反映させると、図１６の破線で示すゲイン曲線が得られる。上記式（４９）から明らかなように、値Ｔｄ・Δｌが０ｄＢ未満であることが制御系の安定条件となるので、相補感度関数Ｔｄのゲインが０ｄＢを下回っている度合が制御系の安定余裕となり、ロバスト性を表すことになる。
【０１５０】
しかし、感度関数Ｓｄと相補感度関数Ｔｄとの間には、前述したように、Ｔｄ（ｚ）＋Ｓｄ（ｚ）＝１の関係が存在しているので、外乱抑制に対する周波数応答特性とロバスト性とを互いに独立して設定することはできない。したがって、外乱抑制に対する周波数応答特性を指定した状態で、遅れ系モデルのモデル化誤差Δｌ（ｚ）に対するロバスト性を向上させるには、このモデル化誤差Δｌ（ｚ）を補償できるような他の制御アルゴリズムが必要となる。
【０１５１】
なお、モデル化誤差Δｌ（ｚ）に対応するために、式（４２）の次数を増やし、感度関数Ｓｄを複雑な形に変更した場合、式（４５）の伝達関数Ｃ（ｚ）において、分母のｚの次数よりも分子のｚの次数が大きくなってしまい、実現不可能なコントローラとなってしまう。また、感度設定パラメータβをトライ＆エラーでチューニングする手法を用いた場合、ＰＩＤ制御のゲインや、最適制御の重み関数Ｑ，Ｒをチューニングする手法と変わらず、外乱抑制の周波数応答特性を直接指定する前述した式（４６）による制御手法の利点を失ってしまう。
【０１５２】
次に、前述した＜問題点２＞について説明する。本実施形態の制御系の場合、空燃比補正係数ＫＡＦと検出当量比ＫＡＣＴとの間にはむだ時間ｄが存在し、その制御対象モデルとして前述した式（２）を用いているが、前述した式（２）の制御対象のｚ領域での伝達関数Ｐ（ｚ）は、下式（５０）となる。
【数５０】
0005400743.tif 000051
【０１５３】
この式（５０）の伝達関数Ｐ（ｚ）においてｄ＝２に設定したときのボード線図は、図１８に示すものとなり、前述した式（４１）のむだ時間のない制御系の伝達関数Ｐ（ｚ）のボード線図は、図１９に示すものとなる。両図を比較すると明らかなように、むだ時間ｄの有無は、ゲイン特性の差異として現れないので、むだ時間を前述したモデル化誤差Δｌとして表すことはできず、前述した式（４６）を用いる制御手法、すなわち感度関数Ｓｄや相補感度関数Ｔｄのゲイン特性によってコントローラのゲインを指定する制御手法では、むだ時間に対するロバスト性を考慮し、これを補償することができない。
【０１５４】
これに対して、むだ時間が制御系に存在する場合、制御系の安定性が著しく低下することは周知であり、そのため、前述した制御手法をむだ時間を有する制御系に適用すると、制御系が発散してしまうおそれがある。
【０１５５】
また、前述した式（４４）に、前述した式（４２）および式（５０）を代入し、コントローラの伝達関数Ｃ（ｚ）を導出すると、下式（５１）が得られる。
【数５１】
0005400743.tif 000052
【０１５６】
この式（５１）を離散時間系の漸化式で表すと、下式（５２）が得られる。
【数５２】
0005400743.tif 000053
【０１５７】
この式（５２）の場合、追従誤差ｅの未来値ｅ（ｋ＋ｄ），ｅ（ｋ＋ｄ−１）が右辺に含まれているので、コントローラの制御アルゴリズムとして、実現不可能なものとなってしまう。
【０１５８】
さらに、本実施形態の制御対象の場合、制御入力としての空燃比補正係数ＫＡＦと制御量としての検出当量比ＫＡＣＴとの間のむだ時間ｄは、前述した図２に示すように、排ガスボリュームＶｅｘに応じて連続的に変化する特性を有しているので、前述した外乱抑制の周波数応答特性を直接的に指定する制御手法の場合、むだ時間を有する制御対象に適用不可能なものである以上、むだ時間ｄが変化するような制御系には当然のごとく適用不可能である。
【０１５９】
以上のように、前述した問題点１，２を解消すべく、前述した感度関数Ｓｄまたは相補感度関数Ｔｄを用いたコントローラすなわち前述した外乱抑制の周波数応答特性を直接的に指定する制御アルゴリズムを用いながら、制御対象モデルのモデルパラメータαの変動やばらつきに対応し、かつ、制御対象のむだ時間ｄが変化する特性に対して対応するできるように、制御アルゴリズムを構築する必要がある。
【０１６０】
この要求を実現するために、本実施形態の制御装置１の場合、まず、オンボード同定器６０で、前述した同定アルゴリズムにより、モデルパラメータαの同定値αｉｄを算出し、状態予測器４０において、前述した予測アルゴリズムにより、むだ時間ｄが経過したタイミングでの検出当量比ＫＡＣＴに相当する予測当量比ＰＲＥ＿ＫＡＣＴを算出する。
【０１６１】
そして、周波数整形コントローラ１３０の制御アルゴリズムとして、検出当量比ＫＡＣＴに代えて予測当量比ＰＲＥ＿ＫＡＣＴを用い、さらに、前述した式（２）のモデルパラメータαを同定値αｉｄに置き換えた下式（５３）を制御対象モデルとして用い、前述した式（４６）を導出した手法と同じ手法によって、前述した式（３４），（３５）が導出される。
【数５３】
0005400743.tif 000054
【０１６２】
この式（５３）は、前述した式（１）のαをαｉｄに置き換えたものである。言い換えれば、前述した制御対象モデルの式（２）のむだ時間特性を除去したもの（むだ時間特性を考慮しないもの）に相当する。
【０１６３】
次に、図２０を参照しながら、ＥＣＵ２によって実行される空燃比制御処理について説明する。この制御処理は、以下に述べるように、燃料噴射弁４から噴射すべき燃料噴射量Ｔｏｕｔを算出するものであり、前述した所定の制御周期ΔＴで実行される。
【０１６４】
この処理では、まず、ステップ１（図では「Ｓ１」と略す。以下同じ）で、エンジン回転数ＮＥおよび吸気圧ＰＢに応じて、図示しないマップを検索することにより、基本噴射量ＴｉＢＳを算出する。
【０１６５】
次いで、ステップ２に進み、ＬＡＦセンサ正常フラグＦ＿ＬＡＦＯＫが「１」であるか否かを判別する。このＬＡＦセンサ正常フラグＦ＿ＬＡＦＯＫは、図示しない判定処理において、ＬＡＦセンサ２３が正常であるときに「１」に、それ以外のときに「０」に設定される。
【０１６６】
ステップ２の判別結果がＮＯで、ＬＡＦセンサ２３が故障しているときには、ステップ１４に進み、燃料噴射量Ｔｏｕｔを基本噴射量ＴｉＢＳに設定した後、本処理を終了する。
【０１６７】
一方、ステップ２の判別結果がＹＥＳで、ＬＡＦセンサ２３が正常であるときには、ステップ３に進み、三元触媒活性化フラグＦ＿ＴＷＣＡＣＴが「１」であるか否かを判別する。この三元触媒活性化フラグＦ＿ＴＷＣＡＣＴは、図示しない判定処理において、２つの三元触媒１１，１２がいずれも活性化しているときに「１」に、それ以外のときに「０」に設定される。
【０１６８】
ステップ３の判別結果がＮＯで、２つの三元触媒１１，１２の少なくとも一方が活性化していないときには、ステップ７に進み、目標当量比ＫＣＭＤを所定のリーン用値ＫＬＥＡＲＮに設定する。この所定のリーン用値ＫＬＥＡＲＮは、エンジン始動直後におけるＨＣの発生を抑制できるような値（例えば０．９）に設定されている。
【０１６９】
一方、ステップ３の判別結果がＹＥＳで、２つの三元触媒１１，１２がいずれも活性化しているときには、ステップ４に進み、ＳＣＲ活性化フラグＦ＿ＳＣＲＡＣＴが「１」であるか否かを判別する。このＳＣＲ活性化フラグＦ＿ＳＣＲＡＣＴは、図示しない判定処理において、２つの選択還元触媒１４，１５の少なくとも一方が活性化しているときに「１」に、それ以外のときに「０」に設定される。
【０１７０】
ステップ４の判別結果がＮＯで、２つの選択還元触媒１４，１５の少なくとも一方が活性化していないときには、ステップ８に進み、目標当量比ＫＣＭＤを所定の理論空燃比用値ＫＳＴＯＩＣに設定する。この理論空燃比用値ＫＳＴＯＩＣは、理論空燃比に相当する当量比の値（＝１）に設定されている。
【０１７１】
一方、ステップ４の判別結果がＹＥＳで、２つの選択還元触媒１４，１５の少なくとも一方が活性化しているときには、ステップ５に進み、エンジン回転数ＮＥおよびアクセル開度ＡＰに応じて、図示しないマップを検索することにより、要求トルクＴＲＱＤＲＶを算出する。
【０１７２】
次いで、ステップ６に進み、エンジン回転数ＮＥおよび要求トルクＴＲＱＤＲＶに応じて、前述した図４のマップを検索することにより、目標当量比ＫＣＭＤを算出する。
【０１７３】
以上のステップ６〜８のいずれかに続くステップ９で、ＬＡＦセンサ活性化フラグＦ＿ＬＡＦＡＣＴが「１」であるか否かを判別する。このＬＡＦセンサ活性化フラグＦ＿ＬＡＦＡＣＴは、図示しない判定処理において、ＬＡＦセンサ２３が活性化しているときに「１」に、それ以外のときに「０」に設定される。
【０１７４】
ステップ９の判別結果がＮＯで、ＬＡＦセンサ２３が活性化していないときには、ステップ１３に進み、燃料噴射量Ｔｏｕｔを目標当量比と基本噴射量の積ＫＣＭＤ・ＴｉＢＳに設定する。その後、本処理を終了する。
【０１７５】
一方、ステップ９の判別結果がＹＥＳで、ＬＡＦセンサ２３が活性化しているときには、ステップ１０に進み、エンジン回転数ＮＥおよび吸気圧ＰＢに応じて、図示しないマップを検索することにより、排ガスボリュームＶｅｘを算出する。
【０１７６】
次いで、ステップ１１に進み、前述した制御アルゴリズムにより、空燃比補正係数ＫＡＦを算出する。具体的には、まず、前述した式（３）〜（７）の予測アルゴリズムと、図６から検索された重み関数値Ｗｄｉとを用いて、予測当量比ＰＲＥ＿ＫＡＣＴを算出する。さらに、前述した式（１７）〜（２９）の同定アルゴリズムと、図１０，１１からそれぞれ検索された基準モデルパラメータαｂｓおよび重み関数値Ｗａｉとを用いて、同定値αｉｄを算出する。そして、算出した予測当量比ＰＲＥ＿ＫＡＣＴおよび同定値αｉｄを用い、前述した式（３４），（３５）により、空燃比補正係数ＫＡＦが最終的に算出される。
【０１７７】
ステップ１１に続くステップ１２で、燃料噴射量Ｔｏｕｔを空燃比補正係数と基本噴射量の積ＫＡＦ・ＴｉＢＳに設定する。その後、本処理を終了する。
【０１７８】
本実施形態の制御装置１は、以上のような空燃比制御処理によって燃料噴射量Ｔｏｕｔを算出し、図示しないが、この燃料噴射量Ｔｏｕｔおよびエンジン回転数ＮＥに応じて、燃料噴射タイミングを算出するとともに、これらの燃料噴射量Ｔｏｕｔおよび燃料噴射タイミングに基づいた制御入力信号で、燃料噴射弁４を駆動することにより、混合気の空燃比を制御する。
【０１７９】
次に、図２１〜２７を参照しながら、本実施形態の制御装置１による空燃比制御のシミュレーション結果（以下「制御結果」という）について説明する。まず、図２１〜２３について説明する。これらの図２１〜図２３はいずれも、シミュレーション条件を、式（２）の制御対象モデルにおいてモデル化誤差が存在しないと設定した場合（具体的には、α＝αｂｓと設定した場合）のものであり、図２１は本実施形態の制御装置１による制御結果を示している。
【０１８０】
また、図２２は、比較のために、シミュレーション条件として、制御装置１における状態予測器４０およびオンボード同定器６０の演算を省略した場合、具体的にはＰＲＥ＿ＫＡＣＴ（ｋ）＝ＫＡＣＴ（ｋ），αｉｄ（ｋ）＝αｂｓ（ｋ）と設定した場合の制御結果（以下「比較例１」という）を示している。さらに、図２３は、比較のために、制御装置１における状態予測器４０およびオンボード同定器６０の演算を省略するとともに、感度設定パラメータβとして、実施形態の設定値の１／６の値を用いた場合の制御結果（以下「比較例２」という）を示している。
【０１８１】
まず、図２２の比較例１を参照すると、排ガスボリュームＶｅｘが小さいときには、予測当量比ＰＲＥ＿ＫＡＣＴすなわち検出当量比ＫＡＣＴが発散し、それに従って、予測追従誤差ＰＲＥ＿ｅや空燃比補正係数ＫＡＦも発散した状態となっていることが判る。すなわち、むだ時間が存在する制御対象を、周波数整形コントローラ１３０のみを用いて制御した場合、相補感度関数Ｔｄによって指定したロバスト性を適切に維持することができず、特に、むだ時間ｄが増大するような、排ガスボリュームＶｅｘが小さい条件下では、制御入力としての空燃比補正係数ＫＡＦが発散してしまうことが判る。
【０１８２】
次に、図２３の比較例２を参照すると、この比較例２の場合、上述した比較例１と比べて、制御系の安定性および制御精度が向上していることが判る。これは、感度関数Ｓｄにおける感度設定パラメータβを比較例１の１／６の値に設定し、フィードバックゲインを低下させたことによるものであり、言い換えれば、外乱抑制能力を低下させたことに起因するものである。この場合の感度設定パラメータβの値は、トライ＆エラー手法によって、制御系の安定性を維持できる限界値に設定されている。したがって、前述した制御条件φを満たすように、感度関数Ｓｄを設定することで、外乱抑制能力を直接的に指定するという本発明の目的を実現できないものとなっている。
【０１８３】
これに対して、図２１の本実施形態の制御結果の場合、状態予測器４０、オンボード同定器６０および周波数整形コントローラ１３０の各種アルゴリズムによって、モデル化誤差が存在しない条件下において、制御の安定性およ制御精度が比較例１，２よりも向上していることが判る。
【０１８４】
次に、図２４〜２７について説明する。これらの図２４〜２７はいずれも、シミュレーション条件を、式（２）の制御対象モデルにおいてモデル化誤差が存在するように設定した場合（具体的には、α＝２・αｂｓとなるように設定した場合）のものであり、図２４は本実施形態の制御装置１による制御結果を示している。
【０１８５】
また、図２５は、比較のために、シミュレーション条件として、制御装置１における状態予測器４０およびオンボード同定器６０の演算を省略した場合の制御結果（以下「比較例３」という）を示している。さらに、図２６は、比較のために、制御装置１における状態予測器４０およびオンボード同定器６０の演算を省略するとともに、感度設定パラメータβとして、実施形態の設定値の１／６の値を用いた場合の制御結果（以下「比較例４」という）を示している。これに加えて、図２７は、比較のために、制御装置１におけるオンボード同定器６０の演算のみを省略した場合、すなわちαｉｄ（ｋ）＝αｂｓ（ｋ）と設定した場合の制御結果（以下「比較例５」という）を示している。
【０１８６】
まず、図２５の比較例３を参照すると、この比較例３の場合、制御対象モデルにおいてモデル化誤差が存在するシミュレーション条件下では、むだ時間に起因する制御系の安定余裕の低下に加えて、モデル化誤差の影響によって、制御系の安定性が損なわれ、排ガスボリュームＶｅｘの全領域に対して、空燃比補正係数ＫＡＦなどの各パラメータがすべて発散した状態となっていることが判る。すなわち、むだ時間が存在する制御対象を周波数整形コントローラ１３０のみを用いて制御した場合、モデル化誤差が存在するシミュレーション条件下では、制御の安定性および制御精度が極めて低下してしまうことが判る。
【０１８７】
次に、図２６の比較例４を参照すると、この比較例４の場合、上述した比較例３のような発散状態は発生しておらず、比較例３と比べて、制御系の安定性が向上していることが判る。これは、前述したように、感度設定パラメータβの設定値に起因するものである。しかし、この比較例４の場合、比較例３と比べて、制御系の安定性が向上しているものの、予測追従誤差ＰＲＥ＿ｅが一時的に過大な値となる状態が発生しており、制御精度が低い状態にあることが判る。これに加えて、前述したように、感度設定パラメータβを実施形態の設定値の１／６の値に設定している以上、制御条件φを満たすように、感度関数Ｓｄを設定することで、外乱抑制能力を直接的に指定するという本発明の目的を実現することができないものとなっている。
【０１８８】
また、図２７の比較例５を参照すると、この比較例５の場合、上述した比較例３と比べて、制御系の安定性が向上していることが判る。これは、排ガスボリュームＶｅｘの変化に伴ってむだ時間ｄが連続的に変化したときでも、状態予測器４０において、そのようなむだ時間ｄの変化を反映させながら予測当量比ＰＲＥ＿ＫＡＣＴが算出されるので、むだ時間ｄの変化の影響を適切に補償でき、制御系の安定性が向上することによる。しかし、この比較例５の場合でも、排ガスボリュームＶｅｘが大きい領域では、モデル化誤差が増大するのに起因して、空燃比補正係数ＫＡＦなどの各パラメータが発散した状態となっていることが判る。
【０１８９】
一方、図２４の本実施形態の制御結果の場合、モデル化誤差が存在するシミュレーション条件下でも、状態予測器４０、オンボード同定器６０および周波数整形コントローラ１３０の各種アルゴリズムによって、制御の安定性および制御精度が比較例３〜５よりも向上していることが判る。例えば、状態予測器４０の予測アルゴリズムによって、予測追従誤差ＰＲＥ＿ｅが小さい値に保持されているとともに、オンボード同定器６０の同定アルゴリズムによって、時間の経過に伴い、同定値αｉｄがモデルパラメータαに収束していることが判る。
【０１９０】
以上のように、第１実施形態の制御装置１によれば、状態予測器４０において、検出当量比ＫＡＣＴと空燃比補正係数ＫＡＦの関係を定義した制御対象モデル（式（２））を用い、むだ時間ｄ＝０〜３がそれぞれ経過したタイミングでの検出当量比ＫＡＣＴとして、第０〜第３予測値ＰＲＥ＿ＫＡＣＴ＿０〜３が算出され、排ガスボリュームＶｅｘに応じて、４つの重み関数値Ｗｄ１〜Ｗｄ４が算出される。そして、予測当量比ＰＲＥ＿ＫＡＣＴが、重み関数値Ｗｄｉと予測値ＰＲＥ＿ＫＡＣＴ＿４−ｉ（ｉ＝１〜４）との乗算値の総和として算出されるので、予測当量比ＰＲＥ＿ＫＡＣＴを、予測値ＰＲＥ＿ＫＡＣＴ＿４−ｉを連続的に結合させた値として算出することができる。それにより、排ガスボリュームＶｅｘの変化に応じてむだ時間ｄが変化した場合でも、予測検出当量比ＫＡＣＴを、そのようなむだ時間ｄの変化を適切に補償しながら精度よく算出することができる。特に、排ガスボリュームＶｅｘの急変に伴ってむだ時間ｄが急変したときでも、予測検出当量比ＫＡＣＴを、むだ時間ｄの急変を適切に補償しながら、段差なく円滑に変化するように算出することができる。
【０１９１】
また、オンボード同定器６０において、前述した同定アルゴリズムにより、同定値αｉｄが算出されるので、前述した同定条件１，２の双方を満たしながら、同定値αｉｄを算出することができる。具体的には、合成信号値Ｗ＿ａｃｔおよび推定合成信号値Ｗ＿ｈａｔが互いに一致するように、同定値αｉｄが算出されるので、同定条件１すなわち拘束条件を満たしながら、同定値αｉｄを算出することができる。さらに、修正後制御入力ＫＡＦ＿ｍｏｄが、４つの重み関数値Ｗｄ４〜１を、むだ時間ｄ＝０〜３前のタイミングでの空燃比補正係数ＫＡＦ（ｋ），ＫＡＦ（ｋ−１），ＫＡＦ（ｋ−２），ＫＡＦ（ｋ−３）にそれぞれ乗算した乗算値の総和として算出されるので、排ガスボリュームＶｅｘの変化に応じてむだ時間ｄが連続的に変化する場合でも、修正後制御入力ＫＡＦ＿ｍｏｄを、そのようなむだ時間ｄの変化を適切に補償しながら精度よく算出することができる。特に、排ガスボリュームＶｅｘの急変に伴ってむだ時間ｄが急変したときでも、修正後制御入力ＫＡＦ＿ｍｏｄを、むだ時間ｄの急変を適切に補償しながら、段差なく円滑に変化するように算出することができる。
【０１９２】
さらに、そのような修正後制御入力ＫＡＦ＿ｍｏｄと検出当量比ＫＡＣＴとの関係を定義した式（３０）のモデルを用いて導出した式（１７）〜（２９）の同定アルゴリズムにより、同定値αｉｄがオンボード同定される。すなわち、２種類の重み関数値Ｗｄｉ，Ｗａｉを用いながら、同定値αｉｄが算出されるので、排ガスボリュームＶｅｘの変化に応じてむだ時間ｄや遅れ特性が変化したときでも、その影響を抑制しながら、同定値αｉｄを精度よく同定することができる。特に、排ガスボリュームＶｅｘの急変に伴ってむだ時間ｄや遅れ特性が急変したときでも、同定値αｉｄを、むだ時間ｄや遅れ特性の急変を適切に補償しながら、段差なく円滑に変化するように算出することができる。そして、そのように算出された同定値αｉｄを用いて、制御入力としての空燃比補正係数ＫＡＦが算出されるので、空燃比制御における制御性や、個体間のばらつきおよび経年変化の影響に対する空燃比制御のロバスト性を飛躍的に向上させることができる。
【０１９３】
これに加えて、周波数整形コントローラ１３０の場合、前述したように、所定の周波数特性が得られるように設定された感度関数Ｓｄに基づいて導出された算出式（３４），（３５）を用いて、空燃比補正係数ＫＡＦが算出されるので、前述した制御条件φを満たしながら、空燃比補正係数ＫＡＦを算出することができる。これに加えて、制御対象モデルのモデルパラメータとして前述した同定値αｉｄを用いているので、むだ時間ｄや遅れ特性の変化を適切に補償しながら、制御装置１における外乱抑制特性やロバスト性を周波数軸上でダイレクトに指定（設定）することができる。それにより、外乱による検出当量比ＫＡＣＴの変動を抑制したい周波数域において、外乱抑制能力やロバスト性を飛躍的に向上させることができる。これに加えて、空燃比補正係数ＫＡＦの算出アルゴリズムにおいて、予測当量比ＰＲＥ＿ＫＡＣＴと目標当量比ＫＣＭＤとの偏差に基づいたフィードバック制御アルゴリズムが用いられているので、むだ時間ｄを補償できることによって、フィードバックゲインを高い状態に維持することができ、高精度かつ高応答性を確保しながら、検出当量比ＫＡＣＴを目標当量比ＫＣＭＤに追従させることができる。
【０１９４】
なお、第１実施形態は、重み関数値として、互いに重なり合う領域における排ガスボリュームＶｅｘのいずれかの値に対応する重み関数値Ｗｄｉの値の和が、各重み関数値Ｗｄｉにおける最大値１と等しくなるように設定されているものを用いた例であるが、本発明の重み関数値はこれに限らず、互いに重なり合う領域における参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が所定値になるように設定されているものであればよい。例えば、互いに重なり合う領域における参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が、重み関数値における絶対値の最大値と等しくなるように設定されたものを用いてもよい。より具体的には、重み関数値として、図６の重み関数値Ｗｄｉの設定値に対して、ｘ軸を対称軸として線対称に設定されたもの、すなわち負値に設定されたものを用いてもよく、その場合には、４つの重み関数値に乗算する値、すなわち４つの予測値ＰＲＥ＿ＫＡＣＴ＿＿４−ｉや、４つの空燃比補正係数ＫＡＦ（ｋ−４＋ｉ）として、負値化した値を用いればよい。
【０１９５】
また、第１実施形態のオンボード同定器６０において、前述した式（１７）〜（２９）の同定アルゴリズムに代えて、以下の式（５４）〜（６７）の同定アルゴリズムを用いて、同定値αｉｄを算出するように構成してもよい。
【０１９６】
【数５４】
0005400743.tif 000055
【数５５】
0005400743.tif 000056
【数５６】
0005400743.tif 000057
【数５７】
0005400743.tif 000058
【数５８】
0005400743.tif 000059
【数５９】
0005400743.tif 000060
【数６０】
0005400743.tif 000061
【数６１】
0005400743.tif 000062
【数６２】
0005400743.tif 000063
【数６３】
0005400743.tif 000064
【数６４】
0005400743.tif 000065
【数６５】
0005400743.tif 000066
【数６６】
0005400743.tif 000067
【数６７】
0005400743.tif 000068
【０１９７】
以上の式（５４）〜（６７）を前述した式（１７）〜（２９）と比較すると明らかなように、式（５４）〜（６４）は式（１７）〜（２７）と同一であり、式（６５）〜（６７）のみが異なっているので、以下、これらの式（６５）〜（６７）についてのみ説明する。まず、上式（６５）を用いて、基準モデルパラメータαｂｓ’が算出される。すなわち、基準モデルパラメータαｂｓ’は、前述した基準モデルパラメータαｂｓを補正係数Ｋαｂｓで補正することによって算出される。
【０１９８】
この補正係数Ｋαｂｓは、エンジン回転数ＮＥおよび検出当量比ＫＡＣＴに応じて、図２８に示すマップを検索することにより算出される。同図において、３つの値ＫＡＣＴ＿Ｒ，ＫＡＣＴ＿Ｓ，ＫＡＣＴ＿Ｌはいずれも検出当量比ＫＡＣＴの所定値であり、ＫＡＣＴ＿Ｓ＝１，ＫＡＣＴ＿Ｌ＜ＫＡＣＴ＿Ｓ＜ＫＡＣＴ＿Ｒが成立するように設定されている。
【０１９９】
このマップでは、補正係数Ｋαｂｓは、値１以下の値に設定されているとともに、エンジン回転数ＮＥが低いほど、より小さい値に設定されている。これは、低回転域では、排ガスボリュームＶｅｘが同じでも、１燃焼サイクルの実行時間が長くなるのに伴い、排ガス成分の周期的変動が大きくなることによって、空燃比補正係数ＫＡＦと検出当量比ＫＡＣＴとの間の応答遅れが大きくなるので、それに対応するためである。
【０２００】
また、補正係数Ｋαｂｓは、検出当量比ＫＡＣＴがリッチ側の値であるほど、より大きい値に設定されている。これは、検出当量比ＫＡＣＴが大きく、排ガス濃度が高い場合には、排ガス中の未燃成分が多くなり、ＬＡＦセンサ２３の検出素子の応答性が高くなることによって、空燃比補正係数ＫＡＦと検出当量比ＫＡＣＴとの間の応答遅れが小さくなるので、それに対応するためである。
【０２０１】
次に、前述した式（６６）を用いて、修正項ｄαｉｊｈ’を算出した後、前述した式（６７）により、同定値αｉｄが最終的に算出される。同式（６６）において、Ｗａｎｊ，Ｗａａｈは重み関数値であり、重み関数値Ｗａｎｊ（ｊ＝１〜４）は、エンジン回転数ＮＥに応じて、図２９に示すマップを検索することにより算出される。同図において、ＮＥ１〜４，ＮＥＭＡＸは、エンジン回転数ＮＥの所定値であり、０＜ＮＥ１＜ＮＥ２＜ＮＥ３＜ＮＥ４＜ＮＥＭＡＸが成立するように設定されている。また、所定値ＮＥＭＡＸは、最大許容回転数に設定されている。
【０２０２】
同図２９に示すように、４つの重み関数値Ｗａｎ１〜４はそれぞれ、エンジン回転数ＮＥが変化し得る領域を０≦ＮＥ≦ＮＥ２，ＮＥ１≦ＮＥ≦ＮＥ３，ＮＥ２≦ＮＥ≦ＮＥ４，ＮＥ３≦ＮＥ≦ＮＥＭＡＸの４つの領域に区分した場合において、これらの４つの領域に対応するように設定されているとともに、対応する領域では、値１以下の正の値に設定され、対応する領域以外では、値０に設定されている。
【０２０３】
具体的には、重み関数値Ｗａｎ１は、これが対応する領域（０≦ＮＥ≦ＮＥ２）では、ＮＥ≦ＮＥ１のときの値１を最大値として、エンジン回転数ＮＥが大きいほど、より小さい正の値に設定されているとともに、それ以外の領域では、値０に設定されている。また、重み関数値Ｗａｎ２は、これが対応する領域（ＮＥ１≦ＮＥ≦ＮＥ３）では、ＮＥ＝ＮＥ２のときの値１を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値０に設定されている。
【０２０４】
さらに、重み関数値Ｗａｎ３は、これが対応する領域（ＮＥ２≦ＮＥ≦ＮＥ４）では、ＮＥ＝ＮＥ３のときの値１を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値０に設定されている。一方、重み関数値Ｗａｎ４は、これが対応する領域（ＮＥ３≦ＮＥ≦ＮＥＭＡＸ）では、ＮＥ４≦ＮＥのときの値１を最大値として、エンジン回転数ＮＥが大きいほど、より大きい正の値に設定されているとともに、それ以外の領域では、値０に設定されている。
【０２０５】
以上に加えて、４つの重み関数値Ｗａｎｊ（ｊ＝１〜４）がそれぞれ対応する４つの領域は、上述したような、隣り合う領域が互いに重なり合うように設定されており、これらの互いに重なり合う領域におけるエンジン回転数ＮＥのいずれかの値に対応する重み関数値Ｗａｎｊの値の和は、各重み関数値Ｗａｎｉにおける最大値１と等しくなるように設定されている。以上のようにエンジン回転数ＮＥに応じて算出される重み関数値Ｗａｎｊを用いた理由は、補正係数Ｋαｂｓの算出で述べた理由とおなじである。
【０２０６】
また、前述した式（６６）の重み関数値Ｗａａｈ（ｈ＝１〜４）は、検出当量比ＫＡＣＴに応じて、図３０に示すマップを検索することにより、４つの重み関数値Ｗａａ１〜４がそれぞれ算出される。同図において、ＫＡＣＴ１〜４，ＫＡＣＴＭＡＸは、検出当量比ＫＡＣＴの所定値であり、０＜ＫＡＣＴ１＜ＫＡＣＴ２＜ＫＡＣＴ３＜ＫＡＣＴ４＜ＫＡＣＴＭＡＸが成立するように設定されている。さらに、所定値ＫＡＣＴＭＡＸは、エンジン３の運転中に検出当量比ＫＡＣＴが変化し得る領域の最大値に設定されている。言い換えれば、検出当量比ＫＡＣＴは、エンジン３の運転中、０〜ＫＡＣＴＭＡＸの領域内で変化する特性を有している。
【０２０７】
同図３０に示すように、４つの重み関数値Ｗａａ１〜４はそれぞれ、検出当量比ＫＡＣＴが変化し得る領域をＫＡＣＴ≦ＫＡＣＴ２，ＫＡＣＴ１≦ＫＡＣＴ≦ＫＡＣＴ３，ＫＡＣＴ２≦ＫＡＣＴ≦ＫＡＣＴ４，ＫＡＣＴ３≦ＫＡＣＴ≦ＫＡＣＴＭＡＸの４つの領域に区分した場合において、これらの４つの領域に対応するように設定されているとともに、対応する領域では、値１以下の正の値に設定され、対応する領域以外では、値０に設定されている。
【０２０８】
具体的には、重み関数値Ｗａａ１は、これが対応する領域（ＫＡＣＴ≦ＫＡＣＴ２）では、ＫＡＣＴ≦ＫＡＣＴ１のときの値１を最大値として、検出当量比ＫＡＣＴが大きいほど、より小さい正の値に設定されているとともに、それ以外の領域では、値０に設定されている。また、重み関数値Ｗａａ２は、これが対応する領域（ＫＡＣＴ１≦ＫＡＣＴ≦ＫＡＣＴ３）では、ＫＡＣＴ＝ＫＡＣＴ２のときの値１を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値０に設定されている。
【０２０９】
さらに、重み関数値Ｗａａ３は、これが対応する領域（ＫＡＣＴ２≦ＫＡＣＴ≦ＫＡＣＴ４）では、ＫＡＣＴ＝ＫＡＣＴ３のときの値１を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値０に設定されている。一方、重み関数値Ｗａａ４は、これが対応する領域（ＫＡＣＴ３≦ＫＡＣＴ≦ＫＡＣＴＭＡＸ）では、ＫＡＣＴ４≦ＫＡＣＴのときの値１を最大値として、検出当量比ＫＡＣＴが大きいほど、より大きい正の値に設定されているとともに、それ以外の領域では、値０に設定されている。
【０２１０】
以上に加えて、４つの重み関数値Ｗａａｈ（ｈ＝１〜４）がそれぞれ対応する４つの領域は、上述したような、隣り合う領域が互いに重なり合うように設定されており、これらの互いに重なり合う領域における検出当量比ＫＡＣＴのいずれかの値に対応する重み関数値Ｗａａｈの値の和は、各重み関数値Ｗａａｈにおける最大値１と等しくなるように設定されている。以上のように、検出当量比ＫＡＣＴに応じて算出される重み関数値Ｗａａｈを用いた理由は、補正係数Ｋαｂｓの算出で述べた理由とおなじである。
【０２１１】
以上の同定アルゴリズムを用いて同定値αｉｄを算出した場合、排ガスボリュームＶｅｘに加えて、エンジン回転数ＮＥおよび検出当量比ＫＡＣＴの変化に伴うむだ時間ｄや遅れ特性の変化を反映させながら、同定値αｉｄを算出することができる。すなわち、３つのパラメータＶｅｘ，ＮＥ，ＫＡＣＴの変化に起因するむだ時間ｄや遅れ特性の変化を補償しながら、同定値αｉｄを算出することができ、同定値αｉｄの同定精度（すなわち算出精度）をさらに向上させることができる。それにより、空燃比制御の制御性およびロバスト性を、第１実施形態のオンボード同定器６０を用いた場合よりも向上させることができる。
【０２１２】
次に、図３１を参照しながら、本発明の第２実施形態に係る制御装置１Ａについて説明する。この制御装置１Ａは、前述した制御装置１と同様に、空燃比補正係数ＫＡＦなどを算出することによって、空燃比を制御するものである。第２実施形態の場合、制御対象モデルとして、第１実施形態で用いた式（２）に代えて、下式（６８）を用いる。
【数６８】
0005400743.tif 000069
【０２１３】
上式（６８）のδはモデルパラメータである。この式（６８）は、式（２）の「１−α」を「δ」に置き換えたものであり、それにより、前述した２つのモデルパラメータ１−α，α間の拘束条件を取り除いたものに相当する。
【０２１４】
図３１に示すように、制御装置１Ａは、目標当量比算出部２３０、可変むだ時間状態予測器（以下「状態予測器」という）２４０、オンボード・スケジュールド・モデルパラメータ同定器（以下「オンボード同定器」という）２６０、周波数整形コントローラ３３０を備えており、これらはいずれもＥＣＵ２によって構成されている。
【０２１５】
この目標当量比算出部２３０では、前述した目標当量比算出部３０と同じ手法によって、目標当量比ＫＣＭＤが算出される。また、状態予測器２４０では、後述する予測アルゴリズムを用いて、予測当量比ＰＲＥ＿ＫＡＣＴが算出され、オンボード同定器２６０では、後述する同定アルゴリズムを用いて、２つのモデルパラメータδ，αを要素とするモデルパラメータベクトルθが算出される。さらに、周波数整形コントローラ３３０では、後述する制御アルゴリズムを用いて、制御入力としての空燃比補正係数ＫＡＦが算出される。
【０２１６】
なお、本実施形態では、目標当量比算出部２３０が目標制御量算出手段に相当し、目標当量比ＫＣＭＤが目標制御量に相当する。また、状態予測器２４０が、予測値算出手段、重み関数値算出手段および予測制御量設定手段に相当し、予測当量比ＰＲＥ＿ＫＡＣＴが予測制御量に相当する。さらに、オンボード同定器２６０が、修正後制御入力設定手段、同定手段および重み関数値算出手段に相当し、周波数整形コントローラが制御入力算出手段に相当する。
【０２１７】
次に、図３２を参照しながら、前述した状態予測器２４０について説明する。同図に示すように、この状態予測器２４０は、前述した図５の状態予測器４０と比較すると、前述した第１〜第３予測値算出部４５〜４７に代えて、第１〜第３予測値算出部２４５〜２４７を備えている点のみが異なっており、それ以外は状態予測器４０と同一の構成を備えているので、以下、異なる点を中心に説明するとともに、状態予測器４０と同一の構成については、同じ符号を付し、その説明は適宜、省略する。
【０２１８】
まず、増幅器４４では、前述した式（３）と同じ下式（６９）により、予測当量比ＰＲＥ＿ＫＡＣＴ０が算出される。
【数６９】
0005400743.tif 000070
【０２１９】
また、第１予測値算出部２４５では、遅延要素４１で１制御サイクル分遅延された空燃比補正係数の値ＫＡＦ（ｋ−１）を用い、下式（７０）により、第１予測値ＰＲＥ＿ＫＡＣＴ＿１が算出される。この式（７０）のモデルパラメータδ，αは、オンボード同定器２６０で同定される。
【数７０】
0005400743.tif 000071
【０２２０】
さらに、第２予測値算出部２４６では、２つの遅延要素４１，４２でそれぞれ１，２制御サイクル分ずつ、遅延された空燃比補正係数の値ＫＡＦ（ｋ−１），ＫＡＦ（ｋ−２）を用い、下式（７１）により、第２予測値ＰＲＥ＿ＫＡＣＴ＿２が算出される。
【数７１】
0005400743.tif 000072
【０２２１】
一方、第３予測値算出部２４７では、３つの遅延要素４１〜４３でそれぞれ１〜３制御サイクル分ずつ、遅延された空燃比補正係数の値ＫＡＦ（ｋ−１），ＫＡＦ（ｋ−２），ＫＡＦ（ｋ−３）を用い、下式（７２）により、第３予測値ＰＲＥ＿ＫＡＣＴ＿３が算出される。
【数７２】
0005400743.tif 000073
【０２２２】
なお、以上の式（７０）〜（７２）は、前述した制御対象モデルの式（６８）に基づき、前述した式（４）〜（６）と同じ手法により導出される。
【０２２３】
さらに、前述したように、４つの重み関数値算出部４８〜５１では、４つの重み関数値Ｗｄ１〜Ｗｄ４がそれぞれ算出され、４つの乗算器５２〜５５では、４つの乗算値Ｗｄ４・ＰＲＥ＿ＫＡＣＴ＿０，Ｗｄ３・ＰＲＥ＿ＫＡＣＴ＿１，Ｗｄ２・ＰＲＥ＿ＫＡＣＴ＿２，Ｗｄ１・ＰＲＥ＿ＫＡＣＴ＿３がそれぞれ算出される。
【０２２４】
そして、加算器５６で、前述した式（７）と同じ下式（７３）により、予測当量比ＰＲＥ＿ＫＡＣＴが算出される。
【数７３】
0005400743.tif 000074
【０２２５】
以上の手法により予測当量比ＰＲＥ＿ＫＡＣＴを算出した場合でも、状態予測器４０と同様の作用効果を得ることができる。すなわち、排ガスボリュームＶｅｘの変化に応じて、むだ時間ｄが値０から値３の間で連続的に変化したときでも、そのようなむだ時間の変化を適切に反映させながら、予測当量比ＰＲＥ＿ＫＡＣＴを算出することができる。
【０２２６】
次に、前述したオンボード同定器２６０について説明する。このオンボード同定器２６０では、以下に述べるスケジュールド補正型の同定アルゴリズムを用いて、モデルパラメータベクトルθが算出される。この同定アルゴリズムは、前述した式（６８）の右辺のＫＡＦ（ｋ−ｄ）を修正後制御入力ＫＡＦ＿ｍｏｄ（ｋ）に置き換えた修正後モデルに基づいて導出される。
【０２２７】
このオンボード同定器２６０は、図３３に示すように、修正後制御入力算出部２７０、３つの遅延要素２６１〜２６３、推定検出当量比算出部２６５、同定ゲインベクトル算出部２６６、減算器２６７、乗算器２６８およびモデルパラメータベクトル算出部２９０を備えている。
【０２２８】
まず、修正後制御入力算出部２７０では、前述した修正後制御入力算出部７０と同一の手法により、修正後制御入力ＫＡＦ＿ｍｏｄが算出される。
【０２２９】
また、推定検出当量比算出部２６５では、遅延要素２６１〜２６３でそれぞれ１制御サイクル分遅延された３つの値ＫＡＣＴ（ｋ−１），ＫＡＦ＿ｍｏｄ（ｋ−１），θ（ｋ−１）を用い、下式（７４）〜（７６）により、推定検出当量比ＫＡＣＴ＿ｈａｔが算出される。
【数７４】
0005400743.tif 000075
【数７５】
0005400743.tif 000076
【数７６】
0005400743.tif 000077
【０２３０】
この式（７６）は、前述した式（６８）の各パラメータを１制御サイクル分過去側にずらした式において、左辺のＫＡＣＴをＫＡＣＴ＿ｈａｔに、右辺のＫＡＦをＫＡＦ＿ｍｏｄにそれぞれ置き換えることによって、導出される。
【０２３１】
さらに、減算器２６７では、下式（７７）により、同定誤差ｅｉｄが算出される。
【数７７】
0005400743.tif 000078
【０２３２】
一方、同定ゲインベクトル算出部２６６では、下式（７８），（７９）により、同定ゲインベクトルＫｐが算出される。この同定ゲインベクトルＫｐは、モデルパラメータベクトルθにおける要素δ，αの修正方向（正負）および修正量を規定するものである。
【数７８】
0005400743.tif 000079
【数７９】
0005400743.tif 000080
【０２３３】
上式（７８）のＩは２次の単位行列であり、Ｐは、その初期値が下式（８０）のように定義される２次の正方行列である。
【数８０】
0005400743.tif 000081
【０２３４】
さらに、上式（７８）の場合、前述したように、重みパラメータλ１、λ２を下記のように設定することにより、同定アルゴリズムとして、３つのアルゴリズムのいずれかを選択することができる。
λ１＝１，λ２＝０  ；固定ゲインアルゴリズム
λ１＝１，λ２＝１  ；最小２乗法アルゴリズム 
λ１＝λ，λ２＝１  ；重み付き最小２乗法アルゴリズム
ここで、λは０＜λ＜１に設定される所定値である。本実施形態の場合、同定精度および制御精度を適切に確保するために、重み付き最小２乗法アルゴリズムを用いている。
【０２３５】
さらに、乗算器２６８で、同定誤差ｅｉｄと同定ゲインベクトルＫｐの乗算値ｅｉｄ・Ｋｐが算出される。
【０２３６】
次いで、モデルパラメータベクトル算出部２９０では、上記の乗算値ｅｉｄ・Ｋｐと排ガスボリュームＶｅｘを用いて、以下に述べるようにモデルパラメータベクトルθが算出される。このモデルパラメータベクトル算出部２９０は、図３４に示すように、基準モデルパラメータ算出部２９１、基準モデルパラメータベクトル算出部３１７、４つの重み関数値算出部２９２〜２９５、８つの乗算器２９６〜３０３、５つの加算器３０４〜３０８、４つの遅延要素３０９〜３１２および４つの増幅器３１３〜３１６を備えている。
【０２３７】
まず、基準モデルパラメータ算出部２９１では、前述した図９の基準モデルパラメータ算出部９１と同一の手法により、基準モデルパラメータαｂｓが算出される。次いで、基準モデルパラメータベクトル算出部３１７では、下式（８１）により、基準モデルパラメータδｂｓを算出した後、下式（８２）により、基準モデルパラメータベクトルθｂｓが算出される。
【数８１】
0005400743.tif 000082
【数８２】
0005400743.tif 000083
【０２３８】
また、４つの重み関数値算出部２９２〜２９５では、前述した図９の重み関数値算出部９２〜９５と同一の手法により、４つの重み関数値Ｗａ１〜Ｗａ４がそれぞれ算出される。また、乗算器２９６で、重み関数値Ｗａ１に値Ｋｐ・ｅｉｄを乗算することにより、乗算値Ｗａ１・Ｋｐ・ｅｉｄが算出され、遅延要素３０９で遅延された修正項ベクトルｄθ１（ｋ−１）に、増幅器３１３で忘却行列ηを乗算することにより、値η・ｄθ１（ｋ−１）が算出される。この忘却行列ηについては後述する。そして、加算器３０４で、値Ｗａ１・Ｋｐ・ｅｉｄに値η・ｄθ１（ｋ−１）を加算することにより、修正項ベクトルｄθ１が算出される。この修正項ベクトルｄθ１は、後述する式（８４）に示すように、２つの修正項ｄδ１，ｄα１を要素とするものである。
【０２３９】
さらに、乗算器２９７で、重み関数値Ｗａ２に値Ｋｐ・ｅｉｄを乗算することにより、乗算値Ｗａ２・Ｋｐ・ｅｉｄが算出され、遅延要素３１０で遅延された修正項ベクトルｄθ２（ｋ−１）に、増幅器３１４で忘却行列ηを乗算することにより、値η・ｄθ２（ｋ−１）が算出される。そして、加算器３０５で、値Ｗａ２・Ｋｐ・ｅｉｄに値η・ｄθ２（ｋ−１）を加算することにより、修正項ベクトルｄθ２が算出される。この修正項ベクトルｄθ２は、後述する式（８４）に示すように、２つの修正項ｄδ２，ｄα２を要素とするものである。
【０２４０】
また、乗算器２９８で、重み関数値Ｗａ３に値Ｋｐ・ｅｉｄを乗算することにより、乗算値Ｗａ３・Ｋｐ・ｅｉｄが算出され、遅延要素３１１で遅延された修正項ベクトルｄθ３（ｋ−１）に、増幅器３１５で忘却行列ηを乗算することにより、値η・ｄθ３（ｋ−１）が算出される。そして、加算器３０６で、値Ｗａ３・Ｋｐ・ｅｉｄに値η・ｄθ３（ｋ−１）を加算することにより、修正項ベクトルｄθ３が算出される。この修正項ベクトルｄθ３は、後述する式（８４）に示すように、３つの修正項ｄδ３，ｄα３を要素とするものである。
【０２４１】
さらに、乗算器２９９で、重み関数値Ｗａ４に値Ｋｐ・ｅｉｄを乗算することにより、乗算値Ｗａ４・Ｋｐ・ｅｉｄが算出され、遅延要素３１２で遅延された修正項ベクトルｄθ４（ｋ−１）に、増幅器３１６で忘却行列ηを乗算することにより、値η・ｄθ４（ｋ−１）が算出される。そして、加算器３０７で、値Ｗａ４・Ｋｐ・ｅｉｄに値η・ｄθ４（ｋ−１）を加算することにより、修正項ベクトルｄθ４が算出される。この修正項ベクトルｄθ４は、後述する式（８４）に示すように、４つの修正項ｄδ４，ｄα４を要素とするものである。
【０２４２】
上述した増幅器３１３〜３１６で用いられる忘却行列ηは、下式（８３）のように定義される。
【数８３】
0005400743.tif 000084
【０２４３】
上式（８３）のη１，η２は、忘却係数であり、０＜η１≦１，０＜η２≦１が成立するように設定される。修正項ベクトルｄθｉ（ｉ＝１〜４）の算出において、この忘却行列ηを用いた理由は、エンジン３の定常運転状態が長時間継続した場合、修正項ベクトルｄθｉが増大化し、不適切な値となるおそれがあるので、それを回避するためである。また、この忘却行列ηにおいて、２つの忘却係数η１，η２のうちの一方を値１に設定した場合、同定誤差ｅｉｄの定常的な発生を抑制することと、制御系の安定性の確保とを両立させることができる。
【０２４４】
また、以上の４つの加算器３０４〜３０７における演算式は、下式（８４），（８５）で表される。
【数８４】
0005400743.tif 000085
【数８５】
0005400743.tif 000086
【０２４５】
さらに、乗算器３００〜３０３ではそれぞれ、４つの重み関数値Ｗａｉを４つの修正項ベクトルｄθｉに乗算することにより、４つのベクトルＷａｉ・ｄθｉが算出される。
【０２４６】
そして、加算器３０８で、下式（８６）により、モデルパラメータベクトルθが最終的に算出される。
【数８６】
0005400743.tif 000087
【０２４７】
オンボード同定器２６０において、以上の同定アルゴリズムを用いたのは、前述した同定条件１，２を満たすためである。すなわち、前述したように、最小２乗法などの一般的な同定アルゴリズムをそのまま用いた場合、同定条件１を満たすことは困難である。そのため、このオンボード同定器２６０の場合、同定条件１を満たしながら、モデルパラメータを同定するために、制御対象モデルの式（６８）のモデルパラメータδ，αの同定演算において、２つのモデルパラメータの基準値δｂｓ，αｂを、両者の間に拘束条件（δｂｓ＝１−αｂｓ）を設定しながら算出するとともに、修正項ベクトルｄθｉの方は一般的な逐次型最小２乗法アルゴリズムを用いて算出する手法を採用している。さらに、同定条件２を満たすために、前述したオンボード同定器６０と同じように、重み関数値Ｗａｉを用いて、修正項ベクトルｄθｉおよびモデルパラメータベクトルθを算出する手法を採用している。
【０２４８】
次に、前述した周波数整形コントローラ３３０について説明する。この周波数整形コントローラ３３０は、予測当量比ＰＲＥ＿ＫＡＣＴが目標当量比ＫＣＭＤに収束するように、空燃比補正係数ＫＡＦを算出するもの、言い換えれば、検出当量比ＫＡＣＴが目標当量比ＫＣＭＤに収束するように、空燃比補正係数ＫＡＦを算出するものである。この周波数整形コントローラ３３０では、まず、前述した式（３４）と同一の下式（８７）により、予測追従誤差ＰＲＥ＿ｅが算出される。
【数８７】
0005400743.tif 000088
【０２４９】
次いで、下式（８８）により、制御入力としての空燃比補正係数ＫＡＦが算出される。
【数８８】
0005400743.tif 000089
【０２５０】
以上の周波数整形コントローラ３３０の制御アルゴリズムは、前述した周波数整形コントローラ１３０の制御アルゴリズムの導出手法と同じ手法により導出される。
【０２５１】
以上のように構成された第２実施形態の制御装置１Ａによれば、第１実施形態の制御装置１と同一の状態予測器４０を用いたことによって、予測検出当量比ＫＡＣＴを、むだ時間ｄの変化を適切に補償しながら精度よく算出することができる。また、周波数整形コントローラ３３０により、前述した周波数整形コントローラ１３０と同様に、制御条件φを満たしかつむだ時間ｄの変化を適切に補償しながら、空燃比補正係数ＫＡＦを算出することができる。すなわち、むだ時間ｄの変化を適切に補償しながら、制御装置１における外乱抑制特性やロバスト性を周波数軸上でダイレクトに指定（設定）することができ、それにより、外乱による検出当量比ＫＡＣＴの変動を抑制したい周波数域において、外乱抑制能力やロバスト性を飛躍的に向上させることができる。
【０２５２】
また、オンボード同定器２６０において、前述したように、修正後制御入力ＫＡＦ＿ｍｏｄおよび重み関数値Ｗａｉを用いる同定アルゴリズムによって、２つのモデルパラメータδ，αが同定されるので、前述した同定条件２を満たしながら、２つのモデルパラメータδ，αを算出することができる。これに加えて、モデルパラメータδ，αの場合、基準モデルパラメータδｂｓ，αｂｓ同士が、同定条件１（すなわち拘束条件）を満たすように設定されているとともに、これらを要素とする基準モデルパラメータベクトルθｂｓを、重み関数値Ｗａｉと修正項ベクトルｄθｉの乗算値の総和で修正することによって、モデルパラメータδ，αを要素とするモデルパラメータベクトルθが算出されるので、２つのモデルパラメータδ，αを、同定条件１を満たす値の近傍の値として同定することができる。
【０２５３】
なお、以上のオンボード同定器２６０の同定アルゴリズムを、前述したオンボード同定器６０の同定アルゴリズムと比較した場合、オンボード同定器６０の同定アルゴリズムでは、同定値αｉｄを同定条件１を完全に満たすように算出できるので、同定条件１を満たすようにモデルパラメータを同定するという観点からは、オンボード同定器６０の同定アルゴリズムの方が優れている。
【０２５４】
次に、図３５を参照しながら、本発明の第３実施形態に係る制御装置１Ｂについて説明する。同図に示すように、この制御装置１Ｂは、前述した第１実施形態の図３に示す制御装置１と比較すると、前述した周波数整形コントローラ１３０に代えて、２自由度応答指定型コントローラ３５０を備えている点のみが異なっており、それ以外は制御装置１と同一の構成を備えているので、以下、２自由度応答指定型コントローラ３５０（制御入力算出手段）についてのみ説明する。
【０２５５】
この２自由度応答指定型コントローラ３５０では、以下の２自由度応答指定型制御アルゴリズムを用いて、空燃比補正係数ＫＡＦが算出される。具体的には、まず、下式（８９）により、目標当量比のフィルタリング値ＫＣＭＤ＿ｆを算出する。
【数８９】
0005400743.tif 000090
ここで、ＰＯＬＥ＿ｆは目標値フィルタ設定パラメータであり、−１＜ＰＯＬＥ＿ｆ＜０の関係が成立するように設定される。
【０２５６】
次いで、下式（９０）により、予測追従誤差ＰＲＥ＿ｅ＿ｆを算出する。
【数９０】
0005400743.tif 000091
【０２５７】
さらに、下式（９１）により、切換関数σ＿ｆを算出する。
【数９１】
0005400743.tif 000092
ここで、ＰＯＬＥは切換関数設定パラメータであり、−１＜ＰＯＬＥ＜０の関係が成立するように設定される。
【０２５８】
次いで、下式（９２）により、等価制御入力Ｕｅｑ＿ｆを算出する。
【数９２】
0005400743.tif 000093
【０２５９】
また、下式（９３）により、到達則入力Ｕｒｃｈ＿ｆを算出する。
【数９３】
0005400743.tif 000094
ここで、Ｋｒｃｈは所定のフィードバックゲインである。
【０２６０】
さらに、下式（９４）により、適応則入力Ｕａｄｐ＿ｆを算出する。
【数９４】
0005400743.tif 000095
ここで、Ｋａｄｐは所定のフィードバックゲインである。
【０２６１】
そして、最終的に、下式（９５）により、空燃比補正係数ＫＡＦが算出される。
【数９５】
0005400743.tif 000096
【０２６２】
なお、以上の式（８９）〜（９５）の２自由度応答指定型アルゴリズムは、前述した式（５３）のＫＡＣＴをＰＲＥ＿ＫＡＣＴに置き換えたモデルに基づいて導出される。
【０２６３】
以上の第３実施形態に係る制御装置１Ｂによれば、第１実施形態の制御装置１と同一の状態予測器４０およびオンボード同定器６０を備えているので、これらによって、第１実施形態の制御装置１と同じ作用効果を得ることができる。これに加えて、２自由度応答指定型コントローラ３５０において、以上の制御アルゴリズムにより、空燃比補正係数ＫＡＦが算出されるので、外乱に起因する、目標当量比ＫＣＭＤと検出当量比ＫＡＣＴとの偏差の値０への時系列的収束挙動と、目標当量比ＫＣＭＤの変化に対する検出当量比ＫＡＣＴの追従特性とを別個に直接的に指定することができる。
【０２６４】
次に、図３６を参照しながら、本発明の第４実施形態に係る制御装置１Ｃについて説明する。同図に示すように、この制御装置１Ｃは、前述した第１実施形態の図３に示す制御装置１と比較した場合、適応外乱オブザーバ３７０（外乱推定値算出手段）を備えている点と、前述した周波数整形コントローラ１３０に代えて２自由度応答指定型コントローラ３８０（制御入力算出手段）を備えている点のみが異なっているので、以下、これらの異なる点についてのみ説明する。
【０２６５】
この適応外乱オブザーバ３７０では、以下に述べる制御アルゴリズムにより、外乱推定値εが算出される。まず、下式（９６）により、外乱推定用の推定検出当量比ＫＡＣＴ＿ａｄｖを算出する。
【数９６】
0005400743.tif 000097
【０２６６】
この式（９６）は、前述した式（２）において、ＫＡＦ（ｋ＋１）をＫＡＣＴ＿ａｄｖ（ｋ）に、αをαｉｄ（ｋ）に、ＫＡＦ（ｋ−ｄ）をＫＡＦ＿ｍｏｄ（ｋ）にそれぞれ置き換えるとともに、右辺に外乱推定値εを加えたもの、すなわち外乱推定用のモデルに相当する。
【０２６７】
次いで、下式（９７）により、外乱推定用の追従誤差ｅ＿ａｄｖを算出する。
【数９７】
0005400743.tif 000098
【０２６８】
そして、最終的に、下式（９８）により、外乱推定値εが算出される。
【数９８】
0005400743.tif 000099
この式（９８）のπは外乱推定ゲインであり、π＞０が成立するように設定される。
【０２６９】
次に、２自由度応答指定型コントローラ３８０について説明する。この２自由度応答指定型コントローラ３８０では、下式（９９）〜（１０４）に示す目標値フィルタ型の２自由度応答指定型制御アルゴリズムよって、空燃比補正係数ＫＡＦが算出される。
【０２７０】
【数９９】
0005400743.tif 000100
【数１００】
0005400743.tif 000101
【数１０１】
0005400743.tif 000102
【数１０２】
0005400743.tif 000103
【数１０３】
0005400743.tif 000104
【数１０４】
0005400743.tif 000105
【０２７１】
以上の式（９９）〜（１０４）は、前述した式（８９）〜（９５）において、外乱推定値εを等価制御入力Ｕｅｑの算出式に付加するとともに、適応則入力Ｕａｄｐを省略したものに相当する。
【０２７２】
以上の第４実施形態に係る制御装置１Ｃによれば、第１実施形態の制御装置１と同一の状態予測器４０およびオンボード同定器６０を備えているので、これらによって、第１実施形態の制御装置１と同じ作用効果を得ることができる。これに加えて、適応外乱オブザーバ３７０において、以上の制御アルゴリズムにより外乱推定値εが算出され、２自由度応答指定型コントローラ３８０において、この外乱推定値εを用いて空燃比補正係数ＫＡＦが算出されるので、空燃比制御における外乱抑制能力すなわちロバスト性を向上させることができる。
【０２７３】
また、適応外乱オブザーバ３７０を備えているので、例えば、π＞Ｐ０に設定し、オンボード同定器６０の同定速度を落とすことによって、制御の安定性を向上させることが可能になる。さらに、同じ理由によって、制御系の共振を防止したり、同定値αｉｄの演算結果を適用した制御対象モデルのゲイン特性が過小になるのを防止するために、同定値αｉｄや同定アルゴリズムで用いる入出力データにフィルタリングを施すことができ、より高い制御性を確保することができる。
【０２７４】
次に、図３７を参照しながら、本発明の第５実施形態に係る制御装置１Ｄについて説明する。なお、以下の説明では、第１実施形態と同じ構成については、同じ符号を付すとともに、その説明は省略する。この制御装置１Ｄは、後述する制御アルゴリズムによって、車両駆動系の自動変速機４００におけるクラッチ４１０の接続・遮断動作などを制御するものである。
【０２７５】
エンジン３は、自動変速機４００および差動ギヤ機構４６０を介して、駆動輪ＷＨ，ＷＨに機械的に連結されており、それにより、エンジン３のトルクは、自動変速機４００および差動ギヤ機構４６０によって変速されながら、駆動輪ＷＨ，ＷＨに伝達される。
【０２７６】
図３７に示すように、自動変速機４００は、クラッチ４１０、主軸４０１、副軸４０２、前進１，２速ギヤ列４２０，４３０、１−２速用同期噛合機構４４０および駆動ギヤ４５０などを備えている。なお、図３７では、前進１，２速ギヤ列４２０，４３０および１−２速用同期噛合機構４４０以外のギヤ列および同期噛合機構は省略されている。
【０２７７】
クラッチ４１０（伝達トルク調整機構）は、乾式クラッチタイプのものであり、エンジン３のクランクシャフト３ａに連結されたクラッチ板４１１と、このクラッチ板４１１と対をなす主軸４０１に連結されたクラッチ板４１２と、クラッチ板４１１をエンジン側に付勢するダイヤフラムスプリング（図示せず）と、クラッチ板４１１をクラッチ板４１２側に駆動するクラッチアクチュエータ４１３などを備えている。
【０２７８】
このクラッチアクチュエータ４１３は、油圧駆動式のものであり、クラッチ電磁弁および油圧アクチュエータなどを組み合わせて構成されている。このクラッチ電磁弁は、ＥＣＵ２に電気的に接続されており、ＥＣＵ２から制御入力信号が供給されるのに伴い、油圧アクチュエータへの供給油圧を変化させる。それにより、クラッチアクチュエータ４１３によるクラッチ板４１１のクラッチ板４１２側への駆動状態が変更されることで、クラッチ４１０の接続・遮断状態が変更される。
【０２７９】
前進１，２速ギヤ列４２０，４３０は、主軸４０１上に回動自在に設けられた１，２速主軸ギヤ４２１，４３１と、副軸４０２上に固定され、１，２速主軸ギヤ４２１，４３１とそれぞれ常に噛み合う１，２速副軸ギヤ４２２，４３２で構成されている。
【０２８０】
また、１，２速主軸ギヤ４２１，４３１の間には、１−２速用同期噛合機構４４０が設けられている。この１−２速用同期噛合機構４４０は、油圧駆動式のものであり、シンクロ電磁弁および油圧アクチュエータなどを組み合わせて構成されている。このシンクロ電磁弁は、ＥＣＵ２に電気的に接続されており、ＥＣＵ２から制御入力信号が供給されるのに伴い、油圧アクチュエータへの供給油圧を変化させる。それにより、１−２速用同期噛合機構４４０は、１速主軸ギヤ４２１または２速主軸ギヤ４３１を主軸４０１に同期させながら噛み合いにより連結する。それにより、前進１速段または前進２速段への変速動作が実行される。
【０２８１】
一方、駆動ギヤ４５０は、差動ギヤ機構４６０の被駆動ギヤ４６１と常に噛み合っており、それにより、副軸４０２の回転に伴い、差動ギヤ機構４６０を介して駆動輪ＷＨ，ＷＨが駆動される。
【０２８２】
さらに、制御装置１Ｄは、ＥＣＵ２を備えており、このＥＣＵ２には、前述したクランク角センサ２０およびアクセル開度センサ２１に加えて、油温センサ２６、４つの車輪速度センサ２７（１つのみ図示）および主軸回転数センサ２８が電気的に接続されている。
【０２８３】
この油温センサ２６は、サーミスタなどで構成され、前述した油圧アクチュエータなどに供給される作動油の温度である油温Ｔｏｉｌを検出して、それを表す検出信号をＥＣＵ２に出力する。ＥＣＵ２は、この油温センサ２６の検出信号に基づいて、油温Ｔｏｉｌを算出する。なお、本実施形態では、油温センサ２６が参照パラメータ検出手段に相当し、油温Ｔｏｉｌが参照パラメータに相当する。
【０２８４】
また、４つの車輪速度センサ２７の各々は、対応する車輪の回転速度を検出して、それを表す検出信号をＥＣＵ２に出力する。ＥＣＵ２は、これらの車輪速度センサ２７の検出信号に基づいて、車速ＶＰなどを算出する。
【０２８５】
さらに、主軸回転数センサ２８は、クランク角センサ２０と同様に、マグネットロータおよびＭＲＥピックアップで構成されており、主軸４０１の回転に伴い、主軸４０１の回転数を表すパルス信号をＥＣＵ２に出力する。ＥＣＵ２は、この主軸回転数センサ２８の検出信号に基づいて、主軸４０１の回転数（以下「主軸回転数」という）ＮＭを算出する。なお、本実施形態では、主軸回転数ＮＭが制御量および出力回転数を表す値に相当する。
【０２８６】
次に、本実施形態の制御装置１Ｄにおけるクラッチ制御の原理について説明する。本実施形態のクラッチ４１０の場合、クラッチアクチュエータ４１３への制御入力Ｕａｃｔと主軸回転数ＮＭとの関係は、下式（１０５）に示すように、一次遅れ系の制御対象モデルとしてモデリングすることができる。
【数１０５】
0005400743.tif 000106
この式（１０５）において、α”はモデルパラメータを、ｄ”はむだ時間をそれぞれ表している。
【０２８７】
また、クラッチ４１０の場合、そのすべり率（またはクランクシャフト３ａと主軸４０１との間の回転差）によって、駆動輪ＷＨ側への伝達トルクが決まるとともに、このすべり率は、クラッチアクチュエータ４１３によるクラッチ板４１１の駆動状態によって調整されるという特性を有している。
【０２８８】
このクラッチアクチュエータ４１３は、前述したように、油圧駆動式のものであるので、油温Ｔｏｉｌの変化に伴って、応答性が変化してしまう特性を有しており、そのため、クラッチ４１０のすべり率すなわちクラッチ４１０のトルク伝達特性は、作動油の温度変化の影響を受けやすいという特性を有している。これに加えて、クラッチ４１０のすべり率は、クラッチ板４１１，４１２の表面温度や各部品の経年変化の影響を受けやすいという特性も有している。
【０２８９】
以上の理由により、前述した式（１０５）のむだ時間ｄ”は、油温Ｔｏｉｌおよびクラッチ板４１１，４１２の表面温度の変化や、各部品の経年変化の影響を受けやすいので、これらに対するロバスト性を確保する必要がある。このむだ時間ｄ”と油温Ｔｏｉｌとの関係をモデリングすると、図３８に示すモデル（マップ）が得られる。同図において、Ｔｏｉｌ１〜Ｔｏｉｌ４，ＴｏｉｌＭＡＸは、油温Ｔｏｉｌの所定値であり、０＜Ｔｏｉｌ１＜Ｔｏｉｌ２＜Ｔｏｉｌ３＜Ｔｏｉｌ４＜ＴｏｉｌＭＡＸが成立するように設定されている。また、所定値ＴｏｉｌＭＡＸは、エンジン３の運転中に油温Ｔｏｉｌが変化し得る領域の最大値に設定されている。言い換えれば、油温Ｔｏｉｌは、エンジン３の運転中、０〜ＴｏｉｌＭＡＸの領域内で変化する特性を有しているので、上記のロバスト性を確保するには、油温Ｔｏｉｌの変化に伴うむだ時間ｄ”の変化を反映させながら、制御入力Ｕａｃｔを算出する必要がある。
【０２９０】
また、一般に、クラッチでは、動作中に「ジャダ」と呼ばれる高周波数の振動的挙動が生じやすく、これが発生した場合、駆動力が振動的に変化し、運転性が低下してしまう。このような問題は、本実施形態のクラッチ４１０のような乾式クラッチではより顕著に発生しやすく、この問題を解消するには、前述した制御条件φを満たす制御アルゴリズムを用いる必要がある。
【０２９１】
以上の理由により、本実施形態では、むだ時間ｄ”を含む前述した式（１０５）の制御対象モデルと、前述した周波数整形コントローラ１３０における制御アルゴリズムと同様の制御アルゴリズムとを用いて、制御入力Ｕａｃｔが算出される。
【０２９２】
以下、本実施形態の制御装置１Ｄの構成および制御アルゴリズムについて説明する。なお、以下に述べる制御アルゴリズムは、車両において、ギヤ段が前進１速段に設定されかつ低速走行中のとき、またはギヤ段が前進１速段に設定されかつ発進中であるときに用いられるアルゴリズムであり、以下の説明では、このようなギヤ段の設定条件および走行条件をまとめて「クラッチ制御条件」という。
【０２９３】
この制御装置１Ｄは、図３９に示すクラッチコントローラ５００と、図４４に示すスロットル弁コントローラ６００とを有しており、これらのコントローラ５００，６００はいずれも、具体的にはＥＣＵ２によって構成されている。
【０２９４】
まず、図３９を参照しながら、クラッチコントローラ５００について説明する。このクラッチコントローラ５００は、上述したクラッチ制御条件が成立しているときに、クラッチ４１０の接続・遮断状態を制御するものであり、同図に示すように、目標主軸回転数算出部５１０、可変むだ時間状態予測器（以下「状態予測器」という）５２０、オンボード・スケジュールド・モデルパラメータ同定器（以下「オンボード同定器」という）５３０および周波数整形コントローラ５４０を備えている。
【０２９５】
この目標主軸回転数算出部５１０では、以下に述べる手法により、目標主軸回転数ＮＭ＿ｃｍｄが算出される。まず、アクセル開度ＡＰおよび車速ＶＰに応じて、図４０に示すマップを検索することにより、目標クラッチすべり率Ｒｓｌｉｐ＿ｃｍｄを算出する。この目標クラッチすべり率Ｒｓｌｉｐ＿ｃｍｄは、クラッチすべり率（クラッチ４１０の入力側回転数と出力側回転数との比ＮＥ／ＮＭ）の目標となる値である。同図において、ＡＰ１〜ＡＰ４は、ＡＰ１＜ＡＰ２＜ＡＰ３＜ＡＰ４が成立するように設定されるアクセル開度ＡＰの所定値であり、特に、ＡＰ１はアクセルペダルが全閉状態のときの値に、ＡＰ４はアクセルペダルが全開状態のときの値に設定されている。また、同図のＶＰ１は、所定車速である。
【０２９６】
同図に示すように、目標クラッチすべり率Ｒｓｌｉｐ＿ｃｍｄは、ＶＰ≦ＶＰ１でかつＡＰ＞ＡＰ１の領域では、アクセル開度ＡＰが大きいほど、または車速ＶＰが高いほど、より小さい値に設定されている。これは、アクセル開度ＡＰが大きいほど、または車速ＶＰが高いほど、クラッチ４１０におけるトルクの伝達効率を高める必要があることによる。
【０２９７】
次いで、以上のように算出した目標クラッチすべり率Ｒｓｌｉｐ＿ｃｍｄを用い、下式（１０６）により、目標主軸回転数ＮＭ＿ｃｍｄが算出される。
【数１０６】
0005400743.tif 000107
【０２９８】
なお、本実施形態では、目標主軸回転数算出部５１０が目標制御量算出手段に相当し、目標主軸回転数ＮＭ＿ｃｍｄが目標制御量に相当する。
【０２９９】
次に、前述した状態予測器５２０について説明する。この状態予測器５２０では、前述した図３８のむだ時間ｄ”の特性を考慮し、前述した第１実施形態の状態予測器４０と同様の予測アルゴリズムを用いて、予測主軸回転数ＰＲＥ＿ＮＭが算出される。なお、本実施形態では、状態予測器５２０が、予測値算出手段、重み関数値算出手段および予測制御量設定手段に相当し、予測主軸回転数ＰＲＥ＿ＮＭが予測制御量に相当する。
【０３００】
この予測主軸回転数ＰＲＥ＿ＮＭは、むだ時間ｄ”が経過したタイミングにおける主軸回転数ＮＭを予測した値に相当するものであり、具体的には、以下の式（１０７）〜（１１１）に示す予測アルゴリズムにより算出される。また、この予測アルゴリズムは、第１実施形態の状態予測器４０の予測アルゴリズムと同じ手法により導出される。
【０３０１】
まず、下式（１０７）により、第０予測値ＰＲＥ＿ＮＭ＿０を算出する。
【数１０７】
0005400743.tif 000108
【０３０２】
また、下式（１０８）により、第１予測値ＰＲＥ＿ＮＭ＿１を算出する。
【数１０８】
0005400743.tif 000109
この式（１０８）のαｉｄ”は、モデルパラメータα”の同定値であり、後述するように、オンボード同定器５３０で算出される。
【０３０３】
さらに、下式（１０９）により、第２予測値ＰＲＥ＿ＮＭ＿２を算出する。
【数１０９】
0005400743.tif 000110
【０３０４】
次いで、下式（１１０）により、第３予測値ＰＲＥ＿ＮＭ＿３を算出する。
【数１１０】
0005400743.tif 000111
【０３０５】
そして、最終的に、下式（１１１）により、予測主軸回転数ＰＲＥ＿ＮＭが算出される。
【数１１１】
0005400743.tif 000112
【０３０６】
上式（１１１）のＷｄｉ”（ｉ＝１〜４）は、重み関数値であり、この重み関数値Ｗｄｉ”（ｉ＝１〜４）は、油温Ｔｏｉｌに応じて、図４１に示すマップを検索することにより算出される。同図に示すように、４つの重み関数値Ｗｄ１”〜Ｗｄ４”はそれぞれ、油温Ｔｏｉｌが変化し得る領域をＴｏｉｌ≦Ｔｏｉｌ２，Ｔｏｉｌ１≦Ｔｏｉｌ≦Ｔｏｉｌ３，Ｔｏｉｌ２≦Ｔｏｉｌ≦Ｔｏｉｌ４，Ｔｏｉｌ３≦Ｔｏｉｌ≦ＴｏｉｌＭＡＸの４つの領域に区分した場合において、これらの４つの領域に対応するように設定されているとともに、対応する領域では、値１以下の正の値に設定され、対応する領域以外では、値０に設定されている。
【０３０７】
具体的には、重み関数値Ｗｄ１”は、これが対応する領域（Ｔｏｉｌ≦Ｔｏｉｌ２）では、Ｔｏｉｌ≦Ｔｏｉｌ１のときの値１を最大値として、油温Ｔｏｉｌが大きいほど、より小さい正の値に設定されているとともに、それ以外の領域では、値０に設定されている。また、重み関数値Ｗｄ２”は、これが対応する領域（Ｔｏｉｌ１≦Ｔｏｉｌ≦Ｔｏｉｌ３）では、Ｔｏｉｌ＝Ｔｏｉｌ２のときの値１を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値０に設定されている。
【０３０８】
さらに、重み関数値Ｗｄ３”は、これが対応する領域（Ｔｏｉｌ２≦Ｔｏｉｌ≦Ｔｏｉｌ４）では、Ｔｏｉｌ＝Ｔｏｉｌ３のときの値１を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値０に設定されている。一方、重み関数値Ｗｄ４”は、これが対応する領域（Ｔｏｉｌ３≦Ｔｏｉｌ≦ＴｏｉｌＭＡＸ）では、Ｔｏｉｌ４≦Ｔｏｉｌのときの値１を最大値として、油温Ｔｏｉｌが大きいほど、より大きい正の値に設定されているとともに、それ以外の領域では、値０に設定されている。
【０３０９】
以上に加えて、４つの重み関数値Ｗｄｉ”（ｉ＝１〜４）がそれぞれ対応する４つの領域は、上述したような、隣り合う領域が互いに重なり合うように設定されており、これらの互いに重なり合う領域における油温Ｔｏｉｌのいずれかの値に対応する重み関数値Ｗｄｉ”の値の和は、各重み関数値Ｗｄｉ”における最大値１と等しくなるように設定されている。
【０３１０】
さらに、図４１と前述した図３８を比較すると明らかなように、互いに重なり合う３つの領域は、むだ時間ｄ”の勾配が一定状態に保持される３つの領域に対応するように設定されている。これに加えて、重み関数値Ｗｄ１”はむだ時間ｄ”＝３に対して、重み関数値Ｗｄ２”はむだ時間ｄ”＝２に対して、重み関数値Ｗｄ３”はむだ時間ｄ”＝１に対して、重み関数値Ｗｄ４”はむだ時間ｄ”＝０に対して重みが最も大きくなるようにそれぞれ設定されている。
【０３１１】
したがって、このように設定されている４つの重み関数値Ｗｄｉ”を、４つの予測値ＰＲＥ＿ＮＭ＿４−ｉにそれぞれ乗算した値の総和として、予測主軸回転数ＰＲＥ＿ＮＭが算出されるので、油温Ｔｏｉｌの変化に応じて、むだ時間ｄ”が図３８に示すように値０から値３の間で連続的に変化したときでも、そのようなむだ時間ｄ”の変化を適切に反映させながら、予測主軸回転数ＰＲＥ＿ＮＭを円滑に変化するような値として算出することができる。
【０３１２】
次に、前述したオンボード同定器５３０について説明する。なお、本実施形態では、オンボード同定器５３０が、修正後制御入力設定手段、同定手段および重み関数値算出手段に相当する。このオンボード同定器５３０では、下式（１１２）〜（１２４）に示す拘束条件付きスケジュールド補正型の同定アルゴリズムを用いて、同定値αｉｄ”が算出される。この同定アルゴリズムは、前述したオンボード同定器６０の同定アルゴリズムと同じ手法によって導出される。
【０３１３】
まず、下式（１１２）により、修正後制御入力Ｕａｃｔ＿ｍｏｄを算出する。
【数１１２】
0005400743.tif 000113
【０３１４】
また、下式（１１３）により、合成信号値Ｗ＿ａｃｔ”を算出する。
【数１１３】
0005400743.tif 000114
【０３１５】
さらに、下式（１１４），（１１５）により、推定合成信号値Ｗ＿ｈａｔ”を算出する。
【数１１４】
0005400743.tif 000115
【数１１５】
0005400743.tif 000116
【０３１６】
次いで、下式（１１６）により、同定誤差ｅｉｄ”を算出する。
【数１１６】
0005400743.tif 000117
【０３１７】
また、下式（１１７），（１１８）により、同定ゲインＫｐ”を算出する。
【数１１７】
0005400743.tif 000118
【数１１８】
0005400743.tif 000119
【０３１８】
上式（１１７）のゲインＰ”の初期値Ｐ”（０）は、下式（１１９）のように定義される。
【数１１９】
0005400743.tif 000120
ここで、Ｐ０”は所定値に設定されている。
【０３１９】
また、上式（１１７）のλ１，λ２は重みパラメータであり、前述したように、これらの値λ１、λ２を下記のように設定することにより、同定アルゴリズムとして、３つのアルゴリズムのいずれかを選択することができる。
λ１＝１，λ２＝０  ；固定ゲインアルゴリズム
λ１＝１，λ２＝１  ；最小２乗法アルゴリズム 
λ１＝λ，λ２＝１  ；重み付き最小２乗法アルゴリズム
ここで、λは０＜λ＜１に設定される所定値である。本実施形態の場合、同定精度および制御精度を適切に確保するために、重み付き最小２乗法アルゴリズムを用いている。
【０３２０】
次いで、下式（１２０）〜（１２２）により、ゲイン係数Ｈ”を算出する。
【数１２０】
0005400743.tif 000121
【数１２１】
0005400743.tif 000122
【数１２２】
0005400743.tif 000123
【０３２１】
上式（１２０）〜（１２２）において、α＿Ｌ”は所定の下限値であり、α＿Ｈ”は所定の上限値である。また、η”は忘却係数であり、０＜η”≦１が成立するように設定される。同定値αｉｄ”の算出において、この忘却係数η”を用いた理由は、定常状態が長時間継続した場合、同定値αｉｄ”が増大化し、不適切な値となるおそれがあるので、それを回避するためである。
【０３２２】
さらに、下式（１２３）により、４つの修正項ｄαｉ”（ｉ＝１〜４）を算出する。
【数１２３】
0005400743.tif 000124
【０３２３】
上式（１２３）のＷａｉ”は、重み関数値であり、油温Ｔｏｉｌに応じて、図４２に示すマップを検索することにより算出される。同図において、Ｔｏｉｌ５〜Ｔｏｉｌ８は、油温Ｔｏｉｌの所定値であり、Ｔｏｉｌ５≦Ｔｏｉｌ６≦Ｔｏｉｌ７≦Ｔｏｉｌ８≦ＴｏｉｌＭＡＸが成立するように設定されている。同図に示すように、４つの重み関数値Ｗａ１”〜Ｗａ４”はそれぞれ、油温Ｔｏｉｌが変化し得る領域をＴｏｉｌ≦Ｔｏｉｌ６，Ｔｏｉｌ５≦Ｔｏｉｌ≦Ｔｏｉｌ７，Ｔｏｉｌ６≦Ｔｏｉｌ≦Ｔｏｉｌ８，Ｔｏｉｌ７≦Ｔｏｉｌ≦ＴｏｉｌＭＡＸの４つの領域に区分した場合において、これらの４つの領域に対応するように設定されているとともに、対応する領域では、値１以下の正の値に設定され、対応する領域以外では、値０に設定されている。
【０３２４】
また、重み関数値Ｗａ１”は、これが対応する領域（Ｔｏｉｌ≦Ｔｏｉｌ６）では、Ｔｏｉｌ≦Ｔｏｉｌ５のときの値１を最大値として、油温Ｔｏｉｌが大きいほど、より小さい正の値に設定されているとともに、それ以外の領域では、値０に設定されている。また、重み関数値Ｗａ２”は、これが対応する領域（Ｔｏｉｌ５≦Ｔｏｉｌ≦Ｔｏｉｌ７）では、Ｔｏｉｌ＝Ｔｏｉｌ６のときの値１を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値０に設定されている。
【０３２５】
さらに、重み関数値Ｗａ３”は、これが対応する領域（Ｔｏｉｌ６≦Ｔｏｉｌ≦Ｔｏｉｌ８）では、Ｔｏｉｌ＝Ｔｏｉｌ７のときの値１を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値０に設定されている。一方、重み関数値Ｗａ４”は、これが対応する領域（Ｔｏｉｌ７≦Ｔｏｉｌ≦ＴｏｉｌＭＡＸ）では、Ｔｏｉｌ８≦Ｔｏｉｌのときの値１を最大値として、油温Ｔｏｉｌが大きいほど、より大きい正の値に設定されているとともに、それ以外の領域では、値０に設定されている。
【０３２６】
以上に加えて、４つの重み関数値Ｗａｉ”（ｉ＝１〜４）がそれぞれ対応する４つの領域は、上述したような、隣り合う領域が互いに重なり合うように設定されており、これらの互いに重なり合う領域における油温Ｔｏｉｌのいずれかの値に対応する重み関数値Ｗａｉ”の値の和は、各重み関数値Ｗａｉ”における最大値１と等しくなるように設定されている。この図４２を後述する図４３と比較すると明らかなように、これらの互いに重なり合う３つの領域は、基準モデルパラメータαｂｓ”の勾配が一定状態に保持される３つの領域に対応するように設定されている。
【０３２７】
そして、下式（１２４）により、同定値αｉｄ”が最終的に算出される。
【数１２４】
0005400743.tif 000125
【０３２８】
上式（１２４）のαｂｓ”は、基準モデルパラメータであり、油温Ｔｏｉｌに応じて、図４３に示すマップを検索することにより算出される。このマップでは、基準モデルパラメータαｂｓ”は、油温Ｔｏｉｌが高いほど、より大きい値に設定されている。これは、油温Ｔｏｉｌが高いほど、クラッチアクチュエータの応答性が高くなり、応答遅れが小さくなることで、制御入力Ｕａｃｔが主軸回転数ＮＭに対して及ぼす影響の度合がより大きくなるので、それに対応するためである。
【０３２９】
次に、前述した周波数整形コントローラ５４０（制御入力算出手段）について説明する。この周波数整形コントローラ５４０では、目標主軸回転数ＮＭ＿ｃｍｄ、予測主軸回転数ＰＲＥ＿ＮＭおよび同定値αｉｄ”を用いて、下式（１２５），（１２６）により、制御入力Ｕａｃｔが算出される。なお、下式（１２５），（１２６）の制御アルゴリズムは、前述した周波数整形コントローラ１３０の制御アルゴリズムと同じ原理によって導出される。
【数１２５】
0005400743.tif 000126
【数１２６】
0005400743.tif 000127
【０３３０】
上式（１２５）のＰＲＥ＿ｅ”は予測追従誤差である。上式（１２６）のβ”は感度設定パラメータであり、前述した制御条件φを満たすように設定されている。
【０３３１】
周波数整形コントローラ５４０では、以上のように制御入力Ｕａｃｔが算出される。そして、ＥＣＵ２によって、この制御入力Ｕａｃｔに対応する制御入力信号がクラッチアクチュエータ４１３に供給されることにより、主軸回転数ＮＭが目標主軸回転数ＮＭ＿ｃｍｄに収束するようにフィードバック制御される。
【０３３２】
次に、図４４を参照しながら、前述したスロットル弁コントローラ６００について説明する。このスロットル弁コントローラ６００は、スロットル弁６ａの開度を制御するものであり、同図に示すように、目標エンジントルク算出部６１０、目標ＴＨ開度算出部６２０およびＴＨコントローラ６３０を備えている。
【０３３３】
この目標エンジントルク算出部６１０では、アクセル開度ＡＰおよび車速ＶＰに応じて、図４５に示すマップを検索することにより、目標エンジントルクＴＲＱ＿ＥＮＧ＿ｃｍｄが算出される。同図において、ＴＲＱ＿ＭＡＸは、エンジン３が発生可能なトルクの最大値である。また、図中のハッチングで示す領域は、アクセルペダルが全閉状態にあり（ＡＰ＝ＡＰ１）、かつ車両が走行状態にある（ＶＰ＞ＶＰ１）ことで、フューエルカット運転を実行すべき領域であり、そのため、この領域では、目標エンジントルクＴＲＱ＿ＥＮＧ＿ｃｍｄが負値に設定されている。
【０３３４】
また、目標ＴＨ開度算出部６２０では、目標エンジントルクＴＲＱ＿ＥＮＧ＿ｃｍｄおよびエンジン回転数ＮＥに応じて、図４６に示すマップを検索することにより、目標ＴＨ開度ＴＨ＿ｃｍｄが算出される。同図において、ＮＥ５〜７は、エンジン回転数ＮＥの所定値であり、０＜ＮＥ５＜ＮＥ６＜ＮＥ７＜ＮＥＭＡＸが成立するように設定されている。このマップの場合、目標ＴＨ開度ＴＨ＿ｃｍｄは、高回転域では、目標エンジントルクＴＲＱ＿ＥＮＧ＿ｃｍｄが大きいほど、それを実現可能な吸入空気量を確保するために、より大きい値に設定されている。また、目標ＴＨ開度ＴＨ＿ｃｍｄは、エンジン回転数ＮＥが高いほど、それを実現可能な吸入空気量を確保するために、より大きい値に設定されている。
【０３３５】
次に、ＴＨコントローラ６３０では、目標ＴＨ開度ＴＨ＿ｃｍｄに応じて、図示しないマップを検索することにより、制御入力Ｕｔｈが算出される。そして、ＥＣＵ２によって、この制御入力Ｕｔｈに対応する制御入力信号がＴＨアクチュエータ６ｂに供給されることにより、スロットル弁６ａの開度が目標ＴＨ開度ＴＨ＿ｃｍｄになるように制御される。
【０３３６】
次に、図４７を参照しながら、第５実施形態の制御装置１Ｄによるクラッチ制御のシミュレーション結果（以下「制御結果」という）について説明する。同図において、Ｄｓｌｉｐは、すべり率偏差であり、実際のクラッチすべり率Ｒｓｌｉｐ（＝ＮＥ／ＮＭ）と、目標クラッチすべり率Ｒｓｌｉｐ＿ｃｍｄとの偏差（＝Ｒｓｌｉｐ−Ｒｓｌｉｐ＿ｃｍｄ）を表している。
【０３３７】
同図に示すように、時刻ｔ１で、アクセルペダルが踏まれ、アクセル開度ＡＰがＡＰ１（＝０）から上昇すると、その直後、実際のクラッチすべり率Ｒｓｌｉｐが目標クラッチすべり率Ｒｓｌｉｐ＿ｃｍｄに対してオーバーシュートすることで、すべり率偏差Ｄｓｌｉｐが一時的に急増する。しかし、制御の進行に伴って、すべり率偏差Ｄｓｌｉｐが減少するとともに、時刻ｔ２〜ｔ３の間では、すべり率偏差Ｄｓｌｉｐが値０近傍に保持されており、高い制御精度が確保されていることが判る。
【０３３８】
そして、時刻ｔ３で、アクセルペダルが開放された以降、実際のクラッチすべり率Ｒｓｌｉｐが目標クラッチすべり率Ｒｓｌｉｐ＿ｃｍｄに対してアンダーシュートすることで、すべり率偏差Ｄｓｌｉｐが一時的に急減する。しかし、制御の進行に伴って、すべり率偏差Ｄｓｌｉｐが値０に向かって上昇するとともに、時刻ｔ４〜ｔ５の間では、すべり率偏差Ｄｓｌｉｐが値０近傍に保持されており、高い制御精度が確保されていることが判る。
【０３３９】
そして、時刻ｔ５で、アクセルペダルが再度、踏まれると、その直後、実際のクラッチすべり率Ｒｓｌｉｐが目標クラッチすべり率Ｒｓｌｉｐ＿ｃｍｄに対してオーバーシュートすることで、すべり率偏差Ｄｓｌｉｐが一時的に急増する。その後、制御の進行に伴って、すべり率偏差Ｄｓｌｉｐが減少する。そして、時刻ｔ６以降、クラッチ４１０が直結状態になることで、すべり率偏差Ｄｓｌｉｐが値０に保持される。
【０３４０】
以上のように、第５実施形態の制御装置１Ｄによれば、状態予測器５２０において、主軸回転数ＮＭと制御入力Ｕａｃｔの関係を定義した制御対象モデル（式（１０５））を用いて、むだ時間ｄ”＝０〜３がそれぞれ経過したタイミングでの主軸回転数ＮＭとして、第０〜第３予測値ＰＲＥ＿ＮＭ＿０〜３が算出され、油温Ｔｏｉｌに応じて、４つの重み関数値Ｗｄ１”〜Ｗｄ４”が算出される。そして、予測当量比ＰＲＥ＿ＮＭが、重み関数値Ｗｄｉ”と予測値ＰＲＥ＿ＮＭ＿４−ｉ（ｉ＝１〜４）との乗算値の総和として算出されるので、予測当量比ＰＲＥ＿ＮＭを、予測値ＰＲＥ＿ＮＭ＿４−ｉを連続的に結合させた値として算出することができる。それにより、油温Ｔｏｉｌの変化に応じてむだ時間ｄ”が変化した場合でも、予測主軸回転数ＮＭを、そのようなむだ時間ｄ”の変化を適切に補償しながら精度よく算出することができる。
【０３４１】
また、オンボード同定器５３０において、前述した同定アルゴリズムにより、同定値αｉｄ”が算出されるので、前述した同定条件１，２の双方を満たしながら、同定値αｉｄ”を算出することができる。具体的には、合成信号値Ｗ＿ａｃｔ”および推定合成信号値Ｗ＿ｈａｔ”が互いに一致するように、同定値αｉｄ”が算出されるので、同定条件１すなわち拘束条件を満たしながら、同定値αｉｄ”を算出することができる。さらに、修正後制御入力Ｕａｃｔ＿ｍｏｄが、４つの重み関数値Ｗｄ４”〜Ｗｄ１”を、むだ時間ｄ”＝０〜３前のタイミングでの制御入力Ｕａｃｔ（ｋ），Ｕａｃｔ（ｋ−１），Ｕａｃｔ（ｋ−２），Ｕａｃｔ（ｋ−３）にそれぞれ乗算した乗算値の総和として算出されるので、油温Ｔｏｉｌの変化に応じてむだ時間ｄ”が連続的に変化する場合でも、修正後制御入力Ｕａｃｔ＿ｍｏｄを、そのようなむだ時間ｄ”の変化を適切に補償しながら精度よく算出することができる。
【０３４２】
さらに、そのような修正後制御入力Ｕａｃｔ＿ｍｏｄを用いた式（１７）〜（２９）の同定アルゴリズムにより、同定値αｉｄ”がオンボード同定されるので、油温Ｔｏｉｌの変化に応じてむだ時間ｄ”が変化したときでも、その影響を抑制しながら、同定値αｉｄ”を精度よく同定することができる。特に、油温Ｔｏｉｌの急変に伴ってむだ時間ｄ”が急変したときでも、同定値αｉｄ”を、むだ時間ｄ”の急変を適切に補償しながら、段差なく円滑に変化するように算出することができる。そして、そのように算出された同定値αｉｄ”を用いて、制御入力Ｕａｃｔが算出されるので、クラッチ制御における制御性や、個体間のばらつきおよび経年変化の影響に対するクラッチ制御のロバスト性を飛躍的に向上させることができる。
【０３４３】
これに加えて、周波数整形コントローラ５４０の場合、第１実施形態の周波数整形コントローラ１３０と同じ手法によって導出された式（１２５），（１２６）を用いて、制御入力Ｕａｃｔが算出されるので、前述した制御条件φを満たしながら、制御入力Ｕａｃｔを算出することができる。これに加えて、制御対象モデルのモデルパラメータとして前述した同定値αｉｄ”を用いているので、むだ時間ｄ”の変化を適切に補償しながら、制御装置１Ｄにおける外乱抑制特性やロバスト性を周波数軸上でダイレクトに指定（設定）することができる。それにより、外乱による主軸回転数ＮＭの変動を抑制したい周波数域において、外乱抑制能力やロバスト性を飛躍的に向上させることができる。これに加えて、制御入力Ｕａｃｔの算出アルゴリズムとして、フィードバック制御アルゴリズムが用いられているので、フィードバックゲインを高い状態に維持することができ、高精度かつ高応答性を確保しながら、主軸回転数ＮＭを目標主軸回転数ＮＭに追従させることができる。
【０３４４】
なお、第５実施形態は、重み関数値として、互いに重なり合う領域における油温Ｔｏｉｌのいずれかの値に対応する重み関数値Ｗｄｉ”の値の和が、各重み関数値Ｗｄｉ”における最大値１と等しくなるように設定されているものを用いた例であるが、本発明の重み関数値はこれに限らず、互いに重なり合う領域における参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が所定値になるように設定されているものであればよい。例えば、互いに重なり合う領域における参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が、重み関数値における絶対値の最大値と等しくなるように設定されたものを用いてもよい。より具体的には、重み関数値として、図４１の重み関数値Ｗｄｉ”の設定値に対して、ｘ軸を対称軸として線対称に設定されたもの、すなわち負値に設定されたものを用いてもよく、その場合には、４つの重み関数値に乗算する値、すなわち４つの予測値ＰＲＥ＿ＮＭ＿４−ｉや、４つの制御入力Ｕａｃｔ（ｋ−４＋ｉ）として、負値化した値を用いればよい。
【０３４５】
次に、図４８を参照しながら、本発明の第６実施形態に係る制御装置１Ｅについて説明する。この制御装置１Ｅは、第５実施形態の制御装置１Ｄと同様に、自動変速機４００におけるクラッチの接続・遮断動作などを制御するものである。この第６実施形態の場合、その機械的な構成は、クラッチとして、乾式クラッチタイプのクラッチ４１０に代えて、湿式クラッチ（図示せず）を用いている点以外は、第５実施形態のものと同一であるので、以下の説明では、第５実施形態と同じ構成については、同じ符号を付すとともに、その説明は省略する。
【０３４６】
ここで、一般に、湿式クラッチの場合、構造上の理由から、乾式クラッチと比べてジャダが発生しにくいという特性を有している。そのため、前述した制御条件φを考慮することなく、クラッチの上流側回転数ＮＥと下流側回転数ＮＭとの間の回転差が時系列的に円滑に値０に収束するように、クラッチを制御すればよいことになる。以上の理由により、本実施形態の制御装置１Ｅでは、以下に述べる制御アルゴリズムを用いて、制御入力Ｕａｃｔが算出される。
【０３４７】
この制御装置１Ｅは、図４８に示すように、クラッチコントローラ７００を備えている。このクラッチコントローラ７００の場合、前述した図３９のクラッチコントローラ５００と比較すると、適応外乱オブザーバ７４０（外乱推定値算出手段）を備えている点と、前述した周波数整形コントローラ５４０に代えて２自由度応答指定型コントローラ７５０（制御入力算出手段）を備えている点のみが異なっているので、以下、これらの異なる点についてのみ説明する。
【０３４８】
最初に、適応外乱オブザーバ７４０について説明する。この適応外乱オブザーバ７４０では、以下に述べる制御アルゴリズムにより、外乱推定値ε”が算出される。まず、下式（１２７）により、外乱推定用の推定主軸回転数ＮＭ＿ａｄｖ（推定制御量）を算出する。
【数１２７】
0005400743.tif 000128
【０３４９】
この式（１２７）は、前述した式（１０５）において、ＮＭ（ｋ＋１）をＮＭ＿ａｄｖ（ｋ）に、α”をαｉｄ”（ｋ）に、Ｕａｃｔ（ｋ−ｄ”）をＵａｃｔ＿ｍｏｄ（ｋ）にそれぞれ置き換えるとともに、右辺に外乱推定値ε”を加えたものに相当する。
【０３５０】
次いで、下式（１２８）により、外乱推定用の追従誤差ｅ＿ａｄｖ”を算出する。
【数１２８】
0005400743.tif 000129
【０３５１】
そして、最終的に、下式（１２９）により、外乱推定値ε”が算出される。
【数１２９】
0005400743.tif 000130
この式（１２９）のπ”は外乱推定ゲインであり、π”＞０が成立するように設定される。
【０３５２】
次に、前述した２自由度応答指定型コントローラ７５０について説明する。この２自由度応答指定型コントローラ７５０では、以下に述べるように、上記外乱推定値ε”を加味した応答指定型制御アルゴリズムよって、制御入力Ｕａｃｔが算出される。
【０３５３】
具体的には、まず、下式（１３０）により、目標主軸回転数のフィルタリング値ＮＭ＿ｃｍｄ＿ｆを算出する。
【数１３０】
0005400743.tif 000131
ここで、ＰＯＬＥ＿ｆ”は目標値フィルタ設定パラメータであり、−１＜ＰＯＬＥ＿ｆ”＜０の関係が成立するように設定される。
【０３５４】
次いで、下式（１３１）により、予測追従誤差ＰＲＥ＿ｅ＿ｆ”を算出する。
【数１３１】
0005400743.tif 000132
【０３５５】
さらに、下式（１３２）により、切換関数σ＿ｆ”を算出する。
【数１３２】
0005400743.tif 000133
ここで、ＰＯＬＥ”は切換関数設定パラメータであり、−１＜ＰＯＬＥ”＜０の関係が成立するように設定される。
【０３５６】
次いで、下式（１３３）により、等価制御入力Ｕｅｑ＿ｆ”を算出する。
【数１３３】
0005400743.tif 000134
【０３５７】
また、下式（１３４）により、到達則入力Ｕｒｃｈ＿ｆ”を算出する。
【数１３４】
0005400743.tif 000135
ここで、Ｋｒｃｈ”は所定のフィードバックゲインである。
【０３５８】
そして、最終的に、下式（１３５）により、制御入力Ｕａｃｔが算出される。
【数１３５】
0005400743.tif 000136
【０３５９】
以上の第６実施形態の制御装置１Ｅによれば、第５実施形態の制御装置１Ｄと同一の状態予測器５２０およびオンボード同定器５３０を備えているので、これらによって、第５実施形態の制御装置１Ｄと同じ作用効果を得ることができる。これに加えて、適応外乱オブザーバ７４０において、以上の制御アルゴリズムにより外乱推定値ε”が算出され、２自由度応答指定型コントローラ７５０において、この外乱推定値ε”を用いて制御入力Ｕａｃｔが算出されるので、クラッチ制御における外乱抑制能力すなわちロバスト性を向上させることができる。
【０３６０】
また、適応外乱オブザーバ７４０を備えているので、例えば、π”＞Ｐ０”に設定し、オンボード同定器５３０の同定速度を落とすことによって、制御の安定性を向上させることが可能になる。さらに、同じ理由によって、制御系の共振を防止したり、同定値αｉｄ”の演算結果を適用した制御対象モデルのゲイン特性が過小になるのを防止するために、同定値αｉｄ”や同定アルゴリズムで用いる入出力データにフィルタリングを施すことができ、より高い制御性を確保することができる。
【０３６１】
なお、第１〜第４実施形態は、本発明の制御装置を制御対象としてのエンジン３の空燃比を制御するものに適用し、第５，６実施形態は、本発明の制御装置を制御対象としてのクラッチ４１０を制御するものに適用した例であるが、本発明の制御装置はこれらに限らず、参照パラメータに応じてむだ時間を含む動特性が変化する特性を備えた制御対象を制御するものに適用可能である。例えば、本発明の制御装置を、制御対象としてのロボットの動作を制御するものに適用してもよい。
【０３６２】
また、各実施形態は、本発明の制御装置をむだ時間が４個の整数値（値０〜３）の間で変化する特性の制御対象に適用した例であるが、本発明の制御装置はこれに限らず、むだ時間がＭ個の整数値の間で変化する特性の制御対象に適用可能である。例えば、本発明を制御装置を、むだ時間が３個以下や５個以上の整数値の間で変化する制御対象に適用してもよい。
【符号の説明】
【０３６３】
１  制御装置
1A〜1E  制御装置
２  ＥＣＵ（目標制御量設定手段、参照パラメータ検出手段、予測値算出手段、重み        関数値算出手段、予測制御量設定手段、制御入力算出手段、修正後制御入力設定        手段、同定手段、外乱推定値算出手段）
３  内燃機関
２０  クランク角センサ（参照パラメータ検出手段）
２２  吸気圧センサ（参照パラメータ検出手段）
２６  油温センサ（参照パラメータ検出手段）
３０  目標当量比算出部（目標制御量設定手段）
４０  可変むだ時間状態予測器（予測値算出手段、重み関数値算出手段、予測制御量設        定手段）
６０  オンボード同定器（修正後制御入力設定手段、同定手段、重み関数値算出手段）
１３０  周波数整形コントローラ（制御入力算出手段）
２３０  目標当量比算出部（目標制御量設定手段）
２４０  可変むだ時間状態予測器（予測値算出手段、重み関数値算出手段、予測制御量設        定手段）
２６０  オンボード同定器（修正後制御入力設定手段、同定手段、重み関数値算出手段）
３３０  周波数整形コントローラ（制御入力算出手段）
３５０  ２自由度応答指定型コントローラ（制御入力算出手段）
３７０  適応外乱オブザーバ（外乱推定値算出手段）
３８０  ２自由度応答指定型コントローラ（制御入力算出手段）
４００  自動変速機
４１０  クラッチ（伝達トルク調整機構）
５１０  目標主軸回転数算出部（目標制御量設定手段）
５２０  可変むだ時間状態予測器（予測値算出手段、重み関数値算出手段、予測制御量設        定手段）
５３０  オンボード同定器（修正後制御入力設定手段、同定手段、重み関数値算出手段）
５４０  周波数整形コントローラ（制御入力算出手段）
７４０  適応外乱オブザーバ（外乱推定値算出手段）
７５０  ２自由度応答指定型コントローラ（制御入力算出手段）
ｄ  むだ時間
Ｖｅｘ  排ガスボリューム（参照パラメータ）
ＫＡＣＴ  検出当量比（制御量、空燃比を表す値）
ＫＡＦ  空燃比補正係数（制御入力、補正係数）
ＫＣＭＤ  目標当量比（目標制御量）
ＰＲＥ＿ＫＡＣＴ＿０〜３  第０〜第３予測値（Ｍ個の予測値）
Ｗｄ１〜Ｗｄ４  重み関数値
ＰＲＥ＿ＫＡＣＴ  予測当量比（予測制御量）
ＫＡＦ＿ｍｏｄ  修正後制御入力
αｉｄ  同定値（モデルパラメータ）
α  モデルパラメータ
δ  モデルパラメータ
Ｓｄ  感度関数
Ｔｄ  相補感度関数
Ｃ  伝達関数
ε  外乱推定値
ＫＡＣＴ＿ａｄｖ  外乱推定用の推定検出当量比（推定制御量）
ｄ”  むだ時間
Ｔｏｉｌ  油温（参照パラメータ）
ＮＭ  主軸回転数（制御量、出力回転数を表す値）
Ｕａｃｔ  制御入力
ＮＭ＿ｃｍｄ  目標主軸回転数（目標制御量）
ＰＲＥ＿ＮＭ＿０〜３  第０〜第３予測値（Ｍ個の予測値）
Ｗｄ１”〜Ｗｄ４”  重み関数値
ＰＲＥ＿ＮＭ  予測主軸回転数（予測制御量）
Ｕａｃｔ＿ｍｏｄ  修正後制御入力
αｉｄ”  同定値（モデルパラメータ）
ＮＭ＿ａｄｖ  外乱推定用の推定主軸回転数（推定制御量）
【図１】
0005400743.tif 000137
【図２】
0005400743.tif 000138
【図３】
0005400743.tif 000139
【図４】
0005400743.tif 000140
【図５】
0005400743.tif 000141
【図６】
0005400743.tif 000142
【図７】
0005400743.tif 000143
【図８】
0005400743.tif 000144
【図９】
0005400743.tif 000145
【図１０】
0005400743.tif 000146
【図１１】
0005400743.tif 000147
【図１２】
0005400743.tif 000148
【図１３】
0005400743.tif 000149
【図１４】
0005400743.tif 000150
【図１５】
0005400743.tif 000151
【図１６】
0005400743.tif 000152
【図１７】
0005400743.tif 000153
【図１８】
0005400743.tif 000154
【図１９】
0005400743.tif 000155
【図２０】
0005400743.tif 000156
【図２１】
0005400743.tif 000157
【図２２】
0005400743.tif 000158
【図２３】
0005400743.tif 000159
【図２４】
0005400743.tif 000160
【図２５】
0005400743.tif 000161
【図２６】
0005400743.tif 000162
【図２７】
0005400743.tif 000163
【図２８】
0005400743.tif 000164
【図２９】
0005400743.tif 000165
【図３０】
0005400743.tif 000166
【図３１】
0005400743.tif 000167
【図３２】
0005400743.tif 000168
【図３３】
0005400743.tif 000169
【図３４】
0005400743.tif 000170
【図３５】
0005400743.tif 000171
【図３６】
0005400743.tif 000172
【図３７】
0005400743.tif 000173
【図３８】
0005400743.tif 000174
【図３９】
0005400743.tif 000175
【図４０】
0005400743.tif 000176
【図４１】
0005400743.tif 000177
【図４２】
0005400743.tif 000178
【図４３】
0005400743.tif 000179
【図４４】
0005400743.tif 000180
【図４５】
0005400743.tif 000181
【図４６】
0005400743.tif 000182
【図４７】
0005400743.tif 000183
【図４８】
0005400743.tif 000184