JP2024044583A - PLANT RESPONSE ESTIMATION APPARATUS, PLANT RESPONSE ESTIMATION METHOD, AND PROGRAM - Google Patents
PLANT RESPONSE ESTIMATION APPARATUS, PLANT RESPONSE ESTIMATION METHOD, AND PROGRAM Download PDFInfo
- Publication number
- JP2024044583A JP2024044583A JP2022150202A JP2022150202A JP2024044583A JP 2024044583 A JP2024044583 A JP 2024044583A JP 2022150202 A JP2022150202 A JP 2022150202A JP 2022150202 A JP2022150202 A JP 2022150202A JP 2024044583 A JP2024044583 A JP 2024044583A
- Authority
- JP
- Japan
- Prior art keywords
- model
- plant
- response
- parameters
- calculation unit
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 230000004044 response Effects 0.000 title claims abstract description 119
- 238000000034 method Methods 0.000 title claims description 30
- PWPJGUXAGUPAHP-UHFFFAOYSA-N lufenuron Chemical compound C1=C(Cl)C(OC(F)(F)C(C(F)(F)F)F)=CC(Cl)=C1NC(=O)NC(=O)C1=C(F)C=CC=C1F PWPJGUXAGUPAHP-UHFFFAOYSA-N 0.000 title 1
- 238000004364 calculation method Methods 0.000 claims abstract description 116
- 230000006641 stabilisation Effects 0.000 claims abstract description 61
- 238000011105 stabilization Methods 0.000 claims abstract description 61
- 230000000087 stabilizing effect Effects 0.000 claims abstract description 20
- 239000013598 vector Substances 0.000 claims description 58
- 241001123248 Arma Species 0.000 claims description 14
- 230000005484 gravity Effects 0.000 claims description 10
- 238000002156 mixing Methods 0.000 claims description 4
- 238000005516 engineering process Methods 0.000 abstract description 2
- 239000000872 buffer Substances 0.000 description 37
- 238000010586 diagram Methods 0.000 description 16
- 238000005070 sampling Methods 0.000 description 16
- 238000006243 chemical reaction Methods 0.000 description 14
- 238000012952 Resampling Methods 0.000 description 10
- 230000008569 process Effects 0.000 description 9
- 239000011159 matrix material Substances 0.000 description 7
- 238000013459 approach Methods 0.000 description 4
- 238000004891 communication Methods 0.000 description 4
- 230000003190 augmentative effect Effects 0.000 description 3
- 238000005316 response function Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 2
- 230000006870 function Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 101000802640 Homo sapiens Lactosylceramide 4-alpha-galactosyltransferase Proteins 0.000 description 1
- 102100035838 Lactosylceramide 4-alpha-galactosyltransferase Human genes 0.000 description 1
- 229910000831 Steel Inorganic materials 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 230000014509 gene expression Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000010248 power generation Methods 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 239000010959 steel Substances 0.000 description 1
Images
Landscapes
- Feedback Control In General (AREA)
Abstract
【課題】操業中のプラントの運転データから応答モデルのパラメータを安定的に推定することができる技術を提供すること。【解決手段】本開示の一態様によるプラント応答推定装置は、制御対象プラントの運転データに基づいて、前記制御対象プラントのプラント応答モデルを表す関数のパラメータであるモデルパラメータの推定値を逐次的に計算するように構成されている推定計算部と、前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算するように構成されている安定化計算部と、前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算するように構成されている応答計算部と、を有し、前記推定計算部は、前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算部によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算するように構成されている。【選択図】図2[Problem] To provide a technology capable of stably estimating parameters of a response model from operational data of a plant in operation. [Solution] A plant response estimation device according to one aspect of the present disclosure includes an estimation calculation unit configured to sequentially calculate, based on operational data of a controlled plant, estimates of model parameters, which are parameters of a function representing a plant response model of the controlled plant, a stabilization calculation unit configured to calculate stabilized model parameters obtained by stabilizing the estimates of the model parameters, and a response calculation unit configured to calculate a response of the controlled plant by the function representing the plant response model based on the estimates of the model parameters, and the estimation calculation unit is configured to calculate the estimates of the model parameters based on the operational data and the stabilized model parameters obtained by stabilizing the estimates of the model parameters calculated previously by the stabilization calculation unit. [Selected Figure] Figure 2
Description
本開示は、プラント応答推定装置、プラント応答推定方法、及びプログラムに関する。 The present disclosure relates to a plant response estimation device, a plant response estimation method, and a program.
温調制御装置やPLC(Programmable Logic Controller)、DCS(Distributed Control System)等の制御装置、パーソナルコンピュータや組み込み制御機器上で実装される制御装置等が産業上広く利用されている。 Controllers such as temperature control devices, programmable logic controllers (PLCs), distributed control systems (DCSs), and controllers implemented on personal computers and embedded control devices are widely used in industry.
また、制御対象の制御量を目標値に追従させることを目的とする制御方式として、PID(Proportional-Integral-Differential)制御、モデル予測制御、内部モデル制御、LQG(Linear-Quadratic-Gaussian)制御、H2制御、H∞制御等の各種の制御方式が知られている。 In addition, various control methods are known that aim to make the control amount of the controlled object follow a target value, such as PID (Proportional-Integral-Differential) control, model predictive control, internal model control, LQG (Linear-Quadratic-Gaussian) control, H2 control, and H∞ control.
モデル予測制御は、制御対象の状態空間モデルや将来の時間応答モデルを用いた最適化計算を逐次的に行うことで望ましい応答を得る方式であり、産業界で広く用いられている(例えば、特許文献1、非特許文献1)。このようなモデル予測制御等で用いられる線形予測モデルに関してそのモデルパラメータを推定(同定)する手法として、逐次最小2乗法(RLS法:Recursive Least Squares法)やカルマンフィルタ法等が知られている(例えば、非特許文献2、非特許文献3)。
Model predictive control is a method to obtain a desired response by sequentially performing optimization calculations using a state space model of the controlled object and a future time response model, and is widely used in industry (for example, patented
なお、制御対象プラントの定常状態を表すオフセット項が含まれるモデルパラメータを推定する技術も知られている(例えば、特許文献3)。 In addition, there is also a known technique for estimating model parameters that include an offset term that represents the steady state of the controlled plant (for example, Patent Document 3).
逐次最小2乗法やカルマンフィルタ法は、操業中(運転中)のプラントのデータから逐次的にモデルパラメータを推定できる点で大変有益であるが、モデルパラメータの推定値が逐次的に更新されるため以下の2つの課題がある。 The iterative least squares method and the Kalman filter method are very useful in that they allow model parameters to be estimated sequentially from data of an operating (operating) plant, but since the estimated values of the model parameters are updated sequentially, the following There are two issues.
1つ目の課題は、真値に収束するまでの途中の段階ではモデルパラメータの推定値の信頼性が担保されないという点である。2つ目の課題は、真のモデルパラメータの次数が事前に未知であるため、推定対象のモデルパラメータが過剰な次数を持つ場合には多重共線性が生じ得るという点である。 The first problem is that the reliability of estimated values of model parameters is not guaranteed during the intermediate stage until convergence to the true value. The second problem is that since the degree of the true model parameter is unknown in advance, multicollinearity may occur if the model parameter to be estimated has an excessive degree.
上記の2つの課題により、例えば、本来安定であるはずのプラントに対して不安定極を持つようなモデルパラメータが途中で計算されたり、観測ノイズの影響によってモデルパラメータが大きく変動したり、多重共線性によって複数の局所最適解が生じたりする等といった不安定性が生じる恐れがある。 Due to the above two issues, for example, model parameters that have unstable poles for a plant that is supposed to be stable may be calculated midway through, model parameters may fluctuate greatly due to the influence of observation noise, or multiple There is a risk that instability may occur due to linearity, such as multiple locally optimal solutions.
オンラインでモデルパラメータの更新と制御対象の制御とを実現する場合、モデルパラメータの不安定性が制御に重要な影響を与えかねないため、収束途中におけるモデルパラメータの安定性を確保することは重要である。 When updating model parameters and controlling the controlled object online, it is important to ensure the stability of the model parameters during convergence, since instability of the model parameters can have a significant impact on the control.
本開示は、上記の点に鑑みてなされたもので、操業中のプラントの運転データから応答モデルのパラメータを安定的に推定することができる技術を提供する。 The present disclosure has been made in view of the above points, and provides a technique that can stably estimate parameters of a response model from operating data of a plant in operation.
本開示の一態様によるプラント応答推定装置は、制御対象プラントの運転データに基づいて、前記制御対象プラントのプラント応答モデルを表す関数のパラメータであるモデルパラメータの推定値を逐次的に計算するように構成されている推定計算部と、前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算するように構成されている安定化計算部と、前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算するように構成されている応答計算部と、を有し、前記推定計算部は、前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算部によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算するように構成されている。 A plant response estimation device according to one aspect of the present disclosure includes an estimation calculation unit configured to sequentially calculate, based on operation data of a controlled plant, estimates of model parameters, which are parameters of a function representing a plant response model of the controlled plant; a stabilization calculation unit configured to calculate stabilized model parameters obtained by stabilizing the estimates of the model parameters; and a response calculation unit configured to calculate a response of the controlled plant using a function representing the plant response model based on the estimates of the model parameters, and the estimation calculation unit is configured to calculate an estimate of the model parameters based on the operation data and the stabilized model parameters obtained by stabilizing the estimates of the previously calculated model parameters by the stabilization calculation unit.
操業中のプラントの運転データから応答モデルのパラメータを安定的に推定することができる技術が提供される。 A technology is provided that can stably estimate response model parameters from operational data of an operating plant.
以下、本発明の一実施形態について説明する。以下の各実施形態では、操業中のプラントの運転データから応答モデルのパラメータ(モデルパラメータ)を安定的に推定すると共に、そのパラメータによりプラント応答を推定することができるプラント応答推定装置10について説明する。なお、操業中とはプラントが正常に稼働し、通常の運転を行っている状態のことであり、例えば、オンライン中、運転中、運用中等と呼ばれてもよい。また、プラントとは1以上の機械、機器、装置等で構成される産業設備のことであり、プラント応答モデルを用いたモデル予測制御によって制御される制御対象である。プラントの具体例としては、例えば、石油化学プラント、食品プラント、鉄鋼プラント、発電プラント等といったものが挙げられるが、これらは一例であって、これらに限られるものではない。
An embodiment of the present invention will be described below. In each of the following embodiments, a plant
[第一の実施形態]
以下、第一の実施形態について説明する。
[First embodiment]
The first embodiment will be described below.
<プラント応答推定装置10のハードウェア構成例>
本実施形態に係るプラント応答推定装置10のハードウェア構成例を図1に示す。図1に示すように、本実施形態に係るプラント応答推定装置10は、入力装置11と、表示装置12と、外部I/F13と、通信I/F14と、プロセッサ15と、メモリ装置16とを有する。これらの各ハードウェアは、それぞれがバス17を介して通信可能に接続される。
<Hardware configuration example of the plant
An example of a hardware configuration of a plant
入力装置11は、例えば、キーボード、マウス、タッチパネル、各種物理ボタン等である。表示装置12は、例えば、ディスプレイ、表示パネル等である。なお、プラント応答推定装置10は、例えば、入力装置11及び表示装置12のうちの少なくとも一方を有していなくてもよい。
The
外部I/F13は、記録媒体13a等の外部装置とのインタフェースである。記録媒体13aとしては、例えば、CD(Compact Disc)、DVD(Digital Versatile Disk)、SDメモリカード(Secure Digital memory card)、USB(Universal Serial Bus)メモリカード等が挙げられる。
The external I/
通信I/F14は、プラント応答推定装置10を通信ネットワークに接続するためのインタフェースである。プロセッサ15は、例えば、CPU(Central Processing Unit)等の各種演算装置である。メモリ装置16は、例えば、HDD(Hard Disk Drive)、SSD(Solid State Drive)、RAM(Random Access Memory)、ROM(Read Only Memory)、フラッシュメモリ等の各種記憶装置である。
Communication I/
なお、図1に示すハードウェア構成は一例であって、プラント応答推定装置10はこれに限られるものではなく、他のハードウェア構成であってもよい。例えば、プラント応答推定装置10は、複数のプロセッサ15や複数のメモリ装置16を有していてもよいし、図示したハードウェアのうちの一部を有していなくてもよいし、図示したハードウェア以外の種々のハードウェアを有していてもよい。
Note that the hardware configuration shown in FIG. 1 is an example, and the plant
<プラント応答推定装置10の機能構成例>
本実施形態に係るプラント応答推定装置10の機能構成例を図2に示す。図2に示すように、本実施形態に係るプラント応答推定装置10は、モデルパラメータ推定部101と、ステップ応答計算部102とを有する。これら各部は、例えば、プラント応答推定装置10にインストールされた1以上のプログラムが、プロセッサ15等に実行させる処理により実現される。
<Example of functional configuration of plant
An example of a functional configuration of the plant
モデルパラメータ推定部101は、サンプリング周期Δ毎に、制御対象のプラント(又はその運転状態を計測するセンサ等の機器)から運転データの観測値(制御量y、操作量u及び外乱v)を受信し、そのサンプリング周期Δに応じて、プラント応答モデルのモデルパラメータを推定する。なお、運転データは、例えば、計測データや観測データ等と呼ばれてもよい。
The model
ここで、時刻をtとすれば制御量y、操作量u及び外乱vはそれぞれy(t)、u(t)及びv(t)と表され、各サンプリング時刻tk(kはサンプリング時刻を表すインデックス)に関してtk+1-tk=Δが成り立つ。なお、kは0以上の整数である。以下では、特に断らない限り、各サンプリング時刻tkをそのインデックスkと同一視し、y(k)、u(k)及びv(k)とも表すことにする。 Here, if time is t, the controlled variable y, manipulated variable u, and disturbance v are represented as y(t), u(t), and v(t), respectively, and t k+1 - t k = Δ holds for each sampling time t k (k is an index representing the sampling time), where k is an integer equal to or greater than 0. In the following, unless otherwise specified, each sampling time t k will be regarded as the same as its index k, and will also be represented as y(k), u(k), and v(k).
また、以下では、プラント応答モデルは、パラメータθを持つプラント応答関数Sθ(t)で表されるものとして、パラメータθをモデルパラメータということにする。プラント応答モデルとしては様々なプラント応答関数Sθで表現されるモデルを採用することが可能であるが、以下では、主に、ARMAモデル(Autoregressive Moving Average Model)といった多項式モデルを想定する。プラント応答モデルがARMAモデル等といった多項式モデルである場合、モデルパラメータθはその多項式の係数を要素として持つベクトルとなる。なお、プラント応答モデルとしてARMAモデル以外の様々なモデルが採用可能であることは言うまでもない。 In the following, the plant response model is represented by a plant response function S θ (t) having a parameter θ, and the parameter θ is referred to as a model parameter. As the plant response model, various models represented by a plant response function S θ can be adopted, but in the following, a polynomial model such as an ARMA model (Autoregressive Moving Average Model) is mainly assumed. When the plant response model is a polynomial model such as an ARMA model, the model parameter θ is a vector having coefficients of the polynomial as elements. It goes without saying that various models other than the ARMA model can be adopted as the plant response model.
ステップ応答計算部102は、モデルパラメータ推定部101によって推定されたモデルパラメータ(以下、モデルパラメータ推定値ともいう。)θを用いて、与えられた時刻tにおけるプラントのステップ応答の推定値(以下、ステップ応答推定値ともいう。)Sθ(t)を計算する。なお、ステップ応答とは、操作量としてステップ信号がプラントに印加されたときの制御量のことである。
The step
ここで、モデルパラメータ推定部101には、バッファ部111と、状態ベクトル変換部112と、逐次推定計算部113と、安定化モデルパラメータ計算部114とが含まれる。
Here, the model
バッファ部111は、或る所定の期間における制御量y(k)、操作量u(k)及び外乱v(k)をメモリ装置16に蓄積(バッファ)する。状態ベクトル変換部112は、メモリ装置16にバッファされている制御量y(k)、操作量u(k)及び外乱v(k)の再サンプリングを行って、状態ベクトルと呼ぶベクトルを作成する。逐次推定計算部113は、状態ベクトルを用いて、RLS法により、プラント応答モデルのモデルパラメータを推定する。安定化モデルパラメータ計算部114は、モデルパラメータから安定化モデルパラメータを計算する。安定化モデルパラメータとは、モデルパラメータを或る統計量(例えば、自己回帰係数の平均値、移動平均係数の平均値等)により安定化したベクトルのことである。以下では、安定化モデルパラメータをθSと表し、またサンプリング時刻tkにおける安定化モデルパラメータをθS(k)と表す。
The
<バッファ部111の動作例>
サンプリング時刻tkにおける制御量バッファをY(k)、操作量バッファをU(k)、外乱バッファをV(k)として、これらの各バッファは以下のベクトルで表されるものとする。
<Example of operation of
Assume that the control amount buffer at sampling time t k is Y(k), the operation amount buffer is U(k), and the disturbance buffer is V(k), and each of these buffers is represented by the following vector.
すなわち、制御量バッファY(k)にはk-B1からkまでの制御量y、操作量バッファU(k)にはk-B2からkまでの操作量u、外乱バッファV(k)にはk-B3からkまでの外乱vがそれぞれ格納されているものとする。ここで、B1、B2及びB3はそれぞれ0以上の整数であり、制御量バッファ、操作量バッファ及び外乱バッファの大きさを決めるパラメータである。これらのB1、B2及びB3はメモリ装置16のサイズ等に応じて、適宜、その値が決定される。
That is, the control amount buffer Y(k) stores the control amount y from k- B1 to k, the manipulated variable buffer U(k) stores the manipulated variable u from k- B2 to k, and the disturbance buffer V(k) stores the disturbance v from k- B3 to k. Here, B1 , B2 , and B3 are each integers equal to or greater than 0, and are parameters that determine the sizes of the control amount buffer, the manipulated variable buffer, and the disturbance buffer. The values of B1 , B2 , and B3 are determined appropriately depending on the size of the
このとき、バッファ部111は、制御量y(k)、操作量u(k)及び外乱v(k)を受信すると、図3に示すように、制御量y(k)と制御量バッファY(k-1)から制御量バッファY(k)、操作量u(k)と操作量バッファU(k-1)から操作量バッファU(k)、外乱v(k)と外乱バッファV(k-1)から外乱バッファV(k)にそれぞれ更新する。
At this time, when the
具体的には、バッファ部111は、制御量バッファY(k-1)に格納されているy(k-B1-1)を削除した上で、新たに観測された制御量y(k)を格納することで、y(k-B1)からy(k)までの制御量が格納された制御量バッファY(k)に更新する。同様に、バッファ部111は、操作量バッファU(k-1)に格納されている操作量u(k-B2-1)を削除した上で、新たに観測された操作量u(k)を格納することで、u(k-B2)からu(k)までの操作量が格納された操作量バッファU(k)に更新する。同様に、バッファ部111は、外乱バッファV(k-1)に格納されている外乱v(k-B3-1)を削除した上で、新たに観測された外乱v(k)を格納することで、v(k-B3)からv(k)までの外乱が格納された外乱バッファV(k)に更新する。
Specifically, the
<状態ベクトル変換部112の動作例>
再サンプリング周期をDとする。このとき、状態ベクトル変換部112は、制御量バッファY(k)、操作量バッファU(k)及び外乱バッファV(k)が更新されると、これらの各バッファY(k)、U(k)及びV(k)から再サンプリング周期Dで再サンプリングを行って、以下の再サンプリング制御量ベクトルYD(k)、再サンプリング操作量ベクトルUD(k)及び再サンプリング外乱ベクトルVD(k)をそれぞれ作成する。
<Example of operation of state
Let D be the resampling period. At this time, when the control amount buffer Y(k), the operation amount buffer U(k), and the disturbance buffer V(k) are updated, the state
なお、一般に、N、M及びLの値は大きい方が多様な表現が可能で、プラント応答の高精度な予測が期待できるが、モデルパラメータθの推定のために多くの計算資源やメモリ量が必要となる。 In general, larger values of N, M, and L allow for more diverse expressions and are expected to enable more accurate prediction of plant response, but they require more computational resources and memory to estimate the model parameter θ.
そして、状態ベクトル変換部112は、再サンプリング制御量ベクトルYD(k)、再サンプリング操作量ベクトルUD(k)及び再サンプリング外乱ベクトルVD(k)から以下により状態ベクトルξ(k)を作成する。
Then, the state
<逐次推定計算部113の動作例>
逐次推定計算部113は、状態ベクトルξ(k)が作成されると、この状態ベクトルξ(k)を用いて、モデルパラメータθの値を推定する。すなわち、逐次推定計算部113は、サンプリング時刻tk毎に、逐次的にモデルパラメータθの値を推定する。
<Example of Operation of Sequential
When the state vector ξ(k) is generated, the recursive
或るkに関してモデルパラメータθの値を推定する処理について、図4を参照しながら説明する。なお、以下では、サンプリング時刻tkにおけるモデルパラメータθ又はその推定値をθ(k)とも表す。 The process of estimating the value of the model parameter θ with respect to a certain k will be explained with reference to FIG. Note that, hereinafter, the model parameter θ or its estimated value at the sampling time tk will also be expressed as θ(k).
ステップS101:まず、逐次推定計算部113は、モデルパラメータθと共分散行列Pとを初期化するか否かを判定する。ここで、初期化すると判定される場合としては、例えば、初回計算時(つまり、k=0のとき)、ユーザ等により初期化指示が行われたとき等が挙げられる。
Step S101: First, the sequential
モデルパラメータθと共分散行列Pとを初期化すると判定した場合(ステップS101でYES)、逐次推定計算部113は、ステップS102に進む。一方で、モデルパラメータθと共分散行列Pとを初期化すると判定しなかった場合(ステップS101でNO)、逐次推定計算部113は、ステップS103に進む。
If it is determined that the model parameters θ and the covariance matrix P are to be initialized (YES in step S101), the sequential
ステップS102:逐次推定計算部113は、サンプリング時刻tkのインデックスkをk=0に初期化すると共に、θ(0)=θ0及びP(0)=Iと初期化する。ここで、θ0は予め設定されたモデルパラメータの初期値、Iは予め設定された任意の行列(例えば、単位行列等)である。
Step S102: The recursive
ステップS103:逐次推定計算部113は、安定化モデルパラメータ計算部114を呼び出して安定化モデルパラメータθS(k-1)を計算する。なお、安定化モデルパラメータθS(k-1)の計算方法については後述する。
Step S103: The sequential
ステップS104:逐次推定計算部113は、安定化モデルパラメータθS(k-1)をモデルパラメータ推定値θ(k-1)に設定する。すなわち、逐次推定計算部113は、θ(k-1)←θS(k-1)とする。これは、モデルパラメータ推定値θ(k-1)を安定化モデルパラメータθS(k-1)に変換している、ということもできる。
Step S104: The sequential
ステップS105:逐次推定計算部113は、モデルパラメータ推定値θ(k-1)と状態ベクトルξ(k)と制御量y(k)とを用いて、予測誤差ε(k)を計算する。予測誤差ε(k)は、例えば、ε(k)=y(k)-ξ(k)Τθ(k-1)により計算される。なお、Τは転置を表す。
Step S105: The sequential
ステップS106:逐次推定計算部113は、共分散行列P(k-1)を以下により更新して共分散行列P(k)を得る。
Step S106: The sequential
ステップS107:そして、逐次推定計算部113は、モデルパラメータ推定値θ(k-1)を以下により更新してモデルパラメータ推定値θ(k)を得る。
Step S107: Then, the sequential
<安定化モデルパラメータ計算部114の動作例>
安定化モデルパラメータ計算部114は、逐次推定計算部113から呼び出されると、安定化モデルパラメータθS(k-1)を計算する。すなわち、安定化モデルパラメータ計算部114は、サンプリング時刻tkで逐次推定計算部113から呼び出される毎に、モデルパラメータ推定値θ(k-1)から安定化モデルパラメータθS(k-1)を計算する。
<Example of Operation of Stabilization Model
The stabilization model
或るサンプリング時刻tkで逐次推定計算部113から呼び出されたときに、安定化モデルパラメータ計算部114が安定化モデルパラメータを計算する処理について、図5を参照しながら説明する。ただし、以下では、一例として、プラント応答モデルとしては以下のARMAモデルを想定する。
The process by which the stabilized model
y(k)=a1y(k-1)+・・・+aNy(k-N)+b0u(k)+・・・+bMu(k-M)
また、このとき、モデルパラメータ推定値θ(k-1)は以下で表される。
y(k)=a 1 y(k-1)+...+a N y(k-N)+b 0 u(k)+...+b M u(k-M)
Further, at this time, the model parameter estimated value θ(k-1) is expressed as follows.
θ(k-1)=[a1,・・・,aN,b0,・・・,bM]Τ
ここで、Nは制御量yに関する次数、Mは操作量uに関する次数である。すなわち、上記のプラント応答モデルは(N,M)次のARMAモデルである。なお、a1,・・・,aNは自己回帰係数、Nは自己回帰係数の次数とも呼ばれる。また、b0,・・・,bMは移動平均係数、Mは移動平均係数の次数とも呼ばれる。
θ(k−1)=[a 1 , . . . , a N , b 0 , . . . , b M ] T
Here, N is the order of the controlled variable y, and M is the order of the manipulated variable u. That is, the above plant response model is an (N, M)-th order ARMA model. Note that a 1 , ..., a N are also called autoregressive coefficients, and N is also called the order of the autoregressive coefficients. Also, b 0 , ..., b M are also called moving average coefficients, and M is also called the order of the moving average coefficients.
ステップS201:まず、安定化モデルパラメータ計算部114は、重心化パラメータθcを計算する。安定化モデルパラメータ計算部114は、自己回帰係数a1,・・・,aNの重心acと移動平均係数b0,・・・,bMの重心bcとをそれぞれ計算した上で、モデルパラメータ推定値θ(k-1)の各自己回帰係数を重心ac、各移動平均係数を重心bcにそれぞれ置き換えたものを重心化パラメータθcとする。すなわち、重心化パラメータθcは以下で表される。
Step S201: First, the stabilization model
θc=[ac,・・・,ac,bc,・・・,bc]Τ
ここで、ac=(a1+・・・+aN)/N、bc=(b0,・・・,bM)/(M+1)である。
θ c = [ ac , ..., a c , b c , ..., b c ] T
Here, a c =(a 1 +...+a N )/N, and b c =(b 0 ,..., b M )/(M+1).
ステップS202:次に、安定化モデルパラメータ計算部114は、安定化モデルパラメータの候補となるパラメータ(以下、候補パラメータともいう。)計算する。安定化モデルパラメータ計算部114は、以下により候補パラメータθrを計算する。
Step S202: Next, the stabilized model
θr=(1-β)×θ(k-1)+β×θc
ここで、βは0<β≦1を取るブレンド係数であり、予め設定された値である。ブレンド係数βが1に近いほど候補パラメータθrは重心化パラメータθcに近づき、0に近いほど候補パラメータθrはモデルパラメータ推定値θ(k-1)に近づく。
θr = (1-β) × θ(k-1) + β × θc
Here, β is a blending coefficient that satisfies 0<β≦1 and is a preset value. As the blending coefficient β approaches 1, the candidate parameter θr approaches the centroid parameter θc , and as the blending coefficient β approaches 0, the candidate parameter θr approaches the model parameter estimate θ(k−1).
ステップS203:次に、安定化モデルパラメータ計算部114は、2種類の予測誤差を計算する。すなわち、安定化モデルパラメータ計算部114は、モデルパラメータ推定値θ(k-1)を用いたときの予測誤差e0と、候補パラメータθrを用いたときの予測誤差erとをそれぞれ以下により計算する。
Step S203: Next, the stabilized model
e0=y0-ξ(k)Τθ(k-1)
er=y0-ξ(k)Τθr
ここで、y0は制御量の現在値、すなわちy0=y(k)である。また、ξ(k)は状態ベクトルである。すなわち、e0は、モデルパラメータ推定値θ(k-1)を用いたときの制御量予測値と、制御量現在値との誤差を表している。同様に、erは、候補パラメータθrを用いたときの制御量予測値と、制御量現在値との誤差を表している。
e 0 =y 0 -ξ(k) Τ θ(k-1)
e r =y 0 -ξ(k) Τ θ r
Here, y 0 is the current value of the controlled variable, ie, y 0 =y(k). Moreover, ξ(k) is a state vector. That is, e 0 represents the error between the predicted value of the controlled variable when using the model parameter estimated value θ(k-1) and the current controlled variable value. Similarly, e r represents the error between the predicted value of the control amount when using the candidate parameter θ r and the current value of the control amount.
ステップS204:次に、安定化モデルパラメータ計算部114は、予測誤差e0の絶対値と予測誤差erの絶対値とを比較する。安定化モデルパラメータ計算部114は、例えば、予測誤差e0の絶対値と予測誤差erの絶対値とを比較し、以下を満たすか否かを判定する。
Step S204: Next, the stabilization model
|er|<α|e0|
ここで、αは0<α≦1を取る比較係数であり、予め設定された値である。比較係数αは、候補パラメータθrを用いたときの予測誤差erをどの程度厳しく評価するかを制御しており、0に近いほど、予測誤差e0に対して予測誤差erの改善が顕著であることが要求される。
|e r |<α|e 0 |
Here, α is a comparison coefficient satisfying 0<α≦1, and is a preset value. The comparison coefficient α controls how strictly the prediction error e r is evaluated when using the candidate parameter θ r , and the closer it is to 0, the better the prediction error e r is improved relative to the prediction error e 0. It is required to be conspicuous.
上記の比較の結果、|er|<α|e0|を満たす場合、安定化モデルパラメータ計算部114は、ステップS205に進む。一方で、|er|<α|e0|を満たさない場合(つまり、|er|≧α|e0|を満たす場合)、安定化モデルパラメータ計算部114は、ステップS206に進む。
As a result of the above comparison, if |e r |<α|e 0 | is satisfied, the stabilization model
ステップS205:|er|<α|e0|を満たす場合、安定化モデルパラメータ計算部114は、候補パラメータθrを安定化モデルパラメータθS(k-1)とする。すなわち、安定化モデルパラメータ計算部114は、θS(k-1)←θrとする。
Step S205: If |e r |<α|e 0 | is satisfied, the stabilized model
ステップS206:|er|<α|e0|を満たさない場合、安定化モデルパラメータ計算部114は、モデルパラメータ推定値θ(k-1)を安定化モデルパラメータθS(k-1)とする。すなわち、安定化モデルパラメータ計算部114は、θS(k-1)←θ(k-1)とする。
Step S206: If |e r |<α|e 0 | is not satisfied, the stabilization model
なお、上記では、一例として、外乱を考慮していないが、外乱が存在する場合であっても同様に適用できることはいうまでもない。外乱を考慮する場合、例えば、ARMAモデルには外乱に関する移動平均モデルが追加され、その結果、モデルパラメータには外乱に関する自己回帰係数が追加される。 Note that although the above description does not take into account disturbances as an example, it goes without saying that the same application is possible even when disturbances exist. When considering disturbance, for example, a moving average model regarding the disturbance is added to the ARMA model, and as a result, an autoregressive coefficient regarding the disturbance is added to the model parameters.
また、上記では、一例として、自己回帰係数a1,・・・,aNの平均値をac、移動平均係数b0,・・・,bMの平均値をbcとしたが、これに限られるものではない。acとして、自己回帰係数a1,・・・,aNから計算可能な任意の統計量を用いることが可能である。同様に、bcとして、移動平均係数b0,・・・,bMから計算可能な任意の統計量を用いることが可能である。より一般には、次元毎に、その次元方向の係数から計算可能な任意の統計量を用いることが可能である。このような統計量としては、平均の他、例えば、加重平均、中央値、最頻値等といったものが挙げられる。 In the above, as an example, the average value of the autoregressive coefficients a 1 , ..., a N is a c , and the average value of the moving average coefficients b 0 , ..., b M is b c , but this is not limited to the above. Any statistic that can be calculated from the autoregressive coefficients a 1 , ..., a N can be used as a c . Similarly, any statistic that can be calculated from the moving average coefficients b 0 , ..., b M can be used as b c . More generally, any statistic that can be calculated from the coefficients in the dimensional direction can be used for each dimension. Examples of such statistics include the average, as well as the weighted average, the median, the mode, and the like.
<ステップ応答計算部102の動作例>
以下、ステップ応答の推定に用いられるモデルパラメータ推定値(つまり、プラント応答関数Sθに設定されるモデルパラメータ推定値)θ=θ(k)のことを「モデルパラメータ設定値θ」ともいう。また、以下では、簡単のため、ステップ信号として単位ステップ信号を想定する。
<Example of Operation of Step
Hereinafter, the model parameter estimate value θ=θ(k) used for estimating the step response (i.e., the model parameter estimate value set in the plant response function Sθ ) is also referred to as the “model parameter set value θ.” In addition, in the following, for simplicity, a unit step signal is assumed as the step signal.
ステップ応答計算部102は、図6に示すように、モデルパラメータ設定値θと時刻tとが与えられると、初期時刻0から時間t経過後の時刻tにおける単位ステップ応答Sθ(t)を計算する。なお、単位ステップ応答とは、操作量uとして単位ステップ信号を印加した場合における応答(つまり、制御対象プラントのプラント応答モデルの出力)のことである。
6 , when a model parameter setting value θ and a time t are given, the step
或る時刻tと或るモデルパラメータ設定値θとが与えられたときに、ステップ応答Sθ(t)を計算する処理について、図7を参照しながら説明する。ただし、以下では、一例として、プラント応答モデルとしては以下のARMAモデルを想定する。 A process for calculating a step response S θ (t) when a certain time t and a certain model parameter setting value θ are given will be described with reference to Fig. 7. In the following, the following ARMA model is assumed as a plant response model as an example.
y(k')=a1y(k'-1)+・・・+aNy(k'-N)+b0u(k')+・・・+bMu(k'-M)
また、このとき、インデックスk'における状態ベクトルφ(k')を以下で表すものとする。
y(k') = a 1 y(k'-1) + ... + a N y(k'-N) + b 0 u(k') + ... + b M u(k'-M)
In addition, in this case, the state vector φ(k′) at index k′ is expressed as follows.
φ(k')=[y(k'-1),y(k'-2),・・・,y(k'-N),u(k'),u(k'-1),・・・,u(k'-M)]Τ
なお、インデックスk'は、サンプリング時刻tkのインデックスkと同じ値を取り得る変数であるが、本処理の中でのみ利用され、インデックスkとは独立に値が更新されることに留意されたい。
φ(k') = [y(k'-1), y(k'-2), ..., y(k'-N), u(k'), u(k'-1), . ..., u(k'-M)] T
Note that index k' is a variable that can take the same value as index k at sampling time tk , but it is used only in this process and its value is updated independently of index k. .
ステップS301:ステップ応答計算部102は、本処理の中でのみ利用する時刻を表すインデックスをτとして、τ=0、k'=0と初期化すると共に、状態ベクトルφ(0)を以下のように初期化する。
Step S301: The step
φ(0)=[0,0,・・・,0,1,0,・・・,0]Τ
すなわち、u(0)のみ1、それ以外の要素は0と状態ベクトルφ(0)を初期化する。これは、単位ステップ応答を模擬する際の初期値を表している。
φ(0)=[0,0,...,0,1,0,...,0] Τ
That is, the state vector φ(0) is initialized so that only u(0) is 1 and the other elements are 0. This represents the initial value when simulating a unit step response.
ステップS302:次に、ステップ応答計算部102は、y(k')=φ(k')Τθにより制御量予測値y(k')を計算する。
Step S302: Next, the step
ステップS303:次に、ステップ応答計算部102は、制御量予測値y(k')を用いて、状態ベクトルφ(k')を、次のインデックスk'+1における状態ベクトルφ(k'+1)に更新する。このとき、ステップ応答計算部102は、状態ベクトルφ(k'+1)のy(k')には上記のステップS302で計算した制御量予測値y(k')を設定し、y(k'-N+1)~y(k'-1)には状態ベクトルφ(k')と同じ値を設定する。また、状態ベクトルφ(k'+1)のu(k'+1)には1を設定し、u(k'-M+1)~u(k')には状態ベクトルφ(k')と同じ値を設定する。
Step S303: Next, the step
ステップS304:次に、ステップ応答計算部102は、時刻τをτ+Δに更新すると共に、インデックスk'をk'+1に更新する。
Step S304: Next, the step
ステップS305:次に、ステップ応答計算部102は、τ≧tであるか否かを判定する。そして、τ≧tであると判定されなかった場合(ステップS305でNO)、ステップ応答計算部102は、ステップS302に戻る。これにより、τ≧tとなるまで、ステップS302~ステップS304が繰り返し実行される。
Step S305: Next, the step
一方で、τ≧tであると判定された場合(ステップS305でYES)、ステップ応答計算部102は、処理を終了する。これにより、最終的に計算された制御量予測値y(k')が単位ステップ応答Sθ(t)として得られる(つまり、Sθ(t)=y(k)が、プラント応答モデルの単位ステップ応答として得られる。)。
On the other hand, if it is determined that τ≧t (YES in step S305), the step
[第二の実施形態]
以下、第二の実施形態について説明する。第二の実施形態では、制御対象のプラントの定常状態を表すオフセット項が含まれるモデルパラメータを対象とする。なお、プラントの定常状態を表すオフセット項が含まれるモデルパラメータとそれに関連する技術については、例えば、上記の特許文献3等を参照されたい。
[Second embodiment]
The second embodiment will be described below. In the second embodiment, model parameters including an offset term representing the steady state of a plant to be controlled are targeted. Note that for model parameters including an offset term representing the steady state of the plant and related techniques, please refer to, for example, the above-mentioned
具体的には、第二の実施形態では、プラント応答モデルとして以下のオフセット項が含まれるARMAXモデル(Autoregressive Moving Average Model with Exogenous Variables)を想定する。 Specifically, in the second embodiment, an ARMAX model (Autoregressive Moving Average Model with Exogenous Variables) including the following offset terms is assumed as a plant response model.
y(k)=a1y(k-1)+・・・+aNy(k-N)+b0u(k)+・・・+bMu(k-M)+c0
ここで、c0がオフセット項である。
y(k) = a 1 y(k-1) + ... + a N y(k-N) + b 0 u(k) + ... + b M u(k-M) + c 0
where c 0 is the offset term.
このとき、モデルパラメータ(及びその推定値)θは、以下で表されるモデルパラメータに拡大することができる。 At this time, the model parameter (and its estimated value) θ can be expanded to the model parameter expressed below.
θe=[a1,・・・,aN,b0,・・・,bM,c0]Τ
このモデルパラメータ(及びその推定値)θeを「拡大モデルパラメータ」とも呼ぶ。
θ e = [a 1 , ..., a N , b 0 , ..., b M , c 0 ] T
This model parameter (and its estimated value) θ e is also called the “augmented model parameter”.
なお、第二の実施形態では、主に、第一の実施形態と相違点について説明し、第一の実施形態と同様としてよい箇所についてはその説明を省略する。 Note that in the second embodiment, the differences from the first embodiment will be mainly described, and descriptions of the parts that may be the same as in the first embodiment will be omitted.
<プラント応答推定装置10の機能構成例>
本実施形態に係るプラント応答推定装置10の機能構成例を図8に示す。図8に示すように、本実施形態に係るプラント応答推定装置10は、パラメータ変換部115がモデルパラメータ推定部101に含まれている点が第一の実施形態と異なる。
<Example of functional configuration of plant
FIG. 8 shows an example of the functional configuration of the plant
パラメータ変換部115は、拡大モデルパラメータθeをモデルパラメータθに変換する。
The
<状態ベクトル変換部112の動作例>
状態ベクトル変換部112は、状態ベクトルを拡大した拡大状態ベクトルと呼ぶベクトルを作成する点が第一の実施形態と異なる。その他の点は第一の実施形態と同様としてよい。
<Example of operation of state
The state
拡大状態ベクトルとは、状態ベクトルξ(k)をオフセット項に対応する要素を加えて拡大したものであり、以下で表される。 The expanded state vector is obtained by expanding the state vector ξ(k) by adding an element corresponding to the offset term, and is expressed as follows.
<逐次推定計算部113の動作例>
逐次推定計算部113は、モデルパラメータθの代わりに拡大モデルパラメータθeを用いる点、状態ベクトルξ(k)の代わりに拡大状態ベクトルξe(k)を用いる点が第一の実施形態と異なる。その他の点は第一の実施形態と同様としてよい。
<Example of Operation of Sequential
The recursive
<安定化モデルパラメータ計算部114の動作例>
安定化モデルパラメータ計算部114は、モデルパラメータθの代わりに拡大モデルパラメータθeを用いる点、状態ベクトルξ(k)の代わりに拡大状態ベクトルξe(k)を用いる点が第一の実施形態と異なる。その他の点は第一の実施形態と同様としてよい。
<Example of operation of stabilization model
The stabilized model
なお、上記のステップS201で重心化パラメータθcを計算する際に、オフセット項は置き換えの対象外であることに留意されたい。すなわち、モデルパラメータ推定値θ(k-1)がθe(k-1)=[a1,・・・,aN,b0,・・・,bM,c0]Τである場合、重心化パラメータθcはθc=[ac,・・・,ac,bc,・・・,bc,c0]Τとなる。 It should be noted that the offset term is not subject to replacement when calculating the centroidization parameter θ c in step S201 above. That is, when the model parameter estimate θ(k-1) is θ e (k-1)=[a 1 ,..., a N , b 0 ,..., b M , c 0 ] T , The centroiding parameter θ c becomes θ c =[ ac , . . . , ac , b c , . . . , b c , c 0 ] T.
<パラメータ変換部115の動作例>
以下、一例として、拡大モデルパラメータ推定値θe(k)は以下で表されるものとする。
<Example of Operation of
Hereinafter, as an example, the augmented model parameter estimate θ e (k) is expressed as follows:
θe(k)=[θY(k),θU(k),θV(k),θC(k)]Τ
ここで、θY(k)はサンプリング時刻tkにおけるARMAXモデルの制御量yに関する項の係数を要素とするN次元ベクトル、θU(k)は操作量uに関する項の係数を要素とするM次元ベクトル、θV(k)は外乱vに関する項の係数を要素とするL次元ベクトル、θC(k)はオフセット項を表すスカラー値である。
θe (k)=[ θY (k), θU (k), θV (k), θC (k)] T
Here, θY (k) is an N-dimensional vector whose elements are the coefficients of the term related to the controlled variable y of the ARMAX model at sampling time tk , θU (k) is an M-dimensional vector whose elements are the coefficients of the term related to the manipulated variable u, θV (k) is an L-dimensional vector whose elements are the coefficients of the term related to the disturbance v, and θC (k) is a scalar value representing the offset term.
このとき、パラメータ変換部115は、拡大モデルパラメータ推定値θe(k)が得られると、図9に示すように、θC(k)を除いたベクトルに変換することで、モデルパラメータ推定値θ(k)を得る。これにより、以下のモデルパラメータ推定値が得られる。
In this case, when the
θ(k)=[θY(k),θU(k),θV(k)]Τ
このように、制御対象プラントのステップ応答を推定する際には、拡大モデルパラメータ推定値θe(k)からθC(k)を除くことで、モデルパラメータ推定値θ(k)が得られる。この理由は、上記の特許文献1に記載されている通り、制御量yの変化量に対してオフセット項は影響しないためである。
θ(k)=[θ Y (k), θ U (k), θ V (k)] Τ
In this way, when estimating the step response of the controlled plant, the model parameter estimate θ(k) is obtained by removing θ C (k) from the expanded model parameter estimate θ e (k). The reason for this is that, as described in
[実施例]
以下、一実施例について説明する。本実施例では、第二の実施形態に係るプラント応答推定装置10によりモデルパラメータとステップ応答を推定した。
[Example]
An example will be described below. In this example, model parameters and a step response are estimated by the plant
制御対象プラントのプラント応答モデルとしては、以下のARMAXモデルを採用した。 The following ARMAX model was adopted as the plant response model for the plant to be controlled.
y(k)=a1y(k-1)+a2y(k-2)+b0u(k)+b1u(k-1)+b2u(k-2)+c0
すなわち、制御量yに関する2次元の要素と、操作量uに関する3次元の要素と、オフセット項とを有するARMAXモデルを採用した。
y(k)=a 1 y(k-1)+a 2 y(k-2)+b 0 u(k)+b 1 u(k-1)+b 2 u(k-2)+c 0
That is, an ARMAX model having a two-dimensional element related to the controlled amount y, a three-dimensional element related to the manipulated variable u, and an offset term was adopted.
このとき、拡大状態ベクトルξe(k)は以下で表される。 In this case, the extended state vector ξ e (k) is expressed as follows.
ξe(k)=[y(k-1),y(k-2),u(k),u(k-1),u(k-2),1]Τ
また、拡大モデルパラメータθeは以下で表される。
ξ e (k) = [y (k-1), y (k-2), u (k), u (k-1), u (k-2), 1] T
Moreover, the enlarged model parameter θ e is expressed as follows.
θe(k)=[a1(k),a2(k),b0(k),b1(k),b2(k),c0(k)]Τ
制御量の予測値(以下、これをyestとも表す。)は以下で表される。
θ e (k) = [a 1 (k), a 2 (k), b 0 (k), b 1 (k), b 2 (k), c 0 (k)] Τ
The predicted value of the control amount (hereinafter also referred to as y est ) is expressed as follows.
yest(k)=ξe(k)Τθe(k-1)
ここで、本実施例における制御対象プラントのステップ応答は、図10に示すようなものであるとする。なお、rは目標値、yは制御量、uは操作量、vは外乱を表す。また、本実施例では、再サンプリング周期D=2、忘却係数λ=0.999、比較係数α=0.9とした。
y est (k)=ξ e (k) Τ θ e (k-1)
Here, it is assumed that the step response of the plant to be controlled in this embodiment is as shown in FIG. Note that r represents the target value, y represents the controlled amount, u represents the manipulated amount, and v represents the disturbance. Further, in this embodiment, the resampling period D=2, the forgetting coefficient λ=0.999, and the comparison coefficient α=0.9.
このとき、本実施例では、2つのデータ(後述するデータ1、データ2)をそれぞれ用いて、モデルパラメータとステップ応答を推定した。
In this embodiment, two pieces of data (
・データ1を用いた場合
制御対象プラントをPID制御で制御したデータ1を図11に示す。このとき、安定化モデルパラメータ計算部114を利用しないでモデルパラメータを推定した結果とモデル予測制御の制御結果とをそれぞれ図12及び図13に示す。また、第二の実施形態に係るプラント応答推定装置10によりモデルパラメータを推定した結果とモデル予測制御の制御結果とをそれぞれ図14及び図15に示す。図12及び図14では、モデルパラメータのa1、a2、b0、b1、b2を時系列で示している。図12(安定化モデルパラメータ計算部114なし)と図14(安定化モデルパラメータ計算部114あり)とを比較すると、自己回帰係数a1、a2は安定化モデルパラメータ計算部114ありの方が重心に近く、移動平均係数b0、b1、b2は安定化モデルパラメータ計算部114ありの方が重心に近いことがわかる。一方で、図13(安定化モデルパラメータ計算部114なし)と図15(安定化モデルパラメータ計算部114あり)とを比較すると、大きな違いはなく、安定化モデルパラメータ計算部114によって制御性能の劣化は生じていないことがわかる。
When
・データ2を用いた場合
制御対象プラントをPID制御で制御したデータ2を図16に示す。データ2では、制御量yに観測ノイズが加わっており、その結果、操作量uにもノイズが加わっている。なお、ノイズの影響以外はデータ1とデータ2とで条件は同一である。このとき、安定化モデルパラメータ計算部114を利用しないでモデルパラメータを推定した結果とモデル予測制御の制御結果とをそれぞれ図17及び図18に示す。また、第二の実施形態に係るプラント応答推定装置10によりモデルパラメータを推定した結果とモデル予測制御の制御結果とをそれぞれ図19及び図20に示す。図17及び図19では、モデルパラメータのa1、a2、b0、b1、b2を時系列で示している。図17(安定化モデルパラメータ計算部114なし)と図19(安定化モデルパラメータ計算部114あり)とを比較すると、自己回帰係数a1、a2は安定化モデルパラメータ計算部114ありの方が重心に近く、移動平均係数b0、b1、b2は安定化モデルパラメータ計算部114ありの方が重心に近いことがわかる。また、安定化モデルパラメータ計算部114なしの場合(図17)では、自己回帰係数a1、a2が徐々に互いに離れた値となっていくが、安定化モデルパラメータ計算部114ありの場合(図19)では、その動きが抑制できていることがわかる。すなわち、ノイズのあるデータ2ではノイズによって自己回帰係数及び移動平均係数が重心から離れた値で推定されるが、安定化モデルパラメータ計算部114によって重心から離れる動きが抑制できている。
- When using
一方で、図18(安定化モデルパラメータ計算部114なし)と図20(安定化モデルパラメータ計算部114あり)とを比較すると、図18では操作量uが比較的鋭く変化している。これに対して、図20では、上記のような鋭い変化は少なく、ノイズのないデータ1での制御結果(図15)の応答波形と比較的良く似た結果となっている。したがって、安定化モデルパラメータ計算部114により、ノイズの影響を抑制しており、モデル予測制御の安定化に寄与しているといえる。
On the other hand, comparing FIG. 18 (without the stabilization model parameter calculation unit 114) with FIG. 20 (with the stabilization model parameter calculation unit 114), the manipulated variable u changes relatively sharply in FIG. 18. In contrast, in FIG. 20, there are fewer sharp changes like those described above, and the result is relatively similar to the response waveform of the control result for noise-free data 1 (FIG. 15). Therefore, it can be said that the stabilization model
[まとめ]
以上のように、第一及び第二の実施形態に係るプラント応答推定装置10は、操業中のプラントの運転データからプラント応答モデルのモデルパラメータを安定的に推定することができる。すなわち、重心化パラメータによってモデルパラメータの推定には重心方向に推定されるような力が働くため、例えば、図17のように、a1、a2が徐々に互いに離れた値となっていく、というような動き(つまり、モデルパラメータが発散するような動き)を抑制することが可能となり、その結果、モデルパラメータを安定的に推定することが可能となる。これにより、真値に収束するまでの途中の段階でもモデルパラメータ推定値の信頼性が担保され、上記の1つ目の課題が解決される。このため、第一及び第二の実施形態に係るプラント応答推定装置10によれば、推定したモデルパラメータを制御(例えば、モデル予測制御)に容易に用いることが可能となる。
[summary]
As described above, the plant
また、第一及び第二の実施形態に係るプラント応答推定装置10では、モデルパラメータの次数が未知であり、必要以上に高い次数が設定された場合であっても、重心方向に推定されるような力が働くため、重心に近いモデルパラメータが推定されやすくなり、多重共線性等の高次数に由来する問題を回避することができるようになる。これにより、上記の2つ目の課題が解決される。
In addition, in the plant
更に、第一及び第二の実施形態に係るプラント応答推定装置10では、安定化モデルパラメータ計算部114にて現状のモデルパラメータを用いた場合の予測誤差と候補パラメータを用いた場合の予測誤差とを比較し、予測精度が改善する場合にのみ候補パラメータを安定化モデルパラメータとして採用する。このため、予測精度の改善も期待できる。
Furthermore, in the plant
加えて、第一及び第二の実施形態に係るプラント応答推定装置10では、観測値に対するノイズに対してもモデルパラメータ推定値が安定するため、ノイズの影響を抑制した推定が可能となる。
In addition, in the plant
本発明は、具体的に開示された上記の実施形態に限定されるものではなく、特許請求の範囲の記載から逸脱することなく、種々の変形や変更、既知の技術との組み合わせ等が可能である。 The present invention is not limited to the above-described specifically disclosed embodiments, and various modifications and changes, combinations with known techniques, etc. can be made without departing from the scope of the claims. be.
10 プラント応答推定装置
11 入力装置
12 表示装置
13 外部I/F
13a 記録媒体
14 通信I/F
15 プロセッサ
16 メモリ装置
17 バス
101 モデルパラメータ推定部
102 ステップ応答計算部
111 バッファ部
112 状態ベクトル変換部
113 逐次推定計算部
114 安定化モデルパラメータ計算部
115 パラメータ変換部
10 Plant
13a Recording medium 14 Communication I/F
15
Claims (8)
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算するように構成されている安定化計算部と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算するように構成されている応答計算部と、
を有し、
前記推定計算部は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算部によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算するように構成されている、プラント応答推定装置。 an estimation calculation unit configured to sequentially calculate, based on operation data of a controlled plant, estimates of model parameters, which are parameters of a function representing a plant response model of the controlled plant;
a stabilization calculation unit configured to calculate stabilized model parameters by stabilizing the estimates of the model parameters;
a response calculation unit configured to calculate a response of the controlled plant using a function representing the plant response model based on the estimated values of the model parameters;
having
The estimation calculation unit is
a plant response estimation device configured to calculate an estimate of the model parameter based on the operation data and a stabilized model parameter obtained by stabilizing a previously calculated estimate of the model parameter by the stabilization calculation unit.
前記モデルパラメータの推定値を表すベクトルの各要素から所定の統計量を計算し、前記統計量により前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算するように構成されている、請求項1に記載のプラント応答推定装置。 The stabilization calculation unit
2. The plant response estimation device according to claim 1, configured to calculate a predetermined statistic from each element of a vector representing the estimated value of the model parameter, and to calculate a stabilized model parameter by stabilizing the estimated value of the model parameter using the statistic.
前記安定化計算部は、
前記多項式モデルの次元毎に、前記モデルパラメータの推定値を表すベクトルの各要素のうち、前記次元方向の要素から前記統計量を計算し、
前記次元毎に、前記次元方向の要素を前記統計量で置換した第1のパラメータを計算し、
前記第1のパラメータと前記モデルパラメータとを所定の比率で混合した第2のパラメータを計算し、
前記第2のパラメータを用いて前記制御対象プラントの応答を予測したときの第1の予測誤差が、前記モデルパラメータを用いて前記制御対象プラントの応答を予測したときの第2の予測誤差よりも所定以上改善している場合は前記第2のパラメータを前記安定化モデルパラメータとし、前記第1の予測誤差が前記第2の予測誤差よりも所定以上改善していない場合は前記モデルパラメータを前記安定化モデルパラメータとする、ように構成されている、請求項2に記載のプラント応答推定装置。 the plant response model is a polynomial model;
The stabilization calculation unit
calculating the statistics from elements in the dimension direction among the elements of a vector representing the estimated values of the model parameters for each dimension of the polynomial model;
calculating a first parameter by replacing an element in the dimension direction with the statistic for each dimension;
Calculating second parameters by mixing the first parameters and the model parameters in a predetermined ratio;
3. The plant response estimation device according to claim 2, configured to: when a first prediction error when predicting a response of the controlled plant using the second parameter is improved by a predetermined amount or more than a second prediction error when predicting a response of the controlled plant using the model parameter, the second parameter is set as the stabilization model parameter; and when the first prediction error is not improved by the predetermined amount or more than the second prediction error, the model parameter is set as the stabilization model parameter.
前記モデルパラメータは、前記ARMAモデル又はARMAXモデルの係数を要素とするベクトルで表される、請求項3に記載のプラント応答推定装置。 the polynomial model is an ARMA model or an ARMAX model;
4. The plant response estimation device according to claim 3, wherein the model parameters are expressed as vectors having coefficients of the ARMA model or the ARMAX model as elements.
前記ARMAモデル又はARMAXモデルの自己回帰係数の平均値を表す第1の重心と、前記ARMAモデル又はARMAXモデルの移動平均係数の平均値を表す第2の重心とを前記統計量として計算し、
前記自己回帰係数を前記第1の重心、前記移動平均係数を前記第2の重心でそれぞれ置換することで、前記第2のパラメータを計算する、ように構成されている、請求項4に記載のプラント応答推定装置。 The stabilization calculation section includes:
Calculating a first centroid representing the average value of the autoregressive coefficient of the ARMA model or ARMAX model and a second centroid representing the average value of the moving average coefficient of the ARMA model or ARMAX model as the statistic,
5. The second parameter is calculated by replacing the autoregressive coefficient with the first center of gravity and the moving average coefficient with the second center of gravity, respectively. Plant response estimation device.
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算する安定化計算手順と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算する応答計算手順と、
をコンピュータが実行し、
前記推定計算手順は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算手順によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算する、プラント応答推定方法。 an estimation calculation procedure for sequentially calculating estimated values of model parameters, which are parameters of a function representing a plant response model of a controlled plant, based on operation data of the controlled plant;
a stabilization calculation step for calculating stabilized model parameters by stabilizing the estimated values of the model parameters;
a response calculation step of calculating a response of the controlled plant using a function representing the plant response model based on the estimated values of the model parameters;
The computer executes
The estimation calculation procedure includes:
the plant response estimation method further comprising: calculating an estimated value of the model parameter based on the operating data and a stabilized model parameter obtained by stabilizing an estimated value of a previously calculated model parameter through the stabilization calculation procedure;
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算する安定化計算手順と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算する応答計算手順と、
をコンピュータに実行させ、
前記推定計算手順は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算手順によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算する、プログラム。 an estimation calculation procedure for sequentially calculating estimated values of model parameters, which are parameters of a function representing a plant response model of a controlled plant, based on operation data of the controlled plant;
a stabilization calculation step for calculating stabilized model parameters by stabilizing the estimated values of the model parameters;
a response calculation step of calculating a response of the controlled plant using a function representing the plant response model based on the estimated values of the model parameters;
Run the following on your computer:
The estimation calculation procedure includes:
a program for calculating an estimated value of the model parameter based on the operating data and a stabilized model parameter obtained by stabilizing an estimated value of a previously calculated model parameter through the stabilization calculation procedure;
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2022150202A JP7231102B1 (en) | 2022-09-21 | 2022-09-21 | PLANT RESPONSE ESTIMATION DEVICE, PLANT RESPONSE ESTIMATION METHOD, AND PROGRAM |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2022150202A JP7231102B1 (en) | 2022-09-21 | 2022-09-21 | PLANT RESPONSE ESTIMATION DEVICE, PLANT RESPONSE ESTIMATION METHOD, AND PROGRAM |
Publications (2)
Publication Number | Publication Date |
---|---|
JP7231102B1 JP7231102B1 (en) | 2023-03-01 |
JP2024044583A true JP2024044583A (en) | 2024-04-02 |
Family
ID=85380661
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2022150202A Active JP7231102B1 (en) | 2022-09-21 | 2022-09-21 | PLANT RESPONSE ESTIMATION DEVICE, PLANT RESPONSE ESTIMATION METHOD, AND PROGRAM |
Country Status (1)
Country | Link |
---|---|
JP (1) | JP7231102B1 (en) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP7428288B1 (en) | 2023-04-25 | 2024-02-06 | 富士電機株式会社 | Plant response estimation device, plant response estimation method, and program |
Family Cites Families (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2000072096A1 (en) * | 1999-05-25 | 2000-11-30 | Siemens Aktiengesellschaft | Method, arrangement and computer program for designing a technical system |
JP2008282159A (en) * | 2007-05-09 | 2008-11-20 | Toyota Motor Corp | Multi-input/multi-output-system control device |
JP4815391B2 (en) * | 2007-05-15 | 2011-11-16 | 株式会社神戸製鋼所 | Model parameter estimation calculation apparatus and method, model parameter estimation calculation processing program, and recording medium recording the same |
JP7047966B1 (en) * | 2021-10-01 | 2022-04-05 | 富士電機株式会社 | Plant response estimator, plant response estimation method, and program |
-
2022
- 2022-09-21 JP JP2022150202A patent/JP7231102B1/en active Active
Also Published As
Publication number | Publication date |
---|---|
JP7231102B1 (en) | 2023-03-01 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP6513015B2 (en) | Method for controlling machine operation, and control system for repetitively controlling machine operation | |
JP6359182B2 (en) | Method and system for controlling the operation of a machine | |
JP6901037B1 (en) | Control devices, control methods and programs | |
JP2004503000A (en) | Multivariate matrix process control | |
JP7044077B2 (en) | Model estimation system, method and program | |
JP6299759B2 (en) | Prediction function creation device, prediction function creation method, and program | |
JP2023053573A (en) | Plant response estimation device, plant response estimation method and program | |
JP7231102B1 (en) | PLANT RESPONSE ESTIMATION DEVICE, PLANT RESPONSE ESTIMATION METHOD, AND PROGRAM | |
JP7283065B2 (en) | Estimation device, optimization device, estimation method, optimization method, and program | |
JP7077667B2 (en) | Control method, control device and program | |
KR20070099330A (en) | Robot and method for localization of the robot using calculated covariance | |
JP6958808B2 (en) | Policy improvement programs, policy improvement methods, and policy improvement devices | |
JP7115654B1 (en) | Control device, control method and program | |
Woodley et al. | Subspace based direct adaptive ℋ︁∞ control | |
JP7428288B1 (en) | Plant response estimation device, plant response estimation method, and program | |
JPWO2019142728A1 (en) | Controls, control methods and programs | |
US20200356622A1 (en) | Arithmetic processing apparatus, arithmetic processing method, and non-transitory computer-readable storage medium for storing arithmetic processing program | |
JP7275492B2 (en) | Control device, control method and program | |
JP7115656B1 (en) | Control device, control method and program | |
JP7472998B2 (en) | Parameter estimation device, secret parameter estimation system, secure computing device, methods thereof, and programs | |
JP6904473B1 (en) | Model creation support device, model creation support method and program | |
JP7505695B2 (en) | Motor Control Device | |
JP2020179438A (en) | Computing system and machine learning method | |
JP7163977B2 (en) | Estimation device, learning device, method thereof, and program | |
JP2017207987A (en) | Objective variable prediction device, method and program |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20220921 |
|
A871 | Explanation of circumstances concerning accelerated examination |
Free format text: JAPANESE INTERMEDIATE CODE: A871 Effective date: 20220921 |
|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20221129 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20221220 |
|
TRDD | Decision of grant or rejection written | ||
A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20230117 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20230130 |
|
R150 | Certificate of patent or registration of utility model |
Ref document number: 7231102 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |