JP3602811B2 - Plant control equipment - Google Patents
Plant control equipment Download PDFInfo
- Publication number
- JP3602811B2 JP3602811B2 JP2001184540A JP2001184540A JP3602811B2 JP 3602811 B2 JP3602811 B2 JP 3602811B2 JP 2001184540 A JP2001184540 A JP 2001184540A JP 2001184540 A JP2001184540 A JP 2001184540A JP 3602811 B2 JP3602811 B2 JP 3602811B2
- Authority
- JP
- Japan
- Prior art keywords
- model parameter
- vector
- plant
- value
- throttle valve
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Lifetime
Links
Images
Landscapes
- Control Of Throttle Valves Provided In The Intake System Or In The Exhaust System (AREA)
- Electrical Control Of Air Or Fuel Supplied To Internal-Combustion Engine (AREA)
- Combined Controls Of Internal Combustion Engines (AREA)
- Feedback Control In General (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、プラントの制御装置に関し、特に制御対象であるプラントをモデル化して得られる制御対象モデルの特性を示すモデルパラメータを、リアルタイムで同定しつつ、そのモデルパラメータを用いて制御を行うものに関する。
【０００２】
【従来の技術】
制御対象であるプラントをモデル化し、その制御対象モデルのモデルパラメータをパラメータ調整機構により算出し、このモデルパラメータを用いる適応制御器により、内燃機関に供給する混合気の空燃比を目標空燃比にフィードバック制御する空燃比制御装置が従来より知られている（例えば特開平１１−７３２０６号公報）。
【０００３】
この制御装置では、モデルパラメータの初期値に、モデルパラメータの同定誤差に応じて算出される更新成分を加算することにより、モデルパラメータが算出される。モデルパラメータの更新成分は、同定誤差の過去値の影響を減少させるように修正され、モデルパラメータが外乱の影響でドリフトし、制御が不安定化することが防止される。
【０００４】
【発明が解決しようとする課題】
上記公報に記載された装置では、外乱の影響によるモデルパラメータのドリフトは防止される。しかし、パラメータ調整機構により算出されたモデルパラメータが、そのまま適応制御器で使用されるため、適応制御器の安定性を維持する上で改善の余地があった。
【０００５】
本発明はこの点に着目してなされたものであり、制御対象のプラントをモデル化した制御対象モデルのモデルパラメータを同定しつつ、そのモデルパラメータを使用する制御の安定性をより向上させることができるプラントの制御装置を提供することを目的とする。
【０００６】
【課題を解決するための手段】
上記目的を達成するため請求項１に記載の発明は、プラントをモデル化した制御対象モデルのモデルパラメータベクトルを、前記プラントの入力及び出力に基づいて同定する同定手段と、該同定手段により同定されたモデルパラメータベクトルを用いて前記プラントを制御する制御手段とを備えたプラントの制御装置において、前記同定手段は、前記モデルパラメータベクトルの同定誤差（ｉｄｅ）を算出する同定誤差算出手段と、該同定誤差に応じて更新ベクトル（ｄθ）を算出する更新ベクトル算出手段と、前記同定誤差の過去値の影響を減少させるように前記更新ベクトルを修正する更新ベクトル修正手段と、前記プラントの動特性変化を示すパラメータに応じて前記モデルパラメータの基準ベクトルを算出する基準ベクトル算出手段と、前記基準ベクトル（θｂａｓｅ）に修正された更新ベクトルを加算することにより、前記モデルパラメータベクトルを算出するモデルパラメータベクトル算出手段と、前記モデルパラメータベクトル算出手段により算出されるモデルパラメータベクトルの要素（ａ１，ａ２，ｂ１，ｃ１）の値を、所定リミット範囲内に制限する制限手段とを有することを特徴とする。
【０００７】
この構成によれば、モデルパラメータベクトルの同定誤差に応じて更新ベクトルが算出され、同定誤差の過去値の影響を減少させるように更新ベクトルが修正され、モデルパラメータの基準ベクトルに修正された更新ベクトルを加算することにより、モデルパラメータベクトルが算出され、さらにモデルパラメータベクトルの要素の値が所定リミット範囲内に制限される。したがって、モデルパラメータのドリフトを防止しつつ、制御の安定性をより一層向上させることができる。またモデルパラメータの基準ベクトルは、前記プラントの動特性変化を示すパラメータに応じて算出されるので、プラントの動特性変化に応じた適切な基準ベクトルが得られる。その結果、特にプラントに非線形要素が含まれる場合でもモデルパラメータを迅速に収束させることができる。
【０００８】
請求項２に記載の発明は、請求項１に記載のプラントの制御装置において、前記更新ベクトル算出手段は、固定ゲインアルゴリズムを用いて前記更新ベクトルの算出を行うことを特徴とする。
この構成によれば、固定ゲインアルゴリズムを用いて更新ベクトルが算出されるので、演算量を低減することができる。
【０００９】
請求項３に記載の発明は、請求項１または２に記載のプラントの制御装置において、前記更新ベクトル修正手段は、前記更新ベクトルの少なくとも１つの要素の過去値に０より大きく１より小さい所定値（ＤＥＬＴＡｉ、ＥＰＳｉ）を乗算することにより、前記更新ベクトルの修正を行うことを特徴とする。
【００１０】
この構成によれば、更新ベクトルの少なくとも１つの要素の過去値に０より大きく１より小さい所定値を乗算することにより、更新ベクトルの修正が行われるので、同定誤差の過去値の影響が減少し、モデルパラメータベクトルのドリフトを防止することができる。
【００１１】
請求項４に記載の発明は、請求項３に記載のプラントの制御装置において、前記更新ベクトル修正手段は、前記更新ベクトルの、前記プラントの入力に関わる要素（ｂ１の演算に関わる要素）または前記プラントの入出力に関わらない要素（ｃ１の演算に関わる要素）については、前記所定値（ＤＥＬＴＡｉ、ＥＰＳｉ）を乗算しないことを特徴とする。
【００１２】
この構成によれば、更新ベクトルの、プラントの入力に関わる要素またはプラントの入出力に関わらない要素については、０より大きく１より小さい所定値が乗算されないので、更新ベクトルの修正による定常偏差の発生を防止することができる。
【００１５】
請求項５に記載の発明は、プラントをモデル化した制御対象モデルのモデルパラメータベクトルを、前記プラントの入力及び出力に基づいて同定する同定手段と、該同定手段により同定されたモデルパラメータベクトルを用いて、前記プラントの出力が制御目標値と一致するように前記プラントを制御する制御手段とを備えたプラントの制御装置において、前記同定手段は、前記モデルパラメータベクトルの同定誤差（ｉｄｅ）を算出する同定誤差算出手段と、前記プラントの出力または制御目標値に応じて不感帯を設定する不感帯設定手段と、該同定誤差が前記不感帯内（−ＥＩＤＮＲＬＭＴ≦ｉｄｅ≦ＥＩＤＮＲＬＭＴ）にあるときは、前記同定誤差を減少方向に修正する同定誤差修正手段と、該同定誤差修正手段により修正された同定誤差（ｉｄｅｎｌ）を用いて前記モデルパラメータベクトルを算出するモデルパラメータベクトル算出手段と、前記モデルパラメータベクトル算出手段により算出されるモデルパラメータベクトルの要素の値を、所定リミット範囲内に制限する制限手段とを有することを特徴とする。
【００１６】
この構成によれば、モデルパラメータの同定誤差が不感帯内にあるときは、同定誤差が減少方向に修正され、該修正された同定誤差を用いてモデルパラメータベクトルが算出され、さらにモデルパラメータベクトルの要素の値が所定リミット範囲内に制限される。したがって、モデルパラメータのドリフトを防止しつつ、制御の安定性をより一層向上させることができる。また前記不感帯は、プラントの出力または制御目標値に応じて設定されるので、モデルパラメータの値に反映させるべき同定誤差を、残留同定誤差として無視してしまうことを防止することができる。
【００１７】
請求項６に記載の発明は、請求項１から５の何れか１項に記載のプラントの制御装置において、前記制限手段は、前記モデルパラメータベクトルの複数の要素（ａ１，ａ２）が所定の関係（図１８）を満たすように、前記複数の要素の値を制限することを特徴とする。
【００１８】
この構成によれば、モデルパラメータベクトルの複数の要素が所定の関係を満たすように、複数の要素の値が制限されるので、そのモデルパラメータベクトルを用いる制御の安定性を向上させることができる。
請求項７に記載の発明は、請求項５に記載のプラントの制御装置において、前記同定誤差修正手段は、前記同定誤差（ｉｄｅ）が前記不感帯内（−ＥＩＤＮＲＬＭＴ≦ｉｄｅ≦ＥＩＤＮＲＬＭＴ）にあるときは、前記同定誤差（ｉｄｅ）を０とすることを特徴とする。
【００１９】
この構成によれば、同定誤差が前記不感帯内にあるときは、同定誤差が０とされるので、モデルパラメータの値に反映されるべきでない同定誤差の影響を無くし、モデルパラメータのドリフト防止効果を高めることができる。
請求項８に記載の発明は、請求項１から７の何れか１項に記載のプラントの制御装置において、前記プラントは、内燃機関のスロットル弁と、該スロットル弁を駆動する駆動手段とを有するスロットル弁駆動装置を含み、前記制御手段は、前記スロットル弁の開度を目標開度に一致させるように、前記スロットル弁駆動装置への制御入力を決定するパラメータを算出することを特徴とする。
【００２０】
この構成によれば、同定手段により同定されたモデルパラメータを用いて、スロットル弁開度を目標開度に一致させる制御が行われるので、スロットル弁開度の目標開度への制御性を向上させるとともに、スロットル弁開度制御の安定性を向上させることができる。
【００２１】
【発明の実施の形態】
以下本発明の実施の形態を図面を参照して説明する。
（第１の実施形態）
図１は本発明の第１の実施形態にかかるスロットル弁制御装置の構成を示す図である。内燃機関（以下「エンジン」という）１の吸気通路２には、スロットル弁３が設けられている。スロットル弁３には、該スロットル弁３を閉弁方向に付勢する第１付勢手段としてのリターンスプリング４と、該スロットル弁３を開弁方向に付勢する第２付勢手段としての弾性部材５とが取り付けられている。またスロットル弁３は、駆動手段としてのモータ６によりギヤ（図示せず）を介して駆動できるように構成されている。モータ６による駆動力がスロットル弁３に加えられない状態では、スロットル弁３の開度ＴＨは、リターンスプリング４の付勢力と、弾性部材５の付勢力とが釣り合うデフォルト開度ＴＨＤＥＦ（例えば５度）に保持される。
【００２２】
モータ６は、電子制御ユニット（以下「ＥＣＵ」という）７に接続されており、その作動がＥＣＵ７により制御される。スロットル弁３には、スロットル弁開度ＴＨを検出するスロットル弁開度センサ８が設けられており、その検出信号は、ＥＣＵ７に供給される。
【００２３】
またＥＣＵ７には、エンジン１が搭載された車両の運転者の要求出力を検出するアクセルペダルの踏み込み量ＡＣＣを検出するアクセルセンサ９が接続されており、その検出信号がＥＣＵ７に供給される。
ＥＣＵ７は、スロットル弁開度センサ８及びアクセルセンサ９の検出信号が供給される入力回路、入力信号をディジタル信号に変換するＡＤ変換回路、各種演算処理を実行する中央演算ユニット（ＣＰＵ）、ＣＰＵが実行するプログラムやプログラムで参照されるマップやテーブルなどを格納するメモリ、及びモータ６に駆動電流を供給する出力回路を備えている。ＥＣＵ７は、アクセルペダルの踏み込み量ＡＣＣに応じてスロットル弁３の目標開度ＴＨＲを決定し、検出したスロットル弁開度ＴＨが目標開度ＴＨＲと一致するようにモータ６の制御量ＤＵＴを決定し、制御量ＤＵＴに応じた電気信号をモータ６に供給する。
【００２４】
本実施形態では、スロットル弁３、リターンスプリング４、弾性部材５及びモータ６からなるスロットル弁駆動装置１０を制御対象とし、該制御対象に対する入力をモータ６に印加する電気信号のデューティ比ＤＵＴとし、制御対象の出力をスロットル弁開度センサ８により検出されるスロットル弁開度ＴＨとする。
【００２５】
スロットル弁駆動装置１０の応答周波数特性を実測すると、図２に実線で示すゲイン特性及び位相特性が得られる。そこで、下記式（１）で定義されるモデルを制御対象モデルとして設定した。このモデルの応答周波数特性は、図２に破線で示すようになり、スロットル弁駆動装置１０の特性に近似していることが確認されている。
ＤＴＨ（ｋ）＝ＴＨ（ｋ）−ＴＨＤＥＦ （２）
ここで、ＴＨは検出したスロットル弁開度、ＴＨＤＥＦは前記デフォルト開度である。
また式（１）のａ１，ａ２，ｂ１，ｃ１は、制御対象モデルの特性を決めるモデルパラメータであり、ｄはむだ時間である。
【００２６】
上記式（１）で定義されるモデルは、適応制御の適用を容易にするために採用した、離散時間系のＤＡＲＸモデル（delayed autoregressive model with exogeneous input:外部入力を持つ自己回帰モデル）である。
式（１）においては、出力の偏差量ＤＴＨに関わるモデルパラメータａ１，ａ２、入力のデューティ比ＤＵＴに関わるモデルパラメータｂ１の他に、入出力に関わらないモデルパラメータｃ１が設定されている。このモデルパラメータｃ１は、デフォルト開度ＴＨＤＥＦのずれやスロットル弁駆動装置に加わる外乱を示すパラメータである。すなわち、モデルパラメータ同定器により、モデルパラメータａ１，ａ２，ｂ１と同時にモデルパラメータｃ１を同定することにより、デフォルト開度ずれや外乱を同定できるようにしている。
【００２７】
図３は、ＥＣＵ７により実現されるスロットル弁制御装置の機能ブロック図であり、この制御装置は、適応スライディングモードコントローラ２１と、モデルパラメータ同定器２２と、むだ時間ｄが経過した後の予測スロットル弁開度偏差量（以下「予測偏差量」という）ＰＲＥＤＴＨ（ｋ）（＝ＤＴＨ（ｋ＋ｄ））を算出する状態予測器２３と、アクセルペダル踏み込み量ＡＣＣに応じてスロットル弁３の目標開度ＴＨＲを設定する目標開度設定部２４とからなる。
【００２８】
適応スライディングモードコントローラ２１は、検出したスロットル弁開度ＴＨが目標開度ＴＨＲと一致するように、適応スライディングモード制御によりデューティ比ＤＵＴを算出し、該算出したデューティ比ＤＵＴを出力する。
適応スライディングモードコントローラ２１を用いることにより、スロットル弁開度ＴＨの目標開度ＴＨＲへの追従応答特性を、所定のパラメータ（ＶＰＯＬＥ）を用いて適宜変更することが可能となり、その結果スロットル弁３を開弁位置から全閉位置に移動させる際の衝撃（スロットル全閉ストッパへの衝突）の回避、及びアクセル操作に対するエンジンレスポンスの可変化が可能となる。また、モデルパラメータの誤差に対する安定性を確保することが可能となる。
【００２９】
モデルパラメータ同定器２２は、修正モデルパラメータベクトルθＬ（θＬＴ＝［ａ１，ａ２，ｂ１，ｃ１］）を算出し、適応スライディングモードコントローラ２１に供給する。より具体的には、モデルパラメータ同定器２２は、スロットル弁開度ＴＨ及びデューティ比ＤＵＴに基づいて、モデルパラメータベクトルθを算出する。さらに、そのモデルパラメータベクトルθに対してリミット処理を行うことにより修正モデルパラメータベクトルθＬを算出し、該修正モデルパラメータベクトルθＬを適応スライディングモードコントローラ２１に供給する。このようにしてスロットル弁開度ＴＨを目標開度ＴＨＲに追従させるために最適なモデルパラメータａ１，ａ２，ｂ１が得られ、さらに外乱及びデフォルト開度ＴＨＤＥＦのずれを示すモデルパラメータｃ１が得られる。
【００３０】
リアルタイムでモデルパラメータを同定するモデルパラメータ同定器２２を用いることにより、エンジン運転条件の変化への適応、ハードウエアの特性ばらつきの補償、電源電圧変動の補償、及びハードウエア特性の経年変化への適応が可能となる。
【００３１】
状態予測器２３は、スロットル弁開度ＴＨ及びデューティ比ＤＵＴに基づいて、むだ時間ｄ後のスロットル弁開度ＴＨ（予測値）、より具体的には予測偏差量ＰＲＥＤＴＨを算出し、適応スライディングモードコントローラ２１に供給する。予測偏差量ＰＲＥＤＴＨを用いることにより、制御対象のむだ時間に対する制御系のロバスト性を確保し、特にむだ時間が大きいデフォルト開度ＴＨＤＥＦ近傍での制御性を向上させることができる。
【００３２】
次に適応スライディングモードコントローラ２１の動作原理を説明する。
先ず下記式（３）により、目標値ＤＴＨＲ（ｋ）を目標開度ＴＨＲ（ｋ）とデフォルト開度ＴＨＤＥＦとの偏差量として定義する。
ＤＴＨＲ（ｋ）＝ＴＨＲ（ｋ）−ＴＨＤＥＦ （３）
ここで、スロットル弁開度偏差量ＤＴＨと、目標値ＤＴＨＲとの偏差ｅ（ｋ）を下記式（４）で定義すると、適応スライディングモードコントローラの切換関数値σ（ｋ）は、下記式（５）にように設定される。
【００３３】
縦軸を偏差ｅ（ｋ）とし、横軸を前回偏差ｅ（ｋ−１）として定義される位相平面上では、σ（ｋ）＝０を満たす偏差ｅ（ｋ）と、前回偏差ｅ（ｋ−１）との組み合わせは、直線となるので、この直線は一般に切換直線と呼ばれる。スライディングモード制御は、この切換直線上の偏差ｅ（ｋ）の振る舞いに着目した制御であり、切換関数値σ（ｋ）が０となるように、すなわち偏差ｅ（ｋ）と前回偏差ｅ（ｋ−１）の組み合わせが位相平面上の切換直線上に載るように制御を行い、外乱やモデル化誤差（実際のプラントの特性と、モデル化した制御対象モデルの特性との差）に対してロバストな制御を実現し、スロットル弁開度偏差量ＤＴＨを目標値ＤＴＨＲに追従させるものである。
【００３４】
また式（５）の切換関数設定パラメータＶＰＯＬＥの値を変更することにより、図４に示すように、偏差ｅ（ｋ）の減衰特性、すなわちスロットル弁開度偏差量ＤＴＨの目標値ＤＴＨＲへの追従特性を変更することができる。具体的には、ＶＰＯＬＥ＝−１とすると、全く追従しない特性となり、切換関数設定パラメータＶＰＯＬＥの絶対値を小さくするほど、追従速度を速めることができる。
【００３５】
スロットル弁制御装置においては、下記要求Ａ１及びＡ２が満たされることが求められる。
Ａ１）スロットル弁３を全閉位置に移動させる際にスロットル全閉ストッパへの衝突を回避すること
Ａ２）デフォルト開度ＴＨＤＥＦ近傍における非線形特性（リターンスプリング４の付勢力と弾性部材５の付勢力とが釣り合うことに起因する弾性特性の変化、モータ６とスロットル弁３と間に介装されたギヤのバックラッシ、デューティ比ＤＵＴの変化してもスロットル弁開度が変化しない不感帯）に対する制御性を向上させること
そのため、スロットル弁の全閉位置近傍では、偏差ｅ（ｋ）の収束速度を低下させ、またデフォルト開度ＴＨＤＥＦの近傍では、収束速度を高める必要がある。
【００３６】
スライディングモード制御によれば、切換関数設定パラメータＶＰＯＬＥを変更することにより、容易に収束速度を変更できるので、本実施形態では、スロットル弁開度ＴＨ及び目標値ＤＴＨＲの変化量ＤＤＴＨＲ（＝ＤＴＨＲ（ｋ）−ＤＴＨＲ（ｋ−１））に応じて、切換関数設定パラメータＶＰＯＬＥを設定するようにした。これにより、上記要求Ａ１及びＡ２を満たすことができる。
【００３７】
上述したようにスライディングモード制御では、偏差ｅ（ｋ）と前回偏差ｅ（ｋ−１）の組み合わせ（以下「偏差状態量」という）を切換直線上に拘束することにより、偏差ｅ（ｋ）を指定した収束速度で、かつ外乱やモデル化誤差に対してロバストに、０に収束させる。したがって、スライディングモード制御では、如何にして偏差状態量を切換直線に載せ、そこに拘束するかが重要となる。
【００３８】
そのような観点から、制御対象への入力（コントローラの出力）ＤＵＴ（ｋ）（Ｕｓｌ（ｋ）とも表記する）は、下記式（６）に示すように、等価制御入力Ｕｅｑ（ｋ）、到達則入力Ｕｒｃｈ（ｋ）及び適応則入力Ｕａｄｐ（ｋ）の和として構成される。
等価制御入力Ｕｅｑ（ｋ）は、偏差状態量を切換直線上に拘束するための入力であり、到達則入力Ｕｒｃｈ（ｋ）は、偏差状態量を切換直線上へ載せるための入力であり、適応則入力Ｕａｄｐ（ｋ）は、モデル化誤差や外乱の影響を抑制し、偏差状態量を切換直線へ載せるための入力である。以下各入力Ｕｅｑ（ｋ），Ｕｒｃｈ（ｋ）及びＵａｄｐ（ｋ）の算出方法を説明する。
【００４０】
等価制御入力Ｕｅｑ（ｋ）は、偏差状態量を切換直線上に拘束するための入力であるから、満たすべき条件は下記式（７）で与えられる。
σ（ｋ）＝σ（ｋ＋１） （７）
式（１）並びに式（４）及び（５）を用いて式（７）を満たすデューティ比ＤＵＴ（ｋ）を求めると、下記式（９）が得られ、これが等価制御入力Ｕｅｑ（ｋ）となる。さらに、到達則入力Ｕｒｃｈ（ｋ）及び適応則入力Ｕａｄｐ（ｋ）を、それぞれ下記式（１０）及び（１１）により定義する。
【数１】
ここで、Ｆ及びＧは、それぞれ到達則制御ゲイン及び適応則制御ゲインであり、以下に述べるように設定される。またΔＴは、制御周期である。
上記式（９）の演算には、むだ時間ｄ経過後のスロットル弁開度偏差量ＤＴＨ（ｋ＋ｄ）及び対応する目標値ＤＴＨＲ（ｋ＋ｄ＋１）が必要である。そこで、むだ時間ｄ経過後のスロットル弁開度偏差量ＤＴＨ（ｋ＋ｄ）として、状態予測器２３により算出される予測偏差量ＰＲＥＤＴＨ（ｋ）を用い、目標値ＤＴＨＲ（ｋ＋ｄ＋１）として、最新の目標値ＤＴＨＲを用いることとする。
【００４２】
次に到達則入力Ｕｒｃｈ及び適応則入力Ｕａｄｐにより、偏差状態量が安定に切換直線上に載せられるように、到達則制御ゲインＦ及び適応則制御ゲインＧの決定を行う。
具体的には外乱Ｖ（ｋ）を想定し、外乱Ｖ（ｋ）に対して切換関数値σ（ｋ）が安定であるための条件を求めることにより、ゲインＦ及びＧの設定条件を求める。その結果、ゲインＦ及びＧの組み合わせが、下記式（１２）〜（１４）を満たすこと、換言すれば図５にハッチングを付して示す領域内にあることが安定条件として得られた。
【００４３】
Ｆ＞０ （１２）
Ｇ＞０ （１３）
Ｆ＜２−（ΔＴ／２）Ｇ （１４）
以上のように、式（９）〜（１１）により、等価制御入力Ｕｅｑ（ｋ）、到達則入力Ｕｒｃｈ（ｋ）及び適応則入力Ｕａｄｐ（ｋ）を算出し、それらの入力の総和として、デューティ比ＤＵＴ（ｋ）を算出することができる。
【００４４】
モデルパラメータ同定器２２は、前述したように制御対象の入力（ＤＵＴ（ｋ））及び出力（ＴＨ（ｋ））に基づいて、制御対象モデルのモデルパラメータベクトルを算出する。具体的には、モデルパラメータ同定器２２は、下記式（１５）による逐次型同定アルゴリズム（一般化逐次型最小２乗法アルゴリズム）により、モデルパラメータベクトルθ（ｋ）を算出する。
θ（ｋ）＝θ（ｋ−１）＋ＫＰ（ｋ）ｉｄｅ（ｋ） （１５）
θ（ｋ）Ｔ＝［ａ１’，ａ２’，ｂ１’，ｃ１’］ （１６）
【００４５】
ここで、ａ１’，ａ２’，ｂ１’及びＣ１’は、後述するリミット処理を実施する前のモデルパラメータである。またｉｄｅ（ｋ）は、下記式（１７）、（１８）及び（１９）により定義される同定誤差である。ＤＴＨＨＡＴ（ｋ）は、最新のモデルパラメータベクトルθ（ｋ−１）を用いて算出される、スロットル弁開度偏差量ＤＴＨ（ｋ）の推定値（以下「推定スロットル弁開度偏差量」という）である。ＫＰ（ｋ）は、下記式（２０）により定義されるゲイン係数ベクトルである。また、式（２０）のＰ（ｋ）は、下記式（２１）により算出される４次の正方行列である。
【数２】
式（２１）の係数λ１，λ２の設定により、式（１５）〜（２１）による同定アルゴリズムは、以下のような４つの同定アルゴリズムのいずれかになる。
λ１＝１，λ２＝０ 固定ゲインアルゴリズム
λ１＝１，λ２＝１ 最小２乗法アルゴリズム
λ１＝１，λ２＝λ 漸減ゲインアルゴリズム（λは０，１以外の所定値）
λ１＝λ，λ２＝１ 重み付き最小２乗法アルゴリズム（λは０，１以外の所定値）
【００４７】
一方本実施形態では、下記Ｂ１）、Ｂ２）、Ｂ３）の要求を満たすことが求められる。
Ｂ１）準静的動特性変化及びハードウエアの特性ばらつきに対する適応
「準静的動特性変化」とは、例えば電源電圧の変動やハードウエアの経年劣化といった変化速度の遅い特性変化を意味する。
Ｂ２）動的な動特性変化への適応
具体的には、スロットル弁開度ＴＨの変化に対応する動特性変化への適応を意味する。
Ｂ３）モデルパラメータのドリフト防止
モデルパラメータに反映すべきでない制御対象の非線形特性などに起因する同定誤差の影響によって、モデルパラメータの絶対値が増大するような不具合を防止する。
【００４８】
先ず上記Ｂ１）及びＢ２）の要求を満たすために、係数λ１及びλ２をそれぞれ所定値λ及び「０」に設定することにより、重み付き最小２乗法アルゴリズムを採用する。
次にモデルパラメータのドリフトについて説明する。図６に示すように、モデルパラメータがある程度収束した後に、スロットル弁の摩擦特性などの非線形特性によって生じる残留同定誤差が存在したり、平均値がゼロでない外乱が定常的に加わるような場合には、残留同定誤差が蓄積し、モデルパラメータのドリフトを引き起こす。
【００４９】
このような残留同定誤差は、モデルパラメータの値に反映すべきものではないので、図７（ａ）に示すような不感帯関数Ｆｎｌを用いて不感帯処理を行う。具体的には、下記式（２３）により、修正同定誤差ｉｄｅｎｌ（ｋ）を算出し、この修正同定誤差ｉｄｅｎｌ（ｋ）を用いてモデルパラメータベクトルθ（ｋ）の算出を行う。すなわち、上記式（１５）に代えて下記式（１５ａ）を用いる。これにより、上記要求Ｂ３）を満たすことができる。
ｉｄｅｎｌ（ｋ）＝Ｆｎｌ（ｉｄｅ（ｋ）） （２３）
θ（ｋ）＝θ（ｋ−１）＋ＫＰ（ｋ）ｉｄｅｎｌ（ｋ） （１５ａ）
【００５０】
なお、不感帯関数Ｆｎｌは、図７（ａ）に示すものに限るものではなく、例えば同図（ｂ）に示すような不連続不感帯関数、または同図（ｃ）に示すような不完全不感帯関数を用いてもよい。ただし、不完全不感帯関数を用いた場合には、ドリフトを完全に防止することはできない。
【００５１】
また、残留同定誤差は、スロットル弁開度ＴＨの変動量に応じてその振幅が変化する。そこで、本実施形態では、図７に示す不感帯の幅を定義する不感帯幅パラメータＥＩＤＮＲＬＭＴを、下記式（２４）により算出される、目標スロットル弁開度ＴＨＲの変化量の二乗平均値ＤＤＴＨＲＳＱＡに応じて設定する（具体的には、二乗平均値ＤＤＴＨＲＳＱＡが増加するほど、不感帯幅パラメータＥＩＤＮＲＬＭＴが増加するように設定する）ようにしている。これにより、モデルパラメータの値に反映させるべき同定誤差を、残留同定誤差として無視してしまうことを防止することができる。式（２４）のＤＤＴＨＲは、目標スロットル弁開度ＴＨＲの変化量であり、下記式（２５）により算出される。
【数４】
ここで、スロットル弁開度偏差量ＤＴＨは目標値ＤＴＨＲへ適応スライディングモードコントローラ２１により制御されているため、同様に式（２５）の目標値ＤＴＨＲをスロットル弁開度偏差量ＤＴＨへ変更し、スロットル弁開度偏差量ＤＴＨの変化量ＤＤＴＨを算出し、式（２４）のＤＤＴＨＲをＤＤＴＨに代えて得られる二乗平均値ＤＤＴＨＲＳＱＡにより不感帯幅パラメータＥＩＤＮＲＬＭＴを変更することもできる。
【００５３】
また制御系のロバスト性をさらに高めるためには、適応スライディングモードコントローラ２１をより安定化させることが有効である。そこで本実施形態では、前記式（１５）により算出されたモデルパラメータベクトルθ（ｋ）の各要素ａ１’，ａ２’，ｂ１’及びｃ１’についてリミット処理を施し、修正モデルパラメータベクトルθＬ（ｋ）（θＬ（ｋ）Ｔ＝［ａ１，ａ２，ｂ１，ｃ１］）を算出する。そして、適応スライディングモードコントローラ２１は、修正モデルパラメータベクトルθＬ（ｋ）を用いて、スライディングモード制御を実行する。なおリミット処理の詳細については、フローチャートを参照して後述する。
【００５４】
次に状態予測器２３による予測偏差量ＰＲＥＤＴＨの算出方法を説明する。
先ず下記式（２６）〜（２９）により、マトリクスＡ及びＢと、ベクトルＸ（ｋ）及びＵ（ｋ）を定義する。
【数５】
Ｘ（ｋ＋１）＝ＡＸ（ｋ）＋ＢＵ（ｋ−ｄ） （３０）
【００５５】
式（３０）からＸ（ｋ＋ｄ）を求めると、下記式（３１）が得られる。
【数６】
【数７】
予測ベクトルＸＨＡＴ（ｋ＋ｄ）の第１行の要素であるＤＴＨＨＡＴ（ｋ＋ｄ）が、予測偏差量ＰＲＥＤＴＨ（ｋ）であり、下記式（３５）で与えられる。
【００５７】
式（３５）により算出される予測偏差量ＰＲＥＤＴＨ（ｋ）を、前記式（９）に適用し、さらに目標値ＤＴＨＲ（ｋ＋ｄ＋１），ＤＴＨＲ（ｋ＋ｄ），及びＤＴＨＲ（ｋ＋ｄ−１）をそれぞれＤＴＨＲ（ｋ），ＤＴＨＲ（ｋ−１），及びＤＴＨＲ（ｋ−２）に置き換えることにより、下記式（９ａ）が得られる。式（９ａ）により、等価制御入力Ｕｅｑ（ｋ）を算出する。
【数８】
また、式（３５）により算出される予測偏差量ＰＲＥＤＴＨ（ｋ）を用いて、下記式（３６）により予測切替関数値σｐｒｅ（ｋ）を定義し、到達則入力Ｕｒｃｈ（ｋ）及び適応則入力Ｕａｄｐ（ｋ）を、それぞれ下記式（１０ａ）及び（１１ａ）により算出する。
次にモデルパラメータｃ１’は、前述したように、デフォルト開度ＴＨＤＥＦのずれ及び外乱を示すパラメータである。したがって、図８に示すように、外乱によって変動するが、デフォルト開度ずれは比較的短い期間内でみればほぼ一定とみなせる。そこで、本実施形態では、モデルパラメータｃ１’を統計処理し、その変動の中心値をデフォルト開度ずれｔｈｄｅｆａｄｐとして算出し、スロットル弁開度偏差量ＤＴＨ及び目標値ＤＴＨＲの算出に用いることとした。
【００６０】
統計処理の手法には、一般に最小２乗法が知られているが、この最小２乗法による統計処理は、通常、ある一定期間内のデータ、すなわち同定されたモデルパラメータｃ１’をすべてメモリに格納しておき、ある時点で一括演算を行うことによって実行される。ところが、この一括演算法では、すべてのデータを格納するために膨大な容量のメモリが必要となり、さらに逆行列演算が必要となって演算量の増大を招く。
【００６１】
そこで本実施形態では、前記式（１５）〜（２１）で示される適応制御の逐次型最小２乗法アルゴリズムを、統計処理に応用し、モデルパラメータｃ１の最小２乗中心値を、デフォルト開度ずれｔｈｄｅｆａｄｐとして算出するようにしている。
【００６２】
具体的には、前記式（１５）〜（２１）のθ（ｋ）及びθ（ｋ）Ｔをｔｈｄｅｆａｄｐに置換し、ζ（ｋ）及びζ（ｋ）Ｔを「１」に置換し、ｉｄｅ（ｋ）をｅｃ１（ｋ）に置換し、ＫＰ（ｋ）をＫＰＴＨ（ｋ）に置換し、Ｐ（ｋ）をＰＴＨ（ｋ）に置換し、λ１及びλ２をそれぞれλ１’及びλ２’に置換することにより、下記式（３７）〜（４０）を得る。
【数１０】
係数λ１’及びλ２’の設定により、前述した４つのアルゴリズムの何れかを選択可能であるが、式（３９）においては、係数λ１’を０または１以外の所定値に設定し、係数λ２’を１に設定することにより、重み付き最小２乗法を採用した。
【００６４】
上記式（３７）〜（４０）の演算においては、記憶すべき値はｔｈｄｅｆａｄｐ（ｋ＋１）及びＰＴＨ（ｋ＋１）のみであり、また逆行列演算は不要である。したがって、逐次型最小２乗法アルゴリズムを採用することにより、一般的な最小２乗法の欠点を克服しつつ、最小２乗法によるモデルパラメータｃ１の統計処理を行うことができる。
【００６５】
統計処理の結果得られるデフォルト開度ずれｔｈｄｅｆａｄｐは、前記式（２）及び（３）に適用され、式（２）及び（３）に代えて下記式（４１）及び（４２）により、スロットル弁開度偏差量ＤＴＨ（ｋ）及び目標値ＤＴＨＲ（ｋ）が算出される。
ＤＴＨ（ｋ）＝ＴＨ（ｋ）−ＴＨＤＥＦ＋ｔｈｄｅｆａｄｐ （４１）
ＤＴＨＲ（ｋ）＝ＴＨＲ（ｋ）−ＴＨＤＥＦ＋ｔｈｄｅｆａｄｐ （４２）
【００６６】
式（４１）及び（４２）を使用することにより、デフォルト開度ＴＨＤＥＦが、ハードウエアの特性ばらつき、あるいは経時変化により、設計値からずれた場合でも、そのずれを補償して正確な制御を行うことができる。
【００６７】
次に上述した適応スライディングモードコントローラ２１、モデルパラメータ同定器２２及び状態予測器２３の機能を実現するための、ＥＣＵ７のＣＰＵにおける演算処理を説明する。
【００６８】
図９は、スロットル弁開度制御の全体フローチャートであり、この処理は所定時間（例えば２ｍｓｅｃ）毎にＥＣＵ７のＣＰＵで実行される。
ステップＳ１１では、図１０に示す状態変数設定処理を実行する。すなわち、式（４１）及び（４２）の演算を実行し、スロットル弁開度偏差量ＤＴＨ（ｋ）及び目標値ＤＴＨＲ（ｋ）を算出する（図１０，ステップＳ２１及びＳ２２）。なお、今回値であることを示す（ｋ）は、省略して示す場合がある。
【００６９】
ステップＳ１２では、図１１に示すモデルパラメータ同定器の演算、すなわち前記式（１５ａ）によるモデルパラメータベクトルθ（ｋ）の算出処理を実行し、さらにリミット処理を実行して修正モデルパラメータベクトルθＬ（ｋ）を算出する。
【００７０】
続くステップＳ１３では、図２１に示す状態予測器の演算を実行し、予測偏差量ＰＲＥＤＴＨ（ｋ）を算出する。
次いでステップＳ１２で算出した修正モデルパラメータベクトルθＬ（ｋ）を用いて、図２２に示す制御入力Ｕｓｌ（ｋ）の演算処理を実行する（ステップＳ１４）。すなわち、等価制御入力Ｕｅｑ、到達則入力Ｕｒｃｈ（ｋ）及び適応則入力Ｕａｄｐ（ｋ）を算出し、それらの入力の総和として、制御入力Ｕｓｌ（ｋ）（＝デューティ比ＤＵＴ（ｋ））を算出する。
【００７１】
続くステップＳ１６では、図２９に示すスライディングモードコントローラの安定判別処理を実行する。すなわち、リアプノフ関数の微分値に基づく安定判別を行い、安定判別フラグＦＳＭＣＳＴＡＢの設定を行う。この安定判別フラグＦＳＭＣＳＴＡＢは、「１」に設定されると適応スライディングモードコントローラ２１が不安定となっていることを示す。安定判別フラグＦＳＭＣＳＴＡＢが「１」に設定され、適応スライディングモードコントローラ２１が不安定となったときは、切換関数設定パラメータＶＰＯＬＥを安定化所定値ＸＰＯＬＥＳＴＢに設定する（図２４、ステップＳ２３１，Ｓ２３２参照）とともに、等価制御入力Ｕｅｑを「０」とし、到達則入力Ｕｒｃｈ及び適応則入力Ｕａｄｐのみによる制御に切り換えることにより、制御の安定化を図る（図２２、ステップＳ２０６，Ｓ２０８参照）。適応スライディングモードコントローラ２１が不安定となったときは、さらに到達則入力Ｕｒｃｈ及び適応則入力Ｕａｄｐの算出式を変更する。すなわち、到達則制御ゲインＦ及び適応則制御ゲインＧの値を、コントローラ２１を安定化させる値に変更するとともに、モデルパラメータｂ１を使用しないで、到達則入力Ｕｒｃｈ及び適応則入力Ｕａｄｐを算出する（図２７，２８参照）。以上のような安定化処理により、適応スライディングモードコントローラ２１の不安定状態を早期に終息させ、安定な状態に戻すことができる。
ステップＳ１７では、図３０に示すｔｈｄｅｆａｄｐ算出処理を実行し、デフォルト開度ずれｔｈｄｅｆａｄｐを算出する。
【００７２】
図１１は、モデルパラメータ同定器２２の演算処理のフローチャートである。ステップＳ３１では、式（２０）によりゲイン係数ベクトルＫＰ（ｋ）を算出し、次いで式（１８）により推定スロットル弁開度偏差量ＤＴＨＨＡＴ（ｋ）を算出する（ステップＳ３２）。ステップＳ３３では、図１２に示すｉｄｅｎｌ（ｋ）の演算処理を実行し、ステップＳ３２で算出した推定スロットル弁開度偏差量ＤＴＨＨＡＴ（ｋ）を、式（１７）に適用して同定誤差ｉｄｅ（ｋ）を算出するとともに、図７（ａ）に示す関数を用いた不感帯処理を行い、修正同定誤差ｉｄｅｎｌを算出する。
【００７３】
続くステップＳ３４では、式（１５ａ）により、モデルパラメータベクトルθ（ｋ）を算出し、次いでモデルパラメータベクトルθ（ｋ）の安定化処理を実行する（ステップＳ３５）。すなわち各モデルパラメータのリミット処理を行って修正モデルパラメータベクトルθＬ（ｋ）を算出する。
【００７４】
図１２は、図１１のステップＳ３３で実行されるｉｄｅｎｌ（ｋ）演算処理のフローチャートである。
ステップＳ５１では、式（１７）により同定誤差ｉｄｅ（ｋ）を算出する。次いで、ステップＳ５３でインクリメントされるカウンタＣＮＴＩＤＳＴの値が、制御対象のむだ時間ｄに応じて設定される所定値ＸＣＮＴＩＤＳＴ（例えば、むだ時間ｄ＝２に対応して、「３」に設定される）より大きいか否かを判別する（ステップＳ５２）。カウンタＣＮＴＩＤＳＴの初期値は「０」であるので、最初はステップＳ５３に進み、カウンタＣＮＴＩＤＳＴを「１」だけインクリメントし、同定誤差ｉｄｅ（ｋ）を「０」に設定して（ステップＳ５４）、ステップＳ５５に進む。モデルパラメータベクトルθ（ｋ）の同定を開始した直後は、式（１７）による演算で正しい同定誤差が得られないので、ステップＳ５２〜Ｓ５４により、式（１７）による演算結果を用いずに同定誤差ｉｄｅ（ｋ）を「０」に設定するようにしている。
【００７５】
ステップＳ５２の答が肯定（ＹＥＳ）となると、直ちにステップＳ５５に進む。
ステップＳ５５では、同定誤差ｉｄｅ（ｋ）のローパスフィルタ処理を行う。具体的には、ローパス特性を有する制御対象のモデルパラメータを同定する場合、最小２乗同定アルゴリズムの同定誤差ｉｄｅ（ｋ）に対する同定重みは、図１３（ａ）に実線Ｌ１で示すような周波数特性を有するが、これをローパスフィルタ処理により、破線Ｌ２で示すように高周波成分を減衰させた特性とする。これは、以下の理由による。
【００７６】
実際の制御対象及びこれをモデル化した制御対象モデルの周波数特性は、それぞれ図１３（ｂ）に実線Ｌ３及びＬ４で示すようになる。すなわち、ローパス特性（高周波成分が減衰する特性）を有する制御対象について、モデルパラメータ同定器２２によりモデルパラメータを同定すると、同定されたモデルパラメータは高周波域阻止特性に大きく影響されたものとなるため、低周波域での制御対象モデルのゲインが実際の特性より低くなる。その結果、スライディングモードコントローラ２１による制御入力の補正が過補正となる。
【００７７】
そこで、ローパスフィルタ処理により同定アルゴリズムの重みの周波数特性を、図１３（ａ）に破線Ｌ２で示すような特性とすることにより、制御対象モデルの周波数特性を、同図（ｂ）に破線Ｌ５で示すような特性とし、実際の周波数特性と一致させ、あるいは制御対象モデルのゲインが実際のゲインよりやや高くなるように修正することとした。これにより、コントローラ２１による過補正を防止し、制御系のロバスト性を高めて制御系をより安定化させることができる。
【００７８】
なお、ローパスフィルタ処理は、同定誤差の過去値ｉｄｅ（ｋ−ｉ）（例えばｉ＝１〜１０に対応する１０個の過去値）をリングバッファに記憶し、それらの過去値に重み係数を乗算して加算することにより実行する。
さらに、同定誤差ｉｄｅ（ｋ）は、前記式（１７）、（１８）及び（１９）を用いて算出しているため、スロットル弁開度偏差量ＤＴＨ（ｋ）と、推定スロットル弁開度偏差量ＤＴＨＨＡＴ（ｋ）とに同様のローパスフィルタ処理を行うこと、あるいは、スロットル弁開度偏差量ＤＴＨ（ｋ−１）及びＤＴＨ（ｋ−２）と、デューティ比ＤＵＴ（ｋ−ｄ−１）とに同様のローパスフィルタ処理を行うことによっても同様の効果が得られる。
【００７９】
図１２に戻り、続くステップＳ５６では、図１４に示す不感帯処理を実行する。図１４のステップＳ６１では、前記式（２４）において例えばｎ＝５として、目標スロットル弁開度ＴＨＲの変化量の二乗平均値ＤＤＴＨＲＳＱＡを算出し、次いで二乗平均値ＤＤＴＨＲＳＱＡに応じて図１５に示すＥＩＤＮＲＬＭＴテーブルを検索し、不感帯幅パラメータＥＩＤＮＲＬＭＴを算出する（ステップＳ６２）。
【００８０】
ステップＳ６３では、同定誤差ｉｄｅ（ｋ）が不感帯幅パラメータＥＩＤＮＲＬＭＴより大きいか否かを判別し、ｉｄｅ（ｋ）＞ＥＩＤＮＲＬＭＴであるときは、下記式（４３）により、修正同定誤差ｉｄｅｎｌ（ｋ）算出する（ステップＳ６７）。
ｉｄｅｎｌ（ｋ）＝ｉｄｅ（ｋ）−ＥＩＤＮＲＬＭＴ （４３）
【００８１】
ステップＳ６３の答が否定（ＮＯ）であるときは、さらに同定誤差ｉｄｅ（ｋ）が不感帯幅パラメータＥＩＤＮＲＬＭＴに負号を付した値より小さいか否かを判別し（ステップＳ６４）、ｉｄｅ（ｋ）＜−ＥＩＤＮＲＬＭＴであるときは、下記式（４４）により、修正同定誤差ｉｄｅｎｌ（ｋ）算出する（ステップＳ６５）。
ｉｄｅｎｌ（ｋ）＝ｉｄｅ（ｋ）＋ＥＩＤＮＲＬＭＴ （４４）
また同定誤差ｉｄｅ（ｋ）が±ＥＩＤＮＲＬＭＴの範囲内にあるときは、修正同定誤差ｉｄｅｎｌ（ｋ）を「０」とする（ステップＳ６６）。
【００８２】
図１６は、図１１のステップＳ３５で実行されるθ（ｋ）の安定化処理のフローチャートである。
ステップＳ７１では、この処理で使用されるフラグＦＡ１ＳＴＡＢ，ＦＡ２ＳＴＡＢ，ＦＢ１ＬＭＴ及びＦＣ１ＬＭＴをそれぞれ「０」に設定することにより、初期化を行う。そして、ステップＳ７２では、図１７に示すａ１’及びａ２’のリミット処理を実行し、ステップＳ７３では、図１９に示すｂ１’のリミット処理を実行し、ステップＳ７４では、図２０に示すｃ１’のリミット処理を実行する。
【００８３】
図１７は、図１６のステップＳ７２で実行されるａ１’及びａ２’のリミット処理のフローチャートである。図１８は、図１７の処理を説明するための図であり、図１７とともに参照する。
図１８においては、リミット処理が必要なモデルパラメータａ１’とａ２’の組み合わせが「×」で示され、また安定なモデルパラメータａ１’及びａ２’の組み合わせの範囲がハッチングを付した領域（以下「安定領域」という）で示されている。図１７の処理は、安定領域外にあるモデルパラメータａ１’及びａ２’の組み合わせを、安定領域内（「○」で示す位置）に移動させる処理である。
【００８４】
ステップＳ８１では、モデルパラメータａ２’が、所定ａ２下限値ＸＩＤＡ２Ｌ以上か否かを判別する。所定ａ２下限値ＸＩＤＡ２Ｌは、「−１」より大きい負の値に設定される。所定ａ２下限値ＸＩＤＡ２Ｌは、「−１」に設定しても、安定な修正モデルパラメータａ１，ａ２が得られるが、前記式（２６）で定義される行列Ａのｎ乗が不安定となる（これは、ａ１’及びａ２’が発散はしないが振動することを意味する）場合があるので、「−１」より大きな値に設定される。
【００８５】
ステップＳ８１でａ２’＜ＸＩＤＡ２Ｌであるときは、修正モデルパラメータａ２を、この下限値ＸＩＤＡ２Ｌに設定するとともに、ａ２安定化フラグＦＡ２ＳＴＡＢを「１」に設定する。ａ２安定化フラグＦＡ２ＳＴＡＢは「１」に設定されると、修正モデルパラメータａ２を下限値ＸＩＤＡ２Ｌに設定したことを示す。図１８においては、ステップＳ８１及びＳ８２のリミット処理Ｐ１によるモデルパラメータの修正が、「Ｐ１」を付した矢線（矢印を付した線）で示されている。
【００８６】
ステップＳ８１の答が肯定（ＹＥＳ）、すなわちａ２’≧ＸＩＤＡ２Ｌであるときは、修正モデルパラメータａ２はモデルパラメータａ２’に設定される（ステップＳ８３）。
ステップＳ８４及びステップＳ８５では、モデルパラメータａ１’が、所定ａ１下限値ＸＩＤＡ１Ｌと所定ａ１上限値ＸＩＤＡ１Ｈできまる範囲内にあるか否かを判別する。所定ａ１下限値ＸＩＤＡ１Ｌは、−２以上且つ０より小さい値に設定され、所定ａ１上限値ＸＩＤＡ１Ｈは、例えば２に設定される。
【００８７】
ステップＳ８４及びＳ８５の答がいずれも肯定（ＹＥＳ）であるとき、すなわちＸＩＤＡ１Ｌ≦ａ１’≦ＸＩＤＡ１Ｈであるときは、修正モデルパラメータａ１はモデルパラメータａ１’に設定される（ステップＳ８８）。
一方ａ１’＜ＸＩＤＡ１Ｌであるときは、修正モデルパラメータａ１を下限値ＸＩＤＡ１Ｌに設定するとともに、ａ１安定化フラグＦＡ１ＳＴＡＢを「１」に設定する（ステップＳ８４，Ｓ８６）。またａ１’＞ＸＩＤＡ１Ｈであるときは、修正モデルパラメータａ１を上限値ＸＩＤＡ１Ｈに設定するとともに、ａ１安定化フラグＦＡ１ＳＴＡＢを「１」に設定する（ステップＳ８５，Ｓ８７）。ａ１安定化フラグＦＡ１ＳＴＡＢは、「１」に設定されると、修正モデルパラメータａ１を下限値ＸＩＤＡ１Ｌまたは上限値ＸＩＤＡ１Ｈに設定したことを示す。図１８においては、ステップＳ８４〜Ｓ８７のリミット処理Ｐ２によるモデルパラメータの修正が、「Ｐ２」を付した矢線で示されている。
【００８８】
ステップＳ９０では、修正モデルパラメータａ１の絶対値と修正モデルパラメータａ２の和が、所定安定判定値ＸＡ２ＳＴＡＢ以下であるか否かを判別する。所定安定判定値ＸＡ２ＳＴＡＢは、「１」に近く「１」より小さい値（例えば０．９９）に設定される。
【００８９】
図１８に示す直線Ｌ１及びＬ２は、下記式（４５）を満たす直線である。
ａ２＋｜ａ１｜＝ＸＡ２ＳＴＡＢ （４５）
したがって、ステップＳ９０は、修正モデルパラメータａ１及びａ２の組み合わせが、図１８に示す直線Ｌ１及びＬ２の線上またはその下側にあるか否かを判別している。ステップＳ９０の答が肯定（ＹＥＳ）であるときは、修正モデルパラメータａ１及びａ２の組み合わせは、図１８の安定領域内にあるので、直ちに本処理を終了する。
【００９０】
一方ステップＳ９０の答が否定（ＮＯ）であるときは、修正モデルパラメータａ１が、所定安定判定値ＸＡ２ＳＴＡＢから所定ａ２下限値ＸＩＤＡ２Ｌを減算した値（ＸＩＤＡ２Ｌ＜０であるので、ＸＡ２ＳＴＡＢ−ＸＩＤＡ２Ｌ＞ＸＡ２ＳＴＡＢが成立する）以下か否かを判別する（ステップＳ９１）。そして修正モデルパラメータａ１が（ＸＡ２ＳＴＡＢ−ＸＩＤＡ２Ｌ）以下であるときは、修正モデルパラメータａ２を（ＸＡ２ＳＴＡＢ−｜ａ１｜）に設定するとともに、ａ２安定化フラグＦＡ２ＳＴＡＢを「１」に設定する（ステップＳ９２）。
【００９１】
ステップＳ９１で修正モデルパラメータａ１が（ＸＡ２ＳＴＡＢ−ＸＩＤＡ２Ｌ）より大きいときは、修正モデルパラメータａ１を（ＸＡ２ＳＴＡＢ−ＸＩＤＡ２Ｌ）に設定し、修正モデルパラメータａ２を所定ａ２下限値ＸＩＤＡ２Ｌに設定するとともに、ａ１安定化フラグＦＡ１ＳＴＡＢ及びａ２安定化フラグＦＡ２ＳＴＡＢをともに「１」に設定する（ステップＳ９３）。
【００９２】
図１８においては、ステップＳ９１及びＳ９２のリミット処理Ｐ３によるモデルパラメータの修正が、「Ｐ３」を付した矢線で示されており、またステップＳ９１及びＳ９３のリミット処理Ｐ４によるモデルパラメータの修正が、「Ｐ４」を付した矢線で示されている。
以上のように図１７の処理により、モデルパラメータａ１’及びａ２’が図１８に示す安定領域内に入るようにリミット処理が実行され、修正モデルパラメータａ１及びａ２が算出される。
【００９３】
図１９は、図１６のステップＳ７３で実行されるｂ１’のリミット処理のフローチャートである。
ステップＳ１０１及びＳ１０２では、モデルパラメータｂ１’が、所定ｂ１下限値ＸＩＤＢ１Ｌと所定ｂ１上限値ＸＩＤＢ１Ｈできまる範囲内にあるか否かを判別する。所定ｂ１下限値ＸＩＤＢ１Ｌは、正の所定値（例えば０．１）に設定され、所定ｂ１上限値ＸＩＤＢ１Ｈは、例えば「１」に設定される。
【００９４】
ステップＳ１０１及びＳ１０２の答がいずれも肯定（ＹＥＳ）であるとき、すなわちＸＩＤＢ１Ｌ≦ｂ１’≦ＸＩＤＢ１Ｈであるときは、修正モデルパラメータｂ１はモデルパラメータｂ１’に設定される（ステップＳ１０５）。
一方ｂ１’＜ＸＩＤＢ１Ｌであるときは、修正モデルパラメータｂ１を下限値ＸＩＤＢ１Ｌに設定するとともに、ｂ１リミットフラグＦＢ１ＬＭＴを「１」に設定する（ステップＳ１０１，Ｓ１０４）。またｂ１’＞ＸＩＤＢ１Ｈであるときは、修正モデルパラメータｂ１を上限値ＸＩＤＢ１Ｈに設定するとともに、ｂ１リミットフラグＦＢ１ＬＭＴを「１」に設定する（ステップＳ１０２，Ｓ１０３）。ｂ１リミットフラグＦＢ１ＬＭＴは、「１」に設定されると、修正モデルパラメータｂ１を下限値ＸＩＤＢ１Ｌまたは上限値ＸＩＤＢ１Ｈに設定したことを示す。
【００９５】
図２０は、図１６のステップＳ７４で実行されるモデルパラメータｃ１’のリミット処理のフローチャートである。
ステップＳ１１１及びＳ１１２では、モデルパラメータｃ１’が、所定ｃ１下限値ＸＩＤＣ１Ｌと所定ｃ１上限値ＸＩＤＣ１Ｈできまる範囲内にあるか否かを判別する。所定ｃ１下限値ＸＩＤＣ１Ｌは、例えば−６０に設定され、所定ｃ１上限値ＸＩＤＣ１Ｈは、例えば６０に設定される。
【００９６】
ステップＳ１１１及びＳ１１２の答がいずれも肯定（ＹＥＳ）であるとき、すなわちＸＩＤＣ１Ｌ≦ｃ１’≦ＸＩＤＣ１Ｈであるときは、修正モデルパラメータｃ１はモデルパラメータｃ１’に設定される（ステップＳ１１５）。
一方ｃ１’＜ＸＩＤＣ１Ｌであるときは、修正モデルパラメータｃ１を下限値ＸＩＤＣ１Ｌに設定するとともに、ｃ１リミットフラグＦＣ１ＬＭＴを「１」に設定する（ステップＳ１１１，Ｓ１１４）。またｃ１’＞ＸＩＤＣ１Ｈであるときは、修正モデルパラメータｃ１を上限値ＸＩＤＣ１Ｈに設定するとともに、ｃ１リミットフラグＦＣ１ＬＭＴを「１」に設定する（ステップＳ１１２，Ｓ１１３）。ｃ１リミットフラグＦＣ１ＬＭＴは、「１」に設定されると、修正モデルパラメータｃ１を下限値ＸＩＤＣ１Ｌまたは上限値ＸＩＤＣ１Ｈに設定したことを示す。
【００９７】
図２１は、図９のステップＳ１３で実行される状態予測器の演算処理のフローチャートである。
ステップＳ１２１では、マトリクス演算を実行して前記式（３５）の行列要素α１，α２，β１〜β２、及びγ１〜γｄを算出する。
ステップＳ１２２では、式（３５）により、予測偏差量ＰＲＥＤＴＨ（ｋ）を算出する。
【００９８】
図２２は、図９のステップＳ１４で実行される、スロットル弁駆動装置１０への制御入力Ｕｓｌ（＝ＤＵＴ）を算出する処理のフローチャートである。
ステップＳ２０１では、図２３に示す予測切換関数値σｐｒｅの演算処理を実行し、ステップＳ２０２では、図２６に示す予測切換関数値σｐｒｅの積算値の演算処理を実行する。ステップＳ２０３では、前記式（９）により、等価制御入力Ｕｅｑを算出する。ステップＳ２０４では、図２７に示す到達則入力Ｕｒｃｈの演算処理を実行し、ステップＳ２０５では、図２８に示す適応則入力Ｕａｄｐの演算処理を実行する。
【００９９】
ステップＳ２０６では、後述する図２９の処理で設定される安定判別フラグＦＳＭＣＳＴＡＢが「１」であるか否かを判別する。安定判別フラグＦＳＭＣＳＴＡＢは、「１」に設定されると、適応スライディングモードコントローラ２１が不安定となっていることを示す。
【０１００】
ステップＳ２０６でＦＳＭＣＳＴＡＢ＝０であって適応スライディングモードコントローラ２１が安定であるときは、ステップＳ２０３〜Ｓ２０５で算出された制御入力Ｕｅｑ，Ｕｒｃｈ及びＵａｄｐを加算することにより、制御入力Ｕｓｌを算出する（ステップＳ２０７）。
【０１０１】
一方ＦＳＭＣＳＴＡＢ＝１であって適応スライディングモードコントローラ２１が不安定となっているときは、到達則入力Ｕｒｃｈ及び適応則入力Ｕａｄｐの和を、制御入力Ｕｓｌとして算出する。すなわち、等価制御入力Ｕｅｑを、制御入力Ｕｓｌの算出に使用しないようにする。これにより、制御系が不安定化することを防止することができる。
【０１０２】
続くステップＳ２０９及びＳ２１０では、算出した制御入力Ｕｓｌが所定上下限値ＸＵＳＬＨ及びＸＵＳＬＬの範囲内にあるか否かを判別し、制御入力Ｕｓｌが所定上下限値の範囲内にあるときは、直ちに本処理を終了する。一方、制御入力Ｕｓｌが所定下限値ＸＵＳＬＬ以下であるときは、制御入力Ｕｓｌを所定下限値ＸＵＳＬＬに設定し（ステップＳ２０９，Ｓ２１２）、制御入力Ｕｓｌが所定上限値ＸＵＳＬＨ以上であるときは、制御入力Ｕｓｌを所定上限値ＸＵＳＬＨに設定する（ステップＳ２１０，Ｓ２１１）。
【０１０３】
図２３は、図２２のステップＳ２０１で実行される予測切換関数値σｐｒｅの演算処理のフローチャートである。
ステップＳ２２１では、図２４に示す切換関数設定パラメータＶＰＯＬＥの演算処理を実行し、次いで前記式（３６）により、予測切換関数値σｐｒｅ（ｋ）の演算を実行する（ステップＳ２２２）。
【０１０４】
続くステップＳ２２３及びＳ２２４では、算出した予測切換関数値σｐｒｅ（ｋ）が所定上下限値ＸＳＧＭＨ及びＸＳＧＭＬの範囲内にあるか否かを判別し、予測切換関数値σｐｒｅ（ｋ）が所定上下限値の範囲内にあるときは、直ちに本処理を終了する。一方、予測切換関数値σｐｒｅ（ｋ）が所定下限値ＸＳＧＭＬ以下であるときは、予測切換関数値σｐｒｅ（ｋ）を所定下限値ＸＳＧＭＬに設定し（ステップＳ２２３，Ｓ２２５）、予測切換関数値σｐｒｅ（ｋ）が所定上限値ＸＳＧＭＨ以上であるときは、予測切換関数値σｐｒｅ（ｋ）を所定上限値ＸＳＧＭＨに設定する（ステップＳ２２４，Ｓ２２６）。
【０１０５】
図２４は、図２３のステップＳ２２１で実行される切換関数設定パラメータＶＰＯＬＥの演算処理のフローチャートである。
ステップＳ２３１では、安定判別フラグＦＳＭＣＳＴＡＢが「１」であるか否かを判別し、ＦＳＭＣＳＴＡＢ＝１であって適応スライディングモードコントローラ２１が不安定となっているときは、切換関数設定パラメータＶＰＯＬＥを安定化所定値ＸＰＯＬＥＳＴＢに設定する（ステップＳ２３２）。安定化所定値ＸＰＯＬＥＳＴＢは、「−１」より大きく「−１」に非常に近い値（例えば−０．９９９）に設定される。
【０１０６】
ＦＳＭＣＳＴＡＢ＝０であって適応スライディングモードコントローラ２１が安定であるときは、下記式（４６）により目標値ＤＴＨＲ（ｋ）の変化量ＤＤＴＨＲ（ｋ）を算出する（ステップＳ２３３）。
ＤＤＴＨＲ（ｋ）＝ＤＴＨＲ（ｋ）−ＤＴＨＲ（ｋ−１） （４６）
ステップＳ２３４では、スロットル弁開度偏差量ＤＴＨ及びステップＳ２３３で算出される目標値の変化量ＤＤＴＨＲに応じてＶＰＯＬＥマップを検索し、切換関数設定パラメータＶＰＯＬＥを算出する。ＶＰＯＬＥマップは、図２５（ａ）に示すように、スロットル弁開度偏差量ＤＴＨが０近傍の値をとるとき（スロットル弁開度ＴＨがデフォルト開度ＴＨＤＥＦ近傍の値をとるとき）増加し、０近傍以外の値ではスロットル弁開度偏差量ＤＴＨの変化に対してはほぼ一定の値となるように設定されている。またＶＰＯＬＥマップは、同図（ｂ）に実線で示すように、目標値の変化量ＤＤＴＨＲが増加するほど、ＶＰＯＬＥ値が増加するように設定されているが、スロットル弁開度偏差量ＤＴＨが０近傍の値をとるときには、同図に破線で示すように目標値の変化量ＤＤＴＨＲが０近傍の値をとるときに増加するように設定されている。
【０１０７】
すなわち、スロットル弁開度の目標値ＤＴＨＲが減少方向の変化が大きいときには、切換関数設定パラメータＶＰＯＬＥは比較的小さな値に設定される。これにより、スロットル弁３がスロットル全閉ストッパに衝突することを防止することができる。また、デフォルト開度ＴＨＤＥＦ近傍においては、切換関数設定パラメータＶＰＯＬＥが比較的大きな値に設定され、デフォルト開度ＴＨＤＥＦ近傍における制御性を向上させることができる。
【０１０８】
なお、同図（ｃ）に示すように、スロットル弁開度ＴＨが全閉開度近傍または全開開度近傍にあるときは、切換関数設定パラメータＶＰＯＬＥを減少させるように設定してもよい。これにより、スロットル弁開度ＴＨが全閉開度近傍または全開開度近傍にあるときは、目標開度ＴＨＲに対する追従速度が遅くなり、スロットル弁３の全閉ストッパ（全開開度でもストッパとして機能する）への衝突防止をより確実にすることができる。
【０１０９】
続くステップＳ２３５及びＳ２３６では、算出した切換関数設定パラメータＶＰＯＬＥが所定上下限値ＸＰＯＬＥＨ及びＸＰＯＬＥＬの範囲内にあるか否かを判別し、切換関数設定パラメータＶＰＯＬＥが所定上下限値の範囲内にあるときは、直ちに本処理を終了する。一方、切換関数設定パラメータＶＰＯＬＥが所定下限値ＸＰＯＬＥＬ以下であるときは、切換関数設定パラメータＶＰＯＬＥを所定下限値ＸＰＯＬＥＬに設定し（ステップＳ２３６，Ｓ２３８）、切換関数設定パラメータＶＰＯＬＥが所定上限値ＸＰＯＬＥＨ以上であるときは、切換関数設定パラメータＶＰＯＬＥを所定上限値ＸＰＯＬＥＨに設定する（ステップＳ２３５，Ｓ２３７）。
【０１１０】
図２６は、図２２のステップＳ２０２で実行される、予測切換関数値σｐｒｅの積算値ＳＵＭＳＩＧＭＡを算出する処理のフローチャートである。積算値ＳＵＭＳＩＧＭＡは、後述する図２８の処理で適応則入力Ｕａｄｐの算出に使用される（前記式（１１ａ）参照）。
【０１１１】
ステップＳ２４１では、下記式（４７）により、積算値ＳＵＭＳＩＧＭＡを算出する。下記式のΔＴは、演算の実行周期である。
ＳＵＭＳＩＧＭＡ（ｋ）＝ＳＵＭＳＩＧＭＡ（ｋ−１）＋σｐｒｅ×ΔＴ （４７）
続くステップＳ２４２及びＳ２４３では、算出した積算値ＳＵＭＳＩＧＭＡが所定上下限値ＸＳＵＭＳＨ及びＸＳＵＭＳＬの範囲内にあるか否かを判別し、積算値ＳＵＭＳＩＧＭＡが所定上下限値の範囲内にあるときは、直ちに本処理を終了する。一方、積算値ＳＵＭＳＩＧＭＡが所定下限値ＸＳＵＭＳＬ以下であるときは、積算値ＳＵＭＳＩＧＭＡを所定下限値ＸＳＵＭＳＬに設定し（ステップＳ２４２，Ｓ２４４）、積算値ＳＵＭＳＩＧＭＡが所定上限値ＸＳＵＭＳＨ以上であるときは、積算値ＳＵＭＳＩＧＭＡを所定上限値ＸＳＵＭＳＨに設定する（ステップＳ２４３，Ｓ２４５）。
【０１１２】
図２７は、図２２のステップＳ２０４で実行される到達則入力Ｕｒｃｈの演算処理のフローチャートである。
ステップＳ２６１では、安定判別フラグＦＳＭＣＳＴＡＢが「１」であるか否かを判別する。安定判別フラグＦＳＭＣＳＴＡＢが「０」であって適応スライディングモードコントローラ２１が安定であるときは、制御ゲインＦを所定通常ゲインＸＫＲＣＨに設定し（ステップＳ２６２）、下記式（４８）（前記式（１０ａ）と同一の式）により、到達則入力Ｕｒｃｈを算出する（ステップＳ２６３）。
Ｕｒｃｈ＝−Ｆ×σｐｒｅ／ｂ１ （４８）
【０１１３】
一方安定判別フラグＦＳＭＣＳＴＡＢが「１」であって適応スライディングモードコントローラ２１が不安定となったときは、制御ゲインＦを、所定安定化ゲインＸＫＲＣＨＳＴＢに設定し（ステップＳ２６４）、モデルパラメータｂ１を使わない下記式（４９）により到達則入力Ｕｒｃｈを算出する（ステップＳ２６５）。
Ｕｒｃｈ＝−Ｆ×σｐｒｅ （４９）
【０１１４】
続くステップＳ２６６及びＳ２６７では、算出した到達則入力Ｕｒｃｈが所定上下限値ＸＵＲＣＨＨ及びＸＵＲＣＨＬの範囲内にあるか否かを判別し、到達則入力Ｕｒｃｈが所定上下限値の範囲内にあるときは、直ちに本処理を終了する。一方、到達則入力Ｕｒｃｈが所定下限値ＸＵＲＣＨＬ以下であるときは、到達則入力Ｕｒｃｈを所定下限値ＸＵＲＣＨＬに設定し（ステップＳ２６６，Ｓ２６８）、到達則入力Ｕｒｃｈが所定上限値ＸＵＲＣＨＨ以上であるときは、到達則入力Ｕｒｃｈを所定上限値ＸＵＲＣＨＨに設定する（ステップＳ２６７，Ｓ２６９）。
【０１１５】
このように適応スライディングモードコントローラ２１が不安定となったときは、制御ゲインＦを所定安定化ゲインＸＫＲＣＨＳＴＢに設定するとともに、モデルパラメータｂ１を使用しないで到達則入力Ｕｒｃｈを算出することにより、適応スライディングモードコントローラ２１を安定な状態に戻すことができる。モデルパラメータ同定器２２による同定が不安定となった場合に、適応スライディングモードコントローラ２１が不安定となるので、不安定となったモデルパラメータｂ１を使わないことによって、適応スライディングモードコントローラ２１を安定化することができる。
【０１１６】
図２８は、図２２のステップＳ２０５で実行される適応則入力Ｕａｄｐの演算処理のフローチャートである。
ステップＳ２７１では、安定判別フラグＦＳＭＣＳＴＡＢが「１」であるか否かを判別する。安定判別フラグＦＳＭＣＳＴＡＢが「０」であって適応スライディングモードコントローラ２１が安定であるときは、制御ゲインＧを所定通常ゲインＸＫＡＤＰに設定し（ステップＳ２７２）、下記式（５０）（前記式（１１ａ）に対応する式）により、適応則入力Ｕａｄｐを算出する（ステップＳ２７３）。
Ｕａｄｐ＝−Ｇ×ＳＵＭＳＩＧＭＡ／ｂ１ （５０）
【０１１７】
一方安定判別フラグＦＳＭＣＳＴＡＢが「１」であって適応スライディングモードコントローラ２１が不安定となったときは、制御ゲインＧを、所定安定化ゲインＸＫＡＤＰＳＴＢに設定し（ステップＳ２７４）、モデルパラメータｂ１を使わない下記式（５１）により適応則入力Ｕａｄｐを算出する（ステップＳ２７５）。
Ｕａｄｐ＝−Ｇ×ＳＵＭＳＩＧＭＡ （５１）
【０１１８】
このように適応スライディングモードコントローラ２１が不安定となったときは、制御ゲインＧを所定安定化ゲインＸＫＡＤＰＳＴＢに設定するとともに、モデルパラメータｂ１を使用しないで適応則入力Ｕａｄｐを算出することにより、適応スライディングモードコントローラ２１を安定な状態に戻すことができる。
【０１１９】
図２９は、図９のステップＳ１６で実行されるスライディングモードコントローラの安定判別処理のフローチャートである。この処理では、リアプノフ関数の微分項に基づく安定判別を行い、安定判別フラグＦＳＭＣＳＴＡＢの設定を行う。
【０１２０】
ステップＳ２８１では下記式（５２）により、切換関数変化量Ｄσｐｒｅを算出し、次いで下記式（５３）により、安定性判別パラメータＳＧＭＳＴＡＢを算出する（ステップＳ２８２）。
Ｄσｐｒｅ＝σｐｒｅ（ｋ）−σｐｒｅ（ｋ−１） （５２）
ＳＧＭＳＴＡＢ＝Ｄσｐｒｅ×σｐｒｅ（ｋ） （５３）
ステップＳ２８３では、安定性判別パラメータＳＧＭＳＴＡＢが安定性判定閾値ＸＳＧＭＳＴＡＢ以下か否かを判別し、ＳＧＭＳＴＡＢ＞ＸＳＧＭＳＴＡＢであるときは、コントローラ２１が不安定である可能性があると判定して不安定検知カウンタＣＮＴＳＭＣＳＴを「１」だけインクリメントする（ステップＳ２８５）。また、ＳＧＭＳＴＡＢ≦ＸＳＧＭＳＴＡＢであるときは、コントローラ２１が安定であると判定して不安定検知カウンタＣＮＴＳＭＣＳＴのカウント値をインクリメントすることなく保持する（ステップＳ２８４）。
【０１２１】
ステップＳ２８６では、不安定検知カウンタＣＮＴＳＭＣＳＴの値が所定カウント値ＸＳＳＴＡＢ以下か否かを判別する。ＣＮＴＳＭＣＳＴ≦ＸＳＳＴＡＢであるときは、コントローラ２１は安定していると判定し、第１判定フラグＦＳＭＣＳＴＡＢ１を「０」に設定する（ステップＳ２８７）。一方ＣＮＴＳＭＣＳＴ＞ＸＳＳＴＡＢであるときは、コントローラ２１は不安定となっていると判定し、第１判定フラグＦＳＭＣＳＴＡＢ１を「１」に設定する（ステップＳ２８８）。なお、不安定検知カウンタＣＮＴＳＭＣＳＴは、イグニッションスイッチオン時にそのカウント値が「０」に初期化される。
【０１２２】
続くステップＳ２８９では、安定判別期間カウンタＣＮＴＪＵＤＳＴを「１」だけデクリメントし、次いでその安定判別期間カウンタＣＮＴＪＵＤＳＴの値が「０」であるか否かを判別する（ステップＳ２９０）。安定判別期間カウンタＣＮＴＪＵＤＳＴは、イグニッションスイッチオン時に所定判別カウント値ＸＣＪＵＤＳＴに初期化される。したがって、最初はステップＳ２９０の答は否定（ＮＯ）となり、直ちにステップＳ２９５に進む。
【０１２３】
その後安定判別期間カウンタＣＮＴＪＵＤＳＴが「０」となると、ステップＳ２９０からステップＳ２９１に進み、第１判定フラグＦＳＭＣＳＴＡＢ１が「１」であるか否かを判別する。そして、第１判定フラグＦＳＭＣＳＴＡＢ１が「０」であるときは、第２判定フラグＦＳＭＣＳＴＡＢ２を「０」に設定し（ステップＳ２９３）、第１判定フラグＦＳＭＣＳＴＡＢ１が「１」であるときは、第２判定フラグＦＳＭＣＳＴＡＢ２を「１」に設定する（ステップＳ２９２）。
【０１２４】
続くステップＳ２９４では、安定判別期間カウンタＣＮＴＪＵＤＳＴの値を所定判別カウント値ＸＣＪＵＤＳＴに設定するとともに、不安定検知カウンタＣＮＴＳＭＣＳＴの値を「０」に設定し、ステップＳ２９５に進む。
ステップＳ２９５では、安定判別フラグＦＳＭＣＳＴＡＢを、第１判定フラグＦＳＭＣＳＴＡＢ１と第２判定フラグＦＳＭＣＳＴＡＢ２の論理和に設定する。第２判定フラグＦＳＭＣＳＴＡＢ２は、ステップＳ２８６の答が肯定（ＹＥＳ）となり、第１判定フラグＦＳＭＣＳＴＡＢ１が「０」に設定されても、安定判別期間カウンタＣＮＴＪＵＤＳＴの値が「０」となるまでは、「１」に維持される。したがって、安定判別フラグＦＳＭＣＳＴＡＢも、安定判別期間カウンタＣＮＴＪＵＤＳＴの値が「０」となるまでは、「１」に維持される。
【０１２５】
図３０は、図９のステップＳ１７で実行されるデフォルト開度ずれｔｈｄｅｆａｄｐの算出処理のフローチャートである。
ステップＳ２５１では、下記式（５４）により、ゲイン係数ＫＰＴＨ（ｋ）を算出する。
ＫＰＴＨ（ｋ）＝ＰＴＨ（ｋ−１）／（１＋ＰＴＨ（ｋ−１）） （５４）
【０１２６】
ここでＰＴＨ(k-1)は、本処理の前回実行時にステップＳ２５３で算出されたゲインパラメータである。
ステップＳ２５２では、図１１に示すモデルパラメータ同定器演算処理で算出されるモデルパラメータｃ１’及びステップＳ２５１で算出したゲイン係数ＫＰＴＨ(k)を下記式（５５）に適用し、デフォルト開度ずれｔｈｄｅｆａｄｐ(k)を算出する。
ステップＳ２５３では、下記式（５６）によりゲインパラメータＰＴＨ（ｋ）を算出する。
【０１２８】
図３０の処理により、モデルパラメータｃ１’が逐次型重み付き最小２乗法により統計処理され、デフォルト開度ずれｔｈｄｅｆａｄｐが算出される。
本実施形態では、スロットル弁駆動装置１０及びＥＣＵ７の一部（モータ６に駆動電流を供給する出力回路）がプラントに相当し、図２２の処理が制御手段に相当し、図１１の処理が同定手段に相当し、図１２の処理が同定誤差算出手段に相当し、図１４の処理が同定誤差修正手段に相当し、図１１のステップＳ３４がモデルパラメータベクトル算出手段に相当し、図１６の処理が制限手段に相当する。
【０１２９】
（第２の実施形態）
上述した第１の実施形態では、制御対象モデルをむだ時間ｄを含む式（１）を用いて定義し、状態予測器２３を用いて、むだ時間ｄ経過後の予測偏差量ＰＲＥＤＴＨを算出することにより、むだ時間を含む制御対象モデルの制御を行っている。そのため、状態予測器２３に対応した演算をＣＰＵで実行する必要があり、ＣＰＵの演算量が大きくなる。そこで、本実施形態では、ＣＰＵに加わる演算負荷の軽減を図るべく、むだ時間ｄを「０」とした下記式（１ａ）により制御対象モデルを定義し、むだ時間ｄを「０」とすることに起因するモデル化誤差は、適応スライディングモード制御のロバスト性によって補償している。
ＣＰＵの演算負荷をさらに軽減するために、モデルパラメータの同定アルゴリズムとして、固定ゲインアルゴリズムを採用している。
また、より一層の制御の安定化を図るべく、モデルパラメータのドリフトを防止する手法として、不感帯処理に代わる他の手法を採用している。
以下本実施形態を、第１の実施形態と異なる点を中心として詳細に説明する。以下に述べる点以外は、第１の実施形態と同一である。
【０１３１】
図３１は、ＥＣＵ７により実現されるスロットル弁制御装置の機能ブロック図であり、この制御装置は、適応スライディングモードコントローラ２１ａと、モデルパラメータ同定器２２ａと、モデルパラメータスケジューラ２５と、アクセルペダル踏み込み量ＡＣＣに応じてスロットル弁３の目標開度ＴＨＲを設定する目標開度設定部２４とからなる。
【０１３２】
適応スライディングモードコントローラ２１ａには、予測偏差量ＰＲＥＤＴＨではなく、検出したスロットル弁開度ＴＨが入力され、このスロットル弁開度ＴＨが目標開度ＴＨＲと一致するように、適応スライディングモード制御によりデューティ比ＤＵＴが算出される。
【０１３３】
適応スライディングモードコントローラ２１ａを用いることにより、第１の実施形態で説明した効果と同様の効果が得られ、また制御対象のむだ時間に対する制御系のロバスト性を確保することができる。したがって、むだ時間ｄを「０」とすることに起因するモデル化誤差を補償することができる。
【０１３４】
モデルパラメータ同定器２２ａは、第１の実施形態とは異なる手法で、修正モデルパラメータベクトルθＬ（θＬＴ＝［ａ１，ａ２，ｂ１，ｃ１］）を算出し、適応スライディングモードコントローラ２１ａに供給する。より具体的には、モデルパラメータ同定器２２ａは、モデルパラメータスケジューラ２５から供給される基準モデルパラメータベクトルθｂａｓｅを、スロットル弁開度ＴＨ及びデューティ比ＤＵＴに基づいて補正することにより、モデルパラメータベクトルθを算出する。さらに、そのモデルパラメータベクトルθに対してリミット処理を行うことにより修正モデルパラメータベクトルθＬを算出し、該修正モデルパラメータベクトルθＬを適応スライディングモードコントローラ２１ａに供給する。このようにしてスロットル弁開度ＴＨを目標開度ＴＨＲに追従させるために最適なモデルパラメータａ１，ａ２，ｂ１が得られ、さらに外乱及びデフォルト開度ＴＨＤＥＦのずれを示すモデルパラメータｃ１が得られる。
【０１３５】
モデルパラメータスケジューラ２５は、スロットル弁開度ＴＨに基づいて、基準モデルパラメータベクトルθｂａｓｅ（θｂａｓｅＴ＝［ａ１ｂａｓｅ，ａ２ｂａｓｅ，ｂ１ｂａｓｅ，ｃ１ｂａｓｅ］）を算出し、モデルパラメータ同定器２２ａに供給する。
【０１３６】
本実施形態では、前記式（１ａ）により制御対象モデルを定義しているので、適応スライディングモードコントローラ２１ａは、等価制御入力Ｕｅｑ、到達則入力Ｕｒｃｈ及び適応則入力Ｕａｄｐを、前記式（９ａ），（１０ａ），（１１ａ）に代えて、下記式（９ｂ），（１０ｂ），（１１ｂ）により算出する。
【数１１】
式（９ｂ）〜（１１ｂ）は、前記式（９）〜（１１）のむだ時間ｄを「０」とすることにより得られる。
モデルパラメータ同定器２２ａは、前述したように制御対象の入力（ＤＵＴ（ｋ））及び出力（ＴＨ（ｋ））に基づいて、制御対象モデルのモデルパラメータベクトルを算出する。具体的には、モデルパラメータ同定器２２ａは、下記式（１５）（再掲）により、モデルパラメータベクトルθ（ｋ）を算出する。
θ（ｋ）＝θ（ｋ−１）＋ＫＰ（ｋ）ｉｄｅ（ｋ） （１５）
【０１３８】
式（１５）の同定誤差ｉｄｅ（ｋ）は、下記式（１７）（再掲）、（１８）（再掲）及び（１９ａ）により定義される。式（１９ａ）は、前記式（１９）のむだ時間ｄを「０」としたものである。ゲイン係数ベクトルＫＰ（ｋ）は、下記式（２０）（再掲）により定義され、式（２０）の正方行列Ｐ（ｋ）は、下記式（２１）（再掲）により算出される
【数１２】
本実施形態では、第１の実施形態と同様に下記要求Ｂ１〜Ｂ３を満たすことに加えて、さらに下記要求Ｂ４及びＢ５を満たすことが求められる。
Ｂ１）準静的動特性変化及びハードウエアの特性ばらつきに対する適応
「準静的動特性変化」とは、例えば電源電圧の変動やハードウエアの経年劣化といった変化速度の遅い特性変化である。
Ｂ２）動的な動特性変化への適応
具体的には、スロットル弁開度ＴＨの変化に対応する動特性変化への適応を意味する。
Ｂ３）モデルパラメータのドリフト防止
モデルパラメータに反映すべきでない制御対象の非線形特性などに起因する同定誤差の影響によって、モデルパラメータの絶対値が増大するような不具合を防止する。
Ｂ４）ＥＣＵの演算能力とマッチング
具体的には、演算量をより低減させることが求められる。
Ｂ５）モデルパラメータ（制御性能）の安定化
具体的は、同定されるモデルパラメータのばらつきを極力抑制することが求められる。
【０１４０】
先ず要求Ｂ４を満たすために、係数λ１及びλ２をそれぞれ１，０に設定することにより、固定ゲインアルゴリズムを採用する。これによって、正方行列Ｐ（ｋ）は一定となるため、式（２１）の演算を省略することができ、演算量を大幅に低減できる。
【０１４１】
すなわち固定ゲインアルゴリズムを採用すると、式（２０）は、下記式（２０ａ）のように簡略化される。式（２０ａ）においてＰは、定数を対角要素とする正方行列である。
【数１４】
【０１４２】
そこで本実施形態では、このようなモデルパラメータのドリフトを防止するために、モデルパラメータベクトルθ（ｋ）を上記式（１５ｂ）に代えて、下記式（１５ｃ）により、算出するようにした。
ＤＥＬＴＡ＝［ＤＥＬＴＡ１，ＤＥＬＴＡ２，ＤＥＬＴＡ３，ＤＥＬＴＡ４］
【０１４３】
忘却係数ＤＥＬＴＡｉは、０から１の間の値に設定され（０＜ＤＥＬＴＡｉ＜１）、過去の同定誤差の影響を徐々に減少させる機能を有する。ただし、モデルパラメータｂ１の演算に係る係数ＤＥＬＴＡ３またはモデルパラメータｃ１の演算にかかる係数ＤＥＬＴＡ４の何れか一方は、「１」として、実質的に忘却係数が乗算されないようにする。このように、忘却係数ベクトルＤＥＬＴＡの要素の一部を「１」とすることにより、目標値ＤＴＨＲと、スロットル弁開度偏差量ＤＴＨとの定常偏差が発生することを防止することができる。なお、係数ＤＥＬＴＡ３及びＤＥＬＴＡ４をともに「１」とすると、モデルパラメータのドリフト防止効果が不十分となるので、何れか一方のみを「１」とすることが望ましい。
【０１４４】
式（１５ｃ）を漸化式形式に書き直すと、下記式（１５ｄ）（１５ｅ）が得られる。前記式（１５）に代えて下記式（１５ｄ）及び（１５ｅ）を用いてモデルパラメータベクトルθ（ｋ）を算出する手法を、以下δ修正法といい、式（１５ｅ）で定義されるｄθ（ｋ）を「更新ベクトル」という。
θ（ｋ）＝θ（０）＋ｄθ（ｋ） （１５ｄ）
ｄθ（ｋ）＝ＤＥＬＴＡ×ｄθ（ｋ−１）＋ＫＰ（ｋ）ｉｄｅ（ｋ） （１５ｅ）
【０１４５】
δ修正法を用いたアルゴリズムによれば、上記要求Ｂ３を満たすドリフト防止効果とともに、上記要求Ｂ５を満たすモデルパラメータの安定化効果も得られる。すなわち、初期値ベクトルθ（０）が常に保存され、更新ベクトルｄθ（ｋ）も忘却係数ベクトルＤＥＬＴＡの働きにより、その要素のとりうる値が制限されるので、各モデルパラメータを初期値近傍に安定させることができる。
【０１４６】
さらに実際の制御対象の入出力データに基づいた同定により更新ベクトルｄθ（ｋ）を調整しつつモデルパラメータを算出するので、実際の制御対象に適合したモデルパラメータを算出でき、上記要求Ｂ１も満たされる。
次に要求Ｂ２を満たすべく、本実施形態では上記式（１５ｄ）の初期値ベクトルθ（０）に代えて、基準モデルパラメータベクトルθｂａｓｅを用いる下記式（１５ｆ）により、モデルパラメータベクトルθ（ｋ）を算出することとした。
θ（ｋ）＝θｂａｓｅ＋ｄθ（ｋ） （１５ｆ）
【０１４７】
基準モデルパラメータベクトルθｂａｓｅは、モデルパラメータスケジューラ２５によりスロットル弁開度偏差量ＤＴＨ（＝ＴＨ−ＴＨＤＥＦ）に応じて設定されるので、スロットル弁開度ＴＨの変化に対応する動特性の変化に適応させることができ、上記要求Ｂ２を満たすことができる。
【０１４８】
以上のように本実施形態では、固定ゲインアルゴリズムを採用することにより、ＥＣＵの演算量の低減を図り（要求Ｂ４）、δ修正法を用いたアルゴリズムを採用することにより、準静的動特性変化及びハードウエアの特性ばらつきに対する適応（要求Ｂ１）、モデルパラメータ（制御性能）の安定化（要求Ｂ５）、及びモデルパラメータのドリフト防止（要求Ｂ３）を実現し、モデルパラメータスケジューラ２５を採用することにより、スロットル弁開度ＴＨの変化に対応した動特性変化への適応（要求Ｂ２）を実現している。
【０１４９】
なお、式（１５ｆ）により算出されるモデルパラメータベクトルθ（ｋ）の各要素ａ１’，ａ２’，ｂ１’及びｃ１’についてリミット処理を施し、修正モデルパラメータベクトルθＬ（ｋ）（θＬ（ｋ）Ｔ＝［ａ１，ａ２，ｂ１，ｃ１］）を算出する点は、第１の実施形態と同様である。
【０１５０】
また、モデルパラメータｃ１’を統計処理し、その変動の中心値をデフォルト開度ずれｔｈｄｅｆａｄｐとして算出し、下記式（４１）（４２）（再掲）によりスロットル弁開度偏差量ＤＴＨ及び目標値ＤＴＨＲを算出する点も、第１の実施形態と同様である。
ＤＴＨ（ｋ）＝ＴＨ（ｋ）−ＴＨＤＥＦ＋ｔｈｄｅｆａｄｐ （４１）
ＤＴＨＲ（ｋ）＝ＴＨＲ（ｋ）−ＴＨＤＥＦ＋ｔｈｄｅｆａｄｐ （４２）
【０１５１】
次に上述した適応スライディングモードコントローラ２１ａ、モデルパラメータ同定器２２ａ及びモデルパラメータスケジューラ２５の機能を実現するための、ＥＣＵ７のＣＰＵにおける演算処理を説明する。
【０１５２】
図３２は、スロットル弁開度制御の全体フローチャートである。この処理は、図９に示すスロットル弁開度制御処理のステップＳ１３（状態予測器の演算）を削除し、ステップＳ１２，Ｓ１４及びＳ１６を、それぞれステップＳ１２ａ，１４ａ及び１６ａに変更したものである。
【０１５３】
ステップＳ１２ａでは、図３３に示すモデルパラメータ同定器の演算、すなわち前記式（１５ｆ）によるモデルパラメータベクトルθ（ｋ）の算出処理を実行し、さらにリミット処理を実行して修正モデルパラメータベクトルθＬ（ｋ）を算出する。
【０１５４】
ステップＳ１４ａでは、修正モデルパラメータベクトルθＬ（ｋ）を用いて、図３６に示す制御入力Ｕｓｌ（ｋ）の演算処理を実行する。すなわち、前記式（９ｂ）（１０ｂ）（１１ｂ）により、等価制御入力Ｕｅｑ、到達則入力Ｕｒｃｈ（ｋ）及び適応則入力Ｕａｄｐ（ｋ）を算出し、それらの入力の総和として、制御入力Ｕｓｌ（ｋ）（＝デューティ比ＤＵＴ（ｋ））を算出する。
【０１５５】
ステップＳ１６ａでは、図４１に示すスライディングモードコントローラの安定判別処理を実行する。すなわち、予測切換関数値σｐｒｅに代えて、切換関数値σを用いてスライディングモードコントローラの安定性判別を行い、安定判別フラグＦＳＭＣＳＴＡＢの設定を行う。この安定判別フラグＦＳＭＣＳＴＡＢは、「１」に設定されたときの処理は、第１の実施形態と同様である。
【０１５６】
図３３は、モデルパラメータ同定器２２ａの演算処理のフローチャートである。この処理は、図１１に示すモデルパラメータ同定器の演算処理のステップＳ３１〜Ｓ３４をそれぞれ、ステップＳ３１ａ〜Ｓ３４ａに変更し、さらにステップＳ３３ｂ及びＳ３３ｃを追加したものである。
【０１５７】
ステップＳ３１ａでは、式（２０ａ）によりゲイン係数ベクトルＫＰ（ｋ）を算出し、次いで式（１８）及び（１９ａ）により推定スロットル弁開度偏差量ＤＴＨＨＡＴ（ｋ）を算出する（ステップＳ３２ａ）。
ステップＳ３３ａでは、図３５に示すｉｄｅ（ｋ）の演算処理を実行し、同定誤差ｉｄｅ（ｋ）を算出する。ステップＳ３３ｂでは、式（１５ｅ）により更新ベクトルｄθ（ｋ）を算出し、次いでスロットル弁開度偏差量ＤＴＨに応じて図３４に示すθｂａｓｅテーブルを検索し、基準モデルパラメータベクトルθｂａｓｅを算出する（ステップＳ３３ｃ）。θｂａｓｅテーブルには、基準モデルパラメータａ１ｂａｓｅ，ａ２ｂａｓｅ及びｂ１ｂａｓｅが設定されている。スロットル弁開度偏差量ＤＴＨが「０」近傍の値をとる（スロットル弁開度ＴＨが、デフォルト開度ＴＨＤＥＦ近傍である）とき、基準モデルパラメータａ１ｂａｓｅ及びｂ１ｂａｓｅは減少し、基準モデルパラメータａ２ｂａｓｅは増加するように設定されている。また、基準モデルパラメータｃ１ｂａｓｅは、「０」に設定される。
【０１５８】
ステップＳ３４ａでは、式（１５ｆ）によりモデルパラメータベクトルθ（ｋ）を算出し、次いで第１の実施形態と同様に、モデルパラメータベクトルθ（ｋ）の安定化処理を実行する（ステップＳ３５）。すなわち各モデルパラメータのリミット処理を行って修正モデルパラメータベクトルθＬ（ｋ）を算出する。
【０１５９】
図３５は、図３３のステップＳ３３ａで実行されるｉｄｅ（ｋ）演算処理のフローチャートである。この処理は、図１２のｉｄｅ（ｋ）演算処理のステップＳ５６（不感帯処理）を削除し、ステップＳ５１をステップＳ５１ａに変更したものである。すなわち本実施形態では、δ修正法により、モデルパラメータのドリフトが防止されるので、不感帯処理は実行しない。
【０１６０】
また、ステップＳ５１ａでは、推定スロットル弁開度偏差量ＤＴＨＨＡＴ（ｋ）を式（１８）及び（１９ａ）により算出し、この推定スロットル弁開度偏差量ＤＴＨＨＡＴ（ｋ）を用いて同定誤差ｉｄｅ（ｋ）を算出する。
本実施形態では、制御対象モデルのむだ時間ｄを「０」としているので、ステップＳ５２の所定値ＸＣＮＴＩＤＳＴは、例えば「２」に設定される。
【０１６１】
図３６は、図３２のステップＳ１４ａで実行される、スロットル弁駆動装置１０への制御入力Ｕｓｌ（＝ＤＵＴ）を算出する処理のフローチャートである。この処理は、図２２に示すＵｓｌ演算処理のステップＳ２０１〜Ｓ２０５を、それぞれステップＳ２０１ａ〜Ｓ２０５ａに変更したものである。
【０１６２】
ステップＳ２０１ａでは、図３７に示す切換関数値σの演算処理を実行し、ステップＳ２０２ａでは、図３８に示す切換関数値σの積算値の演算処理を実行する。ステップＳ２０３ａでは、前記式（９ｂ）により、等価制御入力Ｕｅｑを算出する。ステップＳ２０４ａでは、図３９に示す到達則入力Ｕｒｃｈの演算処理を実行し、ステップＳ２０５ａでは、図４０に示す適応則入力Ｕａｄｐの演算処理を実行する。
【０１６３】
図３７は、図３６のステップＳ２０１ａで実行される切換関数値σの演算処理のフローチャートである。この処理は、図２３に示す予測切換関数値σｐｒｅの演算処理のステップＳ２２２〜Ｓ２２６を、それぞれステップＳ２２２ａ〜２２６ａに変更したものである。
【０１６４】
ステップＳ２２２ａでは、前記式（５）により、切換関数値σ（ｋ）を算出する。続くステップＳ２２３ａ〜Ｓ２２６ａは、図２３のステップＳ２２３〜Ｓ２２６の「σｐｒｅ」を「σ」に置き換えたものであり、切換関数値σ（ｋ）に対して、図２３の処理と同様のリミット処理を行う。
【０１６５】
図３８は、図３６のステップＳ２０２ａで実行される、切換関数値σの積算値ＳＵＭＳＩＧＭＡａを算出する処理のフローチャートである。この処理は、図２６に示すσｐｒｅの積算値演算処理のステップＳ２４１〜Ｓ２４５を、それぞれステップＳ２４１ａ〜Ｓ２４５ａに変更したものである。積算値ＳＵＭＳＩＧＭＡａは、後述する図４０の処理で適応則入力Ｕａｄｐの算出に使用される（前記式（１１ｂ）参照）。
【０１６６】
ステップＳ２４１ａでは、下記式（４７ａ）により、積算値ＳＵＭＳＩＧＭＡａを算出する。
ＳＵＭＳＩＧＭＡａ（ｋ）＝ＳＵＭＳＩＧＭＡａ（ｋ−１）＋σ×ΔＴ （４７ａ）
続くステップＳ２４２ａ〜Ｓ２４５ａでは、算出した積算値ＳＵＭＳＩＧＭＡａに対して、図２６の処理と同様のリミット処理を行う。
【０１６７】
図３９は、図３６のステップＳ２０４ａで実行される到達則入力Ｕｒｃｈの演算処理のフローチャートである。この処理は、図２７に示す到達則入力Ｕｒｃｈ演算処理のステップＳ２６３及びＳ２６５をそれぞれステップＳ２６３ａ及びＳ２６５ａに変更したものである。
【０１６８】
すなわち、本実施形態では、予測切換関数値σｐｒｅではなく、切換関数値σを用いて、適応スライディングモードコントローラ２１ａが安定してるときの到達則入力Ｕｒｃｈの算出（ステップＳ２６３ａ）、及び適応スライディングモードコントローラ２１ａが不安定であるときの到達則入力Ｕｒｃｈの算出（ステップＳ２６５ａ）を実行する。
【０１６９】
図４０は、図３６のステップＳ２０５ａで実行される適応則入力Ｕａｄｐの演算処理のフローチャートである。この処理は、図２８に示す適応則入力Ｕａｄｐ演算処理のステップＳ２７３及びＳ２７５をそれぞれステップＳ２７３ａ及びＳ２７５ａに変更したものである。
【０１７０】
すなわち、本実施形態では、切換関数値σの積算値ＳＵＭＳＩＧＵＭＡａを用いて、適応スライディングモードコントローラ２１ａが安定してるときの適応則入力Ｕａｄｐの算出（ステップＳ２７３ａ）、及び適応スライディングモードコントローラ２１ａが不安定であるときの適応則入力Ｕａｄｐの算出（ステップＳ２７５ａ）を実行する。
【０１７１】
図４１は、図３２のステップＳ１６ａで実行されるスライディングモードコントローラの安定判別処理のフローチャートである。この処理は、図２９のステップＳ２８１及びＳ２８１をそれぞれステップＳ２８１ａ及びＳ２８２ａに変更したものである。
【０１７２】
ステップＳ２８１ａでは下記式（５２ａ）により、切換関数変化量Ｄσを算出し、ステップＳ２８２ａでは、下記式（５３ａ）により、安定性判別パラメータＳＧＭＳＴＡＢを算出する。すなわち、予測切換関数値σｐｒｅではなく、切換関数値σに基づいて安定性判別を行う。
Ｄσ＝σ（ｋ）−σ（ｋ−１） （５２ａ）
ＳＧＭＳＴＡＢ＝Ｄσ×σ（ｋ） （５３ａ）
【０１７３】
本実施形態では、スロットル弁駆動装置１０及びＥＣＵ７の一部（モータ６に駆動電流を供給する出力回路）がプラントに相当し、図３６の処理が制御手段に相当し、図３３の処理が同定手段に相当し、図３５の処理が同定誤差算出手段に相当し、図３３のステップＳ３３ｂが更新ベクトル算出手段及び更新ベクトル修正手段に相当し、図３３のステップＳ３４ａがモデルパラメータベクトル算出手段に相当し、図３３のステップＳ３５が制限手段に相当する。
【０１７４】
（第３の実施形態）
図４２は、本発明の第３の実施形態にかかる制御系の構成を示すブロック図である。
この制御系は、制御対象であるプラント１０１と、プラントの出力である混合液のｐＨ（ペーハー）を検出するｐＨセンサ１０２と、ｐＨセンサ出力Ｖ１ＯＵＴから第１基準値Ｖ１ＢＡＳＥを減算する減算器１０３と、制御目標値Ｖ１ＴＡＲＧＥＴを生成する目標値生成部１０４と、第１操作量Ｕ１を決定する操作量決定部１０５と、第１操作量Ｕ１と第２基準値Ｖ２ＢＡＳＥとを加算し、第２操作量Ｕ２を出力する加算器１０６とからなる。
【０１７５】
減算器１０３、目標値生成部１０４、操作量決定部１０５及び加算器１０６は、具体的にはＣＰＵ、メモリ、入出力回路などからなる電子コントロールユニットにより構成される。
プラント１０１は、第２操作量Ｕ２に応じてアルカリ液の流量を制御する流量制御弁１１１と、流量制御弁１１１を介して供給されるアルカリ液と、酸性液とを攪拌する攪拌器１１２とからなる。プラント１０１は、アルカリ液と酸性液とを攪拌することにより、所望のｐＨ値の混合液を出力するものである。
【０１７６】
操作量決定部１０５は、プラント１０１をモデル化した制御対象モデルのモデルパラメータベクトルを同定する同定器１２１と、適応スライディングモードコントローラ１２２と、予測器１２３とからなる。同定器１２１、適応スライディングモードコントローラ１２２、及び予測器１２３は、それぞれ第１の実施形態におけるモデルパラメータ同定器２２、適応スライディングモードコントローラ２１、及び状態予測器２３に対応し、これらと同様の機能を有する。
【０１７７】
以下本実施形態における構成要素及びパラメータと、第１の実施形態における構成要素及びパラメータとの対応関係を説明する。
ｐＨセンサ１０２は、スロットル弁開度センサ８に相当し、ｐＨセンサ１０２の出力Ｖ１ＯＵＴは、スロットル弁開度ＴＨに相当する。第１目標値Ｖ１ＢＡＳＥは、デフォルト開度ＴＨＤＥＦに相当するものであり、本実施形態では、例えば中性に対応するｐＨ値とする。したがって、偏差量ＤＶ１がスロットル弁開度偏差量ＤＴＨに相当する。また目標値生成部１０４が、目標開度設定部２４に相当し、制御目標値Ｖ１ＴＡＲＧＥＴがスロットル弁開度偏差量の目標値ＤＴＨＲに相当する。なお、第１の実施形態においては、減算器１０３の機能は、モデルパラメータ同定器２２及び状態予測器２３に含まれている。
【０１７８】
第２基準値Ｖ２ＢＡＳＥは、適応スライディングモードコントローラ１２２の出力である第１操作量Ｕ１の中心値をバイアスするために加算されるものである。第１の実施形態では加算器１０６に相当する構成要素はなく、したがって第２基準値Ｖ２ＢＡＳＥは実質的に「０」とされている（すなわち、Ｕ１＝Ｕ２＝Ｕｓｌである）。本実施形態では、第２基準値Ｖ２ＢＡＳＥは、例えば流量制御弁１１１の開度が５０％となるような値に設定される。
【０１７９】
流量制御弁１１１は、デューティ比ＤＵＴのパルス信号でオンオフ制御されるスイッチング素子（ＥＣＵ７の出力回路に含まれ、図示及び説明を省略している）に相当し、アルカリ液は電源電圧に相当する。また流量制御弁１１１の出力流量Ｖ２は、モータ６の駆動電流に相当し、攪拌器１１２は、モータ６及びスロットル弁３の弁体に相当し、酸性液はスロットル弁３の弁体に加わる吸気管負圧や、リターンスプリング４及び弾性部材５の付勢力に相当する。攪拌器１１２から出力されるの混合液のｐＨ値Ｖ１が、実際のスロットル弁開度に相当する。
【０１８０】
以上のような対応関係があるので、プラント１０１を第１の実施形態と同様にモデル化し、同様の制御手法を適用することができる。すなわち、同定器１２１は、第１操作量Ｕ１及び偏差量ＤＶ１に基づいて、第１の実施形態と同様の演算処理により、修正モデルパラメータパラメータベクトルθＬを算出し、予測器１２３は、第１操作量Ｕ１、偏差量ＤＶ１及び修正モデルパラメータベクトルθＬに基づいて、第１の実施形態と同様の演算処理により、予測偏差量ＰＲＥＤＶ１を算出し、適応スライディングモードコントローラ１２２は、予測偏差量ＰＲＥＤＶ１及び修正モデルパラメータベクトルθＬに基づいて、第１の実施形態と同様の演算処理により、予測偏差量ＰＲＥＤＶ１を制御目標値Ｖ１ＴＡＲＧＥＴに一致させるように、第１操作量Ｕ１を算出する。したがって、制御目標値Ｖ１ＴＡＲＧＥＴとして、所望の相対ｐＨ値（第１基準値Ｖ１ＢＡＳＥとの偏差量）を設定することにより、プラントの出力Ｖ１を所望のｐＨ値に一致させることができる。
本実施形態では、同定器１２１が同定手段に相当し、同定誤差算出手段、同定誤差修正手段、モデルパラメータベクトル算出手段及び制限手段を含む。
【０１８１】
（第３の実施形態の変形例）
図４３は、図４２に示す構成の変形例を示す図である。この変形例では、図４２のプラント１０１ではなく、プラント１０１ａが制御対象とされる。プラント１０１ａは、プラント１０１に、流量制御弁１１１の出力流量Ｖ２を検出する流量センサ１１３と、流量センサ出力Ｖ２ＯＵＴが第２操作量Ｕ２に対応する流量値と一致するように、流量制御弁１１１を制御するフィードバック制御器１１４とを追加して構成されている。
【０１８２】
このようにローカルフィードバックループを含むプラントに対しても、第３の実施形態と同様のモデル化及び同様の制御手法の適用が可能である。
なお第１の実施形態では、モータの駆動回路は公知のものであるため、詳細な説明を行っていないが、オンオフ制御されるスイッチング素子の出力電流を検出する電流センサを設け、検出電流値ＩＤが、操作量Ｕｓｌに対応する電流値ＩＲと一致するようにフィードバック制御を行うようにしてもよく、本変形例は第１の実施形態においてそのような回路構成を採用した場合に相当する。
【０１８３】
（第４の実施形態）
図４４は、本発明の第４の実施形態にかかる制御系の構成を示すブロック図である。この制御系は、図４２の操作量決定部１０５を操作量決定部１０５ａに代えたものであり、第２の実施形態として示した制御系に対応する。なお、以下に説明する以外は第３の実施形態と同一である。
【０１８４】
操作量決定部１０５ａは、同定器１２１ａと、適応スライディングモードコントローラ１２２ａと、パラメータスケジューラ１２４とからなる。
同定器１２１ａ、適応スライディングモードコントローラ１２２ａ及びパラメータスケジューラ１２４は、それぞれ第２の実施形態におけるモデルパラメータ同定器２２ａ、適応スライディングモードコントローラ２１ａ及びモデルパラメータスケジューラ２５に対応し、これらと同様の機能を有する。
【０１８５】
すなわち、パラメータスケジューラ１２４は、偏差量ＤＶ１に基づいて第２の実施形態と同様の演算処理により、基準モデルパラメータベクトルθｂａｓｅを算出し、同定器１２１ａは、第１操作量Ｕ１、偏差量ＤＶ１及び基準モデルパラメータベクトルθｂａｓｅに基づいて第２の実施形態と同様の演算処理により、修正モデルパラメータベクトルθＬを算出し、適応スライディングモードコントローラ１２２ａは、偏差量ＤＶ１及び修正モデルパラメータベクトルθＬに基づいて、第２の実施形態と同様の演算処理により偏差量ＤＶ１を制御目標値Ｖ１ＴＡＲＧＥＴに一致させるように、第１操作量Ｕ１を算出する。したがって、制御目標値Ｖ１ＴＡＲＧＥＴとして、所望の相対ｐＨ値（第１基準値Ｖ１ＢＡＳＥとの偏差量）を設定することにより、プラントの出力Ｖ１を所望のｐＨ値に一致させることができる。
【０１８６】
本実施形態では、同定器１２１ａが同定手段に相当し、同定誤差算出手段、更新ベクトル算出手段、更新ベクトル修正手段、モデルパラメータベクトル算出手段及び制限手段を含む。
【０１８７】
（第４の実施形態の変形例）
図４５は、図４４に示す構成の変形例を示す図である。この変形例では、図４４のプラント１０１ではなく、プラント１０１ａが制御対象とされる。プラント１０１ａは、図４３のプラント１０１ａと同一である。
このようにローカルフィードバックループを含むプラントに対しても、第４の実施形態と同様のモデル化及び同様の制御手法の適用が可能である。
【０１８８】
（その他の実施形態）
モデルパラメータの同定誤差ｉｄｅ（ｋ）の算出手法として、δ修正法に代えて、以下に述べるε修正法を採用してもよい。すなわち、前記式（１５ｃ）に代えて、下記式（１５ｇ）により、モデルパラメータベクトルθ（ｋ）を算出するようにしてもよい。
ＥＰＳ＝［ＥＰＳ１，ＥＰＳ２，ＥＰＳ３，ＥＰＳ４］
【０１８９】
忘却係数ＥＰＳ１，ＥＰＳ２及びＥＰＳ４は、前記忘却係数ＤＥＬＴＡｉと同様に、０から１の間の値に設定され（０＜ＥＰＳｉ＜１）、過去の同定誤差の影響を徐々に減少させる機能を有する。
【０１９０】
ただしε修正法の場合、モデルパラメータｂ１の演算に係る係数ＥＰＳ３は、必ず「１」にする必要がある。これは、以下の理由による。ε修正法の場合、同定誤差ｉｄｅ（ｋ）が小さくなると、モデルパラメータはすべてゼロ近傍の値となる。ところが、モデルパラメータｂ１は、式（９ｂ）、（１０ｂ）、（１１ｂ）の分母に適用されるため、モデルパラメータｂ１が「０」に近づくと、制御対象への入力Ｕｓｌが発散するからである。
【０１９１】
式（１５ｇ）は、初期値ベクトルθ（０）にも忘却係数ベクトルＥＰＳが乗算されている点で、式（１５ｃ）と異なる。
式（１５ｇ）を漸化式形式に書き直すと、下記式（１５ｈ）が得られる。前記式（１５）に代えて下記式（１５ｈ）を用いてモデルパラメータベクトルθ（ｋ）を算出する手法を、ε修正法と呼ぶ。
θ（ｋ）＝ＥＰＳ×θ（ｋ−１）＋ＫＰ（ｋ）ｉｄｅ（ｋ） （１５ｈ）
【０１９２】
ε修正法によっても、過去の同定誤差ｅｉｄの影響が低減されるので、モデルパラメータのドリフトを防止することができる。
また第２の実施形態では、δ修正法によりモデルパラメータのドリフトを防止するようにしたが、第１の実施形態と同様に、不感帯処理（図１４）により修正同定誤差ｉｄｅｎｌ（ｋ）を算出し、これを用いてモデルパラメータベクトルθ（ｋ）を算出するようにしてもよい。
【０１９３】
また第１の実施形態において、不感帯処理に代えてδ修正法またはε修正法を採用してもよい。さらに、第１の実施形態においてδ修正法を採用する場合には、第２の実施形態と同様にモデルパラメータスケジューラを導入し、モデルパラメータスケジューラにより算出される基準モデルパラメータベクトルθｂａｓｅに更新ベクトルを加算する形式でモデルパラメータベクトルθを算出することが望ましい。
【０１９４】
【発明の効果】
以上詳述したように請求項１に記載の発明によれば、モデルパラメータベクトルの同定誤差に応じて更新ベクトルが算出され、同定誤差の過去値の影響を減少させるように更新ベクトルが修正され、モデルパラメータの基準ベクトルに修正された更新ベクトルを加算することにより、モデルパラメータベクトルが算出され、さらにモデルパラメータベクトルの要素の値が所定リミット範囲内に制限される。したがって、モデルパラメータのドリフトを防止しつつ、制御の安定性をより一層向上させることができる。またモデルパラメータの基準ベクトルは、前記プラントの動特性変化を示すパラメータに応じて算出されるので、プラントの動特性変化に応じた適切な基準ベクトルが得られる。その結果、特にプラントに非線形要素が含まれる場合でもモデルパラメータを迅速に収束させることができる。
【０１９５】
請求項２に記載の発明によれば、固定ゲインアルゴリズムを用いて更新ベクトルが算出されるので、演算量を低減することができる。
【０１９６】
請求項３に記載の発明によれば、更新ベクトルの少なくとも１つの要素の過去値に０より大きく１より小さい所定値を乗算することにより、更新ベクトルの修正が行われるので、記同定誤差の過去値の影響が減少し、モデルパラメータベクトルのドリフトを防止することができる。
【０１９７】
請求項４に記載の発明によれば、更新ベクトルの、プラントの入力に関わる要素またはプラントの入出力に関わらない要素については、０より大きく１より小さい所定値が乗算されないので、更新ベクトルの修正による定常偏差の発生を防止することができる。
【０２００】
請求項５に記載の発明によれば、モデルパラメータの同定誤差が不感帯内にあるときは、同定誤差が減少方向に修正され、該修正された同定誤差を用いてモデルパラメータベクトルが算出され、さらにモデルパラメータベクトルの要素の値が所定リミット範囲内に制限される。したがって、モデルパラメータのドリフトを防止しつつ、制御の安定性をより一層向上させることができる。また前記不感帯は、プラントの出力または制御目標値に応じて設定されるので、モデルパラメータの値に反映させるべき同定誤差を、残留同定誤差として無視してしまうことを防止することができる。
【０２０１】
請求項６に記載の発明によれば、モデルパラメータベクトルの複数の要素が所定の関係を満たすように、複数の要素の値が制限されるので、そのモデルパラメータベクトルを用いる制御の安定性を向上させることができる。
【０２０２】
請求項７に記載の発明によれば、同定誤差が前記不感帯内にあるときは、同定誤差が０とされるので、モデルパラメータの値に反映されるべきでない同定誤差の影響を無くし、モデルパラメータのドリフト防止効果を高めることができる。
【０２０３】
請求項８に記載の発明によれば、同定手段により同定されたモデルパラメータを用いて、スロットル弁開度を目標開度に一致させる制御が行われるので、スロットル弁開度の目標開度への制御性を向上させるとともに、スロットル弁開度制御の安定性を向上させることができる。
【図面の簡単な説明】
【図１】本発明の一実施形態にかかる内燃機関のスロットル弁駆動装置と、その制御装置を示す図である。
【図２】図１に示すスロットル弁駆動装置の周波数特性を示す図である。
【図３】図１の電子制御ユニット（ＥＣＵ）により実現される機能を示す機能ブロック図である。
【図４】スライディングモードコントローラの制御特性と、切換関数設定パラメータ（ＶＰＯＬＥ）の値との関係を示す図である。
【図５】スライディングモードコントローラの制御ゲイン（Ｆ，Ｇ）の設定範囲を示す図である。
【図６】モデルパラメータのドリフトを説明するための図である。
【図７】同定誤差を修正する関数を示す図である。
【図８】スロットル弁のデフォルト開度ずれがモデルパラメータ（ｃ１’）に反映されることを説明するための図である。
【図９】スロットル弁開度制御処理のフローチャートである。
【図１０】図９の処理において状態変数の設定を行う処理のフローチャートである。
【図１１】図９の処理においてモデルパラメータ同定器の演算を実行する処理のフローチャートである。
【図１２】図１１の処理において同定誤差（ｉｄｅ）の演算を実行する処理のフローチャートである。
【図１３】同定誤差（ｉｄｅ）のローパスフィルタ処理を説明するための図である。
【図１４】図１２の処理における不感帯処理のフローチャートである。
【図１５】図１４の処理で使用されるテーブルを示す図である。
【図１６】図１１の処理におけるモデルパラメータベクトル（θ）の安定化処理のフローチャートである。
【図１７】図１６の処理におけるモデルパラメータ（ａ１’，ａ２’）のリミット処理のフローチャートである。
【図１８】図１６の処理によるモデルパラメータの値の変化を説明するための図である。
【図１９】図１６の処理におけるモデルパラメータ（ｂ１’）のリミット処理のフローチャートである。
【図２０】図１６の処理におけるモデルパラメータ（ｃ１’）のリミット処理のフローチャートである。
【図２１】図９の処理において状態予測器の演算を実行する処理のフローチャートである。
【図２２】図９の処理において制御入力（Ｕｓｌ）の演算を実行する処理のフローチャートである。
【図２３】図２２の処理において予測切換関数値（σｐｒｅ）の演算を実行する処理のフローチャートである。
【図２４】図２３の処理において切換関数設定パラメータ（ＶＰＯＬＥ）の演算を実行する処理のフローチャートである。
【図２５】図２４の処理で使用するマップを示す図である。
【図２６】図２２の処理において予測切換関数値（σｐｒｅ）の積算値の演算を実行する処理のフローチャートである。
【図２７】図２２の処理において到達則入力（Ｕｒｃｈ）の演算を実行する処理のフローチャートである。
【図２８】図２２の処理において適応則入力（Ｕａｄｐ）の演算を実行する処理のフローチャートである。
【図２９】図９の処理においてスライディングモードコントローラの安定判別を実行する処理のフローチャートである。
【図３０】図９の処理においてデフォルト開度ずれ（ｔｈｄｅｆａｄｐ）の演算を実行する処理のフローチャートである。
【図３１】図１の電子制御ユニット（ＥＣＵ）により実現される機能を示す機能ブロック図である（第２の実施形態）。
【図３２】スロットル弁開度制御処理のフローチャートである（第２の実施形態）。
【図３３】図３２の処理においてモデルパラメータ同定器の演算を実行する処理のフローチャートである。
【図３４】図３３の処理で使用するテーブルを示す図である。
【図３５】図３３の処理において同定誤差（ｉｄｅ）の演算を実行する処理のフローチャートである。
【図３６】図３２の処理において制御入力（Ｕｓｌ）の演算を実行する処理のフローチャートである。
【図３７】図３６の処理において切換関数値（σ）の演算を実行する処理のフローチャートである。
【図３８】図３６の処理において切換関数値（σ）の積算値の演算を実行する処理のフローチャートである。
【図３９】図３６の処理において到達則入力（Ｕｒｃｈ）の演算を実行する処理のフローチャートである。
【図４０】図３６の処理において適応則入力（Ｕａｄｐ）の演算を実行する処理のフローチャートである。
【図４１】図３２の処理においてスライディングモードコントローラの安定判別を実行する処理のフローチャートである。
【図４２】本発明の第３の実施形態に係る制御系の構成を示すブロック図である。
【図４３】図４２に示す構成の変形例を示すブロック図である。
【図４４】本発明の第４の実施形態に係る制御系の構成を示すブロック図である。
【図４５】図４４に示す構成の変形例を示すブロック図である。
【符号の説明】
１ 内燃機関
３ スロットル弁
７ 電子制御ユニット
１０ スロットル弁駆動装置
２１，２１ａ 適応スライディングモードコントローラ（制御手段）
２２，２２ａ モデルパラメータ同定器（同定手段、同定誤差算出手段、更新ベクトル算出手段、更新ベクトル修正手段、同定誤差修正手段、モデルパラメータベクトル算出手段、制限手段）
２４ 目標開度設定部
２５ モデルパラメータスケジューラ[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention relates to a control apparatus for a plant, and more particularly to a control apparatus that performs control using the model parameters while real-time identifying model parameters indicating characteristics of a control target model obtained by modeling a plant to be controlled. .
[0002]
[Prior art]
The plant to be controlled is modeled, the model parameters of the model to be controlled are calculated by a parameter adjustment mechanism, and the air-fuel ratio of the air-fuel mixture supplied to the internal combustion engine is fed back to the target air-fuel ratio by an adaptive controller using the model parameters. 2. Description of the Related Art An air-fuel ratio control device to be controlled is conventionally known (for example, JP-A-11-73206).
[0003]
In this control device, the model parameter is calculated by adding an update component calculated according to the identification error of the model parameter to the initial value of the model parameter. The updated component of the model parameter is corrected so as to reduce the influence of the past value of the identification error, thereby preventing the model parameter from drifting under the influence of disturbance and destabilizing the control.
[0004]
[Problems to be solved by the invention]
In the device described in the above publication, drift of model parameters due to the influence of disturbance is prevented. However, since the model parameters calculated by the parameter adjustment mechanism are directly used by the adaptive controller, there is room for improvement in maintaining the stability of the adaptive controller.
[0005]
The present invention has been made by paying attention to this point, and it is possible to further improve the stability of control using the model parameters while identifying the model parameters of the control target model obtained by modeling the plant to be controlled. It is an object of the present invention to provide a control device for a plant that can be used.
[0006]
[Means for Solving the Problems]
In order to achieve the above object, the invention according to
[0007]
According to this configuration, the update vector is calculated according to the identification error of the model parameter vector, the update vector is corrected so as to reduce the influence of the past value of the identification error, and the update vector corrected to the reference vector of the model parameter Is added, the model parameter vector is calculated, and the values of the elements of the model parameter vector are further restricted within a predetermined limit range. Therefore, the stability of the control can be further improved while preventing the drift of the model parameters.In addition, since the reference vector of the model parameter is calculated according to the parameter indicating the dynamic characteristic change of the plant, an appropriate reference vector according to the dynamic characteristic change of the plant can be obtained. As a result, it is possible to quickly converge the model parameters even when a non-linear element is included in the plant.
[0008]
According to a second aspect of the present invention, in the plant control device according to the first aspect, the update vector calculation means calculates the update vector using a fixed gain algorithm.
According to this configuration, since the update vector is calculated using the fixed gain algorithm, the amount of calculation can be reduced.
[0009]
According to a third aspect of the present invention, in the control device for a plant according to the first or second aspect, the update vector correction unit determines that a past value of at least one element of the update vector is greater than 0 and less than 1 The correction of the update vector is performed by multiplying (DELTAIi, EPSi).
[0010]
According to this configuration, the update vector is corrected by multiplying the past value of at least one element of the update vector by a predetermined value larger than 0 and smaller than 1, so that the influence of the past value of the identification error is reduced. , The drift of the model parameter vector can be prevented.
[0011]
According to a fourth aspect of the present invention, in the plant control device according to the third aspect, the update vector correcting unit includes an element (an element relating to an operation of b1) of the update vector relating to an input of the plant. Elements not related to the input / output of the plant (elements related to the calculation of c1) are not multiplied by the predetermined value (DELTAi, EPSi).
[0012]
According to this configuration, the elements related to the input of the plant or the elements not related to the input / output of the plant of the update vector are not multiplied by the predetermined value larger than 0 and smaller than 1. Can be prevented.
[0015]
Claim5The invention described in the above, an identification means for identifying a model parameter vector of a control target model obtained by modeling a plant based on the input and output of the plant, using a model parameter vector identified by the identification meansSo that the output of the plant matches the control target value.In a plant control device provided with control means for controlling the plant, the identification means calculates an identification error (ide) of the model parameter vector,Dead zone setting means for setting a dead zone according to the output or control target value of the plant,The identification error isThe dead zoneWhen (−EIDNRLMT ≦ ide ≦ EIIDNLMT), the model parameter vector is obtained by using an identification error correction unit that corrects the identification error in a decreasing direction, and the identification error (idenl) corrected by the identification error correction unit. And a limiter that limits the values of the elements of the model parameter vector calculated by the model parameter vector calculator within a predetermined limit range.
[0016]
According to this configuration, the identification error of the model parameter is reduced.Dead zoneWhen it is within the range, the identification error is corrected in the decreasing direction, the model parameter vector is calculated using the corrected identification error, and the values of the elements of the model parameter vector are limited to within a predetermined limit range. Therefore, the stability of the control can be further improved while preventing the drift of the model parameters.Further, since the dead zone is set according to the output of the plant or the control target value, it is possible to prevent the identification error to be reflected in the value of the model parameter from being ignored as the residual identification error.
[0017]
Claim6The invention described in claim 15Any of1 itemIn the plant control device according to the above, the restricting means restricts values of the plurality of elements (a1, a2) of the model parameter vector such that the plurality of elements (a1, a2) satisfy a predetermined relationship (FIG. 18). It is characterized by.
[0018]
According to this configuration, the values of the plurality of elements are limited so that the plurality of elements of the model parameter vector satisfy the predetermined relationship, so that the stability of control using the model parameter vector can be improved.
Claim7The invention described in claim53. The control device for a plant according to
[0019]
According to this configuration, the identification error isDead zoneWhen the value is within the range, the identification error is set to 0, so that the effect of the identification error that should not be reflected in the value of the model parameter can be eliminated, and the drift prevention effect of the model parameter can be enhanced.
Claim8The invention described in claim 17Any of1 itemIn the control device of the plant described in the above, the plant includes a throttle valve driving device having a throttle valve of the internal combustion engine, and driving means for driving the throttle valve, the control means, the opening degree of the throttle valve A parameter for determining a control input to the throttle valve driving device is calculated so as to match the target opening.
[0020]
According to this configuration, since the control for matching the throttle valve opening to the target opening is performed using the model parameters identified by the identification unit, the controllability of the throttle valve opening to the target opening is improved. At the same time, the stability of the throttle valve opening control can be improved.
[0021]
BEST MODE FOR CARRYING OUT THE INVENTION
Hereinafter, embodiments of the present invention will be described with reference to the drawings.
(1st Embodiment)
FIG. 1 is a diagram showing a configuration of a throttle valve control device according to a first embodiment of the present invention. A throttle valve 3 is provided in an
[0022]
The
[0023]
The ECU 7 is connected to an accelerator sensor 9 for detecting an accelerator pedal depression amount ACC for detecting a required output of a driver of a vehicle equipped with the
The ECU 7 includes an input circuit to which the detection signals of the throttle
[0024]
In the present embodiment, a throttle
[0025]
When the response frequency characteristics of the throttle
DTH (k) = TH (k) -THDEF (2)
Here, TH is the detected throttle valve opening, and THDEF is the default opening.
Further, a1, a2, b1, and c1 in the equation (1) are model parameters that determine the characteristics of the control target model, and d is a dead time.
[0026]
The model defined by the above equation (1) is a discrete-time DARX model (delayed autoregressive model with exogeneous input: an autoregressive model having an external input) adopted to facilitate the application of adaptive control.
In the equation (1), in addition to the model parameters a1 and a2 related to the output deviation DTH and the model parameter b1 related to the input duty ratio DUT, a model parameter c1 not related to input / output is set. The model parameter c1 is a parameter indicating a deviation of the default opening THDEF and a disturbance applied to the throttle valve driving device. That is, the model parameters a1, a2,b1At the same time, by identifying the model parameter c1, the default opening deviation and disturbance can be identified.
[0027]
FIG. 3 is a functional block diagram of a throttle valve control device realized by the ECU 7. This control device includes an adaptive sliding
[0028]
The adaptive sliding
By using the adaptive sliding
[0029]
The
[0030]
By using the
[0031]
The
[0032]
Next, the operation principle of the adaptive sliding
First, the target value DTHR (k) is defined as a deviation amount between the target opening THR (k) and the default opening THDEF by the following equation (3).
DTHR (k) = THR (k) -THDEF (3)
Here, if the deviation e (k) between the throttle valve opening deviation amount DTH and the target value DTHR is defined by the following equation (4), the switching function value σ (k) of the adaptive sliding mode controller is expressed by the following equation (5). ).
[0033]
On a phase plane defined by the vertical axis as the deviation e (k) and the horizontal axis as the previous deviation e (k−1), the deviation e (k) satisfying σ (k) = 0 and the previous deviation e (k) Since the combination with -1) is a straight line, this straight line is generally called a switching straight line. The sliding mode control is a control focusing on the behavior of the deviation e (k) on the switching line, so that the switching function value σ (k) becomes 0, that is, the deviation e (k) and the previous deviation e (k). -1) is controlled so that the combination is placed on the switching line on the phase plane, and robust against disturbances and modeling errors (differences between the characteristics of the actual plant and the characteristics of the modeled control target model). In this case, the throttle valve opening deviation amount DTH follows the target value DTHR.
[0034]
Also, by changing the value of the switching function setting parameter VPOLE in the equation (5), as shown in FIG. 4, the attenuation characteristic of the deviation e (k), that is, the follow-up of the throttle valve opening deviation amount DTH to the target value DTHR. Properties can be changed. Specifically, if VPOLE = -1, the characteristic does not follow at all, and the following speed can be increased as the absolute value of the switching function setting parameter VPOLE decreases.
[0035]
The throttle valve control device is required to satisfy the following requirements A1 and A2.
A1) Avoiding collision with the throttle fully closed stopper when moving the throttle valve 3 to the fully closed position.
A2) Non-linear characteristics in the vicinity of the default opening THDEF (change in elastic characteristics due to the balance between the urging force of the
Therefore, the convergence speed of the deviation e (k) needs to be reduced near the fully closed position of the throttle valve, and the convergence speed needs to be increased near the default opening THDEF.
[0036]
According to the sliding mode control, the convergence speed can be easily changed by changing the switching function setting parameter VPOLE. Therefore, in this embodiment, the change amount DDTHR (= DTHR (k) of the throttle valve opening TH and the target value DTHR is changed. ) -DTHR (k-1)), the switching function setting parameter VPOLE is set. As a result, the above requirements A1 and A2 can be satisfied.
[0037]
As described above, in the sliding mode control, the combination of the deviation e (k) and the previous deviation e (k−1) (hereinafter referred to as “deviation state quantity”) is constrained on the switching straight line, so that the deviation e (k) is reduced. It converges to 0 at the specified convergence speed and robustly against disturbances and modeling errors. Therefore, in the sliding mode control, it is important how to put the deviation state quantity on the switching straight line and restrict it there.
[0038]
From such a viewpoint, the input (output of the controller) DUT (k) (also referred to as Usl (k)) to the control target is, as shown in the following equation (6), the equivalent control input Ueq (k) It is configured as the sum of the law input Urch (k) and the adaptive law input Uadp (k).
The equivalent control input Ueq (k) is an input for constraining the deviation state quantity on the switching straight line, and the reaching law input Urch (k) is an input for placing the deviation state quantity on the switching straight line. The law input Uadp (k) is an input for suppressing the effects of modeling errors and disturbances and placing the deviation state quantity on the switching straight line. Hereinafter, a method of calculating each input Ueq (k), Urch (k) and Uadp (k) will be described.
[0040]
Since the equivalent control input Ueq (k) is an input for restricting the deviation state quantity on the switching straight line, the condition to be satisfied is given by the following equation (7).
σ (k) = σ (k + 1) (7)
When the duty ratio DUT (k) satisfying the expression (7) is obtained using the expressions (1) and (4) and (5), the following expression (9) is obtained, which is equivalent to the equivalent control input Ueq (k). Become. Further, the reaching law input Urch (k) and the adaptive law input Uadp (k) are defined by the following equations (10) and (11), respectively.
(Equation 1)
Here, F and G are a reaching law control gain and an adaptive law control gain, respectively, and are set as described below. ΔT is a control cycle.
The calculation of equation (9) requires the throttle valve opening deviation amount DTH (k + d) after the elapse of the dead time d and the corresponding target value DTHR (k + d + 1). Therefore, the predicted deviation amount PREDTH (k) calculated by the
[0042]
Next, the reaching law control gain F and the adaptive law control gain G are determined by the reaching law input Urch and the adaptive law input Uadp so that the deviation state quantity is stably placed on the switching line.
Specifically, a condition for the switching function value σ (k) to be stable with respect to the disturbance V (k) is obtained by assuming the disturbance V (k), thereby obtaining the setting conditions of the gains F and G. As a result, it was obtained that the combination of the gains F and G satisfies the following equations (12) to (14), in other words, that the combination is within the area indicated by hatching in FIG.
[0043]
F> 0 (12)
G> 0 (13)
F <2- (ΔT / 2) G (14)
As described above, the equivalent control input Ueq (k), the reaching law input Urch (k), and the adaptive law input Uadp (k) are calculated by the equations (9) to (11), and the duty is calculated as the sum of those inputs. The ratio DUT (k) can be calculated.
[0044]
The
θ (k) = θ (k−1) + KP (k) ide (k) (15)
θ (k)T= [A1 ', a2', b1 ', c1'] (16)
[0045]
Here, a1 ', a2', b1 ', and C1' are model parameters before performing a limit process described later. Ide (k) is an identification error defined by the following equations (17), (18) and (19). DTHHAT (k) is an estimated value of the throttle valve opening deviation amount DTH (k) calculated using the latest model parameter vector θ (k-1) (hereinafter, referred to as “estimated throttle valve opening deviation amount”). It is. KP (k) is a gain coefficient vector defined by the following equation (20). P (k) in Expression (20) is a fourth-order square matrix calculated by Expression (21) below.
(Equation 2)
By setting the coefficients λ1 and λ2 in Expression (21), the identification algorithm according to Expressions (15) to (21) is any of the following four identification algorithms.
λ1 = 1, λ2 = 0 Fixed gain algorithm
λ1 = 1, λ2 = 1 Least squares algorithm
λ1 = 1, λ2 = λ Gradual gain algorithm (λ is a predetermined value other than 0, 1)
λ1 = λ, λ2 = 1 Weighted least squares algorithm (λ is a predetermined value other than 0 and 1)
[0047]
On the other hand, in the present embodiment, it is required to satisfy the following requirements B1), B2), and B3).
B1) Adaptation to quasi-static dynamic characteristic change and hardware characteristic variation
The “quasi-static dynamic characteristic change” means a characteristic change with a slow change speed, such as a change in power supply voltage or aging of hardware.
B2) Adaptation to dynamic changes in dynamic characteristics
Specifically, it means adaptation to a dynamic characteristic change corresponding to a change in the throttle valve opening TH.
B3) Drift prevention of model parameters
A problem that an absolute value of a model parameter increases due to an influence of an identification error caused by a non-linear characteristic of a control target that should not be reflected in the model parameter is prevented.
[0048]
First, in order to satisfy the requirements of B1) and B2), the weighted least squares algorithm is adopted by setting the coefficients λ1 and λ2 to predetermined values λ and “0”, respectively.
Next, the drift of the model parameters will be described. As shown in FIG. 6, after the model parameters have converged to some extent, there is a residual identification error caused by non-linear characteristics such as the friction characteristics of the throttle valve, or when a disturbance having a non-zero average value is constantly applied. , Residual identification errors accumulate, causing model parameter drift.
[0049]
Since such a residual identification error should not be reflected in the value of the model parameter, a dead zone process is performed using a dead zone function Fnl as shown in FIG. Specifically, the corrected identification error idenl (k) is calculated by the following equation (23), and the model parameter vector θ (k) is calculated using the corrected identification error idenl (k). That is, the following equation (15a) is used instead of the above equation (15). Thereby, the above requirement B3) can be satisfied.
idenl (k) = Fnl (ide (k)) (23)
θ (k) = θ (k−1) + KP (k) idenl (k) (15a)
[0050]
The dead zone function Fnl is not limited to the one shown in FIG. 7A, but is, for example, a discontinuous dead zone function as shown in FIG. 7B or an incomplete dead zone function as shown in FIG. May be used. However, when an incomplete dead zone function is used, drift cannot be completely prevented.
[0051]
The amplitude of the residual identification error changes in accordance with the amount of change in the throttle valve opening TH. Therefore, in the present embodiment, the dead zone width parameter EIDNRLMT defining the dead zone width shown in FIG. 7 is changed according to the root mean square value DDTHRSQA of the amount of change in the target throttle valve opening THR calculated by the following equation (24). (Specifically, as the root mean square value DDTHRSQA increases, the dead band width parameter EIDNRLMT is set to increase). As a result, it is possible to prevent the identification error to be reflected in the value of the model parameter from being ignored as the residual identification error. DDTHR in Expression (24) is a change amount of the target throttle valve opening THR, and is calculated by Expression (25) below.
(Equation 4)
Here, since the throttle valve opening deviation DTH is controlled by the adaptive sliding
[0053]
In order to further enhance the robustness of the control system, it is effective to further stabilize the adaptive sliding
[0054]
Next, a method of calculating the predicted deviation amount PREDTH by the
First, matrices A and B and vectors X (k) and U (k) are defined by the following equations (26) to (29).
(Equation 5)
X (k + 1) = AX (k) + BU (k−d) (30)
[0055]
When X (k + d) is obtained from the equation (30), the following equation (31) is obtained.
(Equation 6)
(Equation 7)
DTHHAT (k + d), which is the element of the first row of the prediction vector XHAT (k + d), is the prediction deviation amount PREDTH (k) and is given by the following equation (35).
[0057]
The predicted deviation PREDTH (k) calculated by the equation (35) is applied to the equation (9), and the target values DTHR (k + d + 1), DTHR (k + d), and DTHR (k + d-1) are respectively converted to DTHR (k + d-1). By substituting k), DTHR (k-1) and DTHR (k-2), the following equation (9a) is obtained. The equivalent control input Ueq (k) is calculated by the equation (9a).
(Equation 8)
Further, using the predicted deviation amount PREDTH (k) calculated by the equation (35), the prediction switching function value σpre (k) is defined by the following equation (36), and the reaching law input Urch (k) and the adaptive law input Uadp (k) is calculated by the following equations (10a) and (11a), respectively.
Next, as described above, the model parameter c1 'is a parameter indicating the deviation of the default opening THDEF and a disturbance. Therefore, as shown in FIG. 8, the default opening deviation can be considered to be substantially constant within a relatively short period, though it fluctuates due to disturbance. Therefore, in the present embodiment, the model parameter c1 'is statistically processed, the center value of the variation is calculated as the default opening deviation thdefadp, and used for calculating the throttle valve opening deviation DTH and the target value DTHR.
[0060]
The least squares method is generally known as a statistical processing method. Generally, the statistical processing by the least squares method stores all data within a certain period, that is, all identified model parameters c1 ′ in a memory. This is performed by performing a batch operation at a certain point in time. However, this batch operation method requires an enormous capacity of memory to store all data, and further requires an inverse matrix operation, which leads to an increase in the amount of operation.
[0061]
Therefore, in the present embodiment, the adaptive least-squares algorithm of adaptive control represented by the above equations (15) to (21) is applied to statistical processing, and the least square center value of the model parameter c1 is set to the default opening deviation. It is calculated as thdefadp.
[0062]
Specifically, θ (k) and θ (k) in the above equations (15) to (21)TTo thdefadp, and ζ (k) and ζ (k)TIs replaced by “1”, ide (k) is replaced by ec1 (k), KP (k) is replaced by KPTH (k), P (k) is replaced by PTH (k), λ1 and λ2 Are respectively substituted with λ1 ′ and λ2 ′ to obtain the following formulas (37) to (40).
(Equation 10)
One of the above four algorithms can be selected by setting the coefficients λ1 ′ and λ2 ′. In the equation (39), the coefficient λ1 ′ is set to a predetermined value other than 0 or 1, and the coefficient λ2 ′ is set. Is set to 1 to adopt the weighted least squares method.
[0064]
In the calculations of the above equations (37) to (40), only the values to be stored are thdefadp (k + 1) and PTH (k + 1), and the inverse matrix calculation is unnecessary. Therefore, by employing the sequential least squares algorithm, it is possible to perform statistical processing of the model parameter c1 by the least squares method while overcoming the drawbacks of the general least squares method.
[0065]
The default opening deviation thdefadp obtained as a result of the statistical processing is applied to the equations (2) and (3), and the throttle valve is calculated by the following equations (41) and (42) instead of the equations (2) and (3). An opening deviation amount DTH (k) and a target value DTHR (k) are calculated.
DTH (k) = TH (k) -THDEF + thdefadp (41)
DTHR (k) = THR (k) -THDEF + thdefadp (42)
[0066]
By using equations (41) and (42), even if the default opening THDEF deviates from the design value due to hardware characteristic variations or aging, accurate control is performed by compensating for the deviation. be able to.
[0067]
Next, calculation processing in the CPU of the ECU 7 for realizing the functions of the above-described adaptive sliding
[0068]
FIG. 9 is an overall flowchart of the throttle valve opening control, and this process is executed by the CPU of the ECU 7 every predetermined time (for example, 2 msec).
In step S11, a state variable setting process shown in FIG. 10 is executed. That is, the calculations of the equations (41) and (42) are executed to calculate the throttle valve opening deviation amount DTH (k) and the target value DTHR (k) (FIG. 10, steps S21 and S22). Note that (k) indicating the current value may be omitted.
[0069]
In step S12, the calculation of the model parameter identifier shown in FIG. 11, that is, the calculation processing of the model parameter vector θ (k) according to the equation (15a) is performed, and the limit processing is further performed to execute the correction model parameter vector θL (k ) Is calculated.
[0070]
In the following step S13, the calculation of the state predictor shown in FIG. 21 is executed to calculate the predicted deviation amount PREDTH (k).
Next, using the corrected model parameter vector θL (k) calculated in step S12, the arithmetic processing of the control input Usl (k) shown in FIG. 22 is executed (step S14). That is, the equivalent control input Ueq, the reaching law input Urch (k), and the adaptive law input Uadp (k) are calculated, and the control input Usl (k) (= duty ratio DUT (k)) is calculated as the sum of those inputs. I do.
[0071]
In a succeeding step S16, a stability determination process of the sliding mode controller shown in FIG. 29 is executed. That is, the stability determination is performed based on the differential value of the Lyapunov function, and the stability determination flag FSMCSTAB is set. When the stability determination flag FSMCSTAB is set to “1”, it indicates that the adaptive sliding
In step S17, a thdefadp calculation process shown in FIG. 30 is executed to calculate a default opening deviation thdefadp.
[0072]
FIG. 11 is a flowchart of the calculation processing of the
[0073]
In the following step S34, the model parameter vector θ (k) is calculated by the equation (15a), and then the model parameter vector θ (k) is stabilized (step S35). That is, a limit process of each model parameter is performed to calculate a corrected model parameter vector θL (k).
[0074]
FIG. 12 is a flowchart of the idlen (k) calculation process executed in step S33 of FIG.
In step S51, the identification error ide (k) is calculated by the equation (17). Next, the value of the counter CNTIDST incremented in step S53 is set to a predetermined value XCNTIDST set according to the dead time d of the control target (for example, set to “3” corresponding to the dead time d = 2). It is determined whether or not the value is greater than (step S52). Since the initial value of the counter CNTIDST is "0", the process first proceeds to step S53, where the counter CNTIDST is incremented by "1" and the identification error ide (k) is set to "0" (step S54). Proceed to S55. Immediately after the identification of the model parameter vector θ (k) is started, a correct identification error cannot be obtained by the calculation according to the equation (17). Therefore, in steps S52 to S54, the identification error is calculated without using the calculation result according to the equation (17). ide (k) is set to “0”.
[0075]
When the answer to step S52 is affirmative (YES), the process immediately proceeds to step S55.
In step S55, low-pass filter processing of the identification error ide (k) is performed. Specifically, when identifying a model parameter of a control target having a low-pass characteristic, the identification weight for the identification error ide (k) of the least-squares identification algorithm is represented by a frequency characteristic as shown by a solid line L1 in FIG. However, this is a characteristic in which high-frequency components are attenuated by low-pass filter processing as shown by a broken line L2. This is for the following reason.
[0076]
The frequency characteristics of the actual control target and the control target model obtained by modeling the control target are as shown by solid lines L3 and L4 in FIG. 13B, respectively. That is, when the model parameter is identified by the
[0077]
Therefore, the frequency characteristic of the weight of the identification algorithm is set to the characteristic indicated by the broken line L2 in FIG. 13A by the low-pass filter processing, so that the frequency characteristic of the control target model is indicated by the broken line L5 in FIG. The characteristics are as shown in the drawings, and the actual frequency characteristics are made to coincide with each other, or the gain of the control target model is corrected to be slightly higher than the actual gain. Thus, overcorrection by the
[0078]
In the low-pass filter processing, the past values ide (ki) of the identification error (for example, ten past values corresponding to i = 1 to 10) are stored in a ring buffer, and the past values are multiplied by a weight coefficient. And add them.
Further, since the identification error ide (k) is calculated using the above equations (17), (18) and (19), the throttle valve opening deviation amount DTH (k) and the estimated throttle valve opening deviation The same low-pass filter processing is performed on the amount DTHHAT (k), or the throttle valve opening deviation amounts DTH (k-1) and DTH (k-2), and the duty ratio DUT (k-d-1) The same effect can be obtained by performing the same low-pass filter processing as described above.
[0079]
Returning to FIG. 12, in the following step S56, the dead zone process shown in FIG. 14 is executed. In step S61 of FIG. 14, the root mean square value DDTHRSQA of the amount of change in the target throttle valve opening THR is calculated by, for example, n = 5 in the equation (24), and then EIDNRLMT shown in FIG. 15 according to the root mean square value DDTHRSQA. The table is searched to calculate a dead zone width parameter EIDNLLMT (step S62).
[0080]
In step S63, it is determined whether or not the identification error ide (k) is larger than the dead band width parameter EIDNRLMT. If ide (k)> EIDNRLMT, the modified identification error idenl (k) is calculated by the following equation (43). (Step S67).
idlen (k) = ide (k) -EIDNRLMT (43)
[0081]
If the answer to step S63 is negative (NO), it is further determined whether or not the identification error ide (k) is smaller than a value obtained by adding a minus sign to the dead zone width parameter EIDNRLMT (step S64), and ide (k) If <-EIIDNLLMT, the corrected identification error idenl (k) is calculated by the following equation (44) (step S65).
idlen (k) = ide (k) + EIDNRLMT (44)
If the identification error ide (k) is within the range of ± EIDNRLMT, the modified identification error idenl (k) is set to “0” (step S66).
[0082]
FIG. 16 is a flowchart of the stabilization processing of θ (k) executed in step S35 of FIG.
In step S71, initialization is performed by setting flags FA1STAB, FA2STAB, FB1LMT, and FC1LMT used in this process to “0”, respectively. Then, in step S72, limit processing of a1 'and a2' shown in FIG. 17 is executed, in step S73, limit processing of b1 'shown in FIG. 19 is executed, and in step S74, c1' shown in FIG. Execute limit processing.
[0083]
FIG. 17 is a flowchart of the limit processing of a1 'and a2' executed in step S72 of FIG. FIG. 18 is a diagram for explaining the process of FIG. 17 and is referred to together with FIG.
In FIG. 18, the combination of the model parameters a1 ′ and a2 ′ requiring the limit processing is indicated by “x”, and the range of the stable combination of the model parameters a1 ′ and a2 ′ is indicated by a hatched area (hereinafter referred to as “×”). Stability region). The process of FIG. 17 is a process of moving the combination of the model parameters a1 'and a2' outside the stable region to the inside of the stable region (the position indicated by "○").
[0084]
In step S81, it is determined whether or not the model parameter a2 'is equal to or greater than a predetermined a2 lower limit value XIDA2L. The predetermined a2 lower limit value XIDA2L is set to a negative value larger than “−1”. Even if the predetermined a2 lower limit value XIDA2L is set to “−1”, stable modified model parameters a1 and a2 can be obtained, but the nth power of the matrix A defined by the equation (26) becomes unstable ( This means that a1 'and a2' do not diverge but oscillate), so that the value is set to a value larger than "-1".
[0085]
If a2 '<XIDA2L in step S81, the modified model parameter a2 is set to this lower limit value XIDA2L, and the a2 stabilization flag FA2STAB is set to "1". When the a2 stabilization flag FA2STAB is set to “1”, it indicates that the modified model parameter a2 has been set to the lower limit value XIDA2L. In FIG. 18, the correction of the model parameters by the limit processing P1 in steps S81 and S82 is indicated by an arrow line with “P1” (a line with an arrow).
[0086]
If the answer to step S81 is affirmative (YES), that is, if a2 '≧ XIDA2L, the corrected model parameter a2 is set to the model parameter a2' (step S83).
In steps S84 and S85, it is determined whether or not the model parameter a1 'is within a range defined by a predetermined a1 lower limit value XIDA1L and a predetermined a1 upper limit value XIDA1H. The predetermined al lower limit value XIDA1L is set to a value equal to or greater than -2 and smaller than 0, and the predetermined al upper limit value XIDA1H is set to, for example, 2.
[0087]
If both the answers of steps S84 and S85 are affirmative (YES), that is, if XIDA1L ≦ a1 ′ ≦ XIDA1H, the modified model parameter a1 is set to the model parameter a1 ′ (step S88).
On the other hand, if a1 '<XIDA1L, the correction model parameter a1 is set to the lower limit value XIDA1L, and the a1 stabilization flag FA1STAB is set to "1" (steps S84 and S86). When a1 '> XIDA1H, the modified model parameter a1 is set to the upper limit value XIDA1H, and the a1 stabilization flag FA1STAB is set to "1" (steps S85 and S87). When the a1 stabilization flag FA1STAB is set to “1”, it indicates that the corrected model parameter a1 has been set to the lower limit value XIDA1L or the upper limit value XIDA1H. In FIG. 18, the correction of the model parameters by the limit processing P2 in steps S84 to S87 is indicated by an arrow with “P2”.
[0088]
In step S90, it is determined whether or not the sum of the absolute value of the corrected model parameter a1 and the corrected model parameter a2 is equal to or smaller than a predetermined stability determination value XA2STAB. The predetermined stability determination value XA2STAB is set to a value close to “1” and smaller than “1” (for example, 0.99).
[0089]
The straight lines L1 and L2 shown in FIG. 18 are straight lines satisfying the following equation (45).
a2 + | a1 | = XA2STAB (45)
Therefore, in step S90, it is determined whether or not the combination of the modified model parameters a1 and a2 is on or below the straight lines L1 and L2 shown in FIG. If the answer to step S90 is affirmative (YES), the combination of the modified model parameters a1 and a2 is within the stable region of FIG.
[0090]
On the other hand, when the answer to step S90 is negative (NO), the modified model parameter a1 is a value obtained by subtracting the predetermined a2 lower limit value XIDA2L from the predetermined stability determination value XA2STAB (XIDA2L <0, so that XA2STAB−XIDA2L> XA2STAB is satisfied). (Established) or not (step S91). If the modified model parameter a1 is equal to or less than (XA2STAB-XIDA2L), the modified model parameter a2 is set to (XA2STAB- | a1 |) and the a2 stabilization flag FA2STAB is set to "1" (step S92). .
[0091]
If the corrected model parameter a1 is larger than (XA2STAB-XIDA2L) in step S91, the corrected model parameter a1 is set to (XA2STAB-XIDA2L), the corrected model parameter a2 is set to the predetermined a2 lower limit value XIDA2L, and a1 is stabilized. The flag FA1STAB and the a2 stabilization flag FA2STAB are both set to "1" (step S93).
[0092]
In FIG. 18, the correction of the model parameters by the limit processing P3 of steps S91 and S92 is indicated by an arrow with “P3”, and the correction of the model parameters by the limit processing P4 of steps S91 and S93 is: This is indicated by an arrow with “P4”.
As described above, the limit processing is executed by the processing of FIG. 17 so that the model parameters a1 'and a2' fall within the stable region shown in FIG. 18, and the corrected model parameters a1 and a2 are calculated.
[0093]
FIG. 19 is a flowchart of the limit process of b1 'executed in step S73 of FIG.
In steps S101 and S102, it is determined whether or not the model parameter b1 'is within a range defined by a predetermined b1 lower limit XIDB1L and a predetermined b1 upper limit XIDB1H. The predetermined b1 lower limit value XIDB1L is set to a positive predetermined value (for example, 0.1), and the predetermined b1 upper limit value XIDB1H is set to, for example, "1".
[0094]
When the answers in steps S101 and S102 are both affirmative (YES), that is, when XIDB1L ≦ b1 ′ ≦ XIDB1H, the modified model parameter b1 is set to the model parameter b1 ′ (step S105).
On the other hand, if b1 '<XIDB1L, the correction model parameter b1 is set to the lower limit value XIDB1L, and the b1 limit flag FB1LMT is set to "1" (steps S101 and S104). If b1 '> XIDB1H, the correction model parameter b1 is set to the upper limit value XIDB1H and the b1 limit flag FB1LMT is set to "1" (steps S102 and S103). When the b1 limit flag FB1LMT is set to “1”, it indicates that the modified model parameter b1 has been set to the lower limit value XIDB1L or the upper limit value XIDB1H.
[0095]
FIG. 20 is a flowchart of the limit process of the model parameter c1 'executed in step S74 of FIG.
In steps S111 and S112, it is determined whether or not the model parameter c1 'is within a range defined by a predetermined c1 lower limit value XIDC1L and a predetermined c1 upper limit value XIDC1H. The predetermined c1 lower limit value XIDC1L is set to, for example, -60, and the predetermined c1 upper limit value XIDC1H is set to, for example, 60.
[0096]
If the answers of steps S111 and S112 are both affirmative (YES), that is, if XIDC1L ≦ c1 ′ ≦ XIDC1H, the corrected model parameter c1 is set to the model parameter c1 ′ (step S115).
On the other hand, when c1 '<XIDC1L, the modified model parameter c1 is set to the lower limit value XIDC1L, and the c1 limit flag FC1LMT is set to "1" (steps S111, S114). If c1 '> XIDC1H, the correction model parameter c1 is set to the upper limit value XIDC1H, and the c1 limit flag FC1LMT is set to "1" (steps S112 and S113). When the c1 limit flag FC1LMT is set to “1”, it indicates that the corrected model parameter c1 has been set to the lower limit value XIDC1L or the upper limit value XIDC1H.
[0097]
FIG. 21 is a flowchart of the calculation processing of the state predictor executed in step S13 of FIG.
In step S121, a matrix operation is performed to calculate the matrix elements α1, α2, β1 to β2, and γ1 to γd of the equation (35).
In step S122, the predicted deviation amount PREDTH (k) is calculated by equation (35).
[0098]
FIG. 22 is a flowchart of the process of calculating the control input Usl (= DUT) to the throttle
In step S201, a calculation process of the prediction switching function value σpre shown in FIG. 23 is executed, and in step S202, a calculation process of an integrated value of the prediction switching function value σpre shown in FIG. 26 is executed. In step S203, the equivalent control input Ueq is calculated by the equation (9). In step S204, the arithmetic processing of the reaching law input Urch shown in FIG. 27 is executed, and in step S205, the arithmetic processing of the adaptive law input Uadp shown in FIG. 28 is executed.
[0099]
In step S206, it is determined whether or not a stability determination flag FSMCSTAB set in a process of FIG. 29 described later is “1”. When the stability determination flag FSMCSTAB is set to “1”, it indicates that the adaptive sliding
[0100]
If FSMCSTAB = 0 in step S206 and the adaptive sliding
[0101]
On the other hand, when FSMCSTAB = 1 and the adaptive sliding
[0102]
In the following steps S209 and S210, it is determined whether or not the calculated control input Usl is within the range of the predetermined upper and lower limit values XUSLH and XUSLL. The process ends. On the other hand, when the control input Usl is equal to or less than the predetermined lower limit XUSLL, the control input Usl is set to the predetermined lower limit XUSLL (steps S209 and S212). When the control input Usl is equal to or greater than the predetermined upper limit XUSLH, the control input Usl is set. Usl is set to a predetermined upper limit value XUSLH (steps S210, S211).
[0103]
FIG. 23 is a flowchart of the calculation process of the predicted switching function value σpre performed in step S201 of FIG.
In step S221, the calculation processing of the switching function setting parameter VPOLE shown in FIG. 24 is executed, and then the calculation of the predicted switching function value σpre (k) is executed by the equation (36) (step S222).
[0104]
In the following steps S223 and S224, it is determined whether or not the calculated predicted switching function value σpre (k) is within the range of the predetermined upper and lower limit values XSGMH and XSGML, and the predicted switching function value σpre (k) is determined as the predetermined upper and lower limit value. If it is within the range, this process is immediately terminated. On the other hand, when the predicted switching function value σpre (k) is equal to or smaller than the predetermined lower limit value XSGML, the predicted switching function value σpre (k) is set to the predetermined lower limit value XSGML (steps S223 and S225), and the predicted switching function value σpre (k) is set. If k) is equal to or larger than the predetermined upper limit value XSGMH, the prediction switching function value σpre (k) is set to the predetermined upper limit value XSGMH (steps S224 and S226).
[0105]
FIG. 24 is a flowchart of the calculation process of the switching function setting parameter VPOLE executed in step S221 of FIG.
In step S231, it is determined whether or not the stability determination flag FSMCSTAB is "1". If FSMCSTAB = 1 and the adaptive sliding
[0106]
If FSMCSTAB = 0 and the adaptive sliding
DDTHR (k) = DTHR (k) -DTHR (k-1) (46)
In step S234, a VPOLE map is searched according to the throttle valve opening deviation amount DTH and the target value change amount DDTHR calculated in step S233, and a switching function setting parameter VPOLE is calculated. As shown in FIG. 25A, the VPOLE map increases when the throttle valve opening deviation amount DTH takes a value near 0 (when the throttle valve opening TH takes a value near the default opening THDEF), At a value other than 0, the value is set to be substantially constant with respect to the change of the throttle valve opening deviation amount DTH. The VPOLE map is set so that the VPOLE value increases as the target value change amount DDTHR increases, as indicated by the solid line in FIG. When a near value is taken, the target value change amount DDTHR is set to increase when a near-zero value is taken, as shown by the broken line in FIG.
[0107]
That is, when the target value DTHR of the throttle valve opening greatly changes in the decreasing direction, the switching function setting parameter VPOLE is set to a relatively small value. Thus, it is possible to prevent the throttle valve 3 from colliding with the throttle fully closed stopper. In the vicinity of the default opening THDEF, the switching function setting parameter VPOLE is set to a relatively large value, and controllability in the vicinity of the default opening THDEF can be improved.
[0108]
When the throttle valve opening TH is near the fully closed position or near the fully opened position as shown in FIG. 3C, the switching function setting parameter VPOLE may be set to be decreased. Accordingly, when the throttle valve opening TH is near the full-closed opening or near the fully-opened position, the follow-up speed to the target opening THR becomes slow, and the fully-closed stopper of the throttle valve 3 (even at the fully-opened position functions as a stopper). ) Can be more reliably prevented.
[0109]
In the following steps S235 and S236, it is determined whether or not the calculated switching function setting parameter VPOLE is within the range of predetermined upper and lower limit values XPOLEH and XPOLEL, and when the switching function setting parameter VPOLE is within the range of the predetermined upper and lower limit value. Ends this process immediately. On the other hand, when the switching function setting parameter VPOLE is equal to or smaller than the predetermined lower limit XPOLEL, the switching function setting parameter VPOLE is set to the predetermined lower limit XPOLEL (step S23).6, S238If the switching function setting parameter VPOLE is equal to or larger than the predetermined upper limit value XPOLEH, the switching function setting parameter VPOLE is set to the predetermined upper limit value XPOLEH (step S23).5, S237).
[0110]
FIG. 26 is a flowchart of the process of calculating the integrated value SUMSIGMA of the predicted switching function value σpre, which is performed in step S202 of FIG. The integrated value SUMSIGMA is used to calculate the adaptive law input Uadp in the processing of FIG. 28 described below (see the above equation (11a)).
[0111]
In step S241, an integrated value SUMSIGMA is calculated by the following equation (47). ΔT in the following equation is the execution cycle of the calculation.
SUMSIGMA (k) = SUMSIGMA (k-1) + σpre × ΔT (47)
In the following steps S242 and S243, it is determined whether or not the calculated integrated value SUMSIGMA is within the range of the predetermined upper and lower limit values XSUMSH and XSUMSL. If the integrated value SUMSIGMA is within the range of the predetermined upper and lower limit value, the process proceeds immediately. The process ends. On the other hand, when the integrated value SUMSIGMA is equal to or less than the predetermined lower limit value XSUMSL, the integrated value SUMSIGMA is set to the predetermined lower limit value XSUMSL (steps S242 and S244). SUMSIGMA is set to a predetermined upper limit value XSUMSH (steps S243 and S245).
[0112]
FIG. 27 is a flowchart of the reaching law input Urch calculation process executed in step S204 of FIG.
In step S261, it is determined whether or not the stability determination flag FSMCSTAB is “1”. When the stability determination flag FSMCSTAB is “0” and the adaptive sliding
Urch = −F × σpre / b1 (48)
[0113]
On the other hand, when the stability determination flag FSMCSTAB is "1" and the adaptive sliding
Urch = −F × σpre (49)
[0114]
In subsequent steps S266 and S267, it is determined whether or not the calculated reaching law input Urch is within the range of predetermined upper and lower limits XURCHH and XURCHL. When the reaching law input Urch is within the range of the predetermined upper and lower limits, This processing ends immediately. On the other hand, when the reaching law input Urch is equal to or smaller than the predetermined lower limit value XURCHL, the reaching law input Urch is set to the predetermined lower limit value XURCHL (steps S266 and S268). Then, the reaching law input Urch is set to the predetermined upper limit value XURCHH (steps S267 and S269).
[0115]
When the adaptive sliding
[0116]
FIG. 28 is a flowchart of the calculation process of the adaptive law input Uadp executed in step S205 of FIG.
In the step S271, it is determined whether or not the stability determination flag FSMCSTAB is “1”. When the stability determination flag FSMCSTAB is “0” and the adaptive sliding
Uadp = −G × SUMSIGMA / b1 (50)
[0117]
On the other hand, when the stability determination flag FSMCSTAB is "1" and the adaptive sliding
Uadp = −G × SUMSIGMA (51)
[0118]
When the adaptive sliding
[0119]
FIG. 29 is a flowchart of the stability determination process of the sliding mode controller executed in step S16 of FIG. In this process, stability determination based on the differential term of the Lyapunov function is performed, and a stability determination flag FSMCSTAB is set.
[0120]
In step S281, the switching function change amount Dσpre is calculated by the following equation (52), and then the stability determination parameter SGMSTAB is calculated by the following equation (53) (step S282).
Dσpre = σpre (k) −σpre (k−1) (52)
SGMSTAB = Dσpre × σpre (k) (53)
In step S283, it is determined whether or not the stability determination parameter SGMSTAB is equal to or smaller than a stability determination threshold value XSGSTAB. If SGSTAB> XSGMSTAB, it is determined that the
[0121]
In step S286, it is determined whether or not the value of the instability detection counter CNTSMCST is equal to or smaller than a predetermined count value XSSTAB. If CNTSMCST ≦ XSSTAB, the
[0122]
In the following step S289, the stability determination period counter CNTJUDST is decremented by "1", and then it is determined whether or not the value of the stability determination period counter CNTJUDST is "0" (step S290). The stability determination period counter CNTJUDST is initialized to a predetermined determination count value XCJUDST when the ignition switch is turned on. Therefore, initially, the answer to step S290 is negative (NO), and the process immediately proceeds to step S295.
[0123]
Thereafter, when the stability determination period counter CNTJUDST becomes “0”, the process proceeds from step S290 to step S291, and it is determined whether the first determination flag FSMCSTAB1 is “1”. When the first determination flag FSMCSTAB1 is "0", the second determination flag FSMCSTAB2 is set to "0" (step S293), and when the first determination flag FSMCSTAB1 is "1", the second determination is made. The flag FSMCSTAB2 is set to "1" (step S292).
[0124]
In a succeeding step S294, the value of the stability determination period counter CNTJUDST is set to the predetermined determination count value XCJUDST, and the value of the instability detection counter CNTSMCST is set to "0", and the process proceeds to step S295.
In step S295, the stability determination flag FSMCSTAB is set to the logical sum of the first determination flag FSMCSTAB1 and the second determination flag FSMCSTAB2. Regarding the second determination flag FSMCSTAB2, even if the answer to step S286 is affirmative (YES) and the first determination flag FSMCSTAB1 is set to “0”, the value of the stability determination period counter CNTJUDST becomes “0” until “0”. 1 "is maintained. Therefore, the stability determination flag FSMCSTAB is also maintained at “1” until the value of the stability determination period counter CNTJUDST becomes “0”.
[0125]
FIG. 30 is a flowchart of the process of calculating the default opening deviation thdefadp executed in step S17 of FIG.
In step S251, a gain coefficient KPTH (k) is calculated by the following equation (54).
KPTH (k) = PTH (k-1) / (1 + PTH (k-1)) (54)
[0126]
Here, PTH (k-1) is a gain parameter calculated in step S253 during the previous execution of this process.
In step S252, the model parameters calculated in the model parameter identifier calculation process shown in FIG.c1 'and the gain coefficient KPTH (k) calculated in step S251 are applied to the following equation (55) to calculate a default opening deviation thdefadp (k).
In step S253, a gain parameter PTH (k) is calculated by the following equation (56).
[0128]
By the processing in FIG. 30, the model parameter c1 'is statistically processed by the successive weighted least squares method, and the default opening deviation thdefadp is calculated.
In the present embodiment, a part of the throttle
[0129]
(Second embodiment)
In the above-described first embodiment, the control target model is defined using the equation (1) including the dead time d, and the
In order to further reduce the calculation load on the CPU, a fixed gain algorithm is adopted as an algorithm for identifying model parameters.
In addition, in order to further stabilize control, as a method of preventing model parameter drift,TeensOther techniques are used.
Hereinafter, the present embodiment will be described in detail focusing on differences from the first embodiment. Except for the points described below, the configuration is the same as that of the first embodiment.
[0131]
FIG. 31 is a functional block diagram of a throttle valve control device realized by the ECU 7. This control device includes an adaptive sliding
[0132]
Not the predicted deviation amount PREDTH but the detected throttle valve opening TH is input to the adaptive sliding
[0133]
By using the adaptive sliding
[0134]
The
[0135]
The
[0136]
In the present embodiment, since the control target model is defined by the equation (1a), the adaptive sliding
[Equation 11]
The equations (9b) to (11b) are obtained by setting the dead time d in the equations (9) to (11) to “0”.
The
θ (k) = θ (k−1) + KP (k) ide (k) (15)
[0138]
The identification error ide (k) in the equation (15) is defined by the following equations (17) (re-listed), (18) (re-listed), and (19a). Equation (19a) is obtained by setting the dead time d in Equation (19) to “0”. The gain coefficient vector KP (k) is defined by the following equation (20) (represented), and the square matrix P (k) of the equation (20) is calculated by the following equation (21) (represented).
(Equation 12)
In the present embodiment, in addition to satisfying the following requirements B1 to B3 as in the first embodiment, it is required to further satisfy the following requirements B4 and B5.
B1) Adaptation to quasi-static dynamic characteristic change and hardware characteristic variation
The “quasi-static dynamic characteristic change” is a characteristic change having a slow change speed, such as a change in a power supply voltage or aging of hardware.
B2) Adaptation to dynamic changes in dynamic characteristics
Specifically, it means adaptation to a dynamic characteristic change corresponding to a change in the throttle valve opening TH.
B3) Drift prevention of model parameters
A problem that an absolute value of a model parameter increases due to an influence of an identification error caused by a non-linear characteristic of a control target that should not be reflected in the model parameter is prevented.
B4) Matching with ECU's calculation ability
Specifically, it is required to further reduce the amount of calculation.
B5) Stabilization of model parameters (control performance)
Specifically, it is required to minimize variations in the model parameters to be identified.
[0140]
First, to satisfy the requirement B4, a fixed gain algorithm is adopted by setting the coefficients λ1 and λ2 to 1 and 0, respectively. As a result, the square matrix P (k) becomes constant, so that the calculation of Expression (21) can be omitted, and the amount of calculation can be significantly reduced.
[0141]
That is, when the fixed gain algorithm is adopted, Expression (20) is simplified as Expression (20a) below. In Expression (20a), P is a square matrix having constants as diagonal elements.
[Equation 14]
[0142]
Therefore, in the present embodiment, in order to prevent such a drift of the model parameter, the model parameter vector θ (k) is calculated by the following equation (15c) instead of the above equation (15b).
DELTA = [DELTA1, DELTA2, DELTA3, DELTA4]
[0143]
The forgetting factor DELTAi is set to a value between 0 and 1 (0 <DELTAi <1), and has a function of gradually reducing the influence of past identification errors. However, one of the coefficient DELTA3 related to the calculation of the model parameter b1 and the coefficient DELTA4 related to the calculation of the model parameter c1 is set to “1” so that the forgetting coefficient is not substantially multiplied. In this way, by setting some of the elements of the forgetting coefficient vector DELTA to "1", it is possible to prevent the occurrence of a steady-state deviation between the target value DTHR and the throttle valve opening deviation amount DTH. If both the coefficients DELTA3 and DELTA4 are set to "1", the effect of preventing the drift of the model parameters becomes insufficient. Therefore, it is preferable that only one of the coefficients is set to "1".
[0144]
When the equation (15c) is rewritten into a recurrence formula, the following equations (15d) and (15e) are obtained. A method of calculating the model parameter vector θ (k) using the following expressions (15d) and (15e) instead of the expression (15) is hereinafter referred to as a δ correction method, and dθ (d) defined by the expression (15e) k) is called an “update vector”.
θ (k) = θ (0) + dθ (k) (15d)
dθ (k) = DELTA × dθ (k−1) + KP (k) ide (k) (15e)
[0145]
According to the algorithm using the δ correction method, not only the drift prevention effect that satisfies the requirement B3 but also the model parameter stabilization effect that satisfies the requirement B5 can be obtained. In other words, the initial value vector θ (0) is always stored, and the update vector dθ (k) is also limited in the possible values of its elements by the function of the forgetting coefficient vector DELTA, so that each model parameter is stabilized near the initial value. Can be done.
[0146]
Further, since the model parameters are calculated while adjusting the update vector dθ (k) by the identification based on the input / output data of the actual control target, the model parameters suitable for the actual control target can be calculated, and the above requirement B1 is satisfied. .
Next, in order to satisfy the requirement B2, in the present embodiment, the model parameter vector θ (k) is calculated by the following equation (15f) using the reference model parameter vector θbase instead of the initial value vector θ (0) of the equation (15d). Was calculated.
θ (k) = θbase + dθ (k) (15f)
[0147]
Since the reference model parameter vector θbase is set by the
[0148]
As described above, in the present embodiment, the fixed gain algorithm is used to reduce the amount of computation of the ECU (request B4), and the algorithm using the δ correction method is used to change the quasi-static dynamic characteristics. By adopting the
[0149]
Note that a limit process is performed on each of the elements a1 ', a2', b1 ', and c1' of the model parameter vector θ (k) calculated by the equation (15f), and the corrected model parameter vector θL (k) (θL (k)T= [A1, a2, b1, c1]) is the same as in the first embodiment.
[0150]
Further, the model parameter c1 'is statistically processed, the center value of the variation is calculated as the default opening deviation thdefadp, and the throttle valve opening deviation amount DTH and the target value DTHR are calculated by the following equations (41) and (42). The calculation is also the same as in the first embodiment.
DTH (k) = TH (k) -THDEF + thdefadp (41)
DTHR (k) = THR (k) -THDEF + thdefadp (42)
[0151]
Next, a description will be given of arithmetic processing in the CPU of the ECU 7 for implementing the functions of the above-described adaptive sliding
[0152]
FIG. 32 is an overall flowchart of the throttle valve opening control. In this process, step S13 (operation of the state predictor) of the throttle valve opening control process shown in FIG. 9 is deleted, and steps S12, S14, and S16 are changed to steps S12a, 14a, and 16a, respectively.
[0153]
In step S12a, the calculation of the model parameter identifier shown in FIG. 33, that is, the calculation processing of the model parameter vector θ (k) by the equation (15f) is executed, and the limit processing is further executed to execute the correction model parameter vector θL (k ) Is calculated.
[0154]
In step S14a, the processing of calculating the control input Usl (k) shown in FIG. 36 is executed using the corrected model parameter vector θL (k). That is, the equivalent control input Ueq, the reaching law input Urch (k), and the adaptive law input Uadp (k) are calculated by the equations (9b), (10b), and (11b), and the control input Usl ( k) (= duty ratio DUT (k)) is calculated.
[0155]
In step S16a, a stability determination process of the sliding mode controller shown in FIG. 41 is executed. That is, the stability of the sliding mode controller is determined using the switching function value σ instead of the predicted switching function value σpre, and the stability determination flag FSMCSTAB is set. The processing when the stability determination flag FSMCSTAB is set to "1" is the same as in the first embodiment.
[0156]
FIG. 33 is a flowchart of the calculation process of the
[0157]
In step S31a, the gain coefficient vector KP (k) is calculated by equation (20a), and then the estimated throttle valve opening deviation amount DTHHAT (k) is calculated by equations (18) and (19a) (step S32a).
In step S33a, the calculation process of ide (k) shown in FIG. 35 is executed to calculate the identification error ide (k). In step S33b, the update vector dθ (k) is calculated by the equation (15e), and then the θbase table shown in FIG. 34 is searched according to the throttle valve opening deviation amount DTH to calculate the reference model parameter vector θbase (step S33b). S33c). Reference model parameters a1base, a2base and b1base are set in the θbase table. When the throttle valve opening deviation amount DTH takes a value near “0” (the throttle valve opening TH is near the default opening THDEF), the reference model parameters a1base and b1base decrease, and the reference model parameter a2base increases. Is set to The reference model parameter c1base is set to “0”.
[0158]
In step S34a, the model parameter vector θ (k) is calculated by the equation (15f), and then, as in the first embodiment, the model parameter vector θ (k) is stabilized (step S35). That is, a limit process of each model parameter is performed to calculate a corrected model parameter vector θL (k).
[0159]
FIG. 35 is a flowchart of the ide (k) calculation process executed in step S33a of FIG. This processing is obtained by deleting step S56 (dead zone processing) of the ide (k) calculation processing in FIG. 12 and replacing step S51 with step S51a. That is, in the present embodiment, since the drift of the model parameters is prevented by the δ correction method, the dead zone processing is not executed.
[0160]
In step S51a, the estimated throttle valve opening deviation amount DTHHAT (k) is calculated by the equations (18) and (19a), and the identification error ide (k) is calculated using the estimated throttle valve opening deviation amount DTHHAT (k). ) Is calculated.
In the present embodiment, since the dead time d of the control target model is “0”, the predetermined value XCNTIDST in step S52 is set to, for example, “2”.
[0161]
FIG. 36 is a flowchart of the process of calculating the control input Usl (= DUT) to the throttle
[0162]
In step S201a, the process of calculating the switching function value σ shown in FIG. 37 is executed, and in step S202a, the process of calculating the integrated value of the switching function value σ shown in FIG. 38 is executed. In step S203a, the equivalent control input Ueq is calculated by the equation (9b). In step S204a, the arithmetic processing of the reaching law input Urch shown in FIG. 39 is executed, and in step S205a, the arithmetic processing of the adaptive law input Uadp shown in FIG. 40 is executed.
[0163]
FIG. 37 is a flowchart of the calculation process of the switching function value σ executed in step S201a of FIG. In this processing, steps S222 to S226 of the calculation processing of the prediction switching function value σpre shown in FIG. 23 are changed to steps S222a to 226a, respectively.
[0164]
In step S222a, the switching function value σ (k) is calculated by the equation (5). In subsequent steps S223a to S226a, “σpre” in steps S223 to S226 in FIG. 23 is replaced with “σ”, and a limit process similar to the process in FIG. 23 is performed on the switching function value σ (k). Do.
[0165]
FIG. 38 is a flowchart of the process of calculating the integrated value SUMSIGMAa of the switching function value σ, which is performed in step S202a of FIG. In this processing, steps S241 to S245 of the integrated value calculation processing of σpre shown in FIG. 26 are changed to steps S241a to S245a, respectively. The integrated value SUMSIGMAa is used to calculate the adaptive law input Uadp in the processing of FIG. 40 described later (see the above equation (11b)).
[0166]
In step S241a, the integrated value SUMSIGMAa is calculated by the following equation (47a).
SUMSIGMAa (k) = SUMSIGMAa (k-1) + σ × ΔT (47a)
In subsequent steps S242a to S245a, a limit process similar to the process in FIG. 26 is performed on the calculated integrated value SUMSIGMAa.
[0167]
FIG. 39 is a flowchart of the reaching law input Urch calculation process executed in step S204a of FIG. In this process, steps S263 and S265 of the reaching law input Urch calculation process shown in FIG. 27 are changed to steps S263a and S265a, respectively.
[0168]
That is, in the present embodiment, the reaching law input Urch when the adaptive sliding
[0169]
FIG. 40 is a flowchart of the calculation process of the adaptive law input Uadp executed in step S205a of FIG. This processing is obtained by changing steps S273 and S275 of the adaptive law input Uadp calculation processing shown in FIG. 28 to steps S273a and S275a, respectively.
[0170]
That is, in the present embodiment, the adaptive law input Uadp when the adaptive sliding
[0171]
FIG. 41 is a flowchart of the stability determination processing of the sliding mode controller executed in step S16a of FIG. This processing is obtained by changing steps S281 and S281 in FIG. 29 to steps S281a and S282a, respectively.
[0172]
In step S281a, the switching function change amount Dσ is calculated by the following equation (52a). In step S282a, the stability determination parameter SGMSTAB is calculated by the following equation (53a). That is, the stability determination is performed based on the switching function value σ instead of the predicted switching function value σpre.
Dσ = σ (k) −σ (k−1) (52a)
SGSTAB = Dσ × σ (k) (53a)
[0173]
In the present embodiment, a part of the throttle
[0174]
(Third embodiment)
FIG. 42 is a block diagram illustrating a configuration of a control system according to the third embodiment of the present invention.
This control system includes a
[0175]
The
The
[0176]
The manipulated
[0177]
Hereinafter, the correspondence between the components and parameters in the present embodiment and the components and parameters in the first embodiment will be described.
pH sensor 102Corresponds to the throttle
[0178]
The second reference value V2BASE is added to bias the center value of the first manipulated variable U1, which is the output of the adaptive sliding
[0179]
The
[0180]
Because of the above correspondence, the
In the present embodiment, the
[0181]
(Modification of Third Embodiment)
FIG. 43 is a diagram showing a modification of the configuration shown in FIG. In this modification, a
[0182]
Thus, the same modeling and the same control method as in the third embodiment can be applied to a plant including a local feedback loop.
In the first embodiment, since the motor drive circuit is a known circuit, a detailed description is not given. However, a current sensor for detecting the output current of the switching element that is controlled to be on / off is provided, and the detected current value ID However, the feedback control may be performed so as to match the current value IR corresponding to the operation amount Usl. This modification corresponds to a case where such a circuit configuration is employed in the first embodiment.
[0183]
(Fourth embodiment)
FIG. 44 is a block diagram illustrating a configuration of a control system according to the fourth embodiment of the present invention. This control system is obtained by replacing the operation
[0184]
The operation
The
[0185]
That is, the
[0186]
In the present embodiment, the
[0187]
(Modification of Fourth Embodiment)
FIG. 45 is a diagram showing a modification of the configuration shown in FIG. In this modification, a
Thus, the same modeling and control method as in the fourth embodiment can be applied to a plant including a local feedback loop.
[0188]
(Other embodiments)
As a method for calculating the identification error ide (k) of the model parameter, the ε correction method described below may be employed instead of the δ correction method. That is, the model parameter vector θ (k) may be calculated by the following equation (15g) instead of the equation (15c).
EPS = [EPS1, EPS2, EPS3, EPS4]
[0189]
The forgetting factors EPS1, EPS2, and EPS4 are set to values between 0 and 1 (0 <EPSi <1), similarly to the forgetting factor DELTAI, and have a function of gradually reducing the influence of past identification errors.
[0190]
However, in the case of the ε-correction method, the coefficient EPS3 related to the calculation of the model parameter b1 must be set to “1”. This is for the following reason. In the case of the ε-correction method, when the identification error ide (k) decreases, all the model parameters take values near zero. However, since the model parameter b1 is applied to the denominator of the equations (9b), (10b), and (11b), when the model parameter b1 approaches “0”, the input Usl to the control target diverges. .
[0191]
Expression (15g) differs from Expression (15c) in that the initial value vector θ (0) is also multiplied by the forgetting coefficient vector EPS.
Rewriting the equation (15g) into the recurrence form gives the following equation (15h). A method of calculating the model parameter vector θ (k) using the following equation (15h) instead of the equation (15) is referred to as an ε correction method.
θ (k) = EPS × θ (k−1) + KP (k) ide (k) (15h)
[0192]
The influence of the past identification error eid is also reduced by the ε-correction method, so that the drift of the model parameters can be prevented.
Further, in the second embodiment, the drift of the model parameters is prevented by the δ correction method. However, similarly to the first embodiment, the corrected identification error idenl (k) is calculated by the dead zone processing (FIG. 14). , The model parameter vector θ (k) may be calculated.
[0193]
In the first embodiment, a δ correction method or an ε correction method may be employed instead of the dead zone processing. Further, when the δ correction method is employed in the first embodiment, a model parameter scheduler is introduced as in the second embodiment, and an update vector is added to the reference model parameter vector θbase calculated by the model parameter scheduler. It is desirable to calculate the model parameter vector θ in such a format as follows.
[0194]
【The invention's effect】
As described in detail above, according to the first aspect of the present invention, an update vector is calculated according to the identification error of the model parameter vector, and the update vector is corrected so as to reduce the influence of the past value of the identification error, By adding the modified update vector to the reference vector of the model parameter, the model parameter vector is calculated, and the values of the elements of the model parameter vector are further restricted within a predetermined limit range. Therefore, the stability of the control can be further improved while preventing the drift of the model parameters.In addition, since the reference vector of the model parameter is calculated according to the parameter indicating the dynamic characteristic change of the plant, an appropriate reference vector according to the dynamic characteristic change of the plant can be obtained. As a result, it is possible to quickly converge the model parameters even when a non-linear element is included in the plant.
[0195]
According to the second aspect of the present invention, since the update vector is calculated using the fixed gain algorithm, the amount of calculation can be reduced.
[0196]
According to the third aspect of the present invention, the update vector is corrected by multiplying the past value of at least one element of the update vector by a predetermined value greater than 0 and smaller than 1. The influence of the value is reduced, and the drift of the model parameter vector can be prevented.
[0197]
According to the fourth aspect of the present invention, the update vector is not multiplied by a predetermined value larger than 0 and smaller than 1 for elements related to plant input or elements not related to plant input / output. Can prevent the occurrence of a steady-state error.
[0200]
Claim5According to the invention described in (1), the identification error of the model parameter isDead zoneWhen it is within the range, the identification error is corrected in the decreasing direction, the model parameter vector is calculated using the corrected identification error, and the values of the elements of the model parameter vector are limited to within a predetermined limit range. Therefore, the stability of the control can be further improved while preventing the drift of the model parameters.Further, since the dead zone is set according to the output of the plant or the control target value, it is possible to prevent the identification error to be reflected in the value of the model parameter from being ignored as the residual identification error.
[0201]
Claim6According to the invention described in the above, since the values of the plurality of elements are limited so that the plurality of elements of the model parameter vector satisfy the predetermined relationship, it is possible to improve the stability of control using the model parameter vector. it can.
[0202]
Claim7According to the invention described in the above, the identification error isDead zoneWhen the value is within the range, the identification error is set to 0, so that the effect of the identification error that should not be reflected in the value of the model parameter can be eliminated, and the drift prevention effect of the model parameter can be enhanced.
[0203]
Claim8According to the invention described in the above, since the control for matching the throttle valve opening to the target opening is performed using the model parameters identified by the identification means, the controllability of the throttle valve opening to the target opening is improved. The stability of throttle valve opening control can be improved.
[Brief description of the drawings]
FIG. 1 is a diagram showing a throttle valve driving device for an internal combustion engine and a control device thereof according to an embodiment of the present invention.
FIG. 2 is a diagram showing a frequency characteristic of the throttle valve driving device shown in FIG.
FIG. 3 is a functional block diagram showing functions realized by an electronic control unit (ECU) of FIG.
FIG. 4 is a diagram showing a relationship between control characteristics of a sliding mode controller and values of a switching function setting parameter (VPOLE).
FIG. 5 is a diagram showing a setting range of a control gain (F, G) of a sliding mode controller.
FIG. 6 is a diagram for explaining drift of model parameters.
FIG. 7 is a diagram showing a function for correcting an identification error.
FIG. 8 is a diagram for explaining that a default opening deviation of a throttle valve is reflected on a model parameter (c1 ′).
FIG. 9 is a flowchart of a throttle valve opening control process.
FIG. 10 is a flowchart of a process for setting a state variable in the process of FIG. 9;
FIG. 11 is a flowchart of a process of executing a calculation of a model parameter identifier in the process of FIG. 9;
FIG. 12 is a flowchart of a process of executing an operation of an identification error (ide) in the process of FIG. 11;
FIG. 13 is a diagram for explaining low-pass filtering of an identification error (ide).
FIG. 14 is a flowchart of a dead zone process in the process of FIG.
FIG. 15 is a diagram showing a table used in the processing of FIG. 14;
FIG. 16 is a flowchart of a process of stabilizing a model parameter vector (θ) in the process of FIG. 11;
17 is a flowchart of a limit process of model parameters (a1 ', a2') in the process of FIG.
FIG. 18 is a diagram for explaining a change in a value of a model parameter due to the processing of FIG. 16;
FIG. 19 is a flowchart of a limit process of a model parameter (b1 ′) in the process of FIG. 16;
20 is a flowchart of a limit process of a model parameter (c1 ') in the process of FIG.
FIG. 21 is a flowchart of a process of executing a calculation of a state predictor in the process of FIG. 9;
FIG. 22 is a flowchart of a process for executing a calculation of a control input (Usl) in the process of FIG. 9;
FIG. 23 is a flowchart of a process of executing a calculation of a predictive switching function value (σpre) in the process of FIG. 22;
FIG. 24 is a flowchart of a process of executing a calculation of a switching function setting parameter (VPOLE) in the process of FIG. 23;
FIG. 25 is a diagram showing a map used in the processing of FIG. 24;
FIG. 26 is a flowchart of a process of executing an integrated value of a prediction switching function value (σpre) in the process of FIG. 22;
FIG. 27 is a flowchart of a process of executing a calculation of a reaching law input (Urch) in the process of FIG. 22;
FIG. 28 is a flowchart of a process of executing an adaptive law input (Uadp) operation in the process of FIG. 22;
FIG. 29 is a flowchart of a process for executing a stability determination of the sliding mode controller in the process of FIG. 9;
FIG. 30 is a flowchart of a process of executing a calculation of a default opening deviation (thdefadp) in the process of FIG. 9;
FIG. 31 is a functional block diagram showing functions realized by an electronic control unit (ECU) in FIG. 1 (second embodiment).
FIG. 32 is a flowchart of a throttle valve opening control process (second embodiment).
FIG. 33 is a flowchart of a process of executing a calculation of a model parameter identifier in the process of FIG. 32;
FIG. 34 is a diagram showing a table used in the processing in FIG. 33;
FIG. 35 is a flowchart of a process of executing an operation of an identification error (ide) in the process of FIG. 33;
FIG. 36 is a flowchart of a process for executing a calculation of a control input (Usl) in the process of FIG. 32;
FIG. 37 is a flowchart of a process of executing a calculation of a switching function value (σ) in the process of FIG. 36;
FIG. 38 is a flowchart of a process of executing the calculation of the integrated value of the switching function value (σ) in the process of FIG. 36;
39 is a flowchart of a process of executing a calculation of a reaching law input (Urch) in the process of FIG. 36.
40 is a flowchart of a process of executing an adaptive law input (Uadp) operation in the process of FIG. 36.
FIG. 41 is a flowchart of a process for performing a stability determination of the sliding mode controller in the process of FIG. 32;
FIG. 42 is a block diagram illustrating a configuration of a control system according to a third embodiment of the present invention.
FIG. 43 is a block diagram showing a modification of the configuration shown in FIG. 42;
FIG. 44 is a block diagram illustrating a configuration of a control system according to a fourth embodiment of the present invention.
FIG. 45 is a block diagram showing a modification of the configuration shown in FIG. 44.
[Explanation of symbols]
1 Internal combustion engine
3 Throttle valve
7 Electronic control unit
10 Throttle valve drive
21, 21a Adaptive sliding mode controller (control means)
22, 22a Model parameter identifier (identification means, identification error calculation means, update vector calculation means, update vector correction means, identification error correction means, model parameter vector calculation means, limiting means)
24 Target opening setting section
25 Model Parameter Scheduler
Claims (8)
前記同定手段は、前記モデルパラメータベクトルの同定誤差を算出する同定誤差算出手段と、該同定誤差に応じて更新ベクトルを算出する更新ベクトル算出手段と、前記同定誤差の過去値の影響を減少させるように前記更新ベクトルを修正する更新ベクトル修正手段と、前記プラントの動特性変化を示すパラメータに応じて前記モデルパラメータの基準ベクトルを算出する基準ベクトル算出手段と、前記基準ベクトルに修正された更新ベクトルを加算することにより、前記モデルパラメータベクトルを算出するモデルパラメータベクトル算出手段と、前記モデルパラメータベクトル算出手段により算出されるモデルパラメータベクトルの要素の値を、所定リミット範囲内に制限する制限手段とを有することを特徴とするプラントの制御装置。Identification means for identifying a model parameter vector of a control target model obtained by modeling a plant based on the input and output of the plant; and control means for controlling the plant using the model parameter vector identified by the identification means. In the control device of the plant equipped with
The identification unit includes an identification error calculation unit that calculates an identification error of the model parameter vector, an update vector calculation unit that calculates an update vector according to the identification error, and reduces an influence of a past value of the identification error. wherein the updating vector correction means for correcting the update vector, and the reference vector calculating means for calculating the reference vector of the model parameters according to the parameter indicating the dynamic characteristic changes of the plant, were modified before Kimoto quasi vector updating the Model parameter vector calculating means for calculating the model parameter vector by adding the vectors, and limiting means for limiting the values of the elements of the model parameter vector calculated by the model parameter vector calculating means to within a predetermined limit range. A control device for a plant, comprising:
前記同定手段は、前記モデルパラメータベクトルの同定誤差を算出する同定誤差算出手段と、前記プラントの出力または制御目標値に応じて不感帯を設定する不感帯設定手段と、該同定誤差が前記不感帯内にあるときは、前記同定誤差を減少方向に修正する同定誤差修正手段と、該同定誤差修正手段により修正された同定誤差を用いて前記モデルパラメータベクトルを算出するモデルパラメータベクトル算出手段と、前記モデルパラメータベクトル算出手段により算出されるモデルパラメータベクトルの要素の値を、所定リミット範囲内に制限する制限手段とを有することを特徴とするプラントの制御装置。Identification means for identifying a model parameter vector of a control target model obtained by modeling the plant based on the input and output of the plant, and using the model parameter vector identified by the identification means, the output of the plant is controlled by the control target. Control device for controlling the plant so as to match the value, in a plant control device,
The identification means, an identification error calculation means for calculating an identification error of the model parameter vector, a dead zone setting means for setting a dead zone according to an output of the plant or a control target value, and the identification error is within the dead zone . An identification error correction unit for correcting the identification error in a decreasing direction; a model parameter vector calculation unit for calculating the model parameter vector using the identification error corrected by the identification error correction unit; and the model parameter vector. A plant control device comprising: a limiter that limits a value of an element of a model parameter vector calculated by a calculator within a predetermined limit range.
Priority Applications (15)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2001184540A JP3602811B2 (en) | 2001-06-19 | 2001-06-19 | Plant control equipment |
| EP10181301A EP2261759B8 (en) | 2001-04-20 | 2002-04-19 | Plant control apparatus |
| EP10181400A EP2261760B1 (en) | 2001-04-20 | 2002-04-19 | Plant control apparatus |
| CA002411520A CA2411520C (en) | 2001-04-20 | 2002-04-19 | Control system for plant |
| AT10181400T ATE538418T1 (en) | 2001-04-20 | 2002-04-19 | ARRANGEMENT FOR CONTROLLING APPARATUS |
| PCT/JP2002/003895 WO2002086630A1 (en) | 2001-04-20 | 2002-04-19 | Plant control apparatus |
| BR122012016734-7A BR122012016734B1 (en) | 2001-04-20 | 2002-04-19 | PLANT CONTROL SYSTEM |
| EP10181437A EP2264302B1 (en) | 2001-04-20 | 2002-04-19 | Plant control apparatus |
| CNB028010086A CN1275108C (en) | 2001-04-20 | 2002-04-19 | Plant control apparatus |
| AT10181437T ATE554277T1 (en) | 2001-04-20 | 2002-04-19 | ARRANGEMENT FOR CONTROLLING APPARATUS |
| BR0205088-9A BR0205088A (en) | 2001-04-20 | 2002-04-19 | Plant control system |
| US10/312,108 US7050864B2 (en) | 2001-04-20 | 2002-04-19 | Control system for a plant using identified model parameters |
| KR1020027016490A KR100552579B1 (en) | 2001-04-20 | 2002-04-19 | Plant control apparatus |
| EP02720505.3A EP1293851B1 (en) | 2001-04-20 | 2002-04-19 | Plant control apparatus |
| US11/345,994 US7216006B2 (en) | 2001-04-20 | 2006-02-02 | Control system for a plant including a slide mode controller |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2001184540A JP3602811B2 (en) | 2001-06-19 | 2001-06-19 | Plant control equipment |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2003005804A JP2003005804A (en) | 2003-01-08 |
| JP3602811B2 true JP3602811B2 (en) | 2004-12-15 |
Family
ID=19024286
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2001184540A Expired - Lifetime JP3602811B2 (en) | 2001-04-20 | 2001-06-19 | Plant control equipment |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP3602811B2 (en) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4391789B2 (en) * | 2003-10-03 | 2009-12-24 | 本田技研工業株式会社 | Control device for controlling a plant with an identifier for partially identifying model parameters |
| JP4326386B2 (en) * | 2004-03-26 | 2009-09-02 | 本田技研工業株式会社 | Control device |
| JP4282572B2 (en) | 2004-08-30 | 2009-06-24 | 本田技研工業株式会社 | Control device for controlling the plant |
| JP4952025B2 (en) * | 2006-03-31 | 2012-06-13 | 株式会社日立製作所 | Operation control method, operation control apparatus, and operation control system |
| JP4848320B2 (en) * | 2007-06-25 | 2011-12-28 | 株式会社豊田中央研究所 | Parameter identification apparatus and program thereof |
-
2001
- 2001-06-19 JP JP2001184540A patent/JP3602811B2/en not_active Expired - Lifetime
Also Published As
| Publication number | Publication date |
|---|---|
| JP2003005804A (en) | 2003-01-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| KR100552579B1 (en) | Plant control apparatus | |
| JP4223489B2 (en) | Plant control device | |
| US6668200B2 (en) | Control system for plant | |
| JP3977199B2 (en) | Control device for throttle valve drive device | |
| JP4064159B2 (en) | Plant control equipment | |
| JP3863069B2 (en) | Plant control equipment | |
| JP3602811B2 (en) | Plant control equipment | |
| JP4002460B2 (en) | Control device for throttle valve drive device | |
| JP4549569B2 (en) | Plant control equipment | |
| JP3995899B2 (en) | Plant control equipment | |
| JP3902504B2 (en) | Plant control equipment | |
| JP3602809B2 (en) | Plant control equipment | |
| JP3857169B2 (en) | Plant control equipment | |
| JP3686383B2 (en) | Plant control device | |
| JP3819258B2 (en) | Control device for throttle valve drive device | |
| JP3701600B2 (en) | Plant control device | |
| JP3819257B2 (en) | Control device for throttle valve drive device | |
| JP3819259B2 (en) | Control device for throttle valve drive device | |
| JP3949470B2 (en) | Plant control equipment | |
| JP4464010B2 (en) | Control device for throttle valve drive device | |
| JP4435467B2 (en) | Control device for throttle valve drive device | |
| JP3964284B2 (en) | Control device for throttle valve drive device | |
| EP1331532A2 (en) | Control system for plant | |
| BRPI0205088B1 (en) | PLANT CONTROL SYSTEM |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20040518 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20040707 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20040907 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20040924 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 Ref document number: 3602811 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20071001 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20081001 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20081001 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20091001 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20091001 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20101001 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20101001 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111001 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111001 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121001 Year of fee payment: 8 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20131001 Year of fee payment: 9 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| EXPY | Cancellation because of completion of term |