(57)【要約】【課題】むだ時間の変化を逐次推定して高精度な制御を実施可能にしつつ、実際のむだ時間が急変した場合であってもその急変をむだ時間推定値に迅速に反映できる制御装置を提供する。【解決手段】所定のサンプリング周期ｄｔで離散化された数式により制御対象を表したプラントモデルのパラメータを逐次同定し、同定されたパラメータを用いてプラントモデルに含まれている余むだ時間を推定し、推定した余むだ時間に基づきむだサンプリング回数の値を繰り上げ／繰り下げ更新し、更新されたむだサンプリング回数に基づきむだ時間を推定する。そして、余むだ時間を＋ｄｔと比較して前記更新を実施する処理と、＋２ｄｔと比較して前記更新を実施する処理とを同時に並行して実行する。【選択図】      図１４
2010019107.tif 000002
【特許請求の範囲】
【請求項１】
所定のサンプリング周期（ｄｔ）で離散化された数式により制御対象を表したプラントモデルを対象とし、そのプラントモデルのパラメータを逐次同定する同定手段と、
前記制御対象のむだ時間（Ｌ）を前記サンプリング周期で除算した商をむだサンプリング回数（ｄ）、余りを余むだ時間（Ｌ１）としてモデル化した場合に、前記同定手段により同定された同定パラメータを用いて前記プラントモデルに含まれている前記余むだ時間を推定する余むだ時間推定手段と、
前記余むだ時間推定手段により推定された推定余むだ時間が、所定の増加条件を満たす程度に増加した場合には前記むだサンプリング回数の推定値を一つ繰り上げるよう更新し、所定の減少条件を満たす程度に減少した場合には前記むだサンプリング回数の推定値を一つ繰り下げるよう更新する更新手段と、
前記むだサンプリング回数更新手段により更新された前記むだサンプリング回数の推定値に基づき前記むだ時間を推定するむだ時間推定手段と、
を備え、
前記更新手段は、
推定された前記余むだ時間を、前記サンプリング周期の１倍の長さ（＋ｄｔ）との比較に基づき前記増加条件を満たすか否かを判定する第１増加判定手段と、
推定された前記余むだ時間を前記サンプリング周期の２倍の長さ（＋２ｄｔ）との比較に基づき前記増加条件を満たすか否かを判定する第２増加判定手段と、を有するとともに、
前記第１増加判定手段による判定及び前記第２増加判定手段による判定を同時に並行して実行することを特徴とする制御装置。
【請求項２】
前記第１増加判定手段は、推定された前記余むだ時間の変動推移のうち前記サンプリング周期の１倍長さ（＋ｄｔ）及び２倍長さ（＋２ｄｔ）の範囲で囲まれた部分の面積に相当する積算値が所定の判定閾値を超えた時点で前記増加条件を満たすと判定し、
前記第２増加判定手段は、推定された前記余むだ時間の変動推移のうち前記サンプリング周期の２倍長さ（＋２ｄｔ）及び３倍長さ（＋３ｄｔ）の範囲で囲まれた部分の面積に相当する積算値が所定の判定閾値を超えた時点で前記増加条件を満たすと判定することを特徴とする請求項１に記載の制御装置。
【請求項３】
前記第２増加判定手段にて用いられる判定閾値を、前記第１増加判定手段にて用いられる判定閾値よりも絶対値が大きい値に設定することを特徴とする請求項２に記載の制御装置。
【請求項４】
前記積算値は、前記推定余むだ時間を前記サンプリング時間で除算して無次元化した値の積算値であることを特徴とする請求項２又は３に記載の制御装置。
【請求項５】
前記積算値は、前記推定余むだ時間に相当する値を奇数乗して得られた値を積算した値であることを特徴とする請求項２〜４のいずれか１つに記載の制御装置。
【請求項６】
前記第１増加判定手段における前記積算値は、前記推定余むだ時間を前記サンプリング時間で除算し１を減算した値であって、０以上で１より小さい所定下限値と該下限値より大きく１以下の所定上限値により制限した値であることを特徴とする請求項２〜５のいずれか１つに記載の制御装置。
【請求項７】
前記第２増加判定手段における前記積算値は、前記推定余むだ時間を前記サンプリング時間で除算し２を減算した値であって、０以上で１より小さい所定下限値と該下限値より大きく１以下の所定上限値により制限した値であることを特徴とする請求項２〜６のいずれか１つに記載の制御装置。
【請求項８】
前記更新手段は、推定された前記余むだ時間を前記サンプリング周期の３倍の長さ（＋３ｄｔ）との比較に基づき前記増加条件を満たすか否かを判定する第３増加判定手段を有するとともに、前記第１及び第２増加判定手段による判定及び前記第３増加判定手段による判定を同時に並行して実行することを特徴とする請求項１〜７のいずれか１つに記載の制御装置。
【請求項９】
前記更新手段は、
推定された前記余むだ時間をゼロとの比較に基づき前記減少条件を満たすか否かを判定する第１減少判定手段と、
推定された前記余むだ時間を前記サンプリング周期の−１倍の長さ（−ｄｔ）との比較に基づき前記減少条件を満たすか否かを判定する第２減少判定手段と、を有するとともに、
前記第１減少判定手段による判定及び前記第２減少判定手段による判定を同時に並行して実行することを特徴とする請求項１〜８のいずれか１つに記載の制御装置。
【請求項１０】
所定のサンプリング周期（ｄｔ）で離散化された数式により制御対象を表したプラントモデルを対象とし、そのプラントモデルのパラメータを逐次同定する同定手段と、
前記制御対象のむだ時間（Ｌ）を前記サンプリング周期で除算した商をむだサンプリング回数（ｄ）、余りを余むだ時間（Ｌ１）としてモデル化した場合に、前記同定手段により同定された同定パラメータを用いて前記プラントモデルに含まれている前記余むだ時間を推定する余むだ時間推定手段と、
前記余むだ時間推定手段により推定された推定余むだ時間が、所定の増加条件を満たす程度に増加した場合には前記むだサンプリング回数の推定値を一つ繰り上げるよう更新し、所定の減少条件を満たす程度に減少した場合には前記むだサンプリング回数の推定値を一つ繰り下げるよう更新する更新手段と、
前記むだサンプリング回数更新手段により更新された前記むだサンプリング回数の推定値に基づき前記むだ時間を推定するむだ時間推定手段と、
を備え、
前記更新手段は、
推定された前記余むだ時間をゼロとの比較に基づき前記減少条件を満たすか否かを判定する第１減少判定手段と、
推定された前記余むだ時間を前記サンプリング周期の−１倍の長さ（−ｄｔ）との比較に基づき前記減少条件を満たすか否かを判定する第２減少判定手段と、を有するとともに、
前記第１減少判定手段による判定及び前記第２減少判定手段による判定を同時に並行して実行することを特徴とする制御装置。
【請求項１１】
前記第１減少判定手段は、推定された前記余むだ時間の変動推移のうち前記サンプリング周期の−１倍長さ（−ｄｔ）及びゼロの範囲で囲まれた部分の面積に相当する積算値が所定の判定閾値を超えた時点で前記減少条件を満たすと判定し、
前記第２減少判定手段は、推定された前記余むだ時間の変動推移のうち前記サンプリング周期の−２倍長さ（−２ｄｔ）及び−１倍長さ（−ｄｔ）の範囲で囲まれた部分の面積に相当する積算値が所定の判定閾値を超えた時点で前記減少条件を満たすと判定することを特徴とする請求項９又は１０に記載の制御装置。
【請求項１２】
前記第２減少判定手段にて用いられる判定閾値を、前記第１減少判定手段にて用いられる判定閾値よりも絶対値が大きい値に設定することを特徴とする請求項１１に記載の制御装置。
【請求項１３】
前記積算値は、前記推定余むだ時間を前記サンプリング時間で除算して無次元化した値の積算値であることを特徴とする請求項１１又は１２に記載の制御装置。
【請求項１４】
前記積算値は、前記推定余むだ時間に相当する値を奇数乗して得られた値を積算した値であることを特徴とする請求項１１〜１３のいずれか１つに記載の制御装置。
【請求項１５】
前記第１減少判定手段における前記積算値は、前記推定余むだ時間を前記サンプリング時間で除算した値であって、０以下で−１より大きい所定上限値と該上限値より小さく−１以上の所定下限値により制限した値であることを特徴とする請求項１１〜１４のいずれか１つに記載の制御装置。
【請求項１６】
前記第２減少判定手段における前記積算値は、前記推定余むだ時間を前記サンプリング時間で除算し１を加算した値であって、０以下で−１より大きい所定上限値と該上限値より小さく−１以上の所定下限値により制限した値であることを特徴とする請求項１１〜１５のいずれか１つに記載の制御装置。
【請求項１７】
前記更新手段は、推定された前記余むだ時間を前記サンプリング周期の−２倍の長さ（−２ｄｔ）との比較に基づき前記減少条件を満たすか否かを判定する第３減少判定手段を有するとともに、前記第１及び第２減少判定手段による判定及び前記第３減少判定手段による判定を同時に並行して実行することを特徴とする請求項９〜１６のいずれか１つに記載の制御装置。
【請求項１８】
前記制御対象は、内燃機関の排気通路に配設された排ガス浄化用の触媒と、前記排気通路において前記触媒の上流側で特定排ガス濃度を検出する上流側排ガスセンサと、前記排気通路において前記触媒の下流側で特定排ガス濃度を検出する下流側排ガスセンサと、を有する空燃比制御システムであり、
前記上流側排ガスセンサで検出した空燃比を前記触媒の上流側目標空燃比に一致させるようにフィードバックして、供給空燃比を前記制御対象への制御入力として演算するメインフィードバック制御手段と、
前記下流側排ガスセンサの検出した空燃比を前記触媒の下流側目標空燃比に一致させるようにフィードバックして、前記上流側目標空燃比を前記制御入力の算出に用いる値として演算するサブフィードバック制御手段と、
を備えることを特徴とする請求項１〜１７のいずれか１つに記載の制御装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、内燃機関等の制御対象を数式で表したプラントモデルを用い、該プラントモデルのパラメータを逐次同定することで制御対象の動特性を自動的に適応させるようにした制御装置に関する発明である。
【背景技術】
【０００２】
従来から、制御対象の動特性を逐次同定する制御装置が提案されている。この制御装置は、制御対象を数式で表したプラントモデルを持ち、実際の制御入力をプラントモデルに入力した時のプラントモデル出力と実際の制御対象の出力との誤差をゼロに近づけるようにプラントモデルのパラメータを逐次推定するものであるが、制御対象がむだ時間を含む系の場合、同定に用いるプラントモデルは、実際の制御対象が持つむだ時間を考慮する必要がある。そこで特許文献１，２では、以下の同定手法によりむだ時間の推定を可能にしている。
【０００３】
すなわち、所定のサンプリング周期（ｄｔ）で離散化された数式により制御対象を表したプラントモデルについて、そのプラントモデルのむだ時間（Ｌ）をサンプリング周期で除算した商をむだサンプリング回数（ｄ）、余りを余むだ時間（Ｌ１）とする。そして、逐次同定されたパラメータを用いれば、プラントモデルに含まれている余むだ時間を推定することができることを利用して、以下に説明する図８のタイムチャートに例示される如くむだ時間を推定する。
【０００４】
図８（ａ）に示すように、ある制御入力に対するノミナルむだ時間と実むだ時間との間に図示のような誤差が生じている場合を想定する。なお、サンプリングタイミング（ｄ・ｄｔ）を超えた部分が余むだ時間Ｌ１であり、ノミナルむだ時間は、予め設定されたむだ時間のノミナル値である。
【０００５】
この場合、同定が逐次行われることにより、図８（ｂ）に示すように、余むだ時間の推定値Ｌ１＿hat（以下、ｘ＿hatはｘの同定値もしくは推定値を表す）が変化し、ｄｔ近傍まで増加する。そして、余むだ時間の推定値Ｌ１＿hatがｄｔ近傍で所定時間以上とどまっている場合には、図８（ｃ）に示すように、むだサンプリング回数の推定値ｄ＿hatが１つ繰り上げられ（ｄ＿hat←ｄ＋１）、更に、余むだ時間の推定値Ｌ１＿hatがゼロ近傍の値に更新される。その後、同定が逐次行われることにより、図８（ｄ）に示すように、余むだ時間の推定値Ｌ１＿hatが増加し、それに伴いむだ時間の推定値が実むだ時間に収束する。
【０００６】
このようにして更新されたむだサンプリング回数ｄ＿hat及び余むだ時間Ｌ１＿hatを用いれば、Ｌ（むだ時間）＝ｄ（むだサンプリング回数）×ｄｔ（サンプリング周期）＋Ｌ１（余むだ時間）、との算出式に基づきむだ時間Ｌを推定することができる。
【０００７】
なお、図８では、余むだ時間の推定値Ｌ１＿hatが増加する場合を例示し、それに伴い、むだサンプリング回数の推定値ｄ＿hatが１つ繰り上げられる様子を説明したが、これとは逆に、余むだ時間の推定値Ｌ１＿hatが減少してゼロ近傍で所定時間以上とどまっている場合には、むだサンプリング回数の推定値ｄ＿hatが１つ繰り下げられると共に、余むだ時間の推定値Ｌ１＿hatがｄｔ近傍の値に更新される。
【特許文献１】特開２００６−１１８４２９号公報
【特許文献２】特開２００７−２８５１６７号公報
【発明の開示】
【発明が解決しようとする課題】
【０００８】
しかしながら、制御対象の実際の特性が急激に変化した場合や同定開始時等には、図１７（ｂ）に例示されるように余むだ時間の推定値Ｌ１＿hatが短時間で急激に変化することとなる。しかしながら上記従来の制御では、余むだ時間の推定値Ｌ１＿hatがｄｔ近傍又はゼロ近傍で所定時間Ｔ１０以上とどまっていることを条件としてむだサンプリング回数の推定値ｄ＿hatが更新されるので、急激に変化する余むだ時間に応じたサンプリング回数の更新が速やかに為されず、その更新速度が遅い。
【０００９】
本発明は、上記課題を解決するためになされたものであり、その目的は、むだ時間の変化を逐次推定して高精度な制御を実施可能にしつつ、実際のむだ時間が急変した場合であってもその急変をむだ時間推定値に迅速に反映できる制御装置を提供することにある。
【課題を解決するための手段】
【００１０】
以下、上記課題を解決するための手段、及びその作用効果について記載する。なお、以下に説明する請求項１〜請求項８記載の発明は、実際のむだ時間が増加する向きに急変した場合に対する発明であり、請求項９〜請求項１７記載の発明は、実際のむだ時間が減少する向きに急変した場合に対する発明である。
【００１１】
請求項１記載の発明では、
所定のサンプリング周期（ｄｔ）で離散化された数式により制御対象を表したプラントモデルを対象とし、そのプラントモデルのパラメータを逐次同定する同定手段と、
前記制御対象のむだ時間（Ｌ）を前記サンプリング周期で除算した商をむだサンプリング回数（ｄ）、余りを余むだ時間（Ｌ１）としてモデル化した場合に、前記同定手段により同定された同定パラメータを用いて前記プラントモデルに含まれている前記余むだ時間を推定する余むだ時間推定手段と、
前記余むだ時間推定手段により推定された推定余むだ時間が、所定の増加条件を満たす程度に増加した場合には前記むだサンプリング回数の推定値を一つ繰り上げるよう更新し、所定の減少条件を満たす程度に減少した場合には前記むだサンプリング回数の推定値を一つ繰り下げるよう更新する更新手段と、
前記むだサンプリング回数更新手段により更新された前記むだサンプリング回数の推定値に基づき前記むだ時間を推定するむだ時間推定手段と、
を備え、
前記更新手段は、
推定された前記余むだ時間を、前記サンプリング周期の１倍の長さ（＋ｄｔ）との比較に基づき前記増加条件を満たすか否かを判定する第１増加判定手段と、
推定された前記余むだ時間を前記サンプリング周期の２倍の長さ（＋２ｄｔ）との比較に基づき前記増加条件を満たすか否かを判定する第２増加判定手段と、を有するとともに、
前記第１増加判定手段による判定及び前記第２増加判定手段による判定を同時に並行して実行することを特徴とする。
【００１２】
これによれば、離散化されたプラントモデルのパラメータを逐次同定し（同定手段）、同定されたパラメータを用いてプラントモデルに含まれている余むだ時間を推定し（余むだ時間推定手段）、推定された余むだ時間がサンプリング周期近傍へ変化した場合にはむだサンプリング回数の推定値を一つ繰り上げるよう更新し、ゼロ近傍へ変化した場合には一つ繰り下げるよう更新し（むだサンプリング回数更新手段）、更新されたむだサンプリング回数の推定値に基づきむだ時間を推定する（むだ時間推定手段）。よって、むだ時間の変化を逐次推定して高精度な制御を実施できる。
【００１３】
しかも、第１増加判定手段による判定と第２増加判定手段による判定とを同時に並行して実行するので、第１増加判定手段による判定が終了することを待たずして第２増加判定手段による判定が開始されることとなる。そのため、例えば、推定された余むだ時間が短時間で急激に増加した場合、つまり＋ｄｔとの比較に基づき増加条件を満たすとともに、＋２ｄｔとの比較に基づいても増加条件を満たすほどに急増した場合において、第２増加判定手段による増加条件を満たすとの判定結果を早期に取得できる。
【００１４】
よって、更新手段は、第１増加判定手段による増加判定結果に基づきむだサンプリング回数の推定値を一つ繰り上げた後即座に（或いはその繰り上げと同時に）、第２増加判定手段による増加判定結果に基づきむだサンプリング回数の推定値をさらに一つ繰り上げることができる。よって、急激に変化する余むだ時間に応じたむだサンプリング回数推定値の更新を速やかにできる。したがって、むだ時間の変化を逐次推定して高精度な制御を実施可能にしつつ、実際のむだ時間が急変した場合であってもその急変をむだ時間推定値に迅速に反映できる。
【００１５】
請求項２記載の発明では、前記第１増加判定手段は、推定された前記余むだ時間の変動推移のうち前記サンプリング周期の１倍長さ（＋ｄｔ）及び２倍長さ（＋２ｄｔ）の範囲で囲まれた部分の面積に相当する積算値が所定の判定閾値を超えた時点で前記増加条件を満たすと判定し、前記第２増加判定手段は、推定された前記余むだ時間の変動推移のうち前記サンプリング周期の２倍長さ（＋２ｄｔ）及び３倍長さ（＋３ｄｔ）の範囲で囲まれた部分の面積に相当する積算値が所定の判定閾値を超えた時点で前記増加条件を満たすと判定することを特徴とする。
【００１６】
これによれば、余むだ時間の推定値がノイズ等により瞬時的に大きく変動した場合であっても、当該変動により判定結果が受ける影響を小さくできるので、第１及び第２増加判定手段がノイズ等の影響を受けて誤判定することを抑制できる。
【００１７】
ちなみに、上記請求項１記載の発明を実施するにあたり、第１増加判定手段は、推定された余むだ時間が＋ｄｔ近傍でとどまった状態又は＋ｄｔを超えた状態が所定時間以上継続した場合に、増加条件を満たすと判定するようにしてもよい。第２増加判定手段についても同様にして、推定された余むだ時間が＋２ｄｔ近傍でとどまった状態又は＋２ｄｔを超えた状態が所定時間以上継続した場合に、増加条件を満たすと判定するようにしてもよい。但し、上記請求項２記載の発明に比べれば、余むだ時間推定値の変動による影響を受け易い。つまり、請求項２記載の発明の方が判定精度の点で勝る。
【００１８】
ここで、第２増加判定手段にて用いられる判定閾値を第１増加判定手段にて用いられる判定閾値と同じ値にすると、余むだ時間推定値変動により第２増加判定手段の判定結果が受ける影響は、第１増加判定手段の判定結果が受ける影響に比べて大きい。
【００１９】
これに対し請求項３記載の発明では、前記第２増加判定手段にて用いられる判定閾値を、前記第１増加判定手段にて用いられる判定閾値よりも絶対値が大きい値に設定することを特徴とする。
【００２０】
よって、余むだ時間推定値のノイズ等による変動から第２増加判定手段の判定結果が受ける影響を好適に小さくできる。換言すれば、第１増加判定手段にて用いられる判定閾値を、前記第２増加判定手段にて用いられる判定閾値よりも小さい値に設定することで、第１増加判定手段による増加条件を満たすとの判定結果を取得できる時期を早くできる。つまり、第１増加判定手段による判定結果を早期に取得できるとともに、第２増加判定手段がノイズ等の影響を受けて誤判定することを抑制できる。
【００２１】
請求項４記載の発明では、前記積算値は、前記推定余むだ時間を前記サンプリング時間で除算して無次元化した値の積算値であることを特徴とする。
【００２２】
これによれば、前記推定余りむだ時間は前記サンプリング時間で除算することにより無次元化されるため、前記サンプリング時間が変動する場合においても、前記積算値を好適に算出することができる。
【００２３】
請求項５記載の発明では、前記積算値は、前記推定余むだ時間に相当する値を奇数乗して得られた値を積算した値であることを特徴とする。
【００２４】
これによれば、積算する前に奇数乗（例えば３乗）することによって符号を保持しつつ前記余りむだ時間が大きいほど相対的に積算値が大きくなるような重み付けができ、前記第１および第２増加判定精度を向上させることができる。
【００２５】
請求項６記載の発明では、前記第１増加判定手段における前記積算値は、前記推定余むだ時間を前記サンプリング時間で除算し１を減算した値であって、０以上で１より小さい所定下限値と該下限値より大きく１以下の所定上限値により制限した値であることを特徴とする。また、請求項７記載の発明では、前記第２増加判定手段における前記積算値は、前記推定余むだ時間を前記サンプリング時間で除算し２を減算した値であって、０以上で１より小さい所定下限値と該下限値より大きく１以下の所定上限値により制限した値であることを特徴とする。
【００２６】
これらによれば、上限値及び下限値により制限がかけられた値について積算するので、ノイズ等の影響を受けて誤判定することを抑制できる。
【００２７】
請求項８記載の発明では、前記更新手段は、推定された前記余むだ時間を前記サンプリング周期の３倍の長さ（＋３ｄｔ）との比較に基づき前記増加条件を満たすか否かを判定する第３増加判定手段を有するとともに、前記第１及び第２増加判定手段による判定及び前記第３増加判定手段による判定を同時に並行して実行することを特徴とする。
【００２８】
これによれば、第１及び第２増加判定手段に加え第３増加判定手段による判定をも同時に並行して実行するので、第１及び第２増加判定手段による判定が終了することを待たずして第３増加判定手段による判定が開始されることとなる。そのため、余むだ時間推定値が短時間で急激に増加した場合において、第３増加判定手段による増加条件を満たすとの判定結果を早期に取得できる。よって、更新手段は、第１及び第２増加判定手段による増加判定結果に基づきむだサンプリング回数の推定値を繰り上げた後即座に（或いはその繰り上げと同時に）、第３増加判定手段による増加判定結果に基づきむだサンプリング回数の推定値をさらに一つ繰り上げることができる。よって、急激に変化する余むだ時間に応じたサンプリング回数の更新を速やかにできる。
【００２９】
ちなみに、上記請求項３記載の思想を適用して、第３増加判定手段にて用いられる判定閾値を、第２増加判定手段にて用いられる判定閾値よりも大きい値に設定することが望ましい。
【００３０】
請求項９記載の発明では、
前記更新手段は、
推定された前記余むだ時間をゼロとの比較に基づき前記減少条件を満たすか否かを判定する第１減少判定手段と、
推定された前記余むだ時間を前記サンプリング周期の−１倍の長さ（−ｄｔ）との比較に基づき前記減少条件を満たすか否かを判定する第２減少判定手段と、を有するとともに、
前記第１減少判定手段による判定及び前記第２減少判定手段による判定を同時に並行して実行することを特徴とする。
【００３１】
これによれば、第１減少判定手段による判定と第２減少判定手段による判定とを同時に並行して実行するので、第１減少判定手段による判定が終了することを待たずして第２減少判定手段による判定が開始されることとなる。そのため、例えば、余むだ時間推定値が短時間で急激に増加した場合、つまりゼロとの比較に基づき減少条件を満たすとともに、−ｄｔとの比較に基づいても減少条件を満たすほどに急減した場合において、第２減少判定手段による減少条件を満たすとの判定結果を早期に取得できる。
【００３２】
よって、更新手段は、第１減少判定手段による減少判定結果に基づきむだサンプリング回数の推定値を一つ繰り下げた後即座に（或いはその繰り下げと同時に）、第２減少判定手段による減少判定結果に基づきむだサンプリング回数の推定値をさらに一つ繰り下げることができる。よって、急激に減少変化する場合においても、余むだ時間に応じたむだサンプリング回数推定値の更新を速やかにできる。
【００３３】
請求項１０記載の発明では、上記請求項１と同様の同定手段、余むだ時間推定手段、更新手段及びむだ時間推定手段を備えるので、請求項１と同様にしてむだ時間の変化を逐次推定して高精度な制御を実施できる。また、請求項１０記載の発明では、前記更新手段は、推定された前記余むだ時間をゼロとの比較に基づき前記減少条件を満たすか否かを判定する第１減少判定手段と、推定された前記余むだ時間を前記サンプリング周期の−１倍の長さ（−ｄｔ）との比較に基づき前記減少条件を満たすか否かを判定する第２減少判定手段と、を有するとともに、前記第１減少判定手段による判定及び前記第２減少判定手段による判定を同時に並行して実行することを特徴とする。
【００３４】
これによれば、第１減少判定手段による判定と第２減少判定手段による判定とを同時に並行して実行するので、第１減少判定手段による判定が終了することを待たずして第２減少判定手段による判定が開始されることとなる。そのため、例えば、余むだ時間推定値が短時間で急激に増加した場合、つまりゼロとの比較に基づき減少条件を満たすとともに、−ｄｔとの比較に基づいても減少条件を満たすほどに急減した場合において、第２減少判定手段による減少条件を満たすとの判定結果を早期に取得できる。
【００３５】
よって、更新手段は、第１減少判定手段による減少判定結果に基づきむだサンプリング回数の推定値値を一つ繰り下げた後即座に（或いはその繰り下げと同時に）、第２減少判定手段による減少判定結果に基づきむだサンプリング回数の推定値をさらに一つ繰り下げることができる。よって、急激に変化する余むだ時間に応じたむだサンプリング回数推定値の更新を速やかにできる。したがって、むだ時間の変化を逐次推定して高精度な制御を実施可能にしつつ、実際のむだ時間が急変した場合であってもその急変をむだ時間推定値に迅速に反映できる。
【００３６】
請求項１１記載の発明では、前記第１減少判定手段は、推定された前記余むだ時間の変動推移のうち前記サンプリング周期の−１倍長さ（−ｄｔ）及びゼロの範囲で囲まれた部分の面積に相当する積算値が所定の判定閾値を超えた時点で前記減少条件を満たすと判定し、前記第２減少判定手段は、推定された前記余むだ時間の変動推移のうち前記サンプリング周期の−２倍長さ（−２ｄｔ）及び−１倍長さ（−ｄｔ）の範囲で囲まれた部分の面積に相当する積算値が所定の判定閾値を超えた時点で前記減少条件を満たすと判定することを特徴とする。
【００３７】
これによれば、余むだ時間の推定値がノイズ等により瞬時的に大きく変動した場合であっても、当該変動により判定結果が受ける影響を小さくできるので、第１及び第２減少判定手段がノイズ等の影響を受けて誤判定することを抑制できる。
【００３８】
ちなみに、上記請求項１０記載の発明を実施するにあたり、第１減少判定手段は、推定された余むだ時間がゼロ近傍でとどまった状態又はゼロを下回った状態が所定時間以上継続した場合に、減少条件を満たすと判定するようにしてもよい。第２減少判定手段についても同様にして、推定された余むだ時間が−ｄｔ近傍でとどまった状態又は−ｄｔを下回った状態が所定時間以上継続した場合に、減少条件を満たすと判定するようにしてもよい。但し、上記請求項１１記載の発明に比べれば、余むだ時間推定値の変動による影響を受け易い。つまり、請求項１１記載の発明の方が判定精度の点で勝る。
【００３９】
ここで、第２減少判定手段にて用いられる判定閾値を第１減少判定手段にて用いられる判定閾値と同じ値にすると、余むだ時間推定値変動により第２減少判定手段の判定結果が受ける影響は、第１減少判定手段の判定結果が受ける影響に比べて大きい。
【００４０】
これに対し請求項１２記載の発明では、前記第２減少判定手段にて用いられる判定閾値を、前記第１減少判定手段にて用いられる判定閾値よりも絶対値が大きい値に設定することを特徴とする。
【００４１】
よって、余むだ時間推定値のノイズ等による変動から第２減少判定手段の判定結果が受ける影響を好適に小さくできる。換言すれば、第１減少判定手段にて用いられる判定閾値を、前記第２減少判定手段にて用いられる判定閾値よりも絶対値が小さい値に設定することで、第１減少判定手段による減少条件を満たすとの判定結果を取得できる時期を早くできる。つまり、第１減少判定手段による判定結果を早期に取得できるとともに、第２減少判定手段がノイズ等の影響を受けて誤判定することを抑制できる。
【００４２】
請求項１３記載の発明では、前記積算値は、前記推定余むだ時間を前記サンプリング時間で除算して無次元化した値の積算値であることを特徴とする。
【００４３】
これによれば、前記推定余りむだ時間は前記サンプリング時間で除算することにより無次元化されるため、前記サンプリング時間が変動する場合においても、前記積算値を好適に算出することができる。
【００４４】
請求項１４記載の発明では、前記積算値は、前記推定余むだ時間に相当する値を奇数乗して得られた値を積算した値であることを特徴とする。
【００４５】
これによれば、積算する前に奇数乗（例えば３乗）することによって符号を保持しつつ前記余りむだ時間が小さい（負で絶対値が大きい）ほど相対的に小さく（負で絶対値が大きく）なるような重み付けができ、前記第１および第２減少判定精度を向上させることができる。
【００４６】
請求項１５記載の発明では、前記第１減少判定手段における前記積算値は、前記推定余むだ時間を前記サンプリング時間で除算した値であって、０以下で−１より大きい所定上限値と該上限値より小さく−１以上の所定下限値により制限した値であることを特徴とする。また、請求項１６記載の発明では、前記第２減少判定手段における前記積算値は、前記推定余むだ時間を前記サンプリング時間で除算し１を加算した値であって、０以下で−１より大きい所定上限値と該上限値より小さく−１以上の所定下限値により制限した値であることを特徴とする。
【００４７】
これらによれば、上限値及び下限値により制限がかけられた値について積算するので、ノイズ等の影響を受けて誤判定することを抑制できる。  請求項１７記載の発明では、前記更新手段は、推定された前記余むだ時間を前記サンプリング周期の−２倍の長さ（−２ｄｔ）との比較に基づき前記減少条件を満たすか否かを判定する第３減少判定手段を有するとともに、前記第１及び第２減少判定手段による判定及び前記第３減少判定手段による判定を同時に並行して実行することを特徴とする。
【００４８】
これによれば、第１及び第２減少判定手段に加え第３減少判定手段による判定をも同時に並行して実行するので、第１及び第２減少判定手段による判定が終了することを待たずして第３減少判定手段による判定が開始されることとなる。そのため、余むだ時間推定値が短時間で急激に減少した場合において、第３減少判定手段による減少条件を満たすとの判定結果を早期に取得できる。よって、更新手段は、第１及び第２減少判定手段による減少判定結果に基づきむだサンプリング回数の値を繰り下げた後即座に（或いはその繰り下げと同時に）、第３減少判定手段による減少判定結果に基づきむだサンプリング回数の値をさらに一つ繰り上げることができる。よって、急激に変化する余むだ時間に応じたサンプリング回数の更新を速やかにできる。
【００４９】
ちなみに、上記請求項９記載の思想を適用して、第３減少判定手段にて用いられる判定閾値を、第２減少判定手段にて用いられる判定閾値よりも大きい値に設定することが望ましい。
【００５０】
また、請求項１８記載の如く、前記制御対象を、以下に説明する空燃比制御システムに用いて好適である。すなわち、
前記制御対象は、内燃機関の排気通路に配設された排ガス浄化用の触媒と、前記排気通路において前記触媒の上流側で特定排ガス濃度を検出する上流側排ガスセンサと、前記排気通路において前記触媒の下流側で特定排ガス濃度を検出する下流側排ガスセンサと、を有する空燃比制御システムであり、
前記上流側排ガスセンサで検出した空燃比を前記触媒の上流側目標空燃比に一致させるようにフィードバックして、供給空燃比を前記制御対象への制御入力として演算するメインフィードバック制御手段と、前記下流側排ガスセンサの検出した空燃比を前記触媒の下流側目標空燃比に一致させるようにフィードバックして、前記上流側目標空燃比を前記制御入力の算出に用いる値として演算するサブフィードバック制御手段と、を備えることを特徴とする。
【発明を実施するための最良の形態】
【００５１】
以下、本発明を内燃機関の空燃比制御システムに適用して具体化した一実施形態について図面を用いて説明する。
【００５２】
まず、図１に基づいてエンジン制御システム全体の概略構成を説明する。内燃機関であるエンジン１１の吸気管１２の最上流部には、エアクリーナ１３が設けられ、このエアクリーナ１３の下流側に、吸入空気量を検出するエアフローメータ１４が設けられている。このエアフローメータ１４の下流側には、モータ１５によって開度調節されるスロットルバルブ１６と、このスロットルバルブ１６の開度（スロットル開度）を検出するスロットル開度センサ１７とが設けられている。
【００５３】
更に、スロットルバルブ１６の下流側には、サージタンク１８が設けられ、このサージタンク１８には、吸気管圧力を検出する吸気管圧力センサ１９が設けられている。また、サージタンク１８には、エンジン１１の各気筒に空気を導入する吸気マニホールド２０が設けられ、各気筒の吸気マニホールド２０の吸気ポート近傍に、それぞれ燃料を噴射する燃料噴射弁２１が取り付けられている。また、エンジン１１のシリンダヘッドには、各気筒毎に点火プラグ２２が取り付けられ、各点火プラグ２２の火花放電によって各気筒内の混合気に着火される。
【００５４】
また、エンジン１１のシリンダブロックには、冷却水温を検出する冷却水温センサ２６や、エンジン１１のクランク軸２７が所定クランク角回転する毎にクランク角信号（パルス信号）を出力するクランク角センサ２８が取り付けられている。このクランク角センサ２８のクランク角信号に基づいてクランク角やエンジン回転速度が検出される。
【００５５】
一方、エンジン１１の排気管２３には、排気ガスを浄化するための２つの触媒２５，３０が直列に設けられている。各触媒２５，３０は、例えば、三元触媒、ＮＯｘ吸蔵型三元触媒等により構成され、上流側の触媒２５の上流側と下流側には、それぞれ特定排気ガス濃度（例えば酸素濃度、空燃比等）を検出する上流側排気ガスセンサ３１と下流側排気ガスセンサ３２とが設けられている。本実施形態では、上流側排気ガスセンサ３１として空燃比センサを用い、下流側排気ガスセンサ３２として酸素センサ（Ｏ２センサ）を用いているが、この構成に限定されないことは言うまでもない。
【００５６】
これら各種センサの出力は、エンジン制御回路（以下「ＥＣＵ」と表記する）２９に入力される。このＥＣＵ２９は、マイクロコンピュータを主体として構成され、内蔵されたＲＯＭ（記憶媒体）に記憶された各種のエンジン制御プログラムを実行することで、エンジン運転状態に応じて燃料噴射弁２１の燃料噴射量や点火プラグ２２の点火時期を制御する。
【００５７】
また、ＥＣＵ２９は、図２に示すように、上流側排気ガスセンサ３１で検出した空燃比を上流側触媒２５上流側の目標空燃比に一致させるように供給空燃比（燃料噴射量）をフィードバック補正するメインフィードバック制御（以下の説明では「フィードバック」を「Ｆ／Ｂ」で表記する）を行うメインＦ／ＢコントローラＭ１１と、下流側排気ガスセンサ３２の検出電圧（検出空燃比）を目標電圧（上流側触媒２５下流側の目標空燃比）に一致させるように上流側触媒２５上流側の目標空燃比をＦ／Ｂ補正するサブＦ／Ｂ制御を行うサブＦ／ＢコントローラＭ２１として機能する。
【００５８】
メインＦ／ＢコントローラＭ１１は、事前にモデル化して適合したプラントモデル（１次遅れ系＋むだ時間のモデル）を基に設計されたものであり、このメインＦ／ＢコントローラＭ１１により上流側排気ガスセンサ３１の検出空燃比と上流側触媒２５上流側の目標空燃比との偏差を小さくするように空燃比補正係数が演算され、最適な空燃比Ｆ／Ｂ制御が実現される。
【００５９】
一方、サブＦ／ＢコントローラＭ２１は、下流側排気ガスセンサ３２の検出電圧（検出空燃比）を目標電圧（上流側触媒２５下流側の目標空燃比）に一致させるように上流側触媒２５上流側の目標空燃比をＦ／Ｂ補正するサブＦ／Ｂ制御を行うように、事前にモデル化して適合したプラントモデル（２次遅れ系＋むだ時間のモデル）を基に設計されている。
【００６０】
しかしながら、このサブＦ／ＢコントローラＭ２１の実際の制御対象（上流側触媒２５、下流側排気ガスセンサ３２）の個体差や劣化等によってサブＦ／Ｂ制御誤差が生じる。そこで、本実施形態では、適応制御と称される制御方式を用い、サブＦ／ＢコントローラＭ２１におけるＦ／Ｂゲインを制御対象（プラント）の現時点の動特性に自動的に適応させ、制御系の性能を常に最良の状態に保持するようにしている。すなわち、制御対象を数式モデルで表したプラントモデルを用い、制御対象の入力（上流側触媒２５の流入過不足酸素量）を該プラントモデルに入力した時のプラントモデル出力と制御対象の出力（下流側排気ガスセンサ３２の検出電圧）との誤差ｅをゼロに近づけるようにプラントモデルの可変パラメータを逐次同定する（この機能が特許請求の範囲でいう同定手段に相当する）。
【００６１】
ここで、プラントモデルは、離散時間で表した離散プラントモデルであり、むだ時間を持つ２次遅れ系で近似すると共に、図３に示すように、２次遅れ系は、２つの１次遅れ系に分割して離散化した後結合して近似するようにしている。
【００６２】
プラントモデル（上流側触媒２５＋下流側排気ガスセンサ３２のモデル）の離散化は、オフラインで以下のように行う。本実施形態では、プラントモデルをむだ時間を持つ２次遅れ系で近似する。
【００６３】
連続プラントモデルは、出力をＹ（ｓ）、入力をＵ（ｓ）として次式で表される。
【００６４】
【数１】
2010019107.tif 000003
ここで、ω、ζ、Ｋ、Ｌは連続モデルパラメータで、ωは固有角振動数、ζは減衰係数、Ｋは定常ゲイン、Ｌはむだ時間である。
【００６５】
また、むだ時間Ｌは、次式で表される。
【００６６】
Ｌ＝ｄ×ｄｔ＋Ｌ１
ここで、ｄはむだサンプリング回数、ｄｔはサンプリング周期、Ｌ１は余むだ時間である。つまり、むだ時間Ｌを離散化のサンプリング周期ｄｔで除算した商ｄをむだサンプリング回数とし、余りＬ１を余むだ時間とするものである。むだ時間Ｌをサンプリング周期ｄｔで割り切れない場合は、余むだ時間Ｌ１は、後述するむだ時間更新処理により０〜ｄｔの正の値となる。
【００６７】
上記［数１］式で表される連続プラントモデルの伝達関数Ｇ（ｓ）は次式で表される。
【００６８】
【数２】
2010019107.tif 000004
この連続プラントモデルを離散化（拡張ｚ変換）する際に、部分分数展開を利用して、次のように１次遅れ系の和に分解する。ここで、Z［＊］は＊の拡張ｚ変換を表す。拡張ｚ変換は、例えば「ディジタル制御システム−解析と設計（日刊工業新聞社）」において公知である。
【００６９】
【数３】
2010019107.tif 000005
従って、離散プラントモデルは、次式で表される。
【００７０】
【数４】
2010019107.tif 000006
ここで、ｐ1，ｐ2，q1，q2，q3は離散モデルパラメータであり，q1，q2，q3の中に余むだ時間Ｌ１が含まれる。
【００７１】
以上説明した離散化式をまとめると、図４のように表される。
【００７２】
この離散モデルパラメータを連続化すると、図５に示すように、連続モデルパラメータω、ζ、Ｋ、Ｌ１が導き出される。図５中のＺｒｅ，Ｚｉｍは図６に示す算出式により求められる。
【００７３】
上記各式に関して、実際はオンボードで演算可能にするため近似演算を用いてもよい。例えば三角関数、指数関数、対数関数は、２次までテイラー展開し、誤差大の範囲はテーブルでもつものとする。
【００７４】
さて、図４に示すように、拡張ｚ変換により離散化すると、離散モデルパラメータq1，q2，q3に余むだ時間Ｌ１の情報（これが「むだ時間の一部の情報」に相当）が含まれる。このことが拡張ｚ変換の特徴であり、一般的なｚ変換による離散化では現れないものである。同定システムにより同定するのは、離散モデルパラメータθ（θ＝［ｐ1，ｐ2，q1，q2，q3］Ｔ:上付きＴは転置を表す）であり、同定値θ＿hat（ｘ＿hatはｘの同定値もしくは推定値を表すものとする、以下同様）を用いて、図５の連続化式により連続モデルパラメータの推定値ω＿hat，ζ＿hat，Ｋ＿hat，Ｌ１＿hatを算出できるようになっている。この機能が特許請求の範囲でいう余むだ時間推定手段に相当する。そして、この余むだ時間の推定値Ｌ１＿hatを用いることで、むだ時間Ｌ（むだサンプリング回数ｄ）を推定することが可能となる。この機能が特許請求の範囲でいうむだ時間推定手段に相当する。
【００７５】
以下、離散モデルパラメータを同定する同定処理について説明する。本同定処理では、同定誤差ｅ（制御対象の実出力と離散プラントモデル出力との偏差）をゼロにするように離散モデルパラメータを逐次最小二乗法により推定する。以下、図２を基に詳細に説明する。
【００７６】
離散プラントモデルは上記［数４］式で表される数式モデルであり、離散モデルパラメータθ（θ＝［ｐ1，ｐ2，q1，q2，q3］Ｔ）は適応機構により逐次推定される。
【００７７】
図２に示すように、適応機構Ｍ３２は、離散プラントモデルの予測出力と制御対象の出力（下流側排気ガスセンサ３２の検出電圧）との同定誤差ｅをゼロに近づけるように離散モデルパラメータθを推定するものである。同定誤差ｅは、離散プラントモデルに実際の制御入力ｕを入力した時のモデル出力ｙ＿hatと制御対象の実出力ｙとの誤差である。尚、離散プラントモデルに入力される制御入力ｕや制御対象の実出力ｙ（下流側排気ガスセンサ３２の検出電圧）には、不要な直流成分やノイズ成分などが含まれる。それ故に、制御入力ｕは、低周波成分除去手段としてのＨＰＦ（ハイパスフィルタ）により直流成分が除去される。
【００７８】
また、制御対象の実出力ｙにおいて、Ｆ／Ｂ制御用の出力ｙは、高周波成分除去手段としてのＬＰＦ（ローパスフィルタ）によりノイズ成分が除去され、パラメータ同定用の出力ｙは、ＨＰＦとＬＰＦとを組み合わせたＢＰＦ（バンドパスフィルタ）により直流成分とノイズ成分とが除去される。このとき、空燃比Ｆ／Ｂ用の出力ｙと同定用の出力ｙとには各々別のＬＰＦを設定できるようになっている。
【００７９】
また、同定誤差ｅには必要に応じてフィルタ処理（ＬＰＦ）や不感帯処理が施される。このフィルタ処理により、同定パラメータの振動が抑制され、不感帯処理により過同定が抑制される。尚、このフィルタ処理や不感帯処理は無くても良い。
【００８０】
ここで、前記［数４］式を、入力にむだ時間を有する自己回帰モデルとして表現しなおすと、次のように表される。
【００８１】
【数５】
2010019107.tif 000007
上式を線形パラメトリックな自己回帰モデルとして表現しなおすと、次のように表される。
【００８２】
【数６】
2010019107.tif 000008
ここで、θ＿hat(k-1)＝［ｐ1＿hat(k-1)，ｐ2＿hat(k-1)，ｑ１＿hat(k-1)，ｑ２＿hat(k-1)，ｑ３＿hat(k-1)］Ｔ，ζ(k-1)＝［ｙ(k-1)，ｙ(k-2)，ｕ(k-d)，ｕ(k-d-1)，ｕ(k-d-2)］Ｔである。同定誤差ｅは次式で求められる。
【００８３】
【数７】
2010019107.tif 000009
ここで，同定誤差eに対してノイズ除去や不感帯処理を施すことにより，同定パラメータの振動や過同定を避ける構成としてもよい．
次に、適応機構Ｍ３２において、同定誤差ｅをゼロにするように離散モデルパラメータの推定値θ＿hatを算出するパラメータ調整則は、本実施形態では重み付き最小二乗法の原理に基づいて導出される。次式に示す同定誤差ｅの２乗和を評価関数として考える。
【００８４】
【数８】
2010019107.tif 000010
ここで、λは忘却係数とも呼ばれる重み係数である。上記評価関数が最小となるようなパラメータ調整則は次式のように与えられる。
【００８５】
【数９】
2010019107.tif 000011
上式において、同定誤差ｅに乗算する行列ゲインΓ及びスカラゲインγは次式のようになる。
【００８６】
【数１０】
2010019107.tif 000012
そして、上式で求められた離散モデルパラメータの推定値が離散プラントモデルに反映される。
【００８７】
尚、パラメータ調整則は上記以外に、例えば、固定ゲイン則、漸減ゲイン則、固定忘却ゲイン則，可変ゲイン則、固定トレースゲイン則，上限トレースゲイン則，下限トレースゲイン則，両限トレースゲイン則等を用いても良い。
【００８８】
適応機構Ｍ３２で同定された離散モデルパラメータは、学習判定処理により同定誤差ｅが所定範囲内であるか否かが判定され、同定誤差ｅが所定範囲内であれば、むだ時間Ｌの学習が許可される。
【００８９】
同定された離散モデルパラメータｐ1，ｐ2，ｑ1，ｑ2，ｑ3は、図５に示す連続化式により連続化され、連続モデルパラメータω、ζ、Ｋ、Ｌ１に変換される。そして、後述する方法で、余むだ時間Ｌ１からむだ時間Ｌが学習される。この学習値は、ＥＣＵ２９のバックアップＲＡＭ等の書き換え可能な不揮発性メモリに記憶され、電源オフ中でも、学習値の記憶を保持できるようになっている。
【００９０】
サブＦ／ＢコントローラＭ２１に対しては、連続モデルパラメータω、ζ、Ｋ、Ｌ１の学習値が離散化されて入力される。サブＦ／ＢコントローラＭ２１は、学習値から算出された離散モデルパラメータを用いて制御対象の入力（入力酸素量）を算出する。この際、学習値（同定値）は逐次反映しても良いし、学習値の収束後に反映しても良い。また、学習値の範囲によって、予め算出した離散モデルパラメータ又は制御ゲインそのものを切り替えるようにしても良い。
【００９１】
次に、むだ時間更新処理について説明する。むだ時間Ｌを持つ連続系は、次式で表される。ただし，説明の簡略化のため１次遅れ系を例にとり説明するが，高次の系においても考え方は同様である．
【００９２】
【数１１】
2010019107.tif 000013
そして、上式を拡張ｚ変換により離散化することで、次式が得られる。
【００９３】
【数１２】
2010019107.tif 000014
但し、上記［数１１］式、［数１２］式において、ｙ（ｋ）＝ｘ（ｋ）である。
【００９４】
ここで、図７は、ｋｄｔ≦ｔ＜（ｋ＋１）ｄｔにおける出力ｙ（ｋ＋１）に影響を与える制御入力ｕ（むだ時間分だけ正の方向にずらした制御入力ｕ）の変化を表すタイムチャートである。図７において、むだ時間Ｌがサンプリング周期ｄｔで割り切れない場合、入力ｕ（τ−Ｌ）は、τ＝ｋｄｔ〜（ｋ＋１）ｄｔの期間内で一度だけ値が変わり、当該期間では、ｄ回前の入力ｕの影響とｄ＋１回前の入力ｕの影響を受ける。なお、上記［数１２］式では、右辺の第２項によりｄ回前の入力ｕの影響を表し、同第３項によりｄ＋１回前の入力ｕの影響を表している。
【００９５】
この場合、制御対象の実むだ時間は未知であるが、拡張ｚ変換により離散化した離散モデルパラメータには余むだ時間Ｌ１の情報が含まれており、この余むだ時間Ｌ１を用いることによりむだ時間Ｌの推定が可能となる。本実施形態では、離散モデルパラメータの推定値θ＿hatの連続化により算出される余むだ時間の推定値Ｌ１＿hatの変化に基づき、離散プラントモデルのむだ時間を制御対象の実むだ時間に近づけるように、むだ時間（むだサンプリング回数ｄ）を更新し、更新したむだ時間（むだサンプリング回数ｄ）を離散プラントモデルでのプラントモデル出力の算出に反映させるようにしている。
【００９６】
より具体的には、余むだ時間の推定値Ｌ１＿hatが、サンプリング周期ｄｔで規定される上限値近傍又は下限値近傍にあるか否かを判定する。ここで、推定値Ｌ１＿hatの上限値はｄｔ、下限値は０であり、微小な正の定数ε2により上限値近傍を「ｄｔ−ε2〜ｄｔの範囲」、下限値近傍を「０〜０＋ε2の範囲」として設定する。この場合、余むだ時間の推定値Ｌ１＿hatが上限値近傍（ｄｔ−ε2〜ｄｔ）に変化してその状態が所定時間継続した場合に、むだサンプリング回数ｄを１つ繰り上げると共に、該余むだ時間の推定値Ｌ１＿hatを下限値近傍の所定値（但し下限値＝０以上）とする。また、余むだ時間の推定値Ｌ１＿hatが下限値近傍（０〜０＋ε2）に変化してその状態が所定時間継続した場合に、むだサンプリング回数ｄを１つ繰り下げると共に、該余むだ時間の推定値Ｌ１＿hatを上限値近傍の所定値（但し上限値＝ｄｔ以下）とする。
【００９７】
その概要を図８のタイムチャートを用いて説明する。図８（ａ）に示すように、ある入力に対するノミナルむだ時間（むだ時間のノミナル値）と実むだ時間との間に図示のような誤差が生じている場合を想定する。ノミナルむだ時間は、予め設定したノミナルモデルスケジューラを用いて算出される。サンプリングタイミング（ｄ・ｄｔ）を超えた部分が余むだ時間Ｌ１である。
【００９８】
この場合、同定が逐次行われることにより、図８（ｂ）に示すように、余むだ時間の推定値Ｌ１＿hatが変化し、ｄｔ近傍（上述したｄｔ−ε2〜ｄｔの範囲）まで増加する。そして、余むだ時間の推定値Ｌ１＿hatがｄｔ近傍で所定時間以上とどまっている場合には、図８（ｃ）に示すように、むだサンプリング回数の推定値ｄ＿hatが１つ繰り上げられ（ｄ＿hat←ｄ＋１）、更に、余むだ時間の推定値Ｌ１＿hatがゼロ近傍の値に更新される。
【００９９】
その後、同定が逐次行われることにより、図８（ｄ）に示すように、余むだ時間の推定値Ｌ１＿hatが増加し、それに伴いむだ時間の推定値が実むだ時間に収束する。また、実むだ時間に収束後、劣化等により実むだ時間の変動がある場合にも、上記更新を実行することにより常にむだ時間の推定値を実むだ時間に収束させることができる。
【０１００】
なお、図８では、余むだ時間の推定値Ｌ１＿hatが増加する場合を例示し、それに伴い、むだサンプリング回数の推定値ｄ＿hatが１つ繰り上げられる様子を説明したが、これとは逆に、余むだ時間の推定値Ｌ１＿hatが減少する場合には、それに伴いむだサンプリング回数の推定値ｄ＿hatが１つ繰り下げられると共に、余むだ時間の推定値Ｌ１＿hatがｄｔ近傍の値に更新される。
【０１０１】
以上説明した適応制御処理等は、ＥＣＵ２９によって図９乃至図１６の各プログラムに従って実行される。以下、図９乃至図１６の各プログラムの処理内容を説明する。
［燃料噴射制御プログラム］
図９の燃料噴射制御プログラムは、エンジン運転中に所定周期（例えば３０℃Ａ周期）で実行される。本プログラムが起動されると、まずステップ１０１で、例えば基本噴射量マップ等を用いて、その時点のエンジン回転速度や負荷等の運転状態パラメータに基づいて基本噴射量ＴＰを算出する。
【０１０２】
この後、ステップ１０２に進み、空燃比Ｆ／Ｂ制御実行条件が成立しているか否かを判定する。ここで、空燃比Ｆ／Ｂ制御実行条件は、例えば、次の条件(1)〜(4)によって判定される。
【０１０３】
(1)上流側排気ガスセンサ３１と下流側排気ガスセンサ３２が活性化していること
(2)空燃比制御システムが正常に動作すること（フェイル判定が行われていないこと）
(3)エンジン冷却水温が所定温度以上（例えば７０℃以上）であること、つまりエンジン１１の暖機後であること
(4)エンジン運転状態がＦ／Ｂ実行領域であること
これらの条件(1)〜(4)の中で、１つでも満たさない条件があれば、空燃比Ｆ／Ｂ制御実行条件が不成立となり、ステップ１０３に進み、空燃比補正係数ＦＡＦを「１．０」に設定する。この場合は、空燃比Ｆ／Ｂ制御は行われないこととなる。
【０１０４】
これに対して、上記条件(1)〜(4)を全て満たせば、空燃比Ｆ／Ｂ制御実行条件が成立して、ステップ１０４に進み、後述する図１０の目標空燃比算出プログラムを実行して、上流側触媒２５の上流側の目標空燃比λＴＧを算出した後、ステップ１０５に進み、上流側排気ガスセンサ３１で検出した空燃比と目標空燃比λＴＧとの偏差を小さくするように空燃比補正係数ＦＡＦを算出する。
【０１０５】
以上のようにして、ステップ１０３又は１０５で、空燃比補正係数ＦＡＦを算出した後、ステップ１０６に進み、空燃比補正係数ＦＡＦやその他各種の補正係数（例えば冷却水温補正係数、学習補正係数、加減速時の補正係数等）により基本噴射量ＴＰを補正して、要求燃料噴射量ＴＡＵを算出する。
【０１０６】
ＴＡＵ＝ＴＰ×ＦＡＦ×各種補正係数
［目標空燃比算出プログラム］
図１０の目標空燃比算出プログラムは、上記図９の燃料噴射制御プログラムのステップ１０４で実行されるサブルーチンである。本プログラムが起動されると、まずステップ２０１で、後述する図１１の同定実行条件判定プログラムを実行して、同定実行条件の成否を判定し、その判定結果に応じて同定実行フラグをセット／リセット（ＯＮ／ＯＦＦ）する。
【０１０７】
その後、ステップ２０２に進み、図１１の同定実行条件判定プログラムの実行結果に基づいて同定実行フラグがＯＮであるか否かを判定し、同定実行フラグがＯＮであれば、同定実行条件が成立していると判断して、ステップ２０３に進み、下流側排気ガスセンサ３２の目標電圧を設定する。この際、元々の目標電圧に例えば所定振幅を加えた電圧を目標電圧として設定してもよい。
【０１０８】
この後、ステップ２０４に進み、後述する図１３の同定処理プログラムを実行して、制御対象の入力（上流側触媒２５の流入過不足酸素量）を離散プラントモデルに入力した時のプラントモデル出力と制御対象の出力（下流側排気ガスセンサ３２の検出電圧）との誤差ｅをゼロに近づけるように離散プラントモデルの可変パラメータを同定する。この後、ステップ２０５の処理に進む。
【０１０９】
一方、同定実行フラグがＯＦＦであれば、同定実行条件が不成立であると判断して、上記同定処理（ステップ２０３、２０４）を行わずに、ステップ２０５に進む。
【０１１０】
このステップ２０５では、サブＦ／Ｂ制御入力算出処理を実行する。このサブＦ／Ｂ制御入力算出処理では、同定処理により学習したパラメータに基づいて、下流側排気ガスセンサ３２の検出電圧を目標電圧に近付けるように、モデルベース制御則（例えば現代制御やロバスト制御）によりサブＦ／Ｂ制御入力となる入力酸素量（上流側触媒２５の流入過不足酸素量）を算出する。この際、ＰＩ制御又はＰＩＤ制御のような古典制御を用いてサブＦ／Ｂ制御入力を算出するようにしても良い。
【０１１１】
本実施形態では、入力酸素量＝目標酸素量と見なすようにしているが、入力酸素量＝検出酸素量としても良いし、空燃比補正係数から入力酸素量を算出するようにしても良い。
【０１１２】
この後、ステップ２０６に進み、目標酸素量（入力酸素量）より目標空燃比λＴＧを次のようにして算出する。まず、今回演算時の目標燃料過剰率φTGを次式により算出する。
【０１１３】
φTG＝Ｏ(k)／（Ｋ・Ｗ・ｄｔ）
Ｏ(k)：今回の目標酸素量
Ｋ：標準空気酸素質量割合
Ｗ：上流側触媒２５の流入排気ガス流量［ｇ／ｓ］
ｄｔ：演算周期
ここで、Ｗは、上流側触媒２５の流入排気ガス流量の代わりに、エンジン吸入空気量を用いても良い。また、演算周期ｄｔは、時間同期でも良いし、エンジン回転速度同期でも良い。エンジン回転速度同期の場合は、演算周期ｄｔはエンジン回転速度Ｎeの関数となる。また、ここでφTGは正確には演算周期当たりの目標燃料過剰量を表すが、単位時間当たりの目標燃料過剰量を用いてもよい。その場合、φTGは次式となる。
【０１１４】
φTG＝Ｏ(k)／（Ｋ・Ｗ）
この目標燃料過剰率φＴＧの逆数が目標空燃比λＴＧとなる。
【０１１５】
λＴＧ＝１／φＴＧ
この後、ステップ２０７に進み、目標空燃比λＴＧにベース目標空燃比を加算（又は乗算）して最終的な目標空燃比λＴＧを求める。この際、ベース目標空燃比は、例えばエンジン運転条件をパラメータとするマップ等により算出される。
［同定実行条件判定プログラム］
図１１の同定実行条件判定プログラムは、図１０の目標空燃比算出プログラムのステップ２０１で実行されるサブルーチンである。本プログラムが起動されると、まずステップ３０１で、触媒劣化検出時であるか否かを判定し、触媒劣化検出時であれば、ステップ３０５に進み、同定実行フラグをＯＮにセットする。
【０１１６】
一方、上記ステップ３０１で、触媒劣化検出時でないと判定されれば、ステップ３０２に進み、後述する図１２の運転状態判定プログラムを実行して、現在のエンジン運転状態が定常か過渡かを判定する。この後、ステップ３０３に進み、図１２の運転状態判定プログラムの実行結果に基づいて現在のエンジン運転状態が過渡であるか否かを判定し、過渡であれば、ステップ３０４に進み、同定実行フラグをＯＦＦにリセットし、定常であれば、ステップ３０５に進み、同定実行フラグをＯＮにセットする。
【０１１７】
尚、同定実行条件は、触媒劣化検出時であることのみとしたり、逆に触媒劣化検出時であることを同定実行条件から外したりすることも可能である。また、空燃比Ｆ／Ｂ制御が実行される全運転領域で同定実行フラグをＯＮにセットするようにしても良い。
［運転状態判定プログラム］
図１２の運転状態判定プログラムは、図１１の同定実行条件判定プログラムのステップ３０２で実行されるサブルーチンであり、運転状態判定手段としての役割を果たす。本プログラムが起動されると、まずステップ３０１で、エンジン回転速度変化量（変化率の絶対値）が所定値以下であるか否かを判定し、エンジン回転速度変化量が所定値以下でなければ、ステップ４０４に進み、現在のエンジン運転状態が過渡であると判定する。
【０１１８】
一方、上記ステップ４０１で、エンジン回転速度変化量が所定値以下と判定されれば、ステップ４０２に進み、負荷変化量（負荷変化率の絶対値）が所定値以下であるか否かを判定し、負荷変化量が所定値以下でなければ、ステップ４０４に進み、現在のエンジン運転状態が過渡であると判定する。
【０１１９】
上記ステップ４０１と４０２でいずれも「Ｙｅｓ」と判定された場合、つまりエンジン回転速度変化量が所定値以下で且つ負荷変化量が所定値以下である場合には、ステップ４０３に進み、現在のエンジン運転状態が定常であると判定する。
【０１２０】
尚、エンジン運転中に、所定時間以上の期間、定常と判定されない場合は、同定の実行頻度を確保するために、エンジン回転速度変化量と負荷変化量に対する定常判定値（ステップ４０１と４０２の所定値）を大きい値に変更するようにしても良い。
［同定処理プログラム］
図１３の同定処理プログラムは、図１０の目標空燃比算出プログラムのステップ２０４で実行されるサブルーチンであり、特許請求の範囲でいう同定手段としての役割を果たす。本プログラムが起動されると、まずステップ５０１で、入力酸素量（上流側触媒２５の流入過不足酸素量）を次の(1)、(2)のいずれかの方法で算出する。
【０１２１】
(1)目標酸素量をベース目標空燃比変化分に応じて補正して入力酸素量を算出する。
【０１２２】
(2)上流側触媒２５の流入ガス空燃比を、上流側排気ガスセンサ３１の検出空燃比又は燃料噴射補正係数に基づく値とし、次式により入力酸素量を算出する。
【０１２３】
入力酸素量＝Ｋ・Ｗ・（理論当量比−触媒流入ガス当量比）・ｄｔ
Ｋ：標準空気酸素質量割合
Ｗ：触媒流入排気ガス流量［ｇ／ｓ］（エンジン吸入空気量でもよい）
ｄｔ：演算周期
ここで、Ｗは、上流側触媒２５の流入排気ガス流量の代わりに、エンジン吸入空気量を用いても良い。また、演算周期ｄｔは、時間同期でも良いし、エンジン回転速度同期でも良い。エンジン回転速度同期の場合は、演算周期ｄｔはエンジン回転速度Ｎeの関数となる。また、ここで入力酸素量は正確には演算周期当たりの入力酸素量を表すが、単位時間当たりの入力酸素量を用いてもよい。その場合、入力酸素量は次式により算出される。
【０１２４】
入力酸素量＝Ｋ・Ｗ・（理論当量比−触媒流入ガス当量比）
入力酸素量の算出後、ステップ５０２に進み、離散プラントモデルの出力を算出する。この際、離散プラントモデルの出力と検出出力に同一の所定信号処理を施すものとする。この所定信号処理は、例えば、ＬＰＦ、ＨＰＦ、ＢＰＦのいずれかであるが、他の処理であっても良い。
【０１２５】
この後、ステップ５０３に進み、制御対象の実出力と離散プラントモデルの出力とから同定誤差ｅを算出する。そして、次のステップ５０４で、前述したパラメータ調整則に従ってパラメータ適応処理を実行する。これにより、制御対象の実出力と離散プラントモデルの出力との誤差をゼロに近づけるように離散モデルパラメータの推定値θ＿hatを算出する。
【０１２６】
この後、ステップ５０５に進み、学習判定処理により同定誤差ｅが所定範囲内であるか否かを判定し、同定誤差ｅが所定範囲内であれば、むだ時間Ｌの学習が許可される。この後、ステップ５０６に進み、前記算出した離散モデルパラメータの推定値θ＿hatに対して図５の連続化式を用いて連続化の処理を実行し、連続モデルパラメータ（余むだ時間の推定値Ｌ１＿hat）を算出する。尚、上記ステップ５０５の学習判定処理は、ステップ５０６の連続化処理後に行うようにしても良い。
【０１２７】
そして、次のステップ５０７（更新手段）で、後述する図１４及び図１５のむだ時間更新処理プログラムを実行して、余むだ時間Ｌ１を学習する。この後、ステップ５０８に進み、図４の離散化式を用いて連続モデルパラメータを離散化処理して離散プラントモデルのパラメータを更新する。
【０１２８】
尚、パラメータ適応処理のアルゴリズムは、本実施形態に限定されず、他のアルゴリズムを用いても良い。
［むだ時間更新処理１〜３のプログラム］
図１４，図１５，図１６のむだ時間更新処理１〜３プログラムは、図１３の同定処理プログラムのステップ５０７で、当該プログラムの演算周期ｄｔ（例えば１８０℃Ａ毎又は所定時間毎）で繰り返し実行されるサブルーチンであり、むだ時間推定手段及び学習手段としての役割を果たす。これらのプログラムは同時に並行して実行される。
【０１２９】
なお、図１７は、余むだ時間の推定値Ｌ１＿hatの時間経過に伴う変化を示す図であり、実むだ時間が短時間で急激に上昇したことに伴い推定値Ｌ１＿hatが急上昇する場合の一例を示す。ちなみに、同定の開始直後や、プラントモデルのモデル対象（実プラント）の特性が急変した場合に、余むだ時間の推定値Ｌ１＿hatは図１７の如く急上昇したり急降下したりする。
【０１３０】
［むだ時間更新処理１］
図１４のプログラムが起動されると、まずステップ６０１で、前記図１３のステップ５０６で算出した余むだ時間の推定値Ｌ１＿hatが演算周期ｄｔ以上か否かを判別する。そして、Ｌ１＿hat≧ｄｔであればステップ６０２に進む。つまり、Ｌ１＿hat≧ｄｔとなる程度に実むだ時間が増加する向きに急変した場合に、以降のステップ６０２〜６０６の処理を実行することとなる。
【０１３１】
ステップ６０２では、時間経過とともに変動する推定値Ｌ１＿hatのうちｄｔ〜２ｄｔの範囲についての積算値Ｓ１（図１７（ａ）中の斜線部分面積Ｓ１に相当）を演算する。つまり、推定値Ｌ１＿hatがｄｔ〜２ｄｔの範囲で変動する期間においては、推定値Ｌ１＿hatからｄｔを減算した値を積算し、推定値Ｌ１＿hatが２ｄｔよりも大きい期間においては２ｄｔからｄｔを減算した一定値（つまりｄｔ）を積算する。このようにして算出された値を、現在のＳ１の値に加算してＳ１を更新する。
【０１３２】
その後、ステップ６０３に進み、推定値Ｌ１＿hat＜ｄｔの状態が所定時間継続したか否かを判定する。このステップ６０３で「Ｙｅｓ」と判定されれば、ステップ６０７に進み、ステップ６０２で算出した積算値Ｓ１をゼロにリセットする。
【０１３３】
また、上記ステップ６０３で「Ｎｏ」と判定されればステップ６０４（第１増加判定手段）に進み、ステップ６０２で算出した積算値Ｓ１が予め設定された閾値ｔｈ１０以上か否かを判定する。このステップ６０４で「Ｙｅｓ」と判定されれば、ステップ６０５に進み、余むだ時間Ｌ１＿ａｄｐ、むだサンプリング回数変化量Δｄ＿ａｄｐをそれぞれ更新する。
【０１３４】
Ｌ１＿ａｄｐ＝０
Δｄ＿ａｄｐ＝＋１
その後、ステップ６０６に進み、後述する図１５のステップ７０２（むだ時間更新処理２）で算出された積算値Ｓ２を、むだ時間更新処理１で用いられる積算値Ｓ１に置き換える。
【０１３５】
一方、上記ステップ６０１で「Ｎｏ」と判定されればステップ６０８に進み、余むだ時間の推定値Ｌ１＿hatがゼロ以下か否かを判別する。そして、Ｌ１＿hat≦０であればステップ６０９に進む。つまり、Ｌ１＿hat≦０となる程度に実むだ時間が減少する向きに急変した場合に、以降のステップ６０９〜６１３の処理を実行することとなる。また、このステップ６０８にて「Ｎｏ」と判定されれば処理は６０１に戻る。つまり、実むだ時間が増加及び減少のいずれの向きにも急変せず、余むだ時間の推定値Ｌ１＿hatが０〜ｄｔの範囲で変動している状況では、ステップ６０２〜６０６及びステップ６０９〜６１３のいずれも実行されないこととなる。
【０１３６】
ステップ６０９では、時間経過とともに変動する推定値Ｌ１＿hatのうち−ｄｔ〜０の範囲についての積算値Ｒ１を演算する。つまり、推定値Ｌ１＿hatが−ｄｔ〜０の範囲で変動する期間においては、推定値Ｌ１＿hatをそのまま積算し、推定値Ｌ１＿hatが−ｄｔよりも小さい期間においては−２ｄｔから−ｄｔを減算した一定値（つまり−ｄｔ）を積算する。このようにして算出された値を現在のＲ１の値に加算して、Ｒ１を更新する。なお、この積算値Ｒ１は、正確には絶対値に変換されて正の値となるよう記憶更新される。
【０１３７】
その後、ステップ６１０に進み、推定値Ｌ１＿hat＞０の状態が所定時間継続したか否かを判定する。このステップ６１０で「Ｙｅｓ」と判定されれば、ステップ６１４に進み、ステップ６０９で算出した積算値Ｒ１をゼロにリセットする。
【０１３８】
また、上記ステップ６１０で「Ｎｏ」と判定されればステップ６１１（第１減少判定手段）に進み、ステップ６０９で算出した積算値Ｒ１が予め設定された閾値ｔｈ１１以上か否かを判定する。このステップ６１１で「Ｙｅｓ」と判定されれば、ステップ６１２に進み、余むだ時間Ｌ１＿ａｄｐ、むだサンプリング回数変化量Δｄ＿ａｄｐをそれぞれ更新する。
【０１３９】
Ｌ１＿ａｄｐ＝０
Δｄ＿ａｄｐ＝−１
その後、ステップ６１３に進み、後述する図１５のステップ７０９（むだ時間更新処理２）で算出された積算値Ｒ２を、図１４のむだ時間更新処理１で用いられる積算値Ｒ１に置き換える。
【０１４０】
以上のむだ時間更新処理１によれば、Ｓ１≧ｔｈ１０となる程度に推定値Ｌ１＿hatが上昇すると、図８（ｃ）に示すように、むだサンプリング回数の推定値ｄ＿hatが１つ繰り上げられ（ｄ＿hat←ｄ＋１）、余むだ時間の推定値Ｌ１＿hatがゼロに更新されることとなる。一方、Ｒ１≧ｔｈ１１となる程度に推定値Ｌ１＿hatが下降すると、むだサンプリング回数の推定値ｄ＿hatが１つ繰り下げられ（ｄ＿hat←ｄ−１）、余むだ時間の推定値Ｌ１＿hatがゼロに更新されることとなる。
【０１４１】
［むだ時間更新処理２］
上記図１４の更新処理１では、むだサンプリング回数推定値ｄ＿hatの繰り上げ、繰り下げの判定を、ｄｔ〜２ｄｔの範囲についての積算値Ｓ１、及び−ｄｔ〜０の範囲についての積算値Ｒ１に基づき実施している。これに対し図１５の更新処理２では、２ｄｔ〜３ｄｔの範囲についての積算値Ｓ２、及び−２ｄｔ〜−ｄｔの範囲についての積算値Ｒ２に基づき実施している。以下、主にむだ時間更新処理１との相違点についてむだ時間更新処理２の手順を説明する。
【０１４２】
まずステップ７０１で、推定値Ｌ１＿hatが演算周期２ｄｔ以上か否かを判別する。そして、Ｌ１＿hat≧２ｄｔであればステップ７０２に進む。つまり、Ｌ１＿hat≧２ｄｔとなる程度に実むだ時間が増加する向きに急変した場合に、以降のステップ７０２〜７０６の処理を実行することとなる。
【０１４３】
ステップ７０２では、変動する推定値Ｌ１＿hatのうち２ｄｔ〜３ｄｔの範囲についての積算値Ｓ２（図１７（ａ）中の斜線部分面積Ｓ２に相当）を演算する。つまり、推定値Ｌ１＿hatが２ｄｔ〜３ｄｔの範囲で変動する期間においては、推定値Ｌ１＿hatから２ｄｔを減算した値を積算し、推定値Ｌ１＿hatが３ｄｔよりも大きい期間においては３ｄｔから２ｄｔを減算した一定値（つまりｄｔ）を積算する。このようにして算出された値を、現在のＳ２の値に加算してＳ２を更新する。
【０１４４】
その後、続くステップ７０３において、推定値Ｌ１＿hat＜２ｄｔの状態が所定時間継続したと判定（７０３：Ｙｅｓ）されれば、ステップ７０７に進み積算値Ｓ２をゼロにリセットする。また、上記ステップ７０３で「Ｎｏ」と判定されればステップ７０４（第２増加判定手段）に進み、ステップ７０２で算出した積算値Ｓ２が予め設定された閾値ｔｈ２０以上であると判定（７０４：Ｙｅｓ）されれば、ステップ７０５において、余むだ時間Ｌ１＿ａｄｐ、むだサンプリング回数変化量Δｄ＿ａｄｐをそれぞれ更新する。
【０１４５】
Ｌ１＿ａｄｐ＝０
Δｄ＿ａｄｐ＝＋１
その後、ステップ７０６に進み、後述する図１６のステップ８０２（むだ時間更新処理３）で算出された積算値Ｓ３を、図１５のむだ時間更新処理２で用いられる積算値Ｓ２に置き換える。
【０１４６】
一方、上記ステップ７０１で「Ｎｏ」と判定されればステップ７０８に進み、余むだ時間の推定値Ｌ１＿hatが−ｄｔ以下か否かを判別する。そして、Ｌ１＿hat≦−ｄｔであればステップ７０９に進む。つまり、Ｌ１＿hat≦−ｄｔとなる程度に実むだ時間が減少する向きに急変した場合に、以降のステップ７０９〜７１３の処理を実行することとなる。また、このステップ７０８にて「Ｎｏ」と判定されれば処理は７０１に戻る。つまり、余むだ時間の推定値Ｌ１＿hatが−ｄｔ〜２ｄｔの範囲で変動している状況では、ステップ７０２〜７０６及びステップ７０９〜７１３のいずれも実行されないこととなる。
【０１４７】
ステップ７０９では、変動する推定値Ｌ１＿hatのうち−２ｄｔ〜−ｄｔの範囲についての積算値Ｒ２を演算する。つまり、推定値Ｌ１＿hatが−２ｄｔ〜−ｄｔの範囲で変動する期間においては、推定値Ｌ１＿hatから−ｄｔを減算した値を積算し、推定値Ｌ１＿hatが−２ｄｔよりも小さい期間においては−３ｄｔから−２ｄｔを減算した一定値（つまり−ｄｔ）を積算する。このようにして算出された値を現在のＲ２の値に加算して、Ｒ２を更新する。なお、この積算値Ｒ２は、正確には絶対値に変換されて正の値となるよう記憶更新される。
【０１４８】
その後、続くステップ７１０において、推定値Ｌ１＿hat＞-ｄｔの状態が所定時間継続したと判定（７１０：Ｙｅｓ）されれば、ステップ７１４に進み積算値Ｒ２をゼロにリセットする。また、上記ステップ７１０で「Ｎｏ」と判定されればステップ７１１（第２減少判定手段）に進み、ステップ７０９で算出した積算値Ｒ２が予め設定された閾値ｔｈ２１以上であると判定（７１１：Ｙｅｓ）されれば、ステップ７１２において、余むだ時間Ｌ１＿ａｄｐ、むだサンプリング回数変化量Δｄ＿ａｄｐをそれぞれ更新する。
【０１４９】
Ｌ１＿ａｄｐ＝０
Δｄ＿ａｄｐ＝−１
その後、ステップ７１３に進み、後述する図１６のステップ８０９（むだ時間更新処理３）で算出された積算値Ｒ３を、図１５のむだ時間更新処理２で用いられる積算値Ｒ２に置き換える。
【０１５０】
以上のむだ時間更新処理２によれば、Ｓ２≧ｔｈ２０となる程度に推定値Ｌ１＿hatが上昇すると、図８（ｃ）に示すように、むだサンプリング回数の推定値ｄ＿hatが１つ繰り上げられ（ｄ＿hat←ｄ＋１）、余むだ時間の推定値Ｌ１＿hatがゼロに更新されることとなる。一方、Ｒ２≧ｔｈ２１となる程度に推定値Ｌ１＿hatが下降すると、むだサンプリング回数の推定値ｄ＿hatが１つ繰り下げられ（ｄ＿hat←ｄ−１）、余むだ時間の推定値Ｌ１＿hatがゼロに更新されることとなる。
【０１５１】
［むだ時間更新処理３］
図１６のむだ時間更新処理３では、むだサンプリング回数推定値ｄ＿hatの繰り上げ、繰り下げの判定を、３ｄｔ〜４ｄｔの範囲についての積算値Ｓ３、及び−３ｄｔ〜−２ｄｔの範囲についての積算値Ｒ３に基づき実施している。以下、主に更新処理１との相違点について更新処理３の手順を説明する。
【０１５２】
まずステップ８０１で、推定値Ｌ１＿hatが演算周期３ｄｔ以上か否かを判別する。そして、Ｌ１＿hat≧３ｄｔであればステップ８０２に進む。つまり、Ｌ１＿hat≧３ｄｔとなる程度に実むだ時間が増加する向きに急変した場合に、以降のステップ８０２〜８０６の処理を実行することとなる。
【０１５３】
ステップ８０２では、変動する推定値Ｌ１＿hatのうち３ｄｔ〜４ｄｔの範囲についての積算値Ｓ３（図１７（ａ）中の斜線部分面積Ｓ３に相当）を演算する。つまり、推定値Ｌ１＿hatが３ｄｔ〜４ｄｔの範囲で変動する期間においては、推定値Ｌ１＿hatから３ｄｔを減算した値を積算し、推定値Ｌ１＿hatが４ｄｔよりも大きい期間においては４ｄｔから３ｄｔを減算した一定値（つまりｄｔ）を積算する。このようにして算出された値を、現在のＳ３の値に加算してＳ３を更新する。
【０１５４】
その後、続くステップ８０３において、推定値Ｌ１＿hat＜３ｄｔの状態が所定時間継続したと判定（８０３：Ｙｅｓ）されれば、ステップ８０７に進み積算値Ｓ３をゼロにリセットする。また、上記ステップ８０３で「Ｎｏ」と判定されればステップ８０４（第３増加判定手段）に進み、ステップ８０２で算出した積算値Ｓ３が予め設定された閾値ｔｈ３０以上であると判定（８０４：Ｙｅｓ）されれば、ステップ８０５において、余むだ時間Ｌ１＿ａｄｐ、むだサンプリング回数変化量Δｄ＿ａｄｐをそれぞれ更新する。
【０１５５】
Ｌ１＿ａｄｐ＝０
Δｄ＿ａｄｐ＝＋１
その後、ステップ８０６に進み、図示しない別の更新処理にて算出された積算値Ｓ４を、図１６のむだ時間更新処理３で用いられる積算値Ｓ３に置き換える。或いは、積算値Ｓ３をゼロにリセットする。
【０１５６】
一方、上記ステップ８０１で「Ｎｏ」と判定されればステップ８０８に進み、余むだ時間の推定値Ｌ１＿hatが−２ｄｔ以下か否かを判別する。そして、Ｌ１＿hat≦−２ｄｔであればステップ８０９に進む。つまり、Ｌ１＿hat≦−２ｄｔとなる程度に実むだ時間が減少する向きに急変した場合に、以降のステップ８０９〜８１３の処理を実行することとなる。また、このステップ８０８にて「Ｎｏ」と判定されれば処理は８０１に戻る。つまり、余むだ時間の推定値Ｌ１＿hatが−２ｄｔ〜３ｄｔの範囲で変動している状況では、ステップ８０２〜８０６及びステップ８０９〜８１３のいずれも実行されないこととなる。
【０１５７】
ステップ８０９では、変動する推定値Ｌ１＿hatのうち−３ｄｔ〜−２ｄｔの範囲についての積算値Ｒ３を演算する。つまり、推定値Ｌ１＿hatが−３ｄｔ〜−２ｄｔの範囲で変動する期間においては、推定値Ｌ１＿hatから−２ｄｔを減算した値を積算し、推定値Ｌ１＿hatが−３ｄｔよりも小さい期間においては−４ｄｔから−３ｄｔを減算した一定値（つまり−ｄｔ）を積算する。このようにして算出された値を現在のＲ３の値に加算して、Ｒ３を更新する。なお、この積算値Ｒ３は、正確には絶対値に変換されて正の値となるよう記憶更新される。
【０１５８】
その後、続くステップ８１０において、推定値Ｌ１＿hat＞-２ｄｔの状態が所定時間継続したと判定（８１０：Ｙｅｓ）されれば、ステップ８１４に進み積算値Ｒ３をゼロにリセットする。また、上記ステップ８１０で「Ｎｏ」と判定されればステップ８１１（第３減少判定手段）に進み、ステップ８０９で算出した積算値Ｒ３が予め設定された閾値ｔｈ３１以上であると判定（８１１：Ｙｅｓ）されれば、ステップ８１２において、余むだ時間Ｌ１＿ａｄｐ、むだサンプリング回数変化量Δｄ＿ａｄｐをそれぞれ更新する。
【０１５９】
Ｌ１＿ａｄｐ＝０
Δｄ＿ａｄｐ＝−１
その後、ステップ８１３に進み、図示しない別の更新処理にて算出された積算値Ｒ４を、図１６のむだ時間更新処理３で用いられる積算値Ｒ３に置き換える。或いは、積算値Ｒ３をゼロにリセットする。
【０１６０】
以上のむだ時間更新処理３によれば、Ｓ３≧ｔｈ３０となる程度に推定値Ｌ１＿hatが上昇すると、図８（ｃ）に示すように、むだサンプリング回数の推定値ｄ＿hatが１つ繰り上げられ（ｄ＿hat←ｄ＋１）、余むだ時間の推定値Ｌ１＿hatがゼロに更新されることとなる。一方、Ｒ３≧ｔｈ３１となる程度に推定値Ｌ１＿hatが下降すると、むだサンプリング回数の推定値ｄ＿hatが１つ繰り下げられ（ｄ＿hat←ｄ−１）、余むだ時間の推定値Ｌ１＿hatがゼロに更新されることとなる。
【０１６１】
［むだ時間更新処理１〜３のまとめ］
上述した更新処理１〜３は、それぞれ同時に並行して実施される。そして図１７（ａ）の例では、ｔ１時点において、更新処理１により積算値Ｓ１≧ｔｈ１０と判定され、むだサンプリング回数の推定値ｄ＿hatが１つ繰り上げられている。また、ｔ１時点後のｔ２時点において、更新処理２により積算値Ｓ２≧ｔｈ２０と判定され、むだサンプリング回数の推定値ｄ＿hatがさらに１つ繰り上げられている。なお、更新処理３による積算値Ｓ３についてはｔｈ３０を超えたとは判定されていない状態である。
【０１６２】
また、本実施形態では、むだ時間更新処理３にて繰り上げ更新の判定に用いられる閾値ｔｈ３０は、更新処理２で用いられる閾値ｔｈ２０よりも大きい値に設定され、さらに閾値ｔｈ２０は、更新処理１で用いられる閾値ｔｈ１０よりも大きい値に設定されている。
【０１６３】
ｔｈ３０＞ｔｈ２０＞ｔｈ１０
同様にして、むだ時間更新処理３にて繰り下げ更新の判定に用いられる閾値ｔｈ３１は、更新処理２で用いられる閾値ｔｈ２１よりも大きい値に設定され、さらに閾値ｔｈ２１は、更新処理１で用いられる閾値ｔｈ１１よりも大きい値に設定されている。
【０１６４】
ｔｈ３１＞ｔｈ２１＞ｔｈ１１
以上詳述した本実施形態によれば、以下の効果が得られるようになる。
【０１６５】
（１）先述した特許文献１，２に記載の更新処理では、余むだ時間推定値Ｌ１＿hatがｄｔ近傍又はゼロ近傍で所定時間Ｔ１０（図１７（ｂ））以上とどまっていることを条件として、むだサンプリング回数の推定値ｄ＿hatの繰り上げ又は繰り下げの更新が為される。そして、所定時間Ｔ１０内に余むだ時間推定値Ｌ１＿hatがｄｔを大きく超えて急激に増加した図１７（ｂ）の場合では、余むだ時間推定値Ｌ１＿hat≧ｄｔの状態が所定時間Ｔ１０経過したｔ１時点で、推定値ｄ＿hatを１つ繰り上げる。この繰り上げ更新の後、そのｔ１時点からさらに所定時間Ｔ１０経過したｔ３時点で推定値ｄ＿hatを１つ繰り上げる。要するに、所定時間Ｔ１０内に余むだ時間推定値Ｌ１＿hatがｄｔを大きく超えて急激に増加した場合において、推定値ｄ＿hatを２つ分繰り上げるのに要する時間は所定時間Ｔ１０の２倍となる。
【０１６６】
これに対し、以上説明した本実施形態によれば、複数の更新処理１〜３を同時に並行して実施するので、更新処理１による更新処理が終了することを待たずして更新処理２，３による更新処理が実施される。そのため、例えば、余むだ時間推定値Ｌ１＿hatがｄｔを大きく超えて急激に変化した場合において、推定値ｄ＿hatを２つ分繰り上げるのに要する時間は所定時間Ｔ１０の２倍よりも短くできる。
【０１６７】
（２）さらに、本実施形態によれば、ｔｈ３０＞ｔｈ２０＞ｔｈ１０となるよう繰り上げ更新の判定に用いられる閾値を設定し、ｔｈ３１＞ｔｈ２１＞ｔｈ１１となるよう繰り下げ更新の判定に用いられる閾値を設定している。そして、余むだ時間推定値Ｌ１＿hatが変動することにより更新処理１における更新の判定結果が受ける影響は、更新処理２における判定結果が受ける影響に比べて大きい。また、更新処理２における影響は更新処理３に比べて大きい。
【０１６８】
よって、余むだ時間推定値Ｌ１＿hatのノイズ等による変動から更新処理２，３におけるの更新判定結果が受ける影響を好適に小さくできる。換言すれば、更新処理１による更新判定結果を取得できる時期を早くできる。つまり、更新処理１による更新判定結果を早期に取得できるとともに、更新処理２，３がノイズ等の影響を受けて誤判定することを抑制できる。
【０１６９】
（３）本実施形態では、拡張ｚ変換により離散化を行って離散モデルパラメータに余むだ時間Ｌ１の情報を含ませると共に、該離散モデルパラメータから算出される余むだ時間Ｌ１の推定値を基に、離散プラントモデルのむだ時間を制御対象の実むだ時間に近づけるように推定し、さらに、離散化のむだサンプリング回数ｄを増加又は減少させることで、推定むだ時間を実むだ時間の変化に追従させるようにしたので、応答遅れだけでなくむだ時間も逐次推定できることとなり、制御対象のばらつきや経時変化によらず制御対象の動特性を十分把握することができ、高い制御精度を維持することができる。また、むだ時間の情報を種々に活用することもできる。
【０１７０】
（４）本実施形態では、同定誤差に基づいて逐次推定された離散モデルパラメータを連続モデルパラメータに変換し、その連続モデルパラメータから余むだ時間の推定値Ｌ１＿hatを算出するようにしたので、離散モデルパラメータから演算を行う場合に比べて、余むだ時間の推定値Ｌ１＿hatの演算を容易に行うことができる。
【０１７１】
（５）推定むだ時間を学習する構成としたため、制御対象の個体差や経時変化等により実むだ時間のずれが生じている場合にも、それに追従した制御を継続的に実施できる。このとき、プラントモデル出力と制御対象の実出力との誤差が所定範囲内である場合に学習値を更新するようにしているので、適正な同定が行われていることを条件に学習処理が行われ、学習値の信頼性が高められる利点もある。
【０１７２】
（他の実施形態）
上記実施形態は、以下のように変更して実施してもよい。また、本発明は上記実施形態の記載内容に限定されず、以下に説明する各実施形態の特徴的構成をそれぞれ任意に組み合わせるようにしてもよい。
【０１７３】
・上記実施形態では、むだサンプリング回数ｄの繰り上げ／繰り下げを判定するにあたり、積算値Ｓ１〜Ｓ３，Ｒ１〜Ｒ３に基づき判定しているが、本発明はこれに限られるものではなく、例えば、更新処理１においては、変動する余むだ時間の推定値Ｌ１＿hatが、ｄｔ以上である状態が所定時間以上経過した場合にむだサンプリング回数ｄを繰り上げるよう判定し、ゼロ以下である状態が所定時間以上経過した場合にむだサンプリング回数ｄを繰り下げるよう判定してもよい。
【０１７４】
・更新処理２についても同様にして、Ｌ１＿hat≧２ｄｔ又はＬ１＿hat≦−ｄｔが所定時間以上経過した場合に繰り上げ又は繰り下げるよう判定し、更新処理３についても同様にして、Ｌ１＿hat≧３ｄｔ又はＬ１＿hat≦−２ｄｔが所定時間以上経過した場合に繰り上げ又は繰り下げるよう判定してもよい。但し、複数の更新処理を同時に並行して実施することを要する。
【０１７５】
・上記実施形態では、３つの更新処理１〜３を同時に並行して実施しているが、本発明は少なくとも２つの更新処理を並行して実施していればよく、例えば４つ以上の更新処理を並行して実施してもよい。
【０１７６】
・上記実施形態では、むだサンプリング回数ｄの繰り上げ／繰り下げを判定するにあたり、積算値Ｓ１〜Ｓ３，Ｒ１〜Ｒ３に基づき判定しているが、本発明はこれに限られるものではなく、例えば、更新処理１において、L1_hatをｄｔで除算した後、０〜１を不感帯として不感帯処理（L1_hat＜０ならL1_hat−０を出力，０≦L1_hat≦１なら０を出力，１＜L1_hatならL1_hat−１を出力）を施した値を積算し、該積算値が正の所定閾値以上である場合にむだサンプリング回数ｄを繰り上げるよう判定し、負の所定閾値以下である場合にむだサンプリング回数ｄを繰り下げるよう判定してもよい。また、積算の前に−１以上０以下の所定下限値と０以上１以下の所定上限値で制限してから、積算するようにすると，ノイズ外乱による誤判定を軽減できる。また、積算の前に奇数乗することにより符号を保持しつつ重みを付けてから積算するようにすると、ノイズ外乱による誤判定を軽減しつつ判定精度を向上させることができる。
【０１７７】
・更新処理２についても同様にして、L1_hatをｄｔで除算した後、−１〜２を不感帯として不感帯処理（L1_hat＜−１ならL1_hat−（−１）を出力，−１≦L1_hat≦２なら０を出力，２＜L1_hatならL1_hat−２を出力）を施した値を積算し、該積算値が正の所定閾値以上である場合にむだサンプリング回数ｄを繰り上げるよう判定し、負の所定閾値以下である場合にむだサンプリング回数ｄを繰り下げるよう判定してもよい。また、積算の前に−１以上０以下の所定下限値と０以上１以下の所定上限値で制限してから、積算するようにすると，ノイズ外乱による誤判定を軽減できる。また、積算の前に奇数乗することにより符号を保持しつつ重みを付けてから積算するようにすると、ノイズ外乱による誤判定を軽減しつつ判定精度を向上させることができる。
【０１７８】
・更新処理３についても同様にして、L1_hatをｄｔで除算した後、−２〜３を不感帯として不感帯処理（L1_hat＜−２ならL1_hat−（−２）を出力，−２≦L1_hat≦３なら０を出力，３＜L1_hatならL1_hat−３を出力）を施した値を積算し、該積算値が正の所定閾値以上である場合にむだサンプリング回数ｄを繰り上げるよう判定し、負の所定閾値以下である場合にむだサンプリング回数ｄを繰り下げるよう判定してもよい。また、積算の前に−１以上０以下の所定下限値と０以上１以下の所定上限値で制限してから、積算するようにすると，ノイズ外乱による誤判定を軽減できる。また、積算の前に奇数乗することにより符号を保持しつつ重みを付けてから積算するようにすると、ノイズ外乱による誤判定を軽減しつつ判定精度を向上させることができる。
【０１７９】
・上記実施形態では、離散モデルパラメータの推定値θ＿hatを連続モデルパラメータに変換し、その連続モデルパラメータに含まれる余むだ時間の推定値Ｌ１＿hatを用いてむだ時間Ｌ（むだサンプリング回数ｄ）を推定する構成としたが、これに代えて、離散パラメータの連続化を行わず、余むだ時間Ｌ１を含む離散パラメータの推定値を用いてむだ時間Ｌ（むだサンプリング回数ｄ）を推定することも可能である。
【０１８０】
・上記実施形態では、制御対象（上流側触媒２５と下流側排気ガスセンサ３２）をむだ時間を含む２次遅れ系で近似するようにしたが、本発明は、制御対象をむだ時間を含む３次以上の高次遅れ系で近似するようにしても良い。この場合、３次以上の高次遅れ系は、３つ以上の１次遅れ系に分割して離散化した後結合して近似するようにすれば良い。
【０１８１】
・本発明の適用範囲は、上流側触媒２５と下流側排気ガスセンサ３２を含むサブＦ／Ｂ系（空燃比制御システム）に限定されず、エンジン１１と上流側排気ガスセンサ３１を含むメインＦ／Ｂ系（内燃機関制御システム）に本発明を適用しても良い。この場合、同定手段は、エンジン運転条件に応じて設定されたノミナルモデルパラメータ（モデルパラメータのノミナル値）からの偏差を同定し、エンジン運転条件を表す信号の高周波振動成分を抑制するためのフィルタ手段を設けるようにすると良い。ノミナルモデルパラメータは、予め設定したノミナルモデルスケジューラを用いて算出される。
【０１８２】
・上記実施形態では、排気管２３に２つの触媒２５，３０を設けたが、触媒を１つのみ設けたシステムにも本発明を適用して実施できる。その他、本発明は、サブＦ／Ｂ系、メインＦ／Ｂ系に限定されず、他の制御対象に適用しても良い。
【図面の簡単な説明】
【０１８３】
【図１】本発明の一実施形態におけるエンジン制御システムの概略を示す構成図である。
【図２】ＥＣＵで実現するメインＦ／Ｂ制御機能とサブＦ／Ｂ制御機能の概要を示すブロック図である。
【図３】２次遅れ＋むだ時間で近似する連続プラントモデルを、１次遅れ系に分割して離散化した後結合して離散プラントモデルを作成する手法を説明する図である。
【図４】連続モデルパラメータを離散化して離散モデルパラメータを演算する方法を説明する図である。
【図５】離散モデルパラメータ（同定値）を連続化して連続モデルパラメータを演算する方法を説明する図である。
【図６】図５中のＺｒｅ，Ｚｉｍを算出する方法を説明する図である。
【図７】制御入力ｕの変化を表すタイムチャートである。
【図８】むだ時間更新処理の概要を具体的に示すタイムチャートである。
【図９】燃料噴射制御プログラムの処理の流れを示すフローチャートである。
【図１０】目標空燃比算出プログラムの処理の流れを示すフローチャートである。
【図１１】同定実行条件判定プログラムの処理の流れを示すフローチャートである。
【図１２】運転状態判定プログラムの処理の流れを示すフローチャートである。
【図１３】同定処理プログラムの処理の流れを示すフローチャートである。
【図１４】むだ時間更新処理プログラムの処理の流れを示すフローチャートである（処理１）。
【図１５】むだ時間更新処理プログラムの処理の流れを示すフローチャートである（処理２）。
【図１６】むだ時間更新処理プログラムの処理の流れを示すフローチャートである（処理３）。
【図１７】余りむだ時間推定値Ｌ１＿hatの変化に対するサンプリング回数ｄの更新状態を示すタイムチャートであり、（ａ）は図１４〜図１６の処理による更新状態、（ｂ）は従来の処理による更新状態を示す。
【符号の説明】
【０１８４】
１１…エンジン（内燃機関）、２３…排気管、２５…上流側の触媒、２９…ＥＣＵ（同定手段，余むだ時間推定手段，むだ時間推定手段）、３０…下流側の触媒、３１…上流側排気ガスセンサ、３２…下流側排気ガスセンサ、Ｍ１１…メインＦ／Ｂコントローラ、Ｍ２１…サブＦ／Ｂコントローラ、Ｍ３２…適応機構、５０７…更新手段、６０４…第１増加判定手段、７０４…第２増加判定手段、８０４…第３増加判定手段、６１１…第１減少判定手段、７１１…第２減少判定手段、８１１…第３減少判定手段。
【図１】
2010019107.tif 000015
【図２】
2010019107.tif 000016
【図３】
2010019107.tif 000017
【図４】
2010019107.tif 000018
【図５】
2010019107.tif 000019
【図６】
2010019107.tif 000020
【図７】
2010019107.tif 000021
【図８】
2010019107.tif 000022
【図９】
2010019107.tif 000023
【図１０】
2010019107.tif 000024
【図１１】
2010019107.tif 000025
【図１２】
2010019107.tif 000026
【図１３】
2010019107.tif 000027
【図１４】
2010019107.tif 000028
【図１５】
2010019107.tif 000029
【図１６】
2010019107.tif 000030
【図１７】
2010019107.tif 000031