知財求人 - 知財ポータルサイト「IP Force」
▶ 株式会社東芝の特許一覧 ▶ 東芝電機サービス株式会社の特許一覧
特開2024-17830最適制御装置、最適制御方法およびコンピュータプログラム
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2024017830
(43)【公開日】2024-02-08
(54)【発明の名称】最適制御装置、最適制御方法およびコンピュータプログラム
(51)【国際特許分類】
G05B 13/02 20060101AFI20240201BHJP
【ＦＩ】
G05B13/02 C
G05B13/02 J
【審査請求】未請求
【請求項の数】11
【出願形態】ＯＬ
(21)【出願番号】P 2022120736
(22)【出願日】2022-07-28
(71)【出願人】
【識別番号】000003078
【氏名又は名称】株式会社東芝
(71)【出願人】
【識別番号】598076591
【氏名又は名称】東芝インフラシステムズ株式会社
(74)【代理人】
【識別番号】110003708
【氏名又は名称】弁理士法人鈴榮特許綜合事務所
(72)【発明者】
【氏名】山中 理
(72)【発明者】
【氏名】大西 祐太
(72)【発明者】
【氏名】西室 勇岐
【テーマコード（参考）】
5H004
【Ｆターム（参考）】
5H004GB01
5H004KC33
5H004KC44
5H004KC50
(57)【要約】
【課題】 安定性を維持しながら、制御性能を極力高める最適値の探索を実現する最適制御装置を提供する。
【解決手段】 実施形態による最適制御装置は、制御対象プロセスにおいて取得された計測値を用いて評価関数の評価値を算出するプロセス評価値算出部１２０と、操作量および評価値を入力として、操作量から評価値までの位相遅れの推定値を算出するプロセス位相遅れ推定部４００と、位相遅れの推定値に関する情報と評価値の情報とを用いて、評価値の操作量に対する変化率の推定値を算出する評価関数勾配推定部１４０と、変化率の推定値を積分することにより、操作量の動くべき方向と量とを決める極値探索部１７０と、極値探索部１７０で決定された操作量の動くべき方向と量との情報に基づく操作量を、制御対象プロセスへ出力する操作量出力部１９０と、を備える。
【選択図】図３
【特許請求の範囲】
【請求項1】
制御対象プロセスの操作量と、前記操作量に応じて変化する制御量に基づく評価関数の評価値とに基づいて、前記操作量をリアルタイムに操作して、前記評価値の最適値を探索する最適制御装置であって、
前記制御対象プロセスにおいて取得された計測値を用いて前記評価関数の前記評価値を算出するプロセス評価値算出部と、
前記操作量および前記評価値を用いて、前記操作量から前記評価値までの位相遅れの推定値を算出するプロセス位相遅れ推定部と、
前記位相遅れの推定値に関する情報と前記評価値の情報とを用いて、前記評価値の前記操作量に対する変化率の推定値を算出する評価関数勾配推定部と、
前記変化率の推定値を積分することにより、前記操作量の動くべき方向と量とを決める極値探索部と、
前記極値探索部で決定された前記操作量の動くべき方向と量との情報に基づく前記操作量を、前記制御対象プロセスへ出力する操作量出力部と、を備える最適制御装置。
前記制御対象プロセスにおいて取得された計測値を用いて前記評価関数の前記評価値を算出するプロセス評価値算出部と、
前記操作量および前記評価値を用いて、前記操作量から前記評価値までの位相遅れの推定値を算出するプロセス位相遅れ推定部と、
前記位相遅れの推定値に関する情報と前記評価値の情報とを用いて、前記評価値の前記操作量に対する変化率の推定値を算出する評価関数勾配推定部と、
前記変化率の推定値を積分することにより、前記操作量の動くべき方向と量とを決める極値探索部と、
前記極値探索部で決定された前記操作量の動くべき方向と量との情報に基づく前記操作量を、前記制御対象プロセスへ出力する操作量出力部と、を備える最適制御装置。
【請求項2】
前記プロセス位相遅れ推定部は、前記制御対象プロセスの前記操作量から前記評価値までの伝達関数を、静的な非線形要素を除く線形伝達関数モデルで表現し、線形伝達関数モデルの位相計算により前記位相遅れの推定値を算出する、請求項１記載の最適制御装置。
【請求項3】
前記プロセス位相遅れ推定部は、前記線形伝達関数モデルのパラメータを、オンラインで同定する位相遅れパラメータ推定部を備える、請求項２記載の最適制御装置。
【請求項4】
前記プロセス位相遅れ推定部は、前記制御対象プロセスの評価量から前記評価関数の前記評価値までのモデルを、前記制御対象プロセスのシミュレーションモデルと評価関数モデルとから構成し、前記モデルのシミュレーションにより、前記位相遅れの推定値を算出する、請求項１記載の最適制御装置。
【請求項5】
所定の周期を持つ第１周期信号を生成する復調用ディザー信号生成部と、
前記第１周期信号と同じ周期を持つ第２周期信号を生成する変調用ディザー信号生成部と、を有し、
前記操作量出力部は、前記極値探索部で決定された前記操作量の動くべき方向と量との情報に前記第２周期信号を加えた信号を出力し、
前記評価関数勾配推定部は、前記第２周期信号と前記評価値の情報とを用いて、前記操作量に対する前記評価値の前記変化率の推定値を算出し、
前記第１周期信号は、前記第２周期信号に前記位相遅れの推定値を加えた信号である、請求項２又は請求項４記載の最適制御装置。
前記第１周期信号と同じ周期を持つ第２周期信号を生成する変調用ディザー信号生成部と、を有し、
前記操作量出力部は、前記極値探索部で決定された前記操作量の動くべき方向と量との情報に前記第２周期信号を加えた信号を出力し、
前記評価関数勾配推定部は、前記第２周期信号と前記評価値の情報とを用いて、前記操作量に対する前記評価値の前記変化率の推定値を算出し、
前記第１周期信号は、前記第２周期信号に前記位相遅れの推定値を加えた信号である、請求項２又は請求項４記載の最適制御装置。
【請求項6】
前記変化率の推定値を正規化し、前記変化率の推定値の方向の情報を前記極値探索部へ出力する正規化部を更に備え、
前記極値探索部は、前記変化率の推定値の符号を積分することにより、前記操作量の動くべき方向と量とを決める、請求項５記載の最適制御装置。
前記極値探索部は、前記変化率の推定値の符号を積分することにより、前記操作量の動くべき方向と量とを決める、請求項５記載の最適制御装置。
【請求項7】
前記プロセス位相遅れ推定部は、制御対象プロセスモデルの操作量近傍における線形近似モデルを線形伝達関数で表現した線形伝達関数モデルで表現し、線形伝達関数モデルの位相計算により前記位相遅れの推定値を算出する、請求項４記載の最適制御装置。
【請求項8】
前記プロセス位相遅れ推定部は、前記操作量に対して前記線形伝達関数モデルを通した位相遅れ操作量を前記位相遅れの推定値として出力し、
前記評価関数勾配推定部は、前記位相遅れ操作量の値と前記評価値との、所定の期間の時系列データを用いて、前記位相遅れ操作量を入力とし、前記評価値を出力とする直線を回帰によって算出し、前記直線の１次の項の係数を前記変化率の推定値とする、請求項２又は請求項７記載の最適制御装置。
前記評価関数勾配推定部は、前記位相遅れ操作量の値と前記評価値との、所定の期間の時系列データを用いて、前記位相遅れ操作量を入力とし、前記評価値を出力とする直線を回帰によって算出し、前記直線の１次の項の係数を前記変化率の推定値とする、請求項２又は請求項７記載の最適制御装置。
【請求項9】
前記プロセス位相遅れ推定部は、前記操作量に対して前記線形伝達関数モデルを通した位相遅れ操作量を前記位相遅れの推定値として出力し、
前記評価関数勾配推定部は、前記位相遅れ操作量の値と前記評価値との相関係数を前記変化率の推定値として算出する、請求項２又は請求項７記載の最適制御装置。
前記評価関数勾配推定部は、前記位相遅れ操作量の値と前記評価値との相関係数を前記変化率の推定値として算出する、請求項２又は請求項７記載の最適制御装置。
【請求項10】
制御対象プロセスの操作量と、前記操作量に応じて変化する制御量に基づく評価関数の評価値とに基づいて、前記操作量をリアルタイムに操作して、前記評価値の最適値を探索する最適制御方法であって、
前記制御対象プロセスにおいて取得された計測値を用いて前記評価関数の前記評価値を算出し、
前記操作量および前記評価値を用いて、前記操作量から前記評価値までの位相遅れの推定値を算出し、
前記位相遅れの推定値の情報と前記評価値の情報とを用いて、前記評価値の前記操作量に対する変化率の推定値を算出し、
前記変化率の推定値を積分することにより、前記操作量の動くべき方向と量とを決定し、
決定された前記操作量の動くべき方向と量との情報に基づく前記操作量を、前記制御対象プロセスへ出力する最適制御方法。
前記制御対象プロセスにおいて取得された計測値を用いて前記評価関数の前記評価値を算出し、
前記操作量および前記評価値を用いて、前記操作量から前記評価値までの位相遅れの推定値を算出し、
前記位相遅れの推定値の情報と前記評価値の情報とを用いて、前記評価値の前記操作量に対する変化率の推定値を算出し、
前記変化率の推定値を積分することにより、前記操作量の動くべき方向と量とを決定し、
決定された前記操作量の動くべき方向と量との情報に基づく前記操作量を、前記制御対象プロセスへ出力する最適制御方法。
【請求項11】
コンピュータに、請求項１０に記載の方法を実行させるコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明の実施形態は、最適制御装置、最適制御方法およびコンピュータプログラムに関する。
【背景技術】
【0002】
水処理プラントや環境プラントなどのプラント制御では、外乱や環境条件によってプロセスの状態が変化し、外乱の影響を緩和する様にプロセスの操作量を最適化したり、環境条件に応じて、その環境下で最適なプロセスの操作量を決定したりしなければならない場合が多い。
【0003】
たとえば、下水処理プロセスや排水処理プロセスでは、下水処理場に流入する流入量や流入負荷（＝流入量×流入水質濃度）等の外乱に応じて、ブロワ風量、ポンプ流量、薬品投入量を調整し、放流水質を維持する必要がある。さらに、生物学的排水処理などでは、特に季節に応じた年間の水温変動によって最適な操作量が変化する。
また、上水処理における凝集剤注入制御や塩素注入制御では、取水源の濁度や日射量などの外乱によって、必要な水質を維持するための薬品注入量が変化する。また、凝集剤や塩素の注入量も年間の水温変動の影響を大きく受けることが多い。
また、上水処理における凝集剤注入制御や塩素注入制御では、取水源の濁度や日射量などの外乱によって、必要な水質を維持するための薬品注入量が変化する。また、凝集剤や塩素の注入量も年間の水温変動の影響を大きく受けることが多い。
【0004】
また、火力発電プラントなどにおける脱硝制御では、発電プラントから排出される窒素酸化物（ＮＯｘ）を外乱として、ＮＯｘの排出量に応じてＮＯｘを処理するためのアンモニア注入量を調整する必要がある。
石油化学プラントや鉄鋼プラントでは、原料（重油や鉄鉱石など）の質や量に応じて石油精製や製鉄に必要な操作量（温度、圧力、流量など）を最適値に制御する必要がある。
石油化学プラントや鉄鋼プラントでは、原料（重油や鉄鉱石など）の質や量に応じて石油精製や製鉄に必要な操作量（温度、圧力、流量など）を最適値に制御する必要がある。
【0005】
上記のように、多くの産業におけるプラント制御では、外乱の影響や環境条件によってプラントの状態が変化する。これらのプラントでは、外乱に対処するために各種の操作量を変化させてプラントを最適な状態に保つ様な外乱抑制型の制御や、環境条件に応じて制御を切り替えるスケジューリング型の制御が行われている場合が多く、また、環境条件に適応していく適応型の制御が採用されている場合もある。
【0006】
プラントを最適な状態に保つ代表的な制御手法として、極値制御(ＥＳＣ：Extremum Seeking Control)と呼ばれる技術が注目を浴び、産業界でも利用され始めている。極値制御は、例えば、プラントの運転コストや発電量など、最小化や最大化を図りたい最適性の評価関数を設定し、その評価値の極値（＝局所最小値あるいは局所最大値）をリアルタイムに計測しながら、操作量を常時変化させることによって、極値を探索する制御手法である。極値制御は、数式によるモデル化の難しい複雑な現象を伴うプロセス（例：下水処理プロセス、燃焼プロセス、石油化学プロセスなど）に対する実用的な最適制御技術として注目を集めている。
【0007】
極値制御は、最適化したい評価関数の値を直接計測できるオンラインセンサーにより取得される情報から計算し、評価関数値を最適値（最小値もしくは最大値）に維持する様に操作量を変化させながら適応的に探索するものである。すなわち、極値制御は、複雑な数式モデルを使わずに操作量を変化させながら最適値を探索するため、ＰＬＣ（Programmable Logic Controller）などのコントローラへの実装が容易であるという点で魅力的な方法である。
【0008】
一方で、極値制御は、安定性（極値の探索が可能なこと）を維持しながら、制御性能（極値への収束の速さや最終値と極値の誤差など）を向上させることは一般には難しく、安定性と制御性能（収束速度＋最適性（最適値との誤差））とを両立させることは、極値制御を現実の制御対象（プラント）に適用していく上で欠かせない重要な課題である。
【先行技術文献】
【特許文献】
【0009】
【特許文献1】国際公開第2020/241657号
【特許文献2】特開2019－83030号公報
【非特許文献】
【0010】
【非特許文献1】Yamanaka, Osamu, et al. "Extremum Seeking Based on Approximated Sign of Gradient of Unknown Plant Maps." 2020 59th Annual Conference of the Society of Instrument and Control Engineers of Japan (SICE)., 2020.
【非特許文献2】B. G. B. Hunnekens, M. A. M. Haring, N. van de Wouw, andH. Nijmeijer, “A dither-free Extremum-seeking control approach using1st-order least-squares fits for gradient estimation,” in 53rd Conference on Decision and Control, 2014, pp. 2679-2684.
【非特許文献3】Zengin, Nursefa, and Baris Fidan. "Adaptive extremum seeking using recursive least squares." arXiv preprint arXiv:2003.03891 (2020).
【非特許文献4】Onishi, Yuta, et al. "Extremum Seeking Control for Wastewater Treatment Plant with Prioritized Output Constraints." 2020 59th Annual Conference of the Society of Instrument and Control Engineers of Japan (SICE)., 2020.
【非特許文献5】Y. Tan, D. Nesic and I.M.Y. Mareels, “On non-local stability properties of extremum seeking control”, Automatica, vol. 42 (2006), pp. 889-903.
【発明の概要】
【発明が解決しようとする課題】
【0011】
本発明の実施形態は、上記の事情に鑑みて為されたものであり、安定性を維持しながら、制御性能を極力高める最適値の探索を実現する最適制御装置、最適制御方法、および、コンピュータプログラムを提供することを目的とする。
【課題を解決するための手段】
【0012】
実施形態による最適制御装置は、制御対象プロセスの操作量と、前記操作量に応じて変化する制御量に基づく評価関数の評価値とに基づいて、前記操作量をリアルタイムに操作して、前記評価値の最適値を探索する最適制御装置であって、前記制御対象プロセスにおいて取得された計測値を用いて前記評価関数の前記評価値を算出するプロセス評価値算出部と、前記操作量および前記評価値を用いて、前記操作量から前記評価値までの位相遅れの推定値を算出するプロセス位相遅れ推定部と、前記位相遅れの推定値に関する情報と前記評価値の情報とを用いて、前記評価値の前記操作量に対する変化率の推定値を算出する評価関数勾配推定部と、前記変化率の推定値を積分することにより、前記操作量の動くべき方向と量とを決める極値探索部と、前記極値探索部で決定された前記操作量の動くべき方向と量との情報に基づく前記操作量を、前記制御対象プロセスへ出力する操作量出力部と、を備える。
【図面の簡単な説明】
【0013】
【発明を実施するための形態】
【0014】
以下、実施形態の最適制御装置、最適制御方法、および、コンピュータプログラムについて、図面を参照して詳細に説明する。
図１は、一実施形態の最適制御装置、最適制御方法、および、コンピュータプログラムの制御対象の一例を説明するための図である。
本実施形態の最適制御装置は、制御対象プロセスの操作量と、操作量に応じて変化する制御量に基づく評価関数の評価値とに基づいて、操作量をリアルタイムに操作して、評価値の最適値を探索する装置である。
ここでは、水処理プラント（例えば浄水場）の凝集剤（ＰＡＣ）注入プロセスを、本実施形態の最適制御装置の制御対象の一例として説明する。なお、実施形態の最適制御装置の制御対象は任意のプロセスでよく、図１に示すプロセスに限定されない。
図１は、一実施形態の最適制御装置、最適制御方法、および、コンピュータプログラムの制御対象の一例を説明するための図である。
本実施形態の最適制御装置は、制御対象プロセスの操作量と、操作量に応じて変化する制御量に基づく評価関数の評価値とに基づいて、操作量をリアルタイムに操作して、評価値の最適値を探索する装置である。
ここでは、水処理プラント（例えば浄水場）の凝集剤（ＰＡＣ）注入プロセスを、本実施形態の最適制御装置の制御対象の一例として説明する。なお、実施形態の最適制御装置の制御対象は任意のプロセスでよく、図１に示すプロセスに限定されない。
【0015】
水処理プラントは、例えば、混和池と、フロック形成池と、沈殿池と、排泥池と、砂ろ過池と、浄水池とを含む。
水処理プラントに流入した原水には混和池で凝集剤が注入され、フロック形成池で撹拌されてフロックと呼ばれる固形物が形成される。フロックを含む処理水の一部は排泥池に排出され、一部は沈澱池に沈殿する。フロック形成池から排出された処理水は砂ろ過池へ送水され、ろ過された処理水が浄水池へ排出される。
水処理プラントに流入した原水には混和池で凝集剤が注入され、フロック形成池で撹拌されてフロックと呼ばれる固形物が形成される。フロックを含む処理水の一部は排泥池に排出され、一部は沈澱池に沈殿する。フロック形成池から排出された処理水は砂ろ過池へ送水され、ろ過された処理水が浄水池へ排出される。
【0016】
上記水処理プラントの凝集剤注入プロセスにおいて、凝集剤（ＰＡＣ）の注入によってフロックが形成される。形成されたフロックは、電界をかけることにより移動する。本実施形態の最適制御装置は、フロックの移動速度（ＰＶ）を画像処理によって計測しながら、移動速度がゼロ付近の目標値（ＳＶ）に追従する様にＰＩ制御で注入率（ＭＶ）を調整している。
【0017】
なお、移動速度ＰＶの目標値ＳＶは、ゼロ付近にすることが良いが、その最適値は、水温や流入水質などによって変動すると考えられる。このため、例えば本実施形態の最適制御装置が凝集剤注入プロセスを制御対象とするときには、フロックの移動速度の目標値ＳＶを極値制御により調整する。
図１に示す制御系は、カスケード制御と呼ばれる構成になっており、ＰＩ制御の目標値ＳＶが極値制御の操作量ＭＶ値になる様に、２段のカスケード構成となっている。
図１に示す制御系は、カスケード制御と呼ばれる構成になっており、ＰＩ制御の目標値ＳＶが極値制御の操作量ＭＶ値になる様に、２段のカスケード構成となっている。
【0018】
このとき、本実施形態の最適制御装置における評価関数は、例えば、ＰＡＣの薬品コスト、汚泥処分コスト、および、洗浄コストを考慮してその総和を評価関数（運用コスト）とし、沈殿水とろ過水との濁度に制約条件を設けている。制約条件は、例えば非特許文献５に示されている方法や、ペナルティ関数と呼ばれる方法により、評価関数に換算することができる。以下では簡単のため、制約条件はペナルティ関数で評価関数に組み込まれているものとし、上記運用コストにこの水質制約（制約条件）を組み込んだ評価関数を、総コストと呼ぶ。
【0019】
図２は、極値制御における極値（局所最適値）探索の原理を説明するための図である。
図２において、横軸は極値制御によって操作を行う操作量を意味し、本実施形態では移動速度の目標値（ＳＶ）である。また、図２において、縦軸は極値制御により最適化（最小化）したい評価関数（評価量）であり、本実施形態では総コストである。
図２において、横軸は極値制御によって操作を行う操作量を意味し、本実施形態では移動速度の目標値（ＳＶ）である。また、図２において、縦軸は極値制御により最適化（最小化）したい評価関数（評価量）であり、本実施形態では総コストである。
【0020】
操作量と評価関数との間には何等かの関係性があり、極値（局所最小値）を持つことが予め仮定されている。図２では、操作量と評価関数との間に、下に凸形状の関数関係があることを想定している。
【0021】
実際に極値制御を行っている際は、制御実施時のその時々の操作量に対する評価量の値はリアルタイムで取得できるが、評価関数の全体形状はリアルタイムでは把握できず、未知の状態である。このような状況において、極値制御は、図中の極値（局所最適値、図中の例は単峰性の下に凸な関数であるので大域的最小値）を探索するための制御アルゴリズムを提供するものである。
【0022】
図２によれば、操作量が正弦波などの周期的なディザー信号で駆動されている時、極値の右側（操作量が増加する側）に操作量がある場合には、ディザー信号で駆動された操作量の動きとリアルタイムで取得する評価量（評価値）の動きとは同期して動く、すなわち、同位相で動く。
【0023】
一方、操作量が正弦波などの周期的なディザー信号で駆動されている時、極値の左側（操作量が減少する側）に操作量がある場合には、ディザー信号で駆動された操作量の動きとリアルタイムで取得する評価量（評価値）の動きとは逆になる。すなわち、極値の左側では、ディザー信号で駆動された操作量とリアルタイムで取得する評価量（評価値）とが逆位相で動く。
【0024】
最適制御装置は、この情報を利用することで、制御動作中の操作量が極値に対して左右どちら側にあるかを判断できる。最適制御装置は、制御動作中の操作量が極値に対して右にある場合には操作量を減少させ、制御動作中の操作量が極値に対して左にある場合は操作量を増加させることで、極値探索が可能になる。
【0025】
なお、ディザー信号駆動型の極値制御アルゴリズムの数式を用いて分析すると、最適制御装置は、実際には、例えば図２に示す点Ａのような動作点における平均的な勾配の推定を行っており、この勾配の情報を用いて、制御動作中の操作量を減少させるか増加させるか（操作量の動くべき向き）と、操作量を減少又は増加させる大きさ（操作量の動くべき量）とを判断するものになっている。
本実施形態の最適制御装置は、上記極値制御の安定性および制御性能を改善するものである。
本実施形態の最適制御装置は、上記極値制御の安定性および制御性能を改善するものである。
【0026】
図３は、第１実施形態の最適制御装置が適用されたリアルタイムプロセス最適制御システムの一構成例を概略的に示すブロック図である。
図３に示すリアルタイムプロセス最適制御システムは、最適制御装置と、制御対象２００と、を備えている。
図３に示すリアルタイムプロセス最適制御システムは、最適制御装置と、制御対象２００と、を備えている。
【0027】
制御対象２００は、操作量入力Ｕとプロセス出力Ｙとを持つ任意のプラントにおける任意のプロセスである。本実施形態では、制御対象２００は、例えば図１に示す水処理プラントの凝集剤注入プロセスであり、操作量入力Ｕは凝集剤注入プロセスにおける移動速度目標値（ＳＶ）である。凝集剤注入プロセスにおける実際の操作量はＰＡＣ注入率であるが、先に述べた通り、この制御系は２段のカスケード構成となっているため、極値制御によって調整する操作量は移動速度の目標値である。
また、水処理プラントは、出力Ｙの値を検出する種々のセンサを備えている。凝集剤注入プロセスの出力Ｙは、例えば、凝集剤（ＰＡＣ）の注入率、排泥池の汚泥濃度、汚泥流量、沈砂池の洗浄頻度、沈殿池出口の汚泥濃度、および、浄水池出口の汚泥濃度、である。
また、水処理プラントは、出力Ｙの値を検出する種々のセンサを備えている。凝集剤注入プロセスの出力Ｙは、例えば、凝集剤（ＰＡＣ）の注入率、排泥池の汚泥濃度、汚泥流量、沈砂池の洗浄頻度、沈殿池出口の汚泥濃度、および、浄水池出口の汚泥濃度、である。
【0028】
本実施形態の最適制御装置は、例えば、少なくとも１つのプロセッサと、プロセッサにより実行されるプログラムが記録されたメモリとを備えた演算装置である。なお、最適制御装置のメモリには、プログラムを実行中に一時的にデータ（例えば計測値や評価値の時系列データ）を記憶するための記憶領域（記憶部）が設けられ得る。最適制御装置は、ソフトウエアにより、又は、ソフトウエアとハードウエアとの組み合わせにより、以下に説明する種々の機能を実現することができる。
【0029】
最適制御装置は、プロセス位相遅れ推定部４００と、プロセス計測値取得部１１０と、プロセス評価値算出部１２０と、復調用ディザー信号生成部１３０と、評価関数勾配推定部１４０と、正規化信号発生部１５０と、勾配推定量正規化部１６０と、極値探索部（最適操作量適応調整部）１７０と、変調用ディザー信号生成部１８０と、操作量出力部１９０と、を備えている。プロセス位相遅れ推定部４００は、位相遅れパラメータ推定部３００と、位相遅れ推定部１１００とを含む。
【0030】
プロセス計測値取得部１１０は、評価量や制約条件を算出するために必要となる計測情報を、制御対象２００から取得する。
一例として、薬品費と汚泥処分費とろ過池洗浄費とを合計した運用コストを評価関数として設定する場合には、プロセス計測値取得部１１０は、例えば、ＰＡＣの注入率、排泥池の汚泥濃度、汚泥流量、沈砂池の洗浄頻度、などを計測情報として取得する。また、沈殿池出口や浄水池出口の汚泥濃度に制約を設ける場合は、プロセス計測値取得部１１０は、制約を設けた汚泥濃度の値も所定周期で計測し、所定のフォーマットで時系列データとして保存しておく。また、プロセス計測値取得部１１０は、ろ過池の洗浄履歴データを所定のフォーマットで取得し、保存可能である。
一例として、薬品費と汚泥処分費とろ過池洗浄費とを合計した運用コストを評価関数として設定する場合には、プロセス計測値取得部１１０は、例えば、ＰＡＣの注入率、排泥池の汚泥濃度、汚泥流量、沈砂池の洗浄頻度、などを計測情報として取得する。また、沈殿池出口や浄水池出口の汚泥濃度に制約を設ける場合は、プロセス計測値取得部１１０は、制約を設けた汚泥濃度の値も所定周期で計測し、所定のフォーマットで時系列データとして保存しておく。また、プロセス計測値取得部１１０は、ろ過池の洗浄履歴データを所定のフォーマットで取得し、保存可能である。
【0031】
プロセス評価値算出部１２０は、プロセス計測値取得部１１０で取得した情報（制御対象プロセスの計測値）を用いて、予め設定した評価関数の評価値をリアルタイムに計算する。
運用コストは、薬品費（ＰＡＣ費）と汚泥処分費とろ過池洗浄費との総和で定義されている。
薬品費は、ＰＡＣの注入量の時系列データに薬品単価や希釈率などの係数をかけることで算出でき、時系列データとして得ることができる。
汚泥処分費は、排泥池の汚泥濃度と汚泥流量との積で発生汚泥量を算出し、発生汚泥量に汚泥処分単価をかけることで算出でき、時系列データを得ることができる。
運用コストは、薬品費（ＰＡＣ費）と汚泥処分費とろ過池洗浄費との総和で定義されている。
薬品費は、ＰＡＣの注入量の時系列データに薬品単価や希釈率などの係数をかけることで算出でき、時系列データとして得ることができる。
汚泥処分費は、排泥池の汚泥濃度と汚泥流量との積で発生汚泥量を算出し、発生汚泥量に汚泥処分単価をかけることで算出でき、時系列データを得ることができる。
【0032】
また、ろ過池の洗浄は、通常ろ抗（ろ過抵抗）が所定のしきい値を超過すると行われるため、ろ過抵抗の時系列データや洗浄履歴データから洗浄が行われたタイミングを知ることができ、洗浄が行われたタイミングの洗浄費は履歴データから得られる。洗浄は前回の洗浄実施時からの次に洗浄を行うまでの期間の費用と考えられるので、この洗浄費を当該期間に均等に分配することで洗浄費を時系列データに換算することができる。
【0033】
ただし、洗浄費の時系列データへの換算は、過去の実績データに対しては容易に実施できるが、リアルタイムでは、次回洗浄を行うタイミングは未確定であるため、ろ抗の変化率を監視しながら、次回の洗浄タイミングを推定し、そこから、リアルタイムでの洗浄費の時系列データに換算してもよい。
【0034】
運用コストは、時系列データに換算された各費用の総和である。したがって、プロセス評価値算出部１２０は、上記の方法により得られた各費用を用いて、リアルタイムで運用コストを算出することが可能である。
一方、水質制約は、図１の沈殿池出口および浄水池出口の濁度を時系列データとして取得し、これに、濁度の制約値、例えば、Ｔｌｉｍ（濁度上限値）＝０．８度以下などの制約条件を組み込む。
一方、水質制約は、図１の沈殿池出口および浄水池出口の濁度を時系列データとして取得し、これに、濁度の制約値、例えば、Ｔｌｉｍ（濁度上限値）＝０．８度以下などの制約条件を組み込む。
【0035】
極値制御では、制約条件を直接扱うことはできないため、例えば、非特許文献５の方法などで評価関数として扱える様に変換することができる。ここでは、最適化分野で良く知られたペナルティ関数の考え方で評価関数として扱う。すなわち、例えば、以下の（１）式の様に水質コストを表す。
【0036】
Ｗｃｏｓｔ＝ｍａｘ（０，ａ×（ｅｘｐ（Ｔ－Ｔｌｉｍ）－１））……………（１）
ここで、Ｔは濁度計測値、Ｔｌｉｍは濁度上限値、ａは設計パラメータであり、パラメータａはゼロより大きい（ａ＞０）。
コストＷｃｏｓｔは換算された評価関数である。本実施形態の最適制御装置では、沈殿池出口の濁度と浄水池出口の濁度とについて（１）式の変換を行い、沈殿池出口の濁度についてのコストＷｃｏｓｔと、浄水池出口の濁度についてのコストＷｃｏｓｔとの和を水質コストとする。
ここで、Ｔは濁度計測値、Ｔｌｉｍは濁度上限値、ａは設計パラメータであり、パラメータａはゼロより大きい（ａ＞０）。
コストＷｃｏｓｔは換算された評価関数である。本実施形態の最適制御装置では、沈殿池出口の濁度と浄水池出口の濁度とについて（１）式の変換を行い、沈殿池出口の濁度についてのコストＷｃｏｓｔと、浄水池出口の濁度についてのコストＷｃｏｓｔとの和を水質コストとする。
【0037】
（１）式は、濁度計測値Ｔが濁度上限値Ｔｌｉｍ以下（濁度計測値Ｔ≦濁度上限値Ｔｌｉｍ）のときにはコストＷｃｏｓｔは０となり、濁度計測値Ｔが濁度上限値Ｔｌｉｍを超過（濁度計測値Ｔ＞濁度上限値Ｔｌｉｍ）するとＷｃｏｓｔが急激（指数関数的）に上昇することを意味しており、いわゆるペナルティ関数の一種である。
【0038】
先に述べた運用コストに上記水質コストを加えて、総コストＪを以下の（２）式に示す様に定義することができる。
総コストＪ＝薬品コスト＋汚泥処分コスト＋ろ過池洗浄コスト＋水質コスト…（２）
なお、本実施形態の最適制御装置における凝集剤注入制御では、上述のような評価関数を設定したが、問題設定によっては、対象プラント（制御対象）２００とプロセス評価値算出部１２０とを分離できない場合がある。
総コストＪ＝薬品コスト＋汚泥処分コスト＋ろ過池洗浄コスト＋水質コスト…（２）
なお、本実施形態の最適制御装置における凝集剤注入制御では、上述のような評価関数を設定したが、問題設定によっては、対象プラント（制御対象）２００とプロセス評価値算出部１２０とを分離できない場合がある。
【0039】
例えば、風力発電プラントにおいて、風車のブレードの向きを風向に併せて動かして発電量を最大化する制御に極値制御を適用する場合、評価関数Ｊは発電量であり、操作量Ｕは風車のブレードの回転角である。この例では、プラントの出力Ｙを特に定義する必要は無い。このような問題設定においては、出力Ｙと評価関数Ｊとを区別しないケースもあり、必ずしも出力Ｙが評価関数Ｊと分離した状態で定義されて計測されているとは限らない。
【0040】
一方、対象プラント（制御対象）２００とプロセス評価値算出部１２０とを分離できない場合であっても、本実施形態の最適制御装置と同様に、最適化評価関数を個別に定義することによって、極値制御を適用できる対象になる場合がある。
いずれの場合も評価関数が適切に設定されると、評価関数の値をリアルタイムに所定の制御周期で計測・算出することにより、時々刻々と変化する評価量を取得することができる。
いずれの場合も評価関数が適切に設定されると、評価関数の値をリアルタイムに所定の制御周期で計測・算出することにより、時々刻々と変化する評価量を取得することができる。
【0041】
なお、後述する位相遅れパラメータ推定部３００で導かれる下記（６）式や（７）式の伝達関数モデルの出力Ｙ（ｓ）に対応する信号としてプロセス評価値を用いる場合には、パラメータの推定値を演算するよりも前に、プロセス評価値算出部１２０によりプロセス評価値を取得しておく必要がある。この場合には、プロセス評価値算出部１２０は、位相遅れパラメータ推定部３００にプロセス評価値を供給する。
【0042】
評価関数勾配推定部１４０は、プロセス評価値の操作量に対する勾配（変化率）の推定値を演算する。
図４は、第１実施形態の最適制御装置の一部の構成例を概略的に示す図である。なお、図４では以下の説明に必要な構成を示し、その他の構成の詳細は省略されている。
図４は、第１実施形態の最適制御装置の一部の構成例を概略的に示す図である。なお、図４では以下の説明に必要な構成を示し、その他の構成の詳細は省略されている。
【0043】
評価関数勾配推定部１４０は、ハイパスフィルタ（ＨＰＦ）１４１と、ローパスフィルタ（ＬＰＦ）１４３と、加算器１４２と、を含む。
ハイパスフィルタ（ＨＰＦ）１４１には、プロセス評価値算出部１２０から出力された評価値（時系列データ）が入力される。ハイパスフィルタ１４１は、入力された評価値の所定の周波数以上の成分を出力する。ハイパスフィルタ１４１は、未知の極値（局所最適値）が変化しない、あるいは、変化したとしても非常に緩やかに変化する、と仮定した場合に、極値を強制的に（近似的に）０（ゼロ）にするために導入される。評価関数勾配推定部１４０がハイパスフィルタ１４１を備えることにより制御性能を向上させることができる。なお、ハイパスフィルタ１４１は、極値制御の理論的な構成上は必須のものではなく省略されても構わない。
ハイパスフィルタ（ＨＰＦ）１４１には、プロセス評価値算出部１２０から出力された評価値（時系列データ）が入力される。ハイパスフィルタ１４１は、入力された評価値の所定の周波数以上の成分を出力する。ハイパスフィルタ１４１は、未知の極値（局所最適値）が変化しない、あるいは、変化したとしても非常に緩やかに変化する、と仮定した場合に、極値を強制的に（近似的に）０（ゼロ）にするために導入される。評価関数勾配推定部１４０がハイパスフィルタ１４１を備えることにより制御性能を向上させることができる。なお、ハイパスフィルタ１４１は、極値制御の理論的な構成上は必須のものではなく省略されても構わない。
【0044】
加算器１４２には、ハイパスフィルタ１４１から出力された値と、後述する復調用ディザー信号（位相遅れの推定値に関する情報）とが入力される。加算器１４２は、入力された値を加算した和をローパスフィルタ１４３へ出力する。
ローパスフィルタ１４３は、加算器１４２から供給された信号の所定の周波数よりも小さい成分を出力する。ローパスフィルタ１４３から出力される信号の周期的な平均値は、評価関数の勾配に比例する信号となることが理論的な解析結果から知られているため、ローパスフィルタ１４３の出力信号は、評価関数値の勾配（正確には勾配に比例する信号）と見なすことができる。例えば、勾配推定値は勾配の方向（符号（１又は－１））と大きさとの情報を含む。
ローパスフィルタ１４３は、加算器１４２から供給された信号の所定の周波数よりも小さい成分を出力する。ローパスフィルタ１４３から出力される信号の周期的な平均値は、評価関数の勾配に比例する信号となることが理論的な解析結果から知られているため、ローパスフィルタ１４３の出力信号は、評価関数値の勾配（正確には勾配に比例する信号）と見なすことができる。例えば、勾配推定値は勾配の方向（符号（１又は－１））と大きさとの情報を含む。
【0045】
正規化信号発生部１５０は、正規化信号を発生させる。正規化信号発生部１５０の作用は、例えば特許文献１に記載された正則化信号発生部の作用と同一であり、様々な方法をとることができる。本実施形態の最適制御装置では、評価関数の勾配推定値Ｇ（ｔ）を、以下の性質［１］－［４］を満たす信号Ｇｎ（ｔ）に変換するための信号を「正規化信号」とする。
【0046】
［１］Ｇ（ｔ）＝０⇔Ｇｎ（ｔ）＝０（Ｇ（ｔ）が０の時に限りＧｎ（ｔ）も０となる）
［２］Ｇ（ｔ）が正（負）⇔Ｇｎ（ｔ）が正（負）（Ｇ（ｔ）とＣｎ（ｔ）との符号は同じ）
［３］Ｇ（ｔ）＜∞⇔Ｇｎ（ｔ）＜∞（Ｇ（ｔ）が有限の時はＧｎ（ｔ）も有限に留まり、ゼロ割などが起こらない）
［４］Ｇｎ（ｔ）は、Ｇ（ｔ）→∞の時、Ｇｎ（ｔ）→∞とはならずある正の有限値に近づく、すなわち、ある０＜ｋ＜∞が存在して、Ｇｎ（ｔ）→ｋとなる。
［２］Ｇ（ｔ）が正（負）⇔Ｇｎ（ｔ）が正（負）（Ｇ（ｔ）とＣｎ（ｔ）との符号は同じ）
［３］Ｇ（ｔ）＜∞⇔Ｇｎ（ｔ）＜∞（Ｇ（ｔ）が有限の時はＧｎ（ｔ）も有限に留まり、ゼロ割などが起こらない）
［４］Ｇｎ（ｔ）は、Ｇ（ｔ）→∞の時、Ｇｎ（ｔ）→∞とはならずある正の有限値に近づく、すなわち、ある０＜ｋ＜∞が存在して、Ｇｎ（ｔ）→ｋとなる。
【0047】
正規化信号発生部１５０は、既知の信号発生方法を採用し、上記の性質を満たす任意の信号を正規化信号として採用することができる。
例えば、評価関数の勾配を極値制御に用いるのではなく、勾配の方向（１又は－１）を極値制御に用いる方が、極値制御の制御性能を高める（速く収束させる）ための各種パラメータの調整が容易になる。
例えば、評価関数の勾配を極値制御に用いるのではなく、勾配の方向（１又は－１）を極値制御に用いる方が、極値制御の制御性能を高める（速く収束させる）ための各種パラメータの調整が容易になる。
【0048】
操作量と評価量の関係を表す評価関数形状は未知であり、実際には例えば図２に示すような下に凸の関数になっているとは限らず、動作点によって勾配が極めて急峻になったり、極めて緩やかになったりしている可能性がある。制御実行時には、評価関数の形状の情報を予め知ることができないため、極値制御では評価関数の勾配をリアルタイムで推定している。
【0049】
例えば、推定された勾配が緩やかな場合は、制御を強く効かせることができるが、推定された勾配が急峻な場合には、制御を強く効かせると安定性を損なうリスクがある。したがって、推定された勾配の情報に基づいて極値制御を行うときには、想定される勾配の中で最も急峻と思われる勾配を推測・予想して、極値制御のパラメータの調整を行う必要が出てくる。このとき、制御パラメータの設定は安定性を損なわない様に保守的な調整をせざるを得なくなり、結果として制御性能を向上させる（収束速度を上げる）ことが困難になる。
【0050】
しかしながら、極値探索を行うためには、必ずしも勾配そのものを推定して知る必要は無く、勾配の方向、すなわち、符号（操作量を上げるべきか下げるべきか）の情報だけ得ることができれば、原理的には極値探索を行えるはずである。この観点でみると、正規化信号として次式を用いることが最も適切であると考えられる。
ＲＳ＝１／｜Ｇ（ｔ）｜……………………………………………………………（３）
ここで、ＲＳは正規化信号を意味し、｜Ｇ（ｔ）｜は評価関数の勾配推定値の絶対値である。これを正規化信号として用いると、Ｇ（ｔ）を正規化した信号Ｇｎ（ｔ）は勾配の符号関数となる。
ＲＳ＝１／｜Ｇ（ｔ）｜……………………………………………………………（３）
ここで、ＲＳは正規化信号を意味し、｜Ｇ（ｔ）｜は評価関数の勾配推定値の絶対値である。これを正規化信号として用いると、Ｇ（ｔ）を正規化した信号Ｇｎ（ｔ）は勾配の符号関数となる。
【0051】
ただし、この関数をそのまま用いると、Ｇ（ｔ）とＧｎ（ｔ）との関係を表す平面の原点で不連続となり、上記［１］の条件を満たさなくなる。そのため、上記（３）式は、若干滑らかにする方が好ましい。例えば、次式の様な正則化定数δ＞０を導入することで、以下の様な正規化信号が得られる。
ＲＳ＝１／（δ＋｜Ｇ（ｔ）｜）…………………………………………………（４）
正規化信号発生部１５０は、上記（４）式により得られる正規化信号を勾配推定量正規化部１６０へ出力する。
ＲＳ＝１／（δ＋｜Ｇ（ｔ）｜）…………………………………………………（４）
正規化信号発生部１５０は、上記（４）式により得られる正規化信号を勾配推定量正規化部１６０へ出力する。
【0052】
勾配推定量正規化部１６０は、正規化信号発生部１５０で発生させた正規化信号を用いて、勾配推定量を正規化する。すなわち、勾配推定量正規化部１６０は、評価関数勾配推定部１４０で推定した勾配推定値Ｇ（ｔ）に対して、例えば（４）式の正規化信号を作用させて正規化後の勾配推定値Ｇｎ（ｔ）を生成する。正規化された勾配推定値は勾配の方向（符号）の情報を含む。勾配推定量正規化部１６０は、正規化後の勾配推定値Ｇｎ（ｔ）を極値探索部１７０へ出力する。
【0053】
極値探索部１７０は、積分器１７１と、ゲイン乗算部１７２と、を備え、勾配（変化率の推定値）の情報を積分することにより、操作量の動くべき方向と量とを決める。
積分器１７１は、勾配推定値あるいは正規化された勾配推定値を積分した値を出力する。
ゲイン乗算部１７２は、積分器１７１から出力された値に積分ゲインＫを乗じて出力する。積分ゲインＫは、極値探索の収束速度を決める重要なパラメータであるが、後述する位相遅れ補償を導入することにより、従来の極値制御よりもその値を大きくすることができ、結果的に制御性能を高めることができる。
積分器１７１は、勾配推定値あるいは正規化された勾配推定値を積分した値を出力する。
ゲイン乗算部１７２は、積分器１７１から出力された値に積分ゲインＫを乗じて出力する。積分ゲインＫは、極値探索の収束速度を決める重要なパラメータであるが、後述する位相遅れ補償を導入することにより、従来の極値制御よりもその値を大きくすることができ、結果的に制御性能を高めることができる。
【0054】
なお、積分器１７１およびゲイン乗算部１７２を用いた制御アルゴリズムは、勾配法と呼ばれる最適化アルゴリズムとほぼ同じものであり、積分ゲインＫは勾配法で呼ばれる学習率と呼ばれるパラメータに相当している。
【0055】
変調用ディザー信号生成部１８０は、変調用のディザー信号（第２周期信号）を生成する。
変調用ディザー信号生成部１８０で生成される変調用ディザー信号は、任意の周期的信号で良いが、後述する復調用ディザー信号生成部１３０の信号と同じ波形の信号である必要がある。また、パラメータ同定の観点からは、ディザー信号は、周期信号であっても単一周波数の正弦波よりも複数の周波数の合成となる矩形波などのディザー信号の方が好ましい。後述する復調用ディザー信号（（９）式）に対応する正弦波信号を用いた場合には、変調用ディザー信号生成部１８０では次式の変調用ディザー信号を用いる。
Ｍ（ｔ）＝Ａｓｉｎ（ωｔ）…………………………………………………………（５）
ここで、Ａはディザー信号の振幅であり、調整可能なパラメータである。変調用ディザー信号は、ディザー信号駆動型の極値探索制御では必須の要素である。
変調用ディザー信号生成部１８０で生成される変調用ディザー信号は、任意の周期的信号で良いが、後述する復調用ディザー信号生成部１３０の信号と同じ波形の信号である必要がある。また、パラメータ同定の観点からは、ディザー信号は、周期信号であっても単一周波数の正弦波よりも複数の周波数の合成となる矩形波などのディザー信号の方が好ましい。後述する復調用ディザー信号（（９）式）に対応する正弦波信号を用いた場合には、変調用ディザー信号生成部１８０では次式の変調用ディザー信号を用いる。
Ｍ（ｔ）＝Ａｓｉｎ（ωｔ）…………………………………………………………（５）
ここで、Ａはディザー信号の振幅であり、調整可能なパラメータである。変調用ディザー信号は、ディザー信号駆動型の極値探索制御では必須の要素である。
【0056】
操作量出力部１９０は加算器を備え、極値探索部１７０の出力と変調用ディザー信号とを足し合わせた信号（操作量指令信号）を、制御対象２００に対して出力する。
本実施形態の最適制御装置において、上記のように変調用ディザー信号を足し合わせた信号を操作量指令信号とすることにより、変調用ディザー信号により操作量を周期的に駆動するとともに、極値の探索を行うことができる。
本実施形態の最適制御装置において、上記のように変調用ディザー信号を足し合わせた信号を操作量指令信号とすることにより、変調用ディザー信号により操作量を周期的に駆動するとともに、極値の探索を行うことができる。
【0057】
ここで、実際には、制御対象にダイナミクスがある場合には、操作量の時系列データと評価値の時系列データとの間には、プロセスのダイナミクスによる時間の遅れ（操作量が正弦波で駆動される場合には位相の遅れ）が必ず存在する。
【0058】
図５は、制御対象にダイナミクスがあるときの操作量の時系列データと評価量の時系列データとの関係の一例を概略的に示す図である。
制御対象にダイナミクスがある場合、例えば図５において、制御動作中の操作量が極値の右側にある場合でも、操作量と評価量とが同位相で動くわけではなく、操作量に対して位相の遅れを伴って評価量が変化する。同様に、制御動作中の操作量が極値の左側にある場合でも、操作量と評価量とが逆位相で動くわけではなく、操作量に対して位相の遅れを伴って評価量が変化する。
制御対象にダイナミクスがある場合、例えば図５において、制御動作中の操作量が極値の右側にある場合でも、操作量と評価量とが同位相で動くわけではなく、操作量に対して位相の遅れを伴って評価量が変化する。同様に、制御動作中の操作量が極値の左側にある場合でも、操作量と評価量とが逆位相で動くわけではなく、操作量に対して位相の遅れを伴って評価量が変化する。
【0059】
例えば、操作量に対する評価量の位相の遅れが９０度になると、制御動作中の操作量が極値の右側にある場合に、操作量と評価量とが逆位相になり、制御動作中の操作量が極値の左側にある場合には、操作量と評価量とが同位相になる。
上記のように操作量に対する評価量の位相遅れがある状態で、極値制御による最適化を行うと、操作量が収束するまでの時間を要する可能性があった。
上記のように操作量に対する評価量の位相遅れがある状態で、極値制御による最適化を行うと、操作量が収束するまでの時間を要する可能性があった。
【0060】
操作量に対する評価量の位相遅れを極力小さくするためには、制御対象がほとんど静的なプロセスと見なせるぐらい速いプロセスである必要がある。制御対象が十分に早く応答し、近似的に静的なプロセスと見なせるならば、操作量に対する評価量の位相遅れを小さくすることができる。一方で、実際の制御対象が静的なプロセスと見なせるほど速く応答するプロセスであるとは限らない。
【0061】
応答が速い、あるいは、応答が遅い、という概念は相対的なものである点に着目すると、ディザー信号の周波数を調整することによって、実際の制御対象を近似的に静的と見なすことが可能である。すなわち、実際の制御対象の応答速度は制御対象固有の特性であるため、その応答速度が十分に速い、と判断できるように、ディザー信号の周波数を遅く（小さく）設計する必要がある。そうすることで、制御対象による応答の遅れが所与のものであっても、それを、位相遅れで見た場合に十分小さいくすることができる。これにより、極値制御における「安定性」と「制御性能」とのトレードオフが生じる。
【0062】
そこで、本実施形態の最適制御装置では、以下のように位相遅れ補償を行っている。
位相遅れパラメータ推定部３００は、操作量（＝移動速度目標値）に対するプロセス評価値の位相遅れ情報を取得するために、位相遅れに関するパラメータの推定値を演算する。なお、位相遅れパラメータ推定部３００は、パラメータ推定をオフラインで行ってもよく、オンラインで行ってもよい。
位相遅れパラメータ推定部３００は、操作量（＝移動速度目標値）に対するプロセス評価値の位相遅れ情報を取得するために、位相遅れに関するパラメータの推定値を演算する。なお、位相遅れパラメータ推定部３００は、パラメータ推定をオフラインで行ってもよく、オンラインで行ってもよい。
【0063】
本実施形態の最適制御装置では、位相遅れパラメータ推定部３００において位相遅れを定義するために適切なモデルとして、例えば下記（６）式のような線形伝達関数モデルを想定することができる。
ここで、ｓはラプラス演算子を意味し、ａｉ,ｉ＝０,１,２,…ｎ、ｂｉ,ｉ＝０,１,２,…ｍは、各々伝達関数Ｇ（ｓ）の分母多項式と分子多項式との係数、ｅｘｐ（－Ｌｓ）は、むだ時間（遅れ時間）の伝達関数であり、Ｌは遅れ時間を表す。また、Ｕ（ｓ）は操作量をラプラス変換した値であり、Ｙ（ｓ）はプロセス計測値をラプラス変換した値である。
【数1】
【0064】
プロセス計測値は複数の項目の値を含み得るため、位相遅れパラメータ推定部３００は、（６）式に相当する線形伝達関数モデルをプロセス計測値毎に用いても良い。一般的に、制御系の応答は最も応答の遅いものが律速となるため、例えば、この場合は、Ｙ（ｓ）として浄水池濁度をラプラス変換した値などを代表値としても良い。あるいは、プロセス計測値から、後述するプロセス評価値に換算した上で、プロセス評価値のラプラス変換した値をＹ（ｓ）とすることができる。
【0065】
Ｕ（ｓ）とＹ（ｓ）とを定義することにより、（６）式の伝達関数モデルのパラメータ、ａｉ,ｉ＝０,１,２,…ｎ、ｂｉ,ｉ＝０,１,２,…ｍ、および、Ｌを同定することができる。なお、（６）式は、一般的な連続時間系の線形伝達関数モデルであるが、位相遅れパラメータ推定部３００では、プロセスの静的な情報は不要であるため、予め静的なゲインを１としておいてよい。この場合、Ｕ（ｓ）からＹ（ｓ）の（静的な）ゲインが１になる様に、予めＵ（ｓ）とＹ（ｓ）とに対応するデータを正規化しておく必要がある。予めＵ（ｓ）とＹ（ｓ）とに対応するデータを正規化することでゲインを１に規格化できる。また、（６）式でゲインを１にするために、予めａ０＝ｂ０としておく必要がある。
【0066】
位相遅れパラメータ推定部３００が、オフラインでパラメータ推定を行う場合には、予め蓄積しておいたＵ（ｓ）とＹ（ｓ）とに対応する時系列データを用いて、制御理論分野で良く知られているシステム同定法を流用することで、これらのパラメータを同定することができる。
【0067】
位相遅れパラメータ推定部３００が、オンラインでパラメータ推定を行う場合には、例えば、適応制御で知られている適応オブザーバなどの考え方を流用することで、これらのパラメータを同定することができる。例えば、制御対象２００が外乱によって大きく励起（駆動）されているときには、そのまま適応オブザーバを用いることも可能である。
【0068】
なお、システム同定法や適応オブザーバなどを適用する場合には、（６）式を直接用いるのではなく、（６）式を離散化した離散伝達関数モデルを用いる方が一般的である。また、これらの理論を適用する際には、パラメータが正しく同定できるための条件（可同定条件）が知られているので、可同定を満たす様なデータを用いる必要がある。
具体的には、オフラインの場合には、操作量（＝移動速度目標値ＳＶ）をＭ系列（ＭＬＳ）と呼ばれる信号で変化させたデータを予め取っておくことが好ましい。
上記のように、位相遅れパラメータ推定部３００は、システム同定理論や適応制御理論に基づいて、システマティックに伝達関数モデルの係数を同定することができる。
具体的には、オフラインの場合には、操作量（＝移動速度目標値ＳＶ）をＭ系列（ＭＬＳ）と呼ばれる信号で変化させたデータを予め取っておくことが好ましい。
上記のように、位相遅れパラメータ推定部３００は、システム同定理論や適応制御理論に基づいて、システマティックに伝達関数モデルの係数を同定することができる。
【0069】
なお、システム同定や適応制御分野では、同定すべきパラメータ数が増加するとパラメータを正しく同定できる可同定条件が厳しくなることも広く知られている。一般に、凝集剤注入制御などのプロセス制御分野では、あまり大きなかく乱をプラントに印加することは好まれない場合が多く、実際には、ステップ応答試験ぐらいしか実施できない場合も多い。これは、実際の制御対象を相手にする場合には無視できない強い制約となり、このような現実的な制約のため、実際の現場ではシステム同定や適応制御を用いずに、ステップ応答試験およびＰＩＤ制御の様なシンプルな制御が好まれているという側面がある。
【0070】
位相遅れパラメータ推定部３００においても、（６）式の形を一般系で用いるのは現実的には得策でない場合も多く、より実用的な方法として例えば（７）式のモデルを採用することができる。
上記（７）式のモデルは（６）式の特殊なケースでああり、１次遅れ＋むだ時間モデルである。（７）式は、遅れ時間Ｌと時定数Ｔという２つのパラメータしか含んでおらず、これらの２つのパラメータは、ステップ応答試験を実施できれば、直ちに取得することができるため、実際の取得も容易である。
【数2】
【0071】
本実施形態の最適制御装置では、操作量の時系列データに上記（６）式や（７）式の伝達関数を乗じることで、操作量の時系列データの位相を遅れさせて、操作量に位相遅れ情報を組み込んだ位相遅れ操作量の時系列データと、評価値の時系列データとの間の位相差を解消している。
【0072】
なお、（７）式は、ＰＩＤ制御の調整に用いられる典型的なモデルの特殊系である。ＰＩＤ制御の調整では、ｅｘｐ（－Ｌｓ）/（Ｔｓ＋１）ではなく、Ｋｅｘｐ（－Ｌｓ）/（Ｔｓ＋１）というモデルを用い、プロセスゲインと呼ばれるもう一つのパラメータＫも同定することが通常である。（７）式でＫを組み込んでいないのは、後述するように、位相遅れには、プラントの静的な情報は関係しないためである。
【0073】
また、位相遅れパラメータ推定部３００がオンラインで、遅れ時間Ｌと時定数Ｔとの推定値を演算する場合には、先に述べた適応オブザーバの考え方を流用しても良い。
遅れ時間Ｌと時定数Ｔとの２つのパラメータの推定値を演算するときには、位相遅れパラメータ推定部３００は、ディザー信号の情報を利用することができる。例えば、ディザー信号として正弦波ではなく矩形波を用いている場合には、矩形波の変化前後を一つのステップと考えることができるため、位相遅れパラメータ推定部３００は、ステップ応答の考え方を直接利用して遅れ時間Ｌと時定数Ｔとを同定することができる。
遅れ時間Ｌと時定数Ｔとの２つのパラメータの推定値を演算するときには、位相遅れパラメータ推定部３００は、ディザー信号の情報を利用することができる。例えば、ディザー信号として正弦波ではなく矩形波を用いている場合には、矩形波の変化前後を一つのステップと考えることができるため、位相遅れパラメータ推定部３００は、ステップ応答の考え方を直接利用して遅れ時間Ｌと時定数Ｔとを同定することができる。
【0074】
例えば、ディザー信号として正弦波を用いている場合には、可同定性が劣化するので、直接求めることが困難となる可能性もあるが、位相遅れパラメータ推定部３００がディザー信号の正弦波の周波数を２段階（若しくは２段階以上）変化させ、Ｕ（ｓ）とＹ（ｓ）とに対応する２つの時系列データの位相差を直接推定することで、遅れ時間Ｌと時定数Ｔを逆算して求めることができる。
【0075】
１次遅れむだ時間の位相は、遅れ時間Ｌと時定数Ｔとの関数により表されるため、位相情報が２つ得られれば逆算により遅れ時間Ｌと時定数Ｔとを求められる。例えば３つ以上の周波数を用いた場合は、誤差最小化基準などを設けてパラメータ最適化を行うことで、遅れ時間Ｌと時定数Ｔとを求められる。
【0076】
一方、制御対象プロセスのプロセスモデル（制御対象プロセスモデル）やプロセスシミュレータが既にある場合は、必ずしも、（６）式の様な線形伝達関数、あるいは、その特殊例である（７）式の様な１次遅れ＋むだ時間モデルを同定する必要はなく、プロセスシミュレータを直接用いることもできる。
【0077】
本実施形態の最適制御装置では、浄水場のＰＡＣ注入とＰＡＣ注入とにより移動速度や沈殿池・浄水池出口の濁度が変化する様な凝集プロセスモデル・凝集プロセスシミュレータがある場合には、プロセスシミュレータに含まれる各種パラメータを同定することが位相遅れパラメータ推定処理に相当する。
【0078】
プロセスシミュレータに含まれる各種パラメータの推定方法としては、様々な方法が考えられる。プロセスシミュレータが非線形モデルである場合には、実際のプロセスの入出力とプロセスシミュレータの入出力とが極力一致する様に、同じ入力に対する実際のプロセスの出力とシミュレータ出力との誤差に対して、例えば、２乗誤差最小化などを評価基準として、誤差を最小化する最適化問題を解いてもよく、あるいは、データ同化と呼ばれる方法を適用してもよい。このことにより、位相遅れパラメータ推定部３００は、制御対象プロセスモデルのパラメータを同定することができる。
【0079】
上記のように、位相遅れパラメータ推定部３００は、位相遅れパラメータをオフラインあるいはオンラインで同定することができる。位相遅れパラメータ推定部３００は、同定したパラメータの推定値を、位相遅れ推定部１１００に供給する。
【0080】
位相遅れ推定部１１００は、位相遅れパラメータ推定部３００から供給されたパラメータの推定値を用いて、操作量から評価量までの位相遅れを求めることができる。
本実施形態の最適制御装置において、位相遅れ推定部１１００は、復調用ディザー信号に位相遅れ情報を組み込むため、例えば（６）式あるいは（７）式の伝達関数モデルを用いて、伝達関数モデルの位相遅れを計算する。位相遅れ推定部１１００は、線形伝達関数モデルのボード線図作成時などに用いられる広く一般的に知られた位相計算方法を直接適用することにより、伝達関数モデルの位相遅れを計算できる。ここでは、簡単のため、（７）式の伝達関数モデルの場合の具体的な位相遅れの計算式を（８）式に示す。
∠Ｇ（ｓ）＝－ωＬ－ｔａｎ－１（ωＴ）………………………………………（８）
ここで、ωは周波数である。本実施形態の最適制御装置では、ディザー信号を用いて強制的に操作量を動かすため、例えばディザー信号として正弦波を用いる場合は、ディザー信号の周波数をωの値とすることができる。
本実施形態の最適制御装置において、位相遅れ推定部１１００は、復調用ディザー信号に位相遅れ情報を組み込むため、例えば（６）式あるいは（７）式の伝達関数モデルを用いて、伝達関数モデルの位相遅れを計算する。位相遅れ推定部１１００は、線形伝達関数モデルのボード線図作成時などに用いられる広く一般的に知られた位相計算方法を直接適用することにより、伝達関数モデルの位相遅れを計算できる。ここでは、簡単のため、（７）式の伝達関数モデルの場合の具体的な位相遅れの計算式を（８）式に示す。
∠Ｇ（ｓ）＝－ωＬ－ｔａｎ－１（ωＴ）………………………………………（８）
ここで、ωは周波数である。本実施形態の最適制御装置では、ディザー信号を用いて強制的に操作量を動かすため、例えばディザー信号として正弦波を用いる場合は、ディザー信号の周波数をωの値とすることができる。
【0081】
なお、位相遅れ推定部１１００において（６）式の伝達関数モデルを用いて位相遅れを演算する場合には、（７）式のモデルを用いる場合よりも複雑な計算が必要になるものの、周波数ωと線形伝達関数モデルのパラメータ値とから位相遅れを計算することが可能である。
一方で、線形伝達関数モデルで表す事が適切でない場合は、予め、プロセスシミュレータ（制御対象プロセスを模擬する物理化学式で記述されたシミュレーションモデル）を直接用いて位相推定をすることができる。この際、位相推定に極値制御を適用するためには、出力として何等かの評価関数が必要となるため、いわゆるプロセスシミュレータ（制御対象プロセスのシミュレーションモデル）と評価関数モデルとを組み合わせたシミュレーションモデル、が必要である。
例えば、浄水プロセスや下水処理プロセスではその処理を記述するプロセスシミュレータを作成可能であり、例えば、その出力は水質などである。極値制御を適用する際には、水質だけでなく、水質の良好さとそれにかかる費用などの全体のバランスを見る必要があるため、本実施形態では、最適性を図る評価関数モデルとして、例えば、総合的なコストなどを定義する。この様にプロセスシミュレータと評価関数モデルとから構成されるモデルを直接利用して、極値制御による位相推定を行うことができる。
一方で、線形伝達関数モデルで表す事が適切でない場合は、予め、プロセスシミュレータ（制御対象プロセスを模擬する物理化学式で記述されたシミュレーションモデル）を直接用いて位相推定をすることができる。この際、位相推定に極値制御を適用するためには、出力として何等かの評価関数が必要となるため、いわゆるプロセスシミュレータ（制御対象プロセスのシミュレーションモデル）と評価関数モデルとを組み合わせたシミュレーションモデル、が必要である。
例えば、浄水プロセスや下水処理プロセスではその処理を記述するプロセスシミュレータを作成可能であり、例えば、その出力は水質などである。極値制御を適用する際には、水質だけでなく、水質の良好さとそれにかかる費用などの全体のバランスを見る必要があるため、本実施形態では、最適性を図る評価関数モデルとして、例えば、総合的なコストなどを定義する。この様にプロセスシミュレータと評価関数モデルとから構成されるモデルを直接利用して、極値制御による位相推定を行うことができる。
【0082】
また、プロセスシミュレータを利用する場合には、（８）式のような周波数と位相との関係を数式で陽に表現することは一般にはできない。このため、多数の異なる周波数の正弦波をプロセスシミュレータと評価関数モデルとから構成されたモデルに入力し、シミュレーションの入力（操作量）と出力（評価関数値）との波形のずれを評価して、位相遅れを評価する。各々の入力正弦波に対する出力波形の位相遅れを測定し、位相遅れ推定部１１００において（８）式に相当する周波数ωの関数として位相を表す近似関数を作成しておけばよい。
すなわち、位相という概念は、周波数ωに対して定義されるため、入力の周波数ωを色々と変化させて、周波数ω毎に出力の位相がどのぐらい遅れるかを、シミュレーションにより評価し、周波数ωを横軸とし位相遅れを縦軸としてプロットしたときの関数を表すテーブルを予め作成することができる。位相遅れ推定部１１００は、上記関数（テーブル）を用いて位相遅れの推定値を算出できる。
すなわち、位相という概念は、周波数ωに対して定義されるため、入力の周波数ωを色々と変化させて、周波数ω毎に出力の位相がどのぐらい遅れるかを、シミュレーションにより評価し、周波数ωを横軸とし位相遅れを縦軸としてプロットしたときの関数を表すテーブルを予め作成することができる。位相遅れ推定部１１００は、上記関数（テーブル）を用いて位相遅れの推定値を算出できる。
【0083】
あるいは、位相遅れ推定部１１００は、制御対象プロセスモデルから実際の操作量近傍での線形近似モデルを解析的に導出し、線形近似モデルを線形伝達関数で表現した線形伝達関数モデルの位相計算式を適用したシミュレーションにより位相遅れを算出してもよい。
例えば、下水処理プロセスシミュレータなどは非線形の微分方程式で記述されている。この非線形の微分方程式は、動作点（着目する時点の水質の状態や操作量の状態）の近傍において線形近似モデルを導出することができる。線形近似モデルは、元の非線形の微分方程式を近似した線形の微分方程式であり、本実施形態において線形近似モデルは線形の微分方程式の意味を意味する。線形の微分方程式は、例えば（６）式の線形伝達関数の形でも記述可能である。線形近似モデルを線形伝達関数の形で記述することができれば、位相遅れは（８）式に対応する位相遅れの計算式を解析的に求めることができる。位相遅れ推定部１１００は、上記のように線形近似モデルを線形伝達関数の形で表すことにより得られた位相遅れの計算式で、位相遅れの推定値を算出できる。
位相遅れ推定部１１００は、算出した位相遅れの値を復調用ディザー信号生成部１３０に供給する。
例えば、下水処理プロセスシミュレータなどは非線形の微分方程式で記述されている。この非線形の微分方程式は、動作点（着目する時点の水質の状態や操作量の状態）の近傍において線形近似モデルを導出することができる。線形近似モデルは、元の非線形の微分方程式を近似した線形の微分方程式であり、本実施形態において線形近似モデルは線形の微分方程式の意味を意味する。線形の微分方程式は、例えば（６）式の線形伝達関数の形でも記述可能である。線形近似モデルを線形伝達関数の形で記述することができれば、位相遅れは（８）式に対応する位相遅れの計算式を解析的に求めることができる。位相遅れ推定部１１００は、上記のように線形近似モデルを線形伝達関数の形で表すことにより得られた位相遅れの計算式で、位相遅れの推定値を算出できる。
位相遅れ推定部１１００は、算出した位相遅れの値を復調用ディザー信号生成部１３０に供給する。
【0084】
復調用ディザー信号生成部１３０は、極値制御で、評価関数の勾配（１階微分、ヤコビアン）の推定処理をするために必要となる第１周期信号（復調用ディザー信号）を生成する。
ディザー信号としては、典型的には正弦波を用いるが、必ずしも正弦波である必要はなく、周期的信号であれば矩形波や三角波であっても構わない。復調用ディザー信号は、変調用ディザー信号と同じ周期で同じ形状の波形である。
（９）式に、ディザー信号として正弦波を用いた場合の復調用ディザー信号を示す。
Ｄ（ｔ）＝ｓｉｎ（ωｔ＋φ）…………………………………………………（９）
ディザー信号としては、典型的には正弦波を用いるが、必ずしも正弦波である必要はなく、周期的信号であれば矩形波や三角波であっても構わない。復調用ディザー信号は、変調用ディザー信号と同じ周期で同じ形状の波形である。
（９）式に、ディザー信号として正弦波を用いた場合の復調用ディザー信号を示す。
Ｄ（ｔ）＝ｓｉｎ（ωｔ＋φ）…………………………………………………（９）
【0085】
ここで、ωはディザー信号の周波数であり、φは位相遅れ推定部１１００で演算した位相遅れの推定値である。すなわち、位相遅れパラメータ推定部３００で（７）式のモデルを用いた場合には、（８）式のものに相当し、位相遅れφは下記（１０）式で表すことができる。
φ:＝－ωＬ－ｔａｎ－１（ωＴ）…………………………………………………（１０）
本実施形態の最適制御装置では、上記のように復調用ディザー信号生成部１３０で生成される復調用ディザー信号の中に、位相遅れ情報が組み込まれる。
復調用ディザー信号生成部１３０は、位相遅れφを含む復調用ディザー信号を評価関数勾配推定部１４０に供給する。
φ:＝－ωＬ－ｔａｎ－１（ωＴ）…………………………………………………（１０）
本実施形態の最適制御装置では、上記のように復調用ディザー信号生成部１３０で生成される復調用ディザー信号の中に、位相遅れ情報が組み込まれる。
復調用ディザー信号生成部１３０は、位相遅れφを含む復調用ディザー信号を評価関数勾配推定部１４０に供給する。
【0086】
図６は、位相遅れ補償の前後における操作量の時系列データと評価値の時系列データとの一例を概略的に示す図である。
上記のように位相遅れφを含む復調用ディザー信号を用いて、評価関数勾配推定部１４０にて評価関数の勾配を演算することにより、位相遅れ操作量の時系列データｕ´（ｔ）と評価値の時系列データｙ´（ｔ）との間に生じる位相差を略ゼロとすることができる。
上記のように位相遅れφを含む復調用ディザー信号を用いて、評価関数勾配推定部１４０にて評価関数の勾配を演算することにより、位相遅れ操作量の時系列データｕ´（ｔ）と評価値の時系列データｙ´（ｔ）との間に生じる位相差を略ゼロとすることができる。
【0087】
なお、本実施形態の最適制御装置では、ディザー信号の一例として正弦波を用いたが、正弦波ではなく矩形波や三角波などの周期信号をディザー信号として用いることも可能である。矩形波や三角波などの周期信号をディザー信号として用いる場合には、（１０）式で得られる値をそのまま位相遅れφとすることはできない。（９）式による位相遅れ補償は、（９）式の復調用ディザー信号と、（５）式の変調用ディザー信号との間の位相差をゼロにするようにするための補正処理であるため、位相遅れφがディザー信号の周期に対してどのくらいの割合であるかを求め、変調用ディザー信号に対して、その分だけ周期を補正した信号を復調用のディザー信号とすれば良い。
このような位相遅れ補償を機械的に実施するためには、復調用ディザー信号生成部１３０は、例えば、矩形波などの周期信号をフーリエ級数展開して正弦波の級数として表し、その中に、（１０）式に相当する位相遅れを代入していくという方法を採用しても良い。
このような位相遅れ補償を機械的に実施するためには、復調用ディザー信号生成部１３０は、例えば、矩形波などの周期信号をフーリエ級数展開して正弦波の級数として表し、その中に、（１０）式に相当する位相遅れを代入していくという方法を採用しても良い。
【0088】
次に、本実施形態の最適制御装置において、上記位相遅れ補償を行うことにより得られる効果について説明する。
従来、極値制御の安定性と制御性能との関係について、例えば、極値制御の安定解析を行うことにより、安定性を保証するためにパラメータ設定を適切に行う必要があることを指摘し、定性的には制御パラメータ（ディザー信号の振幅、ディザー信号の周波数、積分器のゲイン（積分ゲイン）、ローパスフィルタのカットオフ周波数、ハイパスフィルタのカットオフ周波数など）の値を十分に小さくしないと、安定性が崩れることが示されている。
従来、極値制御の安定性と制御性能との関係について、例えば、極値制御の安定解析を行うことにより、安定性を保証するためにパラメータ設定を適切に行う必要があることを指摘し、定性的には制御パラメータ（ディザー信号の振幅、ディザー信号の周波数、積分器のゲイン（積分ゲイン）、ローパスフィルタのカットオフ周波数、ハイパスフィルタのカットオフ周波数など）の値を十分に小さくしないと、安定性が崩れることが示されている。
【0089】
制御パラメータを「十分に小さく」設定することは、制御動作を「ゆっくり」、「弱く」働かせることを意味している。すなわち、「安定性」と「制御性能」とにはトレードオフの関係があり、安定性を維持するためには、制御性能を落とすことが本質的に必要であり、逆に制御性能を上げようとしてパラメータを大きく設定すると、安定性が崩れ、結果として制御が失敗するリスクが高くなることを示している。
【0090】
本実施形態の最適制御装置による極値制御では、（５）式の変調用ディザー信号と（９）式の復調用ディザー信号との間に位相差φがあり、操作量の時系列データに対する評価値の時系列データの位相遅れ分が補償されている。この位相遅れ補償により、ディザー信号の周波数ωを大きく（すなわちディザー信号の周期を小さく）設定したり、積分ゲインＫを大きく設定したりすることが可能になり、極値探索の制御性能（収束速度）を改善することができる。
【0091】
図７乃至図１０は、極値制御に用いられる周期信号の周波数を変更したときの、制御の安定性と制御性能との関係の一例を説明するための図である。
凝集剤注入プロセスでは、対象プロセスから沈殿池出口濁度や浄水池出口濁度までの応答時間として数時間オーダーの遅れがある。ここでは、応答に数時間オーダーの遅れがある対象プロセスに対して、上記位相遅れ補償を行わずに、周期が１５００分、１０００分および５００分のディザー信号で操作量を駆動したときの目標値、沈殿池濁度、ＰＡＣ注入量および総コストの時間変化のシミュレーション結果を示している。
凝集剤注入プロセスでは、対象プロセスから沈殿池出口濁度や浄水池出口濁度までの応答時間として数時間オーダーの遅れがある。ここでは、応答に数時間オーダーの遅れがある対象プロセスに対して、上記位相遅れ補償を行わずに、周期が１５００分、１０００分および５００分のディザー信号で操作量を駆動したときの目標値、沈殿池濁度、ＰＡＣ注入量および総コストの時間変化のシミュレーション結果を示している。
【0092】
図７乃至図１０によれば、ディザー信号の周期が１５００分のときには、いずれの評価値も収束しているものの、数値が最適値近傍に収束するまでに長時間要した。ディザー信号の周期が１０００分のときには、ディザー信号の周期が１５００分のときよりも最適値近傍に収束する時間が短くなったが、ディザー信号の周期を５００分とすると、評価値が最適値に収束することがなく発散してしまった。
【0093】
図１１乃至図１３は、実施形態の最適制御装置による効果の一例を説明するための図である。
ここでは、本実施形態の最適制御装置において、プラントの時定数Ｔ＝０．６７ｓ、プラントのむだ時間Ｌ＝０．６７ｓ、操作量の最適値Ｕ＊＝２、ディザー信号の周期をプラントの時定数とプラントのむだ時間との和（Ｔ＋Ｌ）の２倍～５倍として、操作量の時間変換をシミュレーションした結果の一例を示している。
また、図１１乃至図１３には、第１比較例の最適制御装置による操作量の時間変化のシミュレーション結果の一例を併せて示している。第１比較例の最適制御装置は、位相遅れ補償を適用していないこと以外は本実施形態の最適制御装置と同じ条件である。
ここでは、本実施形態の最適制御装置において、プラントの時定数Ｔ＝０．６７ｓ、プラントのむだ時間Ｌ＝０．６７ｓ、操作量の最適値Ｕ＊＝２、ディザー信号の周期をプラントの時定数とプラントのむだ時間との和（Ｔ＋Ｌ）の２倍～５倍として、操作量の時間変換をシミュレーションした結果の一例を示している。
また、図１１乃至図１３には、第１比較例の最適制御装置による操作量の時間変化のシミュレーション結果の一例を併せて示している。第１比較例の最適制御装置は、位相遅れ補償を適用していないこと以外は本実施形態の最適制御装置と同じ条件である。
【0094】
シミュレーション結果によれば、本実施形態の最適制御装置によれば、ディザー信号の周期が小さくなるほど操作量が最適値に収束するまでの時間が短くなった。一方で、第１比較例の最適制御装置については、ディザー信号の周期を小さくすると操作量が最適値に収束するまでの時間が長くなる傾向があり、ディザー信号の周期を２．５１ｓ（（Ｔ＋Ｌ）×２）としたときには操作量が発散して最適値に収束しなかった。
【0095】
上記のように、本実施形態の最適制御装置によれば、操作量とプロセス評価値の間に存在する位相差を補償する機能を導入することにより、極値制御の性能を高める（収束速度を上げる）ことが可能になる。これにより、収束速度が遅いことによる極値探索性能の劣化を改善できる。
【0096】
図１４乃至図１７は、実施形態の最適制御装置による効果の他の例を説明するための図である。
ここでは、第３比較例の最適制御装置と本実施形態の最適制御装置とについて、外乱（流入量の日単位の変動）がある場合のシミュレーション結果の一例を示している。
図１４および図１５には、第３比較例の最適制御装置による極値制御のシミュレーション結果の一例と、目標値ＳＶを一定にしたときのシミュレーション結果の一例とを示している。なお、第３比較例の最適制御装置は、ディザー信号の周期を１５００分（≒１４４０分（一日）に相当）とし、位相遅れ補償を適用せずに極値制御により目標値を演算している。
第３比較例の最適制御装置によるシミュレーション結果によれば、流入量の日単位変動の影響により制御が不安定になっていることが分かる。
ここでは、第３比較例の最適制御装置と本実施形態の最適制御装置とについて、外乱（流入量の日単位の変動）がある場合のシミュレーション結果の一例を示している。
図１４および図１５には、第３比較例の最適制御装置による極値制御のシミュレーション結果の一例と、目標値ＳＶを一定にしたときのシミュレーション結果の一例とを示している。なお、第３比較例の最適制御装置は、ディザー信号の周期を１５００分（≒１４４０分（一日）に相当）とし、位相遅れ補償を適用せずに極値制御により目標値を演算している。
第３比較例の最適制御装置によるシミュレーション結果によれば、流入量の日単位変動の影響により制御が不安定になっていることが分かる。
【0097】
図１６および図１７には、本実施形態の最適制御装置による極値制御のシミュレーション結果の一例と、目標値ＳＶを一定にしたときのシミュレーション結果の一例とを示している。なお、本実施形態の最適制御装置は、ディザー信号の周期を９００分とし、位相遅れ補償を適用した極値制御により目標値を演算している。
第３比較例の最適制御装置に対し、本実施形態の最適制御装置のシミュレーション結果によれば、流入量(外乱)の日単位変動の影響が抑制され、制御が安定している。
第３比較例の最適制御装置に対し、本実施形態の最適制御装置のシミュレーション結果によれば、流入量(外乱)の日単位変動の影響が抑制され、制御が安定している。
【0098】
上記のように、本実施形態によれば、極値制御の「安定性」と「制御性能（収束速度）」との両立という課題を、位相補償という操作により克服することで、安定性を維持しながら、制御性能を極力高める（局所）最適値の探索を実現する最適制御装置、最適制御方法、および、コンピュータプログラムを提供することができる。
【0099】
また、本実施形態の最適制御装置によれば、産業界で広く利用されているＰＩＤ制御の調整時に用いられるステップ応答試験などの簡単な応答試験を実施するだけで、位相遅れ補償を行う手段を構成することでき、容易に実務に利用できる極値制御系を構成可能である。
また、本実施形態の最適制御装置において、勾配推定量を正規化する構成を備えることにより、より収束速度の速い極値制御系を構成することが可能になる。
さらに、位相遅れパラメータ推定部３００においてオンラインでパラメータの推定値を算出することにより、ディザー信号の印加に少しの加工を行うことで、ステップ応答などの試験も不要とする収束速度の速い極値制御系を構成することが可能になる。
また、本実施形態の最適制御装置において、勾配推定量を正規化する構成を備えることにより、より収束速度の速い極値制御系を構成することが可能になる。
さらに、位相遅れパラメータ推定部３００においてオンラインでパラメータの推定値を算出することにより、ディザー信号の印加に少しの加工を行うことで、ステップ応答などの試験も不要とする収束速度の速い極値制御系を構成することが可能になる。
【0100】
図１８は、第１実施形態の最適制御装置の一部の他の構成例を概略的に示す図である。
図１８に示すように、正規化信号発生部１５０および勾配推定量正規化部１６０は必須の構成ではなく、省略されても構わない。その場合には、評価関数勾配推定部１４０の出力値が極値探索部１７０に供給される。
このように、正規化信号発生部１５０および勾配推定量正規化部１６０が省略された構成であっても、安定性を維持しながら、制御性能を極力高める（局所）最適値の探索を実現する最適制御装置、最適制御方法、および、コンピュータプログラムを提供することができる。
図１８に示すように、正規化信号発生部１５０および勾配推定量正規化部１６０は必須の構成ではなく、省略されても構わない。その場合には、評価関数勾配推定部１４０の出力値が極値探索部１７０に供給される。
このように、正規化信号発生部１５０および勾配推定量正規化部１６０が省略された構成であっても、安定性を維持しながら、制御性能を極力高める（局所）最適値の探索を実現する最適制御装置、最適制御方法、および、コンピュータプログラムを提供することができる。
【0101】
次に、第２実施形態の最適制御装置、最適制御方法およびコンピュータプログラムについて図面を参照して詳細に説明する。
なお、以下の説明において上述の第１実施形態と同様の構成については同一の符号を付して説明を省略する。
図１９は、第２実施形態の最適制御装置の一構成例を概略的に示す図である。
本実施形態の最適制御装置は、復調用ディザー信号を用いずに、勾配推定時に位相遅れの情報を組み込んでいる点において上述の第１実施形態と異なっている。
なお、以下の説明において上述の第１実施形態と同様の構成については同一の符号を付して説明を省略する。
図１９は、第２実施形態の最適制御装置の一構成例を概略的に示す図である。
本実施形態の最適制御装置は、復調用ディザー信号を用いずに、勾配推定時に位相遅れの情報を組み込んでいる点において上述の第１実施形態と異なっている。
【0102】
本実施形態の最適制御装置は、プロセス計測値取得部１１０と、プロセス評価値算出部１２０と、評価関数勾配推定部１４０と、極値探索部（最適操作量適応調整部）１７０と、操作量出力部１９０と、プロセス位相遅れ推定部４００と、を備えている。プロセス位相遅れ推定部４００は、位相遅れパラメータ推定部３００と、位相遅れ推定部１１００と、を含む。
【0103】
図２０は、第２実施形態の最適制御装置の一部の構成例を概略的に示す図である。
位相遅れ推定部１１００は、位相遅れを（８）式の様な形で陽に推定するのではなく、（６）式や（７）式のモデルに対して、操作量のデータをリアルタイムで入力した出力の時系列データ（以下、位相遅れ補正操作量データと呼ぶ）を位相遅れ推定部１１００の出力とする。
プロセスシミュレータを用いる場合も同様に、位相遅れ推定部１１００は、操作量データをリアルタイムでプロセスシミュレータに入力した出力の時系列データを、出力とすることができる。
位相遅れ推定部１１００は、位相遅れを（８）式の様な形で陽に推定するのではなく、（６）式や（７）式のモデルに対して、操作量のデータをリアルタイムで入力した出力の時系列データ（以下、位相遅れ補正操作量データと呼ぶ）を位相遅れ推定部１１００の出力とする。
プロセスシミュレータを用いる場合も同様に、位相遅れ推定部１１００は、操作量データをリアルタイムでプロセスシミュレータに入力した出力の時系列データを、出力とすることができる。
【0104】
これにより、操作量と評価値との間の位相のずれを補正できる。すなわち、操作量と評価値との間には、ダイナミクスを持つ制御対象２００を通すことによる位相遅れが必ず存在する。本実施形態の最適制御装置では、操作量に替えて位相遅れ補正操作量を用いることで、位相遅れ補正操作量と評価値との位相遅れを極めて小さくしている。原理的には、線形伝達関数モデルと実際の制御対象との誤差が無ければ、位相遅れ補正操作量と評価値との間の位相遅れをなくすことができる。
【0105】
本実施形態では、評価関数勾配推定部１４０は、復調用ディザー信号を用いずに、動作点における評価関数の傾き（勾配）を直接算出するする。
先に述べたように、例えば図２に示した評価関数の全体の形状は、実際に極値制御を実施している際には不明であるが、動作点近傍の値は制御実施時にリアルタイムに取得できている。従って、動作点近傍のデータ、すなわち、制御実施時の直近の操作量と評価量との時系列データを用いることにより、図２に示すような評価関数の勾配を表す直線の傾きを推定することが可能である。
先に述べたように、例えば図２に示した評価関数の全体の形状は、実際に極値制御を実施している際には不明であるが、動作点近傍の値は制御実施時にリアルタイムに取得できている。従って、動作点近傍のデータ、すなわち、制御実施時の直近の操作量と評価量との時系列データを用いることにより、図２に示すような評価関数の勾配を表す直線の傾きを推定することが可能である。
【0106】
例えば、制御実施時に過去の操作量時系列データと評価量時系列データとを、所定期間分、蓄積し、蓄積された時系列データ利用して最小２乗法や逐次最小２乗法などの回帰アルゴリズムにより次式（１１）の直線の係数を求めることができる。
Ｊ＝ｃＵ＋ｄ………………………………………………………………………（１１）
ここで、Ｊは図２の評価量（プロセス評価値）であり、Ｕは図２の操作量である。ｃとｄとは直線回帰の係数である。これらの係数ｃ、ｄは、最小２乗法などのアルゴリズムを用いて求めることができる。
（１１）式の直線の係数が求まると、直線の傾き＝評価関数の勾配と見なすことができるので、傾きに対応する係数ｃを評価関数の勾配情報（勾配推定値）として用いることができる。
Ｊ＝ｃＵ＋ｄ………………………………………………………………………（１１）
ここで、Ｊは図２の評価量（プロセス評価値）であり、Ｕは図２の操作量である。ｃとｄとは直線回帰の係数である。これらの係数ｃ、ｄは、最小２乗法などのアルゴリズムを用いて求めることができる。
（１１）式の直線の係数が求まると、直線の傾き＝評価関数の勾配と見なすことができるので、傾きに対応する係数ｃを評価関数の勾配情報（勾配推定値）として用いることができる。
【0107】
本実施形態では、評価関数勾配推定部１４０は、制御実施時に過去の操作量時系列データと評価量時系列データとを、所定期間分、蓄積するデータ蓄積部（図示せず）を備え、蓄積されたデータを用いて、（１１）式に替えて以下の（１２）式で係数ｃを勾配推定値として求めている。
Ｊ＝ｃＵｆ＋ｄ……………………………………………………………………（１２）
ここで、Ｕｆは実際の操作量ではなく、操作量に対して位相遅れの推定値の情報を加えた位相遅れ操作量であり、操作量を（６）式や（７）式に示した線形伝達関数を通したものである。すなわち、位相遅れ推定部１１００の線形伝達関数モデルの出力信号の値である。
Ｊ＝ｃＵｆ＋ｄ……………………………………………………………………（１２）
ここで、Ｕｆは実際の操作量ではなく、操作量に対して位相遅れの推定値の情報を加えた位相遅れ操作量であり、操作量を（６）式や（７）式に示した線形伝達関数を通したものである。すなわち、位相遅れ推定部１１００の線形伝達関数モデルの出力信号の値である。
【0108】
上記のように操作量そのものではなく、操作量に対して線形伝達関数モデルを通した値を用いることにより、図２に示すように操作量の位相遅れがない状態で評価関数の勾配推定を行うことができ、勾配推定の精度の向上が見込める。
また、理想的には、（６）式や（７）式の線形伝達関数モデルにより、制御対象２００における位相遅れが正確に推定されているならば、操作量の位相遅れをゼロにすることができる。このため、ディザー信号の周波数をいくらでも速くすることができることから、「安定性」と「制御性能」とのトレードオフ問題を解消することができ、制御の安定性を維持したまま、制御性能を向上させる（収束速度を上げる）ことが可能になる。
また、理想的には、（６）式や（７）式の線形伝達関数モデルにより、制御対象２００における位相遅れが正確に推定されているならば、操作量の位相遅れをゼロにすることができる。このため、ディザー信号の周波数をいくらでも速くすることができることから、「安定性」と「制御性能」とのトレードオフ問題を解消することができ、制御の安定性を維持したまま、制御性能を向上させる（収束速度を上げる）ことが可能になる。
【0109】
なお、線形伝達関数モデルにより、制御対象における位相遅れを正確に同定することは実際には不可能であるため、制御性能（収束速度）の限界はあるが、操作量の時系列データと評価値の時系列データとの位相差を解消することにより、安定性を損なわずに、制御性能を格段に向上させることができる。
【0110】
本実施形態において、勾配推定値を算出する方法は上記に限定されるものではない。例えば、操作量と評価値との「相関係数」を勾配推定値と見なしてもよい。これは、操作量と評価とを表す関数の勾配の傾きの正負は、操作量と評価値との相関係数の正負と一致するとともに、相関係数は勾配の大きさに依存しない、勾配の代用となるからである。すなわち、変数Ｘおよび変数Ｙの直線回帰の回帰係数と、変数Ｘおよび変数Ｙの相関係数との間には、「変数Ｘと変数Ｙとを各々の平均と標準偏差で正規化した上で、回帰直線を求めると、回帰直線の係数ｃは相関係数ｒと一致し、バイアスｂ＝０となる」という事実に基づいている。そのため、回帰係数ｃに替えて相関係数ｒを用いることは、操作量時系列データと評価値時系列データとを各々正規化した上で、回帰係数を求めていることに相当する。
【0111】
本実施形態の最適制御装置においても、評価関数勾配推定部１４０は、勾配推定値を評価量Ｊと位相遅れ操作量Ｕｆとの相関係数により算出することができる。本実施形態の最適制御装置において回帰係数に替えて相関係数を勾配情報として用いることで、第１実施形態の最適制御装置における正規化操作を行うことと同様の効果を得ることができる。
【0112】
なお、図１８における、勾配を直線近似で推定する極値制御の構成は、目標値追従型のいわゆる通常のフィードバック制御において広く知られている、内部モデル制御（Internal Model Control）やスミス補償制御（内部モデル制御の一種）と類似した構造を持っている。
目標値追従型のフィードバック制御において、内部モデル制御は、制御対象と、制御対象を表す制御対象プロセスモデルとを並列に並べ、制御対象の出力と制御対象を表すプロセスモデルの出力との誤差をフィードバックすることで、制御性能を高める制御であり、予測制御の一種ととらえることもでき、特にむだ時間が長いプロセスに対して有効であることが広く知られている。
目標値追従型のフィードバック制御において、内部モデル制御は、制御対象と、制御対象を表す制御対象プロセスモデルとを並列に並べ、制御対象の出力と制御対象を表すプロセスモデルの出力との誤差をフィードバックすることで、制御性能を高める制御であり、予測制御の一種ととらえることもでき、特にむだ時間が長いプロセスに対して有効であることが広く知られている。
【0113】
本実施形態の最適制御装置では、極値探索型のフィードバック制御において、内部モデル制御と同様の考え方を導入したものと考えることができる。すなわち、遅れ時間が長い場合に、それを予測するモデルを制御対象と並列に並べ、一種の予測を行い、制御対象における遅れを補償した操作量に基づいて最適化（勾配推定）を行っている。
本実施形態の最適制御装置における極値探索型のフィードバック制御と、内部モデル制御やスミス補償制御との類似点は、制御対象のモデルを用いて予測を行っている点である。
本実施形態の最適制御装置における極値探索型のフィードバック制御と、内部モデル制御やスミス補償制御との類似点は、制御対象のモデルを用いて予測を行っている点である。
【0114】
一方、相違点は、（１）内部モデル制御では、制御対象プロセスモデルを利用するが、本実施形態においては、制御対象の静的な特性は不要であるため、（制御対象プロセスモデルを利用しても良いが）制御対象プロセスモデルの中の静的な特性（静的な非線形要素）を除く位相特性をあらわすモデル利用するだけで良い点と、（２）内部モデル制御では、制御対象の出力値と制御対象プロセスモデルの予測値との誤差をフィードバックするのに対し、本実施形態では、プロセスの位相特性の予測値と制御対象の出力値（＝評価関数値）との関係を表す勾配を求めて勾配を推定する点、である。
内部モデル制御については既にその効果が明らかになっており、上記のような本実施形態の最適制御装置におけるフィードバック制御と、内部モデル制御との類似性の視点からも本実施形態の効果を類推することができる。
内部モデル制御については既にその効果が明らかになっており、上記のような本実施形態の最適制御装置におけるフィードバック制御と、内部モデル制御との類似性の視点からも本実施形態の効果を類推することができる。
【0115】
図１１乃至図１３に、第２実施形態と第２比較例の最適制御装置による操作量の時間変化のシミュレーション結果の一例を示している。
ここでは、本実施形態の最適制御装置において、プラントの時定数Ｔ＝０．６７ｓ、プラントのむだ時間Ｌ＝０．６７ｓ、操作量の最適値Ｕ＊＝２、ディザー信号の周期をプラントの時定数とプラントのむだ時間との和（Ｔ＋Ｌ）の２倍～５倍として、操作量の時間変換をシミュレーションした結果の一例を示している。
第２比較例の最適制御装置は、位相遅れ補償を適用していないこと以外は本実施形態の最適制御装置と同じ条件である。
ここでは、本実施形態の最適制御装置において、プラントの時定数Ｔ＝０．６７ｓ、プラントのむだ時間Ｌ＝０．６７ｓ、操作量の最適値Ｕ＊＝２、ディザー信号の周期をプラントの時定数とプラントのむだ時間との和（Ｔ＋Ｌ）の２倍～５倍として、操作量の時間変換をシミュレーションした結果の一例を示している。
第２比較例の最適制御装置は、位相遅れ補償を適用していないこと以外は本実施形態の最適制御装置と同じ条件である。
【0116】
シミュレーション結果によれば、本実施形態の最適制御装置および第２比較例の最適制御装置によれば、ディザー信号の周期が小さくなるほど操作量が最適値に収束するまでの時間が短くなる傾向がみられた。また、ディザー信号がいずれの周期の場合にも、本実施形態の最適制御装置によれば操作量が最適値に収束するまでの時間が第２比較例の最適制御装置よりも短かった。
【0117】
上記のように、本実施形態の最適制御装置によれば、操作量とプロセス評価値の間に存在する位相差を補償する機能を導入することにより、極値制御の性能を高める（収束を速くする）ことが可能になる。これにより、収束速度が遅いことによる極値探索性能の劣化を改善できる。
【0118】
すなわち、本実施形態によれば、強制的なディザー信号を印加せずに、安定性を維持しながら、制御性能を極力高める（局所）最適値の探索を実現する最適制御装置、最適制御方法、および、コンピュータプログラムを提供することができる。
【0119】
図２１は、第２実施形態の最適制御装置の構成の変形例を概略的に示す図である。
この例では、最適制御装置が変調用ディザー信号生成部１８０を備えている点において、図１７および図１８に示す例と相違している。
図２２は、第２実施形態の最適制御装置の一部の構成の変形例を概略的に示す図である。
本実施形態では、変調用ディザー信号を操作量に組み込むことにより、回帰推定や相関推定に用いるデータに強制的に変動を与えている。例えば、制御対象２００が外乱などによって十分に励起（駆動）されていない場合には、回帰パラメータや相関係数の推定値の信頼性が劣化する可能性があるため、変調用ディザー信号を操作量に加えておく方が好ましい。
この例では、最適制御装置が変調用ディザー信号生成部１８０を備えている点において、図１７および図１８に示す例と相違している。
図２２は、第２実施形態の最適制御装置の一部の構成の変形例を概略的に示す図である。
本実施形態では、変調用ディザー信号を操作量に組み込むことにより、回帰推定や相関推定に用いるデータに強制的に変動を与えている。例えば、制御対象２００が外乱などによって十分に励起（駆動）されていない場合には、回帰パラメータや相関係数の推定値の信頼性が劣化する可能性があるため、変調用ディザー信号を操作量に加えておく方が好ましい。
【0120】
ただし、この場合、変調用ディザー信号は必ずしも周期的な信号である必要もなく正弦波である必要もない。例えば、乱数の様な信号を変調用ディザー信号とし、制御対象に印加する方法でも構わない。
制御対象が外乱などによって十分に駆動されている場合には、上述の様に変調用ディザー信号生成部１８０を備えない極値制御系を構成することも可能である。
制御対象が外乱などによって十分に駆動されている場合には、上述の様に変調用ディザー信号生成部１８０を備えない極値制御系を構成することも可能である。
【0121】
本変形例は、上記以外の構成は図１９および図２０の最適制御装置と同様の構成である。したがって、本変形例によれば、必要最小限の外部入力の印加により、安定性を維持しながら、制御性能を極力高める（局所）最適値の探索を実現する最適制御装置、最適制御方法、および、コンピュータプログラムを提供することができる。
【0122】
本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれる。
【符号の説明】
【0123】
１００…制御システム、１１０…プロセス計測値取得部、１２０…プロセス評価値算出部、１３０…復調用ディザー信号生成部、１４０…評価関数勾配推定部、１５０…正規化信号発生部、１６０…勾配推定量正規化部、１７０…極値探索部、１７１…積分器、１７２…ゲイン乗算部、１８０…変調用ディザー信号生成部、１９０…操作量出力部、２００…制御対象（対象プラント）、３００…位相遅れパラメータ推定部、４００…プロセス位相遅れ推定部、１１００…位相遅れ推定部
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【図21】
【図22】
