応答曲面モデル構築装置及びプログラム
【課題】エンジンの特性値と制御パラメータの関係を示す測定データに基づきその関係を精度良く再現する応答曲面モデルを効率的に構築する応答曲面モデル構築装置の提供。
【解決手段】応答曲面モデルに使用する制御パラメータの変数増加手段１３による追加処理と変数減少手段１４による削除処理をステップワイズ法による反復処理で行うにあたり、追加処理において、応答曲面モデルに取り込まれていない説明変数群の中から寄与率が最小の説明変数を選択してこれを追加候補説明変数とし、この追加候補説明変数と説明変数群の他の説明変数との間の相関係数を算出する。そして、この相関係数がすべて所定の相関係数閾値より小さい場合に、この追加候補説明変数を応答曲面モデルに取り込む。これにより、多重共線性を正確に考慮し、精度のよい応答曲面モデルの構築を行うことができる。
【解決手段】応答曲面モデルに使用する制御パラメータの変数増加手段１３による追加処理と変数減少手段１４による削除処理をステップワイズ法による反復処理で行うにあたり、追加処理において、応答曲面モデルに取り込まれていない説明変数群の中から寄与率が最小の説明変数を選択してこれを追加候補説明変数とし、この追加候補説明変数と説明変数群の他の説明変数との間の相関係数を算出する。そして、この相関係数がすべて所定の相関係数閾値より小さい場合に、この追加候補説明変数を応答曲面モデルに取り込む。これにより、多重共線性を正確に考慮し、精度のよい応答曲面モデルの構築を行うことができる。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ディーゼルエンジン等の制御パラメータの応答曲面モデルを構築する応答曲面モデル構築技術に関する。
【背景技術】
【０００２】
自動車のディーゼルエンジンは、コモンレール方式による多段噴射装置、ＥＧＲ（排気再循環）装置、過給器などの多数の制御装置を装備している。従って、これらの多数の制御装置のパラメータを運転条件、燃費、多数の排出ガス成分を考慮して最適に設定し制御するには、実験データから多数の制御パラメータと特性値との関係を精度良く表現する応答曲面モデルの構築が必要となる。
【０００３】
かかる応答曲面モデルの構築方法としては、従来から、実験データに基づいて多項式モデルによって近似する方法が広く用いられている（例えば、非特許文献１参照）。多項式モデルの作成においては、制御パラメータ等の目的変数に対して、その予測に有効な説明変数を特性値から選択する必要がある。この選択の方法としては、ステップワイズ法（変数増減法）等が提案されている（非特許文献２，３，特許文献１）。以下、特許文献１に記載のエンジン制御パラメータの適合方法について簡単に概説する。
【０００４】
まず、エンジンの電子制御スロットルの開度、点火時期、燃料噴射時期、ＥＧＲバルブの開度、給気バルブの開弁時期などの制御パラメータを（ｘ１，ｘ２，…，ｘ５）とおく。説明変数を、これら５つの制御パラメータ（ｘ１，ｘ２，…，ｘ５）、及びこれらの制御パラメータを自乗し又はこれらの制御パラメータの２つを掛け合わせて得られる１５個の２次項（ｘ１２，ｘ２２，…ｘ５２，ｘ１ｘ２，ｘ１ｘ３，…，ｘ４ｘ５）とし、これら２０個の説明変数（ｘ１，ｘ２，…，ｘ５，ｘ１２，ｘ２２，…ｘ５２，ｘ１ｘ２，ｘ１ｘ３，…，ｘ４ｘ５）を（ｘ１，ｘ２，…，ｘ２０）と記す。また、目的変数を、ＮＯｘ排出量や燃料消費量等のエンジンの各種特性値とし、これをｙとおく。説明変数ｘｉに対する偏回帰係数をａｉとし、重回帰式を次式で表す。
【０００５】
【数１】
【０００６】
各制御パラメータ（ｘ１，ｘ２，…，ｘ５）を振りながらｎ回の計測を行う。各説明変数及び目的変数の計測値をｘｉ（ｍ）（ｉ＝１，…，２０），ｙ（ｍ）（ｍ＝１，…，ｎ）とする。これらの計測値より、偏差積和Ｓｉｊ，Ｓｙｙを次式のように計算する。
【０００７】
【数２】
【０００８】
さらに、算出された偏差積和Ｓｉｊ，Ｓｙｙを用いて相関係数ｒｉｊを次式で算出する。
【０００９】
【数３】
【００１０】
以上の計算を行った後、まず、算出された相関係数ｒｉｊを用いて、初期の偏回帰係数ａｉ（ｉ＝１，…，５）を次式により算出する。
【００１１】
【数４】
【００１２】
次に、以下のような説明変数の追加処理及び説明変数の削除処理を繰り返すことによって、最適な説明変数の組み合わせによる重回帰式を求め、応答曲面モデルを構築する。
【００１３】
（説明変数の追加処理）
現在の重回帰式にｐ個（ｐ＜２０）の説明変数が含まれているとし、ｐ＋１個目の説明変数ｘｋを追加することを試みるものとする。まず、説明変数がｐ個の重回帰式の決定係数Ｒｐを次式により算出する。ここで、Ｙｐ（ｍ）はｍ番目の目的変数の計測値に対するｐ個の説明変数の重回帰式による予測値である。
【００１４】
【数５】
【００１５】
次に、（ｐ＋１）個目の説明変数ｘｋを追加した（ｐ＋１）個の重回帰式の決定係数Ｒｐ＋１を次式により算出する。
【００１６】
【数６】
【００１７】
これを、追加することが可能なすべての説明変数ｘｋについて行い、決定係数Ｒｐ＋１が最大となる説明変数ｘｋを、新たに追加する説明変数の候補として選出する。
【００１８】
次に、決定係数Ｒｐ＋１が最大となる説明変数ｘｋを追加した重回帰式において、追加した説明変数ｘｋの重回帰式に於ける寄与率を示すＦ値を次式より算出する。
【００１９】
【数７】
【００２０】
そして、このＦｐ＋１を所定の判定値Ｆｉｎと比較し、Ｆｐ＋１＞Ｆｉｎの場合には候補の説明変数ｘｋを追加して、次の説明変数の削除処理に移る。一方、Ｆｐ＋１≦Ｆｉｎの場合には、現在の重回帰式を応答曲面モデルとして出力して終了する。
【００２１】
（説明変数の削除処理）
現在の重回帰式にｐ＋１個（ｐ＋１＞５）の説明変数が含まれているとし、１つの説明変数ｘｋ（１５個の２次項ｘ６〜ｘ２０のうちの何れか）を削除することを試みるものとする。
【００２２】
まず、説明変数としてｘ１，ｘ２，…，ｘｋ−１，ｘｋ＋１，…，ｘｐ＋１のｐ個の変数に対して説明変数ｘｋを追加する場合のＦ値Ｆｐ＋１を式（７）により算出する。そして、このＦｐ＋１が所定の判定値Ｆｏｕｔ（Ｆｏｕｔ＜Ｆｉｎ）よりも小さい場合、この説明変数ｘｋを削除し、そうでなければ説明変数ｘｋを残す。この削除判定を、現在の重回帰式に含まれるすべての２次項の説明変数に対して行い、その後、先の説明変数の追加処理に移る。
【００２３】
以上のようにして、最適なエンジン制御パラメータを含む重回帰式により表される応答曲面モデルを構築する。
【先行技術文献】
【特許文献】
【００２４】
【特許文献１】特開２００４−６８７２９
【非特許文献】
【００２５】
【非特許文献１】山崎光悦，“応答曲面法の展開”，日本機械学会誌，日本国，社団法人日本機械学会，２００６年５月，第１０９巻，第１０５０号，pp.377-379．
【非特許文献２】管民郎，「多変量解析の実践（上）」，初版，上巻，日本国，株式会社現代数学社，1993年12月，pp.70-74．
【非特許文献３】河口至商，「多変量解析入門Ｉ」，初版，日本国，森北出版株式会社，1973年6月，pp.3-33．
【発明の概要】
【発明が解決しようとする課題】
【００２６】
しかしながら、上記従来のエンジン制御パラメータの適合方法においては、変数選択に於いて回帰係数のＦ検定による有意性の判断を行っているものの、厳密には回帰係数の検定のみでは多重共線性が考慮されているとはいえず、説明変数間の相関係数が高い説明変数の組み合わせが選択される場合がある。従って、実際のエンジンの多項式モデルの構築にあたっては、単純にステップワイズ法を各制御パラメータの説明変数の選択に直接適用するだけでは推定精度や実用上の観点から不十分である。
【００２７】
図１０は、実際のエンジンの多項式モデルの構築において、各制御パラメータの説明変数の選択にステップワイズ法を直接適用して重回帰モデルを求めた結果の一例である。図１０において、横軸は測定番号、縦軸はＮＯｘ値を表す。測定は７６回行い、これら７６回の測定から得られる７６組の制御パラメータ及びＮＯｘの測定値の組から、上記従来のステップワイズ法を用いて最適な重回帰式を求めた。重相関係数は０．９４９１であった。黒四角の点が実測値、白丸の点が重回帰式からの推定値である。図１０より、ステップワイズ法を直接適用した場合、ＮＯｘ値の変動が大きい点において誤差が大きく、特に、予測値が負値となる異常な予測点が複数出現しており、従来のエンジン制御パラメータの適合方法は、なお予測精度の改善が必要であることが分かる。そのためには、推定精度向上のために制御パラメータの非線形項を組み合わせた説明変数を複数作成し多項式モデルに導入することや、異常な予測値を出現させないために各説明変数を選択する際に多重共線性を厳密に考慮させることが必要である。また、１次、２次、及び高次の多項式モデルを生成し比較することにより、できる限り少ない説明変数で予測精度の高い多項式モデルを効率的に構築する方法が必要である。
【００２８】
そこで、本発明の目的は、エンジンの特性値と制御パラメータとの関係を示す測定データに基づいてその関係を精度良く再現する応答曲面モデルを効率的に構築する応答曲面モデル構築装置を提供することにある。
【課題を解決するための手段】
【００２９】
本発明に係る応答曲面モデル構築装置の第１の構成は、エンジンの各種制御パラメータである１次の説明変数及びそれに対するエンジンの特性値である目的変数の測定値の組である測定値セットを測定回数だけ格納する測定データ記憶手段と、
各測定値セットについて、各制御パラメータの累乗項を含むｋ個（ｋ＝２，…，Ｍ；Ｍは重回帰式に取り込む説明変数の次数の最大値）の制御パラメータの積項からなる高次の説明変数を算出し、前記測定データ記憶手段に保存する高次説明変数生成手段と、
前記各説明変数から前記目的変数を予測する重回帰式の偏回帰係数及び切片である重回帰係数を記憶する重回帰係数記憶手段と、
前記測定データ記憶手段に格納されたすべての説明変数の中からから重回帰式による目的変数の予測に有効な説明変数を選択するとともに、選択した各説明変数による重回帰式の重回帰係数を算出し前記重回帰係数記憶手段に格納する説明変数選択手段と、を備え、
前記説明変数選択手段は、
各説明変数について、当該説明変数が、重回帰式に取り込まれない予定の未追加変数群に属するのか、重回帰式に取り込む予定の追加変数群に属するのかを示す変数群属性を登録する変数群属性テーブルが格納された変数群属性記憶手段と、
前記変数群属性テーブルの各説明変数の変数群属性を、すべて未追加変数群に初期化する変数群属性初期化手段と、
変数群属性が未追加変数群である説明変数のうち、予測に有効な説明変数の変数群属性を追加変数群に変更する変数取込処理を行う変数増加手段と、
変数群属性が追加変数群である説明変数のうち、予測に有効でない説明変数の変数群属性を未追加変数群に変更する変数除去処理を行う変数減少手段と、
前記変数取込処理及び前記変数除去処理を交互に繰り返し実行し、前記変数群属性テーブルにおいて変数群属性が追加変数群の説明変数が未追加変数群の説明変数よりも少なくなるか、又は前記変数取込処理において何れの説明変数の変数群属性も追加変数群に変更されなかった場合に繰り返し処理を終了する変数選択処理制御手段と、を備え、
前記変数増加手段は、
変数群属性が未追加変数群の説明変数を順次１つずつ選択し、選択した説明変数，変数群属性が追加変数群の説明変数，及び目的変数の測定値セットを前記測定データ記憶手段からすべて読み出し、読み出した測定値セットから重回帰係数を算出し、読み出した測定値セット及び当該重回帰係数から、選択した説明変数の寄与率を算出する追加時寄与率演算手段と、
前記変数群属性テーブルにおいて変数群属性が未追加変数群の説明変数のうち、前記追加時寄与率演算手段により算出された寄与率が最大の説明変数を初期の追加候補説明変数に選択する最大寄与率選定手段と、
前記追加候補説明変数の寄与率が、重回帰式に取り込む取込基準値より大きいか否かを判定し、当該取込基準値より大きくない場合、前記変数取込処理を終了する追加基準判定手段と、
前記追加基準判定手段が、当該追加候補説明変数の寄与率が取込基準値より大きいと判定した場合、当該追加候補説明変数と、前記変数群属性テーブルにおいて変数群属性が未追加変数群の他の各説明変数との相関係数を順次算出する相関係数演算手段と、
前記相関係数演算手段が算出する相関係数がすべて所定の相関係数閾値より小さい場合、前記変数群属性テーブルにおける当該追加候補説明変数の変数群属性を追加変数群に変更し、前記変数取込処理を終了する追加判定手段と、
前記相関係数演算手段が算出する相関係数の何れかが前記相関係数閾値より大きい場合、変数群属性が未追加変数群の説明変数のうち現在の追加候補説明変数の次に寄与率が大きい次候補説明変数がなければ前記変数取込処理を終了し、次候補説明変数があるときは、現在の追加候補説明変数を当該次候補説明変数に変更する追加候補説明変数変更手段と、
前記追加基準判定手段、前記相関係数演算手段、前記追加判定手段、及び前記追加候補説明変数変更手段による処理を繰り返し実行する変数取込処理制御手段と、を備え、
前記変数減少手段は、
前記変数群属性テーブルにおける変数群属性が追加変数群の説明変数を順次１つずつ選択し、当該選択した説明変数を除く追加変数群の説明変数に、当該選択した説明変数を追加する場合の寄与率を算出する除去時寄与率演算手段と、
前記除去時寄与率演算手段により算出される寄与率が最小の説明変数を、除去候補説明変数に選択する最小寄与率選択手段と、
当該除去候補説明変数の寄与率が、重回帰式から除去する除去基準値より小さい場合、前記変数群属性テーブルにおける当該除去候補説明変数の変数群属性を未追加変数群に変更し、それ以外の場合には前記変数除去処理を終了する除去基準判定手段と、
前記最小寄与率選択手段及び前記除去基準判定手段による処理を、前記除去基準判定手段が前記変数除去処理を終了するまで反復実行する制御を行う変数除去処理制御手段と、を備えたことを特徴とする。
【００３０】
この構成によれば、変数増加手段における処理において、追加候補説明変数ｘαｑ＋１の重回帰式への取り込みを行う際に、未追加変数群Ｐの他の各説明変数ｘｊとの相関係数ｒαｑ＋１，ｊを調べ、これらすべての相関係数ｒαｑ＋１，ｊが相関係数閾値ｒｉｎより小さい場合にのみ追加候補説明変数ｘαｑ＋１の重回帰式への取り込みを行うことで、多重共線性を正確に考慮し、精度のよい応答曲面モデルの構築を行うことができる。
【００３１】
また、本発明に係る応答曲面モデル構築装置の第２の構成は、前記第１の構成において、前記寄与率は、Ｆ値又はＴ統計量であることを特徴とする。
【００３２】
また、本発明に係る応答曲面モデル構築装置の第３の構成は、前記第１又は２の構成において、前記説明変数選択手段は、高次の説明変数の最大次数ｈｋが１からｈｍａｘのそれぞれの場合について重回帰式の重回帰係数を算出し、前記重回帰係数記憶手段に格納するものであり、
前記重回帰係数記憶手段に格納された最大次数ｈｋ（ｋ＝１，…，ｈｍａｘ）の場合の重回帰式のそれぞれに対して、根平均自乗誤差ｒｋを算出し、算出した各根平均自乗誤差ｒｋの最大値ｒｍａｘを抽出するとともに、各根平均自乗誤差ｒｋ及びその最大値ｒｍａｘに基づいて、各最大次数ｈｋにおけるコストＪｋを算出し、コストＪｋが最大となる最大次数ｈｋに対する重回帰式を選択し出力する重回帰モデル評価手段を備えたことを特徴とする。
【００３３】
この構成により、重回帰モデル評価手段によって、最も適した最大次数ｈｋの重回帰式による重回帰モデルを選定でき、不必要に説明変数が増大することを防止することができる。
【００３４】
また、本発明に係るプログラムは、コンピュータに読み込んで実行させることにより、前記第１乃至３のいずれかの構成の応答曲面モデル構築装置として機能させることを特徴とする。
【発明の効果】
【００３５】
以上のように、本発明によれば、変数増加手段における処理において、追加候補説明変数ｘαｑ＋１の重回帰式への取り込みを行う際に、未追加変数群Ｐの他の各説明変数ｘｊとの相関係数ｒαｑ＋１，ｊを調べ、これらすべての相関係数ｒαｑ＋１，ｊが相関係数閾値ｒｉｎより小さい場合にのみ追加候補説明変数ｘαｑ＋１の重回帰式への取り込みを行うことで、多重共線性を正確に考慮し、精度のよい応答曲面モデルの構築を行うことができる。
【００３６】
また、さらに重回帰モデル評価手段によって、最も適した最大次数ｈｋの重回帰式による重回帰モデルを選定すれば、不必要に説明変数が増大することを防止することができる。
【図面の簡単な説明】
【００３７】
【図１】本発明の実施例１に係る応答曲面モデル構築装置の構成を表すブロック図である。
【図２】変数増加手段１３の構成を表すブロック図である。
【図３】変数減少手段１４の構成を表すブロック図である。
【図４】本実施例に係る応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理を表すフローチャートである。
【図５】本実施例に係る応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理を表すフローチャートである。
【図６】実際のエンジンの多項式モデルの構築において、本実施例の応答曲面モデル構築装置１により重回帰モデルを求めた結果の一例である。
【図７】本発明の実施例２に係る応答曲面モデル構築装置１の構成を表すブロック図である。
【図８】実施例２の応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理を表すフローチャートである。
【図９】実施例２の応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理を表すフローチャートである。
【図１０】実際のエンジンの多項式モデルの構築において、各制御パラメータの説明変数の選択にステップワイズ法を直接適用して重回帰モデルを求めた結果の一例である。
【発明を実施するための形態】
【００３８】
以下、本発明を実施するための形態について、図面を参照しながら説明する。
【実施例１】
【００３９】
（１）応答曲面モデル構築装置の構成
（１−１）全体構成
図１は、本発明の実施例１に係る応答曲面モデル構築装置の構成を表すブロック図である。
【００４０】
図１において、応答曲面モデル構築装置１は、エンジン１０１のＮＯｘ排出量や燃料消費量等の各種特性値を、所定の目標値の範囲内となるように、エンジン１０１の電子制御スロットルの開度、点火時期、燃料噴射時期、ＥＧＲバルブの開度、吸気バルブの開弁時期などの各種制御パラメータを適合させる応答曲面モデルの構築を行う装置である。すなわち、応答曲面モデル構築装置１は、エンジン１０１の各種制御パラメータを説明変数、エンジン１０１の種特性値を目標変数とする重回帰式の算出を行う。ここで、応答曲面モデルの重回帰式の算出を行うため、各種制御パラメータを人為的に変動させたときの各種特性値の測定データを予め用意しておく必要がある。そこで、制御パラメータ設定手段１０２は、エンジン１０１に対して予め定められた各種制御パラメータの組み合わせに基づいて、順次各種制御パラメータを設定する。そして、制御パラメータ測定手段１０３及び特性値測定手段１０４は、エンジン１０１の各種制御パラメータの設定が変更される毎に、エンジン１０１の各種制御パラメータ及び特性値を測定し、測定データ記憶手段２に格納する。
【００４１】
応答曲面モデル構築装置１は、測定データ記憶手段２、高次説明変数生成手段３、重回帰係数記憶手段４、説明変数選択手段５、及び重回帰モデル出力手段７を備えている。尚、この応答曲面モデル構築装置１は、マイコンやＦＰＧＡなどによりハードウェア的に構成することもできるが、プログラムとして用意してパソコンなどの汎用コンピュータにこのプログラムをロードすることによってソフトウェア的に構成することもできる。
【００４２】
測定データ記憶手段２は、上述したように、制御パラメータ測定手段１０３及び特性値測定手段１０４により測定される、エンジン１０１の各種制御パラメータである１次の説明変数（ｘ１，ｘ２，…，ｘＮ）及びそれに対するエンジン１０１の特性値である目的変数ｙの測定値の組である測定値セット｛（ｘ１，ｘ２，…，ｘＮ），ｙ｝を測定回数Γ分だけ格納するためのものである。ここで、Ｎは制御パラメータの数である。
【００４３】
なお、測定データ記憶手段２において、各測定回に対する測定値セット｛（ｘ１，ｘ２，…，ｘＮ），ｙ｝は、表１に示すような標本値テーブルとして格納される。
【００４４】
【表１】
【００４５】
表１において、「測定番号」のフィールドは、各測定回の番号が格納されるフィールドである。「制御パラメータ」のｘ１〜ｘＮまでのフィールドは、各測定回の１次の説明変数が格納されるフィールドである。また、ｘＮ＋１〜ｘＮ＋Ｌまでのフィールドは、各測定回の高次の説明変数が格納されるフィールドである。「特性値」のフィールドは、各測定回の特性値ｙが格納されるフィールドである。ここで、「高次の説明変数」とは、各制御パラメータの累乗項（ｘｉ２，ｘｉ３等）を含むｋ個（ｋ＝２，…，Ｍ；Ｍは重回帰式に取り込む説明変数の次数の最大値）の制御パラメータの積項からなる説明変数をいう。Ｌは高次の説明変数の数を表す。重回帰式に取り込む説明変数の次数の最大値をＭとすると、Ｌ＝ＮＣＭ＝Ｎ！／（Ｍ！・（Ｎ−Ｍ）！）となる。以下、ｎ番目の測定回に於ける説明変数（ｘ１，ｘ２，…，ｘＮ＋Ｌ）及び目的変数ｙを、それぞれ（ｘ１（ｎ），ｘ２（ｎ），…，ｘＮ＋Ｌ（ｎ））及びｙ（ｎ）と記す。また、高次の説明変数については、どの説明変数を掛け合わせた積項かを明示する必要がある場合、２次の説明変数ｘｉｘｊをｘｉｊ、３次の説明変数ｘｉｘｊｘｋをｘｉｊｋのように記すこととし、一般に、ｎ次の説明変数ｘα１…ｘαｎをｘα１…αｎのように記す。
【００４６】
高次説明変数生成手段３は、測定データ記憶手段２に格納された各測定値セット｛（ｘ１（ｎ），ｘ２（ｎ），…，ｘＮ（ｎ）），ｙ（ｎ）｝（ｎ＝１，…，Γ）について、それぞれの高次の説明変数（ｘＮ＋１（ｎ），ｘＮ＋２（ｎ），…，ｘＮ＋Ｌ（ｎ））を算出し、前記測定データ記憶手段２に保存する。
【００４７】
重回帰係数記憶手段４は、各説明変数ｘα１…αｎ（ｎ＝１，…，Ｍ）から前記目的変数ｙを予測する重回帰式の重回帰係数を記憶する。ここで、「重回帰係数」とは、重回帰式の回帰係数βα１…αｎ及び切片β０をいう。尚、偏回帰係数βα１…αｎは、重回帰式における説明変数ｘα１…αｎの偏回帰係数を表す。
【００４８】
説明変数選択手段５は、測定データ記憶手段２に格納されたすべての説明変数の中からから重回帰式による目的変数の予測に有効な説明変数を選択するとともに、選択した各説明変数による重回帰式の重回帰係数を算出し重回帰係数記憶手段４に格納する。この説明変数選択手段５の構成の詳細については、後述する。
【００４９】
重回帰モデル出力手段７は、最終的に選択された説明変数名と、算出された各重回帰係数を出力する。
【００５０】
（１−２）説明変数選択手段５の構成
説明変数選択手段５は、図１に示すように、変数群属性記憶手段１１、変数群属性初期化手段１２、変数増加手段１３、変数減少手段１４、及び変数選択処理制御手段１５を備えている。
【００５１】
変数群属性記憶手段１１は、上記Ｎ＋Ｌ個の各説明変数について、変数群属性を登録する変数群属性テーブルを記憶する。ここで、「変数群属性」とは、説明変数が、未追加変数群Ｐに属するのか、追加済変数群Ｑに属するのかを示す属性をいう。また、「未追加変数群」とは、重回帰式に取り込まれない予定の変数の集合をいい、「追加済変数群」とは、重回帰式に取り込む予定の変数の集合をいう。表２に、変数群属性テーブルを示す。表２において、ｇｉ（ｉ＝１，２，…，Ｎ＋Ｌ）は、目的変数ｘｉに対する変数群属性であり、「未追加変数群」又は「追加済変数群」の何れかの値をとる。
【００５２】
【表２】
【００５３】
変数群属性初期化手段１２は、変数群属性記憶手段１１に格納された変数群属性テーブルの各説明変数の変数群属性（ｇ１，…ｇＮ＋Ｌ）を、すべて未追加変数群Ｐに初期化する。
【００５４】
変数増加手段１３は、変数群属性が未追加変数群Ｐである説明変数のうち、予測に有効な説明変数の変数群属性を追加済変数群Ｑに変更する変数取込処理を行う。この変数増加手段１３の構成の詳細については、後述する。
【００５５】
変数減少手段１４は、変数群属性が追加済変数群Ｑである説明変数のうち、予測に有効でない説明変数の変数群属性を未追加変数群Ｐに変更する変数除去処理を行う。この変数減少手段１の構成の詳細については、後述する。
【００５６】
変数選択処理制御手段１５は、変数群属性初期化手段１２、変数増加手段１３及び変数減少手段１４を制御して、最初に変数群属性テーブルを初期化した後、変数取込処理及び前記変数除去処理を交互に繰り返し実行し、変数群属性テーブルにおいて変数群属性が追加済変数群Ｑの説明変数が未追加変数群Ｐの説明変数よりも少なくなるか、又は前記変数取込処理において何れの説明変数の変数群属性も追加済変数群Ｑに変更されなかった場合に繰り返し処理を終了する処理を行う。
【００５７】
（１−３）変数増加手段１３の構成
図２は、変数増加手段１３の構成を表すブロック図である。図２に示すように、変数増加手段１３は、追加時寄与率記憶手段２０、追加時寄与率演算手段２１、最大寄与率選定手段２２、追加基準判定手段２３、相関係数演算手段３４、追加判定手段２５、追加候補説明変数変更手段２６、及び変数取込処理制御手段２７を備えている。
【００５８】
追加時寄与率演算手段２１は、変数群属性テーブルを参照しながら、変数群属性が未追加変数群Ｐの説明変数ｘαｑ＋１（ｋ）（ｋ＝１，…，Ｎ＋Ｌ−ｑ）を順次１つずつ選択する。ここで、ｑは変数群属性が追加済変数群Ｑである説明変数の現在の個数である。そして、選択した説明変数（以下「選択説明変数」という。）ｘαｑ＋１（ｋ），現在の変数群属性が追加済変数群Ｑの説明変数（ｘα１，ｘα２，…，ｘαｑ），及び目的変数ｙの測定値セット｛（ｘα１（ｎ），ｘα２（ｎ），…，ｘαｑ（ｎ），ｘαｑ＋１（ｋ）（ｎ）），ｙ（ｎ）｝（ｎ＝１，…，Ｐ）を測定データ記憶手段２からＰ組分すべて読み出し、読み出した測定値セットから重回帰係数Ｆ（βαｑ＋１（ｋ））を算出し、重回帰係数記憶手段４に、表３のような一時回帰係数テーブルとして格納する。ここで、β０（ｋ）は重回帰式の切片、βαｉ（ｋ）（ｉ＝１，…，ｑ）は説明変数ｘαｋに対する偏回帰係数、βαｑ＋１（ｋ）は説明変数ｘαｑ＋１（ｋ）に対する偏回帰係数である。さらに、追加時寄与率演算手段２１は、読み出した測定値セット及び算出した重回帰係数から、選択説明変数ｘαｑ＋１（ｋ）の寄与率Ｆ（βαｑ＋１（ｋ））を算出し、追加時寄与率記憶手段２０に、表４のような寄与率テーブルとして格納する。ここで、Ｆ（βαｑ＋１（ｋ））は、選択説明変数ｘαｑ＋１（ｋ）を追加済変数群Ｑに追加した場合の当該選択説明変数ｘαｑ＋１（ｋ）の寄与率である。
【００５９】
【表３】
【００６０】
【表４】
【００６１】
最大寄与率選定手段２２は、追加時寄与率記憶手段２０から、追加時寄与率演算手段２１により算出された寄与率を読み出し、それらのうち、最大の寄与率に対応する選択説明変数ｘαｑ＋１（ｋ）を、初期の追加候補説明変数ｘαｑ＋１とする。
【００６２】
追加基準判定手段２３は、追加候補説明変数ｘαｑ＋１の寄与率Ｆ（βαｑ＋１）が、重回帰式に取り込む取込基準値Ｆｉｎより大きいか否かを判定し、当該取込基準値Ｆｉｎより大きくない場合には、変数増加手段１３における変数取込処理を終了する。
【００６３】
相関係数演算手段２４は、追加基準判定手段２３が、追加候補説明変数ｘαｑ＋１の寄与率Ｆ（βαｑ＋１）が取込基準値Ｆｉｎより大きいと判定した場合、当該追加候補説明変数ｘαｑ＋１と、変数群属性テーブルにおいて変数群属性が未追加変数群Ｐの他の各説明変数との相関係数を順次算出する。
【００６４】
追加判定手段２５は、相関係数演算手段２４が算出する相関係数がすべて所定の相関係数閾値ｒｉｎより小さい場合、変数群属性テーブルにおける当該追加候補説明変数ｘαｑ＋１の変数群属性を追加済変数群Ｑに変更し、変数増加手段１３における変数取込処理を終了する。
【００６５】
追加候補説明変数変更手段２６は、相関係数演算手段３４が算出する相関係数の何れかが相関係数閾値ｒｉｎ以上の場合、追加時寄与率記憶手段２０を参照し、変数群属性が未追加変数群Ｐの説明変数のうち現在の追加候補説明変数の次に寄与率の大きい説明変数（以下「次候補説明変数」という。）があるか否かを検索する。そして、次候補説明変数がなければ、変数増加手段１３における変数取込処理を終了する。また、次候補説明変数があるときは、現在の追加候補説明変数ｘαｑ＋１を当該次候補説明変数に変更する。
【００６６】
変数取込処理制御手段２７は、変数増加手段１３における変数取込処理の全体制御を行う。すなわち、追加時寄与率演算手段２１及び最大寄与率選定手段２２による処理を最初に行った後、追加基準判定手段２３、相関係数演算手段２４、追加判定手段２５、及び追加候補説明変数変更手段２６による処理を繰り返し実行する制御を行う。
【００６７】
（１−４）変数減少手段１４の構成
図３は、変数減少手段１４の構成を表すブロック図である。変数減少手段１４は、図３に示すように、除去時寄与率記憶手段３０、除去時寄与率演算手段３１、最小寄与率選択手段３２、除去基準判定手段３３、及び変数除去処理制御手段３４を備えている。
【００６８】
除去時寄与率演算手段３１は、変数群属性記憶手段１１に記憶されている変数群属性テーブルを参照し、変数群属性が追加済変数群Ｑである説明変数ｘαｋを順次１つずつ選択し、当該選択した説明変数（以下「選択説明変数」という。）ｘαｋを除く追加済変数群Ｑの説明変数（ｘα１，…，ｘαｋ−１，ｘαｋ＋１，…，ｘαｑ）に、当該選択した説明変数ｘαｋを追加する場合の寄与率Ｆ（βαｋ）を算出し、除去時寄与率記憶手段３０に寄与率テーブルとして格納する。ここで、ｑは変数群属性が追加済変数群Ｑである説明変数の現在の個数である。また、寄与率テーブルは表４に示したものと同様である。
【００６９】
最小寄与率選択手段３２は、除去時寄与率記憶手段３０格納された寄与率テーブルを参照し、除去時寄与率演算手段３１により算出される寄与率が最小の説明変数を、除去候補説明変数ｘminに選択する。
【００７０】
除去基準判定手段３３は、最小寄与率選択手段３２により選択された除去候補説明変数ｘminの寄与率Ｆminが、重回帰式から除去する除去基準値Ｆｏｕｔより小さい場合、変数群属性記憶手段１１に格納された変数群属性テーブルにおける当該除去候補説明変数ｘminの変数群属性を未追加変数群Ｐに変更する。それ以外の場合には、変数減少手段１４における変数除去処理を終了する。
【００７１】
変数除去処理制御手段３４は、変数減少手段１４における変数除去処理の全体制御を行う。すなわち、除去時寄与率演算手段３１による寄与率テーブルの生成を実行した後、最小寄与率選択手段３２及び除去基準判定手段３３による処理を、除去基準判定手段３３が変数除去処理を終了させるまで反復実行する制御を行う。
【００７２】
（２）応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理
以上のように構成された本実施例に係る応答曲面モデル構築装置１について、次に、応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理の詳細について説明する。
【００７３】
図４，図５は、本実施例に係る応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理を表すフローチャートである。
【００７４】
まず、ステップＳ１において、高次説明変数生成手段３は、測定データ記憶手段２に格納された標本値テーブルから、Γ組の測定値セット（ｘ１（ｎ），ｘ２（ｎ），…，ｘＮ（ｎ））（ｎ＝１，…，Γ）を順次読み出し、それぞれ、Ｌ個の高次の説明変数を算出し標本値テーブルに格納する。
【００７５】
次に、ステップＳ２において、変数群属性初期化手段１２は、標本値テーブルに格納したＮ＋Ｌ個の説明変数について、変数群属性記憶手段１１に表２のような変数群属性テーブルを生成し、変数群属性テーブルの各説明変数ｘｉ（ｉ＝１，…，Ｎ＋Ｌ）の変数群属性ｇｉを未追加変数群Ｐに初期化し、未追加変数群Ｐの変数の個数ｐを０，追加済変数群Ｑの変数の個数ｑをＮ＋Ｌとする。
【００７６】
次に、変数取込処理制御手段２７は、以下の変数取込処理（Ｓ３〜Ｓ１５）と変数除去処理（Ｓ１６〜Ｓ２２）を反復実行し、多重共線性を考慮して有効な説明変数を選択するとともに、選択した説明変数による重回帰式を算出する。
【００７７】
（変数取込処理）
変数取込処理では、まず、ステップＳ３において、変数取込処理制御手段２７は、未追加変数群Ｐの変数の個数ｐと追加済変数群Ｑの変数の個数ｑを比較し、ｑ＜ｐの場合には、変数取込処理を打ち切って、ステップＳ２３に移行する。ｑ≧ｐの場合は、次のステップＳ４に進む。
【００７８】
次に、ステップＳ４において、追加時寄与率演算手段２１は、変数属性テーブルの変数属性が未追加変数群Ｐである説明変数ｘαｑ＋１（ｋ）（ｋ＝１，…，ｐ）を順次１つずつ選択する。ここで、ｐ＝Ｎ＋Ｌ−ｑである。
【００７９】
次に、ステップＳ５において、追加時寄与率演算手段２１は、追加済変数群Ｑの説明変数（ｘα１，ｘα２，…，ｘαｑ）と選択説明変数ｘαｑ＋１（ｋ）との重回帰式の重回帰係数βｉ（ｋ）（ｉ＝０，１，…，ｑ，ｑ＋１）を算出し、表３に示した一時回帰係数テーブルに格納する。ここで、選択説明変数ｘαｑ＋１（ｋ）に対する重回帰式は次式（８）のように表される。
【００８０】
【数８】
【００８１】
測定データ記憶手段２の標本値テーブルに格納された測定値セット｛（ｘα１（ｎ），ｘα２（ｎ），…，ｘαｑ（ｎ），ｘαｑ＋１（ｋ）（ｎ）），ｙ（ｎ）｝（ｎ＝１，…，Γ）から重回帰係数βｉ（ｋ）（ｉ＝０，１，…，ｑ，ｑ＋１）を算出する方法については、周知の線形回帰の計算方法が使用される（例えば、非特許文献３参照）。具体的には、目的変数ベクトルＹ、説明変数行列Ｘ、重回帰係数ベクトルＢを次式（９）のように置く。
【００８２】
【数９】
【００８３】
このとき、重回帰係数ベクトルＢは、次式（１０）により算出することができる。但し、Ｘｔは行列Ｘの転置行列を表す。
【００８４】
【数１０】
【００８５】
次に、ステップＳ６において、追加時寄与率演算手段２１は、選択説明変数ｘαｑ＋１（ｋ）に対する寄与率Ｆ（βαｑ＋１（ｋ））を算出し、追加時寄与率記憶手段２０に、表４に示した寄与率テーブルを格納する。具体的には、次のような演算が行われる。
【００８６】
まず、次式のような（ｑ＋１）個の説明変数（ｘα１，ｘα２，…，ｘαｑ，ｘαｑ＋１（ｋ））の分散共分散行列Σの逆行列Σ−１を算出する。ここで、Ｅ［ ］は平均を表す。
【００８７】
【数１１】
【００８８】
また、次式により目的変数ｙの予測誤差二乗和Ｓｅ及び回帰式の誤差分散Ｖｅを算出する。ここで、Γは標本値テーブルに格納された測定値セットの組数（測定回数）である。
【００８９】
【数１２】
【００９０】
そして、これらの算出値から、選択説明変数ｘαｑ＋１（ｋ）に対する寄与率Ｆ（βαｑ＋１（ｋ））を次式により算出する。ここで、Ｓ（βαｑ＋１（ｋ））２は不偏分散、Ｔ（βαｑ＋１（ｋ））はＴ値（Ｔ統計量）である。
【００９１】
【数１３】
【００９２】
ここでは寄与率としてＦ値を使用した。Ｆ値は補修団の重回帰係数βｊがβｊ＝０という仮説のもとで、自由度（１，Ｎ−ｑ−１）のＦ分布に従うとする。Ｆ値がＦ分布の限界値と比較してＦ（βｊ）＞Ｆ（１，Ｎ−ｑ−１，α）（右辺のＦ（ ）はＦ関数）ならば仮説を棄却し、係数βｊが有意であるとする。ここで、αは有意水準である。すなわち、ここではＦ値がＦｉｎの基準値よりも大きい偏回帰係数を仮採用し、それ以外は棄却する処理を行っている。
【００９３】
尚、寄与率としては上式（１３ｂ）のＴ値を使用してもよい。Ｔ値を使用する場合、母集団の偏回帰係数βｊがβｊ＝０という仮説のもとで自由度（Ｎ−ｑ−１）のｔ分布に従うとする。この場合、Ｔ（βｊ）＞ｔ（Ｎ−ｑ−１，α／２）ならば仮説を棄却し、係数βｊは有意であるとする。
【００９４】
次に、ステップＳ７において、追加時寄与率演算手段２１は、未追加変数群Ｐし属するｐ個の説明変数ｘαｑ＋１（ｋ）（ｋ＝１，…，ｐ）についてすべて寄与率Ｆ（βαｑ＋１（ｋ））を算出したか否かを判定し、すべての説明変数ｘαｑ＋１（ｋ）についての算出が終了していないときはステップＳ４に戻り、すべての説明変数ｘαｑ＋１（ｋ）についての算出が終了したときは、次のステップＳ８に移行する。
【００９５】
次に、ステップＳ８において、最大寄与率選定手段２２は、追加時寄与率記憶手段２０に格納された寄与率テーブルを参照し、寄与率Ｆ（βαｑ＋１（ｋ））が最大の説明変数ｘαｑ＋１（ｋ）を索出し、これを初期の追加候補説明変数ｘαｑ＋１に設定する。
【００９６】
次に、ステップＳ９において、追加基準判定手段２３は、追加候補説明変数ｘαｑ＋１の寄与率Ｆ（βαｑ＋１）が、重回帰式に取り込む取込基準値Ｆｉｎより大きいか否かを判定する。当該取込基準値Ｆｉｎより大きい場合には、次のステップＳ１０に移行し、当該取込基準値Ｆｉｎより大きくない場合には、変数取込処理を終了して、ステップＳ２３に移行する。
【００９７】
次に、ステップＳ１０において、相関係数演算手段２４は、変数群属性テーブルにおいて変数群属性が未追加変数群Ｐの他の各説明変数ｘｊを１つ選択し、追加候補説明変数ｘαｑ＋１と、未追加変数群Ｐの他の各説明変数ｘｊとの相関係数ｒαｑ＋１，ｊを順次算出する。
【００９８】
次に、ステップＳ１１において、追加判定手段２５は、相関係数演算手段２４が算出する相関係数ｒαｑ＋１，ｊが所定の相関係数閾値ｒｉｎより小さいか否かを判定する。
【００９９】
相関係数閾値ｒｉｎより小さい場合は、ステップＳ１３において、全ての未追加説明変数群Ｐの説明変数ｘｊが選択されたか否かを判定し、まだ、選択されていない説明変数ｘｊがあれば、ステップＳ１０に戻り、すべての説明変数ｘｊが選択されていれば、ステップＳ１３において、変数群属性テーブルにおける当該追加候補説明変数ｘαｑ＋１の変数群属性を追加済変数群Ｑに変更し、変数増加手段１３における変数取込処理を終了し、次の変数除去処理（ステップＳ１６）に移行する。
【０１００】
一方、相関係数閾値ｒｉｎより小さくない場合は、ステップＳ１４において、追加判定手段２５は、追加時寄与率記憶手段２０を参照し、変数群属性が未追加変数群Ｐの説明変数のうち現在の追加候補説明変数の次に寄与率の大きい次候補説明変数があるか否かを検索する。そして、次候補説明変数がなければ、変数増加手段１３における変数取込処理を終了し、ステップＳ２３に移行する。また、次候補説明変数があるときは、現在の追加候補説明変数ｘαｑ＋１を当該次候補説明変数に変更し、ステップＳ９に戻る。
【０１０１】
（変数除去処理）
変数除去処理については、従来のステップワイズ法と同様である。
【０１０２】
まず、ステップＳ１６において、除去時寄与率演算手段３１は、変数群属性記憶手段１１に記憶されている変数群属性テーブルを参照し、変数群属性が追加済変数群Ｑである説明変数ｘαｋを順次１つずつ選択する。
【０１０３】
次に、ステップＳ１７において、除去時寄与率演算手段３１は、当該選択した説明変数（以下「選択説明変数」という。）ｘαｋを除く追加済変数群Ｑの説明変数（ｘα１，…，ｘαｋ−１，ｘαｋ＋１，…，ｘαｑ）に、当該選択した説明変数ｘαｋを追加する場合の寄与率Ｆ（βαｋ）を算出し、除去時寄与率記憶手段３０に寄与率テーブルとして格納する。寄与率Ｆ（βαｋ）の算出は、変数取込処理における場合と同様である。
【０１０４】
次に、ステップＳにおいて、除去時寄与率演算手段３１は、変数群属性が追加済変数群Ｑである説明変数ｘαｋをすべて選択したかを判定し、まだ選択していない説明変数ｘαｋがあれば、ステップＳ１６に戻り、すべての説明変数ｘαｋについて寄与率Ｆ（βαｋ）を算出した場合には、次のステップＳ１９に移行する。
【０１０５】
次に、ステップＳ１９において、最小寄与率選択手段３２は、除去時寄与率記憶手段３０格納された寄与率テーブルを参照し、除去時寄与率演算手段３１により算出される寄与率が最小の説明変数を、除去候補説明変数ｘminに選択する。
【０１０６】
次に、ステップＳ２０において、除去基準判定手段３３は、最小寄与率選択手段３２により選択された除去候補説明変数ｘminの寄与率Ｆminが、重回帰式から除去する除去基準値Ｆｏｕｔより小さいか否かを判定する。
【０１０７】
寄与率Ｆminが除去基準値Ｆｏｕｔより小さい場合、ステップＳ２１において、変数群属性記憶手段１１に格納された変数群属性テーブルにおける当該除去候補説明変数ｘminの変数群属性を未追加変数群Ｐに変更し、追加済変数群Ｑの個数ｑを１だけ減少させ、未追加変数群Ｐの個数ｐを１だけ増加させる。そして、ｑ＝０となった場合には変数減少手段１４における変数除去処理を終了して、次のステップＳ２３に移行する。ｑ＞０の場合にはステップＳ１６に戻って変数除去処理を繰り返し行う。
【０１０８】
一方、寄与率Ｆminが除去基準値Ｆｏｕｔ以上の場合、他に追加すべき説明変数がないか調べるため、再び変数取込処理（ステップＳ３）に戻る。
【０１０９】
最後に、ステップＳ２３において、重回帰モデル出力手段７は、すべての選択済説明変数群Ｐの説明変数による重回帰係数を出力し、終了する。
【０１１０】
以上のように、本実施例の応答曲面モデル構築装置１によれば、ステップＳ１０〜Ｓ１５の処理において、追加候補説明変数ｘαｑ＋１の重回帰式への取り込みを行う際に、未追加変数群Ｐの他の各説明変数ｘｊとの相関係数ｒαｑ＋１，ｊを調べ、これらすべての相関係数ｒαｑ＋１，ｊが相関係数閾値ｒｉｎより小さい場合にのみ追加候補説明変数ｘαｑ＋１の重回帰式への取り込みを行うことで、多重共線性を正確に考慮し、精度のよい応答曲面モデルの構築を行うことができる。
【０１１１】
図６は、実際のエンジンの多項式モデルの構築において、本実施例の応答曲面モデル構築装置１により重回帰モデルを求めた結果の一例である。
【０１１２】
ステップワイズ法を直接適用した場合（図１０）と比較すると、全体的に回帰誤差が縮小し、また、目的変数が負値となるような異常値の出現が防止されていることがわかる。
【実施例２】
【０１１３】
図７は、本発明の実施例２に係る応答曲面モデル構築装置１の構成を表すブロック図である。本実施例では、図１の応答曲面モデル構築装置１の構成に加えて、重回帰モデル評価手段６が新たに追加されている点が異なる。他の構成は、実施例１と同様である。
【０１１４】
図８，図９は、実施例２の応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理を表すフローチャートである。
【０１１５】
図８，図９において、図４，図５と同様の処理を行うステップについては同符号を付して説明は省略する。
【０１１６】
本実施例においては、ステップＳ１の前に、ステップＳ０において、高次説明変数生成
手段３は、高次の説明変数の最大次数ｈｋを１に初期化する。そして、ステップＳ１におおては、高次説明変数生成手段３は、最大次数ｈｋまでの高次の説明変数を算出し標本値テーブルに格納する。尚、ｈｋ＝１の場合には、高次の説明変数はないため、ステップＳ１においては高次の説明変数は算出されないことになる。
【０１１７】
そして、ステップＳ１〜Ｓ２２までの処理は、実施例１と同様の処理が行われる。
【０１１８】
変数取込処理と説明変数除去処理（ステップＳ３〜Ｓ２２）がすべて終了すると、次に、ステップＳ２２ａにおいて、高次説明変数生成手段３は、最大次数ｈｋが予め定められた最終値ｈｍａｘに達したか否かを判定し、達していない場合は、ステップＳ２２ｂにおいてｈｋを１だけ増加させて、ステップＳ１に戻る。最大次数ｈｋが最終値ｈｍａｘに達した場合には、次のステップＳ２２ｃに移行する。尚、最終値ｈｍａｘは、必要となる可能性のある最大次数ｈｋがすべてカバーできる程度の値に適宜設定すればよく、通常は、ｈｍａｘ＝３〜８の値に設定される。
【０１１９】
ステップＳ２２ｃにおいて、重回帰モデル評価手段６は、最大次数ｈｋが１〜ｈｍａｘの場合において算出された各重回帰式の評価・選定を行う。ここでは、根平均自乗誤差（ＲＭＳＥ）による分析処理により評価を行う。
【０１２０】
まず、重回帰モデル評価手段６は、次式により、各最大次数ｈｋにおける根平均自乗誤差ｒｋを算出する。ここで、ｙ＾（ｎ）は、式（１２ｂ）に示した目的変数の推定値である。
【０１２１】
【数１４】
【０１２２】
そして、重回帰モデル評価手段６は、各最大次数ｈｋ（ｋ＝１，…，ｈｍａｘ）の最大値ｒｍａｘを求める。次いで、根平均自乗誤差ｒｋに基づいて、次式により各最大次数ｈｋにおけるコストＪｋを算出する。
【０１２３】
【数１５】
【０１２４】
ここで、ｗｐは累乗の大きさに対する重み係数で、０〜１の予め所定の値に定めておく。また、ｗｒはＲＭＳＥに対する重み係数でありｗｐ＋ｗｒ＝１である。
【０１２５】
そして、重回帰モデル評価手段６は、コストＪｋが最も小さくなる最大次数ｈｋを選択する。
【０１２６】
最後に、ステップＳ２３において、重回帰モデル出力手段７は、重回帰モデル評価手段６が選択した最大次数ｈｋに対する重回帰式の重回帰係数を出力し、終了する。
【０１２７】
このように、本実施例の応答曲面モデル構築装置１によれば、重回帰モデル評価手段６によって、最も適した最大次数ｈｋの重回帰式による重回帰モデルを選定でき、不必要に説明変数が増大することを防止することができる。
【符号の説明】
【０１２８】
応答曲面モデル構築装置１
測定データ記憶手段２
高次説明変数生成手段３
重回帰係数記憶手段４
説明変数選択手段５
重回帰モデル評価手段６
重回帰モデル出力手段７
変数群属性記憶手段１１
変数群属性初期化手段１２
変数増加手段１３
変数減少手段１４
変数選択処理制御手段１５
追加時寄与率記憶手段２０
追加時寄与率演算手段２１
最大寄与率選定手段２２
追加基準判定手段２３
相関係数演算手段２４
追加判定手段２５
追加候補説明変数変更手段２６
変数取込処理制御手段２７
除去時寄与率記憶手段３０
除去時寄与率演算手段３１
最小寄与率選択手段３２
除去基準判定手段３３
変数除去処理制御手段３４
エンジン１０１
制御パラメータ設定手段１０２
制御パラメータ測定手段１０３
特性値測定手段１０４
【技術分野】
【０００１】
本発明は、ディーゼルエンジン等の制御パラメータの応答曲面モデルを構築する応答曲面モデル構築技術に関する。
【背景技術】
【０００２】
自動車のディーゼルエンジンは、コモンレール方式による多段噴射装置、ＥＧＲ（排気再循環）装置、過給器などの多数の制御装置を装備している。従って、これらの多数の制御装置のパラメータを運転条件、燃費、多数の排出ガス成分を考慮して最適に設定し制御するには、実験データから多数の制御パラメータと特性値との関係を精度良く表現する応答曲面モデルの構築が必要となる。
【０００３】
かかる応答曲面モデルの構築方法としては、従来から、実験データに基づいて多項式モデルによって近似する方法が広く用いられている（例えば、非特許文献１参照）。多項式モデルの作成においては、制御パラメータ等の目的変数に対して、その予測に有効な説明変数を特性値から選択する必要がある。この選択の方法としては、ステップワイズ法（変数増減法）等が提案されている（非特許文献２，３，特許文献１）。以下、特許文献１に記載のエンジン制御パラメータの適合方法について簡単に概説する。
【０００４】
まず、エンジンの電子制御スロットルの開度、点火時期、燃料噴射時期、ＥＧＲバルブの開度、給気バルブの開弁時期などの制御パラメータを（ｘ１，ｘ２，…，ｘ５）とおく。説明変数を、これら５つの制御パラメータ（ｘ１，ｘ２，…，ｘ５）、及びこれらの制御パラメータを自乗し又はこれらの制御パラメータの２つを掛け合わせて得られる１５個の２次項（ｘ１２，ｘ２２，…ｘ５２，ｘ１ｘ２，ｘ１ｘ３，…，ｘ４ｘ５）とし、これら２０個の説明変数（ｘ１，ｘ２，…，ｘ５，ｘ１２，ｘ２２，…ｘ５２，ｘ１ｘ２，ｘ１ｘ３，…，ｘ４ｘ５）を（ｘ１，ｘ２，…，ｘ２０）と記す。また、目的変数を、ＮＯｘ排出量や燃料消費量等のエンジンの各種特性値とし、これをｙとおく。説明変数ｘｉに対する偏回帰係数をａｉとし、重回帰式を次式で表す。
【０００５】
【数１】
【０００６】
各制御パラメータ（ｘ１，ｘ２，…，ｘ５）を振りながらｎ回の計測を行う。各説明変数及び目的変数の計測値をｘｉ（ｍ）（ｉ＝１，…，２０），ｙ（ｍ）（ｍ＝１，…，ｎ）とする。これらの計測値より、偏差積和Ｓｉｊ，Ｓｙｙを次式のように計算する。
【０００７】
【数２】
【０００８】
さらに、算出された偏差積和Ｓｉｊ，Ｓｙｙを用いて相関係数ｒｉｊを次式で算出する。
【０００９】
【数３】
【００１０】
以上の計算を行った後、まず、算出された相関係数ｒｉｊを用いて、初期の偏回帰係数ａｉ（ｉ＝１，…，５）を次式により算出する。
【００１１】
【数４】
【００１２】
次に、以下のような説明変数の追加処理及び説明変数の削除処理を繰り返すことによって、最適な説明変数の組み合わせによる重回帰式を求め、応答曲面モデルを構築する。
【００１３】
（説明変数の追加処理）
現在の重回帰式にｐ個（ｐ＜２０）の説明変数が含まれているとし、ｐ＋１個目の説明変数ｘｋを追加することを試みるものとする。まず、説明変数がｐ個の重回帰式の決定係数Ｒｐを次式により算出する。ここで、Ｙｐ（ｍ）はｍ番目の目的変数の計測値に対するｐ個の説明変数の重回帰式による予測値である。
【００１４】
【数５】
【００１５】
次に、（ｐ＋１）個目の説明変数ｘｋを追加した（ｐ＋１）個の重回帰式の決定係数Ｒｐ＋１を次式により算出する。
【００１６】
【数６】
【００１７】
これを、追加することが可能なすべての説明変数ｘｋについて行い、決定係数Ｒｐ＋１が最大となる説明変数ｘｋを、新たに追加する説明変数の候補として選出する。
【００１８】
次に、決定係数Ｒｐ＋１が最大となる説明変数ｘｋを追加した重回帰式において、追加した説明変数ｘｋの重回帰式に於ける寄与率を示すＦ値を次式より算出する。
【００１９】
【数７】
【００２０】
そして、このＦｐ＋１を所定の判定値Ｆｉｎと比較し、Ｆｐ＋１＞Ｆｉｎの場合には候補の説明変数ｘｋを追加して、次の説明変数の削除処理に移る。一方、Ｆｐ＋１≦Ｆｉｎの場合には、現在の重回帰式を応答曲面モデルとして出力して終了する。
【００２１】
（説明変数の削除処理）
現在の重回帰式にｐ＋１個（ｐ＋１＞５）の説明変数が含まれているとし、１つの説明変数ｘｋ（１５個の２次項ｘ６〜ｘ２０のうちの何れか）を削除することを試みるものとする。
【００２２】
まず、説明変数としてｘ１，ｘ２，…，ｘｋ−１，ｘｋ＋１，…，ｘｐ＋１のｐ個の変数に対して説明変数ｘｋを追加する場合のＦ値Ｆｐ＋１を式（７）により算出する。そして、このＦｐ＋１が所定の判定値Ｆｏｕｔ（Ｆｏｕｔ＜Ｆｉｎ）よりも小さい場合、この説明変数ｘｋを削除し、そうでなければ説明変数ｘｋを残す。この削除判定を、現在の重回帰式に含まれるすべての２次項の説明変数に対して行い、その後、先の説明変数の追加処理に移る。
【００２３】
以上のようにして、最適なエンジン制御パラメータを含む重回帰式により表される応答曲面モデルを構築する。
【先行技術文献】
【特許文献】
【００２４】
【特許文献１】特開２００４−６８７２９
【非特許文献】
【００２５】
【非特許文献１】山崎光悦，“応答曲面法の展開”，日本機械学会誌，日本国，社団法人日本機械学会，２００６年５月，第１０９巻，第１０５０号，pp.377-379．
【非特許文献２】管民郎，「多変量解析の実践（上）」，初版，上巻，日本国，株式会社現代数学社，1993年12月，pp.70-74．
【非特許文献３】河口至商，「多変量解析入門Ｉ」，初版，日本国，森北出版株式会社，1973年6月，pp.3-33．
【発明の概要】
【発明が解決しようとする課題】
【００２６】
しかしながら、上記従来のエンジン制御パラメータの適合方法においては、変数選択に於いて回帰係数のＦ検定による有意性の判断を行っているものの、厳密には回帰係数の検定のみでは多重共線性が考慮されているとはいえず、説明変数間の相関係数が高い説明変数の組み合わせが選択される場合がある。従って、実際のエンジンの多項式モデルの構築にあたっては、単純にステップワイズ法を各制御パラメータの説明変数の選択に直接適用するだけでは推定精度や実用上の観点から不十分である。
【００２７】
図１０は、実際のエンジンの多項式モデルの構築において、各制御パラメータの説明変数の選択にステップワイズ法を直接適用して重回帰モデルを求めた結果の一例である。図１０において、横軸は測定番号、縦軸はＮＯｘ値を表す。測定は７６回行い、これら７６回の測定から得られる７６組の制御パラメータ及びＮＯｘの測定値の組から、上記従来のステップワイズ法を用いて最適な重回帰式を求めた。重相関係数は０．９４９１であった。黒四角の点が実測値、白丸の点が重回帰式からの推定値である。図１０より、ステップワイズ法を直接適用した場合、ＮＯｘ値の変動が大きい点において誤差が大きく、特に、予測値が負値となる異常な予測点が複数出現しており、従来のエンジン制御パラメータの適合方法は、なお予測精度の改善が必要であることが分かる。そのためには、推定精度向上のために制御パラメータの非線形項を組み合わせた説明変数を複数作成し多項式モデルに導入することや、異常な予測値を出現させないために各説明変数を選択する際に多重共線性を厳密に考慮させることが必要である。また、１次、２次、及び高次の多項式モデルを生成し比較することにより、できる限り少ない説明変数で予測精度の高い多項式モデルを効率的に構築する方法が必要である。
【００２８】
そこで、本発明の目的は、エンジンの特性値と制御パラメータとの関係を示す測定データに基づいてその関係を精度良く再現する応答曲面モデルを効率的に構築する応答曲面モデル構築装置を提供することにある。
【課題を解決するための手段】
【００２９】
本発明に係る応答曲面モデル構築装置の第１の構成は、エンジンの各種制御パラメータである１次の説明変数及びそれに対するエンジンの特性値である目的変数の測定値の組である測定値セットを測定回数だけ格納する測定データ記憶手段と、
各測定値セットについて、各制御パラメータの累乗項を含むｋ個（ｋ＝２，…，Ｍ；Ｍは重回帰式に取り込む説明変数の次数の最大値）の制御パラメータの積項からなる高次の説明変数を算出し、前記測定データ記憶手段に保存する高次説明変数生成手段と、
前記各説明変数から前記目的変数を予測する重回帰式の偏回帰係数及び切片である重回帰係数を記憶する重回帰係数記憶手段と、
前記測定データ記憶手段に格納されたすべての説明変数の中からから重回帰式による目的変数の予測に有効な説明変数を選択するとともに、選択した各説明変数による重回帰式の重回帰係数を算出し前記重回帰係数記憶手段に格納する説明変数選択手段と、を備え、
前記説明変数選択手段は、
各説明変数について、当該説明変数が、重回帰式に取り込まれない予定の未追加変数群に属するのか、重回帰式に取り込む予定の追加変数群に属するのかを示す変数群属性を登録する変数群属性テーブルが格納された変数群属性記憶手段と、
前記変数群属性テーブルの各説明変数の変数群属性を、すべて未追加変数群に初期化する変数群属性初期化手段と、
変数群属性が未追加変数群である説明変数のうち、予測に有効な説明変数の変数群属性を追加変数群に変更する変数取込処理を行う変数増加手段と、
変数群属性が追加変数群である説明変数のうち、予測に有効でない説明変数の変数群属性を未追加変数群に変更する変数除去処理を行う変数減少手段と、
前記変数取込処理及び前記変数除去処理を交互に繰り返し実行し、前記変数群属性テーブルにおいて変数群属性が追加変数群の説明変数が未追加変数群の説明変数よりも少なくなるか、又は前記変数取込処理において何れの説明変数の変数群属性も追加変数群に変更されなかった場合に繰り返し処理を終了する変数選択処理制御手段と、を備え、
前記変数増加手段は、
変数群属性が未追加変数群の説明変数を順次１つずつ選択し、選択した説明変数，変数群属性が追加変数群の説明変数，及び目的変数の測定値セットを前記測定データ記憶手段からすべて読み出し、読み出した測定値セットから重回帰係数を算出し、読み出した測定値セット及び当該重回帰係数から、選択した説明変数の寄与率を算出する追加時寄与率演算手段と、
前記変数群属性テーブルにおいて変数群属性が未追加変数群の説明変数のうち、前記追加時寄与率演算手段により算出された寄与率が最大の説明変数を初期の追加候補説明変数に選択する最大寄与率選定手段と、
前記追加候補説明変数の寄与率が、重回帰式に取り込む取込基準値より大きいか否かを判定し、当該取込基準値より大きくない場合、前記変数取込処理を終了する追加基準判定手段と、
前記追加基準判定手段が、当該追加候補説明変数の寄与率が取込基準値より大きいと判定した場合、当該追加候補説明変数と、前記変数群属性テーブルにおいて変数群属性が未追加変数群の他の各説明変数との相関係数を順次算出する相関係数演算手段と、
前記相関係数演算手段が算出する相関係数がすべて所定の相関係数閾値より小さい場合、前記変数群属性テーブルにおける当該追加候補説明変数の変数群属性を追加変数群に変更し、前記変数取込処理を終了する追加判定手段と、
前記相関係数演算手段が算出する相関係数の何れかが前記相関係数閾値より大きい場合、変数群属性が未追加変数群の説明変数のうち現在の追加候補説明変数の次に寄与率が大きい次候補説明変数がなければ前記変数取込処理を終了し、次候補説明変数があるときは、現在の追加候補説明変数を当該次候補説明変数に変更する追加候補説明変数変更手段と、
前記追加基準判定手段、前記相関係数演算手段、前記追加判定手段、及び前記追加候補説明変数変更手段による処理を繰り返し実行する変数取込処理制御手段と、を備え、
前記変数減少手段は、
前記変数群属性テーブルにおける変数群属性が追加変数群の説明変数を順次１つずつ選択し、当該選択した説明変数を除く追加変数群の説明変数に、当該選択した説明変数を追加する場合の寄与率を算出する除去時寄与率演算手段と、
前記除去時寄与率演算手段により算出される寄与率が最小の説明変数を、除去候補説明変数に選択する最小寄与率選択手段と、
当該除去候補説明変数の寄与率が、重回帰式から除去する除去基準値より小さい場合、前記変数群属性テーブルにおける当該除去候補説明変数の変数群属性を未追加変数群に変更し、それ以外の場合には前記変数除去処理を終了する除去基準判定手段と、
前記最小寄与率選択手段及び前記除去基準判定手段による処理を、前記除去基準判定手段が前記変数除去処理を終了するまで反復実行する制御を行う変数除去処理制御手段と、を備えたことを特徴とする。
【００３０】
この構成によれば、変数増加手段における処理において、追加候補説明変数ｘαｑ＋１の重回帰式への取り込みを行う際に、未追加変数群Ｐの他の各説明変数ｘｊとの相関係数ｒαｑ＋１，ｊを調べ、これらすべての相関係数ｒαｑ＋１，ｊが相関係数閾値ｒｉｎより小さい場合にのみ追加候補説明変数ｘαｑ＋１の重回帰式への取り込みを行うことで、多重共線性を正確に考慮し、精度のよい応答曲面モデルの構築を行うことができる。
【００３１】
また、本発明に係る応答曲面モデル構築装置の第２の構成は、前記第１の構成において、前記寄与率は、Ｆ値又はＴ統計量であることを特徴とする。
【００３２】
また、本発明に係る応答曲面モデル構築装置の第３の構成は、前記第１又は２の構成において、前記説明変数選択手段は、高次の説明変数の最大次数ｈｋが１からｈｍａｘのそれぞれの場合について重回帰式の重回帰係数を算出し、前記重回帰係数記憶手段に格納するものであり、
前記重回帰係数記憶手段に格納された最大次数ｈｋ（ｋ＝１，…，ｈｍａｘ）の場合の重回帰式のそれぞれに対して、根平均自乗誤差ｒｋを算出し、算出した各根平均自乗誤差ｒｋの最大値ｒｍａｘを抽出するとともに、各根平均自乗誤差ｒｋ及びその最大値ｒｍａｘに基づいて、各最大次数ｈｋにおけるコストＪｋを算出し、コストＪｋが最大となる最大次数ｈｋに対する重回帰式を選択し出力する重回帰モデル評価手段を備えたことを特徴とする。
【００３３】
この構成により、重回帰モデル評価手段によって、最も適した最大次数ｈｋの重回帰式による重回帰モデルを選定でき、不必要に説明変数が増大することを防止することができる。
【００３４】
また、本発明に係るプログラムは、コンピュータに読み込んで実行させることにより、前記第１乃至３のいずれかの構成の応答曲面モデル構築装置として機能させることを特徴とする。
【発明の効果】
【００３５】
以上のように、本発明によれば、変数増加手段における処理において、追加候補説明変数ｘαｑ＋１の重回帰式への取り込みを行う際に、未追加変数群Ｐの他の各説明変数ｘｊとの相関係数ｒαｑ＋１，ｊを調べ、これらすべての相関係数ｒαｑ＋１，ｊが相関係数閾値ｒｉｎより小さい場合にのみ追加候補説明変数ｘαｑ＋１の重回帰式への取り込みを行うことで、多重共線性を正確に考慮し、精度のよい応答曲面モデルの構築を行うことができる。
【００３６】
また、さらに重回帰モデル評価手段によって、最も適した最大次数ｈｋの重回帰式による重回帰モデルを選定すれば、不必要に説明変数が増大することを防止することができる。
【図面の簡単な説明】
【００３７】
【図１】本発明の実施例１に係る応答曲面モデル構築装置の構成を表すブロック図である。
【図２】変数増加手段１３の構成を表すブロック図である。
【図３】変数減少手段１４の構成を表すブロック図である。
【図４】本実施例に係る応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理を表すフローチャートである。
【図５】本実施例に係る応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理を表すフローチャートである。
【図６】実際のエンジンの多項式モデルの構築において、本実施例の応答曲面モデル構築装置１により重回帰モデルを求めた結果の一例である。
【図７】本発明の実施例２に係る応答曲面モデル構築装置１の構成を表すブロック図である。
【図８】実施例２の応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理を表すフローチャートである。
【図９】実施例２の応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理を表すフローチャートである。
【図１０】実際のエンジンの多項式モデルの構築において、各制御パラメータの説明変数の選択にステップワイズ法を直接適用して重回帰モデルを求めた結果の一例である。
【発明を実施するための形態】
【００３８】
以下、本発明を実施するための形態について、図面を参照しながら説明する。
【実施例１】
【００３９】
（１）応答曲面モデル構築装置の構成
（１−１）全体構成
図１は、本発明の実施例１に係る応答曲面モデル構築装置の構成を表すブロック図である。
【００４０】
図１において、応答曲面モデル構築装置１は、エンジン１０１のＮＯｘ排出量や燃料消費量等の各種特性値を、所定の目標値の範囲内となるように、エンジン１０１の電子制御スロットルの開度、点火時期、燃料噴射時期、ＥＧＲバルブの開度、吸気バルブの開弁時期などの各種制御パラメータを適合させる応答曲面モデルの構築を行う装置である。すなわち、応答曲面モデル構築装置１は、エンジン１０１の各種制御パラメータを説明変数、エンジン１０１の種特性値を目標変数とする重回帰式の算出を行う。ここで、応答曲面モデルの重回帰式の算出を行うため、各種制御パラメータを人為的に変動させたときの各種特性値の測定データを予め用意しておく必要がある。そこで、制御パラメータ設定手段１０２は、エンジン１０１に対して予め定められた各種制御パラメータの組み合わせに基づいて、順次各種制御パラメータを設定する。そして、制御パラメータ測定手段１０３及び特性値測定手段１０４は、エンジン１０１の各種制御パラメータの設定が変更される毎に、エンジン１０１の各種制御パラメータ及び特性値を測定し、測定データ記憶手段２に格納する。
【００４１】
応答曲面モデル構築装置１は、測定データ記憶手段２、高次説明変数生成手段３、重回帰係数記憶手段４、説明変数選択手段５、及び重回帰モデル出力手段７を備えている。尚、この応答曲面モデル構築装置１は、マイコンやＦＰＧＡなどによりハードウェア的に構成することもできるが、プログラムとして用意してパソコンなどの汎用コンピュータにこのプログラムをロードすることによってソフトウェア的に構成することもできる。
【００４２】
測定データ記憶手段２は、上述したように、制御パラメータ測定手段１０３及び特性値測定手段１０４により測定される、エンジン１０１の各種制御パラメータである１次の説明変数（ｘ１，ｘ２，…，ｘＮ）及びそれに対するエンジン１０１の特性値である目的変数ｙの測定値の組である測定値セット｛（ｘ１，ｘ２，…，ｘＮ），ｙ｝を測定回数Γ分だけ格納するためのものである。ここで、Ｎは制御パラメータの数である。
【００４３】
なお、測定データ記憶手段２において、各測定回に対する測定値セット｛（ｘ１，ｘ２，…，ｘＮ），ｙ｝は、表１に示すような標本値テーブルとして格納される。
【００４４】
【表１】
【００４５】
表１において、「測定番号」のフィールドは、各測定回の番号が格納されるフィールドである。「制御パラメータ」のｘ１〜ｘＮまでのフィールドは、各測定回の１次の説明変数が格納されるフィールドである。また、ｘＮ＋１〜ｘＮ＋Ｌまでのフィールドは、各測定回の高次の説明変数が格納されるフィールドである。「特性値」のフィールドは、各測定回の特性値ｙが格納されるフィールドである。ここで、「高次の説明変数」とは、各制御パラメータの累乗項（ｘｉ２，ｘｉ３等）を含むｋ個（ｋ＝２，…，Ｍ；Ｍは重回帰式に取り込む説明変数の次数の最大値）の制御パラメータの積項からなる説明変数をいう。Ｌは高次の説明変数の数を表す。重回帰式に取り込む説明変数の次数の最大値をＭとすると、Ｌ＝ＮＣＭ＝Ｎ！／（Ｍ！・（Ｎ−Ｍ）！）となる。以下、ｎ番目の測定回に於ける説明変数（ｘ１，ｘ２，…，ｘＮ＋Ｌ）及び目的変数ｙを、それぞれ（ｘ１（ｎ），ｘ２（ｎ），…，ｘＮ＋Ｌ（ｎ））及びｙ（ｎ）と記す。また、高次の説明変数については、どの説明変数を掛け合わせた積項かを明示する必要がある場合、２次の説明変数ｘｉｘｊをｘｉｊ、３次の説明変数ｘｉｘｊｘｋをｘｉｊｋのように記すこととし、一般に、ｎ次の説明変数ｘα１…ｘαｎをｘα１…αｎのように記す。
【００４６】
高次説明変数生成手段３は、測定データ記憶手段２に格納された各測定値セット｛（ｘ１（ｎ），ｘ２（ｎ），…，ｘＮ（ｎ）），ｙ（ｎ）｝（ｎ＝１，…，Γ）について、それぞれの高次の説明変数（ｘＮ＋１（ｎ），ｘＮ＋２（ｎ），…，ｘＮ＋Ｌ（ｎ））を算出し、前記測定データ記憶手段２に保存する。
【００４７】
重回帰係数記憶手段４は、各説明変数ｘα１…αｎ（ｎ＝１，…，Ｍ）から前記目的変数ｙを予測する重回帰式の重回帰係数を記憶する。ここで、「重回帰係数」とは、重回帰式の回帰係数βα１…αｎ及び切片β０をいう。尚、偏回帰係数βα１…αｎは、重回帰式における説明変数ｘα１…αｎの偏回帰係数を表す。
【００４８】
説明変数選択手段５は、測定データ記憶手段２に格納されたすべての説明変数の中からから重回帰式による目的変数の予測に有効な説明変数を選択するとともに、選択した各説明変数による重回帰式の重回帰係数を算出し重回帰係数記憶手段４に格納する。この説明変数選択手段５の構成の詳細については、後述する。
【００４９】
重回帰モデル出力手段７は、最終的に選択された説明変数名と、算出された各重回帰係数を出力する。
【００５０】
（１−２）説明変数選択手段５の構成
説明変数選択手段５は、図１に示すように、変数群属性記憶手段１１、変数群属性初期化手段１２、変数増加手段１３、変数減少手段１４、及び変数選択処理制御手段１５を備えている。
【００５１】
変数群属性記憶手段１１は、上記Ｎ＋Ｌ個の各説明変数について、変数群属性を登録する変数群属性テーブルを記憶する。ここで、「変数群属性」とは、説明変数が、未追加変数群Ｐに属するのか、追加済変数群Ｑに属するのかを示す属性をいう。また、「未追加変数群」とは、重回帰式に取り込まれない予定の変数の集合をいい、「追加済変数群」とは、重回帰式に取り込む予定の変数の集合をいう。表２に、変数群属性テーブルを示す。表２において、ｇｉ（ｉ＝１，２，…，Ｎ＋Ｌ）は、目的変数ｘｉに対する変数群属性であり、「未追加変数群」又は「追加済変数群」の何れかの値をとる。
【００５２】
【表２】
【００５３】
変数群属性初期化手段１２は、変数群属性記憶手段１１に格納された変数群属性テーブルの各説明変数の変数群属性（ｇ１，…ｇＮ＋Ｌ）を、すべて未追加変数群Ｐに初期化する。
【００５４】
変数増加手段１３は、変数群属性が未追加変数群Ｐである説明変数のうち、予測に有効な説明変数の変数群属性を追加済変数群Ｑに変更する変数取込処理を行う。この変数増加手段１３の構成の詳細については、後述する。
【００５５】
変数減少手段１４は、変数群属性が追加済変数群Ｑである説明変数のうち、予測に有効でない説明変数の変数群属性を未追加変数群Ｐに変更する変数除去処理を行う。この変数減少手段１の構成の詳細については、後述する。
【００５６】
変数選択処理制御手段１５は、変数群属性初期化手段１２、変数増加手段１３及び変数減少手段１４を制御して、最初に変数群属性テーブルを初期化した後、変数取込処理及び前記変数除去処理を交互に繰り返し実行し、変数群属性テーブルにおいて変数群属性が追加済変数群Ｑの説明変数が未追加変数群Ｐの説明変数よりも少なくなるか、又は前記変数取込処理において何れの説明変数の変数群属性も追加済変数群Ｑに変更されなかった場合に繰り返し処理を終了する処理を行う。
【００５７】
（１−３）変数増加手段１３の構成
図２は、変数増加手段１３の構成を表すブロック図である。図２に示すように、変数増加手段１３は、追加時寄与率記憶手段２０、追加時寄与率演算手段２１、最大寄与率選定手段２２、追加基準判定手段２３、相関係数演算手段３４、追加判定手段２５、追加候補説明変数変更手段２６、及び変数取込処理制御手段２７を備えている。
【００５８】
追加時寄与率演算手段２１は、変数群属性テーブルを参照しながら、変数群属性が未追加変数群Ｐの説明変数ｘαｑ＋１（ｋ）（ｋ＝１，…，Ｎ＋Ｌ−ｑ）を順次１つずつ選択する。ここで、ｑは変数群属性が追加済変数群Ｑである説明変数の現在の個数である。そして、選択した説明変数（以下「選択説明変数」という。）ｘαｑ＋１（ｋ），現在の変数群属性が追加済変数群Ｑの説明変数（ｘα１，ｘα２，…，ｘαｑ），及び目的変数ｙの測定値セット｛（ｘα１（ｎ），ｘα２（ｎ），…，ｘαｑ（ｎ），ｘαｑ＋１（ｋ）（ｎ）），ｙ（ｎ）｝（ｎ＝１，…，Ｐ）を測定データ記憶手段２からＰ組分すべて読み出し、読み出した測定値セットから重回帰係数Ｆ（βαｑ＋１（ｋ））を算出し、重回帰係数記憶手段４に、表３のような一時回帰係数テーブルとして格納する。ここで、β０（ｋ）は重回帰式の切片、βαｉ（ｋ）（ｉ＝１，…，ｑ）は説明変数ｘαｋに対する偏回帰係数、βαｑ＋１（ｋ）は説明変数ｘαｑ＋１（ｋ）に対する偏回帰係数である。さらに、追加時寄与率演算手段２１は、読み出した測定値セット及び算出した重回帰係数から、選択説明変数ｘαｑ＋１（ｋ）の寄与率Ｆ（βαｑ＋１（ｋ））を算出し、追加時寄与率記憶手段２０に、表４のような寄与率テーブルとして格納する。ここで、Ｆ（βαｑ＋１（ｋ））は、選択説明変数ｘαｑ＋１（ｋ）を追加済変数群Ｑに追加した場合の当該選択説明変数ｘαｑ＋１（ｋ）の寄与率である。
【００５９】
【表３】
【００６０】
【表４】
【００６１】
最大寄与率選定手段２２は、追加時寄与率記憶手段２０から、追加時寄与率演算手段２１により算出された寄与率を読み出し、それらのうち、最大の寄与率に対応する選択説明変数ｘαｑ＋１（ｋ）を、初期の追加候補説明変数ｘαｑ＋１とする。
【００６２】
追加基準判定手段２３は、追加候補説明変数ｘαｑ＋１の寄与率Ｆ（βαｑ＋１）が、重回帰式に取り込む取込基準値Ｆｉｎより大きいか否かを判定し、当該取込基準値Ｆｉｎより大きくない場合には、変数増加手段１３における変数取込処理を終了する。
【００６３】
相関係数演算手段２４は、追加基準判定手段２３が、追加候補説明変数ｘαｑ＋１の寄与率Ｆ（βαｑ＋１）が取込基準値Ｆｉｎより大きいと判定した場合、当該追加候補説明変数ｘαｑ＋１と、変数群属性テーブルにおいて変数群属性が未追加変数群Ｐの他の各説明変数との相関係数を順次算出する。
【００６４】
追加判定手段２５は、相関係数演算手段２４が算出する相関係数がすべて所定の相関係数閾値ｒｉｎより小さい場合、変数群属性テーブルにおける当該追加候補説明変数ｘαｑ＋１の変数群属性を追加済変数群Ｑに変更し、変数増加手段１３における変数取込処理を終了する。
【００６５】
追加候補説明変数変更手段２６は、相関係数演算手段３４が算出する相関係数の何れかが相関係数閾値ｒｉｎ以上の場合、追加時寄与率記憶手段２０を参照し、変数群属性が未追加変数群Ｐの説明変数のうち現在の追加候補説明変数の次に寄与率の大きい説明変数（以下「次候補説明変数」という。）があるか否かを検索する。そして、次候補説明変数がなければ、変数増加手段１３における変数取込処理を終了する。また、次候補説明変数があるときは、現在の追加候補説明変数ｘαｑ＋１を当該次候補説明変数に変更する。
【００６６】
変数取込処理制御手段２７は、変数増加手段１３における変数取込処理の全体制御を行う。すなわち、追加時寄与率演算手段２１及び最大寄与率選定手段２２による処理を最初に行った後、追加基準判定手段２３、相関係数演算手段２４、追加判定手段２５、及び追加候補説明変数変更手段２６による処理を繰り返し実行する制御を行う。
【００６７】
（１−４）変数減少手段１４の構成
図３は、変数減少手段１４の構成を表すブロック図である。変数減少手段１４は、図３に示すように、除去時寄与率記憶手段３０、除去時寄与率演算手段３１、最小寄与率選択手段３２、除去基準判定手段３３、及び変数除去処理制御手段３４を備えている。
【００６８】
除去時寄与率演算手段３１は、変数群属性記憶手段１１に記憶されている変数群属性テーブルを参照し、変数群属性が追加済変数群Ｑである説明変数ｘαｋを順次１つずつ選択し、当該選択した説明変数（以下「選択説明変数」という。）ｘαｋを除く追加済変数群Ｑの説明変数（ｘα１，…，ｘαｋ−１，ｘαｋ＋１，…，ｘαｑ）に、当該選択した説明変数ｘαｋを追加する場合の寄与率Ｆ（βαｋ）を算出し、除去時寄与率記憶手段３０に寄与率テーブルとして格納する。ここで、ｑは変数群属性が追加済変数群Ｑである説明変数の現在の個数である。また、寄与率テーブルは表４に示したものと同様である。
【００６９】
最小寄与率選択手段３２は、除去時寄与率記憶手段３０格納された寄与率テーブルを参照し、除去時寄与率演算手段３１により算出される寄与率が最小の説明変数を、除去候補説明変数ｘminに選択する。
【００７０】
除去基準判定手段３３は、最小寄与率選択手段３２により選択された除去候補説明変数ｘminの寄与率Ｆminが、重回帰式から除去する除去基準値Ｆｏｕｔより小さい場合、変数群属性記憶手段１１に格納された変数群属性テーブルにおける当該除去候補説明変数ｘminの変数群属性を未追加変数群Ｐに変更する。それ以外の場合には、変数減少手段１４における変数除去処理を終了する。
【００７１】
変数除去処理制御手段３４は、変数減少手段１４における変数除去処理の全体制御を行う。すなわち、除去時寄与率演算手段３１による寄与率テーブルの生成を実行した後、最小寄与率選択手段３２及び除去基準判定手段３３による処理を、除去基準判定手段３３が変数除去処理を終了させるまで反復実行する制御を行う。
【００７２】
（２）応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理
以上のように構成された本実施例に係る応答曲面モデル構築装置１について、次に、応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理の詳細について説明する。
【００７３】
図４，図５は、本実施例に係る応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理を表すフローチャートである。
【００７４】
まず、ステップＳ１において、高次説明変数生成手段３は、測定データ記憶手段２に格納された標本値テーブルから、Γ組の測定値セット（ｘ１（ｎ），ｘ２（ｎ），…，ｘＮ（ｎ））（ｎ＝１，…，Γ）を順次読み出し、それぞれ、Ｌ個の高次の説明変数を算出し標本値テーブルに格納する。
【００７５】
次に、ステップＳ２において、変数群属性初期化手段１２は、標本値テーブルに格納したＮ＋Ｌ個の説明変数について、変数群属性記憶手段１１に表２のような変数群属性テーブルを生成し、変数群属性テーブルの各説明変数ｘｉ（ｉ＝１，…，Ｎ＋Ｌ）の変数群属性ｇｉを未追加変数群Ｐに初期化し、未追加変数群Ｐの変数の個数ｐを０，追加済変数群Ｑの変数の個数ｑをＮ＋Ｌとする。
【００７６】
次に、変数取込処理制御手段２７は、以下の変数取込処理（Ｓ３〜Ｓ１５）と変数除去処理（Ｓ１６〜Ｓ２２）を反復実行し、多重共線性を考慮して有効な説明変数を選択するとともに、選択した説明変数による重回帰式を算出する。
【００７７】
（変数取込処理）
変数取込処理では、まず、ステップＳ３において、変数取込処理制御手段２７は、未追加変数群Ｐの変数の個数ｐと追加済変数群Ｑの変数の個数ｑを比較し、ｑ＜ｐの場合には、変数取込処理を打ち切って、ステップＳ２３に移行する。ｑ≧ｐの場合は、次のステップＳ４に進む。
【００７８】
次に、ステップＳ４において、追加時寄与率演算手段２１は、変数属性テーブルの変数属性が未追加変数群Ｐである説明変数ｘαｑ＋１（ｋ）（ｋ＝１，…，ｐ）を順次１つずつ選択する。ここで、ｐ＝Ｎ＋Ｌ−ｑである。
【００７９】
次に、ステップＳ５において、追加時寄与率演算手段２１は、追加済変数群Ｑの説明変数（ｘα１，ｘα２，…，ｘαｑ）と選択説明変数ｘαｑ＋１（ｋ）との重回帰式の重回帰係数βｉ（ｋ）（ｉ＝０，１，…，ｑ，ｑ＋１）を算出し、表３に示した一時回帰係数テーブルに格納する。ここで、選択説明変数ｘαｑ＋１（ｋ）に対する重回帰式は次式（８）のように表される。
【００８０】
【数８】
【００８１】
測定データ記憶手段２の標本値テーブルに格納された測定値セット｛（ｘα１（ｎ），ｘα２（ｎ），…，ｘαｑ（ｎ），ｘαｑ＋１（ｋ）（ｎ）），ｙ（ｎ）｝（ｎ＝１，…，Γ）から重回帰係数βｉ（ｋ）（ｉ＝０，１，…，ｑ，ｑ＋１）を算出する方法については、周知の線形回帰の計算方法が使用される（例えば、非特許文献３参照）。具体的には、目的変数ベクトルＹ、説明変数行列Ｘ、重回帰係数ベクトルＢを次式（９）のように置く。
【００８２】
【数９】
【００８３】
このとき、重回帰係数ベクトルＢは、次式（１０）により算出することができる。但し、Ｘｔは行列Ｘの転置行列を表す。
【００８４】
【数１０】
【００８５】
次に、ステップＳ６において、追加時寄与率演算手段２１は、選択説明変数ｘαｑ＋１（ｋ）に対する寄与率Ｆ（βαｑ＋１（ｋ））を算出し、追加時寄与率記憶手段２０に、表４に示した寄与率テーブルを格納する。具体的には、次のような演算が行われる。
【００８６】
まず、次式のような（ｑ＋１）個の説明変数（ｘα１，ｘα２，…，ｘαｑ，ｘαｑ＋１（ｋ））の分散共分散行列Σの逆行列Σ−１を算出する。ここで、Ｅ［ ］は平均を表す。
【００８７】
【数１１】
【００８８】
また、次式により目的変数ｙの予測誤差二乗和Ｓｅ及び回帰式の誤差分散Ｖｅを算出する。ここで、Γは標本値テーブルに格納された測定値セットの組数（測定回数）である。
【００８９】
【数１２】
【００９０】
そして、これらの算出値から、選択説明変数ｘαｑ＋１（ｋ）に対する寄与率Ｆ（βαｑ＋１（ｋ））を次式により算出する。ここで、Ｓ（βαｑ＋１（ｋ））２は不偏分散、Ｔ（βαｑ＋１（ｋ））はＴ値（Ｔ統計量）である。
【００９１】
【数１３】
【００９２】
ここでは寄与率としてＦ値を使用した。Ｆ値は補修団の重回帰係数βｊがβｊ＝０という仮説のもとで、自由度（１，Ｎ−ｑ−１）のＦ分布に従うとする。Ｆ値がＦ分布の限界値と比較してＦ（βｊ）＞Ｆ（１，Ｎ−ｑ−１，α）（右辺のＦ（ ）はＦ関数）ならば仮説を棄却し、係数βｊが有意であるとする。ここで、αは有意水準である。すなわち、ここではＦ値がＦｉｎの基準値よりも大きい偏回帰係数を仮採用し、それ以外は棄却する処理を行っている。
【００９３】
尚、寄与率としては上式（１３ｂ）のＴ値を使用してもよい。Ｔ値を使用する場合、母集団の偏回帰係数βｊがβｊ＝０という仮説のもとで自由度（Ｎ−ｑ−１）のｔ分布に従うとする。この場合、Ｔ（βｊ）＞ｔ（Ｎ−ｑ−１，α／２）ならば仮説を棄却し、係数βｊは有意であるとする。
【００９４】
次に、ステップＳ７において、追加時寄与率演算手段２１は、未追加変数群Ｐし属するｐ個の説明変数ｘαｑ＋１（ｋ）（ｋ＝１，…，ｐ）についてすべて寄与率Ｆ（βαｑ＋１（ｋ））を算出したか否かを判定し、すべての説明変数ｘαｑ＋１（ｋ）についての算出が終了していないときはステップＳ４に戻り、すべての説明変数ｘαｑ＋１（ｋ）についての算出が終了したときは、次のステップＳ８に移行する。
【００９５】
次に、ステップＳ８において、最大寄与率選定手段２２は、追加時寄与率記憶手段２０に格納された寄与率テーブルを参照し、寄与率Ｆ（βαｑ＋１（ｋ））が最大の説明変数ｘαｑ＋１（ｋ）を索出し、これを初期の追加候補説明変数ｘαｑ＋１に設定する。
【００９６】
次に、ステップＳ９において、追加基準判定手段２３は、追加候補説明変数ｘαｑ＋１の寄与率Ｆ（βαｑ＋１）が、重回帰式に取り込む取込基準値Ｆｉｎより大きいか否かを判定する。当該取込基準値Ｆｉｎより大きい場合には、次のステップＳ１０に移行し、当該取込基準値Ｆｉｎより大きくない場合には、変数取込処理を終了して、ステップＳ２３に移行する。
【００９７】
次に、ステップＳ１０において、相関係数演算手段２４は、変数群属性テーブルにおいて変数群属性が未追加変数群Ｐの他の各説明変数ｘｊを１つ選択し、追加候補説明変数ｘαｑ＋１と、未追加変数群Ｐの他の各説明変数ｘｊとの相関係数ｒαｑ＋１，ｊを順次算出する。
【００９８】
次に、ステップＳ１１において、追加判定手段２５は、相関係数演算手段２４が算出する相関係数ｒαｑ＋１，ｊが所定の相関係数閾値ｒｉｎより小さいか否かを判定する。
【００９９】
相関係数閾値ｒｉｎより小さい場合は、ステップＳ１３において、全ての未追加説明変数群Ｐの説明変数ｘｊが選択されたか否かを判定し、まだ、選択されていない説明変数ｘｊがあれば、ステップＳ１０に戻り、すべての説明変数ｘｊが選択されていれば、ステップＳ１３において、変数群属性テーブルにおける当該追加候補説明変数ｘαｑ＋１の変数群属性を追加済変数群Ｑに変更し、変数増加手段１３における変数取込処理を終了し、次の変数除去処理（ステップＳ１６）に移行する。
【０１００】
一方、相関係数閾値ｒｉｎより小さくない場合は、ステップＳ１４において、追加判定手段２５は、追加時寄与率記憶手段２０を参照し、変数群属性が未追加変数群Ｐの説明変数のうち現在の追加候補説明変数の次に寄与率の大きい次候補説明変数があるか否かを検索する。そして、次候補説明変数がなければ、変数増加手段１３における変数取込処理を終了し、ステップＳ２３に移行する。また、次候補説明変数があるときは、現在の追加候補説明変数ｘαｑ＋１を当該次候補説明変数に変更し、ステップＳ９に戻る。
【０１０１】
（変数除去処理）
変数除去処理については、従来のステップワイズ法と同様である。
【０１０２】
まず、ステップＳ１６において、除去時寄与率演算手段３１は、変数群属性記憶手段１１に記憶されている変数群属性テーブルを参照し、変数群属性が追加済変数群Ｑである説明変数ｘαｋを順次１つずつ選択する。
【０１０３】
次に、ステップＳ１７において、除去時寄与率演算手段３１は、当該選択した説明変数（以下「選択説明変数」という。）ｘαｋを除く追加済変数群Ｑの説明変数（ｘα１，…，ｘαｋ−１，ｘαｋ＋１，…，ｘαｑ）に、当該選択した説明変数ｘαｋを追加する場合の寄与率Ｆ（βαｋ）を算出し、除去時寄与率記憶手段３０に寄与率テーブルとして格納する。寄与率Ｆ（βαｋ）の算出は、変数取込処理における場合と同様である。
【０１０４】
次に、ステップＳにおいて、除去時寄与率演算手段３１は、変数群属性が追加済変数群Ｑである説明変数ｘαｋをすべて選択したかを判定し、まだ選択していない説明変数ｘαｋがあれば、ステップＳ１６に戻り、すべての説明変数ｘαｋについて寄与率Ｆ（βαｋ）を算出した場合には、次のステップＳ１９に移行する。
【０１０５】
次に、ステップＳ１９において、最小寄与率選択手段３２は、除去時寄与率記憶手段３０格納された寄与率テーブルを参照し、除去時寄与率演算手段３１により算出される寄与率が最小の説明変数を、除去候補説明変数ｘminに選択する。
【０１０６】
次に、ステップＳ２０において、除去基準判定手段３３は、最小寄与率選択手段３２により選択された除去候補説明変数ｘminの寄与率Ｆminが、重回帰式から除去する除去基準値Ｆｏｕｔより小さいか否かを判定する。
【０１０７】
寄与率Ｆminが除去基準値Ｆｏｕｔより小さい場合、ステップＳ２１において、変数群属性記憶手段１１に格納された変数群属性テーブルにおける当該除去候補説明変数ｘminの変数群属性を未追加変数群Ｐに変更し、追加済変数群Ｑの個数ｑを１だけ減少させ、未追加変数群Ｐの個数ｐを１だけ増加させる。そして、ｑ＝０となった場合には変数減少手段１４における変数除去処理を終了して、次のステップＳ２３に移行する。ｑ＞０の場合にはステップＳ１６に戻って変数除去処理を繰り返し行う。
【０１０８】
一方、寄与率Ｆminが除去基準値Ｆｏｕｔ以上の場合、他に追加すべき説明変数がないか調べるため、再び変数取込処理（ステップＳ３）に戻る。
【０１０９】
最後に、ステップＳ２３において、重回帰モデル出力手段７は、すべての選択済説明変数群Ｐの説明変数による重回帰係数を出力し、終了する。
【０１１０】
以上のように、本実施例の応答曲面モデル構築装置１によれば、ステップＳ１０〜Ｓ１５の処理において、追加候補説明変数ｘαｑ＋１の重回帰式への取り込みを行う際に、未追加変数群Ｐの他の各説明変数ｘｊとの相関係数ｒαｑ＋１，ｊを調べ、これらすべての相関係数ｒαｑ＋１，ｊが相関係数閾値ｒｉｎより小さい場合にのみ追加候補説明変数ｘαｑ＋１の重回帰式への取り込みを行うことで、多重共線性を正確に考慮し、精度のよい応答曲面モデルの構築を行うことができる。
【０１１１】
図６は、実際のエンジンの多項式モデルの構築において、本実施例の応答曲面モデル構築装置１により重回帰モデルを求めた結果の一例である。
【０１１２】
ステップワイズ法を直接適用した場合（図１０）と比較すると、全体的に回帰誤差が縮小し、また、目的変数が負値となるような異常値の出現が防止されていることがわかる。
【実施例２】
【０１１３】
図７は、本発明の実施例２に係る応答曲面モデル構築装置１の構成を表すブロック図である。本実施例では、図１の応答曲面モデル構築装置１の構成に加えて、重回帰モデル評価手段６が新たに追加されている点が異なる。他の構成は、実施例１と同様である。
【０１１４】
図８，図９は、実施例２の応答曲面モデル構築装置１によるエンジン制御パラメータ適合処理を表すフローチャートである。
【０１１５】
図８，図９において、図４，図５と同様の処理を行うステップについては同符号を付して説明は省略する。
【０１１６】
本実施例においては、ステップＳ１の前に、ステップＳ０において、高次説明変数生成
手段３は、高次の説明変数の最大次数ｈｋを１に初期化する。そして、ステップＳ１におおては、高次説明変数生成手段３は、最大次数ｈｋまでの高次の説明変数を算出し標本値テーブルに格納する。尚、ｈｋ＝１の場合には、高次の説明変数はないため、ステップＳ１においては高次の説明変数は算出されないことになる。
【０１１７】
そして、ステップＳ１〜Ｓ２２までの処理は、実施例１と同様の処理が行われる。
【０１１８】
変数取込処理と説明変数除去処理（ステップＳ３〜Ｓ２２）がすべて終了すると、次に、ステップＳ２２ａにおいて、高次説明変数生成手段３は、最大次数ｈｋが予め定められた最終値ｈｍａｘに達したか否かを判定し、達していない場合は、ステップＳ２２ｂにおいてｈｋを１だけ増加させて、ステップＳ１に戻る。最大次数ｈｋが最終値ｈｍａｘに達した場合には、次のステップＳ２２ｃに移行する。尚、最終値ｈｍａｘは、必要となる可能性のある最大次数ｈｋがすべてカバーできる程度の値に適宜設定すればよく、通常は、ｈｍａｘ＝３〜８の値に設定される。
【０１１９】
ステップＳ２２ｃにおいて、重回帰モデル評価手段６は、最大次数ｈｋが１〜ｈｍａｘの場合において算出された各重回帰式の評価・選定を行う。ここでは、根平均自乗誤差（ＲＭＳＥ）による分析処理により評価を行う。
【０１２０】
まず、重回帰モデル評価手段６は、次式により、各最大次数ｈｋにおける根平均自乗誤差ｒｋを算出する。ここで、ｙ＾（ｎ）は、式（１２ｂ）に示した目的変数の推定値である。
【０１２１】
【数１４】
【０１２２】
そして、重回帰モデル評価手段６は、各最大次数ｈｋ（ｋ＝１，…，ｈｍａｘ）の最大値ｒｍａｘを求める。次いで、根平均自乗誤差ｒｋに基づいて、次式により各最大次数ｈｋにおけるコストＪｋを算出する。
【０１２３】
【数１５】
【０１２４】
ここで、ｗｐは累乗の大きさに対する重み係数で、０〜１の予め所定の値に定めておく。また、ｗｒはＲＭＳＥに対する重み係数でありｗｐ＋ｗｒ＝１である。
【０１２５】
そして、重回帰モデル評価手段６は、コストＪｋが最も小さくなる最大次数ｈｋを選択する。
【０１２６】
最後に、ステップＳ２３において、重回帰モデル出力手段７は、重回帰モデル評価手段６が選択した最大次数ｈｋに対する重回帰式の重回帰係数を出力し、終了する。
【０１２７】
このように、本実施例の応答曲面モデル構築装置１によれば、重回帰モデル評価手段６によって、最も適した最大次数ｈｋの重回帰式による重回帰モデルを選定でき、不必要に説明変数が増大することを防止することができる。
【符号の説明】
【０１２８】
応答曲面モデル構築装置１
測定データ記憶手段２
高次説明変数生成手段３
重回帰係数記憶手段４
説明変数選択手段５
重回帰モデル評価手段６
重回帰モデル出力手段７
変数群属性記憶手段１１
変数群属性初期化手段１２
変数増加手段１３
変数減少手段１４
変数選択処理制御手段１５
追加時寄与率記憶手段２０
追加時寄与率演算手段２１
最大寄与率選定手段２２
追加基準判定手段２３
相関係数演算手段２４
追加判定手段２５
追加候補説明変数変更手段２６
変数取込処理制御手段２７
除去時寄与率記憶手段３０
除去時寄与率演算手段３１
最小寄与率選択手段３２
除去基準判定手段３３
変数除去処理制御手段３４
エンジン１０１
制御パラメータ設定手段１０２
制御パラメータ測定手段１０３
特性値測定手段１０４
【特許請求の範囲】
【請求項１】
エンジンの各種制御パラメータである１次の説明変数及びそれに対するエンジンの特性値である目的変数の測定値の組である測定値セットを測定回数だけ格納する測定データ記憶手段と、
各測定値セットについて、各制御パラメータの累乗項を含むｋ個（ｋ＝２，…，Ｍ；Ｍは重回帰式に取り込む説明変数の次数の最大値）の制御パラメータの積項からなる高次の説明変数を算出し、前記測定データ記憶手段に保存する高次説明変数生成手段と、
前記各説明変数から前記目的変数を予測する重回帰式の偏回帰係数及び切片である重回帰係数を記憶する重回帰係数記憶手段と、
前記測定データ記憶手段に格納されたすべての説明変数の中からから重回帰式による目的変数の予測に有効な説明変数を選択するとともに、選択した各説明変数による重回帰式の重回帰係数を算出し前記重回帰係数記憶手段に格納する説明変数選択手段と、を備え、
前記説明変数選択手段は、
各説明変数について、当該説明変数が、重回帰式に取り込まれない予定の未追加変数群に属するのか、重回帰式に取り込む予定の追加変数群に属するのかを示す変数群属性を登録する変数群属性テーブルが格納された変数群属性記憶手段と、
前記変数群属性テーブルの各説明変数の変数群属性を、すべて未追加変数群に初期化する変数群属性初期化手段と、
変数群属性が未追加変数群である説明変数のうち、予測に有効な説明変数の変数群属性を追加変数群に変更する変数取込処理を行う変数増加手段と、
変数群属性が追加変数群である説明変数のうち、予測に有効でない説明変数の変数群属性を未追加変数群に変更する変数除去処理を行う変数減少手段と、
前記変数取込処理及び前記変数除去処理を交互に繰り返し実行し、前記変数群属性テーブルにおいて変数群属性が追加変数群の説明変数が未追加変数群の説明変数よりも少なくなるか、又は前記変数取込処理において何れの説明変数の変数群属性も追加変数群に変更されなかった場合に繰り返し処理を終了する変数選択処理制御手段と、を備え、
前記変数増加手段は、
変数群属性が未追加変数群の説明変数を順次１つずつ選択し、選択した説明変数，変数群属性が追加変数群の説明変数，及び目的変数の測定値セットを前記測定データ記憶手段からすべて読み出し、読み出した測定値セットから重回帰係数を算出し、読み出した測定値セット及び当該重回帰係数から、選択した説明変数の寄与率を算出する追加時寄与率演算手段と、
前記変数群属性テーブルにおいて変数群属性が未追加変数群の説明変数のうち、前記追加時寄与率演算手段により算出された寄与率が最大の説明変数を初期の追加候補説明変数に選択する最大寄与率選定手段と、
前記追加候補説明変数の寄与率が、重回帰式に取り込む取込基準値より大きいか否かを判定し、当該取込基準値より大きくない場合、前記変数取込処理を終了する追加基準判定手段と、
前記追加基準判定手段が、当該追加候補説明変数の寄与率が取込基準値より大きいと判定した場合、当該追加候補説明変数と、前記変数群属性テーブルにおいて変数群属性が未追加変数群の他の各説明変数との相関係数を順次算出する相関係数演算手段と、
前記相関係数演算手段が算出する相関係数がすべて所定の相関係数閾値より小さい場合、前記変数群属性テーブルにおける当該追加候補説明変数の変数群属性を追加変数群に変更し、前記変数取込処理を終了する追加判定手段と、
前記相関係数演算手段が算出する相関係数の何れかが前記相関係数閾値より大きい場合、変数群属性が未追加変数群の説明変数のうち現在の追加候補説明変数の次に寄与率が大きい次候補説明変数がなければ前記変数取込処理を終了し、次候補説明変数があるときは、現在の追加候補説明変数を当該次候補説明変数に変更する追加候補説明変数変更手段と、
前記追加基準判定手段、前記相関係数演算手段、前記追加判定手段、及び前記追加候補説明変数変更手段による処理を繰り返し実行する変数取込処理制御手段と、を備え、
前記変数減少手段は、
前記変数群属性テーブルにおける変数群属性が追加変数群の説明変数を順次１つずつ選択し、当該選択した説明変数を除く追加変数群の説明変数に、当該選択した説明変数を追加する場合の寄与率を算出する除去時寄与率演算手段と、
前記除去時寄与率演算手段により算出される寄与率が最小の説明変数を、除去候補説明変数に選択する最小寄与率選択手段と、
当該除去候補説明変数の寄与率が、重回帰式から除去する除去基準値より小さい場合、前記変数群属性テーブルにおける当該除去候補説明変数の変数群属性を未追加変数群に変更し、それ以外の場合には前記変数除去処理を終了する除去基準判定手段と、
前記最小寄与率選択手段及び前記除去基準判定手段による処理を、前記除去基準判定手段が前記変数除去処理を終了するまで反復実行する制御を行う変数除去処理制御手段と、を備えたことを特徴とする応答曲面モデル構築装置。
【請求項２】
前記寄与率は、Ｆ値又はＴ統計量であることを特徴とする請求項１記載の応答曲面モデル構築装置。
【請求項３】
前記説明変数選択手段は、高次の説明変数の最大次数ｈｋが１からｈｍａｘのそれぞれの場合について重回帰式の重回帰係数を算出し、前記重回帰係数記憶手段に格納するものであり、
前記重回帰係数記憶手段に格納された最大次数ｈｋ（ｋ＝１，…，ｈｍａｘ）の場合の重回帰式のそれぞれに対して、根平均自乗誤差ｒｋを算出し、算出した各根平均自乗誤差ｒｋの最大値ｒｍａｘを抽出するとともに、各根平均自乗誤差ｒｋ及びその最大値ｒｍａｘに基づいて、各最大次数ｈｋにおけるコストＪｋを算出し、コストＪｋが最大となる最大次数ｈｋに対する重回帰式を選択し出力する重回帰モデル評価手段を備えたことを特徴とする請求項１又は２に記載の応答曲面モデル構築装置。
【請求項４】
コンピュータに読み込んで実行させることにより、コンピュータを請求項１乃至３の何れかに記載の応答曲面モデル構築装置として機能させることを特徴とするプログラム。
【請求項１】
エンジンの各種制御パラメータである１次の説明変数及びそれに対するエンジンの特性値である目的変数の測定値の組である測定値セットを測定回数だけ格納する測定データ記憶手段と、
各測定値セットについて、各制御パラメータの累乗項を含むｋ個（ｋ＝２，…，Ｍ；Ｍは重回帰式に取り込む説明変数の次数の最大値）の制御パラメータの積項からなる高次の説明変数を算出し、前記測定データ記憶手段に保存する高次説明変数生成手段と、
前記各説明変数から前記目的変数を予測する重回帰式の偏回帰係数及び切片である重回帰係数を記憶する重回帰係数記憶手段と、
前記測定データ記憶手段に格納されたすべての説明変数の中からから重回帰式による目的変数の予測に有効な説明変数を選択するとともに、選択した各説明変数による重回帰式の重回帰係数を算出し前記重回帰係数記憶手段に格納する説明変数選択手段と、を備え、
前記説明変数選択手段は、
各説明変数について、当該説明変数が、重回帰式に取り込まれない予定の未追加変数群に属するのか、重回帰式に取り込む予定の追加変数群に属するのかを示す変数群属性を登録する変数群属性テーブルが格納された変数群属性記憶手段と、
前記変数群属性テーブルの各説明変数の変数群属性を、すべて未追加変数群に初期化する変数群属性初期化手段と、
変数群属性が未追加変数群である説明変数のうち、予測に有効な説明変数の変数群属性を追加変数群に変更する変数取込処理を行う変数増加手段と、
変数群属性が追加変数群である説明変数のうち、予測に有効でない説明変数の変数群属性を未追加変数群に変更する変数除去処理を行う変数減少手段と、
前記変数取込処理及び前記変数除去処理を交互に繰り返し実行し、前記変数群属性テーブルにおいて変数群属性が追加変数群の説明変数が未追加変数群の説明変数よりも少なくなるか、又は前記変数取込処理において何れの説明変数の変数群属性も追加変数群に変更されなかった場合に繰り返し処理を終了する変数選択処理制御手段と、を備え、
前記変数増加手段は、
変数群属性が未追加変数群の説明変数を順次１つずつ選択し、選択した説明変数，変数群属性が追加変数群の説明変数，及び目的変数の測定値セットを前記測定データ記憶手段からすべて読み出し、読み出した測定値セットから重回帰係数を算出し、読み出した測定値セット及び当該重回帰係数から、選択した説明変数の寄与率を算出する追加時寄与率演算手段と、
前記変数群属性テーブルにおいて変数群属性が未追加変数群の説明変数のうち、前記追加時寄与率演算手段により算出された寄与率が最大の説明変数を初期の追加候補説明変数に選択する最大寄与率選定手段と、
前記追加候補説明変数の寄与率が、重回帰式に取り込む取込基準値より大きいか否かを判定し、当該取込基準値より大きくない場合、前記変数取込処理を終了する追加基準判定手段と、
前記追加基準判定手段が、当該追加候補説明変数の寄与率が取込基準値より大きいと判定した場合、当該追加候補説明変数と、前記変数群属性テーブルにおいて変数群属性が未追加変数群の他の各説明変数との相関係数を順次算出する相関係数演算手段と、
前記相関係数演算手段が算出する相関係数がすべて所定の相関係数閾値より小さい場合、前記変数群属性テーブルにおける当該追加候補説明変数の変数群属性を追加変数群に変更し、前記変数取込処理を終了する追加判定手段と、
前記相関係数演算手段が算出する相関係数の何れかが前記相関係数閾値より大きい場合、変数群属性が未追加変数群の説明変数のうち現在の追加候補説明変数の次に寄与率が大きい次候補説明変数がなければ前記変数取込処理を終了し、次候補説明変数があるときは、現在の追加候補説明変数を当該次候補説明変数に変更する追加候補説明変数変更手段と、
前記追加基準判定手段、前記相関係数演算手段、前記追加判定手段、及び前記追加候補説明変数変更手段による処理を繰り返し実行する変数取込処理制御手段と、を備え、
前記変数減少手段は、
前記変数群属性テーブルにおける変数群属性が追加変数群の説明変数を順次１つずつ選択し、当該選択した説明変数を除く追加変数群の説明変数に、当該選択した説明変数を追加する場合の寄与率を算出する除去時寄与率演算手段と、
前記除去時寄与率演算手段により算出される寄与率が最小の説明変数を、除去候補説明変数に選択する最小寄与率選択手段と、
当該除去候補説明変数の寄与率が、重回帰式から除去する除去基準値より小さい場合、前記変数群属性テーブルにおける当該除去候補説明変数の変数群属性を未追加変数群に変更し、それ以外の場合には前記変数除去処理を終了する除去基準判定手段と、
前記最小寄与率選択手段及び前記除去基準判定手段による処理を、前記除去基準判定手段が前記変数除去処理を終了するまで反復実行する制御を行う変数除去処理制御手段と、を備えたことを特徴とする応答曲面モデル構築装置。
【請求項２】
前記寄与率は、Ｆ値又はＴ統計量であることを特徴とする請求項１記載の応答曲面モデル構築装置。
【請求項３】
前記説明変数選択手段は、高次の説明変数の最大次数ｈｋが１からｈｍａｘのそれぞれの場合について重回帰式の重回帰係数を算出し、前記重回帰係数記憶手段に格納するものであり、
前記重回帰係数記憶手段に格納された最大次数ｈｋ（ｋ＝１，…，ｈｍａｘ）の場合の重回帰式のそれぞれに対して、根平均自乗誤差ｒｋを算出し、算出した各根平均自乗誤差ｒｋの最大値ｒｍａｘを抽出するとともに、各根平均自乗誤差ｒｋ及びその最大値ｒｍａｘに基づいて、各最大次数ｈｋにおけるコストＪｋを算出し、コストＪｋが最大となる最大次数ｈｋに対する重回帰式を選択し出力する重回帰モデル評価手段を備えたことを特徴とする請求項１又は２に記載の応答曲面モデル構築装置。
【請求項４】
コンピュータに読み込んで実行させることにより、コンピュータを請求項１乃至３の何れかに記載の応答曲面モデル構築装置として機能させることを特徴とするプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【公開番号】特開２０１２−９２７６０（Ｐ２０１２−９２７６０Ａ）
【公開日】平成２４年５月１７日（２０１２．５．１７）
【国際特許分類】
【出願番号】特願２０１０−２４１５７０（Ｐ２０１０−２４１５７０）
【出願日】平成２２年１０月２７日（２０１０．１０．２７）
【出願人】（８９９００００６８）学校法人早稲田大学 (602)
【Ｆターム（参考）】
【公開日】平成２４年５月１７日（２０１２．５．１７）
【国際特許分類】
【出願日】平成２２年１０月２７日（２０１０．１０．２７）
【出願人】（８９９００００６８）学校法人早稲田大学 (602)
【Ｆターム（参考）】
[ Back to top ]