"内燃機関の気筒に吸入される新気量である気筒吸入空気量を制御する吸入空気量制御手段を備える内燃機関の制御装置において、前記機関のスロットル弁を通過する新気の流量であるスロットル弁通過空気流量を算出または検出するスロットル弁通過空気流量取得手段と、前記機関の吸気圧を検出する吸気圧検出手段と、前記機関に吸入される空気の温度である吸気温を検出する吸気温検出手段と、前記吸気圧及び吸気温に基づいて理論気筒吸入空気量を算出する理論気筒吸入空気量算出手段と、前記気筒吸入空気量の過去値及び前記理論気筒吸入空気量を用いて、前記機関の体積効率を算出する体積効率算出手段と、前記機関の吸気管をモデル化した吸気管モデル式に、前記スロットル弁通過空気流量及び体積効率を適用して、前記気筒吸入空気量を算出する気筒吸入空気量算出手段と、該気筒吸入空気量を用いて該気筒吸入空気量の予測値である予測気筒吸入空気量を算出する気筒吸入空気量予測手段と、目標気筒吸入空気量を算出する目標気筒吸入空気量算出手段と、前記吸気管モデル式の逆モデル式に前記目標気筒吸入空気量及び体積効率を適用して、目標スロットル弁通過空気流量を算出する目標スロットル弁通過空気流量算出手段とを備えることを特徴とする内燃機関の制御装置。前記体積効率算出手段は、前記気筒吸入空気量算出手段により算出された気筒吸入空気量を用いて、前記体積効率を少なくとも1回更新することを特徴とする請求項1に記載の内燃機関の制御装置。前記体積効率算出手段は、前記気筒の容積と前記吸気管の容積との比、及び前記スロットル弁通過空気流量をさらに用いて前記体積効率を算出することを特徴とする請求項1に記載の内燃機関の制御装置。前記スロットル弁の開度を検出するスロットル弁開度検出手段と、前記目標スロットル弁通過空気流量に応じて前記スロットル弁の目標開度を算出する目標開度算出手段とを備え、前記スロットル弁通過空気流量取得手段は、前記スロットル弁の開度と該スロットル弁を通過する空気の流量との関係をモデル化した弁通過空気流量モデル式に前記スロットル弁開度を適用して、前記スロットル弁通過空気流量を算出し、前記目標開度算出手段は、前記弁通過空気流量モデル式の逆モデル式に前記目標スロットル弁通過空気流量を適用して前記目標開度を算出することを特徴とする請求項1から3の何れか1項に記載の内燃機関の制御装置。前記スロットル弁通過空気流量を検出する空気流量検出手段と、前記スロットル弁通過空気流量の検出値に基づいて前記弁通過空気流量モデル式を補正する弁通過空気流量モデル補正手段とを備え、前記スロットル弁通過空気流量取得手段は、補正された弁通過空気流量モデル式を用いて前記スロットル弁通過空気流量の算出を行い、前記目標開度算出手段は、補正された弁通過空気流量モデル式の逆モデル式を用いて前記目標開度の算出を行うことを特徴とする請求項4に記載の内燃機関の制御装置。前記スロットル弁の開度の開き側及び/または閉じ側の限界値を算出する限界値算出手段と、前記スロットル弁の開度を前記限界値に設定したときのスロットル弁通過空気流量に相当する限界スロットル弁通過空気流量を算出する限界スロットル弁通過空気流量算出手段と、前記限界スロットル弁通過空気流量に対応する限界気筒吸入空気量を算出する限界気筒吸入空気量算出手段と、前記目標気筒吸入空気量が前記限界気筒吸入空気量で規定される制限範囲内の値をとるように制限処理を行う制限処理手段とを備え、前記目標スロットル弁通過空気流量算出手段は、前記制限処理後の目標気筒吸入空気量を用いて前記目標スロットル弁通過空気流量の算出を行うことを特徴とする請求項1から5の何れか1項に記載の内燃機関の制御装置。機関の出力トルクを要求トルクに精度良く一致させるためには、気筒吸入空気量に適した量の燃料を供給することが必要であり、実際の気筒吸入空気量を、要求トルクに応じた目標気筒吸入空気量に一致させる吸気量制御と、実際の気筒吸入空気量に応じた燃料供給量制御とを行う必要がある。特許文献1には、燃料供給量制御については詳細な記載が無いが、気筒吸入空気量の予測値を算出する点が示されていないので、算出される気筒吸入空気量が燃料供給量制御に適用されると推測される。しかしながら、そのような燃料供給量制御では、要求トルクが変化する過渡状態では燃料量の過不足が生じ、制御精度が低下する。そこで特許文献1に示された装置に、特許文献2に示された気筒吸入空気量の予測値算出手法を適用することにより、改善は可能であるが、予測値は予測誤差を当然に含むため、気筒吸入空気量の予測値を、前回予測値を用いて算出する特許文献1の手法では、予測値の算出精度の点で改善の余地があった。本発明は、上述した点に着目してなされたものであり、内燃機関の気筒吸入空気量の予測値をより高い精度で算出しつつ、実際の気筒吸入空気量を目標値に精度良く制御することができる内燃機関の制御装置を提供することを目的とする。上記目的を達成するため請求項1に記載の発明は、内燃機関の気筒に吸入される新気量である気筒吸入空気量GAIRCYLNを制御する吸入空気量制御手段を備える内燃機関の制御装置において、前記機関のスロットル弁3を通過する新気の流量であるスロットル弁通過空気流量HGAIRTH,GAIRTHを算出または検出するスロットル弁通過空気流量取得手段と、前記機関の吸気圧PBAを検出する吸気圧検出手段と、前記機関に吸入される空気の温度である吸気温TAを検出する吸気温検出手段と、前記吸気圧PBA及び吸気温TAに基づいて理論気筒吸入空気量GAIRSTDを算出する理論気筒吸入空気量算出手段と、前記気筒吸入空気量の過去値GAIRCYLN及び前記理論気筒吸入空気量GAIRSTDを用いて、前記機関の体積効率ηvを算出する体積効率算出手段と、前記機関の吸気管をモデル化した吸気管モデル式に、前記スロットル弁通過空気流量HGAIRTH及び体積効率ηvを適用して、前記気筒吸入空気量GAIRCYLNを算出する気筒吸入空気量算出手段と、該気筒吸入空気量GAIRCYLNを用いて該気筒吸入空気量の予測値である予測気筒吸入空気量GAIRCYLPを算出する気筒吸入空気量予測手段と、目標気筒吸入空気量GAIRCYLCMDを算出する目標気筒吸入空気量算出手段と、前記吸気管モデル式の逆モデル式に前記目標気筒吸入空気量GAIRCYLCMD及び体積効率ηvを適用して、目標スロットル弁通過空気流量GAIRTHCMDを算出する目標スロットル弁通過空気流量算出手段とを備えることを特徴とする。請求項2に記載の発明は、請求項1に記載の内燃機関の制御装置において、前記体積効率算出手段は、前記気筒吸入空気量算出手段により算出された気筒吸入空気量GAIRCYLNを用いて、前記体積効率ηvを少なくとも1回更新することを特徴とする。請求項3に記載の発明は、請求項1に記載の内燃機関の制御装置において、前記体積効率算出手段は、前記気筒の容積Vcylと前記吸気管の容積Vinとの比、及び前記スロットル弁通過空気流量HGAIRTHをさらに用いて前記体積効率ηvを算出することを特徴とする。請求項4に記載の発明は、請求項1から3の何れか1項に記載の内燃機関の制御装置において、前記スロットル弁の開度THを検出するスロットル弁開度検出手段と、前記目標スロットル弁通過空気流量GAIRTHCMDに応じて前記スロットル弁の目標開度THCMDを算出する目標開度算出手段とを備え、前記スロットル弁通過空気流量取得手段は、前記スロットル弁の開度と該スロットル弁を通過する空気の流量との関係をモデル化した弁通過空気流量モデル式に前記スロットル弁開度THを適用して、前記スロットル弁通過空気流量HGAIRTHを算出し、前記目標開度算出手段は、前記弁通過空気流量モデル式の逆モデル式に前記目標スロットル弁通過空気流量GAIRTHCMDを適用して前記目標開度THCMDを算出することを特徴とする。請求項5に記載の発明は、請求項4に記載の内燃機関の制御装置において、前記スロットル弁通過空気流量GAIRTHを検出する空気流量検出手段と、前記スロットル弁通過空気流量の検出値GAIRTHに基づいて前記弁通過空気流量モデル式を補正する弁通過空気流量モデル式補正手段とを備え、前記スロットル弁通過空気流量取得手段は、補正された弁通過空気流量モデル式を用いて前記スロットル弁通過空気流量HGAIRTHの算出を行い、前記目標開度算出手段は、補正された弁通過空気流量モデル式の逆モデル式を用いて前記目標開度THCMDの算出を行うことを特徴とする。請求項6に記載の発明は、請求項1から5の何れか1項に記載の内燃機関の制御装置において、前記スロットル弁の開度THの開き側及び/または閉じ側の限界値THOFH,THOFLを算出する限界値算出手段と、前記スロットル弁の開度を前記限界値THOFH,THOFLに設定したときのスロットル弁通過空気流量に相当する限界スロットル弁通過空気流量GAIRTHLMH,GAIRTHLMLを算出する限界スロットル弁通過空気流量算出手段と、前記限界スロットル弁通過空気流量GAIRTHLMH,GAIRTHLMLに対応する限界気筒吸入空気量GAIRCYLLMH,GAIRCYLLMLを算出する限界気筒吸入空気量算出手段と、前記目標気筒吸入空気量GAIRCYLCMDが前記限界気筒吸入空気量GAIRCYLLMH,GAIRCYLLMLで規定される制限範囲内の値をとるように制限処理を行う制限処理手段とを備え、前記目標スロットル弁通過空気流量算出手段は、前記制限処理後の目標気筒吸入空気量GAIRCYLCMDを用いて前記目標スロットル弁通過空気流量GAIRTHCMDの算出を行うことを特徴とする。請求項1に記載の発明によれば、スロットル弁通過空気流量が算出または検出され、検出される吸気圧及び吸気温に基づいて理論気筒吸入空気量が算出され、気筒吸入空気量の過去値及び理論気筒吸入空気量を用いて、機関の体積効率が算出される。また機関の吸気管をモデル化した吸気管モデル式に、スロットル弁通過空気流量及び体積効率を適用して、気筒吸入空気量が算出され、その気筒吸入空気量を用いて予測気筒吸入空気量が算出される。すなわち、算出された気筒吸入空気量の最新値を用いて予測気筒吸入空気量が算出されるので、前回予測値を用いる先行例と比較して、より高い精度で予測気筒吸入空気量を算出することができる。また吸気管モデル式の逆モデル式に目標気筒吸入空気量及び体積効率を適用して、目標スロットル弁通過空気流量が算出されるので、気筒吸入空気量の算出に適用される体積効率をそのまま使用して、目標スロットル弁通過空気流量を精度良く算出することができる。したがって、算出された目標スロットル弁通過空気流量に応じてスロットル弁の目標開度を算出し、スロットル弁開度をその目標開度に制御することによって、実際の気筒吸入空気量を目標気筒吸入空気量に精度良く制御することができる。請求項2に記載の発明によれば、算出された気筒吸入空気量を用いて、体積効率が少なくとも1回更新されるので、体積効率の算出精度を高めることができる。請求項3に記載の発明によれば、気筒容積と吸気管容積との比、及びスロットル弁通過空気流量をさらに用いて体積効率が算出される。気筒容積と吸気管容積との比及びスロットル弁通過空気流量をさらに用いることによって、繰り返し更新演算を行うことなく、体積効率の算出精度を高めることができる。請求項4に記載の発明によれば、弁通過空気流量モデル式に検出スロットル弁開度を適用して、スロットル弁通過空気流量が算出されるとともに、弁通過空気流量モデル式の逆モデル式に目標スロットル弁通過空気流量を適用してスロットル弁の目標開度が算出される。これにより、検出スロットル弁開度に応じたスロットル弁通過空気流量の算出、及び目標スロットル弁通過空気流量に応じた目標開度の算出をともに精度良く行うことができる。請求項5に記載の発明によれば、スロットル弁通過空気流量の検出値に基づいて弁通過空気流量モデル式が補正され、補正された弁通過空気流量モデル式を用いてスロットル弁通過空気流量の算出が行われるとともに、補正された弁通過空気流量モデル式の逆モデル式を用いて目標開度の算出が行われる。したがって、スロットル弁の特性ばらつきや経時変化、有効吸気管断面積の減少、運転環境の変化、あるいはセンサの温度特性などに起因する弁通過空気流量モデルのモデル化誤差を抑制し、良好な制御性能を得ることができる。請求項6に記載の発明によれば、スロットル弁開度の開き側及び/または閉じ側の限界値が算出され、スロットル弁の開度を算出される限界値に設定したときのスロットル弁通過空気流量に相当する限界スロットル弁通過空気流量が算出され、さらに限界スロットル弁通過空気流量に対応する限界気筒吸入空気量が算出され、目標気筒吸入空気量が限界気筒吸入空気量で規定される制限範囲内の値をとるように制限処理が行われ、制限処理後の目標気筒吸入空気量を用いて目標スロットル弁通過空気流量の算出が行われる。したがって、実現不可能なスロットル弁開度に対応した目標スロットル弁通過空気流量が算出されることを確実に防止し、特に要求トルクが変化する過渡状態において目標開度を適切に設定することができる。以下本発明の実施の形態を図面を参照して説明する。[第1の実施形態]図1は、本発明の一実施形態にかかる内燃機関とその制御装置の構成を示す図であり、図1において、例えば4気筒を有する内燃機関以下単にエンジンという1は、吸気弁の作動位相を連続的に変更する弁作動特性可変機構40を備えている。エンジン1の吸気管2の途中にはスロットル弁3が配置されている。また、スロットル弁3にはその開度THを検出するスロットル弁開度センサ4が連結されており、スロットル弁開度THに応じた電気信号を出力して電子制御ユニット以下ECUという5に供給する。スロットル弁3には、スロットル弁3を駆動するアクチュエータ7が接続されており、アクチュエータ7は、ECU5によりその作動が制御される。吸気管2には、スロットル弁3を介してエンジン1に吸入される空気新気の流量である吸入空気流量GAIRを検出する吸入空気流量センサ13が設けられ、さらにスロットル弁3の上流側に吸気温TAを検出する吸気温センサ9が設けられている。これらのセンサ13及び9の検出信号は、ECU5に供給される。燃料噴射弁6はエンジン1とスロットル弁3との間かつ吸気管2の図示しない吸気弁の少し上流側に各気筒毎に設けられており、各噴射弁は図示しない燃料ポンプに接続されていると共にECU5に電気的に接続されて当該ECU5からの信号により燃料噴射弁6の開弁時間が制御される。エンジン1の各気筒の点火プラグ12は、ECU5に接続されており、ECU5は点火プラグ12に点火信号を供給し、点火時期制御を行う。スロットル弁3の下流には吸気圧PBAを検出する吸気圧センサ8が取付けられている。またエンジン1の本体には、エンジン冷却水温TWを検出するエンジン冷却水温センサ10が取り付けられている。これらのセンサ8及び10の検出信号は、ECU5に供給される。ECU5には、エンジン1のクランク軸図示せずの回転角度を検出するクランク角度位置センサ11が接続されており、クランク軸の回転角度に応じた信号がECU5に供給される。クランク角度位置センサ11は、エンジン1の特定の気筒の所定クランク角度位置でパルス以下CYLパルスというを出力する気筒判別センサ、各気筒の吸入行程開始時の上死点TDCに関し所定クランク角度前のクランク角度位置で4気筒エンジンではクランク角180度毎にTDCパルスを出力するTDCセンサ及びTDCパルスより短い一定クランク角周期例えば6度周期で1パルス以下CRKパルスというを発生するCRKセンサから成り、CYLパルス、TDCパルス及びCRKパルスがECU5に供給される。これらのパルスは、燃料噴射時期、点火時期等の各種タイミング制御、エンジン回転数エンジン回転速度NEの検出に使用される。ECU5には、エンジン1によって駆動される車両のアクセルペダルの踏み込み量以下アクセルペダル操作量というAPを検出するアクセルセンサ31、エンジン1により駆動される車両の走行速度車速VPを検出する車速センサ32、及び大気圧PAを検出する大気圧センサ33が接続されている。これらのセンサの検出信号は、ECU5に供給される。またエンジン1は排気還流機構図示せずを備えており、エンジン1の排気が吸気管2のスロットル弁3の下流側に還流される。ECU5は各種センサからの入力信号波形を整形し、電圧レベルを所定レベルに修正し、アナログ信号値をデジタル信号値に変換する等の機能を有する入力回路、中央演算処理ユニット以下CPUという、CPUで実行される演算プログラム及び演算結果等を記憶する記憶回路のほか、アクチュエータ7、燃料噴射弁6、点火プラグ12、弁作動特性可変機構40に駆動信号を供給する出力回路等から構成される。ECU5のCPUは、上記センサの検出信号に応じて、スロットル弁3の開度制御、エンジン1に供給する燃料量燃料噴射弁6の開弁時間の制御、点火時期制御、及び吸気弁の作動位相制御を行う。またECU5のCPUは、スロットル弁3を通過する空気流量の推定値以下推定スロットル弁通過空気流量というHGAIRTHを算出するとともに、推定スロットル弁通過空気流量HGAIRTH、吸気圧PBA、及び吸気温TAに基づいてエンジン1の気筒に吸入される新気量である気筒吸入空気量GAIRCYLN[g/TDC]1TDC期間、すなわち4気筒エンジンであればエンジン1のクランク軸が180度回転するのに要する時間当たりの空気量を算出する。さらにECU5のCPUは、スロットル弁3を通過する空気流量の1行程期間後の予測値以下予測スロットル弁通過空気流量というHGAIRTHPを算出するとともに、予測スロットル弁通過空気流量HGAIRTHP及び気筒吸入空気量GAIRCYLNを用いて、気筒吸入空気量GAIRCYLNの1行程期間後の予測値以下予測気筒吸入空気量というGAIRCYLPを算出する。算出した予測気筒吸入空気量GAIRCYLPは、燃料供給量や点火時期の制御に適用される。また、ECU5のCPUは、エンジン1の要求出力に応じた目標気筒吸入空気量GAIRCYLCMDを算出し、実際の気筒吸入空気量を目標気筒吸入空気量GAIRCYLCMDに制御するためのスロットル弁3の目標開度THCMDを算出し、検出されるスロットル弁開度THが目標開度THCMDと一致するように、アクチュエータ7の駆動制御を行う。図2は、予測気筒吸入空気量GAIRCYLP及び目標開度THCMDを算出するモジュールの構成を示すブロック図である。図2に示す各ブロックの機能は、ECU5のCPUによる演算処理により実現される。図2に示す演算モジュールは、スロットル弁通過空気流量算出部51と、予測気筒吸入空気量算出部52と、モデル補正係数算出部53と、目標気筒吸入空気量算出部54と、リミット処理部55と、目標スロットル弁通過空気流量算出部56と、目標開度算出部57とを備えている。スロットル弁通過空気流量算出部51は、検出されるスロットル弁開度THを、弁通過空気流量モデル式に適用して、推定吸入空気流量HGAIR[g/sec]を算出し、後述する長期モデル補正係数KMDLLを用いて推定吸入空気流量HGAIRを補正することにより、第1補正推定吸入空気流量HGAIRLを算出する。弁通過空気流量モデル式は、スロットル弁開度THと、スロットル弁3を通過する空気の流量との関係をモデル化したものである。スロットル弁通過空気流量算出部51は、さらに後述する短期モデル補正係数KMDLSを用いて第1補正推定吸入空気流量HGAIRLを補正することにより、第2補正推定吸入空気流量HGAIRLSを算出し、第2補正推定吸入空気流量HGAIRLS[g/sec]の単位変換を行って、推定スロットル弁通過空気流量HGAIRTH[g/TDC]を算出する。予測気筒吸入空気量算出部52は、吸気管2より具体的には吸気管2のスロットル弁下流側をモデル化した吸気管モデル式に、推定スロットル弁通過空気流量HGAIRTHを適用して、気筒吸入空気量GAIRCYLNを算出するとともに、推定スロットル弁通過空気流量HGAIRTHを用いて予測スロットル弁通過空気流量HGAIRTHPを算出し、気筒吸入空気量GAIRCYLN及び予測スロットル弁通過空気流量HGAIRTHPを用いて、予測気筒吸入空気量GAIRCYLPを算出する。また、予測気筒吸入空気量算出部52は、吸気管モデル式のモデルパラメータの1つであるエンジン1の体積効率ηvを算出し、気筒吸入空気量GAIRCYLNの算出に適用する。算出される体積効率ηvは、目標スロットル弁通過空気流量算出部56における演算にも適用される。モデル補正係数算出部53は、検出される吸入空気流量GAIR、推定吸入空気流量HGAIR、及び第1補正推定吸入空気流量HGAIRLに基づいて、弁通過空気流量モデルのモデル化誤差を補正するための長期モデル補正係数KMDLL及び短期モデル補正係数KMDLSを算出する。長期モデル補正係数KMDLLは、スロットル弁3の特性ばらつき及び経時変化に起因する、弁通過空気流量モデルのモデル化誤差を補正するための補正係数であり、短期モデル補正係数KMDLSは、環境例えば気温や大気圧などの変化に起因する、弁通過空気流量モデルのモデル化誤差を補正するための補正係数である。目標気筒吸入空気量算出部54は、アクセルペダル操作量AP及びエンジン回転数NEに応じて、目標気筒吸入空気量GAIRCYLCMDを算出する。リミット処理部55は、スロットル弁開度の1制御周期で変更可能なスロットル弁開度変化量最大変化量DTHLMTと、検出されるスロットル弁開度TH現在値とに応じて、開き側限界値THOFH及び閉じ側限界値THOFLを算出し、スロットル弁開度を開き側限界値THOFH及び閉じ側限界値THOFLに設定したときのスロットル弁通過空気流量に相当する上限スロットル弁通過空気流量GAIRTHLMH及び下限スロットル弁通過空気流量GAIRTHLMLを算出する。この演算には、弁通過空気流量モデル式が適用され、長期モデル補正係数KMDLL及び短期モデル補正係数KMDLSが使用される。リミット処理部55は、さらに上限スロットル弁通過空気流量GAIRTHLMH及び下限スロットル弁通過空気流量GAIRTHLMLに対応する上限気筒吸入空気量GAIRCYLHMH及び下限気筒吸入空気量GAIRCYLHMLを算出し、目標気筒吸入空気量GAIRCYLCMDが上限気筒吸入空気量GAIRCYLHMH及び下限気筒吸入空気量GAIRCYLHMLで定義される制限範囲内の値をとるようにリミット処理を行う。目標スロットル弁通過空気流量算出部56は、リミット処理後の目標気筒吸入空気量GAIRCYLCMDを吸気管モデル式の逆モデル式に適用し、目標スロットル弁通過空気流量GAIRTHCMDを算出する。目標開度算出部57は、目標スロットル弁通過空気流量GAIRTHCMDを弁通過空気流量モデル式の逆モデル式に適用し、目標開度THCMDを算出する。この演算には、長期モデル補正係数KMDLL及び短期モデル補正係数KMDLSが適用される。次に、スロットル弁通過空気流量算出部51における推定スロットル弁通過空気流量HGAIRTHの算出手法と、予測気筒吸入空気量算出部52における予測スロットル弁通過空気流量HGAIRTHP、気筒吸入空気量GAIRCYLN、及び予測気筒吸入空気量GAIRCYLPの算出手法とを詳細に説明する。図3はエンジン1を模式的に示す図であり、スロットル弁3,吸気弁21、排気弁22、及び気筒1aが示されている。スロットル弁3の開度THと、スロットル弁通過空気流量吸入空気流量GAIRの推定値HGAIR[g/sec]との関係は、弁通過空気流量モデルを定義する下記式1で表すことができる。式1のKCは流量の単位を[g/sec]とするための変換定数であり、KTHTHはスロットル弁開度THに応じて算出される開口面積流量関数であり、ΨRPは、スロットル弁3の上流側圧力である大気圧PAと、下流側圧力である吸気圧PBAとの比率RP=PBA/PAに応じて算出される圧力比流量関数であり、Rは気体定数である。開口面積流量関数KTHTHの値は、予め実験的に求められた図7aに示すKTHテーブルを用いて算出される。また圧力比流量関数Ψは、下記式2で与えられる。式2のκは空気の比熱比である。ただし、空気流速が音速を超えると、圧力比流量関数Ψは圧力比に拘わらず極大値をとるので、実際の演算処理では、圧力比流量関数ΨRPの値も予め設定されたΨRPテーブル図7bを用いて算出される。本実施形態では、弁通過空気流量モデルのモデル化誤差を補正する長期モデル補正係数KMDLL及び短期モデル補正係数KMDLSを下記式3に適用し、第2補正推定吸入空気流量HGAIRLSを算出し、式3により算出される第2補正推定吸入空気流量HGAIRLSを、エンジン回転数NEを用いて1TDC期間当たりの流量に変換することにより、推定スロットル弁通過空気流量HGAIRTH[g/TDC]が算出される。HGAIRLS=HGAIR×KMDLL×KMDLS3一方吸気管2のスロットル弁下流側部分2a内の空気量の変化量DGAIRINは、下記式11で与えられる。式11のVinはスロットル弁下流側部分2aの吸気管容積、TAKは絶対温度に変換した吸気温TA、Rは気体定数、DPBAは吸気圧PBAの変化量PBA-PBA-GAIRCYLN12一方、気筒吸入空気量GAIRCYLNは、下記式13で与えられる。式13のVcylは気筒容積であり、ηvは体積効率である。GAIRCYLN=Vcyl×ηv×PBA/R×TAK13式13を用いると、吸気圧変化量DPBAは、下記式14で与えられる。式14で与えられるDPBA及び式12の関係を式11に適用することにより、下記の式15吸気管モデルの定義式が得られる。したがって、遅れ係数CGACYLを下記式16で定義すると、式15は下記式15aで示され、気筒吸入空気量GAIRCYLNは、推定スロットル弁通過空気流量HGAIRTHを入力とする一次遅れモデルの式を用いて算出することができる。CGACYL=Vcyl×ηv/Vin16GAIRCYLN=1-CGACYL×GAIRCYLNCGACYL×HGAIRTH15a式16により遅れ係数CGACYLを算出するためには、体積効率ηvを算出することが必要である。体積効率ηvは、エンジン運転状態エンジン回転数NE,吸気圧PBA、吸気弁の作動位相、排気還流率などに依存して変化するものであるため、本実施形態では、下記式17により、気筒吸入空気量GAIRCYLNの算出に用いる体積効率ηvを算出するようにしている。ηv=GAIRCYLN/GAIRSTD17式17のGAIRSTDは下記式18により算出される理論気筒吸入空気量である。GAIRSTD=PBA×Vcyl/R×TAK18式17を用いることにより、マップやテーブルを用いることなく体積効率ηvを算出することが可能となり、常に更新されるのでエンジン特性の経時変化の影響を受けることなく最適な値を得ることできる。図4は、本実施形態における予測気筒吸入空気量GAIRCYLPの算出手法を説明するためのタイムチャートである。一点鎖線及び細い実線は、それぞれ推定スロットル弁通過空気流量HGAIRTH及び予測スロットル弁通過空気流量HGAIRTHPの推移を示し、破線及び太い実線は、それぞれ気筒吸入空気量GAIRCYLN及び予測気筒吸入空気量GAIRCYLPの推移を示す。本実施形態では、推定スロットル弁通過空気流量HGAIRTHに基づいて予測スロットル弁通過空気流量HGAIRTHPが算出され、時刻kにおける予測スロットル弁通過空気流量HGAIRTHPと、時刻kにおける気筒吸入空気量GAIRCYLNとを用いて、時刻kにおける予測気筒吸入空気量GAIRCYLPが算出される。予測気筒吸入空気量GAIRCYLPは、時刻k1における気筒吸入空気量の予測値を示す。したがって、予測気筒吸入空気量GAIRCYLPの前回値図の点PPに相当する値を用いて予測気筒吸入空気量GAIRCYLPを算出する従来手法では、前回値GAIRCYLPが予測誤差ERRPを含むのに対し、本実施形態では予測誤差ERRPを含まない気筒吸入空気量GAIRCYLNを用いて予測気筒吸入空気量GAIRCYLPが算出されるので、予測気筒吸入空気量GAIRCYLPの算出精度を向上させることができる。図4に示すように燃料噴射が吸気行程で実行されるため、燃料噴射量を予測気筒吸入空気量GAIRCYLPに基づいて算出することにより、実際の気筒吸入空気量GAIRCYLNとほぼ等しいと推定される空気量に適した量の燃料を供給することができ、過渡運転状態での空燃比制御精度を向上させることができる。時刻kにおいて算出される気筒吸入空気量GAIRCYLNは、その直前の1行程期間に吸入される空気量を示すため、図示される加速状態では、実際の吸入空気量より小さな値となる。次に弁通過空気流量モデルの逆モデル、及び吸気管モデルの逆モデルを定義する式について説明する。弁通過空気流量モデルを定義する式1を変形すると下記式21が得られる。したがって、弁通過空気流量モデルの逆モデルは、下記式22で定義される。一方、吸気管モデルを定義する式15aを変形すると、下記式23が得られる。したがって、吸気管モデルの逆モデルは、式23で定義される。HGAIRTH=GAIRCYLN-GAIRCYLN/CGACYLGAIRCYLN23目標スロットル弁通過空気流量算出部56では、目標気筒吸入空気量GAIRCYLCMDから目標スロットル弁通過空気流量GAIRTHCMDが算出されるので、式23のHGAIRTH及びGAIRCYLNを、それぞれGAIRTHCMD及びGAIRCYLCMDに置換した下記式24を用いて、目標スロットル弁通過空気流量GAIRTHCMDが算出される。GAIRTHCMD=GAIRCYLCMD-GAIRCYLCMD/CGACYLGAIRCYLCMD24また目標開度算出部57では、目標スロットル弁通過空気流量GAIRTHCMD[g/TDC]を単位時間当たりの流量[g/sec]に変換することにより、目標吸入空気流量GAIRCMDが算出され、目標吸入空気流量GAIRCMD、長期モデル補正係数KMDLL、及び短期モデル補正係数KMDLSを下記式25に適用して、補正目標吸入空気流量GAIRCMDCが算出される。GAIRCMDC=GAIRCMD/KMDLL×KMDLS25さらに、補正目標吸入空気流量GAIRCMDCを下記式26に適用して目標開度THCMDが算出される。式26は、式22のHGAIRをGAIRCMDCに置換したものである。なお、式26の演算は、図7aに示すKTHテーブルを逆検索することにより行われる。次に弁通過空気流量モデルのモデル化誤差を補正するための長期モデル補正係数KMDLL及び短期モデル補正係数KMDLSの算出手法の概要を説明する。モデル化誤差は、基本的には吸入空気流量センサ13により検出される吸入空気流量GAIRと、推定吸入空気流量HGAIRとのずれに相当するものであるが、本実施形態では吸入空気流量センサ13の検出遅れを考慮し、推定吸入空気流量HGAIRについて検出遅れ補正を行って、遅れ補正推定吸入空気流量HGAIRAFMを算出し、モデル化誤差を示すパラメータとして、モデル化誤差係数KTHERRを下記式31で定義する。KTHERR=GAIR/HGAIRAFM31モデル化誤差がなければ、モデル化誤差係数KTHERRは1.0となり、1.0との差が大きくなるほど、モデル化誤差が大きいことを示す。そこで、本実施形態では、モデル化誤差係数KTHERRに基づいて長期モデル補正係数KMDLLを算出する。一方短期モデル補正係数KMDLSは、以下のようにして算出する。推定吸入空気流量HGAIRを長期モデル補正係数KMDLLで補正した第1補正推定吸入空気流量HGAIRLについて検出遅れ補正を行って遅れ補正推定吸入空気流量HGAIRLAFMを算出し、これを下記式32に適用して補正モデル化誤差係数KTHERRSを算出する。そして、補正モデル化誤差係数KTHERRSを下記式33に適用してなまし演算を行うことにより、短期モデル補正係数KMDLSを算出する。式33のCSは、例えば0.02に設定されるなまし係数である。jは演算周期TS例えば1TDC期間に設定されるで離散化した離散化時刻である。KTHERRS=GAIR/HGAIRLAFM32KMDLS=CS×1/KTHERRS1-CS×KMDLS33なお、検出吸入空気流量GAIRが脈動しているとき、及び圧力比流量関数Ψに適用される圧力比RPが上限値に到達しているときは、なまし係数CSが0に設定され、短期モデル補正係数KMDLSは前回値に維持される。遅れ補正推定吸入空気流量HGAIRAFMHGAIRLAFMを算出するために、吸入空気流量センサ13をモデル化した吸入空気流量センサモデル以下AFMモデルというが使用される。AFMモデルは、むだ時間要素と一次遅れ要素の結合で近似することができるので、本実施形態では、下記式34により、AFMモデルを定義する。式34のjTNはむだ時間TNを演算周期TSで離散化した離散化むだ時間であり、CTDは下記式35で与えられる遅れ定数である。遅れ定数CTDは0より大きく1より小さい値をとる。式35のTDは一次遅れ要素の遅れ時定数である。遅れ補正推定吸入空気流量HGAIRLAFMは、下記式34aで与えられる。HGAIRAFM=CTD×HGAIR1-CTD×HGAIRAFM34CTD=1-e35HGAIRLAFM=CTD×HGAIRL1-CTD×HGAIRLAFM34a次に長期モデル補正係数KMDLLの算出手法を詳細に説明する。図5は、スロットル弁開度THと、モデル化誤差係数KTHERRとの関係を示す図であり、この図に示す黒丸はモデル化誤差が存在する状態における実験データを示し、実線L1は実験データからスロットル弁開度THとモデル化誤差係数KTHERRとの関係を近似するための曲線を示す。この図から明らかなように、スロットル弁開度THが所定開度THBより大きい範囲では、モデル化誤差係数KTHERRは1.0とほぼ等しくなる。そこで、本実施形態では、曲線L1TH≦THBの範囲を2次曲線で近似することとし、下記式36で誤差モデルを定義する。式36のKTHERRCORを以下誤差パラメータという。KTHERRCOR=KTHERR-1=A×TH-THB236式36は、図5に示す破線L2を示す数式に相当し、Aはモデルパラメータである。式36を用いることにより、誤差パラメータKTHERRCORを算出するために必要なモデルパラメータが1つとなり一般的な2次式では、モデルパラメータは3つである、モデルパラメータの同定演算を簡単化することができる。モデルパラメータAは、最小二乗法を用いる場合、下記式37で与えられる。同定されたモデルパラメータA及びスロットル弁開度THを式36に適用して、誤差パラメータKTHERRCORを算出し、その誤差パラメータKTHERRCORを下記式38に適用することにより、長期モデル補正係数KMDLLが得られる。KMDLL=1/KTHERR=1/KTHERRCOR138このように長期モデル補正係数KMDLLは、短期モデル補正係数KMDLSに比べて、より長い期間における検出データTH,GAIRを統計処理することにより算出されるモデルパラメータAを用いて算出される。弁通過空気流量モデルのモデル化誤差は、スロットル弁の特性ばらつきや経時変化、有効吸気管断面積の減少などに起因するので、モデルを補正するための補正係数の正確な値を得るためには、比較的長い期間で得られる検出データに基づく演算を行う必要がある。したがって、長期モデル補正係数KMDLLを用いることにより、モデル化誤差をより高い精度で補正することが可能となる。図6は、スロットル弁通過空気流量算出部51における演算処理のフローチャートである。ステップS11では、スロットル弁開度THに応じて図7aに示す開口面積流量関数テーブルを検索し、開口面積流量関数値KTHTHを算出する。ステップS12では、圧力比RP=PBA/PAに応じて図7bに示す圧力比流量関数テーブルを検索し、圧力比流量関数値ΨRPを算出する。ステップS13では、開口面積流量関数値KTHTH、圧力比流量関数値ΨRP、及び吸気温TAを前記式1に適用し、推定吸入空気流量HGAIRを算出する。ステップS14では、推定吸入空気流量HGAIR及び長期モデル補正係数KMDLLを下記式41に適用して、第1補正推定吸入空気流量HGAIRLを算出し、さらに第1補正推定吸入空気流量HGAIRL及び短期モデル補正係数KMDLSを下記式42に適用して、第2補正推定吸入空気流量HGAIRLSを算出する。HGAIRL=HGAIR×KMDLL41HGAIRLS=HGAIRL×KMDLS42ステップS15では、第2補正推定吸入空気流量HGAIRLS[g/sec]及びエンジン回転数NEを下記式43に適用し、1TDC期間当たりの吸入空気量である推定スロットル弁通過空気流量HGAIRTH[g/TDC]を算出する。式43のKCVは変換係数である。HGAIRTH=HGAIRLS×KCV/NE43図8は、予測気筒吸入空気量算出部52における演算処理のフローチャートである。ステップS21では前記式17を用いて体積効率ηvを算出する。なお、この算出には、気筒吸入空気量GAIRCYLNの前回値を用いるため、気筒吸入空気量GAIRCYLNの初期値GAIRCYLNINIの設定が必要である。本実施形態では、初期値GAIRCYLNINIは、下記式45により、理論気筒吸入空気量GAIRSTDに設定される。よって体積効率ηvの初期値は1となる式17参照。GAIRCYLNINI=GAIRSTD=PBA×Vcyl/R×TAK45ステップS22では、算出された体積効率ηvを前記式16に適用し、遅れ係数CGACYLを算出する。ステップS23では、算出された遅れ係数CGACYLを前記式15aに適用し、気筒吸入空気量GAIRCYLNを算出する。ステップS24では、下記式46により予測スロットル弁通過空気流量HGAIRTHPを算出する。式46のKPは、実験により設定される所定予測ゲインである。HGAIRTHP=HGAIRTHHGAIRTH-HGAIRTH}×KP46ステップS25では、下記式47により予測気筒吸入空気量GAIRCYLPを算出する。式47は、式15aと同形式の吸気管モデル式であり、式15aのGAIRCYLN及びHGAIRTHを、それぞれGAIRCYLN及びHGAIRTHPに代えたものである。GAIRCYLP=1-CGACYL×GAIRCYLNCGACYL×HGAIRTHP47図9は、図2のモデル補正係数算出部53において、短期モデル補正係数KMDLS及び長期モデル補正係数KMDLLを算出する処理のフローチャートである。この処理は、ECU5のCPUで、演算周期TSで実行される。ステップS101では、前記式1を用いて推定吸入空気流量HGAIRを算出する。ステップS102では、推定吸入空気流量HGAIRに応じて図10aに示すTNテーブルを検索し、無駄時間TNを算出するとともに、推定吸入空気流量HGAIRに応じて図10bに示すTDテーブルを検索し、遅れ時定数TDを算出する。なお、無駄時間TN及び遅れ時定数TDは、検出される吸入空気流量GAIRに応じて算出するようにしてもよい。ステップS103では、スロットル弁開度THを式36に適用して、誤差パラメータKTHERRCORを算出し、その誤差パラメータKTHERRCORを式38に適用することにより、長期モデル補正係数KMDLLを算出する。ステップS104では、長期モデル補正係数KMDLL及び推定吸入空気流量HGAIRを式41に適用し、第1補正推定吸入空気流量HGAIRLを算出する。ステップS105では、第1補正推定吸入空気流量HGAIRLを式34aに適用して、遅れ補正推定吸入空気流量HGAIRLAFMを算出し、ステップS106では、遅れ補正推定吸入空気流量HGAIRLAFMを式32に適用して、補正モデル化誤差係数KTHERRSを算出し、ステップS107では、補正モデル化誤差係数KTHERRSを式33に適用して、短期モデル補正係数KMDLSを算出する。図11は、モデルパラメータAを算出する学習処理のフローチャートである。この処理は、ECU5のCPUで、演算周期TSで実行される。ステップS110では、学習実行フラグFLEARNが1であるか否かを判別する。学習実行フラグFLEARNは、エンジン1の今回の運転開始エンジン始動時点からの走行距離DISTが所定値DLEARN未満であるとき1に設定される。ステップS110の答が否定NOであるときは、直ちに処理を終了する。学習実行フラグFLEARNが1であるときは、ステップS111~S116の学習処理を実行する。ステップS111及びS112における演算は、図9のステップS101及びS102における演算と同一である。ステップS113では、式34により遅れ補正推定吸入空気流量HGAIRAFMを算出し、算出した遅れ補正推定吸入空気流量HGAIRAFM及び検出吸入空気流量GAIRを式31に適用して、モデル化誤差係数KTHERRを算出するステップS114。ステップS115では、スロットル弁開度THが所定開度THBより小さいか否かを判別し、その答が肯定YESであるときは図12に示す同定演算を実行しステップS116、スロットル弁開度THが所定開度THB以上であるときは直ちに処理を終了する。図12は図11のステップS116で実行される同定演算処理のフローチャートである。ステップS140では、モデル化誤差係数KTHERRから1を減算して、誤差パラメータKTHERRCORを算出する。S141~S143では、スロットル弁開度THが、所定開度TH1~TH3TH1<TH2<TH3によって定義される第1~第4領域のいずれにあるかを判別し、スロットル弁開度THが第1領域TH<TH1にあるときは領域パラメータnを1に設定しステップS144、スロットル弁開度THが第2領域TH1≦TH<TH2にあるときは領域パラメータnを2に設定しステップS145、スロットル弁開度THが第3領域TH2≦TH<TH3にあるときは領域パラメータnを3に設定しステップS146、スロットル弁開度THが第4領域TH3≦THにあるときは領域パラメータnを4に設定するステップS147。ステップS148では、下記式61により第n領域の分母積算値XXXXnを算出する。分母積算値XXXXnは、式37の分母に対応するものであり、式61のXXXXn]zは、分母積算値XXXXn]の前回算出値である。XXXXn]=XXXXnzTH-THB461ステップS149では、分子積算値XXYnを下記式62により算出する。分子積算値XXYnは、式37の分子に対応するものであり、式62のXXYnzは、分子積算値XXYnの前回算出値である。XXYn=XXYnzKTHERRCOR×TH-THB262ステップS150では、第n領域のサンプリング数NSAMPLnを1だけインクリメントする前回値NSAMPLnzに1を加算する。ステップS151では、下記式63により、分母加重平均値XXXXTTLを算出する。式63によれば、4つの領域毎に算出される分母平均値を、均等の重み付けで加算することにより、分母加重平均値XXXXTTLが算出される。ステップS152では、下記式64により、分子加重平均値XXYTTLを算出する。式64によれば、4つの領域毎に算出される分子平均値を、均等の重み付けで加算することにより、分子加重平均値XXYTTLが算出される。ステップS153では、分子加重平均値XXYTTLを分母加重平均値XXXXTTLで除算することにより、モデルパラメータAを算出する。図12に示す処理によれば、スロットル弁開度THの変化範囲が4つの領域に分割され、各領域毎の分子平均値を、均等の重み付けで加算することにより、分子加重平均値XXYTTLが算出されるとともに、各領域毎の分母平均値を、均等の重み付けで加算することにより、分母加重平均値XXXXTTLが算出され、分子加重平均値XXYTTLを分母加重平均値XXXXTTLで除算してモデルパラメータAが算出される。これにより、サンプリングされるスロットル弁開度THの値に偏りがある場合でも、高い精度を確保することができる。図12の処理により算出されるモデルパラメータAは、エンジン運転期間エンジン始動時点からイグニッションスイッチがオフされるまでの期間の終了時点でメモリに格納され、次のエンジン運転期間において式36の演算に適用される。図13は、図2の目標気筒吸入空気量算出部54、リミット処理部55、目標スロットル弁通過空気流量算出部56、及び目標開度算出部57における処理のフローチャートである。ステップS31では、アクセルペダル操作量AP及びエンジン回転数NEに応じて目標気筒吸入空気量GAIRCYLCMDを算出する。ステップS32では、図14に示すリミット処理を実行し、目標気筒吸入空気量GAIRCYLCMDを、上限値GAIRCYLLMH及び下限値GAIRCYLLMLの範囲内の値に制限する。ステップS33では、吸気管モデル式の逆モデル式式24にリミット処理後の目標気筒吸入空気量GAIRCYLCMD及びGAIRCYLCMDを適用し、目標スロットル弁通過空気流量GAIRTHCMDを算出する。ステップS34では、目標スロットル弁通過空気流量GAIRTHCMD[g/TDC]を下記式65に適用し、目標吸入空気流量GAIRCMD[g/sec]に変換する。GAIRCMD=GAIRTHCMD×NE/KCV65ステップS35では、目標吸入空気流量GAIRCMDを前記式25に適用し、補正目標吸入空気流量GAIRCMDCを算出する。ステップS36では、式21のHGAIRとして補正目標吸入空気流量GAIRCMDCを適用し、開口面積流量関数値KTHTHを算出する。ステップS37では、ステップS36で算出した開口面積流量関数値KTHTHに応じて、図7aに示す開口面積流量関数テーブルを逆検索し、目標開度THCMDを算出する式26の演算を行う。図14は、図13のステップS32におけるリミット処理のフローチャートである。ステップS41では、スロットル弁3を駆動するアクチュエータ7の電源電圧VBDBWに応じて、スロットル弁開度の単位時間当たりの最大変化量DTHlmt[deg/ms]を算出し、スロットル弁3の制御周期TCTL例えば2ms当たりの最大変化量DTHLMT[deg/TCTL]に変換する。なお、最大変化量DTHLMTは、予め設定した値を用いてもよい。ステップS42では、検出スロットル弁開度TH及び最大変化量DTHLMTを下記式66に適用し、基本開き側限界値THOFHBを算出する。式66は、THMAXとTHDTHLMTの小さい方を選択する演算を行うものであり、THMAXは所定最大開度例えば85degである。THOFHB=MinTHMAX,THDTHLMT66ステップS43では、検出スロットル弁開度TH及び最大変化量DTHLMTを下記式67に適用し、基本閉じ側限界値THOFLBを算出する。式67は、THMINとTH-DTHLMTの大きい方を選択する演算を行うものであり、THMINは所定最小開度例えば0.5degである。THOFLB=MaxTHMIN,TH-DTHLMT67ステップS44では、検出スロットル弁開度TH及び基本開き側限界値THOFHBを下記式68に適用し、開き側限界値THOFHを算出する。式68のCTHADDは、0から1の間の値に設定される開き側なまし係数である。THOFH=THCTHADD×THOFHB-TH68ステップS45では、検出スロットル弁開度TH及び基本閉じ側限界値THOFLBを下記式69に適用し、閉じ側限界値THOFLを算出する。式69のCTHDECは、0から1の間の値に設定される閉じ側なまし係数である。THOFL=THCTHDEC×THOFLB-TH69閉じ側なまし係数CTHDECは、開き側なまし係数CTHADDより大きな値に設定され、かつなまし係数CTHDEC及びCTHADDは、ともに制御周期TCTLが長くなるほどより大きな値をとるように設定される。ステップS46では、開き側限界値THOFH及び閉じ側限界値THOFLを、弁通過空気流量モデル式式1に適用し、対応する上限吸入空気流量GAIRLMH[g/sec]及び下限吸入空気流量GAIRLML[g/sec]を算出する。さらに上限吸入空気流量GAIRLMH及び下限吸入空気流量GAIRLMLを下記式70及び71に適用し、補正上限吸入空気流量GAIRLMHLS及び補正下限吸入空気流量GAIRLMLLSを算出する。GAIRLMHLS=GAIRLMH×KMDLL×KMDLS70GAIRLMLLS=GAIRLML×KMDLL×KMDLS71ステップS47では、補正上限吸入空気流量GAIRLMHLS[g/sec]及び補正下限吸入空気流量GAIRLMLLS[g/sec]にKCV/NEを乗算して、上限スロットル弁通過空気流量GAIRTHLMH[g/TDC]及び下限スロットル弁通過空気流量GAIRTHLML[g/TDC]を算出する。ステップS48では、上限スロットル弁通過空気流量GAIRTHLMH及び下限スロットル弁通過空気流量GAIRTHLMLを吸気管モデル式式15aに適用し、上限気筒吸入空気量GAIRCYLLMH及び下限気筒吸入空気量GAIRCYLLMLを算出する。ステップS49では、目標気筒吸入空気量GAIRCYLCMDが上限気筒吸入空気量GAIRCYLLMHより大きいか否かを判別し、その答が肯定YESであるときは、目標気筒吸入空気量GAIRCYLCMDを上限気筒吸入空気量GAIRCYLLMHに設定しステップS50、ステップS51に進む。ステップS49の答が否定NOであるときは直ちにステップS51に進む。ステップS51では、目標気筒吸入空気量GAIRCYLCMDが下限気筒吸入空気量GAIRCYLLMLより小さいか否かを判別し、その答が肯定YESであるときは、目標気筒吸入空気量GAIRCYLCMDを下限気筒吸入空気量GAIRCYLLMLに設定しステップS52、処理を終了する。ステップS51の答が否定NOであるときは直ちに処理を終了する。以上のように本実施形態では、推定スロットル弁通過空気流量HGAIRTHが算出され、検出される吸気圧PBA及び吸気温TAに基づいて理論気筒吸入空気量GAIRSTDが算出され、気筒吸入空気量の前回値GAIRCYLN及び理論気筒吸入空気量GAIRSTDを用いて、体積効率ηvが算出される。また吸気管2スロットル弁下流側をモデル化した吸気管モデルの定義式15aに、スロットル弁通過空気流量HGAIRTH及び体積効率ηvを適用して、気筒吸入空気量GAIRCYLNが算出され、その気筒吸入空気量GAIRCYLNを用いて予測気筒吸入空気量GAIRCYLPが算出される。すなわち、算出された気筒吸入空気量GAIRCYLNの最新値を用いて予測気筒吸入空気量GAIRCYLPが算出されるので、前回予測値を用いる先行例と比較して、より高い精度で予測気筒吸入空気量GAIRCYLPを算出することができる。また吸気管モデル式の逆モデル式に目標気筒吸入空気量GAIRCYLCMD及び体積効率ηvを適用して、目標スロットル弁通過空気流量GAIRTHCMDが算出されるので、気筒吸入空気量GAIRCYLNの算出に適用される体積効率ηvをそのまま使用して、目標スロットル弁通過空気流量GAIRTHCMDを精度良く算出することができる。したがって、算出された目標スロットル弁通過空気流量GAIRTHCMDに応じてスロットル弁の目標開度THCMDを算出し、スロットル弁開度THをその目標開度THCMDに制御することによって、実際の気筒吸入空気量を目標気筒吸入空気量GAIRCYLCMDに精度良く制御することができる。また弁通過空気流量モデル式に検出スロットル弁開度THを適用して、推定スロットル弁通過空気流量HGAIRTHが算出されるとともに、弁通過空気流量モデル式の逆モデル式に目標スロットル弁通過空気流量GAIRTHCMDを適用して目標開度THCMDが算出される。これにより、検出スロットル弁開度THに応じた推定スロットル弁通過空気流量HGAIRTHの算出、及び目標スロットル弁通過空気流量GAIRTHCMDに応じた目標開度THCMDの算出をともに精度良く行うことができる。またスロットル弁通過空気流量の検出値、すなわち吸入空気流量センサ13により検出される吸入空気流量GAIRに基づいて、長期モデル補正係数KMDLL及び短期モデル補正係数KMDLSが算出され、これらのモデル補正係数KMDLL,KMLDS補正された弁通過空気流量モデル式を用いて、スロットル弁通過空気流量HGARITHの算出が行われるとともに、モデル補正係数KMDLL,KMLDS補正された弁通過空気流量モデル式の逆モデル式を用いて、目標開度THCMDが算出される。したがって、スロットル弁の特性ばらつきや経時変化、有効吸気管断面積の減少、運転環境の変化、あるいはセンサの温度特性などに起因する弁通過空気流量モデルのモデル化誤差を抑制し、良好な制御性能を得ることができる。またスロットル弁開度の開き側限界値THOFH及び閉じ側限界値THOHLが算出され、スロットル弁開度を算出される限界値THOFH,THOFLに設定したときのスロットル弁通過空気流量に相当する上限スロットル弁通過空気流量GAIRTHLMH及び下限スロットル弁通過空気流量GAIRTHLMLが算出される。さらに上限スロットル弁通過空気流量GAIRTHLMH及び下限スロットル弁通過空気流量GAIRTHLMLに対応する上限気筒吸入空気量GAIRCYLLMH及び下限気筒吸入空気量GAIRCYLLMLが算出され、目標気筒吸入空気量GAIRCYLCMDが,上限気筒吸入空気量GAIRCYLLMH及び下限気筒吸入空気量GAIRCYLLMLで定義される制限範囲内の値をとるようにリミット処理が行われ、リミット処理後の目標気筒吸入空気量GAIRCYLCMDを用いて目標スロットル弁通過空気流量GAIRTHCMDの算出が行われる。したがって、実現不可能なスロットル弁開度に対応した目標スロットル弁通過空気流量GAIRTHCMDが算出されることを確実に防止し、特に要求トルクが変化する過渡状態において目標開度THCMDを適切に設定することができる。本実施形態では、スロットル弁開度センサ4、吸気圧センサ8、吸気温センサ9、及び吸入空気流量センサ13がそれぞれスロットル弁開度検出手段、吸気圧検出手段、吸気温検出手段、及び空気流量検出手段に相当し、ECU5が吸入空気量制御手段の一部、スロットル弁通過空気流量取得手段、理論気筒吸入空気量算出手段、体積効率算出手段、気筒吸入空気量算出手段、気筒吸入空気量予測手段、目標気筒吸入空気量算出手段、目標スロットル弁通過空気流量算出手段、目標開度算出手段、弁通過空気流量モデル補正手段、限界値算出手段、限界スロットル弁通過空気流量算出手段、限界気筒吸入空気量算出手段、及び制限処理手段を構成する。具体的には、図6の処理がスロットル弁通過空気流量取得手段に相当し、図8のステップS21が体積効率算出手段に相当し、ステップS23が気筒吸入空気量算出手段に相当し、ステップS24及びS25が気筒吸入空気量予測手段に相当し、図13のステップS31が目標気筒吸入空気量算出手段に相当し、ステップS33が目標スロットル弁通過空気流量算出手段に相当し、ステップS34~S37が目標開度算出手段に相当し、図9の処理が弁通過空気流量モデル式補正手段に相当し、図14のステップS41~S45が限界値算出手段に相当し、ステップS46及びS47が限界スロットル弁通過空気流量算出手段に相当し、ステップS48が限界気筒吸入空気量算出手段に相当し、ステップS49~S52が制限処理手段に相当する。[変形例1]上述した実施形態では、長期モデル補正係数KMDLLの算出に適用するモデル化誤差係数KTHERRを、遅れ補正推定吸入空気流量HGAIRAFMを用いて算出するようにしたが式31、定常的なエンジン運転状態においてモデルパラメータAの学習を行うことにより、吸入空気流量センサ13の検出遅れの影響を排除することができる。したがって、遅れ補正推定吸入空気流量HGAIRAFMに代えて、検出遅れ補正を行わない推定吸入空気流量HGAIRを用いてモデル化誤差係数KTHERRを算出するようにしてもよい。[変形例2]上述した実施形態では、モデルパラメータAの算出に適用される分母加重平均値XXXTTL及び分子加重平均値XXYTTLを、それぞれ式63及び64により算出するようにしたが、これらの数式に代えて、下記式63a及び64aを用いて算出するようにしてもよい。式63a及び64aのKG1~KG4は、重み係数であり下記式72及び73を満たすように設定される。KG1>KG2>KG3>KG472KG1KG2KG3KG4=173弁通過空気流量モデル式のモデル化誤差は、図5に示したように、スロットル弁開度THが低下するほど増加する傾向がある。したがって、式72を満たすように設定される重み係数KG1~KG4を用いることにより、モデルパラメータAの同定演算に上記傾向を適切に反映させ、長期モデル補正係数KMDLLの算出精度を高めることができる。なお、重み係数KG1~KG4は、これらのうちいずれかが同じ値となるように設定してもよい。[変形例3]理論気筒吸入空気量GAIRSTDの算出方法として、公知の特許文献3WO2010/095477号公報の段落0083~0089に開示された以下の手法を用いてもよい。この場合には、吸気圧PBAが基準大気圧PA0に等しくかつ吸気温TAが基準温度TA0に等しい状態における最大気筒吸入空気量GAIRWOTを、エンジン回転数NEに対応して複数の所定エンジン回転数毎に予め求めてテーブルとして記憶しておき、エンジン運転時にはエンジン回転数NEに応じたテーブル検索により最大気筒吸入空気量GAIRWOTを算出し、検出される吸気圧PBA、基準大気圧PA0、及び最大気筒吸入空気量GAIRWOTを特許文献3の式21に適用して基準理論気筒吸入空気量GAIRSTDBを算出し、さらに検出される吸気温TA及びエンジン冷却水温TWに応じて基準理論気筒吸入空気量GAIRSTDBを補正することにより、理論気筒吸入空気量GAIRSTDを算出することができる。[第2の実施形態]本実施形態は、第1の実施形態において離散化時刻kにおける体積効率ηv、遅れ係数CGACYL、及び気筒吸入空気量GAIRCYLNの演算を複数回実行することにより、エンジンの過渡運転状態においてより正確な気筒吸入空気量GAIRCYLNを得られるようにしたものである。以下に説明する点以外は第1の実施形態と同一である。図15は、本実施形態における気筒吸入空気量算出処理のフローチャートである。この処理は、TDCパルスの発生に同期して1行程毎に4気筒エンジンであればクランク軸が180度回転する毎にECU5のCPUで実行される。ステップS81では、前記式18により理論気筒吸入空気量GAIRSTDを算出する。ステップS82では、初期化フラグFINIが1であるか否かを判別する。エンジンの始動直後は、初期化フラグFINIは0であるので、ステップS83に進み、気筒吸入空気量GAIRCYLNを理論気筒吸入空気量GAIRSTDに設定するとともに、体積効率ηvを1.0に設定する。次いで初期化フラグFINIを1に設定するステップS84。初期化フラグFINIが1であるときは、ステップS83からステップS85に進み、更新演算の実行回数を計数するインデクスパラメータiを0に設定する。以下の説明ではインデクスパラメータiを付したGAIRCYLN、ηv、及びCGACYLをそれぞれ更新気筒吸入空気量、更新体積効率、及び更新遅れ係数という。ステップS86では、更新気筒吸入空気量GAIRCYLNi=0を気筒吸入空気量の前回値GAIRCYLNに設定するとともに、更新体積効率ηvi=0を体積効率の前回値ηvに設定する。ステップS87では、インデクスパラメータiを1だけインクリメントし、ステップS88では下記式17aにより、更新体積効率ηvを算出する。ηv=GAIRCYLN/GAIRSTD17aステップS89では、下記式16aにより、更新遅れ係数CGACYLを算出する。CGACYL=Vcyl×ηv/Vin16aステップS90では、下記式15bにより、更新気筒吸入空気量GAIRCYLNを算出する。GAIRCYLN=1-CGACYL×GAIRCYLNCGACYL×HGAIRTH15bステップS91では、インデクスパラメータiが上限値iMAXに達したか否かを判別する。本実施形態では、上限値iMAXは例えばCPUの処理能力演算速度に応じて2以上の値に設定される。最初はステップS91の答は否定NOであるので、ステップS92に進み、下記式74により体積効率変化量Dηvを算出する。Dηv=|ηv-ηv|74ステップS93では、体積効率変化量Dηvが所定閾値DηvLより小さいか否かを判別し、その答が否定NOであるときはステップS87に戻り、ステップS87~S90により、更新体積効率ηv及び更新気筒吸入空気量GAIRCYLNの算出を再度実行する。ステップS91またはS93の答が肯定YESとなると、ステップS94に進み、その時点の体積効率ηv及び気筒吸入空気量GAIRCYLNを、それぞれその時点の更新体積効率ηv及び更新気筒吸入空気量GAIRCYLNに設定する。図16は、図15の処理を説明するためのタイムチャートであり、気筒吸入空気量GAIRCYLNが増加する過渡状態における理論気筒吸入空気量GAIRSTD、気筒吸入空気量GAIRCYLN、及び体積効率ηvの推移が示されている。気筒吸入空気量GAIRCYLN及び体積効率ηvの推移を示す破線は、第1の実施形態の算出手法に対応し、実線が本実施形態の算出手法に対応する。時刻kにおける演算において、細い実線の矢印がi=1の演算を示し、破線の矢印がi=2の演算を示し、一点鎖線の矢印がi=3の演算を示す。この例では、時刻kにおいてインデクスパラメータiが3となるまで更新演算がおこなわれたことが示されており、時刻k1、k2においても同様に更新演算が行われ図示省略、時刻k2において、定常状態に達した気筒吸入空気量GAIRCYLNを得ることができる。このように更新演算を行うことにより、過渡的な運転状態においてより正確な体積効率ηv及び気筒吸入空気量GAIRCYLNを得ることができる。またインデクスパラメータiが上限値iMAXに達する前でも、体積効率変化量Dηvが所定閾値DηvLより小さくなったときは、更新演算を終了するようにしたので、更新演算を適切な時期に終了することができる。本実施形態では、図15のステップS82~S94が体積効率算出手段及び気筒吸入空気量算出手段に相当する。[変形例]図15のステップS92及びS93を削除し、ステップS91の答が否定NOであるときは直ちにステップS87に戻るようにしてもよい。この変形例では、更新演算は常にインデクスパラメータiが上限値iMAXに達するまで実行される。[第3の実施形態]第2の実施形態に示した更新演算を用いる手法は、時刻kにおける繰り返し演算が必要となり、CPUの演算負荷が増大するという課題がある。そこで、本実施形態は、第2の実施形態に示した繰り返し演算を行うことなく、過渡状態においても正確な体積効率ηvが得られるようにしたものである。以下に説明する点以外は、第1の実施形態と同一である。上述した実施形態では、体積効率ηvの算出式として下記式17再掲を用いたが、過渡状態を考慮するとより厳密には体積効率ηvは下記式81を用いて算出することが望ましい。ηv=GAIRCYLN/GAIRSTD17ηv=GAIRCYLN/GAIRSTD81式81を用いると遅れ係数CGACYLは、下記式82で与えられる。ここで演算パラメータλを下記式83で定義すると、遅れ係数CGACYLは、下記式84で与えられる。CGACYL=λ×GAIRCYLN84気筒吸入空気量GAIRCYLNは、下記式15a再掲で与えられるので、これを式84に適用すると、下記式85が得られる。GAIRCYLN=1-CGACYL×GAIRCYLNCGACYL×HGAIRTH15aCGACYL=λ×{1-CGACYL×GAIRCYLNCGACYL×HGAIRTH}=λ×{CGACYL×HGAIRTH-GAIRCYLNGAIRCYLN=λ×HGAIRTH-GAIRCYLN×CGACYLλ×GAIRCYLN85式85を用いると、第2の実施形態における更新遅れ係数CGACYLに相当するパラメータは、下記数式群86で与えられるiはインデクスパラメータ。ここで、演算パラメータΓ及びHを下記式87及び88で定義すると、数式群86は、下記数式群89で表される。Γ=λ×HGAIRTH-HGAIRCYLN87H=λ×GAIRCYLN88第2の実施形態における繰り返し演算の結果得られる遅れ係数CGACYLは、理論的にはインデクスパラメータiを無限大とした極値に相当するので、演算パラメータΓが-1から1の間の値をとるという条件の下で、下記式90で与えられる。式90から明らかなように、遅れ係数CGACYLは、初期値CGACYL0に依存しない。式90の演算パラメータΓ及びH並びにΓ及びHに含まれるλをもとの数式に戻すことにより、下記式91が得られる。したがって、繰り返し演算により得られる体積効率ηvは、式92で与えられる。図17は、要求トルクTRQCMDが急激に増加する過渡状態における動作例のタイムチャートである。図17の破線が要求トルクTRQCMDの推移を示し、実線が実出力トルクTRQの推移を示す。このように、本実施形態に示した算出手法による体積効率ηvを使用することにより、良好なトルク制御性能を得ることができる。以上のように本実施形態では、図8のステップS21において式92を用いて体積効率ηvを算出するようにしている。式92を用いることにより、気筒吸入空気量GAIRCYLN及び理論気筒吸入空気量GAIRSTDとともに、気筒容積Vcylと吸気管容積Vinとの比、及び推定スロットル弁通過空気流量HGAIRTHをさらに用いて体積効率ηvが算出され、繰り返し更新演算を行うことなく、すなわちCPUの負荷の大幅に増加させることなく、体積効率ηvの算出精度を高めることができる。[第4の実施形態]本実施形態は、第1の実施形態における長期モデル補正係数KMDLLの算出手法を変更したものである。以下に説明する点以外は第1の実施形態と同一である。第1の実施形態では、スロットル弁開度THを式36に適用して誤差パラメータKTHERRCORを算出したが、本実施形態では、式36で示される関係をKTHERRCORテーブルとして設定し、テーブル検索により誤差パラメータKTHERRCORを算出するようにしたものである。図18は、KTHERRCORテーブルを表形式で示す図である。KTHERRCORテーブルには、スロットル弁開度THについて1[deg]からTHB-1[deg]に対応して、テーブル設定値KTHERRCOR1,KTHERRCOR2,...,KTHERRCORTHB1KTHERRCORp,p=1~THB-1が設定され、所定開度THBに対応して1が設定されている。図19は、KTHERRCORテーブルを更新する処理のフローチャートである。この処理は、ECU5のCPUで演算周期TSで実行される。ステップS161~S165は、図11のステップS111~S115と同一である。ステップS166では、図20に示す設定値更新処理を実行し、テーブル設定値KTHERRCORpp=1~THB-1の更新を行う。図20のステップS171では、図19のステップS164で算出されたモデル化誤差係数KTHERRから1を減算して誤差パラメータKTHERRCORを算出する。ステップS172では、スロットル弁開度THに応じて更新すべきテーブル設定値KTHERRCORpSを選択する。例えばスロットル弁開度THが0.3[deg]であるときは、テーブル設定値KTHERRCOR1pS=1が選択され、スロットル弁開度THが1.8[deg]であるときは、テーブル設定値KTHERRCOR2pS=2が選択される。ステップS173では、更新すべきテーブル設定値KTHERRCORpSの算出に適用されたデータ数NDATApSを1だけインクリメントする。NDATApSzは、前回更新時のデータ数である。ステップS174では、ステップS171で算出した誤差パラメータKTHERRCORを下記式101に適用し、テーブル設定値KTHERRCORpSの更新を行う。式101のKTHERRCORpSzは、更新前のテーブル設定値である。図20の処理により、テーブル設定値KTHERRCORpの更新を行うことにより、スロットル弁開度THの1[deg]からTHB-1[deg]までの範囲において、テーブルの設定格子点TH=1,2,...,THB-1毎に設定値の更新が行われ、更新回数が増加するほど各格子点におけるテーブル設定値KTHERRCORpの精度を高め、長期モデル補正係数KMDLLの算出精度を高めることができる。[変形例]第4の実施形態では、誤差パラメータKTHERRCORが設定されたKTHERRCORテーブルを使用するようにしたが、誤差パラメータKTHERRCORに代えて長期モデル補正係数KMDLLが設定されたKMDLLテーブルを使用し、KMDLLテーブルの設定値を更新するようにしてもよい。また、KTHERRCORテーブルにおけるスロットル弁開度THの格子点は1[deg]間隔で設定したが、これに限るものではなく、1[deg]より大きい間隔または小さい間隔で設定するようにしてもよい。なお本発明は上述した第1~第4実施形態に限るものではなく、種々の他の変形が可能である。例えば、上述した実施形態では、上限気筒吸入空気量GAIRCYLLMH及び下限気筒吸入空気量GAIRCYLLMLを用いてリミット処理を行うようにしたが、例えばエンジンの加速時は上限気筒吸入空気量GAIRCYLLMHのみを用いてリミット処理を行い、減速時は下限気筒吸入空気量GAIRCYLLMLのみを用いてリミット処理を行うようにしてもよい。また上述した実施形態では、気筒吸入空気量GAIRCYLN及び予測気筒吸入空気量GAIRCYLPの算出には、推定スロットル弁通過空気流量HGAIRTHを使用したが、吸入空気流量センサ13により検出される吸入空気流量GAIR[g/sec]を1TDC期間毎の流量に変換した検出スロットル弁通過空気流量GAIRTHを使用するようにしてもよい。また、上述した実施形態では大気圧センサ33により検出した大気圧PAを用いて推定吸入空気流量HGAIRを算出するようにしたが、公知の大気圧推定手法例えば米国特許第6016460号公報参照を用いて算出した推定大気圧HPAを用いて推定吸入空気流量HGAIRを算出するようにしてもよい。また上述した実施形態では、本発明をガソリン内燃エンジンに適用した例を示したが、本発明はディーゼル内燃エンジンにも適用可能である。また本発明は、クランク軸を鉛直方向とした船外機などのような船舶推進機用エンジンなどにも適用が可能である。"