JP2007260504A - Method of constructing system model - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide a method of constructing a system model which allows utilization of closed-loop data and is stable numerically. <P>SOLUTION: Assuming that the target system is linear, the method identifies the target system by determining the parameters for a state-space model created imaginarily from input and output data of the target system. In order to model the target system, the parameters for the state-space model are determined by carrying out Eigensystem Realization Algorithm (ERA) processing of the matrix that is formed by extending the polynomial for the time of the state-space model created imaginarily by a subspace identification method into the future region and subjecting the resultant factor matrix to singular value decomposition. <P>COPYRIGHT: (C)2008,JPO&INPIT

Description

  The present invention relates to a method for modeling a system having an input and an output, and a method for estimating an output for an input using the model, and more particularly to a method for modeling a system in a chemical plant.

  Conventionally, in order to adjust the production amount to the optimum amount in a chemical plant, in order to grasp the movement and characteristics of the reaction tower, the operable valve is moved to change the behavior of the reaction tower, that is, temperature change and product Conduct plant tests to understand changes in composition. Normally, a valve is automatically controlled by a regulator. However, when performing a plant test, the control by the regulator is disconnected and changed from an automatically controlled mode to a manual mode operated by a human. Then, operate the valve directly to collect data.

  However, such a manual mode is likely to cause operational problems that may cause the plant to become unstable or the composition of the product to deteriorate, and testing is difficult even if the staff actually operates it. There was sometimes. For this reason, it was necessary to be able to perform the preparation in the control mode that is a closed loop for adjusting the input according to the output.

  In order to avoid these problems, a system based on input / output data is used to obtain appropriate input values through simulation without having to perform tests to separate the actual control implemented in the plant. The model of the transfer function is estimated, the parameters of the function are estimated, and a system model for system identification is being constructed. By investigating the response of the output value to the input value using this modeled transfer function, it is not necessary to perform a plant test in which the control implemented in the plant is separated. Further, it is also possible to analyze the contents of the system that is in the black box state as a whole from the parameters of the obtained transfer function.

  However, the system to perform such identification is not only a system that can obtain output from simple input, but it also adjusts the valve that is the operation end at any time according to the production amount and components that are output in the chemical plant. Some have a recycling process in which the output is fed back to the input, such as when the output is converged to a constant by operating, or when a part of the product is returned to the raw material. For this reason, it cannot be simply modeled, and various methods are being studied.

  As a classical method, a prediction error method (PEM) that determines parameters of a linear polynomial using a maximum likelihood method that reduces the variance of a prediction error that occurs when a constructed linear polynomial model is applied to observation data. There is. Although this method can handle closed-loop data and has high accuracy, it requires convergence calculation, so it is not easy for numerical calculations, and it is complicated to handle multiple variables with multiple inputs and outputs. turn into.

  As another classical method, there is an autoregressive method (ARX) in which parameters of the constructed linear polynomial model are calculated using a least square method. This method can use closed-loop data and does not require convergence calculation, so it is easy to calculate, but it becomes complicated to handle multi-variables with multiple inputs and outputs, and the prediction accuracy is lower than that of the PEM. End up.

  In order to solve these complications, by calculating the subspace of the matrix in which the data is matrixed, it is possible to handle multiple variables with multiple inputs and outputs without performing convergence calculation, and necessary calculations A reduced amount of subspace method (SIM) has been developed. However, when this method is used for closed-loop data, it has the disadvantage that it calculates the input / output of the system including the closed-loop and cannot model the system alone.

  Apart from these, a method called Observer / Kalman filter Identification (OKID) has been developed that can construct a state space model using state variables and matrices and can handle multiple variables with multiple inputs and outputs. . This is because closed loop data can be used only in cases where the order is high, and multivariables can be easily handled. However, the accuracy of this state space model is simple, and it is often necessary to obtain accuracy comparable to PEM. Because it is necessary to use these parameters, it is difficult to use.

  In order to solve these problems, Subspace identification method (SSARX) described in Non-Patent Document 1 has been developed as a subspace method expanded so that closed loop data can be used. Specifically, this method is based on the following steps.

（ｃ）前記（２）式をカルマンプレディクター形式に変形する。
（ｄ）変形した式中のカルマンプレディクターマルコフパラメータを最小二乗法により算出する。
（ｅ）部分空間法の準備として状態空間モデルを未来部分へ拡張する。
（ｆ）前記状態空間モデルの係数を、サンプリングデータから最小二乗法を用いて求める。
（ｇ）前記（ｆ）の係数行列に特異値分解を施す。
（ｈ）特異値分解した結果から、時刻ｔにおける状態変数ｘを割り当てる。
（ｉ）時刻（ｔ＋１）における状態変数ｘを、前記（ｃ）乃至（ｇ）と同様に求める。
（ｊ）前記（ｂ）の（２）式の第二式から、パラメータＣを最小二乗法で求める。
（ｋ）前記（ｂ）の（２）式の第一式から、パラメータＡ，Ｂ，Ｋを最小二乗法で求める。 (A) Assume that the target system is a linear system.
(B) Using the output y, the input u, the state variable x, the innovation process e, the matrices A, B, C, and K, the input / output of the target system is temporarily constructed as a state space model of the following equation (2).

(C) The equation (2) is transformed into the Kalman predictor format.
(D) The Kalman predictor Markov parameter in the transformed equation is calculated by the method of least squares.
(E) Extend the state space model to the future part in preparation for the subspace method.
(F) The coefficient of the state space model is obtained from the sampling data using the least square method.
(G) Singular value decomposition is performed on the coefficient matrix of (f).
(H) The state variable x at time t is assigned from the result of singular value decomposition.
(I) The state variable x at time (t + 1) is obtained in the same manner as (c) to (g).
(J) The parameter C is determined by the least square method from the second expression of the expression (2) in (b).
(K) Parameters A, B, and K are obtained from the first equation of equation (2) in (b) by the least square method.

13th IFAC Symposium on System Identification P1625-1630

  However, in the SSARX method, after the coefficient matrix is subjected to singular value decomposition in (g) above, the error e generated using the parameter C obtained by the least square method using the inverse matrix in (j) is obtained, and The remaining parameters are obtained using the least square method, but the inverse square calculation is performed in the least square method. However, the inverse matrix calculation has a problem that the numerical value becomes unstable because of the possibility of division or divergence to zero.

  Accordingly, an object of the present invention is to provide a model construction method that can use closed-loop data and construct a more numerically stable model when constructing a model of a target system.

The present invention provides a system model construction method in which a target system is assumed to be a linear system, and parameters of a state space model modeled hypothetically from input / output data of the target system are obtained to perform system identification of the target system. ,
The state space model is obtained by performing Eigensystem Realization Algorithm (ERA) processing on a matrix obtained by singular value decomposition of a coefficient matrix obtained by extending a polynomial with respect to the time of the state space model assumed using the subspace method to the future part. The above-mentioned problem has been solved by a system model construction method that calculates the parameters of and models the target system.

  That is, in the conventional SSARX, the least square method using an inverse matrix after singular value decomposition was repeated, but in the present invention, by directly calculating a parameter by ERA processing that performs matrix calculation without using an inverse matrix, The model can be constructed stably numerically. Further, unlike a simple subspace method, since a coefficient matrix extended to the future part is used, closed loop data in which an output value is fed back to an input value can be used.

  By using the model construction method according to the present invention, it is possible to stably construct a model of a closed-loop system involving multiple variables without the need for convergence calculation. With the model thus made, it is possible to obtain an accuracy equal to or higher than that of the conventional SSARX model.

  The present invention will be described in detail below. The present invention assumes that the target system is a linear system and is a parameter of a state space model that is modeled hypothetically from input data and output data of the target system (hereinafter collectively referred to as “input / output data”). Is a system model construction method for performing system identification of the target system.

  The target system modeled by the present invention is a system having an output with respect to an input, and both the input and the output can be quantified. For example, mass, energy amount, volume, temperature, pressure, etc. The present invention is used for a system in which a plurality of input / output sampling data obtained by quantifying input / output can be obtained and the data can be assumed to be linear. In the system according to the present invention, even if the output is fed back to the input, and the output is closed-loop data that affects the future output, it can be modeled.

  In the present invention, it is assumed that the target system is a linear system from the input / output sampling data, and a state space model is created. This state space model is a polynomial composed of a product of a state variable x, an input u, an output y, and an error noise e with respect to a time t, and a feedback, by expressing the state variable x at a time (t + 1). Is in a possible format. If a part of the matrix as this parameter is collected, the number of terms can be reduced and transformed into a format that is easy to analyze. This format is a prediction format using a Kalman filter and is called a Kalman predictor format.

  Among the collected parameters, a Markov parameter including a coefficient matrix of the state variable x for the output y, a state variable x at the time t with respect to the state variable at the time (t + 1), and a coefficient matrix of the input u, and the input u and the output y Is calculated by the least square method. The procedure so far is the same procedure as the so-called ARX model and OKID model.

  Next, the above state space model is extended from time t to future time (t + f-1). With respect to a polynomial including a coefficient matrix relating to this future time, the coefficient matrix is obtained by least squares the input / output sampling data, and this coefficient matrix is subjected to singular value decomposition. The procedure so far is the same as SSARX.

The above singular value decomposition is a method of matrix calculation, and an o × q matrix X is converted to an o-order orthogonal matrix U, a q-order orthogonal matrix V, and an o × q-order matrix only in diagonal elements. This is a method of decomposing X = USV ^{T including} a matrix S having values.

Further, by performing the Eigensystem Realization Algorithm process (hereinafter abbreviated as “ERA process”) on the coefficient matrix subjected to the above singular value decomposition, the parameters of the state space model are calculated in advance. The target system can be modeled by obtaining a coefficient matrix of a polynomial in the state space model. Here, the ERA process is a method for obtaining a matrix corresponding to an eigen matrix subjected to singular value decomposition. As a specific example, for the coefficient matrices A, B, and C corresponding to the Markov parameters, a Hankel matrix: H _K whose matrix elements are CA ^K B is defined, and an eigen matrix of this Hankel matrix is obtained. Since the eigenmatrix does not change from the original matrix no matter how many times it is raised, A of the Markov parameters can be obtained first. Further, the remaining matrices B and C can be calculated as the remainder obtained by carving out the eigenmatrix A. As another method, there is a singular value decomposition method (SVD) in which A, B, and C are calculated by obtaining a singular matrix without obtaining an eigenmatrix that must be a square matrix. Hereinafter, description will be made on the assumption of ERA processing using the singular value decomposition method. Note that A, B, and C described here are Markov parameters of the equation (1) described later.

This ERA process is conventionally known. In order to perform singular value decomposition on a matrix, the input / output sampling data matrix is used to decompose the matrix into the above X = USV ^T format. A certain normalization coefficient W ₁ and W ₂ are modified to have at the beginning and end of matrix multiplication, and this normalization coefficient is used as a canonical correlation analysis (hereinafter abbreviated as CCA) weight. By performing the ERA process on the coefficient matrix having the parameters, the parameters can be obtained in correspondence with each other without using the least square method using the inverse matrix.

  By using the above CCA weights, if the target matrix is divided into a standardized matrix that is separate from the left and right in terms of calculation, whether the state that the portion carved on the left side can observe can be observed. Since the portion corresponding to the portion and cut to the right side corresponds to the portion of whether or not the state can be controlled, there is an advantage that the state of the process can be observed and controlled. On the contrary, if the separation is biased to only one of the left and right, only one of observation or control can be dominant.

  With these parameters, a polynomial in which the target system is modeled with the state space model is completed. By introducing the input data u into this polynomial, estimated output data y can be obtained. When operating the target system, using the model identified by the system model construction method according to the present invention and estimating the estimated output data of the target system for the input data in a plurality of input data as described above The input data from which the requested output data can be obtained can be examined. In other words, it becomes possible to introduce input data for obtaining suitable output data into the target system, and the target system can be operated efficiently.

  The estimated output data is compared with the measured output data obtained by actually executing the same input data on the target system, and the accuracy of the obtained state space model and parameters is confirmed. May be. The accuracy is adaptability, and indicates how well the estimated output data for the input data matches the actual output data. When it is determined that the accuracy is insufficient, it is preferable to further improve the accuracy by re-measurement or addition of the input / output sampling data and modeling by the above procedure.

  Examples of the target system to which such a system model construction method is applied include, for example, a multidimensional system that generates a plurality of types of products from a plurality of types of manufacturing materials in a chemical plant. Since such a system has a large number of variables to be analyzed, it is difficult to model with the conventional method, but by using the method according to the present invention, modeling with a stable numerical value can be performed.

  In such a chemical plant system, as input data to be modeled, it is preferable to select an element whose contribution to output is not negligible among values that can be quantified, such as raw material supply amount, component ratio, temperature, and pressure. . Here, the raw material is not limited to a simple manufacturing raw material, but may also include a material used indirectly such as water vapor used as a heat source for adjusting the reaction temperature and the like. This is because any element related to the target system can be considered. On the other hand, the output data includes the production amount of the product, the component ratio, the temperature, the pressure, and the like. Among these, it is preferable to select an element that cannot be ignored due to input. Here, the product includes not only the desired product but also by-products, waste liquids and wastes that are separated and discarded.

  As a method of using the system model construction method according to the present invention, the system model construction method according to the present invention calculates parameters A, B, C, and K of the state space model modeled from the input / output data. System identifying means for performing system identification, input means for inputting the input / output data of the target system, output means for displaying a state space model composed of the calculated parameters A, B, C, and K Can be used as an integrated system identification device for modeling by inputting the input / output data. In addition, the system identification device that performs modeling in this way further inputs estimated input data corresponding to the input data additionally input by using the modeled state space model by inputting the input data additionally. It is preferable that the value can be calculated because it is easier to use in the actual target system.

  Here, as the input means, in addition to an input device such as a keyboard, a numerical value is directly input using a sensor for measuring a physical property value of a manufacturing raw material in the target system or a sensor for measuring a physical property value of a product. You may do it. Examples of the output means include general output devices such as a display and a printer. The system identifying means is a computer that executes a program for realizing the model construction method according to the present invention. In addition, the input data that can be used as the physical property value includes at least one numerical value among the supply amount, the component ratio, the temperature, and the pressure of the manufacturing raw material used in the target system. Furthermore, the output data that can be used as the physical property value includes at least one numerical value among the production amount, component ratio, temperature, and pressure of the product obtained by the target system.

  More specifically, as a method of using the system model construction method according to the present invention, input means for inputting the input / output data of the target system and desired output data as target output data, and the present invention System identification means for performing system identification by calculating parameters A, B, C, and K of the state space model modeled from the input / output data by the model construction method according to the above, and the calculated parameter A , B, C, and K, the output calculation means for calculating the estimated output data for the new input data is compared with the desired output data and the estimated output data. Means for judging suitability of the input data, and the new input data for which the suitability satisfies the standard condition How to use the system operation device having execution means for executing the stem and the like.

  In this case, as the input means, in addition to the input devices and sensors enumerated above, the desired desired output data necessary for the operation of the target system, the operation results and scale of the target system so far And those that are automatically calculated from past information. In addition to the output device, the output means includes a regulator that directly adjusts a production raw material supply amount, a heat amount, a pressure, and the like in the target system. Thereby, it is possible to cause the actual target system to execute input data that satisfies the above-described reference conditions.

  The output calculation means is a computer that executes a program capable of matrix calculation using the state space model. Further, the suitability judging means compares the estimated output data with the desired output data, and whether or not the difference is permissible in the operation of the target system satisfies the reference condition. It is a computer that executes a program that judges whether or not. Here, the reference condition is input from the input means, and may be input directly by the administrator of the target system, or is calculated from the operation results of the target system so far. It may be a thing.

  Furthermore, the system identification means, the output calculation means, and the suitability determination means may be executed by a program that plays the respective roles on a single computer. The system model construction method according to the present invention is useful when the system control program that prepares the input means so as to play these roles and obtain an approximate value of desired output data as actual output data. It can be applied to the target system.

  Such a system operation device can be used in the above-described chemical plant system. As in the above, the input data and output data to be modeled include the supply amount of the manufacturing material and the posture object, the component ratio, the temperature. Among the values that can be quantified, such as pressure, it is preferable to select an element whose contribution to the output cannot be ignored at the input and an element that cannot ignore the contribution from the input at the output.

Hereinafter, a specific parameter identification procedure for constructing a model according to the present invention will be described.
First, a linear time-invariant state space model expressed in an innovation model form such as the following equation (1) is assumed. An innovation model is a model that includes disturbance data that cannot be estimated.

Where y (t) εR ^ny is the system output, u (t) εR ^nu is the system input, ny is the order of y, nu is the order of u, and x (t) εR ⁿ is the order. a state variable that is a vector of n, e (t) εR ^ny is an innovation process in which the mean value is zero and the covariance matrix is unknown, and each matrix is a matrix AεR ⁿ × ⁿ and a matrix Bε R ⁿ × ^nu , matrix CεR ^ny × ⁿ , and matrix KεR ⁿ × ^ny . The innovation process e (t) is disturbance data that enters the system from time to time, and is estimated from input / output data because the covariance matrix is unknown. Here, by focusing on the discrete-time system, the direct term D directly connected to the output y from the input u at the same time is set to zero, and the discussion proceeds with the model of the following equation (2).

The above formula (2) is rearranged into a predictor format like the following formula (3).
Note that the predictor format is a Kalman expressed using Kalman gain K, which will be described later, among observer formats in which the input u (t) and the output y (t) are collectively treated as a new input z (t). This is the Predictor format. In other words, in the following analysis, analysis is performed as an observer that considers the input and output of the system as one. This formula (3) is based on the following formulas (4) to (6).

  Here, the parameter K included in the above equation (4) is a coefficient matrix of the innovation process e in the above equation (2), which is a steady Kalman gain that converges and becomes a constant when sufficient time has elapsed. Since it is a called matrix, (A-) (indicating the left side of the above formula (4). For the notation in the text, the symbol having a bar on the matrix symbol is represented similarly below) is asymptotically stable. Can be considered. Therefore, the state vector x can be expressed by the state variable x, the input u, and the output y at an earlier time point from the above equation (3) by eliminating the unstable innovation process. It can be written as equation (7). Here, p is a positive number.

Here, since (A−) is asymptotically stable, the first term on the right side of the above equation (7) can be approximated to (A−) ^p = 0 for a sufficiently large positive number p. . Therefore, the first term on the right-hand side can be approximated to 0, and the state vector x (t) is the past input u (t), output y (t), and Kalman predictor Markov parameters (A−) and (B− ) From the following equation (8).

  Here, (x∧) (shows the left side of the above equation (8). In the notation in the text, the symbol having の 上 on the matrix symbol is also shown below) indicates the estimated value of the state variable x. . Next, by replacing the above equation (2) with the state variable estimated by the above equation (8), the output y can be approximated as the following equation (9).

  For input / output time series data, the above equation (9) can be summarized as the following equation (10). In addition, the symbol in following formula (10) is as following formula (11) thru | or (14), and N is the number of data in a formula, and is a positive number.

Further, the equation (10) is a approximate expression of y _t, and also with respect to u _t and e _t define the same matrix as the above (13). The Kalman predictor Markov parameter C (L−) _z in the above equation (10) has a least square solution represented by the following equation (15). This calculation procedure is obtained by multiplying the above equation (10) by applying an inverse matrix of Zp on both sides from the right. This investigation is equivalent to the least squares method. Also, the product of the e _t and Z _p are, e _t innovation, namely because of the nature of white, the expected value of the product is 0.
Note that R in the following equation (15) is a correlation matrix for the sampling data of the input / output data, and is represented by the following equations (16) and (17).

  The Kalman predictor Markov parameter element estimated by this equation (15) can be expressed as the following equation (18) from the above equation (5) and the above equation (11). The equation (15) corresponds to a so-called ARX model.

  Here, in the case of the conventional OKID algorithm, a system Markov parameter is derived from the Markov parameter of the above equation (15), a Hankel matrix is created from the Markov parameter, and the Hankel matrix is finally subjected to ERA processing. Calculate the matrix A, B, C, K of the system. However, in the system model construction method used in the present invention, as described above, the subspace method for extending the state space model to the future part by using the Kalman Predictor Markov parameter, which is an observer with input / output as z (t), is used. Using this, the Hankel matrix is obtained and ERA processing is performed. The procedure in appearance is more complicated than the conventional OKID, but the matrices A, B, C, and K can be obtained directly and the accuracy is high, so that the result is stable in numerical analysis.

With respect to the state space model of the above formula (3) estimated in this way, the time is extended from the basic time t to the future time (t + f-1), and the Kalman predictor Markov parameter is analyzed by the subspace method. . The procedure is as follows. First, the above equation (3) is expanded to the following equation (19). The contents of each symbol are as shown in the following formulas (20) to (25).

The Kalman predictor Markov parameter estimated in the above equation (15) is used as an element constituting the matrix (Φ−) of the above equation (22). Further, the state vector X (t) composed of the state function x is estimated from the controllable matrix (L−) _z and the past data Z _p as shown in the following equation (26).

  Using the above equation (26), the above equation (19) can be written as the following equation (27). This expression corresponds to the above expression (10) extended to the future time (t + f-1).

Next, in order to analyze the Kalman predictor Markov parameter, a Hankel matrix composed of the Kalman predictor Markov parameter is defined and estimated. The Hankel matrix is defined by the following formula (28), a matrix of _{fn y × p (n u +} n y).

  In the above equation (28), if k = 0, the above equation (28) becomes the following equation (29).

Next, in the above equation (27), the part contributing to the future input is removed from the part relating to the future output. That is, the term of Z _{f−1 on} the right side is moved to the left side. The part related to this output is defined as (Y˜) f as shown in the following equation (30).

On the other hand, (H−) ₀ can be estimated as the following equation (31) by the least square method using the above equation (17) which is a correlation function for the sampling data. However, (R∧) _fp is a correlation function defined by the following equation (32).

Further, when (R∧) _ff is defined as the following equation (33) as the correlation matrix for the sampling data in the time element f related to the future time, the above equation (17), which is the correlation matrix for the sampling data, and Using equation (27), (H−) ₀ can be transformed into equation (34) below. However, the normalization coefficients W ₁ and W ₂ are the following formula (35), and G is the following formula (36).

Here, the equation (36) is a canonical correlation matrix between the signals (Y˜) _f and Z _p . Canonical correlation analysis (hereinafter abbreviated as “CCA”) for performing singular value decomposition on G is performed. When singular value decomposition is performed assuming that the dimension of the system is n, the following equation (37) is obtained.

Here, U and V are normal orthogonal matrices, respectively, S is a diagonal matrix whose diagonal components are singular values, and the singular values are those whose values decrease in order along the diagonal components. It is. U _n and V _n are matrices each consisting of the first n columns of U and V, respectively. Further, S _n is a diagonal matrix composed of n singular values of S. Thereby, from the above equation (34), (H−∧) ₀ can be expressed as the following equation (38).

For the coefficient matrix subjected to singular value decomposition as described above, W ₁ and W ₂ that are normalization coefficients for normalization are included as inverse matrices corresponding to the denominator when normalizing in one dimension. That is, the following ERA processing is performed in a format including CCA weights.

First, when the matrix (H−) ₀ is constructed so as to be described using a Moore-Penrose pseudo inverse matrix H †, the following equation (39) is obtained.

Here, by using the relationship of ^{_{U n T U n = I n}} , and ^{_{V n T V n = I n}} , when to apply the relation of the equation (39) into the equation (38), the following equation (40) become that way.

Next, in order to select (Γ−) _f and (L−) _z used to describe the matrix (H−) ₀ in the above equation (39), from the idea of parallel realization, (Γ−) ^{_{f = W 1 -1 U n S}} n 1/2 and include a method of determining as ^{_{(L-) z = S n 1/2}} V n T W 2. Moreover, the relationship of following formula (41) is materialized.

Further, I _i is an i-th order unit matrix, and O _i is an i-th order zero matrix. Since n _u is the number of input variables and n _y is the number of output variables, the following matrices (42) and (43) are defined.

  When these definitions are used to perform CCA-weighted ERA processing, the process and result are expressed by the following equation (44).

  In addition, the said Formula (41) is applied to the calculation from the 3rd line of the said formula to the 4th line. From the above equation (44), among the Kalman predictor Markov parameters, (A−), (B−), and C are calculated as the following equations (45) to (47).

However, (A−) includes (H−) ₁ which is a Hankel matrix, and it is necessary to separately obtain (H−) ₁ . From the above equation (28), which is the Hankel matrix definition, (H−) ₁ is as follows.

The second row of the above equation (48) selects the second to p-th column blocks in the controllable matrix of the matrix (H−) _{0 and} starts from the first in the controllable matrix of (H−) 1. Since (A−) ^p can be approximated to 0, the last p block is obtained by assigning a zero matrix to the p−1th column block.

  Accordingly, the Kalman predictor Markov parameters of the state space model defined in the above equation (1) are derived from the defined equations (4) through (6) as in the following equations (49) through (52). be able to. The expression method uses the description method of MATLAB (registered trademark).

  Thus, the system model construction method according to the present invention can be realized by obtaining the Kalman predictor Markov parameters. By using the method according to the present invention, it is possible to reduce the number of times of use of the least square method compared to the conventional SSARX, and to identify a numerically more stable model and thereby estimate the output data.

Examples using the present invention will be described below with specific examples.
The example taken up by Verhagen, 1993 is used as a closed loop data application verification system. This system is a fifth-order system having a parameter of the following formula (53) modeled in the form of the above formula (2).

  From this parameter, the controller described by the state space model can express the parameter in the above equation (2) as the following equation (54). The controller feeds back the output to the input and adjusts the output as desired, and corresponds to K (q) in FIG. 3 described later.

  The above was determined and a unit circle was plotted. The result is shown in FIG. The curve in FIG. 1 is a unit circle, the horizontal axis represents the real number axis, and the vertical axis represents the imaginary number axis. Using the Monte Carlo method for this target process, 100 sets of closed-loop data consisting of 1200 points created as a reference signal r is a Gaussian white signal with variance 1 and an innovation process e is a Gaussian white signal with variance 1/9. Create 100 data sets each.

  This data set was introduced into a model (denoted as “SOKID” in FIG. 1) modeled by the system model construction method according to the present invention as an example, and its eigenvalue was calculated and plotted in FIG. In addition, about this model, as a method of determining the order of the past lag p and the future lag f, it determined by FOE (Final output error criterion) (Zhu, 2001). Note that FOE refers to a method of comparing a measured value with an estimated output value of an identified model and searching for a value with the smallest mean square value of the error. However, although the order n of the above formula (53) is 5, it is necessary to satisfy the condition that n cannot be larger than the smaller value of p or f. As a result of calculation for 100 data sets, the past lag p and the future lag f were set to 20.

  Similarly, as a comparative example, a data set similar to the above was introduced into a model modeled by conventional SSARX and OKID, and its eigenvalue was calculated and plotted in FIG. However, since OKID is not a subspace method, F = 1 was set. In addition, in a model modeled by N4SID, which is a subspace method function of the system identification toolbox, installed in the matrix calculation software Matlab (registered trademark), a similar data set is introduced, and its eigenvalues are calculated. Plotted to 1.

  As a result, it is shown that the variation of the eigenvalue of SOKID using the system model construction method according to the present invention is almost the same as SSARX, and the error is as small as SSARX. Further, it can be seen that the OKID has a large variation and a large error compared to these, and that the N4SID has a small variation but is biased and the identification of the closed loop is not successful. Therefore, it has been found that the modeling by the system model construction method according to the present invention enables simulation with an accuracy equal to or higher than that of the conventional method.

  The system shown in FIG. 2 was modeled using actual data. First, this plant will be described. In this plant, the amount of the production material gas F (FreshFeed) is adjusted by a valve 11, and this valve is prepared by a command from a temperature controller provided after the production material gas passes through the heat exchanger 12 and the heat exchanger 13. And is introduced into the reactor R1 so as to be heated to an appropriate temperature. The high-temperature gas reacted from the bottom of the reactor R1 is separated in the middle, and a part of the heat is exchanged with the gas before being introduced into the reactor R1 by the heat exchanger 13 and cooled, and then introduced into the reactor R2. The remainder is directly introduced from the reactor R1 into the reactor R2. Here, the valve 14 in the line directly introduced from the reactor R1 to the reactor R2 is controlled by the temperature controller TC2 to control the temperature of the mixed gas of the gas that has passed through the heat exchanger 13 and the gas that has not passed through the reactor. Valve 14 is prepared to achieve the desired temperature at R2. Further, the amount of heat of the product gas P (Product) discharged from the reactor R2 is transferred to the production material gas in the heat exchanger 12.

  The exhaust temperatures from the reactors R1 and R2 are measured by the temperature sensors TI1 and TI2, and the inputs to the respective sensors are the disturbance variables DV1 and DV2. The operation signal to the valve 14 is a disturbance variable DV3.

  Further, each temperature controller TC1 and TC2 measures the indicated value of the temperature, and calculates the signal value to be operated to the valves 11 and 14 so that the deviation between the set value given by the operator and the indicated value is reduced. , Output operation to each valve via a signal line. Specifically, a PID controller is used.

  Here, the output y (t) is defined as an instruction value of the temperature of the gas detected by the temperature controller TC1, and the input u (t) is defined as the supply amount of the gas prepared by the valve 11. The operation of this system can be expressed as shown in FIG. The heat exchanger 12 and the heat exchanger 13 are the system G (q), and a disturbance v (t) that combines the disturbance variables DV1 to DV3 is involved, and an output y (t) is made. Further, the output y (t) is fed back, and the controller K (q) controls the input u (t) so that the deviation between the set value r (t) and the output y (t) decreases.

  The correspondence between the input u (t) and the output y (t) at that time is as shown in FIG. The horizontal axis indicates the number of data points measured at 1 minute intervals, and the vertical axis indicates the temperature (° C.) for y (t) and the controller output (%) for u (t). On the other hand, FIG. 4 shows data obtained by performing a plant test with the controller turned off, that is, data obtained in an open loop. The result of analysis based on the obtained response from the plant test data obtained by raising and lowering the process input u (t) and taking the process response is shown as a step response model in the frequency domain. This is indicated by the bold line 5. Treat this as the true value of the process and compare it to the value modeled in the following procedure.

  Next, FIG. 6 shows normal operation using the controller K (q) in the process of FIG. 2, that is, operation data in a closed loop.

  On the other hand, the process was identified using the closed loop data shown in FIG. 6, and the step response model in the frequency domain was obtained. The output value estimated by the system model construction method (SOKID) according to the present invention is indicated by a solid line in FIG. 5 by using the FOE norm as in the first embodiment and setting the order n suitable for minimizing the FOE to 9. Note that the data sampling period is 2500 points per minute, and the same applies to the comparative example described later.

  Further, as a comparative example, as in the first embodiment, a model modeled with SSARX, OKID, and Matlab-N4SID (CVA weight) is used, and the FOE norm is minimized using the FOE norm as in the above SOKID embodiment. The output value with respect to the frequency estimated by selecting and modeling the order was obtained. The result is shown in FIG. In addition, a dash-dot line is SSARX, a dotted line is OKID, and a broken line is Matlab-N4SID.

  Furthermore, FIG. 7 shows the absolute value of the deviation of the estimated value calculated by each modeling for the true process model represented by the solid line in FIG. Similarly, a solid line is a system model construction method (SOKID) according to the present invention, a dash-dot line is SSARX, a dotted line is OKID, and a broken line is Matlab-N4SID. When the mean square error with respect to the frequency was calculated for this absolute value, SOKID was 0.22, SSARX was 0.34, OKID was 0.29, and N4SID was 0.31. That is, in modeling this system, an estimated value obtained by performing modeling using the system model construction method according to the present invention was the closest to the true process value.

The graph which shows the eigenvalue distribution of the true value and estimated value in an Example Schematic of the plant modeled in Example 2 Schematic of closed loop showing mechanism of feedback from output in embodiment 2 The graph which shows the sampling data of the input u in the open loop which turned off the control of the controller in Example 2, and the output y The graph which shows the true value and estimated value of the frequency response of the gain in Example 2 The graph which shows the sampling data containing the disturbance variable in the closed loop data which turned on control of the controller used in Example 2 The graph which shows the absolute value of deviation with the true value in FIG.

Explanation of symbols

r (t) Set value u (t) Input v (t) Disturbance y (t) Output DV1, DV2, DV3 Disturbance variable G (q) System K (q) Controller F Production material gas P Product gas R1, R2 Reactor TC1, TC2 Temperature controller TI1, TI2 Temperature sensor 11, 14 Valve 12, 13 Heat exchanger

Claims (7)

Assuming the target system is a linear system, a system model construction method for performing system identification of the target system by obtaining parameters of the state space model from input data and output data of the target system,
The parameters of the state space model are calculated by performing Eigensystem Realization Algorithm on a matrix obtained by singular value decomposition of a coefficient matrix obtained by extending a polynomial with respect to the time of the state space model assumed using the subspace method to the future part. Then, a system model construction method for modeling the target system. By the following steps (A) to (H),
The system identification method according to claim 1, wherein a polynomial parameter of a state space model modeled on the assumption that the target system is a linear system is obtained and system identification of the target system is performed.
(A) The input data and output data of the process plant to be identified are collected, and the process plant is assumed to be a linear system.
(B) y (t) εR ^ny is the output of the system, u (t) εR ^nu is the input of the system, ny is the order of y, nu is the order of u, and x (t) εR ⁿ is the order n state vector, e (t) εR ^ny is an innovation process with zero mean and unknown covariance matrix, matrix AεR ⁿ × ⁿ , matrix BεR ⁿ × ^nu , matrix Cε A state space model for the time t in the following equation (2), which is R ^ny × ⁿ and the matrix KεR ⁿ × ^ny , is constructed.

(C) The above equation (2) is transformed into the following equation (7) by the following equations (4) to (6).

(D) A Kalman Director Markov parameter for obtaining an approximate value of the state vector x (t) together with the input u (t) and the output y (t) is obtained by the method of least squares.
(E) Extend the state space model to the future part (t + f-1) using the Kalman predictor Markov parameter.
(F) Define a Hankel matrix of the following formula (28) having the Kalman predictor Markov parameter as an element, and estimate the Hankel matrix when k = 0 by a least square method using a correlation function for the sampling data To do.
(G) A canonical correlation analysis is performed in which a singular value decomposition is performed on a matrix including the estimated normalization coefficient. (Note that p represents a positive number equal to or less than t.)

(H) Eigensystem Realization Algorithm processing including a normalization matrix for normalization by canonical correlation analysis is applied to the coefficient matrix subjected to the singular value decomposition, and the above expressions (4) and (5) The parameter represented by the left side of C and the above C are calculated.
(I) The parameters A, B, and K are calculated from the parameters calculated in the above (H) and the above C by the above formulas (4) to (7).   Using the state space model obtained by the system model construction method according to claim 1 or 2 and the parameter, estimated output data for input data is obtained, and measured output data for the input data is measured in the target system. A system model construction method for comparing the estimated output data and the actually measured output data to confirm the accuracy of the state space model and the parameters.   The system output estimation method which calculates | requires the estimated output data with respect to the input data to the said object system using the said space state model and the said parameter which were calculated | required by the system model construction method in any one of Claim 1 thru | or 3. In the system that compares the actual process output data actually obtained in the target system with the desired output data desired, and calculates the input value to the actual process so that the error is reduced,
The input data given to the target system is calculated so that the estimated output data estimated by the system output estimation method according to claim 4 is close to the desired output data, and the input data is input to the target system. , System operation method. The target system is a system that produces one or more products from one or more manufacturing materials in a chemical plant;
The input data is a numerical value of at least one of the supply amount, the component ratio, the temperature, and the pressure of the manufacturing material to be introduced into the chemical plant, and the output data is the product output, the component ratio, the temperature, and the pressure. At least one number,
The system operation method according to claim 5, wherein a production plan of a necessary product in the target system is operated. Input means for inputting the input data and output data of the target system;
System identification means for performing system identification by calculating parameters A, B, C, K of the state space model modeled from the input data and output data by the model construction method according to claim 1 or 3; ,
A system modeling device having output means for outputting a state space model composed of the calculated parameters A, B, C, and K.

Images