(57)【要約】【課題】モデル予測制御において制御対象の状態制約を適切に取り扱う。【解決手段】制御対象の状態制約を導入しない場合に用意するデータに加え、モデル予測制御の制御対象における状態制約に対応するペナルティ関数を状態偏差について１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、上記データに含まれる一次式群から線形化制御入力についての係数を抽出し、ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する。【選択図】図１
2014206802.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とは加算されて操作量が生成され、当該操作量が飽和要素２３０を介してプラント２４０に入力される。このような処理が繰り返される。
【００２０】
このようなエンジン本体１の吸気系システムは、ＥＧＲバルブ及びＶＮＴノズルについて非線形であり、エンジンのモデルとしても、エンジン回転数ω及び燃料噴射量σに応じて非線形に変化する。
【００２１】
ここで、エンジン回転数ω及び燃料噴射量σが固定であるとして、エンジンモデルは、以下のように定式化される。
【数１】
2014206802.tif 000003
【００２２】
ここでＡω(k),σ(k)及びＢω(k),σ(k)は、時刻ｋにおける運転条件（ω（ｋ），σ（ｋ））によって変化する非線形の２×２行列関数である。なお、時刻ｋにおける回転数をω（ｋ）と表し、燃料噴射量をσ（ｋ）と表す。また、時刻ｋにおけるＭＡＦの状態偏差をｘmaf（ｋ）と表し、ＭＡＰの状態偏差をｘmap（ｋ）と表す。また、時刻ｋにおけるＭＡＦの測定値をｙmaf（ｋ）と表し、ＭＡＰの測定値をｙmap（ｋ）と表す。さらに、時刻ｋにおけるＥＧＲバルブの状態をｖegr（ｋ）と表し、ＶＮＴノズルの状態をｖvnt（ｋ）と表す。また、吸気圧ＭＡＰの目標値をＸmaprefと表し、新気量ＭＡＦの目標値をＸmafrefと表す。
【００２３】
従って、状態偏差ベクトルｘ（ｋ）、測定値ベクトルｙ（ｋ）、目標値ベクトルＸref（ｋ）及び状態ベクトルｖ（ｋ）は、以下のように表される。
【数２】
2014206802.tif 000004
【００２４】
ｘ（ｋ）、ｙ（ｋ）及びＸref（ｋ）の関係は以下のようになっている。
【数３】
2014206802.tif 000005
【００２５】
また、非線形関数ｇを以下のように表す。
【数４】
2014206802.tif 000006
【００２６】
ｇ1及びｇ2は、単調増加関数である。実験データからすると、このような仮定は妥当性がある。非線型モデルでは計算速度等から取り扱いが困難なため、以下では、非線形なエンジンモデルを線型モデルに変換する。そこで、上で述べた線形化制御入力ｕegr及びｕvntを導入することで、線形化する。
【数５】
2014206802.tif 000007
【００２７】
なお、線形化制御入力ベクトルｕ（ｋ）を以下のように表す。
【数６】
2014206802.tif 000008
【００２８】
本実施の形態では、典型的な運転条件を選択して、セットポイントと呼ぶことにする。そして、モデル予測制御部２２０は、運転条件（ω（ｋ），σ（ｋ））を取得すると、対応するセットポイントについてのω及びσについての行列Ａω,σ及びＢω,σを選択し、以下で述べる処理を行って線形化制御入力ｕを算出する。また、変換関数２２５は、線形化制御入力ｕを、変換関数ｇ-1により非線形の最適状態ｖに変換する。この最適状態ｖと状態目標値Ｖrefとは加算されて操作量が生成され、エンジンのプラント２４０に出力される。
【００２９】
モデル予測制御部２２０は、セットポイントについて最適な線形化制御入力ｕを、１サンプル時間毎に、以下の線形モデルに従って算出する。
【数７】
2014206802.tif 000009
【００３０】
但し、線形化制御入力ｕについての制約条件は以下のとおりである。
【数８】
2014206802.tif 000010
【００３１】
例えば、ω＝２０００ｒｐｍ、σ＝２０ｍｍ3／ｓｔの場合、以下のような行列が用いられる。
【数９】
2014206802.tif 000011
【数１０】
2014206802.tif 000012
【００３２】
各時刻ｋについての線形化制御入力ｕは、以下に示すような評価関数Ｊを最小化することで決定される。なお、予測ホライズンＨp及び制御ホライズンＨcは、予め設定されるものとする。
【数１１】
2014206802.tif 000013
【００３３】
ここで、ＱTは行列Ｑの転置行列である。また、Ｑ∈Ｍ（ｎ，Ｒ）、Ｒ∈Ｍ（ｎ，Ｒ）及びＳ∈Ｍ（ｎ，Ｒ）は、重み行列であり、予め設定される。なお、Ｍ（ｎ，Ｒ）はｎ×ｎ行列を表す。（９）式において、無限ホライズンについての最適化問題を、有限ホライズンについての最適化問題にしており、このため、Ricatti代数方程式を解く行列として重み行列Ｓを用いている。
【００３４】
本実施の形態では、モデル予測制御の問題は、有限ホライズンの最適制御問題を解くことに変形される。
【数１２】
2014206802.tif 000014
ここで、ｘは、現在の状態偏差ベクトルを表す。
【００３５】
最適解を求める手法として、ラグランジュの未定乗数法が効率的であるが、不等式制約がある場合にはこの手法をそのまま用いることができない。従って、不等式制約の代わりにペナルティ関数を導入して、不等式制約のない評価関数に緩和する。
【００３６】
すなわち、以下のような評価関数Ｊbに緩和する。
【数１３】
2014206802.tif 000015
【００３７】
ここで、ｐ（ｕ（ｋ＋ｉ））はペナルティ関数であり、入力値ｕ（ｋ＋ｉ）が、可能な値域の境界に近づくと、大きな値を評価関数Ｊに加えることになる。
【００３８】
一般的に、不等式制約ｕ≦ｕmaxに対して用いられるペナルティ関数としては以下のようなものがある。
【数１４】
2014206802.tif 000016
ｒは、ペナルティ関数の係数であり、不等式制約に対する重みを表している。
【００３９】
本実施の形態では、以下の演算のためには、（ｃ）のようなシンプルなペナルティ関数が好ましいので、（ｃ）を用いることにする。
【００４０】
ディーゼルエンジンの場合、線形化制御入力ｕには以下のような制約が存在する。
【数１５】
2014206802.tif 000017
【００４１】
（ｃ）のような形のペナルティ関数を採用する場合には、ペナルティ関数は以下のように表される。
【数１６】
2014206802.tif 000018
【００４２】
（１５）式を用いて制約条件を緩和する場合、以下に示すような、有限ホライズンについての最適制御問題を解くことに帰着する。
【数１７】
2014206802.tif 000019
【００４３】
このような最適制御問題を解くためには、ラグランジュの未定乗数法が用いられる。ここで、ラグランジュ乗数λ（ｊ）＝［λmaf（ｊ），λmap（ｊ）］Tを用い、ハミルトニアンＨを以下のように表す。
【数１８】
2014206802.tif 000020
【００４４】
ここで、Х（χの大文字）は、以下のように定義される。
【数１９】
2014206802.tif 000021
【００４５】
∇x,u,λＨ＝０と置けば、以下のような方程式が得られる。
【数２０】
2014206802.tif 000022
なお、（２０）式において∂Ｊb／∂ｘ＝∂Ｊ／∂ｘ＋∂ｐ（ｕ）／∂ｘ＝∂Ｊ／∂ｘと変形されている。
【００４６】
なお、∂／∂ｘ，∂／∂ｕ，∂／∂λは、以下のような演算である。
【数２１】
2014206802.tif 000023
【００４７】
なお、方程式（２０）乃至（２２）は、∂ｐ／∂ｕの部分を除き線形方程式となっている。∂ｐ／∂ｕを、直前のステップの最適線形化制御入力ｕ~（”~”はｕの上に付される）でテイラー展開することで線形化すると、以下のように表される。
【数２２】
2014206802.tif 000024
【００４８】
（２３）式の第１項のＥＧＲ成分は、以下のように表される。なお、ＶＮＴ成分も同様の形で表される。
【数２３】
2014206802.tif 000025
【００４９】
このような線形近似を行えば、（２０）乃至（２２）式は線形となる。そうすると、線形近似された（２０）乃至（２２）式は以下のように表される。
【数２４】
2014206802.tif 000026
ここでＭは、（２Ｈp＋Ｈc）×（２Ｈp＋Ｈc）の行列であり、Ｕ及びｂは、（２Ｈp＋Ｈc）ベクトルである。ベクトルＵは、以下のように表される。
【数２５】
2014206802.tif 000027
【００５０】
また、ベクトルｂは、以下のように表される。なお、ｂmaf、ｂmap及びｂconstは、行列の変換によって算出される（２Ｈp＋Ｈc）行列である。
ｂ＝ｘmafｂmaf＋ｘmapｂmap＋ｂconst    （２７）
【００５１】
最適な線形化制御入力は、この線形方程式を解くことによって得られる。しかしながら、行列Ｍのサイズが大きければ計算時間は長くなってしまうので、行列Ｍのサイズをオフラインの前処理にて削減する。
【００５２】
まず、（２２）式の第１乃至第３の方程式を具体的に書き下すと以下のようになる。
【数２６】
2014206802.tif 000028
【００５３】
（２８−１）式及び（２８−２）式の状態偏差ｘmaf（ｋ）及びｘmap（ｋ）は、現在値である。従って、（２８−１）式によれば、ｘmaf（ｋ＋１）は線形化制御入力ｕ（ｋ）の１次式で表される。（２８−２）式によれば、ｘmap（ｋ＋１）も、線形化制御入力ｕ（ｋ）の１次式で表される。また、（２９）式によれば、ｘmaf（ｋ＋２）は、状態偏差及び線形化制御入力の１次式で表される。しかしながら、（２８−１）式及び（２８−２）式を（２９）式に代入すると、ｘmaf（ｋ＋２）を、線形化制御入力ｕ（ｋ）及びｕ（ｋ＋１）で表すことができるようになる。これをｉを増加させつつ繰り返せば、状態偏差ｘmaf（ｋ＋ｉ）及びｘmap（ｋ＋ｉ）は、線形化制御入力の１次式として以下のように表される。
【数２７】
2014206802.tif 000029
【００５４】
一方、（２０）式の（２Ｈp−１）番目及び（２Ｈp−３）番目の方程式を具体的に書き下すと以下のようになる。
【数２８】
2014206802.tif 000030
ここで、ｓ11、ｓ22、ｒ11は、重み行列Ｓ及びＲの行列成分を表す。
【００５５】
この（３１）式及び（３２）式は、λmaf（Ｈp）及びλmap（Ｈp）の一次式で表されることが分かる。そこで、（３１）式及び（３２）式を、（３３）式に代入すると、λmaf（Ｈp−１）が、目標値Ｘrefmaf、状態偏差ｘ（ｋ＋Ｈp）及びｘ（ｋ＋Ｈp−１）で表すことができる。このような処理をｉを減少させつつ順次行うことで、ラグランジュ乗数λ（ｋ＋ｉ）は、以下のように表される。
【数２９】
2014206802.tif 000031
【００５６】
さらに、（３０）式を（３４）式に代入することで、ラグランジュ乗数λ（ｋ＋ｉ）は、線形化制御入力ｕ（ｋ）の一次式に変形される。
【数３０】
2014206802.tif 000032
【００５７】
最後に、（３０）式及び（３５）式を、（２１）式に代入すると、以下のように表される。
【数３１】
2014206802.tif 000033
【００５８】
このように、（２１）式は、状態偏差の初期値ｘ（ｋ）及び線形化制御入力ｕの一次式で表されることができるようになる。そうすると、（３６）式は、以下の形式に直すことができる。
【数３２】
2014206802.tif 000034
【００５９】
ここでＭ2はＨc×Ｈc行列である。Ｕ2及びｂ2は、以下のような形のベクトルである。
【数３３】
2014206802.tif 000035
【数３４】
2014206802.tif 000036
なお、ｂ2maf、ｂ2map及びｂ2constは、行列の変換によって算出される（２Ｈp＋Ｈc）行列である。
【００６０】
このような処理を数式処理を用いてオフラインで行っておけば、エンジン制御装置１０００では、（３７）式を解く演算を行えばよい。行列Ｍ2のサイズは、行列Ｍよりも小さくなっており、計算量も削減されるため、高速に制御入力を得ることができるようになる。
【００６１】
次に、定常偏差を削減するために、状態制約を導入する場合を考える。説明を簡単にするために、ＭＡＦの状態偏差について以下のような制約を導入する。ＭＡＰについてはここでは取り扱わないが、同様に設定する場合もある。
【数３５】
2014206802.tif 000037
【００６２】
しかしながら、このままではモデル予測制御問題が不等式制約を含む問題となり、現実的には解けないので、例えば、以下に示すようなペナルティ関数を導入することで、不等式制約問題を等式制約問題に緩和する。
【数３６】
2014206802.tif 000038
【００６３】
ここでは、評価関数Ｊbに対して（４２）式のペナルティ関数を加える。そうすると、以下のように表される。
【数３７】
2014206802.tif 000039
【００６４】
そして、（３０）式からすると、ｘ（ｋ＋ｉ）は、以下のように書き下される。
【数３８】
2014206802.tif 000040
【００６５】
そうすると、ｘ（ｋ＋ｉ）は、全てのｉについて、ｕ（０）乃至ｕ（Ｈc−１）の関数であるので、（４２）式のペナルティ関数についてもｕ（０）乃至ｕ（Ｈc−１）の関数と考えることができる。よって、入力制約について（２０）式乃至（２２）式で記述された最適化問題は、以下のように変形される。
【数３９】
2014206802.tif 000041
【００６６】
このような（４６）式乃至（４８）式を解くことになるが、そのためには状態についてのペナルティ関数に関する∂ｑ（ｘ）／∂ｕを線形化する。
【００６７】
まず、∂ｑ（ｘmaf（ｋ＋ｉ））／∂ｕegr（ｋ＋ｊ）は以下のように表される。
【数４０】
2014206802.tif 000042
【００６８】
ここで（４５）式があるので、（４９）式は、以下のように表される。
【数４１】
2014206802.tif 000043
【００６９】
次に、（５１）式中の∂ｑ（ｘmaf（ｋ＋ｉ））／∂ｘmaf（ｋ＋ｉ）を線形化する。その際、以下の最適な線形化制御入力ｕ~で１次のテイラー展開を行うと、∂ｑ（ｘmaf（ｋ＋ｉ））／∂ｘmaf（ｋ＋ｉ）は以下のように表される。
【数４２】
2014206802.tif 000044
【数４３】
2014206802.tif 000045
【００７０】
また、（５１）式と同様に、以下のように表される。
【数４４】
2014206802.tif 000046
【００７１】
従って、以下のような近似式を得られる。
【数４５】
2014206802.tif 000047
【００７２】
よって、（５１）式及び（５４）式を組み合わせると、以下のような近似式が得られる。
【数４６】
2014206802.tif 000048
【００７３】
ここで、ｘ~maf（ｋ＋ｉ）（”~”はｘの上に載っているものとする）を、ｕ~に対応する状態偏差ベースポイントと設定する。そして、短い時間における状態偏差の値の変動は小さいので、ｘbasemaf（ｋ）は以下のように表される。
【数４７】
2014206802.tif 000049
そして、全てのｉについてｘ~maf（ｋ＋ｉ）をｘbasemaf（ｋ）と近似する。ｘ（ｋ）は、上でも述べたように現在の状態偏差の値を表しており、ｕ*（ｋ−１）は、１単位時間前の最適な線形化制御入力である。
【００７４】
そうすると、（４５）式から、（５５）式の定数部分は、以下のように表される。
【数４８】
2014206802.tif 000050
【００７５】
また、ｘ~maf（ｋ）をｘbasemaf（ｋ）と近似する。そうすると、（５６）式から、以下のような近似式が得られる。
【数４９】
2014206802.tif 000051
【００７６】
そこで、
ｄ（ｘbasemaf（ｋ））＝（１  ０）Ｂｕ（ｋ−１）  （１０１）と置く。
【００７７】
さらに、ｘ~maf（ｋ＋ｉ）をｘbasemaf（ｋ）と近似すれば、｛∂2ｑ（ｘmaf（ｋ＋ｉ））／∂ｘmaf（ｋ＋ｉ）2｝（ｕ~）は、全てのｉについてベース状態偏差ｘbasemaf（ｋ）に関して同じ関数になる。同様に、｛∂ｑ（ｘmaf（ｋ＋ｉ））／∂ｘmaf（ｋ＋ｉ）｝（ｕ~）は、全てのｉについてベース状態偏差ｘbasemaf（ｋ）に関して同じ関数になる。
【００７８】
以下、ペナルティ関数のｘ（（４２）式であればｘmaf）についての一階偏微分を∂（ｘ）と表し、二階偏微分を∂2（ｘ）と表し、ｘ＝ｘbasemaf（ｋ）についての∂（ｘbasemaf（ｋ））及び∂2（ｘbasemaf（ｋ））によって、（５５）式は以下のように表される。
【数５０】
2014206802.tif 000052
【００７９】
この（６２）式から、以下のような近似も行われる。
【数５１】
2014206802.tif 000053
【００８０】
さらに、この（６３）式から、以下のような近似も行われる。
【数５２】
2014206802.tif 000054
ここで、Λ及びνは以下のように表される。
【数５３】
2014206802.tif 000055
【数５４】
2014206802.tif 000056
なお、Λにおけるｐは２ｊ＋１で固定で、ｑについては、１からＨcまで増加する。すなわち、Λは、１行Ｈc列の行列である。
【００８１】
従って、（６４）式及び（３７）式から、状態についてのペナルティ関数を伴う最適化問題は、以下のように緩和される。
【数５５】
2014206802.tif 000057
【００８２】
（６５）式をエンジン制御装置１０００で解くためには、Ｍ2、ｂ2、Λ及びνを予めオフラインで用意しておき、ペナルティ関数のｘmafについての一階偏微分∂（ｘ）と、二階偏微分∂2（ｘ）を予めオフラインで算出しておく。そして、エンジン制御装置１０００では、ｘbasemaf（ｋ）を（５６）式で算出するとともに、∂（ｘbasemaf（ｋ））及び∂2（ｘbasemaf（ｋ））を算出し、さらにｄ（ｘbasemaf（ｋ））＝（１  ０）Ｂｕ（ｋ−１）を算出する。そして（６５）式を解くことになる。
【００８３】
なお、上ではｘmafについて状態制約を設ける例を示したが、ｘmapについて状態制約を設けるようにしてもよいし、それら両方に状態制約を設けるようにしてもよい。基本的な処理は上で述べたとおりである。
【００８４】
次に、本実施の形態に係るオフライン処理装置１００の機能ブロック図を図３に示す。本実施の形態では、オフライン処理装置１００は、入力部１０１と、入力データ格納部１０２と、第１変形処理部１０３と、第１データ格納部１０４と、第２変形処理部１０５と、第２データ格納部１０６と、第３変形処理部１０７と、第３データ格納部１０８と、第１行列生成部１０９と、第４データ格納部１１０と、第２行列生成部１１１と、微分処理部１１２とを有する。
【００８５】
次に、オフライン処理装置１００の処理内容を図４を用いて説明する。まず、入力部１０１は、ユーザからの入力データとして、行列Ａ、Ｂ、Ｓ及びＲ、予測ホライズンＨp及び制御ホライズンＨcなどの値の入力を受け付け、入力データ格納部１０２に格納する（図４：ステップＳ１）。なお、入力データ格納部１０２に格納されていない場合には、（４６）乃至（４８）式のデータも入力部１０１を介して入力データ格納部１０２に格納されるものとする。同様に、入力制約に係る項の式のデータ（（２４）式）についても同様に、入力データ格納部１０２に格納されており、格納されていない場合には入力部１０１を介して入力される。さらに、状態制約に係る式のデータについても同様に、入力データ格納部１０２に格納されており、格納されていない場合には入力部１０１を介して入力される。
【００８６】
次に、第１変形処理部１０３は、入力データ格納部１０２に格納されているデータを用いて、状態偏差ｘについての第１変形処理を実行し、処理結果を第１データ格納部１０４に格納する（ステップＳ３）。（２２）式に由来する（２８−１）及び（２８−２）式を（２９）式に代入する処理をｉを増加させながら繰り返して実行し、（３０）式を生成する処理を実行する。なお、（３０）式については、（４５）式のように書き下される。
【００８７】
また、第２変形処理部１０５は、入力データ格納部１０２に格納されているデータを用いて、ラグランジュ乗数λについての第２変形処理を実行し、処理結果を第２データ格納部１０６に格納する（ステップＳ５）。（２０）式に由来する（３１）及び（３２）式を（３３）式に代入する処理をｉを減少させながら繰り返して実行し、（３４）式を生成する処理を実行する。
【００８８】
さらに、第３変形処理部１０７は、第１データ格納部１０４及び第２データ格納部１０６に格納されているデータを用いて、ラグランジュ乗数λについての第３変形処理を実行し、処理結果を第３データ格納部１０８に格納する（ステップＳ７）。（３０）式を（３４）式に代入する処理を実行して、（３５）式を生成する。
【００８９】
そして、第１行列生成部１０９は、（３０）式（第１データ格納部１０４に格納されている）及び（３５）式（第３データ格納部１０８に格納されている）を（２１）式（入力データ格納部１０２に格納されている）に代入して連立方程式（３６）式を生成し、これを（３７）式の形に変形することで、行列Ｍ2及びｂ2を算出し、第４データ格納部１１０に格納する（ステップＳ９）。ここまでの処理は、特開２０１２−１９４９６０号公報に記載の技術とほぼ同様である。
【００９０】
その後、第２行列生成部１１１は、第１データ格納部１０４に格納されている、ｘについての第１変形処理の結果（（３０）式）を用いて、入力データ格納部１０２に格納されている式に従って、係数行列Λ及び係数ベクトルνを生成し、第４データ格納部１１０に格納する（ステップＳ１１）。
【００９１】
さらに、微分処理部１１２は、第１データ格納部１０２に格納されている、状態制約についてのペナルティ関数のｘについての１階偏微分関数∂（ｘ）と、２階偏微分関数∂2（ｘ）を生成し、第４データ格納部１１０に格納する（ステップＳ１３）。１階偏微分関数及び２階偏微分関数の生成については、周知であるからここでは説明を省略する。
【００９２】
以上のような処理を実行することで、エンジン制御装置１０００の処理負荷を大幅に削減でき、処理速度を高速化できる。
【００９３】
次に、図５にエンジン制御装置１０００におけるモデル予測制御部２２０の構成例を示す。モデル予測制御部２２０は、微分値算出部２２０１と、行列算出部２２０２と、出力生成部２２０３と、データ格納部２２０５とを有する。
【００９４】
データ格納部２２０５には、第４データ格納部１１０に格納されているデータが格納される。
【００９５】
次に、図６を用いて、エンジン制御装置１０００におけるモデル予測制御部２２０の処理内容について説明する。
【００９６】
まず、モデル予測制御部２２０は、ｘ（ｋ）（＝Ｘref（ｋ）−ｙ（ｋ））を取得する（図６：ステップＳ２１）。
【００９７】
そうすると、モデル予測制御部２２０は、関数ベクトルｂ2にｘ（ｋ）を代入して関数ベクトルｂ2を数値ベクトル化する（ステップＳ２３）。
【００９８】
また、微分値算出部２２０１は、ｘbase（ｋ）（ｘbasemaf（ｋ）であれば（５６）式）を算出し、データ格納部２２０５に格納されているペナルティ関数の１階偏微分関数∂（ｘ）と、２階偏微分関数∂2（ｘ）にｘbase（ｋ）を代入し、∂（ｘbase（ｋ））及び∂2（ｘbase（ｋ））を求める（ステップＳ２５）。
【００９９】
そして、行列算出部２２０２は、Ｍ3及びｂ3を算出する（ステップＳ２７）。
【０１００】
（６５）式の左辺におけるベクトルＵ2以外の部分である（Ｍ2＋∂2（ｘbase（ｋ））Λ）がＭ3である。また、（６５）式の右辺がｂ3である。すなわち、ｂ3＝ｂ2−（∂（ｘbase（ｋ））−∂2（ｘbase（ｋ））ｄ（ｘbase（ｋ））νである。ｂ3を算出するためには、ｄ（ｘbase（ｋ））を、Ｂｕ（ｋ−１）で算出する。
【０１０１】
そして、出力生成部２２０３は、Ｕ＝Ｍ-13ｂ3によってベクトルＵを算出する（ステップＳ２９）。さらに、出力生成部２２０３は、ベクトルＵから線形化制御入力ｕ（ｋ）（ｉ＝０）を抽出し、出力する（ステップＳ３１）。
【０１０２】
その後、モデル予測制御部２２０は、処理が終了するのかを判断し（ステップＳ３３）、処理終了でなければ処理はステップＳ２１に戻り、処理終了であれば処理は終了する。
【０１０３】
以上のような処理で分かるように、オフライン処理によって算出されているデータを活用でき、Ｍ3のサイズも削減されているので、状態制約についてのペナルティ関数を導入することでステップＳ２５及びＳ２７が追加されていても、高速に線形化制御入力を得ることができるようになる。さらに、ペナルティ関数の近似も適切に行われており、高い精度で適切な制御がなされるようになる。
【０１０４】
なお、ディーゼルエンジンの吸気系制御では、図７（縦軸はＭＡＦを表し、横軸は時間を表す）に示すようにＭＡＦの状態ｙmafが目標値Ｘrefmafを下回ると大量の微粒子物質が発生し、吸気系のつまりの原因となることがある。そのため、上で述べたようにＭＡＦに対して下限値を設定して制御を行うこととする。
【０１０５】
例えば、以下のようなモデルを採用するものとする。
【数５６】
2014206802.tif 000058
【０１０６】
この時に、例えばｘmaf≧−１０という状態制約を設定する。そして、予測ホライズンＨp＝４ステップ、制御ホライズンＨc＝２ステップとする。この際、状態制約についてのペナルティ関数は、以下のような１階偏微分関数及び２階偏微分関数が用いられる。
【数５７】
2014206802.tif 000059
【０１０７】
また、以下に示すようなΛ及びνが得られる。
【数５８】
2014206802.tif 000060
【０１０８】
このようにすれば、大量の微粒子物質の発生を抑制させることができる。
【０１０９】
なお、エンジンの制御だけではなく、例えば自動車の運転支援のための制御アルゴリズムにおいて、ハンドル（ステアリング）の位置とアクセル開度から、自動車の位置と速度を制御する場合にも本実施の形態は適用可能である。例えば、自動車は白線を越えてはいけないため、白線の位置が状態制約となる。自動車のモデルは非線形であるため、区分的線形なモデルに近似することによって、本実施の形態を適用すれば、安全な運転支援が可能となる。
【０１１０】
さらに、空調機で室内の温度を使用電力が最小になるように制御する場合にも、本実施の形態を適用可能である。室内で、ある一定の温度区間外の室温となった時、室内にいる住人は不快に感じるため、それらの境界の温度が状態制約となる。本実施の形態を適用することで、不快温度にすることなく室内の温度の制御することができる。
【０１１１】
［実施の形態２］
ディーゼルエンジンは、ＥＧＲバルブとＶＮＴノズルを操作することでＭＡＦとＭＡＰを最適に制御することで、最適な燃焼状態を保ち、有害物質（例えばＮＯｘ、ＰＭ）の発生を抑制する。
【０１１２】
ＥＧＲ及びＶＮＴによってＭＡＦ及びＭＡＰを制御する際、過渡運転状態においてＭＡＦ状態が最適な状態（目標値）を上回るオーバーシュートを起こすことがあり、ＮＯｘ発生の原因となっている。図８は、縦軸がＭＡＦを表し、横軸は時間を表し、点線が目標値を表し、実線が実際のＭＡＦの値ｘmafを表す。この例でも、オーバーシュートが発生している。また、図８のように、目標値が下がっていく時でも、ＭＡＦの状態が目標値を上回った状態を続けることが多くなると、この点にてもＮＯｘ発生の要因となる。また、図７に示したように目標値を大きく下回ると、ＰＭが過剰に発生する。
【０１１３】
第１の実施の形態のように状態制約を取り扱うモデル予測制御を行えば、改善が見込まれるが、以下のような問題もある。
【０１１４】
（ａ）目標値が上がっていく運転状態（上りの状態）と下がっていく状態（下りの状態）、変わらない状態（定常状態）で状態制約の意味合いが変わってしまう。
すなわち、上りの状態の時は上限の設定が難しく、下りの時は目標値を上回る状態が続くことが多いため、現在の状態が制約を満たさない場合がある。
（ｂ）運転状態は運転手に依存するため，状態制約を切り替えて制御することが難しい。
（ｃ）外乱やノイズの影響で状態制約を超えてしまう可能性があり、最適解を求めることができなくなることがある。
【０１１５】
そこで、本実施の形態では、以下のような状態制約を導入してオーバーシュートやアンダーシュートを抑制する。
【数５９】
2014206802.tif 000061
γ1乃至γ4についてはチューニングパラメータである。また、本実施の形態ではＭＡＦについての状態制約として採用する。
【０１１６】
このペナルティ関数は、図９に示すような関数である。ｘ（ｉ）、ｘmin及びｘmaxは、目標値を基準にした値、すなわちｘ（ｉ）は状態偏差であり、ｘmin及びｘmaxはそれぞれ目標値から下又は上にどこまでずれても対処しないかを表しており、その範囲を出た場合には、その範囲内に戻ろうとする。図１０に示すように、目標値曲線の上側の上限曲線と、下側の下限曲線との間であればペナルティ関数の値は０であるが、それを外れると値が急激に大きくなるペナルティ関数を採用する。
【０１１７】
なお、上で述べたペナルティ関数は一例であり、（ｘ−ｘmax）や（ｘmin−ｘ）をｅのべき数としていたが、これらを基数としてｎ乗（ｎは３以上）するような関数（γmax（ｘ−ｘmax）n、γmin（ｘmin−ｘ）n）を採用してもよい。
【０１１８】
このような状態制約を導入したとしても、基本的に第１の実施の形態で説明したオフライン処理は同じである。但し、上限値を超えた場合の関数と、下限値を超えた場合の関数とでは形が異なるので、それぞれについてｘmafについての１階偏微分関数及び２階偏微分関数を算出する。すなわち、ステップＳ１３では、下限側の１階偏微分関数∂min（ｘ）及び２階偏微分関数∂2min（ｘ）と、上限側の１階偏微分関数∂max（ｘ）及び２階偏微分関数∂2max（ｘ）とを生成する。
【０１１９】
そして、ｘmaf（ｋ）の値に応じて、下限側の１階偏微分関数∂min（ｘ）及び２階偏微分関数∂2min（ｘ）と、上限側の１階偏微分関数∂max（ｘ）及び２階偏微分関数∂2max（ｘ）と、制約なしとを切り替えることになる。
【０１２０】
具体的に、図１１を用いてモデル予測制御部２２０の処理内容を説明する。
【０１２１】
まず、モデル予測制御部２２０は、ｘ（ｋ）（＝Ｘref（ｋ）−ｙ（ｋ））を取得する（図１１：ステップＳ４１）。
【０１２２】
そうすると、モデル予測制御部２２０は、関数ベクトルｂ2にｘ（ｋ）を代入して関数ベクトルｂ2を数値ベクトル化する（ステップＳ４３）。
【０１２３】
また、微分値算出部２２０１は、ｘ（ｋ）がペナルティ関数の適用範囲のうち下限値未満についての適用範囲０であるか、上限値を超える範囲についての適用範囲１であるか、それ以外の範囲である適用範囲２であるかを判断する（ステップＳ４５）。
【０１２４】
適用範囲１であれば、微分値算出部２２０１は、ｘbase（ｋ）（ｘbasemaf（ｋ）であれば（５６）式）を算出し、データ格納部２２０５に格納されているペナルティ関数の１階偏微分関数∂min（ｘ）と、２階偏微分関数∂2min（ｘ）にｘbase（ｋ）を代入し、∂min（ｘbase（ｋ））及び∂2min（ｘbase（ｋ））を求める（ステップＳ４７）。
【０１２５】
そして、行列算出部２２０２は、範囲１についてのＭ3及びｂ3を算出する（ステップＳ４９）。
【０１２６】
（６５）式の左辺におけるベクトルＵ2以外の部分である（Ｍ2＋∂2min（ｘbase（ｋ））Λ）がＭ3である。また、（６５）式の右辺がｂ3である。すなわち、ｂ3＝ｂ2−（（∂min（ｘbase（ｋ））−∂2min（ｘbase（ｋ））ｄ（ｘbase（ｋ）））νである。ｂ3を算出するためには、ｄ（ｘbase（ｋ））を、Ｂｕ（ｋ−１）で算出する。そして処理はステップＳ５７に移行する。
【０１２７】
一方、適用範囲０であれば、微分値算出部２２０１は、ｘbase（ｋ）（ｘbasemaf（ｋ）であれば（５６）式）を算出し、データ格納部２２０５に格納されているペナルティ関数の１階偏微分関数∂max（ｘ）と、２階偏微分関数∂2max（ｘ）にｘbase（ｋ）を代入し、∂max（ｘbase（ｋ））及び∂2max（ｘbase（ｋ））を求める（ステップＳ５１）。
【０１２８】
そして、行列算出部２２０２は、範囲０についてのＭ3及びｂ3を算出する（ステップＳ５３）。
【０１２９】
（６５）式の左辺におけるベクトルＵ2以外の部分である（Ｍ2＋∂2max（ｘbase（ｋ））Λ）がＭ3である。また、（６５）式の右辺がｂ3である。すなわち、ｂ3＝ｂ2−（（∂max（ｘbase（ｋ））−∂2max（ｘbase（ｋ））ｄ（ｘbase（ｋ）））νである。ｂ3を算出するためには、ｄ（ｘbase（ｋ））を、Ｂｕ（ｋ−１）で算出する。そして処理はステップＳ５７に移行する。
【０１３０】
さらに、適用範囲２であれば、微分値算出部２２０１は、Ｍ3＝Ｍ2及びｂ3＝ｂ2と設定する（ステップＳ５５）。状態制約についてのペナルティ関数は適用されない。その後処理はステップＳ５７に移行する。
【０１３１】
そして、出力生成部２２０３は、Ｕ＝Ｍ-13ｂ3によってベクトルＵを算出する（ステップＳ５７）。さらに、出力生成部２２０３は、ベクトルＵから線形化制御入力ｕ（ｋ）（ｉ＝０）を抽出し、出力する（ステップＳ５９）。
【０１３２】
その後、モデル予測制御部２２０は、処理が終了するのかを判断し（ステップＳ６１）、処理終了でなければ処理はステップＳ４１に戻り、処理終了であれば処理は終了する。
【０１３３】
このような処理を実行すれば、適切なペナルティ関数が適用されるようになり、オーバーシュートやアンダーシュートを効果的に抑制できるようになる。
【０１３４】
［実施の形態３］
一般にＭＡＦとＭＡＰの追従性にはトレードオフの関係があるため、ＭＡＦの追従性を上げるためにＭＡＰは悪化するという現象が発生する。そこで、以下のような制御を行うと効果的である。
【０１３５】
（Ａ）ＭＡＰがある程度目標値に近い場合は、できるだけＭＡＦを近づけたいので、ＭＡＦのペナルティ関数を加える。
（Ｂ）ＭＡＦはある程度目標値に近いがＭＡＰが目標値から大きく離れている場合、ＭＡＦにペナルティ関数を加えて得られる効果よりＭＡＰの悪化が大きくなることが多いので、ペナルティ関数を加えない。
（Ｃ）ＭＡＦ及びＭＡＰともに目標値から離れている場合には、ＭＡＦ重視であってもＭＡＰもある程度追従していなければＮＯｘ及びＰＭは減らないし、ＭＡＦ及びＭＡＰともに大きく変化させようとする場合、入力のバルブが全開もしくは全閉になることが多く、ペナルティ関数を入れたとしても入力制約の影響で効果が出ないことが多い。従って、計算量の観点もあり、ペナルティ関数を加えない。
【０１３６】
そこで、図１２に示すような領域分割を行って、ＭＡＦ及びＭＡＰがいずれの領域に入っているか否かを判断することで、制御を切り替えるものとする。図１２の例では、縦軸はＭＡＰの偏差ｘmapを表し、横軸はＭＡＦの偏差ｘmafを表す。そして、ＭＡＰの基準値Ａ及びＢを設け、ｘmaxとＭＡＰの基準値Ａ及びＢで囲まれる範囲を領域０と設定してＭＡＦの状態制約についての上限側のペナルティ関数を適用する。また、ｘminとＭＡＰの基準値Ａ及びＢで囲まれる範囲を領域１として設定してＭＡＦの状態制約についての下限側のペナルティ関数を適用する。それ以外の領域については、領域２として設定し、ペナルティ関数を適用しない。
【０１３７】
このような領域設定を行う場合には、第２の実施の形態における処理フローのステップＳ４５において、ＭＡＦ及びＭＡＰの状態偏差ｘ（ｋ）に基づき、いずれの領域にｘ（ｋ）が入っているかを判定するようにする。後の処理については図１１と同様である。
【０１３８】
このような制御の切り替えを行った結果の効果を、図１３Ａ及び図１３Ｂに示す。図１３Ａは、ＭＡＦの値の時間変化を表しており、図１３Ｂは、ＭＡＰの値の時間変化を表している。ＭＡＰについてはあまり変化していないが、ＭＡＦについては、部分Ｐではオーバーシュートが低減されており、部分Ｑでは目標値への追従性が向上している。これによって、ＮＯｘ及びＰＭ排出量のさらなる低減がなされる。
【０１３９】
以上本技術の実施の形態を説明したが、本技術はこれに限定されるものではない。例えば、上で述べた機能ブロック構成については、プログラムモジュール構成とは一致しない場合もある。また、処理結果が変わらない限り、処理フローについても、ステップの順番を入れ替えたり、複数のステップを並列実行するようにしてもよい。
【０１４０】
なお、上で述べたオフライン処理装置１００は、コンピュータ装置であって、図１４に示すように、メモリ２５０１とＣＰＵ２５０３とハードディスク・ドライブ（ＨＤＤ）２５０５と表示装置２５０９に接続される表示制御部２５０７とリムーバブル・ディスク２５１１用のドライブ装置２５１３と入力装置２５１５とネットワークに接続するための通信制御部２５１７とがバス２５１９で接続されている。オペレーティング・システム（ＯＳ：Operating System）及び本実施例における処理を実施するためのアプリケーション・プログラムは、ＨＤＤ２５０５に格納されており、ＣＰＵ２５０３により実行される際にはＨＤＤ２５０５からメモリ２５０１に読み出される。ＣＰＵ２５０３は、アプリケーション・プログラムの処理内容に応じて表示制御部２５０７、通信制御部２５１７、ドライブ装置２５１３を制御して、所定の動作を行わせる。また、処理途中のデータについては、主としてメモリ２５０１に格納されるが、ＨＤＤ２５０５に格納されるようにしてもよい。本技術の実施例では、上で述べた処理を実施するためのアプリケーション・プログラムはコンピュータ読み取り可能なリムーバブル・ディスク２５１１に格納されて頒布され、ドライブ装置２５１３からＨＤＤ２５０５にインストールされる。インターネットなどのネットワーク及び通信制御部２５１７を経由して、ＨＤＤ２５０５にインストールされる場合もある。このようなコンピュータ装置は、上で述べたＣＰＵ２５０３、メモリ２５０１などのハードウエアとＯＳ及びアプリケーション・プログラムなどのプログラムとが有機的に協働することにより、上で述べたような各種機能を実現する。
【０１４１】
また、上で述べたエンジン制御装置１０００は、コンピュータ装置であって、図１５に示すように、ＲＡＭ（Random Access Memory）４５０１とプロセッサ４５０３とＲＯＭ（Read Only Memory）４５０７とセンサ群４５１５とがバス４５１９で接続されている。本実施の形態における処理を実施するための制御プログラム及び存在している場合にはオペレーティング・システム（ＯＳ：Operating System））は、ＲＯＭ４５０７に格納されており、プロセッサ４５０３により実行される際にはＲＯＭ４５０７からＲＡＭ４５０１に読み出される。プロセッサ４５０３は、センサ群４５１５（例えば、吸気圧センサ及び新気量センサ。場合によっては燃料噴射量測定部及びエンジン回転数測定部など。）を制御して、測定値を取得する。また、処理途中のデータについては、ＲＡＭ４５０１に格納される。なお、プロセッサ４５０３は、ＲＯＭ４５０７を含む場合もあり、さらに、ＲＡＭ４５０１を含む場合もある。本技術の実施の形態では、上で述べた処理を実施するための制御プログラムは、コンピュータ読み取り可能なリムーバブル・ディスクに格納されて頒布され、ＲＯＭライタによってＲＯＭ４５０７に書き込まれる場合もある。このようなコンピュータ装置は、上で述べたプロセッサ４５０３、ＲＡＭ４５０１、ＲＯＭ４５０７などのハードウエアと制御プログラム（場合によってはＯＳも）とが有機的に協働することにより、上で述べたような各種機能を実現する。
【０１４２】
以上述べた本実施の形態をまとめると以下のようになる。
【０１４３】
本実施の形態の第１の態様に係るオフライン処理方法は、（Ａ）（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ制御対象の状態方程式に関する第１の関数との和を含む第２の関数をラグランジュ乗数について偏微分することで得られ且つ制御対象の状態偏差及び制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）第２の関数を制御対象の状態偏差について偏微分することで得られ且つ制御対象の状態偏差及びラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）第２の関数を制御対象の線形化制御入力について偏微分することで得られ且つラグランジュ乗数及び制御対象の状態偏差を変数として含む第３の一次式群とを格納するデータ格納部から第１の一次式群を読み出して、状態偏差の現在値を除き状態偏差の変数を消し込むように順次代入を行って、第１の一次式群を、状態偏差の現在値と線形化制御入力とを変数として含む第４の一次式群に変形し、（Ｂ）データ格納部から第２の一次式群を読み出して、ラグランジュ乗数の変数を消し込むように順次代入を行って、第２の一次式群を、状態偏差を変数として含む第５の一次式群に変形し、（Ｃ）データ格納部から第３の一次式群を読み出して、第４の一次式群及び第５の一次式群を代入することで、第３の一次式群を、状態偏差の現在値と線形化制御入力とを変数として含む第６の一次式群に変形し、（Ｄ）第６の一次式群を、第１の係数行列と線形化制御入力についてのベクトルとの積が、状態偏差の現在値の関数ベクトルと等式で表されるように変形し、（Ｅ）モデル予測制御の制御対象における状態制約に対応するペナルティ関数を状態偏差について１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、（Ｆ）第４の一次式群から線形化制御入力についての係数を抽出し、ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する処理を含む。
【０１４４】
このような処理をオフラインで実行しておくことにより、状態制約を導入したモデル予測制御を行うことができるようになる。
【０１４５】
上で述べたペナルティ関数が、状態偏差が第１の値以上第２の値以下で０となり、状態偏差が第１の値より小さい場合及び第２の値より大きい場合状態偏差の値に応じて指数関数的に増加する関数である場合もある。制御対象の特性に応じた適切な関数が用いられる。
【０１４６】
また、上で述べたペナルティ関数が、状態偏差が第１の値以上第２の値以下で０となり、状態偏差が第１の値より小さい場合に状態偏差の値に応じて指数関数的に増加する第１の関数と、第２の値より大きい場合状態偏差の値に応じて指数関数的に増加する第２の関数とで表される場合もある。この場合、第１及び第２の微分関数を、第１の関数及び第２の関数の各々について生成する場合もある。このようにすれば、制御実行時に第１の関数と第２の関数とを切り替えて用いることができるようになる。
【０１４７】
本実施の形態の第２の態様に係る制御方法は、（Ａ）モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、制御期間内における制御対象に対する線形化制御入力のベクトルとの積が関数ベクトルと等しくなるように予め算出されている第１の係数行列と、モデル予測制御の制御対象における状態制約に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、ペナルティ関数を状態偏差について１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から第１の式及び第２の式を読み出して、状態偏差の現在値から得られる基点の値を代入することで、第１の値及び第２の値を算出し、（Ｂ）データ格納部から第１の係数行列及び第２の係数行列を読み出して、第２の値と第２の係数行列との積と第１の係数行列との和である第１の行列を算出し、（Ｃ）データ格納部から関数ベクトルを読み出して、関数ベクトルに状態偏差の現在値を代入することで得られる第１のベクトルと、第１の値と第２の値との差に線形化制御入力の１単位時間前の値から得られる第３の値を乗じた結果と係数ベクトルとの積との和である第２のベクトルを算出し、（Ｄ）第１の行列の逆行列と第２のベクトルとの積により、線形化制御入力のベクトルを生成し、（Ｅ）線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する処理を含む。
【０１４８】
このような処理を実行することで、高速に状態制約を伴うモデル予測制御を実行できるようになる。
【０１４９】
上で述べた制御方法は、（Ｆ）状態偏差の現在値が所定範囲内であるか判断する処理をさらに含むようにしてもよい。この場合、状態偏差の現在値が所定範囲内であれば、第１の値及び第２の値を算出する処理と、第１の行列を算出する処理と、第２のベクトルを算出する処理と、線形化制御入力のベクトルを生成する処理と、線形化制御入力を抽出する処理とを実行するようにしてもよい。制御対象によっては、状態制約に基づく処理を行うことが好ましくない場合もあるためである。
【０１５０】
さらに、上で述べたデータ格納部が、状態偏差の現在値の値域に応じた第１の式及び第２の式を格納している場合もある。この場合、（ａ１）上記第１の値及び第２の値を算出する処理が、状態偏差の現在値が属する値域を特定し、当該特定された値域に応じた第１の式及び第２の式を読み出す処理を含むようにしてもよい。適切な式を状況に応じて切り分けて用いることで、適切な制御を行うことができるようになる。
【０１５１】
また、本制御方法が、上で述べた状態偏差の現在値が所定範囲外であれば、（Ｇ）データ格納部から関数ベクトル及び第１の係数行列を読み出して、関数ベクトルに状態偏差の現在値を代入することで得られる第１のベクトルと、第１の係数行列の逆行列との積により、線形化制御入力の第２のベクトルを生成し、（Ｈ）線形化制御入力の第２のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する処理をさらに含むようにしてもよい。状態制約を考慮しないモデル予測制御を高速に実行できる。
【０１５２】
なお、上記方法による処理をプロセッサ（又はコンピュータ）に行わせるためのプログラムを作成することができ、当該プログラムは、例えばフレキシブルディスク、ＣＤ−ＲＯＭ、光磁気ディスク、半導体メモリ、ハードディスク等のコンピュータ読み取り可能な記憶媒体又は記憶装置に格納される。尚、中間的な処理結果はメインメモリ等の記憶装置に一時保管される。
【０１５３】
以上の実施例を含む実施形態に関し、さらに以下の付記を開示する。
【０１５４】
（付記１）
（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ前記制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形し、
前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む第５の一次式群に変形し、
前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を代入することで、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形し、
前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形し、
前記モデル予測制御の制御対象についての制約条件に対応するペナルティ関数を前記状態偏差について１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、
前記第４の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する
処理を、コンピュータに実行させるためのプログラム。
【０１５５】
（付記２）
前記ペナルティ関数が、前記状態偏差が第１の値以上第２の値以下で０となり、前記状態偏差が前記第１の値より小さい場合及び前記第２の値より大きい場合前記状態偏差の値に応じて指数関数的に増加する関数である
付記１記載のプログラム。
【０１５６】
（付記３）
前記ペナルティ関数が、前記状態偏差が第１の値以上第２の値以下で０となり、前記状態偏差が前記第１の値より小さい場合に前記状態偏差の値に応じて指数関数的に増加する第１の関数と、前記第２の値より大きい場合前記状態偏差の値に応じて指数関数的に増加する第２の関数とで表され、
前記第１及び第２の微分関数を、前記第１の関数及び前記第２の関数の各々について生成する
付記１記載のプログラム。
【０１５７】
（付記４）
モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、前記モデル予測制御の制御対象における状態制約に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、前記ペナルティ関数を前記状態偏差について１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び第２の式を読み出して、前記状態偏差の現在値から得られる基点の値を代入することで、第１の値及び第２の値を算出し、
前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出し、
前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルと、前記第１の値と前記第２の値との差に前記線形化制御入力の１単位時間前の値から得られる第３の値を乗じた結果と前記係数ベクトルとの積との和である第２のベクトルを算出し、
前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成し、
前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
処理をプロセッサに実行させるためのプログラム。
【０１５８】
（付記５）
前記状態偏差の現在値が所定範囲内であるか判断する処理をさらに含み、
前記状態偏差の現在値が前記所定範囲内であれば、前記第１の値及び第２の値を算出する処理と、前記第１の行列を算出する処理と、前記第２のベクトルを算出する処理と、前記線形化制御入力のベクトルを生成する処理と、前記線形化制御入力を抽出する処理とを実行する
付記４記載のプログラム。
【０１５９】
（付記６）
前記データ格納部が、
前記状態偏差の現在値の値域に応じた前記第１の式及び第２の式を格納しており、
前記第１の値及び前記第２の値を算出する処理が、
前記状態偏差の現在値が属する値域を特定し、当該特定された値域に応じた前記第１の式及び前記第２の式を読み出す処理
を含む付記４又は５記載のプログラム。
【０１６０】
（付記７）
前記状態偏差の現在値が前記所定範囲外であれば、
前記データ格納部から前記関数ベクトル及び前記第１の係数行列を読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルと、前記第１の係数行列の逆行列との積により、前記線形化制御入力の第２のベクトルを生成し、
前記線形化制御入力の第２のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
処理をさらに含む付記３乃至５のいずれか１つ記載のプログラム。
【０１６１】
（付記８）
（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ前記制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形し、
前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む第５の一次式群に変形し、
前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を代入することで、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形し、
前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形し、
前記モデル予測制御の制御対象についての制約条件に対応するペナルティ関数を前記状態偏差について１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、
前記第４の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する
処理を含み、コンピュータにより実行されるオフライン処理方法。
【０１６２】
（付記９）
モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、前記モデル予測制御の制御対象における状態制約に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、前記ペナルティ関数を前記状態偏差について１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び第２の式を読み出して、前記状態偏差の現在値から得られる基点の値を代入することで、第１の値及び第２の値を算出し、
前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出し、
前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルと、前記第１の値と前記第２の値との差に前記線形化制御入力の１単位時間前の値から得られる第３の値を乗じた結果と前記係数ベクトルとの積との和である第２のベクトルを算出し、
前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成し、
前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
処理を含み、プロセッサにより実行される制御方法。
【０１６３】
（付記１０）
（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ前記制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形する手段と、
前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む第５の一次式群に変形する手段と、
前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を代入することで、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形する手段と、
前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形する手段と、
前記モデル予測制御の制御対象についての制約条件に対応するペナルティ関数を前記状態偏差について１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成する手段と、
前記第４の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する手段と、
を有するオフライン処理装置。
【０１６４】
（付記１１）
モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、前記モデル予測制御の制御対象における状態制約に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、前記ペナルティ関数を前記状態偏差について１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び第２の式を読み出して、前記状態偏差の現在値から得られる基点の値を代入することで、第１の値及び第２の値を算出する第１の算出部と、
前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出し、前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルと、前記第１の値と前記第２の値との差に前記線形化制御入力の１単位時間前の値から得られる第３の値を乗じた結果と前記係数ベクトルとの積との和である第２のベクトルを算出する第２の算出部と、
前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成する生成部と、
前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を出力する出力部と、
を有する制御装置。
【符号の説明】
【０１６５】
１０１  入力部
１０２  入力データ格納部
１０３  第１変形処理部
１０４  第１データ格納部
１０５  第２変形処理部
１０６  第２データ格納部
１０７  第３変形処理部
１０８  第３データ格納部
１０９  第１行列生成部
１１０  第４データ格納部
１１１  第２行列生成部
１１２  微分処理部
【図１】
2014206802.tif 000062
【図２】
2014206802.tif 000063
【図３】
2014206802.tif 000064
【図４】
2014206802.tif 000065
【図５】
2014206802.tif 000066
【図６】
2014206802.tif 000067
【図７】
2014206802.tif 000068
【図８】
2014206802.tif 000069
【図９】
2014206802.tif 000070
【図１０】
2014206802.tif 000071
【図１１】
2014206802.tif 000072
【図１２】
2014206802.tif 000073
【図１３Ａ】
2014206802.tif 000074
【図１３Ｂ】
2014206802.tif 000075
【図１４】
2014206802.tif 000076
【図１５】
2014206802.tif 000077