(57)【特許請求の範囲】
【請求項１】
内燃機関の制御対象を運転状態に応じて模擬した制御モデルを用いて該制御対象を操作するアクチュエータの操作量をフィードバック制御する内燃機関の制御装置において、
前記制御対象の状態を検出する状態検出手段と、
前記アクチュエータの現在及び過去の操作量並びに前記状態検出手段で検出した現在及び過去の状態検出値を前記制御モデルの内部状態を代表する状態変数量として出力する状態変数量出力手段と、
前記運転状態に応じて模擬した制御モデルのモデルパラメータを用いて制御パラメータを演算する制御パラメータ演算手段と、
前記制御パラメータ演算手段で前記制御モデルのモデルパラメータを用いて演算した制御パラメータと、前記状態変数量出力手段から出力される状態変数量の差分値と、前記状態検出手段で検出した検出値と制御目標値との偏差と、過去の操作量補正値とに基づいて前記アクチュエータの操作量を補正する今回の操作量補正値を演算する補正値演算手段と、
前記補正値演算手段で演算した前記今回の操作量補正値を前回の操作量に加算して今回の操作量を求める操作量演算手段と
を備えていることを特徴とする内燃機関の制御装置。
【請求項２】
前記制御パラメータ演算手段は、前記モデルパラメータを運転条件に応じて切り換えることを特徴とする請求項１に記載の内燃機関の制御装置。
【請求項３】
前記制御パラメータ演算手段は、極配置法に基づいて前記制御モデルの特性多項式の係数を演算し、この特性多項式の係数と前記モデルパラメータとから前記制御パラメータを演算することを特徴とする請求項１又は２に記載の内燃機関の制御装置。
【請求項４】
前記制御パラメータ演算手段は、目標とする極を運転条件に応じて切り換えることを特徴とする請求項３に記載の内燃機関の制御装置。
【請求項５】
前記運転条件として、内燃機関の回転速度、吸入空気量、負荷、冷却水温、始動後経過時間のいずれかを用いることを特徴とする請求項２又は４に記載の内燃機関の制御装置。
【請求項６】
前記制御パラメータ演算手段は、高空気量運転時は目標とする極を高応答に設定し、低空気量運転時は目標とする極を低応答に設定することを特徴とする請求項３乃至５のいずれかに記載の内燃機関の制御装置。
【請求項７】
前記操作量演算手段は、内燃機関の空燃比フィードバック制御系を模擬した制御モデルを用いて空燃比補正係数を演算することを特徴とする請求項１乃至６のいずれかに記載の内燃機関の制御装置。
【請求項８】
前記操作量演算手段は、内燃機関のアイドル回転速度制御系を模擬した制御モデルを用いて、アイドル回転速度を操作するアクチュエータの操作量を演算することを特徴とする請求項１乃至６のいずれかに記載の内燃機関の制御装置。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、燃料噴射弁から空燃比検出手段までの制御対象を模擬した制御モデルを用いて燃料噴射量（空燃比）を制御する内燃機関の制御装置に関するものである。
【０００２】
【従来の技術】
近年の自動車は、排気管に排出ガス浄化用の三元触媒を設置すると共に、この三元触媒の上流側に空燃比センサを設置し、この空燃比センサの出力に基づいて排出ガスの空燃比を触媒の浄化ウインドウ（理論空燃比付近）に制御するように状態フィードバックを実行して燃料噴射量を制御することで、排出ガスを効率良く浄化するようにしている。このような空燃比制御は、一般に、特開平７−１１９９５号公報に示すように、燃料噴射弁から空燃比センサまでの制御対象をモデル化して、最適レギュレータにより状態フィードバックのフィードバックゲインを算出すると共に、このフィードバックゲインを用いて空燃比補正係数を算出し、エンジン運転状態に応じて求めた基本噴射量を上記空燃比補正係数等で補正して燃料噴射量を求めるようにしている。
【０００３】
【発明が解決しようとする課題】
上記従来の空燃比制御は、エンジン運転状態に応じてフィードバックゲインを連続的に変化させることができず、制御系の安定化のために低いフィードバックゲインで制御せざるを得ず、その分、空燃比制御精度が低下するという欠点があった。
【０００４】
本発明の第１の目的は、内燃機関のフィードバック制御系の制御パラメータ（フィードバックゲイン）を運転状態に応じて連続的に変化させることができ、制御精度を向上することができる内燃機関の制御装置を提供することにある。
【０００５】
更に、第２の目的は、内燃機関のフィードバック制御系の制御パラメータ（フィードバックゲイン）をオンラインリアルタイム処理で算出することができる内燃機関の制御装置を提供することにある。
【０００６】
【課題を解決するための手段】
一般に、状態フィードバックの制御パラメータ（フィードバックゲイン）Ｆ1 〜Ｆd+1 ，Ｆ０に基づいて空燃比補正係数ＦＡＦ(i) を算出する場合は、次式が用いられることが多い。
0003758134.tif 000002
ここで、λ(i) は現在の空燃比（空気過剰率）、ＦＡＦ(i-1) 〜ＦＡＦ(i-d) は過去の空燃比補正係数、λref は目標空燃比（目標空気過剰率）である。
【０００７】
しかし、この空燃比補正係数の演算方法では、制御パラメータＦ1 〜Ｆd+1 ，Ｆ０を運転条件等に応じて切り換えると、その瞬間に、空燃比補正係数ＦＡＦが一時的に乱れ、その結果、空燃比λが一時的に乱れる現象が発生する可能性がある。
【０００８】
そこで、制御パラメータ演算手段で演算した制御パラメータと、空燃比検出手段で検出した空燃比の変化量と、該空燃比と目標空燃比との偏差と、過去の空燃比補正係数補正値とに基づいて今回の空燃比補正係数補正値ΔＦＡＦ(i) を補正値演算手段により演算し、今回の空燃比補正係数補正値ΔＦＡＦ(i) を前回の空燃比補正係数ＦＡＦ(i-1) に加算して今回の空燃比補正係数ＦＡＦ(i) を求めるようにすると良い。
ＦＡＦ(i) ＝ＦＡＦ(i-1) ＋ΔＦＡＦ(i)
【０００９】
このようにすれば、制御パラメータを運転条件等に応じて切り換えても、空燃比補正係数が乱れることがなくなり、空燃比が乱れる現象が発生しなくなる。これにより、制御パラメータを運転条件等に応じて切り換えながら、安定した空燃比制御を行うことが可能となる。
【００１０】
ところで、内燃機関のフィードバック制御系は、空燃比フィードバック制御系の他に、例えば、アイドル回転速度制御系等がある。
これらのフィードバック制御系に関しては、請求項１のように、制御対象の状態を検出する状態検出手段と、アクチュエータの現在及び過去の操作量並びに前記状態検出手段で検出した現在及び過去の状態検出値を制御モデルの内部状態を代表する状態変数量として出力する状態変数量出力手段と、運転状態に応じて模擬した制御モデルのモデルパラメータを用いて制御パラメータ（フィードバックゲイン）を演算する制御パラメータ演算手段とを備え、前記制御パラメータ演算手段で前記制御モデルのモデルパラメータを用いて演算した制御パラメータと、前記状態変数量出力手段から出力される状態変数量の差分値と、前記状態検出手段で検出した検出値と制御目標値との偏差と、過去の操作量補正値とに基づいて補正値演算手段によってアクチュエータの操作量を補正する今回の操作量補正値を演算すると共に、操作量演算手段によって、この今回の操作量補正値を前回の操作量に加算して今回の操作量を求めるようにすると良い。このようにすれば、制御対象の制御パラメータを運転状態等に応じて切り換えても、制御対象の状態が乱れることがなくなり、制御パラメータを運転条件等に応じて切り換えながら、制御対象を安定して制御することが可能となる。
【００１１】
この場合、請求項２のように、制御モデルのモデルパラメータを運転条件に応じて切り換えるようにすると良い。このようにすれば、モデルパラメータから算出する制御パラメータを、運転条件に応じて切り換えることができる。
【００１２】
また、請求項３のように、極配置法に基づいて制御モデルの特性多項式の係数を演算し、この特性多項式の係数とモデルパラメータとから制御パラメータを演算するようにしても良い。このようにすれば、オンラインリアルタイム処理で特性多項式の係数を演算して制御パラメータを更新することができるため、制御パラメータを運転条件に応じて連続的に変化させることができ、運転条件の変化に対する制御特性を向上することができる。
【００１３】
更に、請求項４のように、目標とする極を運転条件に応じて切り換えるようにしても良い。このようにすれば、制御モデルの特性多項式の係数を運転条件に応じて変化させて、制御パラメータを運転条件に応じて変化させることが可能となる。
【００１４】
この場合、請求項５のように、運転条件として、内燃機関の回転速度、吸入空気量、負荷、冷却水温、始動後経過時間のいずれかを用いるようにしても良い。これら内燃機関の回転速度、吸入空気量、負荷、冷却水温、始動後経過時間は、いずれも内燃機関のフィードバック制御系に影響を与える主要な運転パラメータである。
【００１５】
更に、請求項６のように、高空気量運転時は目標とする極を高応答に設定し、低空気量運転時は目標とする極を低応答に設定しても良い。このようにすれば、フィードバック制御系の応答性と安定性とを両立させることができる。
【００１６】
以上説明した請求項１〜６に係る発明は、請求項７，８のように、空燃比フィードバック制御系、アイドル回転速度制御系のいずれに適用しても良く、要は、内燃機関の制御対象をフィードバック制御する制御系に適用して実施できる。
【００１７】
【発明の実施の形態】
《実施形態（１）》
以下、本発明の実施形態（１）を図１乃至図１３に基づいて説明する。
【００１８】
まず、図１に基づいてエンジン制御システム全体の概略構成を説明する。内燃機関であるエンジン１１の吸気管１２の最上流部には、エアクリーナ１３が設けられ、このエアクリーナ１３の下流側には、吸入空気量を検出するエアフローメータ１４が設けられている。このエアフローメータ１４の下流側には、スロットルバルブ１５とスロットル開度を検出するスロットル開度センサ１６とが設けられている。
【００１９】
更に、スロットルバルブ１５の下流側には、サージタンク１７が設けられ、このサージタンク１７に、吸気管圧力を検出する吸気管圧力センサ１８が設けられている。また、サージタンク１７には、エンジン１１の各気筒に空気を導入する吸気マニホールド１９が設けられ、各気筒の吸気マニホールド１９の吸気ポート近傍に、それぞれ燃料を噴射する燃料噴射弁２０が取り付けられている。
【００２０】
一方、エンジン１１の排気管２１の途中には、排出ガス中の有害成分（ＣＯ，ＨＣ，ＮＯｘ等）を低減させる三元触媒等の触媒２２が設置されている。この触媒２２の上流側には、排出ガスの空燃比を検出するリニアＡ／Ｆセンサ等の空燃比センサ２３（空燃比検出手段，状態検出手段）が設けられている。また、エンジン１１のシリンダブロックには、冷却水温を検出する冷却水温センサ２４や、エンジン回転速度を検出するクランク角センサ２５が取り付けられている。
【００２１】
これら各種のセンサ出力は、エンジン制御回路（以下「ＥＣＵ」と表記する）２６に入力される。このＥＣＵ２６は、マイクロコンピュータを主体として構成され、内蔵されたＲＯＭ（記憶媒体）に記憶された後述する各プログラムを実行することで、空燃比補正係数ＦＡＦを演算して燃料噴射弁２０の燃料噴射量を制御する。
【００２２】
一般に、状態フィードバックの制御パラメータ（フィードバックゲイン）Ｆ１〜Ｆd+1 ，Ｆ０に基づいて空燃比補正係数ＦＡＦ(i) を算出する場合は、次式が用いられることが多い。
0003758134.tif 000003
は過去の空燃比補正係数、λref は目標空燃比（目標空気過剰率）である。
【００２３】
しかし、この空燃比補正係数の演算方法では、制御パラメータＦ１〜Ｆd+1 ，Ｆ０を運転条件等に応じて切り換えると、その瞬間に、空燃比補正係数ＦＡＦが一時的に乱れ、その結果、空燃比λが一時的に乱れる現象が発生する可能性がある。
【００２４】
そこで、本実施形態（１）では、今回の空燃比補正係数補正値ΔＦＡＦ(i) を算出し、今回の空燃比補正係数補正値ΔＦＡＦ(i) を前回の空燃比補正係数ＦＡＦ(i-1) に加算して今回の空燃比補正係数ＦＡＦ(i) を求める。
ＦＡＦ(i) ＝ＦＡＦ(i-1) ＋ΔＦＡＦ(i)
この場合、空燃比補正係数ＦＡＦが特許請求の範囲でいう操作量に相当し、空燃比補正係数補正値ΔＦＡＦが特許請求の範囲でいう操作量補正値に相当する。
【００２５】
今回の空燃比補正係数補正値ΔＦＡＦ(i) は、次式により算出する。
ΔＦＡＦ(i) ＝Ｆ１・Δφ(i) ＋Ｆ２・ΔＦＡＦ(i-1) ＋…
…＋Ｆd+1 ・ΔＦＡＦ(i-d) ＋Ｆd+2 ・ΔＦＡＦ(i-d-1)
＋Ｆ０・（φref −φ(i) ）
ここで、Δφ(i) は燃料過剰率の変化量、つまりΔφ(i) ＝φ(i) −φ(i-1) であり、このΔφ(i) が特許請求の範囲でいう状態変数量の差分値に相当する。また、ΔＦＡＦ(i-1) 〜ΔＦＡＦ(i-d-1) は過去の空燃比補正係数補正値であり、φref は目標燃料過剰率である。尚、上式では、空燃比の代用情報として燃料過剰率φを用いたが、空気過剰率λを用いても良いことは言うまでもない。
【００２６】
上式を用いて空燃比補正係数ＦＡＦを演算すれば、状態フィードバックの制御パラメータＦ１〜Ｆd+2 ，Ｆ０を運転条件等に応じて切り換えても、空燃比補正係数ＦＡＦが乱れることがなくなり、空燃比が乱れる現象が発生しなくなる。これにより、制御パラメータＦ１〜Ｆd+2 ，Ｆ０を運転条件等に応じて切り換えながら、安定した空燃比制御を行うことが可能となる。
【００２７】
図２は、上式を用いて空燃比補正係数ＦＡＦを演算する空燃比フィードバック制御系の各部の機能を表す機能ブロック図である。この空燃比フィードバック制御系の各機能は、ＥＣＵ２６が実行する図３乃至図１０の各プログラムによって実現される。以下、これら各プログラムの処理内容を説明する。
【００２８】
［燃料噴射量演算］
図３の燃料噴射量演算プログラムは、各気筒の噴射タイミングに同期して起動され、次のようにして燃料噴射量ＴＡＵを算出する。まず、ステップ３０１で、現在のエンジン運転状態に応じてマップ等により基本噴射量Ｔｐを算出する。この後、ステップ３０２で、基本噴射量Ｔｐに対する各種の補正係数ＦＡＬＬ（例えば冷却水温による補正係数、加減速時の補正係数等）を算出し、次のステップ３０３で、空燃比フィードバック条件が成立しているか否かを判定する。もし、空燃比フィードバック条件が成立していなければ、空燃比補正係数ＦＡＦを「１」にセットして、オープンループ制御により空燃比を制御する。
【００２９】
一方、空燃比フィードバック条件が成立していれば、ステップ３０５に進み、排出ガスの空燃比を触媒２２の浄化ウインドウ（理論空燃比付近）に収めるように目標燃料過剰率φref を設定し、次のステップ３０６で、後述する図１０のＦＡＦ演算プログラムを実行して空燃比補正係数ＦＡＦを算出する。
【００３０】
以上のようにして、ステップ３０４又は３０６で空燃比補正係数ＦＡＦを設定した後、ステップ３０７に進み、基本噴射量Ｔｐに空燃比補正係数ＦＡＦと各種補正係数ＦＡＬＬを乗算して燃料噴射量ＴＡＵを求める。これにより、排出ガスの空燃比を触媒２２の浄化ウインドウ内に制御する。
【００３１】
［制御対象特性値演算］
図４の制御対象特性値演算プログラムは、各気筒の噴射タイミングに同期して起動され、制御対象の特性値であるモデル時定数Ｔとむだ時間Ｌを次のようにして算出する。まず、ステップ３０１で、吸入空気量Ｑａを読み込み、次のステップ４０２で、基本モデル時定数Ｔsen と基本むだ時間Ｌsen を、それぞれ吸入空気量Ｑａをパラメータとするマップ等により算出する。
【００３２】
この後、ステップ４０３に進み、負荷（吸入空気量／エンジン回転速度）と冷却水温ＴＨＷを読み込んだ後、ステップ４０４に進み、時定数補正係数α1 とむだ時間補正係数α2 を、それぞれ負荷と冷却水温ＴＨＷをパラメータとするマップにより算出する。尚、この補正係数α1 ，α2 の算出マップに用いる運転パラメータは、負荷と冷却水温ＴＨＷの他に、エンジン回転速度や始動後経過時間を用いるようにしても良い。
【００３３】
各補正係数α1 ，α2 を算出した後、ステップ４０５に進み、基本モデル時定数Ｔsen 、基本むだ時間Ｌsen と、それぞれの補正係数α1 ，α2 を用いて、次式によりモデル時定数Ｔとむだ時間Ｌを算出して、本プログラムを終了する。
Ｔ＝（１＋α1 ）・Ｔsen
Ｌ＝（１＋α1 ）・Ｌsen
【００３４】
［噴射間隔演算］
図５の噴射間隔演算プログラムは、各気筒の噴射タイミングに同期して起動され、噴射間隔ｄｔを次のようにして算出する。まず、ステップ４１１で、エンジン回転速度Ｎｅ（ｒｐｍ）を読み込んだ後、ステップ４１２に進み、噴射間隔ｄｔを次式により算出して、本プログラムを終了する。
ｄｔ＝３０／Ｎｅ×気筒数
【００３５】
［減衰係数ζ、固有角周波数ω演算］
図６の減衰係数ζ、固有角周波数ω演算プログラムは、各気筒の噴射タイミングに同期して起動され、極配置法の演算に用いる減衰係数ζと固有角周波数ωを次のようにして算出する。まず、ステップ４２１で、吸入空気量Ｑａを読み込み、次のステップ４２２で、基本減衰係数ζsen と基本固有角周波数ωsen を、それぞれ吸入空気量Ｑａをパラメータとするマップにより算出する。
【００３６】
この後、ステップ４２３に進み、負荷（吸入空気量／エンジン回転速度）と冷却水温ＴＨＷを読み込んだ後、ステップ４２４に進み、減衰係数補正係数α3 と固有角周波数補正係数α4 を、それぞれ負荷と冷却水温ＴＨＷをパラメータとするマップにより算出する。尚、この補正係数α3 ，α4 の算出マップに用いる運転パラメータは、負荷と冷却水温ＴＨＷの他に、エンジン回転速度や始動後経過時間を用いるようにしても良い。
【００３７】
各補正係数α3 ，α4 の算出後、ステップ４２５に進み、基本減衰係数ζsen 、基本固有角周波数ωsen と、それぞれの補正係数α3 ，α4 を用いて、次式により減衰係数ζと固有角周波数ωを算出して本プログラムを終了する。
ζ＝（１＋α3 ）・ζsen
ω＝（１＋α4 ）・ωsen
【００３８】
この場合、減衰係数ζと固有角周波数ωは、特許請求の範囲でいう目標とする極に相当する。本実施形態（１）では、減衰係数ζと固有角周波数ωは、高空気量運転時には高応答に設定し、低空気量運転時には低応答に設定する。これにより、空燃比フィードバック制御系の応答性と安定性とを両立させることが可能となる。
【００３９】
［モデルパラメータ演算］
図７のモデルパラメータ演算プログラムは、各気筒の噴射タイミングに同期して起動されて、モデルパラメータａ，ｂ1 ，ｂ2 を次のようにして算出する。まず、ステップ４３１で、モデル時定数Ｔ、むだ時間Ｌ、噴射間隔ｄｔを読み込み、次のステップ４３２で、噴射間隔ｄｔ（演算間隔）を基準にして換算したむだ時間ｄ（＝Ｌ／ｄｔ）を小数点以下を切り捨てて求めると共に、その切り捨て誤差Ｌ1 （＝Ｌ−ｄ・ｄｔ）を算出する。
【００４０】
この後、ステップ４３３に進み、モデル時定数Ｔと噴射間隔ｄｔを用いてモデルパラメータａを次式により算出する。
ａ＝ｅｘｐ（−ｄｔ／Ｔ）
【００４１】
この演算は、高性能のＣＰＵを必要とするため、現在の車載コンピュータのＣＰＵの演算能力では、ｅｘｐ（−ｄｔ／Ｔ）の演算を高速で行うことは困難であると思われる。そこで、本実施形態（１）では、演算負荷を軽減するために、ｄｔ／Ｔが例えば０．３５以下の時は、ｅｘｐ（−ｄｔ／Ｔ）を次式により近似し、この近似式によりモデルパラメータａを算出する。
ａ＝１−ｄｔ／Ｔ＋０．５（ｄｔ／Ｔ）2 
【００４２】
この近似式は、ｄｔ／Ｔが大きくなるに従って、演算誤差が大きくなるため、ｄｔ／Ｔが例えば０．３５よりも大きい領域では、予め、ｄｔ／Ｔとモデルパラメータａとの関係をテーブル化してＲＯＭに記憶しておき、このテーブルを検索して、現在のｄｔ／Ｔに応じたモデルパラメータａを求める。尚、ｄｔ／Ｔが、０．３５以下の時にも、予め設定したテーブルからモデルパラメータａを求めるようにしても良い。
【００４３】
この後、ステップ４３４に進み、モデルパラメータｂ1 ，ｂ2 の算出に用いる変数βを次式により算出する。
β＝ｅｘｐ｛−（ｄｔ−Ｌ1 ）／Ｔ｝
この変数βを算出する際も、演算負荷を軽減するために、（ｄｔ−Ｌ1 ）／Ｔが例えば０．３５以下の時は、ｅｘｐ｛−（ｄｔ−Ｌ1 ）／Ｔ｝を次式により近似し、この近似式により変数βを算出する。
β＝１−（ｄｔ−Ｌ1 ）／Ｔ＋０．５｛（ｄｔ−Ｌ1 ）／Ｔ｝2 
【００４４】
この近似式は、（ｄｔ−Ｌ1 ）／Ｔが大きくなるに従って、演算誤差が大きくなるため、（ｄｔ−Ｌ1 ）／Ｔが例えば０．３５よりも大きい領域では、予め、（ｄｔ−Ｌ1 ）／Ｔと変数βとの関係をテーブル化してＲＯＭに記憶しておき、このテーブルを検索して現在の（ｄｔ−Ｌ1 ）／Ｔに応じた変数βを求める。尚、（ｄｔ−Ｌ1 ）／Ｔが０．３５以下の時にも、予め設定したテーブルから変数βを求めるようにしても良い。
【００４５】
この後、ステップ４３５に進み、変数βとモデルパラメータａを用いてモデルパラメータｂ1 ，ｂ2 を次式により算出する。
ｂ1 ＝１−β
ｂ2 ＝１−ｂ1 −ａ
【００４６】
［特性多項式の係数演算］
図８の特性多項式の係数演算プログラムは、各気筒の噴射タイミングに同期して起動され、制御モデルのむだ時間ｄ分の根を０とする極配置法に基づいて特性多項式の係数Ａ１，Ａ２を次のようにして算出する。尚、極配置法に関しては、本出願人が先に出願した特願２０００−１８９７３４号の明細書に詳細に記載されている。
【００４７】
本プログラムが起動されると、まず、ステップ４４１で、減衰係数ζ、固有角周波数ω、噴射間隔ｄｔを読み込み、次のステップ４４２で、ω・ｄｔを上限ガード値（例えば０．６２８３）でガード処理する。つまり、ω・ｄｔが上限ガード値よりも大きい場合は、ω・ｄｔ＝上限ガード値にセットし、ω・ｄｔが上限ガード値以下の場合は、その時のω・ｄｔの値をそのまま用いる。このように、ω・ｄｔを上限ガード値でガード処理する理由は、ω・ｄｔの値が大きくなり過ぎると、制御精度が低下するためである。
【００４８】
ω・ｄｔのガード処理後、ステップ４４３に進み、特性多項式の係数Ａ１，Ａ２の算出に用いる変数ｅｚｗｄｔを次式により算出する。
ｅｚｗｄｔ＝ｅｘｐ（−ζ・ω・ｄｔ）
【００４９】
この変数ｅｚｗｄｔを算出する場合も、ＣＰＵの演算負荷を軽減するために、ζ・ω・ｄｔが例えば０．３５以下の時は、ｅｘｐ（−ζ・ω・ｄｔ）を次式により近似し、この近似式により変数ｅｚｗｄｔを算出する。
ｅｚｗｄｔ＝１−ζ・ω・ｄｔ＋０．５（ζ・ω・ｄｔ）2 
【００５０】
この近似式は、ζ・ω・ｄｔが大きくなるに従って、演算誤差が大きくなるため、ζ・ω・ｄｔが例えば０．３５よりも大きい領域では、予め、ζ・ω・ｄｔと変数ｅｚｗｄｔとの関係をテーブル化してＲＯＭに記憶しておき、このテーブルを検索して現在のζ・ω・ｄｔに応じた変数ｅｚｗｄｔを求める。尚、ζ・ω・ｄｔが０．３５以下の時にも、予め設定したテーブルから変数ｅｚｗｄｔを求めるようにしても良い。
【００５１】
この後、ステップ４４４に進み、特性多項式の係数Ａ１，Ａ２の算出に用いる他の変数ｃｏｓｚｗｔを次式により算出する。
ｃｏｓｚｗｔ＝ｃｏｓ｛（１−ζ2 ）0.5 ・ω・ｄｔ｝
【００５２】
この変数ｃｏｓｚｗｔを算出する場合も、ＣＰＵの演算負荷を軽減するために次の近似式を用いる。
ｃｏｓｚｗｔ＝１−０．５（１−ζ2 ）（ω・ｄｔ）2 
【００５３】
この後、ステップ４４５に進み、変数ｅｚｗｄｔ，ｃｏｓｚｗｔを用いて、特性多項式の係数Ａ１，Ａ２を次式により算出する。
Ａ１＝−２・ｅｚｗｄｔ・ｃｏｓｚｗｔ
Ａ２＝（ｅｚｗｄｔ）2 
【００５４】
［制御パラメータ演算］
図９の制御パラメータ演算プログラムは、特許請求の範囲でいう制御パラメータ演算手段としての役割を果たし、各気筒の噴射タイミングに同期して起動され、状態フィードバックの制御パラメータＦ０〜Ｆ８（ｄ＝６の場合）を次のようにして算出する。まず、ステップ４５１で、制御モデルのモデルパラメータａ，ｂ1 ，ｂ2 を読み込み、次のステップ４５２で、特性多項式の係数Ａ１，Ａ２を読み込む。
【００５５】
この後、ステップ４５３に進み、モデルパラメータａ，ｂ1 ，ｂ2 と係数Ａ１，Ａ２を用いて、次式により制御パラメータＦ０〜Ｆ８（ｄ＝６の場合）を順番に算出する。
Ｆ０＝（１＋Ａ１＋Ａ２）／（ｂ1 ＋ｂ2 ）
Ｆ２＝−１−ａ−Ａ１
Ｆ３＝ａ−Ａ２＋（１＋ａ）・Ｆ２
Ｆ４＝（１＋ａ）・Ｆ３−ａ・Ｆ２
Ｆ５＝（１＋ａ）・Ｆ４−ａ・Ｆ３
Ｆ６＝（１＋ａ）・Ｆ５−ａ・Ｆ４
Ｆ７＝（１＋ａ）・Ｆ６−ａ・Ｆ５
Ｆ１＝ａ／（ａ・ｂ1 ＋ｂ2 ）・（ａ・Ｆ７−ｂ1 ・Ｆ０）
Ｆ８＝ｂ2 ／ａ・Ｆ１
【００５６】
［ＦＡＦ演算］
図１０のＦＡＦ演算プログラムは、前述した図３の燃料噴射量演算プログラムのステップ３０６で起動され、次のようにして空燃比補正係数ＦＡＦを算出する。まず、ステップ４６１で、エンジン始動後の演算回数ｋをカウントするカウンタをクリアし、次のステップ４６２で、現在の燃料過剰率φ(i) 、目標燃料過剰率φref 、制御パラメータＦ０〜Ｆ８（ｄ＝６の場合）を読み込む。この後、ステップ４６３に進み、目標燃料過剰率φref と実際の燃料過剰率φ(i) との偏差ｅ(i) を算出する。
ｅ(i) ＝φref −φ(i)
【００５７】
この後、ステップ４６４に進み、演算回数ｋ＝０であるか否か（エンジン始動直後の最初の演算タイミングであるか否か）を判定し、演算回数ｋ＝０（エンジン始動直後の最初の演算タイミング）である場合は、ステップ４６５に進み、初期化処理を実行し、前回の燃料過剰率φ(i-1) ＝今回の燃料過剰率φ(i) にセットすると共に、過去の空燃比補正係数補正値ΔＦＡＦ(i-7) 〜ΔＦＡＦ(i-1) の記憶値を全て０にクリアする。この後、ステップ４６６に進み、前回から今回までの燃料過剰率の変化量Δφ(i) を算出する。
Δφ(i) ＝φ(i) −φ(i-1)
【００５８】
一方、エンジン始動後の２回目以降の演算タイミング（ｋ≠０）の場合は、上記ステップ４６５の初期化処理を行わずに、ステップ４６６に進み、燃料過剰率変化量Δφ(i) を算出する。
【００５９】
この後、ステップ４６７に進み、今回の空燃比補正係数補正値ΔＦＡＦ(i) を次式により算出する。
ΔＦＡＦ(i) ＝Ｆ１・Δφ(i) ＋Ｆ２・ΔＦＡＦ(i-1) ＋…
…＋Ｆ７・ΔＦＡＦ(i-6) ＋Ｆ８・ΔＦＡＦ(i-5) ＋Ｆ０・ｅ(i)
この場合、上記ステップ４６３〜４６７の処理が特許請求の範囲でいう補正値演算手段としての役割を果たす。
【００６０】
この後、ステップ４６８に進み、前回の空燃比補正係数ＦＡＦ(i-1) に今回の空燃比補正係数補正値ΔＦＡＦ(i) を加算して今回の空燃比補正係数ＦＡＦ(i) を求める。
ＦＡＦ(i) ＝ＦＡＦ(i-1) ＋ΔＦＡＦ(i)
このステップ４６８の処理が特許請求の範囲でいう操作量演算手段としての役割を果たす。
【００６１】
この後、ステップ４６９に進み、次回の空燃比補正係数ＦＡＦの演算に備えてφ(i-1) 、ΔＦＡＦ(i-7) 〜ΔＦＡＦ(i-1) の記憶データを更新する。
φ(i-1) ＝φ(i)
ＦＡＦ(i-7) ＝ΔＦＡＦ(i-6)
ＦＡＦ(i-6) ＝ΔＦＡＦ(i-5)
ＦＡＦ(i-5) ＝ΔＦＡＦ(i-4)
ＦＡＦ(i-4) ＝ΔＦＡＦ(i-3)
ＦＡＦ(i-3) ＝ΔＦＡＦ(i-2)
ＦＡＦ(i-2) ＝ΔＦＡＦ(i-1)
ＦＡＦ(i-1) ＝ΔＦＡＦ(i)
このステップ４６９と前記ステップ４６２の処理が特許請求の範囲でいう状態変数量出力手段としての役割を果たす。
【００６２】
この後、ステップ４７０に進み、エンジン始動後の演算回数ｋをカウントするカウンタをカウントアップしてステップ４６２に戻る。以後、ステップ４６２〜４７０の処理を噴射間隔で繰り返し、各気筒の噴射タイミングに同期して空燃比補正係数ＦＡＦを算出する。
【００６３】
以上説明した本実施形態（１）の効果を、図１１乃至図１３を用いて従来仕様と比較して説明する。従来仕様は、制御パラメータＦ０〜Ｆd+1 を用いて、次式により空燃比補正係数ＦＡＦ(i) を算出する。
0003758134.tif 000004
【００６４】
図１１は、制御パラメータを切り換えた時の空燃比補正係数ＦＡＦと燃料過剰率φの挙動を示すタイムチャートである。従来仕様では、制御パラメータを切り換えると、その瞬間に、空燃比補正係数ＦＡＦが一時的に乱れ、その結果、燃料過剰率φ（空燃比）も一時的に乱れる。
【００６５】
これに対し、本実施形態（１）では、制御パラメータＦ０〜Ｆ８と、燃料過剰率変化量Δφ(i) と、目標燃料過剰率φref と実際の燃料過剰率φ(i) との偏差ｅ(i) と、過去の空燃比補正係数補正値ΔＦＡＦ(i-1) 〜ΔＦＡＦ(i-7) とに基づいて今回の空燃比補正係数補正値ΔＦＡＦ(i) を演算し、今回の空燃比補正係数補正値ΔＦＡＦ(i) を前回の空燃比補正係数ＦＡＦ(i-1) に加算して今回の空燃比補正係数ＦＡＦ(i) を求める。このため、制御パラメータＦ０〜Ｆ８を運転条件等に応じて切り換えても、空燃比補正係数ＦＡＦと燃料過剰率φ（空燃比）が乱れることがなくなる。これにより、制御パラメータＦ０〜Ｆ８を運転条件等に応じて切り換えながら、安定した空燃比制御を実施することが可能となる。
【００６６】
図１２は、外乱が発生した時の空燃比補正係数ＦＡＦの上限ガード値からの復帰性能を示すタイムチャートである。従来仕様では、外乱が発生して、空燃比補正係数ＦＡＦが上限ガード値に張り付いた状態になると、空燃比補正係数ＦＡＦの上限ガード値が無い場合に空燃比補正係数ＦＡＦが上限ガード値よりも低下するまで、空燃比補正係数ＦＡＦが上限ガード値に張り付いた状態に維持される。このため、空燃比補正係数ＦＡＦが１．０に復帰するのが遅れ、外乱による燃料過剰率φ（空燃比）の乱れが収まるのが遅れる傾向がある。
【００６７】
これに対し、本実施形態（１）では、空燃比補正係数ＦＡＦが上限ガード値に張り付いた状態に維持される時間が従来仕様よりも短くなり、空燃比補正係数ＦＡＦが上限ガード値から従来仕様よりも早期に下がり始めるようになる。その結果、空燃比補正係数ＦＡＦが従来仕様よりも早期に１．０に復帰するようになり、外乱による燃料過剰率φ（空燃比）の乱れが収まるのが早くなる。
【００６８】
図１３は、モデル時定数、むだ時間、制御パラメータを変化させながら、外乱を入力した時の空燃比補正係数ＦＡＦと燃料過剰率φ（空燃比）の挙動を示すタイムチャートである。従来仕様では、制御パラメータが固定されているため、外乱入力時の空燃比補正係数ＦＡＦと燃料過剰率φ（空燃比）の乱れが大きくなる傾向がある。
【００６９】
これに対し、本実施形態（１）では、制御パラメータが運転条件に応じて変化するため、外乱入力時の空燃比補正係数ＦＡＦと燃料過剰率φ（空燃比）の乱れが従来仕様よりも小さくなり、外乱に対しても比較的安定した空燃比制御が可能となる。
【００７０】
《実施形態（２）》
上記実施形態（１）は、本発明を空燃比フィードバック制御系に適用した実施形態であるが、内燃機関の制御対象をフィードバック制御する制御系であれば、本発明を適用できる。図１４乃至図２２は、本発明をアイドル回転速度制御系に適用した実施形態（２）である。以下、本実施形態（２）で実行する各プログラムの処理内容を説明する。本実施例（２）では、ＩＳＣフィードバック補正量ＤＦＢが特許請求の範囲でいう操作量に相当し、ＩＳＣフィードバック補正量補正値ΔＤＦＢが特許請求の範囲でいう操作量補正値に相当し、前回から今回までのエンジン回転速度変化量ΔＮｅが特許請求の範囲でいう状態変数量の差分値に相当する。
【００７１】
［ＩＳＣＶ開度演算］
図１４のＩＳＣＶ開度演算プログラムは、所定時間毎又は所定クランク角毎に起動され、次のようにしてＩＳＣＶ開度ＤＯＰを算出する。ＩＳＣＶ開度ＤＯＰは、アイドル回転速度制御バルブを備えたシステムでは、アイドル回転速度制御バルブ（ＩＳＣＶ）の開度であり、スロットルバルブの開度でアイドル回転速度を制御する電子スロットルシステムでは、アイドル運転時のスロットル開度がＩＳＣＶ開度ＤＯＰとなる。
【００７２】
本プログラムが起動されると、まず、ステップ５０１で、現在のエンジン運転状態に応じてマップ等により基本開度Ｄbaseを算出する。この後、ステップ５０２で、基本開度Ｄbaseに対する各種の補正量ＤＡＬＬ（例えば冷却水温による補正量等）を算出し、次のステップ５０３で、アイドル回転速度制御（ＩＳＣ）のフィードバック条件が成立しているか否かを判定する。もし、ＩＳＣフィードバック条件が成立していなければ、ＩＳＣフィードバック補正量ＤＦＢを「０」にセットする。
【００７３】
一方、ＩＳＣフィードバック条件が成立していれば、ステップ５０５に進み、冷却水温ＴＨＷ、エアコンのＯＮ／ＯＦＦ信号、トルクコンバータ負荷信号等に応じてマップ等により目標アイドル回転速度を設定し、次のステップ６０６で、後述する図２３のＩＳＣフィードバック補正量演算プログラムを実行し、ＩＳＣフィードバック補正量ＤＦＢを算出する。
【００７４】
以上のようにして、ステップ５０４又は５０６でＩＳＣフィードバック補正量ＤＦＢを設定した後、ステップ５０７に進み、基本開度Ｄbase、各種補正量ＤＡＬＬ、ＩＳＣフィードバック補正量ＤＦＢを加算してＩＳＣＶ開度ＤＯＰを求める。
ＤＯＰ＝Ｄbase＋ＤＡＬＬ＋ＤＦＢ
【００７５】
［制御対象特性値演算］
図１５の制御対象特性値演算プログラムは、所定時間毎又は所定クランク角毎に起動され、制御対象の特性値であるモデルパラメータａ1 ，ａ2 ，ｂ1 ，ｂ2 を次のようにして算出する。まず、ステップ６０１で、冷却水温ＴＨＷを読み込み、次のステップ６０２で、冷却水温ＴＨＷに応じてマップ等によりモデルパラメータａ1 ，ａ2 ，ｂ1 ，ｂ2 を算出する。
【００７６】
ここで、モデルパラメータａ1 ，ａ2 ，ｂ1 ，ｂ2 をそれぞれ冷却水温ＴＨＷのみに基づいて算出する理由は、アイドル回転速度制御の実行中は、エンジン回転速度等の運転条件の変動が少なく、冷却水温ＴＨＷ等に応じて目標アイドル回転速度を設定するためである。尚、モデルパラメータａ1 ，ａ2 ，ｂ1 ，ｂ2 は、冷却水温ＴＨＷの他に、エアコンのＯＮ／ＯＦＦ信号、トルクコンバータ負荷信号等に応じて設定しても良い。
【００７７】
［減衰係数ζ、固有角周波数ω演算］
図１６の減衰係数ζ、固有角周波数ω演算プログラムは、所定時間毎又は所定クランク角毎に起動され、極配置法の演算に用いる減衰係数ζと固有角周波数ωを次のようにして算出する。まず、ステップ６１１で、冷却水温ＴＨＷを読み込み、次のステップ６１２で、減衰係数ζと固有角周波数ωを、それぞれ冷却水温ＴＨＷをパラメータとするマップにより算出する。
【００７８】
ここで、減衰係数ζと固有角周波数ωを、それぞれ冷却水温ＴＨＷのみに基づいて算出する理由は、アイドル回転速度制御の実行中は、エンジン回転速度等の運転条件の変動が少なく、冷却水温ＴＨＷ等に応じて目標アイドル回転速度を設定するためである。尚、減衰係数ζと固有角周波数ωは、冷却水温ＴＨＷの他に、エアコンのＯＮ／ＯＦＦ信号、トルクコンバータ負荷信号等に応じて設定しても良い。
【００７９】
尚、特性多項式の係数Ａ１，Ａ２は、前記実施形態（１）で説明した図８の特性多項式の係数演算プログラムと同じプログラムによって算出される。
【００８０】
［制御パラメータ演算］
図１７の制御パラメータ演算プログラムは、所定時間毎又は所定クランク角毎に起動され、状態フィードバックの制御パラメータＦ０〜Ｆ４（ｄ＝１の場合）を次のようにして算出する。まず、ステップ６２１で、制御モデルのモデルパラメータａ1 ，ａ2 ，ｂ1 ，ｂ2 を読み込み、次のステップ６２２で、特性多項式の係数Ａ１，Ａ２を読み込む。
【００８１】
この後、ステップ６２３に進み、モデルパラメータａ1 ，ａ2 ，ｂ1 ，ｂ2 と係数Ａ１，Ａ２を用いて、次式により制御パラメータＦ０〜Ｆ４（ｄ＝１の場合）を順番に算出する。
Ｆ０＝（１＋Ａ１＋Ａ２）／（ｂ1 ＋ｂ2 ）
Ｆ３＝−１−ａ1 −Ａ１
Ｆ１＝｛ａ2 （ａ2 ｂ1 −ａ1 ｂ2 ）
＋（−ａ1 ａ2 ｂ1 ＋ａ1 2 ｂ2 ＋ａ2 ｂ2 ）Ｆ３
＋ｂ2 （ａ2 ｂ1 −ａ1 ｂ2 ）Ｆ０｝
／（ａ1 ｂ1 ｂ2 −ａ2 ｂ1 2 ＋ｂ2 2 ）
Ｆ４＝−ａ2 −ｂ1 Ｆ１＋ａ1 Ｆ３−ｂ2 Ｆ０
Ｆ２＝（−ｂ2 Ｆ１＋ａ2 Ｆ３＋ａ1 Ｆ４）／ｂ1
【００８２】
［ＩＳＣフィードバック補正量演算］
図１８のＩＳＣフィードバック補正量演算プログラムは、前述した図１４のＩＳＣＶ開度演算プログラムのステップ５０６で起動され、次のようにしてＩＳＣフィードバック補正量ＤＦＢを算出する。まず、ステップ６３１で、エンジン始動後の演算回数ｋをカウントするカウンタをクリアし、次のステップ６３２で、現在のエンジン回転速度Ｎｅ、目標アイドル回転速度Ｎｔ、制御パラメータＦ０〜Ｆ４（ｄ＝１の場合）を読み込む。この後、ステップ６３３に進み、目標アイドル回転速度Ｎｔと現在のエンジン回転速度Ｎｅとの偏差ｅ(i) を算出する。
ｅ(i) ＝Ｎｔ−Ｎｅ
【００８３】
この後、ステップ６３４に進み、演算回数ｋ＝０であるか否か（エンジン始動直後の最初の演算タイミングであるか否か）を判定し、演算回数ｋ＝０（エンジン始動直後の最初の演算タイミング）である場合は、ステップ６３５に進み、初期化処理を実行し、前回のエンジン回転速度Ｎｅ(i-1) ＝今回のエンジン回転速度Ｎｅ(i) にセットし、前回のエンジン回転速度変化量ΔＮｅ(i-1) を０にセットすると共に、過去のＩＳＣフィードバック補正量補正値ΔＤＦＢ(i-2) ，ΔＤＦＢ(i-1) を０にセットする。
【００８４】
この後、ステップ６３６に進み、前回から今回までのエンジン回転速度変化量ΔＮｅ(i) を算出する。
ΔＮｅ(i) ＝Ｎｅ(i) −Ｎｅ(i-1)
【００８５】
この後、ステップ６３７に進み、今回のＩＳＣフィードバック補正量補正値ΔＤＦＢ(i) を次式により算出する。
0003758134.tif 000005
【００８６】
この後、ステップ４６８に進み、前回のＩＳＣフィードバック補正量ＤＦＢ(i-1) に今回のＩＳＣフィードバック補正量補正値ΔＤＦＢ(i) を加算して今回のＩＳＣフィードバック補正量ＤＦＢ(i）を求める。
ＤＦＢ(i) ＝ＤＦＢ(i-1) ＋ΔＤＦＢ(i)
【００８７】
この後、ステップ６３９に進み、次回のＩＳＣフィードバック補正量ＤＦＢの演算に備えて、ΔＮｅ(i-1) 、ΔＤＦＢ(i-2) 、ΔＤＦＢ(i-1) の記憶データを更新する。
ΔＮｅ(i-1) ＝ΔＮｅ(i)
ΔＤＦＢ(i-2) ＝ΔＤＦＢ(i-1)
ΔＤＦＢ(i-1) ＝ΔＤＦＢ(i)
【００８８】
この後、ステップ６４０に進み、エンジン始動後の演算回数ｋをカウントするカウンタをカウントアップしてステップ６３２に戻る。以後、ステップ６３２〜６４０の処理を所定時間毎又は所定クランク角毎に繰り返して、ＩＳＣフィードバック補正量ＤＦＢ(i）を算出する。
【００８９】
以上説明した本実施形態（２）では、制御パラメータＦ０〜Ｆ４と、エンジン回転速度変化量ΔＮｅ(i) と、目標アイドル回転速度Ｎｔと現在のエンジン回転速度Ｎｅとの偏差ｅ(i) と、過去のＩＳＣフィードバック補正量補正値ΔＤＦＢ(i-2) ，ΔＤＦＢ(i-1) とに基づいて今回のＩＳＣフィードバック補正量補正値ΔＤＦＢ(i) を演算し、今回のＩＳＣフィードバック補正量補正値ΔＤＦＢ(i）を前回のＩＳＣフィードバック補正量ＤＦＢ(i-1) に加算して今回のＩＳＣフィードバック補正量ＤＦＢ(i）を求める。このため、制御パラメータＦ０〜Ｆ４を運転条件等に応じて切り換えても、ＩＳＣフィードバック補正量ＤＦＢ(i）とエンジン回転速度Ｎｅが乱れることがなくなる。これにより、制御パラメータＦ０〜Ｆ４を運転条件等に応じて切り換えながら、安定したアイドル回転速度制御を行うことが可能となる。
【００９０】
《その他の実施形態》
前記実施形態（１），（２）では、制御パラメータを極配置法で算出するようにしたが、最適レギュレータで算出するようにしても良い。
【００９１】
また、モデルパラメータは、操作量（ＦＡＦ、ＩＳＣフィードバック補正量）と制御量（空燃比、エンジン回転速度）との関係から、例えば、システム同定手法を用いて、オンボードで算出しても良い。
【図面の簡単な説明】
【図１】本発明の実施形態（１）を示すエンジン制御システム全体の概略構成図
【図２】実施形態（１）の空燃比フィードバック制御系の各部の機能を表す機能ブロック図
【図３】実施形態（１）の燃料噴射量演算プログラムの処理の流れを示すフローチャート
【図４】実施形態（１）の制御対象特性値演算プログラムの処理の流れを示すフローチャート
【図５】実施形態（１）の噴射間隔演算プログラムの処理の流れを示すフローチャート
【図６】実施形態（１）の減衰係数ζ、固有角周波数ω演算プログラムの処理の流れを示すフローチャート
【図７】実施形態（１）のモデルパラメータ演算プログラムの処理の流れを示すフローチャート
【図８】実施形態（１）の特性多項式の係数演算プログラムの処理の流れを示すフローチャート
【図９】実施形態（１）の制御パラメータ演算プログラムの処理の流れを示すフローチャート
【図１０】実施形態（１）のＦＡＦ演算プログラムの処理の流れを示すフローチャート
【図１１】制御パラメータを切り換えた時の空燃比補正係数ＦＡＦと燃料過剰率φの挙動を示すタイムチャート
【図１２】外乱が発生した時の空燃比補正係数ＦＡＦの上限ガード値からの復帰性能を示すタイムチャート
【図１３】モデル時定数、むだ時間、制御パラメータを変化させながら、外乱を入力した時の空燃比補正係数ＦＡＦと燃料過剰率φ（空燃比）の挙動を示すタイムチャート
【図１４】実施形態（２）のＩＳＣＶ開度演算プログラムの処理の流れを示すフローチャート
【図１５】実施形態（２）の制御対象特性値演算プログラムの処理の流れを示すフローチャート
【図１６】実施形態（２）の減衰係数ζ、固有角周波数ω演算プログラムの処理の流れを示すフローチャート
【図１７】実施形態（２）の制御パラメータ演算プログラムの処理の流れを示すフローチャート
【図１８】実施形態（２）のＩＳＣフィードバック補正量演算プログラムの処理の流れを示すフローチャート
【符号の説明】
１１…エンジン（内燃機関）、１２…吸気管、１４…エアフローメータ（空気量検出手段）、２０…燃料噴射弁、２１…排気管、２２…触媒、２３…空燃比センサ（空燃比検出手段，状態検出手段）、２６…ＥＣＵ（制御パラメータ演算手段，状態変数量出力手段，操作量演算手段）。
【図１】
0003758134.tif 000006
【図２】
0003758134.tif 000007
【図３】
0003758134.tif 000008
【図４】
0003758134.tif 000009
【図５】
0003758134.tif 000010
【図６】
0003758134.tif 000011
【図７】
0003758134.tif 000012
【図８】
0003758134.tif 000013
【図９】
0003758134.tif 000014
【図１０】
0003758134.tif 000015
【図１１】
0003758134.tif 000016
【図１２】
0003758134.tif 000017
【図１３】
0003758134.tif 000018
【図１４】
0003758134.tif 000019
【図１５】
0003758134.tif 000020
【図１６】
0003758134.tif 000021
【図１７】
0003758134.tif 000022
【図１８】
0003758134.tif 000023