以下、図面を参照しながら、本発明の第1実施形態に係る制御装置について説明する。本実施形態の図1に示す制御装置1は、後述する制御アルゴリズムによって、内燃機関(以下「エンジン」という)3の混合気の空燃比を制御するものであり、ECU2を備えている。
エンジン3は、図示しない車両に搭載された直噴式のガソリンエンジンであり、このエンジン3には、燃料噴射弁4(1つのみ図示)が気筒ごとに取り付けられている。この燃料噴射弁4は、ECU2に電気的に接続されており、その開弁時間および開弁タイミングがECU2によって制御され、それにより、燃料噴射制御が実行される。この場合、通常運転状態では、空燃比がリーン側の値になるように、燃料噴射制御が実行され、それによって、エンジン3はリーンバーン運転される。
ECU2には、クランク角センサ20およびアクセル開度センサ21がそれぞれ接続されている。このクランク角センサ20(参照パラメータ検出手段)は、マグネットロータおよびMREピックアップで構成されており、クランクシャフト(図示せず)の回転に伴い、いずれもパルス信号であるCRK信号およびTDC信号をECU2に出力する。
このCRK信号は、所定クランク角(例えば1゜)毎に1パルスが出力され、ECU2は、このCRK信号に基づき、エンジン3の機関回転数(以下「エンジン回転数」という)NEを算出する。また、TDC信号は、各気筒のピストン(いずれも図示せず)が吸気行程のTDC位置よりも若干、手前の所定のクランク角位置にあることを表す信号であり、所定クランク角毎に1パルスが出力される。
また、アクセル開度センサ21は、車両の図示しないアクセルペダルの踏み込み量(以下「アクセル開度」という)APを検出して、それを表す検出信号をECU2に出力する。
一方、エンジン3の吸気通路5には、上流側から順に、スロットル弁機構6および吸気圧センサ22が設けられている。このスロットル弁機構6は、スロットル弁6aおよびこれを開閉駆動するTHアクチュエータ6bなどを備えている。スロットル弁6aは、吸気通路5の途中に回動自在に設けられており、当該回動に伴う開度の変化によりスロットル弁6aを通過する空気量を変化させる。THアクチュエータ6bは、ECU2に接続されたモータにギヤ機構(いずれも図示せず)を組み合わせたものであり、ECU2からの制御入力信号によって制御されることにより、スロットル弁6aの開度を変化させる。
さらに、吸気圧センサ22(参照パラメータ検出手段)は、吸気通路5のスロットル弁6aよりも下流側のサージタンクの部分に配置されており、吸気通路5内の圧力(以下「吸気圧」という)を検出して、それを表す検出信号をECU2に出力する。ECU2は、この吸気圧センサ22の検出信号に基づいて、吸気圧PBを算出する。なお、この吸気圧PBは絶対圧として算出される。
一方、エンジン3の排気通路10には、上流側から順に、LAFセンサ23、上流側三元触媒11、酸素濃度センサ24、下流側三元触媒12、ユリア噴射弁13、上流側選択還元触媒14、NH3濃度センサ25および下流側選択還元触媒15が設けられている。
このLAFセンサ23は、ジルコニアおよび白金電極などで構成され、理論空燃比よりもリッチなリッチ領域から極リーン領域までの広範囲な空燃比の領域において、排気通路10内を流れる排ガス中の酸素濃度をリニアに検出し、それを表す検出信号をECU2に出力する。ECU2は、このLAFセンサ23の検出信号の値に基づき、排ガス中の当量比を表す検出当量比KACTを算出する。なお、本実施形態では、検出当量比KACTが制御量および空燃比を表す値に相当する。
また、上流側三元触媒11は、その温度が所定の活性化温度よりも高い領域にあるときに活性化し、排ガス中の有害な未燃成分を浄化する。さらに、下流側三元触媒12は、上流側三元触媒11と同じタイプのものであり、上流側選択還元触媒14に流入する排ガス成分をNOx浄化に最適な状態に調整し、上流側選択還元触媒14において高いNOx浄化率を確保することを目的として、上流側選択還元触媒14の上流側に設けられている。なお、下流側三元触媒12として、上流側三元触媒11と異なるタイプのもの、例えば、リーンバーン運転中のHCやCOの酸化能力をより高めたものや、NOを酸化してNO2化する能力をより高めたものを用いてもよい。
さらに、酸素濃度センサ24は、ジルコニアおよび白金電極などで構成され、上流側三元触媒11を通過した排ガス中の酸素濃度に基づいた出力をECU2に送る。この酸素濃度センサ24の出力は、理論空燃比よりもリッチな混合気が燃焼したときには、ハイレベルの電圧値(例えば0.8V)となり、混合気がリーンのときには、ローレベルの電圧値(例えば0.2V)となるとともに、混合気が理論空燃比付近のときには、ハイレベルとローレベルの間の所定の目標値(例えば0.6V)となる。
一方、ユリア噴射弁13は、ECU2に電気的に接続されており、ECU2からの制御入力信号によって開弁状態に駆動されると、ユリアタンク(図示せず)からの尿素水を排気通路10内に噴射する。この場合、ユリア噴射弁13から噴射された尿素水の尿素は、その一部が、排ガスの熱および選択還元触媒14との接触によって、アンモニアに変化する。
また、上流側選択還元触媒14は、還元剤としての尿素(Urea)が存在する雰囲気下で、排ガス中の窒素酸化物(NOx)を選択的に還元する。上流側選択還元触媒14では、そのNOxの還元作用において、尿素水の噴射時に尿素から変化したアンモニアも、尿素と一緒に消費されるとともに、消費されなかった分のアンモニアは、上流側選択還元触媒14内に貯蔵される。
さらに、下流側選択還元触媒15は、上流側選択還元触媒14と同じタイプのものであり、排ガス中のNOxを浄化することに加えて、上流側選択還元触媒15を通過したアンモニアを捕捉することを目的として、上流側選択還元触媒14の下流側に設けられている。本実施形態の場合、以上のユリア噴射弁13および選択還元触媒14,15によって、尿素SCR装置が構成されている。なお、下流側選択還元触媒15として、上流側選択還元触媒14と比べて低温でのNOx浄化性能をより高めたもの、例えば、Cuゼオライトタイプのものや、後側に酸化触媒をゾーンコート化したものを用いてもよい。
さらに、NH3濃度センサ25は、上流側選択還元触媒14を通過した排ガス中のアンモニア濃度を検出し、それを表す検出信号をECU2に出力する。ECU2は、このNH3濃度センサ25の検出信号に基づき、ユリア噴射弁13を介してユリア噴射量を制御し、それにより、尿素SCR装置によるNOx浄化率やNOx浄化量を制御する。
一方、ECU2は、CPU、RAM、ROM、I/Oインターフェースおよび駆動回路(いずれも図示せず)などからなるマイクロコンピュータで構成されており、前述した各種のセンサ20〜25の検出信号などに応じて、エンジン3の運転状態を判別するとともに、後述する空燃比制御処理などを実行する。
なお、本実施形態では、ECU2が、目標制御量設定手段、参照パラメータ検出手段、予測値算出手段、重み関数値算出手段、予測制御量設定手段、制御入力算出手段、修正後制御入力設定手段、同定手段、および外乱推定値算出手段に相当する。
次に、本実施形態の制御装置1について説明する。まず、本実施形態の制御装置1で用いる制御対象モデルについて説明する。エンジン3の燃料噴射弁4からLAFセンサ23までの系を、空燃比補正係数KAFを制御入力とし、検出当量比KACTを制御量とする一次遅れ系の制御対象と見なしてモデリングすると、下式(1)が得られる。この場合、空燃比補正係数KAFは、後述する制御アルゴリズムを用いて、当量比と同じ次元の値として算出される。
この式(1)のαはモデルパラメータである。また、式(1)において、記号(k)付きの各離散データは、所定の制御周期ΔT(本実施形態ではTDC信号の発生周期)でサンプリングまたは算出されたデータであることを示しており、記号k(kは正の整数)は各離散データのサンプリングまたは算出サイクルの順番を表している。この点は、以下の離散データにおいても同様である。なお、以下の説明では、各離散データにおける記号(k)を適宜省略する。
上式(1)の場合、空燃比補正係数KAFと検出当量比KACTとの間に存在するむだ時間dが考慮されていないので、このむだ時間dを上式(1)に反映させると、下式(2)が得られる。この式(2)を制御対象モデルとして用いた理由については後述する。
この場合、上式(2)のむだ時間dは、エンジン3の運転状態に応じて変化するものであり、このむだ時間dと排ガスボリュームVexとの関係をモデリング(マッピング)すると、図2に示すモデル(マップ)が得られる。この排ガスボリュームVex(参照パラメータ)は、排ガスの空間速度に相当する値であり、具体的には、エンジン回転数NEおよび吸気圧PBに応じて、図示しないマップを検索することにより算出される。
図2において、Vex1〜Vex4,VexMAXは、排ガスボリュームVexの所定値であり、0<Vex1<Vex2<Vex3<Vex4<VexMAXが成立するように設定されている。また、所定値VexMAXは、エンジン3の運転中に排ガスボリュームVexが変化し得る領域の最大値に設定されている。言い換えれば、排ガスボリュームVexは、エンジン3の運転中、0〜VexMAXの領域内で変化する特性を有している。
本実施形態の制御装置1では、以上のむだ時間dを含む上式(2)の制御対象モデルを用いて、以下に述べるように、空燃比補正係数KAFなどの各種の演算値が算出される。図3に示すように、制御装置1は、目標当量比算出部30、可変むだ時間状態予測器(以下「状態予測器」という)40、オンボード・スケジュールド・モデルパラメータ同定器(以下「オンボード同定器」という)60、および周波数整形コントローラ130を備えており、これらはいずれもECU2によって構成されている。
この目標当量比算出部30では、前述した検出当量比KACTの目標となる値として、目標当量比KCMDが算出される。具体的には、エンジン回転数NEおよびアクセル開度APに応じて、図示しないマップを検索することにより、要求トルクTRQDRVを算出し、この要求トルクTRQDRVおよびエンジン回転数NEに応じて、図4に示すマップを検索することにより、目標当量比KCMDが算出される。同図において、KCMD1〜KCMD4は、目標当量比KCMDの所定値であり、KCMD1=1,KCMD1>KCMD2>KCMD3>KCMD4が成立するように設定される。
また、状態予測器40では、後述する予測アルゴリズムを用いて、検出当量比KACTの予測値として、予測当量比PRE_KACTが算出され、オンボード同定器60では、後述する同定アルゴリズムを用いて、前述したモデルパラメータαをオンボード同定した値として、同定値αidが算出される。さらに、周波数整形コントローラ130では、後述する制御アルゴリズムを用いて、制御入力としての空燃比補正係数KAFが算出される。
なお、本実施形態では、目標当量比算出部30が目標制御量設定手段に相当し、目標当量比KCMDが目標制御量に相当する。また、状態予測器40が、予測値算出手段、重み関数値算出手段および予測制御量設定手段に相当し、予測当量比PRE_KACTが予測制御量に相当する。さらに、オンボード同定器60が、修正後制御入力設定手段、同定手段および重み関数値算出手段に相当し、周波数整形コントローラ130が制御入力算出手段に相当する。
次に、前述した状態予測器40について説明する。この状態予測器40は、以下に述べる予測アルゴリズムにより、予測当量比PRE_KACTを算出するものであり、この予測当量比PRE_KACTは、現在の制御系におけるむだ時間dが経過した後の制御タイミングにおける検出当量比KACTを予測した値に相当する。
図5に示すように、状態予測器40は、3つの遅延要素41〜43と、増幅器44と、3つの予測値算出部45〜47と、4つの重み関数値算出部48〜51と、4つの乗算器52〜55と、加算器56とを備えている。
まず、増幅器44で、下式(3)により、第0予測値PRE_KACT_0が算出される。すなわち、第0予測値PRE_KACT_0は、むだ時間d=0のときの検出当量比KACT(k)として算出される。
また、第1予測値算出部45では、遅延要素41で1制御サイクル分遅延された空燃比補正係数の値KAF(k−1)を用い、下式(4)により、第1予測値PRE_KACT_1が算出される。
この第1予測値PRE_KACT_1は、d=1のときのむだ時間dが経過したときの検出当量比KACTを予測した値に相当する。なお、上式(4)の導出手法については後述する。
さらに、第2予測値算出部46では、2つの遅延要素41,42でそれぞれ1,2制御サイクル分ずつ、遅延された空燃比補正係数の値KAF(k−1),KAF(k−2)を用い、下式(5)により、第2予測値PRE_KACT_2が算出される。
この第2予測値PRE_KACT_2は、d=2のときのむだ時間dが経過したときの検出当量比KACTを予測した値に相当する。なお、上式(5)の導出手法については後述する。
また、第3予測値算出部47では、3つの遅延要素41〜43でそれぞれ1〜3制御サイクル分ずつ、遅延された空燃比補正係数の値KAF(k−1),KAF(k−2),KAF(k−3)を用い、下式(6)により、第3予測値PRE_KACT_3が算出される。
この第3予測値PRE_KACT_3は、d=3のときのむだ時間dが経過したときの検出当量比KACTを予測した値に相当する。なお、上式(6)の導出手法については後述する。
さらに、4つの重み関数値算出部48〜51では、排ガスボリュームVexに応じて、図6に示すマップを検索することにより、4つの重み関数値Wd1〜Wd4がそれぞれ算出される。同図に示すように、4つの重み関数値Wd1〜Wd4はそれぞれ、排ガスボリュームVexが変化し得る領域を0≦Vex≦Vex2,Vex1≦Vex≦Vex3,Vex2≦Vex≦Vex4,Vex3≦Vex≦VexMAXの4つの領域に区分した場合において、これらの4つの領域に対応するように設定されているとともに、対応する領域では、値1以下の正の値に設定され、対応する領域以外では、値0に設定されている。
具体的には、重み関数値Wd1は、これが対応する領域(0≦Vex≦Vex2)では、Vex≦Vex1のときの値1を最大値として、Vex1<Vexの領域で排ガスボリュームVexが大きいほど、より小さい正の値に設定されているとともに、それ以外の領域では、値0に設定されている。また、重み関数値Wd2は、これが対応する領域(Vex1≦Vex≦Vex3)では、Vex=Vex2のときの値1を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値0に設定されている。
さらに、重み関数値Wd3は、これが対応する領域(Vex2≦Vex≦Vex4)では、Vex=Vex3のときの値1を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値0に設定されている。一方、重み関数値Wd4は、これが対応する領域(Vex3≦Vex≦VexMAX)では、Vex4≦Vexのときの値1を最大値として、排ガスボリュームVexが大きいほど、より大きい正の値に設定されているとともに、それ以外の領域では、値0に設定されている。
以上に加えて、4つの重み関数値Wdi(i=1〜4)がそれぞれ対応する4つの領域は、上述したような、隣り合う領域が互いに重なり合うように設定されており、これらの互いに重なり合う領域における排ガスボリュームVexのいずれかの値に対応する重み関数値Wdiの値の和は、各重み関数値Wdiにおける最大値1と等しくなるように設定されている。
さらに、図6と前述した図2を比較すると明らかなように、互いに重なり合う3つの領域は、むだ時間dの勾配が一定状態に保持される3つの領域に対応するように設定されている。これに加えて、重み関数値Wd1はむだ時間d=3に対して、重み関数値Wd2はむだ時間d=2に対して、重み関数値Wd3はむだ時間d=1に対して、重み関数値Wd4はむだ時間d=0に対して重みが最も大きくなるようにそれぞれ設定されている。
また、乗算器52では、重み関数値Wd4に第0予測値PRE_KACT_0を乗算することにより、乗算値Wd4・PRE_KACT_0が算出され、乗算器53では、重み関数値Wd3に第1予測値PRE_KACT_1を乗算することにより、乗算値Wd3・PRE_KACT_1が算出される。さらに、乗算器53では、重み関数値Wd2に第2予測値PRE_KACT_2を乗算することにより、乗算値Wd2・PRE_KACT_2が算出され、乗算器54では、重み関数値Wd1に第3予測値PRE_KACT_3を乗算することにより、乗算値Wd1・PRE_KACT_3が算出される。
そして、加算器56で、以上の4つの乗算値を互いに加算することにより、予測当量比PRE_KACTが算出される。すなわち、下式(7)により、予測当量比PRE_KACが算出される。
以上のように、前述した4つの重み関数値Wdiを、4つの予測値PRE_KACT_4−iにそれぞれ乗算した値の総和として、予測当量比PRE_KACTが算出されるので、排ガスボリュームVexの変化に応じて、むだ時間dが図2に示すように値0から値3の間で連続的に変化したときでも、そのようなむだ時間dの変化を適切に反映させながら、予測当量比PRE_KACTを円滑に段差なく変化するような値として算出することができる。
なお、前述した第1〜3予測値PRE_KACT_1〜3の算出式(4)〜(6)は、以下に述べるように導出される。まず、前述した式(2)において、d=1とすると、下式(8)が得られる。
上式(8)において、右辺のKACT(k+1)をPRE_KACT_1(k)に、左辺のαをαid(k)にそれぞれ置き換えると、前述した式(4)が得られる。
また、前述した式(2)において、d=2とすると、下式(9)が得られる。
上式(9)の各変数を1制御サイクル分、未来側にシフトさせると、下式(10)が得られる。
上式(10)に上式(9)を代入すると、下式(11)が得られる。
上式(11)において、右辺のKACT(k+2)をPRE_KACT_2(k)に、左辺のαをαid(k)にそれぞれ置き換えると、前述した式(5)が得られる。
また、前述した式(2)において、d=3とすると、下式(12)が得られる。
上式(12)の各変数を1制御サイクル分、未来側にシフトさせると、下式(13)が得られる。
上式(13)に上式(12)を代入すると、下式(14)が得られる。
さらに、上式(13)の各変数を1制御サイクル分、未来側にシフトさせると、下式(15)が得られる。
この式(15)に上式(14)を代入すると、下式(16)が得られる。
この式(16)において、右辺のKACT(k+3)をPRE_KACT_3(k)に、左辺のαをαid(k)にそれぞれ置き換えると、前述した式(6)が得られる。
次に、前述したオンボード同定器60について説明する。このオンボード同定器60の場合、本実施形態の制御対象のように、むだ時間dが排ガスボリュームVexに応じて連続的に変化する場合において、以下に述べる拘束条件付きスケジュールド補正型の同定アルゴリズムを用いて、そのようなむだ時間dの変化を反映させながら、同定値αidが算出される。なお、このオンボード同定器60の同定アルゴリズムは、後述するように、前述した式(2)の右辺の値KAF(k−d)を後述する修正後制御入力KAF_mod(k)に置き換えた修正後モデル(後述する式(30))に基づいて導出される。
このオンボード同定器60は、図7に示すように、修正後制御入力算出部70、3つの遅延要素61〜63、合成信号値算出部64、推定合成信号値算出部65、同定ゲイン算出部66、減算器67、乗算器68および同定値算出部90を備えている。
まず、修正後制御入力算出部70について説明する。この修正後制御入力算出部70は、修正後制御入力KAF_modを算出するものであり、図8に示すように、修正後制御入力算出部70は、3つの遅延要素71〜73、4つの重み関数値算出部74〜77、4つの乗算器78〜81および加算器82を備えている。
まず、4つの重み関数値算出部74〜77で、前述した4つの重み関数値算出部48〜51と同様に、排ガスボリュームVexに応じて、前述した図6に示すマップを検索することにより、4つの重み関数値Wd1〜Wd4がそれぞれ算出される。
また、乗算器78では、空燃比補正係数の今回値KAF(k)を重み関数値Wd4(k)に乗算することにより、乗算値Wd4(k)・KAF(k)が算出され、乗算器79では、遅延要素71で1制御サイクル分遅延された空燃比補正係数の値KAF(k−1)を重み関数値Wd3(k)に乗算することにより、乗算値Wd3(k)・KAF(k−1)が算出される。
さらに、乗算器80では、2つの遅延要素71,72によって2制御サイクル分遅延された空燃比補正係数の値KAF(k−2)を重み関数値Wd2(k)に乗算することにより、乗算値Wd2(k)・KAF(k−2)が算出され、乗算器81では、3つの遅延要素71〜73によって3制御サイクル分遅延された空燃比補正係数の値KAF(k−3)を重み関数値Wd1(k)に乗算することにより、乗算値Wd1(k)・KAF(k−3)が算出される。
そして、加算器82で、以上の4つの乗算値を用い、下式(17)により、修正後制御入力KAF_modが算出される。
図7に戻り、合成信号値算出部64では、検出当量比KACTと、遅延要素61によって1制御サイクル分遅延された検出当量比の値KACT(k−1)を用い、下式(18)により、合成信号値W_actが算出される。
また、推定合成信号値算出部65では、遅延要素61によって1制御サイクル分遅延された検出当量比の値KACT(k−1)と、遅延要素62によって1制御サイクル分遅延された修正後制御入力の値KAF_mod(k−1)とを用いて、下式(19)により、偏差ζ’(k−1)を算出した後、これと遅延要素63によって1制御サイクル分遅延された同定値の値αid(k−1)を用い、下式(20)により、推定合成信号値W_hatが算出される。
さらに、減算器67では、下式(21)により、同定誤差eid’が算出される。
一方、同定ゲイン算出部68では、下式(22),(23)により、同定ゲインKp’が算出される。この同定ゲインKp’は、同定値αidの修正方向(正負)および修正量を規定するものである。
上式(22)のゲインP’(k)の初期値P’(0)は、下式(24)のように定義される。
ここで、P0は所定値に設定されている。
また、上式(22)において、λ1,λ2は重みパラメータであり、これらの値λ1、λ2を下記のように設定することにより、同定アルゴリズムとして、3つのアルゴリズムのいずれかを選択することができる。
λ1=1,λ2=0 ;固定ゲインアルゴリズム
λ1=1,λ2=1 ;最小2乗法アルゴリズム
λ1=λ,λ2=1 ;重み付き最小2乗法アルゴリズム
ここで、λは0<λ<1に設定される所定値である。本実施形態の場合、同定精度および制御精度を適切に確保するために、重み付き最小2乗法アルゴリズムを用いている。
さらに、乗算器68で、同定ゲインKp’と同定誤差eid’の乗算値Kp’・eid’が算出される。
次いで、同定値算出部90では、上記の乗算値Kp’・eid’と排ガスボリュームVexを用いて、以下に述べるように同定値αidが算出される。図9に示すように、この同定値算出部90は、基準モデルパラメータ算出部91、4つの重み関数値算出部92〜95、8つの乗算器96〜103、5つの加算器104〜108、4つの遅延要素109〜112および4つの増幅器113〜116を備えている。
まず、基準モデルパラメータ算出部91では、排ガスボリュームVexに応じて、図10に示すマップを検索することにより、基準モデルパラメータαbsが算出される。同図において、Vex5〜Vex8は、排ガスボリュームVexの所定値であり、0<Vex5<Vex6<Vex7<Vex8<VexMAXが成立するように設定されている。このマップでは、基準モデルパラメータαbsは、排ガスボリュームVexが大きいほど、より大きい値に設定されている。これは、排ガスボリュームVexが大きいほど、LAFセンサ23におけるセンサカバーの孔を介した排ガスの交換が促進され、LAFセンサ23の遅れ特性が小さくなることで、空燃比補正係数KAFが検出当量比KACTに及ぼす影響の度合がより大きくなることによる。
また、4つの重み関数値算出部92〜95では、排ガスボリュームVexに応じて、図11に示すマップを検索することにより、4つの重み関数値Wa1〜Wa4がそれぞれ算出される。同図に示すように、4つの重み関数値Wa1〜Wa4はそれぞれ、排ガスボリュームVexが変化し得る領域を0≦Vex≦Vex6,Vex5≦Vex≦Vex7,Vex6≦Vex≦Vex8,Vex7≦Vex≦VexMAXの4つの領域に区分した場合において、これらの4つの領域に対応するように設定されているとともに、対応する領域では、値1以下の正の値に設定され、対応する領域以外では、値0に設定されている。
すなわち、重み関数値Wa1は、これが対応する領域(0≦Vex≦Vex6)では、Vex≦Vex5のときの値1を最大値として、排ガスボリュームVexが大きいほど、より小さい正の値に設定されているとともに、それ以外の領域では、値0に設定されている。また、重み関数値Wa2は、これが対応する領域(Vex5≦Vex≦Vex7)では、Vex=Vex6のときの値1を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値0に設定されている。
さらに、重み関数値Wa3は、これが対応する領域(Vex6≦Vex≦Vex8)では、Vex=Vex7のときの値1を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値0に設定されている。一方、重み関数値Wa4は、これが対応する領域(Vex7≦Vex≦VexMAX)では、Vex8≦Vexのときの値1を最大値として、排ガスボリュームVexが大きいほど、より大きい正の値に設定されているとともに、それ以外の領域では、値0に設定されている。
以上に加えて、4つの重み関数値Wai(i=1〜4)がそれぞれ対応する4つの領域は、上述したような、隣り合う領域が互いに重なり合うように設定されており、これらの互いに重なり合う領域における排ガスボリュームVexのいずれかの値に対応する重み関数値Waiの値の和は、各重み関数値Waiにおける最大値1と等しくなるように設定されている。図11と前述した図10を比較すると明らかなように、これらの互いに重なり合う3つの領域は、基準モデルパラメータαbsの勾配が一定状態に保持される3つの領域に対応するように設定されている。
また、乗算器96で、重み関数値Wa1に値Kp’・eid’を乗算することにより、乗算値Wa1・Kp’・eid’が算出され、遅延要素109で1制御サイクル分遅延された修正項dα1(k−1)に、増幅器113でゲイン係数H(k)を乗算することにより、乗算値H(k)・dα1(k−1)が算出される。なお、このゲイン係数Hについては後述する。そして、加算器104で、値Wa1・Kp’・eid’に値H(k)・dα1(k−1)を加算することにより、修正項dα1が算出される。
さらに、乗算器97で、重み関数値Wa2に値Kp’・eid’を乗算することにより、乗算値Wa2・Kp’・eid’が算出され、遅延要素110で1制御サイクル分遅延された修正項dα2(k−1)に、増幅器114でゲイン係数H(k)を乗算することにより、乗算値H(k)・dα2(k−1)が算出される。そして、加算器105で、値Wa2・Kp’・eid’に値H(k)・dα2(k−1)を加算することにより、修正項dα2が算出される。
また、乗算器98で、重み関数値Wa3に値Kp’・eid’を乗算することにより、乗算値Wa3・Kp’・eid’が算出され、遅延要素111で1制御サイクル分遅延された修正項dα3(k−1)に、増幅器115でゲイン係数H(k)を乗算することにより、乗算値H(k)・dα3(k−1)が算出される。そして、加算器106で、値Wa3・Kp’・eid’に値H(k)・dα3(k−1)を加算することにより、修正項dα3が算出される。
さらに、乗算器99で、重み関数値Wa4に値Kp’・eid’を乗算することにより、乗算値Wa4・Kp’・eid’が算出され、遅延要素112で1制御サイクル分遅延された修正項dα4(k−1)に、増幅器116でゲイン係数H(k)を乗算することにより、乗算値H(k)・dα4(k−1)が算出される。そして、加算器107で、値Wa4・Kp’・eid’に値H(k)・dα4(k−1)を加算することにより、修正項dα4が算出される。
以上の増幅器113〜116では、ゲイン係数Hが以下の式(25)〜(27)に示すように算出される。
上式(25)〜(27)において、α_Lは所定の下限値であり、α_Hは所定の上限値である。また、η’は忘却係数であり、0<η’≦1が成立するように設定される。同定値αidの算出において、この忘却係数η’を用いた理由は、エンジン3の定常運転状態が長時間継続した場合、同定値αidが増大化し、不適切な値となるおそれがあるので、それを回避するためである。さらに、上式(26)に示すように、同定値αidが下限値α_Lと上限値α_Hとの間の領域にあるときに、忘却係数η’による忘却効果を停止している理由は、このオンボード同定器60で用いている同定アルゴリズムの場合、後述する同定条件1(拘束条件)を満たすように、同定値αidを常に同定することができるので、拘束条件を満たすことを目的として、同定値αidを後述する基準モデルパラメータαbsの近傍に無理に拘束する必要がないことによる。
また、以上の4つの加算器104〜107における演算は、下式(28)で表される。
さらに、乗算器100〜103ではそれぞれ、4つの重み関数値Waiを4つの修正項dαiに乗算することにより、4つの乗算値Wai・dαiが算出される。
そして、加算器108で、下式(29)により、同定値αidが最終的に算出される。
以上のように、オンボード同定器60において、修正後制御入力KAF_modが、4つの制御タイミングでの検出当量比KACTに4つの重み関数値Wdiをそれぞれ乗算した値の総和として算出され、4つの修正項dαiが、この修正後制御入力KAF_modを用いて算出した同定誤差eid’と同定ゲインKp’との積Kp’・eid’に4つの重み関数値Waiをそれぞれ乗算した値の総和として算出される。そして、その総和を基準モデルパラメータαbsに加算することにより、同定値αidが算出されるので、排ガスボリュームVexの変化に応じて遅れ特性やむだ時間dが連続的に変化したときでも、2種類の重み関数値Wdi,Waiの効果によって、その影響を抑制しながら、同定値αidを円滑に変化するような値として同定することができる。
同定値αidの算出において、以上の式(17)〜(29)の同定アルゴリズムを用いた理由は以下による。まず、本実施形態の制御装置1における制御系は、空燃比補正係数KAFを制御入力とし、検出当量比KACTを制御量とするものであり、外乱がない状態では、定常偏差を生じない系となる。それに起因して、前述した式(2)の制御対象モデルの場合、入力と出力との間に定常偏差を生じないようにするために、入力項および出力項の乗算係数すなわちモデルパラメータα,1−αは、両者の和が値1になるように設定されている。
この場合、2つのモデルパラメータα,1−αは、互いに独立した値を取ることができず、一方の増大に従って他方が減少するという互いに拘束し合う関係にある。したがって、これらのモデルパラメータα,1−αを同定する場合、一方の増大に従って他方が減少するという互いに拘束し合う条件(以下「拘束条件」という)を満たすように、2つのモデルパラメータα,1−αを同定する必要があり、以下、これを「同定条件1」という。ここで、最小2乗法などの一般的な同定アルゴリズムをそのまま用いた場合、この同定条件1を満たすのは困難である。
これに加えて、前述したように、遅れ特性やむだ時間dが排ガスボリュームVexに応じて変化する特性を有しているので、最小2乗法などの一般的な同定アルゴリズムをそのまま用いた場合、遅れ特性やむだ時間dの変化を反映させながら、2つのモデルパラメータα,1−αを同定することができず、結果的に、2つのモデルパラメータα,1−αの同定精度が低下してしまう。したがって、同定精度を高めるために、遅れ特性やむだ時間dが変化した場合でも、それを適切に反映させながら、モデルパラメータα,1−αを同定する必要があり、以下、これを「同定条件2」という。
まず、上記の同定条件2を満たすために、前述した式(2)に代えて、下式(30)を制御対象モデルとして用いる。
この式(30)は、前述した式(2)の右辺の値KAF(k−d)を値KAF_mod(k)に置き換えたものに相当する。この修正後制御入力KAF_modは、前述した式(17)に示すように、4つの重み関数値Wdiと4つの空燃比補正係数KAFとの乗算値の和として算出され、4つの重み関数値Wdiは、前述した手法により算出されるので、修正後制御入力KAF_modを、むだ時間dが変化したときでも、それを適切に反映させながら算出することができる。これに加えて、重み関数値Waiを用いることによって、遅れ特性の変化を反映させながら、4つの修正項dαiを算出することができる。すなわち、上記の同定条件2を満たすことができる。
上式(30)を変形すると、下式(31)が得られる。
上式(31)の左辺および右辺をそれぞれ、下式(32),(33)に示すように、合成信号値W_actおよび推定合成信号値W_hatとして定義する。
このように定義した場合、上記の同定条件1を満たすには、合成信号値W_actおよび推定合成信号値W_hatが互いに一致するように、制御対象モデルのモデルパラメータを同定すればよいことになる。すなわち、前述した同定誤差eid’が値0になるように、同定値αidを同定(算出)すればよいことになる。以上の理由により、前述した式(17)〜(29)の同定アルゴリズムを用いて、同定値αidが算出される。
また、制御対象モデルのモデルパラメータαと排ガスボリュームVexとが前述した図10の関係を有している場合、最小2乗法などの一般的な同定アルゴリズムでは、図10の関係を基準としてモデルパラメータを同定することができない。これに対して、本実施形態のオンボード同定器60の場合、排ガスボリュームVexに応じて図10のマップを検索することによって、基準モデルパラメータαbsを算出し、この基準モデルパラメータαbsを、前述した重み関数値Waiと修正項dαiとの積の総和で修正することにより、同定値αidが算出されるので、高い同定精度を確保することができる。
次に、前述した周波数整形コントローラ130について説明する。この周波数整形コントローラ130は、予測当量比PRE_KACTが目標当量比KCMDに収束するように、空燃比補正係数KAFを算出するもの、言い換えれば、検出当量比KACTが目標当量比KCMDに収束するように、空燃比補正係数KAFを算出するものである。この周波数整形コントローラ130では、まず、下式(34)に示すように、前述した予測当量比PRE_KACTから目標当量比KCMDを減算することにより、予測追従誤差PRE_eが算出される。
次いで、下式(35)により、制御入力としての空燃比補正係数KAFが算出される。
この式(35)において、βは感度設定パラメータであり、以下に述べる手法により所定値(例えば値0.6)に設定されている。
次に、上記周波数整形コントローラ130の制御アルゴリズムの導出原理について説明する。本実施形態の場合、良好な排ガス特性と良好な燃費とを両立させることを目的として、ガソリンエンジンタイプのエンジン3において、その空燃比をリーン側に制御し、リーンバーン運転するとともに、尿素SCR装置によって排ガス中のNOxを浄化するように構成されている。
このように構成した場合、ガソリンエンジンでは、リーンバーン運転中の燃焼安定性が低く、燃焼可能な混合気の空燃比が所定範囲内に制限されることになるので、空燃比が一時的に過度にリーン化する現象を抑制する必要がある。この現象は、過渡運転状条件下で特に発生しやすい。これに加えて、リーンバーン運転中は、燃焼変動に起因して、サージング現象が発生しやすくなるので、それを抑制するために、燃料量が過度に変動しないように制御する必要がある。これらの要求を満たすには、低周波数の外乱に対する抑制能力が低く、かつ高周波数の外乱に対する抑制能力が高くなるように、空燃比を制御することが必要となる。以下、これを「制御条件φ」という。
ここで、本実施形態の制御装置1のようなフィードバック制御系の構成、すなわち制御入力としての空燃比補正係数KAFを制御対象に入力することにより、検出当量比KACTを目標当量比KCMDに収束させるようにフィードバック制御する系の構成を、z領域のブロック図で表すと、図12のようになる。同図において、C(z)はコントローラの伝達関数を、P(z)は制御対象の伝達関数を、D(z)は外乱をそれぞれ表している。なお、以下の説明では、各データにおける記号(z)を適宜省略する。
この制御系の場合、外乱Dと検出当量比KACTとの間の伝達関数すなわち感度関数Sは、下式(36)のようになる。
この場合、上記制御条件φを満たすには、感度関数Sにおいて、そのゲイン特性(すなわち周波数応答特性)を表すゲイン曲線が図13に示すものになることが必要になる。以下、制御条件φを満たす図13のゲイン曲線が得られる感度関数を「最適感度関数Sopt」という。同図において、FQ1は、所定周波数を表しており、この所定周波数FQ1は予め実験により設定される。同図に示すように、この最適感度関数Soptの場合、所定周波数FQ1以上の、外乱を抑制する必要性が高い周波数領域(以下「外乱抑制領域」という)では、高ゲインに設定されている一方、所定周波数FQ1未満の、外乱を抑制する必要性が低い周波数域(以下「外乱非抑制領域」という)では、外乱抑制領域よりも低いゲインに設定されている。より具体的には、最適感度関数Soptは、外乱抑制領域では、高ゲインかつフラットな状態のゲイン特性を示すとともに、外乱非抑制領域では、周波数が低いほど、ゲインが連続的に急減するような特性を示すように設定されている。なお、本実施形態では、制御条件φを満たすように設定された最適感度関数Soptが、所定の周波数特性が得られるように設定された感度関数に相当する。
ここで、前述した特許文献1におけるスライディングモード制御アルゴリズムを図12の制御系に適用した場合、以下のようになる。すなわち、スライディングモード制御アルゴリズムの場合、追従誤差eおよび切換関数σが以下の式(37),(38)のように定義される。
ここで、POLE_Eは、−1<POLE_E<0が成立するように設定される切換関数設定パラメータである。
スライディングモード制御アルゴリズムの場合、制御対象の動特性をσ=0になるように拘束する制御手法であり、σ=0を上式(38)に適用すると、下式(39)が得られるとともに、この式(39)を整理すると、下式(40)が得られる。
上式(40)は、入力の存在しない一次遅れ系を表している。すなわち、スライディングモード制御アルゴリズムは、制御対象の動特性を入力の存在しない一次遅れ系に拘束する制御アルゴリズムであり、そのような一次遅れ系の感度関数Ssldのゲイン曲線は、図14に実線で示すものとなる。同図を参照すると明らかなように、この感度関数Ssldの場合、そのゲイン曲線は同図に破線で示す最適感度関数Soptのゲイン曲線とかなり近似しており、前述した制御条件φを満たすものとなることが判る。
ここで、スライディングモード制御アルゴリズムの場合、追従誤差eが切換直線上の値となる(すなわちσ=0となる)までの到達モードと、追従誤差eが切換直線上の値となった(すなわち、制御対象の動特性が入力の存在しない一次遅れ系に拘束された)以降のスライディングモードとが存在している。そのため、スライディングモードにあるときには、制御条件φを満たすことができるものの、到達モードにあるときには、制御条件φを満たすことができないことになる。すなわち、スライディングモード制御アルゴリズムの場合、制御条件φを常に満たすことができない。
そこで、本実施形態では、制御条件φを常に満たす制御アルゴリズムとして、以下に述べるように、感度関数Sdを制御条件φを満たすように予め設定する制御アルゴリズムを採用した。まず、当量比の次元を備える空燃比補正係数KAFを制御入力とし、検出当量比KACTを制御量とする系を一次遅れ系と想定すると、その制御対象モデルは前述した式(1)となるとともに、その制御対象モデルのz領域の伝達関数Pは、下式(41)となる。
一方、制御条件φを満たす感度関数Sdを下式(42)に示すように定義する。
上式(42)のβは、感度関数設定パラメータであり、0<β<1となる所定値に設定される。上式(42)の感度関数Sdにおいて、β=0.6と設定したときのゲイン曲線は、図15に示す実線で示すものとなる。同図を参照すると明らかなように、この感度関数Sdの場合、そのゲイン曲線は、同図に破線で示す最適感度関数Soptのゲイン曲線とかなり近似しており、前述した制御条件φを満たすものとなることが判る。
この感度関数Sdと、コントローラの伝達関数Cおよび制御対象の伝達関数Pの関係は、下式(43)に示すものとなる。
上式(43)を変形し、感度関数Sdの定義式をコントローラCについて解くと、下式(44)が得られる。
この式(44)に式(42)を代入すると、下式(45)が得られる。
この式(45)を離散時間系の漸化式で表すと、下式(46)が得られる。
この式(46)を参照すると、制御対象モデルのモデルパラメータαと、感度関数Sdの周波数応答特性(ゲイン特性)を決定する感度関数設定パラメータβとによって、コントローラのフィードバックゲインを指定(設定)できることが判る。
一方、前述した図12の制御系の場合、相補感度関数Tは下式(47)で表される。
ここで、相補感度関数Tと感度関数Sの関係は下式(48)で表されることが知られている。
以上の式(47),(48)から明らかなように、前述した式(46)を導出した手法は、外乱Dと検出当量比KACTとの間の周波数応答特性(ゲイン特性)を決定すると同時に、目標当量比KCMDと検出当量比KACTとの間の周波数応答特性(ゲイン特性)を決定していることになる。
ここで、図15の感度関数Sdに対応する相補感度関数をTdとした場合、この相補感度関数Tdのゲイン曲線は、図16に示すものとなる。なお、同図において、破線で示す曲線は、後述するモデル化誤差Δlを相補感度関数Tdに反映させた場合のゲイン曲線である。
以上のように、式(46)の制御アルゴリズムは、制御条件φを満たす感度関数Sdを用いて導出されたものであるが、この式(46)をそのまま用いて制御を実行しようとした場合、以下に述べる問題点1,2が生じてしまう。
<問題点1>:制御対象モデルのモデルパラメータαの変動やばらつきに対応することができず、高いロバスト性を確保することができない。例えば、従来のPID制御アルゴリズムや最適制御アルゴリズムと同様のロバスト性しか確保することができない。
<問題点2>:制御対象がむだ時間特性を備えている場合、それに対応することができず、制御精度の低下を招いてしまう。
まず、<問題点1>について具体的に説明する。式(1)の制御対象モデルと、実際の制御対象とのモデル式誤差をΔl(z)とした場合、制御系の安定条件として、下式(49)の成立が必要であることが知られている。
ここで、式(1)の一次遅れ系モデルのような遅れ系モデルの場合、そのモデル化誤差Δlは、図17に示すように、周波数域が高くなるほど、より増大するという特性を備えているので、このモデル化誤差Δlを前述した相補感度関数Tdに反映させると、図16の破線で示すゲイン曲線が得られる。上記式(49)から明らかなように、値Td・Δlが0dB未満であることが制御系の安定条件となるので、相補感度関数Tdのゲインが0dBを下回っている度合が制御系の安定余裕となり、ロバスト性を表すことになる。
しかし、感度関数Sdと相補感度関数Tdとの間には、前述したように、Td(z)+Sd(z)=1の関係が存在しているので、外乱抑制に対する周波数応答特性とロバスト性とを互いに独立して設定することはできない。したがって、外乱抑制に対する周波数応答特性を指定した状態で、遅れ系モデルのモデル化誤差Δl(z)に対するロバスト性を向上させるには、このモデル化誤差Δl(z)を補償できるような他の制御アルゴリズムが必要となる。
なお、モデル化誤差Δl(z)に対応するために、式(42)の次数を増やし、感度関数Sdを複雑な形に変更した場合、式(45)の伝達関数C(z)において、分母のzの次数よりも分子のzの次数が大きくなってしまい、実現不可能なコントローラとなってしまう。また、感度設定パラメータβをトライ&エラーでチューニングする手法を用いた場合、PID制御のゲインや、最適制御の重み関数Q,Rをチューニングする手法と変わらず、外乱抑制の周波数応答特性を直接指定する前述した式(46)による制御手法の利点を失ってしまう。
次に、前述した<問題点2>について説明する。本実施形態の制御系の場合、空燃比補正係数KAFと検出当量比KACTとの間にはむだ時間dが存在し、その制御対象モデルとして前述した式(2)を用いているが、前述した式(2)の制御対象のz領域での伝達関数P(z)は、下式(50)となる。
この式(50)の伝達関数P(z)においてd=2に設定したときのボード線図は、図18に示すものとなり、前述した式(41)のむだ時間のない制御系の伝達関数P(z)のボード線図は、図19に示すものとなる。両図を比較すると明らかなように、むだ時間dの有無は、ゲイン特性の差異として現れないので、むだ時間を前述したモデル化誤差Δlとして表すことはできず、前述した式(46)を用いる制御手法、すなわち感度関数Sdや相補感度関数Tdのゲイン特性によってコントローラのゲインを指定する制御手法では、むだ時間に対するロバスト性を考慮し、これを補償することができない。
これに対して、むだ時間が制御系に存在する場合、制御系の安定性が著しく低下することは周知であり、そのため、前述した制御手法をむだ時間を有する制御系に適用すると、制御系が発散してしまうおそれがある。
また、前述した式(44)に、前述した式(42)および式(50)を代入し、コントローラの伝達関数C(z)を導出すると、下式(51)が得られる。
この式(51)を離散時間系の漸化式で表すと、下式(52)が得られる。
この式(52)の場合、追従誤差eの未来値e(k+d),e(k+d−1)が右辺に含まれているので、コントローラの制御アルゴリズムとして、実現不可能なものとなってしまう。
さらに、本実施形態の制御対象の場合、制御入力としての空燃比補正係数KAFと制御量としての検出当量比KACTとの間のむだ時間dは、前述した図2に示すように、排ガスボリュームVexに応じて連続的に変化する特性を有しているので、前述した外乱抑制の周波数応答特性を直接的に指定する制御手法の場合、むだ時間を有する制御対象に適用不可能なものである以上、むだ時間dが変化するような制御系には当然のごとく適用不可能である。
以上のように、前述した問題点1,2を解消すべく、前述した感度関数Sdまたは相補感度関数Tdを用いたコントローラすなわち前述した外乱抑制の周波数応答特性を直接的に指定する制御アルゴリズムを用いながら、制御対象モデルのモデルパラメータαの変動やばらつきに対応し、かつ、制御対象のむだ時間dが変化する特性に対して対応するできるように、制御アルゴリズムを構築する必要がある。
この要求を実現するために、本実施形態の制御装置1の場合、まず、オンボード同定器60で、前述した同定アルゴリズムにより、モデルパラメータαの同定値αidを算出し、状態予測器40において、前述した予測アルゴリズムにより、むだ時間dが経過したタイミングでの検出当量比KACTに相当する予測当量比PRE_KACTを算出する。
そして、周波数整形コントローラ130の制御アルゴリズムとして、検出当量比KACTに代えて予測当量比PRE_KACTを用い、さらに、前述した式(2)のモデルパラメータαを同定値αidに置き換えた下式(53)を制御対象モデルとして用い、前述した式(46)を導出した手法と同じ手法によって、前述した式(34),(35)が導出される。
この式(53)は、前述した式(1)のαをαidに置き換えたものである。言い換えれば、前述した制御対象モデルの式(2)のむだ時間特性を除去したもの(むだ時間特性を考慮しないもの)に相当する。
次に、図20を参照しながら、ECU2によって実行される空燃比制御処理について説明する。この制御処理は、以下に述べるように、燃料噴射弁4から噴射すべき燃料噴射量Toutを算出するものであり、前述した所定の制御周期ΔTで実行される。
この処理では、まず、ステップ1(図では「S1」と略す。以下同じ)で、エンジン回転数NEおよび吸気圧PBに応じて、図示しないマップを検索することにより、基本噴射量TiBSを算出する。
次いで、ステップ2に進み、LAFセンサ正常フラグF_LAFOKが「1」であるか否かを判別する。このLAFセンサ正常フラグF_LAFOKは、図示しない判定処理において、LAFセンサ23が正常であるときに「1」に、それ以外のときに「0」に設定される。
ステップ2の判別結果がNOで、LAFセンサ23が故障しているときには、ステップ14に進み、燃料噴射量Toutを基本噴射量TiBSに設定した後、本処理を終了する。
一方、ステップ2の判別結果がYESで、LAFセンサ23が正常であるときには、ステップ3に進み、三元触媒活性化フラグF_TWCACTが「1」であるか否かを判別する。この三元触媒活性化フラグF_TWCACTは、図示しない判定処理において、2つの三元触媒11,12がいずれも活性化しているときに「1」に、それ以外のときに「0」に設定される。
ステップ3の判別結果がNOで、2つの三元触媒11,12の少なくとも一方が活性化していないときには、ステップ7に進み、目標当量比KCMDを所定のリーン用値KLEARNに設定する。この所定のリーン用値KLEARNは、エンジン始動直後におけるHCの発生を抑制できるような値(例えば0.9)に設定されている。
一方、ステップ3の判別結果がYESで、2つの三元触媒11,12がいずれも活性化しているときには、ステップ4に進み、SCR活性化フラグF_SCRACTが「1」であるか否かを判別する。このSCR活性化フラグF_SCRACTは、図示しない判定処理において、2つの選択還元触媒14,15の少なくとも一方が活性化しているときに「1」に、それ以外のときに「0」に設定される。
ステップ4の判別結果がNOで、2つの選択還元触媒14,15の少なくとも一方が活性化していないときには、ステップ8に進み、目標当量比KCMDを所定の理論空燃比用値KSTOICに設定する。この理論空燃比用値KSTOICは、理論空燃比に相当する当量比の値(=1)に設定されている。
一方、ステップ4の判別結果がYESで、2つの選択還元触媒14,15の少なくとも一方が活性化しているときには、ステップ5に進み、エンジン回転数NEおよびアクセル開度APに応じて、図示しないマップを検索することにより、要求トルクTRQDRVを算出する。
次いで、ステップ6に進み、エンジン回転数NEおよび要求トルクTRQDRVに応じて、前述した図4のマップを検索することにより、目標当量比KCMDを算出する。
以上のステップ6〜8のいずれかに続くステップ9で、LAFセンサ活性化フラグF_LAFACTが「1」であるか否かを判別する。このLAFセンサ活性化フラグF_LAFACTは、図示しない判定処理において、LAFセンサ23が活性化しているときに「1」に、それ以外のときに「0」に設定される。
ステップ9の判別結果がNOで、LAFセンサ23が活性化していないときには、ステップ13に進み、燃料噴射量Toutを目標当量比と基本噴射量の積KCMD・TiBSに設定する。その後、本処理を終了する。
一方、ステップ9の判別結果がYESで、LAFセンサ23が活性化しているときには、ステップ10に進み、エンジン回転数NEおよび吸気圧PBに応じて、図示しないマップを検索することにより、排ガスボリュームVexを算出する。
次いで、ステップ11に進み、前述した制御アルゴリズムにより、空燃比補正係数KAFを算出する。具体的には、まず、前述した式(3)〜(7)の予測アルゴリズムと、図6から検索された重み関数値Wdiとを用いて、予測当量比PRE_KACTを算出する。さらに、前述した式(17)〜(29)の同定アルゴリズムと、図10,11からそれぞれ検索された基準モデルパラメータαbsおよび重み関数値Waiとを用いて、同定値αidを算出する。そして、算出した予測当量比PRE_KACTおよび同定値αidを用い、前述した式(34),(35)により、空燃比補正係数KAFが最終的に算出される。
ステップ11に続くステップ12で、燃料噴射量Toutを空燃比補正係数と基本噴射量の積KAF・TiBSに設定する。その後、本処理を終了する。
本実施形態の制御装置1は、以上のような空燃比制御処理によって燃料噴射量Toutを算出し、図示しないが、この燃料噴射量Toutおよびエンジン回転数NEに応じて、燃料噴射タイミングを算出するとともに、これらの燃料噴射量Toutおよび燃料噴射タイミングに基づいた制御入力信号で、燃料噴射弁4を駆動することにより、混合気の空燃比を制御する。
次に、図21〜27を参照しながら、本実施形態の制御装置1による空燃比制御のシミュレーション結果(以下「制御結果」という)について説明する。まず、図21〜23について説明する。これらの図21〜図23はいずれも、シミュレーション条件を、式(2)の制御対象モデルにおいてモデル化誤差が存在しないと設定した場合(具体的には、α=αbsと設定した場合)のものであり、図21は本実施形態の制御装置1による制御結果を示している。
また、図22は、比較のために、シミュレーション条件として、制御装置1における状態予測器40およびオンボード同定器60の演算を省略した場合、具体的にはPRE_KACT(k)=KACT(k),αid(k)=αbs(k)と設定した場合の制御結果(以下「比較例1」という)を示している。さらに、図23は、比較のために、制御装置1における状態予測器40およびオンボード同定器60の演算を省略するとともに、感度設定パラメータβとして、実施形態の設定値の1/6の値を用いた場合の制御結果(以下「比較例2」という)を示している。
まず、図22の比較例1を参照すると、排ガスボリュームVexが小さいときには、予測当量比PRE_KACTすなわち検出当量比KACTが発散し、それに従って、予測追従誤差PRE_eや空燃比補正係数KAFも発散した状態となっていることが判る。すなわち、むだ時間が存在する制御対象を、周波数整形コントローラ130のみを用いて制御した場合、相補感度関数Tdによって指定したロバスト性を適切に維持することができず、特に、むだ時間dが増大するような、排ガスボリュームVexが小さい条件下では、制御入力としての空燃比補正係数KAFが発散してしまうことが判る。
次に、図23の比較例2を参照すると、この比較例2の場合、上述した比較例1と比べて、制御系の安定性および制御精度が向上していることが判る。これは、感度関数Sdにおける感度設定パラメータβを比較例1の1/6の値に設定し、フィードバックゲインを低下させたことによるものであり、言い換えれば、外乱抑制能力を低下させたことに起因するものである。この場合の感度設定パラメータβの値は、トライ&エラー手法によって、制御系の安定性を維持できる限界値に設定されている。したがって、前述した制御条件φを満たすように、感度関数Sdを設定することで、外乱抑制能力を直接的に指定するという本発明の目的を実現できないものとなっている。
これに対して、図21の本実施形態の制御結果の場合、状態予測器40、オンボード同定器60および周波数整形コントローラ130の各種アルゴリズムによって、モデル化誤差が存在しない条件下において、制御の安定性およ制御精度が比較例1,2よりも向上していることが判る。
次に、図24〜27について説明する。これらの図24〜27はいずれも、シミュレーション条件を、式(2)の制御対象モデルにおいてモデル化誤差が存在するように設定した場合(具体的には、α=2・αbsとなるように設定した場合)のものであり、図24は本実施形態の制御装置1による制御結果を示している。
また、図25は、比較のために、シミュレーション条件として、制御装置1における状態予測器40およびオンボード同定器60の演算を省略した場合の制御結果(以下「比較例3」という)を示している。さらに、図26は、比較のために、制御装置1における状態予測器40およびオンボード同定器60の演算を省略するとともに、感度設定パラメータβとして、実施形態の設定値の1/6の値を用いた場合の制御結果(以下「比較例4」という)を示している。これに加えて、図27は、比較のために、制御装置1におけるオンボード同定器60の演算のみを省略した場合、すなわちαid(k)=αbs(k)と設定した場合の制御結果(以下「比較例5」という)を示している。
まず、図25の比較例3を参照すると、この比較例3の場合、制御対象モデルにおいてモデル化誤差が存在するシミュレーション条件下では、むだ時間に起因する制御系の安定余裕の低下に加えて、モデル化誤差の影響によって、制御系の安定性が損なわれ、排ガスボリュームVexの全領域に対して、空燃比補正係数KAFなどの各パラメータがすべて発散した状態となっていることが判る。すなわち、むだ時間が存在する制御対象を周波数整形コントローラ130のみを用いて制御した場合、モデル化誤差が存在するシミュレーション条件下では、制御の安定性および制御精度が極めて低下してしまうことが判る。
次に、図26の比較例4を参照すると、この比較例4の場合、上述した比較例3のような発散状態は発生しておらず、比較例3と比べて、制御系の安定性が向上していることが判る。これは、前述したように、感度設定パラメータβの設定値に起因するものである。しかし、この比較例4の場合、比較例3と比べて、制御系の安定性が向上しているものの、予測追従誤差PRE_eが一時的に過大な値となる状態が発生しており、制御精度が低い状態にあることが判る。これに加えて、前述したように、感度設定パラメータβを実施形態の設定値の1/6の値に設定している以上、制御条件φを満たすように、感度関数Sdを設定することで、外乱抑制能力を直接的に指定するという本発明の目的を実現することができないものとなっている。
また、図27の比較例5を参照すると、この比較例5の場合、上述した比較例3と比べて、制御系の安定性が向上していることが判る。これは、排ガスボリュームVexの変化に伴ってむだ時間dが連続的に変化したときでも、状態予測器40において、そのようなむだ時間dの変化を反映させながら予測当量比PRE_KACTが算出されるので、むだ時間dの変化の影響を適切に補償でき、制御系の安定性が向上することによる。しかし、この比較例5の場合でも、排ガスボリュームVexが大きい領域では、モデル化誤差が増大するのに起因して、空燃比補正係数KAFなどの各パラメータが発散した状態となっていることが判る。
一方、図24の本実施形態の制御結果の場合、モデル化誤差が存在するシミュレーション条件下でも、状態予測器40、オンボード同定器60および周波数整形コントローラ130の各種アルゴリズムによって、制御の安定性および制御精度が比較例3〜5よりも向上していることが判る。例えば、状態予測器40の予測アルゴリズムによって、予測追従誤差PRE_eが小さい値に保持されているとともに、オンボード同定器60の同定アルゴリズムによって、時間の経過に伴い、同定値αidがモデルパラメータαに収束していることが判る。
以上のように、第1実施形態の制御装置1によれば、状態予測器40において、検出当量比KACTと空燃比補正係数KAFの関係を定義した制御対象モデル(式(2))を用い、むだ時間d=0〜3がそれぞれ経過したタイミングでの検出当量比KACTとして、第0〜第3予測値PRE_KACT_0〜3が算出され、排ガスボリュームVexに応じて、4つの重み関数値Wd1〜Wd4が算出される。そして、予測当量比PRE_KACTが、重み関数値Wdiと予測値PRE_KACT_4−i(i=1〜4)との乗算値の総和として算出されるので、予測当量比PRE_KACTを、予測値PRE_KACT_4−iを連続的に結合させた値として算出することができる。それにより、排ガスボリュームVexの変化に応じてむだ時間dが変化した場合でも、予測検出当量比KACTを、そのようなむだ時間dの変化を適切に補償しながら精度よく算出することができる。特に、排ガスボリュームVexの急変に伴ってむだ時間dが急変したときでも、予測検出当量比KACTを、むだ時間dの急変を適切に補償しながら、段差なく円滑に変化するように算出することができる。
また、オンボード同定器60において、前述した同定アルゴリズムにより、同定値αidが算出されるので、前述した同定条件1,2の双方を満たしながら、同定値αidを算出することができる。具体的には、合成信号値W_actおよび推定合成信号値W_hatが互いに一致するように、同定値αidが算出されるので、同定条件1すなわち拘束条件を満たしながら、同定値αidを算出することができる。さらに、修正後制御入力KAF_modが、4つの重み関数値Wd4〜1を、むだ時間d=0〜3前のタイミングでの空燃比補正係数KAF(k),KAF(k−1),KAF(k−2),KAF(k−3)にそれぞれ乗算した乗算値の総和として算出されるので、排ガスボリュームVexの変化に応じてむだ時間dが連続的に変化する場合でも、修正後制御入力KAF_modを、そのようなむだ時間dの変化を適切に補償しながら精度よく算出することができる。特に、排ガスボリュームVexの急変に伴ってむだ時間dが急変したときでも、修正後制御入力KAF_modを、むだ時間dの急変を適切に補償しながら、段差なく円滑に変化するように算出することができる。
さらに、そのような修正後制御入力KAF_modと検出当量比KACTとの関係を定義した式(30)のモデルを用いて導出した式(17)〜(29)の同定アルゴリズムにより、同定値αidがオンボード同定される。すなわち、2種類の重み関数値Wdi,Waiを用いながら、同定値αidが算出されるので、排ガスボリュームVexの変化に応じてむだ時間dや遅れ特性が変化したときでも、その影響を抑制しながら、同定値αidを精度よく同定することができる。特に、排ガスボリュームVexの急変に伴ってむだ時間dや遅れ特性が急変したときでも、同定値αidを、むだ時間dや遅れ特性の急変を適切に補償しながら、段差なく円滑に変化するように算出することができる。そして、そのように算出された同定値αidを用いて、制御入力としての空燃比補正係数KAFが算出されるので、空燃比制御における制御性や、個体間のばらつきおよび経年変化の影響に対する空燃比制御のロバスト性を飛躍的に向上させることができる。
これに加えて、周波数整形コントローラ130の場合、前述したように、所定の周波数特性が得られるように設定された感度関数Sdに基づいて導出された算出式(34),(35)を用いて、空燃比補正係数KAFが算出されるので、前述した制御条件φを満たしながら、空燃比補正係数KAFを算出することができる。これに加えて、制御対象モデルのモデルパラメータとして前述した同定値αidを用いているので、むだ時間dや遅れ特性の変化を適切に補償しながら、制御装置1における外乱抑制特性やロバスト性を周波数軸上でダイレクトに指定(設定)することができる。それにより、外乱による検出当量比KACTの変動を抑制したい周波数域において、外乱抑制能力やロバスト性を飛躍的に向上させることができる。これに加えて、空燃比補正係数KAFの算出アルゴリズムにおいて、予測当量比PRE_KACTと目標当量比KCMDとの偏差に基づいたフィードバック制御アルゴリズムが用いられているので、むだ時間dを補償できることによって、フィードバックゲインを高い状態に維持することができ、高精度かつ高応答性を確保しながら、検出当量比KACTを目標当量比KCMDに追従させることができる。
なお、第1実施形態は、重み関数値として、互いに重なり合う領域における排ガスボリュームVexのいずれかの値に対応する重み関数値Wdiの値の和が、各重み関数値Wdiにおける最大値1と等しくなるように設定されているものを用いた例であるが、本発明の重み関数値はこれに限らず、互いに重なり合う領域における参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が所定値になるように設定されているものであればよい。例えば、互いに重なり合う領域における参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が、重み関数値における絶対値の最大値と等しくなるように設定されたものを用いてもよい。より具体的には、重み関数値として、図6の重み関数値Wdiの設定値に対して、x軸を対称軸として線対称に設定されたもの、すなわち負値に設定されたものを用いてもよく、その場合には、4つの重み関数値に乗算する値、すなわち4つの予測値PRE_KACT__4−iや、4つの空燃比補正係数KAF(k−4+i)として、負値化した値を用いればよい。
また、第1実施形態のオンボード同定器60において、前述した式(17)〜(29)の同定アルゴリズムに代えて、以下の式(54)〜(67)の同定アルゴリズムを用いて、同定値αidを算出するように構成してもよい。
以上の式(54)〜(67)を前述した式(17)〜(29)と比較すると明らかなように、式(54)〜(64)は式(17)〜(27)と同一であり、式(65)〜(67)のみが異なっているので、以下、これらの式(65)〜(67)についてのみ説明する。まず、上式(65)を用いて、基準モデルパラメータαbs’が算出される。すなわち、基準モデルパラメータαbs’は、前述した基準モデルパラメータαbsを補正係数Kαbsで補正することによって算出される。
この補正係数Kαbsは、エンジン回転数NEおよび検出当量比KACTに応じて、図28に示すマップを検索することにより算出される。同図において、3つの値KACT_R,KACT_S,KACT_Lはいずれも検出当量比KACTの所定値であり、KACT_S=1,KACT_L<KACT_S<KACT_Rが成立するように設定されている。
このマップでは、補正係数Kαbsは、値1以下の値に設定されているとともに、エンジン回転数NEが低いほど、より小さい値に設定されている。これは、低回転域では、排ガスボリュームVexが同じでも、1燃焼サイクルの実行時間が長くなるのに伴い、排ガス成分の周期的変動が大きくなることによって、空燃比補正係数KAFと検出当量比KACTとの間の応答遅れが大きくなるので、それに対応するためである。
また、補正係数Kαbsは、検出当量比KACTがリッチ側の値であるほど、より大きい値に設定されている。これは、検出当量比KACTが大きく、排ガス濃度が高い場合には、排ガス中の未燃成分が多くなり、LAFセンサ23の検出素子の応答性が高くなることによって、空燃比補正係数KAFと検出当量比KACTとの間の応答遅れが小さくなるので、それに対応するためである。
次に、前述した式(66)を用いて、修正項dαijh’を算出した後、前述した式(67)により、同定値αidが最終的に算出される。同式(66)において、Wanj,Waahは重み関数値であり、重み関数値Wanj(j=1〜4)は、エンジン回転数NEに応じて、図29に示すマップを検索することにより算出される。同図において、NE1〜4,NEMAXは、エンジン回転数NEの所定値であり、0<NE1<NE2<NE3<NE4<NEMAXが成立するように設定されている。また、所定値NEMAXは、最大許容回転数に設定されている。
同図29に示すように、4つの重み関数値Wan1〜4はそれぞれ、エンジン回転数NEが変化し得る領域を0≦NE≦NE2,NE1≦NE≦NE3,NE2≦NE≦NE4,NE3≦NE≦NEMAXの4つの領域に区分した場合において、これらの4つの領域に対応するように設定されているとともに、対応する領域では、値1以下の正の値に設定され、対応する領域以外では、値0に設定されている。
具体的には、重み関数値Wan1は、これが対応する領域(0≦NE≦NE2)では、NE≦NE1のときの値1を最大値として、エンジン回転数NEが大きいほど、より小さい正の値に設定されているとともに、それ以外の領域では、値0に設定されている。また、重み関数値Wan2は、これが対応する領域(NE1≦NE≦NE3)では、NE=NE2のときの値1を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値0に設定されている。
さらに、重み関数値Wan3は、これが対応する領域(NE2≦NE≦NE4)では、NE=NE3のときの値1を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値0に設定されている。一方、重み関数値Wan4は、これが対応する領域(NE3≦NE≦NEMAX)では、NE4≦NEのときの値1を最大値として、エンジン回転数NEが大きいほど、より大きい正の値に設定されているとともに、それ以外の領域では、値0に設定されている。
以上に加えて、4つの重み関数値Wanj(j=1〜4)がそれぞれ対応する4つの領域は、上述したような、隣り合う領域が互いに重なり合うように設定されており、これらの互いに重なり合う領域におけるエンジン回転数NEのいずれかの値に対応する重み関数値Wanjの値の和は、各重み関数値Waniにおける最大値1と等しくなるように設定されている。以上のようにエンジン回転数NEに応じて算出される重み関数値Wanjを用いた理由は、補正係数Kαbsの算出で述べた理由とおなじである。
また、前述した式(66)の重み関数値Waah(h=1〜4)は、検出当量比KACTに応じて、図30に示すマップを検索することにより、4つの重み関数値Waa1〜4がそれぞれ算出される。同図において、KACT1〜4,KACTMAXは、検出当量比KACTの所定値であり、0<KACT1<KACT2<KACT3<KACT4<KACTMAXが成立するように設定されている。さらに、所定値KACTMAXは、エンジン3の運転中に検出当量比KACTが変化し得る領域の最大値に設定されている。言い換えれば、検出当量比KACTは、エンジン3の運転中、0〜KACTMAXの領域内で変化する特性を有している。
同図30に示すように、4つの重み関数値Waa1〜4はそれぞれ、検出当量比KACTが変化し得る領域をKACT≦KACT2,KACT1≦KACT≦KACT3,KACT2≦KACT≦KACT4,KACT3≦KACT≦KACTMAXの4つの領域に区分した場合において、これらの4つの領域に対応するように設定されているとともに、対応する領域では、値1以下の正の値に設定され、対応する領域以外では、値0に設定されている。
具体的には、重み関数値Waa1は、これが対応する領域(KACT≦KACT2)では、KACT≦KACT1のときの値1を最大値として、検出当量比KACTが大きいほど、より小さい正の値に設定されているとともに、それ以外の領域では、値0に設定されている。また、重み関数値Waa2は、これが対応する領域(KACT1≦KACT≦KACT3)では、KACT=KACT2のときの値1を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値0に設定されている。
さらに、重み関数値Waa3は、これが対応する領域(KACT2≦KACT≦KACT4)では、KACT=KACT3のときの値1を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値0に設定されている。一方、重み関数値Waa4は、これが対応する領域(KACT3≦KACT≦KACTMAX)では、KACT4≦KACTのときの値1を最大値として、検出当量比KACTが大きいほど、より大きい正の値に設定されているとともに、それ以外の領域では、値0に設定されている。
以上に加えて、4つの重み関数値Waah(h=1〜4)がそれぞれ対応する4つの領域は、上述したような、隣り合う領域が互いに重なり合うように設定されており、これらの互いに重なり合う領域における検出当量比KACTのいずれかの値に対応する重み関数値Waahの値の和は、各重み関数値Waahにおける最大値1と等しくなるように設定されている。以上のように、検出当量比KACTに応じて算出される重み関数値Waahを用いた理由は、補正係数Kαbsの算出で述べた理由とおなじである。
以上の同定アルゴリズムを用いて同定値αidを算出した場合、排ガスボリュームVexに加えて、エンジン回転数NEおよび検出当量比KACTの変化に伴うむだ時間dや遅れ特性の変化を反映させながら、同定値αidを算出することができる。すなわち、3つのパラメータVex,NE,KACTの変化に起因するむだ時間dや遅れ特性の変化を補償しながら、同定値αidを算出することができ、同定値αidの同定精度(すなわち算出精度)をさらに向上させることができる。それにより、空燃比制御の制御性およびロバスト性を、第1実施形態のオンボード同定器60を用いた場合よりも向上させることができる。
次に、図31を参照しながら、本発明の第2実施形態に係る制御装置1Aについて説明する。この制御装置1Aは、前述した制御装置1と同様に、空燃比補正係数KAFなどを算出することによって、空燃比を制御するものである。第2実施形態の場合、制御対象モデルとして、第1実施形態で用いた式(2)に代えて、下式(68)を用いる。
上式(68)のδはモデルパラメータである。この式(68)は、式(2)の「1−α」を「δ」に置き換えたものであり、それにより、前述した2つのモデルパラメータ1−α,α間の拘束条件を取り除いたものに相当する。
図31に示すように、制御装置1Aは、目標当量比算出部230、可変むだ時間状態予測器(以下「状態予測器」という)240、オンボード・スケジュールド・モデルパラメータ同定器(以下「オンボード同定器」という)260、周波数整形コントローラ330を備えており、これらはいずれもECU2によって構成されている。
この目標当量比算出部230では、前述した目標当量比算出部30と同じ手法によって、目標当量比KCMDが算出される。また、状態予測器240では、後述する予測アルゴリズムを用いて、予測当量比PRE_KACTが算出され、オンボード同定器260では、後述する同定アルゴリズムを用いて、2つのモデルパラメータδ,αを要素とするモデルパラメータベクトルθが算出される。さらに、周波数整形コントローラ330では、後述する制御アルゴリズムを用いて、制御入力としての空燃比補正係数KAFが算出される。
なお、本実施形態では、目標当量比算出部230が目標制御量算出手段に相当し、目標当量比KCMDが目標制御量に相当する。また、状態予測器240が、予測値算出手段、重み関数値算出手段および予測制御量設定手段に相当し、予測当量比PRE_KACTが予測制御量に相当する。さらに、オンボード同定器260が、修正後制御入力設定手段、同定手段および重み関数値算出手段に相当し、周波数整形コントローラが制御入力算出手段に相当する。
次に、図32を参照しながら、前述した状態予測器240について説明する。同図に示すように、この状態予測器240は、前述した図5の状態予測器40と比較すると、前述した第1〜第3予測値算出部45〜47に代えて、第1〜第3予測値算出部245〜247を備えている点のみが異なっており、それ以外は状態予測器40と同一の構成を備えているので、以下、異なる点を中心に説明するとともに、状態予測器40と同一の構成については、同じ符号を付し、その説明は適宜、省略する。
まず、増幅器44では、前述した式(3)と同じ下式(69)により、予測当量比PRE_KACT0が算出される。
また、第1予測値算出部245では、遅延要素41で1制御サイクル分遅延された空燃比補正係数の値KAF(k−1)を用い、下式(70)により、第1予測値PRE_KACT_1が算出される。この式(70)のモデルパラメータδ,αは、オンボード同定器260で同定される。
さらに、第2予測値算出部246では、2つの遅延要素41,42でそれぞれ1,2制御サイクル分ずつ、遅延された空燃比補正係数の値KAF(k−1),KAF(k−2)を用い、下式(71)により、第2予測値PRE_KACT_2が算出される。
一方、第3予測値算出部247では、3つの遅延要素41〜43でそれぞれ1〜3制御サイクル分ずつ、遅延された空燃比補正係数の値KAF(k−1),KAF(k−2),KAF(k−3)を用い、下式(72)により、第3予測値PRE_KACT_3が算出される。
なお、以上の式(70)〜(72)は、前述した制御対象モデルの式(68)に基づき、前述した式(4)〜(6)と同じ手法により導出される。
さらに、前述したように、4つの重み関数値算出部48〜51では、4つの重み関数値Wd1〜Wd4がそれぞれ算出され、4つの乗算器52〜55では、4つの乗算値Wd4・PRE_KACT_0,Wd3・PRE_KACT_1,Wd2・PRE_KACT_2,Wd1・PRE_KACT_3がそれぞれ算出される。
そして、加算器56で、前述した式(7)と同じ下式(73)により、予測当量比PRE_KACTが算出される。
以上の手法により予測当量比PRE_KACTを算出した場合でも、状態予測器40と同様の作用効果を得ることができる。すなわち、排ガスボリュームVexの変化に応じて、むだ時間dが値0から値3の間で連続的に変化したときでも、そのようなむだ時間の変化を適切に反映させながら、予測当量比PRE_KACTを算出することができる。
次に、前述したオンボード同定器260について説明する。このオンボード同定器260では、以下に述べるスケジュールド補正型の同定アルゴリズムを用いて、モデルパラメータベクトルθが算出される。この同定アルゴリズムは、前述した式(68)の右辺のKAF(k−d)を修正後制御入力KAF_mod(k)に置き換えた修正後モデルに基づいて導出される。
このオンボード同定器260は、図33に示すように、修正後制御入力算出部270、3つの遅延要素261〜263、推定検出当量比算出部265、同定ゲインベクトル算出部266、減算器267、乗算器268およびモデルパラメータベクトル算出部290を備えている。
まず、修正後制御入力算出部270では、前述した修正後制御入力算出部70と同一の手法により、修正後制御入力KAF_modが算出される。
また、推定検出当量比算出部265では、遅延要素261〜263でそれぞれ1制御サイクル分遅延された3つの値KACT(k−1),KAF_mod(k−1),θ(k−1)を用い、下式(74)〜(76)により、推定検出当量比KACT_hatが算出される。
この式(76)は、前述した式(68)の各パラメータを1制御サイクル分過去側にずらした式において、左辺のKACTをKACT_hatに、右辺のKAFをKAF_modにそれぞれ置き換えることによって、導出される。
さらに、減算器267では、下式(77)により、同定誤差eidが算出される。
一方、同定ゲインベクトル算出部266では、下式(78),(79)により、同定ゲインベクトルKpが算出される。この同定ゲインベクトルKpは、モデルパラメータベクトルθにおける要素δ,αの修正方向(正負)および修正量を規定するものである。
上式(78)のIは2次の単位行列であり、Pは、その初期値が下式(80)のように定義される2次の正方行列である。
さらに、上式(78)の場合、前述したように、重みパラメータλ1、λ2を下記のように設定することにより、同定アルゴリズムとして、3つのアルゴリズムのいずれかを選択することができる。
λ1=1,λ2=0 ;固定ゲインアルゴリズム
λ1=1,λ2=1 ;最小2乗法アルゴリズム
λ1=λ,λ2=1 ;重み付き最小2乗法アルゴリズム
ここで、λは0<λ<1に設定される所定値である。本実施形態の場合、同定精度および制御精度を適切に確保するために、重み付き最小2乗法アルゴリズムを用いている。
さらに、乗算器268で、同定誤差eidと同定ゲインベクトルKpの乗算値eid・Kpが算出される。
次いで、モデルパラメータベクトル算出部290では、上記の乗算値eid・Kpと排ガスボリュームVexを用いて、以下に述べるようにモデルパラメータベクトルθが算出される。このモデルパラメータベクトル算出部290は、図34に示すように、基準モデルパラメータ算出部291、基準モデルパラメータベクトル算出部317、4つの重み関数値算出部292〜295、8つの乗算器296〜303、5つの加算器304〜308、4つの遅延要素309〜312および4つの増幅器313〜316を備えている。
まず、基準モデルパラメータ算出部291では、前述した図9の基準モデルパラメータ算出部91と同一の手法により、基準モデルパラメータαbsが算出される。次いで、基準モデルパラメータベクトル算出部317では、下式(81)により、基準モデルパラメータδbsを算出した後、下式(82)により、基準モデルパラメータベクトルθbsが算出される。
また、4つの重み関数値算出部292〜295では、前述した図9の重み関数値算出部92〜95と同一の手法により、4つの重み関数値Wa1〜Wa4がそれぞれ算出される。また、乗算器296で、重み関数値Wa1に値Kp・eidを乗算することにより、乗算値Wa1・Kp・eidが算出され、遅延要素309で遅延された修正項ベクトルdθ1(k−1)に、増幅器313で忘却行列ηを乗算することにより、値η・dθ1(k−1)が算出される。この忘却行列ηについては後述する。そして、加算器304で、値Wa1・Kp・eidに値η・dθ1(k−1)を加算することにより、修正項ベクトルdθ1が算出される。この修正項ベクトルdθ1は、後述する式(84)に示すように、2つの修正項dδ1,dα1を要素とするものである。
さらに、乗算器297で、重み関数値Wa2に値Kp・eidを乗算することにより、乗算値Wa2・Kp・eidが算出され、遅延要素310で遅延された修正項ベクトルdθ2(k−1)に、増幅器314で忘却行列ηを乗算することにより、値η・dθ2(k−1)が算出される。そして、加算器305で、値Wa2・Kp・eidに値η・dθ2(k−1)を加算することにより、修正項ベクトルdθ2が算出される。この修正項ベクトルdθ2は、後述する式(84)に示すように、2つの修正項dδ2,dα2を要素とするものである。
また、乗算器298で、重み関数値Wa3に値Kp・eidを乗算することにより、乗算値Wa3・Kp・eidが算出され、遅延要素311で遅延された修正項ベクトルdθ3(k−1)に、増幅器315で忘却行列ηを乗算することにより、値η・dθ3(k−1)が算出される。そして、加算器306で、値Wa3・Kp・eidに値η・dθ3(k−1)を加算することにより、修正項ベクトルdθ3が算出される。この修正項ベクトルdθ3は、後述する式(84)に示すように、3つの修正項dδ3,dα3を要素とするものである。
さらに、乗算器299で、重み関数値Wa4に値Kp・eidを乗算することにより、乗算値Wa4・Kp・eidが算出され、遅延要素312で遅延された修正項ベクトルdθ4(k−1)に、増幅器316で忘却行列ηを乗算することにより、値η・dθ4(k−1)が算出される。そして、加算器307で、値Wa4・Kp・eidに値η・dθ4(k−1)を加算することにより、修正項ベクトルdθ4が算出される。この修正項ベクトルdθ4は、後述する式(84)に示すように、4つの修正項dδ4,dα4を要素とするものである。
上述した増幅器313〜316で用いられる忘却行列ηは、下式(83)のように定義される。
上式(83)のη1,η2は、忘却係数であり、0<η1≦1,0<η2≦1が成立するように設定される。修正項ベクトルdθi(i=1〜4)の算出において、この忘却行列ηを用いた理由は、エンジン3の定常運転状態が長時間継続した場合、修正項ベクトルdθiが増大化し、不適切な値となるおそれがあるので、それを回避するためである。また、この忘却行列ηにおいて、2つの忘却係数η1,η2のうちの一方を値1に設定した場合、同定誤差eidの定常的な発生を抑制することと、制御系の安定性の確保とを両立させることができる。
また、以上の4つの加算器304〜307における演算式は、下式(84),(85)で表される。
さらに、乗算器300〜303ではそれぞれ、4つの重み関数値Waiを4つの修正項ベクトルdθiに乗算することにより、4つのベクトルWai・dθiが算出される。
そして、加算器308で、下式(86)により、モデルパラメータベクトルθが最終的に算出される。
オンボード同定器260において、以上の同定アルゴリズムを用いたのは、前述した同定条件1,2を満たすためである。すなわち、前述したように、最小2乗法などの一般的な同定アルゴリズムをそのまま用いた場合、同定条件1を満たすことは困難である。そのため、このオンボード同定器260の場合、同定条件1を満たしながら、モデルパラメータを同定するために、制御対象モデルの式(68)のモデルパラメータδ,αの同定演算において、2つのモデルパラメータの基準値δbs,αbを、両者の間に拘束条件(δbs=1−αbs)を設定しながら算出するとともに、修正項ベクトルdθiの方は一般的な逐次型最小2乗法アルゴリズムを用いて算出する手法を採用している。さらに、同定条件2を満たすために、前述したオンボード同定器60と同じように、重み関数値Waiを用いて、修正項ベクトルdθiおよびモデルパラメータベクトルθを算出する手法を採用している。
次に、前述した周波数整形コントローラ330について説明する。この周波数整形コントローラ330は、予測当量比PRE_KACTが目標当量比KCMDに収束するように、空燃比補正係数KAFを算出するもの、言い換えれば、検出当量比KACTが目標当量比KCMDに収束するように、空燃比補正係数KAFを算出するものである。この周波数整形コントローラ330では、まず、前述した式(34)と同一の下式(87)により、予測追従誤差PRE_eが算出される。
次いで、下式(88)により、制御入力としての空燃比補正係数KAFが算出される。
以上の周波数整形コントローラ330の制御アルゴリズムは、前述した周波数整形コントローラ130の制御アルゴリズムの導出手法と同じ手法により導出される。
以上のように構成された第2実施形態の制御装置1Aによれば、第1実施形態の制御装置1と同一の状態予測器40を用いたことによって、予測検出当量比KACTを、むだ時間dの変化を適切に補償しながら精度よく算出することができる。また、周波数整形コントローラ330により、前述した周波数整形コントローラ130と同様に、制御条件φを満たしかつむだ時間dの変化を適切に補償しながら、空燃比補正係数KAFを算出することができる。すなわち、むだ時間dの変化を適切に補償しながら、制御装置1における外乱抑制特性やロバスト性を周波数軸上でダイレクトに指定(設定)することができ、それにより、外乱による検出当量比KACTの変動を抑制したい周波数域において、外乱抑制能力やロバスト性を飛躍的に向上させることができる。
また、オンボード同定器260において、前述したように、修正後制御入力KAF_modおよび重み関数値Waiを用いる同定アルゴリズムによって、2つのモデルパラメータδ,αが同定されるので、前述した同定条件2を満たしながら、2つのモデルパラメータδ,αを算出することができる。これに加えて、モデルパラメータδ,αの場合、基準モデルパラメータδbs,αbs同士が、同定条件1(すなわち拘束条件)を満たすように設定されているとともに、これらを要素とする基準モデルパラメータベクトルθbsを、重み関数値Waiと修正項ベクトルdθiの乗算値の総和で修正することによって、モデルパラメータδ,αを要素とするモデルパラメータベクトルθが算出されるので、2つのモデルパラメータδ,αを、同定条件1を満たす値の近傍の値として同定することができる。
なお、以上のオンボード同定器260の同定アルゴリズムを、前述したオンボード同定器60の同定アルゴリズムと比較した場合、オンボード同定器60の同定アルゴリズムでは、同定値αidを同定条件1を完全に満たすように算出できるので、同定条件1を満たすようにモデルパラメータを同定するという観点からは、オンボード同定器60の同定アルゴリズムの方が優れている。
次に、図35を参照しながら、本発明の第3実施形態に係る制御装置1Bについて説明する。同図に示すように、この制御装置1Bは、前述した第1実施形態の図3に示す制御装置1と比較すると、前述した周波数整形コントローラ130に代えて、2自由度応答指定型コントローラ350を備えている点のみが異なっており、それ以外は制御装置1と同一の構成を備えているので、以下、2自由度応答指定型コントローラ350(制御入力算出手段)についてのみ説明する。
この2自由度応答指定型コントローラ350では、以下の2自由度応答指定型制御アルゴリズムを用いて、空燃比補正係数KAFが算出される。具体的には、まず、下式(89)により、目標当量比のフィルタリング値KCMD_fを算出する。
ここで、POLE_fは目標値フィルタ設定パラメータであり、−1<POLE_f<0の関係が成立するように設定される。
次いで、下式(90)により、予測追従誤差PRE_e_fを算出する。
さらに、下式(91)により、切換関数σ_fを算出する。
ここで、POLEは切換関数設定パラメータであり、−1<POLE<0の関係が成立するように設定される。
次いで、下式(92)により、等価制御入力Ueq_fを算出する。
また、下式(93)により、到達則入力Urch_fを算出する。
ここで、Krchは所定のフィードバックゲインである。
さらに、下式(94)により、適応則入力Uadp_fを算出する。
ここで、Kadpは所定のフィードバックゲインである。
そして、最終的に、下式(95)により、空燃比補正係数KAFが算出される。
なお、以上の式(89)〜(95)の2自由度応答指定型アルゴリズムは、前述した式(53)のKACTをPRE_KACTに置き換えたモデルに基づいて導出される。
以上の第3実施形態に係る制御装置1Bによれば、第1実施形態の制御装置1と同一の状態予測器40およびオンボード同定器60を備えているので、これらによって、第1実施形態の制御装置1と同じ作用効果を得ることができる。これに加えて、2自由度応答指定型コントローラ350において、以上の制御アルゴリズムにより、空燃比補正係数KAFが算出されるので、外乱に起因する、目標当量比KCMDと検出当量比KACTとの偏差の値0への時系列的収束挙動と、目標当量比KCMDの変化に対する検出当量比KACTの追従特性とを別個に直接的に指定することができる。
次に、図36を参照しながら、本発明の第4実施形態に係る制御装置1Cについて説明する。同図に示すように、この制御装置1Cは、前述した第1実施形態の図3に示す制御装置1と比較した場合、適応外乱オブザーバ370(外乱推定値算出手段)を備えている点と、前述した周波数整形コントローラ130に代えて2自由度応答指定型コントローラ380(制御入力算出手段)を備えている点のみが異なっているので、以下、これらの異なる点についてのみ説明する。
この適応外乱オブザーバ370では、以下に述べる制御アルゴリズムにより、外乱推定値εが算出される。まず、下式(96)により、外乱推定用の推定検出当量比KACT_advを算出する。
この式(96)は、前述した式(2)において、KAF(k+1)をKACT_adv(k)に、αをαid(k)に、KAF(k−d)をKAF_mod(k)にそれぞれ置き換えるとともに、右辺に外乱推定値εを加えたもの、すなわち外乱推定用のモデルに相当する。
次いで、下式(97)により、外乱推定用の追従誤差e_advを算出する。
そして、最終的に、下式(98)により、外乱推定値εが算出される。
この式(98)のπは外乱推定ゲインであり、π>0が成立するように設定される。
次に、2自由度応答指定型コントローラ380について説明する。この2自由度応答指定型コントローラ380では、下式(99)〜(104)に示す目標値フィルタ型の2自由度応答指定型制御アルゴリズムよって、空燃比補正係数KAFが算出される。
以上の式(99)〜(104)は、前述した式(89)〜(95)において、外乱推定値εを等価制御入力Ueqの算出式に付加するとともに、適応則入力Uadpを省略したものに相当する。
以上の第4実施形態に係る制御装置1Cによれば、第1実施形態の制御装置1と同一の状態予測器40およびオンボード同定器60を備えているので、これらによって、第1実施形態の制御装置1と同じ作用効果を得ることができる。これに加えて、適応外乱オブザーバ370において、以上の制御アルゴリズムにより外乱推定値εが算出され、2自由度応答指定型コントローラ380において、この外乱推定値εを用いて空燃比補正係数KAFが算出されるので、空燃比制御における外乱抑制能力すなわちロバスト性を向上させることができる。
また、適応外乱オブザーバ370を備えているので、例えば、π>P0に設定し、オンボード同定器60の同定速度を落とすことによって、制御の安定性を向上させることが可能になる。さらに、同じ理由によって、制御系の共振を防止したり、同定値αidの演算結果を適用した制御対象モデルのゲイン特性が過小になるのを防止するために、同定値αidや同定アルゴリズムで用いる入出力データにフィルタリングを施すことができ、より高い制御性を確保することができる。
次に、図37を参照しながら、本発明の第5実施形態に係る制御装置1Dについて説明する。なお、以下の説明では、第1実施形態と同じ構成については、同じ符号を付すとともに、その説明は省略する。この制御装置1Dは、後述する制御アルゴリズムによって、車両駆動系の自動変速機400におけるクラッチ410の接続・遮断動作などを制御するものである。
エンジン3は、自動変速機400および差動ギヤ機構460を介して、駆動輪WH,WHに機械的に連結されており、それにより、エンジン3のトルクは、自動変速機400および差動ギヤ機構460によって変速されながら、駆動輪WH,WHに伝達される。
図37に示すように、自動変速機400は、クラッチ410、主軸401、副軸402、前進1,2速ギヤ列420,430、1−2速用同期噛合機構440および駆動ギヤ450などを備えている。なお、図37では、前進1,2速ギヤ列420,430および1−2速用同期噛合機構440以外のギヤ列および同期噛合機構は省略されている。
クラッチ410(伝達トルク調整機構)は、乾式クラッチタイプのものであり、エンジン3のクランクシャフト3aに連結されたクラッチ板411と、このクラッチ板411と対をなす主軸401に連結されたクラッチ板412と、クラッチ板411をエンジン側に付勢するダイヤフラムスプリング(図示せず)と、クラッチ板411をクラッチ板412側に駆動するクラッチアクチュエータ413などを備えている。
このクラッチアクチュエータ413は、油圧駆動式のものであり、クラッチ電磁弁および油圧アクチュエータなどを組み合わせて構成されている。このクラッチ電磁弁は、ECU2に電気的に接続されており、ECU2から制御入力信号が供給されるのに伴い、油圧アクチュエータへの供給油圧を変化させる。それにより、クラッチアクチュエータ413によるクラッチ板411のクラッチ板412側への駆動状態が変更されることで、クラッチ410の接続・遮断状態が変更される。
前進1,2速ギヤ列420,430は、主軸401上に回動自在に設けられた1,2速主軸ギヤ421,431と、副軸402上に固定され、1,2速主軸ギヤ421,431とそれぞれ常に噛み合う1,2速副軸ギヤ422,432で構成されている。
また、1,2速主軸ギヤ421,431の間には、1−2速用同期噛合機構440が設けられている。この1−2速用同期噛合機構440は、油圧駆動式のものであり、シンクロ電磁弁および油圧アクチュエータなどを組み合わせて構成されている。このシンクロ電磁弁は、ECU2に電気的に接続されており、ECU2から制御入力信号が供給されるのに伴い、油圧アクチュエータへの供給油圧を変化させる。それにより、1−2速用同期噛合機構440は、1速主軸ギヤ421または2速主軸ギヤ431を主軸401に同期させながら噛み合いにより連結する。それにより、前進1速段または前進2速段への変速動作が実行される。
一方、駆動ギヤ450は、差動ギヤ機構460の被駆動ギヤ461と常に噛み合っており、それにより、副軸402の回転に伴い、差動ギヤ機構460を介して駆動輪WH,WHが駆動される。
さらに、制御装置1Dは、ECU2を備えており、このECU2には、前述したクランク角センサ20およびアクセル開度センサ21に加えて、油温センサ26、4つの車輪速度センサ27(1つのみ図示)および主軸回転数センサ28が電気的に接続されている。
この油温センサ26は、サーミスタなどで構成され、前述した油圧アクチュエータなどに供給される作動油の温度である油温Toilを検出して、それを表す検出信号をECU2に出力する。ECU2は、この油温センサ26の検出信号に基づいて、油温Toilを算出する。なお、本実施形態では、油温センサ26が参照パラメータ検出手段に相当し、油温Toilが参照パラメータに相当する。
また、4つの車輪速度センサ27の各々は、対応する車輪の回転速度を検出して、それを表す検出信号をECU2に出力する。ECU2は、これらの車輪速度センサ27の検出信号に基づいて、車速VPなどを算出する。
さらに、主軸回転数センサ28は、クランク角センサ20と同様に、マグネットロータおよびMREピックアップで構成されており、主軸401の回転に伴い、主軸401の回転数を表すパルス信号をECU2に出力する。ECU2は、この主軸回転数センサ28の検出信号に基づいて、主軸401の回転数(以下「主軸回転数」という)NMを算出する。なお、本実施形態では、主軸回転数NMが制御量および出力回転数を表す値に相当する。
次に、本実施形態の制御装置1Dにおけるクラッチ制御の原理について説明する。本実施形態のクラッチ410の場合、クラッチアクチュエータ413への制御入力Uactと主軸回転数NMとの関係は、下式(105)に示すように、一次遅れ系の制御対象モデルとしてモデリングすることができる。
この式(105)において、α”はモデルパラメータを、d”はむだ時間をそれぞれ表している。
また、クラッチ410の場合、そのすべり率(またはクランクシャフト3aと主軸401との間の回転差)によって、駆動輪WH側への伝達トルクが決まるとともに、このすべり率は、クラッチアクチュエータ413によるクラッチ板411の駆動状態によって調整されるという特性を有している。
このクラッチアクチュエータ413は、前述したように、油圧駆動式のものであるので、油温Toilの変化に伴って、応答性が変化してしまう特性を有しており、そのため、クラッチ410のすべり率すなわちクラッチ410のトルク伝達特性は、作動油の温度変化の影響を受けやすいという特性を有している。これに加えて、クラッチ410のすべり率は、クラッチ板411,412の表面温度や各部品の経年変化の影響を受けやすいという特性も有している。
以上の理由により、前述した式(105)のむだ時間d”は、油温Toilおよびクラッチ板411,412の表面温度の変化や、各部品の経年変化の影響を受けやすいので、これらに対するロバスト性を確保する必要がある。このむだ時間d”と油温Toilとの関係をモデリングすると、図38に示すモデル(マップ)が得られる。同図において、Toil1〜Toil4,ToilMAXは、油温Toilの所定値であり、0<Toil1<Toil2<Toil3<Toil4<ToilMAXが成立するように設定されている。また、所定値ToilMAXは、エンジン3の運転中に油温Toilが変化し得る領域の最大値に設定されている。言い換えれば、油温Toilは、エンジン3の運転中、0〜ToilMAXの領域内で変化する特性を有しているので、上記のロバスト性を確保するには、油温Toilの変化に伴うむだ時間d”の変化を反映させながら、制御入力Uactを算出する必要がある。
また、一般に、クラッチでは、動作中に「ジャダ」と呼ばれる高周波数の振動的挙動が生じやすく、これが発生した場合、駆動力が振動的に変化し、運転性が低下してしまう。このような問題は、本実施形態のクラッチ410のような乾式クラッチではより顕著に発生しやすく、この問題を解消するには、前述した制御条件φを満たす制御アルゴリズムを用いる必要がある。
以上の理由により、本実施形態では、むだ時間d”を含む前述した式(105)の制御対象モデルと、前述した周波数整形コントローラ130における制御アルゴリズムと同様の制御アルゴリズムとを用いて、制御入力Uactが算出される。
以下、本実施形態の制御装置1Dの構成および制御アルゴリズムについて説明する。なお、以下に述べる制御アルゴリズムは、車両において、ギヤ段が前進1速段に設定されかつ低速走行中のとき、またはギヤ段が前進1速段に設定されかつ発進中であるときに用いられるアルゴリズムであり、以下の説明では、このようなギヤ段の設定条件および走行条件をまとめて「クラッチ制御条件」という。
この制御装置1Dは、図39に示すクラッチコントローラ500と、図44に示すスロットル弁コントローラ600とを有しており、これらのコントローラ500,600はいずれも、具体的にはECU2によって構成されている。
まず、図39を参照しながら、クラッチコントローラ500について説明する。このクラッチコントローラ500は、上述したクラッチ制御条件が成立しているときに、クラッチ410の接続・遮断状態を制御するものであり、同図に示すように、目標主軸回転数算出部510、可変むだ時間状態予測器(以下「状態予測器」という)520、オンボード・スケジュールド・モデルパラメータ同定器(以下「オンボード同定器」という)530および周波数整形コントローラ540を備えている。
この目標主軸回転数算出部510では、以下に述べる手法により、目標主軸回転数NM_cmdが算出される。まず、アクセル開度APおよび車速VPに応じて、図40に示すマップを検索することにより、目標クラッチすべり率Rslip_cmdを算出する。この目標クラッチすべり率Rslip_cmdは、クラッチすべり率(クラッチ410の入力側回転数と出力側回転数との比NE/NM)の目標となる値である。同図において、AP1〜AP4は、AP1<AP2<AP3<AP4が成立するように設定されるアクセル開度APの所定値であり、特に、AP1はアクセルペダルが全閉状態のときの値に、AP4はアクセルペダルが全開状態のときの値に設定されている。また、同図のVP1は、所定車速である。
同図に示すように、目標クラッチすべり率Rslip_cmdは、VP≦VP1でかつAP>AP1の領域では、アクセル開度APが大きいほど、または車速VPが高いほど、より小さい値に設定されている。これは、アクセル開度APが大きいほど、または車速VPが高いほど、クラッチ410におけるトルクの伝達効率を高める必要があることによる。
次いで、以上のように算出した目標クラッチすべり率Rslip_cmdを用い、下式(106)により、目標主軸回転数NM_cmdが算出される。
なお、本実施形態では、目標主軸回転数算出部510が目標制御量算出手段に相当し、目標主軸回転数NM_cmdが目標制御量に相当する。
次に、前述した状態予測器520について説明する。この状態予測器520では、前述した図38のむだ時間d”の特性を考慮し、前述した第1実施形態の状態予測器40と同様の予測アルゴリズムを用いて、予測主軸回転数PRE_NMが算出される。なお、本実施形態では、状態予測器520が、予測値算出手段、重み関数値算出手段および予測制御量設定手段に相当し、予測主軸回転数PRE_NMが予測制御量に相当する。
この予測主軸回転数PRE_NMは、むだ時間d”が経過したタイミングにおける主軸回転数NMを予測した値に相当するものであり、具体的には、以下の式(107)〜(111)に示す予測アルゴリズムにより算出される。また、この予測アルゴリズムは、第1実施形態の状態予測器40の予測アルゴリズムと同じ手法により導出される。
まず、下式(107)により、第0予測値PRE_NM_0を算出する。
また、下式(108)により、第1予測値PRE_NM_1を算出する。
この式(108)のαid”は、モデルパラメータα”の同定値であり、後述するように、オンボード同定器530で算出される。
さらに、下式(109)により、第2予測値PRE_NM_2を算出する。
次いで、下式(110)により、第3予測値PRE_NM_3を算出する。
そして、最終的に、下式(111)により、予測主軸回転数PRE_NMが算出される。
上式(111)のWdi”(i=1〜4)は、重み関数値であり、この重み関数値Wdi”(i=1〜4)は、油温Toilに応じて、図41に示すマップを検索することにより算出される。同図に示すように、4つの重み関数値Wd1”〜Wd4”はそれぞれ、油温Toilが変化し得る領域をToil≦Toil2,Toil1≦Toil≦Toil3,Toil2≦Toil≦Toil4,Toil3≦Toil≦ToilMAXの4つの領域に区分した場合において、これらの4つの領域に対応するように設定されているとともに、対応する領域では、値1以下の正の値に設定され、対応する領域以外では、値0に設定されている。
具体的には、重み関数値Wd1”は、これが対応する領域(Toil≦Toil2)では、Toil≦Toil1のときの値1を最大値として、油温Toilが大きいほど、より小さい正の値に設定されているとともに、それ以外の領域では、値0に設定されている。また、重み関数値Wd2”は、これが対応する領域(Toil1≦Toil≦Toil3)では、Toil=Toil2のときの値1を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値0に設定されている。
さらに、重み関数値Wd3”は、これが対応する領域(Toil2≦Toil≦Toil4)では、Toil=Toil3のときの値1を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値0に設定されている。一方、重み関数値Wd4”は、これが対応する領域(Toil3≦Toil≦ToilMAX)では、Toil4≦Toilのときの値1を最大値として、油温Toilが大きいほど、より大きい正の値に設定されているとともに、それ以外の領域では、値0に設定されている。
以上に加えて、4つの重み関数値Wdi”(i=1〜4)がそれぞれ対応する4つの領域は、上述したような、隣り合う領域が互いに重なり合うように設定されており、これらの互いに重なり合う領域における油温Toilのいずれかの値に対応する重み関数値Wdi”の値の和は、各重み関数値Wdi”における最大値1と等しくなるように設定されている。
さらに、図41と前述した図38を比較すると明らかなように、互いに重なり合う3つの領域は、むだ時間d”の勾配が一定状態に保持される3つの領域に対応するように設定されている。これに加えて、重み関数値Wd1”はむだ時間d”=3に対して、重み関数値Wd2”はむだ時間d”=2に対して、重み関数値Wd3”はむだ時間d”=1に対して、重み関数値Wd4”はむだ時間d”=0に対して重みが最も大きくなるようにそれぞれ設定されている。
したがって、このように設定されている4つの重み関数値Wdi”を、4つの予測値PRE_NM_4−iにそれぞれ乗算した値の総和として、予測主軸回転数PRE_NMが算出されるので、油温Toilの変化に応じて、むだ時間d”が図38に示すように値0から値3の間で連続的に変化したときでも、そのようなむだ時間d”の変化を適切に反映させながら、予測主軸回転数PRE_NMを円滑に変化するような値として算出することができる。
次に、前述したオンボード同定器530について説明する。なお、本実施形態では、オンボード同定器530が、修正後制御入力設定手段、同定手段および重み関数値算出手段に相当する。このオンボード同定器530では、下式(112)〜(124)に示す拘束条件付きスケジュールド補正型の同定アルゴリズムを用いて、同定値αid”が算出される。この同定アルゴリズムは、前述したオンボード同定器60の同定アルゴリズムと同じ手法によって導出される。
まず、下式(112)により、修正後制御入力Uact_modを算出する。
また、下式(113)により、合成信号値W_act”を算出する。
さらに、下式(114),(115)により、推定合成信号値W_hat”を算出する。
次いで、下式(116)により、同定誤差eid”を算出する。
また、下式(117),(118)により、同定ゲインKp”を算出する。
上式(117)のゲインP”の初期値P”(0)は、下式(119)のように定義される。
ここで、P0”は所定値に設定されている。
また、上式(117)のλ1,λ2は重みパラメータであり、前述したように、これらの値λ1、λ2を下記のように設定することにより、同定アルゴリズムとして、3つのアルゴリズムのいずれかを選択することができる。
λ1=1,λ2=0 ;固定ゲインアルゴリズム
λ1=1,λ2=1 ;最小2乗法アルゴリズム
λ1=λ,λ2=1 ;重み付き最小2乗法アルゴリズム
ここで、λは0<λ<1に設定される所定値である。本実施形態の場合、同定精度および制御精度を適切に確保するために、重み付き最小2乗法アルゴリズムを用いている。
次いで、下式(120)〜(122)により、ゲイン係数H”を算出する。
上式(120)〜(122)において、α_L”は所定の下限値であり、α_H”は所定の上限値である。また、η”は忘却係数であり、0<η”≦1が成立するように設定される。同定値αid”の算出において、この忘却係数η”を用いた理由は、定常状態が長時間継続した場合、同定値αid”が増大化し、不適切な値となるおそれがあるので、それを回避するためである。
さらに、下式(123)により、4つの修正項dαi”(i=1〜4)を算出する。
上式(123)のWai”は、重み関数値であり、油温Toilに応じて、図42に示すマップを検索することにより算出される。同図において、Toil5〜Toil8は、油温Toilの所定値であり、Toil5≦Toil6≦Toil7≦Toil8≦ToilMAXが成立するように設定されている。同図に示すように、4つの重み関数値Wa1”〜Wa4”はそれぞれ、油温Toilが変化し得る領域をToil≦Toil6,Toil5≦Toil≦Toil7,Toil6≦Toil≦Toil8,Toil7≦Toil≦ToilMAXの4つの領域に区分した場合において、これらの4つの領域に対応するように設定されているとともに、対応する領域では、値1以下の正の値に設定され、対応する領域以外では、値0に設定されている。
また、重み関数値Wa1”は、これが対応する領域(Toil≦Toil6)では、Toil≦Toil5のときの値1を最大値として、油温Toilが大きいほど、より小さい正の値に設定されているとともに、それ以外の領域では、値0に設定されている。また、重み関数値Wa2”は、これが対応する領域(Toil5≦Toil≦Toil7)では、Toil=Toil6のときの値1を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値0に設定されている。
さらに、重み関数値Wa3”は、これが対応する領域(Toil6≦Toil≦Toil8)では、Toil=Toil7のときの値1を最大値として、三角形の斜辺状に変化する値に設定されているとともに、それ以外の領域では、値0に設定されている。一方、重み関数値Wa4”は、これが対応する領域(Toil7≦Toil≦ToilMAX)では、Toil8≦Toilのときの値1を最大値として、油温Toilが大きいほど、より大きい正の値に設定されているとともに、それ以外の領域では、値0に設定されている。
以上に加えて、4つの重み関数値Wai”(i=1〜4)がそれぞれ対応する4つの領域は、上述したような、隣り合う領域が互いに重なり合うように設定されており、これらの互いに重なり合う領域における油温Toilのいずれかの値に対応する重み関数値Wai”の値の和は、各重み関数値Wai”における最大値1と等しくなるように設定されている。この図42を後述する図43と比較すると明らかなように、これらの互いに重なり合う3つの領域は、基準モデルパラメータαbs”の勾配が一定状態に保持される3つの領域に対応するように設定されている。
そして、下式(124)により、同定値αid”が最終的に算出される。
上式(124)のαbs”は、基準モデルパラメータであり、油温Toilに応じて、図43に示すマップを検索することにより算出される。このマップでは、基準モデルパラメータαbs”は、油温Toilが高いほど、より大きい値に設定されている。これは、油温Toilが高いほど、クラッチアクチュエータの応答性が高くなり、応答遅れが小さくなることで、制御入力Uactが主軸回転数NMに対して及ぼす影響の度合がより大きくなるので、それに対応するためである。
次に、前述した周波数整形コントローラ540(制御入力算出手段)について説明する。この周波数整形コントローラ540では、目標主軸回転数NM_cmd、予測主軸回転数PRE_NMおよび同定値αid”を用いて、下式(125),(126)により、制御入力Uactが算出される。なお、下式(125),(126)の制御アルゴリズムは、前述した周波数整形コントローラ130の制御アルゴリズムと同じ原理によって導出される。
上式(125)のPRE_e”は予測追従誤差である。上式(126)のβ”は感度設定パラメータであり、前述した制御条件φを満たすように設定されている。
周波数整形コントローラ540では、以上のように制御入力Uactが算出される。そして、ECU2によって、この制御入力Uactに対応する制御入力信号がクラッチアクチュエータ413に供給されることにより、主軸回転数NMが目標主軸回転数NM_cmdに収束するようにフィードバック制御される。
次に、図44を参照しながら、前述したスロットル弁コントローラ600について説明する。このスロットル弁コントローラ600は、スロットル弁6aの開度を制御するものであり、同図に示すように、目標エンジントルク算出部610、目標TH開度算出部620およびTHコントローラ630を備えている。
この目標エンジントルク算出部610では、アクセル開度APおよび車速VPに応じて、図45に示すマップを検索することにより、目標エンジントルクTRQ_ENG_cmdが算出される。同図において、TRQ_MAXは、エンジン3が発生可能なトルクの最大値である。また、図中のハッチングで示す領域は、アクセルペダルが全閉状態にあり(AP=AP1)、かつ車両が走行状態にある(VP>VP1)ことで、フューエルカット運転を実行すべき領域であり、そのため、この領域では、目標エンジントルクTRQ_ENG_cmdが負値に設定されている。
また、目標TH開度算出部620では、目標エンジントルクTRQ_ENG_cmdおよびエンジン回転数NEに応じて、図46に示すマップを検索することにより、目標TH開度TH_cmdが算出される。同図において、NE5〜7は、エンジン回転数NEの所定値であり、0<NE5<NE6<NE7<NEMAXが成立するように設定されている。このマップの場合、目標TH開度TH_cmdは、高回転域では、目標エンジントルクTRQ_ENG_cmdが大きいほど、それを実現可能な吸入空気量を確保するために、より大きい値に設定されている。また、目標TH開度TH_cmdは、エンジン回転数NEが高いほど、それを実現可能な吸入空気量を確保するために、より大きい値に設定されている。
次に、THコントローラ630では、目標TH開度TH_cmdに応じて、図示しないマップを検索することにより、制御入力Uthが算出される。そして、ECU2によって、この制御入力Uthに対応する制御入力信号がTHアクチュエータ6bに供給されることにより、スロットル弁6aの開度が目標TH開度TH_cmdになるように制御される。
次に、図47を参照しながら、第5実施形態の制御装置1Dによるクラッチ制御のシミュレーション結果(以下「制御結果」という)について説明する。同図において、Dslipは、すべり率偏差であり、実際のクラッチすべり率Rslip(=NE/NM)と、目標クラッチすべり率Rslip_cmdとの偏差(=Rslip−Rslip_cmd)を表している。
同図に示すように、時刻t1で、アクセルペダルが踏まれ、アクセル開度APがAP1(=0)から上昇すると、その直後、実際のクラッチすべり率Rslipが目標クラッチすべり率Rslip_cmdに対してオーバーシュートすることで、すべり率偏差Dslipが一時的に急増する。しかし、制御の進行に伴って、すべり率偏差Dslipが減少するとともに、時刻t2〜t3の間では、すべり率偏差Dslipが値0近傍に保持されており、高い制御精度が確保されていることが判る。
そして、時刻t3で、アクセルペダルが開放された以降、実際のクラッチすべり率Rslipが目標クラッチすべり率Rslip_cmdに対してアンダーシュートすることで、すべり率偏差Dslipが一時的に急減する。しかし、制御の進行に伴って、すべり率偏差Dslipが値0に向かって上昇するとともに、時刻t4〜t5の間では、すべり率偏差Dslipが値0近傍に保持されており、高い制御精度が確保されていることが判る。
そして、時刻t5で、アクセルペダルが再度、踏まれると、その直後、実際のクラッチすべり率Rslipが目標クラッチすべり率Rslip_cmdに対してオーバーシュートすることで、すべり率偏差Dslipが一時的に急増する。その後、制御の進行に伴って、すべり率偏差Dslipが減少する。そして、時刻t6以降、クラッチ410が直結状態になることで、すべり率偏差Dslipが値0に保持される。
以上のように、第5実施形態の制御装置1Dによれば、状態予測器520において、主軸回転数NMと制御入力Uactの関係を定義した制御対象モデル(式(105))を用いて、むだ時間d”=0〜3がそれぞれ経過したタイミングでの主軸回転数NMとして、第0〜第3予測値PRE_NM_0〜3が算出され、油温Toilに応じて、4つの重み関数値Wd1”〜Wd4”が算出される。そして、予測当量比PRE_NMが、重み関数値Wdi”と予測値PRE_NM_4−i(i=1〜4)との乗算値の総和として算出されるので、予測当量比PRE_NMを、予測値PRE_NM_4−iを連続的に結合させた値として算出することができる。それにより、油温Toilの変化に応じてむだ時間d”が変化した場合でも、予測主軸回転数NMを、そのようなむだ時間d”の変化を適切に補償しながら精度よく算出することができる。
また、オンボード同定器530において、前述した同定アルゴリズムにより、同定値αid”が算出されるので、前述した同定条件1,2の双方を満たしながら、同定値αid”を算出することができる。具体的には、合成信号値W_act”および推定合成信号値W_hat”が互いに一致するように、同定値αid”が算出されるので、同定条件1すなわち拘束条件を満たしながら、同定値αid”を算出することができる。さらに、修正後制御入力Uact_modが、4つの重み関数値Wd4”〜Wd1”を、むだ時間d”=0〜3前のタイミングでの制御入力Uact(k),Uact(k−1),Uact(k−2),Uact(k−3)にそれぞれ乗算した乗算値の総和として算出されるので、油温Toilの変化に応じてむだ時間d”が連続的に変化する場合でも、修正後制御入力Uact_modを、そのようなむだ時間d”の変化を適切に補償しながら精度よく算出することができる。
さらに、そのような修正後制御入力Uact_modを用いた式(17)〜(29)の同定アルゴリズムにより、同定値αid”がオンボード同定されるので、油温Toilの変化に応じてむだ時間d”が変化したときでも、その影響を抑制しながら、同定値αid”を精度よく同定することができる。特に、油温Toilの急変に伴ってむだ時間d”が急変したときでも、同定値αid”を、むだ時間d”の急変を適切に補償しながら、段差なく円滑に変化するように算出することができる。そして、そのように算出された同定値αid”を用いて、制御入力Uactが算出されるので、クラッチ制御における制御性や、個体間のばらつきおよび経年変化の影響に対するクラッチ制御のロバスト性を飛躍的に向上させることができる。
これに加えて、周波数整形コントローラ540の場合、第1実施形態の周波数整形コントローラ130と同じ手法によって導出された式(125),(126)を用いて、制御入力Uactが算出されるので、前述した制御条件φを満たしながら、制御入力Uactを算出することができる。これに加えて、制御対象モデルのモデルパラメータとして前述した同定値αid”を用いているので、むだ時間d”の変化を適切に補償しながら、制御装置1Dにおける外乱抑制特性やロバスト性を周波数軸上でダイレクトに指定(設定)することができる。それにより、外乱による主軸回転数NMの変動を抑制したい周波数域において、外乱抑制能力やロバスト性を飛躍的に向上させることができる。これに加えて、制御入力Uactの算出アルゴリズムとして、フィードバック制御アルゴリズムが用いられているので、フィードバックゲインを高い状態に維持することができ、高精度かつ高応答性を確保しながら、主軸回転数NMを目標主軸回転数NMに追従させることができる。
なお、第5実施形態は、重み関数値として、互いに重なり合う領域における油温Toilのいずれかの値に対応する重み関数値Wdi”の値の和が、各重み関数値Wdi”における最大値1と等しくなるように設定されているものを用いた例であるが、本発明の重み関数値はこれに限らず、互いに重なり合う領域における参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が所定値になるように設定されているものであればよい。例えば、互いに重なり合う領域における参照パラメータのいずれかの値に対応する重み関数値の総和の絶対値が、重み関数値における絶対値の最大値と等しくなるように設定されたものを用いてもよい。より具体的には、重み関数値として、図41の重み関数値Wdi”の設定値に対して、x軸を対称軸として線対称に設定されたもの、すなわち負値に設定されたものを用いてもよく、その場合には、4つの重み関数値に乗算する値、すなわち4つの予測値PRE_NM_4−iや、4つの制御入力Uact(k−4+i)として、負値化した値を用いればよい。
次に、図48を参照しながら、本発明の第6実施形態に係る制御装置1Eについて説明する。この制御装置1Eは、第5実施形態の制御装置1Dと同様に、自動変速機400におけるクラッチの接続・遮断動作などを制御するものである。この第6実施形態の場合、その機械的な構成は、クラッチとして、乾式クラッチタイプのクラッチ410に代えて、湿式クラッチ(図示せず)を用いている点以外は、第5実施形態のものと同一であるので、以下の説明では、第5実施形態と同じ構成については、同じ符号を付すとともに、その説明は省略する。
ここで、一般に、湿式クラッチの場合、構造上の理由から、乾式クラッチと比べてジャダが発生しにくいという特性を有している。そのため、前述した制御条件φを考慮することなく、クラッチの上流側回転数NEと下流側回転数NMとの間の回転差が時系列的に円滑に値0に収束するように、クラッチを制御すればよいことになる。以上の理由により、本実施形態の制御装置1Eでは、以下に述べる制御アルゴリズムを用いて、制御入力Uactが算出される。
この制御装置1Eは、図48に示すように、クラッチコントローラ700を備えている。このクラッチコントローラ700の場合、前述した図39のクラッチコントローラ500と比較すると、適応外乱オブザーバ740(外乱推定値算出手段)を備えている点と、前述した周波数整形コントローラ540に代えて2自由度応答指定型コントローラ750(制御入力算出手段)を備えている点のみが異なっているので、以下、これらの異なる点についてのみ説明する。
最初に、適応外乱オブザーバ740について説明する。この適応外乱オブザーバ740では、以下に述べる制御アルゴリズムにより、外乱推定値ε”が算出される。まず、下式(127)により、外乱推定用の推定主軸回転数NM_adv(推定制御量)を算出する。
この式(127)は、前述した式(105)において、NM(k+1)をNM_adv(k)に、α”をαid”(k)に、Uact(k−d”)をUact_mod(k)にそれぞれ置き換えるとともに、右辺に外乱推定値ε”を加えたものに相当する。
次いで、下式(128)により、外乱推定用の追従誤差e_adv”を算出する。
そして、最終的に、下式(129)により、外乱推定値ε”が算出される。
この式(129)のπ”は外乱推定ゲインであり、π”>0が成立するように設定される。
次に、前述した2自由度応答指定型コントローラ750について説明する。この2自由度応答指定型コントローラ750では、以下に述べるように、上記外乱推定値ε”を加味した応答指定型制御アルゴリズムよって、制御入力Uactが算出される。
具体的には、まず、下式(130)により、目標主軸回転数のフィルタリング値NM_cmd_fを算出する。
ここで、POLE_f”は目標値フィルタ設定パラメータであり、−1<POLE_f”<0の関係が成立するように設定される。
次いで、下式(131)により、予測追従誤差PRE_e_f”を算出する。
さらに、下式(132)により、切換関数σ_f”を算出する。
ここで、POLE”は切換関数設定パラメータであり、−1<POLE”<0の関係が成立するように設定される。
次いで、下式(133)により、等価制御入力Ueq_f”を算出する。
また、下式(134)により、到達則入力Urch_f”を算出する。
ここで、Krch”は所定のフィードバックゲインである。
そして、最終的に、下式(135)により、制御入力Uactが算出される。
以上の第6実施形態の制御装置1Eによれば、第5実施形態の制御装置1Dと同一の状態予測器520およびオンボード同定器530を備えているので、これらによって、第5実施形態の制御装置1Dと同じ作用効果を得ることができる。これに加えて、適応外乱オブザーバ740において、以上の制御アルゴリズムにより外乱推定値ε”が算出され、2自由度応答指定型コントローラ750において、この外乱推定値ε”を用いて制御入力Uactが算出されるので、クラッチ制御における外乱抑制能力すなわちロバスト性を向上させることができる。
また、適応外乱オブザーバ740を備えているので、例えば、π”>P0”に設定し、オンボード同定器530の同定速度を落とすことによって、制御の安定性を向上させることが可能になる。さらに、同じ理由によって、制御系の共振を防止したり、同定値αid”の演算結果を適用した制御対象モデルのゲイン特性が過小になるのを防止するために、同定値αid”や同定アルゴリズムで用いる入出力データにフィルタリングを施すことができ、より高い制御性を確保することができる。
なお、第1〜第4実施形態は、本発明の制御装置を制御対象としてのエンジン3の空燃比を制御するものに適用し、第5,6実施形態は、本発明の制御装置を制御対象としてのクラッチ410を制御するものに適用した例であるが、本発明の制御装置はこれらに限らず、参照パラメータに応じてむだ時間を含む動特性が変化する特性を備えた制御対象を制御するものに適用可能である。例えば、本発明の制御装置を、制御対象としてのロボットの動作を制御するものに適用してもよい。
また、各実施形態は、本発明の制御装置をむだ時間が4個の整数値(値0〜3)の間で変化する特性の制御対象に適用した例であるが、本発明の制御装置はこれに限らず、むだ時間がM個の整数値の間で変化する特性の制御対象に適用可能である。例えば、本発明を制御装置を、むだ時間が3個以下や5個以上の整数値の間で変化する制御対象に適用してもよい。