JP3960106B2 - エンジンの空燃比制御装置 - Google Patents

Description

【０００１】
【発明の属する技術分野】
本発明はエンジンの空燃比制御装置、特にスライディングモード制御を空燃比のフィードバック制御に応用するものに関する。
【０００２】
【従来の技術】
エンジンの排気通路に設けた空燃比センサにより検出される検出空燃比が目標空燃比と一致するように空燃比のフィードバック制御を行う場合、古典制御のうち代表的なＰＩＤ制御が広く用いられている。
【０００３】
【発明が解決しようとする課題】
ところで、上記の古典制御であるＰＩＤ制御ではエンジンコントローラの入力に実際に観測される値を用いるため制御応答性及びロバスト性に限界があることから、最近になって現代制御の一手法であるスライディングモード制御をフィードバック制御に用いる方法が脚光を浴びている。
【０００４】
これを簡単に説明すると、これは図５に示したようにシステム（系）の誤差の状態量を非線形入力を与えてまず切換関数上に落ち着かせ、その後に線形入力により切換関数上を滑らせながら（スライディング状態）位相平面の原点へと収束させようとするものである。誤差の状態量が原点に収束したとき状態変数が目標値に収束する。スライディングモード制御の特徴として誤差の状態量が切換関数上にのってしまえばＰＩＤ制御を用いる場合よりもロバスト性が非常にあることが知られている。
【０００５】
しかしながら、スライディングモード制御についての公知文献（例えば野波健蔵他著『スライディングモード制御』コロナ社発行参照）をみても理論的な説明が主であり、空燃比のフィードバック制御に適用した具体例を記載してはいない。
【０００６】
従ってスライディングモード制御を空燃比フィードバック制御に適用できれば、ＰＩＤ制御を用いる場合よりもロバスト性を向上できることになる。
【０００７】
そこで本願の発明者はこの文献を参考にして空燃比フィードバック制御への適用を試みたわけであるが、この場合に本願の発明者は先にエンジンの排気系及び前記空燃比センサの動特性を離散系２次の伝達関数で記述できることを提案している（特開２０００−２９１４８４号公報参照）。
【０００８】
そこで本発明はこの離散系２次の伝達関数に基づいて離散系２次の状態方程式を作成し、この離散系２次の状態方程式に基づいてスライディングモード制御を行い空燃比操作量を演算することにより、スライディングモード制御を空燃比フィードバック制御に適用することを可能としてＰＩＤ制御を用いる場合よりも応答性とロバスト性をともに向上することを目的とする。
【０００９】
なお、スライディングモード制御を空燃比フィードバック制御に適用した先行技術として特開平８−２３２７１３号公報がある。この先行技術では制御対象となるプラント（エンジン）を連続系の物理モデルで表し、その連続系の物理モデルから出発してより離散化した状態方程式を求めている。これはいわば物理モデルに従って最初は連続系の方程式を立て、その連続系の方程式を離散系の方程式へと変更するものであり、最初から離散系で考えている本願発明とは技術的思想が異なっている。このため、この先行技術のように連続系の物理モデルより離散化した状態方程式を用いるのでは、物理モデルの精度によってはプラントを正確にモデル化するのが難しい。特にエンジンのような非線形性の強いものに対しては運転の各領域において正確に物理モデルを記述することは非常に困難である。具体的には先行技術ではプラントの遅れを連続系の二次遅れで近似するものに過ぎない。
【００１０】
また、目標空燃比が一定でなく変化する場合には目標値追従型であるサーボ系となるが、この目標値追従型であるサーボ系においては目標値への収束性をよくするため及び定常偏差を吸収するために入力を積分する機能が必要である。一般的なサーボ系問題では積分器を含んだ状態変数を追加し切換関数を設計する。この場合に切換関数によって算出される非線形入力も積分要素を含んでしまうため、応答性が悪化することが考えられる。しかしながら、目標空燃比が変化する場合を扱っている先行技術にはどの部分にも積分器が設けられていない。
【００１１】
【課題を解決するための手段】
第１の発明は、排気の空燃比を検出する空燃比検出手段と、エンジンの排気系及び前記空燃比検出手段の動特性を離散系２次の伝達関数で近似しシリンダ吸入空燃比のステップ応答より同定することによりこの離散系２次の伝達関数を設定する伝達関数設定手段と、この離散系２次の伝達関数を離散系２次の状態方程式に変換する手段と、この離散系２次の状態方程式に基づいて前記空燃比検出手段により検出される排気の空燃比が目標空燃比と一致するようにスライディングモード制御を用いた空燃比のフィードバック制御を行う空燃比フィードバック制御手段とを備え、前記離散系２次の状態方程式が
mat ｘ ( ｎ＋１ ) ＝ mat Ａ mat ｘ ( ｎ ) ＋ mat Ｂｕ ( ｎ )
mat ｙ ( ｎ ) ＝ mat Ｃ mat ｘ ( ｎ )
ただし、 mat ｘ ( ｎ＋１ ) ：時刻ｎ＋１の状態変数、
mat ｘ ( ｎ ) ：時刻ｎの状態変数、
mat ｙ ( ｎ ) ：時刻ｎの出力、
ｕ ( ｎ ) ：時刻ｎの入力、
mat Ａ ：システム行列、
mat Ｂ ：入力行列、
mat Ｃ ：出力行列、
である。
【００１２】
第２の発明は、第１の発明おいて前記同定を最小２乗法を用いて行う。
【００１３】
第３の発明は、第１の発明おいて前記同定をＡＲＸモデルを用いて行う。
【００１４】
第４の発明は、第１から第３までのいずれか一つの発明において前記空燃比のフィードバック制御が燃料供給量のフィードバック制御である。
【００１５】
第５の発明は、第４の発明において前記空燃比フィードバック制御手段が、前記離散系２次の状態方程式の状態変数matｘ(ｎ)と目標値の状態量matθ(ｎ)の差を誤差の状態量matｅ(ｎ)として演算する手段と、この誤差の状態量matｅ(ｎ)に切換関数ゲインmatＳを乗算して状態量σ(ｎ)を演算する手段と、この切換関数ゲインの乗算された状態量σ(ｎ)に基づいて非線形入力ｕnl(ｎ)を演算する手段と、同じくこの切換関数ゲインの乗算された状態量σ(ｎ)に基づいて線形入力ｕeq(ｎ)を演算する手段と、これら非線形入力と線形入力の和を総制御入力ｕsl(ｎ)として算出する手段と、この総制御入力に基づいて前記燃料供給量のフィードバック制御値ＡＬＰＨＡを演算する手段とからなる。
【００１６】
第６の発明では、第５の発明において目標空燃比が変化する場合に前記線形入力ｕeq(ｎ)に代えてこれを積分した値を用いる。
【００１７】
第７の発明では、第６の発明において排気通路に酸素ストレージ機能を有する触媒と、この触媒の前後にあって排気中の酸素濃度を検出する一対のセンサ（例えば触媒上流側が広域空燃比センサ、触媒下流側がＯ₂センサ）とを備える場合に、この一対のセンサの出力に基づいて前記触媒の酸素ストレージ量を最適化する空燃比を目標空燃比として演算する。
【００１８】
第８の発明では、第６の発明において前記積分する際の積分ゲインμと前記切換関数ゲインmatＳとの関係を前記線形入力の式から求め、前記目標空燃比への応答の時定数が最適となるように前記積分ゲインμを決定し、この決定した積分ゲインμより前記関係を用いて切換関数ゲインmatＳを算出する。
【００１９】
第９の発明では、第８の発明おいて前記積分ゲインμを回転速度に応じて定める（例えば回転速度が高くなるほど積分ゲインを大きくする）。
【００２１】
第１０の発明では、第１の発明において前記システム行列matＡの状態（例えば要素ａ₀の符号）により前記切換関数ゲインmatＳ（例えば要素Ｓ₂）の符号を変化させる。
【００２２】
第１１の発明では、第１の発明において前記離散系２次の伝達関数が
Ｇeng(ｑ)＝（ｂ₁ｑ＋ｂ₀）／（ｑ²＋ａ₁ｑ＋ａ₀）
（無駄時間の項を除く）
ただし、Ｇeng(ｑ)：プラントの伝達関数、
ｑ ：離散系シフトオペレータ、
ａ₁、ａ₀ ：微分係数、
ｂ₁、ｂ₀ ：微分係数、
である場合に、前記システム行列matＡ、入力行列matＢ、出力行列matＣが前記数１式により定義される値である。
【００２３】
第１２の発明では、第１の発明において前記離散系２次の伝達関数が
Ｇeng(ｑ)＝（ｂ₁ｑ＋ｂ₀）／（ｑ²＋ａ₁ｑ＋ａ₀）
（無駄時間の項を除く）
ただし、Ｇeng(ｑ)：プラントの伝達関数、
ｑ ：離散系シフトオペレータ、
ａ₁、ａ₀ ：微分係数、
ｂ₁、ｂ₀ ：微分係数、
である場合に、前記システム行列matＡ、入力行列matＢ、出力行列matＣが前記数２式により定義される値である。
【００２４】
第１３の発明では、第１２の発明において前記切換関数ゲインの乗算された状態量σ(ｎ)が
σ(ｎ)＝Ｓ（ｘ₁(ｎ)−θ₁(ｎ) ）＋（ｘ₁(ｎ)−ｘ₂(ｎ) ）
ただし、ｘ₁(ｎ)、ｘ₂(ｎ) ：状態変数matｘ(ｎ) の要素、
θ₁(ｎ) ：ｘ₁(ｎ)の目標値、
である。
【００２５】
【発明の効果】
第１、第２、第３、第４、第５、第１１の発明によれば、非線形性の強いエンジンのようなプラントであっても運転領域毎に空燃比ステップ応答を測定することにより運転領域毎の伝達関数を同定することが可能であり、この離散系２次の伝達関数を離散系２次の状態方程式に変換し、この離散系２次の状態方程式に基づいてスライディングモード制御を用いた空燃比のフィードバック制御を行うので、上記の先行技術のように連続系の物理モデルからスライディングモード制御に用いる状態方程式を算出する場合よりも精度良くプラントモデルを構築することができるとともに、ＰＩＤ制御を用いる場合より応答性とロバスト性が共に向上する。かつステップ応答測定より伝達関数が一意的に決まるのでスライディングモード制御の切換関数及び線形入力の設計が容易となる。
【００２６】
変化する目標値への追従性を良くし及び定常偏差を吸収するため目標値追従問題であるサーボ系では状態方程式に積分項を追加した拡大系を用いるのが一般的であるのでスライディングモード制御においても拡大系を用いたサーボ系としたとき、切換関数及び非線形入力に積分項を持つことになり応答性が低下するのであるが、第６、第７の発明によれば応答性を確保しつつ定常偏差の吸収が可能となる。
【００２７】
目標空燃比への収束を早めようと過大な空燃比操作量を与えたのでは運転性能や排気性能が悪化する可能性があるため各運転領域毎に最適な切換関数ゲイン及び積分ゲインを決定する必要がありこの双方を最適化するのは非常に工数がかかるのであるが、第８の発明によれば積分ゲインを決定しさえすれば、この積分ゲインより一定の関係式を用いて切換関数ゲインを算出することができるので、切換関数の設計が容易となる。
【００２８】
フィードバック制御を固定時間間隔で演算する場合に回転速度が上昇すると単位時間当たりの燃料噴射イベントが増加し見かけ上フィードバックゲインが落ちてしまい高回転速度域での目標値への追従性が悪化するのであるが、第９の発明によればフィードバック制御を固定周期で行った場合でも見かけ上のゲインの低下を防ぐことが可能となり、全運転領域で最適な積分ゲインを得ることができる。
【００２９】
第１０の発明によればシステム行列の状態によって切換関数ゲインの符号を切換えることでシステム行列の状態が変化しても空燃比操作量を正しい方向に付加することが可能となる。
【００３０】
第１２の発明によれば、状態量ｘ₁(ｎ)、ｘ₂(ｎ)が
ｘ₂(ｎ)＝ｘ₁(ｎ−１)
の関係を満たすことから、状態量ｘ₂(ｎ)についても物理的な意味が明確となった。すなわち、第１２の発明によれば、状態量ｘ₂(ｎ)は状態量ｘ₁(ｎ)の前回値を表す。
【００３１】
第１３の発明によれば、状態量ｘ₁(ｎ)の目標値θ₁(ｎ)との差分及び状態量ｘ₁(ｎ)の変化量がゼロとなるような切換関数を設定できるので、目標値θ₁(ｎ)の変化に対しても応答良く空燃比操作量を変化させることができる。
【００３２】
【発明の実施の形態】
以下、添付図面に基づき本発明の実施の形態について説明する。
【００３３】
図１はガソリンエンジンのシステム構成図である。同図において１はエンジン本体、２は吸気通路、３は排気通路、４はスロットル弁、５は燃料噴射弁、６は点火プラグである。
【００３４】
前記燃料噴射弁５からの燃料噴射量、燃料噴射時期、点火プラグ６による火花点火の時期を運転条件に応じて制御するためエンジンコントローラ１１を備える。エンジンコントローラ１１にはクランク角センサ（ポジションセンサ１２と位相センサ１３からなる）からの信号、エアフローメータ１４からの吸入空気流量の信号が、水温センサ１５からの信号などと共に入力し、これら信号に基づいて燃料噴射弁５からの燃料噴射量、噴射時期を制御し、また点火プラグ６による火花点火の時期を制御する。
【００３５】
エンジンの排気通路３には酸素ストレージ機能を有し触媒雰囲気が理論空燃比付近のときにＮＯｘ、ＨＣ、ＣＯを同時に浄化する三元触媒７を備える。この三元触媒７の酸素ストレージ量が飽和量に達したり酸素を全く保持しない状態となったりしてしまうとＨＣ、ＣＯ、ＮＯｘを効率よく浄化できなくなるため、エンジンコントローラ１１では触媒７の酸素ストレージ量を演算し、触媒７の転換効率を最大に保つべく触媒７の酸素ストレージ量が一定となるように触媒７の上流側に設けた広域空燃比センサ１６の出力と触媒７の下流側に設けたＯ₂センサ１７の出力とに基づいてエンジンの空燃比制御を行う。
【００３６】
この制御は概略次のようなものである（詳細は特願２０００−３８６８８号参照）。すなわちエンジンコントローラ１１は触媒７に流入する排気の空燃比とエンジンの吸入空気量に基づき触媒７の酸素ストレージ量を推定演算するが 、このとき酸素ストレージ量の演算を高速成分ＨＯ２と低速成分ＬＯ２とで分けて行う。具体的には酸素吸収時は高速成分ＨＯ２が優先して吸収し、高速成分ＨＯ２が吸収しきれない状態となったら低速成分ＬＯ２が吸収し始めるとして演算を行い、また酸素放出時は低速成分ＬＯ２と高速成分ＨＯ２の比（ＬＯ２／ＨＯ２）が一定割合ＡＲ以下の場合は高速成分ＨＯ２から優先して酸素が放出されるとし、比ＬＯ２／ＨＯ２が一定割合になったらその比ＬＯ２／ＨＯ２を保つように低速成分ＬＯ２と高速成分ＨＯ２の両方から酸素が放出されるとして酸素ストレージ量の演算を行う。
【００３７】
そして、演算された酸素ストレージ量の高速成分ＨＯ２が目標値よりも多いときは、エンジンの空燃比をリッチ側に制御して高速成分ＨＯ２を減少させ、目標値よりも少ないときは空燃比をリーン側に制御して高速成分ＨＯ２を増大させる。
【００３８】
この結果、酸素ストレージ量の高速成分ＨＯ２が目標とする値に保たれるので、触媒７に流入する排気の空燃比が理論空燃比からずれたとしても、応答性の高い高速成分ＨＯ２から直ちに酸素が吸収あるいは放出されて触媒雰囲気が理論空燃比方向に修正され、触媒７の転換効率が最大に保たれる。
【００３９】
さらに、演算誤差が累積すると演算される酸素ストレージ量が実際の酸素ストレージ量とずれてくるが、触媒７下流がリッチあるいはリーンになったタイミングで酸素ストレージ量（高速成分ＨＯ２及び低速成分ＬＯ２）のリセットを行うことで演算値と実際の酸素ストレージ量とのずれを修正している。
【００４０】
このように酸素ストレージ量に応じて目標空燃比を変化させ、こうして変化する目標空燃比が得られるように空燃比のフィードバック制御を行うのであるが、本発明ではこの空燃比フィードバック制御にスライディングモード制御を適用する。
【００４１】
これを図２により説明すると、同図はエンジンコントローラ１１の制御ブロック図である。同図において目標空燃比演算部２１では広域空燃比センサ１６及びＯ₂センサ１７の出力に基づいて前述のように触媒７の酸素ストレージ量を最適化する空燃比を目標空燃比ＴＧＡＢＦとして演算する。
【００４２】
理論空燃比の運転時にこの目標空燃比が得られるように従来はＰＩＤ制御により空燃比フィードバック補正値ＡＬＰＨＡを演算していたが、これに代えてスライディングモードコントローラ（スライディングモード制御部）２２を設けている。このスライディングモードコントローラ２２が後述するようにして空燃比フィードバック補正値ＡＬＰＨＡを演算する。そして燃料噴射量演算部３１では基本噴射パルス幅ＴＰ０から得られるシリンダ吸気量パルス幅ＴＰに対してこの空燃比フィードバック補正値ＡＬＰＨＡとこれ以外の各種補正を行って次の式により燃料噴射パルス幅ＣＴＩを演算する。そしてこの燃料噴射パルス幅ＣＴＩを用いて燃料噴射弁５を間欠的に駆動する。
【００４３】
ＣＴＩ＝（ＴＰ×ＴＦＢＹＡ＋ＫＡＴＯＨＯＳ）
×（ＡＬＰＨＡ＋ＫＢＬＲＣ−１）＋ＴＳ＋ＣＨＯＳ…（補１）
ただし、ＴＦＢＹＡ ：目標当量比、
ＫＡＴＯＨＯＳ：燃料フィードフォワード補正値、
ＡＬＰＨＡ ：空燃比フィードバック補正値、
ＫＢＬＲＣ ：空燃比学習値、
ＴＳ ：無効噴射パルス幅、
ＣＨＯＳ ：気筒別燃料フィードフォワード補正値、
スライディングモードコントローラ２２は切換関数演算部２３、非線形入力演算部２４、線形入力演算部２５、積分器２６、加算器２７、換算部２８及び補正制限部２９からなっている。スライディングモードコントローラ２２により実行されるこれらの制御内容は後述するフローチャート（図７参照）によりまとめて説明する。従ってここでは図２に示す個々のブロック２３〜２９の働きの説明は省略する。
【００４４】
次にスライディングモード制御を用いた空燃比フィードバック制御の考え方を詳述する。これは本願の発明者が新たに創作したものである。
【００４５】
図３はプラント（エンジン）の物理モデルである。シリンダ内の空燃比は燃料量のフィードフォワード制御により（例えば上記（補１）式のＫＡＴＯＨＯＳ、ＣＨＯＳによる）目標空燃比となることが補償されているため排気のダイナミクス（ガス混合ダイナミクスを含む）及び空燃比センサ１６のダイナミクスを統合し離散系２次で記述する。プラント（エンジン）の入力をシリンダ内空燃比、出力を検出空燃比とすると、エンジンの排気系及びセンサ１６の動特性は次のように離散系２次の伝達関数Ｇeng(ｑ)で表すことができる。
【００４６】
Ｇeng(ｑ)＝（ｂ₁ｑ＋ｂ₀）／（ｑ²＋ａ₁ｑ＋ａ₀） …（１）
（無駄時間の項を除く）
ただし、Ｇeng(ｑ)：プラントの伝達関数、
ｑ ：離散系シフトオペレータ、
ａ₁、ａ₀ ：微分係数、
ｂ₁、ｂ₀ ：微分係数、
（１）式の伝達関数Ｇeng(ｑ)は前記プラントの図４のような入出力信号ステップ応答よりＡＲＸモデル（或いは最小２乗法）を用いて同定することができる。ＡＲＸモデルによる同定手法は一般的に知られており（例えば特開２０００−２９１４８４号公報参照）、ここでは説明を省略する。なお、非線形性の強いエンジンのようなプラントにあっても運転領域毎に空燃比ステップ応答を測定することにより運転領域毎の伝達関数を同定することが可能である。
【００４７】
離散系２次の伝達関数Ｇeng(ｑ)を次のように離散系２次の状態空間表現（状態方程式）で記述する。
【００４８】
【数３】

【００４９】
ここで、入力ｕ(ｎ)は燃料噴射パルス幅より算出したシリンダ吸入空燃比、出力ｙ(ｎ)は広域空燃比センサ１６により検出される空燃比である。
【００５０】
なお、（２ａ）、（２ｂ）の状態方程式に使われる状態変数などは行列であるため記号を太字で記載するのが一般的である。このため例えばｘ(ｎ)が行列（matrix）であるときは「matｘ(ｎ)」のように「mat」を記号の前に付けて行列であることを表現している。またｎは現時刻を、ｎ＋１は次の制御周期の時刻を表している。
【００５１】
一般的なシステムの状態方程式は入力の数をｍ、出力の数をｌ（ｍ、ｌは任意の正の整数）として記載されるが、ここでは物理量として把握しやすいよう（２ａ）、（２ｂ）のように入力を２つ（ｍ＝２）、出力を２つ（ｌ＝２）として状態方程式を記述している。このとき状態変数matｘ(ｎ)はｘ₁(ｎ)、ｘ₂(ｎ)を２要素とする列ベクトル、また出力matｙ(ｎ)はｙ₁(ｎ)、ｙ₂(ｎ)を２要素とする列ベクトルでもある。実質的には１出力となるように出力行列matＣを定めている。すなわち（２ｂ）式よりｙ₁(ｎ)＝ｘ₁(ｎ)であるのに対してｙ₂(ｎ)＝０である。ここでｙ₁(ｎ)には検出空燃比ＡＦＳＡＦを採用する。従って、まとめると次のようになる。
【００５２】
ｙ₁(ｎ)＝ｘ₁(ｎ)＝ＡＦＳＡＦ …（補２）
ｙ₂(ｎ)＝０ …（補３）
また、（２ａ）式においてシステム行列matＡの要素であるａ₁、ａ₀、入力行列matＢの要素であるｂ₁、ｂ₀は（１）式の伝達関数を決定する際に定まっている。 さて、離散系で同定されたシステムに対して、スライディングモード制御を用いて状態変数matｘ(ｎ)を状態空間内のある軌道に追従させるサーボ系問題を考える。この離散系モデルの誤差の状態量matｅ(ｎ)を次のように状態変数matｘ(ｎ)と軌道である目標値の状態量matθ(ｎ)の差で与える。
【００５３】
【数４】

【００５４】
ここで（３）式の目標値θ₁(ｎ)はｘ₁(ｎ)（＝ＡＦＳＡＦ）に対する目標値であるから、θ₁(ｎ)は目標空燃比ＴＧＡＢＦである。
【００５５】
また同離散系モデルの切換関数を次式で定義する。
【００５６】
【数５】

【００５７】
図５はスライディングモード制御における誤差の状態量matｅ(ｎ)の変化を位相平面上に示したものである。誤差の状態量matｅ(ｎ)は、図中の切換関数までの距離に応じた非線形入力の分で切換関数上に拘束された後、線形入力である等価制御入力により切換関数上を滑る（スライディング状態）ようにして位相平面の原点へと収束する。これにより状態変数matｘ(ｎ)を目標値の状態量matθ(ｎ)へと収束させることができる。なお、Ｓ₂が固定の場合、切換関数ゲインmatＳの要素であるＳ₁を大きくすれば直線の傾きが急になる（制御ゲインが大きくなる）。
【００５８】
なお、状態変数matｘ(ｎ)の要素が２つであるため切換関数が位相平面上で直線として表され物理的に把握しやすいものとなっているが、状態変数matｘ(ｎ)の要素が３つになると切換関数は位相空間上で平面となり、さらに要素の数が増すと切換関数は超平面という概念になる。
【００５９】
次に線形入力及び非線形入力の具体的な操作量の演算方法を説明する。
【００６０】
システムがスライディングモードにある条件(連続系ではσ＝０、σドット＝０、ただし「σドット」はσの微分値である。）は離散系の場合、
σ(ｎ)＝０ …（５ａ）
Δσ(ｎ＋１)＝σ(ｎ＋１)−σ(ｎ)＝０ …（５ｂ）
の条件を満たす。すると（５ｂ）式は次のようになる。
【００６１】
Δσ(ｎ＋１)
＝σ(ｎ＋１)−σ(ｎ)
＝matＳ｛matｘ(ｎ＋１)−matθ(ｎ＋１)｝−σ(ｎ)
＝matＳ｛matＡmatｘ(ｎ)＋matＢｕeq(ｎ)−matθ(ｎ＋１)｝
−σ(ｎ) …（６）
従ってこの（６）式においてdet(matＳmatＢ)^-1≠０であれば線形入力（等価制御入力）は次の式となる。
【００６２】
matｕeq(ｎ)＝−（matＳmatＢ）^-1matＳ
×｛matＡmatｘ(ｎ)−matθ(ｎ＋１)｝−σ(ｎ）…（７）
（２ａ）式、（４）式よりスカラーで記述すると、このシステムの線形入力は次のようになる。
【００６３】
ｕeq(ｎ）＝（Ｓ₁ｂ₁＋Ｓ₂ｂ₀）^-1（ａ₁Ｓ₁＋ａ₀Ｓ₂＋Ｓ₁ ²／Ｓ₂）
×｛ｘ₁(ｎ）−θ₁(ｎ)｝ …（８）
次に非線形入力は任意に与えることができるので、ここでは次のように定義する。すなわち外乱に対する分もまとめて１つの関数とすることもできるのであるが、ここでは外乱に対する分をα(ｎ)とは別にβ(ｎ)として導入している。
【００６４】
ｕnl(ｎ）＝｛α(ｎ)＋β(ｎ)｝sgn［σ(ｎ)］ …（９）
ただし、β(ｎ）≧Ｈｍａｘ、
α(ｎ）、β(ｎ)：正のスカラー値関数、
Ｈｍａｘ ：外乱の最大推定値、
ここで（９）式の符号関数sgn［σ(ｎ)］は状態量σ(ｎ)＞０のときsgn［σ(ｎ)］＝１、σ(ｎ)＜０のときsgn［σ(ｎ)］＝−１となる関数である。この符号関数は例えば誤差の状態量が図５において切換関数を横切って下側にきたとき次回の制御時に切換関数の上側へと誤差の状態量を戻すために必要となるものである。
【００６５】
ここでチャタリングを防止するための非線形入力の条件は次のようになる。
【００６６】
（i）σ(ｎ）＞０のとき：
０≦σ(ｎ＋１)≦σ(ｎ) …（１０ａ）
（ii）σ(ｎ）＜０のとき：
σ(ｎ)≦σ(ｎ＋１)≦０ …（１０ｂ）
これは（１０ａ）の条件式についていえば現時刻ｎのときよりその次の時刻であるｎ＋１でのほうが誤差の状態量が収束点である０に近い位置にあることを表したものである。
【００６７】
（１０ａ）の条件を満足するスカラー値関数α(ｎ）を次にようにして算出す る。（４）式から、
σ(ｎ＋１)＝matＳmatｅ(ｎ＋１) …（１１）
とする。（２ａ）、（２ｂ）の状態方程式で表されたシステムに対して、状態量σ(ｎ）＝matＳmatｅ(ｎ)とするとき、σ(ｎ）＝σ(ｎ＋１)＝σ(ｎ＋２)…を満たす入力は（８）式の線形入力であり、このとき（１１）式は次のようになる。
【００６８】
σ(ｎ＋１)
＝matＳ｛matｘ(ｎ＋１)−matθ(ｎ＋１)｝
＝matＳ｛matＡmatｘ(ｎ）＋matＢ｛ｕnl(ｎ）＋ｕeq(ｎ)｝
−matＡmatθ(ｎ)｝
…
＝σ(ｎ)＋matＳmatＢｕnl(ｎ) …（１２）
（１２）式より（１０ａ）式を満たす条件は次のようになる。
【００６９】
‖matＳmatＢ‖α(ｎ)≦‖σ(ｎ）‖ …（１３）
これより（１０ａ）式を満たすスカラー値関数α(ｎ）は、
α(ｎ)＝η‖matＳmatＢ‖^-1‖σ(ｎ)‖ …（１４）
ただし、η：非線形ゲイン（０＜η＜１）、
σ(ｎ＋１)σ(ｎ)＞０、
よって非線形入力ｕnl(ｎ)は、
ｕnl(ｎ)
＝｛η‖matＳmatＢ‖^-1‖σ(ｎ)‖＋β(ｎ)｝sgn［σ(ｎ)］
＝｛η｜Ｓ₁ｂ₁＋Ｓ₂ｂ₀｜^-1｜σ(ｎ)｜＋β(ｎ)｝sgn［σ(ｎ)］…（１５）
ただし、０＜η＜１、σ(ｎ＋１)σ(ｎ)＞０、
となり、この（１５）式によればチャタリングを防止しかつ外乱に対してもロバストとなる。
【００７０】
このようにして（８）式により線形入力ｕeq(ｎ)を、また（１５）式により非線形入力ｕnl(ｎ）が演算したので総制御入力である空燃比操作量ｕsl(ｎ)を次式により算出する。
【００７１】
ｕsl(ｎ)＝κ∫ｕeq(ｎ)ｄｔ＋ｕnl(ｎ） …（１６）
ただし、κ：積分ゲイン、
（１６）式右辺第１項では線形入力ｕeq(ｎ)を積分しており、そのための積分器２６を備えるが（図２参照）、これは可変値としての目標空燃比への収束性の向上及び定常偏差の吸収のため導入したものである。
【００７２】
図６のフローチャートはスライディングモード制御を設計するためのもので、これは開発段階で行う。
【００７３】
ステップ１でシリンダ内空燃比のステップ応答を測定し（図４参照）、ステップ２では離散系２次の伝達関数Ｇeng(ｑ)をＡＲＸモデル（あるいは最小２乗法）を用いて同定する。ステップ３ではこの伝達関数Ｇeng(ｑ)を（２ａ）、（２ｂ）式で示した離散系２次の状態方程式に変換する。ステップ４では切換関数σ(ｎ)＝０を設計する。ステップ５、６、７では線形入力ｕeq(ｎ)の演算式である（８）式、非線形入力ｕnl(ｎ)の演算式である（１５）式、空燃比操作量ｕsl(ｎ)の演算式である（１６）式をそれぞれ設計する。
【００７４】
次に一般的なスライディングモード制御のサーボ系問題を記述する。
【００７５】
仮に
【００７６】
【数６】

【００７７】
のような正準系のシステムの場合、（１７）式の状態変数に目標値θ_rと出力θ_yの差を積分した値であるｚを付加した次の拡大系を用いる。
【００７８】
【数７】

【００７９】
拡大系では切換関数を次式で定義する。
【００８０】
σ(ｎ）＝matＳmatｘ＝Ｓ₁ｚ＋Ｓ₂ｘ₁＋ｘ₂＝０ …（１９）
拡大系を用いた場合、（１９）式中に積分項ｚがあるため、状態が変化した場合に本来応答性を確保するための非線形入力（（１５）式）の操作が遅くなり、外乱ロバスト性も低下する。しかしながら、（１６）式のように積分項を線形入力に持たせることで、拡大系を用いた場合であっても応答性を低下させることなく目標空燃比へと収束させることができる。
【００８１】
図７のフローチャートは空燃比フィードバック補正値ＡＬＰＨＡを演算するためのもので、理論空燃比での運転中に一定時間毎またはクランク角の基準位置信号（Ｒｅｆ）の入力毎に実行する。これは図２のスライディングモードコントローラ２２の各ブロック２３〜２９をフローチャートで再構成したものである。
【００８２】
ステップ１１、１２では目標空燃比ＴＧＡＢＦ、実空燃比ＡＦＳＡＦを読み込み、ステップ１３で状態変数matｘの要素である状態量ｘ₁、ｘ₂を算出する。すなわち（２ｂ）式よりｙ₁＝ｘ₁、ｙ₂＝０である。ここでｙ₁は検出空燃比ＡＦＳＡＦであるからｘ₁＝ＡＦＳＡＦである。また（１７）式は変数が２つ（ｘ₁、ｘ₂）の連立２元方程式であるから、これをｘ₁、ｘ₂について解けばｘ₂がｘ₁の関数として得られるので、これにｘ₁＝ＡＦＳＡＦを代入してやればｘ₂が得られる。
【００８３】
ステップ１４では状態変数matｘ(ｎ)、目標値の状態量matθ及び切換関数ゲインmatＳを用いて現時刻ｎにおける切換関数σ(ｎ)＝０を算出する。ステップ１４は図２の切換関数演算部２３に相当する。
【００８４】
ステップ１５では状態量σ(ｎ)、非線形ゲインη、切換関数ゲインmatＳ、入力行列matＢ、スカラー値関数β(ｎ)、符号関数sgn［σ(ｎ)］を用いて（１５）式により現時刻ｎにおける非線形入力ｕnl(ｎ)を算出する。ステップ５は図２の非線形入力演算部２４に相当する。
【００８５】
ステップ１６では状態量ｘ₁(ｎ)、その目標値θ₁(ｎ)、切換関数ゲインmatＳ、入力行列matＢ、システム行列matＡを用いて（８）式により現時刻ｎにおける線形入力ｕeq(ｎ)を算出する。ステップ６は図２の線形入力演算部２５に相当する。
【００８６】
ステップ１７ではこの線形入力ｕeq(ｎ)を積分した後、ステップ１８で非線形入力ｕnl(ｎ)とこの積分した線形入力とを加算して現時刻ｎにおける空燃比操作量ｕsl(ｎ)を算出する。ステップ１７、１８は図２の積分器２６、加算器２７に相当する。
【００８７】
ステップ１９では空燃比操作量ｕsl(ｎ)を
ＡＬＰＨＡ＝ＣＹＬＡＦ／（ＣＹＬＡＦ＋ｕsl(ｎ)）×１００…（補４）
ただし、ＣＹＬＡＦ：シリンダ吸入空燃比、
の式により空燃比フィードバック補正値ＡＬＰＨＡに換算し（単位が空燃比から%に変換される）、ステップ２０でリミッタ処理を実行する。ステップ１９、２０は図２の換算部２８、補正制限部２９に相当する。
【００８８】
なお、（補４）式のシリンダ吸入空燃比ＣＹＬＡＦは次の式により演算される値である。
【００８９】
ＣＹＬＡＦ＝１４．７×ＴＰ
／｛ＴＰ×ＴＦＢＹＡ×（ＡＬＰＨＡ＋ＫＢＬＲＣ−１）｝…（補５）
次にフィードバックを行う場合の最適な切換関数ゲインmatＳの要素Ｓ₁を算出するための積分ゲインの設定方法を説明する。なお、切換関数ゲインのもう一つの要素Ｓ₂については後述する考察によりＳ₂＝１または−１となるので、一方の要素Ｓ₁についてだけ考えておけばよいことになる。
【００９０】
空燃比フィードバック制御の場合、シリンダ内の空燃比が急変すると運転性や排気性能に悪い影響を及ぼすので空燃比操作量が大きくなりすぎないように限界値で制限する必要がある。その一方で空燃比操作量を小さくすると目標空燃比への収束性が悪くなる。従って目標値への収束性と排気、運転性能とはトレードオフの関係を有する。そこで目標値への応答の時定数が所望の値となるときの積分ゲインをμとすると線形入力の（８）式及び空燃比操作量の（１６）式右辺の積分項からμは次のように置くことができる。
【００９１】
μ＝κ（Ｓ₁ｂ₁＋Ｓ₂ｂ₀）^-1（ａ₁Ｓ₁＋ａ₀Ｓ₂＋Ｓ₁ ²／Ｓ₂） …（２０）
ここで、線形入力の積分ゲインκを１とすると図５において切換関数である直線の傾きを表すＳ₁は次式となり、所望の応答を実現するＳ₁を演算できる。
【００９２】
Ｓ₁ ＝
｛−Ｓ₂（ａ₁−μｂ₁）±Ｓ₂｛（ａ₁−μｂ₁）²−（ａ₀−μｂ₀）｝^1/2｝
／２ …（２１）
なお、（２１）式のＳ₁は２値を採りうるので安定性判別を行って安定性のよいほうを選択する。
【００９３】
図８のフローチャートは最適な切換関数ゲイン（Ｓ₁）を設定するためのものである。なおこれは図６のスライディングモード制御設計フローチャートにおけるステップ４の内容の一部をなすものである。
【００９４】
ステップ２１で所望の応答が得られる空燃比操作量となる積分ゲインμを設定し（μは経験により適当な値がわかっている）、ステップ２２ではこの積分ゲインμに基づき上記の（２１）式より切換関数ゲインmatＳの要素Ｓ₁を算出する。なお、（２１）式によりＳ₁を算出する際にＳ₂が定まっている必要があるが、これは後述するようにａ₀＞０のとき１に、またａ₀＜０のとき−１になる値である。
【００９５】
次に積分ゲインμを一定として空燃比のフィードバック制御を行った場合のＡＬＰＨＡの変化を図９、図１０に示すと、図９は低回転速度時の、これに対して図１０は高回転速度時のものである。ＡＬＰＨＡを固定時間間隔（例えば１０ｍｓ毎のジョブ）で演算する場合、回転速度の上昇に伴い単位時間当たりの噴射回数が増加するため１噴射当りのＡＬＰＨＡの変化代が例えば５％から２％へと減少する。これにより低回転速度においてマッチングした積分ゲインμを高回転速度域においてもそのまま使用すると、ＡＬＰＨＡの変化代が減少して応答性が悪化する。これを回避するためには積分ゲインμを回転速度に応じて変化させることが必要であり、このとき図１１のように高回転速度域でも低回転速度時と同じ応答性を確保することが可能となる（１噴射当りのＡＬＰＨＡの変化代が５％へと回復している）。
【００９６】
図１２のフローチャートは積分ゲインμを演算するためのものである。
【００９７】
ステップ３１では回転速度を読み込み、ステップ３２でこの回転速度に応じて積分ゲインの回転速度補正係数ＫＮＥを算出する。ステップ３３では積分ゲイン基本値μ０にこの補正係数ＫＮＥを乗算した値を積分ゲインμとして算出する。ここで積分ゲイン基本値μ０は最低の回転速度時にマッチングした値である。このようにして演算された積分ゲインμを用いて上記の（２１）式によりＳ₁が演算され、このＳ₁が図７の空燃比フィードバック補正値演算フローチャートにおけるステップ１５、１６で非線形入力、線形入力の算出に用いられる。
【００９８】
次に（４）式の状態量σ(ｎ)を（３）式によって変形する。
【００９９】
σ(ｎ)＝matＳmatｅ(ｎ)
＝Ｓ₁（ｘ₁(ｎ)−θ₁）＋Ｓ₂（ｘ₂(ｎ)−θ₂）
＝Ｓ₁（ｘ₁(ｎ)−θ₁）−Ｓ₂ａ₀（ｘ₁(ｎ−１)−θ₁）…（２２）
誤差の状態量matｅ(ｎ)が切換関数をよぎる前に誤差の状態量の符号が変化するためには（２２）式右辺の第１項と第２項の符号は逆である必要があるので次の条件を仮定する。
【０１００】
（i）ａ₀＞０のとき ：Ｓ₁＝Ｓ＞０、Ｓ₂＝１ …（補６）
（ii）ａ₀＜０のとき：Ｓ₁＝Ｓ＞０、Ｓ₂＝−１ …（補７）
このときの状態量σ(ｎ)を成分で記述すると次のようになる。
【０１０１】
（i）α₀＞０のとき：
σ(ｎ)＝Ｓ｛ｘ₁(ｎ)−θ₁｝＋｛ｘ₂(ｎ)−θ₂｝ …（２３）
（ii）α₀＜０のとき：
σ(ｎ)＝Ｓ｛ｘ₁(ｎ)−θ₁｝−｛ｘ₂(ｎ）−θ₂｝ …（２４）
また（補６）、（補７）式を用いると、（８）式の線形入力、（１５）式の非線形入力及び（１６）式の総制御入力（空燃比操作量）は次の（２５）式〜（３０）式のように書き直される。
〈１〉線形入力（等価制御入力）：
（i）ａ₀＞０のとき：
ｕeq(ｎ)＝（Ｓｂ₁＋ｂ₀）^-1（ａ₁Ｓ＋ａ₀＋１）×｛ｘ₁(ｎ)−θ₁(ｎ)｝…（２５）
（ii）ａ₀＜０のとき：
ｕeq(ｎ)＝（Ｓｂ₁−ｂ₀）^-1（ａ₁Ｓ−ａ₀−１）×｛ｘ₁(ｎ)−θ₁(ｎ)｝…（２６）
〈２〉非線形入力：
（i）ａ₀＞０のとき：
ｕnl(ｎ)＝｛η｜Ｓｂ₁＋ｂ₀｜^-1｜σ(ｎ)｜＋β(ｎ)｝×sgn［σ(ｎ)］…（２７）
（ii）ａ₀＜０のとき：
ｕnl(ｎ)＝｛η｜Ｓｂ₁−ｂ₀｜^-1｜σ(ｎ)｜＋β(ｎ)｝×sgn［σ(ｎ)］…（２８）
ただし、０＜η＜１、σ(ｎ＋１)σ(ｎ)＞０、
〈２〉総制御入力（空燃比操作量）：
（i）ａ₀＞０のとき：

（ii）ａ₀＜０のとき：

図１３のフローチャートはシステム行列の要素ａ₀の値に応じた総制御入力を演算するためのものである。これは図７の空燃比フィードバック補正値演算フローチャートにおけるステップ１４からステップ１８の具体例をなすものである。
【０１０２】
ステップ４１ではシステム行列matＡの要素ａ₀の符号をみる。ａ₀＞０となったときにはステップ４２〜４５で切換関数、線形入力、非線形入力、総制御入力の各演算を行う。ステップ１でａ₀＜０となったときにはステップ４６〜４９で同様の演算を行う。
【０１０３】
ここで、本実施形態（第１実施形態）の作用を説明する。
【０１０４】
本実施形態によれば非線形性の強いエンジンのようなプラントであっても運転領域毎に空燃比ステップ応答を測定することにより運転領域毎の伝達関数Ｇeng(ｑ)を同定することが可能である。従ってこの離散系２次の伝達関数Ｇeng(ｑ)を（２ａ）、（２ｂ）式に示す離散系２次の状態方程式に変換し、この離散系２次の状態方程式に基づいて空燃比センサ１６により検出される排気の空燃比ＡＦＳＡＦが目標空燃比ＴＧＡＢＦと一致するようにスライディングモード制御を用いた空燃比のフィードバック制御を行うことで、連続系の物理モデルからスライディングモード制御に用いる状態方程式を算出する先行技術の場合よりも精度良くプラントモデルを構築することができると共にＰＩＤ制御のような古典制御を用いる場合より応答性とロバスト性が向上する。かつステップ応答測定より一意的に切換関数が決まるのでスライディングモード制御の切換関数及び線形入力の設計が容易となる。
【０１０５】
また、変化する目標値（目標空燃比）への追従性を良くし及び定常偏差を吸収するため目標値追従問題であるサーボ系では状態方程式に積分項を追加した拡大系を用いるのが一般的であるのでスライディングモード制御においても拡大系を用いたサーボ系としたとき、切換関数及び非線形入力に積分項を持つことになり応答性が低下するのであるが、本実施形態によれば積分器２６を追加して設け線形入力ｕeq(ｎ)に代えてこれを積分した値を用いるようにしたので、応答性を確保しつつ定常偏差の吸収が可能となる。
【０１０６】
また、目標空燃比への収束を早めようと過大な空燃比操作量を与えたのでは運転性能や排気性能が悪化する可能性があるため各領域毎に最適な切換関数ゲインmatＳ及び積分ゲインμを決定する必要がありこの双方を最適化するのは非常に工数がかかるのであるが、本実施形態によれば積分ゲインμを決定しさえすれば、この積分ゲインμより（２１）に示した一定の関係式を用いて切換関数ゲインmatＳを算出することができるので、切換関数の設計が容易となる。
【０１０７】
また、フィードバック制御を固定時間間隔で演算する場合に回転速度が上昇すると単位時間当たりの燃料噴射イベントが増加し見かけ上フィードバックゲインが落ちてしまい高回転速度域での目標値への追従性が悪化するのであるが、本実施形態によれば積分ゲインμを回転速度に応じて演算するので、フィードバック制御を固定周期で行った場合でも見かけ上のゲインの低下を防ぐことが可能となり、全運転領域で最適な積分ゲインを得ることができる。
【０１０８】
また、本実施形態によればシステム行列matＡの要素であるａ₀の符号（システム行列の状態）によって切換関数ゲインmatＳの要素であるＳ₂の符号を切換えることでシステム行列の状態が変化しても空燃比操作量を正しい方向に付加することが可能となっている。
【０１０９】
次に本願の第２実施形態を説明する。
【０１１０】
ここで、第２実施形態と第１実施形態とでは創作した時期にずれがあり、本願の第２実施形態は、本願の基礎出願（特願２００１−２１１５７１号）とほぼ同時期に提案したもう一つの先願装置（特願２００１−２２４１６２号）にもヒントを得てなされているので、次にはそのもう一つの先願装置（このもう一つの先願装置を、以下「先願装置２」という。）の実施形態を説明し、その後に本願の第２実施形態を説明する。
【０１１１】
本願の基礎出願ではプラントに無駄時間がないシンプルなシステムで考えた。しかしながら実際のプラントでは無駄時間を無視できないので、次には無駄時間を有するプラントについて考える。
【０１１２】
図１４は無駄時間を有するプラントを対象とするシステムである場合のエンジンコントローラ１１の制御ブロック図で、図２に対して状態予測部３１を追加している。この状態予測部３１はシリンダ吸入空燃比演算部３２、モデルスケジュール部３３、プラントモデル状態変数演算部３４、予測状態変数演算部３５からなる。状態予測部３１により実行される制御内容は後述するフローチャート（図１７参照）によりまとめて説明する。従ってここでは図１４に示す個々のブロック３２〜３５の働きの説明は省略する。
【０１１３】
次に無駄時間を有するプラントを対象とするシステムについての考え方を詳述する。
【０１１４】
制御対象である無駄時間を有するプラント（以下「制御対象プラント」という）は次のように記述できる。
【０１１５】
matｘ(ｎ＋１)＝matＡmatｘ(ｎ)＋matＢｕ(ｎ−ｄ) …（３１ａ）
matｙ(ｎ)＝matＣmatｘ(ｎ) …（３１ｂ）
ただし、matｘ(ｎ) ：制御対象プラントの状態変数、
matｙ(ｎ) ：制御対象プラントの出力、
ｕ(ｎ−ｄ)：制御対象プラントの入力、
matＡ ：制御対象プラントのシステム行列、
matＢ ：制御対象プラントの入力行列、
matＣ ：制御対象プラントの出力行列、
ｄ ：制御対象プラントの無駄時間、
これを上記の（２ａ）、（２ｂ）式と比較すれば（３１ａ）のｕ(ｎ−ｄ)が違うだけである。
【０１１６】
また上記の（２ａ）、（２ｂ）式は無駄時間の無いプラント（以下「プラントモデル」という。）を対象とするシステムの状態方程式であり、これを制御対象プラントの場合と区別して扱うためｘの後に「Ｍ」の添え字を付け改めて次のように記載する。
【０１１７】
matｘ_M(ｎ＋１)＝matＡmatｘ_M(ｎ)＋matＢｕ(ｎ) …（３２ａ）
matｙ_M(ｎ)＝matＣmatｘ_M(ｎ) …（３２ｂ）
ただし、matｘ_M(ｎ) ：プラントモデルの状態変数、
matｙ_M(ｎ) ：プラントモデルの出力、
ｕ(ｎ) ：プラントモデルの入力、
matＡ ：プラントモデルのシステム行列、
matＢ ：プラントモデルの入力行列、
matＣ ：プラントモデルの出力行列、
ここで、プラントモデルの伝達関数は上記の（１）式であるが、制御対象プラントの伝達関数は次の通りである。
【０１１８】
Ｇeng(ｚ)＝（ｂ₁ｚ^-1＋ｂ₀ｚ^-2）／｛１＋（ｚ^-1−ｄ）ａ₁＋（ｚ^-1−ｄ）
²ａ₀｝ …（補８）
ただし、Ｇeng(ｑ)：制御対象プラントの伝達関数、
ａ₁、ａ₀ ：微分係数、
ｂ₁、ｂ₀ ：微分係数、
ｄ ：無駄時間、
さて制御対象プラントの有する無駄時間が空燃比に与える影響を排除するために図１５に示すようなスミス法を用いて補償する手法が知られている。このスミス法によればモデル無駄時間ｄ_M経過後の制御対象プラントの状態変数（この状態変数を以下「予測状態変数」という。）は一般に
【０１１９】
【数８】

【０１２０】
のように表される。ここで（３３）式左辺のmatｘに付けた山型（ハット）の印は予測値を表すので、この表示を以下「matｘ_hat」とも表記する。
【０１２１】
（３３）式を（３２ａ）、（３２ｂ）式を用いて変形すると次式が得られる。
【０１２２】
matｘ_hat(ｎ)
＝matｘ(ｎ＋ｄ_M)
＝matＡ^dmatｘ(ｎ)＋matｘ_M(ｎ)−matＡ^dmatｘ_M(ｎ−ｄ_M)…（３４）
なお（３４）式右辺のmatｘ_M(ｎ−ｄ_M)はモデル無駄時間ｄ_M経過後のプラントモデルの状態変数である。
【０１２３】
しかしながらスミス法は目標値入力に対してはよい応答を与えることができるが、外乱に対する過渡応答は必ずしも良くないことも知られている。そのために図１６に示すような外乱補償器を追加する手法が提案されている。外乱補償器の設計手法には様々なものが提案されているが（『むだ時間システムの制御』コロナ社発行参照）、本願発明では図１６のブロック図から実際の動特性に基づいた検討を行い、その検討結果に基づいて外乱補償器の伝達関数Ｍ(ｚ)の設計を行った。この検討結果を結論から先に述べると、状態変数matｘ(ｎ)が可観測であることを考慮して予測状態変数matｘ_hat(ｎ)を次のように表す。
【０１２４】
matｘ_hat(ｎ)
＝matｘ(ｎ＋ｄ_M)
＝matＭ(ｎ)｛matｘ(ｎ)−matｘ_M(ｎ−ｄ_M)｝＋matｘ_M(ｎ)…（３５ａ）
matＭ(ｎ)＝１ …（３５ｂ）
または
matｘ_hat(ｎ)
＝matｘ(ｎ＋ｄ_M)
＝matｘ(ｎ)−matｘ_M(ｎ−ｄ_M)＋matｘ_M(ｎ) …（３６）
ただし、matｘ_hat(ｎ) ：モデル無駄時間ｄ_M後の制御対象プラントの状態変数、
matｘ(ｎ) ：制御対象プラントの状態変数、
matｘ_M(ｎ) ：プラントモデルの状態変数、
matｘ_M(ｎ−ｄ_M) ：モデル無駄時間ｄ_M経過後のプラントモデルの状態変数、
matＭ(ｎ) ：外乱補償器のシステム行列、
なお、（３６）式をスカラー表示で表すと予測状態変数matｘ_hat(ｎ)の状態量（この状態量を以下「予測状態量」という。）ｘ_1hat(ｎ)、ｘ_2hat(ｎ)は次のようになる。
【０１２５】
ｘ_1hat(ｎ)＝ｘ₁(ｎ＋ｄ_M)＝ｘ₁(ｎ)−ｘ_1M(ｎ−ｄ_M)＋ｘ_1M(ｎ)…（３７ａ）
ｘ_2hat(ｎ)＝ｘ₂(ｎ＋ｄ_M)＝ｘ₂(ｎ)−ｘ_2M(ｎ−ｄ_M)＋ｘ_2M(ｎ)…（３７ｂ）
上記（３６）式の意味するところを外乱がない場合と外乱がある場合について分けて述べる。
〈１〉外乱がない場合：
（３６）式右辺の第２項以降の部分に着目すると、外乱が無くプラントモデルが正確に制御対象プラントを記述しているとすればモデル無駄時間ｄ_M経過後のプラントモデルの状態変数matｘ_M(ｎ−ｄ_M)は実際に観測される状態変数matｘ(ｎ)と等しくなるため（matｘ_M(ｎ−ｄ_M)＝matｘ(ｎ)）（３６）式右辺はプラントモデルの状態変数matｘ_M(ｎ)のみが残る。後に残ったこのプラントモデルの状態変数matｘ_M(ｎ)は実際の無駄時間ｄ経過後の制御対象プラントで検出される状態変数matｘ(ｎ)を予測しており、従ってこれを予測状態変数matｘ_hat(ｎ)とすることができる。この場合にはフィードバックコントローラは予測状態変数matｘ_hat(ｎ)をパラメータにフィードバック制御する形となる。
〈２〉外乱がある場合：
一方、外乱ｇが付加される場合またはプラントモデルが制御対象プラントを正確に記述できない場合にはモデル無駄時間ｄ_M経過後のプラントモデルの状態変数matｘ_M(ｎ−ｄ_M)は実際に観測される制御対象プラントの状態変数matｘ(ｎ)と異なる結果となるため、このままでは（３６）式右辺第２項以降のいずれも消去することができない。
【０１２６】
ここで外乱ｇがステップ入力であるかまたは制御対象プラントの状態変数matｘ(ｎ)の変化がステップ応答であると仮定すると、ステップ変化からモデル無駄時間ｄ_M経過後にプラントモデルの状態変数matｘ_M(ｎ)とプラント無駄時間ｄ_M経過後のプラントモデルの状態変数matｘ_M(ｎ−ｄ_M)とが等しい値となる（matｘ_M(ｎ−ｄ_M)＝matｘ_M(ｎ) ）。そこで外乱補償器のシステム行列matＭ(ｎ)＝１とすると、（３５）式右辺は実際に観測される制御対象プラントの状態変数matｘ(ｎ)のみとなる。この場合、フィードバックコントローラは制御対象プラントの状態変数matｘ(ｎ)をパラメータにフィードバックする形となり、状態予測フィードバックではないものの外乱の影響を受けずにフィードバックすることが可能となる。
【０１２７】
このように図１６において外乱補償器の伝達関数Ｍ(ｚ)＝１とすることでスミス法に対して外乱補償器を用いた無駄時間を補償する式が得られている。
【０１２８】
図１７のフローチャートは制御対象プラントの予測状態変数matｘ_hat(ｎ)を演算するためのものである。
【０１２９】
ステップ５１では制御対象プラントの状態量ｘ₁(ｎ)、ｘ₂(ｎ)を算出する。ここでｘ₁(ｎ)＝ＡＦＳＡＦである。ｘ₂(ｎ)は後述するようにして求める。これら２つの状態量ｘ₁(ｎ)とｘ₂(ｎ)から制御対象プラントの状態変数matｘ(ｎ)が定まる。
【０１３０】
ステップ５２ではシリンダ吸入空気量Ｑｃを算出する。このシリンダ吸入空気量Ｑｃは例えば次式のようにシリンダ吸気量パルス幅ＴＰ［ｍｓ／ｃｙｃｌｅ］より逆算して算出している。
【０１３１】
Ｑｃ［ｋｇ／ｃｙｃｌｅ］
＝ＴＰ［ｍｓ／ｃｙｃｌｅ］／ＫＣＯＮＳＴ＃［ｍｓ／ｋｇ］…（補９）
ステップ５３ではこのシリンダ吸入空気量Ｑｃよりプラントモデルのシステム行列matＡ、入力行列matＢ、モデル無駄時間ｄ_Mを算出する。これらの算出については図１８のフローにより説明する。図１８のフローチャートは図１７のステップ５３のサブルーチンである。
【０１３２】
ステップ６１ではシリンダ吸入空気量Ｑｃを読み込む。ステップ６２ではこのシリンダ吸入空気量Ｑｃに応じてシステム行列matＡの要素ａ₁、ａ₀と入力行列matＢの要素ｂ₁、ｂ₀を選択する。これは、エンジンの空燃比の動特性が運転領域に対して強い非線形を有し一つだけの動特性モデルで全運転領域を表すことが不可能であることから（運転領域毎にシステム行列matＡ、入力行列matＢが変化しプラントモデルの離散系２次の状態方程式が違ったものになる）、シリンダ吸入空気量Ｑｃにより運転領域を例えば８つに分けて各運転領域毎に異なる値の要素ａ₁、ａ₀、ｂ₁、ｂ₀を入れておくことにより運転領域毎にプラントモデルの離散系２次の状態方程式を切換えるようにするためのものである。従ってシリンダ吸入空気量Ｑｃが今Ｑｃ１とＱｃ２の間にあればモデル＃２に格納されている要素ａ₁、ａ₀、ｂ₁、ｂ₀を取り出す。なお、各運転領域毎の要素ａ₁、ａ₀、ｂ₁、ｂ₀は各運転領域毎の伝達関数を求めることにより得られる。
【０１３３】
ステップ６３ではシリンダ吸入空気量Ｑｃより所定のテーブルを検索することによりモデル無駄時間ｄ_Mを演算する。この値はシリンダ吸入空気量Ｑｃが大きくなるほど小さくなる値である。これは燃料噴射から空燃比センサ１６が実際の空燃比を検出するまでの無駄時間がシリンダ吸入空気の流速に依存し低流速域で無駄時間が大きくなるため、全運転領域を単一の無駄時間モデルで記述することはできない。そこでモデルシリンダ吸入空気量Ｑｃに対応して変化する実際の無駄時間に対応させて無駄時間ｄ_Mを算出するようにしたものである。
【０１３４】
図１７に戻りステップ５４ではシリンダ吸入空燃比ＣＹＬＡＦを算出する。すなわちシリンダ内燃料噴射パルス幅（壁流補正なし）ＴＩＰＳＣＹＬを
ＴＩＰＳＣＹＬ＝ＴＰ×ＴＦＢＹＡ×（ＡＬＰＨＡ＋ＫＢＬＲＣ−１）…（補１０）
の式により算出し、これから上記の（補９）式と同様にしてシリンダ吸入燃料量Ｆｃを次のように算出する。
【０１３５】
Ｆｃ［ｋｇ／ｃｙｃｌｅ］
＝ＴＩＰＳＣＹＬ［ｍｓ／ｃｙｃｌｅ］／ＫＣＯＮＳＴ＃［ｍｓ／ｋｇ］…（補１１）
ここで（補１０）式のようにシリンダ内燃料噴射パルス幅ＴＩＰＳＣＹＬを空燃比フィードバック補正値ＡＬＰＨＡと空燃比学習値ＫＢＬＲＣにより修正した値に基づいて算出するようにしたのは次の理由による。制御対象プラントの空燃比センサ１６や燃料噴射弁５には個体差があるので、空燃比フィードバック補正値ＡＬＰＨＡと空燃比学習値ＫＢＬＲＣにより修正する前の値に基づいてシリンダ内燃料噴射パルス幅を算出したのでは、空燃比センサ１６や燃料噴射弁５の個体差によってプラントモデルの状態方程式に誤差が生じてしまう。そこで空燃比センサ１６や燃料噴射弁５に個体差にがあってもプラントモデルの状態方程式に誤差が生じないように空燃比フィードバック補正値ＡＬＰＨＡと空燃比学習値ＫＢＬＲＣにより修正した値に基づいてシリンダ内燃料噴射パルス幅ＴＩＰＳＣＹＬを算出するようにしたものである。
【０１３６】
（補９）、（補１１）式より空気過剰率λｃｙｌはλｃｙｌ＝Ｑｃ／Ｆｃであるから、シリンダ吸入空燃比ＣＹＬＡＦは
ＣＹＬＡＦ＝λｃｙｌ×１４．７ …（補１２）
の式により求めることができる。（補９）、（補１１）式を（補１２）式に代入すると、上記の（補５）式が得られる。
【０１３７】
ステップ５５ではこのシリンダ吸入空燃比ＣＹＬＡＦを入力ｕ(ｎ)とする。
【０１３８】
ステップ５６では入力ｕ(ｎ)、システム行列matＡ、入力行列matＢを用いて上記（３２ａ）式の漸化式を用いてプラントモデルの状態量ｘ_1M(ｎ)とｘ_2M(ｎ)を算出する。これら２つの状態量ｘ_1M(ｎ)とｘ_2M(ｎ)からプラントモデルの状態変数matｘ_M(ｎ)が定まる。
【０１３９】
ステップ５７ではこのようにして求めた制御対象プラントの状態変数matｘ(ｎ)、プラントモデルの状態変数matｘ_M(ｎ)、モデル無駄時間ｄ_Mから予測状態量ｘ_1hat(ｎ)とｘ_2hat(ｎ)を上記の（３７ａ）、（３７ｂ）式を用いて算出する。これら２つの予測状態量ｘ_1hat(ｎ)とｘ_2hat(ｎ)から予測状態変数matｘ_hat(ｎ)が定まる。
【０１４０】
図１７は図１４の状態予測部３１をフローチャートで構成したものでもある。両者を対応づけると、図１７のステップ５３が図１４のモデルスケジュール部３３に、ステップ５４が図１４のシリンダ吸入空燃比演算部３２に、ステップ５６が図１４のプラントモデル状態変数演算部３４に、ステップ５７が図１４の予測状態変数演算部３５に相当する。
【０１４１】
ここで図１４の状態予測部３１で行われるところを物理的な実体に即した状態量に基づいて再度述べると次のようになる。プラントモデル状態変数演算部３４では現時刻ｎでのプラント入力ｕ(ｎ)であるシリンダ吸入空燃比ＣＹＬＡＦに基づいてそのときの運転領域におけるプラントモデルの現時刻ｎでの空燃比（モデル空燃比）を演算し、予測状態変数演算部３５では図１６に示した外乱補償型スミス法を用いて現時刻ｎでのモデル空燃比と実際に検出される現時刻ｎでの空燃比とから制御対象プラントの現時刻ｎよりモデル無駄時間ｄ_Mのちの触媒７上流の空燃比を予測（推定）する。
【０１４２】
このようにして演算される予測空燃比がスライディングモードコントローラ２２に渡されると、切換関数演算部２３では予測空燃比及び目標空燃比ＴＧＡＢＦに基づいて現時刻ｎにおける切換関数を上記（４）式に代え次の式により演算する。
【０１４３】
σ(ｎ)＝Ｓ₁（ｘ_1hat(ｎ)−θ₁(ｎ)）＋Ｓ₂（ｘ_2hat(ｎ)−θ₂(ｎ)）＝０…（３８）
ただし、σ(ｎ)：状態量、
また線形入力演算部２５では線形入力を上記（７）式に代えて次式により演算する。
【０１４４】
ｕeq(ｎ)＝−（matＳmatＢ）^-1matＳ
×｛matＡmatｘ_hat(ｎ)−matθ(ｎ＋１)｝−σ(ｎ）…（３９）
（３６）式、（３８）式よりスカラーで記述すると、（３９）式は次のようになる。
【０１４５】
ｕeq(ｎ）＝（Ｓ₁ｂ₁＋Ｓ₂ｂ₀）^-1（ａ₁Ｓ₁＋ａ₀Ｓ₂＋Ｓ₁ ²／Ｓ₂）
×｛ｘ_1hat(ｎ）−θ₁(ｎ)｝ …（４０）
ところで、モデル無駄時間ｄ_Mが小の運転領域から大の運転領域へと切換えられる場合に小領域で噴射した燃料による検出空燃比が大領域で過剰に無駄時間分遅らされ（切換時の予測状態量ｘ_2hat(ｎ)の誤差が大きくなる）、これに起因して空燃比が変動することが実験により判明している。これについて図１９を用いて具体的に説明すると、制御対象プラントの状態量ｘ₂(ｎ)に対する目標値θ₂(ｎ)、状態量ｘ₂(ｎ)、プラントモデルの状態量ｘ_2M(ｎ)、予測状態量ｘ_2hat(ｎ)の演算式は次の通りである。
【０１４６】
θ₂(ｎ)＝−ａ₀×θ₁(ｎ)＋ｂ₀×ｕ(ｎ−ｄ_M−１) …（４１）
ｘ₂(ｎ)＝−ａ₀×ｘ₁(ｎ−１)＋ｂ₀×ｕ(ｎ−ｄ_M−１) …（４２）
ｘ_2M(ｎ)＝−ａ₀×ｘ_1M(ｎ−１)＋ｂ₀×ｕ(ｎ−１) …（４３）
ｘ_2hat(ｎ)＝ｘ₂(ｎ＋ｄ_M)
＝ｘ₂(ｎ)−ｘ_2M(ｎ−ｄ_M)＋ｘ_2M(ｎ) …（４４）
ここで（４１）、（４２）、（４３）式は次のようにして得られるものである。まず（３１ａ）式よりｙ(ｎ)＝ｘ₁(ｎ)であるので、出力目標値θ₁はそのままｘ₁(ｎ)の目標値となる。（３１ａ）式を展開すると次のようになる。
【０１４７】
ｘ₁(ｎ＋１)＝−ａ₁ｘ₁(ｎ)＋ｘ₂(ｎ)＋ｂ₁ｕ(ｎ−ｄ) …（補１３）
ｘ₂(ｎ＋１)＝−ａ₀ｘ₁(ｎ)＋ｂ₀ｕ(ｎ−ｄ) …（補１４）
（補１４）式においてｎ＋１に代えてｎをおくと
ｘ₂(ｎ)＝−ａ₀ｘ₁(ｎ−１)＋ｂ₀ｕ(ｎ−１−ｄ) …（補１５）
となり（４２）式が得られる。ｘ₂(ｎ)とｘ_2M(ｎ)の間にはｘ₂(ｎ)に無駄時間がありｘ_2M(ｎ)に無駄時間がないとの違いがあるだけなので、（補１５）式のｕ(ｎ−１−ｄ) よりｄを除くとｘ_2M(ｎ)の式となる。すなわち
ｘ_2M(ｎ)＝−ａ₀ｘ₁(ｎ−１)＋ｂ₀ｕ(ｎ−１) …（補１６）
となり（４３）式が得られる。
【０１４８】
状態量ｘ₂(ｎ)の目標値は（補１４）式より
θ₂(ｎ＋１)＝−ａ₀θ₁(ｎ)＋ｂ₀ｕ(ｎ−ｄ) …（補１７）
である。これより
θ₂(ｎ)＝−ａ₀θ₁(ｎ−１)＋ｂ₀ｕ(ｎ−ｄ−１) …（補１８）
となる。θ₁(ｎ)＝θ₁(ｎ−１)であるから、
θ₂(ｎ)＝−ａ₀θ₁(ｎ)＋ｂ₀ｕ(ｎ−ｄ−１) …（補１９）
となり（４１）式が得られる。
【０１４９】
さて図１９に示したようにシリンダ吸入空気量Ｑｃの減少で（モデル入力ｕ(ｎ)は変わらないと仮定する）モデル無駄時間ｄ_Mが２から３へと変化しこれに対して目標値θ₂(ｎ)が減少すると仮定したとき、この目標値θ₂(ｎ)を追っかけるように状態量ｘ₂(ｎ)、ｘ_2M(ｎ)がいずれも減少する。予測状態量ｘ_2hat(ｎ)はモデル無駄時間ｄ_Mが３になったとき
ｘ_2hat(ｎ)＝ｘ₂(ｎ)−ｘ_2M(ｎ−３)＋ｘ_2M(ｎ) …（補２０）
となり、この値は時刻ｎにおいてｎ＋３での状態量ｘ₂(ｎ)を予測するものである。値を簡単にするため状態量ｘ₂(ｎ)、ｘ_2M(ｎ)がいずれも３から１へと減少したとすると、時刻ｎでの状態量ｘ₂(ｎ)、ｘ_2M(ｎ)はともに１であり、ｘ_2M(ｎ−３)はｎ−３でのｘ_2M(ｎ)の値であるから３である。故に時刻ｎでの（補２０）式は
ｘ_2hat(ｎ)＝１−３＋１＝−１
となりこれはｘ₂(ｎ)の値より小さい。そして状態量ｘ_2hat(ｎ)が−１となる状態はｎ＋３の直前まで（切換タイミングより切換後の領域でのモデル無駄時間が経過する直前まで）続き、ｎ＋３のタイミングでｘ₂(ｎ)と同じ１になる（図１９の実線参照）。しかしながら、予測値としてのｘ_2hat(ｎ)の値は制御上は同図においてｎからｎ＋３の直前までの間も一点鎖線で示したように１であることが望ましい。従ってｎからｎ＋３の直前までの間も上記の（４４）式により演算するｘ_2hat(ｎ)の値に基づくのではモデル誤差が生じる。
【０１５０】
そこでモデル無駄時間ｄ_Mが大きくなる側への切換時には切換直後より切換後の領域に対するモデル無駄時間が経過するまでｘ_2hat(ｎ＋ｄ_M)を上記の（４４）式により演算するのではなくｘ₂(ｎ)の値を用いる。つまりｎからｎ＋３の直前までの間は
ｘ_2hat(ｎ)＝ｘ₂(ｎ) …（補２１）
として外乱補償形スミス法による予測を停止する。
【０１５１】
図２０のフローチャートは図１７のステップ５７の内容をなすものである。
【０１５２】
ステップ７１でシリンダ吸入空気量Ｑｃを読み込み、ステップ７２でそのＱｃに基づいてモデル無駄時間ｄ_Mが大きくなる側に運転領域が切換わったかどうかをみる。モデル無駄時間ｄ_Mが大きくなる側へと切換わったタイミングにあれば図１６に示した外乱補償形スミス法による予測を停止するためステップ７３に進み予測状態変数matｘ_hat(ｎ)＝matｘ(ｎ)とする。そうでないときにはステップ７２よりステップ７４に進みモデル無駄時間ｄ_Mが大きくなる側への切換タイミングより切換後の運転領域に対応する無駄時間ｄ_Mが経過したかどうかをみる。モデル無駄時間ｄ_Mが大きくなる側への切換タイミングより切換後の運転領域に対応する無駄時間ｄ_Mが経過していなければステップ７３の処理を実行し、モデル無駄時間ｄ_Mが大きくなる側へのモデル切換タイミングより切換後の運転領域に対応する無駄時間ｄ_Mが経過したときステップ７５に進み図１６に示した外乱補償形スミス法による予測を行う。
【０１５３】
次に、無駄時間を有するプラントを対象とするシステムである場合の先願装置２の実施形態の作用を説明する。
【０１５４】
先願装置２の実施形態では燃料噴射（燃料供給）により定まる空燃比が空燃比センサ１６（空燃比検出手段）により検出されるまでの無駄時間を、スミス法に外乱補償器を加えることにより補償する無駄時間補償手段を設け、この無駄時間補償手段により無駄時間の補償された空燃比が目標空燃比と一致するようにスライディングモード制御を用いた空燃比のフィードバック制御を燃料供給量の増減で行う。詳細には状態予測部３１が、制御対象プラント（実際のエンジン）の排気系及び空燃比センサ１６の動特性を無駄時間を有する離散系２次の伝達関数で近似しシリンダ吸入空燃比のステップ応答よりＡＲＸモデルで同定することによりこの無駄時間を有する離散系２次の伝達関数を設定し、この無駄時間を有する離散系２次の伝達関数を制御対象プラントについての無駄時間を有する離散系２次の状態程式に変換し、プラントモデル（無駄時間のない仮想のエンジン）の排気系及び空燃比センサ１６の動特性を無駄時間のない離散系２次の伝達関数で近似しシリンダ吸入空燃比のステップ応答よりＡＲＸモデルで同定することによりこの無駄時間のない離散系２次の伝達関数を設定し、この無駄時間のない離散系２次の伝達関数をプラントモデルについての無駄時間のない離散系２次の状態方程式に変換し、スミス法に外乱補償器を加えて無駄時間補償を行った状態方程式を制御対象プラントについての無駄時間を有する離散系２次の状態方程式の状態変数matｘ(ｎ)とプラントモデルについての無駄時間のない離散系２次の状態方程式の状態変数matｘ_M(ｎ)とモデル無駄時間ｄ_Mとに基づいて演算する。そしてスライディングコントローラ２２が、このスミス法に外乱補償器を加えて無駄時間補償を行った状態方程式の状態変数matｘhat(ｎ)と目標値の状態量matθ(ｎ)の差を誤差の状態量matｅ(ｎ)として演算し、この誤差の状態量matｅ(ｎ)に切換関数ゲインmatＳを乗算して状態量σ(ｎ)を演算し、この切換関数ゲインの乗算された状態量σ(ｎ)に基づいて非線形入力ｕnl(ｎ)を演算し、同じくこの切換関数ゲインの乗算された状態量σ(ｎ)に基づいて線形入力ｕeq(ｎ)を演算し、これら非線形入力と線形入力の和を総制御入力ｕsl(ｎ)として算出し、この総制御入力に基づいて燃料供給量のフィードバック制御を行っている。
【０１５５】
このように構成することで、無駄時間を有するプラントを制御対象としてフィードバックゲインを大きくとっても空燃比が振動的になることがなく、かつ無駄時間の影響が特に大きいアイドル状態付近の低空気流量域においても先行技術と相違してフィードバックゲインを大きくできる。
【０１５６】
また、制御対象プラントには空燃比センサ１６や燃料噴射弁５（燃料供給手段）の個体差により誤差が生じるが、こうした個体差による誤差はシリンダ吸気量パルス幅ＴＰ（燃料供給量）を空燃比フィードバック補正値ＡＬＰＨＡや空燃比学習値ＫＢＬＲＣで修正する手段により解消される。従って空燃比フィードバック補正値ＡＬＰＨＡまたは空燃比学習値ＫＢＬＲＣで修正された燃料噴射パルス幅であるシリンダ内燃料噴射パルス幅ＴＩＰＳＣＹＬ（燃料供給量）に基づいてシリンダ吸入空燃比ＣＹＬＡＦを演算するようにした先願装置２の実施形態によれば制御対象エンジンの個体差による誤差をなくすことができる。
【０１５７】
また、制御対象プラント（実際のエンジン）の有する空燃比の動特性はエンジン負荷及びエンジン回転速度によって変化するため１つの状態方程式だと全ての運転領域に最適とならないのであるが、プラントモデルについての無駄時間のない離散系２次の状態方程式を負荷と回転速度によって定まるシリンダ吸入空気量Ｑｃにより多段階に切換えるようにした先願装置２の実施形態によれば、負荷及び回転速度から定まるどのような運転条件でもプラントモデルについての無駄時間のない離散系２次の状態方程式を最適に与えることができる。またこれにより空燃比フィードバック制御を固定時間間隔で実行した場合でも、単位時間当たりの燃焼回数の影響を排除できる。
【０１５８】
また、先願装置２の実施形態ではモデル無駄時間ｄ_Mをシリンダ吸入空気量Ｑｃに応じて演算するので制御対象プラントの実際の無駄時間ｄに近い値を推定できる。
【０１５９】
また、モデル無駄時間ｄ_Mが大きくなる側へと切換わる場合にもスミス法に外乱補償器を加えての無駄時間補償を継続するとフィードバックの過補正による空燃比変動が生じるのであるが、先願装置２の実施形態によればモデル無駄時間ｄ_Mが大きくなる側へと切換わったタイミングより切換後の運転領域に対応する無駄時間ｄ_Mが経過する直前まで過補正の原因となるスミス法に外乱補償器を加えての無駄時間補償を行わないことにしたので、モデル無駄時間ｄ_Mが大きくなる側へと切換わる場合に発生していたフィードバックの過補正による空燃比変動を抑えることができる。
【０１６０】
これで、先願装置２の実施形態の説明を終える。
【０１６１】
さて、本願の基礎出願後の実験により本願の上記第１実施形態には次の問題があり、なお改良の余地があることがわかった。
【０１６２】
問題点１：上記（補１４）式によれば、状態量ｘ₁(ｎ)とｘ₂(ｎ)の関係が複雑であり、入力ｕ(ｎ)が決まらないと、状態量ｘ₂(ｎ)を計算できない。このため切換関数の構成を上記（４）式のようにせざるを得ず、状態変数ｘ₂(ｎ)の物理的な意味が不明確になっていた。
【０１６３】
問題点２：上記（４）式の状態量σ(ｎ)は次式のようにも記載できる。
【０１６４】
σ(ｎ)＝Ｓ₁（ｘ₁(ｎ)−θ₁(ｎ)）−（ｘ₁(ｎ−１)−θ₁(ｎ−１)）…（５１）
この（５１）式は次のようにして導いたものである。すなわち、上記の（補１５）式、（補１８）式を上記（４）式に代入する。
【０１６５】

ここで、（５２）式のＳ₁／ａ₀Ｓ₂を改めてＳ₁と定義すれば、（５１）式が得られる。
【０１６６】
（５１）式により状態量σ(ｎ)を定義した場合、目標値θ₁(ｎ)の変化時には（５１）式右辺第１項、第２項とも実値と目標値に偏差を生じさせるためにσ(ｎ)の値が変化するが、目標値θ₁(ｎ)への到達以前に状態量σ(ｎ)の符号を変化させるために、（５１）式右辺は差分式になっており、結果的にσ(ｎ)の値は偏差量に対して小さなものとなる。一方、（５１）式の切換関数のゲインＳ₁は応答性能によって決定されるため必要以上に大きな値とすることはできない。このため小さな値にしかならないσ(ｎ)の値に応じて空燃比操作量ｕsl(ｎ)が決定されるとなると、目標値θ₁(ｎ)の変化当初の空燃比フィードバック補正値ＡＬＰＨＡの変化が図２１のように小さなものとなる場合がある（図２１下段の実線参照）。言い換えると、本願の第１実施形態は出力matｙ(ｎ)と状態量ｘ₁(ｎ)が一致するように構成したため（（補２）式、（補３）式参照）、状態量ｘ₂(ｎ)に対しては物理的な意味を持たすことができず、従って状態量ｘ₂(ｎ)の動作を物理的に説明するのが難しい構成であった。このため、切換関数の設計が最適な形となっておらず、目標値θ₁(ｎ)への追従性が十分でなかったのである（図２１上段の実線参照）。
【０１６７】
そこで、第２実施形態では、
ｘ₂(ｎ)＝ｘ₁(ｎ−１)…（５３）
の式が成立するように無駄時間を有する制御プラントと、無駄時間のないプラントモデルとについて、システム行列、入力行列、出力行列の要素（ａ₁、ａ₀、ｂ₁、ｂ₂、１、０）を組み換えてみたところ、無駄時間を有する制御プラントについての離散系２次の状態空間表現（状態方程式）を、上記（３１ａ）式、（３１ｂ）式に代えて、次の（５４ａ）式、（５４ｂ）式で、また無駄時間のないプラントモデルとについての離散系２次の状態空間表現（状態方程式）を、上記（３２ａ）式、（３２ｂ）式に代えて、次の（５５ａ）式、（５５ｂ）の式で記述すれば、（５３）式が成立することを見出した。
【０１６８】
ア）無駄時間を有する制御プラントについての状態方程式：
matｘ(ｎ＋１)＝matＤmatｘ(ｎ)＋matＥｕ(ｎ−ｄ) …（５４ａ）
matｙ(ｎ)＝matＦmatｘ(ｎ) …（５４ｂ）
ただし、matｘ(ｎ) ：制御対象プラントの状態変数、
matｙ(ｎ) ：制御対象プラントの出力、
ｕ(ｎ−ｄ)：制御対象プラントの入力、
matＤ ：制御対象プラントのシステム行列、
matＥ ：制御対象プラントの入力行列、
matＦ ：制御対象プラントの出力行列、
ｄ ：制御対象プラントの無駄時間、
イ）プラントモデルについての状態方程式：
matｘ_M(ｎ＋１)＝matＤmatｘ_M(ｎ)＋matＥｕ(ｎ) …（５５ａ）
matｙ_M(ｎ)＝matＦmatｘ_M(ｎ) …（５５ｂ）
ただし、matｘ_M(ｎ) ：プラントモデルの状態変数、
matｙ_M(ｎ) ：プラントモデルの出力、
ｕ(ｎ) ：プラントモデルの入力、
matＤ ：プラントモデルのシステム行列、
matＥ ：プラントモデルの入力行列、
matＦ ：プラントモデルの出力行列、
ここで、（５４ａ）、（５４ｂ）、（５５ａ）、（５５ｂ）の各式においてシステム行列matＤ、入力行列matＥ、出力行列matＦは、次の式で定義される値である。
【０１６９】
【数９】

【０１７０】
また、切換関数を、上記（４）式に代えて次式で定義する。
【０１７１】
σ(ｎ)＝matＳmatｅ(ｎ)
＝［Ｓ １］matｅ(ｎ)
＝Ｓ（ｘ₁(ｎ)−θ₁(ｎ)）＋（ｘ₁(ｎ)−ｘ₂(ｎ)）＝０ …（５９）
ただし、σ(ｎ) ：状態量、
matＳ ：切換関数ゲイン、
Ｓ、１ ：切換関数ゲインmatＳの要素、
θ₁(ｎ) ：ｘ₁(ｎ)の目標値、
このときのスライディングモード制御における誤差の状態量の変化を、図５に対応させて図２２に示す。
【０１７２】
上記の（５３）式をこの（５９）式右辺のｘ₂(ｎ)に代入してみると、次式が得られる。
【０１７３】
σ(ｎ)＝Ｓ（ｘ₁(ｎ)−θ₁(ｎ)）＋（ｘ₁(ｎ)−ｘ₁(ｎ−１)）…（６０）
ここで、（６０）式右辺の第１項は状態量ｘ₁(ｎ)とその目標値θ₁(ｎ)との差分を、第２項は状態量ｘ₁(ｎ)の微分値（制御周期当たりの変化量）を表す。従って、σ(ｎ)＝０とすることは、差分をゼロ、微分値をゼロすることであり、差分をゼロにすることは目標値に到達させることを、しかも微分値をゼロにすることはその目標値の位置に静止させることを意味する。すなわち、（５９）式によれば、切換関数を、状態量ｘ₁(ｎ)を目標値θ₁(ｎ)へと収束させるという所望の動作特性が得られるために必要な構成とすることができている。言い換えると切換関数は、本来、状態を最適な形に持っていくために設計するのであるが、上記（５３）式により状態量ｘ₂(ｎ)についても物理的な意味が明確になったため、状態量ｘ₁(ｎ)の目標値θ₁(ｎ)との差分及び状態量ｘ₁(ｎ)の微分値がゼロとなるような切換関数を設定できた。
【０１７４】
また、上記（６）式に代えて本願の第２実施形態では次のようになる。
【０１７５】
【数１０】

【０１７６】
従って、この（６１）式においてdet(matＳmatＥ)^-1≠０であればこのシステムの線形入力（等価制御入力）はスカラーで記述すると、次のようになる。
【０１７７】
ｕeq(ｎ）＝（ｂ₁＋ｂ₀）｛ａ₁ｘ₁(ｎ)＋ａ₀ｘ₂(ｎ)
−（ａ₀＋ａ₁）θ₁(ｎ)＋（Ｓ＋１）^-1（ｘ₁(ｎ)−θ₁(ｎ)））｝…（６２）
次に、非線形入力ｕnl(ｎ)は状態量を切換線に向けて移動することを目的としており、チャタリングを最小に抑えることが要求される。切換関数σ(ｎ)＝matＳmatｅ(ｎ)＝０と状態量との距離がチャタリングに密接に関係して影響を及ぼすため、図２３のような平滑化関数（縦軸の値）を導入する。すなわち、１入力の場合の非線形入力は、
ｕnl(ｎ)＝−η・sgn［σ(ｎ)］＝−η・σ(ｎ)／｜σ(ｎ)｜…（６３）
の式で定義されることから、チャタリングを除去するため平滑化関数を導入すると、非線形入力は（６３）式に代えて次のようになる。
【０１７８】
ｕnl(ｎ)＝−η・σ(ｎ)／（｜σ(ｎ)｜＋δ）…（６４）
ただし、η：非線形ゲイン、
δ（＞０）：平滑化定数、
ここで、本願の第１実施形態の上記（１５）式の非線形入力ｕnl(ｎ)は、周知の到達則を用いて与えたものであったが、この（１５）式の非線形入力ｕnl(ｎ)によれば、チャタリングは出ないものの応答性が悪いという特性がある。一方、上記（６３）式によれば応答性はよいがチャタリングが出る。これに対して上記（６４）式によれば応答性がよい上にチャタリングも出ないのである。
【０１７９】
また、（６４）式によれば単純な正負で補正する量を与えるため、本願の第１実施形態の上記（１５）式よりも式そのものがシンプルとなった。
【０１８０】
このように、本願の第２実施形態によれば、状態量ｘ₁(ｎ)、ｘ₂(ｎ)が（５３）式の関係を満たすように、上記（５４ａ）、（５４ｂ）、（５５ａ）、（５５ｂ）の各式で示した状態方程式のシステム行列matＤ、入力行列matＥ、出力行列matＦを見直し、それらの行列を上記（５６）式、（５７）式、（５８）式により定義される値としたので、状態量ｘ₂(ｎ)についても物理的な意味が明確となった。すなわち、本願の第２実施形態によれば、上記（５３）式で示したように状態量ｘ₂(ｎ)は、状態量ｘ₁(ｎ)の前回値を表すこととなる。
【０１８１】
また、切換関数についても、上記（５９）式に示したように切換関数を、状態量ｘ₁(ｎ)について所望の動作特性が得られるために必要な構成としたので、目標値θ₁(ｎ)の変化に対しても十分な空燃比補正値を得ることができる。
【０１８２】
例えば、目標値θ₁(ｎ)が理論空燃比（ほぼ１４．７）からリーン空燃比（例えば１６）へとステップ変化した場合の空燃比応答を図２１に示すと、本願の第２実施形態では状態量ｘ₁(ｎ)の応答が本願の第１実施形態より良くなるため（図２１上段の破線参照）、空燃比フィードバック補正値ＡＬＰＨＡの応答も本願の第１実施形態より良くなり（図２１下段の破線参照）、空燃比制御性能が向上している。
【０１８３】
実施形態では２入力、２出力の状態方程式の場合で説明したが、これに限られるものでない。
【０１８４】
実施形態では目標空燃比が変化する場合で説明したが、目標空燃比が一定の場合にも適用がある。例えば理論空燃比を目標とする空燃比フィードバック制御にも適用することができる。
【図面の簡単な説明】
【図１】ガソリンエンジンのシステム構成図。
【図２】スライディングモードコントローラの制御ブロック図。
【図３】無駄時間のない離散系２次のプラントモデル。
【図４】プラントモデルを同定するのに用いるシリンダ内空燃比のステップ応答を示す波形図。
【図５】スライディングモード制御における誤差の状態量の変化を位相平面上に示した図。
【図６】スライディングモード制御設計を説明するためのフローチャート。
【図７】空燃比フィードバック補正値の演算を説明するためのフローチャート。
【図８】切換関数ゲインの設定を説明するためのフローチャート。
【図９】固定積分ゲインを用いたときの低回転速度時のＡＬＰＨＡの変化波形図。
【図１０】固定積分ゲインを用いたときの高回転速度時のＡＬＰＨＡの変化波形図。
【図１１】可変積分ゲインを用いたときの高回転速度時のＡＬＰＨＡの変化波形図。
【図１２】積分ゲインの演算を説明するためのフローチャート。
【図１３】システム行列の要素ａ₀に応じた総制御入力の演算を説明するためのフローチャート。
【図１４】先願装置２の実施形態のスライディングモードコントローラの制御ブロック図。
【図１５】スミス法を伝達関数で表示したブロック図。
【図１６】外乱補償形スミス法を伝達関数で表示したブロック図。
【図１７】先願装置２の実施形態の予測状態変数の演算を説明するためのフローチャート。
【図１８】先願装置２の実施形態のシステム行列、入力行列、モデル無駄時間の演算を説明するためのフローチャート。
【図１９】先願装置２の実施形態のモデル無駄時間の切換時の解析を示すための波形図。
【図２０】先願装置２の実施形態のモデル無駄時間切換時の予測状態変数の演算を説明するためのフローチャート。
【図２１】本願の第２実施形態の作用を説明するための波形図。
【図２２】本願の第２実施形態のスライディングモード制御における誤差の状態量の変化を位相平面上に示した図。
【図２３】本願の第２実施形態の平滑化関数の特性図。
【符号の説明】
５ 燃料噴射弁
７ 三元触媒
１１ エンジンコントローラ
１６ 上流側広域空燃比センサ（空燃比検出手段）
２１ 目標空燃比演算部
２２ スライディングモードコントローラ

Claims (13)

1. 排気の空燃比を検出する空燃比検出手段と、
   エンジンの排気系及び前記空燃比検出手段の動特性を離散系２次の伝達関数で近似しシリンダ吸入空燃比のステップ応答より同定することによりこの離散系２次の伝達関数を設定する伝達関数設定手段と、
   この離散系２次の伝達関数を離散系２次の状態方程式に変換する手段と、
   この離散系２次の状態方程式に基づいて前記空燃比検出手段により検出される排気の空燃比が目標空燃比と一致するようにスライディングモード制御を用いた空燃比のフィードバック制御を行う空燃比フィードバック制御手段と
   を備え、
   前記離散系２次の状態方程式は、
   mat ｘ ( ｎ＋１ ) ＝ mat Ａ mat ｘ ( ｎ ) ＋ mat Ｂｕ ( ｎ )
   mat ｙ ( ｎ ) ＝ mat Ｃ mat ｘ ( ｎ )
   ただし、 mat ｘ ( ｎ＋１ ) ：時刻ｎ＋１の状態変数、
   mat ｘ ( ｎ ) ：時刻ｎの状態変数、
   mat ｙ ( ｎ ) ：時刻ｎの出力、
   ｕ ( ｎ ) ：時刻ｎの入力、
   mat Ａ ：システム行列、
   mat Ｂ ：入力行列、
   mat Ｃ ：出力行列、
   であることを特徴とするエンジンの空燃比制御装置。
2. 前記同定を最小２乗法を用いて行うことを特徴とする請求項１に記載のエンジンの空燃比制御装置。
3. 前記同定をＡＲＸモデルを用いて行うことを特徴とする請求項１に記載のエンジンの空燃比制御装置。
4. 前記空燃比のフィードバック制御は燃料供給量のフィードバック制御であることを特徴とする請求項１から３までのいずれか一つに記載のエンジンの空燃比制御装置。
5. 前記空燃比フィードバック制御手段は、前記離散系２次の状態方程式の状態変数と目標値の状態量との差を誤差の状態量として演算する手段と、この誤差の状態量に切換関数ゲインを乗算して状態量を演算する手段と、この切換関数ゲインの乗算された状態量に基づいて非線形入力を演算する手段と、同じくこの切換関数ゲインの乗算された状態量に基づいて線形入力を演算する手段と、これら非線形入力と線形入力の和を総制御入力として算出する手段と、この総制御入力に基づいて前記燃料供給量のフィードバック制御値を演算する手段とからなることを特徴とする請求項４に記載のエンジンの空燃比制御装置。
6. 目標空燃比が変化する場合に前記線形入力に代えてこれを積分した値を用いることを特徴とする請求項５に記載のエンジンの空燃比制御装置。
7. 排気通路に酸素ストレージ機能を有する触媒と、この触媒の前後にあって排気中の酸素濃度を検出する一対のセンサとを備える場合に、この一対のセンサの出力に基づいて前記触媒の酸素ストレージ量を最適化する空燃比を目標空燃比として演算することを特徴とする請求項６に記載のエンジンの空燃比制御装置。
8. 前記積分する際の積分ゲインと前記切換関数ゲインとの関係を前記線形入力の式から求め、前記目標空燃比への応答の時定数が最適となるように前記積分ゲインを決定し、この決定した積分ゲインより前記関係を用いて切換関数ゲインを算出することを特徴とする請求項６に記載のエンジンの空燃比制御装置。
9. 前記積分ゲインを回転速度に応じて定めることを特徴とする請求項８に記載のエンジンの空燃比制御装置。
10. 前記システム行列matＡの状態により前記切換関数ゲインmatＳの符号を変化させることを特徴とする請求項１に記載のエンジンの空燃比制御装置。
11. 前記離散系２次の伝達関数が
    Ｇeng(ｑ)＝（ｂ₁ｑ＋ｂ₀）／（ｑ²＋ａ₁ｑ＋ａ₀）
    （無駄時間の項を除く）
    ただし、Ｇeng(ｑ)：プラントの伝達関数、
    ｑ ：離散系シフトオペレータ、
    ａ₁、ａ₀ ：微分係数、
    ｂ₁、ｂ₀ ：微分係数、
    である場合に、前記システム行列matＡ、入力行列matＢ、出力行列matＣは次式により定義される値であることを特徴とする請求項１に記載のエンジンの空燃比制御装置。
    

    
12. 前記離散系２次の伝達関数が
    Ｇeng(ｑ)＝（ｂ₁ｑ＋ｂ₀）／（ｑ²＋ａ₁ｑ＋ａ₀）
    （無駄時間の項を除く）
    ただし、Ｇeng(ｑ)：プラントの伝達関数、
    ｑ ：離散系シフトオペレータ、
    ａ₁、ａ₀ ：微分係数、
    ｂ₁、ｂ₀ ：微分係数、
    である場合に、前記システム行列matＡ、入力行列matＢ、出力行列matＣは次式により定義される値であることを特徴とする請求項１に記載のエンジンの空燃比制御装置。
    

    
13. 前記切換関数ゲインの乗算された状態量σ(ｎ)は、
    σ(ｎ)＝Ｓ（ｘ₁(ｎ)−θ₁(ｎ) ）＋（ｘ₁(ｎ)−ｘ₂(ｎ) ）
    ただし、ｘ₁(ｎ)、ｘ₂(ｎ) ：状態変数matｘ(ｎ) の要素、
    θ₁(ｎ) ：ｘ₁(ｎ)の目標値、
    であることを特徴とする請求項１２に記載のエンジンの空燃比制御装置。

Images