(57)【要約】【課題】ガウス過程を用いたモデルによる出力を用いて内燃機関を適切に制御する制御装置を提供する。【解決手段】制御装置は、運転パラメータの値に基づいて制御対象である制御パラメータを制御する。制御装置は、運転パラメータの現在の値を取得し、取得された運転パラメータの現在の値に基づいて制御パラメータの値に対する出力パラメータの発生確率の確率分布をモデルを用いて算出し、算出した出力パラメータの発生確率の確率分布に基づいて、出力パラメータの値が目標値以上となる確率が目標確率に最も近づくように、制御パラメータの目標値を設定する、ように構成される。制御パラメータと運転パラメータと出力パラメータとは互いに異なるパラメータである。モデルは、運転パラメータの値及び制御パラメータの値が入力されると前力パラメータの発生確率の確率分布を出力するガウス過程を用いたモデルである。【選択図】図５
2019157652.tif 000002
【特許請求の範囲】
【請求項１】
内燃機関の運転に関する複数の運転パラメータの値に基づいて制御対象である制御パラメータを制御する、内燃機関の制御装置であって、
前記制御装置は、
前記運転パラメータの現在の値を取得し、
前記取得された運転パラメータの現在の値に基づいて前記制御パラメータの値に対する出力パラメータの発生確率の確率分布をモデルを用いて算出し、
前記算出した出力パラメータの発生確率の確率分布に基づいて、前記出力パラメータの値が基準値以上又は基準値以下となる確率が目標確率に最も近づくように、前記制御パラメータの目標値を設定する、ように構成され、
前記制御パラメータと前記運転パラメータと前記出力パラメータとは互いに異なるパラメータであり、
前記モデルは、前記運転パラメータの値及び前記制御パラメータの値が入力されると前記出力パラメータの発生確率の確率分布を出力するガウス過程を用いたモデルである、内燃機関の制御装置。
【請求項２】
前記内燃機関は燃焼室内の混合気に点火する点火プラグを備え、
前記制御パラメータは点火時期であり、前記出力パラメータはノック強度である、請求項１に記載の内燃機関の制御装置。
【請求項３】
内燃機関の運転に関する複数の運転パラメータの値に基づいて制御対象である制御パラメータを制御する、内燃機関の制御装置であって、
前記制御装置は、
前記運転パラメータの現在の値を取得し、
前記取得された運転パラメータの現在の値に基づいて前記制御パラメータの値に対する出力パラメータの発生確率の確率分布をモデルを用いて算出し、
前記算出した出力パラメータの発生確率の確率分布に基づいて、前記出力パラメータの値が目標値となる確率が最大となるように、前記制御パラメータの目標値を設定する、ように構成され、
前記制御パラメータと前記運転パラメータと前記出力パラメータとは互いに異なるパラメータであり、
前記モデルは、前記運転パラメータの値及び前記制御パラメータの値が入力されると前記出力パラメータの発生確率の確率分布を出力するガウス過程を用いたモデルである、内燃機関の制御装置。
【請求項４】
前記内燃機関は燃焼室に燃料を供給する燃料噴射弁を備え、
前記制御パラメータは前記燃料噴射弁からの燃料噴射量であり、前記出力パラメータは排気ガスの空燃比である、請求項３に記載の内燃機関の制御装置。
【請求項５】
前記制御装置は、前記内燃機関の運転中にオンボードで前記モデルを更新するように構成され、
前記モデルの更新は、前記内燃機関の運転中に取得された前記運転パラメータの値と前記制御パラメータの値とに基づいて、前記モデルを表すハイパーパラメータを更新することなく、逐次ガウス過程によって行われる、請求項１〜４のいずれか１項に記載の内燃機関の制御装置。
【請求項６】
前記モデルは、前記運転パラメータの値及び前記制御パラメータの値に応じて分散が変化する異分散ガウス過程を用いたモデルである、請求項１〜５のいずれか１項に記載の内燃機関の制御装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、内燃機関の制御装置に関する。
【背景技術】
【０００２】
従来から、内燃機関のデータに基づいて関数モデルを作成し、この関数モデルを用いて入力に対する出力の値を算出することが知られている。また、このような関数モデルを作成するにあたって、ガウス過程を用いることも知られている（例えば、特許文献１）。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２０１０−７０１７１号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
ところで、ガウス過程を用いたモデルでは、その出力は所定のパラメータの発生確率の確率分布の形となる。このため、ガウス過程を用いたモデルを内燃機関の制御に用いる場合であっても、そのモデルをそのまま内燃機関の制御に用いることはできない。したがって、斯かるモデルを用いて内燃機関の制御を行うには、斯かるモデルによって出力された確率分布に対する処理が必要である。
【０００５】
本発明は、上記課題に鑑みてなされたものであって、その目的は、ガウス過程を用いたモデルによる出力を用いて内燃機関を適切に制御する制御装置を提供することにある。
【課題を解決するための手段】
【０００６】
本発明は、上記課題を解決するためになされたものであり、その要旨は以下のとおりである。
【０００７】
（１）内燃機関の運転に関する複数の運転パラメータの値に基づいて制御対象である制御パラメータを制御する、内燃機関の制御装置であって、前記制御装置は、前記運転パラメータの現在の値を取得し、前記取得された運転パラメータの現在の値に基づいて前記制御パラメータの値に対する出力パラメータの発生確率の確率分布をモデルを用いて算出し、前記算出した出力パラメータの発生確率の確率分布に基づいて、前記出力パラメータの値が基準値以上となる確率が目標確率に最も近づくように、前記制御パラメータの目標値を設定する、ように構成され、前記制御パラメータと前記運転パラメータと前記出力パラメータとは互いに異なるパラメータであり、前記モデルは、前記運転パラメータの値及び制御パラメータの値が入力されると前記出力パラメータの発生確率の確率分布を出力するガウス過程を用いたモデルである、内燃機関の制御装置。
【０００８】
（２）前記内燃機関は燃焼室内の混合気に点火する点火プラグを備え、前記制御パラメータは点火時期であり、前記出力パラメータはノック強度である、上記（１）に記載の内燃機関の制御装置。
【０００９】
（３）内燃機関の運転に関する複数の運転パラメータの値に基づいて制御対象である制御パラメータを制御する、内燃機関の制御装置であって、前記制御装置は、前記運転パラメータの現在の値を取得し、前記取得された運転パラメータの現在の値に基づいて前記制御パラメータの値に対する出力パラメータの発生確率の確率分布をモデルを用いて算出し、前記算出した出力パラメータの発生確率の確率分布に基づいて、前記出力パラメータの値が目標値となる確率が最大となるように、前記制御パラメータの目標値を設定する、ように構成され、前記制御パラメータと前記運転パラメータと前記出力パラメータとは互いに異なるパラメータであり、前記モデルは、前記運転パラメータの値及び前記制御パラメータの値が入力されると前記出力パラメータの発生確率の確率分布を出力するガウス過程を用いたモデルである、内燃機関の制御装置。
【００１０】
（４）前記内燃機関は燃焼室に燃料を供給する燃料噴射弁を備え、前記制御パラメータは燃料噴射弁からの燃料噴射量であり、前記出力パラメータは排気ガスの空燃比である、上記（３）に記載の内燃機関の制御装置。
【００１１】
（５）前記制御装置は、内燃機関の運転中にオンボードで前記モデルを更新するように構成され、前記モデルの更新は、内燃機関の運転中に取得された前記運転パラメータの値と前記制御パラメータの値とに基づいて、前記モデルを表すハイパーパラメータを更新することなく、逐次ガウス過程によって行われる、上記（１）〜（４）のいずれか１つに記載の内燃機関の制御装置。
【００１２】
（６）前記モデルは、前記運転パラメータの値及び前記制御パラメータの値に応じて分散が変化する異分散ガウス過程を用いたモデルである、上記（１）〜（５）のいずれか１つに記載の内燃機関の制御装置。
【発明の効果】
【００１３】
本発明によれば、ガウス過程を用いたモデルによる出力を用いて内燃機関を適切に制御する制御装置が提供される。
【図面の簡単な説明】
【００１４】
【図１】図１は、制御装置が用いられる内燃機関を概略的に示す図である。
【図２】図２は、内燃機関の制御装置の機能ブロック図である。
【図３】図３は、ノック強度モデルによって算出されるノック強度の発生確率の確率分布を示している。
【図４】図４は、図３に示した確率分布のうち所定の点火時期におけるノック強度の対数値と発生確率との関係を示している。
【図５】図５は、基本点火時期算出部における基本点火時期の算出制御の制御ルーチンを示すフローチャートである。
【図６】図６は、内燃機関の制御装置の機能ブロック図である。
【図７】図７は、空燃比モデルによって算出される排気空燃比の発生確率の確率分布を示している。
【発明を実施するための形態】
【００１５】
以下、図面を参照して本発明の実施形態について詳細に説明する。なお、以下の説明では、同様な構成要素には同一の参照番号を付す。
【００１６】
なお、本明細書において、基本的に、小文字のみの文字列で表されるパラメータ（例えば、ｅｓａ）はスカラーを表し、Ｍを含まない大文字を含む文字列で表されるパラメータ（例えば、Ｘ）はベクトルを表し、Ｍを含む大文字を含む文字列で表されるパラメータ（例えば、ＭＸ）は行列を表す。
【００１７】
＜第一実施形態＞
≪内燃機関全体の説明≫
図１は、第一実施形態に係る制御装置が用いられる内燃機関を概略的に示す図である。図１に示したように、内燃機関１は、機関本体２、シリンダブロック３、シリンダブロック３内で往復動するピストン４、シリンダブロック３上に固定されたシリンダヘッド５、吸気弁６、吸気ポート７、排気弁８、排気ポート９を備える。ピストン４とシリンダヘッド５との間には燃焼室１０が形成される。吸気弁６は吸気ポート７を開閉し、排気弁８は排気ポート９を開閉する。また、機関本体２には、吸気弁６のバルブタイミングを制御する可変バルブタイミング機構２８が設けられる。なお、機関本体２には、排気弁８のバルブタイミングを制御する可変バルブタイミング機構が設けられてもよい。
【００１８】
図１に示したようにシリンダヘッド５の内壁面の中央部には点火プラグ１１が配置され、シリンダヘッド５の内壁面周辺部には燃料噴射弁１２が配置される。点火プラグ１１は、点火信号に応じて火花を発生させるように構成される。また、燃料噴射弁１２は、噴射信号に応じて、所定量の燃料を燃焼室１０内に噴射する。なお、燃料噴射弁１２は、吸気ポート７内に燃料を噴射するように配置されてもよい。
【００１９】
各気筒の吸気ポート７はそれぞれ対応する吸気枝管１３を介してサージタンク１４に連結され、サージタンク１４は吸気管１５を介してエアクリーナ１６に連結される。吸気ポート７、吸気枝管１３、サージタンク１４、吸気管１５は吸気通路を形成する。また、吸気管１５内にはスロットル弁駆動アクチュエータ１７によって駆動されるスロットル弁１８が配置される。
【００２０】
一方、各気筒の排気ポート９は排気マニホルド１９に連結され、排気マニホルド１９は排気浄化触媒２０を内蔵したケーシング２１に連結される。ケーシング２１は排気管２２に連結される。排気ポート９、排気マニホルド１９、ケーシング２１及び排気管２２は、排気通路を形成する。
【００２１】
排気マニホルド１９とサージタンク１４とはＥＧＲ管２４によって互いに連通せしめられる。ＥＧＲ管２４には、排気マニホルド１９からサージタンク１４へとＥＧＲ管２４内を流れるＥＧＲガスを冷却するＥＧＲクーラ２５が設けられる。加えて、ＥＧＲ管２４には、サージタンク１４へと供給されるＥＧＲガスの流量を制御するＥＧＲ制御弁２６が設けられる。ＥＧＲ管２４、ＥＧＲクーラ２５及びＥＧＲ制御弁２６は、排気ガスの一部を吸気通路に供給するＥＧＲ機構を構成する。
【００２２】
また、内燃機関１は、電子制御ユニット（ＥＣＵ）３１を備える。ＥＣＵ３１は、デジタルコンピュータからなり、双方向性バス３２を介して相互に接続されたＲＡＭ（ランダムアクセスメモリ）３３、ＲＯＭ（リードオンリメモリ）３４、ＣＰＵ（マイクロプロセッサ）３５、入力ポート３６及び出力ポート３７を具備する。
【００２３】
吸気管１５には、吸気管１５内を流れる空気流量を検出するためのエアフロメータ３９が設けられ、スロットル弁１８には、スロットル弁１８の開度を検出するスロットル開度センサ４０が設けられる。加えて、シリンダブロック３には、ノック強度を検出するためのノックセンサ４１が設けられ、排気マニホルド１９には、排気マニホルド１９内を流れる排気ガスの空燃比（以下、「排気空燃比」ともいう）を検出する空燃比センサ４２が設けられる。これらエアフロメータ３９、スロットル開度センサ４０、ノックセンサ４１及び空燃比センサ４２の出力は、対応するＡＤ変換器３８を介して入力ポート３６に入力される。なお、本実施形態では、ノックセンサ４１によってノック強度を検出しているが、燃焼室１０内の圧力を検出する筒内圧センサをシリンダヘッド５に設けると共にこの筒内圧センサの出力に基づいてノック強度を算出するようにしてもよい。
【００２４】
また、アクセルペダル４３にはアクセルペダル４３の踏込み量に比例した出力電圧を発生する負荷センサ４４が接続され、負荷センサ４４の出力電圧は対応するＡＤ変換器３８を介して入力ポート３６に入力される。クランク角センサ４５は例えばクランクシャフトが１５度回転する毎に出力パルスを発生し、この出力パルスが入力ポート３６に入力される。ＣＰＵ３５ではこのクランク角センサ４５の出力パルスから機関回転速度が計算される。
【００２５】
一方、出力ポート３７は対応する駆動回路４６を介して点火プラグ１１、燃料噴射弁１２及びスロットル弁駆動アクチュエータ１７に接続される。したがって、ＥＣＵ３１は、点火プラグ１１による点火時期、燃料噴射弁１２からの燃料噴射時期や燃料噴射量、スロットル弁１８の開度等を制御する制御装置として機能する。
【００２６】
≪点火時期制御≫
次に、図２を参照して、本実施形態における点火プラグ１１による燃焼室１０内の混合気への点火時期の目標値の算出方法について説明する。図２は、本実施形態に係るＥＣＵ３１の機能ブロック図である。
【００２７】
図２からわかるように、ＥＣＵ３１は、制御対象たる制御パラメータである点火時期の算出にあたり、大きく分けて二つの機能ブロックを有する。具体的には、ＥＣＵは、内燃機関１の運転に関する各種パラメータ（以下、「運転パラメータ」ともいう）の値に基づいてノック強度モデルを用いて基本点火時期を算出するモデル利用部Ａと、ノックセンサ４１によって検出されたノック強度に基づいて点火時期をフィードバック制御するＦＢ制御部Ｂとを備える。したがって、モデル利用部Ａは、各種運転パラメータの値に基づいて基本点火時期を算出するフィードフォワード制御を行っており、ＦＢ制御部Ｂは、検出されたノック強度に基づいて点火時期の目標値を算出するフィードバック制御を行っている。
【００２８】
モデル利用部Ａは、基本点火時期算出部Ａ１と、モデル更新部Ａ２とを備える。基本点火時期算出部Ａ１では、各種運転パラメータの現在の値に基づいて、基本点火時期ｅｓａｂａｓｅが算出される。基本点火時期算出部Ａ１に入力される運転パラメータとしては、具体的には、例えば、スロットル弁１８の開度θｔ、機関回転速度ｎｅ、燃焼室１０へ吸入された空気量（吸入空気量）ｍｃ、吸気弁６のバルブタイミングｉｖｔ、及び制御弁２６の開度ｄｅｇｒ等が含まれる（なお、本実施形態では、運転パラメータには点火時期、ノック強度は含まれない）。
【００２９】
また、基本点火時期算出部Ａ１には、モデル更新部Ａ２によって更新されたノック強度モデルを表す各パラメータ（以下、「モデルパラメータ」ともいう）の値がＲＡＭ３３から読み込まれる。ノック強度モデルは、上述した各種運転パラメータの値に対するノック強度の発生確率の確率分布を表すモデルである。換言すると、本実施形態におけるモデルは、運転パラメータの値に対する出力パラメータの発生確率の確率分布を表すモデルであるということができる。基本点火時期算出部Ａ１は、各種運転パラメータの現在の値に基づいて基本点火時期ｅｓａｂａｓｅを算出するにあたって、ノック強度モデルを用いる。基本点火時期算出部Ａ１における具体的な点火時期の算出方法については後述する。
【００３０】
モデル更新部Ａ２には、上述した内燃機関１の運転状態に関する各種運転パラメータに加えて、点火プラグ１１における点火時期ｅｓａと、その点火時期ｅｓａにおいて点火プラグ１１によって混合気へ点火されたときのノック強度ｋｉとが入力される。モデル更新部Ａ２では、このようにして入力された運転パラメータの値、点火時期ｅｓａ、ノック強度ｋｉを学習データとして、ノック強度モデルが更新される。モデル更新部Ａ２は、このようにして更新されたノック強度モデルを表す各モデルパラメータの値をＲＡＭ３３へ書き込む。ノック強度モデルの具体的な更新方法については後述する。
【００３１】
ＦＢ制御部Ｂは、点火時期算出部Ｂ１と、ノック発生判定部Ｂ２と、ＦＢ補正量算出部Ｂ３とを備える。点火時期算出部Ｂ１は、基本点火時期算出部Ａ１から出力された基本点火時期ｅｓａｂａｓｅに、ＦＢ補正量算出部によって算出されたＦＢ補正量Δｅｓａを加算することによって、点火時期ｅｓａを算出する（ｅｓａ＝ｅｓａｂａｓｅ＋Δｅｓａ）。算出された点火時期ｅｓａは制御信号として点火プラグ１１に送信され、点火プラグ１１はこの点火時期ｅｓａにおいて混合気へ点火する。
【００３２】
ノック発生判定部Ｂ２は、ノックセンサ４１によって検出されたノック強度ｋｉからノック基準強度ｋｉｒｅｆを減算したノック強度差Δｋｉを算出する（Δｋｉ＝ｋｉｒｅｆ−ｋｉ）。本実施形態では、ノック強度がノック基準強度ｋｉｒｅｆ以上である場合にノッキングが発生していると判定する。したがって、ノック発生判定部Ｂ２において算出されたノック強度差Δｋｉが負の値である場合にはノッキングが発生していると判定されていることを意味し、逆にノック強度差Δｋｉが正の値である場合にはノッキングは発生していないと判定されていることを意味する。
【００３３】
ＦＢ補正量算出部Ｂ３は、ノック強度差Δｋｉに基づいて、ＦＢ補正量Δｅｓａを算出する。具体的には、ＦＢ補正量Δｅｓａは、下記式（１）に基づいて算出される。
Δesak=Δesak-1+a・Δki      …（１）
上記式（１）において、Δｅｓａkは今回算出されるＦＢ補正量を、Δｅｓａk-1はＦＢ補正量算出部Ｂ３において前回算出されたＦＢ補正量をそれぞれ意味する。また、ａは予め定められた所定の正の定数である。式（１）からわかるように、ノッキングが発生していてノック強度差Δｋｉが負の値であるときにはＦＢ補正量Δｅｓａは小さくなる。逆に、ノッキングが発生していなくてノック強度差Δｋｉが正の値であるときにはＦＢ補正量Δｅｓａは大きくなる。
【００３４】
ＦＢ補正量算出部Ｂ３によって算出されたＦＢ補正量Δｅｓａは、上述したように点火時期算出部Ｂ１において、基本点火時期ｅｓａｂａｓｅに加算される。ここで、本実施形態における点火時期は、圧縮上死点からの進角度合い（°ＢＴＤＣ）で合わされるため、点火時期ｅｓａの値が大きくなるほど点火時期は進角される。ノッキングが発生している場合には、ＦＢ補正量Δｅｓａが小さくなるため、ＦＢ制御部Ｂにおけるフィードバック制御により点火時期が遅角されることになる。一方、ノッキングが発生していない場合には、ＦＢ補正量Δｅｓａが大きくなるため、ＦＢ制御部Ｂにおけるフィードバック制御により点火時期が進角されることになる。
【００３５】
なお、上述したＦＢ制御部Ｂにおけるフィードバック制御は一例であり、ＦＢ制御部Ｂでは、ＰＩＤ制御やＰＩ制御等、様々なフィードバック制御を用いることができる。また、ＥＣＵ３１での計算負荷を低減させるという観点からは、ＦＢ制御部Ｂにおけるフィードバック制御は行われなくてもよい。この場合には、モデル利用部Ａによるフィードフォワード制御のみが実行され、基本点火時期算出部Ａ１によって算出された基本点火時期ｅｓａｂａｓｅが制御信号として点火プラグ１１へ送信される。
【００３６】
≪基本点火時期の算出≫
次に、図３及び図４を参照して、基本点火時期算出部Ａ１における基本点火時期の算出方法について説明する。図３は、ノック強度モデルによって算出されるノック強度の発生確率の確率分布を示している。図４は、図３に示した確率分布のうち所定の点火時期におけるノック強度の対数値と発生確率との関係を示している。
【００３７】
ここで、ノック強度は、内燃機関１の運転状態が同一でも必ずしも同一の値になるとは限らず、確率的に発生することが知られている。特に、ノック強度の発生確率の確率分布は、対数正規分布で近似される。したがって、内燃機関１の運転状態をＸとし、各ノック強度の発生確率をｙとすると、ノック強度モデルにおけるＸとｙとの関係は下記式（２）によって表される。なお、Ｘは、点火時期ｅｓａ並びにスロットル弁の開度θｔ及び機関回転速度ｎｅ等の各種運転パラメータをパラメータとして有するベクトルを示している（Ｘ＝［ｅｓａ，θｔ，ｎｅ，．．．］）。
y|X〜N(f(X),σ2)      …（２）
【００３８】
上記式（２）においてｆ（Ｘ）は平均値を表し、σ2は分散を表しており、またＮ（μ，σ2）は平均値がμであって分散がσ2の正規分布を表している。したがって、上記式（２）は、ノック強度モデルでは、ノック強度の発生確率ｙが、平均値がｆ（Ｘ）であって分散がσ（Ｘ）の正規分布に従うことを表している。
【００３９】
点火時期以外の内燃機関１の運転状態が固定されている場合、ノック強度モデルによって算出される各ノック強度の発生確率ｙは、点火時期に応じて変化する。この様子を、図３に示す。図３は、点火時期以外の内燃機関１の運転状態が固定されている状態で、ノック強度モデルにおいて算出される点火時期と、ノック強度の対数値と、各ノック強度の発生確率との関係の一例を示す。
【００４０】
図４は、図３に示した確率分布のうち或る点火時期（例えば、１０°ＢＴＤＣ）におけるノック強度の対数値とその発生確率との関係を示す図である。図４は、点火時期も固定されている場合の確率分布を示していることから、図４は任意一つの運転状態Ｘにおけるノック強度の発生確率ｙの確率分布を示しているといえる。図４に示したように、本実施形態では、或る運転状態Ｘにおける各ノック強度の発生確率ｙは正規分布に従うものとして近似されている。
【００４１】
ここで、本実施形態では、ノック強度ｋｉが所定の基準値ｋｉｒｅｆ以上であるときに内燃機関１にノッキングが発生していると判定される。したがって、或る運転状態Ｘにおいてノック強度ｋｉが基準値ｋｉｒｅｆ未満である領域における発生確率ｙの積分値（図４中のα）は、その運転状態Ｘにおいてノッキングが発生しない確率ｐｎｔを表している。一方、或る運転状態Ｘにおいてノック強度ｋｉが基準値ｋｉｒｅｆ以上である領域における発生確率ｙの積分値（図４中のβ）は、その運転状態Ｘにおいてノッキングが発生する確率（以下、「ノック発生確率」ともいう）ｐｋｎを表している。
【００４２】
そして、本実施形態では、ノック発生確率ｐｋｎが目標ノック発生確率ｐｔｒｇとなる点火時期が基準点火時期ｅｓａｂａｓｅとして算出される。ノック発生確率ｐｋｎが目標ノック発生確率ｐｔｒｇとなる点火時期は基本的に一意に決まるが、複数の点火時期においてノック発生確率ｐｋｎが目標ノック発生確率ｐｔｒｇとなるような場合には、これら複数の点火時期のうち最も進角側の点火時期が基準点火時期ｅｓａｂａｓｅとして算出される。
【００４３】
すなわち、本実施形態では、出力パラメータ（ノック強度）の発生確率の確率分布に基づいて、出力パラメータの値が基準値以上となる確率（ノック発生確率）が目標確率（目標ノック発生確率）に最も近づくように、制御パラメータ（点火時期）の目標値を設定しているといえる。
【００４４】
ただし、上述したノック強度ｋｉの算出は、例えば点火時期を所定角度（例えば、０．１°）毎にずらして入力して行われる。したがって、ノック発生確率ｐｋｎは所定角度の点火時期毎にしか算出することができず、よって点火時期に対するノック発生確率ｐｋｎを連続的に算出することはできない。このため、必ずしも目標ノック発生確率ｐｔｒｇに一致する点火時期を算出することはできない。そこで、本実施形態では、離散的に入力した点火時期のうちノック発生確率ｐｋｎが目標ノック発生確率ｐｔｒｇに最も近い値となる点火時期を基準点火時期ｅｓａｂａｓｅとして算出するようにしてもよい。或いは、離散的に入力した点火時期のうちノック発生確率ｐｋｎが目標ノック発生確率ｐｔｒｇ以下であって目標ノック発生確率ｐｔｒｇに最も近い値となる点火時期を基準点火時期ｅｓａｂａｓｅとして算出するようにしてもよい。
【００４５】
なお、図３からわかるように、ノック強度の平均値（各点火時期において発生確率がピークをとるノック強度）は、基本的に、点火時期が進角されるほど、すなわち図３における点火時期の角度が大きくなるほど、大きくなる。したがって、基本的には、ノック発生確率ｐｋｎも、点火時期が進角されるほど大きくなる。このため、ノック発生確率ｐｋｎが目標ノック発生確率ｐｔｒｇとなる点火時期は上述したように一意に決まることになる。
【００４６】
また、ノック発生確率はｐｋｎは点火時期が進角されるほど大きくなる。したがって、ノック発生確率ｐｋｎが目標ノック発生確率ｐｔｒｇ又はそれに最も近い値となるように基本点火時期を決めることは、実質的にノック発生確率ｐｋｎが目標ノック発生確率ｐｔｒｇ以下である点火時期のうち最も進角側の点火時期を基準点火時期ｅｓａｂａｓｅに設定しているといえる。
【００４７】
さらに、上記実施形態では、ノック発生確率ｐｋｎが目標ノック発生確率ｐｔｒｇに最も近づくように点火時期の目標値を設定している。しかしながら、ノッキングが発生しない確率ｐｎｔ、すなわち出力パラメータ（ノック強度）の値が基準値以下となる確率が目標確率に最も近づくように、制御パラメータ（点火時期）の目標値を設定してもよい。
【００４８】
ここで、点火時期が遅角されると、基本的に、燃焼室１０内における混合気の燃焼に伴う熱発生が全体的に遅角側にシフトされることになると共に、混合気の燃焼が緩慢になる。このため、点火時期が遅角されると、基本的には、熱効率が悪化し、よって燃費や機関出力の悪化を招く。したがって、本実施形態では、ノック発生確率ｐｋｎを目標ノック発生確率ｐｔｒｇ以下に維持しつつ、できる限り燃費や機関出力が高くなるように点火時期を設定しているといえる。
【００４９】
図５は、基本点火時期算出部Ａ１における基本点火時期の算出制御の制御ルーチンを示すフローチャートである。図示した制御ルーチンは一定時間間隔毎に実行される。
【００５０】
図５に示したように、まず、ステップＳ１１において各種運転パラメータの現在の値が取得される。具体的には、斯かる運転パラメータには、例えば、スロットル弁１８の開度θｔ、機関回転速度ｎｅ、吸入空気量ｍｃ、吸気弁６のバルブタイミングｉｖｔ、及びＥＧＲ制御弁２６の開度ｄｅｇｒ等の少なくともいずれか一つが含まれる。
【００５１】
スロットル弁１８の開度θｔはスロットル開度センサ４０によって検出され、機関回転速度ｎｅはクランク角センサ４５の出力に基づいて算出され、吸入空気量ｍｃはエアフロメータ３９の出力に基づいて算出される。吸気弁６のバルブタイミングｉｖｔは、吸気弁のバルブタイミングを検出するセンサ（図示せず）によって検出されてもよいし、可変バルブタイミング機構２８への制御信号に基づいて算出されてもよい。また、ＥＧＲ制御弁２６の開度ｄｅｇｒは、ＥＧＲ制御弁２６の開度を検出するセンサ（図示せず）によって検出されてもよいし、ＥＧＲ制御弁２６への制御信号に基づいて算出されてもよい。
【００５２】
次いで、ステップＳ１２では、モデル更新部Ａ２によって算出されたノック強度モデルを表す各モデルパラメータがＲＡＭ３３から取得される。モデル更新部Ａ２では、ノック強度モデルを表す各種モデルパラメータの一部の値が学習により更新されるため、ステップＳ１２では具体的には更新された各種パラメータの値が取得されることになる。
【００５３】
次いで、ステップＳ１３では、ステップＳ１１で取得された内燃機関１の運転状態に関するパラメータの現在の値に基づいて、ステップＳ１２で取得されたノック強度モデルを用いて、図３に示したような点火時期に対するノック強度の発生確率の確率分布が算出される。
【００５４】
次いで、ステップＳ１４では、ステップＳ１３で算出された点火時期に対するノック強度の発生確率の確率分布に基づいて、各点火時期におけるノック発生確率ｐｋｎが算出される。そして、算出されたノック発生確率ｐｋｎが目標ノック発生確率ｐｔｒｇに最も近い値となる点火時期が基本点火時期ｅｓａｂａｓｅとして算出される。
【００５５】
≪ノック強度モデル≫
次に、ノック強度モデルの作成方法及び更新方法について説明する。上述したように、ノック現象は同じ運転状態においても確率的に発生し、特にノック強度の対数値の確率分布は正規分布でよく近似されることが知られている。そこで、本実施形態では、ノック強度モデルとして、ガウス過程（Gaussian Process、ＧＰ）モデルを用いる。このようにノック強度モデルとしてＧＰモデルを用いることにより、少量の学習データからモデルを構築することができるようになる。
【００５６】
≪ノック強度モデルの作成≫
まず、ノック強度モデルの作成方法について説明する。ノック強度モデルの作成は、ＧＰモデルであるノック強度モデルを表す各モデルパラメータの値を設定することを意味する。ノック強度モデルの作成は、例えば、内燃機関１を搭載した車両の出荷前に行われる。ノック強度モデルの作成にあたっては、複数の組の学習データが利用される。
【００５７】
ここで、ノック強度モデルを作成するのにｎ組の学習データを用いることを考えたとき、ノック強度モデルへの入力の学習データをＭＸ＝［Ｘ1，Ｘ2，．．．，Ｘn］、ノック強度モデルからの出力の学習データをＹ＝［ｙ1，ｙ2，．．．，ｙn］T、学習データをＤ＝（ＭＸ、Ｙ）とする。入力の各学習データＸnは、内燃機関の運転状態を表す各種運転パラメータ（スロットル弁の開度θｔn、機関回転速度ｎｅn等）及び点火時期ｅｓａnを含む。また、出力の各学習データは、ノックセンサ４１によって検出されるノック強度ｋｉを含む。
【００５８】
任意のカーネル関数をｋ（・，・）として表すと（・にはベクトルや行列が代入される）、ＧＰの事前分布をｆ（Ｘ）〜ＧＰ（０，ｋ（Ｘ、Ｘ’））、観測ノイズをσ2としたとき、すなわちｙ｜Ｘ〜Ｎ（ｆ（Ｘ）、σ2）としたとき、予測分布は下記式（３）のように表される。
yt|X*,Θ,D〜N(μf*,σ2f*)      …（３）
【００５９】
ここで、式（３）におけるＸ*はノック強度モデルによって実際にノック強度の発生確率の確率分布を算出する際の任意の入力データを表し、ｙ*はこの入力データに対応する出力データ（すなわち、ノック強度の発生確率の確率分布）を表す。また、Θは、ノック強度モデルを表すモデルパラメータを表している。
【００６０】
加えて、式（３）における平均値μf*及び分散σ2f*は、それぞれ下記式（４）及び（５）で表される。
μf*=k(X*,MX)(MK+σ2MI)-1Y      …（４）
σf*=k(X*,X*)-k(X*,MX)(MK+σ2MI)-1k(MX,X*)+σ2MI      …（５）
【００６１】
式（４）及び式（５）において行列ＭＩは単位行列を表している。また、行列ＭＫ＝ｋ（ＭＸ，ＭＸ）であり、行列Ｘが与えられたときにカーネル関数を表す行列は下記式（６）及び式（７）で定義される。
k(MX,MX)=[kij], kij=k(Xi,Xj)      …（６）
k(X*,MX)=[k(X*,X1),...,k(X*,Xn)]
=k(MX,X*)T                …（７）
【００６２】
ＧＰは主にカーネル関数ｋ（・，・）によってその性質が定まる。本実施形態ではカーネルとして、ガウスカーネルを拡張したＡＲＤカーネルを用いる。したがって、本実施形態におけるカーネル関数は下記式（８）のように表される。
【数１】
2019157652.tif 000003
【００６３】
式（８）において、MΛ＝ｄｉａｇ（ｌ12，ｌ22，．．．，ｌd2）であり、ベクトルＸの各要素間の関係又はベクトルＸの各要素のノック強度に対する影響度を特徴付けるスケールである。また、λ2は、潜在関数の分散を表すパラメータである。これら全てのパラメータΘ＝［ｌ12，ｌ22，．．．，ｌd2，λ2，σ2］はハイパーパラメータと呼ばれ、ノック強度モデルを表すモデルパラメータの一部を構成する。
【００６４】
これらパラメータΘについては、例えば、ＥＭ法を用いて、下記式（９）に示した周辺尤度最大化により最適な値が求められる。また、式（９）におけるｌｏｇ（ｐ（Ｙ｜ＭＸ，Θ））は下記式（１０）によって表される。
【数２】
2019157652.tif 000004
【００６５】
上述した式（３）〜式（１０）を用いることによって、ｎ組の学習データ（ＭＸ及びＹ）から、ノック強度モデルを作成することができる。具体的には、式（３）〜式（１０）により、ｎ組の学習データに基づいて、ノック強度モデルにおける各モデルパラメータの値が算出される。
【００６６】
このようにして作成されたノック強度モデルでは、入力データＸ*が入力されると、上記式（４）を用いて平均値μf*を算出することができ、上記式（５）を用いて分散σf2*を算出することができる。すなわち、各種運転パラメータ及び点火時期ｅｓａが入力されると、その運転状態におけるノック強度の発生確率の確率分布を平均値がμf*で分散がσf2*である図４に示したような正規分布として算出することができる。
【００６７】
なお、上記実施形態では、カーネルとして、ＡＲＤカーネルを用いている。ＡＲＤカーネルは学習モデルが連続且つ滑らかである場合に良い性能を示すことから、本実施形態でも比較的高い精度でノック強度の発生確率の確率分布を算出することができる。しかしながら、カーネルとしてはガウスカーネルやSpectral Mixture（ＳＭ）カーネル、ニューラルネットワークカーネル等、様々なカーネルを用いることができる。
【００６８】
この場合、ガウスカーネルを用いると、学習計算に伴う計算負荷を低減することができるもののＡＲＤカーネルに比べて表現力が低下する。また、ＳＭカーネルでは、学習モデルが複数の高周波数成分を持つ場合に良い性能を示す可能性があるものの、学習計算に伴う計算負荷が増大する。
【００６９】
≪ノック強度モデルの更新≫
ところで、内燃機関１の各運転状態に対するノック強度は必ずしも一定ではなく、内燃機関１の運転時間が長くなるにつれて変化していく。これは、例えば、燃焼室１０内にカーボン等が付着して燃焼室１０内での混合気の燃焼状態が変化することによって生じる。このため、ノック強度モデルによるノック強度の発生確率の確率分布の推定精度を高く維持するためには、ノック強度モデルを所定の間隔で更新することが必要になる。
【００７０】
ところが、上述したノック強度モデルの作成方法と同様な手法によりノック強度モデルの更新を行うと、ノック強度モデルを更新する毎に、上述した全ての計算を行うことが必要になる。ノック強度モデルの更新は基本的に内燃機関の運転中にオンボードで行う必要があるため、このようにしてノック強度モデルの更新を行うと、ＥＣＵ３１での計算負荷が極めて高いものとなってしまう。
【００７１】
ここでノック境界近傍のみを考える局所モデルでは、観測ノイズがスカラー値であるＧＰモデルで十分な近似が可能である。すなわち、上述したようにしてノック強度モデルを作成することによりノック境界を大まかに推定することができれば、詳細なノック境界の変化は計算負荷の低いＧＰによって求めることができる。そこで、本実施形態では、ノック強度モデルの更新を行うにあたっては逐次ガウス過程（Recursive Gaussian Process、ＲＧＰ）が用いられる。
【００７２】
具体的には、以下の手法によりノック強度モデルの更新が行われる。まず、ＧＰと同様に、学習データＤ＝（ＭＸ，Ｙ）とし、Ｆ＝ｆ（ＭＸ）として定義する。ここで、Ｆは、初期分布ｐ（Ｆ）＝Ｎ（Ｆ｜μ0f，ＭＣ0f）を持つＧＰモデルであると仮定する（なお、μ0fはベクトルを表す。以下、μについて同様）。ここで、μ0fは上述したノック強度モデルの作成時に算出され、ＭＣ0f＝ｋ（ＭＸ，ＭＸ）である。ＧＰモデルでは、一度定義された事前分布は変わらないが、ＲＧＰモデルにおける事前分布は、新たな入力の学習データＸkとそれに対応する出力の学習データｙkが与えられると、この学習データによってオンボードで更新される。ノック強度モデルの更新は、カルマンフィルターの更新則と同様に、以下の計算式により行われる。
【００７３】
まず、ステップｋ−１の事後分布を用いて、ステップｋでの予測分布ｐ（Ｙk｜Ｙ1:k-1）＝Ｎ（Ｙk｜μkp，ＭＣkp＋σ2ＭＩ）が下記式（１１）及び（１２）により計算される。また、式（１１）及び式（１２）におけるＭＪk及び式（１２）におけるＭＢkはそれぞれ下記式（１３）及び式（１４）により算出される。
【数３】
2019157652.tif 000005
【００７４】
次いで、新たな出力の学習データｙkを用いて、ｆの事後分布は下記式（１５）及び（１６）により算出される。また、式（１５）及び（１６）におけるＭＧkは下記式（１７）により算出される。
【数４】
2019157652.tif 000006
【００７５】
上記式（１１）〜（１７）からわかるように、本実施形態におけるノック強度モデルの更新においては、ハイパーパラメータの更新は行われない。加えて、本実施形態におけるノック強度モデルの更新においては、新たに追加された学習データのみについて計算が行われ、過去の学習データについては計算が行われない。このため、ノック強度モデルの更新に伴うＥＣＵ３１の計算負荷を低減することができる。
【００７６】
＜第二実施形態＞
次に、第二実施形態に係る制御装置について説明する。第二実施形態に係る制御装置における構成及び制御は基本的に第一実施形態に係る制御装置における構成及び制御と同様である。したがって、以下では、第一実施形態に係る制御装置とは異なる部分を中心に説明する。
【００７７】
ところで、上記第一実施形態では、式（３）の予測分布を求めるにあたり、ＧＰの観測ノイズσ2を入力値に異存しないスカラー値としている。このため、上記第一実施形態のノック強度モデルでは、観測ノイズσ2が入力値に依存した分散を持つモデルとしては表現されていない。しかしながら、ノック強度の発生確率の確率分布における分散は運転パラメータに応じて変化すると考えられることから、上記第一実施形態におけるノック強度モデルでは必ずしも高い精度でノック強度の発生確率の確率分布を推定できない可能性がある。
【００７８】
そこで、本実施形態では、入力値、すなわち運転パラメータの値に依存する分散を表現することができるように、下記式（１８）に示すノイズモデルを追加した異分散ガウス過程（Heteroscedastic Gaussian Process、ＨＧＰ）を考える。
【数５】
2019157652.tif 000007
なお、式（１８）及び式（１９）において、σn2（Ｘ）は運転パラメータの値に依存する分散を示している。
【００７９】
ｖも正規分布に従うことから、下記式（２０）のように表すことができる。また、本実施形態ではｖについてもＡＲＤカーネルを用いることからカーネル関数は下記式（２１）のように表される。
【数６】
2019157652.tif 000008
ここで、MΛn＝ｄｉａｇ（ｍ12，ｍ22，．．．，ｍd2）であり、ベクトルＸの各要素間の関係又はベクトルＸの各要素の分散に対する影響度を特徴付けるスケールである。また、λn2は、潜在関数の分散を表すパラメータである。これらパラメータもハイパーパラメータであり、したがって本実施形態におけるノック強度モデルで用いられるハイパーパラメータはΦ＝［ｌ12，ｌ22，．．．，ｌd2，λ2，ｍ12，ｍ22，．．．，ｍd2，λn2，σ2］のように表される。
【００８０】
ここでの学習は、例えば、期待値伝搬法やＥＭ法を応用することによって行われる。この学習では、ｖの事後分布ｐ（ｖ｜Ｄ）をガウス分布ｑ（ｖ｜Ｄ）に近似し、周辺尤度最大化によりハイパーパラメータΦの最適な値が求められる。入力データＸ*を与えたときの出力データの予測値ｙ*は、ガウス分布で近似したｑ（ｖ*｜Ｘ*，Ｄ）＝Ｎ（μv*，σv2*）を用いて下記式（２２）により算出される。式（２２）における平均値μ*及び分散σ2*は、それぞれ下記式（２３）及び（２４）で表される。
【数７】
2019157652.tif 000009
【００８１】
このようにして作成されたノック強度モデルでも、入力データＸ*が入力されると、上記式（４）及び式（２３）を用いて平均値μ*を算出することができ、上記式（５）及び式（２４）を用いて分散σ2*を算出することができる。すなわち、各種運転パラメータ及び点火時期ｅｓａが入力されると、その運転状態におけるノック強度の発生確率の確率分布を平均値がμf*で分散がσf2*である図４に示したような正規分布として算出することができる。
【００８２】
本実施形態によれば、ノック強度モデルにおける分散を入力データに応じて変化するものとしてノック強度の発生確率の確率分布を算出している。したがって、より高い精度で、ノック強度の発生確率の確率分布を求めることができる。
【００８３】
加えて、本実施形態の手法に基づいてノック強度モデルを作成した場合であっても、ノック強度モデルの更新には逐次ガウス過程を用いることができる。この場合には、上記第一実施形態と同様にハイパーパラメータΦの更新は行われず、よって本実施形態においてもノック強度モデルの更新に伴うＥＣＵ３１の計算負荷を低減することができる。
【００８４】
＜第三実施形態＞
次に、第三実施形態に係る制御装置について説明する。第三実施形態に係る制御装置における構成及び制御は基本的に第一実施形態及び第二実施形態に係る制御装置における構成及び制御と同様である。したがって、以下では、第一実施形態及び第二実施形態に係る制御装置とは異なる部分を中心に説明する。
【００８５】
ところで、上記第一実施形態では、ノック強度に基づいた点火時期の制御を行っている。これに対して本実施形態では、排気空燃比に基づいた燃料噴射弁１２からの燃料噴射量の制御が行われる。
【００８６】
≪燃料噴射量制御≫
図６を参照して、本実施形態における燃料噴射弁１２からの燃料噴射量の算出方法について説明する。図６は、本実施形態に係るＥＣＵ３１の機能ブロック図である。
【００８７】
図６からわかるように、ＥＣＵ３１は、制御対象たる制御パラメータである燃料噴射量の算出にあたり、大きく分けて二つの機能ブロックを有する。具体的には、ＥＣＵは、運転パラメータの値に基づいて空燃比モデルを用いて基本燃料噴射量を算出するモデル利用部Ａと、空燃比センサ４２に基づいて燃料噴射量をフィードバック制御するＦＢ制御部Ｂとを備える。したがって、モデル利用部Ａは、各種運転パラメータの値に基づいて基本噴射量を算出するフィードフォワード制御を行っており、ＦＢ制御部Ｂは、検出された排気空燃比に基づいて燃料噴射量を算出するフィードバック制御を行っている。
【００８８】
モデル利用部Ａは、基本噴射量算出部Ａ１と、モデル更新部Ａ２とを備える。基本噴射量算出部Ａ１では、各種運転パラメータの現在の値に基づいて、基本燃料噴射量ｑｂａｓｅが算出される。なお、本実施形態では、運転パラメータには燃料噴射量、排気空燃比は含まれないものとする。
【００８９】
また、基本点火時期算出部Ａ１には、モデル更新部Ａ２によって更新された空燃比モデルを表す各モデルパラメータの値がＲＡＭ３３から読み込まれる。空燃比モデルは、上述した各種運転パラメータの値に対する排気空燃比の発生確率の確率分布を表すモデルである。基本噴射量算出部Ａ１は、各種運転パラメータの現在の値に基づいて基本噴射量ｑｂａｓｅを算出するにあたって、空燃比モデルを用いる。基本噴射量算出部Ａ１における具体的な点火時期の算出方法については後述する。
【００９０】
モデル更新部Ａ２には、上述した内燃機関１の運転状態に関する各種運転パラメータに加えて、燃料噴射量ｑと、その燃料噴射量ｑの燃料噴射が行われたときの空燃比ａｆとが入力される。モデル更新部Ａ２では、このようにして入力された運転パラメータの値、燃料噴射量ｑ、空燃比ａｆを学習データとして、空燃比モデルが更新される。モデル更新部Ａ２は、このようにして更新された空燃比モデルを表す各モデルパラメータの値をＲＡＭ３３へ書き込む。空燃比モデルの具体的な更新方法については後述する。
【００９１】
ＦＢ制御部Ｂは、噴射量算出部Ｂ１と、空燃比差算出部Ｂ２と、ＦＢ補正量算出部Ｂ３とを備える。点火時期算出部Ｂ１は、基本噴射量算出部Ａ１から出力された基本燃料噴射量ｑｂａｓｅに、ＦＢ補正量算出部によって算出されたＦＢ補正量Δｑを加算することによって、燃料噴射量ｑを算出する（ｑ＝ｑｂａｓｅ＋Δｑ）。算出された燃料噴射量ｑは制御信号として燃料噴射弁１２に送信され、燃料噴射弁１２はこの燃料噴射量ｑにて燃料噴射を行う。
【００９２】
空燃比差算出部Ｂ２は、空燃比センサ４２によって検出された排気空燃比ａｆから目標空燃比ａｆｔｇｔを減算した空燃比差Δａｆを算出する（Δａｆ＝ａｆ−ａｆｔｇｔ）。ＦＢ補正量算出部Ｂ３は、空燃比差Δａｆに基づいて、ＦＢ補正量Δｑを算出する。具体的には、ＦＢ補正量Δｑは、下記式（２５）に基づいて算出される。
Δqk=Δqk-1+b・Δq      …（２５）
上記式（２５）において、Δｑkは今回算出されるＦＢ補正量を、Δｑk-1はＦＢ補正量算出部Ｂ３において前回算出されたＦＢ補正量をそれぞれ意味する。また、ｂは予め定められた所定の正の定数である。
【００９３】
なお、本実施形態においても、ＦＢ制御部Ｂでは様々なフィードバック制御を用いることができる。また、ＦＢ制御部Ｂにおけるフィードバック制御は行われなくてもよい。
【００９４】
≪基本燃料噴射量の算出≫
次に、図７を参照して、基本噴射量算出部Ａ１における基本燃料噴射量の算出方法について説明する。図７は、空燃比モデルによって算出される排気空燃比の発生確率の確率分布を示している。
【００９５】
ここで、排気空燃比は、内燃機関１の運転状態が同一でも必ずしも同一の値になるとは限らず、確率的に発生する。特に、排気空燃比の発生確率の確率分布は、対数正規分布で近似される。したがって、内燃機関１の運転状態をＸとし、各空燃比の発生確率をｙとすると、空燃比モデルにおけるＸとｙとの関係は、上記式（２）と同様に、下記式（２６）によって表される。なお、Ｘは、燃料噴射量ｑ並びにスロットル弁の開度θｔ及び機関回転速度ｎｅ等の各種運転パラメータをパラメータとして有するベクトルを示している（Ｘ＝［ｑ，θｔ，ｎｅ，．．．］）。
y|X〜N(f(X),σ2)      …（２６）
【００９６】
燃料噴射量以外の内燃機関１の運転状態が固定されている場合、空燃比モデルによって算出される各空燃比の発生確率ｙは、燃料噴射量に応じて変化する。この様子を、図７に示す。図７は、燃料噴射量以外の内燃機関１の運転状態が固定されている状態で、空燃比モデルにおいて算出される燃料噴射量と、排気空燃比の対数値と、各空燃比の発生確率との関係の一例を示す。
【００９７】
図７に示したような排気空燃比の発生確率の確率分布を求めることができれば、空燃比が或る特定の目標空燃比となる確率が最も高くなる燃料噴射量を算出することができる。そこで、本実施形態では、出力パラメータである空燃比が目標空燃比となる確率が最大となるような燃料噴射量が基本燃料噴射量ｑｂａｓｅとして算出される。すなわち、本実施形態では、出力パラメータ（空燃比）の発生確率の確率分布に基づいて、出力パラメータの値が目標値（目標空燃比）となる確率が最大となるように、制御パラメータ（燃料噴射量）の目標値を設定しているといえる。
【００９８】
なお、本実施形態における空燃比モデルも、上記第一実施形態及び第二実施形態におけるノック強度モデルと同様に、ガウス過程や異分散ガウス過程を用いて作成される。加えて、本実施形態における空燃比モデルも、上記第一実施形態及び第二実施形態におけるノック強度モデルと同様に、逐次ガウス過程を用いて更新される。
【００９９】
なお、本実施形態では、排気空燃比に基づいた燃料噴射量の制御を行っているが、本実施形態における制御と同様な制御は他の制御にも適用することができる。例えば、本実施形態における制御と同様な制御を、燃焼室１０へのＥＧＲガスの供給量に基づいたＥＧＲ弁の開度の制御や、燃焼室１０へのＥＧＲガスの供給量に基づいた吸気弁６のバルブタイミング又は排気弁８のバルブタイミングの制御に用いることも可能である。
【符号の説明】
【０１００】
１    内燃機関
６    吸気弁
８    排気弁
１１  点火プラグ
１２  燃料噴射弁
３１  ＥＣＵ
３９  エアフロメータ
４０  スロットル開度センサ
４１  ノックセンサ
【図１】
2019157652.tif 000010
【図２】
2019157652.tif 000011
【図３】
2019157652.tif 000012
【図４】
2019157652.tif 000013
【図５】
2019157652.tif 000014
【図６】
2019157652.tif 000015
【図７】
2019157652.tif 000016