JPH083743B2 - Method for identifying process parameters in a reactor - Google Patents

Description

【発明の詳細な説明】 ［産業上の利用分野］ 本発明は、化学反応器などのように非線形性の大きい
プロセスにおけるパラメータを、運転時のデータにもと
づき、逐次型最尤推定法によって同定する方法に関す
る。DETAILED DESCRIPTION OF THE INVENTION [Industrial field of application] The present invention identifies parameters in a process having a large non-linearity such as a chemical reactor by a sequential maximum likelihood estimation method based on operating data. Regarding the method.

［従来の技術］ 石油工業、化学工業あるいは製紙工業等のプロセス工
業においては各種のプラント等の制御をプロセス制御に
よって行なっている。[Prior Art] In the process industries such as the petroleum industry, the chemical industry, and the paper manufacturing industry, various plants are controlled by process control.

このうち、非線形性の大きい反応器などの非線形プロ
セスにおいては、反応器設計時の反応速度式，伝熱係
数，触媒量，触媒活性などが、プロセスを実際に運転す
ると変化する。したがって、これらの変化による反応速
度，原料濃度，製品濃度などのパラメータの変化を推定
し、このパラメータ変化に対応した適応制御や異常診断
が行なわれている。Among them, in a non-linear process such as a reactor having a large non-linearity, the reaction rate formula, heat transfer coefficient, catalyst amount, catalyst activity, etc. at the time of designing the reactor change when the process is actually operated. Therefore, changes in parameters such as reaction rate, raw material concentration, and product concentration due to these changes are estimated, and adaptive control and abnormality diagnosis are performed in accordance with these parameter changes.

従来、非線形性の大きいプロセスのパラメータ推定を
行なう方法としては、非線形モデルを線形化して拡張カ
ルマンフィルタ法によってパラメータを推定する方法
や、理論プロセス・モデルにおける複数のパラメータの
うち、いくつかのパラメータを固定しておいて他のパラ
メータを推定する方法などが採用されていた。Conventionally, as a method for parameter estimation of a process having a large non-linearity, a parameter is estimated by the extended Kalman filter method by linearizing a non-linear model, or some parameters among a plurality of parameters in a theoretical process model are fixed. Then, the method of estimating other parameters was adopted.

［解決すべき問題点］ 上述した従来の方法のうち、拡張カルマンフィルタ法
を用いる方法は、非線形モデルを線形化する際の誤差が
大きく、実際のプロセスに適用するのは困難であった。
また、理論プロセス・モデルにおける複数のパラメータ
のうち、いくつかのパラメータを固定しておいて他のパ
ラメータを推定する方法は、それぞれのパラメータを対
応させて決定することが困難であり、この方法も実際の
プロセスに適用することは難しかった。[Problems to be Solved] Among the above-mentioned conventional methods, the method using the extended Kalman filter method has a large error in linearizing the nonlinear model and is difficult to apply to an actual process.
Further, in a method of fixing some parameters among a plurality of parameters in the theoretical process model and estimating other parameters, it is difficult to determine each parameter in correspondence, and this method also It was difficult to apply to the actual process.

なお、特開昭60-150109号公報に、入力情報と出力情
報との関係にもとづいて規定される残差の二乗の荷重和
が最小になるようにプロセス・パラメータを決定する
際、過去のデータに対して重み付けをすることにより過
去の影響を無くするようにしたパラメータの同定方法が
記載されている。It should be noted that when determining process parameters so that the weighted sum of the squares of residuals, which is specified in JP-A-60-150109, based on the relationship between input information and output information, is minimized, past data There is described a method of identifying a parameter in which the influence of the past is eliminated by weighting.

このパラメータ同定方法は非線形性が不明な場合にも
適用できると思われるが、プロセス・パラメータを変化
に対応して必ずしも正確に同定することができないとい
う問題点を有している。This parameter identification method is considered to be applicable even when the nonlinearity is unknown, but it has a problem that it is not always possible to accurately identify process parameters in response to changes.

本発明は上記の問題点にかんがみてなされたもので、
非線形プロセスの理論モデルにもとづいて、逐次型最尤
推定法によりパラメータを正確に同定できるようにした
プロセス・パラメータの同定方法の提供を目的とする。The present invention has been made in view of the above problems,
It is an object of the present invention to provide a process parameter identification method capable of accurately identifying parameters by the iterative maximum likelihood estimation method based on a theoretical model of a nonlinear process.

［問題点の解決手段］ 本発明の反応器におけるプロセス・パラメータの同定
方法は、上記目的を達成するため、 現在の観測値から１サンプル後のプロセス出力を推定
し、１サンプル後における観測値と推定値の誤差が小さ
くなるように頻度因子、熱伝達係数等のプロセス・パラ
メータを修正する逐次型最尤推定法によって、非線形プ
ロセスの理論モデルにおけるパラメータの変化を同定す
る方法としてある。[Means for Solving Problems] In order to achieve the above-mentioned object, the method for identifying process parameters in a reactor according to the present invention estimates the process output after one sample from the current observed value, and estimates the process output after one sample. This is a method for identifying parameter changes in a theoretical model of a non-linear process by a sequential maximum likelihood estimation method that modifies process parameters such as a frequency factor and a heat transfer coefficient so that an error in an estimated value becomes small.

さらに詳述すれば、本発明のプロセス・パラメータの
同定方法は、 イ．反応器の非線形プロセスモデルを、プロセス・パラ
メータに関する線形差分方程式に変換し、 ロ．反応器におけるプロセス・パラメータである頻度因
子，熱伝達係数等についてのスケーリング・ファクタを
求めて推定式を設定し、 ハ．反応器からデータをサンプリングし、このサンプリ
ング周期ごとに、上記推定式により上記頻度因子，熱伝
達係数等を推定し、 ニ．誤差修正ゲインを計算し、 ホ．反応器出力の推定値を計算し、かつ、この推定値と
観測値より推定誤差を計算し、 ヘ．上記推定誤差が収束しているか否かを判断し、 ト．収束していない場合には、上記誤差修正ゲインおよ
び推定誤差にもとづいて推定した頻度因子，熱伝達係数
等を修正するとともに、誤差共分散を修正し、上記ハ．
〜ホ．を繰返して同定を継続し、 チ．収束している場合には同定を終了する 方法としてある。More specifically, the process parameter identification method of the present invention comprises: Convert the nonlinear process model of the reactor into a linear difference equation for the process parameters, b. Set the estimation formulas by obtaining the scaling factors for the process parameters in the reactor, such as frequency factors and heat transfer coefficients, and c. Data is sampled from the reactor, and the frequency factor, heat transfer coefficient, etc. are estimated by the above estimation formula for each sampling period, and d. Calculate the error correction gain, and e. Calculate the estimated value of the reactor output, and the estimated error from this estimated value and the observed value. F. It is judged whether the above estimation error has converged, and If it has not converged, the frequency factor, the heat transfer coefficient, and the like estimated based on the error correction gain and the estimation error are corrected, and the error covariance is corrected.
~ E. Repeat the above to continue the identification, and If it has converged, this is a method of ending the identification.

［実施例］ 以下、本発明の実施例について図面を参照して説明す
る。[Embodiment] An embodiment of the present invention will be described below with reference to the drawings.

第１図は、制御系の構成要素に顕著な非線形性を有す
るプロセスの一種である反応器のプロセスフローチャー
トであり、ヘプタン（Ｃ₇Ｈ₁₆）を原料としトルエン
（Ｃ₇Ｈ₈）を製品として取り出す場合の反応例を示して
いる。FIG. 1 is a process flow chart of a reactor, which is a kind of process having a significant non-linearity in the components of the control system, and uses heptane (C ₇ H ₁₆ ) as a raw material and toluene (C ₇ H ₈ ) as a product. The example of reaction when taking out is shown.

第２図は、第１図に示す反応器のプロセス・パラメー
タの同定を行なう場合のブロック図である。同図におい
て、１はプロセス（反応器）、２はプロセス１の入出力
データをサンプリングするデータ・サンプリング部、３
は非線形プロセスモデルを求める非線形プロセスモデル
部、４は離散化した差分方程モデルを求める差分方程式
モデル部、５は求めたいパラメータに関する線形方程式
に変換する線形方程式モデル部,6は推定したいパラメー
タの大きさからスケーリング・ファクタを決めて推定式
を求める推定式決定部である。FIG. 2 is a block diagram for identifying process parameters of the reactor shown in FIG. In the figure, 1 is a process (reactor), 2 is a data sampling unit for sampling input / output data of the process 1, and 3
Is a non-linear process model part for obtaining a non-linear process model, 4 is a difference equation model part for obtaining a discretized difference method model, 5 is a linear equation model part for converting into a linear equation for a desired parameter, and 6 is a size of a parameter to be estimated. This is an estimation formula determination unit that determines a scaling factor from Eq.

また、10は逐次型最尤推定演算部であり、パラメータ
係数行列計算部11、誤差修正ゲイン計算部12、推定値計
算部13、推定誤差計算部14、推定誤差収束判定部15、推
定パラメータ修正部16、誤差共分散修正部17とからなっ
ている。そして、これら各部において、現在（ｋ時点）
の観測値から１サンプル後（ｋ＋１時点）の反応器の出
力を推定し、ｋ＋１時点における観測値と推定値の誤差
を求め（第４図参照）、この誤差が収束しているときに
は同定を終了し、収束していないときには、誤差にもと
づいて推定プロセス・パラメータを修正し、かつ、誤差
共分散を修正して同定を継続させる。Further, 10 is a sequential maximum likelihood estimation calculation unit, which includes a parameter coefficient matrix calculation unit 11, an error correction gain calculation unit 12, an estimated value calculation unit 13, an estimated error calculation unit 14, an estimated error convergence determination unit 15, and an estimated parameter correction. It comprises a unit 16 and an error covariance correction unit 17. And in each of these parts, now (at time k)
Estimate the output of the reactor one sample after (k + 1 time point) from the observed value, and obtain the error between the observed value and the estimated value at k + 1 time point (see Fig. 4). If this error converges, the identification is completed. If not, the estimation process parameter is corrected based on the error, and the error covariance is corrected to continue the identification.

次に、本実施例における反応器のプロセス・パラメー
タの同定方法を、第３図に示すフローチャートを参照し
て説明する。Next, the method for identifying the process parameters of the reactor in this example will be described with reference to the flowchart shown in FIG.

非線形プロセス・モデル部３において、非線形プロ
セス・モデルを求める（301）。The non-linear process model unit 3 obtains a non-linear process model (301).

○微分方程式モデル 第１図に示す反応器における熱収支と物質収支は次の
通りである。○ Differential equation model The heat balance and mass balance in the reactor shown in Fig. 1 are as follows.

イ．熱収支 dx₁/dt＝q/v・｛Ｔ_i−ｘ₁（ｔ）｝ −〔ΔH/ρｃ_p〕・ｋ（ｘ₁）ｘ₄ ＋〔ah/ρｃ_pｖ〕・（Ｔ_h−ｘ₁） ……［１］ ロ．物資収支 dx₂/dt＝−〔q/v〕・ｘ₂（ｔ） ＋ｋ（ｘ₁）・ｘ₄（ｔ） ……［２］ dx₃/dt＝−〔q/v〕・ｘ₃（ｔ） ＋ｋ（ｘ₁）・ｘ₄（ｔ） ……［３］ dx₄/dt＝−〔q/v〕・（ｘ₄）（ｔ） −4k（ｘ₁）・ｘ₄（ｔ） ＋〔q/v〕・ｄ₁（ｔ） ……［４］ ここで、 Ｔ_i：原料の入口温度 Ｔ_h：熱媒の入口温度 T:出口温度 a:伝熱面積 h:熱伝達係数 ｘ₁：反応温度 ｘ₂：トルエン出口温度 ｘ₃：水素出口濃度 ｘ₄：ヘプタン出口濃度 ρ：内部流体平均密度 Ｃ_p：内部流体平均比熱 ΔH:反応熱 ｋ（ｘ）：反応速度定数 Ｋ（ｘ₁）＝ｋ₀ e^-Ea/RT ｋ₀：頻度因子 ｄ₁：原料入口濃度 である。I. Heat balance _{dx 1 / dt = q / v} · {T i -x 1 (t)} - [[Delta] H / rho] c _{p] · k (x 1) x 4} + [ah / ρc _p v] · (T _h -x ₁ ) ...... [1] b. Material balance dx ₂ / dt = − [q / v] ・ x ₂ (t) + k (x ₁ ) ・ x ₄ (t) …… [2] dx ₃ / dt = − [q / v] ・ x ₃ ( _{t) + k (x 1)} · x 4 (t) ...... [3] dx 4 / dt = - [q / v] _{· (x 4) (t)} -4k (x 1) · x 4 (t) + [Q / v] · d ₁ (t) ··· [4] where T _i : raw material inlet temperature T _h : heat medium inlet temperature T: outlet temperature a: heat transfer area h: heat transfer coefficient x ₁ : Reaction temperature x ₂ : Toluene outlet temperature x ₃ : Hydrogen outlet concentration x ₄ : Heptane outlet concentration ρ: Internal fluid average density C _p : Internal fluid average specific heat ΔH: Reaction heat k (x): Reaction rate constant K (x ₁ ) = K ₀ e ^{−Ea / RT} k ₀ : frequency factor d ₁ : raw material inlet concentration.

次いで、差分方程式モデル部４において、離散化し
て差分方程式モデルを求める（302）。Next, the difference equation model unit 4 discretizes the difference equation model (302).

ここで、Δｔは時間刻みで、θはパラメータを示す。 Here, Δt is a time step, and θ is a parameter.

上記［１］，［２］，［３］，［４］式をそれぞれ
［５］式に代入する。The above expressions [1], [2], [3], and [4] are substituted into the expression [5].

ｘ₁（ｋ＋１）＝ｘ₁（ｋ）＋Δｔ〔q/v・｛Ｔ_i−ｘ
₁（ｋ）｝ −ΔH/ρＣ_p・｛ｋ（ｘ₁（ｋ））ｘ₄（ｋ）｝ ＋ah/ρｃ_pｖ・｛Ｔ_h−ｘ₁（ｋ）｝〕 ……［６］ ｘ₂（ｋ＋１）＝ｘ₂（ｋ）＋Δｔ〔−q/v・ｘ₂（ｋ） ＋ｋ｛ｘ₁（ｋ）｝ｘ₄（ｋ）〕 ……［７］ ｘ₃（ｋ＋１）＝ｘ₃（ｋ）＋Δｔ〔−q/v・ｘ₃（ｋ） ＋ｋ｛ｘ₁（ｋ）｝ｘ₄（ｋ）〕 ……［８］ ｘ₄（ｋ＋１）＝ｘ₄（ｋ）＋Δｔ〔−q/v・ｘ₄（ｋ） −4k｛ｘ₁（ｋ）｝ｘ₄（ｋ）＋q/v・ｄ₁（ｋ）〕 ……［９］ 次に、線形方程式モデル部５において、差分方程式
モデルを、求めたいパラメータに関する線形方程式に変
換する（303）。x ₁ (k + 1) = x ₁ (k) + Δt [q / v · {T _i −x
_{1 (k)} -ΔH / ρC} p · {k (x 1 (k)) x 4 (k)} + ah / ρc p v · {T h -x 1 (k)} ] ...... [6] x _{2 (k + 1) = x 2} (k) + Δt _{[-q / v · x 2 (k} ) + k {x 1 (k)} x 4 (k) ] _{...... [7] x 3 (k} + 1) = x 3 (k ) + Δt [−q / v · x ₃ (k) + k {x ₁ (k)} x ₄ (k)] …… [8] x ₄ (k + 1) = x ₄ (k) + Δt [−q / v · _{x 4 (k) -4k {x} 1 (k)} x 4 (k) + q / v · d 1 (k) ] ...... [9] then, in the linear equation model section 5, the difference equation model, determined Convert to a linear equation for the desired parameter (303).

上記［６］，［７］，［８］，［９］式から ｘ₁（ｋ＋１）＝ｘ₁（ｋ）＋Δｔ〔｛−ΔH/ρｃ_p・ exp（−Ｅ_a/RX₁（ｋ））ｘ₄（ｋ）｝・ｋ₀ ＋a/ρｃ_pｖ・｛Ｔ_h−ｘ₁（ｋ）｝・ｈ ＋q/v・｛Ｔ_i−ｘ₁（ｋ）｝〕 ……［10］ ｘ₂（ｋ＋１）＝ｘ₂（ｋ）＋Δｔ｛−q/v・ｘ₂（ｋ） ＋exp（−Ｅ_a/RX₁（ｋ））ｘ₄（ｋ）・ｋ₀｝ ……［11］ ｘ₃（ｋ＋１）＝ｘ₃（ｋ）＋Δｔ｛−q/v・ｘ₃（ｋ） ＋exp（−Ｅ_a/RX₁（ｋ））ｘ₄（ｋ）・ｋ₀｝ ……［12］ ｘ₄（ｋ＋１）＝ｘ₄（ｋ）＋Δｔ｛−q/v・ｘ₄（ｋ） −4 exp（−Ｅ_a/RX₁（ｋ））ｘ₄（ｋ）・ｋ₀ ＋q/v・ｄ₁（ｋ）｝ ……［13］ 次いで、推定式決定部６において、推定したいパラ
メータの大きさからスケーリング・ファクタを決めて、
推定式を求める（304）。 The [6], [7], [8], x 1 from [9] Formula _{(k + 1) = x 1} (k) + Δt _{[{- ΔH / ρc p · exp} (-E a / RX 1 (k)) _{x 4 (k)} · k} 0 + a / ρc p v · {T h -x 1 (k)} · h + q / v · {T i -x 1 (k)} ] ...... [10] x ₂ ( _{k + 1) = x 2 (} k) + Δt {-q / v · x 2 (k) + exp (-E a / RX 1 (k)) x 4 (k) · k 0} ...... [11] x 3 (k + 1 _{) = x 3 (k) +} Δt {-q / v · x 3 (k) + exp (-E a / RX 1 (k)) x 4 (k) · k 0} ...... [12] x 4 (k + 1) = X ₄ (k) + Δt {−q / v · x ₄ (k) −4 exp (−E _a / RX ₁ (k)) x ₄ (k) · k ₀ + q / v · d ₁ (k)} ... [13] Next, in the estimation formula determination unit 6, the scaling factor is determined from the size of the parameter to be estimated,
Obtain an estimation formula (304).

ここで、スケーリング・ファクタとは、各パラメータ
に対する各変数からの寄与が同程度になるように設定し
た係数である。Here, the scaling factor is a coefficient set such that the contributions from the variables to the parameters are similar.

一例として、反応温度の推定値に関する線形方程式
［10］について説明する。As an example, the linear equation [10] regarding the estimated value of the reaction temperature will be described.

ｘ₁（ｋ＋１）＝ｘ₁（ｋ）＋Δｔ〔｛−ΔH/ρｃ_p・ exp（−Ｅ_a/Rx₁（ｋ））ｘ₄（Ｋ）｝・ｋ₀ ＋a/ρｃ_pｖ・｛Ｔ_h−ｘ₁（ｋ）｝・ｈ ＋q/v・｛Ｔ_i−ｘ₁（ｋ）｝〕 ……［10］ ［10］式においてｋ₀＝5.01×10⁸ ｈ＝6.069×10⁵（J/g・mol・ｋ） であり、ｋ₀の値はｈの約10³倍である。 _{x 1 (k + 1) =} x 1 (k) + Δt _{[{- ΔH / ρc p · exp} (-E a / Rx 1 (k)) x 4 (K)} · k 0 + a / ρc p v · {T h −x ₁ (k)} · h + q / v · {T _i −x ₁ (k)}] ... [10] In the formula [10], k ₀ = 5.01 × 10 ⁸ h = 6.069 × 10 ⁵ (J / g · mol · k), and the value of k ₀ is about 10 ³ times h.

ここで、Ｆ（1,1）＝Δｔ｛−ΔH/ρｃ_p・ exp（−Ｅ_a/Rx₁（ｋ））・ｘ₄（ｋ） ……［14］ Ｆ（1,2）＝Δｔ〔a/ρｃ_pｖ・｛Ｔ_h −ｘ₁（ｋ）｝〕 ……［15］ とおく。Here, F (1,1) = Δt { -ΔH / ρc p · exp (-E a / Rx 1 (k)) · x 4 (k) ...... [14] F (1,2) = Δt [ _{a / ρc p v · {T} h -x 1 (k)} ] is denoted by ...... [15].

Δｘ₁＝ｘ₁（ｋ＋１）−ｘ₁（ｋ）−Δｔ・q/v・ ｛Ｔ_i−ｘ₁（ｋ）｝＝Ｆ（1,1）×ｋ₀−Ｆ（1,2） ×ｈ ……［16］ ここで、Δｘ₁からｋ₀,hを推定する場合、Ｆ（1,1）
が非常に小さな値になる。また、ｋ₀を少し修正すると
大きく影響が出ることになる。Δx ₁ = x ₁ (k + 1) −x ₁ (k) −Δt · q / v · {T _i −x ₁ (k)} = F (1,1) × k ₀ −F (1,2) × h ... [16] Here, when estimating k ₀ , h from Δx ₁ , F (1,1)
Becomes a very small value. Also, if k ₀ is modified a little, it will have a great effect.

そこで、推定したいパラメータの大きさを等しくする
ために［16］式を次のように置き換える必要がある。Therefore, it is necessary to replace equation [16] as follows in order to equalize the size of the parameters to be estimated.

Δｘ₁＝Ｆ（1,1）×10⁸×ｋ₀′＋Ｆ（1,2） ×10⁵×ｈ′ ここで、第１項の10⁸、第２項の10⁵がスケーリング・
ファクタである。Δx ₁ = F (1,1) × 10 ⁸ × k ₀ ′ + F (1,2) × 10 ⁵ × h ′ where 10 ⁸ of the first term and 10 ^{5 of} the second term are scaling
Is a factor.

Δｘ₁＝Ｆ（1,1）′×ｋ₀′＋Ｆ（1,2）′×ｈ′ ……
［17］ ［15］式が求める推定式である。Δx ₁ = F (1,1) ′ × k ₀ ′ + F (1,2) ′ × h ′ ...
[17] [15] is the estimation formula obtained.

ｋ₀′,h′の初期値を適宜定めて代入することによ
り、パラメータｋ₀,hを推定する。The parameters k ₀ and h are estimated by appropriately setting and substituting initial values of k ₀ ′ and h ′.

続いて、逐次型最尤推定演算部10において、次の推定
演算処理を行なう。すなわち、 パラメータ係数行列部11で、パラメータ係数行列を
計算し、プロセス・パラメータを推定する（305）。Then, the successive maximum likelihood estimation calculation unit 10 performs the following estimation calculation processing. That is, the parameter coefficient matrix unit 11 calculates a parameter coefficient matrix and estimates process parameters (305).

前述した反応温度におけるパラメータ推定例を示す。 An example of parameter estimation at the above-mentioned reaction temperature will be shown.

［10］，［14］，［15］式より ｘ₁（ｋ＋１）＝ｘ₁（ｋ）＋Ｆ（1,1）′×ｋ₀′＋Ｆ
（1,2）′ ×ｈ′＋Ｂ（ｋ） マトリックス表現すると ｘ₁（ｋ＋１）＝ｘ₁（Ｋ）＋Ｆ（ｋ）・Ａ（ｋ）＋Ｂ
（ｋ） 誤差修正ゲイン部12で、誤差修正ゲインを計算する
（306）。From the expressions [10], [14], and [15], x ₁ (k + 1) = x ₁ (k) + F (1,1) ′ × k ₀ ′ + F
(1,2) ′ × h ′ + B (k) When expressed in a matrix, x ₁ (k + 1) = x ₁ (K) + F (k) · A (k) + B
(K) The error correction gain unit 12 calculates the error correction gain (306).

［18］式より 誤差共分散行列Ｐ（ｋ）の初期値 誤差修正ゲイン Ｋ（ｋ＋１）＝Ｐ（ｋ）・Ｆ（ｋ）・〔１＋Ｆ
^t（ｋ） ・Ｐ（ｋ）・Ｆ（ｋ）〕^-1 ……［19］ において、ｋの初期値を０とし、計算を開始する。From the formula [18] Initial value of error covariance matrix P (k) Error correction gain K (k + 1) = P (k) · F (k) · [1 + F
^t (k) · P (k) · F (k)] ⁻¹ ... [19], the initial value of k is set to 0, and calculation is started.

次に、推定値計算部13で、推定値を求める（30
7）。Next, the estimated value calculation unit 13 obtains an estimated value (30
7).

推定誤差計算部14で、推定値と観測値より推定誤差
を計算する（308）。 The estimation error calculator 14 calculates an estimation error from the estimated value and the observed value (308).

Ｅ（ｋ）＝ｘ₁（ｋ＋１）−ｘ₁（ｋ＋１）……［21］ 推定誤差収束判定部15で、推定誤差が収束したかど
うかを判定する（309）。E (k) = x ₁ (k + 1) −x ₁ (k + 1) ... [21] The estimation error convergence determination unit 15 determines whether the estimation error has converged (309).

が閾値εより小さいかどうかを判断する。（ｋ時点にお
いて、ｎサンプル前（ｋ−ｎ時点）からｋ時点までの収
束状況を判断する）。 Is smaller than the threshold ε. (At the time point k, the convergence status from n samples before (the time point k-n) to the time point k is determined).

この結果、Ｊ≦εのときはパラメータを決定して（31
0）、同定を終了し、Ｊ≧εのときは、次のようにして
同定を続ける。 As a result, when J ≦ ε, the parameters are determined (31
0), identification is completed, and when J ≧ ε, identification is continued as follows.

推定パラメータ修正部16で、推定パラメータを修正
する（311）。The estimation parameter correction unit 16 corrects the estimation parameter (311).

［19］，［21］式より Ａ（ｋ＋１）＝Ａ_(k)＋Ｅ（ｋ）・Ｋ（ｋ＋１） ……
［22］ 誤差共分散修正部17で、誤差共分散を修正する（31
2）。From equations [19] and [21], A (k + 1) = A _(k) + E (k) · K (k + 1)
[22] The error covariance correction unit 17 corrects the error covariance (31
2).

Ｐ（ｋ＋１）＝Ｐ（ｋ）−Ｐ（ｋ）・Ｆ（ｋ）・〔１
＋Ｆ^t（ｋ） ・Ｐ（ｋ）・Ｆ（ｋ）〕^-1・Ｆ^t（ｋ）・Ｐ（ｋ）
……［23］ データサンプリングの時点を１ステップずらして
（313）、１サンプル後のデータサンプリングを行な
い、上述した［５］以下の計算，修正を行なう。P (k + 1) = P (k) −P (k) · F (k) · [1
+ F ^t (k) · P (k) · F (k)] ⁻¹ · F ^t (k) · P (k)
... [23] The data sampling time is shifted by one step (313), data sampling is performed one sample later, and the calculation and correction of [5] and above are performed.

このように本発明は、サンプル周期ごとにパラメータ
を修正しながら計算している。例えば、第４図に示すよ
うに１分ごと計算する場合、その計算ごとにパラメータ
を修正しながら所定の回数計算するので、何分か後に反
応状態が変化したとき（第４図参照）には、その近辺で
パラメータが変化することから異常を診断したり、予測
制御のモデルパラメータの修正などの処置を早期に行な
うことができる。As described above, the present invention calculates while correcting the parameter for each sample period. For example, as shown in FIG. 4, when the calculation is performed every 1 minute, the calculation is performed a predetermined number of times while correcting the parameters for each calculation, so when the reaction state changes after a few minutes (see FIG. 4). Since the parameters change in the vicinity thereof, it is possible to diagnose an abnormality or take an action such as correction of the model parameter of the predictive control at an early stage.

［発明の効果］ 以上のように本発明によれば、反応器のプロセス・パ
ラメータを正確に同定することにより、パラメータの変
化から早期にプロセスの変動を検出できるとともに、毎
回更新されるパラメータにより、反応器の変化に適応し
た制御を行なうことができる。[Effects of the Invention] As described above, according to the present invention, by accurately identifying the process parameters of the reactor, it is possible to detect the process variation early from the change of the parameters, and the parameters updated every time The control adapted to the change of the reactor can be performed.

【図面の簡単な説明】[Brief description of drawings]

第１図は反応器のプロセスフローチャートを示し、第２
図は第１図の反応器をプロセスとしてプロセスの同定を
行なう場合のブロック図、第３図は本発明実施例方法を
説明するためのフローチャート、第４図は逐次型最尤推
定法を説明するための曲線図を示す。 1:プロセス、2:データ・サンプリング部 3:非線形プロセス・モデル部 10:逐次型最尤推定演算部FIG. 1 shows a process flow chart of the reactor,
FIG. 1 is a block diagram for identifying a process by using the reactor of FIG. 1 as a process, FIG. 3 is a flow chart for explaining the method of the present invention, and FIG. 4 is a sequential maximum likelihood estimation method. A curve diagram for is shown. 1: Process, 2: Data sampling unit 3: Non-linear process model unit 10: Recursive maximum likelihood estimation operation unit

Claims (1)

【特許請求の範囲】[Claims] 【請求項１】イ．反応器の非線形プロセスモデルを、プ
ロセス・パラメータに関する線形差分方程式に変換し、 ロ．反応器におけるプロセス・パラメータである頻度因
子，熱伝達係数等についてのスケーリング・ファクタを
求めて推定式を設定し、 ハ．反応器からデータをサンプリングし、このサンプリ
ング周期ごとに、上記推定式により上記頻度因子，熱伝
達係数等を推定し、 ニ．誤差修正ゲインを計算し、 ホ．反応器出力の推定値を計算し、かつ、この推定値と
観測値より推定誤差を計算し、 ヘ．上記推定誤差が収束しているか否かを判断し、 ト．収束していない場合には、上記誤差修正ゲインおよ
び推定誤差にもとづいて推定した頻度因子，熱伝達係数
等を修正するとともに、誤差共分散を修正し、上記ハ．
〜ホ．を繰返して同定を継続し、 チ．収束している場合には同定を終了する、 ことを特徴とした反応器におけるプロセス・パラメータ
の同定方法。1. A. Convert the nonlinear process model of the reactor into a linear difference equation for the process parameters, b. Set the estimation formulas by obtaining the scaling factors for the process parameters in the reactor, such as frequency factors and heat transfer coefficients, and c. Data is sampled from the reactor, and the frequency factor, heat transfer coefficient, etc. are estimated by the above estimation formula for each sampling period, and d. Calculate the error correction gain, and e. Calculate the estimated value of the reactor output, and the estimated error from this estimated value and the observed value. F. It is judged whether the above estimation error has converged, and If it has not converged, the frequency factor, the heat transfer coefficient, and the like estimated based on the error correction gain and the estimation error are corrected, and the error covariance is corrected.
~ E. Repeat the above to continue the identification, and A method for identifying process parameters in a reactor, characterized in that the identification is terminated when the values converge.