(57)【特許請求の範囲】
【請求項１】
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算に関する第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和に関するハミルトニアンの前記入力についての偏微分であって前記第１のラグランジュ乗数と前記第２のラグランジュ乗数と前記解入力との関係を表す第１の条件式と、
（４）前記第１のラグランジュ乗数と前記ハミルトニアンの前記出力状態についての偏微分との関係を表す第２の条件式と、
（５）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
を格納するデータ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第１の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部から前記第２の条件式及び前記関係式を読み出して、前記第２の条件式及び前記関係式を、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解出力状態とで表す第２の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第１の一次式群を前記第２の一次式群に代入することで、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解入力と前記出力状態とで表す第３の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部から前記第１の条件式を読み出して、前記第２のラグランジュ乗数の項を除外した第１の条件式に、前記データ格納部に格納されている前記第３の一次式群を代入することで、前記目標値と前記解入力と前記出力状態との各時刻における関係を表す第４の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第４の一次式群を、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等式で表されるように変形して、前記係数行列と前記関数ベクトルとを前記データ格納部に格納する処理と、
を、コンピュータに実行させるためのプログラム。
【請求項２】
前記制約条件に含まれる前記入力の下限値に関連する前記第２のラグランジュ乗数を前記解入力の第１の一次式で近似した場合における当該第１の一次式の定数項であって前記制約条件に含まれる前記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む定数項を、各時刻について含む第１の制約ベクトルと、
前記第１の一次式における前記解入力の係数であって前記制約条件に含まれる前記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第２の制約ベクトルと、
前記制約条件に含まれる前記入力の上限値に関連する前記第２のラグランジュ乗数を前記解入力の第２の一次式で近似した場合における当該第２の一次式の定数項であって前記制約条件に含まれる前記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む定数項を、各時刻について含む第３の制約ベクトルと、
前記第２の一次式における前記解入力の係数であって前記制約条件に含まれる前記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第４の制約ベクトルと、
を生成し、前記データ格納部に格納する処理、
をさらに前記コンピュータに実行させるための請求項１記載のプログラム。
【請求項３】
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算に関する第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和に関するハミルトニアンの前記出力状態についての偏微分と前記第１のラグランジュ乗数との関係を表す条件式と、
（４）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
（５）前記制約条件に含まれる前記入力の下限値に関連する前記第２のラグランジュ乗数を前記解入力の第１の一次式で近似した場合における当該第１の一次式群と、前記制約条件に含まれる前記入力の上限値に関連する前記第２のラグランジュ乗数を前記解入力の第２の一次式で近似した場合における当該第２の一次式群と、
を格納するデータ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第３の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部から前記条件式及び前記関係式を読み出して、前記条件式及び前記関係式を、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解出力状態とで表す第４の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第３の一次式群を前記第４の一次式群に代入することで、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解入力と前記出力状態とで表す第５の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第３の一次式群に対して、前記データ格納部に格納されている前記第５の一次式群、前記第１の一次式群及び前記第２の一次式群を代入することで、前記目標値と前記解入力と前記出力状態との各時刻における関係を表す第６の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第６の一次式群を、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等式で表されるように変形して、前記係数行列と前記関数ベクトルとを前記データ格納部に格納する処理と、
を、コンピュータに実行させるためのプログラム。
【請求項４】
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算に関する第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和に関するハミルトニアンの前記入力についての偏微分であって前記第１のラグランジュ乗数と前記第２のラグランジュ乗数と前記解入力との関係を表す第１の条件式と、
（４）前記第１のラグランジュ乗数と前記ハミルトニアンの前記出力状態についての偏微分との関係を表す第２の条件式と、
（５）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
を格納するデータ格納部と、
前記データ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第１の一次式群に変換し、前記データ格納部に格納する第１変形処理部と、
前記データ格納部から前記第２の条件式及び前記関係式を読み出して、前記第２の条件式及び前記関係式を、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解出力状態とで表す第２の一次式群に変換し、前記データ格納部に格納する第２変形処理部と、
前記データ格納部に格納されている前記第１の一次式群を前記第２の一次式群に代入することで、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解入力と前記出力状態とで表す第３の一次式群に変換し、前記データ格納部に格納する第３変形処理部と、
前記データ格納部から前記第１の条件式を読み出して、前記第２のラグランジュ乗数の項を除外した第１の条件式に、前記データ格納部に格納されている前記第３の一次式群を代入することで、前記目標値と前記解入力と前記出力状態との各時刻における関係を表す第４の一次式群に変換し、前記データ格納部に格納する式生成部と、
前記データ格納部に格納されている前記第４の一次式群を、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等式で表されるように変形して、前記係数行列と前記関数ベクトルとを前記データ格納部に格納する行列生成部と、
を有する行列生成装置。
【請求項５】
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算に関する第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和に関するハミルトニアンの前記出力状態についての偏微分と前記第１のラグランジュ乗数との関係を表す条件式と、
（４）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
（５）前記制約条件に含まれる前記入力の下限値に関連する前記第２のラグランジュ乗数を前記解入力の第１の一次式で近似した場合における当該第１の一次式群と、前記制約条件に含まれる前記入力の上限値に関連する前記第２のラグランジュ乗数を前記解入力の第２の一次式で近似した場合における当該第２の一次式群と、
を格納するデータ格納部と、
前記データ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第３の一次式群に変換し、前記データ格納部に格納する第１変形処理部と、
前記データ格納部から前記条件式及び前記関係式を読み出して、前記条件式及び前記関係式を、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解出力状態とで表す第４の一次式群に変換し、前記データ格納部に格納する第２変形処理部と、
前記データ格納部に格納されている前記第３の一次式群を前記第４の一次式群に代入することで、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解入力と前記出力状態とで表す第５の一次式群に変換し、前記データ格納部に格納する第３変形処理部と、
前記データ格納部に格納されている前記第３の一次式群に対して、前記データ格納部に格納されている前記第５の一次式群、前記第１の一次式群及び前記第２の一次式群を代入することで、前記目標値と前記解入力と前記出力状態との各時刻における関係を表す第６の一次式群に変換し、前記データ格納部に格納する式生成部と、
前記データ格納部に格納されている前記第６の一次式群を、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等式で表されるように変形して、前記係数行列と前記関数ベクトルとを前記データ格納部に格納する行列生成部と、
を有する行列生成装置。
【請求項６】
制御対象であるプラントの現在の出力状態と現在の目標値とを取得する処理と、
（１）前記プラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算を表す第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和で表されるハミルトニアンの前記入力についての偏微分であって前記第１のラグランジュ乗数と前記第２のラグランジュ乗数と前記解入力との関係を表す第１の条件式と、
（４）前記第１のラグランジュ乗数と前記ハミルトニアンの前記出力状態についての偏微分との関係を表す第２の条件式と、
（５）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
等価な一次方程式群であって、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等しくなるように表されている前記一次方程式群に対して、前記現在の出力状態と前記現在の目標値とを代入して解くことによって、前記制御期間内における解入力を算出する算出処理と、
前記制御期間内における解入力のうち最初の時刻のための解入力を、前記制御対象のプラントに出力する処理と、
を、プロセッサに実行させ、
前記係数行列のサイズが、前記解入力の次元と前記制御期間とに依存する
プログラム。
【請求項７】
前記算出処理が、
前記制御対象のプラントの解入力に対して上限値及び下限値を含む制約条件が設定されている場合に、
前記制約条件に含まれる前記入力の下限値に関連する前記第２のラグランジュ乗数を前記解入力の第１の一次式で近似した場合における当該第１の一次式における前記解入力の係数であって前記制約条件に含まれる前記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第１の制約ベクトルに対して、前記１単位時間前に計算された対応時刻の解入力を設定し、対角行列化することで得られる第１の行列と、
前記制約条件に含まれる前記入力の上限値に関連する前記第２のラグランジュ乗数を前記解入力の第２の一次式で近似した場合における当該第２の一次式における前記解入力の係数であって前記制約条件に含まれる前記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第２の制約ベクトルに対して、前記１単位時間前に計算された対応時刻の解入力を設定し、対角行列化することで得られる第２の行列と、
を生成し、前記係数行列から前記第１の行列を差し引き、前記第２の行列を加算する処理と、
前記制約条件に含まれる前記入力の下限値に関連する前記第２のラグランジュ乗数を前記解入力の第１の一次式で近似した場合における当該第１の一次式における定数項であって前記制約条件に含まれる前記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む定数項に対して前記１単位時間前に計算された対応時刻の解入力を設定して各時刻について並べた第１のベクトルと、
前記制約条件に含まれる前記入力の上限値に関連する前記第２のラグランジュ乗数を前記解入力の第２の一次式で近似した場合における当該第２の一次式における定数項であって前記制約条件に含まれる前記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む定数項に対して前記１単位時間前に計算された対応時刻の解入力を設定して各時刻について並べた第２のベクトルと、
を生成し、前記関数ベクトルに前記第１ベクトルを加算して前記第２ベクトルを差し引く処理と、
を含む請求項６記載のプログラム。
【請求項８】
制御対象であるプラントの現在の出力状態と現在の目標値とを取得する取得部と、
（１）前記プラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算を表す第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和で表されるハミルトニアンの前記入力についての偏微分であって前記第１のラグランジュ乗数と前記第２のラグランジュ乗数と前記解入力との関係を表す第１の条件式と、
（４）前記第１のラグランジュ乗数と前記ハミルトニアンの前記出力状態についての偏微分との関係を表す第２の条件式と、
（５）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
等価な一次方程式群であって、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等しくなるように表されている前記一次方程式群を格納するデータ格納部と、
前記データ格納部から前記一次方程式群を読み出して、前記一次方程式群に対して、前記現在の出力状態と前記現在の目標値とを代入して解くことによって、前記制御期間内における解入力を算出する演算部と、
前記制御期間内における解入力のうち最初の時刻のための解入力を、前記制御対象のプラントに出力する出力部と、
を有し、
前記係数行列のサイズが、前記解入力の次元と前記制御期間とに依存する
プラント制御装置。
【請求項９】
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記解出力状態との誤差の二乗を第１の重み行列で重み付けして初期時刻から１単位時間後の時点から予測区間後の１単位時間前の時点まで加算する第１の項と、前記出力状態の目標値と前記予測区間後における前記解出力状態との誤差の二乗を第２の重み行列で重み付けして加算する第２の項と、前記解入力の二乗を第３の重み行列で重み付けして初期時刻から前記予測区間後の１単位時間前の時点まで加算する第３の項と、制約条件についての第４の項と、を含む評価関数のうち前記第１の項と前記第２の項とを、初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られる、前記解出力状態についての第１の一次式群と、
（４）前記第３の項を前記初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られ且つ前記第３の重み行列の要素値、前記制御区間の値及び前記予測区間の値とで決定される行列要素を有する第１の行列と、
を格納するデータ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第２の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第２の一次式群及び前記第２の一次式群における前記解入力の係数値を、前記データ格納部に格納されている前記第１の一次式群に代入することで第３の一次式群を生成し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第３の一次式群における各前記解入力の項の係数を抽出して第２の行列を生成し、前記第３の一次式群における前記解入力の項以外の項を抽出して、前記制御対象であるプラントの解入力を算出するために用いられるベクトルを生成し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第１の行列と前記第２の行列とを加算して、前記制御対象であるプラントの解入力を算出するために用いられる第３の行列を生成し、前記データ格納部に格納する処理と、
を、コンピュータに実行させるためのプログラム。
【請求項１０】
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記解出力状態との誤差の二乗を第１の重み行列で重み付けして初期時刻から１単位時間後の時点から予測区間後の１単位時間前の時点まで加算する第１の項と、前記出力状態の目標値と前記予測区間後における前記解出力状態との誤差の二乗を第２の重み行列で重み付けして加算する第２の項と、前記解入力の二乗を第３の重み行列で重み付けして初期時刻から前記予測区間後の１単位時間前の時点まで加算する第３の項と、制約条件についての第４の項と、を含む評価関数のうち前記第１の項及び前記第２の項を、係数として乗じられている、前記第１の重み行列及び前記第２の重み行列の要素毎に分類することで得られ且つ当該要素が除外された複数の式の各々を、初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られる、前記解出力状態についての第１の一次式群と、
を格納するデータ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第２の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第２の一次式群及び前記第２の一次式群における前記解入力の係数値を、前記データ格納部に格納されている前記第１の一次式群に代入することで、前記複数の式の各々について第３の一次式群を生成し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第３の一次式群における各前記解入力の項の係数を前記複数の式の各々について抽出して前記複数の式の各々について前記制御対象であるプラントの解入力を算出するために用いられる行列を生成し、前記第３の一次式群における前記解入力の項以外の項を前記複数の式の各々について抽出して前記複数の式の各々について前記制御対象であるプラントの解入力を算出するために用いられるベクトルを生成し、前記データ格納部に格納する処理と、
を、コンピュータに実行させるためのプログラム。
【請求項１１】
前記制約条件に含まれる入力の下限値に関連し且つ前記第１の一次式群に出現するラグランジュ乗数を前記解入力の第１の一次式で近似した場合における当該第１の一次式の定数項であって前記制約条件に含まれる前記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む定数項を、各時刻について含む第１の制約ベクトルと、
前記第１の一次式における前記解入力の係数であって前記制約条件に含まれる前記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第２の制約ベクトルと、
前記制約条件に含まれる前記入力の上限値に関連する前記ラグランジュ乗数を前記解入力の第２の一次式で近似した場合における当該第２の一次式の定数項であって前記制約条件に含まれる前記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む定数項を、各時刻について含む第３の制約ベクトルと、
前記第２の一次式における前記解入力の係数であって前記制約条件に含まれる前記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第４の制約ベクトルと、
を生成し、前記データ格納部に格納する処理、
をさらに前記コンピュータに実行させるための請求項９又は１０記載のプログラム。
【請求項１２】
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記解出力状態との誤差の二乗を第１の重み行列で重み付けして初期時刻から１単位時間後の時点から予測区間後の１単位時間前の時点まで加算する第１の項と、前記出力状態の目標値と前記予測区間後における前記解出力状態との誤差の二乗を第２の重み行列で重み付けして加算する第２の項と、前記解入力の二乗を第３の重み行列で重み付けして初期時刻から前記予測区間後の１単位時間前の時点まで加算する第３の項と、制約条件についての第４の項と、を含む評価関数のうち前記第１の項と前記第２の項とを、初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られる、前記解出力状態についての第１の一次式群と、
（４）前記第３の項を前記初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られ且つ前記第３の重み行列の要素値、前記制御区間の値及び前記予測区間の値とで決定される行列要素を有する第１の行列と、
を格納するデータ格納部と、
前記データ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第２の一次式群に変換し、前記データ格納部に格納する変形部と、
前記データ格納部に格納されている前記第２の一次式群及び前記第２の一次式群における前記解入力の係数値を、前記データ格納部に格納されている前記第１の一次式群に代入することで第３の一次式群を生成し、前記データ格納部に格納し、前記データ格納部に格納されている前記第３の一次式群における各前記解入力の項の係数を抽出して第２の行列を生成し、前記第３の一次式群における前記解入力の項以外の項を抽出して、前記制御対象であるプラントの解入力を算出するために用いられるベクトルを生成し、前記データ格納部に格納し、前記データ格納部に格納されている前記第１の行列と前記第２の行列とを加算して、前記制御対象であるプラントの解入力を算出するために用いられる第３の行列を生成し、前記データ格納部に格納する生成部と、
を有する行列生成装置。
【請求項１３】
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記解出力状態との誤差の二乗を第１の重み行列で重み付けして初期時刻から１単位時間後の時点から予測区間後の１単位時間前の時点まで加算する第１の項と、前記出力状態の目標値と前記予測区間後における前記解出力状態との誤差の二乗を第２の重み行列で重み付けして加算する第２の項と、前記解入力の二乗を第３の重み行列で重み付けして初期時刻から前記予測区間後の１単位時間前の時点まで加算する第３の項と、制約条件についての第４の項と、を含む評価関数のうち前記第１の項及び前記第２の項を、係数として乗じられている、前記第１の重み行列及び前記第２の重み行列の要素毎に分類することで得られ且つ当該要素が除外された複数の式の各々を、初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られる、前記解出力状態についての第１の一次式群と、
を格納するデータ格納部と、
前記データ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第２の一次式群に変換し、前記データ格納部に格納する変形部と、
前記データ格納部に格納されている前記第２の一次式群及び前記第２の一次式群における前記解入力の係数値を、前記データ格納部に格納されている前記第１の一次式群に代入することで、前記複数の式の各々について第３の一次式群を生成し、前記データ格納部に格納し、前記データ格納部に格納されている前記第３の一次式群における各前記解入力の項の係数を前記複数の式の各々について抽出して前記複数の式の各々について前記制御対象であるプラントの解入力を算出するために用いられる行列を生成し、前記第３の一次式群における前記解入力の項以外の項を前記複数の式の各々について抽出して前記複数の式の各々について前記制御対象であるプラントの解入力を算出するために用いられるベクトルを生成し、前記データ格納部に格納する生成部と、
を有する行列生成装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本技術は、プラントの制御技術に関する。
【背景技術】
【０００２】
モデル予測制御とは、各時刻において有限時間未来までの応答を最適化することによって制御対象への入力を決定する制御方法である。非定常目標値への追従制御を実現できる制御技術として注目されおり、飽和などの制約条件も扱えるといった長所もある。
【０００３】
従来は主として線形系を対象にした研究が行われてきたが、近年は非線形系への適用も盛んに研究されている。しかし現在のところ、現場では非線形な問題も線形な問題に簡単化して適用されることが多いが、その線形系に関しても高速に変化するダイナミクスを持つ場合に実適用するのは容易ではなく計算時間の高速化が問題となっている。
【０００４】
非線形なプラントモデルを持つ制御対象も含めて、モデル予測制御についての問題を実時間で実行するための方法としては，大きく分けて２つに分類される。１つはオフラインで状態フィードバック制御則をプラントの状態ｘ（ｋ）の関数として表しておき、オンラインでは代入計算のみを行う方法であり、もう１つはオンラインですべて計算する方法である。
【０００５】
オフラインでの計算は、ハミルトン・ヤコビ・ベルマンの解法に帰着するなどして、現在のプラントの状態ｘ（ｋ）の空間を場合分けし、各領域に対してフィードバック制御則を与える方法などが考えられている。しかし、これらの手法では場合分けの領域の数が非常に多くなることがある。また、状態の次元が大きくなると分割数の急激な増大が起こってしまう。
【０００６】
一方、オンラインで計算する方法は、勾配法や内点法などの最適化問題のさまざまな数値解法が適用できる。また、ラグランジュ乗数に関する双対問題ととらえることによって効率的に解く方法や、ヘッシアンなどの構造を用いて効率化を図ったアルゴリズムも提唱されている。しかし、数値計算の手法として通常の反復計算を用いると、求める解が収束するまで計算を繰り返さなければならず、計算量が大きくなってしまう。
【０００７】
そこで、時刻をパラメータとする連続変形法としてとらえて反復計算を減らす工夫も考えられている。それでもプラントの状態ｘ（ｋ）の変化が大きいときは反復回数が大きくなり、求める時間内に計算が終わらないステップが出てくる可能性がある。
【０００８】
また数値解法では、効率的に解を求めるため、最適問題の式の近似を多く行うことが多く、本来求めたい値から大きくずれる可能性を含んでしまう。
【先行技術文献】
【特許文献】
【０００９】
【特許文献１】ＷＯ９９／５７６１６号公報
【特許文献２】特開２００４−８６９０３号公報
【発明の概要】
【発明が解決しようとする課題】
【００１０】
従って、本技術の目的は、モデル予測制御においてプラントに対する入力を高速に算出するための技術を提供することである。
【課題を解決するための手段】
【００１１】
本技術の第１の態様に係る行列生成方法は、（Ａ）（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、（２）出力状態の初期条件式と、（３）出力状態の目標値と出力状態との誤差の二乗及び入力の二乗の重み付け加算に関する第１の関数と、第１のラグランジュ乗数と状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和に関するハミルトニアンの入力についての偏微分であって第１のラグランジュ乗数と第２のラグランジュ乗数と解入力との関係を表す第１の条件式と、（４）第１のラグランジュ乗数とハミルトニアンの出力状態についての偏微分との関係を表す第２の条件式と、（５）第１のラグランジュ乗数の予測期間後の値と、出力状態の目標値と予測期間後の解出力状態との差に関する関数との関係式と、を格納するデータ格納部から、状態方程式及び初期条件式とを読み出して、状態方程式及び初期条件式を、各時刻における解出力状態を出力状態と解入力とで表す第１の一次式群に変換し、データ格納部に格納する処理と、（Ｂ）データ格納部から第２の条件式及び関係式を読み出して、第２の条件式及び関係式を、各時刻における第１のラグランジュ乗数を目標値と解出力状態とで表す第２の一次式群に変換し、データ格納部に格納する処理と、（Ｃ）データ格納部に格納されている第１の一次式群を第２の一次式群に代入することで、各時刻における第１のラグランジュ乗数を目標値と解入力と出力状態とで表す第３の一次式群に変換し、データ格納部に格納する処理と、（Ｄ）データ格納部から第１の条件式を読み出して、第２のラグランジュ乗数の項を除外した第１の条件式に、データ格納部に格納されている第３の一次式群を代入することで、目標値と解入力と出力状態との各時刻における関係を表す第４の一次式群に変換し、データ格納部に格納する処理と、（Ｅ）データ格納部に格納されている第４の一次式群を、係数行列と、制御期間内における解入力についてのベクトルとの積が、目標値と出力状態とについての関数ベクトルと等式で表されるように変形して、係数行列と関数ベクトルとをデータ格納部に格納する処理とを含む。
【００１２】
本技術の第２の態様に係る行列生成方法は、（Ａ）（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、（２）出力状態の初期条件式と、（３）出力状態の目標値と出力状態との誤差の二乗及び入力の二乗の重み付け加算に関する第１の関数と、第１のラグランジュ乗数と状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和に関するハミルトニアンの前記入力についての偏微分であって第１のラグランジュ乗数と第２のラグランジュ乗数と解入力との関係を表す第１の条件式と、（４）第１のラグランジュ乗数とハミルトニアンの出力状態についての偏微分との関係を表す第２の条件式と、（５）第１のラグランジュ乗数の予測期間後の値と、出力状態の目標値と予測期間後の解出力状態との差に関する関数との関係式と、（６）制約条件に含まれる上記入力の下限値に関連する第２のラグランジュ乗数を解入力の第１の一次式で近似した場合における当該第１の一次式群と、制約条件に含まれる上記入力の上限値に関連する第２のラグランジュ乗数を解入力の第２の一次式で近似した場合における当該第２の一次式群と、を格納するデータ格納部から、状態方程式及び初期条件式とを読み出して、状態方程式及び初期条件式を、各時刻における解出力状態を出力状態と解入力とで表す第３の一次式群に変換し、データ格納部に格納する処理と、（Ｂ）データ格納部から第２の条件式及び関係式を読み出して、第２の条件式及び関係式を、各時刻における第１のラグランジュ乗数を前記目標値と前記解出力状態とで表す第４の一次式群に変換し、データ格納部に格納する処理と、（Ｃ）データ格納部に格納されている第３の一次式群を第４の一次式群に代入することで、各時刻における第１のラグランジュ乗数を目標値と解入力と出力状態とで表す第５の一次式群に変換し、データ格納部に格納する処理と、（Ｄ）データ格納部に格納されている第３の一次式群に対して、データ格納部に格納されている第５の一次式群、第１の一次式群及び第２の一次式群を代入することで、目標値と解入力と出力状態との各時刻における関係を表す第６の一次式群に変換し、データ格納部に格納する処理と、（Ｅ）データ格納部に格納されている第６の一次式群を、係数行列と、制御期間内における解入力についてのベクトルとの積が、目標値と出力状態とについての関数ベクトルと等式で表されるように変形して、係数行列と関数ベクトルとをデータ格納部に格納する処理とを含む。
【００１３】
本技術の第３の態様に係るプラント制御方法は、（Ａ）制御対象であるプラントの現在の出力状態と現在の目標値とを取得する処理と、（Ｂ）（１）プラントの解出力状態と解入力との関係を表す状態方程式と、（２）出力状態の初期条件式と、（３）出力状態の目標値と出力状態との誤差の二乗及び入力の二乗の重み付け加算を表す第１の関数と、第１のラグランジュ乗数と状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和で表されるハミルトニアンの入力についての偏微分であって第１のラグランジュ乗数と第２のラグランジュ乗数と解入力との関係を表す第１の条件式と、（４）第１のラグランジュ乗数とハミルトニアンの出力状態についての偏微分との関係を表す第２の条件式と、（５）第１のラグランジュ乗数の予測期間後の値と、出力状態の目標値と予測期間後の解出力状態との差に関する関数との関係式と、等価な一次方程式群であって、係数行列と、制御期間内における解入力についてのベクトルとの積が、目標値と出力状態とについての関数ベクトルと等しくなるように表されている一次方程式群に対して、現在の出力状態と現在の目標値とを代入して解くことによって、制御期間内における解入力を算出する算出処理と、（Ｃ）制御期間内における解入力のうち最初の時刻のための解入力を、制御対象のプラントに出力する処理とを含む。
【００１４】
本技術の第４の態様に係る行列生成方法は、（Ａ）（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、（２）出力状態の初期条件式と、（３）出力状態の目標値と解出力状態との誤差の二乗を第１の重み行列で重み付けして初期時刻から１単位時間後の時点から予測区間後の１単位時間前の時点まで加算する第１の項と、出力状態の目標値と予測区間後における解出力状態との誤差の二乗を第２の重み行列で重み付けして加算する第２の項と、解入力の二乗を第３の重み行列で重み付けして初期時刻から予測区間後の１単位時間前の時点まで加算する第３の項と、制約条件についての第４の項と、を含む評価関数のうち第１の項と第２の項とを、初期時刻における解入力から制御区間後の１単位時間前の時点における解入力までの各解入力で偏微分することで得られる、解出力状態についての第１の一次式群と、（４）第３の項を初期時刻における解入力から制御区間後の１単位時間前の時点における解入力までの各解入力で偏微分することで得られ且つ第３の重み行列の要素値、制御区間の値及び予測区間の値とで決定される行列要素を有する第１の行列と、を格納するデータ格納部から、状態方程式及び初期条件式とを読み出して、状態方程式及び初期条件式を、各時刻における解出力状態を出力状態と解入力とで表す第２の一次式群に変換し、データ格納部に格納する処理と、（Ｂ）データ格納部に格納されている第２の一次式群及び第２の一次式群における解入力の係数値を、データ格納部に格納されている第１の一次式群に代入することで第３の一次式群を生成し、データ格納部に格納する処理と、（Ｃ）データ格納部に格納されている第３の一次式群における各解入力の項の係数を抽出して第２の行列を生成し、第３の一次式群における解入力の項以外の項を抽出してベクトルを生成し、データ格納部に格納する処理と、（Ｄ）データ格納部に格納されている第１の行列と第２の行列とを加算して、第３の行列を生成し、データ格納部に格納する処理とを含む。
【００１５】
本技術の第５の態様に係る行列生成方法は、（Ａ）（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、（２）出力状態の初期条件式と、（３）出力状態の目標値と解出力状態との誤差の二乗を第１の重み行列で重み付けして初期時刻から１単位時間後の時点から予測区間後の１単位時間前の時点まで加算する第１の項と、出力状態の目標値と予測区間後における解出力状態との誤差の二乗を第２の重み行列で重み付けして加算する第２の項と、解入力の二乗を第３の重み行列で重み付けして初期時刻から予測区間後の１単位時間前の時点まで加算する第３の項と、制約条件についての第４の項と、を含む評価関数のうち第１の項及び第２の項を、係数として乗じられている、第１の重み行列及び第２の重み行列の要素毎に分類することで得られ且つ当該要素が除外された複数の式の各々を、初期時刻における解入力から制御区間後の１単位時間前の時点における解入力までの各解入力で偏微分することで得られる、解出力状態についての第１の一次式群と、を格納するデータ格納部から、状態方程式及び初期条件式とを読み出して、状態方程式及び初期条件式を、各時刻における解出力状態を出力状態と解入力とで表す第２の一次式群に変換し、データ格納部に格納する処理と、（Ｂ）データ格納部に格納されている第２の一次式群及び第２の一次式群における解入力の係数値を、データ格納部に格納されている第１の一次式群に代入することで、複数の式の各々について第３の一次式群を生成し、データ格納部に格納する処理と、（Ｃ）データ格納部に格納されている第３の一次式群における各解入力の項の係数を複数の式の各々について抽出して複数の式の各々について行列を生成し、第３の一次式群における解入力の項以外の項を複数の式の各々について抽出して複数の式の各々についてベクトルを生成し、データ格納部に格納する処理とを含む。
【発明の効果】
【００１６】
モデル予測制御においてプラントに対する入力を高速に算出できるようになる。
【図面の簡単な説明】
【００１７】
【図１】図１は、プラントの一例であるエンジンの模式図である。
【図２】図２は、第１のオフライン行列生成装置の機能ブロック図である。
【図３】図３は、第１のオフライン行列生成処理の処理フローを示す図である。
【図４】図４は、プラントの一例であるエンジンのエンジン制御装置を説明するための図である。
【図５】図５は、第１のエンジン制御装置の簡略化したブロック線図である。
【図６】図６は、第１モデル予測制御部の処理フローを示す図である。
【図７】図７は、第２のオフライン行列生成装置の機能ブロック図である。
【図８】図８は、第２のオフライン行列生成処理の処理フローを示す図である。
【図９】図９は、第２のエンジン制御装置の簡略化したブロック線図である。
【図１０】図１０は、第２モデル予測制御部の処理フローを示す図である。
【図１１】図１１は、オフライン行列生成処理前の連立方程式の状態を示す図である。
【図１２】図１２は、オフライン行列生成処理後の連立方程式の状態を示す図である。
【図１３】図１３は、第３のオフライン行列生成装置の機能ブロック図である。
【図１４】図１４は、第３のオフライン行列生成処理の処理フローを示す図である。
【図１５】図１５は、行列Ｍ１及びベクトルｂの生成について説明するための模式図である。
【図１６】図１６は、第２の実施の形態に係る第１のオフライン行列生成装置の機能ブロック図である。
【図１７】図１７は、第２の実施の形態に係る第１のオフライン行列生成処理の処理フローを示す図である。
【図１８】図１８は、第２の実施の形態に係る第２のオフライン行列生成装置の機能ブロック図である。
【図１９】図１９は、第２の実施の形態に係る第２のオフライン行列生成処理の処理フローを示す図である。
【図２０】図２０は、第３の実施の形態に係る第１のオフライン行列生成装置の機能ブロック図である。
【図２１】図２１は、第３の実施の形態に係る第１のオフライン行列生成処理の処理フローを示す図である。
【図２２】図２２は、第３の実施の形態に係る第２のオフライン行列生成装置の機能ブロック図である。
【図２３】図２３は、第３の実施の形態に係る第２のオフライン行列生成処理の処理フローを示す図である。
【図２４】図２４は、コンピュータの機能ブロック図である。
【図２５】図２５は、エンジン制御装置の機能ブロック図である。
【図２６】図２６は、プラント制御装置の機能ブロック図である。
【発明を実施するための形態】
【００１８】
［実施の形態１］
本実施の形態では、線形システムに対するモデル予測制御（以下、線形モデル予測制御と呼ぶ）を取り扱うことにする。
【００１９】
最初に、この線形モデル予測制御の定式化について説明するが、詳細については、「大塚敏之：非線形Receding Horizon 制御の計算方法について，計測と制御, 41 第5号, 366/371 (2002)」等を参照のこと。
【００２０】
ｍ入力ｎ出力の制御対象は以下の状態方程式で与えられるとする。
【数１】
0005834759.tif 000002
【００２１】
ここで、ｔは時間を表し、ｕ（ｋ）∈Ｒmはプラントへの制御入力ベクトル、ｘ（ｋ）∈Ｒnはプラントの状態ベクトルを表す。また、Ａはｎ×ｎ行列、Ｂはｎ×ｍ行列である。制御入力関数は時間ｔに対して区分的連続とする。
【００２２】
また、入力に関して以下のような制約が加えられているとする。
【数２】
0005834759.tif 000003
【００２３】
このような不等式制約を等式制約として表すために、ダミー変数を用いて以下のような等式制約と表し直す。但し、これらをＣ（ｕ）＝０と表す。
【数３】
0005834759.tif 000004
【００２４】
このような前提の下、まず現在の時間をｔ＝ｋ、時刻ｔ＝ｋ＋Ｔまでを評価区間とし、状態ｘの目標値をＸrefとする。そして以下のような評価関数Ｖを考える。
【数４】
0005834759.tif 000005
【００２５】
ここでＳ，Ｑ及びＲは、それぞれｎ×ｎ、ｎ×ｎ、ｍ×ｍの重み行列である。また、ダミー変数ｕα及びｕβは（３）式のように２次の項に現れるため、符号が正負どちらでも最適解になり得る。そのため、正負どちらが望ましいかを反映させるために、ｐtopbar（上にバーが付されたｐ）及びｐunderbar（下にバーが付されたｐ）を微小な正の値として−ｐTtopbarｕα−ｐTunderbarｕβという項を評価関数の被積分項に付け加えている。この評価関数Ｖを最小にする入力ｕ*を求めることが、線形モデル予測制御の核となる問題である。なお、ｐtopbarは、ｕiの最大値に対応し、ｐunderbarは、ｕiの最小値に対応する。
【００２６】
ラグランジュ乗数をλ及びμとし、ハミルトニアンＨを以下のように定義する。
【数５】
0005834759.tif 000006
【００２７】
評価関数Ｖを最小にする最適制御とその時の軌道が満たすべき必要条件は以下のとおりである。
【数６】
0005834759.tif 000007
【００２８】
従って、Ｖを最小にするｕ*を求めることは、（５）乃至（７）式を満たすｕ*を求めることに帰着する。
【００２９】
ここまでは連続系において定式化を行ったが、実際に計算機で解を求めるためには、評価区間を分割し離散近似した問題を各時刻で解くことになる。通常、最適制御問題の評価区間をプラントの入力周期を刻み幅として分割し、離散近似する。以下、各時刻で解くべき問題を説明する。
【００３０】
入力の周期Δτを分割の幅とし、予測期間を入力周期の倍数Ｔ＝ＨpΔτと設定する。このＨp∈Ｒを予測ホライズンと呼ぶ。以下、簡単化のため、連続系における関数ｆ（ｋ＋ｉ×Δτ）を離散系ではｆ（ｉ）と表すことにする。また、制御入力を変えることができる回数Ｈc∈Ｒを制御ホライズンと呼ぶ。このとき、以下のように表される。
【数７】
0005834759.tif 000008
【００３１】
そして、離散近似された問題に対する最適性の必要条件は、以下の２点境界値問題になる。
【数８】
0005834759.tif 000009
【００３２】
ここでｘ*（ｉ）、ｕ*（ｉ）、λ*（ｉ）及びμ*（ｉ）はそれぞれｉステップ目の最適な出力、入力、ラグランジュ乗数を表す。ΦxTは、ΦTの、ｘについての偏微分を表す。また、Ｈu及びＨxはそれぞれハミルトニアンのｕ及びｘによる偏微分を表す。すなわち、以下のとおりである。
【数９】
0005834759.tif 000010
【００３３】
そして、（８）乃至（１３）式をまとめると、以下の非線形方程式系とみなすことができる。
【数１０】
0005834759.tif 000011
【００３４】
従って、線形モデル予測制御では、各時刻において連立１次方程式（１４）を解くことになる。この（１４）式の計算効率と解の精度がモデル予測制御の解法の速度及び精度を決定づける。
【００３５】
なお、近年のディーゼルエンジンにおいては、エミッションの低減と燃費の向上を目的として、新気量（ＭＡＦ：Mass Air Flow）及び吸気圧（ＭＡＰ：Manifold Air Pressure）が吸気系制御器により最適に制御されている。
【００３６】
一般的に、図１に示すように、ディーゼルエンジンの吸気制御系は、吸気圧制御系と新気量制御系を含み、吸気圧と新気量は、互いに独立に制御されている。吸気圧制御系は、排気中のスス（ＰＭ：Particulate Matter）を低減するために、可変ノズルターボＶＮＴ（Variable Nozzle Turbo）のノズル径を制御して吸気圧を吸気圧目標値に追従するようにコントロールしている。一方、新気量制御系は、排気中の窒素酸化物（ＮＯｘ）を低減するために、排気をシリンダ内に再循環させる排気循環器ＥＧＲ（Exhaust Gas Recirculation）のバルブ開度を制御して新気量を新気量目標値に追従するようにコントロールしている。
【００３７】
このようにディーゼルエンジンでは、２入力２出力制御系となっており、本実施の形態でも、２入力２出力系を検討する。
【００３８】
２入力２出力系の場合には、制御出力ｘ（ｋ）、制御入力ｕ（ｋ）、ダミー変数ｕα及びｕβ、ラグランジュ乗数λ（ｋ）、制御出力の目標値Ｘrefを以下のように表す。
【数１１】
0005834759.tif 000012
なお、ｄtopbar（上にバーが付されたｄ）及びｄunderbar（下にバーが付されたｄ）は変数である。
【００３９】
さらに、制御対象の状態方程式における係数行列を以下のように表す。
【数１２】
0005834759.tif 000013
【００４０】
また、評価関数Ｖの重み行列を以下のように表すものとする。
【数１３】
0005834759.tif 000014
【００４１】
重み行列Ｓ、Ｑ及びＲに関しては一般的には対角行列であるとは限らないが、実際上しばしば用いられる対角行列として説明する。なお、対角行列でない場合においても、特段問題なく以下で述べる処理を行うことができる。
【００４２】
このような前提の下、（１４）式は、以下のように展開される。
【数１４】
0005834759.tif 000015
【数１５】
0005834759.tif 000016
【数１６】
0005834759.tif 000017
【数１７】
0005834759.tif 000018
なお、ｘｒ1及びｘｒ2は、Ｘrefの成分値を表す。
【００４３】
上で述べたようなディーゼルエンジンの吸気系制御のような応答が速い対象に対しては、より高速に連立方程式を解くことが望まれる。さらに、できるだけ誤差がなく正確な解を安定して得るために直接法を採用することが好ましいが、直接法は一般に反復法に比べて変数の増加に対応する計算時間の増加が大きいという問題もある。よって、許容される時間内で計算可能とするためには、連立１次方程式のサイズをできるだけ小さくすることが好ましい。このような背景から、本実施の形態では、オフラインで代数的な簡略手法を用いて連立方程式の変数の数を減らし行列サイズを小さくすることによって、オンラインでの計算量を減らし、直接法の適用も実用的な範囲で可能にする。
【００４４】
上で述べたように、線形モデル予測制御の最適制御問題は、制約のない場合には以下の変数に対して（ｍ＋２ｎ）×Ｈp元連立方程式（３４）となっている。
【数１８】
0005834759.tif 000019
【数１９】
0005834759.tif 000020
【００４５】
［制約条件を考慮しない場合］
制約条件を考慮しない線形モデル予測制御の場合、この連立方程式を解くことで最適な入力を得ることができる。
【００４６】
なお、最終的には、（３６）式を変数とするｍ×Ｈp元連立方程式（３５）を求める。
【数２０】
0005834759.tif 000021
【数２１】
0005834759.tif 000022
【００４７】
本実施の形態では、以下、（３３）式及び（３４）式を以下のように簡略化する処理を実施するオフライン行列生成装置を導入する。
【００４８】
図２に本実施の形態に係るオフライン行列生成装置１００の機能ブロック図を示す。なお、図２では制約なしの場合を説明する。オフライン行列生成装置１００は、入力部１０１と、入力データ格納部１０２と、第１変形処理部１０３と、第１データ格納部１０４と、第２変形処理部１０５と、第２データ格納部１０６と、第３変形処理部１０７と、第３データ格納部１０８と、連立方程式生成部１０９と、第４データ格納部１１０と、行列生成部１１１と、第５データ格納部１１２とを有する。
【００４９】
入力部１０１は、例えばユーザからの入力データである行列Ａ、Ｂ、Ｓ、Ｑ、Ｒ、Ｈp、Ｈc、Δτの入力を受け付け、入力データ格納部１０２に格納する。入力データ格納部１０２は、これらのデータに加えて解くべき連立方程式を表すデータをも格納している。解くべき連立方程式のデータについても入力部１０１から入力されることもある。さらに、入力部１０１は、行列Ａ、Ｂ、Ｓ、Ｑ、Ｒ、Ｈp、Ｈc、Δτのデータを、解くべき連立方程式に設定する処理をも行う。第１変形処理部１０３は、入力データ格納部１０２に格納されたデータを用いて第１変形処理を実施し、処理結果を第１データ格納部１０４に格納する。第２変形処理部１０５は、入力データ格納部１０２に格納されたデータを用いて第２変形処理を実施し、処理結果を第２データ格納部１０６に格納する。第３変形処理部１０７は、第１データ格納部１０４及び第２データ格納部１０６に格納されているデータを用いて第３変形処理を実施し、第３データ格納部１０８に格納する。
【００５０】
連立方程式生成部１０９は、第１データ格納部１０４及び第３データ格納部１０８に格納されているデータを用いて処理を行い、処理結果を第４データ格納部１１０に格納する。行列生成部１１１は、第４データ格納部１１０に格納されているデータを用いて処理を行い、第５データ格納部１１２に格納する。
【００５１】
次に、図３を用いて、オフライン行列生成装置１００の処理内容を説明する。まず、入力部１０１は、ユーザからの入力データである行列Ａ、Ｂ、Ｓ、Ｑ、Ｒ、Ｈp、Ｈc、Δτの入力を受け付け、入力データ格納部１０２に格納する（図３：ステップＳ１）。解くべき連立方程式のデータについても入力部１０１から入力され、入力データ格納部１０２に格納されることもある。さらに、入力部１０１は、行列Ａ、Ｂ、Ｓ、Ｑ、Ｒ、Ｈp、Ｈc、Δτのデータを、解くべき連立方程式に設定する処理をも行って、処理結果を入力データ格納部１０２に格納する。
【００５２】
そして、第１変形処理部１０３は、入力データ格納部１０２に格納されたデータを用いて、最適出力状態ｘ*（ｉ）についての第１変形処理を実施し、処理結果を第１データ格納部１０４に格納する（ステップＳ３）。
【００５３】
具体的には以下のような処理を実施する。ｉ＝０として、（９）式を（８）式に代入する。２入力２出力系では、（１７）式を（１５）式に代入し、（１８）式を（１５）式に代入する。但し、説明を簡略化するために、以下では（８）式及び（９）式で説明する。
【００５４】
このような処理を実施すると、ｘ（ｋ）及びｕ*（０）を変数とする以下の式が得られる。
【数２２】
0005834759.tif 000023
【００５５】
同様に、ｉ＝１とすると、（８）式からｘ*（２）はｘ*（１）とｕ*（１）との一次式で表されるため、上で求めたｘ*（１）の式を代入すると、ｘ（ｋ）、ｕ*（０）及びｕ*（１）を変数とする以下の式が得られる。
【数２３】
0005834759.tif 000024
【００５６】
ｘ*（３）以降も同様に処理を繰り返すことによって、ｘ（ｋ）とｕ*（０），．．，ｕ*（ｉ−１）を変数とする以下の式が得られる。
【数２４】
0005834759.tif 000025
なお、ｉ＝０，１，．．．Ｈp−１である。
【００５７】
このようにｘ*（ｉ）の式を生成する。このような処理結果は第１データ格納部１０４に格納される。
【００５８】
また、第２変形処理部１０５は、入力データ格納部１０２に格納されたデータを用いて、最適ラグランジュ乗数λ*（ｉ）についての第２変形処理を実施し、処理結果を第２データ格納部１０６に格納する（ステップＳ５）。
【００５９】
具体的には以下のような処理を実施する。すなわち、ｉ＝Ｈpとして、（１３）式（２入力２出力系では（３１）式及び（３２）式）をそのまま採用する。但し、ΦxT（ｘ*（Ｈp））は、目標値Ｘref及びｘ*（Ｈp）が変数となっているので、以下のように表される。
【数２５】
0005834759.tif 000026
【００６０】
そして、ｉ＝Ｈp−１として、（１３）式を（１２）式に代入する。２入力２出力系では、（３１）式を（２９）式に代入し、（３２）式を（３０）式に代入する。但し、説明を簡略化するために、以下では（１３）式及び（１２）式で説明する。
【数２６】
0005834759.tif 000027
【００６１】
ここで、ΦxT（ｘ*（Ｈp））は目標値Ｘrefとｘ*（Ｈp）の関数であり、第２項はｘ*（Ｈp−１）とΦxT（ｘ*（Ｈp））の関数であるから、結局λ*（Ｈp−１）は、Ｘrefとｘ*（Ｈp）とｘ*（Ｈp−１）との関数となり、以下のように表される。
【数２７】
0005834759.tif 000028
【００６２】
その後、ｉ＝Ｈp−２として、（３９）式を（１２）式に代入する。このような処理を繰り返すと、λ*（Ｈp−ｊ）はＸrefとｘ*（Ｈp）、ｘ*（Ｈp−１），．．．ｘ*（Ｈp−ｊ）の一次式になる。
【数２８】
0005834759.tif 000029
なお、ｊ＝０，１，．．．，Ｈp−１である。
【００６３】
このようにλ*（Ｈp−ｊ）の式を生成する。このような処理結果は第２データ格納部１０６に格納される。
【００６４】
さらに、第３変形処理部１０７は、第１データ格納部１０４及び第２データ格納部１０６に格納されているデータを用いて、最適ラグランジュ乗数λ*（ｉ）についての第３変形処理を実施し、処理結果を第３データ格納部１０８に格納する（ステップＳ７）。
【００６５】
具体的には以下のような処理を実施する。（４０）式の各λ*（Ｈp−ｊ）に対して各（３７）式を代入する。これによって、各λ*（Ｈp−１）は、以下に表すように、Ｘrefとｕ*（０），．．．，ｕ*（Ｈp−１），ｘ（ｋ）の１次式となる。
【数２９】
0005834759.tif 000030
【００６６】
このようにして得られた各λ*（Ｈp−ｊ）のデータは、第３データ格納部１０８に格納される。
【００６７】
その後、連立方程式生成部１０９は、第１データ格納部１０４及び第３データ格納部１０８に格納されているデータを用いて連立方程式生成処理を実施し、処理結果を第４データ格納部１１０に格納する（ステップＳ９）。
【００６８】
制約がない場合には、μ*＝０とみなすことができ、（４１）式を（１０）式（２入力２出力系では（１９）式及び（２０）式）に代入する。そうすると、（１０）式は、Ｘref、ｕ*（０），．．．，ｕ*（Ｈp−１），ｘ（ｋ）を変数とする一次式の組として以下のように表される。
【数３０】
0005834759.tif 000031
【００６９】
以上から制約条件のない場合における離散化された線形モデル予測制御の最適制御解を求める連立方程式は、ｕ*（０），．．．，ｕ*（Ｈp−１）についての一次式となる。
【００７０】
このような連立方程式のデータは、第４データ格納部１１０に格納される。
【００７１】
そして、行列生成部１１１は、第４データ格納部１１０に格納されたデータを用いて行列Ｍ及びｂを生成して、第５データ格納部１１２に格納する（ステップＳ１１）。
【００７２】
具体的には、以下のような処理を実施する。すなわち、Ｈc＜Ｈpのとき、以下のようになる。
【数３１】
0005834759.tif 000032
【００７３】
従って（４２）式にこの関係式を代入する。そうすると以下のように表される。
【数３２】
0005834759.tif 000033
【００７４】
さらに、ζHc-1＝０，．．．，ζHp-1＝０は同等の式になるので除外する。そうすると、最適制御解を求める連立方程式は、ｕ*（０），．．．，ｕ*（Ｈc−１）を変数としたＨc本の連立方程式と同等となる。
【数３３】
0005834759.tif 000034
【００７５】
この方程式はｍ×Ｈc変数の連立一次方程式となる。ｍは入力の数である。この（４４）式を纏めて、以下の形式に変換する。
【数３４】
0005834759.tif 000035
【００７６】
ここでＭはｍ×Ｈc次の正方行列となり、ｂはＸrefとｘ（ｋ）の１次式を成分として有するｍ×Ｈc次ベクトルである。
【００７７】
この行列Ｍとベクトル関数ｂ（Ｘref，ｘ（ｋ））が、オフライン行列生成装置１００の出力データとなる。
【００７８】
この手法を用いれば、主たる計算部分である連立１次方程式の計算が制御入力の次元と制御ホライズンにのみ依存するため、出力状態の次元を大きくした場合でも制御入力の数が変わらなければ、全体の計算量はほとんど変わらなくなる。
【００７９】
ここで、２入力２出力系における具体例（例１と呼ぶ）について示しておく。制御対象の状態方程式の係数行列を以下のように設定したとする。
【数３５】
0005834759.tif 000036
【００８０】
また、評価関数Ｖの重み行列を設定したとする。
【数３６】
0005834759.tif 000037
【００８１】
さらに、目標値をＸref＝（１，１０）Tとし、Ｈp＝４、Ｈc＝２とすると、制約条件がない場合には、以下のような行列が得られる。
【数３７】
0005834759.tif 000038
【数３８】
0005834759.tif 000039
【数３９】
0005834759.tif 000040
【００８２】
上で述べた処理を行う前の連立１次方程式の行列のサイズは２４×２４であり、行列のサイズが１／６になっている。連立１次方程式を直接法で解く場合、計算時間は行列サイズの３乗に比例するため、連立方程式を解く部分の計算時間は１／２１６となっている。上で述べた処理によって入力の次元と制御ホライズンの設定によっては１回のステップで解くべき連立方程式の行列サイズはきわめて小さくすることができ、それによって反復法ではなく直接法を採用することができ、安定且つ高速に解くことが可能になる。また、上で述べた処理によって得られる連立方程式は制御入力の次元と制御ホライズンＨcのみに依存するため、出力状態ｘ（ｋ）の次元が大きくなったとしても計算量は大幅に増加することはない。例えば、プラントのむだ時間をパディ近似などによって考慮する場合には、その分状態の数が増加する。例えば２次のパディ近似を適用すると２入力６出力系となるが、状態ｘ（ｉ）の要素が６つになっても、上で述べたように最適出力状態ｘ*（ｉ）の項をｕ*（ｉ）の項に置換してゆくので、計算量はｕ*（ｉ）の次数に依存するためほとんど影響がない。
【００８３】
図４に、本技術の実施の形態に係るエンジンの一例としてディーゼルエンジンを示す。エンジン本体１には、エンジン本体１からの排ガスを供給する排気循環器ＥＧＲと、排ガスの圧力にてタービンを回して新気（Fresh Air）を圧縮してエンジン本体１に供給する可変ノズルターボＶＮＴとが接続されている。可変ノズルターボＶＮＴのノズル開度を調整することによって、可変ノズルターボＶＮＴのタービンの回転が調整され、吸気圧（ＭＡＰ）センサで測定される吸気圧（ＭＡＰ）が調整される。一方、排気循環器ＥＧＲに設けられているＥＧＲバルブのバルブ開度を調整することによって、新気量（ＭＡＦ）センサで測定される新気量（ＭＡＦ）が調整される。
【００８４】
本実施の形態に係るエンジン制御装置１０００には、ＭＡＰセンサからの吸気圧測定値と、ＭＡＦセンサからの新気量測定値と、外部から与えられる燃料噴射量の設定値と、同じく外部から与えられるエンジン回転数の設定値とが入力されるようになっている。また、エンジン制御装置１０００からは、ＥＧＲバルブのバルブ開度の操作量（制御入力又は単に入力とも呼ぶ）がＥＧＲバルブに出力され、ＶＮＴノズルのノズル開度の操作量（制御入力又は単に入力とも呼ぶ）がＶＮＴノズルに出力されるようになっている。
【００８５】
本実施の形態に係るエンジン制御装置１０００のブロック線図を図５に示す。エンジン制御装置１０００は、計画器２１０と、第１モデル予測制御部２２０とを有する。プラント２３０は、例えばエンジン特性に相当する。第１モデル予測制御部２２０は、直接法によって連立一次方程式を解く演算部である連立方程式演算部２２１と、データ格納部２２２とを含む。データ格納部２２２には、Ｍ及びｂについてのデータが格納されているものとする。
【００８６】
すなわち、計画器２１０には、燃料噴射量の設定値とエンジン回転数の設定値とが入力されて、燃料噴射量の値及びエンジン回転数の値に対応付けてＥＧＲバルブ開度の目標値及びＶＮＴノズル開度の目標値Ｕref、燃料噴射量の値及びエンジン回転数の値に対応付けて吸気圧ＭＡＰの目標値及び新気量ＭＡＦの目標値Ｘrefが登録されている。本実施の形態では、計画器２１０から、ＥＧＲバルブ開度の目標値及びＶＮＴノズル開度の目標値Ｕref、及び吸気圧ＭＡＰの目標値及び新気量ＭＡＦの目標値Ｘrefを出力する。
【００８７】
この目標値Ｘrefとプラント２３０からの出力ｘ（ｋ）（例えばＭＡＦセンサ出力及びＭＡＰセンサ出力）とが、第１モデル予測制御部２２０に入力される。第１モデル予測制御部２２０は、以下で述べるモデル予測制御の処理を実施して、入力ｕ*（ｉ）をプラント２３０に出力する。このような処理が繰り返される。
【００８８】
次に、第１モデル予測制御部２２０の処理について図６を用いて説明する。第１モデル予測制御部２２０は、目標値Ｘrefと出力ｘ（ｋ）を取得し、ｂに目標値Ｘref及び出力ｘ（ｋ）を代入してｂを数値ベクトル化する（図６：ステップＳ２１）。そして、第１モデル予測制御部２２０の連立方程式演算部２２１は、ＭＵ＝ｂという連立一次方程式を、よく知られている直接法で解き、解Ｕを導出する（ステップＳ２３）。解法自体は従来と同じであるが、用いている行列自体が全く異なっている。上でも述べたように、Ｍもｂも全くサイズが非常に小さくなっているだけではなく、その各成分（要素とも呼ぶ）も異なっている。
【００８９】
すなわち、本実施の形態のオフライン行列生成処理を実施しない場合、（８）式乃至（１３）式をまとめると、その左辺は、プラントの離散状態方程式の係数行列Ａbar及びＢbar（Ａbar＝Ｉ＋ΔτＡ、Ｂbar＝Ｉ＋ΔτＢ）と評価関数Ｖの重み行列Ｓ、Ｑ、Ｒの成分を含む係数行列と、予測期間Ｈp個×入力数ｍの最適入力ｕ*（ｉ）と予測期間Ｈp個×出力数ｎの最適出力状態ｘ*（ｉ）と予測期間Ｈp個×出力数ｎの最適ラグランジュ乗数λ*（ｉ）とを含むベクトルとの積になっている。一方、右辺は、オフライン行列生成処理時には未定のｘ（ｋ）及び目標値Ｘrefについての成分を含むベクトルである。
【００９０】
そしてオフライン行列生成処理を実施することによって、最適出力状態ｘ*（ｉ）及び最適ラグランジュ乗数λ*（ｉ）を最適入力ｕ*（ｉ）の式に変換して代入することで考慮すべき要素を入力ｕ*（ｉ）に限定する。そして、最適入力ｕ*についても最適入力ｕ*（Ｈc）からｕ*（Ｈp−１）をも除外して解くべき要素数をも削減している。このようなオフライン行列生成処理を実施した後であれば、ＭＵ＝ｂにおけるベクトルＵは、制御期間Ｈc個の入力ｕ*（ｉ）のみを成分として含む。また、関数ベクトルｂは、目標値Ｘref及びｘ（ｋ）の関数となっており、制御期間Ｈc個の要素を含む。係数行列Ｍは、ＭＵ＝ｂの形でオフライン行列生成処理前の方程式との等価性を担保するための係数行列である。
【００９１】
従って、連立一次方程式を解く解法については同一でも、全く異なる形の連立一次方程式を解くことになっており、処理負荷が大幅に削減されている。
【００９２】
その後、第１モデル予測制御部２２０は、解Ｕから最適入力値ｕ*を抽出し、プラント２３０に出力する（ステップＳ２５）。最適入力値ｕ*は、解Ｕにおける成分ｕ*（０）である。以上の処理を処理終了になるまで実施する（ステップＳ２７）。
【００９３】
以上のような処理を実施することによって、オンライン処理では高速に最適入力値ｕ*を得ることができるようになる。なお、ｕ*に計画器２１０から出力された、ＥＧＲバルブ開度の目標値及びＶＮＴノズル開度の目標値Ｕrefが重畳されて、プラント２３０に入力される。
【００９４】
［制約条件を考慮する場合］
次に、制約条件を考慮する場合について説明する。制約条件を考慮しない場合には代数的簡略化の手法を用いて行列サイズを小さくしていた。しかしながら、（１１）式の制約条件Ｃ（ｕ）＝０（２入力２出力系であれば（２５）式乃至（２８）式）は線形式ではないため、そのままでは取り扱うことができない。そこで、代数的処理と１次近似を巧く行うことによって、制約条件を、制約条件を考慮しない場合の連立１次方程式に組み込む。
【００９５】
制約条件を考慮する場合のオフライン行列生成処理装置５００の構成を図７に示す。制約条件を考慮しない場合と同じ機能を有するものについては、同じ参照番号が付されている。
【００９６】
オフライン行列生成装置５００は、入力部１０１と、入力データ格納部１０２と、第１変形処理部１０３と、第１データ格納部１０４と、第２変形処理部１０５と、第２データ格納部１０６と、第３変形処理部１０７と、第３データ格納部１０８と、連立方程式生成部１０９と、第４データ格納部１１０と、行列生成部１１１と、第５データ格納部１１２と、第４変形処理部５０１とを有する。このように第４変形処理部５０１が新たに導入されている。
【００９７】
入力部１０１は、例えばユーザからの入力データである行列Ａ、Ｂ、Ｓ、Ｑ、Ｒ、Ｈp、Ｈc、Δτ、ｐtopbar及びｐunderbarの入力を受け付け、入力データ格納部１０２に格納する。入力データ格納部１０２は、これらのデータに加えて解くべき連立方程式を表すデータをも格納している。解くべき連立方程式のデータについても入力部１０１から入力されることもある。さらに、入力部１０１は、行列Ａ、Ｂ、Ｓ、Ｑ、Ｒ、Ｈp、Ｈc、Δτ、ｐtopbar及びｐunderbarのデータを、解くべき連立方程式に設定する処理をも行う。制約条件も入力され、入力データ格納部１０２に格納される場合がある。第１変形処理部１０３は、入力データ格納部１０２に格納されたデータを用いて第１変形処理を実施し、処理結果を第１データ格納部１０４に格納する。第２変形処理部１０５は、入力データ格納部１０２に格納されたデータを用いて第２変形処理を実施し、処理結果を第２データ格納部１０６に格納する。第３変形処理部１０７は、第１データ格納部１０４及び第２データ格納部１０６に格納されているデータを用いて第３変形処理を実施し、第３データ格納部１０８に格納する。
【００９８】
連立方程式生成部１０９は、第１データ格納部１０４及び第３データ格納部１０８に格納されているデータを用いて処理を行い、処理結果を第４データ格納部１１０に格納する。行列生成部１１１は、第４データ格納部１１０に格納されているデータを用いて処理を行い、第５データ格納部１１２に格納する。
【００９９】
第４変形処理部５０１は、入力データ格納部１０２に格納されているデータを用いて処理を行い、処理結果を第５データ格納部１１２に格納する。
【０１００】
次に、図８を用いて、オフライン行列生成装置５００の処理内容を説明する。まず、入力部１０１は、ユーザからの入力データである行列Ａ、Ｂ、Ｓ、Ｑ、Ｒ、Ｈp、Ｈc、Δτ、ｐtopbar及びｐunderbarの入力を受け付け、入力データ格納部１０２に格納する（図８：ステップＳ３１）。制約条件ｕminj及びｕmaxjについても入力する場合がある。解くべき連立方程式のデータについても入力部１０１から入力され、入力データ格納部１０２に格納される場合もある。さらに、入力部１０１は、行列Ａ、Ｂ、Ｓ、Ｑ、Ｒ、Ｈp、Ｈc、Δτ、ｐtopbar及びｐunderbarのデータを、解くべき連立方程式に設定する処理をも行って、処理結果を入力データ格納部１０２に格納する。
【０１０１】
そして、第１変形処理部１０３は、入力データ格納部１０２に格納されたデータを用いて、最適出力状態ｘ*（ｉ）についての第１変形処理を実施し、処理結果を第１データ格納部１０４に格納する（ステップＳ３３）。本ステップは、ステップＳ３と同じであるから説明を省略する。
【０１０２】
また、第２変形処理部１０５は、入力データ格納部１０２に格納されたデータを用いて、最適ラグランジュ乗数λ*（ｉ）についての第２変形処理を実施し、処理結果を第２データ格納部１０６に格納する（ステップＳ３５）。本ステップは、ステップＳ５と同じであるから説明を省略する。
【０１０３】
さらに、第３変形処理部１０７は、第１データ格納部１０４及び第２データ格納部１０６に格納されているデータを用いて、最適ラグランジュ乗数λ*（ｉ）についての第３変形処理を実施し、処理結果を第３データ格納部１０８に格納する（ステップＳ３７）。本ステップは、ステップＳ７と同じであるから、詳細な説明を省略する。
【０１０４】
そして、第４変形処理部５０１は、入力データ格納部１０２及び第３データ格納部１０８に格納されているデータを用いて、制約一次近似演算を実施し、処理結果を第５データ格納部１１２に格納する（ステップＳ３９）。本処理が、制約条件を考慮しない場合と異なる処理である。但し、他の処理への影響はないので、図８の処理の説明の後に詳細に説明する。
【０１０５】
その後、連立方程式生成部１０９は、入力データ格納部１０２に格納されているデータを用いて連立方程式生成処理を実施し、処理結果を第４データ格納部１１０に格納する（ステップＳ４１）。本ステップは、ステップＳ９と同じであるから、詳細な説明を省略する。
【０１０６】
そして、行列生成部１１１は、第４データ格納部１１０に格納されたデータを用いて行列Ｍ及びｂを生成して、第５データ格納部１１２に格納する（ステップＳ４３）。本ステップは、ステップＳ１１と同様であるから、詳細な説明を省略する。これにより、第５データ格納部１１２には、オフライン処理で使用するデータが格納されるようになる。
【０１０７】
次に、制約一次近似演算について、詳細に説明する。ここでは、（１０）式（２入力２出力の場合には（１９）式乃至（２４）式）におけるλ*及びμ*を除去する。すなわち、ｕ*（ｉ）及び１単位時間前のｕ*（ｉ）の値であるｕ*bar（ｉ）（ｕの上にバーが付いたものを、ｕbarと表すものとする。）とで表される１次式を生成する。
【０１０８】
２入力２出力の場合を詳しく説明する。既に、ステップＳ２３において（１５）式乃至（１８）式によりｘ*（ｉ）がｕ*（ｉ）で表されており、その結果と（２９）式乃至（３２）式によりλ*（ｉ）がｕ*で表されているものとする。そしてｉ＝０から、以下の処理を実施する。
【０１０９】
すなわち、制約条件の下限についての（２１）式は、以下のように変形される。
【数４０】
0005834759.tif 000041
【０１１０】
また、（２５）式は、以下のように変形される。
【数４１】
0005834759.tif 000042
【０１１１】
ここで、（４７）式を（４６）式に代入すると、以下のようになる。
【数４２】
0005834759.tif 000043
【０１１２】
これによって、μ*underbar1（０）（下にバーが付されたμ1*（０））はｕ*1（０）の関数の形にすることができる。しかし、一次関数ではないのでこのままでは連立一次方程式として取り扱うことができない。従って、１次近似を用いる。
【０１１３】
１次近似を行う場合には、近似の基点の選び方が重要であるが、本実施の形態では、ｔ＝ｋ−１において最適化計算により得られた同じ時間の入力値ｕ1bar*（１）（すなわちプラント２３０には出力されないが計算される１単位時間後の出力）を基点に選ぶ。そうすると（４８）式は以下のように変形される。
【数４３】
0005834759.tif 000044
【０１１４】
これによって連立一次方程式に組み込むことができるようになる。なお、分かりやすくするために、ｗunderbar及びｖunderbar（underbarは、下線を表す。）を導入する。
【数４４】
0005834759.tif 000045
【０１１５】
従って（４９）式は以下のように表される。
【数４５】
0005834759.tif 000046
【０１１６】
μunderbar2（０）についても、同様に変形を実施した上で１次近似を行う。さらに、ｉを１インクリメントした後に、さらに同様に変形を実施した上で、ｔ＝ｋ−１において最適計算によって得られた（ｉ＋１）単位時間後の計算結果ｕ1bar*（２）及びｕ2bar*（２）をフィードバックして、それを基に１次近似を実施する。このような処理をｉ＝Ｈc−１になるまで繰り返す。そうすると、一般的には以下のように表される。
【数４６】
0005834759.tif 000047
【０１１７】
次に、制約条件の上限についての（２２）式は、ｉ＝０で以下のように変形される。
【数４７】
0005834759.tif 000048
【０１１８】
また、（２６）式は、以下のように変形される。
【数４８】
0005834759.tif 000049
【０１１９】
ここで、（５５）式を（５４）式に代入すると、以下のようになる。
【数４９】
0005834759.tif 000050
【０１２０】
（５６）式を見れば、上限についてのμ*topbar1（０）（上にバーが付されたμ1*（０））についてもｕ*1（０）の関数になることが分かる。従って、同様の１次近似を行う。すなわち、本実施の形態では、ｔ＝ｋ−１において最適化計算により得られた同じ時間の入力値ｕ1bar*（１）（すなわちプラント２３０には出力されないが計算される１単位時間後の出力）を基点に選ぶ。そうすると（５６）式は以下のように変形される。
【数５０】
0005834759.tif 000051
【０１２１】
これによって連立一次方程式に組み込むことができるようになる。なお、分かりやすくするために、ｗtopbar及びｖtopbar（topbarは、上線を表す。）を導入する。
【数５１】
0005834759.tif 000052
【０１２２】
このような処理をｉ＝Ｈc−１になるまで繰り返す。そうすると、一般的には以下のように表される。
【数５２】
0005834759.tif 000053
【０１２３】
μtopbar2（０）についても、同様に変形を実施した上で１次近似を行う。さらに、ｉを１インクリメントした後に、さらに同様に変形を実施した上で、ｔ＝ｋ−１において最適計算によって得られた（ｉ＋１）単位時間後の計算結果ｕ1bar*（２）及びｕ2bar*（２）をフィードバックして、それを基に１次近似を実施する。このような処理をｉ＝Ｈc−１になるまで繰り返す。そうすると、一般的には以下のように表される。
【数５３】
0005834759.tif 000054
【０１２４】
一般的には、上で述べたような変形を行った後、これらと第３変形処理部１０７の処理結果とで連立一次方程式を解くことになるのだが、実際には、全体の連立方程式は、以下のように表されることが分かっている。
【数５４】
0005834759.tif 000055
【数５５】
0005834759.tif 000056
なお、diag（Ｖ）は、ベクトルＶの各成分を対角成分とする対角行列である。
【０１２５】
すなわち、実際に第４変形処理部５０１は、（６３）乃至（６６）式で表されるベクトルを生成することになる。
【数５６】
0005834759.tif 000057
【０１２６】
このようにＭに対しては、上限値と前回計算結果における該当時刻における最適入力値とで表される成分を入力数×制御期間Ｈc分含む第１の制約ベクトルを対角行列に変換した結果を加算し、下限値と前回計算結果における該当時刻における最適入力値とで表される成分を入力数×制御期間Ｈc分含む第２の制約ベクトルを対角行列に変換した結果を差し引くようになる。また、ｂに対しても、上限値と前回計算結果における該当時刻における最適入力値とで表される成分を入力数×制御期間Ｈc分含む第３の制約ベクトルを差し引き、下限値と前回計算結果における該当時刻における最適入力値とで表される成分を入力数×制御期間Ｈc分含む第４の制約ベクトルを加算するようになっている。
【０１２７】
次に、制約条件を考慮する場合のエンジン制御装置１０００のブロック線図を図９に示す。
【０１２８】
エンジン制御装置１０００は、計画器２１０と、第２モデル予測制御部２４０とを有する。プラント２３０は、エンジン特性に相当する。第２モデル予測制御部２４０は、連立方程式演算部２４１と、制約条件反映部２４２と、Ｕ格納部２４３と、データ格納部２４４とを有する。データ格納部２４４は、Ｍ、ｂ、Ｖunderbar、Ｖtopbar、Ｗunderbar及びＷtopbarについてのデータを格納しているものとする。
【０１２９】
すなわち、計画器２１０には、燃料噴射量の設定値とエンジン回転数の設定値とが入力されて、燃料噴射量の値及びエンジン回転数の値に対応付けてＥＧＲバルブ開度の目標値及びＶＮＴノズル開度の目標値Ｕref、燃料噴射量の値及びエンジン回転数の値に対応付けて吸気圧ＭＡＰの目標値及び新気量ＭＡＦの目標値Ｘrefが登録されている。本実施の形態では、計画器２１０から、ＥＧＲバルブ開度の目標値及びＶＮＴノズル開度の目標値Ｕref、及び吸気圧ＭＡＰの目標値及び新気量ＭＡＦの目標値Ｘrefを出力する。
【０１３０】
この目標値Ｘrefとプラント２３０からの出力ｘ（ｋ）（ＭＡＦセンサ出力及びＭＡＰセンサ出力）と、１単位時間前に計算された最適入力値Ｕ*（ｉ−１）と、制約条件とが、第２モデル予測制御部２４０に入力される。制約条件については不変であるならば、オフライン行列生成処理においてＷunderbar、Ｖunderbar、Ｗtopbar及びＶtopbarに埋め込んでしまう場合もある。第２モデル予測制御部２４０は、以下で述べるモデル予測制御の処理を実施して、入力ｕ*（ｉ）をプラント２３０に出力する。
【０１３１】
次に、第２モデル予測制御部２４０の処理について図１０を用いて説明する。第２モデル予測制御部２４０は、目標値Ｘref及びｘ（ｋ）を取得し、ｂに目標値Ｘref及び出力ｘ（ｋ）を代入してｂを数値ベクトル化する（図１０：ステップＳ５１）。そして、第２モデル予測制御部２４０は、制約条件が有効であるか判断する（ステップＳ５３）。固定の設定である場合もあれば、有効無効が外部から設定される場合もある。
【０１３２】
制約条件が有効であれば、制約条件反映部２４２は、Ｕ格納部２４３から１単位時間前の解Ｕ*（ｉ−１）を読み出し、上で述べたようにＶunderbar、Ｖtopbar、Ｗunderbar及びＷtopbarを算出する（ステップＳ５５）。なお、制約条件が、その都度設定される場合には、当該制約条件のデータを取得して、ステップＳ５５で使用する。そして、制約条件反映部２４２は、Ｍ＝Ｍ−diag（Ｖunderbar）＋diag（Ｖtopbar）を算出すると共に（ステップＳ５７）、ｂ＝ｂ＋Ｗunderbar−Ｗtopbarを算出する（ステップＳ５９）。このようにしてＭ及びｂに制約条件が反映される。その後処理はステップＳ６１に移行する。
【０１３３】
制約条件が有効でない場合又はステップＳ５９の後に、連立方程式演算部２４１は、ＭＵ＝ｂという連立一次方程式を、よく知られている直接法で解き、解Ｕを導出し、解ＵをＵ格納部２４３に格納する（ステップＳ６１）。解法自体は従来と同じであるが、用いている行列自体が全く異なっている。
【０１３４】
すなわち、本実施の形態のオフライン行列生成処理を実施しない場合、図１１に示すような行列である。具体的には、左辺は、プラントの離散状態方程式の係数行列Ａbar及びＢbar（Ａbar＝Ｉ＋ΔτＡ。Ｂbar＝Ｉ＋ΔτＢ。）と評価関数Ｖの重み行列Ｓ、Ｑ、Ｒの要素を含む係数行列と、予測期間Ｈp個×入力数ｍの最適入力ｕ*（ｉ）と予測期間Ｈp個×出力数ｎの最適出力状態ｘ*（ｉ）と予測期間Ｈp個×出力数ｎの最適ラグランジュ乗数λ*（ｉ）とを含むベクトルとの積になっている。一方、右辺は、予測期間Ｈp個×入力数ｍ×出力数ｎの最適ラグランジュ乗数μ*（ｉ）と、オフライン行列生成処理時には未定のｘ（ｋ）及び目標値Ｘrefについての成分とを含む。
【０１３５】
そしてオフライン行列生成処理を実施することによって、最適出力状態ｘ*（ｉ）及び最適ラグランジュ乗数λ*（ｉ）を最適入力ｕ*（ｉ）の式に変換して代入することで考慮すべき要素を最適入力ｕ*（ｉ）に限定する。そして、最適入力ｕ*についても入力ｕ*（Ｈc）からｕ*（Ｈp−１）をも除外して解くべき要素数をも削減している。このようなオフライン行列生成処理を実施した後であれば、図１２に示すようなベクトルｂが得られる。すなわち、ＭＵ＝ｂにおけるベクトルＵは、制御期間Ｈc個の最適入力ｕ*（ｉ）のみを成分として含む。また、関数ベクトルｂは、目標値Ｘref及びｘ（ｋ）の関数となっており、制御期間Ｈc個の要素を含む。係数行列Ｍは、ＭＵ＝ｂの形でオフライン行列生成処理前の式との等価性を担保するための係数行列である。
【０１３６】
制約条件がある場合でも、制約条件の値及び１単位時間前の計算結果Ｕ*（ｉ−１）の値が設定されてしまえば定数として扱うことができるので、行列及びベクトルの形は同じである。
【０１３７】
従って、連立一次方程式を解く解法については同一でも、全く異なる形の連立一次方程式を解くことになっており、処理負荷が大幅に削減されている。
【０１３８】
その後、第２モデル予測制御部２４０は、解Ｕから最適入力値ｕ*を抽出し、プラント２３０に出力する（ステップＳ６３）。最適入力値ｕ*は、解Ｕにおけるｕ*（０）である。以上の処理を処理終了になるまで実施する（ステップＳ６５）。
【０１３９】
以上のような処理を実施することによって、オンライン処理では高速に最適入力値ｕ*を得ることができるようになる。このような処理が繰り返される。なお、ｕ*に計画器２１０から出力された、ＥＧＲバルブ開度の目標値及びＶＮＴノズル開度の目標値Ｕrefが重畳されて、プラント２３０に入力される。
【０１４０】
ここで上で述べた例１の設定において、以下のような制約条件が設定された場合を考える。
【数５７】
0005834759.tif 000058
【０１４１】
そうすると、以下のような連立一次方程式が得られる。
【数５８】
0005834759.tif 000059
【数５９】
0005834759.tif 000060
【数６０】
0005834759.tif 000061
【数６１】
0005834759.tif 000062
【数６２】
0005834759.tif 000063
【０１４２】
以上のように、線形モデル予測制御に対して、オフラインで代数的手法を用いて解くべき問題を簡略化し、オンラインにおける計算量を小さくし、そのうえで連立方程式の直接解法を採用する。これにより、速いダイナミクスを持つ機械系に対しても有効なほど高速に計算可能で、且つ安定的に解を求めることができるようになる。さらに、プラントの状態ｘ（ｋ）の次元が仮に大きくなるような場合においても、大幅な計算量の増大を抑えることが可能となる。従って、プラントの実時間制御が可能となる。
【０１４３】
また、１次近似は制約条件を考慮する部分のみに限定しており、精度の低下も最低限に抑えている。
【０１４４】
なお、（５３−１）式及び（５３−２）式並びに（６０）式及び（６１）式が得られれば、これらを（１０）式（２入力２出力系であれば（１９）式及び（２０）式）に代入すれば、μ*（ｉ）の項が除去される。このように変形した後の式に、（４１）式をさらに代入することも可能である。
【０１４５】
具体的には、図１３のようなオフライン行列生成装置６００を導入する。図７との差は、第４変形処理部５０１を除去した点と、連立方程式生成部１０９の代わりに第２連立方程式生成部５０２を導入した点と、第２連立方程式生成部５０２が入力データ格納部１０２に格納されているデータを用いる点とである。また、入力データ格納部１０２には、（５３−１）式及び（５３−２）式並びに（６０）式及び（６１）式のデータも格納されているものとする。
【０１４６】
そして、図１４に示すように、処理フローは、ステップＳ３１乃至Ｓ３７については同じである。ステップＳ３７の後には、第２連立方程式生成部５０２は、上で述べたよう代入処理を実施して連立方程式を生成し、第４データ格納部１１０に格納する（ステップＳ５１）。
【０１４７】
その後、行列生成部１１１は、第４データ格納部１１０に格納されたデータを用いてＭ2Ｕ＝ｂ2における行列Ｍ2及びｂ2を生成して、第５データ格納部１１２に格納する（ステップＳ５３）。本ステップは、ステップＳ１１と基本的には同様であるから、詳細な説明を省略する。但し、（５３−１）式及び（５３−２）式並びに（６０）式及び（６１）式には定数でない項も含まれているので、行列Ｍ及びｂとは中身が異なる。
【０１４８】
［実施の形態２］
本実施の形態では、２入力２出力系について議論するが、同様な方式でｍ入力ｎ出力系に変形できる。
【０１４９】
予測区間Ｈpと制御区間Ｈcとから、第１の実施の形態で述べた評価関数Ｖを離散化すると、以下のように表される。
【数６３】
0005834759.tif 000064
但し、ｕ*（Ｈc）＝ｕ*（Ｈc＋１）＝・・・＝ｕ*（Ｈp−１）
【０１５０】
また、ここで重み行列などを以下のように定義する。
【数６４】
0005834759.tif 000065
【０１５１】
この離散化された評価関数Ｊを最小化するという最適化問題は、以下のような連立方程式を解く問題に帰着する。
【数６５】
0005834759.tif 000066
【０１５２】
一方、この離散化された評価関数Ｊを具体的に書き下すと以下のように表される。
【数６６】
0005834759.tif 000067
【０１５３】
ここで評価関数Ｊを以下のように分けて考える。
【数６７】
0005834759.tif 000068
【０１５４】
すなわち、Ｊ＝Ｊ１＋Ｊ２＋ＣＯＮＳと表すことができ、ＣＯＮＳについては、制約条件についての項であり、制約条件がなければ０として取り扱うことができ、０でなければ第１の実施の形態と同様に処理できる。従って、Ｊ’＝Ｊ１＋Ｊ２について処理することを考える。そうすると、（６８）式乃至（７１）式は、以下の計算を行うことから求めることができる。
【数６８】
0005834759.tif 000069
【０１５５】
ここで、（３７）式で示したように、ｘ1*（ｉ）及びｘ2*（ｉ）は、ｕ1*（０），ｕ2*（０），．．．，ｕ1*（Ｈc−１），ｕ2*（Ｈc−１）の一次式として表される。従って、∂Ｊ１／∂ｕ1*（０）は以下のように表される。
【数６９】
0005834759.tif 000070
【０１５６】
このように書き下せば、ｘ1*（ｉ）及びｘ2*（ｉ）と、ｘ1*（０）におけるｕ1*（０）の係数等とを代入することで、（６８）式の具体的な式を得ることができる。（６８）式乃至（７１）式のうち（６８）式以外についても同様の形に変形でき、同様に式及び係数を代入することで、具体的な式を得ることができる。すなわち、（７２）式のような式を（６８）式乃至（７１）式のそれぞれについて予め用意しておく。
【０１５７】
そして、導入後の式をｕ1*（０），ｕ2*（０），．．．，ｕ1*（Ｈc−１）及びｕ2*（Ｈc−１）の１次式に整理する。そうすると、以下のような形で表される。
【数７０】
0005834759.tif 000071
【０１５８】
なお、ｂ1（Ｘref,x(k)）等は、目標値Ｘref及び初期状態ｘ1（ｋ）及びｘ2（ｋ）の一次式の形にした定数項（ｕ1*（０），ｕ2*（０），．．．，ｕ1*（Ｈc−１）及びｕ2*（Ｈc−１）を変数としてみた場合の定数項）である。
【０１５９】
そうすると、図１５に示すように、ｕ1*（０），ｕ2*（０），．．．，ｕ1*（Ｈc−１）及びｕ2*（Ｈc−１）の係数を抽出して並べることで行列Ｍ１を生成する。さらに、ｂベクトルについては、ｕ1*（０），ｕ2*（０），．．．，ｕ1*（Ｈc−１）及びｕ2*（Ｈc−１）以外の項を、抽出して符号を反転させて生成される。
【０１６０】
また、Ｊ２についてはｕ1*（０），ｕ2*（０），．．．，ｕ1*（Ｈc−１）及びｕ2*（Ｈc−１）の各々について偏微分すると、以下のようになる。
【数７１】
0005834759.tif 000072
【０１６１】
このように重み行列Ｒと予測区間Ｈp及び制御区間Ｈcとが用意されれば、値が得られるようになる。従って、以下のような行列Ｍ２を用意しておき、重み行列Ｒと予測区間Ｈp及び制御区間Ｈcが設定されれば、Ｍ２が得られる。但し、本実施の形態では、Ｈp及びＨcについては予め決定された上で（７２）式等が用意されるので、Ｍ２においてもＨc及びＨpについては設定されているものとする。
【数７２】
0005834759.tif 000073
【０１６２】
このようにして得られたＭ１及びＭ２から、Ｍ＝Ｍ１＋Ｍ２を実行する。この結果として第１の実施の形態の処理結果と同じ解を算出するＭ，ｂを得ることができる。制約については、第１の実施の形態と同様に取り扱えばよい。
【０１６３】
例えば、以下のような行列等が設定されているとする。
【数７３】
0005834759.tif 000074
【０１６４】
そして、（３７）式として以下のような一次式群が得られたものとする。
【数７４】
0005834759.tif 000075
【０１６５】
そうすると、（７２）式等の式に、（８３）式乃至（９０）式及びｕj*（ｉ）の係数を代入すると、以下のような式が得られる。
【数７５】
0005834759.tif 000076
【０１６６】
このような（９１）式乃至（９４）式のｕ1*（０），ｕ2*（０），ｕ1*（１），ｕ2*（１）の係数を抽出すると、以下のような行列Ｍ１が得られる。
【数７６】
0005834759.tif 000077
【０１６７】
さらに、ｕ1*（０），ｕ2*（０），ｕ1*（１），ｕ2*（１）以外の項を抽出して符号を反転させると、以下のようなベクトルｂが得られる。
【数７７】
0005834759.tif 000078
【０１６８】
また、Ｍ２行列に行列Ｒの要素値とＨc及びＨpとを代入すると、以下のような行列が得られる。
【数７８】
0005834759.tif 000079
【０１６９】
そうすると、以下のような行列Ｍが得られる。
【数７９】
0005834759.tif 000080
【０１７０】
このように、演算量は第１の実施の形態よりも削減されており、高速に行列Ｍ及びベクトルｂを生成することができるようになる。
【０１７１】
［制約条件を考慮しない場合］
次に、このような行列Ｍ及びベクトルｂを生成するためのオフライン行列生成装置７００の構成例（制約条件無し）を図１６に示す。オフライン行列生成装置７００は、入力部１０１と、入力データ格納部１０２と、第１変形処理部１０３と、第１データ格納部１０４と、第２行列生成部７０１と、第９データ格納部７０２と、代入処理部７０３と、第６データ格納部７０４と、第１行列生成部７０５と、第７データ格納部７０６と、第３行列生成部７０７と、第８データ格納部７０８とを有する。入力部１０１、入力データ格納部１０２、第１変形処理部１０３及び第１データ格納部１０４については、第１の実施の形態で示したものとほぼ同じである。
【０１７２】
入力部１０１は、例えばユーザからの入力データである行列Ａ、Ｂ、Ｓ、Ｑ、Ｒ、Δτの入力を受け付け、入力データ格納部１０２に格納する。入力データ格納部１０２は、これらのデータに加えて解くべき連立方程式を表すデータと、（７２）式等の評価関数の部分Ｊ１を偏微分した式及び行列Ｍ２のデータをも格納している。解くべき連立方程式のデータ等（（８）式及び（９）式、２入力２出力系では（１５）式乃至（１８）式）についても入力部１０１から入力されることもある。また、第１変形処理部１０３は、入力データ格納部１０２に格納されたデータを用いて第１変形処理を実施し、処理結果を第１データ格納部１０４に格納する。なお、行列Ａ及びＢなどの要素値などを、関係する連立方程式に代入する処理をも実施する。
【０１７３】
第２行列生成部７０１は、入力データ格納部１０２に格納されているデータ（すなわち行列Ｍ２のデータ及び重み行列Ｒ）から、行列Ｍ２を生成して、第９データ格納部７０２に格納する。代入処理部７０３は、入力データ格納部１０２と第１データ格納部１０４とに格納されているデータを用いて、上で述べたような（３７）式とその係数を、（７２）式等に代入する処理を実施して、処理結果を第６データ格納部７０４に格納する。また、行列Ｓ及びＱなどの要素値を、（７２）式等に代入する処理をも実施する。
【０１７４】
第１行列生成部７０５は、第６データ格納部７０４に格納されているデータを用いて、行列Ｍ１及びベクトルｂを生成する処理を実施し、行列Ｍ１について第７データ格納部７０６に格納し、ベクトルｂについては第８データ格納部７０８に格納する。第３行列生成部７０７は、第９データ格納部７０２及び第７データ格納部７０６に格納されている行列Ｍ１と行列Ｍ２を加算して、加算結果である行列Ｍを第８データ格納部７０８に格納する。これによって、第８データ格納部７０８には、行列Ｍ及びベクトルｂが格納されるようになる。
【０１７５】
次に、図１６に示したオフライン行列生成装置７００の処理フローを図１７に示す。まず、入力部１０１は、ユーザからの入力データである行列Ａ、Ｂ、Ｓ、Ｑ、Ｒ、Δτの入力を受け付け、入力データ格納部１０２に格納する（図１７：ステップＳ６１）。解くべき連立方程式のデータについても入力部１０１から入力され、入力データ格納部１０２に格納される場合もある。
【０１７６】
そして、第１変形処理部１０３は、入力データ格納部１０２に格納されたデータを用いて、最適出力状態ｘ*（ｉ）についての第１変形処理を実施し、処理結果である一次式群（すなわち（３７）式）を第１データ格納部１０４に格納する（ステップＳ６３）。本ステップは、ステップＳ３とほぼ同様であるから説明を省略する。但し、行列Ａ及びＢ等のデータについても、例えば本ステップで代入した上で処理する。
【０１７７】
そして、代入処理部７０３は、入力データ格納部１０２から、行列Ｑ及びＳのデータ及び（７２）式等のデータ（すなわち評価関数Ｊの部分Ｊ１についての偏微分の式）を読み出し、第１データ格納部１０４から第１変形処理の処理結果である一次式群のデータを読み出して、行列Ｑ及びＳの要素値等と、第１変形処理の処理結果である一次式群及びその係数を、（７２）式等に代入し、整理する代入処理を実施し、処理結果（（７３）乃至（７５）式）を第６データ格納部７０４に格納する（ステップＳ６５）。
【０１７８】
さらに、第１行列生成部７０５は、第６データ格納部７０４に格納されている代入処理の処理結果から、ｕ1*（０），ｕ2*（０），．．．，ｕ1*（Ｈc−１）及びｕ2*（Ｈc−１）についての係数値を抽出して行列Ｍ１を生成し、第７データ格納部７０６に格納すると共に、定数項を抽出して符号を反転させてベクトルｂを生成し、第８データ格納部７０８に格納する（ステップＳ６７）。
【０１７９】
一方、第２行列生成部７０１は、入力データ格納部１０２から、重み行列Ｒの要素値と行列Ｍ２の式（（８２）式）とを読み出して、行列Ｍ２の式に重み行列Ｒの要素値を代入して、行列Ｍ２を生成し、第９データ格納部７０２に格納する（ステップＳ６９）。
【０１８０】
その後、第３行列生成部７０７は、第７データ格納部７０６に格納されている行列Ｍ１と、第９データ格納部７０２に格納されている行列Ｍ２とを加算することで行列Ｍを生成し、第８データ格納部７０８に格納する（ステップＳ７１）。
【０１８１】
このような処理を実施することで、行列Ｍ及びベクトルｂが生成されて、図５のエンジン制御装置で処理される。
【０１８２】
［制約条件を考慮する場合］
制約条件ありの場合には、図１８に示すようなオフライン行列生成装置８００が用いられる。図１６に示したオフライン行列生成装置７００との差は、第１の実施の形態において制約条件を考慮する場合に導入された第４変形処理部５０１を導入した点である。すなわち、第４変形処理部５０１は、入力データ格納部１０２に格納されているデータから、Ｖtopbar、Ｖunderbar、Ｗtopbar及びＷunderbarを生成し、第８データ格納部７０８に格納する。
【０１８３】
従って制約条件を考慮する場合の処理フローは、図１９に示すように、第４変形処理部５０１が実施する制約１次近似演算を実施するステップであるステップＳ７３が追加された点が、図１７と異なる。ステップＳ７３は、図８のステップＳ３９と同様であり、（６３）乃至（６６）式で表されるベクトルを生成し、第８データ格納部７０８に格納する。
【０１８４】
このような処理を実施することで、行列Ｍ及びベクトルｂと制約条件に関係するＶtopbar、Ｖunderbar、Ｗtopbar及びＷunderbarが生成されて、図９のエンジン制御装置で処理される。
【０１８５】
［実施の形態３］
第１及び第２の実施の形態では、重み行列Ｑ、Ｒ及びＳの要素値を埋め込んだ上で行列Ｍ及びベクトルｂを生成していた。従って、実際のエンジン制御装置の代わりにシミュレーションを行って当該重み行列の要素値の適否を判断する際、当該重み行列の要素値が適切ではないと判断されると、行列Ｍ及びベクトルｂの生成から処理を行うことになる。これでは適切な行列Ｍ及びベクトルｂを生成する際にかかる時間が長くなってしまう。そこで、以下のような工夫を行うことで、繰り返し行列及びベクトルを生成することを回避でき、シミュレーションにおいて重み行列の要素値のみを変更すれば良くなる。
【０１８６】
具体的には、評価関数Ｊの一部分Ｊ’を以下のように分割する。
【数８０】
0005834759.tif 000081
【０１８７】
すなわち、基本的にはＪ１を重み行列の要素毎に項をグループ化する。具体的には、Ｊ３は重み行列Ｑの要素ｑ1が乗じられている項を含み、Ｊ４は、重み行列Ｑの要素ｑ2が乗じられている項を含み、Ｊ５は、重み行列Ｓの要素ｓ1が乗じられている項を含み、Ｊ６は、重み行列Ｓの要素ｓ2が乗じられている項を含む。２入力２出力系以外であっても同様に分割する。なお、Ｊ２については第２の実施の形態と同様である。なお、制約条件についても第２の実施の形態と同様である。
【０１８８】
従って、評価関数Ｊの最小化は、第２の実施の形態と同様に、以下の偏微分を算出することに分解される。
【数８１】
0005834759.tif 000082
【０１８９】
但し、重み行列の要素値については変数扱いにするため、以下のように例えばＪ３であればｑ1を除外した部分をＪ３１として取り扱い、偏微分の式を以下のように用意する。
【数８２】
0005834759.tif 000083
【数８３】
0005834759.tif 000084
【０１９０】
（９５）式のような偏微分の式をｕ1*（０），ｕ2*（０），．．．，ｕ1*（Ｈc−１）及びｕ2*（Ｈc−１）の各々について用意しておき、（３７）式及び（３７）式における係数を代入し、さらに整理する。そして、第２の実施の形態と同様に、ｕ1*（０），ｕ2*（０），．．．，ｕ1*（Ｈc−１）及びｕ2*（Ｈc−１）についての係数を抽出して行列Ｍ３を生成し、定数項を抽出してベクトルｂ３を生成する。
【０１９１】
Ｊ４についても重み行列の要素を除外した部分Ｊ４１についても、（９５）式のような偏微分の式をｕ1*（０），ｕ2*（０），．．．，ｕ1*（Ｈc−１）及びｕ2*（Ｈc−１）の各々について用意しておき、（３７）式及び（３７）式における係数を代入し、さらに整理する。そして、第２の実施の形態と同様に、ｕ1*（０），ｕ2*（０），．．．，ｕ1*（Ｈc−１）及びｕ2*（Ｈc−１）についての係数を抽出して行列Ｍ４を生成し、定数項を抽出してベクトルｂ４を生成する。
【０１９２】
Ｊ５及びＪ６についても同様にして、行列Ｍ５及びベクトルｂ５、行列Ｍ６及びベクトルｂ６を生成する。
【０１９３】
最終的に、生成すべき行列Ｍ及びベクトルｂは、以下のように表される。
Ｍ＝ｑ1Ｍ３＋ｑ2Ｍ４＋ｓ1Ｍ５＋ｓ2Ｍ６＋Ｍ２  （９３）
ｂ＝ｑ1ｂ３＋ｑ2ｂ４＋ｓ1ｂ５＋ｓ2ｂ６        （９４）
【０１９４】
但し、このように加算せず、個々の行列及びベクトルを、重み行列の要素と対応付けて保存しておけばよい。そうすれば、シミュレーション実施前に、行列Ｓ、Ｑ及びＲを設定して、行列Ｍ及びベクトルｂを（９３）式及び（９４）式で生成してから、シミュレーションを実施すればよい。
【０１９５】
例えば、行列Ａ及びＢが以下のように設定され、Ｈp＝４、Ｈc＝２と設定され、行列Ｑ，Ｒ及びＳが以下のような２行２列の対角行列であるものとする。
【数８４】
0005834759.tif 000085
【数８５】
0005834759.tif 000086
【０１９６】
そして（３７）式が、（８３）式乃至（９０）式であったとすると、Ｊ３１の偏微分の結果は以下のようになる。
【数８６】
0005834759.tif 000087
【０１９７】
従って、行列Ｍ３及びベクトルｂ３については以下のように表される。
【数８７】
0005834759.tif 000088
【数８８】
0005834759.tif 000089
【０１９８】
また、Ｊ４１の偏微分の結果は以下のようになる。
【数８９】
0005834759.tif 000090
【０１９９】
従って、行列Ｍ４及びベクトルｂ４については以下のように表される。
【数９０】
0005834759.tif 000091
【数９１】
0005834759.tif 000092
【０２００】
さらに、Ｊ５１の偏微分の結果は以下のようになる。
【数９２】
0005834759.tif 000093
【０２０１】
従って、行列Ｍ５及びベクトルｂ５については以下のように表される。
【数９３】
0005834759.tif 000094
【数９４】
0005834759.tif 000095
【０２０２】
また、Ｊ６１の偏微分の結果は以下のようになる。
【数９５】
0005834759.tif 000096
【０２０３】
従って、行列Ｍ６及びベクトルｂ６については以下のように表される。
【数９６】
0005834759.tif 000097
【数９７】
0005834759.tif 000098
【０２０４】
［制約条件を考慮しない場合］
次に、このような行列Ｍ２乃至Ｍ６及びベクトルｂ３乃至ｂ６を生成するためのオフライン行列生成装置９００の構成例（制約条件無し）を図２０に示す。オフライン行列生成装置９００は、入力部１０１と、入力データ格納部１０２と、第１変形処理部１０３と、第１データ格納部１０４と、第２代入処理部９０１と、第１１データ格納部９０２と、第４行列生成部９０３と、第１２データ格納部９０４と、第５行列生成部９０５と、第１３データ格納部９０６とを有する。入力部１０１、入力データ格納部１０２、第１変形処理部１０３及び第１データ格納部１０４については、第１の実施の形態で示したものとほぼ同じである。
【０２０５】
入力部１０１は、例えばユーザからの入力データである行列Ａ、Ｂ、Δτの入力を受け付け、入力データ格納部１０２に格納する。入力データ格納部１０２は、これらのデータに加えて解くべき連立方程式を表すデータ（（８）式及び（９）式、２入力２出力系であれば（１５）乃至（１８）式）と、（９５）式等の、評価関数の部分Ｊ３１乃至Ｊ６１を偏微分した式及び行列Ｍ２のデータをも格納している。行列Ｍ２のデータについては、重み行列Ｒが変数扱いで、Ｈp及びＨcが確定しているので、行列Ｍ２のデータについては処理せずとも良い。また、解くべき連立方程式のデータ等についても入力部１０１から入力されることもある。さらに、第１変形処理部１０３は、入力データ格納部１０２に格納されたデータを用いて第１変形処理を実施し、処理結果を第１データ格納部１０４に格納する。なお、行列Ａ及びＢなどの要素値などを、関係する連立方程式に代入する処理をも実施する。
【０２０６】
第２代入処理部９０１は、入力データ格納部１０２と第１データ格納部１０４とに格納されているデータを用いて、上で述べたような（３７）式とその係数を、（９５）式等に代入する処理を実施して、処理結果を第１１データ格納部９０２に格納する。
【０２０７】
第４行列生成部９０３は、第１１データ格納部９０２に格納されているデータを用いて、行列Ｍ３、Ｍ４、Ｍ５及びＭ６と、ベクトルｂ３、ｂ４、ｂ５及びｂ６を生成する処理を実施し、第１２データ格納部９０４に格納する。第５行列生成部９０５は、入力データ格納部１０２及び第１２データ格納部９０４に格納されているデータを用いて、（９３）式及び（９４）式の演算を実施して、演算結果を第１３データ格納部９０６に格納する。このようにすれば、重み行列Ｑ，Ｓ及びＲの要素が変数化された状態の行列Ｍ及びベクトルｂが、第１３データ格納部９０６に格納されており、シミュレーションなどに用いることができるようになる。
【０２０８】
次に、図２０に示したオフライン行列生成装置９００の処理フローを図２１に示す。まず、入力部１０１は、ユーザからの入力データである行列Ａ、Ｂ、Δτの入力を受け付け、入力データ格納部１０２に格納する（図２１：ステップＳ８１）。解くべき連立方程式のデータについても入力部１０１から入力され、入力データ格納部１０２に格納される場合もある。
【０２０９】
そして、第１変形処理部１０３は、入力データ格納部１０２に格納されたデータを用いて、最適出力状態ｘ*（ｉ）についての第１変形処理を実施し、処理結果である一次式群（すなわち（３７）式）を第１データ格納部１０４に格納する（ステップＳ８３）。本ステップは、ステップＳ３とほぼ同様であるから説明を省略する。但し、行列Ａ及びＢ等のデータについても、例えば本ステップで代入した上で処理する。これより前に代入しても良い。
【０２１０】
そして、第２代入処理部９０１は、入力データ格納部１０２から、（９２）式等のデータ（すなわち評価関数Ｊの部分Ｊ３１乃至Ｊ６１についての偏微分の式）を読み出し、第１データ格納部１０４から第１変形処理の処理結果である一次式群のデータを読み出して、第１変形処理の処理結果である一次式群及びその係数を、（９２）式等に代入し、整理する代入処理を実施し、処理結果を第１１データ格納部９０２に格納する（ステップＳ８５）。
【０２１１】
さらに、第４行列生成部９０３は、第１１データ格納部９０２に格納されている代入処理の処理結果から、ｕ1*（０），ｕ2*（０），．．．，ｕ1*（Ｈc−１）及びｕ2*（Ｈc−１）についての係数値を抽出して行列Ｍ３乃至Ｍ６を生成し、さらに定数項を抽出して符号を反転させてベクトルｂ３乃至ｂ６を生成し、第１２データ格納部９０４に格納する（ステップＳ８７）。
【０２１２】
その後、第５行列生成部９０５は、第１２データ格納部９０４に格納されている行列Ｍ３乃至Ｍ６に対して、対応する重み行列Ｑ及びＳの要素を乗じると共に、入力データ格納部１０２に格納されている行列Ｍ２を加算して行列Ｍを生成し、第１３データ格納部９０６に格納する（ステップＳ８９）。すなわち、ｑ1Ｍ３＋ｑ2Ｍ４＋ｓ1Ｍ５＋ｓ2Ｍ６＋Ｍ２により行列Ｍを生成する。
【０２１３】
さらに、第５行列生成部９０５は、第１２データ格納部９０４に格納されているベクトルｂ３乃至ｂ６に対して、対応する重み行列Ｑ及びＳの要素を乗じることでベクトルｂを生成し、第１３データ格納部９０６に格納する（ステップＳ９１）。すなわち、ｑ1ｂ３＋ｑ2ｂ４＋ｓ1ｂ５＋ｓ2ｂ６によりベクトルｂを生成する。
【０２１４】
このような処理を実施することで、重み行列Ｑ、Ｓ及びＲの要素が変数化された行列Ｍ及びベクトルｂが生成されて、シミュレーションにおいて重み行列の調整を実施することが容易になる。
【０２１５】
なお、ステップＳ８９及びＳ９１については、シミュレーションにおいて実施するようにしても良い。
【０２１６】
［制約条件を考慮する場合］
制約条件ありの場合には、図２２に示すようなオフライン行列生成装置１１００が用いられる。図２０に示したオフライン行列生成装置９００との差は、第１の実施の形態において制約条件を考慮する場合に導入された第４変形処理部５０１を導入した点である。すなわち、第４変形処理部５０１は、入力データ格納部１０２に格納されているデータから、Ｖtopbar、Ｖunderbar、Ｗtopbar及びＷunderbarを生成し、第１３データ格納部９０６に格納する。
【０２１７】
従って制約条件を考慮する場合の処理フローは、図２３に示すように、第４変形処理部５０１が実施する制約１次近似演算を実施するステップであるステップＳ９３が追加された点が、図２１と異なる。ステップＳ９３は、図８のステップＳ３９と同様であり、（６３）乃至（６６）式で表されるベクトルを生成し、第１３データ格納部９０６に格納する。
【０２１８】
このような処理を実施することで、重み行列Ｑ、Ｒ及びＳの要素が変数化された行列Ｍ及びベクトルｂと、制約条件に関係するＶtopbar、Ｖunderbar、Ｗtopbar及びＷunderbarとが生成されて、シミュレーションにおいて用いられる。
【０２１９】
以上本技術の実施の形態を説明したが、本技術はこれに限定されるものではない。例えば、機能ブロック図及びブロック線図は一例であって、必ずしも実際のプログラムモジュール構成と一致しない場合もある。また、処理フローについても、処理結果が変わらない限り、並列実施したり、処理順番を入れ替えるなどの変形が可能である。さらに、上ではエンジン制御装置への適用を述べたが、他のプラントの制御装置への適用も可能である。
【０２２０】
なお、上で述べたオフライン行列生成装置１００、５００乃至１１００は、コンピュータ装置であって、図２４に示すように、メモリ２５０１とＣＰＵ２５０３とハードディスク・ドライブ（ＨＤＤ）２５０５と表示装置２５０９に接続される表示制御部２５０７とリムーバブル・ディスク２５１１用のドライブ装置２５１３と入力装置２５１５とネットワークに接続するための通信制御部２５１７とがバス２５１９で接続されている。オペレーティング・システム（ＯＳ：Operating System）及び本実施例における処理を実施するためのアプリケーション・プログラムは、ＨＤＤ２５０５に格納されており、ＣＰＵ２５０３により実行される際にはＨＤＤ２５０５からメモリ２５０１に読み出される。ＣＰＵ２５０３は、アプリケーション・プログラムの処理内容に応じて表示制御部２５０７、通信制御部２５１７、ドライブ装置２５１３を制御して、所定の動作を行わせる。また、処理途中のデータについては、主としてメモリ２５０１に格納されるが、ＨＤＤ２５０５に格納されるようにしてもよい。本技術の実施例では、上で述べた処理を実施するためのアプリケーション・プログラムはコンピュータ読み取り可能なリムーバブル・ディスク２５１１に格納されて頒布され、ドライブ装置２５１３からＨＤＤ２５０５にインストールされる。インターネットなどのネットワーク及び通信制御部２５１７を経由して、ＨＤＤ２５０５にインストールされる場合もある。このようなコンピュータ装置は、上で述べたＣＰＵ２５０３、メモリ２５０１などのハードウエアとＯＳ及びアプリケーション・プログラムなどのプログラムとが有機的に協働することにより、上で述べたような各種機能を実現する。
【０２２１】
また、上で述べたエンジン制御装置１０００は、コンピュータ装置であって、図２５に示すように、ＲＡＭ（Random Access Memory）４５０１とプロセッサ４５０３とＲＯＭ（Read Only Memory）４５０７とセンサ群４５１５とがバス４５１９で接続されている。本実施の形態における処理を実施するための制御プログラム及び存在している場合にはオペレーティング・システム（ＯＳ：Operating System））は、ＲＯＭ４５０７に格納されており、プロセッサ４５０３により実行される際にはＲＯＭ４５０７からＲＡＭ４５０１に読み出される。プロセッサ４５０３は、センサ群４５１５（例えば、吸気圧センサ及び新気量センサ。場合によっては燃料噴射量測定部及びエンジン回転数測定部など。）を制御して、測定値を取得する。また、処理途中のデータについては、ＲＡＭ４５０１に格納される。なお、プロセッサ４５０３は、ＲＯＭ４５０７を含む場合もあり、さらに、ＲＡＭ４５０１を含む場合もある。本技術の実施の形態では、上で述べた処理を実施するための制御プログラムは、コンピュータ読み取り可能なリムーバブル・ディスクに格納されて頒布され、ＲＯＭライタによってＲＯＭ４５０７に書き込まれる場合もある。このようなコンピュータ装置は、上で述べたプロセッサ４５０３、ＲＡＭ４５０１、ＲＯＭ４５０７などのハードウエアと制御プログラム（場合によってはＯＳも）とが有機的に協働することにより、上で述べたような各種機能を実現する。
【０２２２】
以上述べた本実施の形態をまとめると以下のようになる。
【０２２３】
本実施の形態の第１の態様に係る行列生成方法は、（Ａ）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）出力状態の目標値と出力状態との誤差の二乗及び入力の二乗の重み付け加算に関する第１の関数と、第１のラグランジュ乗数と状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和に関するハミルトニアンの前記入力についての偏微分であって第１のラグランジュ乗数と第２のラグランジュ乗数と解入力との関係を表す第１の条件式と、
（４）第１のラグランジュ乗数とハミルトニアンの出力状態についての偏微分との関係を表す第２の条件式と、
（５）第１のラグランジュ乗数の予測期間後の値と、出力状態の目標値と予測期間後の解出力状態との差に関する関数との関係式と、
を格納するデータ格納部から、状態方程式及び初期条件式とを読み出して、状態方程式及び初期条件式を、各時刻における解出力状態を出力状態と解入力とで表す第１の一次式群に変換し、データ格納部に格納する処理と、（Ｂ）データ格納部から第２の条件式及び関係式を読み出して、第２の条件式及び関係式を、各時刻における第１のラグランジュ乗数を目標値と解出力状態とで表す第２の一次式群に変換し、データ格納部に格納する処理と、（Ｃ）データ格納部に格納されている第１の一次式群を第２の一次式群に代入することで、各時刻における第１のラグランジュ乗数を目標値と解入力と出力状態とで表す第３の一次式群に変換し、データ格納部に格納する処理と、（Ｄ）データ格納部から第１の条件式を読み出して、第２のラグランジュ乗数の項を除外した第１の条件式に、データ格納部に格納されている第３の一次式群を代入することで、目標値と解入力と出力状態との各時刻における関係を表す第４の一次式群に変換し、データ格納部に格納する処理と、（Ｅ）データ格納部に格納されている第４の一次式群を、係数行列と、制御期間内における解入力についてのベクトルとの積が、目標値と出力状態とについての関数ベクトルと等式で表されるように変形して、係数行列と関数ベクトルとをデータ格納部に格納する処理とを含む。
【０２２４】
このようにすることで係数行列及び関数ベクトルが大幅にコンパクト化されており、係数行列と制御期間内における解入力についてのベクトルとの積が関数ベクトルと等式で結ぶことによって表される連立一次方程式を高速に解くことができるようになる。
【０２２５】
また、実施の形態の第１の態様に係る行列生成方法は、（Ｆ）
制約条件に含まれる上記入力の下限値に関連する第２のラグランジュ乗数を解入力の第１の一次式で近似した場合における当該第１の一次式の定数項であって制約条件に含まれる上記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む定数項を、各時刻について含む第１の制約ベクトルと、
第１の一次式における解入力の係数であって制約条件に含まれる上記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第２の制約ベクトルと、
制約条件に含まれる上記入力の上限値に関連する第２のラグランジュ乗数を解入力の第２の一次式で近似した場合における当該第２の一次式の定数項であって制約条件に含まれる上記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む定数項を、各時刻について含む第３の制約ベクトルと、
第２の一次式における解入力の係数であって制約条件に含まれる上記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第４の制約ベクトルと、
を生成し、データ格納部に格納する処理をさらに含むようにしても良い。
【０２２６】
このようにすれば、簡単に制約条件を考慮した演算を実施できるようになる。
【０２２７】
実施の形態の第２の態様に係る行列生成方法は、（Ａ）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）出力状態の目標値と出力状態との誤差の二乗及び入力の二乗の重み付け加算に関する第１の関数と、第１のラグランジュ乗数と状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和に関するハミルトニアンの前記入力についての偏微分であって第１のラグランジュ乗数と第２のラグランジュ乗数と解入力との関係を表す第１の条件式と、
（４）第１のラグランジュ乗数とハミルトニアンの出力状態についての偏微分との関係を表す第２の条件式と、
（５）第１のラグランジュ乗数の予測期間後の値と、出力状態の目標値と予測期間後の解出力状態との差に関する関数との関係式と、
（６）制約条件に含まれる上記入力の下限値に関連する第２のラグランジュ乗数を解入力の第１の一次式で近似した場合における当該第１の一次式群と、制約条件に含まれる上記入力の上限値に関連する第２のラグランジュ乗数を解入力の第２の一次式で近似した場合における当該第２の一次式群と、
を格納するデータ格納部から、状態方程式及び初期条件式とを読み出して、状態方程式及び初期条件式を、各時刻における解出力状態を出力状態と解入力とで表す第３の一次式群に変換し、データ格納部に格納する処理と、（Ｂ）データ格納部から第２の条件式及び関係式を読み出して、第２の条件式及び関係式を、各時刻における第１のラグランジュ乗数を前記目標値と前記解出力状態とで表す第４の一次式群に変換し、データ格納部に格納する処理と、（Ｃ）データ格納部に格納されている第３の一次式群を第４の一次式群に代入することで、各時刻における第１のラグランジュ乗数を目標値と解入力と出力状態とで表す第５の一次式群に変換し、データ格納部に格納する処理と、（Ｄ）データ格納部に格納されている第３の一次式群に対して、データ格納部に格納されている第５の一次式群、第１の一次式群及び第２の一次式群を代入することで、目標値と解入力と出力状態との各時刻における関係を表す第６の一次式群に変換し、データ格納部に格納する処理と、（Ｅ）データ格納部に格納されている第６の一次式群を、係数行列と、制御期間内における解入力についてのベクトルとの積が、目標値と出力状態とについての関数ベクトルと等式で表されるように変形して、係数行列と関数ベクトルとをデータ格納部に格納する処理とを含む。
【０２２８】
制約条件を考慮する場合には、このような形で（１）乃至（６）の式を変形しても、係数行列及び関数ベクトルとが大幅にコンパクト化されるようになる。
【０２２９】
本実施の形態の第３の態様に係るプラント制御方法は、（Ａ）制御対象であるプラントの現在の出力状態と現在の目標値とを取得する処理と、（Ｂ）
（１）プラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）出力状態の目標値と出力状態との誤差の二乗及び入力の二乗の重み付け加算を表す第１の関数と、第１のラグランジュ乗数と状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和で表されるハミルトニアンの前記入力についての偏微分であって第１のラグランジュ乗数と第２のラグランジュ乗数と解入力との関係を表す第１の条件式と、
（４）第１のラグランジュ乗数とハミルトニアンの出力状態についての偏微分との関係を表す第２の条件式と、
（５）第１のラグランジュ乗数の予測期間後の値と、出力状態の目標値と予測期間後の解出力状態との差に関する関数との関係式と、
等価な一次方程式群であって、係数行列と、制御期間内における解入力についてのベクトルとの積が、目標値と出力状態とについての関数ベクトルと等しくなるように表されている一次方程式群に対して、現在の出力状態と現在の目標値とを代入して解くことによって、制御期間内における解入力を算出する算出処理と、（Ｃ）制御期間内における解入力のうち最初の時刻のための解入力を、制御対象のプラントに出力する処理とを含む。
【０２３０】
一次方程式群は、元の（１）乃至（５）の式と比して大幅にコンパクト化されており、通常どおり演算を行って一次方程式群を解いたとしても短時間で解入力を得ることができるようになる。
【０２３１】
上で述べた算出処理が、（ｂ１）制御対象のプラントの解入力に対して上限値及び下限値を含む制約条件が設定されている場合に、
（６）制約条件に含まれる入力の下限値に関連する第２のラグランジュ乗数を解入力の第１の一次式で近似した場合における当該第１の一次式における解入力の係数であって制約条件に含まれる上記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第１の制約ベクトルに対して、１単位時間前に計算された対応時刻の解入力を設定し、対角行列化することで得られる第１の行列と、
（７）制約条件に含まれる上記入力の上限値に関連する第２のラグランジュ乗数を解入力の第２の一次式で近似した場合における当該第２の一次式における解入力の係数であって制約条件に含まれる上記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第２の制約ベクトルに対して、１単位時間前に計算された対応時刻の解入力を設定し、対角行列化することで得られる第２の行列と、
を生成し、係数行列から第１の行列を差し引き、第２の行列を加算する処理と、（ｂ２）
（８）制約条件に含まれる上記入力の下限値に関連する第２のラグランジュ乗数を解入力の第１の一次式で近似した場合における当該第１の一次式における定数項であって制約条件に含まれる上記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む定数項に対して１単位時間前に計算された対応時刻の解入力を設定して各時刻について並べた第１のベクトルと、
（９）制約条件に含まれる上記入力の上限値に関連する第２のラグランジュ乗数を解入力の第２の一次式で近似した場合における当該第２の一次式における定数項であって制約条件に含まれる上記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む定数項に対して１単位時間前に計算された対応時刻の解入力を設定して各時刻について並べた第２のベクトルと、
を生成し、関数ベクトルに第１ベクトルを加算して第２ベクトルを差し引く処理とを含むようにしても良い。
【０２３２】
このようにすれば、少ない演算で制約条件を考慮した形で解入力を出力できるようになる。
【０２３３】
本実施の形態の第４の態様に係るプラント制御装置（図２６）は、（Ａ）制御対象であるプラントの現在の出力状態と現在の目標値とを取得する取得部（図２６：３１００）と、（Ｂ）
（１）プラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）出力状態の目標値と出力状態との誤差の二乗及び入力の二乗の重み付け加算を表す第１の関数と、第１のラグランジュ乗数と状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和で表されるハミルトニアンの上記入力についての偏微分であって第１のラグランジュ乗数と第２のラグランジュ乗数と解入力との関係を表す第１の条件式と、
（４）前記第１のラグランジュ乗数と前記ハミルトニアンの前記出力状態についての偏微分との関係を表す第２の条件式と、
（５）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
等価な一次方程式群であって、係数行列と、制御期間内における解入力についてのベクトルとの積が、目標値と出力状態とについての関数ベクトルと等しくなるように表されている一次方程式群を格納するデータ格納部（図２６：３４００）と、（Ｃ）データ格納部から一次方程式群を読み出して、一次方程式群に対して、現在の出力状態と現在の目標値とを代入して解くことによって、制御期間内における解入力を算出する演算部（図２６：３２００）と、（Ｄ）制御期間内における解入力のうち最初の時刻のための解入力を、制御対象のプラントに出力する出力部（図２６：３３００）とを有する。
【０２３４】
本実施の形態の第５の態様に係る行列生成方法は、（Ａ）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）出力状態の目標値と解出力状態との誤差の二乗を第１の重み行列で重み付けして初期時刻から１単位時間後の時点から予測区間後の１単位時間前の時点まで加算する第１の項と、出力状態の目標値と予測区間後における解出力状態との誤差の二乗を第２の重み行列で重み付けして加算する第２の項と、解入力の二乗を第３の重み行列で重み付けして初期時刻から予測区間後の１単位時間前の時点まで加算する第３の項と、制約条件についての第４の項と、を含む評価関数のうち第１の項と第２の項とを、初期時刻における解入力から制御区間後の１単位時間前の時点における解入力までの各解入力で偏微分することで得られる、解出力状態についての第１の一次式群と、
（４）第３の項を初期時刻における解入力から制御区間後の１単位時間前の時点における解入力までの各解入力で偏微分することで得られ且つ第３の重み行列の要素値、制御区間の値及び予測区間の値とで決定される行列要素を有する第１の行列と、
を格納するデータ格納部から、状態方程式及び初期条件式とを読み出して、状態方程式及び初期条件式を、各時刻における解出力状態を出力状態と解入力とで表す第２の一次式群に変換し、データ格納部に格納する処理と、（Ｂ）データ格納部に格納されている第２の一次式群及び第２の一次式群における解入力の係数値を、データ格納部に格納されている第１の一次式群に代入することで第３の一次式群を生成し、データ格納部に格納する処理と、（Ｃ）データ格納部に格納されている第３の一次式群における各解入力の項の係数を抽出して第２の行列を生成し、第３の一次式群における解入力の項以外の項を抽出してベクトルを生成し、データ格納部に格納する処理と、（Ｄ）データ格納部に格納されている第１の行列と第２の行列とを加算して、第３の行列を生成し、データ格納部に格納する処理とを含む。
【０２３５】
このような処理を行えば演算量が削減されているので、高速に第３の行列及びベクトルを生成することができるようになる。
【０２３６】
本実施の形態の第６の態様に係る行列生成方法は、（Ａ）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）出力状態の目標値と解出力状態との誤差の二乗を第１の重み行列で重み付けして初期時刻から１単位時間後の時点から予測区間後の１単位時間前の時点まで加算する第１の項と、出力状態の目標値と予測区間後における解出力状態との誤差の二乗を第２の重み行列で重み付けして加算する第２の項と、解入力の二乗を第３の重み行列で重み付けして初期時刻から予測区間後の１単位時間前の時点まで加算する第３の項と、制約条件についての第４の項と、を含む評価関数のうち第１の項及び第２の項を、係数として乗じられている、第１の重み行列及び第２の重み行列の要素毎に分類することで得られ且つ当該要素が除外された複数の式の各々を、初期時刻における解入力から制御区間後の１単位時間前の時点における解入力までの各解入力で偏微分することで得られる、解出力状態についての第１の一次式群と、
を格納するデータ格納部から、状態方程式及び初期条件式とを読み出して、状態方程式及び初期条件式を、各時刻における解出力状態を出力状態と解入力とで表す第２の一次式群に変換し、データ格納部に格納する処理と、（Ｂ）データ格納部に格納されている第２の一次式群及び第２の一次式群における解入力の係数値を、データ格納部に格納されている第１の一次式群に代入することで、複数の式の各々について第３の一次式群を生成し、データ格納部に格納する処理と、（Ｃ）データ格納部に格納されている第３の一次式群における各解入力の項の係数を複数の式の各々について抽出して複数の式の各々について行列を生成し、第３の一次式群における解入力の項以外の項を複数の式の各々について抽出して複数の式の各々についてベクトルを生成し、データ格納部に格納する処理とを含む。
【０２３７】
このようにすれば、行列及びベクトルを生成し直すことなく、後に実施するシミュレーションなどにおいて重み行列の要素値を調整することができるようになる。なお、制約条件については、本実施の形態の第５及び第６の態様においては、本実施の形態の第１の態様と同様に処理することができる。
【０２３８】
なお、上記方法による処理をプロセッサに行わせるためのプログラムを作成することができ、当該プログラムは、例えばフレキシブルディスク、ＣＤ−ＲＯＭ、光磁気ディスク、半導体メモリ、ハードディスク等のコンピュータ読み取り可能な記憶媒体又は記憶装置に格納される。尚、中間的な処理結果はメインメモリ等の記憶装置に一時保管される。
【０２３９】
以上の実施例を含む実施形態に関し、さらに以下の付記を開示する。
【０２４０】
（付記１）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算に関する第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和に関するハミルトニアンの前記入力についての偏微分であって前記第１のラグランジュ乗数と前記第２のラグランジュ乗数と前記解入力との関係を表す第１の条件式と、
（４）前記第１のラグランジュ乗数と前記ハミルトニアンの前記出力状態についての偏微分との関係を表す第２の条件式と、
（５）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
を格納するデータ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第１の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部から前記第２の条件式及び前記関係式を読み出して、前記第２の条件式及び前記関係式を、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解出力状態とで表す第２の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第１の一次式群を前記第２の一次式群に代入することで、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解入力と前記出力状態とで表す第３の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部から前記第１の条件式を読み出して、前記第２のラグランジュ乗数の項を除外した第１の条件式に、前記データ格納部に格納されている前記第３の一次式群を代入することで、前記目標値と前記解入力と前記出力状態との各時刻における関係を表す第４の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第４の一次式群を、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等式で表されるように変形して、前記係数行列と前記関数ベクトルとを前記データ格納部に格納する処理と、
を、コンピュータに実行させるためのプログラム。
【０２４１】
（付記２）
前記制約条件に含まれる前記入力の下限値に関連する前記第２のラグランジュ乗数を前記解入力の第１の一次式で近似した場合における当該第１の一次式の定数項であって前記制約条件に含まれる前記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む定数項を、各時刻について含む第１の制約ベクトルと、
前記第１の一次式における前記解入力の係数であって前記制約条件に含まれる前記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第２の制約ベクトルと、
前記制約条件に含まれる前記入力の上限値に関連する前記第２のラグランジュ乗数を前記解入力の第２の一次式で近似した場合における当該第２の一次式の定数項であって前記制約条件に含まれる前記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む定数項を、各時刻について含む第３の制約ベクトルと、
前記第２の一次式における前記解入力の係数であって前記制約条件に含まれる前記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第４の制約ベクトルと、
を生成し、前記データ格納部に格納する処理、
をさらに前記コンピュータに実行させるための付記１記載のプログラム。
【０２４２】
（付記３）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算に関する第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和に関するハミルトニアンの前記入力についての偏微分であって前記第１のラグランジュ乗数と前記第２のラグランジュ乗数と前記解入力との関係を表す第１の条件式と、
（４）前記第１のラグランジュ乗数と前記ハミルトニアンの前記出力状態についての偏微分との関係を表す第２の条件式と、
（５）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
（６）前記制約条件に含まれる前記入力の下限値に関連する前記第２のラグランジュ乗数を前記解入力の第１の一次式で近似した場合における当該第１の一次式群と、前記制約条件に含まれる前記入力の上限値に関連する前記第２のラグランジュ乗数を前記解入力の第２の一次式で近似した場合における当該第２の一次式群と、
を格納するデータ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第３の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部から前記第２の条件式及び前記関係式を読み出して、前記第２の条件式及び前記関係式を、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解出力状態とで表す第４の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第３の一次式群を前記第４の一次式群に代入することで、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解入力と前記出力状態とで表す第５の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第３の一次式群に対して、前記データ格納部に格納されている前記第５の一次式群、前記第１の一次式群及び前記第２の一次式群を代入することで、前記目標値と前記解入力と前記出力状態との各時刻における関係を表す第６の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第６の一次式群を、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等式で表されるように変形して、前記係数行列と前記関数ベクトルとを前記データ格納部に格納する処理と、
を、コンピュータに実行させるためのプログラム。
【０２４３】
（付記４）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算に関する第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和に関するハミルトニアンの前記入力についての偏微分であって前記第１のラグランジュ乗数と前記第２のラグランジュ乗数と前記解入力との関係を表す第１の条件式と、
（４）前記第１のラグランジュ乗数と前記ハミルトニアンの前記出力状態についての偏微分との関係を表す第２の条件式と、
（５）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
を格納するデータ格納部と、
前記データ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第１の一次式群に変換し、前記データ格納部に格納する第１変形処理部と、
前記データ格納部から前記第２の条件式及び前記関係式を読み出して、前記第２の条件式及び前記関係式を、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解出力状態とで表す第２の一次式群に変換し、前記データ格納部に格納する第２変形処理部と、
前記データ格納部に格納されている前記第１の一次式群を前記第２の一次式群に代入することで、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解入力と前記出力状態とで表す第３の一次式群に変換し、前記データ格納部に格納する第３変形処理部と、
前記データ格納部から前記第１の条件式を読み出して、前記第２のラグランジュ乗数の項を除外した第１の条件式に、前記データ格納部に格納されている前記第３の一次式群を代入することで、前記目標値と前記解入力と前記出力状態との各時刻における関係を表す第４の一次式群に変換し、前記データ格納部に格納する式生成部と、
前記データ格納部に格納されている前記第４の一次式群を、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等式で表されるように変形して、前記係数行列と前記関数ベクトルとを前記データ格納部に格納する行列生成部と、
を有する行列生成装置。
【０２４４】
（付記５）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算に関する第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和に関するハミルトニアンの前記入力についての偏微分であって前記第１のラグランジュ乗数と前記第２のラグランジュ乗数と前記解入力との関係を表す第１の条件式と、
（４）前記第１のラグランジュ乗数と前記ハミルトニアンの前記出力状態についての偏微分との関係を表す第２の条件式と、
（５）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
（６）前記制約条件に含まれる前記入力の下限値に関連する前記第２のラグランジュ乗数を前記解入力の第１の一次式で近似した場合における当該第１の一次式群と、前記制約条件に含まれる前記入力の上限値に関連する前記第２のラグランジュ乗数を前記解入力の第２の一次式で近似した場合における当該第２の一次式群と、
を格納するデータ格納部と、
前記データ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第３の一次式群に変換し、前記データ格納部に格納する第１変形処理部と、
前記データ格納部から前記第２の条件式及び前記関係式を読み出して、前記第２の条件式及び前記関係式を、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解出力状態とで表す第４の一次式群に変換し、前記データ格納部に格納する第２変形処理部と、
前記データ格納部に格納されている前記第３の一次式群を前記第４の一次式群に代入することで、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解入力と前記出力状態とで表す第５の一次式群に変換し、前記データ格納部に格納する第３変形処理部と、
前記データ格納部に格納されている前記第３の一次式群に対して、前記データ格納部に格納されている前記第５の一次式群、前記第１の一次式群及び前記第２の一次式群を代入することで、前記目標値と前記解入力と前記出力状態との各時刻における関係を表す第６の一次式群に変換し、前記データ格納部に格納する式生成部と、
前記データ格納部に格納されている前記第６の一次式群を、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等式で表されるように変形して、前記係数行列と前記関数ベクトルとを前記データ格納部に格納する行列生成部と、
を有する行列生成装置。
【０２４５】
（付記６）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算に関する第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和に関するハミルトニアンの前記入力についての偏微分であって前記第１のラグランジュ乗数と前記第２のラグランジュ乗数と前記解入力との関係を表す第１の条件式と、
（４）前記第１のラグランジュ乗数と前記ハミルトニアンの前記出力状態についての偏微分との関係を表す第２の条件式と、
（５）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
を格納するデータ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第１の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部から前記第２の条件式及び前記関係式を読み出して、前記第２の条件式及び前記関係式を、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解出力状態とで表す第２の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第１の一次式群を前記第２の一次式群に代入することで、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解入力と前記出力状態とで表す第３の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部から前記第１の条件式を読み出して、前記第２のラグランジュ乗数の項を除外した第１の条件式に、前記データ格納部に格納されている前記第３の一次式群を代入することで、前記目標値と前記解入力と前記出力状態との各時刻における関係を表す第４の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第４の一次式群を、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等式で表されるように変形して、前記係数行列と前記関数ベクトルとを前記データ格納部に格納する処理と、
を、コンピュータが実行する行列生成方法。
【０２４６】
（付記７）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算に関する第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和に関するハミルトニアンの前記入力についての偏微分であって前記第１のラグランジュ乗数と前記第２のラグランジュ乗数と前記解入力との関係を表す第１の条件式と、
（４）前記第１のラグランジュ乗数と前記ハミルトニアンの前記出力状態についての偏微分との関係を表す第２の条件式と、
（５）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
（６）前記制約条件に含まれる前記入力の下限値に関連する前記第２のラグランジュ乗数を前記解入力の第１の一次式で近似した場合における当該第１の一次式群と、前記制約条件に含まれる前記入力の上限値に関連する前記第２のラグランジュ乗数を前記解入力の第２の一次式で近似した場合における当該第２の一次式群と、
を格納するデータ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第３の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部から前記第２の条件式及び前記関係式を読み出して、前記第２の条件式及び前記関係式を、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解出力状態とで表す第４の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第３の一次式群を前記第４の一次式群に代入することで、各時刻における前記第１のラグランジュ乗数を前記目標値と前記解入力と前記出力状態とで表す第５の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第３の一次式群に対して、前記データ格納部に格納されている前記第５の一次式群、前記第１の一次式群及び前記第２の一次式群を代入することで、前記目標値と前記解入力と前記出力状態との各時刻における関係を表す第６の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第６の一次式群を、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等式で表されるように変形して、前記係数行列と前記関数ベクトルとを前記データ格納部に格納する処理と、
を、コンピュータが実行する行列生成方法。
【０２４７】
（付記８）
制御対象であるプラントの現在の出力状態と現在の目標値とを取得する処理と、
（１）前記プラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算を表す第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和で表されるハミルトニアンの前記入力についての偏微分であって前記第１のラグランジュ乗数と前記第２のラグランジュ乗数と前記解入力との関係を表す第１の条件式と、
（４）前記第１のラグランジュ乗数と前記ハミルトニアンの前記出力状態についての偏微分との関係を表す第２の条件式と、
（５）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
等価な一次方程式群であって、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等しくなるように表されている前記一次方程式群に対して、前記現在の出力状態と前記現在の目標値とを代入して解くことによって、前記制御期間内における解入力を算出する算出処理と、
前記制御期間内における解入力のうち最初の時刻のための解入力を、前記制御対象のプラントに出力する処理と、
を、プロセッサに実行させるためのプログラム。
【０２４８】
（付記９）
前記算出処理が、
前記制御対象のプラントの解入力に対して上限値及び下限値を含む制約条件が設定されている場合に、
前記制約条件に含まれる前記入力の下限値に関連する前記第２のラグランジュ乗数を前記解入力の第１の一次式で近似した場合における当該第１の一次式における前記解入力の係数であって前記制約条件に含まれる前記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第１の制約ベクトルに対して、前記１単位時間前に計算された対応時刻の解入力を設定し、対角行列化することで得られる第１の行列と、
前記制約条件に含まれる前記入力の上限値に関連する前記第２のラグランジュ乗数を前記解入力の第２の一次式で近似した場合における当該第２の一次式における前記解入力の係数であって前記制約条件に含まれる前記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第２の制約ベクトルに対して、前記１単位時間前に計算された対応時刻の解入力を設定し、対角行列化することで得られる第２の行列と、
を生成し、前記係数行列から前記第１の行列を差し引き、前記第２の行列を加算する処理と、
前記制約条件に含まれる前記入力の下限値に関連する前記第２のラグランジュ乗数を前記解入力の第１の一次式で近似した場合における当該第１の一次式における定数項であって前記制約条件に含まれる前記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む定数項に対して前記１単位時間前に計算された対応時刻の解入力を設定して各時刻について並べた第１のベクトルと、
前記制約条件に含まれる前記入力の上限値に関連する前記第２のラグランジュ乗数を前記解入力の第２の一次式で近似した場合における当該第２の一次式における定数項であって前記制約条件に含まれる前記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む定数項に対して前記１単位時間前に計算された対応時刻の解入力を設定して各時刻について並べた第２のベクトルと、
を生成し、前記関数ベクトルに前記第１ベクトルを加算して前記第２ベクトルを差し引く処理と、
を含む付記８記載のプログラム。
【０２４９】
（付記１０）
制御対象であるプラントの現在の出力状態と現在の目標値とを取得する取得部と、
（１）前記プラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算を表す第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和で表されるハミルトニアンの前記入力についての偏微分であって前記第１のラグランジュ乗数と前記第２のラグランジュ乗数と前記解入力との関係を表す第１の条件式と、
（４）前記第１のラグランジュ乗数と前記ハミルトニアンの前記出力状態についての偏微分との関係を表す第２の条件式と、
（５）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
等価な一次方程式群であって、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等しくなるように表されている前記一次方程式群を格納するデータ格納部と、
前記データ格納部から前記一次方程式群を読み出して、前記一次方程式群に対して、前記現在の出力状態と前記現在の目標値とを代入して解くことによって、前記制御期間内における解入力を算出する演算部と、
前記制御期間内における解入力のうち最初の時刻のための解入力を、前記制御対象のプラントに出力する出力部と、
を有するプラント制御装置。
【０２５０】
（付記１１）
制御対象であるプラントの現在の出力状態と現在の目標値とを取得する処理と、
（１）前記プラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記出力状態との誤差の二乗及び入力の二乗の重み付け加算を表す第１の関数と、第１のラグランジュ乗数と前記状態方程式との積と、第２のラグランジュ乗数と制約条件に関する第２の関数との積と、の和で表されるハミルトニアンの前記入力についての偏微分であって前記第１のラグランジュ乗数と前記第２のラグランジュ乗数と前記解入力との関係を表す第１の条件式と、
（４）前記第１のラグランジュ乗数と前記ハミルトニアンの前記出力状態についての偏微分との関係を表す第２の条件式と、
（５）前記第１のラグランジュ乗数の予測期間後の値と、前記出力状態の目標値と前記予測期間後の解出力状態との差に関する関数との関係式と、
等価な一次方程式群であって、係数行列と、制御期間内における前記解入力についてのベクトルとの積が、前記目標値と前記出力状態とについての関数ベクトルと等しくなるように表されている前記一次方程式群に対して、前記現在の出力状態と前記現在の目標値とを代入して解くことによって、前記制御期間内における解入力を算出する算出処理と、
前記制御期間内における解入力のうち最初の時刻のための解入力を、前記制御対象のプラントに出力する処理と、
を、プロセッサが実行する制御方法。
【０２５１】
（付記１２）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記解出力状態との誤差の二乗を第１の重み行列で重み付けして初期時刻から１単位時間後の時点から予測区間後の１単位時間前の時点まで加算する第１の項と、前記出力状態の目標値と前記予測区間後における前記解出力状態との誤差の二乗を第２の重み行列で重み付けして加算する第２の項と、前記解入力の二乗を第３の重み行列で重み付けして初期時刻から前記予測区間後の１単位時間前の時点まで加算する第３の項と、制約条件についての第４の項と、を含む評価関数のうち前記第１の項と前記第２の項とを、初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られる、前記解出力状態についての第１の一次式群と、
（４）前記第３の項を前記初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られ且つ前記第３の重み行列の要素値、前記制御区間の値及び前記予測区間の値とで決定される行列要素を有する第１の行列と、
を格納するデータ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第２の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第２の一次式群及び前記第２の一次式群における前記解入力の係数値を、前記データ格納部に格納されている前記第１の一次式群に代入することで第３の一次式群を生成し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第３の一次式群における各前記解入力の項の係数を抽出して第２の行列を生成し、前記第３の一次式群における前記解入力の項以外の項を抽出してベクトルを生成し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第１の行列と前記第２の行列とを加算して、第３の行列を生成し、前記データ格納部に格納する処理と、
を、コンピュータに実行させるためのプログラム。
【０２５２】
（付記１３）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記解出力状態との誤差の二乗を第１の重み行列で重み付けして初期時刻から１単位時間後の時点から予測区間後の１単位時間前の時点まで加算する第１の項と、前記出力状態の目標値と前記予測区間後における前記解出力状態との誤差の二乗を第２の重み行列で重み付けして加算する第２の項と、前記解入力の二乗を第３の重み行列で重み付けして初期時刻から前記予測区間後の１単位時間前の時点まで加算する第３の項と、制約条件についての第４の項と、を含む評価関数のうち前記第１の項及び前記第２の項を、係数として乗じられている、前記第１の重み行列及び前記第２の重み行列の要素毎に分類することで得られ且つ当該要素が除外された複数の式の各々を、初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られる、前記解出力状態についての第１の一次式群と、
を格納するデータ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第２の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第２の一次式群及び前記第２の一次式群における前記解入力の係数値を、前記データ格納部に格納されている前記第１の一次式群に代入することで、前記複数の式の各々について第３の一次式群を生成し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第３の一次式群における各前記解入力の項の係数を前記複数の式の各々について抽出して前記複数の式の各々について行列を生成し、前記第３の一次式群における前記解入力の項以外の項を前記複数の式の各々について抽出して前記複数の式の各々についてベクトルを生成し、前記データ格納部に格納する処理と、
を、コンピュータに実行させるためのプログラム。
【０２５３】
（付記１４）
前記制約条件に含まれる前記入力の下限値に関連する前記第２のラグランジュ乗数を前記解入力の第１の一次式で近似した場合における当該第１の一次式の定数項であって前記制約条件に含まれる前記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む定数項を、各時刻について含む第１の制約ベクトルと、
前記第１の一次式における前記解入力の係数であって前記制約条件に含まれる前記入力の下限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第２の制約ベクトルと、
前記制約条件に含まれる前記入力の上限値に関連する前記第２のラグランジュ乗数を前記解入力の第２の一次式で近似した場合における当該第２の一次式の定数項であって前記制約条件に含まれる前記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む定数項を、各時刻について含む第３の制約ベクトルと、
前記第２の一次式における前記解入力の係数であって前記制約条件に含まれる前記入力の上限値と１単位時間前に計算された対応時刻の解入力とを含む係数を、各時刻について含む第４の制約ベクトルと、
を生成し、前記データ格納部に格納する処理、
をさらに前記コンピュータに実行させるための付記１２又は１３記載のプログラム。
【０２５４】
（付記１５）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記解出力状態との誤差の二乗を第１の重み行列で重み付けして初期時刻から１単位時間後の時点から予測区間後の１単位時間前の時点まで加算する第１の項と、前記出力状態の目標値と前記予測区間後における前記解出力状態との誤差の二乗を第２の重み行列で重み付けして加算する第２の項と、前記解入力の二乗を第３の重み行列で重み付けして初期時刻から前記予測区間後の１単位時間前の時点まで加算する第３の項と、制約条件についての第４の項と、を含む評価関数のうち前記第１の項と前記第２の項とを、初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られる、前記解出力状態についての第１の一次式群と、
（４）前記第３の項を前記初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られ且つ前記第３の重み行列の要素値、前記制御区間の値及び前記予測区間の値とで決定される行列要素を有する第１の行列と、
を格納するデータ格納部と、
前記データ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第２の一次式群に変換し、前記データ格納部に格納する変形部と、
前記データ格納部に格納されている前記第２の一次式群及び前記第２の一次式群における前記解入力の係数値を、前記データ格納部に格納されている前記第１の一次式群に代入することで第３の一次式群を生成し、前記データ格納部に格納し、前記データ格納部に格納されている前記第３の一次式群における各前記解入力の項の係数を抽出して第２の行列を生成し、前記第３の一次式群における前記解入力の項以外の項を抽出してベクトルを生成し、前記データ格納部に格納し、前記データ格納部に格納されている前記第１の行列と前記第２の行列とを加算して、第３の行列を生成し、前記データ格納部に格納する生成部と、
を有する行列生成装置。
【０２５５】
（付記１６）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記解出力状態との誤差の二乗を第１の重み行列で重み付けして初期時刻から１単位時間後の時点から予測区間後の１単位時間前の時点まで加算する第１の項と、前記出力状態の目標値と前記予測区間後における前記解出力状態との誤差の二乗を第２の重み行列で重み付けして加算する第２の項と、前記解入力の二乗を第３の重み行列で重み付けして初期時刻から前記予測区間後の１単位時間前の時点まで加算する第３の項と、制約条件についての第４の項と、を含む評価関数のうち前記第１の項及び前記第２の項を、係数として乗じられている、前記第１の重み行列及び前記第２の重み行列の要素毎に分類することで得られ且つ当該要素が除外された複数の式の各々を、初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られる、前記解出力状態についての第１の一次式群と、
を格納するデータ格納部と、
前記データ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第２の一次式群に変換し、前記データ格納部に格納する変形部と、
前記データ格納部に格納されている前記第２の一次式群及び前記第２の一次式群における前記解入力の係数値を、前記データ格納部に格納されている前記第１の一次式群に代入することで、前記複数の式の各々について第３の一次式群を生成し、前記データ格納部に格納し、前記データ格納部に格納されている前記第３の一次式群における各前記解入力の項の係数を前記複数の式の各々について抽出して前記複数の式の各々について行列を生成し、前記第３の一次式群における前記解入力の項以外の項を前記複数の式の各々について抽出して前記複数の式の各々についてベクトルを生成し、前記データ格納部に格納する生成部と、
を有する行列生成装置。
【０２５６】
（付記１７）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記解出力状態との誤差の二乗を第１の重み行列で重み付けして初期時刻から１単位時間後の時点から予測区間後の１単位時間前の時点まで加算する第１の項と、前記出力状態の目標値と前記予測区間後における前記解出力状態との誤差の二乗を第２の重み行列で重み付けして加算する第２の項と、前記解入力の二乗を第３の重み行列で重み付けして初期時刻から前記予測区間後の１単位時間前の時点まで加算する第３の項と、制約条件についての第４の項と、を含む評価関数のうち前記第１の項と前記第２の項とを、初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られる、前記解出力状態についての第１の一次式群と、
（４）前記第３の項を前記初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られ且つ前記第３の重み行列の要素値、前記制御区間の値及び前記予測区間の値とで決定される行列要素を有する第１の行列と、
を格納するデータ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第２の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第２の一次式群及び前記第２の一次式群における前記解入力の係数値を、前記データ格納部に格納されている前記第１の一次式群に代入することで第３の一次式群を生成し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第３の一次式群における各前記解入力の項の係数を抽出して第２の行列を生成し、前記第３の一次式群における前記解入力の項以外の項を抽出してベクトルを生成し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第１の行列と前記第２の行列とを加算して、第３の行列を生成し、前記データ格納部に格納する処理と、
を、コンピュータが実行する行列生成方法。
【０２５７】
（付記１８）
（１）制御対象であるプラントの解出力状態と解入力との関係を表す状態方程式と、
（２）出力状態の初期条件式と、
（３）前記出力状態の目標値と前記解出力状態との誤差の二乗を第１の重み行列で重み付けして初期時刻から１単位時間後の時点から予測区間後の１単位時間前の時点まで加算する第１の項と、前記出力状態の目標値と前記予測区間後における前記解出力状態との誤差の二乗を第２の重み行列で重み付けして加算する第２の項と、前記解入力の二乗を第３の重み行列で重み付けして初期時刻から前記予測区間後の１単位時間前の時点まで加算する第３の項と、制約条件についての第４の項と、を含む評価関数のうち前記第１の項及び前記第２の項を、係数として乗じられている、前記第１の重み行列及び前記第２の重み行列の要素毎に分類することで得られ且つ当該要素が除外された複数の式の各々を、初期時刻における前記解入力から前記制御区間後の１単位時間前の時点における前記解入力までの各前記解入力で偏微分することで得られる、前記解出力状態についての第１の一次式群と、
を格納するデータ格納部から、前記状態方程式及び前記初期条件式とを読み出して、前記状態方程式及び前記初期条件式を、各時刻における前記解出力状態を前記出力状態と前記解入力とで表す第２の一次式群に変換し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第２の一次式群及び前記第２の一次式群における前記解入力の係数値を、前記データ格納部に格納されている前記第１の一次式群に代入することで、前記複数の式の各々について第３の一次式群を生成し、前記データ格納部に格納する処理と、
前記データ格納部に格納されている前記第３の一次式群における各前記解入力の項の係数を前記複数の式の各々について抽出して前記複数の式の各々について行列を生成し、前記第３の一次式群における前記解入力の項以外の項を前記複数の式の各々について抽出して前記複数の式の各々についてベクトルを生成し、前記データ格納部に格納する処理と、
を、コンピュータが実行する行列生成方法。
【符号の説明】
【０２５８】
１００  オフライン行列生成装置
１０１  入力部
１０２  入力データ格納部
１０３  第１変形処理部
１０４  第１データ格納部
１０５  第２変形処理部
１０６  第２データ格納部
１０７  第３変形処理部
１０８  第３データ格納部
１０９  連立方程式生成部
１１０  第４データ格納部
１１１  行列生成部
１１２  第５データ格納部
５０１  第４変形処理部
２１０  計画器
２２０  第１モデル予測制御部
２３０  プラント
２４０  第２モデル予測制御部
５０２  第２連立方程式生成部
７０１  第２行列生成部
７０２  第９データ格納部
７０３  代入処理部
７０４  第６データ格納部
７０５  第１行列生成部
７０６  第７データ格納部
７０７  第３行列生成部
７０８  第８データ格納部
９０１  第２代入処理部
９０２  第１１データ格納部
９０３  第４行列生成部
９０４  第１２データ格納部
９０５  第５行列生成部
９０６  第１３データ格納部
【図１】
0005834759.tif 000099
【図２】
0005834759.tif 000100
【図３】
0005834759.tif 000101
【図４】
0005834759.tif 000102
【図５】
0005834759.tif 000103
【図６】
0005834759.tif 000104
【図７】
0005834759.tif 000105
【図８】
0005834759.tif 000106
【図９】
0005834759.tif 000107
【図１０】
0005834759.tif 000108
【図１１】
0005834759.tif 000109
【図１２】
0005834759.tif 000110
【図１３】
0005834759.tif 000111
【図１４】
0005834759.tif 000112
【図１５】
0005834759.tif 000113
【図１６】
0005834759.tif 000114
【図１７】
0005834759.tif 000115
【図１８】
0005834759.tif 000116
【図１９】
0005834759.tif 000117
【図２０】
0005834759.tif 000118
【図２１】
0005834759.tif 000119
【図２２】
0005834759.tif 000120
【図２３】
0005834759.tif 000121
【図２４】
0005834759.tif 000122
【図２５】
0005834759.tif 000123
【図２６】
0005834759.tif 000124