"排気循環器及び可変ノズルターボを有するエンジンに対する燃料噴射量の設定値、エンジン回転数の設定値、前記エンジンの吸気圧の測定値及び新気量の測定値を取得するステップと、前記吸気圧の測定値及び前記新気量の測定値と前記エンジンが有するむだ時間特性をパディ近似により表すことで追加される前記吸気圧に関する内部状態及び前記新気量に関する内部状態とを表す状態ベクトルと第1の行列との積と、前記排気循環器のバルブ開度と前記可変ノズルターボのノズル開度とを表す操作量ベクトルと第2の行列との積との和により1単位時間後の状態ベクトルが表される合成状態方程式における前記第1の行列と前記第2の行列とを用いて前記吸気圧に関する内部状態及び前記新気量に関する内部状態を状態オブザーバによって推定し、前記エンジンの吸気圧の測定値及び新気量の測定値と併せて状態ベクトルを生成するステップと、前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じた前記新気量の目標値及び前記吸気圧の目標値と、生成した前記状態ベクトルとから、前記第1の行列と前記第2の行列とから得られる行列を用いてモデル予測制御を実施することで、前記排気循環器のバルブ開度の操作量と前記可変ノズルターボのノズル開度の操作量とを算出し、前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じた、前記排気循環器のバルブ開度の基準値及び前記可変ノズルターボのノズル開度の基準値を加算することで、前記排気循環器のバルブ開度の指令値及び前記可変ノズルターボのノズル開度の指令値を算出するステップと、をプロセッサに実行させ、前記第1の行列により、前記状態ベクトルに含まれる要素間の相互作用を表し、前記第2の行列が、前記排気循環器のバルブ開度が前記吸気圧に影響を与える干渉要素と、前記可変ノズルターボのノズル開度が前記新気量に影響を与える干渉要素とを含む、前記エンジンの状態方程式における行列と、前記排気循環器のバルブ開度が前記吸気圧に関する内部状態に影響を与える干渉要素と、前記可変ノズルターボのノズル開度が前記新気量に関する内部状態に影響を与える干渉要素とを含む、前記むだ時間特性をパディ近似で表すことで規定される状態方程式における行列とを含む、エンジン制御プログラム。前記第1の行列が、前記吸気圧が前記新気量に影響を与える干渉要素と前記新気量が前記吸気圧に影響を与える干渉要素とを含む、前記エンジンの状態方程式における行列と、前記吸気圧に関する内部状態が前記新気量に関する内部状態に影響を与える干渉要素と前記新気量に関する内部状態が前記吸気圧に関する内部状態に影響を与える干渉要素とを含む、前記むだ時間特性をパディ近似で表すことで規定される状態方程式における行列と、前記吸気圧に関する内部状態及び前記新気量に関する内部状態が前記新気量に影響を与える干渉要素、及び前記吸気圧に関する内部状態及び前記新気量に関する内部状態が前記吸気圧に影響を与える干渉要素とを含む請求項1記載のエンジン制御プログラム。前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じて前記第1の行列及び前記第2の行列を更新するステップをさらに前記プロセッサに実行させるための請求項1又は2記載のエンジン制御プログラム。前記第1の行列及び前記第2の行列が、前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じて変化するむだ時間の補償量に応じて算出される請求項3記載のエンジン制御プログラム。前記モデル予測制御の処理が、前記排気循環器のバルブ開度の操作量の時系列値及び前記可変ノズルターボのノズル開度の操作量の時系列値が、前記排気循環器のバルブ開度及び前記可変ノイズターボのノズル開度の飽和による制約条件を満たすか判断するステップと、前記制約条件が満たされないと判断された場合には、前記制約条件を満たしつつ、追従誤差に関する値と、操作量に関する値とを重み付け加算した評価関数の値を最小化するように、前記排気循環器のバルブ開度の操作量の時系列値及び前記可変ノズルターボのノズル開度の操作量の時系列値を計算し直すステップと、を含む請求項1記載のエンジン制御プログラム。排気循環器及び可変ノズルターボを有するエンジンに対する燃料噴射量の設定値、エンジン回転数の設定値、前記エンジンの吸気圧の測定値及び新気量の測定値を取得するステップと、前記吸気圧の測定値及び前記新気量の測定値と前記エンジンが有するむだ時間特性をパディ近似により表すことで追加される前記吸気圧に関する内部状態及び前記新気量に関する内部状態とを表す状態ベクトルと第1の行列との積と、前記排気循環器のバルブ開度と前記可変ノズルターボのノズル開度とを表す操作量ベクトルと第2の行列との積との和により1単位時間後の状態ベクトルが表される合成状態方程式における前記第1の行列と前記第2の行列とを用いて前記吸気圧に関する内部状態及び前記新気量に関する内部状態を状態オブザーバによって推定し、前記エンジンの吸気圧の測定値及び新気量の測定値と併せて状態ベクトルを生成するステップと、前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じた前記新気量の目標値及び前記吸気圧の目標値と、生成した前記状態ベクトルとから、前記第1の行列と前記第2の行列とから得られる行列を用いてモデル予測制御を実施することで、前記排気循環器のバルブ開度の操作量と前記可変ノズルターボのノズル開度の操作量とを算出し、前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じた、前記排気循環器のバルブ開度の基準値及び前記可変ノズルターボのノズル開度の基準値を加算することで、前記排気循環器のバルブ開度の指令値及び前記可変ノズルターボのノズル開度の指令値を算出するステップと、を含み、前記第1の行列により、前記状態ベクトルに含まれる要素間の相互作用を表し、前記第2の行列が、前記排気循環器のバルブ開度が前記吸気圧に影響を与える干渉要素と、前記可変ノズルターボのノズル開度が前記新気量に影響を与える干渉要素とを含む、前記エンジンの状態方程式における行列と、前記排気循環器のバルブ開度が前記吸気圧に関する内部状態に影響を与える干渉要素と、前記可変ノズルターボのノズル開度が前記新気量に関する内部状態に影響を与える干渉要素とを含む、前記むだ時間特性をパディ近似で表すことで規定される状態方程式における行列とを含む、エンジン制御方法。排気循環器及び可変ノズルターボを有するエンジンに対する燃料噴射量の設定値、エンジン回転数の設定値、前記エンジンの吸気圧の測定値及び新気量の測定値を取得するデータ取得部と、前記吸気圧の測定値及び前記新気量の測定値と前記エンジンが有するむだ時間特性をパディ近似により表すことで追加される前記吸気圧に関する内部状態及び前記新気量に関する内部状態とを表す状態ベクトルと第1の行列との積と、前記排気循環器のバルブ開度と前記可変ノズルターボのノズル開度とを表す操作量ベクトルと第2の行列との積との和により1単位時間後の状態ベクトルが表される合成状態方程式における前記第1の行列と前記第2の行列とを用いて前記吸気圧に関する内部状態及び前記新気量に関する内部状態を推定し、前記エンジンの吸気圧の測定値及び新気量の測定値と併せて状態ベクトルを生成する状態オブザーバと、前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じた前記新気量の目標値及び前記吸気圧の目標値と、生成した前記状態ベクトルとから、前記第1の行列と前記第2の行列とから得られる行列を用いてモデル予測制御を実施することで、前記排気循環器のバルブ開度の操作量と前記可変ノズルターボのノズル開度の操作量とを算出し、前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じた、前記排気循環器のバルブ開度の基準値及び前記可変ノズルターボのノズル開度の基準値を加算することで、前記排気循環器のバルブ開度の指令値及び前記可変ノズルターボのノズル開度の指令値を算出するモデル予測制御部と、を有し、前記第1の行列により、前記状態ベクトルに含まれる要素間の相互作用を表し、前記第2の行列が、前記排気循環器のバルブ開度が前記吸気圧に影響を与える干渉要素と、前記可変ノズルターボのノズル開度が前記新気量に影響を与える干渉要素とを含む、前記エンジンの状態方程式における行列と、前記排気循環器のバルブ開度が前記吸気圧に関する内部状態に影響を与える干渉要素と、前記可変ノズルターボのノズル開度が前記新気量に関する内部状態に影響を与える干渉要素とを含む、前記むだ時間特性をパディ近似で表すことで規定される状態方程式における行列とを含む、エンジン制御装置。従って、本技術の目的は、吸気圧制御系と新気量制御系とが互いに干渉を生じる制御系において安定的な制御を行うためのモデル予測制御技術を提供することである。本エンジン制御方法は、A排気循環器及び可変ノズルターボを有するエンジンに対する燃料噴射量の設定値、エンジン回転数の設定値、エンジンの吸気圧の測定値及び新気量の測定値を取得するステップと、B吸気圧の測定値及び新気量の測定値とエンジンが有するむだ時間特性をパディ近似により表すことで追加される吸気圧に関する内部状態及び新気量に関する内部状態とを表す状態ベクトルと第1の行列との積と、排気循環器のバルブ開度と可変ノズルターボのノズル開度とを表す操作量ベクトルと第2の行列との積との和により1単位時間後の状態ベクトルが表される合成状態方程式における第1の行列と第2の行列とを用いて吸気圧に関する内部状態及び新気量に関する内部状態を状態オブザーバによって推定し、第2の状態ベクトルを生成するステップと、C燃料噴射量の設定値及びエンジン回転数の設定値に応じた新気量の目標値及び吸気圧の目標値と、生成した第2の状態ベクトルと、第1の行列と第2の行列とから得られる行列とを用いてモデル予測制御を実施することで、排気循環器のバルブ開度の操作量と可変ノズルターボのノズル開度の操作量とを算出し、燃料噴射量の設定値及びエンジン回転数の設定値に応じた、排気循環器のバルブ開度の基準値及び可変ノズルターボのノズル開度の基準値を加算することで、排気循環器のバルブ開度の指令値及び可変ノズルターボのノズル開度の指令値を算出するステップとを含む。そして、第1の行列により、状態ベクトルに含まれる要素間の相互作用を表し、第2の行列が、排気循環器のバルブ開度が吸気圧に影響を与える干渉要素と、可変ノズルターボのノズル開度が新気量に影響を与える干渉要素とを含む、エンジンの状態方程式における行列と、排気循環器のバルブ開度が吸気圧に関する内部状態に影響を与える干渉要素と、可変ノズルターボのノズル開度が新気量に関する内部状態に影響を与える干渉要素とを含む、むだ時間特性をパディ近似で表すことで規定される状態方程式における行列とを含む。吸気圧制御系と新気量制御系とが互いに干渉を生じている制御系にモデル予測制御を適用した場合、安定的な制御を行うことができるようになる。[実施の形態1]一般的なプラント例えばディーゼルエンジンを状態方程式で表現すると以下のようになる。xpkは2次元の状態ベクトルを表し、upkは2次元の操作量ベクトル又は入力ベクトルを表し、ypkは2次元測定出力ベクトルを表す。xmafkは新気量、xmapkは吸気圧、uegrkはEGRバルブ開度、uvntkはVNTノズル開度を表す。さらに、むだ時間が存在する場合を考える。その場合の状態方程式は以下のように表される。但し、lは、むだ時間分の遅延を表している。一方、むだ時間要素esτでは、状態方程式表現ができないため、パディ近似により線形化して表現することにする。そうすると、むだ時間要素esτの2次のパディ近似による伝達関数表現は、以下のようになる。むだ時間は、図5に示したように、新気量制御系のEGRバルブ開度から新気量MAFへの直達項の係数τ11として、吸気圧制御系のVNTノズル開度から吸気圧MAPへの直達項の係数τ22として、EGRバルブ開度から吸気圧MAPへの干渉項の係数τ12として、VNTノズル開度から新気量MAFへの干渉項の係数τ21として現れる。そしてそれらは異なる値を持つため、実際のエンジンの特性に合わせて設定する。なお、τ11乃至τ22をむだ時間補償量と記す。1-3式を、周知の方法で連続系の状態方程式に変換して、さらに離散系の状態方程式に変換すると、以下のように表される。なお、τと以下で表す行列Ad、Bd及びCdとの関係式についてはエンジン制御装置内のメモリなどに保持しておく。本実施の形態では、むだ時間のない系1-1式と、むだ時間要素をパディ近似したモデルの系1-4式とがカスケードに接続されていると考えて、その合成系を考える。そうすると、以下のように表される。ここで1-8式を改めて書き直す。1-11式に含まれる行列A及びBは、吸気圧制御系と新気量制御系とで干渉がある制御系においてむだ時間を適切に取り扱う上で主要な役割を有している。図6に示すように、行列Aにおける要素1は、ディーゼルエンジンの元々の状態方程式に含まれる行列Apであり、現時刻の新気量MAFが1単位時間後の新気量MAFに対して影響を与える要素a11と、現時刻の吸気圧MAPが1単位時間後の新気量MAFに対して影響を与える干渉要素a12と、現時刻の新気量MAFが1単位時間後の吸気圧MAPに対して影響を与える干渉要素a21と、現時刻の吸気圧MAPが1単位時間後の吸気圧MAPに対して影響を与える要素a22とを含む。さらに、行列Aにおける要素2は、むだ時間のパディ近似による状態方程式で導入された内部状態xmaf1が1単位時間後の新気量MAFに影響を与える要素b11c11と、内部状態xmap1が1単位時間後の吸気圧MAPに影響を与える要素b21c11を含む。これらは、むだ時間のパディ近似による状態方程式のBd及びCdに由来している。同様に、行列Aにおける要素3は、むだ時間のパディ近似による状態方程式で導入された内部状態xmaf1が1単位時間後の吸気圧MAPに影響を与える要素b12c11と、内部状態xmap1が1単位時間後の吸気圧MAPに影響を与える要素b22c23を含む。これらは、むだ時間のパディ近似による状態方程式のBd及びCdに由来している。さらに、行列Aにおける要素4乃至7は、むだ時間のパディ近似による状態方程式における行列Adに由来する。そして、行列Aにおける要素4は、現時刻の内部状態xmaf1及びxmaf2が1単位時間後の内部状態xmaf1及びxmaf2に影響を与える要素である。また、行列Aの要素5は、現時刻の内部状態xmap1及びxmap2が1単位時間後の内部状態xmaf1及びxmaf2に影響を与える干渉要素である。また、行列Aの要素6は、現時刻の内部状態xmaf1及びxmaf2が1単位時間後の内部状態xmap1及びxmap2に影響を与える干渉要素である。また、行列Aの要素7は、現時刻の内部状態xmap1及びxmap2が1単位時間後の内部状態xmap1及びxmap2に影響を与える要素である。さらに、行列Bにおける要素8は、エンジンの元々の状態方程式における行列Bpに由来する。そして、行列Bにおける要素8は、現時刻のEGRバルブ開度の指令値uegrが1単位時間後の新気量に影響を与える要素b11と、現時刻のVNTノズル開度の指令値uvntが1単位時間後の新気量に影響を与える干渉要素b12と、現時刻のEGRバルブ開度の指令値uegrが1単位時間後の吸気圧に影響を与える干渉要素b21と、現時刻のVNTノズル開度の指令値uvntが1単位時間後の吸気圧に影響を与える要素b22とを含む。さらに、行列Bにおける要素9乃至12は、むだ時間のパディ近似による状態方程式における行列Bdに由来する。行列Bにおける要素9は、現時刻のEGRバルブ開度の指令値uegrが、むだ時間のパディ近似による状態方程式における内部状態xmaf1及びxmaf2に影響を与える要素dd11及びdd21を含む。また、行列Bにおける要素10は、現時刻のEGRバルブ開度の指令値uegrが、内部状態xmap1及びxmap2に影響を与える干渉要素bd31及びbd41を含む。さらに、行列Bにおける要素11は、現時刻のVNTノズル開度の指令値uvntが、内部状態xmaf1及びxmaf2に影響を与える干渉要素bd21及びbd22を含む。また、行列Bにおける要素12は、現時刻のVNTノズル開度の指令値uvntが、内部状態xmap1及びxmap2に影響を与える要素bd32及びbd42を含む。そして、1-9乃至1-12式は、むだ時間要素の位相遅れを表現しているため、1-1式に対して内部モデルとして位相遅れ特性を組み込んだことに相当する。なお、パディ近似の次数を2次の例を示したが、3次、4次と次数を上げることも可能である。次数を上げると、近似精度を高くできる。なお、むだ時間おおよそ40ms毎に1次次数を上げると良いことが実験的に分かっている。このようにパディ近似で増やした状態xmaf1、xmaf2、xmap1、xmap2については外部から観測できないので、状態オブザーバを導入して推定する。本実施の形態では、以下に示すような同次元オブザーバを採用することにする。ここでオブザーバゲインLは、以下のような行列式が例えば6次のバタワースフィルタの離散パルス伝達関数に一致するように定める。1-11及び1-12式で1-13式を展開すると以下のようになる。この行列式は6次式であり、6次のバタワースフィルタの離散パルス伝達関数の分母は、以下のように表される。このように1-18式も6次式となるので、z5の係数同士、z4の係数同士、z3の係数同士、z2の係数同士、z1の係数同士、z0の係数同士が一致するものとして、連立方程式を構成して解く。そうすると、l11、l12、l21、l22、l31、l32、l61及びl62が算出される。以下でも述べるように、行列A及びBは変更される場合があるので、その場合にはLも計算し直すことになる。なお、今回の場合、xmaf及びxmapは観測できるので、xmaf1、xmaf2、xmap1、xmap2のみについて最小状態推定オブザーバを採用するようにしても良い。次に、1-9式のプラントモデルについてモデル予測制御による最適操作量を導出するものとする。ここでは一般的な評価関数を最小にする操作量を求める有限時間最適制御問題として解くことにする。但し、upk=Δupkupk-1である。Qi及びRiは、重み係数である。重みQとRの設定により、第1項と第2項とはトレードオフ関係となる。具体的には、Rが一定で、重みQを大きくすることは、同じukのエネルギーで、目標値への追従誤差を相対的に小さくするような解が得られ、逆に重みQを小さくすることは、目標値への追従誤差は相対的に大きくてもよいというような解が得られる。一方、重みQが一定で、Rを大きくすることは、少ないエネルギーで一定の目標値追従性を達成することを意味し、逆に小さくすることは、大きなエネルギーで一定の目標値追従性を達成することを意味している。また、Yを2N次元の状態ベクトル、Uを2N次元の操作量ベクトルとして拡大系を定義すると、1-9式は以下のように書き換えられる。さらに、目標値に関しても、同様に2N行の目標ベクトルとして拡大系を表す。なお、Tk=Tskであり、sk=[smaf,smap]Tである。Tは以下のように表される。ここで、Trは参照軌道の1次遅れの時定数である。Trを短くすると即応性のある応答を、長くするとゆっくりした応答となる。1-21乃至1-24式を用いて1-19式を書き換えると以下のようになる。ここで、Ekを以下のように定義する。そうすると、1-26式を以下のように変形できる。1-28式のJkを最適にするUoptkは、以下の条件から1-28式は1-30式のように書くことができる。従って、ΔUoptkは以下のように表される。最終的に、Δuoptkは、以下のように表される。なお、ΔuoptkはΔupkとも表される。但し、xk=[xmaf,xmap,xmaf1,xmaf2,xmap1,xmap2]Tについては一部の状態を観測できないので、上でも述べたように状態オブザーバによって推定された値x^k|k-1x又はxprは、xハットを表すを用いる。すなわち、以下のように表される。本実施の形態では、1-35式に基づく制御を行う。図7に、本技術の実施の形態に係るエンジンの一例としてディーゼルエンジンを示す。エンジン本体1には、エンジン本体1からの排ガスを供給する排気循環器EGRと、排ガスの圧力にてタービンを回して新気FreshAirを圧縮してエンジン本体1に供給する可変ノズルターボVNTとが接続されている。可変ノズルターボVNTのノズル開度を調整することによって、可変ノズルターボVNTのタービンの回転が調整され、吸気圧MAPセンサで測定される吸気圧MAPが調整される。一方、排気循環器EGRに設けられているEGRバルブのバルブ開度を調整することによって、新気量MAFセンサで測定される新気量MAFが調整される。本実施の形態に係るエンジン制御装置1000には、MAPセンサからの吸気圧測定値と、MAFセンサからの新気量測定値と、外部から与えられる燃料噴射量の設定値と、同じく外部から与えられるエンジン回転数の設定値とが入力されるようになっている。また、エンジン制御装置1000からは、EGRバルブのバルブ開度の操作量がEGRバルブに出力され、VNTノズルのノズル開度の操作量がVNTノズルに出力されるようになっている。本実施の形態に係るエンジン制御装置1000のブロック線図を図8に示す。エンジン制御装置1000は、計画器100と、モデル予測制御部110と、状態オブザーバ120とを有する。すなわち、燃料噴射量の設定値とエンジン回転数の設定値とが入力されて、燃料噴射量の値及びエンジン回転数の値に対応付けてEGRバルブ開度の基準値及びVNTノズル開度の基準値の組み合わせと新気量MAFの目標値及び吸気圧MAPの目標値の組み合わせが登録されている計画器100から、燃料噴射量の設定値及びエンジン回転数の設定値に対応する新気量MAFの目標値及び吸気圧MAPの目標値skと、EGRバルブ開度の基準値及びVNTノズル開度の基準値urefkとを読み出す。そして、計画器100が出力した新気量MAFの目標値及び吸気圧MAPの目標値skに、参照軌道行列Tを乗じて2N次元の目標ベクトルTkを算出すると共に、状態オブザーバ120の出力xk|k-1にΨcを乗じて2N次元の応答ベクトルΨcxkを算出する。また、1サンプル前の操作量upk-1に対してΓcを乗じたΓcupk-1も算出する。そうすると、2N次元の誤差ベクトルEk=Tk-Ψcxk-Γcupk-1が得られる。この誤差ベクトルEkは、MPC主制御部111に入力される。MPC主制御部111は、1-35式に従ってM2x2NH1ΘcTQQ又はQbarは、Qバーを表す。をEkに乗じてΔupkを算出して出力する。積分要素112では、Δupkを基に例えばΔupkupk-1、upkを算出する。そして、この操作量upkは、計画器100が出力したフィードフォワード値urefkと加算されてエンジン特性130に出力される。すなわち、upkurefkによって、VNTノズル開度及びEGRバルブ開度が制御される。なお、遅れ要素113は、操作量upkを1サンプル遅らせてupk-1を出力する。このupk-1は、状態オブザーバ120にも出力される。状態オブザーバ120は、1-13式に従った演算を実施する。具体的には、エンジン特性130の測定値ypk=xpkに対して-Cxk|k-1を加算してxpk-Cxk|k-1を算出する。さらにこれにLを乗じて、Lxpk-Cxk|k-1を算出する。一方、モデル予測制御部110の出力upk-1に対してBを乗じてBupk-1を算出し、xk|k-1にAを乗じてAxk|k-1を算出する。そして、xk|k-1、Bupk-1及びLxpk-Cxk|k-1を加算して遅れ要素121に入力して、1サンプル遅延させてから、モデル予測制御部110に次のxk|k-1として出力する。次に、エンジン制御装置1000の詳細な処理内容について図9及び図10を用いて説明する。まず、エンジン制御装置1000は、k=0に設定する図9:ステップS1。また、エンジン制御装置1000は、初期設定を実施するステップS3。具体的には、エンジン制御装置1000は、k=0におけるエンジン回転数及び燃料噴射量の設定値に応じたプラント行列A及びB並びにむだ時間補償量τ=[τ11,τ12,τ21,τ22]を取得すると共に、重み行列Q及びR、予測期間N、初期値x0及びu0、時定数Tr等を取得する。なお、これらのデータについては、予めメモリなどに保持しておき、該当データを読み出す。さらに、エンジン制御装置1000は、初期計算として、参照軌道行列Tを算出し、設定するステップS5。さらに、エンジン制御装置1000は、モデル予測制御部110及び状態オブザーバ120において用いる行列Hc1=ΘcTQΘcR1及びΘcTQ、L、Γc、Ψcなどを算出し、設定するステップS7。さらに、Hc1ΘcTQも算出してしまっても良い。処理は端子Aを介して、図10の処理に移行する。図10の処理の説明に移行して、エンジン制御装置1000は、kを1インクリメントしステップS9、MAPセンサ及びMAFセンサから測定値ypkを取得し、さらに計画器100から現時刻におけるエンジン回転数及び燃料噴射量を取得して当該エンジン回転数及び燃料噴射量に対応する目標値skを取得するステップS11。なお、エンジン制御装置1000は、計画器100から現時刻におけるエンジン回転数及び燃料噴射量の設定値に対応する基準値urefも取得する。そして、エンジン制御装置1000のモデル予測制御部110は、目標値ベクトルTk=Tskを算出するステップS13。また、エンジン制御装置1000の状態オブザーバ120は、1-13式に従ってx^k|k-1を算出するステップS15。そして、モデル予測制御部110のMPC主制御部111は、応答ベクトルΨcxk|k-1を算出すると共に、Γcupk-1を算出して、誤差ベクトルEk=Tk-Ψcxk|k-1-Γcupkを算出するステップS17。そして、モデル予測制御部110は、1-32式に従ってΔUoptkを算出するステップS19。さらに、モデル予測制御部110のMPC主制御部111は、ΔUoptkに対してM2x2Nを乗ずることで、Δuoptkを算出するステップS21。その後、モデル予測制御部110は、ΔuoptkからΔuoptkuoptk-1によりuoptk=upkを算出し、さらに基準値urefkを加算することで最終的な操作量の指令値を算出し、当該指令値をエンジン1のVNTノズル及びEGRバルブに出力するステップS23。また、モデル予測制御部110は、upk-1=uoptkと設定して、時刻k1についての処理で用いるようにメモリに保持しておくステップS25。そして、エンジン制御装置1000は、kが予め設定された最大値kmaxに達したか判断するステップS27。kがkmaxに達していない場合には、処理はステップS9に戻る。一方、kがkmaxに達した場合には、処理を終了する。例えば、ステップS1に戻って、初期設定から処理を開始する。以上のような処理を実施することで、むだ時間が存在する多入力多出力制御系MIMOにおいて、むだ時間要素の位相遅れに起因するシステムの応答遅れや不安定性を補償して安定化することができる。なお、図9のステップS3の内、システム行列A、B、重みQ、R、むだ時間補償量τの読み込み、および、ステップS7の定数行列の計算を、図10のステップS9とステップS11の間に入れることにより、毎回、運転条件が変わってシステム行列が変化する場合に対応することもできる。例えば、図11Aに示すような燃料噴射量qが設定され、図11Bに示すようなエンジン回転数rpmが設定された場合を想定する。なお、図11Aの縦軸は燃料噴射量Qを表し、横軸は時間を表し、図11Bの縦軸はエンジン回転数rpmを表し、横軸は時間を表す。なお、むだ時間が1サンプル遅れ、3サンプル遅れ、6サンプル遅れの各ケースについて本実施の形態を実施した場合を検討する。図11Cは、新気量MAFの時間変化を表す。図11Cにおいて、太点線は目標値を表しており、このケースではむだ時間の補正量が多くなると目標値に近付くことが分かる。同様に、図11Dは、吸気圧MAPの時間変化を表す。図11Dにおいて、太点線は目標値を表しており、このケースではむだ時間の補正量が多くなると目標値に近付くことが分かる。図11Eは、EGRバルブ開度の操作量指令値の時間変化を表しており、図11Fは、VNTノズル開度の操作量指令値の時間変化を表している。むだ時間の補正量に応じて、大きく波形が異なっており、本実施の形態によるモデル予測制御が適切に作用していることが分かる。[実施の形態2]1-19式ではΔupを用いて評価関数を定義していたが、単にupを用いて定義することも可能である。この場合評価関数は以下のように表すことができる。これを拡大系に変形した場合、Ekは以下のように表される。従って、1-33式の代わりに、本実施の形態では、uoptkは以下のように表される。2-3式に従ってエンジンを制御する場合のエンジン制御装置1000は、例えば図12のようなブロック線図で表される。エンジン制御装置1000は、計画器100と、第2モデル予測制御部210と、状態オブザーバ120とを有する。2-3式を採用することにしたので、モデル予測制御部110の代わりに第2モデル予測制御部210が採用されている。具体的には、燃料噴射量の設定値とエンジン回転数の設定値とが入力されて、燃料噴射量の値及びエンジン回転数の値に対応付けてEGRバルブ開度の基準値及びVNTノズル開度の基準値の組み合わせと新気量MAFの目標値及び吸気圧MAPの目標値の組み合わせが登録されている計画器100から、燃料噴射量の設定値及びエンジン回転数の設定値に対応するEGRバルブ開度の基準値及びVNTノズル開度の基準値urefkと、新気量MAFの目標値及び吸気圧MAPの基準値skとを読み出す。そして、計画器100が出力した新気量MAFの目標値及び吸気圧MAPの目標値skに、参照軌道行列Tを乗じて2N次元の目標ベクトルTkを算出すると共に、状態オブザーバ120の出力xk|k-1にΨcを乗じて2N次元の応答ベクトルΨcxk|k-1を算出する。第1の実施の形態とは異なり、第2モデル予測制御部210内でのフィードバックはない。従って、2N次元の誤差ベクトルEk=Tk-Ψcxk|k-1が算出される。この誤差ベクトルEkは、第2MPC主制御部211に入力される。第2MPC主制御部211は、2-3式に従ってM2x2NH1ΘcTQをEkに乗じてupkを算出して出力する。そして、この操作量upkは、計画器100が出力したフィードフォワード値urefkと加算されてエンジン特性130に出力される。すなわち、upkurefkによって、VNTノズル開度及びEGRバルブ開度が制御される。なお、遅れ要素212は、操作量upkを1サンプル遅らせてupk-1を出力する。このupk-1は、状態オブザーバ120に出力される。状態オブザーバ120は、1-13式に従った演算を実施する。具体的には、エンジン特性130の測定値ypk=xpkに対して-Cxk|k-1を加算してxpk-Cxk|k-1を算出する。さらにこれにLを乗じて、Lxpk-Cxk|k-1を算出する。一方、第2モデル予測制御部210の出力upk-1に対してBを乗じてBupk-1を算出し、xk|k-1にAを乗じてAxk|k-1を算出する。そして、Axk|k-1、Bupk-1及びLxpk-Cxk|k-1を加算して遅れ要素121に入力して、1サンプル遅延させてから、第2モデル予測制御部210に次のxk|k-1として出力する。次に、エンジン制御装置1000の詳細な処理内容について図13及び図14を用いて説明する。まず、エンジン制御装置1000は、k=0に設定する図13:ステップS31。また、エンジン制御装置1000は、初期設定を実施するステップS33。具体的には、エンジン制御装置1000は、k=0におけるエンジン回転数及び燃料噴射量の設定値に応じたプラント行列A及びB並びにむだ時間補償量τ=[τ11,τ12,τ21,τ22]を取得すると共に、重み行列Q及びR、予測期間N、初期値x0及びu0、時定数Tr等を取得する。なお、これらのデータについては、予めメモリなどに保持しておき、該当データを読み出す。さらに、エンジン制御装置1000は、初期計算として、参照軌道行列Tを算出し、設定するステップS35。さらに、エンジン制御装置1000は、第2モデル予測制御部210及び状態オブザーバ120において用いる行列Hc1=ΘcTQΘcR1及びΘcTQ、L、Ψcなどを算出し、設定するステップS37。さらに、Hc1ΘcTQも算出してしまっても良い。処理は端子Bを介して、図14の処理に移行する。図14の処理の説明に移行して、エンジン制御装置1000は、kを1インクリメントしステップS41、MAPセンサ及びMAFセンサから測定値ypkを取得し、さらに現時刻におけるエンジン回転数及び燃料噴射量の設定値を取得して当該エンジン回転数及び燃料噴射量の設定値に対応する目標値skを取得するステップS43。なお、エンジン制御装置1000は、計画器100から現時刻におけるエンジン回転数及び燃料噴射量の設定値に対応する基準値urefも取得する。そして、エンジン制御装置1000の第2モデル予測制御部210は、目標値ベクトルTk=Tskを算出するステップS45。また、エンジン制御装置1000の状態オブザーバ120は、1-13式に従ってx^k|k-1を算出するステップS47。そして、第2モデル予測制御部210は、応答ベクトルΨcxk|k-1を算出し、さらに誤差ベクトルEk=Tk-Ψcxk|k-1を算出するステップS49。そして、第2モデル予測制御部210の第2MPC主制御部211は、2-3式に従ってuoptkを算出するステップS51。その後、第2モデル予測制御部210は、uoptkと基準値urefkとを加算することで操作量の指令値を算出し、当該指令値をエンジン1のVNTノズル及びEGRバルブに出力するステップS53。また、第2モデル予測制御部210は、upk-1=uoptkと設定して、時刻k1についての処理で用いるようにメモリに保持しておくステップS55。そして、エンジン制御装置1000は、kが予め設定された最大値kmaxに達したか判断するステップS57。kがkmaxに達していない場合には、処理はステップS41に戻る。一方、kがkmaxに達した場合には、処理を終了する。例えば、ステップS31に戻って、初期設定から処理を開始する。以上のような処理を実施することで、むだ時間が存在する多入力多出力制御系MIMOにおいて、むだ時間要素の位相遅れに起因するシステムの応答遅れや不安定性を補償して安定化することができる。なお、図13のステップS33の内、システム行列A、B、重みQ、R、むだ時間補償量τの読み込み、および、ステップS37の定数行列の計算を、図14のステップS41とステップS43の間に入れることにより、毎回、運転条件が変わってシステム行列が変化する場合に対応することもできる。[実施の形態3]EGRバルブ及びVNTノズルは、開閉されるものであるから、全閉が0、全開が1又は100といったように、制約を有する。従って、このような制約を考慮して上で述べたような評価関数Jkを最小にする解を求めることになる。すなわち、モデル予測制御では、このような制約条件の下での最適化問題を解くことになる。より具体的には、操作量に対する制約条件は以下のように表される。ここで、a、bはEGRの下限値、上限値を、c、dはVNTの下限値、上限値を示す。3-1式を以下のような拡大系で表すことを考える。そうすると3-1式は以下のように変形される。そして3-3式は以下のように表現される。さらに、Ukの全ての要素について上記の制約が課される。すなわち、以下のように表される。これをまとめると以下のように表される。このように3-2式の形に変形される。一方、制御量に対する制約条件は、以下のように表される。なお、上限値及び下限値に同じa,b,c,dを使用しているが、上で述べたa,b,c,dとは異なる。3-9式は、以下のような拡大系で表すことを考える。さらに3-11式は以下のように表すこともできる。さらに、Xkの全ての要素について上記の制約が課される。すなわち、以下のように表される。従って、これをまとめると以下のように表される。3-14式にXk=ΨcxkΘcΔUkを代入すると以下のように表される。以上の制約条件をまとめると以下のように表される。但し、xkはそのままでは観測できない状態があるので、xk|k-1を用いる。すなわち、以下のように表される。3-18式を満たさないようなUkが算出された場合には、以下のような3-18式の制約条件を満たし且つ1-28式のJkを最小化するUkを見つけるという最適化問題を解くことになる。すなわち、最適化問題は以下のように表される。なお、3-18式に含まれる、EGRとVNTに関する制約のa、b、c、dの値は、フィードフォワード値urefkを加算した後の値が飽和しているか否かを判定できるような値に設定されているものとする。3-19式の解法としては、QP問題、ラグランジュ乗数法、アクティブセット法、内点法などが存在しているので、ここでは詳しく述べない。このような制約条件を考慮した場合のブロック線図の一例を図15に示す。具体的には、燃料噴射量の設定値とエンジン回転数の設定値とが入力されて、燃料噴射量の値及びエンジン回転数の値に対応付けてEGRバルブ開度の基準値及びVNTノズル開度の基準値の組み合わせと新気量MAFの目標値及び吸気圧MAPの目標値の組み合わせが登録されている計画器100から、燃料噴射量の設定値及びエンジン回転数の設定値に対応するEGRバルブ開度の基準値及びVNTノズル開度の基準値urefkと、新気量MAFの目標値及び吸気圧MAPの目標値skとを読み出す。そして、計画器100が出力した新気量MAFの目標値及び吸気圧MAPの目標値skに、参照軌道行列Tを乗じて2N次元の目標ベクトルTkを算出すると共に、状態オブザーバ120の出力xk|k-1にΨcを乗じて2N次元の応答ベクトルΨcxk|k-1を算出する。また、1サンプル前の操作量upk-1に対してΓcを乗じたΓcupkも算出する。そうすると、2N次元の誤差ベクトルEk=Tk-Ψcxk|k-1-Γcupkが得られる。この誤差ベクトルEkは、最適化処理部311に入力される。最適化処理部311は、1-32式に従ってΔUk=H1ΘcTQEkを算出する。ここで、最適化処理部311は、状態オブザーバ120からの出力xk|k-1とΔUkとを用いて3-18式が満たされているか判断する。3-18式が満たされていない場合には、最適化処理部311は、3-19式に表される最適化処理を実施して、ΔUoptkを算出する。一方、3-18式が満たされている場合には、既に計算されているΔUkをΔUoptkに設定する。そして、最適化処理部311は、M2x2NをΔUoptkを乗ずることによって、Δupkを算出して出力する。積分要素312は、Δupを基に例えばΔupkupk-1、upkを算出する。そして、この操作量upkは、計画器100が出力したフィードフォワード値urefkと加算されてエンジン特性130に出力される。すなわち、upkurefkによって、VNTノズル開度及びEGRバルブ開度が制御される。なお、遅れ要素313は、操作量upkを1サンプル遅らせてupk-1を出力する。このupk-1は、状態オブザーバ120にも出力される。状態オブザーバ120は、1-13式に従った演算を実施する。具体的には、エンジン特性130の測定値ypk=xpkに対して-Cxk|k-1を加算してxpk-Cxk|k-1を算出する。さらにこれにLを乗じて、Lxpk-Cxk|k-1を算出する。一方、第3モデル予測制御部310の出力upk-1に対してBを乗じてBupk-1を算出し、xk|k-1にAを乗じてAxk|k-1を算出する。そして、Axk|k-1、Bupk-1及びLxpk-Cxk|k-1を加算して遅れ要素121に入力して、1サンプル遅延させてから、第3モデル予測制御部310に次のxk|k-1として出力する。次に、エンジン制御装置1000の詳細な処理内容について図16乃至図18を用いて説明する。まず、エンジン制御装置1000は、k=0に設定する図16:ステップS61。また、エンジン制御装置1000は、初期設定を実施するステップS63。具体的には、エンジン制御装置1000は、k=0におけるエンジン回転数及び燃料噴射量の設定値に応じたプラント行列A及びB並びにむだ時間補償量τ=[τ11,τ12,τ21,τ22]を取得すると共に、重み行列Q及びR、予測期間N、初期値x0及びu0、時定数Tr等を取得する。なお、これらのデータについては、予めメモリなどに保持しておき、該当データを読み出す。また、エンジン制御装置1000は、メモリなどに格納されている制約条件上の説明におけるa,b,c及びdを2セット、具体的には、VNTノズル開度の上限値及び下限値、EGRバルブ開度の上限値及び下限値、xmafの上限値及び下限値、並びにxmapの上限値及び下限値を取得するステップS65。さらに、エンジン制御装置1000は、初期計算として、参照軌道行列Tを算出し、設定するステップS67。また、エンジン制御装置1000は、最適化処理部311で用いる定数行列H1=ΘTQΘR1及びΘTQ等、L、Ψc、Γc、制約拡大行列f、F、g及びGを算出し、設定するステップS69。H1ΘTQも算出してしまっても良い。処理は端子Cを介して、図17の処理に移行する。図17の処理の説明に移行して、エンジン制御装置1000は、kを1インクリメントしステップS71、MAPセンサ及びMAFセンサから測定値ypkを取得し、さらに現時刻におけるエンジン回転数及び燃料噴射量の設定値を取得して当該エンジン回転数及び燃料噴射量の設定値に対応する目標値skを取得するステップS73。なお、エンジン制御装置1000は、計画器100から現時刻におけるエンジン回転数及び燃料噴射量の設定値に対応する基準値urefも取得する。そして、第3モデル予測制御部310は、目標値ベクトルTk=Tskを算出するステップS75。また、エンジン制御装置1000の状態オブザーバ120は、1-13式に従ってx^k|k-1を算出するステップS77。そして、第3モデル予測制御部310は、応答ベクトルΨcxk|k-1を算出すると共に、Γcupk-1を算出して、誤差ベクトルEk=Tk-Ψcxk|k-1-Γcupk-1を算出するステップS79。そして、最適化処理部311は、1-32式に従ってΔUkを算出するステップS81。そして、最適化処理部311は、3-18式の制約条件を満たしているかを判断するステップS83。制約条件を満たしているΔUkが得られた場合には、ΔUkをΔUoptkに設定して、処理は端子Dを介して図18の処理に移行する。一方、制約条件を満たしているΔUkが得られなかった場合には、最適化処理部311は、3-19式に従って最適化処理を実施してΔUoptkを算出するステップS85。そして処理は端子Dを介して図18の処理に移行する。図18の処理の説明に移行して、第3モデル予測制御部310は、ΔUoptkに対してM2x2Nを乗ずることで、Δuoptkを算出するステップS87。その後、第3モデル予測制御部310は、ΔuoptkからΔuoptkuoptk-1によりuoptk=upkを算出し、さらに基準値urefkを加算することで操作量の指令値を算出し、当該指令値をエンジン1のVNTノズル及びEGRバルブに出力するステップS89。また、第3モデル予測制御部310は、upk-1=uoptkと設定して、時刻k1についての処理で用いるようにメモリに保持しておくステップS91。そして、エンジン制御装置1000は、kが予め設定された最大値kmaxに達したか判断するステップS93。kがkmaxに達していない場合には、処理はステップS71に戻る。一方、kがkmaxに達した場合には、処理を終了する。例えば、ステップS61に戻って、初期設定から処理を開始する。以上のような処理を実施することで、むだ時間が存在する多入力多出力制御系MIMOにおいて、制約条件を満たしつつ、むだ時間要素の位相遅れに起因するシステムの応答遅れや不安定性を補償して安定化することができる。なお、図16のステップS63の内、システム行列A、B、重みQ、R、むだ時間補償量τの読み込み、ステップS65の制約条件及びステップS69の定数行列の計算を、図17のステップS71とステップS73の間に入れることにより、毎回、運転条件が変わってシステム行列が変化する場合に対応することもできる。[実施の形態4]第2の実施の形態についても、第3の実施の形態のように、制約条件に従った形で制御を行うように変形することができる。第2の実施の形態の場合には、3-14式にXk=ΨcxkΘcUkを代入すると以下のように表される。以上の制約条件をまとめると以下のように表される。但し、xkはそのままでは観測できない状態があるので、xk|k-1を用いる。すなわち、以下のように表される。4-3式を満たさないようなUkが算出された場合には、以下のような4-3式の制約条件を満たし且つ2-1式の拡大系の評価式を最小化するUkを見つけるという最適化問題を解くことになる。すなわち、最適化問題は以下のように表される。なお、4-3式に含まれる、EGRとVNTに関する制約fのa、b、c、dの値は、フィードフォワード値urefkを加算した後の値が飽和しているか否かを判定できるような値に設定されているものとする。4-4式の解法としては、QP問題、ラグランジュ乗数法、アクティブセット法、内点法などが存在しているので、ここでは詳しく述べない。このような制約条件を考慮した場合のブロック線図の一例を図19に示す。具体的には、燃料噴射量の設定値とエンジン回転数の設定値とが入力されて、燃料噴射量の値及びエンジン回転数の値に対応付けてEGRバルブ開度の基準値及びVNTノズル開度の基準値の組み合わせと新気量MAFの目標値及び吸気圧MAPの目標値の組み合わせが登録されている計画器100から、燃料噴射量の設定値及びエンジン回転数の設定値に対応するEGRバルブ開度の基準値及びVNTノズル開度の基準値urefkと、新気量MAFの目標値及び吸気圧MAPの目標値skとを読み出す。そして、計画器100が出力した新気量MAFの目標値及び吸気圧MAPの目標値skに、参照軌道行列Tを乗じて2N次元の目標ベクトルTkを算出すると共に、状態オブザーバ120の出力xk|k-1にΨcを乗じて2N次元の応答ベクトルΨcxk|k-1を算出する。そうすると、2N次元の誤差ベクトルEk=Tk-Ψcxk|k-1が得られる。この誤差ベクトルEkは、第2最適化処理部411に入力される。第2最適化処理部411は、2-3式に従ってUk=H1ΘTQEkを算出する。ここで、第2最適化処理部411は、状態オブザーバ120からの出力xk|k-1とUkとを用いて4-3式が満たされているか判断する。4-3式が満たされていない場合には、第2最適化処理部411は、4-4式に表される最適化処理を実施して、Uoptkを算出する。一方、4-3式が満たされている場合には、既に計算されているUkをUoptkに設定する。そして、第2最適化処理部411は、M2x2NをUoptkを乗ずることによって、upkを算出して出力する。そして、この操作量upkは、計画器100が出力したフィードフォワード値urefkと加算されてエンジン特性130に出力される。すなわち、upkurefkによって、VNTノズル開度及びEGRバルブ開度が制御される。なお、遅れ要素412は、操作量upkを1サンプル遅らせてupk-1を出力する。このupk-1は、状態オブザーバ120に出力される。状態オブザーバ120は、1-13式に従った演算を実施する。具体的には、エンジン特性130の測定値ypk=xpkに対して-Cxk|k-1を加算してxpk-Cxk|k-1を算出する。さらにこれにLを乗じて、Lxpk-Cxk|k-1を算出する。一方、第4モデル予測制御部410の出力upk-1に対してBを乗じてBupk-1を算出し、xk|k-1にAを乗じてAxk|k-1を算出する。そして、Axk|k-1、Bupk-1及びLxpk-Cxk|k-1を加算して遅れ要素121に入力して、1サンプル遅延させてから、第4モデル予測制御部410に次のxk|k-1として出力する。次に、エンジン制御装置1000の詳細な処理内容について図20及び図21を用いて説明する。エンジン制御装置1000は、k=0に設定する図20:ステップS101。また、エンジン制御装置1000は、初期設定を実施するステップS103。具体的には、エンジン制御装置1000は、k=0におけるエンジン回転数及び燃料噴射量の設定値に応じたプラント行列A及びB並びにむだ時間補償量τ=[τ11,τ12,τ21,τ22]を取得すると共に、重み行列Q及びR、予測期間N、初期値x0及びu0、時定数Tr等を取得する。なお、これらのデータについては、予めメモリなどに保持しておき、該当データを読み出す。また、エンジン制御装置1000は、メモリなどに格納されている制約条件上の説明におけるa,b,c及びdを2セット分、具体的には、VNTノズル開度の上限値及び下限値、EGRバルブ開度の上限値及び下限値、xmafの上限値及び下限値、並びにxmapの上限値及び下限値を取得するステップS105。さらに、エンジン制御装置1000は、初期計算として、参照軌道行列Tを算出し、設定するステップS107。また、エンジン制御装置1000は、第2最適化処理部411で用いる定数行列H1=ΘTQΘR1及びΘTQ等、L、Ψc、制約拡大行列f、F、g及びGを算出し、設定するステップS109。H1ΘTQも算出してしまっても良い。処理は端子Fを介して、図21の処理に移行する。図21の処理の説明に移行して、エンジン制御装置1000は、kを1インクリメントしステップS111、MAPセンサ及びMAFセンサから測定値ypkを取得し、さらに現時刻におけるエンジン回転数及び燃料噴射量の設定値を取得して当該エンジン回転数及び燃料噴射量の設定値に対応する目標値skを取得するステップS113。なお、エンジン制御装置1000は、計画器100から現時刻におけるエンジン回転数及び燃料噴射量の設定値に対応する基準値urefも取得する。そして、エンジン制御装置1000の第4モデル予測制御部410は、目標値ベクトルTk=Tskを算出するステップS115。また、エンジン制御装置1000の状態オブザーバ120は、1-13式に従ってx^k|k-1を算出するステップS117。そして、第4モデル予測制御部410は、応答ベクトルΨcx^k|k-1を算出し、さらに誤差ベクトルEk=Tk-Ψcx^k|k-1を算出するステップS119。そして、第4モデル予測制御部410は、2-3式に従ってUkを算出するステップS121。第4モデル予測制御部410の第2最適化処理部411は、4-3式の制約条件を満たしているかを判断するステップS123。制約条件を満たしているUkが得られた場合には、UkをUoptkに設定して処理はステップS127に移行する。一方、制約条件を満たしているUkが得られなかった場合には、第2最適化処理部411は、4-4式に従って最適化処理を実施してUoptkを算出するステップS125。そして、第2最適化処理部411は、Uoptkに対してM2x2Nを乗ずることで、uoptkを算出し、第4モデル予測制御部410は、uoptkに基準値urefkを加算することで操作量の指令値を算出し、当該指令値をエンジン1のVNTノズル及びEGRバルブに出力するステップS127。また、第4モデル予測制御部310は、upk-1=uoptkと設定して、時刻k1についての処理で用いるようにメモリに保持しておくステップS129。そして、エンジン制御部1000は、kが予め設定された最大値kmaxに達したか判断するステップS131。kがkmaxに達していない場合には、処理はステップS111に戻る。一方、kがkmaxに達した場合には、処理を終了する。例えば、ステップS101に戻って、初期設定から処理を開始する。以上のような処理を実施することで、むだ時間が存在する多入力多出力制御系MIMOにおいて、制約条件を満たしつつ、むだ時間要素の位相遅れに起因するシステムの応答遅れや不安定性を補償して安定化することができる。なお、図20のステップS103の内、システム行列A、B、重みQ、R、むだ時間補償量τの読み込み、ステップS105の制約条件及びステップS109の定数行列の計算を、図21のステップS111とステップS113の間に入れることにより、毎回、運転条件が変わってシステム行列が変化する場合に対応することもできる。[実施の形態5]1-4式におけるむだ時間補償量τ=[τ11,τ12,τ21,τ22]、プラントの行列Ap並びにBp等は運転モードによって異なる。従って、運転モードの切替が検出された場合又は所定時間毎にこれらによって影響を受ける行列を計算し直すと、より効果的である。従って、図22に示すように行列設定部500を導入する。行列設定部500は、例えばエンジン回転数rpm及び燃料噴射量qで特定される運転モード毎にむだ時間補償量τの値が登録されている第1のテーブルと、エンジン回転数rpm及び燃料噴射量qで特定される運転モード毎にエンジン特性の行列Ap及びBpとが登録されている第2のテーブルとを格納するデータ格納部501を有しており、パディ近似表現行列生成部502と、拡大システム行列生成部503と、オブザーバゲイン算出部504と、定数行列生成部505とを有する。パディ近似表現行列生成部502は、エンジン回転数rpm及び燃料噴射量qに基づき第1のテーブルから読み出したむだ時間補償量τを用いて、パディ近似についての状態方程式で用いられる行列Ad及びBdを算出する。拡大システム行列生成部503は、エンジン回転数rpm及び燃料噴射量qに基づき第2のテーブルから読み出したエンジン特性の行列Ap及びBpと、パディ近似表現行列生成部502によって生成された行列Ad及びBdを用いて、拡大系の行列A及びBを生成する。オブザーバゲイン算出部504は、拡大系の行列A及びBを用いてオブザーバゲインLを算出する。定数行列生成部505は、行列A及びB等を用いて、Ψc、Γc及びΘc等を算出する。より具体的には、図23に示すような処理を実施する。まず、行列設定部500は、例えば計画器100からエンジン回転数rpm及び燃料噴射量qの設定値を取得する図23:ステップS141。そして、行列設定部500は、エンジン回転数rpm及び燃料噴射量qの設定値に応じてデータ格納部501における第1のテーブルから、該当するむだ時間補償量τを特定するステップS143。そして、パディ近似表現行列生成部502は、むだ時間補償量τから、むだ時間についての状態方程式における行列Ad及びBdを算出し、例えばメモリに格納するステップS145。さらに、拡大システム行列生成部503は、エンジン回転数rpm及び燃料噴射量qの設定値から、データ格納部501における第2のテーブルにおける、エンジンの状態方程式の行列Ap及びBpを特定するとともに、行列Ap及びBpと行列Ad及びBdとを用いて、合成状態方程式の行列A及び行列Bを算出し、例えばメモリに格納するステップS147。その後、オブザーバゲイン算出部504は、行列A及びBから、オブザーバゲイン行列Lを生成し、例えばメモリに格納するステップS149。そして、定数行列生成部505は、行列A及びBから定数行列Ψc、Θc、Γc及びその他のモデル予測制御に用いる定数行列を算出し、例えばメモリに格納するステップS151。なお、Γcは、第2及び第4の実施の形態では不要である。このようにすれば、モデル予測制御で用いられる行列を運転モードに応じた適切な行列に設定することができるようになる。なお、第1乃至第4の実施の形態との関係を、図24に一例を示す。図24は、行列設定部500を第1の実施の形態に係るエンジン制御装置1000に組み込んだ場合のブロック線図を示している。このように、計画器100からエンジン回転数rpm及び燃料噴射量qの設定値を取得して、定期的に又はこれらの設定値が基準以上に変動した場合など、図23に示したような処理を実施して、モデル予測制御部110のMPC主制御部111、T、Ψc及びΓc、状態オブザーバ120の行列L、A及びBを設定する。第1乃至第4の実施の形態では、kをkmaxに達するまで処理を繰り返すようになっているが、例えばエンジン回転数rpm及び燃料噴射量qの設定値が急激に変化した場合には、割り込みをかけて図23に示したような処理を実施して、行列を新たに設定するようにしても良い。図24と同様に、他の実施の形態のエンジン制御装置1000にも行列設定部500を導入することができる。以上本技術の実施の形態を説明したが本技術はこれに限定されるものではない。例えば、ブロック線図に従ったプログラムモジュール構成ではなくてもよい。また、処理フローについても、処理結果が変わらない限り処理順番を入れ替えたり、並列実施しても良い。なお、上で述べたエンジン制御装置1000は、コンピュータ装置であって、図25に示すように、RAMRandomAccessMemory2501とプロセッサ2503とROMReadOnlyMemory2507とセンサ群2515とがバス2519で接続されている。本実施の形態における処理を実施するための制御プログラム及び存在している場合にはオペレーティングシステムOS:OperatingSystemは、ROM2507に格納されており、プロセッサ2503により実行される際にはROM2507からRAM2501に読み出される。所定のプロセッサ2503は、センサ群吸気圧センサ及び新気量センサ。場合によっては燃料噴射量測定部及びエンジン回転数測定部など。を制御して、所定の測定値を取得する。また、処理途中のデータについては、RAM2501に格納される。なお、プロセッサ2503は、ROM2507を含む場合もあり、さらに、RAM2501を含む場合もある。本技術の実施の形態では、上で述べた処理を実施するための制御プログラムは、コンピュータ読み取り可能なリムーバブルディスクに格納されて頒布され、ROMライタによってROM2507に書き込まれる場合もある。このようなコンピュータ装置は、上で述べたプロセッサ2503、RAM2501、ROM2507などのハードウエアと制御プログラム場合によってはOSもとが有機的に協働することにより、上で述べたような各種機能を実現する。但し、エンジン制御装置全体をハードウエアのみにて実装することも可能である。以上述べた本実施の形態をまとめると以下のようになる。実施の形態に係るエンジン制御方法は、A排気循環器及び可変ノズルターボを有するエンジンに対する燃料噴射量の設定値、エンジン回転数の設定値、エンジンの吸気圧の測定値及び新気量の測定値を取得するステップと、B吸気圧の測定値及び新気量の測定値とエンジンが有するむだ時間特性をパディ近似により表すことで追加される吸気圧に関する内部状態及び新気量に関する内部状態とを表す状態ベクトルと第1の行列との積と、排気循環器のバルブ開度と可変ノズルターボのノズル開度とを表す操作量ベクトルと第2の行列との積との和により1単位時間後の状態ベクトルが表される合成状態方程式における第1の行列と第2の行列とを用いて吸気圧に関する内部状態及び新気量に関する内部状態を状態オブザーバによって推定し、第2の状態ベクトルを生成するステップと、C燃料噴射量の設定値及びエンジン回転数の設定値に応じた新気量の目標値及び吸気圧の目標値と、生成した第2の状態ベクトルと、第1の行列と第2の行列とから得られる行列とを用いてモデル予測制御を実施することで、排気循環器のバルブ開度の操作量と可変ノズルターボのノズル開度の操作量とを算出し、燃料噴射量の設定値及びエンジン回転数の設定値に応じた、排気循環器のバルブ開度の基準値及び可変ノズルターボのノズル開度の基準値を加算することで、排気循環器のバルブ開度の指令値及び可変ノズルターボのノズル開度の指令値を算出するステップとを含む。そして、第1の行列により、状態ベクトルに含まれる要素間の相互作用を表し、第2の行列が、排気循環器のバルブ開度が吸気圧に影響を与える干渉要素と、可変ノズルターボのノズル開度が新気量に影響を与える干渉要素とを含む、エンジンの状態方程式における行列と、排気循環器のバルブ開度が吸気圧に関する内部状態に影響を与える干渉要素と、可変ノズルターボのノズル開度が新気量に関する内部状態に影響を与える干渉要素とを含む、むだ時間特性をパディ近似で表すことで規定される状態方程式における行列とを含む。このように、状態オブザーバを導入すると共に、上記のような第1及び第2の行列を用いることで、吸気圧制御系及び新気量制御系で干渉があってもむだ時間を適切に補償するモデル予測制御を実施することができるようになる。なお、上で述べた第1の行列は、A吸気圧が新気量に影響を与える干渉要素と新気量が吸気圧に影響を与える干渉要素とを含む、エンジンの状態方程式における行列と、B吸気圧に関する内部状態が新気量に関する内部状態に影響を与える干渉要素と新気量に関する内部状態が吸気圧に関する内部状態に影響を与える干渉要素とを含む、むだ時間特性をパディ近似で表すことで規定される状態方程式における行列と、C吸気圧に関する内部状態及び新気量に関する内部状態が新気量に影響を与える干渉要素、及び吸気圧に関する内部状態及び新気量に関する内部状態が吸気圧に影響を与える干渉要素とを含むようにしてもよい。制御系の干渉に対処できるようになる。なお、燃料噴射量の設定値及びエンジン回転数の設定値に応じて第1の行列及び第2の行列を更新するステップをさらに含むようにしても良い。このようにすれば、運転モードの切り替えに適切に対処することができ、安定的な制御が可能となる。なお、第1の行列及び第2の行列が、燃料噴射量の設定値及びエンジン回転数の設定値に応じて変化するむだ時間の補償量に応じて算出されるような場合もある。また、上で述べたモデル予測制御の処理が、A排気循環器のバルブ開度の操作量の時系列値及び可変ノズルターボのノズル開度の操作量の時系列値が、排気循環器のバルブ開度及び可変ノイズターボのノズル開度の飽和による制約条件を満たすか判断するステップと、B制約条件が満たされないと判断された場合には、制約条件を満たしつつ、追従誤差に関する値と、操作量に関する値とを重み付け加算した評価関数の値を最小化するように、排気循環器のバルブ開度の操作量の時系列値及び可変ノズルターボのノズル開度の操作量の時系列値を計算し直すステップとを含む場合がある。このようにすれば、排気循環器のバルブ開度及び可変ノズルターボのノズル開度の飽和による制約条件を満たすように操作量を決定することができるようになる。なお、制約条件については、吸気圧及び新気量の測定値についての制約を含む場合もある。なお、実施の形態に係るエンジン制御装置図26:6000は、A排気循環器及び可変ノズルターボを有するエンジンに対する燃料噴射量の設定値、エンジン回転数の設定値、エンジンの吸気圧の測定値及び新気量の測定値を取得するデータ取得部図26:6010と、B吸気圧の測定値及び新気量の測定値とエンジンが有するむだ時間特性をパディ近似により表すことで追加される吸気圧に関する内部状態及び新気量に関する内部状態とを表す状態ベクトルと上で述べた第1の行列との積と、排気循環器のバルブ開度と可変ノズルターボのノズル開度とを表す操作量ベクトルと上で述べた第2の行列との積との和により1単位時間後の状態ベクトルが表される合成状態方程式における第1の行列と第2の行列とを用いて吸気圧に関する内部状態及び新気量に関する内部状態を推定し、第2の状態ベクトルを生成する状態オブザーバ図26:6020と、C燃料噴射量の設定値及びエンジン回転数の設定値に応じた新気量の目標値及び吸気圧の目標値と、生成した第2状態ベクトルと、第1の行列と第2の行列とから得られる行列とを用いてモデル予測制御を実施することで、排気循環器のバルブ開度の操作量と可変ノズルターボのノズル開度の操作量とを算出し、燃料噴射量の設定値及びエンジン回転数の設定値に応じた、排気循環器のバルブ開度の基準値及び可変ノズルターボのノズル開度の基準値を加算することで、排気循環器のバルブ開度の指令値及び可変ノズルターボのノズル開度の指令値を算出するモデル予測制御部図26:6030とを有する。なお、上記方法による処理をプロセッサに行わせるためのプログラムを作成することができ、当該プログラムは、例えばフレキシブルディスク、CD-ROM、光磁気ディスク、半導体メモリ、ハードディスク等のコンピュータ読み取り可能な記憶媒体又は記憶装置に格納される。尚、中間的な処理結果はメインメモリ等の記憶装置に一時保管される。以上の実施例を含む実施形態に関し、さらに以下の付記を開示する。付記1排気循環器及び可変ノズルターボを有するエンジンに対する燃料噴射量の設定値、エンジン回転数の設定値、前記エンジンの吸気圧の測定値及び新気量の測定値を取得するステップと、前記吸気圧の測定値及び前記新気量の測定値と前記エンジンが有するむだ時間特性をパディ近似により表すことで追加される前記吸気圧に関する内部状態及び前記新気量に関する内部状態とを表す状態ベクトルと第1の行列との積と、前記排気循環器のバルブ開度と前記可変ノズルターボのノズル開度とを表す操作量ベクトルと第2の行列との積との和により1単位時間後の状態ベクトルが表される合成状態方程式における前記第1の行列と前記第2の行列とを用いて前記吸気圧に関する内部状態及び前記新気量に関する内部状態を状態オブザーバによって推定し、前記エンジンの吸気圧の測定値及び新気量の測定値と併せて状態ベクトルを生成するステップと、前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じた前記新気量の目標値及び前記吸気圧の目標値と、生成した前記状態ベクトルとから、前記第1の行列と前記第2の行列とから得られる行列を用いてモデル予測制御を実施することで、前記排気循環器のバルブ開度の操作量と前記可変ノズルターボのノズル開度の操作量とを算出し、前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じた、前記排気循環器のバルブ開度の基準値及び前記可変ノズルターボのノズル開度の基準値を加算することで、前記排気循環器のバルブ開度の指令値及び前記可変ノズルターボのノズル開度の指令値を算出するステップと、をプロセッサに実行させ、前記第1の行列により、前記状態ベクトルに含まれる要素間の相互作用を表し、前記第2の行列が、前記排気循環器のバルブ開度が前記吸気圧に影響を与える干渉要素と、前記可変ノズルターボのノズル開度が前記新気量に影響を与える干渉要素とを含む、前記エンジンの状態方程式における行列と、前記排気循環器のバルブ開度が前記吸気圧に関する内部状態に影響を与える干渉要素と、前記可変ノズルターボのノズル開度が前記新気量に関する内部状態に影響を与える干渉要素とを含む、前記むだ時間特性をパディ近似で表すことで規定される状態方程式における行列とを含む、エンジン制御プログラム。付記2前記第1の行列が、前記吸気圧が前記新気量に影響を与える干渉要素と前記新気量が前記吸気圧に影響を与える干渉要素とを含む、前記エンジンの状態方程式における行列と、前記吸気圧に関する内部状態が前記新気量に関する内部状態に影響を与える干渉要素と前記新気量に関する内部状態が前記吸気圧に関する内部状態に影響を与える干渉要素とを含む、前記むだ時間特性をパディ近似で表すことで規定される状態方程式における行列と、前記吸気圧に関する内部状態及び前記新気量に関する内部状態が前記新気量に影響を与える干渉要素、及び前記吸気圧に関する内部状態及び前記新気量に関する内部状態が前記吸気圧に影響を与える干渉要素とを含む付記1記載のエンジン制御プログラム。付記3前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じて前記第1の行列及び前記第2の行列を更新するステップをさらに前記プロセッサに実行させるための付記1又は2記載のエンジン制御プログラム。付記4前記第1の行列及び前記第2の行列が、前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じて変化するむだ時間の補償量に応じて算出される付記3記載のエンジン制御プログラム。付記5前記モデル予測制御の処理が、前記排気循環器のバルブ開度の操作量の時系列値及び前記可変ノズルターボのノズル開度の操作量の時系列値が、前記排気循環器のバルブ開度及び前記可変ノイズターボのノズル開度の飽和による制約条件を満たすか判断するステップと、前記制約条件が満たされないと判断された場合には、前記制約条件を満たしつつ、追従誤差に関する値と、操作量に関する値とを重み付け加算した評価関数の値を最小化するように、前記排気循環器のバルブ開度の操作量の時系列値及び前記可変ノズルターボのノズル開度の操作量の時系列値を計算し直すステップと、を含む付記1記載のエンジン制御プログラム。付記6排気循環器及び可変ノズルターボを有するエンジンに対する燃料噴射量の設定値、エンジン回転数の設定値、前記エンジンの吸気圧の測定値及び新気量の測定値を取得するステップと、前記吸気圧の測定値及び前記新気量の測定値と前記エンジンが有するむだ時間特性をパディ近似により表すことで追加される前記吸気圧に関する内部状態及び前記新気量に関する内部状態とを表す状態ベクトルと第1の行列との積と、前記排気循環器のバルブ開度と前記可変ノズルターボのノズル開度とを表す操作量ベクトルと第2の行列との積との和により1単位時間後の状態ベクトルが表される合成状態方程式における前記第1の行列と前記第2の行列とを用いて前記吸気圧に関する内部状態及び前記新気量に関する内部状態を状態オブザーバによって推定し、第2の状態ベクトルを生成するステップと、前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じた前記新気量の目標値及び前記吸気圧の目標値と、生成した前記第2の状態ベクトルと、前記第1の行列と前記第2の行列とから得られる行列とを用いてモデル予測制御を実施することで、前記排気循環器のバルブ開度の操作量と前記可変ノズルターボのノズル開度の操作量とを算出し、前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じた、前記排気循環器のバルブ開度の基準値及び前記可変ノズルターボのノズル開度の基準値を加算することで、前記排気循環器のバルブ開度の指令値及び前記可変ノズルターボのノズル開度の指令値を算出するステップと、を含み、前記第1の行列により、前記状態ベクトルに含まれる要素間の相互作用を表し、前記第2の行列が、前記排気循環器のバルブ開度が前記吸気圧に影響を与える干渉要素と、前記可変ノズルターボのノズル開度が前記新気量に影響を与える干渉要素とを含む、前記エンジンの状態方程式における行列と、前記排気循環器のバルブ開度が前記吸気圧に関する内部状態に影響を与える干渉要素と、前記可変ノズルターボのノズル開度が前記新気量に関する内部状態に影響を与える干渉要素とを含む、前記むだ時間特性をパディ近似で表すことで規定される状態方程式における行列とを含む、エンジン制御方法。付記7排気循環器及び可変ノズルターボを有するエンジンに対する燃料噴射量の設定値、エンジン回転数の設定値、前記エンジンの吸気圧の測定値及び新気量の測定値を取得するデータ取得部と、前記吸気圧の測定値及び前記新気量の測定値と前記エンジンが有するむだ時間特性をパディ近似により表すことで追加される前記吸気圧に関する内部状態及び前記新気量に関する内部状態とを表す状態ベクトルと第1の行列との積と、前記排気循環器のバルブ開度と前記可変ノズルターボのノズル開度とを表す操作量ベクトルと第2の行列との積との和により1単位時間後の状態ベクトルが表される合成状態方程式における前記第1の行列と前記第2の行列とを用いて前記吸気圧に関する内部状態及び前記新気量に関する内部状態を推定し、第2の状態ベクトルを生成する状態オブザーバと、前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じた前記新気量の目標値及び前記吸気圧の目標値と、生成した前記第2の状態ベクトルと、前記第1の行列と前記第2の行列とから得られる行列とを用いてモデル予測制御を実施することで、前記排気循環器のバルブ開度の操作量と前記可変ノズルターボのノズル開度の操作量とを算出し、前記燃料噴射量の設定値及び前記エンジン回転数の設定値に応じた、前記排気循環器のバルブ開度の基準値及び前記可変ノズルターボのノズル開度の基準値を加算することで、前記排気循環器のバルブ開度の指令値及び前記可変ノズルターボのノズル開度の指令値を算出するモデル予測制御部と、を有し、前記第1の行列により、前記状態ベクトルに含まれる要素間の相互作用を表し、前記第2の行列が、前記排気循環器のバルブ開度が前記吸気圧に影響を与える干渉要素と、前記可変ノズルターボのノズル開度が前記新気量に影響を与える干渉要素とを含む、前記エンジンの状態方程式における行列と、前記排気循環器のバルブ開度が前記吸気圧に関する内部状態に影響を与える干渉要素と、前記可変ノズルターボのノズル開度が前記新気量に関する内部状態に影響を与える干渉要素とを含む、前記むだ時間特性をパディ近似で表すことで規定される状態方程式における行列とを含む、エンジン制御装置。"