(57)【要約】【課題】  この発明は内燃機関の点火時期制御装置に関し、点火時期を調整しようとする際に、点火時期調整後の運転安定性を事前に予測することを目的とする。【解決手段】  現点火時期SAcよりも進角した想定点火時期SAiを仮想する（ステップ１１８）。想定点火時期SAiの下での筒内圧履歴データを、相当数のサイクル分、生成する（ステップ１２０）。それらの筒内圧履歴データに基づいて、想定点火時期SAiの下でのノックサイクル頻度を予測する（ステップ１３４）。予測されたノックサイクル頻度に基づいて、ノッキングを起こさない範囲でなるべくMBTに近くなるように、点火時期を設定する（ステップ１５０〜１５８）。【選択図】      図９
2007056727.tif 000002
【特許請求の範囲】
【請求項１】
内燃機関の点火時期を現在値から変更したと仮定した場合の、クランク角に対する筒内圧履歴を模擬した模擬筒内圧履歴データを生成する模擬筒内圧履歴データ生成手段と、
多数サイクル分の前記模擬筒内圧履歴データに基づいて、点火時期変更後の運転安定性を事前に予測する運転安定性予測手段と、
前記運転安定性予測手段の予測結果に基づいて、実際の点火時期を設定する点火時期設定手段と、
を備えることを特徴とする内燃機関の点火時期制御装置。
【請求項２】
前記模擬筒内圧履歴データ生成手段は、
所定の燃焼パラメータに基づいて燃焼状態を模擬する燃焼モデル手段と、
点火時期変更後の前記燃焼パラメータの分布の標準偏差を予測する標準偏差予測手段と、
前記標準偏差予測手段の予測値を標準偏差とする正規分布から、前記燃焼パラメータの値をランダムに抽出する燃焼パラメータ抽出手段と、
前記燃焼パラメータ抽出手段により抽出された値を前記燃焼モデル手段に入力することで得られる模擬の燃焼状態に基づいて、前記模擬筒内圧履歴データを算出する筒内圧履歴算出手段と、
を含むことを特徴とする請求項１記載の内燃機関の点火時期制御装置。
【請求項３】
前記燃焼パラメータは、燃焼効率に関する値と、所定の燃焼割合となるクランク角に関する値と、第１の燃焼割合から第２の燃焼割合に変化するまでの燃焼期間に関する値とを含むことを特徴とする請求項２記載の内燃機関の点火時期制御装置。
【請求項４】
前記標準偏差予測手段は、
現在の運転状態における前記燃焼パラメータの現実値の平均値を取得する平均値取得手段と、
前記平均値取得手段により取得された平均値と現在の運転状態を表す所定の運転パラメータとに基づいて前記燃焼パラメータの標準偏差予測値を導出するための演算式を記憶した手段と、
現在の運転状態における前記運転パラメータの現実値を取得する運転パラメータ取得手段と、
を含むことを特徴とする請求項２または３記載の内燃機関の点火時期制御装置。
【請求項５】
前記平均値取得手段は、
実筒内圧を検出する筒内圧センサと、
クランク角に対する前記実筒内圧の履歴に基づいて、現在の運転状態における前記燃焼パラメータの現実値を各サイクルごとに算出する燃焼パラメータ算出手段と、
前記燃焼パラメータ算出手段の算出値の種別ごとの算術平均を算出する算術平均手段と、
を含むことを特徴とする請求項４記載の内燃機関の点火時期制御装置。
【請求項６】
前記燃焼パラメータ算出手段により算出された前記燃焼パラメータの現実値の分布から、前記燃焼パラメータの実標準偏差を算出する実標準偏差算出手段と、
前記標準偏差予測値を導出するための演算式に、前記実標準偏差と、前記運転パラメータの現実値とを入力することにより、前記演算式内の係数を較正する較正手段を更に備えることを特徴とする請求項５記載の内燃機関の点火時期制御装置。
【請求項７】
前記運転安定性予測手段は、
各々の前記模擬筒内圧履歴データごとに、そのサイクルがノックを起こすノックサイクルであるか否かを判定するノックサイクル判定手段と、
前記ノックサイクルの発生頻度が許容レベル以下であるか否かを判別するノッキング判定手段と、
を含むことを特徴とする請求項１乃至６の何れか１項記載の内燃機関の点火時期制御装置。
【請求項８】
前記ノックサイクル判定手段は、
前記模擬筒内圧履歴データに基づいて、そのサイクルの自己着火時点を予測する自己着火予測手段と、
前記予測された自己着火時点での燃焼割合が判定値以下である場合に、そのサイクルがノックサイクルであると判定する燃焼割合算出判定手段と、
を含むことを特徴とする請求項７記載の内燃機関の点火時期制御装置。
【請求項９】
前記運転安定性予測手段は、
各々の前記模擬筒内圧履歴データごとに、そのサイクルの図示トルクを算出する図示トルク算出手段と、
前記図示トルクの分布のバラツキ度合いを算出するバラツキ度合い算出手段と、
前記バラツキ度合いが許容レベル以下であるか否かを判別するトルク変動判定手段と、
を含むことを特徴とする請求項１乃至６の何れか１項記載の内燃機関の点火時期制御装置。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、内燃機関の点火時期制御装置に関する。
【背景技術】
【０００２】
火花点火内燃機関においては、点火時期が運転状態に多大な影響を及ぼす。通常運転かつ定常状態では、一般に、ノッキングが生じない範囲において出力トルクが最大となるような点火時期が、最適な点火時期であると言える。
【０００３】
最適な点火時期は、機関回転数や、負荷率（吸気圧）、冷却水温などのパラメータによって変化する。そこで、それらのパラメータと最適点火時期との関係を定めたマップを作っておき、そのマップを参照して点火時期を制御することが従来から行われている。
【０００４】
最適点火時期は、更に、機差バラツキや、経年変化によっても、違いが生じる。このため、予め用意された上記マップを参照するだけでは、真の最適点火時期は得られない。そこで、ノッキングを検知するノックセンサを用いたノッキングフィードバック制御が従来から行われている。この制御では、ノックセンサでノッキングが検知されたら、ノッキングが検知されなくなるまで点火時期を少しずつ遅角していく。そして、ノッキングが検知されなくなったら、点火時期を少しずつ進角していき、ノッキングが再発した場合は、再度点火時期を遅角していく。
【０００５】
また、特開平９−２５０４３５号公報には、所定クランク角までの実際の燃焼割合を検知し、その燃焼割合が目標値となるように、点火時期をフィードバック制御する技術が開示されている。
【０００６】
【特許文献１】特開平９−２５０４３５号公報
【特許文献２】特開平９−３１７５２２号公報
【特許文献３】特開昭６４−８７８７４号公報
【特許文献４】特開平１０−１１０６３８号公報
【発明の開示】
【発明が解決しようとする課題】
【０００７】
しかしながら、ノッキングフィードバック制御は、上述した通り、ノッキングを検知することを前提としている。つまり、一旦はノッキングを発生させなければ、最適な点火時期を見つけ出すことができない。このため、ノッキングによる内燃機関へのダメージや騒音を完全に防止することができる訳ではない。これと同様の問題は、特開平９−２５０４３５号公報に開示された技術についても、残っている。
【０００８】
この発明は、上記の点に鑑みてなされたものであり、点火時期を調整しようとする際に、点火時期調整後の運転安定性を事前に予測することのできる内燃機関の点火時期制御装置を提供することを目的とする。
【課題を解決するための手段】
【０００９】
第１の発明は、上記の目的を達成するため、内燃機関の点火時期制御装置であって、
内燃機関の点火時期を現在値から変更したと仮定した場合の、クランク角に対する筒内圧履歴を模擬した模擬筒内圧履歴データを生成する模擬筒内圧履歴データ生成手段と、
多数サイクル分の前記模擬筒内圧履歴データに基づいて、点火時期変更後の運転安定性を事前に予測する運転安定性予測手段と、
前記運転安定性予測手段の予測結果に基づいて、実際の点火時期を設定する点火時期設定手段と、
を備えることを特徴とする。
【００１０】
また、第２の発明は、第１の発明において、
前記模擬筒内圧履歴データ生成手段は、
所定の燃焼パラメータに基づいて燃焼状態を模擬する燃焼モデル手段と、
点火時期変更後の前記燃焼パラメータの分布の標準偏差を予測する標準偏差予測手段と、
前記標準偏差予測手段の予測値を標準偏差とする正規分布から、前記燃焼パラメータの値をランダムに抽出する燃焼パラメータ抽出手段と、
前記燃焼パラメータ抽出手段により抽出された値を前記燃焼モデル手段に入力することで得られる模擬の燃焼状態に基づいて、前記模擬筒内圧履歴データを算出する筒内圧履歴算出手段と、
を含むことを特徴とする。
【００１１】
また、第３の発明は、第２の発明において、
前記燃焼パラメータは、燃焼効率に関する値と、所定の燃焼割合となるクランク角に関する値と、第１の燃焼割合から第２の燃焼割合に変化するまでの燃焼期間に関する値とを含むことを特徴とする。
【００１２】
また、第４の発明は、第２または第３の発明において、
前記標準偏差予測手段は、
現在の運転状態における前記燃焼パラメータの現実値の平均値を取得する平均値取得手段と、
前記平均値取得手段により取得された平均値と現在の運転状態を表す所定の運転パラメータとに基づいて前記燃焼パラメータの標準偏差予測値を導出するための演算式を記憶した手段と、
現在の運転状態における前記運転パラメータの現実値を取得する運転パラメータ取得手段と、
を含むことを特徴とする。
【００１３】
また、第５の発明は、第４の発明において、
前記平均値取得手段は、
実筒内圧を検出する筒内圧センサと、
クランク角に対する前記実筒内圧の履歴に基づいて、現在の運転状態における前記燃焼パラメータの現実値を各サイクルごとに算出する燃焼パラメータ算出手段と、
前記燃焼パラメータ算出手段の算出値の種別ごとの算術平均を算出する算術平均手段と、
を含むことを特徴とする。
【００１４】
また、第６の発明は、第５の発明において、
前記燃焼パラメータ算出手段により算出された前記燃焼パラメータの現実値の分布から、前記燃焼パラメータの実標準偏差を算出する実標準偏差算出手段と、
前記標準偏差予測値を導出するための演算式に、前記実標準偏差と、前記運転パラメータの現実値とを入力することにより、前記演算式内の係数を較正する較正手段を更に備えることを特徴とする。
【００１５】
また、第７の発明は、第１乃至第６の発明の何れかにおいて、
前記運転安定性予測手段は、
各々の前記模擬筒内圧履歴データごとに、そのサイクルがノックを起こすノックサイクルであるか否かを判定するノックサイクル判定手段と、
前記ノックサイクルの発生頻度が許容レベル以下であるか否かを判別するノッキング判定手段と、
を含むことを特徴とする。
【００１６】
また、第８の発明は、第７の発明において、
前記ノックサイクル判定手段は、
前記模擬筒内圧履歴データに基づいて、そのサイクルの自己着火時点を予測する自己着火予測手段と、
前記予測された自己着火時点での燃焼割合が判定値以下である場合に、そのサイクルがノックサイクルであると判定する燃焼割合算出判定手段と、
を含むことを特徴とする。
【００１７】
また、第９の発明は、第１乃至第６の発明の何れかにおいて、
前記運転安定性予測手段は、
各々の前記模擬筒内圧履歴データごとに、そのサイクルの図示トルクを算出する図示トルク算出手段と、
前記図示トルクの分布のバラツキ度合いを算出するバラツキ度合い算出手段と、
前記バラツキ度合いが許容レベル以下であるか否かを判別するトルク変動判定手段と、
を含むことを特徴とする。
【発明の効果】
【００１８】
第１の発明によれば、点火時期変更後のクランク角に対する筒内圧履歴を模擬した模擬筒内圧履歴データを、予測的に生成することができる。この模擬筒内圧履歴データを多数サイクル分集めると、点火時期変更後に予測される燃焼変動の様子を表現することができる。運転安定性は、燃焼変動に深く関係している。この発明によれば、点火時期変更後の運転安定性を、燃焼変動の様子に基づいて、精度良く事前に予測することができる。そして、この発明によれば、運転安定性の予測結果に基づいて、運転安定性を損う事態を招くことのないように、点火時期を最適な時期に調整することができる。また、機差バラツキや経年変化の影響を受けることなく、確実に上記効果を奏することができる。
【００１９】
第２の発明によれば、燃焼状態を模擬する燃焼モデル手段に入力すべき燃焼パラメータの、点火時期変更後の標準偏差を予測することができる。そして、その予測値を標準偏差とする正規分布からランダムに抽出した燃焼パラメータの値を燃焼モデル手段に入力することにより、精度の高い模擬筒内圧履歴データを生成することができる。そのため、この発明によれば、点火時期変更後の運転安定性をより良い精度で予測することができる。
【００２０】
第３の発明によれば、燃焼モデル手段は、燃焼効率に関する値と、所定の燃焼割合となるクランク角に関する値と、第１の燃焼割合から第２の燃焼割合に変化するまでの燃焼期間に関する値とに基づいて、燃焼状態を模擬することができる。このため、燃焼状態をより実態に近く模擬することができる。それゆえ、この発明によれば、点火時期変更後の運転安定性をより良い精度で予測することができる。
【００２１】
第４の発明によれば、点火時期変更後の燃焼パラメータの標準偏差をより良い精度で予測することができる。それゆえ、この発明によれば、更に精度良く模擬筒内圧履歴データを生成することができ、ひいては、点火時期変更後の運転安定性を更に良い精度で予測することができる。
【００２２】
第５の発明によれば、実筒内圧を検出する筒内圧センサの出力に基づいて、現在の運転状態における燃焼パラメータの現実値の平均値を実測することができる。そして、その実測平均値を基礎として、点火時期変更後の燃焼パラメータの標準偏差を予測するので、その予測精度を更に向上することができる。
【００２３】
第６の発明によれば、実測された燃焼パラメータの現実値の分布から、燃焼パラメータの実標準偏差を取得することができる。そして、この発明によれば、その実標準偏差の値を用いて、標準偏差予測値を導出するための演算式内の係数を較正することができる。それゆえ、燃焼パラメータの標準偏差予測値を導出するための演算式に、機差バラツキや経年変化の影響を反映させることができる。このため、機差バラツキや経年変化の影響を更に正確に織り込むことができるので、燃焼変動の様子をより良い精度で予測することができる。
【００２４】
第７の発明によれば、運転安定性として、ノッキングが起きるか否かを予測することができる。このため、この発明によれば、ノッキングを一度も生じさせることなく、点火時期を最適な時期に調整することができる。
【００２５】
第８の発明によれば、ノッキングが起きるか否かを、より良い精度で予測することができる。
【００２６】
第９の発明によれば、運転安定性として、トルク変動の大きさが許容レベル以下に収まるか否かを予測することができる。このため、この発明によれば、トルク変動による弊害を一度も生じさせることなく、点火時期を最適な時期に調整することができる。
【発明を実施するための最良の形態】
【００２７】
実施の形態１．
［システム構成の説明］
図１は、本発明の実施の形態１のシステム構成を説明するための図である。図１に示すように、本実施形態のシステムは、火花点火式の内燃機関１０を備えている。内燃機関１０には、クランク角を検出するクランク角センサ１２が組み込まれている。クランク角センサ１２は、クランク軸が所定回転角だけ回転する毎に、Hi出力とLo出力を反転させるセンサである。クランク角センサ１２の出力によれば、クランク角（クランク軸の回転位置）や、機関回転数NEなどを検知することができる。
【００２８】
更に、内燃機関１０には、冷却水温THWを検出する水温センサ１６と、筒内圧センサ１８とが組み込まれている。筒内圧センサ１８によれば、筒内（燃焼室内）に生ずる圧力を検出することができる。
【００２９】
内燃機関１０の吸気通路１９の途中には、サージタンク２０が設けられている。サージタンク２０には、その内部の圧力、すなわち吸気管圧力を検出する吸気圧センサ２１が設置されている。吸気管圧力は、内燃機関１０の１行程当たりの吸入空気量に比例する。従って、吸気圧センサ２１の出力によれば、内燃機関１０の負荷率KL［％］を検出することができる。負荷率KLとは、標準状態で総行程容積を占める新気の質量に対する、筒内に吸入された新気の質量の比率を百分率で表したものである。
【００３０】
また、吸気通路１９には、その内部を流通する吸入空気量GAを検出するエアフロメータ２２が配置されている。エアフロメータ２２の下流には、スロットル弁２４が配置されている。スロットル弁２４の近傍には、スロットル開度TAを検出するスロットル開度センサ２６が組み付けられている。
【００３１】
内燃機関１０の吸気ポートには、その内部にガソリン等の燃料を噴射するための燃料噴射弁２８が配設されている。また、内燃機関１０には、燃焼室内の混合気に点火するための点火プラグ３０が設置されている。更に、内燃機関１０の排気通路３２には、その内部の圧力、すなわち排気管圧力を検出する排気圧センサ３３が設置されている。また、排気通路３２には、排気ガスを浄化するための触媒３４が組み込まれている。
【００３２】
内燃機関１０は、可変動弁機構３６を備えている。可変動弁機構３６によれば、吸気弁３８のバルブタイミングを進角したり遅角したりすることができる。可変動弁機構３６の近傍には、吸気カム軸の回転位置を検出するカム角センサ４０を備えている。カム角センサ４０の出力によれば、吸気弁３８のバルブタイミングの実進角量を検知することができる。以下、吸気弁３８のバルブタイミングの実進角量を「VVT進角量」と称し、符号VT［degCA］で表す。VVT進角量VTは、吸気弁３８のバルブタイミングが最遅角状態のとき、０degCAとされる。
【００３３】
本実施形態のシステムは、ECU(Electronic Control Unit)５０を備えている。ECU５０には、上述した各種のセンサからセンサ信号が供給されている。ECU５０は、それらのセンサ信号に基づいて、燃料噴射弁２８や、点火プラグ３０、可変動弁機構３６などの各種アクチュエータを制御することができる。
【００３４】
［実施の形態１の基本思想］
本実施形態では、ECU５０は、内燃機関１０が通常運転状態かつ定常状態にあるとき、点火時期が最適な時期になるように制御する。この場合の最適な点火時期とは、ノッキングが発生しない範囲において、内燃機関１０のトルクが最大となる点火時期である。なお、通常運転とは、触媒暖気制御等の特別な制御を行っていない運転状態のことを意味する。以下では、通常運転かつ定常状態のことを、単に「定常運転」と称することがある。
【００３５】
一般に、点火時期以外の運転条件が一定の下で、内燃機関１０のトルクを最大にするような点火時期が存在する。その点火時期は、MBT(Minimum spark advance for Best Torque)と称される。一方、ノッキングは、点火時期を進角させるほど、発生し易くなる。
【００３６】
一般に、ノッキングが発生し始めるノック点火時期と、MBTとは、互いに近接している。そして、ノック点火時期やMBTは、運転条件によって変化するほか、経年によっても変化してくる。更に、ノック点火時期やMBTには、機差バラツキ（個体差）も存在する。また、ノック点火時期とMBTとのどちらが進角側に存在するかという順番も、定まってはいない。
【００３７】
従来からあるノッキングフィードバック制御の思想は、ノックセンサでノッキングが検出されるところまで点火時期を実際に徐々に進角させていくことによって、最適な点火時期を得ようとするものであった。しかしながら、この方法では、最適な点火時期に調整される前に、一旦は必然的にノッキングが発生することとなる。このため、ノッキングによる騒音や機関へのダメージを完全には排除することができなかった。
【００３８】
そこで、本実施形態では、点火時期を実際に進角させるのに先立って、点火時期進角後にノッキングが発生するかどうかを事前に予測することとした。そして、点火時期を進角してもノッキングは発生しないとの予測結果が得られた場合にのみ、点火時期を実際に進角させることとした。これにより、点火時期を最適時期に調整するに際して、ノッキングが発生することを未然に防止することができる。
【００３９】
［燃焼変動］
内燃機関１０の運転においては、一般に、燃焼変動が生ずる。燃焼変動とは、筒内での混合気の燃焼の様子がサイクルごとに変化することを言う。燃焼変動は、過渡運転状態だけでなく、定常運転状態においても生ずる。燃焼変動が生ずる原因は、空気と燃料の混合比や、前サイクルの既燃ガスの残留量、筒内での混合気の流動状態（乱れ）などが、サイクルごとに変化するためであると考えられている。
【００４０】
図２は、クランク角に対する筒内圧履歴の実測データを示すグラフである。図２中では、圧縮上死点を０degCAとしている。図２のグラフには、機関回転数NEや負荷率KLなどを一定に保った定常運転中に実測した、多数のサイクルの筒内圧履歴が重ねて描かれている。図２によれば、定常運転状態においても、燃焼期間中の筒内圧履歴（燃焼圧履歴）を示す波形が、燃焼変動のためにサイクルごとに変動していることが分かる。以下では、クランク角に対する筒内圧履歴のデータを「筒内圧履歴データ」と称する。
【００４１】
［ノッキング］
火花点火機関における燃焼は、点火プラグ３０により点火されて生じた火炎が周囲に伝播していくことによって進行する。ノッキングは、火炎が未だ到達していない所の混合気が自己着火することで生じる圧力波によって引き起こされる現象である。
【００４２】
上述した燃焼変動により、燃焼状態はサイクルごとに変化している。このため、ノッキングが発生しているときであっても、サイクルごとに見れば、ノックが起きているサイクルと、そうでないサイクルとがある。そして、ノックが起きているサイクル（以下「ノックサイクル」と称する）の発生頻度が多くなるほど、ノッキングの強度が強まるものと考えられる。
【００４３】
このように、ノッキングが発生するかどうかは、燃焼変動の様子と深く関係している。それゆえ、点火時期を現在値から進角したと仮定した場合の燃焼変動の様子を予測することができれば、点火時期進角後にノッキングが生ずるかどうかを事前に予測する基礎が築かれることとなる。燃焼変動の様子は、図２のように、多数サイクル分の筒内圧履歴データを集めることによって、表現することができる。
【００４４】
そこで、本実施形態では、点火時期を現在値から変更したと仮定した条件下での模擬的な筒内圧履歴データを作成することとした。このデータを以下「模擬筒内圧履歴データ」と称する。この模擬筒内圧履歴データを多数サイクル分集めることにより、点火時期変更後の燃焼変動の様子を事前に予測することができる。模擬筒内圧履歴データは、以下に説明する燃焼変動予測モデルによって、生成することができる。
【００４５】
［燃焼変動予測モデル］
本燃焼変動予測モデルでは、熱発生パターンの近似関数として知られているWiebe関数を用いて、筒内での燃焼状態を模擬する。Wiebe関数は、次式で表される。
【００４６】
【数１】
2007056727.tif 000003
【００４７】
上記（１）式は、燃焼期間中に成立する。この式中の各記号の意味は、下記の通りである。
Q：累積加熱率［J］
θ：クランク角［deg］
k：燃焼効率
Qf：燃焼に用いた燃料（筒内に吸入された燃料）の発熱量［J］
ａ＝−ｌｎ（１−０．９９９）
m：形状パラメータ
θ０：燃焼開始点［degCA］
θ０−１００：０−１００％燃焼期間［degCA］
【００４８】
上記の燃料発熱量Qfは、燃料の種類によって定まる。本モデルにおいては、この燃料発熱量Qfは既定値とされる。また、燃焼効率kは、累積加熱率Qの最大値Qmaxと、燃料発熱量Qfとを用いて、k＝Qmax／Qfで表される。すなわち、燃焼効率kは、そのサイクルの燃焼が完全燃焼にどれだけ近いかを示す値である。０−１００％燃焼期間θ０−１００は、燃焼割合Xbが０％から１００％に変化するまでの期間をクランク角で表したものである。燃焼割合Xbとは、燃焼の進行状態を表す量である。すなわち、０−１００％燃焼期間θ０−１００とは、燃焼開始点θ０から燃焼終了点までの全燃焼期間をクランク角で表したものである。
【００４９】
上記（１）式中の０−１００％燃焼期間θ０−１００および燃焼開始点θ０は、更に、次式で表される。
【００５０】
【数２】
2007056727.tif 000004
【００５１】
上記（２）式中、θ１０−７０は、１０−７０％燃焼期間［degCA］を意味する。１０−７０％燃焼期間θ１０−７０とは、燃焼割合Xbが１０％から７０％に変化するまでの期間をクランク角で表したものである。また、上記（３）式中、θ５０は、５０％燃焼点［degCA］を意味する。５０％燃焼点θ５０とは、燃焼割合Xbが５０％となる時点のクランク角である。
【００５２】
Wiebe関数は、公知の方法、つまり熱力学上のエネルギー保存則と組み合わせる方法により、そのサイクルの燃焼期間中の筒内圧履歴に変換することができる。この方法については、後述する。
【００５３】
上記（２）、（３）式を用いた場合には、Wiebe関数に代入すべきパラメータは、燃焼効率k、形状パラメータm、５０％燃焼点θ５０、および１０−７０％燃焼期間θ１０−７０の４つとなる。燃焼変動は、これら４つのパラメータがサイクルごとにばらつくことに相当するものと考えられる。換言すれば、燃焼変動は、Wiebe関数をサイクルごとに変化させることで表現（模擬）することができると考えられる。
【００５４】
ただし、本発明の発明者らの研究によれば、上記４つのパラメータのうち、形状パラメータmは、燃焼変動に有意な影響を与えないことが明らかとなった。このため、本実施形態では、形状パラメータmは既定値とする。そして、残りの燃焼効率k、５０％燃焼点θ５０、および１０−７０％燃焼期間θ１０−７０の３つのパラメータのサイクルごとのバラツキ具合を考慮することとした。以下では、上記３つのパラメータのことを「燃焼パラメータ」と称する。
【００５５】
図３中の棒グラフは、定常運転状態における、サイクルごとの燃焼効率ｋの分布の実測データを示している。この実測データは、図３中の実線で示される正規分布にほぼ一致している。
また、図４中の棒グラフは、定常運転状態における、サイクルごとの５０％燃焼点θ５０の分布の実測データを示している。この実測データは、図４中の実線で示される正規分布にほぼ一致している。
そして、図５中の棒グラフは、定常運転状態における、サイクルごとの１０−７０％燃焼期間θ１０−７０の分布の実測データを示している。この実測データは、図５中の実線で示される正規分布にほぼ一致している。
【００５６】
図６は、定常運転状態における５０％燃焼点θ５０および１０−７０％燃焼期間θ１０−７０の実測値の相関図である。図７は、定常運転状態における燃焼効率kと５０％燃焼点θ５０の実測値の相関図である。図８は、定常運転状態における燃焼効率kと１０−７０％燃焼期間θ１０−７０の実測値の相関図である。
【００５７】
図６から、５０％燃焼点θ５０と１０−７０％燃焼期間θ１０−７０とは、相関を有することが分かる。これに対し、図７からは、燃焼効率kと５０％燃焼点θ５０との間にはほぼ相関がないことが分かり、図８からは、燃焼効率kと１０−７０％燃焼期間θ１０−７０との間にはほぼ相関がないことが分かる。
【００５８】
多数サイクル分の模擬筒内圧履歴データを作成するには、点火時期を変更した場合の燃焼変動を表現する多数のWiebe関数を生成する必要がある。それには、Wiebe関数に代入すべき３つの燃焼パラメータのサイクルごとのバラツキを予測する必要がある。
【００５９】
一方、図３乃至図８に示される実測データから、３つの燃焼パラメータのサイクルごとのバラツキは、θ５０とθ１０−７０との相関を加味した多変数正規分布にほぼ従っていることが明らかである。そこで、本実施形態では、この多変数正規分布から、燃焼パラメータk、θ５０、θ１０−７０の組をランダムに抽出し、それらの値をWiebe関数に代入することとした。そして、このようなWiebe関数を多数生成し、その各々を筒内圧履歴に変換することで、多数サイクル分の模擬筒内圧履歴データを生成することとした。
【００６０】
本実施形態では、このようにして得られる多数サイクル分の模擬筒内圧履歴データと、後述する自己着火予測式とを組み合わせることにより、点火時期変更後にノッキングが発生するかどうかを事前に予測することとした。
【００６１】
［実施の形態１における具体的処理］
図９は、上記の機能を実現するために本実施形態においてECU５０が実行するルーチンのフローチャートである。なお、本ルーチンは、クランク角に同期してサイクルごとに実行されるものとする。図９に示すルーチンによれば、まず、現在の運転状態を表す各種のパラメータが計測され、保存される（ステップ１００）。具体的には、機関回転数NE、負荷率KL、およびVVT進角量VTが計測され、保存される。
【００６２】
次に、現在の運転状態が、通常運転かつ定常状態か否かが判別される（ステップ１０２）。本実施形態では、触媒暖気制御等の特別な制御を行っていない通常運転をしているときの点火時期を制御することを目的とする。このため、通常運転状態でなかった場合には、上記ステップ１００で保存したデータをリセットし（ステップ１０４）、その後、今回の処理サイクルを速やかに終了する。
【００６３】
また、本実施形態では、内燃機関１０が定常状態にあることも前提とされる。上記ステップ１０２において、定常状態であるか否かは、ステップ１００で取得された機関回転数NE、負荷率KL、およびVVT進角量VTの計測値が、前回の処理サイクルでの計測値と同値であるか否かによって判別される。その判別の結果、定常状態でなかった場合には、上記ステップ１００および後述するステップ１１０で保存されたデータをリセットし（ステップ１０４）、その後、今回の処理サイクルを速やかに終了する。
【００６４】
一方、上記ステップ１０２において、通常運転かつ定常状態であることが認めれた場合には、カウント値cycの値を一つインクリメントする。本実施形態では、連続した５０サイクルにおける、燃焼パラメータk、θ５０、θ１０−７０の平均値を算出する。カウント値cycは、その５０サイクルをカウントするための値であり、その初期値は０とされている。
【００６５】
次に、今サイクルにおける、クランク角に対する実筒内圧の履歴が計測される（ステップ１０８）。具体的には、所定クランク角Δθごとに、筒内圧センサ１８の出力をサンプリングし、得られた実筒内圧値を記憶する。上記Δθは、例えば１〜３degCA程度とされる。また、ここでは、実筒内圧は、７２０degCAの全範囲に渡って計測しなくてもよく、例えば、吸気弁閉弁時から排気弁開弁時までのクランク角範囲で計測すればよい。
【００６６】
次に、上記ステップ１０８で計測された実筒内圧の履歴に基づいて、今サイクルにおける３つの燃焼パラメータk、θ５０、θ１０−７０の現実値を算出する（ステップ１１０）。以下、この算出処理について説明する。この算出処理には、次式が用いられる。
【００６７】
【数３】
2007056727.tif 000005
【００６８】
上記（４）および（５）式中、Pは筒内圧［Pa］を、Vは筒内容積［ｍ３］を、κは比熱比を、それぞれ表す。この両式は、熱力学上のエネルギー保存則を、累積加熱率Q、筒内圧P、筒内容積V、およびクランク角θで表したものに相当する。また、上記（５）式は、上記（４）式の両辺をクランク角θで積分したものに相当する。
【００６９】
上記両式中のVおよびdV／dθは、クランク角θに応じて幾何学的に決定される。つまり、VおよびdV／dθは、クランク角θの関数である。その関係を用いると、上記両式からVを消去することができる。ECU５０は、そのVが消去された状態の上記両式を用いて演算を行う。また、演算の際、比熱比κは一定の既定値とされる。
【００７０】
上記（５）式の右辺の積分は、点火時のクランク角から開始される。そして、あるクランク角θまでこの積分を計算した値を、以下、Q(θ)と表す。ここでのθは、点火時と排気弁閉弁時との間において定義される。
【００７１】
ステップ１１０においては、まず、上記ステップ１０８で計測された実筒内圧の履歴に対するQ(θ)の値が、所定クランク角ごとに、数値計算により近似的に算出される。すなわち、上記（５）式の右辺のPおよびdP／dθに、上記ステップ１０８で計測された各点での筒内圧値およびその変化率を代入した値が、積算される。
【００７２】
このようにして算出される累積加熱率Q(θ)は、燃焼が継続している間は、その燃焼熱により、増大を続ける。燃焼が終了した後は、燃焼室内壁やシリンダ内壁へ熱が逃げることの損失により、累積加熱率Q(θ)は低下に転じる。よって、累積加熱率Q(θ)は、点火時と排気弁閉弁時との間の、あるクランク角において、最大値Qmaxをとる。ステップ１１０においては、所定クランク角ごとに算出されたQ(θ)のうちから、最大値Qmaxが求められる。
【００７３】
燃焼効率kは、既述した通り、k＝Qmax／Qfと表すことができる。この式に従って、今サイクルの燃焼効率kの現実値が算出される。なお、燃料発熱量Qfは、既定値である燃料の特性値と、現在の空燃比および筒内吸入新気量に基づいて算出することができる。
【００７４】
クランク角θにおける燃焼割合Xbは、Xb＝Q(θ)／Qmax×１００として算出することができる。よって、今サイクルの５０％燃焼点θ５０の現実値は、Q(θ)／Qmax＝０．５が成立するようなクランク角として求めることができる。また、１０−７０％燃焼期間θ１０−７０の現実値は、Q(θ)／Qmaxが０．１から０．７に変化するまでのクランク角範囲として求めることができる。
【００７５】
ステップ１１０においては、以上のようにして、３つの燃焼パラメータk、θ５０、θ１０−７０の現実値が算出され、その算出値が保存される。次に、カウント値cycが５０を超えたか否かが判別される（ステップ１１２）。この判別は、上記ステップ１１０で算出された燃焼パラメータの実測値データが、５０サイクル分、蓄積されたか否か判別に相当する。
【００７６】
ステップ１１２の判別の結果、カウント値cycが５０に達していない場合には、今回の処理サイクルが速やかに終了される。この場合、内燃機関１０の次の作動サイクルに同期して、図９に示すルーチンが再度実行される。そのとき、定常運転が継続していれば、上記ステップ１０６〜１０８の処理がまた行われて、燃焼パラメータの実測値データが蓄積される。これに対し、運転条件が変わっていた場合には、前回までに蓄積された燃焼パラメータの実測値データは、上記ステップ１０４においてリセットされ、カウント値cycも０にリセットされる。
【００７７】
定常運転が５０サイクル以上継続すると、ステップ１１２において、cyc＞５０の成立が認められる。この場合には、５０サイクル分の燃焼パラメータの実測値データが蓄積されていることになる。この場合には、次に、その５０サイクル分の燃焼パラメータk、θ５０、θ１０−７０の実測値データの平均値（算術平均）が算出される（ステップ１１４）。ここで算出される燃焼効率kの平均値をkc、５０％燃焼点θ５０の平均値をBPc、１０−７０％燃焼期間θ１０−７０の平均値をBDc、とそれぞれ表記する。
【００７８】
本実施形態では、ノッキングが発生するか否かを予測する対象とする仮想上の点火時期を、想定点火時期SAiと称する。また、想定点火時期SAiが実現された場合の５０％燃焼点θ５０の平均値の予測値を、想定５０％燃焼点平均値BPiと表す。
【００７９】
上記ステップ１１４の処理が終わったら、次に、この想定点火時期SAiおよび想定５０％燃焼点平均値BPiに、それらの現在値が、基準値として代入される（ステップ１１６）。すなわち、想定点火時期SAiには、点火時期の現在値SAcが代入される。点火時期の現在値SAcは、ECU５０が現在設定している点火時期の値である。また、想定５０％燃焼点平均値BPiには、上記ステップ１１４で算出された、５０サイクル分の実測値の平均値BPcが代入される。
【００８０】
次に、想定点火時期SAiおよび想定５０％燃焼点平均値BPiが、それぞれ、２degCA小さい値に更新される（ステップ１１８）。このステップ１１８の処理は、想定点火時期SAiを２degCA進角させることに相当する。本発明者らの知見によれば、点火時期を移動した場合、これに伴って、５０％燃焼点θ５０の平均値もほぼ同じ分だけ移動する。よって、点火時期を２degCA進角させた場合、５０％燃焼点θ５０の平均値も２degCA進角すると予測される。この現象を表現するため、ステップ１１８においては、想定点火時期SAiだけでなく、想定５０％燃焼点平均値BPiも２degCA小さい値に更新される。
【００８１】
ステップ１１８の処理が初めて行われた場合、その処理後には、想定点火時期SAiおよび想定５０％燃焼点平均値BPiは、それぞれ、現在値より２degCA進角側の値を表すこととなる。
【００８２】
次に、想定点火時期SAiにおける筒内圧履歴のサイクル変動を予測する（ステップ１２０）。具体的には、想定点火時期SAiにおける模擬筒内圧履歴データを１００サイクル分生成する。図１０は、ステップ１２０の処理を行うためにECU５０が実行するサブルーチンのフローチャートである。
【００８３】
図１０に示すサブルーチンでは、まず、想定点火時期SAiにおける、各燃焼パラメータk、θ５０、θ１０−７０の分布の標準偏差と、５０％燃焼点θ５０と１０−７０％燃焼期間θ１０−７０との共分散とを予測する（ステップ１２２）。具体的には、次の関数func１〜func４に基づいて、それらの値を算出する。
【００８４】
σ１／μ１＝func１（NE，KL，KL２）      ・・・（６）
σ２／μ３＝func２（NE，KL，KL２，VT２）      ・・・（７）
σ３／μ３＝func３（NE，NE２，KL，KL２，VT）      ・・・（８）
σ２３／μ３２＝func４（NE，NE２，KL，KL２，VT）      ・・・（９）
【００８５】
上記（６）〜（９）式中、σ１は燃焼効率kの標準偏差、σ２は５０％燃焼点θ５０の標準偏差、σ３は１０−７０％燃焼期間θ１０−７０の標準偏差、σ２３は５０％燃焼点θ５０と１０−７０％燃焼期間θ１０−７０との共分散、μ１は燃焼効率kの平均値、μ３は１０−７０％燃焼期間θ１０−７０の平均値、NEは機関回転数、KLは負荷率、VTはVVT進角量、をそれぞれ表す。
【００８６】
上記（６）〜（９）式は、予め行われた実験による計測データから求められた重回帰式とされる。すなわち、内燃機関１０において機関回転数NE、負荷率KLおよびVVT進角量VTを異ならせた種々の運転条件下で、標準偏差σ１、σ２、σ３および共分散σ２３を実測し、その実測データから回帰することにより、上記（６）〜（９）式を得ることができる。上記（６）〜（９）式の右辺は、カッコ内の変数の多項式となっている。
【００８７】
ステップ１２２においては、上記（６）〜（９）式に、NE、KL、VT、μ１、μ３の各々の値を代入することにより、標準偏差σ１、σ２、σ３および共分散σ２３が算出される。このとき、NE、KL、VTには、ステップ１００で計測された機関回転数NE、負荷率KL、VVT進角量VTの値が代入される。μ１には、ステップ１１４で算出されたkcの値が代入される。μ３には、ステップ１１４で算出されたBDcの値が代入される。
【００８８】
想定点火時期SAiの下での５０％燃焼点θ５０の平均値は、既述した通り、想定５０％燃焼点平均値BPiであると予測される。また、想定点火時期SAiの下での５０％燃焼点θ５０のバラツキの標準偏差は、上記ステップ１２２で算出されたσ２であると予測される。従って、想定点火時期SAiの下での５０％燃焼点θ５０のバラツキは、平均値をBPi、標準偏差をσ２とする正規分布、つまり図１１に示す正規分布に従うと予測される。
【００８９】
一方、想定点火時期SAiの下での燃焼効率kおよび１０−７０％燃焼期間θ１０−７０のバラツキは、次のように予測される。本発明者らの知見によれば、点火時期を移動した場合でも、燃焼効率kの平均値は、ほとんど変化しない。また、１０−７０％燃焼期間θ１０−７０の平均値は、点火時期に対して感度を有するが、２degCA程度の点火時期の変更によってはほとんど変化しないとしてよい。従って、想定点火時期SAiの下での燃焼効率kの平均値および１０−７０％燃焼期間θ１０−７０の平均値は、それぞれ、上記ステップ１１４で算出された現在値kcおよびBDcに等しいと予測される。また、想定点火時期SAiの下での燃焼効率kおよび１０−７０％燃焼期間θ１０−７０のバラツキの標準偏差は、それぞれ、上記ステップ１２２で算出されたσ１、σ３であると予測される。
【００９０】
以上より、想定点火時期SAiの下での燃焼効率kのバラツキは、平均値をkc、標準偏差をσ１とする正規分布、つまり図１２に示す正規分布に従うと予測される。また、想定点火時期SAiの下での１０−７０％燃焼期間θ１０−７０のバラツキは、平均値をBDc、標準偏差をσ３とする正規分布、つまり図１３に示す正規分布に従うと予測される。
【００９１】
上記ステップ１２２の次には、上述した図１１乃至図１３に示す正規分布（多変数正規分布）から、３つの燃焼パラメータk、θ５０、θ１０−７０の値が一組、ランダムに抽出される（ステップ１２４）。ただし、既述した通り、５０％燃焼点θ５０と１０−７０％燃焼期間θ１０−７０とは、相関を有している。想定点火時期SAiの下での両者の相関の度合いは、上記ステップ１２２で算出された共分散σ２３により予測される。よって、ステップ１２４における、５０％燃焼点θ５０および１０−７０％燃焼期間θ１０−７０の抽出に際しては、両者の共分散σ２３が織り込まれる。なお、ステップ１２４において、正規分布から値をランダムに抽出する方法は、特に限定されないが、例えば乱数を利用する方法をとることができる。
【００９２】
上記ステップ１２４の処理により、想定点火時期SAiの下で確率的に起こり得る一組のk、θ５０、θ１０−７０の値を生成することができる。次に、これらk、θ５０、θ１０−７０の値が、上記（２）および（３）式を介して、上記（１）式のWiebe関数に代入される（ステップ１２６）。このステップ１２６の処理によって得られたWiebe関数は、想定点火時期SAiの下で確率的に起こり得る、ある一つのサイクルの燃焼を模擬していると言える。
【００９３】
次に、ステップ１２６で得られたWiebe関数に基づいて、燃焼期間中の模擬筒内圧履歴、つまり点火時から燃焼終了時までの模擬筒内圧履歴が算出される（ステップ１２８）。具体的には、以下のような処理が行われる。まず、そのWiebe関数の、上記（１）式の右辺に相当するものと、上記（４）式の右辺とを等しいと置く。これにより、筒内圧Pの微分方程式が得られる。ECU５０が、この微分方程式を数値計算によって近似的に解くことにより、燃焼期間中の模擬筒内圧履歴が算出される。
【００９４】
図１４は、図１０に示すサブルーチンの処理によって作成される模擬筒内圧履歴データの一例をグラフ化して示す図である。図１４中では、圧縮上死点を０degCAとしている。上記ステップ１２８で算出された燃焼期間中の模擬筒内圧履歴は、上述したように、図１４中（３）の期間の模擬筒内圧履歴に相当する。
【００９５】
上記ステップ１２８の次には、図１４中（３）の燃焼期間中以外の模擬筒内圧履歴が算出される（ステップ１３０）。このステップ１３０において、燃焼終了時から排気弁開弁時(EVO)までの期間（図１４中（４）の期間）の模擬筒内圧履歴は、Woschiniの熱伝達モデルによる熱損失の式と、上記（４）式とに基づいて、算出される。この方法は公知であるので、ここではこれ以上の説明を省略する。そして、吸気弁閉弁時(IVC)から点火時までの期間（図１４中（２）の期間）の模擬筒内圧履歴は、ポリトロープ指数をサイクル平均値としたポリトロープ変化で近似して、算出される。この方法は公知であるので、ここではこれ以上の説明を省略する。また、吸気行程期間（図１４中（１）の期間）の模擬筒内圧履歴データは、吸気圧センサ２１で検出されるサージタンク圧の平均値に等しい一定値とされる。また、排気行程期間（図１４中（５）の期間）の模擬筒内圧履歴データは、その時に排気圧センサ３３で検出される排気管圧力の平均値に等しい一定値とされる。
【００９６】
上記ステップ１２４乃至１３０の処理により、図１４に例示するような１サイクル分の模擬筒内圧履歴データが生成される。本実施形態では、このような模擬筒内圧履歴データを多数サイクル分（ここでは１００サイクル分）集めて図２のようなデータを得ることにより、想定点火時期SAiの下での燃焼変動の様子が表現される。
【００９７】
上記ステップ１３０の次には、作成された模擬筒内圧履歴データが１００サイクル分に到達したか否かが判断される（ステップ１３２）。そして、作成された模擬筒内圧履歴データが１００サイクル分に達するまで、上記ステップ１２４乃至１３０の処理が繰り返し行われ、その後、図１０に示すサブルーチンの処理が終了される。
【００９８】
図１０に示すサブルーチンの処理によって１００サイクル分の模擬筒内圧履歴データが作成された後は、その１００サイクルのうち、ノックを起こすノックサイクルがいくつあるかを調査する処理が行われる（図９中のステップ１３４）。図１５は、ステップ１３４の処理を行うためにECU５０が実行するサブルーチンのフローチャートである。
【００９９】
図１５に示すサブルーチンでは、各々の模擬筒内圧履歴データに対して、ループ始端１３６およびループ終端１３８で示されるループ１の処理が行われる。すなわち、ループ１内の処理は、１００回（１００サイクル分）繰り返し行われる。ループ１の処理では、まず、対象とされた模擬筒内圧履歴データについて、その自己着火時点が算出される（ステップ１４０）。このステップ１４０の処理には、次式で表される自己着火予測式が用いられる。
【０１００】
【数４】
2007056727.tif 000006
【０１０１】
上記（１０）式中の各記号の意味は、下記の通りである。
ｔ：時間
P(t)：筒内圧［kgf/cm２］
Tu(t)：未燃ガス温度［Ｋ］
Φ：当量比
ON：燃料のオクタン価
A,B,C,F,Ea：適合定数（既定値）
【０１０２】
上記自己着火予測式によれば、筒内の未燃ガスが自己着火する時点を次のように予測することができる。上記（１０）式の左辺の積分を吸気弁閉弁時から開始し、その値が１になったとき、つまり上記（１０）式が成立したとき、自己着火が生じると予測される。
【０１０３】
ステップ１４０においては、上記（１０）式の左辺の積分を数値計算し、上記（１０）式が成立する時点を算出する。この計算において、上記（１０）式中に代入すべき値は、次のようにして求めることができる。筒内圧P(t)は、時間に対する筒内圧の履歴である。クランク角に対する筒内圧Pの履歴である模擬筒内圧履歴データは、ステップ１００で計測された機関回転数NEを用いて、P(t)に換算することができる。この換算により得られたP(t)が上記（１０）式に代入される。また、未燃ガス温度Tu(t)は、上記P(t)と、クランク角θから定まる筒内容積Vとを、気体の状態方程式に代入することで求められる。また、当量比Φは、ECU５０が現在設定している目標空燃比からの換算により、あるいは図示しない空燃比センサで検出される実空燃比からの換算により、求めることができる。燃料オクタン価ONおよび適合定数A,B,C,F,Eaは、予め定められた定数とされる。
【０１０４】
ここで、本実施形態におけるノックサイクルの判定方法について説明する。図１６は、クランク角に対する燃焼割合Xbの履歴を示すグラフである。本実施形態では、自己着火時点での燃焼割合Xbが、判定値（ここでは９０％とする）以下である場合には、ノックサイクルであると判定される。
【０１０５】
例えば、自己着火時点で、クランク角が図１６中のＡ点であった場合、燃焼割合Xbは９０％に達していない。このような場合には、自己着火時点で未燃ガスがまだ比較的多く残っていることになる。このため、自己着火が起きると、多量の未燃ガスが異常燃焼することよって強い圧力波が生じる。よって、このようなサイクルは、ノックを起こすと予測することができる。
【０１０６】
一方、自己着火時点で、クランク角が図１６中のＢ点であった場合、燃焼割合Xbは９０％を超えている。このような場合には、自己着火時点で未燃ガスはもうほとんど残っていないことになる。このため、自己着火が起きても、生ずる圧力波は弱い。よって、このようなサイクルは、ノックを起こさないと予測することができる。
【０１０７】
上記ステップ１４０で自己着火時点が算出された後は、上述したようなノックサイクル判定を行うべく、次のような処理が実行される。まず、自己着火時点における累積加熱率Qsが算出される（ステップ１４２）。具体的には、まず、自己着火時点の値を、機関回転数NEの値を用いて、自己着火時点のクランク角θsに換算する。次に、上記（５）式の積分を数値計算することにより、クランク角θsまで積分、すなわちQ(θs)を算出する。この算出処理は、上記ステップ１１０の算出処理に類似した処理である。ただし、この場合、上記（５）式中のPおよびdP／dθには、それぞれ、模擬筒内圧履歴データのPおよびその変化率が代入される。このようにして算出されたQ(θs)が、自己着火時点の累積加熱率Qsに相当する。
【０１０８】
次に、模擬筒内圧履歴データに基づいて、累積加熱率の最大値Qmaxを算出する（ステップ１４４）。具体的には、累積加熱率Q(θ)を、クランク角θsを超えて、排気弁閉弁時まで所定クランク角ごとに算出する。この算出処理は、上記ステップ１４２と同様にして計算可能である。このようにして算出されたクランク角ごとのQ(θ)のうちの最大値が、求めるQmaxに相当する。
【０１０９】
次に、Qs／Qmaxの値が０．９以下であるか否かが判別される（ステップ１４６）。Qs／Qmax≦０．９が成立することは、自己着火時点での燃焼割合Xbが９０％に達していないことを意味する。この場合には、前述したノックサイクル判定方法の考え方により、このサイクルがノックサイクルであると判断できる。この場合には、カウント値KNCを一つインクリメントし（ステップ１４８）、その後ループ１を終了する。カウント値KNCは、ノックサイクルの数を表す値である。このカウント値KNCの初期値は、０とされている。
【０１１０】
一方、Qs／Qmax≦０．９が不成立の場合には、自己着火時点での燃焼割合Xbが９０％を超えていることを意味する。従って、このサイクルはノックサイクルでないと判断できる。そこで、この場合には、カウント値KNCのインクリメント処理を飛ばして、ループ１を終了する。
【０１１１】
以上のようなループ１の処理が１００サイクル分の模擬筒内圧履歴データに対して繰り返し行われた後、図１５に示すサブルーチンの処理が終了される。そのとき、カウント値KNCは、その１００サイクルのうちのノックサイクルの数を表す値となっている。すなわち、カウント値KNCは、１００サイクル中のノックサイクルの発生頻度を表す値となっている。
【０１１２】
図１５に示すサブルーチンの処理終了後は、次に、カウント値KNCが許容値α以下であるか否かが判別される（図９中のステップ１５０）。既述した通り、ノックサイクルの頻度が高いほど、ノッキング強度は強まる。一方、ノックサイクルの頻度がある程度以下であれば、ノッキング強度は、機関へのダメージや騒音等の害が実際上生じないレベルに収まる。上記許容値αは、ノッキング強度が実害のない範囲に収まることを判定可能な値として、予め設定されている。
【０１１３】
以下では、説明の便宜上、ノッキング強度が実害を生じるレベルである場合を「ノッキング」と定義する。この定義によれば、カウント値KNCが許容値αを超えている場合は、ノッキングが発生すると予測されていることになり、カウント値KNCが許容値α以下である場合は、ノッキングが発生しないと予測されていることになる。
【０１１４】
上記ステップ１５０において、KNC≦αの不成立が認められた場合には、想定点火時期SAiが実現されるとノッキングが発生する、との予測が成り立つ。この場合には、想定点火時期SAiが２degCA遅角された値に更新される（ステップ１５２）。すなわち、想定点火時期SAiに２degCAが加算される。その後、その更新された想定点火時期SAiが実際の点火時期SAmとして設定される（ステップ１５４）。すなわち、ステップ１５４の処理以後、内燃機関１０は、点火時期SAm（＝SAi）で運転される。
【０１１５】
上記ステップ１５０の判別が行われたのが１回だけだった場合には、想定点火時期SAiは、上記ステップ１５２の処理により、現点火時期SAcに等しい値に戻っている。つまり、この場合には、実際の点火時期は変更されず、現点火時期SAcのままで内燃機関１０の運転が継続されることとなる。よって、この場合、実際にノッキングが発生することはない。
【０１１６】
一方、上記ステップ１５０において、KNC≦αの成立が認められた場合には、次に、想定点火時期SAiがMBTより進角側にあるか否か、つまりSAi＞MBTの成否、が判別される（ステップ１５６）。MBTと、機関回転数NEおよび負荷率KLとの関係は、予め実験で調査されており、既知となっている。ECU５０には、その関係を定めたマップが記憶されている。上記ステップ１５６では、そのマップを参照して算出されたMBTの値が用いられる。
【０１１７】
本実施形態では、実際の点火時期をMBTよりも進角させる実益はない。そこで、上記ステップ１５６において、想定点火時期SAiがMBTより進角側にあると判別された場合には、想定点火時期SAiがMBTの値に更新され（ステップ１５８）、その後、想定点火時期SAiが実際の点火時期SAmとして設定される（ステップ１５４）。すなわち、この場合には、ステップ１５４の処理以後、内燃機関１０は、MBTに等しい点火時期SAmで運転される。この場合には、上記ステップ１５０の処理において、点火時期がMBTより進角側にあってもノッキングが発生しないとの予測が成り立っている。よって、この場合、MBTで内燃機関１０を運転しても、実際上ノッキングが発生することはない。
【０１１８】
一方、上記ステップ１５６において想定点火時期SAiがMBTより遅角側にあると判別された場合には、現時点の想定点火時期SAiよりもさらに点火時期を進角させた方が、トルク向上の観点からは好ましいと言える。そこで、この場合には、現時点の想定点火時期SAiよりもさらに点火時期を進角させた場合にノッキングが発生するかどうかを予測するべく、ステップ１１８以下の処理が再度実行される。具体的には、ステップ１１８において、想定点火時期SAiおよび想定５０％燃焼点平均値BPiがさらに２degCA進角された値に更新される。そして、その更新された想定点火時期SAiの下でノッキングが発生するか否かの予測が再度行われる（ステップ１２０〜１５０）。
【０１１９】
このように、ステップ１５６からステップ１１８へと戻るループが設けられていることにより、ノッキングを起こさないと予測される点火時期のうちでなるべくMBTに近い点火時期を見つけ出すことができる。ステップ１５６からステップ１１８へと戻るループを通った場合には、ステップ１５０の判別が複数回行われることになる。複数回目のステップ１５０においてKNC≦αの不成立が認められた場合、すなわちノッキングが発生するとの予測がなされた場合には、ステップ１５２において想定点火時期SAiが２degCA遅角された値に更新され、その後、更新された想定点火時期SAiがステップ１５４において実際の点火時期SAmとして設定される。この場合に実際の点火時期SAmとして設定された値は、１回前のステップ１５０の処理においてノッキングを起こさないと予測された点火時期の値に等しい。よって、この場合、実際上ノッキングが発生することはない。
【０１２０】
以上説明したように、本実施形態によれば、点火時期を進角させようとする際、点火時期進角後にノッキングが発生するかどうかを事前に予測することができる。そして、ノッキングが発生しないとの予測が得られた場合に、点火時期を実際に進角させる。このため、実際のノッキングを一瞬たりとも起こさせることなく、点火時期を最適値に制御することができる。よって、ノッキングによる内燃機関１０へのダメージや騒音を従来よりも更に低減することができる。
【０１２１】
また、本実施形態では、筒内圧センサ１８を用いて実測された燃焼パラメータに基づいて、模擬筒内圧履歴データを生成する。筒内圧センサ１８を用いて実測された燃焼パラメータは、機差バラツキや経年変化の影響を受けた、内燃機関１０の現在の状態を反映している。このため、模擬筒内圧履歴データにも内燃機関１０の機差バラツキや経年変化の影響が反映される。従って、本実施形態によれば、点火時期変更後の燃焼変動の様子を、機差バラツキや経年変化に応じて、正確に予測することができる。このため、機差バラツキや経年変化によらず、上述した効果を奏することができる。
【０１２２】
ところで、上述した実施の形態１では、燃料をポート内に噴射する内燃機関について説明したが、本実施形態は、燃料を筒内に直接噴射する内燃機関の制御装置にも適用することができる。
【０１２３】
また、上述した実施の形態１では、ステップ１１８において想定点火時期SAiを２degCA進角させた値に更新することとしているが、ここで進角させる値は２degCAに限定されるものではない。更に、ここで進角させる値を固定値とせずに、機関回転数NE等の運転条件に応じて変化する値としてもよい。
【０１２４】
また、上述した実施の形態１では、上記（１０）式を利用して自己着火時点を予測することとしているが、自己着火時点を予測する方法は、これに限定されるものではない。例えば、国際公開第２００２／０７９６２９号パンフレットに記載された方法で自己着火時点を予測することとしてもよい。
【０１２５】
また、上述した実施の形態１では、１サイクルの全クランク角範囲（７２０degCA）に渡る模擬筒内圧履歴データを作成することとしているが、本発明における模擬筒内圧履歴データは、全クランク角範囲に渡るものに限定されない。例えば、模擬筒内圧履歴データは、燃焼期間中のみの筒内圧の履歴を表すものであってもよい。
なお、以上述べたような各変形例は、後述する各実施の形態においても同様に適用可能である。
【０１２６】
また、上述した実施の形態１においては、ノックサイクルの発生頻度が前記第１の発明における「運転安定性」に相当している。また、ECU５０が、上記ステップ１００〜１３２の処理を実行することにより前記第１の発明における「模擬筒内圧履歴データ生成手段」が、上記ステップ１３４〜１５０の処理を実行することにより前記第１の発明における「運転安定性予測手段」が、上記ステップ１５２〜１５８の処理を実行することにより前記第１の発明における「点火時期設定手段」が、それぞれ実現されている。
【０１２７】
また、上述した実施の形態１においては、ECU５０が、上記ステップ１２６および１２８の処理を実行することにより前記第２の発明における「燃焼モデル手段」が、上記ステップ１００〜１１８および１２２の処理を実行することにより前記第２の発明における「標準偏差予測手段」が、上記ステップ１２４の処理を実行することにより前記第２の発明における「燃焼パラメータ抽出手段」が、上記ステップ１２６および１２８の処理を実行することにより前記第２の発明における「筒内圧履歴算出手段」が、それぞれ実現されている。
【０１２８】
また、上述した実施の形態１においては、燃焼効率kが前記第３の発明における「燃焼効率に関する値」に、５０％燃焼点θ５０が前記第３の発明における「所定の燃焼割合となるクランク角に関する値」に、１０−７０％燃焼期間θ１０−７０が前記第３の発明における「第１の燃焼割合から第２の燃焼割合に変化するまでの燃焼期間に関する値」に、それぞれ相当している。
【０１２９】
また、上述した実施の形態１においては、機関回転数NE、負荷率KLおよびVVT進角量VTが前記第４の発明における「所定の運転パラメータ」に、ECU５０が前記第４の発明における「記憶した手段」に、それぞれ相当している。また、ECU５０が、上記ステップ１０８〜１１４の処理を実行することにより前記第４の発明における「平均値取得手段」が、上記ステップ１００の処理を実行することにより前記第４の発明における「運転パラメータ取得手段」が、それぞれ実現されている。
【０１３０】
また、上述した実施の形態１においては、ECU５０が、上記ステップ１１０の処理を実行することにより前記第５の発明における「燃焼パラメータ算出手段」が、上記ステップ１１４の処理を実行することにより前記第５の発明における「算術平均手段」が、それぞれ実現されている。
【０１３１】
また、上述した実施の形態１においては、ECU５０が、上記ステップ１４０〜１４６の処理を実行することにより前記第７の発明における「ノックサイクル判定手段」が、上記ステップ１４８および１５０の処理を実行することにより前記第７の発明における「ノッキング判定手段」が、それぞれ実現されている。
【０１３２】
また、上述した実施の形態１においては、ECU５０が、上記ステップ１４０の処理を実行することにより前記第８の発明における「自己着火予測手段」が、上記ステップ１４２〜１４６の処理を実行することにより前記第８の発明における「燃焼割合算出判定手段」が、それぞれ実現されている。
【０１３３】
実施の形態２．
次に、本発明の実施の形態２について説明するが、上述した実施の形態１との相違点を中心に説明し、同様の事項については、その説明を省略または簡略する。本実施形態のシステムは、図１に示すハードウェア構成を用いて、ECU５０に、後述する図１７および図１８に示すルーチンと、既述した図１０に示すルーチンとを実行させることにより実現することができる。
【０１３４】
［実施の形態２の特徴］
前述した実施の形態１では、通常運転状態において最適な点火時期を得るための制御に本発明を適用した場合について説明した。これに対し、本実施形態は、冷間始動時の触媒暖機運転において最適な点火時期を得るための制御を対象とする。
【０１３５】
内燃機関１０の排気通路３２に設けられた触媒３４は、活性温度まで加熱された状態において、良好な排気浄化作用を発揮する。触媒３４は、排気ガスの熱によって暖められる。排気エミッションを低減するためには、冷間始動後、触媒３４の温度をなるべく早期に活性温度まで昇温させる必要がある。このため、従来より、冷間始動後に点火時期を遅角させる触媒暖機制御が行われている。点火時期を遅角させると、燃焼ガスがピストンに対してする仕事が減少する分、排気温度が上昇する。このため、触媒３４を早期に暖機することができる。
【０１３６】
排気温度は、点火時期を遅角させるほど、高くなる。このため、触媒３４をより早期に暖機する観点からは、点火時期をより遅角側にするのが好ましい。しかしながら、点火時期を遅くするほど、燃焼が不安定となるので、燃焼変動が大きくなる。燃焼変動が大きくなると、トルク変動が大きくなり易い。従来の触媒暖機制御では、点火時期を遅角しすぎたために過大なトルク変動が生じて車両の快適性を損ねたり、あるいは、点火遅角幅が十分でないために触媒暖機に要する時間が長くなったりし易かった。
【０１３７】
以上のような問題に鑑み、本実施形態では、点火時期を現在値よりも遅角したと仮定した場合のトルク変動の大きさを予測することとした。そして、予測されたトルク変動の大きさが許容レベルのものであった場合に、点火時期を実際に遅角させることとした。
【０１３８】
点火時期を遅角した場合のトルク変動の大きさは、多数サイクル分の模擬筒内圧履歴データに基づいて予測することができる。すなわち、各々の模擬筒内圧履歴データについて図示トルクを算出することにより、サイクルごとのトルク変動を予測することができる。
【０１３９】
［実施の形態２における具体的処理］
図１７は、上記の機能を実現するために本実施形態においてECU５０が実行するルーチンのフローチャートである。なお、図１７において、図９に示すステップと同一のステップには、同一の符号を付してその説明を省略または簡略する。
【０１４０】
図１７に示すルーチンによれば、まず、現在の運転状態を表す各種のパラメータが計測され、保存される（ステップ１００）。次に、現在の運転状態が、触媒暖機制御中かつ定常状態か否かが判別される（ステップ１６０）。本実施形態では、内燃機関１０が触媒暖機制御中かつ定常状態にあることを前提とする。このため、ステップ１６０の条件が否定された場合には、データをリセットし（ステップ１０４）、その後、今回の処理サイクルを速やかに終了する。ステップ１０６〜１１６は、実施の形態１と同様の処理であるので、説明を省略する。
【０１４１】
次に、想定点火時期SAiおよび想定５０％燃焼点平均値BPiが、それぞれ、２degCA大きい値に更新される（ステップ１６２）。本実施形態における想定点火時期SAiは、トルク変動が許容範囲内に収まるか否かを予測する対象とする、仮想上の点火時期である。このステップ１６２の処理は、想定点火時期SAiを２degCA遅角させることに相当する。ステップ１６２の処理が初めて行われた場合、その処理後には、想定点火時期SAiおよび想定５０％燃焼点平均値BPiは、それぞれ、現在値より２degCA遅角側の値を表すこととなる。
【０１４２】
次のステップ１２０においては、実施の形態１と同様に、図１０に示すサブルーチンの処理が実行される。その処理により、想定点火時期SAiの下での、１００サイクル分の模擬筒内圧履歴データが作成される。
【０１４３】
次いで、その１００サイクル分の模擬筒内圧履歴データに基づいて、図示トルクの標準偏差が算出される（図１７中のステップ１６４）。図１８は、ステップ１６４の処理を行うためにECU５０が実行するサブルーチンのフローチャートである。
【０１４４】
図１８に示すサブルーチンでは、各々の模擬筒内圧履歴データについて図示トルクを算出するため、図示トルクを算出する処理（ステップ１６６）が１００回（１００サイクル分）繰り返し行われる。ステップ１６６において、図示トルクは、模擬筒内圧履歴データの筒内圧Pを７２０degCAの範囲に渡って積分（積算）することにより、算出することができる。以上の処理により、想定点火時期SAiの下での、サイクルごとの図示トルクが算出される。
【０１４５】
次に、１００サイクル分の図示トルクデータに基づいて、それらのバラツキ度合いを示す標準偏差σＴを算出する処理が行われる（ステップ１６８）。その後、図１８に示すサブルーチンの処理を終了する。
【０１４６】
次に、図示トルクの標準偏差σＴが許容値β以下であるか否かが判別される（図１７中のステップ１７０）。標準偏差σＴが大きいほど、サイクルごとの図示トルクのバラツキ度合いが大きく、よって、トルク変動が大きくなると考えられる。上記許容値βは、トルク変動の大きさが実際上問題とならない範囲に収まることを判定可能な値として、予め設定されている。
【０１４７】
上記ステップ１７０において、σＴ≦βの不成立が認められた場合には、想定点火時期SAiが実現されるとトルク変動が大きくなり過ぎる、との予測が成り立つ。この場合には、想定点火時期SAiが２degCA進角された値に更新される（ステップ１７２）。すなわち、想定点火時期SAiから２degCAが減算される。その後、その更新された想定点火時期SAiが実際の点火時期SAmとして設定される（ステップ１５４）。
【０１４８】
上記ステップ１５０の判別が行われたのが１回だけだった場合には、想定点火時期SAiは、上記ステップ１５２の処理により、現点火時期SAcに等しい値に戻っている。つまり、この場合には、実際の点火時期は変更されず、現点火時期SAcのままで内燃機関１０の運転が継続されることとなる。よって、この場合、実際にトルク変動が大きくなり過ぎることはない。
【０１４９】
一方、上記ステップ１７０において、σＴ≦βの成立が認められた場合には、想定点火時期SAiが実現されたとしても、トルク変動は許容できるレベルに収まる、との予測が成り立つ。しかし、触媒早期暖機の観点からは点火時期は遅いほどよい。そこで、この場合には、現時点の想定点火時期SAiよりもさらに点火時期を遅角させた場合にトルク変動が許容レベルに収まるかどうかを予測するべく、ステップ１６２以下の処理が再度実行される。具体的には、ステップ１６２において、想定点火時期SAiおよび想定５０％燃焼点平均値BPiがさらに２degCA遅角された値に更新される。そして、その更新された想定点火時期SAiの下でトルク変動が許容レベルに収まるか否かの予測が再度行われる（ステップ１２０〜１７０）。
【０１５０】
このように、ステップ１７０からステップ１６２へと戻るループが設けられていることにより、トルク変動が許容レベルに収まると予測される点火時期のうちでなるべく遅い点火時期を見つけ出すことができる。ステップ１７０からステップ１６２へと戻るループを通った場合には、ステップ１７０の判別が複数回行われることになる。複数回目のステップ１７０においてσＴ≦βの不成立が認められた場合、すなわちトルク変動が大きくなり過ぎるとの予測がなされた場合には、ステップ１７２において想定点火時期SAiが２degCA進角された値に更新され、その後、更新された想定点火時期SAiがステップ１５４において実際の点火時期SAmとして設定される。この場合に実際の点火時期SAmとして設定された値は、１回前のステップ１７０の処理においてトルク変動が許容レベルに収まると予測された点火時期の値に等しい。よって、この場合、実際上トルク変動が大きくなり過ぎることはない。
【０１５１】
以上説明したように、本実施形態によれば、点火時期を遅角させようとする際、点火時期遅角後のトルク変動の大きさが許容レベルを超えないかどうかを事前に予測することができる。そして、トルク変動の大きさが許容レベル以下となるとの予測が得られた場合に、点火時期を実際に遅角させる。このため、トルク変動による弊害を一瞬たりとも生じさせることなく、点火時期をなるべく遅い時期に制御することができる。よって、触媒早期暖機と、トルク変動の弊害防止とを両立することができる。
【０１５２】
ところで、上述した実施の形態２では、図示トルクのバラツキ度合いを示す値として、標準偏差σＴを算出しているが、図示トルクのバラツキ度合いを示す値は、これに限定されるものではない。例えば、図示トルクの分散を代わりに用いてもよい。
【０１５３】
なお、上述した実施の形態２においては、ECU５０が、上記ステップ１６６の処理を実行することにより前記第９の発明における「図示トルク算出手段」が、上記ステップ１６８の処理を実行することにより前記第９の発明における「バラツキ度合い算出手段」が、上記ステップ１７０の処理を実行することにより前記第９の発明における「トルク変動判定手段」が、それぞれ実現されている。
【０１５４】
また、上述した実施の形態２においては、トルク変動の大きさが前記第１の発明における「運転安定性」に相当している。また、ECU５０が、上記ステップ１３４〜１７０の処理を実行することにより前記第１の発明における「運転安定性予測手段」が、上記ステップ１７２および１５４の処理を実行することにより前記第１の発明における「点火時期設定手段」が、それぞれ実現されている。
【０１５５】
実施の形態３．
次に、本発明の実施の形態３について説明するが、上述した実施の形態１との相違点を中心に説明し、同様の事項については、その説明を省略または簡略する。本実施形態のシステムは、図２に示すハードウェア構成を用いて、ECU５０に、後述する図１９に示すルーチンを実行させることにより実現することができる。
【０１５６】
［実施の形態３の特徴］
本実施形態は、上述した実施の形態１の処理に加えて、上記（６）〜（９）式を実測データに基づいて較正する処理を更に行うことを特徴とする。既述した通り、実施の形態１では、燃焼パラメータk、θ５０、θ１０−７０の標準偏差σ１、σ２、σ３、およびθ５０とθ１０−７０との共分散σ２３を、上記（６）〜（９）式で表される回帰式に基づいて予測している。
【０１５７】
上記（６）〜（９）式中には、複数の係数が含まれている。本実施形態では、標準偏差σ１、σ２、σ３、および共分散σ２３の予測精度を更に高めるべく、これらの係数を、実測データに基づいて較正することとした。
【０１５８】
［実施の形態３における具体的処理］
図１９は、上記の機能を実現するために本実施形態においてECU５０が実行するルーチンのフローチャートである。このルーチンは、図９に示すルーチン中のステップ１１４の処理が行われるのに同期して、実行される。以下では、説明を簡単にするため、上記（６）式を較正する処理について代表的に説明する。
【０１５９】
上記（６）式の右辺には、４つの係数が含まれている。次式は、その４つの係数を未定係数Ｘ１、Ｘ２、Ｘ３、Ｘ４に置き換えた式である。
σ１／μ１＝Ｘ１NE＋Ｘ２KL＋Ｘ３KL２＋Ｘ４      ・・・（１１）
図１９に示すルーチンは、上記（１１）式中の未定係数Ｘ１、Ｘ２、Ｘ３、Ｘ４を実測データから逆算することにより、上記（６）式を較正するものである。
【０１６０】
図９中のステップ１１４の処理においては、既述した通り、５０サイクル分の燃焼パラメータk、θ５０、θ１０−７０の実測値データに基づいて、それらの平均値kc、BPc、BDcが算出される。このとき、図１９に示すルーチンが実行される。
【０１６１】
図１９に示すルーチンでは、まず、５０サイクル分の燃焼効率kの実測値から、その実標準偏差σｋが算出される（ステップ１７４）。次に、この実標準偏差σｋと、図９中のステップ１１４で算出された燃焼効率kの実測平均値kcと、図９中のステップ１００で計測された機関回転数NEおよび負荷率KLの値とが上記（１１）に代入される（ステップ１７６）。このとき、σｋが上記（１１）式中のσ１に、kcがμ１に、それぞれ代入される。これにより、未定係数Ｘ１、Ｘ２、Ｘ３、Ｘ４の関係を定めた較正式が一つ得られる。
【０１６２】
未定係数Ｘ１、Ｘ２、Ｘ３、Ｘ４の値を求めるには、上記較正式が未定係数Ｘ１、Ｘ２、Ｘ３、Ｘ４の数（ここでは４つ）だけ必要である。そこで、上記ステップ１７６の次には、上記較正式が４つ揃ったか否かが判別される。較正式が４つ揃っていない場合には、そのまま今回の処理サイクルを終了する。
【０１６３】
既述した通り、図９に示すルーチンは、定常運転状態が５０サイクル以上続いた場合のみ、ステップ１１４以降の処理が実行され、完遂される。ステップ１１４の処理が実行されると、図１９に示すルーチンが実行されて、上記ステップ１７６の処理により、較正式を１つ作成することができる。従って、図９に示すルーチンが４回完遂されると、上記較正式が４つ揃うことになる。
【０１６４】
上記較正式が４つ揃ったことが認められた場合には、次に、それら４つの較正式を連立方程式として解くことにより、未定係数Ｘ１、Ｘ２、Ｘ３、Ｘ４の値が算出される（ステップ１８０）。続いて、上記（６）式中の係数を、ステップ１８０で算出されたＸ１、Ｘ２、Ｘ３、Ｘ４の値に更新する（ステップ１８２）。この処理により、上記（６）式の回帰式が較正される。
【０１６５】
以上、上記（６）式を較正する場合について代表して説明したが、上記（７）〜（９）式についても、同様にして較正することができる。ただし、上記（７）式には係数が５個、上記（８）および（９）式には係数がそれぞれ６個ある。各式の較正に際しては、それら係数の数だけ較正式を揃えた上で、行えばよい。
【０１６６】
以上説明した通り、本実施形態では、上記（６）〜（９）式の回帰式を実測データに基づいて較正する処理が行われる。このため、燃焼パラメータk、θ５０、θ１０−７０の標準偏差σ１、σ２、σ３や、θ５０とθ１０−７０との共分散σ２３を予測する上で、内燃機関１０の機差バラツキや経年変化の影響をより正確に反映させることができる。従って、それらの予測精度を更に高めることができる。このため、点火時期変更後の燃焼変動の様子をより高い精度で予測することができる。ひいては、点火時期変更後にノッキングが起きるか否かといった運転安定性をより高い精度で予測することができる。
【０１６７】
なお、上述した実施の形態３においては、ECU５０が、上記ステップ１７４の処理を実行することにより前記第６の発明における「実標準偏差算出手段」が、上記ステップ１７６乃至１８２の処理を実行することにより前記第６の発明における「較正手段」が、それぞれ実現されている。
【０１６８】
また、上述した実施の形態３においては、上記（６）〜（９）式の回帰式を較正する処理を、実施の形態１に組み込むこととしているが、本発明はこれに限定されるものではない。すなわち、本実施形態において特有な上記の処理は、実施の形態２に組み込むこととしてもよい。
【０１６９】
実施の形態４．
次に、本発明の実施の形態４について説明するが、上述した実施の形態１との相違点を中心に説明し、同様の事項については、その説明を省略または簡略する。本実施形態のシステムは、図１に示すハードウェア構成を用いて、ECU５０に、図９に示すルーチンに代えて、後述する図２０に示すルーチンを実行させることにより実現することができる。
【０１７０】
［実施の形態４の特徴］
前述した実施の形態１においては、５０サイクルに渡り、燃焼パラメータk、θ５０、θ１０−７０の値を筒内圧センサ１８の出力に基づいてサイクルごとに実測していた。そして、その実測データに基づいて、現在の運転状態における燃焼パラメータk、θ５０、θ１０−７０の現実値の平均値kc、BPc、BDcを取得していた。これに対し、本実施形態では、数式モデルを用いることよって、燃焼パラメータk、θ５０、θ１０−７０の現実値の平均値kc、BPc、BDcを取得することとした。これにより、本実施形態では、筒内圧センサ１８を不要とすることができる。
【０１７１】
［実施の形態４における具体的処理］
図２０は、上記の機能を実現するために本実施形態においてECU５０が実行するルーチンのフローチャートである。なお、図２０において、図９に示すステップと同一のステップには、同一の符号を付してその説明を省略または簡略する。
【０１７２】
図２０に示すルーチンによれば、まず、現在の機関回転数NE、負荷率KL、およびVVT進角量VTが計測され、保存される（ステップ１００）。次に、現在の運転状態における燃焼パラメータk、θ５０、θ１０−７０の現実値の平均値kc、BPc、BDcが、次式の関数に基づいて算出される（ステップ１８４）。
Kc ＝ｆ１（NE，KL，VT，SAc，A/F）        ・・・（１２）
BPc＝ｆ２（NE，KL，VT，SAc，A/F）        ・・・（１３）
BDc＝ｆ３（NE，KL，VT，SAc，A/F）        ・・・（１４）
【０１７３】
上記（１２）〜（１４）式は、それぞれ、機関回転数NE、負荷率KL、VVT進角量VT、現点火時期SAc、および現空燃比A/Fを変数として、平均値kc、BPc、BDcを求めるための関数である。上記（１２）〜（１４）式は、設計段階で予め行われた実験に基づいて作成された多変数重回帰式である。すなわち、上記（１２）〜（１４）式は、NE，KL，VT，SAc，A/Fの異なる種々の運転条件下で実測された燃焼パラメータk、θ５０、θ１０−７０の平均値データから回帰することによって求められた、所定の関数である。ECU５０には、上記（１２）〜（１４）式が予め記憶されている。
【０１７４】
上記ステップ１８４においては、ステップ１００で計測されたNE、KL、VTの値が上記（１２）〜（１４）式に代入される。また、上記（１２）〜（１４）式中のSAcには、ECU５０が現在設定している点火時期の値が代入され、A/Fには、ECU５０が現在設定している目標空燃比の値、あるいは図示しない空燃比センサで検出される実空燃比の値が代入される。
【０１７５】
上記ステップ１８４の後、ステップ１１６以降は、実施の形態１と同様の処理が行われる。本実施形態によれば、筒内圧センサ１８が不要であるので、筒内圧センサ１８を備えないシステムにも適用することが可能である。また、本実施形態は、内燃機関１０が定常運転状態であることを前提としないので、過渡運転状態における点火時期の制御にも適用することができる。
【０１７６】
なお、上述した実施の形態４においては、ECU５０が、上記ステップ１８４の処理を実行することにより前記第４の発明における「平均値取得手段」が実現されている。
【０１７７】
また、上述した実施の形態４においては、燃焼パラメータk、θ５０、θ１０−７０の現実値の平均値kc、BPc、BDcを数式モデルを用いて算出する処理を、実施の形態１に組み込むこととしているが、本発明はこれに限定されるものではない。すなわち、本実施形態において特有な上記の処理は、実施の形態２に組み込むこととしてもよい。
【図面の簡単な説明】
【０１７８】
【図１】本発明の実施の形態１のシステム構成を説明するための図である。
【図２】クランク角に対する筒内圧履歴の実測データを示すグラフである。
【図３】定常運転状態における、サイクルごとの燃焼効率kの分布の実測データを示す図である。
【図４】定常運転状態における、サイクルごとの５０％燃焼点θ５０の分布の実測データを示す図である。
【図５】定常運転状態における、サイクルごとの１０−７０％燃焼期間θ１０−７０の分布の実測データを示す図である。
【図６】定常運転状態における５０％燃焼点θ５０および１０−７０％燃焼期間θ１０−７０の実測値の相関図である。
【図７】定常運転状態における燃焼効率kと５０％燃焼点θ５０の実測値の相関図である。
【図８】定常運転状態における燃焼効率kと１０−７０％燃焼期間θ１０−７０の実測値の相関図である。
【図９】本発明の実施の形態１において実行されるルーチンのフローチャートである。
【図１０】本発明の実施の形態１において実行されるルーチンのフローチャートである。
【図１１】想定点火時期SAiの下での５０％燃焼点θ５０のバラツキを示す正規分布である。
【図１２】想定点火時期SAiの下での燃焼効率kのバラツキを示す正規分布である。
【図１３】想定点火時期SAiの下での１０−７０％燃焼期間θ１０−７０のバラツキを示す正規分布である。
【図１４】模擬筒内圧履歴データの一例をグラフ化して示す図である。
【図１５】本発明の実施の形態１において実行されるルーチンのフローチャートである。
【図１６】クランク角に対する燃焼割合Xbの履歴を示すグラフである。
【図１７】本発明の実施の形態２において実行されるルーチンのフローチャートである。
【図１８】本発明の実施の形態２において実行されるルーチンのフローチャートである。
【図１９】本発明の実施の形態３において実行されるルーチンのフローチャートである。
【図２０】本発明の実施の形態４において実行されるルーチンのフローチャートである。
【符号の説明】
【０１７９】
１０  内燃機関
１２  クランク角センサ
１８  筒内圧センサ
２１  吸気圧センサ
３０  点火プラグ
３３  排気圧センサ
３４  触媒
３６  可変動弁機構
４０  カム角センサ
５０  ECU(Electronic Control Unit)
【図１】
2007056727.tif 000007
【図２】
2007056727.tif 000008
【図３】
2007056727.tif 000009
【図４】
2007056727.tif 000010
【図５】
2007056727.tif 000011
【図６】
2007056727.tif 000012
【図７】
2007056727.tif 000013
【図８】
2007056727.tif 000014
【図９】
2007056727.tif 000015
【図１０】
2007056727.tif 000016
【図１１】
2007056727.tif 000017
【図１２】
2007056727.tif 000018
【図１３】
2007056727.tif 000019
【図１４】
2007056727.tif 000020
【図１５】
2007056727.tif 000021
【図１６】
2007056727.tif 000022
【図１７】
2007056727.tif 000023
【図１８】
2007056727.tif 000024
【図１９】
2007056727.tif 000025
【図２０】
2007056727.tif 000026