JP6136494B2 - 情報処理方法、装置及びプログラム - Google Patents

Description

本技術は、モデル予測制御に関する。

モデル予測制御は、各時刻において有限時間未来までの応答を最適化することによって制御対象への入力を決定する制御方法である。操作量飽和などの不等式制約条件を陽に扱うことができる制御手法として高い制御性能を期待されている。

このモデル予測制御の適用については様々な分野について考察されているが、ディーゼルエンジンの吸気系制御への適用も考察されている。ディーゼルエンジンの吸気系は、ＥＧＲ（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とは加算されて操作量が生成され、当該操作量が飽和要素２３０を介してプラント２４０に入力される。このような処理が繰り返される。

このようなエンジン本体１の吸気系システムは、ＥＧＲバルブ及びＶＮＴノズルについて非線形であり、エンジンのモデルとしても、エンジン回転数ω及び燃料噴射量σに応じて非線形に変化する。

ここで、エンジン回転数ω及び燃料噴射量σが固定であるとして、エンジンモデルは、以下のように定式化される。

ここでＡ_ω(k),σ(k)及びＢ_ω(k),σ(k)は、時刻ｋにおける運転条件（ω（ｋ），σ（ｋ））によって変化する非線形の２×２行列関数である。なお、時刻ｋにおける回転数をω（ｋ）と表し、燃料噴射量をσ（ｋ）と表す。また、時刻ｋにおけるＭＡＦの状態偏差をｘ_maf（ｋ）と表し、ＭＡＰの状態偏差をｘ_map（ｋ）と表す。また、時刻ｋにおけるＭＡＦの測定値をｙ_maf（ｋ）と表し、ＭＡＰの測定値をｙ_map（ｋ）と表す。さらに、時刻ｋにおけるＥＧＲバルブの状態をｖ_egr（ｋ）と表し、ＶＮＴノズルの状態をｖ_vnt（ｋ）と表す。また、吸気圧ＭＡＰの目標値をＸ_map ^refと表し、新気量ＭＡＦの目標値をＸ_maf ^refと表す。

従って、状態偏差ベクトルｘ（ｋ）、測定値ベクトルｙ（ｋ）、目標値ベクトルＸ^ref（ｋ）及び状態ベクトルｖ（ｋ）は、以下のように表される。

ｘ（ｋ）、ｙ（ｋ）及びＸ^ref（ｋ）の関係は以下のようになっている。

また、非線形関数ｇを以下のように表す。

ｇ₁及びｇ₂は、単調増加関数である。実験データからすると、このような仮定は妥当性がある。非線型モデルでは計算速度等から取り扱いが困難なため、以下では、非線形なエンジンモデルを線型モデルに変換する。そこで、上で述べた線形化制御入力ｕ_egr及びｕ_vntを導入することで、線形化する。

なお、線形化制御入力ベクトルｕ（ｋ）を以下のように表す。

本実施の形態では、典型的な運転条件を選択して、セットポイントと呼ぶことにする。そして、モデル予測制御部２２０は、運転条件（ω（ｋ），σ（ｋ））を取得すると、対応するセットポイントについてのω及びσについての行列Ａ_ω,σ及びＢ_ω,σを選択し、以下で述べる処理を行って線形化制御入力ｕを算出する。また、変換関数２２５は、線形化制御入力ｕを、変換関数ｇ^-1により非線形の最適状態ｖに変換する。この最適状態ｖと状態目標値Ｖrefとは加算されて操作量が生成され、エンジンのプラント２４０に出力される。

モデル予測制御部２２０は、セットポイントについて最適な線形化制御入力ｕを、１サンプル時間毎に、以下の線形モデルに従って算出する。

但し、線形化制御入力ｕについての制約条件は以下のとおりである。

例えば、ω＝２０００ｒｐｍ、σ＝２０ｍｍ³／ｓｔの場合、以下のような行列が用いられる。

各時刻ｋについての線形化制御入力ｕは、以下に示すような評価関数Ｊを最小化することで決定される。なお、予測ホライズンＨp及び制御ホライズンＨcは、予め設定されるものとする。

ここで、Ｑ^Tは行列Ｑの転置行列である。また、Ｑ∈Ｍ（ｎ，Ｒ）、Ｒ∈Ｍ（ｎ，Ｒ）及びＳ∈Ｍ（ｎ，Ｒ）は、重み行列であり、予め設定される。なお、Ｍ（ｎ，Ｒ）はｎ×ｎ行列を表す。（９）式において、無限ホライズンについての最適化問題を、有限ホライズンについての最適化問題にしており、このため、Ricatti代数方程式を解く行列として重み行列Ｓを用いている。

本実施の形態では、モデル予測制御の問題は、有限ホライズンの最適制御問題を解くことに変形される。

ここで、ｘは、現在の状態偏差ベクトルを表す。

最適解を求める手法として、ラグランジュの未定乗数法が効率的であるが、不等式制約がある場合にはこの手法をそのまま用いることができない。従って、不等式制約の代わりにペナルティ関数を導入して、不等式制約のない評価関数に緩和する。

すなわち、以下のような評価関数Ｊbに緩和する。

ここで、ｐ（ｕ（ｋ＋ｉ））はペナルティ関数であり、入力値ｕ（ｋ＋ｉ）が、可能な値域の境界に近づくと、大きな値を評価関数Ｊに加えることになる。

一般的に、不等式制約ｕ≦ｕ^maxに対して用いられるペナルティ関数としては以下のようなものがある。

ｒは、ペナルティ関数の係数であり、不等式制約に対する重みを表している。

本実施の形態では、以下の演算のためには、（ｃ）のようなシンプルなペナルティ関数が好ましいので、（ｃ）を用いることにする。

ディーゼルエンジンの場合、線形化制御入力ｕには以下のような制約が存在する。

（ｃ）のような形のペナルティ関数を採用する場合には、ペナルティ関数は以下のように表される。

（１５）式を用いて制約条件を緩和する場合、以下に示すような、有限ホライズンについての最適制御問題を解くことに帰着する。

このような最適制御問題を解くためには、ラグランジュの未定乗数法が用いられる。ここで、ラグランジュ乗数λ（ｊ）＝［λ_maf（ｊ），λ_map（ｊ）］^Tを用い、ハミルトニアンＨを以下のように表す。

ここで、Х（χの大文字）は、以下のように定義される。

∇_x,u,λＨ＝０と置けば、以下のような方程式が得られる。

なお、（２０）式において∂Ｊb／∂ｘ＝∂Ｊ／∂ｘ＋∂ｐ（ｕ）／∂ｘ＝∂Ｊ／∂ｘと変形されている。

なお、∂／∂ｘ，∂／∂ｕ，∂／∂λは、以下のような演算である。

なお、方程式（２０）乃至（２２）は、∂ｐ／∂ｕの部分を除き線形方程式となっている。∂ｐ／∂ｕを、直前のステップの最適線形化制御入力ｕ~（”~”はｕの上に付される）でテイラー展開することで線形化すると、以下のように表される。

（２３）式の第１項のＥＧＲ成分は、以下のように表される。なお、ＶＮＴ成分も同様の形で表される。

このような線形近似を行えば、（２０）乃至（２２）式は線形となる。そうすると、線形近似された（２０）乃至（２２）式は以下のように表される。

ここでＭは、（２Ｈ_p＋Ｈ_c）×（２Ｈ_p＋Ｈ_c）の行列であり、Ｕ及びｂは、（２Ｈ_p＋Ｈ_c）ベクトルである。ベクトルＵは、以下のように表される。

また、ベクトルｂは、以下のように表される。なお、ｂ_maf、ｂ_map及びｂ_constは、行列の変換によって算出される（２Ｈ_p＋Ｈ_c）行列である。
ｂ＝ｘ_mafｂ_maf＋ｘ_mapｂ_map＋ｂ_const （２７）

最適な線形化制御入力は、この線形方程式を解くことによって得られる。しかしながら、行列Ｍのサイズが大きければ計算時間は長くなってしまうので、行列Ｍのサイズをオフラインの前処理にて削減する。

まず、（２２）式の第１乃至第３の方程式を具体的に書き下すと以下のようになる。

（２８−１）式及び（２８−２）式の状態偏差ｘ_maf（ｋ）及びｘ_map（ｋ）は、現在値である。従って、（２８−１）式によれば、ｘ_maf（ｋ＋１）は線形化制御入力ｕ（ｋ）の１次式で表される。（２８−２）式によれば、ｘ_map（ｋ＋１）も、線形化制御入力ｕ（ｋ）の１次式で表される。また、（２９）式によれば、ｘ_maf（ｋ＋２）は、状態偏差及び線形化制御入力の１次式で表される。しかしながら、（２８−１）式及び（２８−２）式を（２９）式に代入すると、ｘ_maf（ｋ＋２）を、線形化制御入力ｕ（ｋ）及びｕ（ｋ＋１）で表すことができるようになる。これをｉを増加させつつ繰り返せば、状態偏差ｘ_maf（ｋ＋ｉ）及びｘ_map（ｋ＋ｉ）は、線形化制御入力の１次式として以下のように表される。

一方、（２０）式の（２Ｈ_p−１）番目及び（２Ｈ_p−３）番目の方程式を具体的に書き下すと以下のようになる。

ここで、ｓ₁₁、ｓ₂₂、ｒ₁₁は、重み行列Ｓ及びＲの行列成分を表す。

この（３１）式及び（３２）式は、λ_maf（Ｈ_p）及びλ_map（Ｈ_p）の一次式で表されることが分かる。そこで、（３１）式及び（３２）式を、（３３）式に代入すると、λ_maf（Ｈ_p−１）が、目標値Ｘ^ref _maf、状態偏差ｘ（ｋ＋Ｈ_p）及びｘ（ｋ＋Ｈ_p−１）で表すことができる。このような処理をｉを減少させつつ順次行うことで、ラグランジュ乗数λ（ｋ＋ｉ）は、以下のように表される。

さらに、（３０）式を（３４）式に代入することで、ラグランジュ乗数λ（ｋ＋ｉ）は、線形化制御入力ｕ（ｋ）の一次式に変形される。

最後に、（３０）式及び（３５）式を、（２１）式に代入すると、以下のように表される。

このように、（２１）式は、状態偏差の初期値ｘ（ｋ）及び線形化制御入力ｕの一次式で表されることができるようになる。そうすると、（３６）式は、以下の形式に直すことができる。

ここでＭ₂はＨ_c×Ｈ_c行列である。Ｕ₂及びｂ₂は、以下のような形のベクトルである。

なお、ｂ_2maf、ｂ_2map及びｂ_2constは、行列の変換によって算出される（２Ｈ_p＋Ｈ_c）行列である。

このような処理を数式処理を用いてオフラインで行っておけば、エンジン制御装置１０００では、（３７）式を解く演算を行えばよい。行列Ｍ₂のサイズは、行列Ｍよりも小さくなっており、計算量も削減されるため、高速に制御入力を得ることができるようになる。

次に、定常偏差を削減するために、状態制約を導入する場合を考える。説明を簡単にするために、ＭＡＦの状態偏差について以下のような制約を導入する。ＭＡＰについてはここでは取り扱わないが、同様に設定する場合もある。

しかしながら、このままではモデル予測制御問題が不等式制約を含む問題となり、現実的には解けないので、例えば、以下に示すようなペナルティ関数を導入することで、不等式制約問題を等式制約問題に緩和する。

ここでは、評価関数Ｊbに対して（４２）式のペナルティ関数を加える。そうすると、以下のように表される。

そして、（３０）式からすると、ｘ（ｋ＋ｉ）は、以下のように書き下される。

そうすると、ｘ（ｋ＋ｉ）は、全てのｉについて、ｕ（０）乃至ｕ（Ｈ_c−１）の関数であるので、（４２）式のペナルティ関数についてもｕ（０）乃至ｕ（Ｈ_c−１）の関数と考えることができる。よって、入力制約について（２０）式乃至（２２）式で記述された最適化問題は、以下のように変形される。

このような（４６）式乃至（４８）式を解くことになるが、そのためには状態についてのペナルティ関数に関する∂ｑ（ｘ）／∂ｕを線形化する。

まず、∂ｑ（ｘ_maf（ｋ＋ｉ））／∂ｕ_egr（ｋ＋ｊ）は以下のように表される。

ここで（４５）式があるので、（４９）式は、以下のように表される。

次に、（５１）式中の∂ｑ（ｘ_maf（ｋ＋ｉ））／∂ｘ_maf（ｋ＋ｉ）を線形化する。その際、以下の最適な線形化制御入力ｕ~で１次のテイラー展開を行うと、∂ｑ（ｘ_maf（ｋ＋ｉ））／∂ｘ_maf（ｋ＋ｉ）は以下のように表される。

また、（５１）式と同様に、以下のように表される。

従って、以下のような近似式を得られる。

よって、（５１）式及び（５４）式を組み合わせると、以下のような近似式が得られる。

ここで、ｘ~_maf（ｋ＋ｉ）（”~”はｘの上に載っているものとする）を、ｕ~に対応する状態偏差ベースポイントと設定する。そして、短い時間における状態偏差の値の変動は小さいので、ｘ^base _maf（ｋ）は以下のように表される。

そして、全てのｉについてｘ~_maf（ｋ＋ｉ）をｘ^base _maf（ｋ）と近似する。ｘ（ｋ）は、上でも述べたように現在の状態偏差の値を表しており、ｕ*（ｋ−１）は、１単位時間前の最適な線形化制御入力である。

そうすると、（４５）式から、（５５）式の定数部分は、以下のように表される。

また、ｘ~_maf（ｋ）をｘ^base _maf（ｋ）と近似する。そうすると、（５６）式から、以下のような近似式が得られる。

そこで、
ｄ（ｘ^base _maf（ｋ））＝（１ ０）Ｂｕ（ｋ−１） （１０１）と置く。

さらに、ｘ~_maf（ｋ＋ｉ）をｘ^base _maf（ｋ）と近似すれば、｛∂²ｑ（ｘ_maf（ｋ＋ｉ））／∂ｘ_maf（ｋ＋ｉ）²｝（ｕ~）は、全てのｉについてベース状態偏差ｘ^base _maf（ｋ）に関して同じ関数になる。同様に、｛∂ｑ（ｘ_maf（ｋ＋ｉ））／∂ｘ_maf（ｋ＋ｉ）｝（ｕ~）は、全てのｉについてベース状態偏差ｘ^base _maf（ｋ）に関して同じ関数になる。

以下、ペナルティ関数のｘ（（４２）式であればｘ_maf）についての一階偏微分を∂（ｘ）と表し、二階偏微分を∂²（ｘ）と表し、ｘ＝ｘ^base _maf（ｋ）についての∂（ｘ^base _maf（ｋ））及び∂²（ｘ^base _maf（ｋ））によって、（５５）式は以下のように表される。

この（６２）式から、以下のような近似も行われる。

さらに、この（６３）式から、以下のような近似も行われる。

ここで、Λ及びνは以下のように表される。

なお、Λにおけるｐは２ｊ＋１で固定で、ｑについては、１からＨ_cまで増加する。すなわち、Λは、１行Ｈ_c列の行列である。

従って、（６４）式及び（３７）式から、状態についてのペナルティ関数を伴う最適化問題は、以下のように緩和される。

（６５）式をエンジン制御装置１０００で解くためには、Ｍ₂、ｂ₂、Λ及びνを予めオフラインで用意しておき、ペナルティ関数のｘ_mafについての一階偏微分∂（ｘ）と、二階偏微分∂²（ｘ）を予めオフラインで算出しておく。そして、エンジン制御装置１０００では、ｘ^base _maf（ｋ）を（５６）式で算出するとともに、∂（ｘ^base _maf（ｋ））及び∂²（ｘ^base _maf（ｋ））を算出し、さらにｄ（ｘ^base _maf（ｋ））＝（１ ０）Ｂｕ（ｋ−１）を算出する。そして（６５）式を解くことになる。

なお、上ではｘ_mafについて状態制約を設ける例を示したが、ｘ_mapについて状態制約を設けるようにしてもよいし、それら両方に状態制約を設けるようにしてもよい。基本的な処理は上で述べたとおりである。

次に、本実施の形態に係るオフライン処理装置１００の機能ブロック図を図３に示す。本実施の形態では、オフライン処理装置１００は、入力部１０１と、入力データ格納部１０２と、第１変形処理部１０３と、第１データ格納部１０４と、第２変形処理部１０５と、第２データ格納部１０６と、第３変形処理部１０７と、第３データ格納部１０８と、第１行列生成部１０９と、第４データ格納部１１０と、第２行列生成部１１１と、微分処理部１１２とを有する。

次に、オフライン処理装置１００の処理内容を図４を用いて説明する。まず、入力部１０１は、ユーザからの入力データとして、行列Ａ、Ｂ、Ｓ及びＲ、予測ホライズンＨ_p及び制御ホライズンＨ_cなどの値の入力を受け付け、入力データ格納部１０２に格納する（図４：ステップＳ１）。なお、入力データ格納部１０２に格納されていない場合には、（４６）乃至（４８）式のデータも入力部１０１を介して入力データ格納部１０２に格納されるものとする。同様に、入力制約に係る項の式のデータ（（２４）式）についても同様に、入力データ格納部１０２に格納されており、格納されていない場合には入力部１０１を介して入力される。さらに、状態制約に係る式のデータについても同様に、入力データ格納部１０２に格納されており、格納されていない場合には入力部１０１を介して入力される。

次に、第１変形処理部１０３は、入力データ格納部１０２に格納されているデータを用いて、状態偏差ｘについての第１変形処理を実行し、処理結果を第１データ格納部１０４に格納する（ステップＳ３）。（２２）式に由来する（２８−１）及び（２８−２）式を（２９）式に代入する処理をｉを増加させながら繰り返して実行し、（３０）式を生成する処理を実行する。なお、（３０）式については、（４５）式のように書き下される。

また、第２変形処理部１０５は、入力データ格納部１０２に格納されているデータを用いて、ラグランジュ乗数λについての第２変形処理を実行し、処理結果を第２データ格納部１０６に格納する（ステップＳ５）。（２０）式に由来する（３１）及び（３２）式を（３３）式に代入する処理をｉを減少させながら繰り返して実行し、（３４）式を生成する処理を実行する。

さらに、第３変形処理部１０７は、第１データ格納部１０４及び第２データ格納部１０６に格納されているデータを用いて、ラグランジュ乗数λについての第３変形処理を実行し、処理結果を第３データ格納部１０８に格納する（ステップＳ７）。（３０）式を（３４）式に代入する処理を実行して、（３５）式を生成する。

そして、第１行列生成部１０９は、（３０）式（第１データ格納部１０４に格納されている）及び（３５）式（第３データ格納部１０８に格納されている）を（２１）式（入力データ格納部１０２に格納されている）に代入して連立方程式（３６）式を生成し、これを（３７）式の形に変形することで、行列Ｍ₂及びｂ₂を算出し、第４データ格納部１１０に格納する（ステップＳ９）。ここまでの処理は、特開２０１２−１９４９６０号公報に記載の技術とほぼ同様である。

その後、第２行列生成部１１１は、第１データ格納部１０４に格納されている、ｘについての第１変形処理の結果（（３０）式）を用いて、入力データ格納部１０２に格納されている式に従って、係数行列Λ及び係数ベクトルνを生成し、第４データ格納部１１０に格納する（ステップＳ１１）。

さらに、微分処理部１１２は、第１データ格納部１０２に格納されている、状態制約についてのペナルティ関数のｘについての１階偏微分関数∂（ｘ）と、２階偏微分関数∂²（ｘ）を生成し、第４データ格納部１１０に格納する（ステップＳ１３）。１階偏微分関数及び２階偏微分関数の生成については、周知であるからここでは説明を省略する。

以上のような処理を実行することで、エンジン制御装置１０００の処理負荷を大幅に削減でき、処理速度を高速化できる。

次に、図５にエンジン制御装置１０００におけるモデル予測制御部２２０の構成例を示す。モデル予測制御部２２０は、微分値算出部２２０１と、行列算出部２２０２と、出力生成部２２０３と、データ格納部２２０５とを有する。

データ格納部２２０５には、第４データ格納部１１０に格納されているデータが格納される。

次に、図６を用いて、エンジン制御装置１０００におけるモデル予測制御部２２０の処理内容について説明する。

まず、モデル予測制御部２２０は、ｘ（ｋ）（＝Ｘ_ref（ｋ）−ｙ（ｋ））を取得する（図６：ステップＳ２１）。

そうすると、モデル予測制御部２２０は、関数ベクトルｂ₂にｘ（ｋ）を代入して関数ベクトルｂ₂を数値ベクトル化する（ステップＳ２３）。

また、微分値算出部２２０１は、ｘ^base（ｋ）（ｘ^base _maf（ｋ）であれば（５６）式）を算出し、データ格納部２２０５に格納されているペナルティ関数の１階偏微分関数∂（ｘ）と、２階偏微分関数∂²（ｘ）にｘ^base（ｋ）を代入し、∂（ｘ^base（ｋ））及び∂²（ｘ^base（ｋ））を求める（ステップＳ２５）。

そして、行列算出部２２０２は、Ｍ₃及びｂ₃を算出する（ステップＳ２７）。

（６５）式の左辺におけるベクトルＵ₂以外の部分である（Ｍ₂＋∂²（ｘ^base（ｋ））Λ）がＭ₃である。また、（６５）式の右辺がｂ₃である。すなわち、ｂ₃＝ｂ₂−（∂（ｘ^base（ｋ））−∂²（ｘ^base（ｋ））ｄ（ｘ^base（ｋ））νである。ｂ₃を算出するためには、ｄ（ｘ^base（ｋ））を、Ｂｕ（ｋ−１）で算出する。

そして、出力生成部２２０３は、Ｕ＝Ｍ^-1 ₃ｂ₃によってベクトルＵを算出する（ステップＳ２９）。さらに、出力生成部２２０３は、ベクトルＵから線形化制御入力ｕ（ｋ）（ｉ＝０）を抽出し、出力する（ステップＳ３１）。

その後、モデル予測制御部２２０は、処理が終了するのかを判断し（ステップＳ３３）、処理終了でなければ処理はステップＳ２１に戻り、処理終了であれば処理は終了する。

以上のような処理で分かるように、オフライン処理によって算出されているデータを活用でき、Ｍ₃のサイズも削減されているので、状態制約についてのペナルティ関数を導入することでステップＳ２５及びＳ２７が追加されていても、高速に線形化制御入力を得ることができるようになる。さらに、ペナルティ関数の近似も適切に行われており、高い精度で適切な制御がなされるようになる。

なお、ディーゼルエンジンの吸気系制御では、図７（縦軸はＭＡＦを表し、横軸は時間を表す）に示すようにＭＡＦの状態ｙ_mafが目標値Ｘ^ref _mafを下回ると大量の微粒子物質が発生し、吸気系のつまりの原因となることがある。そのため、上で述べたようにＭＡＦに対して下限値を設定して制御を行うこととする。

例えば、以下のようなモデルを採用するものとする。

この時に、例えばｘ_maf≧−１０という状態制約を設定する。そして、予測ホライズンＨ_p＝４ステップ、制御ホライズンＨ_c＝２ステップとする。この際、状態制約についてのペナルティ関数は、以下のような１階偏微分関数及び２階偏微分関数が用いられる。

また、以下に示すようなΛ及びνが得られる。

このようにすれば、大量の微粒子物質の発生を抑制させることができる。

なお、エンジンの制御だけではなく、例えば自動車の運転支援のための制御アルゴリズムにおいて、ハンドル（ステアリング）の位置とアクセル開度から、自動車の位置と速度を制御する場合にも本実施の形態は適用可能である。例えば、自動車は白線を越えてはいけないため、白線の位置が状態制約となる。自動車のモデルは非線形であるため、区分的線形なモデルに近似することによって、本実施の形態を適用すれば、安全な運転支援が可能となる。

さらに、空調機で室内の温度を使用電力が最小になるように制御する場合にも、本実施の形態を適用可能である。室内で、ある一定の温度区間外の室温となった時、室内にいる住人は不快に感じるため、それらの境界の温度が状態制約となる。本実施の形態を適用することで、不快温度にすることなく室内の温度の制御することができる。

［実施の形態２］
ディーゼルエンジンは、ＥＧＲバルブとＶＮＴノズルを操作することでＭＡＦとＭＡＰを最適に制御することで、最適な燃焼状態を保ち、有害物質（例えばＮＯｘ、ＰＭ）の発生を抑制する。

ＥＧＲ及びＶＮＴによってＭＡＦ及びＭＡＰを制御する際、過渡運転状態においてＭＡＦ状態が最適な状態（目標値）を上回るオーバーシュートを起こすことがあり、ＮＯｘ発生の原因となっている。図８は、縦軸がＭＡＦを表し、横軸は時間を表し、点線が目標値を表し、実線が実際のＭＡＦの値ｘ_mafを表す。この例でも、オーバーシュートが発生している。また、図８のように、目標値が下がっていく時でも、ＭＡＦの状態が目標値を上回った状態を続けることが多くなると、この点にてもＮＯｘ発生の要因となる。また、図７に示したように目標値を大きく下回ると、ＰＭが過剰に発生する。

第１の実施の形態のように状態制約を取り扱うモデル予測制御を行えば、改善が見込まれるが、以下のような問題もある。

（ａ）目標値が上がっていく運転状態（上りの状態）と下がっていく状態（下りの状態）、変わらない状態（定常状態）で状態制約の意味合いが変わってしまう。
すなわち、上りの状態の時は上限の設定が難しく、下りの時は目標値を上回る状態が続くことが多いため、現在の状態が制約を満たさない場合がある。
（ｂ）運転状態は運転手に依存するため，状態制約を切り替えて制御することが難しい。
（ｃ）外乱やノイズの影響で状態制約を超えてしまう可能性があり、最適解を求めることができなくなることがある。

そこで、本実施の形態では、以下のような状態制約を導入してオーバーシュートやアンダーシュートを抑制する。

γ₁乃至γ₄についてはチューニングパラメータである。また、本実施の形態ではＭＡＦについての状態制約として採用する。

このペナルティ関数は、図９に示すような関数である。ｘ（ｉ）、ｘ_min及びｘ_maxは、目標値を基準にした値、すなわちｘ（ｉ）は状態偏差であり、ｘ_min及びｘ_maxはそれぞれ目標値から下又は上にどこまでずれても対処しないかを表しており、その範囲を出た場合には、その範囲内に戻ろうとする。図１０に示すように、目標値曲線の上側の上限曲線と、下側の下限曲線との間であればペナルティ関数の値は０であるが、それを外れると値が急激に大きくなるペナルティ関数を採用する。

なお、上で述べたペナルティ関数は一例であり、（ｘ−ｘ_max）や（ｘ_min−ｘ）をｅのべき数としていたが、これらを基数としてｎ乗（ｎは３以上）するような関数（γ_max（ｘ−ｘ_max）ⁿ、γ_min（ｘ_min−ｘ）ⁿ）を採用してもよい。

このような状態制約を導入したとしても、基本的に第１の実施の形態で説明したオフライン処理は同じである。但し、上限値を超えた場合の関数と、下限値を超えた場合の関数とでは形が異なるので、それぞれについてｘ_mafについての１階偏微分関数及び２階偏微分関数を算出する。すなわち、ステップＳ１３では、下限側の１階偏微分関数∂min（ｘ）及び２階偏微分関数∂²min（ｘ）と、上限側の１階偏微分関数∂max（ｘ）及び２階偏微分関数∂²max（ｘ）とを生成する。

そして、ｘ_maf（ｋ）の値に応じて、下限側の１階偏微分関数∂min（ｘ）及び２階偏微分関数∂²min（ｘ）と、上限側の１階偏微分関数∂max（ｘ）及び２階偏微分関数∂²max（ｘ）と、制約なしとを切り替えることになる。

具体的に、図１１を用いてモデル予測制御部２２０の処理内容を説明する。

まず、モデル予測制御部２２０は、ｘ（ｋ）（＝Ｘ_ref（ｋ）−ｙ（ｋ））を取得する（図１１：ステップＳ４１）。

そうすると、モデル予測制御部２２０は、関数ベクトルｂ₂にｘ（ｋ）を代入して関数ベクトルｂ₂を数値ベクトル化する（ステップＳ４３）。

また、微分値算出部２２０１は、ｘ（ｋ）がペナルティ関数の適用範囲のうち下限値未満についての適用範囲０であるか、上限値を超える範囲についての適用範囲１であるか、それ以外の範囲である適用範囲２であるかを判断する（ステップＳ４５）。

適用範囲１であれば、微分値算出部２２０１は、ｘ^base（ｋ）（ｘ^base _maf（ｋ）であれば（５６）式）を算出し、データ格納部２２０５に格納されているペナルティ関数の１階偏微分関数∂min（ｘ）と、２階偏微分関数∂²min（ｘ）にｘ^base（ｋ）を代入し、∂min（ｘ^base（ｋ））及び∂²min（ｘ^base（ｋ））を求める（ステップＳ４７）。

そして、行列算出部２２０２は、範囲１についてのＭ₃及びｂ₃を算出する（ステップＳ４９）。

（６５）式の左辺におけるベクトルＵ₂以外の部分である（Ｍ₂＋∂²min（ｘ^base（ｋ））Λ）がＭ₃である。また、（６５）式の右辺がｂ₃である。すなわち、ｂ₃＝ｂ₂−（（∂min（ｘ^base（ｋ））−∂²min（ｘ^base（ｋ））ｄ（ｘ^base（ｋ）））νである。ｂ₃を算出するためには、ｄ（ｘ^base（ｋ））を、Ｂｕ（ｋ−１）で算出する。そして処理はステップＳ５７に移行する。

一方、適用範囲０であれば、微分値算出部２２０１は、ｘ^base（ｋ）（ｘ^base _maf（ｋ）であれば（５６）式）を算出し、データ格納部２２０５に格納されているペナルティ関数の１階偏微分関数∂max（ｘ）と、２階偏微分関数∂²max（ｘ）にｘ^base（ｋ）を代入し、∂max（ｘ^base（ｋ））及び∂²max（ｘ^base（ｋ））を求める（ステップＳ５１）。

そして、行列算出部２２０２は、範囲０についてのＭ₃及びｂ₃を算出する（ステップＳ５３）。

（６５）式の左辺におけるベクトルＵ₂以外の部分である（Ｍ₂＋∂²max（ｘ^base（ｋ））Λ）がＭ₃である。また、（６５）式の右辺がｂ₃である。すなわち、ｂ₃＝ｂ₂−（（∂max（ｘ^base（ｋ））−∂²max（ｘ^base（ｋ））ｄ（ｘ^base（ｋ）））νである。ｂ₃を算出するためには、ｄ（ｘ^base（ｋ））を、Ｂｕ（ｋ−１）で算出する。そして処理はステップＳ５７に移行する。

さらに、適用範囲２であれば、微分値算出部２２０１は、Ｍ₃＝Ｍ₂及びｂ₃＝ｂ₂と設定する（ステップＳ５５）。状態制約についてのペナルティ関数は適用されない。その後処理はステップＳ５７に移行する。

そして、出力生成部２２０３は、Ｕ＝Ｍ^-1 ₃ｂ₃によってベクトルＵを算出する（ステップＳ５７）。さらに、出力生成部２２０３は、ベクトルＵから線形化制御入力ｕ（ｋ）（ｉ＝０）を抽出し、出力する（ステップＳ５９）。

その後、モデル予測制御部２２０は、処理が終了するのかを判断し（ステップＳ６１）、処理終了でなければ処理はステップＳ４１に戻り、処理終了であれば処理は終了する。

このような処理を実行すれば、適切なペナルティ関数が適用されるようになり、オーバーシュートやアンダーシュートを効果的に抑制できるようになる。

［実施の形態３］
一般にＭＡＦとＭＡＰの追従性にはトレードオフの関係があるため、ＭＡＦの追従性を上げるためにＭＡＰは悪化するという現象が発生する。そこで、以下のような制御を行うと効果的である。

（Ａ）ＭＡＰがある程度目標値に近い場合は、できるだけＭＡＦを近づけたいので、ＭＡＦのペナルティ関数を加える。
（Ｂ）ＭＡＦはある程度目標値に近いがＭＡＰが目標値から大きく離れている場合、ＭＡＦにペナルティ関数を加えて得られる効果よりＭＡＰの悪化が大きくなることが多いので、ペナルティ関数を加えない。
（Ｃ）ＭＡＦ及びＭＡＰともに目標値から離れている場合には、ＭＡＦ重視であってもＭＡＰもある程度追従していなければＮＯｘ及びＰＭは減らないし、ＭＡＦ及びＭＡＰともに大きく変化させようとする場合、入力のバルブが全開もしくは全閉になることが多く、ペナルティ関数を入れたとしても入力制約の影響で効果が出ないことが多い。従って、計算量の観点もあり、ペナルティ関数を加えない。

そこで、図１２に示すような領域分割を行って、ＭＡＦ及びＭＡＰがいずれの領域に入っているか否かを判断することで、制御を切り替えるものとする。図１２の例では、縦軸はＭＡＰの偏差ｘ_mapを表し、横軸はＭＡＦの偏差ｘ_mafを表す。そして、ＭＡＰの基準値Ａ及びＢを設け、ｘ_maxとＭＡＰの基準値Ａ及びＢで囲まれる範囲を領域０と設定してＭＡＦの状態制約についての上限側のペナルティ関数を適用する。また、ｘ_minとＭＡＰの基準値Ａ及びＢで囲まれる範囲を領域１として設定してＭＡＦの状態制約についての下限側のペナルティ関数を適用する。それ以外の領域については、領域２として設定し、ペナルティ関数を適用しない。

このような領域設定を行う場合には、第２の実施の形態における処理フローのステップＳ４５において、ＭＡＦ及びＭＡＰの状態偏差ｘ（ｋ）に基づき、いずれの領域にｘ（ｋ）が入っているかを判定するようにする。後の処理については図１１と同様である。

このような制御の切り替えを行った結果の効果を、図１３Ａ及び図１３Ｂに示す。図１３Ａは、ＭＡＦの値の時間変化を表しており、図１３Ｂは、ＭＡＰの値の時間変化を表している。ＭＡＰについてはあまり変化していないが、ＭＡＦについては、部分Ｐではオーバーシュートが低減されており、部分Ｑでは目標値への追従性が向上している。これによって、ＮＯｘ及びＰＭ排出量のさらなる低減がなされる。

以上本技術の実施の形態を説明したが、本技術はこれに限定されるものではない。例えば、上で述べた機能ブロック構成については、プログラムモジュール構成とは一致しない場合もある。また、処理結果が変わらない限り、処理フローについても、ステップの順番を入れ替えたり、複数のステップを並列実行するようにしてもよい。

なお、上で述べたオフライン処理装置１００は、コンピュータ装置であって、図１４に示すように、メモリ２５０１とＣＰＵ２５０３とハードディスク・ドライブ（ＨＤＤ）２５０５と表示装置２５０９に接続される表示制御部２５０７とリムーバブル・ディスク２５１１用のドライブ装置２５１３と入力装置２５１５とネットワークに接続するための通信制御部２５１７とがバス２５１９で接続されている。オペレーティング・システム（ＯＳ：Operating System）及び本実施例における処理を実施するためのアプリケーション・プログラムは、ＨＤＤ２５０５に格納されており、ＣＰＵ２５０３により実行される際にはＨＤＤ２５０５からメモリ２５０１に読み出される。ＣＰＵ２５０３は、アプリケーション・プログラムの処理内容に応じて表示制御部２５０７、通信制御部２５１７、ドライブ装置２５１３を制御して、所定の動作を行わせる。また、処理途中のデータについては、主としてメモリ２５０１に格納されるが、ＨＤＤ２５０５に格納されるようにしてもよい。本技術の実施例では、上で述べた処理を実施するためのアプリケーション・プログラムはコンピュータ読み取り可能なリムーバブル・ディスク２５１１に格納されて頒布され、ドライブ装置２５１３からＨＤＤ２５０５にインストールされる。インターネットなどのネットワーク及び通信制御部２５１７を経由して、ＨＤＤ２５０５にインストールされる場合もある。このようなコンピュータ装置は、上で述べたＣＰＵ２５０３、メモリ２５０１などのハードウエアとＯＳ及びアプリケーション・プログラムなどのプログラムとが有機的に協働することにより、上で述べたような各種機能を実現する。

また、上で述べたエンジン制御装置１０００は、コンピュータ装置であって、図１５に示すように、ＲＡＭ（Random Access Memory）４５０１とプロセッサ４５０３とＲＯＭ（Read Only Memory）４５０７とセンサ群４５１５とがバス４５１９で接続されている。本実施の形態における処理を実施するための制御プログラム及び存在している場合にはオペレーティング・システム（ＯＳ：Operating System））は、ＲＯＭ４５０７に格納されており、プロセッサ４５０３により実行される際にはＲＯＭ４５０７からＲＡＭ４５０１に読み出される。プロセッサ４５０３は、センサ群４５１５（例えば、吸気圧センサ及び新気量センサ。場合によっては燃料噴射量測定部及びエンジン回転数測定部など。）を制御して、測定値を取得する。また、処理途中のデータについては、ＲＡＭ４５０１に格納される。なお、プロセッサ４５０３は、ＲＯＭ４５０７を含む場合もあり、さらに、ＲＡＭ４５０１を含む場合もある。本技術の実施の形態では、上で述べた処理を実施するための制御プログラムは、コンピュータ読み取り可能なリムーバブル・ディスクに格納されて頒布され、ＲＯＭライタによってＲＯＭ４５０７に書き込まれる場合もある。このようなコンピュータ装置は、上で述べたプロセッサ４５０３、ＲＡＭ４５０１、ＲＯＭ４５０７などのハードウエアと制御プログラム（場合によってはＯＳも）とが有機的に協働することにより、上で述べたような各種機能を実現する。

以上述べた本実施の形態をまとめると以下のようになる。

本実施の形態の第１の態様に係るオフライン処理方法は、（Ａ）（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ制御対象の状態方程式に関する第１の関数との和を含む第２の関数をラグランジュ乗数について偏微分することで得られ且つ制御対象の状態偏差及び制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）第２の関数を制御対象の状態偏差について偏微分することで得られ且つ制御対象の状態偏差及びラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）第２の関数を制御対象の線形化制御入力について偏微分することで得られ且つラグランジュ乗数及び制御対象の状態偏差を変数として含む第３の一次式群とを格納するデータ格納部から第１の一次式群を読み出して、状態偏差の現在値を除き状態偏差の変数を消し込むように順次代入を行って、第１の一次式群を、状態偏差の現在値と線形化制御入力とを変数として含む第４の一次式群に変形し、（Ｂ）データ格納部から第２の一次式群を読み出して、ラグランジュ乗数の変数を消し込むように順次代入を行って、第２の一次式群を、状態偏差を変数として含む第５の一次式群に変形し、（Ｃ）データ格納部から第３の一次式群を読み出して、第４の一次式群及び第５の一次式群を用いて、第３の一次式群を、状態偏差の現在値と線形化制御入力とを変数として含む第６の一次式群に変形し、（Ｄ）第６の一次式群を、第１の係数行列と線形化制御入力についてのベクトルとの積が、状態偏差の現在値の関数ベクトルと等式で表されるように変形し、（Ｅ）モデル予測制御の制御対象における状態制約に対応するペナルティ関数を状態偏差について１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、（Ｆ）第４の一次式群から線形化制御入力についての係数を抽出し、ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する処理を含む。

このような処理をオフラインで実行しておくことにより、状態制約を導入したモデル予測制御を行うことができるようになる。

上で述べたペナルティ関数が、状態偏差が第１の値以上第２の値以下で０となり、状態偏差が第１の値より小さい場合及び第２の値より大きい場合状態偏差の値に応じて指数関数的に増加する関数である場合もある。制御対象の特性に応じた適切な関数が用いられる。

また、上で述べたペナルティ関数が、状態偏差が第１の値以上第２の値以下で０となり、状態偏差が第１の値より小さい場合に状態偏差の値に応じて指数関数的に増加する第１の関数と、第２の値より大きい場合状態偏差の値に応じて指数関数的に増加する第２の関数とで表される場合もある。この場合、第１及び第２の微分関数を、第１の関数及び第２の関数の各々について生成する場合もある。このようにすれば、制御実行時に第１の関数と第２の関数とを切り替えて用いることができるようになる。

本実施の形態の第２の態様に係る制御方法は、（Ａ）モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、制御期間内における制御対象に対する線形化制御入力のベクトルとの積が関数ベクトルと等しくなるように予め算出されている第１の係数行列と、モデル予測制御の制御対象における状態制約に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、ペナルティ関数を状態偏差について１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から第１の式及び第２の式を読み出して、状態偏差の現在値から得られる基点の値を代入することで、第１の値及び第２の値を算出し、（Ｂ）データ格納部から第１の係数行列及び第２の係数行列を読み出して、第２の値と第２の係数行列との積と第１の係数行列との和である第１の行列を算出し、（Ｃ）データ格納部から関数ベクトルを読み出して、関数ベクトルに状態偏差の現在値を代入することで得られる第１のベクトルと、第１の値と第２の値との差に線形化制御入力の１単位時間前の値から得られる第３の値を乗じた結果と係数ベクトルとの積との差である第２のベクトルを算出し、（Ｄ）第１の行列の逆行列と第２のベクトルとの積により、線形化制御入力のベクトルを生成し、（Ｅ）線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する処理を含む。

このような処理を実行することで、高速に状態制約を伴うモデル予測制御を実行できるようになる。

上で述べた制御方法は、（Ｆ）状態偏差の現在値が所定範囲内であるか判断する処理をさらに含むようにしてもよい。この場合、状態偏差の現在値が所定範囲内であれば、第１の値及び第２の値を算出する処理と、第１の行列を算出する処理と、第２のベクトルを算出する処理と、線形化制御入力のベクトルを生成する処理と、線形化制御入力を抽出する処理とを実行するようにしてもよい。制御対象によっては、状態制約に基づく処理を行うことが好ましくない場合もあるためである。

さらに、上で述べたデータ格納部が、状態偏差の現在値の値域に応じた第１の式及び第２の式を格納している場合もある。この場合、（ａ１）上記第１の値及び第２の値を算出する処理が、状態偏差の現在値が属する値域を特定し、当該特定された値域に応じた第１の式及び第２の式を読み出す処理を含むようにしてもよい。適切な式を状況に応じて切り分けて用いることで、適切な制御を行うことができるようになる。

また、本制御方法が、上で述べた状態偏差の現在値が所定範囲外であれば、（Ｇ）データ格納部から関数ベクトル及び第１の係数行列を読み出して、関数ベクトルに状態偏差の現在値を代入することで得られる第１のベクトルと、第１の係数行列の逆行列との積により、線形化制御入力の第２のベクトルを生成し、（Ｈ）線形化制御入力の第２のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する処理をさらに含むようにしてもよい。状態制約を考慮しないモデル予測制御を高速に実行できる。

なお、上記方法による処理をプロセッサ（又はコンピュータ）に行わせるためのプログラムを作成することができ、当該プログラムは、例えばフレキシブルディスク、ＣＤ−ＲＯＭ、光磁気ディスク、半導体メモリ、ハードディスク等のコンピュータ読み取り可能な記憶媒体又は記憶装置に格納される。尚、中間的な処理結果はメインメモリ等の記憶装置に一時保管される。

以上の実施例を含む実施形態に関し、さらに以下の付記を開示する。

（付記１）
（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ前記制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形し、
前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む第５の一次式群に変形し、
前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を代入することで、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形し、
前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形し、
前記モデル予測制御の制御対象についての制約条件に対応するペナルティ関数を前記状態偏差について１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、
前記第４の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する
処理を、コンピュータに実行させるためのプログラム。

（付記２）
前記ペナルティ関数が、前記状態偏差が第１の値以上第２の値以下で０となり、前記状態偏差が前記第１の値より小さい場合及び前記第２の値より大きい場合前記状態偏差の値に応じて指数関数的に増加する関数である
付記１記載のプログラム。

（付記３）
前記ペナルティ関数が、前記状態偏差が第１の値以上第２の値以下で０となり、前記状態偏差が前記第１の値より小さい場合に前記状態偏差の値に応じて指数関数的に増加する第１の関数と、前記第２の値より大きい場合前記状態偏差の値に応じて指数関数的に増加する第２の関数とで表され、
前記第１及び第２の微分関数を、前記第１の関数及び前記第２の関数の各々について生成する
付記１記載のプログラム。

（付記４）
モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、前記モデル予測制御の制御対象における状態制約に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、前記ペナルティ関数を前記状態偏差について１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び第２の式を読み出して、前記状態偏差の現在値から得られる基点の値を代入することで、第１の値及び第２の値を算出し、
前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出し、
前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルと、前記第１の値と前記第２の値との差に前記線形化制御入力の１単位時間前の値から得られる第３の値を乗じた結果と前記係数ベクトルとの積との和である第２のベクトルを算出し、
前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成し、
前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
処理をプロセッサに実行させるためのプログラム。

（付記５）
前記状態偏差の現在値が所定範囲内であるか判断する処理をさらに含み、
前記状態偏差の現在値が前記所定範囲内であれば、前記第１の値及び第２の値を算出する処理と、前記第１の行列を算出する処理と、前記第２のベクトルを算出する処理と、前記線形化制御入力のベクトルを生成する処理と、前記線形化制御入力を抽出する処理とを実行する
付記４記載のプログラム。

（付記６）
前記データ格納部が、
前記状態偏差の現在値の値域に応じた前記第１の式及び第２の式を格納しており、
前記第１の値及び前記第２の値を算出する処理が、
前記状態偏差の現在値が属する値域を特定し、当該特定された値域に応じた前記第１の式及び前記第２の式を読み出す処理
を含む付記４又は５記載のプログラム。

（付記７）
前記状態偏差の現在値が前記所定範囲外であれば、
前記データ格納部から前記関数ベクトル及び前記第１の係数行列を読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルと、前記第１の係数行列の逆行列との積により、前記線形化制御入力の第２のベクトルを生成し、
前記線形化制御入力の第２のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
処理をさらに含む付記３乃至５のいずれか１つ記載のプログラム。

（付記８）
（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ前記制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形し、
前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む第５の一次式群に変形し、
前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を代入することで、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形し、
前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形し、
前記モデル予測制御の制御対象についての制約条件に対応するペナルティ関数を前記状態偏差について１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、
前記第４の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する
処理を含み、コンピュータにより実行されるオフライン処理方法。

（付記９）
モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、前記モデル予測制御の制御対象における状態制約に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、前記ペナルティ関数を前記状態偏差について１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び第２の式を読み出して、前記状態偏差の現在値から得られる基点の値を代入することで、第１の値及び第２の値を算出し、
前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出し、
前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルと、前記第１の値と前記第２の値との差に前記線形化制御入力の１単位時間前の値から得られる第３の値を乗じた結果と前記係数ベクトルとの積との和である第２のベクトルを算出し、
前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成し、
前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
処理を含み、プロセッサにより実行される制御方法。

（付記１０）
（ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ前記制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形する手段と、
前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む第５の一次式群に変形する手段と、
前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を代入することで、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形する手段と、
前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形する手段と、
前記モデル予測制御の制御対象についての制約条件に対応するペナルティ関数を前記状態偏差について１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成する手段と、
前記第４の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する手段と、
を有するオフライン処理装置。

（付記１１）
モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、前記モデル予測制御の制御対象における状態制約に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、前記ペナルティ関数を前記状態偏差について１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び第２の式を読み出して、前記状態偏差の現在値から得られる基点の値を代入することで、第１の値及び第２の値を算出する第１の算出部と、
前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出し、前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルと、前記第１の値と前記第２の値との差に前記線形化制御入力の１単位時間前の値から得られる第３の値を乗じた結果と前記係数ベクトルとの積との和である第２のベクトルを算出する第２の算出部と、
前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成する生成部と、
前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を出力する出力部と、
を有する制御装置。

１０１ 入力部
１０２ 入力データ格納部
１０３ 第１変形処理部
１０４ 第１データ格納部
１０５ 第２変形処理部
１０６ 第２データ格納部
１０７ 第３変形処理部
１０８ 第３データ格納部
１０９ 第１行列生成部
１１０ 第４データ格納部
１１１ 第２行列生成部
１１２ 微分処理部

Claims (11)

1. （ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形し、
   前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む、前記ラグランジュ乗数についての第５の一次式群に変形し、
   前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を用いて、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形し、
   前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形し、
   前記モデル予測制御の制御対象における状態制約に対応するペナルティ関数を前記状態偏差について１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、
   前記第４の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する
   処理を、コンピュータに実行させるためのプログラム。
2. 前記ペナルティ関数が、前記状態偏差が第１の値以上第２の値以下で０となり、前記状態偏差が前記第１の値より小さい場合及び前記第２の値より大きい場合前記状態偏差の値に応じて指数関数的に増加する関数である
   請求項１記載のプログラム。
3. 前記ペナルティ関数が、前記状態偏差が第１の値以上第２の値以下で０となり、前記状態偏差が前記第１の値より小さい場合に前記状態偏差の値に応じて指数関数的に増加する第１の関数と、前記第２の値より大きい場合前記状態偏差の値に応じて指数関数的に増加する第２の関数とで表され、
   前記第１及び第２の微分関数を、前記第１の関数及び前記第２の関数の各々について生成する
   請求項１記載のプログラム。
4. モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、前記モデル予測制御の制御対象における状態制約に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、前記ペナルティ関数を前記状態偏差について１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び第２の式を読み出して、前記状態偏差の現在値から得られる基点の値を代入することで、第１の値及び第２の値を算出し、
   前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出し、
   前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルと、前記第１の値と前記第２の値との差に前記線形化制御入力の１単位時間前の値から得られる第３の値を乗じた結果と前記係数ベクトルとの積との差である第２のベクトルを算出し、
   前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成し、
   前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
   処理をプロセッサに実行させるためのプログラム。
5. 前記状態偏差の現在値が所定範囲内であるか判断する処理をさらに含み、
   前記状態偏差の現在値が前記所定範囲内であれば、前記第１の値及び第２の値を算出する処理と、前記第１の行列を算出する処理と、前記第２のベクトルを算出する処理と、前記線形化制御入力のベクトルを生成する処理と、前記線形化制御入力を抽出する処理とを実行する
   請求項４記載のプログラム。
6. 前記データ格納部が、
   前記状態偏差の現在値の値域に応じた前記第１の式及び第２の式を格納しており、
   前記第１の値及び前記第２の値を算出する処理が、
   前記状態偏差の現在値が属する値域を特定し、当該特定された値域に応じた前記第１の式及び前記第２の式を読み出す処理
   を含む請求項４又は５記載のプログラム。
7. 前記状態偏差の現在値が前記所定範囲外であれば、
   前記データ格納部から前記関数ベクトル及び前記第１の係数行列を読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルと、前記第１の係数行列の逆行列との積により、前記線形化制御入力の第２のベクトルを生成し、
   前記線形化制御入力の第２のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
   処理をさらに含む請求項４乃至６のいずれか１つ記載のプログラム。
8. （ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形し、
   前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む、前記ラグランジュ乗数についての第５の一次式群に変形し、
   前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を用いて、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形し、
   前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形し、
   前記モデル予測制御の制御対象における状態制約に対応するペナルティ関数を前記状態偏差について１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成し、
   前記第４の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する
   処理を含み、コンピュータにより実行されるオフライン処理方法。
9. モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、前記モデル予測制御の制御対象における状態制約に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、前記ペナルティ関数を前記状態偏差について１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び第２の式を読み出して、前記状態偏差の現在値から得られる基点の値を代入することで、第１の値及び第２の値を算出し、
   前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出し、
   前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルと、前記第１の値と前記第２の値との差に前記線形化制御入力の１単位時間前の値から得られる第３の値を乗じた結果と前記係数ベクトルとの積との差である第２のベクトルを算出し、
   前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成し、
   前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を抽出する
   処理を含み、プロセッサにより実行される制御方法。
10. （ａ）モデル予測制御の評価関数とラグランジュ乗数を含み且つ制御対象の状態方程式に関する第１の関数との和を含む第２の関数を前記ラグランジュ乗数について偏微分することで得られ且つ前記制御対象の状態偏差及び前記制御対象に対する線形化制御入力を変数として含む第１の一次式群と、（ｂ）前記第２の関数を前記制御対象の状態偏差について偏微分することで得られ且つ前記制御対象の状態偏差及び前記ラグランジュ乗数を変数として含む第２の一次式群と、（ｃ）前記第２の関数を前記制御対象の線形化制御入力について偏微分することで得られ且つ前記ラグランジュ乗数及び前記制御対象の状態偏差を変数として含む第３の一次式群とを格納するデータ格納部から前記第１の一次式群を読み出して、前記状態偏差の現在値を除き前記状態偏差の変数を消し込むように順次代入を行って、前記第１の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第４の一次式群に変形する手段と、
    前記データ格納部から前記第２の一次式群を読み出して、前記ラグランジュ乗数の変数を消し込むように順次代入を行って、前記第２の一次式群を、前記状態偏差を変数として含む、前記ラグランジュ乗数についての第５の一次式群に変形する手段と、
    前記データ格納部から前記第３の一次式群を読み出して、前記第４の一次式群及び前記第５の一次式群を用いて、前記第３の一次式群を、前記状態偏差の現在値と前記線形化制御入力とを変数として含む第６の一次式群に変形する手段と、
    前記第６の一次式群を、第１の係数行列と前記線形化制御入力についてのベクトルとの積が、前記状態偏差の現在値の関数ベクトルと等式で表されるように変形する手段と、
    前記モデル予測制御の制御対象における状態制約に対応するペナルティ関数を前記状態偏差について１階偏微分した第１の微分関数及び２階偏微分した第２の微分関数を生成する手段と、
    前記第４の一次式群から前記線形化制御入力についての係数を抽出し、前記ペナルティ関数について予め定められている第２の係数行列及び係数ベクトルを算出する手段と、
    を有するオフライン処理装置。
11. モデル予測制御の制御対象についての状態偏差の現在値についての関数ベクトルと、制御期間内における前記制御対象に対する線形化制御入力のベクトルとの積が前記関数ベクトルと等しくなるように予め算出されている第１の係数行列と、前記モデル予測制御の制御対象における状態制約に対応するペナルティ関数について予め算出されている第２の係数行列及び係数ベクトルと、前記ペナルティ関数を前記状態偏差について１階偏微分した第１の式及び２階偏微分した第２の式とを格納するデータ格納部から前記第１の式及び第２の式を読み出して、前記状態偏差の現在値から得られる基点の値を代入することで、第１の値及び第２の値を算出する第１の算出部と、
    前記データ格納部から前記第１の係数行列及び前記第２の係数行列を読み出して、前記第２の値と前記第２の係数行列との積と前記第１の係数行列との和である第１の行列を算出し、前記データ格納部から前記関数ベクトルを読み出して、前記関数ベクトルに前記状態偏差の現在値を代入することで得られる第１のベクトルと、前記第１の値と前記第２の値との差に前記線形化制御入力の１単位時間前の値から得られる第３の値を乗じた結果と前記係数ベクトルとの積との差である第２のベクトルを算出する第２の算出部と、
    前記第１の行列の逆行列と前記第２のベクトルとの積により、前記線形化制御入力のベクトルを生成する生成部と、
    前記線形化制御入力のベクトルの要素のうち最初の時刻のための線形化制御入力を出力する出力部と、
    を有する制御装置。

Images