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