JP6260214B2 - 情報処理方法、装置及びプログラム - 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−１）の関数と考えることができる。よって、入力制約について（２０）式乃至（２２）式で記述された最適化問題は、以下のように変形される。

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

まず、∂ｑ_b（τ（ｋ＋ｉ））／∂ｕ_egr（ｋ＋ｊ）は以下のように表される。また、（４３）式から、２行目も得られる。

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

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

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

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

ここで、近似のためペナルティ関数部分の変数として代入するτのベースポイントを設定する。ｘ~（ｋ＋ｉ）（”~”はｘの上に載っているものとする）を、ｕ~に対応する状態偏差ベースポイントと設定する。そして、短い時間における状態偏差の値の変動は小さいので、ｘ^base（ｋ）は以下のように表される。

従って、（４０−２）式から、τのベースポイントτ^base（ｋ）は、以下のように表される。なお、ｕ*（ｋ−１）は、１単位時間前の最適な線形化制御入力である。

このようなベースポイントτ^base（ｋ）を用いると、（４３）式は以下のように変形される。

ここで（５７）式を、以下のように表す。

さらに、τ（ｋ＋ｉ）をτ^base（ｋ）と近似すれば、｛∂²ｑ_b（τ（ｋ＋ｉ））／∂τ（ｋ＋ｉ）²｝（ｕ~）は、全てのｉについてτ^base（ｋ）に関して同じ関数になる。同様に、｛∂ｑ_b（τ（ｋ＋ｉ））／∂τ（ｋ＋ｉ）｝（ｕ~）は、全てのｉについてτ^base（ｋ）に関して同じ関数になる。

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

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

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

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

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

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

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

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

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

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

本実施の形態では、さらに、第１変形処理部１０３は、入力データ格納部１０２に格納されているデータを用いて、制約されるパラメータである空気過剰率τについての変形処理を実行し、処理結果を第１データ格納部１０４に格納する（ステップＳ４）。（３０）式を（４２）式に代入することで、（４３）式に示すようにτ（ｋ＋ｉ）は書き下される。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

なお、エンジンの制御だけではなく、例えばリチウムイオン電池の充放電制御に適用しても良い。リチウムイオン電池の等価回路を図７に示す。リチウムイオン電池の端子電圧Ｖ（ｋ）は、充電率ＳｏＣ（ｋ）と係数ａの積と、抵抗Ｒ１及び容量Ｃ１が並列に接続された部分の電圧ｖ（ｋ）と、直列等価抵抗Ｒ₀と流れる電流ｉ（ｋ）との積との和で表される。すなわち、以下の式で表される。

また、充電率ＳｏＣ（ｋ）は、以下のような関係式に従って変化する。なお、ＦＣは蓄電池の最大蓄積容量を表す。

さらに、電圧ｖ（ｋ）は、以下のような関係式に従って変化する。

一方、状態偏差ではなく、状態についての一般的なモデル予測制御は、以下のように表される。

そうすると、リチウムイオン電池の場合、電圧ｖ（ｋ）と充電率ＳｏＣ（ｋ）を状態ｘ、端子電圧Ｖ（ｋ）を出力ｙ、電流ｉ（ｋ）を入力ｕ（ｋ）として扱えば、以下のように表される。

ここで、ｙ_min＜ｙ（ｋ）＜ｙ_maxという制約条件を適用する場合には、以下のようなペナルティ関数を導入する。

これは、上で述べた（４１）式と同様の形を有する。

また、（３０）式相当の式は、以下のように表されるので、（１０１）式をそれらの式に代入する。そうすると、（４３）式相当の式が得られるようになる。

以上のようにすれば、上で述べたのと同様の式の変形ができる。

［実施の形態２］
上で述べたようなディーゼルエンジンでは、ＭＡＦ及びＭＡＰを変化させる制御を行うが、耐久性の問題から、ＭＡＦ及びＭＡＰの急激な変化は避けることが好ましい。

本実施の形態では、モデル予測制御を用いた、ディーゼルエンジンのＭＡＦ及びＭＡＰ制御において、ＭＡＦ及びＭＡＰの急激な変化を抑制するための制約を加えるものとする。但し、説明を簡単にするためにＭＡＦの急激な変化を抑制する制約のみを考えるものとする。

このような問題は、以下のように表される。

（２０１−５）式に対してペナルティ関数ｑ_cは、以下のように表される。

そうすると、モデル予測制御の目的関数は、以下のように表される。

上でも述べたように、ＭＡＦのみ取り扱うので、（２０１−２）式から、以下の式が得られる。

そうすると、（３０）式を（２０３）式に代入すると、（４３）式相当の式が得られる。

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

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

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

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

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

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

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

ここで、近似のためペナルティ関数部分の変数として代入するｄｘ_mafのベースポイントを設定する。この場合近似は行わずとも、ベースポイントｄｘ^base _maf（ｋ）は以下のように表される。

ディーゼルエンジンの特性上、状態の変化は入力の影響が支配的になるため、初期値の影響が状態の変化に与える影響は極めて小さいものとして取り扱って近似できる。従って、（２０４）式の定数部分は、（２１１）式から以下のように近似される。

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

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

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

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

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

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

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

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

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

なお、本実施の形態に係るオフライン処理装置１００は、図３に示したものと同様の構成を有する。

次に、オフライン処理装置１００の処理内容も、おおよそ第１の実施の形態と同様であるから、図４を用いて説明する。まず、入力部１０１は、ユーザからの入力データとして、行列Ａ、Ｂ、Ｓ及びＲ、予測ホライズンＨ_p及び制御ホライズンＨ_cなどの値の入力を受け付け、入力データ格納部１０２に格納する（図４：ステップＳ１）。なお、入力データ格納部１０２に格納されていない場合には、（２０５）乃至（２０７）式のデータも入力部１０１を介して入力データ格納部１０２に格納されるものとする。同様に、入力制約に係る項の式のデータ（（２４）式）についても同様に、入力データ格納部１０２に格納されており、格納されていない場合には入力部１０１を介して入力される。状態偏差の変化量ｄｘ（ｋ＋ｉ）（＝ｘ（ｋ＋ｉ）−ｘ（ｋ＋ｉ−１））に係る式のデータについても同様に、入力データ格納部１０２に格納されており、格納されていない場合には入力部１０１を介して入力される。

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

本実施の形態では、さらに、第１変形処理部１０３は、入力データ格納部１０２に格納されているデータを用いて、制約されるパラメータである状態偏差の変化量ｄｘについての変形処理を実行し、処理結果を第１データ格納部１０４に格納する（ステップＳ４）。制約条件が異なるので、このステップは第１の実施の形態とは異なる。（３０）式を（２０３）式に代入することで、（２０４）式に示すようにｄｘ（ｋ＋ｉ）は書き下される。

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

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

そして、第１行列生成部１０９は、（３０）式（第１データ格納部１０４に格納されている）及び（３５）式（第３データ格納部１０８に格納されている）を（２１）式（入力データ格納部１０２に格納されている）に代入して連立方程式（３６）式を生成し、これを（３７）式の形に変形することで、行列Ｍ₂及びｂ₂を算出し、第４データ格納部１１０に格納する（ステップＳ９）。

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

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

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

また、エンジン制御装置１０００におけるモデル予測制御部２２０の構成例は、第１の実施の形態と同様である。

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

次に、エンジン制御装置１０００におけるモデル予測制御部２２０の処理内容についても、第１の実施の形態と同様であるから、図６を用いて説明する。

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

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

また、微分値算出部２２０１は、制約されるパラメータのベースポイントの値であるｄｘ^base _maf（ｋ）（（２１３）式）を算出し、データ格納部２２０５に格納されているペナルティ関数の１階偏微分関数∂（ｄｘ）と、２階偏微分関数∂²（ｄｘ）にｄｘ^base _maf（ｋ）を代入し、∂（ｄｘ^base _maf（ｋ））及び∂²（ｄｘ^base _maf（ｋ））を求める（ステップＳ２５）。ここは第１の実施の形態とは異なる。

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

（２１８）式の左辺におけるベクトルＵ₂以外の部分である（Ｍ₂＋∂²（ｄｘ^base _maf（ｋ））Λ^d）がＭ₃である。また、（２１８）式の右辺がｂ₃である。すなわち、ｂ₃＝ｂ₂−（∂（ｄｘ^base _maf（ｋ））−∂²（ｄｘ^base _maf（ｋ））ｄｘ^base _maf（ｋ））ν^dである。ｂ₃を算出するためには、ｄｘ^base _maf（ｋ）を、（２１３）式で算出する。ここは第１の実施の形態とは異なる。

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

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

［実施の形態３］
実際に制御を行う場合、空気過剰率τが制限値を超えてしまう場合や、多少制限値を超えることを許容する場合もある。この場合でも、できるだけ過剰に制限値を超えることは避けることが好ましい。

この場合には、例えば、以下のようなペナルティ関数を採用する。

このペナルティ関数は、図８に示すように変化する関数である。このように、τがτ_min以下となるｑ_bが指数関数的に増加する。

さらに、設計者によっては空気過剰率τが、所望する目標値周辺で推移するような制御を行う場合もある。

この場合も、過度に目標値から離れることを抑制するために、目標値想定推移区間をτ_minからτ_maxとした場合、以下のようなペナルティ関数を採用する場合もある。

このペナルティ関数は、図９に示すように変化する関数である。すなわち、τ_min乃至τ_maxの間ｑ_bは「０」となるが、その区間を離れるほどｑ_bは指数関数的に増加する。

さらに、空気過剰率τをより目標値τ_refに近づけたい場合には、以下のようなペナルティ関数を採用する場合もある。

このようなペナルティ関数を導入したとしても、基本的に第１の実施の形態で説明したオフライン処理は同じである。但し、空気過剰率τの値域によって関数の形が異なるので、それぞれについてτ（ｋ＋ｉ）についての１階偏微分関数及び２階偏微分関数を算出する。すなわち、ステップＳ１３では、各適用範囲について、１階偏微分関数∂（τ）及び２階偏微分関数∂²（τ）を生成する。

そして、τ（ｋ＋ｉ）の値に応じて、ペナルティ関数の適用範囲を特定し、適用範囲に応じた１階偏微分関数∂（τ）及び２階偏微分関数∂²（τ）又は「制約なし」（ペナルティ関数が「０」の場合）を用いることになる。

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

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

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

また、微分値算出部２２０１は、制約されるパラメータの現在値であるτ（ｋ）（＝Ｃｘ（ｋ−１）＋Ｄｕ（ｋ−１））を算出する（ステップＳ４４）。さらに、微分値算出部２２０１は、τ（ｋ）から、ペナルティ関数の適用範囲を特定する（ステップＳ４５）。

そうすると、微分値算出部２２０１は、特定された適用範囲が、ペナルティ関数が０の範囲であるか否かを判断する（ステップＳ４７）。

特定された適用範囲が、ペナルティ関数が０の範囲ではない場合には、微分値算出部２２０１は、制約されるパラメータのベースポイントにおける値であるτ^base（ｋ）（（５６）式）を算出し、データ格納部２２０５に格納されているペナルティ関数の１階偏微分関数∂（τ）と、２階偏微分関数∂²（τ）にτ^base（ｋ）を代入し、∂（τ^base（ｋ））及び∂²（τ^base（ｋ））を求める（ステップＳ５１）。

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

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

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

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

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

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

［実施の形態４］
状態偏差の変化量に関する制約についても、過度の変化を抑える目的もあるので、制御の方法によっては、制約を超えた変化も許す場合もある。この場合も、過度に制約を超えることを抑制するため、以下のようなペナルティ関数を採用する場合もある。

このペナルティ関数は、図１１に示すが、図９に示したものと同様である。

従って、第３の実施の形態で説明したように、オフラインでは、ｄｘの適用範囲毎に、ペナルティ関数の１階偏微分関数∂（ｄｘ）と、２階偏微分関数∂²（ｄｘ）とを算出しておく。さらに、オンライン処理では、ｄｘ（ｋ）（＝ｘ（ｋ）−ｘ（ｋ−１））を算出し、ｄｘ（ｋ）が含まれる適用範囲についての、１階偏微分関数∂（ｄｘ）及び２階偏微分関数∂²（ｄｘ）を用いて、演算を行えばよい。

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

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

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

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

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

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

なお、上で述べたパラメータが、状態偏差と線形化制御入力とから表されるパラメータである場合もある。このような場合には、ディーゼルエンジンの空気過剰率などに対して制約条件を設定することができるようになる。

さらに、上で述べたパラメータが、状態偏差の変化量である場合もある。このような場合には、状態偏差についての変化が急激に発生しないようにすることができるようになる。

さらに、上で述べたペナルティ関数が、パラメータの値域毎に異なる関数で表される場合がある。この場合、上で述べた第１及び第２の微分関数を、パラメータの値域毎に生成するようにしても良い。このようにすれば、制御実行時に第１の関数と第２の関数とを切り替えて用いることができるようになる。

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

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

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

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

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

なお、上で述べたパラメータが、状態偏差と線形化制御入力とから表されるパラメータである場合には、上で述べた第３の値は、制御対象の状態方程式で用いられる行列と線形化制御入力の１単位時間前の値とから算出される場合もある。例えば、空気過剰率の場合には、このような演算が行われる。

さらに、上で述べたパラメータが、状態偏差の変化量である場合には、第３の値は、状態偏差の現在値と状態偏差の１単位時間前の値との差である場合もある。状態偏差の変化量に制約がある場合には、このような演算が行われる。

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

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

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

（付記２）
前記パラメータが、前記状態偏差と前記線形化制御入力とから表されるパラメータである
付記１記載のプログラム。

（付記３）
前記パラメータが、前記状態偏差の変化量である
付記１記載のプログラム。

（付記４）
前記ペナルティ関数が、前記パラメータの値域毎に異なる関数で表され、
前記第１及び第２の微分関数を、前記パラメータの値域毎に生成する
付記１乃至３のいずれか１つ記載のプログラム。

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

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

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

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

（付記９）
前記パラメータが、前記状態偏差と前記線形化制御入力とから表されるパラメータである場合には、
前記第３の値は、前記制御対象の状態方程式で用いられる行列と前記線形化制御入力の１単位時間前の値とから算出される
付記５乃至８のいずれか１つ記載のプログラム。

（付記１０）
前記パラメータが、前記状態偏差の変化量である場合には、
前記第３の値は、前記状態偏差の現在値と前記状態偏差の１単位時間前の値との差である
付記１記載のプログラム。

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

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

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

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

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

Claims (14)

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

Images