(57)【要約】【課題】モデル予測制御において制御対象の状態制約を適切に取り扱う。【解決手段】制御対象の状態制約を導入しない場合に用意するデータに加え、モデル予測制御の制御対象における状態制約に対応するペナルティ関数を、制約条件において制限されるパラメータについて１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、上記データに含まれる一次式群から線形化制御入力についての係数を抽出し、ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する。【選択図】図４
2015095126.tif 000002
【特許請求の範囲】
【請求項１】
（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ前記制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群と、（ｄ）前記モデル予測制御の制御対象についての制約条件に対応するペナルティ関数と、（ｅ）前記制約条件において制限されるパラメータに関する関係式とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形し、
前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む第５の一次式群に変形し、
前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を代入することで、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形し、
前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形し、
前記データ格納部に格納された前記ペナルティ関数を前記パラメータについて１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、
前記データ格納部から前記関係式を読み出して、前記第４の一次式群を前記関係式に代入することで、前記関係式を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第７の一次式群に変形し、
前記第７の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する
処理を、コンピュータに実行させるためのプログラム。
【請求項２】
前記パラメータが、前記状態偏差と前記線形化制御入力とから表されるパラメータである
請求項１記載のプログラム。
【請求項３】
前記パラメータが、前記状態偏差の変化量である
請求項１記載のプログラム。
【請求項４】
前記ペナルティ関数が、前記パラメータの値域毎に異なる関数で表され、
前記第１及び第２の微分関数を、前記パラメータの値域毎に生成する
請求項１乃至３のいずれか１つ記載のプログラム。
【請求項５】
（ａ）モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、（ｂ）制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、（ｃ）前記モデル予測制御の制御対象における制約条件に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、（ｄ）前記制約条件において制限されるパラメータについて前記ペナルティ関数を１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び前記第２の式を読み出して、前記状態偏差の現在値から得られる、前記パラメータの基点の値を代入することで、第１の値及び第２の値を算出し、
前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出し、
前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルから、（（前記第１の値）−（前記第２の値と前記パラメータの基点に関連し且つ前記状態偏差の現在値又は前記線形化制御入力の１単位時間前の値から得られる第３の値との積））と前記係数ベクトルとの積を差し引いた第２のベクトルを算出し、
前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成し、
前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
処理をプロセッサに実行させるためのプログラム。
【請求項６】
前記状態偏差の現在値から算出される前記パラメータの現在値が所定範囲内であるか判断する処理をさらに含み、
前記パラメータの現在値が前記所定範囲内であれば、前記第１の値及び第２の値を算出する処理と、前記第１の行列を算出する処理と、前記第２のベクトルを算出する処理と、前記線形化制御入力のベクトルを生成する処理と、前記線形化制御入力を抽出する処理とを実行する
請求項５記載のプログラム。
【請求項７】
前記データ格納部が、
前記パラメータの現在値の値域に応じた前記第１の式及び前記第２の式を格納しており、
前記第１の値及び前記第２の値を算出する処理が、
前記パラメータの現在値が属する値域を特定し、当該特定された値域に応じた前記第１の式及び前記第２の式を読み出す処理
を含む請求項５又は６記載のプログラム。
【請求項８】
前記パラメータの現在値が前記所定範囲外であれば、
前記データ格納部から前記関数ベクトル及び前記第１の係数行列を読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルと、前記第１の係数行列の逆行列との積により、前記線形化制御入力の第２のベクトルを生成し、
前記線形化制御入力の第２のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
処理をさらに含む請求項５乃至７のいずれか１つ記載のプログラム。
【請求項９】
前記パラメータが、前記状態偏差と前記線形化制御入力とから表されるパラメータである場合には、
前記第３の値は、前記制御対象の状態方程式で用いられる行列と前記線形化制御入力の１単位時間前の値とから算出される
請求項５乃至８のいずれか１つ記載のプログラム。
【請求項１０】
前記パラメータが、前記状態偏差の変化量である場合には、
前記第３の値は、前記状態偏差の現在値と前記状態偏差の１単位時間前の値との差である
請求項１記載のプログラム。
【請求項１１】
（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ前記制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群と、（ｄ）前記モデル予測制御の制御対象についての制約条件に対応するペナルティ関数と、（ｅ）前記制約条件において制限されるパラメータに関する関係式とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形し、
前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む第５の一次式群に変形し、
前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を代入することで、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形し、
前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形し、
前記データ格納部に格納された前記ペナルティ関数を前記パラメータについて１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、
前記データ格納部から前記関係式を読み出して、前記第４の一次式群を前記関係式に代入することで、前記関係式を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第７の一次式群に変形し、
前記第７の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する
処理を含み、コンピュータにより実行されるオフライン処理方法。
【請求項１２】
（ａ）モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、（ｂ）制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、（ｃ）前記モデル予測制御の制御対象における制約条件に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、（ｄ）前記制約条件において制限されるパラメータについて前記ペナルティ関数を１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び前記第２の式を読み出して、前記状態偏差の現在値から得られる、前記パラメータの基点の値を代入することで、第１の値及び第２の値を算出し、
前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出し、
前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルから、（（前記第１の値）−（前記第２の値と前記パラメータの基点に関連し且つ前記状態偏差の現在値又は前記線形化制御入力の１単位時間前の値から得られる第３の値との積））と前記係数ベクトルとの積を差し引いた第２のベクトルを算出し、
前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成し、
前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
処理を含み、プロセッサにより実行される制御方法。
【請求項１３】
（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ前記制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群と、（ｄ）前記モデル予測制御の制御対象についての制約条件に対応するペナルティ関数と、（ｅ）前記制約条件において制限されるパラメータに関する関係式とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形する手段と、
前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む第５の一次式群に変形する手段と、
前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を代入することで、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形する手段と、
前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形する手段と、
前記データ格納部に格納された前記ペナルティ関数を前記パラメータについて１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成する手段と、
前記データ格納部から前記関係式を読み出して、前記第４の一次式群を前記関係式に代入することで、前記関係式を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第７の一次式群に変形する手段と、
前記第７の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する手段と、
を有するオフライン処理装置。
【請求項１４】
（ａ）モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、（ｂ）制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、（ｃ）前記モデル予測制御の制御対象における制約条件に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、（ｄ）前記制約条件において制限されるパラメータについて前記ペナルティ関数を１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び前記第２の式を読み出して、前記状態偏差の現在値から得られる、前記パラメータの基点の値を代入することで、第１の値及び第２の値を算出する手段と、
前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出する手段と、
前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルから、（（前記第１の値）−（前記第２の値と前記パラメータの基点に関連し且つ前記状態偏差の現在値又は前記線形化制御入力の１単位時間前の値から得られる第３の値との積））と前記係数ベクトルとの積を差し引いた第２のベクトルを算出する手段と、
前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成する手段と、
前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する手段と、
を有する制御装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本技術は、モデル予測制御に関する。
【背景技術】
【０００２】
モデル予測制御は、各時刻において有限時間未来までの応答を最適化することによって制御対象への入力を決定する制御方法である。操作量飽和などの不等式制約条件を陽に扱うことができる制御手法として高い制御性能を期待されている。
【０００３】
このモデル予測制御の適用については様々な分野について考察されているが、ディーゼルエンジンの吸気系制御への適用も考察されている。ディーゼルエンジンの吸気系は、ＥＧＲ（Exhaust Gas Recirculator：排気循環器）とＶＮＴ（Variable Nozzle Turbo：可変ノズルターボ）を用いて、ＭＡＦ（Mass Air Flow：新気量）とＭＡＰ（Manifold Absolute Pressure：吸気圧）を独立に、ＳＩＳＯ（Single Input Single Output）としてＰＩＤ（Proportional, Integral and Differential）制御を行っていることが多い。しかしながら、ＥＧＲとＶＮＴは排気を共有して互いに干渉することから、現行の手法では低環境負荷性能を達成するために十分な目標追従性を確保することが難しくなっている。
【０００４】
このため、ＥＧＲ及びＶＮＴでＭＡＦ及びＭＡＰを制御する２入力２出力のＭＩＭＯ（Multi-Input and Multi-Output）制御が検討されているが、ここでＥＧＲバルブとＶＮＴノズルの操作量の飽和に起因する目標値追従性能の劣化が大きな問題となっている。この問題の解決策として、上で述べたモデル予測制御が期待されている。
【０００５】
しかしながら、吸気系のような応答の速い対象に対してモデル予測制御を適用する場合、大きな課題となるのが、サンプル時間間隔毎に行う最適化計算に要する計算時間である。この問題に対しては様々な取り組みがなされているが、その中で使われる手法が、操作量（制御入力とも呼ぶ）飽和条件付きの最適化問題をペナルティ関数を用いて緩和する方法である。
【０００６】
一方、エンジンなどの制御対象となるプラントに状態制約が存在する場合については、あまり検討がなされていない。そして、どのようにすれば実時間処理における計算時間を抑えることができるかという点については、具体的に考察されていない。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】梅田裕平，丸山次人，志村知洋，穴井宏和，”代数的簡略化を用いたモデル予測制御の効率的解法”，第１１回計測自動制御学会制御部門大会予稿集，１７１−１−３（２０１１）
【非特許文献２】Y. Umeda, T. Maruyama, H. Anai, A. Ejiri and K. Shimotani,"A Fast Model Predictive Control Algorithm for Diesel Engines", Proceedings of 4th IFAC Nonlinear Model Predictive Control Conference International Federation of Automatic Control., Paper SuDR.6, 2012.
【非特許文献３】梅田裕平, 丸山次人, 志村知洋，穴井宏和，下谷 圭司,”ディーゼルエンジンの効率的なモデル予測制御のためのペナルティ関数の設定に関する一考察”，  第１２回計測自動制御学会制御部門大会予稿集, P0066 （２０１２）
【非特許文献４】梅田裕平, 丸山次人, 志村知洋，穴井宏和，下谷 圭司, ”ＧＡを用いたディーゼルエンジンの効率的なモデル予測制御のためのペナルティ関数の設定”，第３回進化計算学会研究会予稿集, P1-5 （２０１２）
【発明の概要】
【発明が解決しようとする課題】
【０００８】
従って、本技術の目的は、一側面として、モデル予測制御において制御対象の状態制約を適切に取り扱うための技術を提供することである。
【課題を解決するための手段】
【０００９】
本発明の第１の態様に係るオフライン処理方法は、（１）（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ制御対象の状態方程式に関する第１の関数との和を含む第２の関数をラグランジュ乗数について偏微分することで得られ且つ制御対象の状態偏差及び制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）第２の関数を制御対象の状態偏差について偏微分することで得られ且つ制御対象の状態偏差及びラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）第２の関数を制御対象の線形化制御入力について偏微分することで得られ且つラグランジュ乗数及び制御対象の状態偏差を変数として含む第３の一次式群と、（ｄ）モデル予測制御の制御対象についての制約条件に対応するペナルティ関数と、（ｅ）制約条件において制限されるパラメータに関する関係式とを格納するデータ格納部から第１の一次式群を読み出して、状態偏差の現在値を除き状態偏差の変数を消し込むように順次代入を行って、第１の一次式群を、状態偏差の現在値と線形化制御入力とを変数として含む第４の一次式群に変形し、（２）データ格納部から第２の一次式群を読み出して、ラグランジュ乗数の変数を消し込むように順次代入を行って、第２の一次式群を、状態偏差を変数として含む第５の一次式群に変形し、（３）データ格納部から第３の一次式群を読み出して、第４の一次式群及び第５の一次式群を代入することで、第３の一次式群を、状態偏差の現在値と線形化制御入力とを変数として含む第６の一次式群に変形し、（４）第６の一次式群を、第１の係数行列と線形化制御入力についてのベクトルとの積が、状態偏差の現在値の関数ベクトルと等式で表されるように変形し、（５）データ格納部に格納されたペナルティ関数をパラメータについて１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、（６）データ格納部から関係式を読み出して、第４の一次式群を関係式に代入することで、関係式を、状態偏差の現在値と線形化制御入力とを変数として含む第７の一次式群に変形し、（７）第７の一次式群から線形化制御入力についての係数を抽出し、ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する処理を含む。
【００１０】
本発明の第２の態様に係る制御方法は、（１）（ａ）モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、（ｂ）制御期間内における制御対象に対する線形化制御入力のベクトルとの積が関数ベクトルと等しくなるように予め算出されている第１の係数行列と、（ｃ）モデル予測制御の制御対象における制約条件に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、（ｄ）制約条件において制限されるパラメータについてペナルティ関数を１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から第１の式及び第２の式を読み出して、状態偏差の現在値から得られる、パラメータの基点の値を代入することで、第１の値及び第２の値を算出し、（２）データ格納部から第１の係数行列及び第２の係数行列を読み出して、第２の値と第２の係数行列との積と第１の係数行列との和である第１の行列を算出し、（３）データ格納部から関数ベクトルを読み出して、関数ベクトルに状態偏差の現在値を代入することで得られる第１のベクトルから、（（第１の値）−（第２の値とパラメータの基点に関連し且つ状態偏差の現在値又は線形化制御入力の１単位時間前の値から得られる第３の値との積））と係数ベクトルとの積を差し引いた第２のベクトルを算出し、（４）第１の行列の逆行列と第２のベクトルとの積により、線形化制御入力のベクトルを生成し、（５）線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する処理を含む。
【発明の効果】
【００１１】
一側面によれば、モデル予測制御において制御対象の状態制約をペナルティ関数で適切に取り扱うことができるようになる。
【図面の簡単な説明】
【００１２】
【図１】図１は、ディーゼルエンジンの吸気システムの一例を示す図である。
【図２】図２は、エンジン制御装置のブロック線図の一例を示す図である。
【図３】図３は、オフライン処理装置の機能ブロック図である。
【図４】図４は、オフライン処理装置により実行される処理の処理フローを示す図である。
【図５】図５は、モデル予測制御部の機能ブロック図を示す図である。
【図６】図６は、モデル予測制御部により実行される処理の処理フローを示す図である。
【図７】図７は、リチウムイオン電池の等価回路を示す図である。
【図８】図８は、第３の実施の形態におけるペナルティ関数の第１の例を示す図である。
【図９】図９は、第３の実施の形態におけるペナルティ関数の第２の例を示す図である。
【図１０】図１０は、第３の実施の形態におけるモデル予測制御部により実行される処理の処理フローを示す図である。
【図１１】図１１は、第４の実施の形態におけるペナルティ関数の一例を示す図である。
【図１２】図１２は、コンピュータの機能ブロック図である。
【図１３】図１３は、制御装置の機能ブロック図である。
【発明を実施するための形態】
【００１３】
［実施の形態１］
図１にディーゼルエンジンの吸気システムを示す。ディーゼルエンジンの吸気制御系は、吸気圧制御系と新気量制御系を含む。吸気圧制御系は、排気中の微粒子物質（ＰＭ：Particulate Matter）を低減するために、可変ノズルターボＶＮＴのノズル径を制御して吸気圧を吸気圧目標値に追従するようにコントロールしている。一方、新気量制御系は、排気中の窒素酸化物（ＮＯｘ）を低減するために、排気をシリンダ内に再循環させる排気循環器ＥＧＲのバルブ開度を制御して新気量を新気量目標値に追従するようにコントロールしている。
【００１４】
すなわち、エンジン本体１には、エンジン本体１からの排ガスを供給する排気循環器ＥＧＲと、排ガスの圧力にてタービンを回して新気（Fresh Air）を圧縮してエンジン本体１に供給する可変ノズルターボＶＮＴとが接続されている。可変ノズルターボＶＮＴのノズル開度を調整することによって、可変ノズルターボＶＮＴのタービンの回転が調整され、吸気圧（ＭＡＰ）センサで測定される吸気圧（ＭＡＰ）が調整される。一方、排気循環器ＥＧＲに設けられているＥＧＲバルブのバルブ開度を調整することによって、新気量（ＭＡＦ）センサで測定される新気量（ＭＡＦ）が調整される。
【００１５】
本実施の形態に係るエンジン制御装置１０００には、ＭＡＰセンサからの吸気圧測定値と、ＭＡＦセンサからの新気量測定値と、外部から与えられる燃料噴射量の設定値と、同じく外部から与えられるエンジン回転数の設定値とが入力されるようになっている。また、エンジン制御装置１０００からは、ＥＧＲバルブのバルブ開度の操作量がＥＧＲバルブに出力され、ＶＮＴノズルのノズル開度の操作量がＶＮＴノズルに出力されるようになっている。
【００１６】
以下、エンジンをプラントの一例として説明するが、本実施の形態の適用範囲は、エンジンに限られない。
【００１７】
本実施の形態に係るエンジン制御装置１０００のブロック線図を図２に示す。エンジン制御装置１０００は、計画器２１０と、モデル予測制御部２２０と、変換部２２５とを有する。プラント２４０は、例えばモデル化されたエンジン特性に相当し、ＥＧＲバルブ及びＶＮＴノズルの開度には物理的に制限があるため、飽和要素２３０も設けられている。すなわち、プラント２４０への操作量は、飽和要素２３０による制限を受ける。
【００１８】
このようなエンジン制御装置１０００の計画器２１０には、燃料噴射量の設定値とエンジン回転数の設定値とが入力されて、燃料噴射量の値及びエンジン回転数の値に対応付けてＥＧＲバルブの状態目標値及びＶＮＴノズルの状態目標値Ｖref、燃料噴射量の値及びエンジン回転数の値に対応付けて吸気圧ＭＡＰの目標値及び新気量ＭＡＦの目標値Ｘrefが登録されている。本実施の形態では、計画器２１０は、ＥＧＲバルブの状態目標値及びＶＮＴノズルの状態目標値Ｖref、及び吸気圧ＭＡＰの目標値及び新気量ＭＡＦの目標値Ｘrefを出力する。
【００１９】
この目標値Ｘrefとプラント２４０からの出力ｙ（例えばＭＡＦセンサ出力及びＭＡＰセンサ出力）との差が、モデル予測制御部２２０に入力される。モデル予測制御部２２０は、以下で述べるモデル予測制御の処理を実施して、線形化制御入力ｕを出力する。この線形化制御入力ｕは、以下で説明する変換部２２５（ｇ-1）で、ＥＧＲバルブの状態及びＶＮＴノズルの状態ｖに変換される。そして、状態ｖと、計画器２１０からの状態目標値Ｖrefとは加算されて操作量が生成され、当該操作量が飽和要素２３０を介してプラント２４０に入力される。このような処理が繰り返される。
【００２０】
このようなエンジン本体１の吸気系システムは、ＥＧＲバルブ及びＶＮＴノズルについて非線形であり、エンジンのモデルとしても、エンジン回転数ω及び燃料噴射量σに応じて非線形に変化する。
【００２１】
ここで、エンジン回転数ω及び燃料噴射量σが固定であるとして、エンジンモデルは、以下のように定式化される。
【００２２】
【数１】
2015095126.tif 000003
【００２３】
ここでＡω(k),σ(k)及びＢω(k),σ(k)は、時刻ｋにおける運転条件（ω（ｋ），σ（ｋ））によって変化する非線形の２×２行列関数である。なお、時刻ｋにおける回転数をω（ｋ）と表し、燃料噴射量をσ（ｋ）と表す。また、時刻ｋにおけるＭＡＦの状態偏差をｘmaf（ｋ）と表し、ＭＡＰの状態偏差をｘmap（ｋ）と表す。また、時刻ｋにおけるＭＡＦの測定値をｙmaf（ｋ）と表し、ＭＡＰの測定値をｙmap（ｋ）と表す。さらに、時刻ｋにおけるＥＧＲバルブの状態をｖegr（ｋ）と表し、ＶＮＴノズルの状態をｖvnt（ｋ）と表す。また、吸気圧ＭＡＰの目標値をＸmaprefと表し、新気量ＭＡＦの目標値をＸmafrefと表す。
【００２４】
従って、状態偏差ベクトルｘ（ｋ）、測定値ベクトルｙ（ｋ）、目標値ベクトルＸref（ｋ）及び状態ベクトルｖ（ｋ）は、以下のように表される。
【００２５】
【数２】
2015095126.tif 000004
【００２６】
ｘ（ｋ）、ｙ（ｋ）及びＸref（ｋ）の関係は以下のようになっている。
【００２７】
【数３】
2015095126.tif 000005
【００２８】
また、非線形関数ｇを以下のように表す。
【００２９】
【数４】
2015095126.tif 000006
【００３０】
ｇ1及びｇ2は、単調増加関数である。実験データからすると、このような仮定は妥当性がある。非線型モデルでは計算速度等から取り扱いが困難なため、以下では、非線形なエンジンモデルを線型モデルに変換する。そこで、上で述べた線形化制御入力ｕegr及びｕvntを導入することで、線形化する。
【００３１】
【数５】
2015095126.tif 000007
【００３２】
なお、線形化制御入力ベクトルｕ（ｋ）を以下のように表す。
【００３３】
【数６】
2015095126.tif 000008
【００３４】
本実施の形態では、典型的な運転条件を選択して、セットポイントと呼ぶことにする。そして、モデル予測制御部２２０は、運転条件（ω（ｋ），σ（ｋ））を取得すると、対応するセットポイントについてのω及びσについての行列Ａω,σ及びＢω,σを選択し、以下で述べる処理を行って線形化制御入力ｕを算出する。また、変換関数２２５は、線形化制御入力ｕを、変換関数ｇ-1により非線形の最適状態ｖに変換する。この最適状態ｖと状態目標値Ｖrefとは加算されて操作量が生成され、エンジンのプラント２４０に出力される。
【００３５】
モデル予測制御部２２０は、セットポイントについて最適な線形化制御入力ｕを、１サンプル時間毎に、以下の線形モデルに従って算出する。
【００３６】
【数７】
2015095126.tif 000009
【００３７】
但し、線形化制御入力ｕについての制約条件は以下のとおりである。
【００３８】
【数８】
2015095126.tif 000010
【００３９】
例えば、ω＝２０００ｒｐｍ、σ＝２０ｍｍ3／ｓｔの場合、以下のような行列が用いられる。
【００４０】
【数９】
2015095126.tif 000011
【００４１】
【数１０】
2015095126.tif 000012
【００４２】
各時刻ｋについての線形化制御入力ｕは、以下に示すような評価関数Ｊを最小化することで決定される。なお、予測ホライズンＨp及び制御ホライズンＨcは、予め設定されるものとする。
【００４３】
【数１１】
2015095126.tif 000013
【００４４】
ここで、ＱTは行列Ｑの転置行列である。また、Ｑ∈Ｍ（ｎ，Ｒ）、Ｒ∈Ｍ（ｎ，Ｒ）及びＳ∈Ｍ（ｎ，Ｒ）は、重み行列であり、予め設定される。なお、Ｍ（ｎ，Ｒ）はｎ×ｎ行列を表す。（９）式において、無限ホライズンについての最適化問題を、有限ホライズンについての最適化問題にしており、このため、Ricatti代数方程式を解く行列として重み行列Ｓを用いている。
【００４５】
本実施の形態では、モデル予測制御の問題は、有限ホライズンの最適制御問題を解くことに変形される。
【００４６】
【数１２】
2015095126.tif 000014
【００４７】
ここで、ｘは、現在の状態偏差ベクトルを表す。
【００４８】
最適解を求める手法として、ラグランジュの未定乗数法が効率的であるが、不等式制約がある場合にはこの手法をそのまま用いることができない。従って、不等式制約の代わりにペナルティ関数を導入して、不等式制約のない評価関数に緩和する。
【００４９】
すなわち、以下のような評価関数Ｊbに緩和する。
【００５０】
【数１３】
2015095126.tif 000015
【００５１】
ここで、ｐ（ｕ（ｋ＋ｉ））はペナルティ関数であり、入力値ｕ（ｋ＋ｉ）が、可能な値域の境界に近づくと、大きな値を評価関数Ｊに加えることになる。
【００５２】
一般的に、不等式制約ｕ≦ｕmaxに対して用いられるペナルティ関数としては以下のようなものがある。
【００５３】
【数１４】
2015095126.tif 000016
【００５４】
ｒは、ペナルティ関数の係数であり、不等式制約に対する重みを表している。
【００５５】
本実施の形態では、以下の演算のためには、（ｃ）のようなシンプルなペナルティ関数が好ましいので、（ｃ）を用いることにする。
【００５６】
ディーゼルエンジンの場合、線形化制御入力ｕには以下のような制約が存在する。
【００５７】
【数１５】
2015095126.tif 000017
【００５８】
（ｃ）のような形のペナルティ関数を採用する場合には、ペナルティ関数は以下のように表される。
【００５９】
【数１６】
2015095126.tif 000018
【００６０】
（１５）式を用いて制約条件を緩和する場合、以下に示すような、有限ホライズンについての最適制御問題を解くことに帰着する。
【００６１】
【数１７】
2015095126.tif 000019
【００６２】
このような最適制御問題を解くためには、ラグランジュの未定乗数法が用いられる。ここで、ラグランジュ乗数λ（ｊ）＝［λmaf（ｊ），λmap（ｊ）］Tを用い、ハミルトニアンＨを以下のように表す。
【００６３】
【数１８】
2015095126.tif 000020
【００６４】
ここで、Х（χの大文字）は、以下のように定義される。
【００６５】
【数１９】
2015095126.tif 000021
【００６６】
∇x,u,λＨ＝０と置けば、以下のような方程式が得られる。
【００６７】
【数２０】
2015095126.tif 000022
【００６８】
なお、（２０）式において∂Ｊb／∂ｘ＝∂Ｊ／∂ｘ＋∂ｐ（ｕ）／∂ｘ＝∂Ｊ／∂ｘと変形されている。
【００６９】
なお、∂／∂ｘ，∂／∂ｕ，∂／∂λは、以下のような演算である。
【００７０】
【数２１】
2015095126.tif 000023
【００７１】
なお、方程式（２０）乃至（２２）は、∂ｐ／∂ｕの部分を除き線形方程式となっている。∂ｐ／∂ｕを、直前のステップの最適線形化制御入力ｕ~（”~”はｕの上に付される）でテイラー展開することで線形化すると、以下のように表される。
【００７２】
【数２２】
2015095126.tif 000024
【００７３】
（２３）式の第１項のＥＧＲ成分は、以下のように表される。なお、ＶＮＴ成分も同様の形で表される。
【００７４】
【数２３】
2015095126.tif 000025
【００７５】
このような線形近似を行えば、（２０）乃至（２２）式は線形となる。そうすると、線形近似された（２０）乃至（２２）式は以下のように表される。
【００７６】
【数２４】
2015095126.tif 000026
【００７７】
ここでＭは、（２Ｈp＋Ｈc）×（２Ｈp＋Ｈc）の行列であり、Ｕ及びｂは、（２Ｈp＋Ｈc）ベクトルである。ベクトルＵは、以下のように表される。
【００７８】
【数２５】
2015095126.tif 000027
【００７９】
また、ベクトルｂは、以下のように表される。なお、ｂmaf、ｂmap及びｂconstは、行列の変換によって算出される（２Ｈp＋Ｈc）行列である。
【００８０】
ｂ＝ｘmafｂmaf＋ｘmapｂmap＋ｂconst    （２７）
【００８１】
最適な線形化制御入力は、この線形方程式を解くことによって得られる。しかしながら、行列Ｍのサイズが大きければ計算時間は長くなってしまうので、行列Ｍのサイズをオフラインの前処理にて削減する。
【００８２】
まず、（２２）式の第１乃至第３の方程式を具体的に書き下すと以下のようになる。
【００８３】
【数２６】
2015095126.tif 000028
【００８４】
（２８−１）式及び（２８−２）式の状態偏差ｘmaf（ｋ）及びｘmap（ｋ）は、現在値である。従って、（２８−１）式によれば、ｘmaf（ｋ＋１）は線形化制御入力ｕ（ｋ）の１次式で表される。（２８−２）式によれば、ｘmap（ｋ＋１）も、線形化制御入力ｕ（ｋ）の１次式で表される。また、（２９）式によれば、ｘmaf（ｋ＋２）は、状態偏差及び線形化制御入力の１次式で表される。しかしながら、（２８−１）式及び（２８−２）式を（２９）式に代入すると、ｘmaf（ｋ＋２）を、線形化制御入力ｕ（ｋ）及びｕ（ｋ＋１）で表すことができるようになる。これをｉを増加させつつ繰り返せば、状態偏差ｘmaf（ｋ＋ｉ）及びｘmap（ｋ＋ｉ）は、線形化制御入力の１次式として以下のように表される。
【００８５】
【数２７】
2015095126.tif 000029
【００８６】
一方、（２０）式の（２Ｈp−１）番目及び（２Ｈp−３）番目の方程式を具体的に書き下すと以下のようになる。
【００８７】
【数２８】
2015095126.tif 000030
【００８８】
ここで、ｓ11、ｓ22、ｒ11は、重み行列Ｓ及びＲの行列成分を表す。
【００８９】
この（３１）式及び（３２）式は、λmaf（Ｈp）及びλmap（Ｈp）の一次式で表されることが分かる。そこで、（３１）式及び（３２）式を、（３３）式に代入すると、λmaf（Ｈp−１）が、目標値Ｘrefmaf、状態偏差ｘ（ｋ＋Ｈp）及びｘ（ｋ＋Ｈp−１）で表すことができる。このような処理をｉを減少させつつ順次行うことで、ラグランジュ乗数λ（ｋ＋ｉ）は、以下のように表される。
【００９０】
【数２９】
2015095126.tif 000031
【００９１】
さらに、（３０）式を（３４）式に代入することで、ラグランジュ乗数λ（ｋ＋ｉ）は、線形化制御入力ｕ（ｋ）の一次式に変形される。
【００９２】
【数３０】
2015095126.tif 000032
【００９３】
最後に、（３０）式及び（３５）式を、（２１）式に代入すると、以下のように表される。
【００９４】
【数３１】
2015095126.tif 000033
【００９５】
このように、（２１）式は、状態偏差の初期値ｘ（ｋ）及び線形化制御入力ｕの一次式で表すことができるようになる。そうすると、（３６）式は、以下の形式に直すことができる。
【００９６】
【数３２】
2015095126.tif 000034
【００９７】
ここでＭ2はＨc×Ｈc行列である。Ｕ2及びｂ2は、以下のような形のベクトルである。
【００９８】
【数３３】
2015095126.tif 000035
【００９９】
【数３４】
2015095126.tif 000036
【０１００】
なお、ｂ2maf、ｂ2map及びｂ2constは、行列の変換によって算出される（２Ｈp＋Ｈc）行列である。
【０１０１】
このような処理を数式処理を用いてオフラインで行っておけば、エンジン制御装置１０００では、（３７）式を解く演算を行えばよい。行列Ｍ2のサイズは、行列Ｍよりも小さくなっており、計算量も削減されるため、高速に制御入力を得ることができるようになる。
【０１０２】
次に、制御入力に加えて、一般的な状態制約を導入する場合を考える。例えばディーゼルエンジンでは、燃費効率を上げるために、運転時にエンジン本体に送る気体の空気過剰率（一般的にはλを用いて表すが、ここでは記号の重複を避けるためにτを用いる）を一定以上に保つことを求められる。
【０１０３】
空気過剰率τは、各運転条件（すなわち回転数及び燃料噴射料）に応じて、各時刻ｋ＋ｉに対し、ＭＡＦ及びＭＡＰの状態偏差ベクトルｘ（ｋ＋ｉ）と線形化入力ベクトルｕ（ｋ＋ｉ）を用いて以下のように近似される。
【数３５】
2015095126.tif 000037
ここで、行列Ｃ及びＤは、１行２列行列である。
【０１０４】
モデル予測制御を用いた、ディーゼルエンジンのＭＡＦ及びＭＡＰ制御には、以下のような空気過剰率τの制約条件を、各時刻に対して適用する。
【数３６】
2015095126.tif 000038
【０１０５】
よって、モデル予測制御の問題は、以下のように表される。
【数３７】
2015095126.tif 000039
τ0は、空気過剰率の初期値である。
【０１０６】
しかしながら、このままではモデル予測制御問題が不等式制約を含む問題となり、現実的には解けないので、例えば、以下に示すようなペナルティ関数を導入することで、不等式制約問題を等式制約問題に緩和する。
【０１０７】
【数３８】
2015095126.tif 000040
【０１０８】
ここでは、評価関数Ｊbに対して（４１）式のペナルティ関数を加える。そうすると、以下のように表される。
【０１０９】
【数３９】
2015095126.tif 000041
【０１１０】
また、（４０−２）式から各時刻のτ（ｋ＋ｉ）は、具体的には以下のように表される。
【０１１１】
【数４０】
2015095126.tif 000042
【０１１２】
そして、（３０）式を（４２）式に代入すれば、τ（ｋ＋ｉ）は、以下のように書き下される。
【０１１３】
【数４１】
2015095126.tif 000043
【０１１４】
そうすると、τ（ｋ＋ｉ）は、全てのｉについて、ｕ（０）乃至ｕ（Ｈc−１）の関数であるので、（４２）式のペナルティ関数についてもｕ（０）乃至ｕ（Ｈc−１）の関数と考えることができる。よって、入力制約について（２０）式乃至（２２）式で記述された最適化問題は、以下のように変形される。
【０１１５】
【数４２】
2015095126.tif 000044
【０１１６】
このような（４６）式乃至（４８）式を解くことになるが、そのためには状態についてのペナルティ関数に関する∂ｑb（ｘ）／∂ｕを線形化する。
【０１１７】
まず、∂ｑb（τ（ｋ＋ｉ））／∂ｕegr（ｋ＋ｊ）は以下のように表される。また、（４３）式から、２行目も得られる。
【０１１８】
【数４３】
2015095126.tif 000045
【０１１９】
次に、（５１）式中の∂ｑb（τ（ｋ＋ｉ））／∂τ（ｋ＋ｉ）を線形化する。その際、以下の最適な線形化制御入力ｕ~で１次のテイラー展開を行うと、∂ｑb（τ（ｋ＋ｉ））／∂τ（ｋ＋ｉ）は以下のように表される。
【０１２０】
【数４４】
2015095126.tif 000046
【０１２１】
【数４５】
2015095126.tif 000047
【０１２２】
また、（５１）式と同様に、以下のように表される。
【０１２３】
【数４６】
2015095126.tif 000048
【０１２４】
従って、以下のような近似式を得られる。
【０１２５】
【数４７】
2015095126.tif 000049
【０１２６】
よって、（５１）式及び（５４）式を組み合わせると、以下のような近似式が得られる。
【０１２７】
【数４８】
2015095126.tif 000050
【０１２８】
ここで、近似のためペナルティ関数部分の変数として代入するτのベースポイントを設定する。ｘ~（ｋ＋ｉ）（”~”はｘの上に載っているものとする）を、ｕ~に対応する状態偏差ベースポイントと設定する。そして、短い時間における状態偏差の値の変動は小さいので、ｘbase（ｋ）は以下のように表される。
【０１２９】
【数４９】
2015095126.tif 000051
【０１３０】
従って、（４０−２）式から、τのベースポイントτbase（ｋ）は、以下のように表される。なお、ｕ*（ｋ−１）は、１単位時間前の最適な線形化制御入力である。
【０１３１】
【数５０】
2015095126.tif 000052
【０１３２】
このようなベースポイントτbase（ｋ）を用いると、（４３）式は以下のように変形される。
【０１３３】
【数５１】
2015095126.tif 000053
【０１３４】
ここで（５７）式を、以下のように表す。
【０１３５】
【数５２】
2015095126.tif 000054
【０１３６】
さらに、τ（ｋ＋ｉ）をτbase（ｋ）と近似すれば、｛∂2ｑb（τ（ｋ＋ｉ））／∂τ（ｋ＋ｉ）2｝（ｕ~）は、全てのｉについてτbase（ｋ）に関して同じ関数になる。同様に、｛∂ｑb（τ（ｋ＋ｉ））／∂τ（ｋ＋ｉ）｝（ｕ~）は、全てのｉについてτbase（ｋ）に関して同じ関数になる。
【０１３７】
以下、ペナルティ関数のτについての一階偏微分を∂（τ）と表し、二階偏微分を∂2（τ）と表し、τ＝τbase（ｋ）についての∂（τbase（ｋ））及び∂2（τbase（ｋ））によって、（５５）式は以下のように表される。
【０１３８】
【数５３】
2015095126.tif 000055
【０１３９】
この（６２）式から、以下のような近似も行われる。
【０１４０】
【数５４】
2015095126.tif 000056
【０１４１】
さらに、この（６３）式から、以下のような近似も行われる。
【０１４２】
【数５５】
2015095126.tif 000057
【０１４３】
ここで、Λl及びνlは以下のように表される。
【０１４４】
【数５６】
2015095126.tif 000058
【０１４５】
【数５７】
2015095126.tif 000059
【０１４６】
なお、Λlにおけるｐは２ｊ＋１で固定で、ｑについては、１からＨcまで増加する。すなわち、Λlは、１行Ｈc列の行列である。
【０１４７】
従って、（６４）式及び（３７）式から、状態についてのペナルティ関数を伴う最適化問題は、以下のように緩和される。
【０１４８】
【数５８】
2015095126.tif 000060
【０１４９】
（６５）式をエンジン制御装置１０００で解くためには、Ｍ2、ｂ2、Λｌ及びνlを予めオフラインで用意しておき、ペナルティ関数のτについての一階偏微分∂（τ）と、二階偏微分∂2（τ）を予めオフラインで算出しておく。そして、エンジン制御装置１０００では、τbase（ｋ）を（５６）式で算出するとともに、∂（τbase（ｋ））及び∂2（τbase（ｋ））を算出し、さらにｄl（τbase（ｋ））を（５８）式で算出する。そして（６５）式を解くことになる。
【０１５０】
次に、本実施の形態に係るオフライン処理装置１００の機能ブロック図を図３に示す。本実施の形態では、オフライン処理装置１００は、入力部１０１と、入力データ格納部１０２と、第１変形処理部１０３と、第１データ格納部１０４と、第２変形処理部１０５と、第２データ格納部１０６と、第３変形処理部１０７と、第３データ格納部１０８と、第１行列生成部１０９と、第４データ格納部１１０と、第２行列生成部１１１と、微分処理部１１２とを有する。
【０１５１】
次に、オフライン処理装置１００の処理内容を図４を用いて説明する。まず、入力部１０１は、ユーザからの入力データとして、行列Ａ、Ｂ、Ｃ、Ｄ、Ｓ及びＲ、予測ホライズンＨp及び制御ホライズンＨcなどの値の入力を受け付け、入力データ格納部１０２に格納する（図４：ステップＳ１）。なお、入力データ格納部１０２に格納されていない場合には、（４６）乃至（４８）式のデータも入力部１０１を介して入力データ格納部１０２に格納されるものとする。同様に、入力制約に係る項の式のデータ（（２４）式）についても同様に、入力データ格納部１０２に格納されており、格納されていない場合には入力部１０１を介して入力される。さらに、状態偏差ベクトルｘ（ｋ＋ｉ）と線形化入力ベクトルｕ（ｋ＋ｉ）を用いて表されるτ（ｋ＋ｉ）に係る式のデータについても同様に、入力データ格納部１０２に格納されており、格納されていない場合には入力部１０１を介して入力される。
【０１５２】
次に、第１変形処理部１０３は、入力データ格納部１０２に格納されているデータを用いて、状態偏差ｘについての第１変形処理を実行し、処理結果を第１データ格納部１０４に格納する（ステップＳ３）。（２２）式に由来する（２８−１）及び（２８−２）式を（２９）式に代入する処理をｉを増加させながら繰り返して実行し、状態偏差ベクトルについての（３０）式を生成する処理を実行する。
【０１５３】
本実施の形態では、さらに、第１変形処理部１０３は、入力データ格納部１０２に格納されているデータを用いて、制約されるパラメータである空気過剰率τについての変形処理を実行し、処理結果を第１データ格納部１０４に格納する（ステップＳ４）。（３０）式を（４２）式に代入することで、（４３）式に示すようにτ（ｋ＋ｉ）は書き下される。
【０１５４】
また、第２変形処理部１０５は、入力データ格納部１０２に格納されているデータを用いて、ラグランジュ乗数λについての第２変形処理を実行し、処理結果を第２データ格納部１０６に格納する（ステップＳ５）。（２０）式に由来する（３１）及び（３２）式を（３３）式に代入する処理をｉを減少させながら繰り返して実行し、（３４）式を生成する処理を実行する。
【０１５５】
さらに、第３変形処理部１０７は、第１データ格納部１０４及び第２データ格納部１０６に格納されているデータを用いて、ラグランジュ乗数λについての第３変形処理を実行し、処理結果を第３データ格納部１０８に格納する（ステップＳ７）。（３０）式を（３４）式に代入する処理を実行して、（３５）式を生成する。
【０１５６】
そして、第１行列生成部１０９は、（３０）式（第１データ格納部１０４に格納されている）及び（３５）式（第３データ格納部１０８に格納されている）を（２１）式（入力データ格納部１０２に格納されている）に代入して連立方程式（３６）式を生成し、これを（３７）式の形に変形することで、行列Ｍ2及びｂ2を算出し、第４データ格納部１１０に格納する（ステップＳ９）。ステップＳ４以外のここまでの処理は、特開２０１２−１９４９６０号公報に記載の技術とほぼ同様である。
【０１５７】
その後、第２行列生成部１１１は、第１データ格納部１０４に格納されているτについての変形処理の結果（（４３）式）を用いて、入力データ格納部１０２に格納されている式に従って、係数行列Λl及び係数ベクトルνlを生成し、第４データ格納部１１０に格納する（ステップＳ１１）。
【０１５８】
さらに、微分処理部１１２は、第１データ格納部１０２に格納されている、状態制約に対するペナルティ関数のτについての１階偏微分関数∂（τ）と、２階偏微分関数∂2（τ）を生成し、第４データ格納部１１０に格納する（ステップＳ１３）。１階偏微分関数及び２階偏微分関数の生成については、周知であるからここでは説明を省略する。
【０１５９】
以上のような処理を実行することで、エンジン制御装置１０００の処理負荷を大幅に削減でき、処理速度を高速化できる。
【０１６０】
次に、図５にエンジン制御装置１０００におけるモデル予測制御部２２０の構成例を示す。モデル予測制御部２２０は、微分値算出部２２０１と、行列算出部２２０２と、出力生成部２２０３と、データ格納部２２０５とを有する。
【０１６１】
データ格納部２２０５には、第４データ格納部１１０に格納されているデータが格納される。
【０１６２】
次に、図６を用いて、エンジン制御装置１０００におけるモデル予測制御部２２０の処理内容について説明する。
【０１６３】
まず、モデル予測制御部２２０は、ｘ（ｋ）（＝Ｘref（ｋ）−ｙ（ｋ））を取得する（図６：ステップＳ２１）。
【０１６４】
そうすると、モデル予測制御部２２０は、関数ベクトルｂ2にｘ（ｋ）を代入して関数ベクトルｂ2を数値ベクトル化する（ステップＳ２３）。
【０１６５】
また、微分値算出部２２０１は、制約されるパラメータのベースポイントの値であるτbase（ｋ）（（５６）式）を算出し、データ格納部２２０５に格納されているペナルティ関数の１階偏微分関数∂（τ）と、２階偏微分関数∂2（τ）にτbase（ｋ）を代入し、∂（τbase（ｋ））及び∂2（τbase（ｋ））を求める（ステップＳ２５）。
【０１６６】
そして、行列算出部２２０２は、Ｍ3及びｂ3を算出する（ステップＳ２７）。
【０１６７】
（６５）式の左辺におけるベクトルＵ2以外の部分である（Ｍ2＋∂2（τbase（ｋ））Λl）がＭ3である。また、（６５）式の右辺がｂ3である。すなわち、ｂ3＝ｂ2−（∂（τbase（ｋ））−∂2（τbase（ｋ））ｄl（τbase（ｋ））νlである。ｂ3を算出するためには、ｄl（τbase（ｋ））を、（５８）式で算出する。
【０１６８】
そして、出力生成部２２０３は、Ｕ＝Ｍ-13ｂ3によってベクトルＵを算出する（ステップＳ２９）。さらに、出力生成部２２０３は、ベクトルＵから線形化制御入力ｕ（ｋ）（ｉ＝０）を抽出し、出力する（ステップＳ３１）。
【０１６９】
その後、モデル予測制御部２２０は、処理が終了するのかを判断し（ステップＳ３３）、処理終了でなければ処理はステップＳ２１に戻り、処理終了であれば処理は終了する。
【０１７０】
以上のような処理で分かるように、オフライン処理によって算出されているデータを活用でき、Ｍ3のサイズも削減されているので、状態制約についてのペナルティ関数を導入することでステップＳ２５及びＳ２７が追加されていても、高速に線形化制御入力を得ることができるようになる。さらに、ペナルティ関数の近似も適切に行われており、高い精度で適切な制御がなされるようになる。
【０１７１】
なお、エンジンの制御だけではなく、例えばリチウムイオン電池の充放電制御に適用しても良い。リチウムイオン電池の等価回路を図７に示す。リチウムイオン電池の端子電圧Ｖ（ｋ）は、充電率ＳｏＣ（ｋ）と係数ａの積と、抵抗Ｒ１及び容量Ｃ１が並列に接続された部分の電圧ｖ（ｋ）と、直列等価抵抗Ｒ0と流れる電流ｉ（ｋ）との積との和で表される。すなわち、以下の式で表される。
【０１７２】
【数５９】
2015095126.tif 000061
【０１７３】
また、充電率ＳｏＣ（ｋ）は、以下のような関係式に従って変化する。なお、ＦＣは蓄電池の最大蓄積容量を表す。
【０１７４】
【数６０】
2015095126.tif 000062
【０１７５】
さらに、電圧ｖ（ｋ）は、以下のような関係式に従って変化する。
【０１７６】
【数６１】
2015095126.tif 000063
【０１７７】
一方、状態偏差ではなく、状態についての一般的なモデル予測制御は、以下のように表される。
【０１７８】
【数６２】
2015095126.tif 000064
【０１７９】
そうすると、リチウムイオン電池の場合、電圧ｖ（ｋ）と充電率ＳｏＣ（ｋ）を状態ｘ、端子電圧Ｖ（ｋ）を出力ｙ、電流ｉ（ｋ）を入力ｕ（ｋ）として扱えば、以下のように表される。
【０１８０】
【数６３】
2015095126.tif 000065
【０１８１】
ここで、ｙmin＜ｙ（ｋ）＜ｙmaxという制約条件を適用する場合には、以下のようなペナルティ関数を導入する。
【０１８２】
【数６４】
2015095126.tif 000066
【０１８３】
これは、上で述べた（４１）式と同様の形を有する。
【０１８４】
また、（３０）式相当の式は、以下のように表されるので、（１０１）式をそれらの式に代入する。そうすると、（４３）式相当の式が得られるようになる。
【０１８５】
【数６５】
2015095126.tif 000067
【０１８６】
【数６６】
2015095126.tif 000068
【０１８７】
以上のようにすれば、上で述べたのと同様の式の変形ができる。
【０１８８】
［実施の形態２］
上で述べたようなディーゼルエンジンでは、ＭＡＦ及びＭＡＰを変化させる制御を行うが、耐久性の問題から、ＭＡＦ及びＭＡＰの急激な変化は避けることが好ましい。
【０１８９】
本実施の形態では、モデル予測制御を用いた、ディーゼルエンジンのＭＡＦ及びＭＡＰ制御において、ＭＡＦ及びＭＡＰの急激な変化を抑制するための制約を加えるものとする。但し、説明を簡単にするためにＭＡＦの急激な変化を抑制する制約のみを考えるものとする。
【０１９０】
このような問題は、以下のように表される。
【０１９１】
【数６７】
2015095126.tif 000069
【０１９２】
（２０１−５）式に対してペナルティ関数ｑcは、以下のように表される。
【０１９３】
【数６８】
2015095126.tif 000070
【０１９４】
そうすると、モデル予測制御の目的関数は、以下のように表される。
【０１９５】
【数６９】
2015095126.tif 000071
【０１９６】
上でも述べたように、ＭＡＦのみ取り扱うので、（２０１−２）式から、以下の式が得られる。
【０１９７】
【数７０】
2015095126.tif 000072
【０１９８】
そうすると、（３０）式を（２０３）式に代入すると、（４３）式相当の式が得られる。
【０１９９】
【数７１】
2015095126.tif 000073
【０２００】
そうすると、ｄｘ（ｋ＋ｉ）は、全てのｉについて、ｕ（０）乃至ｕ（Ｈc−１）の関数であるので、（２０２）式のペナルティ関数についてもｕ（０）乃至ｕ（Ｈc−１）の関数と考えることができる。よって、入力制約について（２０）式乃至（２２）式で記述された最適化問題は、以下のように変形される。
【０２０１】
【数７２】
2015095126.tif 000074
【０２０２】
このような（２０４）式乃至（２０６）式を解くことになるが、そのためには状態についてのペナルティ関数に関する∂ｑc（ｘ）／∂ｕを線形化する。
【０２０３】
まず、∂ｑc（ｄｘmaf（ｋ＋ｉ））／∂ｕegr（ｋ＋ｊ）は以下のように表される。
【０２０４】
【数７３】
2015095126.tif 000075
【０２０５】
次に、（２０８）式中の∂ｑc（ｄｘmaf（ｋ＋ｉ））／∂ｄｘmaf（ｋ＋ｉ）を線形化する。その際、以下の最適な線形化制御入力ｕ~で１次のテイラー展開を行うと、∂ｑc（ｄｘmaf（ｋ＋ｉ））／∂ｄｘmaf（ｋ＋ｉ）は以下のように表される。
【０２０６】
【数７４】
2015095126.tif 000076
【０２０７】
【数７５】
2015095126.tif 000077
【０２０８】
また、（２０８）式と同様に、以下のように表される。
【０２０９】
【数７６】
2015095126.tif 000078
【０２１０】
従って、以下のような近似式を得られる。
【０２１１】
【数７７】
2015095126.tif 000079
【０２１２】
よって、（２０８）式及び（２１１）式を組み合わせると、以下のような近似式が得られる。
【０２１３】
【数７８】
2015095126.tif 000080
【０２１４】
ここで、近似のためペナルティ関数部分の変数として代入するｄｘmafのベースポイントを設定する。この場合近似は行わずとも、ベースポイントｄｘbasemaf（ｋ）は以下のように表される。
【０２１５】
【数７９】
2015095126.tif 000081
【０２１６】
ディーゼルエンジンの特性上、状態の変化は入力の影響が支配的になるため、初期値の影響が状態の変化に与える影響は極めて小さいものとして取り扱って近似できる。従って、（２０４）式の定数部分は、（２１１）式から以下のように近似される。
【０２１７】
【数８０】
2015095126.tif 000082
【０２１８】
さらに、ｄｘ~maf（ｋ＋ｉ）をｄｘbasemaf（ｋ）と近似すれば、｛∂2ｑc（ｄｘmaf（ｋ＋ｉ））／∂ｄｘmaf（ｋ＋ｉ）2｝（ｕ~）は、全てのｉについてベースポイントｄｘbasemaf（ｋ）に関して同じ関数になる。同様に、｛∂ｑc（ｄｘmaf（ｋ＋ｉ））／∂ｄｘmaf（ｋ＋ｉ）｝（ｕ~）は、全てのｉについてベースポイントｄｘbasemaf（ｋ）に関して同じ関数になる。
【０２１９】
以下、ペナルティ関数のｄｘ（ここではｄｘmaf）についての一階偏微分を∂（ｄｘ）と表し、二階偏微分を∂2（ｄｘ）と表し、ｄｘ＝ｄｘbasemaf（ｋ）についての∂（ｄｘbasemaf（ｋ））及び∂2（ｄｘbasemaf（ｋ））によって、（２１２）式は以下のように表される。
【０２２０】
【数８１】
2015095126.tif 000083
【０２２１】
この（２１５）式から、以下のような近似も行われる。
【０２２２】
【数８２】
2015095126.tif 000084
【０２２３】
さらに、この（２１６）式から、以下のような近似も行われる。
【０２２４】
【数８３】
2015095126.tif 000085
【０２２５】
ここで、Λd及びνdは以下のように表される。
【０２２６】
【数８４】
2015095126.tif 000086
【０２２７】
【数８５】
2015095126.tif 000087
【０２２８】
なお、Λdにおけるｐは２ｊ＋１で固定で、ｑについては、１からＨcまで増加する。すなわち、Λdは、１行Ｈc列の行列である。
【０２２９】
従って、（２１７）式及び（３７）式から、状態についてのペナルティ関数を伴う最適化問題は、以下のように緩和される。
【０２３０】
【数８６】
2015095126.tif 000088
【０２３１】
（２１８）式をエンジン制御装置１０００で解くためには、Ｍ2、ｂ2、Λd及びνdを予めオフラインで用意しておき、ペナルティ関数のｄｘmafについての一階偏微分∂（ｄｘ）と、二階偏微分∂2（ｄｘ）を予めオフラインで算出しておく。そして、エンジン制御装置１０００では、ｄｘbasemaf（ｋ）を（２１３）式で算出するとともに、∂（ｄｘbasemaf（ｋ））及び∂2（ｄｘbasemaf（ｋ））を算出する。そして（２１８）式を解くことになる。
【０２３２】
なお、上ではｄｘmafについて状態制約を設ける例を示したが、ｄｘmapについて状態制約を設けるようにしてもよいし、それら両方に状態制約を設けるようにしてもよい。基本的な処理は上で述べたとおりである。
【０２３３】
なお、本実施の形態に係るオフライン処理装置１００は、図３に示したものと同様の構成を有する。
【０２３４】
次に、オフライン処理装置１００の処理内容も、おおよそ第１の実施の形態と同様であるから、図４を用いて説明する。まず、入力部１０１は、ユーザからの入力データとして、行列Ａ、Ｂ、Ｓ及びＲ、予測ホライズンＨp及び制御ホライズンＨcなどの値の入力を受け付け、入力データ格納部１０２に格納する（図４：ステップＳ１）。なお、入力データ格納部１０２に格納されていない場合には、（２０５）乃至（２０７）式のデータも入力部１０１を介して入力データ格納部１０２に格納されるものとする。同様に、入力制約に係る項の式のデータ（（２４）式）についても同様に、入力データ格納部１０２に格納されており、格納されていない場合には入力部１０１を介して入力される。状態偏差の変化量ｄｘ（ｋ＋ｉ）（＝ｘ（ｋ＋ｉ）−ｘ（ｋ＋ｉ−１））に係る式のデータについても同様に、入力データ格納部１０２に格納されており、格納されていない場合には入力部１０１を介して入力される。
【０２３５】
次に、第１変形処理部１０３は、入力データ格納部１０２に格納されているデータを用いて、状態偏差ｘについての第１変形処理を実行し、処理結果を第１データ格納部１０４に格納する（ステップＳ３）。（２２）式に由来する（２８−１）及び（２８−２）式を（２９）式に代入する処理をｉを増加させながら繰り返して実行し、状態偏差ベクトルについての（３０）式を生成する処理を実行する。
【０２３６】
本実施の形態では、さらに、第１変形処理部１０３は、入力データ格納部１０２に格納されているデータを用いて、制約されるパラメータである状態偏差の変化量ｄｘについての変形処理を実行し、処理結果を第１データ格納部１０４に格納する（ステップＳ４）。制約条件が異なるので、このステップは第１の実施の形態とは異なる。（３０）式を（２０３）式に代入することで、（２０４）式に示すようにｄｘ（ｋ＋ｉ）は書き下される。
【０２３７】
また、第２変形処理部１０５は、入力データ格納部１０２に格納されているデータを用いて、ラグランジュ乗数λについての第２変形処理を実行し、処理結果を第２データ格納部１０６に格納する（ステップＳ５）。（２０）式に由来する（３１）及び（３２）式を（３３）式に代入する処理をｉを減少させながら繰り返して実行し、（３４）式を生成する処理を実行する。
【０２３８】
さらに、第３変形処理部１０７は、第１データ格納部１０４及び第２データ格納部１０６に格納されているデータを用いて、ラグランジュ乗数λについての第３変形処理を実行し、処理結果を第３データ格納部１０８に格納する（ステップＳ７）。（３０）式を（３４）式に代入する処理を実行して、（３５）式を生成する。
【０２３９】
そして、第１行列生成部１０９は、（３０）式（第１データ格納部１０４に格納されている）及び（３５）式（第３データ格納部１０８に格納されている）を（２１）式（入力データ格納部１０２に格納されている）に代入して連立方程式（３６）式を生成し、これを（３７）式の形に変形することで、行列Ｍ2及びｂ2を算出し、第４データ格納部１１０に格納する（ステップＳ９）。
【０２４０】
その後、第２行列生成部１１１は、第１データ格納部１０４に格納されているｄｘについての変形処理の結果（（２０４）式）を用いて、入力データ格納部１０２に格納されている式に従って、係数行列Λd及び係数ベクトルνdを生成し、第４データ格納部１１０に格納する（ステップＳ１１）。このステップについても第１の実施の形態とは異なる。
【０２４１】
さらに、微分処理部１１２は、第１データ格納部１０２に格納されている、状態偏差の変化量についてのペナルティ関数におけるτについての１階偏微分関数∂（ｄｘ）と、２階偏微分関数∂2（ｄｘ）を生成し、第４データ格納部１１０に格納する（ステップＳ１３）。１階偏微分関数及び２階偏微分関数の生成については、周知であるからここでは説明を省略する。
【０２４２】
以上のような処理を実行することで、エンジン制御装置１０００の処理負荷を大幅に削減でき、処理速度を高速化できる。
【０２４３】
また、エンジン制御装置１０００におけるモデル予測制御部２２０の構成例は、第１の実施の形態と同様である。
【０２４４】
データ格納部２２０５には、第４データ格納部１１０に格納されているデータが格納される。
【０２４５】
次に、エンジン制御装置１０００におけるモデル予測制御部２２０の処理内容についても、第１の実施の形態と同様であるから、図６を用いて説明する。
【０２４６】
まず、モデル予測制御部２２０は、ｘ（ｋ）（＝Ｘref（ｋ）−ｙ（ｋ））を取得する（図６：ステップＳ２１）。
【０２４７】
そうすると、モデル予測制御部２２０は、関数ベクトルｂ2にｘ（ｋ）を代入して関数ベクトルｂ2を数値ベクトル化する（ステップＳ２３）。
【０２４８】
また、微分値算出部２２０１は、制約されるパラメータのベースポイントの値であるｄｘbasemaf（ｋ）（（２１３）式）を算出し、データ格納部２２０５に格納されているペナルティ関数の１階偏微分関数∂（ｄｘ）と、２階偏微分関数∂2（ｄｘ）にｄｘbasemaf（ｋ）を代入し、∂（ｄｘbasemaf（ｋ））及び∂2（ｄｘbasemaf（ｋ））を求める（ステップＳ２５）。ここは第１の実施の形態とは異なる。
【０２４９】
そして、行列算出部２２０２は、Ｍ3及びｂ3を算出する（ステップＳ２７）。
【０２５０】
（２１８）式の左辺におけるベクトルＵ2以外の部分である（Ｍ2＋∂2（ｄｘbasemaf（ｋ））Λd）がＭ3である。また、（２１８）式の右辺がｂ3である。すなわち、ｂ3＝ｂ2−（∂（ｄｘbasemaf（ｋ））−∂2（ｄｘbasemaf（ｋ））ｄｘbasemaf（ｋ））νdである。ｂ3を算出するためには、ｄｘbasemaf（ｋ）を、（２１３）式で算出する。ここは第１の実施の形態とは異なる。
【０２５１】
そして、出力生成部２２０３は、Ｕ＝Ｍ-13ｂ3によってベクトルＵを算出する（ステップＳ２９）。さらに、出力生成部２２０３は、ベクトルＵから線形化制御入力ｕ（ｋ）（ｉ＝０）を抽出し、出力する（ステップＳ３１）。
【０２５２】
その後、モデル予測制御部２２０は、処理が終了するのかを判断し（ステップＳ３３）、処理終了でなければ処理はステップＳ２１に戻り、処理終了であれば処理は終了する。
【０２５３】
［実施の形態３］
実際に制御を行う場合、空気過剰率τが制限値を超えてしまう場合や、多少制限値を超えることを許容する場合もある。この場合でも、できるだけ過剰に制限値を超えることは避けることが好ましい。
【０２５４】
この場合には、例えば、以下のようなペナルティ関数を採用する。
【０２５５】
【数８７】
2015095126.tif 000089
【０２５６】
このペナルティ関数は、図８に示すように変化する関数である。このように、τがτmin以下となるｑbが指数関数的に増加する。
【０２５７】
さらに、設計者によっては空気過剰率τが、所望する目標値周辺で推移するような制御を行う場合もある。
【０２５８】
この場合も、過度に目標値から離れることを抑制するために、目標値想定推移区間をτminからτmaxとした場合、以下のようなペナルティ関数を採用する場合もある。
【０２５９】
【数８８】
2015095126.tif 000090
【０２６０】
このペナルティ関数は、図９に示すように変化する関数である。すなわち、τmin乃至τmaxの間ｑbは「０」となるが、その区間を離れるほどｑbは指数関数的に増加する。
【０２６１】
さらに、空気過剰率τをより目標値τrefに近づけたい場合には、以下のようなペナルティ関数を採用する場合もある。
【０２６２】
【数８９】
2015095126.tif 000091
【０２６３】
【数９０】
2015095126.tif 000092
【０２６４】
このようなペナルティ関数を導入したとしても、基本的に第１の実施の形態で説明したオフライン処理は同じである。但し、空気過剰率τの値域によって関数の形が異なるので、それぞれについてτ（ｋ＋ｉ）についての１階偏微分関数及び２階偏微分関数を算出する。すなわち、ステップＳ１３では、各適用範囲について、１階偏微分関数∂（τ）及び２階偏微分関数∂2（τ）を生成する。
【０２６５】
そして、τ（ｋ＋ｉ）の値に応じて、ペナルティ関数の適用範囲を特定し、適用範囲に応じた１階偏微分関数∂（τ）及び２階偏微分関数∂2（τ）又は「制約なし」（ペナルティ関数が「０」の場合）を用いることになる。
【０２６６】
具体的に、図１０を用いてモデル予測制御部２２０の処理内容を説明する。
【０２６７】
まず、モデル予測制御部２２０は、ｘ（ｋ）（＝Ｘref（ｋ）−ｙ（ｋ））を取得する（図１０：ステップＳ４１）。
【０２６８】
そうすると、モデル予測制御部２２０は、関数ベクトルｂ2にｘ（ｋ）を代入して関数ベクトルｂ2を数値ベクトル化する（ステップＳ４３）。
【０２６９】
また、微分値算出部２２０１は、制約されるパラメータの現在値であるτ（ｋ）（＝Ｃｘ（ｋ−１）＋Ｄｕ（ｋ−１））を算出する（ステップＳ４４）。さらに、微分値算出部２２０１は、τ（ｋ）から、ペナルティ関数の適用範囲を特定する（ステップＳ４５）。
【０２７０】
そうすると、微分値算出部２２０１は、特定された適用範囲が、ペナルティ関数が０の範囲であるか否かを判断する（ステップＳ４７）。
【０２７１】
特定された適用範囲が、ペナルティ関数が０の範囲ではない場合には、微分値算出部２２０１は、制約されるパラメータのベースポイントにおける値であるτbase（ｋ）（（５６）式）を算出し、データ格納部２２０５に格納されているペナルティ関数の１階偏微分関数∂（τ）と、２階偏微分関数∂2（τ）にτbase（ｋ）を代入し、∂（τbase（ｋ））及び∂2（τbase（ｋ））を求める（ステップＳ５１）。
【０２７２】
そして、行列算出部２２０２は、特定された適用範囲についてのＭ3及びｂ3を算出する（ステップＳ５３）。
【０２７３】
（６５）式の左辺におけるベクトルＵ2以外の部分である（Ｍ2＋∂2（τbase（ｋ））Λl）がＭ3である。また、（６５）式の右辺がｂ3である。すなわち、ｂ3＝ｂ2−（∂（τbase（ｋ））−∂2（τbase（ｋ））ｄl（τbase（ｋ）））νlである。ｂ3を算出するためには、ｄl（τbase（ｋ））を、（５８）式で算出する。そして処理はステップＳ５７に移行する。
【０２７４】
一方、特定された適用範囲が、ペナルティ関数が０の範囲であれば、微分値算出部２２０１は、Ｍ3＝Ｍ2及びｂ3＝ｂ2と設定する（ステップＳ５５）。状態制約についてのペナルティ関数は適用されない。その後処理はステップＳ５７に移行する。
【０２７５】
そして、出力生成部２２０３は、Ｕ＝Ｍ-13ｂ3によってベクトルＵを算出する（ステップＳ５７）。さらに、出力生成部２２０３は、ベクトルＵから線形化制御入力ｕ（ｋ）（ｉ＝０）を抽出し、出力する（ステップＳ５９）。
【０２７６】
その後、モデル予測制御部２２０は、処理が終了するのかを判断し（ステップＳ６１）、処理終了でなければ処理はステップＳ４１に戻り、処理終了であれば処理は終了する。
【０２７７】
このような処理を実行すれば、適切なペナルティ関数が適用されるようになり、オーバーシュートやアンダーシュートを効果的に抑制できるようになる。
【０２７８】
［実施の形態４］  
状態偏差の変化量に関する制約についても、過度の変化を抑える目的もあるので、制御の方法によっては、制約を超えた変化も許す場合もある。この場合も、過度に制約を超えることを抑制するため、以下のようなペナルティ関数を採用する場合もある。
【０２７９】
【数９１】
2015095126.tif 000093
【０２８０】
このペナルティ関数は、図１１に示すが、図９に示したものと同様である。
【０２８１】
従って、第３の実施の形態で説明したように、オフラインでは、ｄｘの適用範囲毎に、ペナルティ関数の１階偏微分関数∂（ｄｘ）と、２階偏微分関数∂2（ｄｘ）とを算出しておく。さらに、オンライン処理では、ｄｘ（ｋ）（＝ｘ（ｋ）−ｘ（ｋ−１））を算出し、ｄｘ（ｋ）が含まれる適用範囲についての、１階偏微分関数∂（ｄｘ）及び２階偏微分関数∂2（ｄｘ）を用いて、演算を行えばよい。
【０２８２】
以上本技術の実施の形態を説明したが、本技術はこれに限定されるものではない。例えば、上で述べた機能ブロック構成については、プログラムモジュール構成とは一致しない場合もある。また、処理結果が変わらない限り、処理フローについても、ステップの順番を入れ替えたり、複数のステップを並列実行するようにしてもよい。
【０２８３】
なお、上で述べたオフライン処理装置１００は、コンピュータ装置であって、図１２に示すように、メモリ２５０１とＣＰＵ２５０３とハードディスク・ドライブ（ＨＤＤ）２５０５と表示装置２５０９に接続される表示制御部２５０７とリムーバブル・ディスク２５１１用のドライブ装置２５１３と入力装置２５１５とネットワークに接続するための通信制御部２５１７とがバス２５１９で接続されている。オペレーティング・システム（ＯＳ：Operating System）及び本実施例における処理を実施するためのアプリケーション・プログラムは、ＨＤＤ２５０５に格納されており、ＣＰＵ２５０３により実行される際にはＨＤＤ２５０５からメモリ２５０１に読み出される。ＣＰＵ２５０３は、アプリケーション・プログラムの処理内容に応じて表示制御部２５０７、通信制御部２５１７、ドライブ装置２５１３を制御して、所定の動作を行わせる。また、処理途中のデータについては、主としてメモリ２５０１に格納されるが、ＨＤＤ２５０５に格納されるようにしてもよい。本技術の実施例では、上で述べた処理を実施するためのアプリケーション・プログラムはコンピュータ読み取り可能なリムーバブル・ディスク２５１１に格納されて頒布され、ドライブ装置２５１３からＨＤＤ２５０５にインストールされる。インターネットなどのネットワーク及び通信制御部２５１７を経由して、ＨＤＤ２５０５にインストールされる場合もある。このようなコンピュータ装置は、上で述べたＣＰＵ２５０３、メモリ２５０１などのハードウエアとＯＳ及びアプリケーション・プログラムなどのプログラムとが有機的に協働することにより、上で述べたような各種機能を実現する。
【０２８４】
また、上で述べたエンジン制御装置１０００は、コンピュータ装置であって、図１３に示すように、ＲＡＭ（Random Access Memory）４５０１とプロセッサ４５０３とＲＯＭ（Read Only Memory）４５０７とセンサ群４５１５とがバス４５１９で接続されている。本実施の形態における処理を実施するための制御プログラム及び存在している場合にはオペレーティング・システム（ＯＳ：Operating System））は、ＲＯＭ４５０７に格納されており、プロセッサ４５０３により実行される際にはＲＯＭ４５０７からＲＡＭ４５０１に読み出される。プロセッサ４５０３は、センサ群４５１５（例えば、吸気圧センサ及び新気量センサ。場合によっては燃料噴射量測定部及びエンジン回転数測定部など。）を制御して、測定値を取得する。また、処理途中のデータについては、ＲＡＭ４５０１に格納される。なお、プロセッサ４５０３は、ＲＯＭ４５０７を含む場合もあり、さらに、ＲＡＭ４５０１を含む場合もある。本技術の実施の形態では、上で述べた処理を実施するための制御プログラムは、コンピュータ読み取り可能なリムーバブル・ディスクに格納されて頒布され、ＲＯＭライタによってＲＯＭ４５０７に書き込まれる場合もある。このようなコンピュータ装置は、上で述べたプロセッサ４５０３、ＲＡＭ４５０１、ＲＯＭ４５０７などのハードウエアと制御プログラム（場合によってはＯＳも）とが有機的に協働することにより、上で述べたような各種機能を実現する。
【０２８５】
以上述べた本実施の形態をまとめると以下のようになる。
【０２８６】
本実施の形態の第１の態様に係るオフライン処理方法は、（１）（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ制御対象の状態方程式に関する第１の関数との和を含む第２の関数をラグランジュ乗数について偏微分することで得られ且つ制御対象の状態偏差及び制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）第２の関数を制御対象の状態偏差について偏微分することで得られ且つ制御対象の状態偏差及びラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）第２の関数を制御対象の線形化制御入力について偏微分することで得られ且つラグランジュ乗数及び制御対象の状態偏差を変数として含む第３の一次式群と、（ｄ）モデル予測制御の制御対象についての制約条件に対応するペナルティ関数と、（ｅ）制約条件において制限されるパラメータに関する関係式とを格納するデータ格納部から第１の一次式群を読み出して、状態偏差の現在値を除き状態偏差の変数を消し込むように順次代入を行って、第１の一次式群を、状態偏差の現在値と線形化制御入力とを変数として含む第４の一次式群に変形し、（２）データ格納部から第２の一次式群を読み出して、ラグランジュ乗数の変数を消し込むように順次代入を行って、第２の一次式群を、状態偏差を変数として含む第５の一次式群に変形し、（３）データ格納部から第３の一次式群を読み出して、第４の一次式群及び第５の一次式群を代入することで、第３の一次式群を、状態偏差の現在値と線形化制御入力とを変数として含む第６の一次式群に変形し、（４）第６の一次式群を、第１の係数行列と線形化制御入力についてのベクトルとの積が、状態偏差の現在値の関数ベクトルと等式で表されるように変形し、（５）データ格納部に格納されたペナルティ関数をパラメータについて１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、（６）データ格納部から関係式を読み出して、第４の一次式群を関係式に代入することで、関係式を、状態偏差の現在値と線形化制御入力とを変数として含む第７の一次式群に変形し、（７）第７の一次式群から線形化制御入力についての係数を抽出し、ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する処理を含む。
【０２８７】
このような処理をオフラインで実行しておくことにより、様々な状態制約を導入したモデル予測制御を行うことができるようになる。
【０２８８】
なお、上で述べたパラメータが、状態偏差と線形化制御入力とから表されるパラメータである場合もある。このような場合には、ディーゼルエンジンの空気過剰率などに対して制約条件を設定することができるようになる。
【０２８９】
さらに、上で述べたパラメータが、状態偏差の変化量である場合もある。このような場合には、状態偏差についての変化が急激に発生しないようにすることができるようになる。
【０２９０】
さらに、上で述べたペナルティ関数が、パラメータの値域毎に異なる関数で表される場合がある。この場合、上で述べた第１及び第２の微分関数を、パラメータの値域毎に生成するようにしても良い。このようにすれば、制御実行時に第１の関数と第２の関数とを切り替えて用いることができるようになる。
【０２９１】
本実施の形態の第２の態様に係る制御方法は、（１）（ａ）モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、（ｂ）制御期間内における制御対象に対する線形化制御入力のベクトルとの積が関数ベクトルと等しくなるように予め算出されている第１の係数行列と、（ｃ）モデル予測制御の制御対象における制約条件に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、（ｄ）制約条件において制限されるパラメータについてペナルティ関数を１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から第１の式及び第２の式を読み出して、状態偏差の現在値から得られる、パラメータの基点の値を代入することで、第１の値及び第２の値を算出し、（２）データ格納部から第１の係数行列及び第２の係数行列を読み出して、第２の値と第２の係数行列との積と第１の係数行列との和である第１の行列を算出し、（３）データ格納部から関数ベクトルを読み出して、関数ベクトルに状態偏差の現在値を代入することで得られる第１のベクトルから、（（第１の値）−（第２の値とパラメータの基点に関連し且つ状態偏差の現在値又は線形化制御入力の１単位時間前の値から得られる第３の値との積））と係数ベクトルとの積を差し引いた第２のベクトルを算出し、（４）第１の行列の逆行列と第２のベクトルとの積により、線形化制御入力のベクトルを生成し、（５）線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する処理を含む。
【０２９２】
このような処理を実行することで、高速に状態制約を伴うモデル予測制御を実行できるようになる。
【０２９３】
また、第２の態様に係る制御方法は、状態偏差の現在値から算出されるパラメータの現在値が所定範囲内であるか判断する処理をさらに含むようにしても良い。この場合、パラメータの現在値が所定範囲内であれば、第１の値及び第２の値を算出する処理と、第１の行列を算出する処理と、第２のベクトルを算出する処理と、線形化制御入力のベクトルを生成する処理と、線形化制御入力を抽出する処理とを実行するようにしても良い。制御対象によっては、状態制約に基づく処理を行うことが好ましくない場合もあるためである。
【０２９４】
さらに、上で述べたデータ格納部が、パラメータの現在値の値域に応じた第１の式及び第２の式を格納している場合もある。この場合、上記第１の値及び第２の値を算出する処理が、パラメータの現在値が属する値域を特定し、当該特定された値域に応じた第１の式及び第２の式を読み出す処理を含むようにしても良い。適切な式を状況に応じて切り分けて用いることで、適切な制御を行うことができるようになる。
【０２９５】
また、上で述べたパラメータの現在値が所定範囲外であれば、データ格納部から関数ベクトル及び第１の係数行列を読み出して、関数ベクトルに状態偏差の現在値を代入することで得られる第１のベクトルと、第１の係数行列の逆行列との積により、線形化制御入力の第２のベクトルを生成し、線形化制御入力の第２のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する処理をさらに含むようにしても良い。状態制約を考慮しないモデル予測制御を高速に実行できる。
【０２９６】
なお、上で述べたパラメータが、状態偏差と線形化制御入力とから表されるパラメータである場合には、上で述べた第３の値は、制御対象の状態方程式で用いられる行列と線形化制御入力の１単位時間前の値とから算出される場合もある。例えば、空気過剰率の場合には、このような演算が行われる。
【０２９７】
さらに、上で述べたパラメータが、状態偏差の変化量である場合には、第３の値は、状態偏差の現在値と状態偏差の１単位時間前の値との差である場合もある。状態偏差の変化量に制約がある場合には、このような演算が行われる。
【０２９８】
なお、上記方法による処理をプロセッサ（又はコンピュータ）に行わせるためのプログラムを作成することができ、当該プログラムは、例えばフレキシブルディスク、ＣＤ−ＲＯＭ、光磁気ディスク、半導体メモリ、ハードディスク等のコンピュータ読み取り可能な記憶媒体又は記憶装置に格納される。尚、中間的な処理結果はメインメモリ等の記憶装置に一時保管される。
【０２９９】
以上の実施例を含む実施形態に関し、さらに以下の付記を開示する。
【０３００】
（付記１）
（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ前記制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群と、（ｄ）前記モデル予測制御の制御対象についての制約条件に対応するペナルティ関数と、（ｅ）前記制約条件において制限されるパラメータに関する関係式とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形し、
前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む第５の一次式群に変形し、
前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を代入することで、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形し、
前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形し、
前記データ格納部に格納された前記ペナルティ関数を前記パラメータについて１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、
前記データ格納部から前記関係式を読み出して、前記第４の一次式群を前記関係式に代入することで、前記関係式を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第７の一次式群に変形し、
前記第７の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する
処理を、コンピュータに実行させるためのプログラム。
【０３０１】
（付記２）
前記パラメータが、前記状態偏差と前記線形化制御入力とから表されるパラメータである
付記１記載のプログラム。
【０３０２】
（付記３）
前記パラメータが、前記状態偏差の変化量である
付記１記載のプログラム。
【０３０３】
（付記４）
前記ペナルティ関数が、前記パラメータの値域毎に異なる関数で表され、
前記第１及び第２の微分関数を、前記パラメータの値域毎に生成する
付記１乃至３のいずれか１つ記載のプログラム。
【０３０４】
（付記５）
（ａ）モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、（ｂ）制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、（ｃ）前記モデル予測制御の制御対象における制約条件に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、（ｄ）前記制約条件において制限されるパラメータについて前記ペナルティ関数を１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び前記第２の式を読み出して、前記状態偏差の現在値から得られる、前記パラメータの基点の値を代入することで、第１の値及び第２の値を算出し、
前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出し、
前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルから、（（前記第１の値）−（前記第２の値と前記パラメータの基点に関連し且つ前記状態偏差の現在値又は前記線形化制御入力の１単位時間前の値から得られる第３の値との積））と前記係数ベクトルとの積を差し引いた第２のベクトルを算出し、
前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成し、
前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
処理をプロセッサに実行させるためのプログラム。
【０３０５】
（付記６）
前記状態偏差の現在値から算出される前記パラメータの現在値が所定範囲内であるか判断する処理をさらに含み、
前記パラメータの現在値が前記所定範囲内であれば、前記第１の値及び第２の値を算出する処理と、前記第１の行列を算出する処理と、前記第２のベクトルを算出する処理と、前記線形化制御入力のベクトルを生成する処理と、前記線形化制御入力を抽出する処理とを実行する
付記５記載のプログラム。
【０３０６】
（付記７）
前記データ格納部が、
前記パラメータの現在値の値域に応じた前記第１の式及び前記第２の式を格納しており、
前記第１の値及び前記第２の値を算出する処理が、
前記パラメータの現在値が属する値域を特定し、当該特定された値域に応じた前記第１の式及び前記第２の式を読み出す処理
を含む付記５又は６記載のプログラム。
【０３０７】
（付記８）
前記パラメータの現在値が前記所定範囲外であれば、
前記データ格納部から前記関数ベクトル及び前記第１の係数行列を読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルと、前記第１の係数行列の逆行列との積により、前記線形化制御入力の第２のベクトルを生成し、
前記線形化制御入力の第２のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
処理をさらに含む付記５乃至７のいずれか１つ記載のプログラム。
【０３０８】
（付記９）
前記パラメータが、前記状態偏差と前記線形化制御入力とから表されるパラメータである場合には、
前記第３の値は、前記制御対象の状態方程式で用いられる行列と前記線形化制御入力の１単位時間前の値とから算出される
付記５乃至８のいずれか１つ記載のプログラム。
【０３０９】
（付記１０）
前記パラメータが、前記状態偏差の変化量である場合には、
前記第３の値は、前記状態偏差の現在値と前記状態偏差の１単位時間前の値との差である
付記１記載のプログラム。
【０３１０】
（付記１１）
（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ前記制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群と、（ｄ）前記モデル予測制御の制御対象についての制約条件に対応するペナルティ関数と、（ｅ）前記制約条件において制限されるパラメータに関する関係式とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形し、
前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む第５の一次式群に変形し、
前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を代入することで、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形し、
前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形し、
前記データ格納部に格納された前記ペナルティ関数を前記パラメータについて１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、
前記データ格納部から前記関係式を読み出して、前記第４の一次式群を前記関係式に代入することで、前記関係式を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第７の一次式群に変形し、
前記第７の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する
処理を含み、コンピュータにより実行されるオフライン処理方法。
【０３１１】
（付記１２）
（ａ）モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、（ｂ）制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、（ｃ）前記モデル予測制御の制御対象における制約条件に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、（ｄ）前記制約条件において制限されるパラメータについて前記ペナルティ関数を１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び前記第２の式を読み出して、前記状態偏差の現在値から得られる、前記パラメータの基点の値を代入することで、第１の値及び第２の値を算出し、
前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出し、
前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルから、（（前記第１の値）−（前記第２の値と前記パラメータの基点に関連し且つ前記状態偏差の現在値又は前記線形化制御入力の１単位時間前の値から得られる第３の値との積））と前記係数ベクトルとの積を差し引いた第２のベクトルを算出し、
前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成し、
前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
処理を含み、プロセッサにより実行される制御方法。
【０３１２】
（付記１３）
（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ前記制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群と、（ｄ）前記モデル予測制御の制御対象についての制約条件に対応するペナルティ関数と、（ｅ）前記制約条件において制限されるパラメータに関する関係式とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形する手段と、
前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む第５の一次式群に変形する手段と、
前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を代入することで、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形する手段と、
前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形する手段と、
前記データ格納部に格納された前記ペナルティ関数を前記パラメータについて１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成する手段と、
前記データ格納部から前記関係式を読み出して、前記第４の一次式群を前記関係式に代入することで、前記関係式を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第７の一次式群に変形する手段と、
前記第７の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する手段と、
を有するオフライン処理装置。
【０３１３】
（付記１４）
（ａ）モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、（ｂ）制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、（ｃ）前記モデル予測制御の制御対象における制約条件に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、（ｄ）前記制約条件において制限されるパラメータについて前記ペナルティ関数を１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び前記第２の式を読み出して、前記状態偏差の現在値から得られる、前記パラメータの基点の値を代入することで、第１の値及び第２の値を算出する手段と、
前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出する手段と、
前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルから、（（前記第１の値）−（前記第２の値と前記パラメータの基点に関連し且つ前記状態偏差の現在値又は前記線形化制御入力の１単位時間前の値から得られる第３の値との積））と前記係数ベクトルとの積を差し引いた第２のベクトルを算出する手段と、
前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成する手段と、
前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する手段と、
を有する制御装置。
【符号の説明】
【０３１４】
１０１  入力部
１０２  入力データ格納部
１０３  第１変形処理部
１０４  第１データ格納部
１０５  第２変形処理部
１０６  第２データ格納部
１０７  第３変形処理部
１０８  第３データ格納部
１０９  第１行列生成部
１１０  第４データ格納部
１１１  第２行列生成部
１１２  微分処理部
【図１】
2015095126.tif 000094
【図２】
2015095126.tif 000095
【図３】
2015095126.tif 000096
【図４】
2015095126.tif 000097
【図５】
2015095126.tif 000098
【図６】
2015095126.tif 000099
【図７】
2015095126.tif 000100
【図８】
2015095126.tif 000101
【図９】
2015095126.tif 000102
【図１０】
2015095126.tif 000103
【図１１】
2015095126.tif 000104
【図１２】
2015095126.tif 000105
【図１３】
2015095126.tif 000106