CN101872181A - Prediction method for monitoring performance of power plant instruments - Google Patents

Abstract

Disclosed is a prediction method for monitoring performance of power plant instruments. The prediction method extracts a principal component of an instrument signal, obtains an optimized constant of a SVR model through a response surface methodology using data for optimization, and trains a model using training data. Therefore, compared to an existing Kernel regression method, accuracy for calculating a prediction value can be improved.

Description

The Forecasting Methodology that is used for monitoring performance of power plant instruments

The cross reference of related application

The application requires to submit on April 22nd, 2009, sequence number is the right of priority of the korean patent application of 10-2009-0035254, and its in full content be incorporated into this by reference.

Technical field

Embodiment relates to the Forecasting Methodology that is used for monitoring performance of power plant instruments, more specifically, relates to a kind of being used in the nuclear power plant running technology of the performance monitoring of generating plant safety monitoring device always.

Background technology

In general, all power generating equipments all are equipped with large number quipments and obtain various signals in real time from described large number quipments, to utilize the various signals that obtained in generating plant supervision and the protection system.Especially, the measurement passage of the nuclear power plant relevant with security system uses many device concept to guarantee the accuracy and the reliability of measuring-signal, and put down in writing by the operative technique guide, check and proofread and correct power plant equipment with the interval (about 18 months) of nuclear fuel cycle.In worldwide, nuclear power plant is developing a kind of method always, is used for that (Condition Based Monitoring, CBM) method prolongs the monitoring and the calibration cycle of inessential equipment calibration task by the monitoring based on condition.

Fig. 1 is the block diagram of legacy equipment performance rule monitoring system.This system is called as (Auto-Associative) model from association.With reference to figure 1, described legacy equipment performance rule monitoring system comprises forecast model, comparison module and decision logic.Described legacy equipment performance rule monitoring system can by measuring-signal is input to described forecast model, the described forecast model of output is input to described decision logic and continues to monitor described equipment about the predicted value of input measurement value, difference with measured value and described predicted value in comparison module, with the drift and the fault of monitoring equipment.

On the other hand, Argonne National Laboratory (Argonne National Laborary) develops a kind of polynary state estimation technology (MSET) and obtained to be used to develop the United States Patent (USP) of MSET.SmartSignal company and expert's micro-system utilize this United States Patent (USP) to make a kind of product that is used for commercial object.Expert's micro-system has been installed a kind of Palo Verde that utilizes the U.S., and Limerick 1/2, TMI, and V.C.Summer, the product that the MSET in Sequoyah1 and Salem 1 field makes, and measure the in-service monitoring of passage.In addition, because SmartSignal company can not use the relevant United States Patent (USP) of MSET, it has developed a kind of technology based on nuclear homing method monitoring equipment performance.

Linear regression method is widely used as a kind of method of computing equipment predicted value.This method selects to have with the signal of equipment to be predicted the signal of other relevant equipment of highly linear, and obtains regression coefficient, so that the error sum of squares of predicted value and measured value is a minimum.This method can be in order to equation 1 expression down.

＜equation 1 〉

∑E ²＝∑(Y-Y′) ²

In case utilize known dependent variable and independent variable to determine regression coefficient, described linear regression method can be predicted the independent variable about unknown dependent variable.

Yet, in existing linear regression method, if dependent variable has big linear relationship, multicollinearity may appear, so that occurring big error for the little noise in being included in dependent variable in the independent variable.

The nuclear homing method is a kind of distribution-free regression procedure, its storage selects measurement data as memory vector (memory vector), from Euclidean distance, obtain the weighted value of kernel about training data group sets of measurement signals, that the memory vector, set, and described weighted value is applied in the described memory vector, to obtain the predicted value of measuring equipment, and need not operation parameter, as regression coefficient, weighted value (it optimizes input/output relation as existing linear regression method) or neural network.As examine homing method as described in distribution-free regression procedure on the model of nonlinear state, have very strong advantage with input/output relation and signal noise.

A kind of existing Kernel regression method has the calculation procedure of following five steps.

The first step: with the matrix representation of training data with equation 2.

＜equation 2 〉

$X=\left[\begin{array}{cccc}{\mathrm{<figure-callout\; id="1"\; label="X"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">X</figure-callout>}}_{\mathrm{1,1}}& {\mathrm{<figure-callout\; id="1"\; label="X"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">\; <figure-callout\; id="2"\; label="X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">X</figure-callout>\; </figure-callout>}}_{\mathrm{1,2}}& \mathrm{\&Lambda;}& {\mathrm{<figure-callout\; id="1"\; label="X"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">X</figure-callout>}}_{1,\mathrm{<figure-callout\; id="2"\; label="m\; X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">m</figure-callout>}}& \\ X_{}{}_{\mathrm{2,1}}& {\mathrm{<figure-callout\; id="2"\; label="X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">\; <figure-callout\; id="2"\; label="X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">X</figure-callout>\; </figure-callout>}}_{\mathrm{2,2}}& \mathrm{\&Lambda;}& {\mathrm{<figure-callout\; id="2"\; label="X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">X</figure-callout>}}_{2,m}\\ M& M& O& M\\ X_{n_{\mathrm{TRAN}},1}& X_{n_{\mathrm{TRAN}},2}& \mathrm{\&Lambda;}& X_{n_{\mathrm{TRAN}},m}\end{array}\right]$

Wherein, X is the training data matrix that is stored in the memory vector, and n is the number of training data, and m is the number (a number of a instrument) of equipment.

Second step: the summation of Euclidean distance that obtains to be used for the training data of the first device signal group by following equation 3.

＜equation 3 〉

$d(x_1,q_1)=\sqrt{\mathrm{\&Sigma;}{(x_{1,j}-{\mathrm{<figure-callout\; id="1"\; label="q"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">q</figure-callout>}}_{1,j})}^2}$

$d(x_2,q_1)=\sqrt{\mathrm{\&Sigma;}{(x_{2,j}-{\mathrm{<figure-callout\; id="1"\; label="q"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">q</figure-callout>}}_{1,j})}^2}$

.

.

.

$d(x_{\mathrm{trn}},q_1)=\sqrt{\underset{j}{\mathrm{\&Sigma;}}{(x_{\mathrm{trn},j}-{\mathrm{<figure-callout\; id="1"\; label="q"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">q</figure-callout>}}_{1,j})}^2}$

Wherein x is a training data, and q is test data (perhaps, data query), and trn is the number of training data, and j is the number of equipment.

The 3rd step: by the following weighted value that comprises equation 4 acquisitions of kernel function about each training data group and given test data set.

＜equation 4 〉

${\mathrm{<figure-callout\; id="1"\; label="w"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">w</figure-callout>}}_1=\sqrt{K\left(d(x_1,q_1)\right)}$

${\mathrm{<figure-callout\; id="2"\; label="w"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">w</figure-callout>}}_2=\sqrt{K\left(d(x_2,q_1)\right)}$

.

.

.

$w_{\mathrm{trn}}=\sqrt{K\left(d(x_{\mathrm{trn}},q_1)\right)}$

In equation 4, the gaussian kernel that is used as weight function is defined as follows:

$K\left(d\right)=e^{-\left(\frac{{\mathrm{<figure-callout\; id="2"\; label="d"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">d</figure-callout>}}^2}{{\mathrm{\&sigma;}}^2}\right)}$

The 4th step: by each training data be multiply by weighted value, then with its result divided by weighted value and obtain the predicted value of test data, as following equation 5.

＜equation 5 〉

The 5th step: repeat second and went on foot for the 4th step, to obtain the predicted value of whole test data.

Existing Kernel regression method has very strong advantage aspect nonlinear model and the signal noise.

Yet, to compare with the linear regression analysis method, existing Kernel regression method for example has because the defective of the low precision that the dispersion degree increase of prediction of output value causes.Because AAKR method storage selective measurement data are as the memory vector, from Euclidean distance, obtain the weighted value of kernel about training data group sets of measurement signals, that the memory vector, set, and use this weighted value to described memory vector obtaining the predicted value of equipment, so dispersion degree increases.

Summary of the invention

Thereby embodiment is at a kind of Forecasting Methodology that is used for the shop equipment performance monitoring, and it has overcome basically because the limitation of correlation technique and one or more problem that defective causes.

Therefore, a feature of embodiment provides a kind of Forecasting Methodology of utilizing principal element analysis and support vector regression (SVR) to be used for the shop equipment performance monitoring, compares with known nuclear homing method, improves the computational accuracy of predicted value.That is to say that this method is used to improve traditional low precision of prediction.This method utilization uses the plant data normalization of the SVR model homing method of response surface optimization, major component to extract, parameter (the wide σ of nucleus band, loss function constant ε, penalties C) is optimized, by above-mentioned to the realization of SVR model and the reverse method for normalizing of prediction of output data, so that modeling factory system, monitoring equipment signal then.

In above-mentioned and other features and the advantage at least one can realize by a kind of Forecasting Methodology that is used for the shop equipment performance monitoring is provided, described method is extracted the major component of generating plant data, make up various instance models about the system that uses the SVR method, utilize the response surface analysis method to optimize three parameters of regression equation, and after utilizing described three described power plant systems of parameter model the monitoring equipment signal.Therefore, compare with the existing nuclear homing method of widespread use, the computational accuracy of predicted value can improve.

Description of drawings

By with reference to the accompanying drawings typical embodiment being described in detail, it is more obvious that above-mentioned and other feature and advantage will become for a person skilled in the art.

Fig. 1 is the block diagram of legacy equipment performance rule monitoring system;

Fig. 2 is the process flow diagram that illustration is used for the Forecasting Methodology of monitoring performance of power plant instruments according to embodiments of the present invention;

Fig. 3 illustration based on the Forecasting Methodology that is used for monitoring performance of power plant instruments of foundation embodiment of the present invention and the SVR model that produces;

Fig. 4 illustration the synoptic diagram of the ORL that produces by SVR;

Fig. 5 illustration when the number of model parameter is 3 the experimental point of center combination design;

Fig. 6 A and Fig. 6 B illustration be used for reactor core power data accuracy of detection, nuclear power plant according to an embodiment of the present invention;

Fig. 7 A and Fig. 7 B illustration be used for pressurizer horizontal data accuracy of detection, nuclear power plant (pressurizer level data) according to an embodiment of the present invention;

Fig. 8 A and Fig. 8 B illustration be used for the steam flow data of steam generator accuracy of detection, nuclear power plant according to an embodiment of the present invention;

Fig. 9 A and Fig. 9 B illustration be used for the horizontal data in a narrow margin of steam generator accuracy of detection, nuclear power plant according to an embodiment of the present invention;

Figure 10 A and Figure 10 B illustration be used for the pressure data of steam generator accuracy of detection, nuclear power plant according to an embodiment of the present invention;

Figure 11 A and Figure 11 B illustration be used for the wide cut horizontal data of steam generator accuracy of detection, nuclear power plant according to an embodiment of the present invention;

Figure 12 A and Figure 12 B illustration be used for the main feed flow data on flows of steam generator accuracy of detection, nuclear power plant according to an embodiment of the present invention;

Figure 13 A and Figure 13 B illustration be used for turbine output data accuracy of detection, nuclear power plant according to an embodiment of the present invention;

Figure 14 A and Figure 14 B illustration be used for first cycle charging flow accuracy of detection, nuclear power plant (loop charging flow) data according to an embodiment of the present invention;

Figure 15 A and Figure 15 B illustration be used for the data on flows that removes waste heat accuracy of detection, nuclear power plant according to an embodiment of the present invention; And

Figure 16 A and Figure 16 B illustration be used for the temperature data of reactor roof liquid coolant accuracy of detection, nuclear power plant according to an embodiment of the present invention.

Embodiment

Describe typical embodiment in detail below in conjunction with accompanying drawing, yet it can be presented as various form, and should not be interpreted as being restricted to listed embodiment herein.On the contrary, what provided makes the described embodiment that the disclosure is complete and complete, covers scope of the present invention to those skilled in the art fully.

Fig. 2 is the process flow diagram that illustration is used for the Forecasting Methodology of monitoring performance of power plant instruments according to embodiments of the present invention.With reference to figure 2, described Forecasting Methodology is included in the total data of display matrix among the operation S100, normalization total data in operation S200, in operation S300 is three parts with data component, as training, optimize and test, in operation S400, extract the major component of each normalization data group, in operation S500, utilize the response surface method to calculate support vector regression (SVR) Model Optimization constant to optimize the predicted value error of data, in operation S600, utilize described optimization constant to produce training pattern, utilize among the S700 normalization test data to obtain kernel function and predict the output of training pattern in operation, and normalization output that will be given in operation S800 oppositely is normalized to the predicted value of initial range with the acquisition variable as input.

According to an embodiment of the present invention, comprise that the Forecasting Methodology of monitoring performance of power plant instruments of said structure is as follows.

The first, the demonstration of matrix has shown the total data in the matrix shown in the following equation 6 among the operation S100.

＜equation 6 〉

$X=\left[\begin{array}{cccc}{\mathrm{<figure-callout\; id="1"\; label="X"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">X</figure-callout>}}_{\mathrm{1,1}}& {\mathrm{<figure-callout\; id="1"\; label="X"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">\; <figure-callout\; id="2"\; label="X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">X</figure-callout>\; </figure-callout>}}_{\mathrm{1,2}}& \mathrm{\&Lambda;}& {\mathrm{<figure-callout\; id="1"\; label="X"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">X</figure-callout>}}_{1,\mathrm{<figure-callout\; id="2"\; label="m\; X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">m</figure-callout>}}& \\ X_{}{}_{\mathrm{2,1}}& {\mathrm{<figure-callout\; id="2"\; label="X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">\; <figure-callout\; id="2"\; label="X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">X</figure-callout>\; </figure-callout>}}_{\mathrm{2,2}}& \mathrm{\&Lambda;}& {\mathrm{<figure-callout\; id="2"\; label="X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">X</figure-callout>}}_{2,m}\\ M& M& O& M\\ X_{3n,1}& {\mathrm{<figure-callout\; id="2"\; label="X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">X</figure-callout>}}_{3n,2}& \mathrm{\&Lambda;}& X_{3n,m}\end{array}\right]=\left[\begin{array}{cccc}{X}_1& {\mathrm{<figure-callout\; id="2"\; label="X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">X</figure-callout>}}_2& \mathrm{\&Lambda;}& X_m\end{array}\right]$

X _ts＝[X _3i+1，1?X _3i+1，2?Λ，X _3i+1，m]＝[X _ts1?X _ts2?Λ?X _tsm]

X _tr＝[X _3i+2，1?X _3i+2，2?Λ??X _3i+2，m]＝[X _tr1?X _tr2?Λ?X _trm]

X _op＝[X _3i+3，1?X _3i+3，2?Λ??X _3i+3，m]＝[X _op1?X _op2?Λ?X _opm]

Wherein X is whole data set.X _Tr, X _OpAnd X _TsBe respectively the data set that is used to train, the data set that is used to optimize and the data set that is used to test.3n is the number of total data, and m is the number of equipment.

Secondly, the normalization operation among the operation S200 utilizes following equation 7 normalization total datas.

＜equation 7 〉

$Z_i=\frac{X_i-\min \left(X_i\right)}{\max \left(X_i\right)-\min \left(X_i\right)}$

Wherein, i=1,2 ..., 3n.

All the data set Z of normalization data can be by following equation 8 expressions.

＜equation 8 〉

$X=\left[\begin{array}{cccc}{\mathrm{<figure-callout\; id="1"\; label="Z"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">Z</figure-callout>}}_{\mathrm{1,1}}& {\mathrm{<figure-callout\; id="1"\; label="Z"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">\; <figure-callout\; id="2"\; label="Z"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">Z</figure-callout>\; </figure-callout>}}_{\mathrm{1,2}}& \mathrm{\&Lambda;}& {\mathrm{<figure-callout\; id="1"\; label="X"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">X</figure-callout>}}_{1,\mathrm{<figure-callout\; id="2"\; label="m\; X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">m</figure-callout>}}& \\ X_{}{}_{\mathrm{2,1}}& {\mathrm{<figure-callout\; id="2"\; label="X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">\; <figure-callout\; id="2"\; label="X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">X</figure-callout>\; </figure-callout>}}_{\mathrm{2,2}}& \mathrm{\&Lambda;}& {\mathrm{<figure-callout\; id="2"\; label="X"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">X</figure-callout>}}_{2,m}\\ M& M& O& M\\ Z_{3n,1}& {\mathrm{<figure-callout\; id="2"\; label="Z"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">Z</figure-callout>}}_{3n,2}& \mathrm{\&Lambda;}& Z_{3n,m}\end{array}\right]=\left[\begin{array}{cccc}{Z}_1& {\mathrm{<figure-callout\; id="2"\; label="Z"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">Z</figure-callout>}}_2& \mathrm{\&Lambda;}& Z_m\end{array}\right]$

The operation that separates among the operation S300 is divided into three parts with described data set Z, as training, optimization and test.Shown in following equation 9, in this embodiment of the present invention, this three partial data is called Z _Tr, Z _OpAnd Z _Ts, its size is n * m.

＜equation 9 〉

Z _ts＝[Z _3i+1，1?Z _3i+1，2?Λ?Z _3i+1，m]

Z _tr＝[Z _3i+2，1?Z _3i+2，2?Λ?Z _3i+2，m]

Z _op＝[Z _3i+3，1?Z _3i+3，2?Λ?Z _3i+3，m]

Wherein, i=0,1,2 ..., n-1.

The leaching process of major component extracts each normalization data group Z among the operation S400 _Tr, Z _OpAnd Z _TsMajor component.The dispersion degree of major component (being the eigenwert of covariance matrix) is arranged according to its size, from largest percentage dispersion degree value, selects about Z _Tr, Z _OpAnd Z _TsMajor component P _Tr, P _OpAnd P _Ts, up to its accumulative total and reach 99.5%.

Obtain the method for major component

Principal component analysis (PCA) (PCA) is with the effective ways of a plurality of input variables by the several variablees of linear transformation boil down to.Variable after the described compression is called major component.Described PCA utilize correlativity between the variable with the data radiation of original dimension for its quadratic sum the low dimension lineoid during for maximum.

The process of extracting major component is described below.

If being tieed up input variable, m is called x ₁, x ₂..., x _m, the new variables that produces by its linear combination is called θ ₁, θ ₂..., θ _m, its relation can be by equation 10 and equation 11 expressions.

＜equation 10 〉

θ ₁＝q ₁₁Z ₁+q ₁₂z ₂+...+q _1mz _m

θ ₂＝q ₂₁z ₁+q ₂₂z ₂+...+q _2mz _m

θ _m＝q _m1z ₁+q _m2z ₂+...+q _mmz _m

＜equation 11 〉

Θ＝QZ

$\mathrm{\&Theta;}=\left[\begin{array}{c}{\mathrm{\&theta;}}_1\\ {\mathrm{\&theta;}}_2\\ M\\ {\mathrm{\&theta;}}_m\end{array}\right]$ $Q=\left[\begin{array}{cccc}{\mathrm{<figure-callout\; id="1"\; label="q"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">q</figure-callout>}}_{\mathrm{1,1}}& {\mathrm{<figure-callout\; id="2"\; label="q"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">\; <figure-callout\; id="1"\; label="q"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">q</figure-callout>\; </figure-callout>}}_{\mathrm{2,1}}& \mathrm{\&Lambda;}& q_{m,1}\\ {\mathrm{<figure-callout\; id="1"\; label="q"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">q</figure-callout>}}_{\mathrm{1,2}}& {\mathrm{<figure-callout\; id="2"\; label="q"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">\; <figure-callout\; id="2"\; label="q"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">q</figure-callout>\; </figure-callout>}}_{\mathrm{2,2}}& \mathrm{\&Lambda;}& q_{m,2}\\ M& M& \mathrm{<figure-callout\; id="1"\; label="O\; M\; q"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">O</figure-callout>}& M\\ q_{}{}_{1,\mathrm{<figure-callout\; id="2"\; label="m\; q"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">m</figure-callout>}}& q_{}{}_{2,m}& \mathrm{\&Lambda;}& q_{n,m}\end{array}\right]=\left[\begin{array}{cccc}{q}_1& {\mathrm{<figure-callout\; id="2"\; label="q"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">q</figure-callout>}}_2& \mathrm{\&Lambda;}& q_m\end{array}\right]$ $Z=\left[\begin{array}{c}{z}_1\\ z_2\\ M\\ z_m\end{array}\right]$

In this case, θ ₁, θ ₂..., θ _mThe major component that is called linear system.For convenience's sake, specifying first principal component is most important major component.The maximum that described most important major component has been described input variable changes, i.e. the major component of maximum dispersion degree.Described linear transformation Q is used for determining to satisfy the condition of following equation 12 and 13.

＜equation 12 〉

${\mathrm{<figure-callout\; id="1"\; label="q\; i"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">q</figure-callout>}}_i^1+{\mathrm{<figure-callout\; id="2"\; label="q\; i"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">q</figure-callout>}}_i^2+...+{\mathrm{<figure-callout\; id="2"\; label="q\; im"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">q</figure-callout>}}_{\mathrm{im}}^{}=1$ for?i＝1，2，…，m

＜equation 13 〉

q _i1q _j1+q _i2q _j2+…+q _imq _im＝0?for??i≠j

The condition of equation 12 remains unchanged after conversion, and the condition of equation 13 is in the relevance of having removed after the conversion between the variable.

Described major component can obtain by following process.

A. from each data set Z _Tr, Z _OpAnd Z _TsIn deduct the mean value of each variable, and be known as matrix A.Utilize described data set Z in this embodiment _TrThe embodiment of book is described as an illustration, and is expressed as following equation 14.

＜equation 14 〉

$A=Z_{\mathrm{tr}}-\stackrel{\&OverBar;}{Z_{\mathrm{tr}}}$

B. utilize equation 15 to obtain A by 17 ^TThe eigenvalue of A and the singular value S of A.Press descending sort by the eigenvalue except that 0 that equation 15 obtains, be called λ ₁, λ ₂..., λ _m

＜equation 15 〉

|A ^TA-λI|＝0

＜equation 16 〉

${\mathrm{<figure-callout\; id="1"\; label="s"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">s</figure-callout>}}_1=\sqrt{{\mathrm{\&lambda;}}_1},{\mathrm{<figure-callout\; id="2"\; label="s"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">s</figure-callout>}}_2=\sqrt{{\mathrm{\&lambda;}}_2},\mathrm{\&Lambda;},s_m=\sqrt{{\mathrm{\&lambda;}}_m},({\mathrm{\&lambda;}}_1\&GreaterEqual;{\mathrm{\&lambda;}}_2\&GreaterEqual;\mathrm{\&Lambda;}\&GreaterEqual;{\mathrm{\&lambda;}}_m)$

＜equation 17 〉

$S=\left[\begin{array}{cccc}{\mathrm{<figure-callout\; id="1"\; label="S"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">S</figure-callout>}}_1& 0& \mathrm{\&Lambda;}& 0\\ 0& {\mathrm{<figure-callout\; id="2"\; label="S"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">S</figure-callout>}}_2& \mathrm{\&Lambda;}& 0\\ \mathrm{<figure-callout\; id="0"\; label="M\; M\; O"\; filenames="H2009101779708E0000071.png,H2009101779708E0000091.png"\; state="\{\{state\}\}">M</figure-callout>}& M& O& 0\\ 0& 0& \mathrm{\&Lambda;}& S_m\end{array}\right]$

C. obtain AA ^TProper vector, promptly n * n rank matrix obtains unitary matrix U then.Utilize equation 18 to obtain eigenvalue, λ is brought into obtains in the equation 19 then corresponding to each eigenwert e ₁N * 1 rank proper vector e ₁, e ₂..., e _m

＜equation 18 〉

|AA ^T-λI|＝0

＜equation 19 〉

(AA ^T-λI)X＝0

D. utilize equation 20 to obtain the dispersion degree of each major component.

＜equation 20 〉

${\mathrm{\&sigma;}}_p={\left(\frac{\left[\begin{array}{cccc}{S}_1& {\mathrm{<figure-callout\; id="2"\; label="S"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">S</figure-callout>}}_2& \mathrm{\&Lambda;}& S_m\end{array}\right]}{\sqrt{n-1}}\right)}^2$

E. utilize equation 21 and 22 with the dispersion degree of each major component divided by the dispersion degree of whole major components and obtain number percent.

＜equation 21 〉

${\mathrm{\&sigma;}}_{p\_\mathrm{tot}}=\mathrm{sum}{\left(\frac{\left[\begin{array}{cccc}{S}_1& {\mathrm{<figure-callout\; id="2"\; label="S"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">S</figure-callout>}}_2& \mathrm{\&Lambda;}& S_m\end{array}\right]}{\sqrt{n-1}}\right)}^2$

＜equation 22 〉

$\%{\mathrm{\&sigma;}}_p=\left(\frac{{\mathrm{\&sigma;}}_p}{{\mathrm{\&sigma;}}_{p\_\mathrm{tot}}}\right)\&times;100$

F. from largest percentage dispersion degree % σ _pBeginning is calculated by carrying out accumulation, and the p quantity of selecting major component is until reaching required number percent dispersion degree (as 99.98%).

G. utilize equation 23 to extract major component.

＜equation 23 〉

Ptr＝[S ₁e ₁?S ₂e ₂?Λ?S _pe _p]

H. utilize above-mentioned identical process to extract about Z _OpAnd Z _TsMajor component.

Table 1 illustration the dispersion degree of the major component of extracting.

Table 1

??No.	??PC?Var	??Cum	??Cum％
				??1	??0.70234	??0.70234	??84.12％
??2	??0.07859	??0.78093	??93.54％

??No.	??PC?Var	??Cum	??Cum％
				??3	??0.0.02905	??0.80999	??97.02％
??4	??0.01818	??0.82816	??99.20％
				??5	??0.00357	??0.83173	??99.62％
??6	??0.00226	??0.83400	??99.89％
				??7	??0.00071	??0.83470	??99.98％
??8	??0.00011	??0.83481	??99.99％
				??9	??0.00004	??0.83485	??100.00％
??10	??0.00002	??0.83486	??100.00％
				??11	??0.00001	??0.83488	??100.00％

According to this embodiment of the present invention, utilize 7 major components.Reference table 1 when using described 7 major components, can be explained 99.9% of whole dispersion degree.Therefore, because to give up the information loss that causes of residue major component only be 0.02%.

The SVR modeling

M is tieed up input variable x ₁, x ₂..., x _mBoil down to p dimension major component θ ₁, θ ₂..., θ _m, be expressed as following equation 24.

θ ₁＝q ₁₁x ₁+q ₁₂x ₂+q _1mx _m

θ ₂＝q ₂₁x ₁+q ₂₂x ₂+q _2mx _m

...

θ _p＝q _p1x ₁+q _p2x ₂+q _pmx _m

Wherein p is the integer that is equal to or less than m.

Can be about the optimization tropic (ORL) that the k time output obtains by following equation 25 expressions as SVR.

＜equation 25 〉

$f_k\left(\mathrm{\&theta;}\right)=w_K^T\mathrm{\&theta;}+b_k$

Wherein k is 1,2 ..., m.

If will be about the k time output variable y ^(k)ε insensitive loss function definition be following equation 26, be used for obtaining about y ^(k)The optimization equation of ORL can be expressed as following equation 27.

＜equation 26 〉

＜equation 27 〉

$\mathrm{Minimize\&Phi;}(w_k,{\mathrm{\&xi;}}_k)=\frac{1}{2}{w}_k^T{w}_k+C_k\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}({\mathrm{\&xi;}}_{k,i}+{\mathrm{\&xi;}}_{k,i}^{*})$

$s.t.y_i^{\left(k\right)}-w_k^T{\mathrm{\&theta;}}_i-b\&le;{\mathrm{\&epsiv;}}_k+{\mathrm{\&xi;}}_{k,i}$

$w_k^T{\mathrm{\&theta;}}_i+b-y_i^{\left(k\right)}\&le;{\mathrm{\&epsiv;}}_k+{\mathrm{\&xi;}}_{k,i}^{*}$

${\mathrm{\&epsiv;}}_k,{\mathrm{\&xi;}}_{k,i},{\mathrm{\&xi;}}_{k,i}^{*}\&GreaterEqual;0,i=\mathrm{1,2},...,n$

Wherein the k in the equation 26 and 27 is 1,2 ..., m.ξ among Fig. 4 _KiAnd ξ _Ki ^*Be slack variable.Fig. 4 illustration the synoptic diagram of the ORL that produces by SVR.θ herein _iFor about x, corresponding to the major component vector of the i time observation vector, not i the component of θ.

Described optimization problem can utilize following equation 28 to be expressed as dual problem.

＜equation 28 〉

$\max {\mathrm{\&lambda;}}_k,{\mathrm{\&lambda;}}_k^{*}\{-\frac{1}{2}{\mathrm{\&Sigma;}}_{i=1}^n{\mathrm{\&Sigma;}}_{j=1}^n({\mathrm{\&lambda;}}_{k,i}-{\mathrm{\&lambda;}}_{k,j}^{*}){\mathrm{\&theta;}}_i^T{\mathrm{\&theta;}}_j+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}[{\mathrm{\&lambda;}}_{k,i}(y_i^{\left(k\right)}-{\mathrm{\&epsiv;}}_k)-{\mathrm{\&lambda;}}_{k,j}^{*}(y_i^{\left(k\right)}-{\mathrm{\&epsiv;}}_k)\left]\right\}$

$s.t.0\&le;{\mathrm{\&lambda;}}_{k,i},{\mathrm{\&lambda;}}_{k,j}^{*}\&le;C_k$ for?i＝1，2，…，n

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}({\mathrm{\&lambda;}}_{k,i}-{\mathrm{\&lambda;}}_{k,j}^{*})=0$

Wherein k is 1,2 ..., m.

Correspondingly, with Lagrange multiplier λ _{K, 1}And λ ^* _KiBring in the equation 29 to determine ORL about the k time output variable associating support vector regression (AASVR) certainly.

＜equation 29 〉

$f_k\left(\mathrm{\&theta;}\right)={w^{*}}_k^T\mathrm{\&theta;}+b_k^{*}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}({\mathrm{\&lambda;}}_{k,i}-{\mathrm{\&lambda;}}_{k,j}^{*}){\mathrm{\&theta;}}_i^T\mathrm{\&theta;}+b_k^{*}$

Foregoing is the step that obtains the equation of linear regression of optimization.If the Nonlinear Mapping result from the raw data to the higher dimensional space is called vectorial Φ (), the function that is defined as the equation 30 of Φ () (being the result of Nonlinear Mapping) inner product is called nuclear.

＜equation 30 〉

When attempting in higher dimensional space, to seek ORL, needn't know Φ (x _i) and Φ (x _j), only need know kernel function K (x _i, x _j) get final product.Use Gaussian radial basis function in this embodiment.

When using kernel function K (θ _i, θ)=Φ (θ _i) ^TDuring Φ (θ), can obtain as optimization Nonlinear regression equation as described in the following equation 31.

＜equation 31 〉

$f_k\left(0\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}({\mathrm{\&lambda;}}_{k,i}-{\mathrm{\&lambda;}}_{k,j}^{*})K({\mathrm{\&theta;}}_i,\mathrm{\&theta;})+b_k^{*}$

Utilize θ by following equation 32 herein, _rAnd θ _s(being any support vector) calculates bias term.

＜equation 32 〉

$b_k^{*}=-\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}({\mathrm{\&lambda;}}_{k,i}-{\mathrm{\&lambda;}}_{k,j}^{*})[K({\mathrm{\&theta;}}_i,{\mathrm{\&theta;}}_r)+K({\mathrm{\&theta;}}_i,{\mathrm{\&theta;}}_s)]$

For the optimization Nonlinear regression equation that is obtained,, can provide the wide σ of nucleus band in advance if the constant ε of loss function, the penalties C and the radial basis function (RBF) of dual objective function are used as nuclear.Can obtain total m number by repeating these processes about the SVR of each output, and the AASVR in can design of graphics 3.

Optimize the calculating of constant and use the response surface method to obtain SVR Model Optimization constant (ε, C and σ) in operation S500, described model has minimized optimization data Z _OpThe predicted value error.

A. with σ, ε and C (being described SVR model parameter) are appointed as v respectively ₁, v ₂And v ₃

B. determine with respect to each v ₁, v ₂And v ₃The hunting zone.Can obtain suitable hunting zone by existing experiment and small-scale preliminary experiment.In this embodiment, v ₁=0.56486～1.63514, v ₂=0.010534～0.010534, v ₃=2.1076～7.9933.

C. the upper and lower bound of specified search range is U respectively ₁, U ₂, U ₃And L ₁, L ₂, L ₃, utilize following equation 33 normalized mode shape parameters.

＜equation 33 〉

$x_1=\frac{v_1-(({\mathrm{<figure-callout\; id="1"\; label="U"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">U</figure-callout>}}_1+L_1)/2)}{(U_1-L_1)/2},$ $x_2=\frac{v_2-(({\mathrm{<figure-callout\; id="2"\; label="U"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">U</figure-callout>}}_2+L_2)/2)}{(U_2-L_2)/2},$ $x_3=\frac{v_3-((U_3+L_3)/2)}{(U_3-L_3)/2}$

D. experimental point, promptly the evaluation point of model performance is by considering normalized mode shape parameter x1, the hunting zone of x2 and x3 is determined.For this reason, can utilize center combination design (being a kind of of planning of statistical experiment) to finish.If the experimental point by center combination design appointment represents that with three dimensions it will be as shown in Figure 5.

E. the experimental point of center combination design comprises point on eight summits, central point and six axles.For the size of estimating experiment error, operation repeated experiments five times.The coordinate of point is defined as α=2 on the axle ^3/4=1.68179.

＜equation 34 〉

α=[number of factor experimental point] 1/4

If operate repeated experiments five times at described central point, the x1 of center combination design, the experimental point of x2 and x3 is as shown in table 2 below.

Table 2

??No	??X 1	??X 2	??X 3
				??1	??-1	??-1	??-1
??2	??1	??-1	??-1
				??3	??-1	??1	??-1
??4	??1	??1	??-1
				??5	??-1	??-1	??1
??6	??1	??-1	??1

??No	??X 1	??X 2	??X 3
				??7	??-1	??1	??1
??8	??1	??1	??1
				??9	??-1.68179	??0	??0
??10	??1.68179	??0	??0
				??11	??0	??-1.68179	??0
??12	??0	??1.68179	??0
				??13	??0	??0	??-1.68179
??14	??0	??0	??1.68179
				??15	??0	??0	??0
??16	??0	??0	??0
				??17	??0	??0	??0
??18	??0	??0	??0
				??19	??0	??0	??0

F. as shown in table 2, determine model parameter v1, the value of v2 and v3 obtains table 3 (v1 of center combination design, the experimental point of v2 and v3).This is the value that will be directly used in the model parameter that obtains model.

Table 3

??No	??V 1	??V 2	??V 3
				??1	??0.564855	??0.010534	??2.1067
??2	??1.635144	??0.010534	??2.1067
				??3	??0.564856	??0.039967	??2.1067

??No	??V 1	??V 2	??V 3
				??4	??1.635144	??0.039967	??2.1067
??5	??0.564856	??0.010534	??7.9933
				??6	??1.635144	??0.010534	??7.9933
??7	??0.564856	??0.039967	??7.9933
				??8	??1.635144	??0.039967	??7.9933
??9	??0.2	??0.02525	??5.05
				??10	??2	??0.02525	??5.05
??11	??1.1	??0.0005	??5.05
				??12	??1.1	??0.05	??5.05
??13	??1.1	??0.02525	??.01
				??14	??1.1	??0.02525	??10
??15	??1.1	??0.02525	??5.05
				??16	??1.1	??0.02525	??5.05
??17	??1.1	??0.02525	??5.05
				??18	??1.1	??0.02525	??5.05
??19	??1.1	??0.02525	??5.05

G. on each experimental point of table 3, utilize Z _TrAnd P _TrObtain the β vector and the deviation constant of SVR model.Utilize data set Z for this reason _TrIn fact, obtain identical model in sequence number 15 to sequence number 19 corresponding to central point.

H. import data set P _OpTo the m number of AASVR, to estimate the degree of accuracy of each model, the normalization predicted value of the data that are optimized then

Thus, can calculate the degree of accuracy of output model (being MSE) by equation 35.Because five central points have identical model, with P _OpBe divided into five parts, obtain independent MSE result of calculation by the subdata group after cutting apart.

＜equation 35 〉

Wherein, z _IjBe P _OpThe j time input data of middle sensor i,

It is the estimated value that obtains by model.Utilize P _TrObtain described model.As shown in table 4 below by the MSE result of calculation (being experimental result) that experiment obtains.

Table 4

??No	??MSE	??Log(MSE)
			??1	??0.000126	??-8.97923
??2	??0.000100	??-9.21034
			??3	??0.000552	??-7.50196
??4	??0.000538	??-7.52765
			??5	??0.000121	??-9.01972
??6	??0.000094	??-9.27222
			??7	??0.000552	??-7.50196
??8	??0.000541	??-7.52209
			??9	??0.001083	??-6.82802
??10	??0.000263	??-8.24336
			??11	??0.000072	??-9.53884
??12	??0.000829	??-7.09529
			??13	??0.000542	??-7.52024
??14	??0.000265	??-8.23578
			??15	??0.000221	??-8.41735
??16	??0.000204	??-8.49739
			??17	??0.000265	??-8.23578

??No	??MSE	??Log(MSE)
			??18	??0.000294	??-8.13193
??19	??0.000194	??-8.54765

I. utilize log (MSE) but not MSE meets with a response face.Given this, assessment models parameter x 1, x2, the response surface between x3 and the log (MSE).Described response surface is assumed to the 2D model as following equation 36.

＜equation 36 〉

$\log \left(\mathrm{MSE}\right)={\mathrm{\&beta;}}_0+{\mathrm{\&beta;}}_1{x}_1+{\mathrm{\&beta;}}_2{x}_2+{\mathrm{\&beta;}}_3{x}_3+{\mathrm{\&beta;}}_{11}{x}_1^2+{\mathrm{\&beta;}}_{22}{x}_2^2+{\mathrm{\&beta;}}_{33}{x}_3^2+$

${\mathrm{\&beta;}}_{12}{x}_1{x}_2+{\mathrm{\&beta;}}_{13}{x}_1{x}_4+{\mathrm{\&beta;}}_{23}{x}_2{x}_3+e$

Wherein e is a stochastic error.

The response surface that described embodiment is estimated is as follows:

$\log \left(\mathrm{MSE}\right)=-8.3492b_0-0.2131x_1-0.7716x_2-0.0952x_3$

$-0.2010{\mathrm{<figure-callout\; id="1"\; label="x"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">x</figure-callout>}}_1^2-0.0753{\mathrm{<figure-callout\; id="2"\; label="x"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">x</figure-callout>}}_2^2+0.0799x_3^2$

J. the response surface of utilization estimation obtains minimizing log (MSE), the v1 of e, the optimal conditions of v2 and v3.Because the two-dimentional response surface of supposition is determined described optimal conditions by local derviation.That is, be met the x1 of following equation 37, x2 and x3.

＜equation 37 〉

$\frac{\&PartialD;\log \left(\mathrm{MSE}\right)}{\&PartialD;x_1}=0,$ $\frac{\&PartialD;\log \left(\mathrm{MSE}\right)}{\&PartialD;x_2}=0,$ $\frac{\&PartialD;\log \left(\mathrm{MSE}\right)}{\&PartialD;x_3}=0$

Obtaining by this embodiment under the situation of described response surface, described optimal conditions is (x*1, x*2, and x*3)=(1.5438 ,-0.6818,1.5929).

K. utilize following equation 38 with described optimal conditions x*1, x*2 and x*3 are converted into initial cell.

＜equation 38 〉

${\mathrm{<figure-callout\; id="1"\; label="v"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">v</figure-callout>}}_1=\left(\frac{U_1-{\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">L</figure-callout>}}_1}{2}\right)x_1+\left(\frac{{\mathrm{<figure-callout\; id="1"\; label="U"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">U</figure-callout>}}_1+{\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">L</figure-callout>}}_1}{2}\right),$ ${\mathrm{<figure-callout\; id="2"\; label="v"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">v</figure-callout>}}_2=\left(\frac{U_2-{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">L</figure-callout>}}_2}{2}\right)x_2+\left(\frac{{\mathrm{<figure-callout\; id="2"\; label="U"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">U</figure-callout>}}_2+{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">L</figure-callout>}}_2}{2}\right),$ $v_3=\left(\frac{U_3-L_3}{2}\right)x_3+\left(\frac{U_3+L_3}{2}\right)$

According to this embodiment, v*1=1.3910=σ *, v*2=0.0005=ε *, v*3=6.7951=C*.Because described prediction log (MSE) is-9.9446 under this condition, if it is with index and is converted into MSE, it becomes 0.000048.

The generation of training pattern receives described three and optimizes constant ε * in operation S600, the major component P of C* and σ *, the training data that obtains from operation S500 _Tr, and the training data (Z of first signal _TrFirst hurdle), utilize quadratic equation to separate the optimization equation then.Then, the generation of training data obtains β among the operation S600 ₁(n * 1) (as, the difference of Lagrange multiplier) and deviation constant b ₁, to produce the SVR of Fig. 3 ₁Model.By identical method, on the 2nd to m sensor, repeat above-mentioned steps and obtain β ₂, β ₃..., β _m, and b ₂, b ₃..., b _m, like this, by SVR ₂Produce SVR _mModel is so that make up SVR model corresponding to as shown in Figure 3 whole sensor.

The prediction of the output of training data utilizes the major component P of training data among the operation S700 _TrMajor component P with test data _TsWith the kernel function matrix K that obtains Gaussian radial basis function (n * n).Then, the prediction of the output of training data utilizes β among the operation S700 ₁(as, the difference of the Lagrange multiplier of the SVR model that obtains among the operation S600) and deviation constant b ₁To produce SVR ₁Output.By identical method, on the 2nd to m sensor, repeat above-mentioned steps to obtain model predication value, promptly from SVR ₂To SVR _mOutput.

It is by following equation 39 expressions.

＜equation 39 〉

＜carry out reverse normalization 〉

The predicted value that reverse normalized execution will operate the normalization experimental data that obtains among the S700 among the operation S800 oppositely is normalized to initial range, so that the predicted value of each sensor of acquisition initial value scope.It is by following equation 40 expressions.

＜equation 40 〉

${\hat{X}}_{\mathrm{<figure-callout\; id="1"\; label="ts"\; filenames="H2009101779708E0000171.png,H2009101779708E0000181.png"\; state="\{\{state\}\}">ts</figure-callout>}1}={\hat{Z}}_{\mathrm{ts}1}\{\max \left(X_{\mathrm{ts}1}\right)-\min \left(X_{\mathrm{ts}1}\right)\}+\min \left(X_{\mathrm{ts}1}\right)$

${\hat{X}}_{\mathrm{<figure-callout\; id="2"\; label="ts"\; filenames="H2009101779708E0000091.png,H2009101779708E0000101.png"\; state="\{\{state\}\}">ts</figure-callout>}2}={\hat{Z}}_{\mathrm{ts}2}\{\max \left(X_{\mathrm{ts}2}\right)-\min \left(X_{\mathrm{ts}2}\right)\}+\min \left(X_{\mathrm{ts}2}\right)$

.?????????.

.?????????.

.?????????.

${\hat{X}}_{\mathrm{tsm}}={\hat{Z}}_{\mathrm{tsm}}\{\max \left(X_{\mathrm{tsm}}\right)-\min \left(X_{\mathrm{tsm}}\right)\}+\min \left(X_{\mathrm{tsm}}\right)$

In order to confirm superiority according to the Forecasting Methodology that is used for shop equipment monitoring performance of embodiment of the present invention, utilize the device signal data that described Forecasting Methodology and existing method are compared, described data by the unloading phase nuclear power plant first and second the circulation measure.The described data that are used to analyze are by the value of 11 sensor measurements altogether.

Table 5 illustration traditional core homing method and being used to according to embodiments of the present invention monitor predicted value between the method for shop equipment ratio of precision.

Table 5

The data that are used to analyze among Fig. 5 following value of 11 sensor measurements altogether of serving as reasons.

1. reactor capability (%)

2. pressurizer level (%)

3. the steam flow of steam generator (kg/hr)

4. the horizontal data in a narrow margin (%) of steam generator

5. pressure data (the kg/cm of steam generator ²)

6. the wide cut horizontal data (%) of steam generator

7. the main feed flow data on flows (Mkg/hr) of steam generator

8. turbine output data (MWe)

9. reactor liquid coolant charging data on flows (m ³/ hr)

10. remove the data on flows (m of waste heat ³/ hr)

11. the temperature data of reactor roof liquid coolant (℃)

If described data are used to draw function of time figure, be Fig. 6 by 16 these figure.

Fig. 6 A and Fig. 6 B illustration the power data of reactor core in the nuclear power plant.Fig. 6 A represents the test input data X of equation 6 _Ts1, Fig. 6 B represents to utilize the predicted data that obtains according to the algorithm of importing the embodiment of data about the test of equation 40

Fig. 7 A and Fig. 7 B illustration be used for pressurizer horizontal data measuring accuracy, nuclear power plant according to an embodiment of the present invention.Fig. 7 A represents the test input data X of equation 6 _Ts2, Fig. 7 B represents to utilize the data estimator that obtains according to the algorithm predicts about the embodiment of the input data of test equation 40

Fig. 8 A and Fig. 8 B illustration be used for the steam flow data of steam generator accuracy of detection, nuclear power plant according to an embodiment of the present invention.Fig. 8 A represents the test input data X of equation 6 _Ts3, the data estimator that Fig. 8 B represents to utilize the algorithm predicts according to the embodiment relevant with the input data of test equation 40 to obtain

Fig. 9 A and Fig. 9 B illustration be used for the horizontal data in a narrow margin of steam generator accuracy of detection, nuclear power plant according to an embodiment of the present invention.Fig. 9 A represents the test input data X of equation 6 _Ts4, the data estimator that Fig. 9 B represents to utilize the algorithm predicts according to the embodiment relevant with the input data of test equation 40 to obtain

Figure 10 A and Figure 10 B illustration be used for the pressure data of steam generator accuracy of detection, nuclear power plant according to an embodiment of the present invention.Figure 10 A represents the test input data X of equation 6 _Ts5, the data estimator that Figure 10 B represents to utilize the algorithm predicts according to the embodiment relevant with the input data of test equation 40 to obtain

Figure 11 A and Figure 11 B illustration be used for the wide cut horizontal data of steam generator accuracy of detection, nuclear power plant according to an embodiment of the present invention.Figure 11 A represents the test input data X of equation 6 _Ts6, the data estimator that Figure 11 B represents to utilize the algorithm predicts according to the embodiment relevant with the input data of test equation 40 to obtain

Figure 12 A and Figure 12 B illustration be used for the main feed flow data on flows of steam generator accuracy of detection, nuclear power plant according to an embodiment of the present invention.Figure 12 A represents the test input data X of equation 6 _Ts7, the data estimator that Figure 12 B represents to utilize the algorithm predicts according to the embodiment relevant with the input data of test equation 40 to obtain

Figure 13 A and Figure 13 B illustration be used for turbine output data accuracy of detection, nuclear power plant according to an embodiment of the present invention.Figure 13 A represents the test input data X of equation 6 _Ts8, the data estimator that Figure 13 B represents to utilize the algorithm predicts according to the embodiment relevant with the input data of test equation 40 to obtain

Figure 14 A and Figure 14 B illustration be used for first cycle charging data on flows accuracy of detection, nuclear power plant according to an embodiment of the present invention.Figure 14 A represents the test input data X of equation 6 _Ts9, the data estimator that Figure 14 B represents to utilize the algorithm predicts according to the embodiment relevant with the input data of test equation 40 to obtain

Figure 15 A and Figure 15 B illustration be used for the data on flows that removes waste heat accuracy of detection, nuclear power plant according to an embodiment of the present invention.Figure 15 A represents the test input data X of equation 6 _Ts10, the data estimator that Figure 15 B represents to utilize the algorithm predicts according to the embodiment relevant with the input data of test equation 40 to obtain

Figure 16 A and Figure 16 B illustration be used for the temperature data of reactor roof liquid coolant accuracy of detection, nuclear power plant according to an embodiment of the present invention.Figure 16 A represents the test input data X of equation 6 _Ts11, the data estimator that Figure 16 B represents to utilize the algorithm predicts according to the embodiment relevant with the input data of test equation 40 to obtain

In table 5, when forecast model was applied to operational monitoring, precision was the most basic barometer.Described precision is represented with the mean square deviation of model predication value and actual measured value usually.

The precision of an equipment of equation 41 expressions.

＜equation 41 〉

Wherein N is the number of experimental data,

Be the estimated value of the model of the i time experimental data, x _iIt is the measured value of the i time experimental data.

The Forecasting Methodology that is used for the performance monitoring of shop equipment is according to embodiments of the present invention extracted the major component of device signal, and the data that are used to optimize obtain SVR Model Optimization constant by the response surface method, and utilize the training data training pattern.Therefore, compare with existing Kernel regression method, the precision of calculating predicted value can improve.

Embodiment of the present invention provide a kind of Forecasting Methodology that is used for the shop equipment performance monitoring, and by the performance of on-line system monitoring shop equipment, described embodiment applies to the safety monitoring passage of nuclear power plant to this method in factory's operation.Therefore, thus the fault of equipment can obtain the reliability that real-time monitoring improves equipment.In addition, the equipment calibration cycle of nuclear power plant increases to maximum 8 years by present 18 months fuel replacement period, so that correction expense and reduce at the correction technician's of radiation area irradiation.Unnecessary number of corrections has been avoided owing to stop transport suddenly in the generating plant that inaccurate correction causes by reducing.Stop transport the cycle by the prevention of shortening factory, improve the utilization factor of factory.

Disclosed typical embodiment herein, although used specific term, its only at large, use on the describing significance and explain, be not purpose in order to limit.Correspondingly, those of ordinary skill in the art will understand, and can carry out the conversion of various forms and details under the prerequisite that does not break away from the listed spirit and scope of the invention of following claim.

Claims (19)

1. Forecasting Methodology that is used for monitoring performance of power plant instruments, this method comprises:

In matrix, show total data;

Described total data is normalized to data set;

Described normalized data set three is divided into three data sets, is used for training, optimizes and test;

Each of the data set of described normalization and three minutes is extracted major component;

Utilize the response surface method to calculate support vector regression (SVR) Model Optimization constant, the predicted value error of the data that are used to optimize with optimization;

Utilize described optimization constant to produce the SVR training pattern;

Use normalized test data to obtain the kernel function matrix, and predict the output valve of described SVR model as input; And

Given normalization output oppositely is normalized to initial range, to obtain the predicted value of variable.

2. Forecasting Methodology as claimed in claim 1, wherein, the demonstration of the described total data in the matrix is represented by following equation:

$X=\left[\begin{array}{cccc}{X}_{\mathrm{1,1}}& X_{\mathrm{1,2}}& \mathrm{\&Lambda;}& X_{1,m}\\ X_{\mathrm{2,1}}& X_{\mathrm{2,2}}& \mathrm{\&Lambda;}& X_{2,m}\\ M& M& O& M\\ X_{3n,1}& X_{3n,2}& \mathrm{\&Lambda;}& X_{3n,m}\end{array}\right]=\left[\begin{array}{cccc}{X}_1& X_2& \mathrm{\&Lambda;}& X_m\end{array}\right]$

X _ts＝[X _3i+1，1?X _3i+1，2?Λ，X _3i+1，m]＝[X _ts1?X _ts2?Λ?X _tsm]

X _tr＝[X _3i+2，1?X _3i+2，2?Λ??X _3i+2，m]＝[X _tr1?X _tr2?Λ?X _trm]

X _op＝[X _3i+3，1?X _3i+3，2?Λ??X _3i+3，m]＝[X _op1?X _op2?Λ?X _opm]

Wherein X is whole data set, the data set that Xtr, Xop and Xts are respectively the data set that is used to train, the data set that is used to optimize and are used to test, and 3n is the aggregate of data, m is the number of equipment.

3. Forecasting Methodology as claimed in claim 1, wherein, the normalization of described total data is undertaken by following equation:

$Z_i=\frac{X_i-\min \left(X_i\right)}{\max \left(X_i\right)-\min \left(X_i\right)}$

I=1 wherein, 2 ..., 3n.

4. Forecasting Methodology as claimed in claim 1, wherein, described normalized data set was undertaken by following n * m dimension equation in three minutes, and described normalization data group is called Z, and the data set that will be used to train is called Z _Tr, the data set that will be used to optimize is called Z _Op, the data set that will be used to test is called Z _Ts,

Z _ts＝[Z _3i+1，1?Z _3i+1，2?Λ?Z _3i+1，m]

Z _tr＝[Z _3i+2，1?Z _3i+2，2?Λ?Z _3i+2，m]

Z _op＝[Z _3i+3，1?Z _3i+3，2?Λ?Z _3i+3，m]

I=0 wherein, 1,2 ..., n-1.

5. Forecasting Methodology as claimed in claim 1, wherein, in the leaching process of described major component, the dispersion degree of described major component (being the eigenwert of covariance matrix) is arranged according to its size, from the major component of the number percent dispersion degree value of maximum, select about Z _Tr, Z _OpAnd Z _TsMajor component P _Tr, P _OpAnd P _TsUp to accumulative total with above 99.5%, to extract each normalization data group Z _Tr, Z _OpAnd Z _TsMajor component.

6. Forecasting Methodology as claimed in claim 5, wherein, from each data set Z _Tr, Z _Op, and Z _TsIn deduct the matrix that the mean value of each variable obtains and be called matrix A, and represent by following equation:

$A=Z_{\mathrm{tr}}-\stackrel{\&OverBar;}{Z_{\mathrm{tr}}}.$

7. Forecasting Methodology as claimed in claim 6, wherein:

Obtain A by following equation ^TThe eigenvalue of A and the singular value S of A, and

Eigenvalue except that 0 is by descending sort, and the eigenvalue after these arrangements is called λ ₁, λ ₂..., λ _m,

|A ^TA-λI|＝0

$s_1=\sqrt{{\mathrm{\&lambda;}}_1},s_2=\sqrt{{\mathrm{\&lambda;}}_2},\mathrm{\&Lambda;},s_m=\sqrt{{\mathrm{\&lambda;}}_m},({\mathrm{\&lambda;}}_1\&GreaterEqual;{\mathrm{\&lambda;}}_2\&GreaterEqual;\mathrm{\&Lambda;}{\&GreaterEqual;\mathrm{\&lambda;}}_m)$

$S=\left[\begin{array}{cccc}{S}_1& 0& \mathrm{\&Lambda;}& 0\\ 0& S_2& \mathrm{\&Lambda;}& 0\\ M& M& O& 0\\ 0& 0& \mathrm{\&Lambda;}& S_m\end{array}\right].$

8. Forecasting Methodology as claimed in claim 7, wherein:

Obtain A ^T(that is, n * n matrix) proper vector obtains unitary matrix U to A then; And

Obtain eigenvalue by following equation, then λ is brought into following equation to obtain n * 1 rank proper vector e about each eigenvalue ₁, e ₂..., e _m,

|AA ^T-λI|＝0

(AA ^T-λI)X＝0。

9. Forecasting Methodology as claimed in claim 8 wherein, obtains the dispersion degree of each major component by following equation:

${\mathrm{\&sigma;}}_p={\left(\frac{\left[\begin{array}{cccc}{S}_1& S_2& \mathrm{\&Lambda;}& S_m\end{array}\right]}{\sqrt{n-1}}\right)}^2.$

10. Forecasting Methodology as claimed in claim 9, wherein, by following equation with the dispersion degree of each major component divided by the dispersion degree of whole major components and obtain number percent:

${\mathrm{\&sigma;}}_{p\_\mathrm{tot}}=\mathrm{sum}{\left(\frac{\left[\begin{array}{cccc}{S}_1& S_2& \mathrm{\&Lambda;}& S_m\end{array}\right]}{\sqrt{n-1}}\right)}^2$

$\%{\mathrm{\&sigma;}}_p=\left(\frac{{\mathrm{\&sigma;}}_p}{{\mathrm{\&sigma;}}_{p\_\mathrm{tot}}}\right)\&times;100.$

11. Forecasting Methodology as claimed in claim 10, wherein, from the number percent dispersion degree % σ of maximum _pBeginning is calculated by carrying out accumulation, and the p quantity of selecting major component is until reaching required number percent dispersion degree (as 99.98%).

12. Forecasting Methodology as claimed in claim 11 wherein, is extracted major component by following equation:

Ptr＝[S ₁e ₁?S ₂e ₂?Λ?S _pe _p]。

13. Forecasting Methodology as claimed in claim 12, wherein, by the input variable x of following equation with the m dimension ₁, x ₂..., x _mThe major component θ of boil down to p dimension ₁, θ ₂..., θ _m:

θ ₁＝q ₁₁x ₁+q ₁₂x ₂+q _1mx _m

θ ₂＝q ₂₁x ₁+q ₂₂x ₂+q _2mx _m

θ _p＝q _p1x ₁+q _p2x ₂+q _pmx _m

Wherein p is the integer that is equal to or less than m.

14. Forecasting Methodology as claimed in claim 13 wherein, is represented by following equation about the optimization tropic (ORL) as SVR that the k time output obtains:

$f_k\left(\mathrm{\&theta;}\right)=w_K^t\mathrm{\&theta;}+b_k$

Wherein k is 1,2 ..., m.

15. Forecasting Methodology as claimed in claim 14, wherein, when defining about the k time output variable y in order to following equation ^(k)ε insensitive loss function the time, be used for obtaining about y ^(k)The optimization equation of ORL represent by following equation:

$\mathrm{Minimize\&Phi;}(w_k,{\mathrm{\&xi;}}_k)=\frac{1}{2}{w}_k^T{w}_k+C_k\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}({\mathrm{\&xi;}}_{k,i}+{\mathrm{\&xi;}}_{k,i}^{*})$

$s.t.y_i^{\left(k\right)}-w_k^T{\mathrm{\&theta;}}_i-b\&le;{\mathrm{\&epsiv;}}_k+{\mathrm{\&xi;}}_{k,i}$

$w_k^T{\mathrm{\&theta;}}_i+b-y_i^{\left(k\right)}\&le;{\mathrm{\&epsiv;}}_k+{\mathrm{\&xi;}}_{k,i}^{*}$

${\mathrm{\&epsiv;}}_k,{\mathrm{\&xi;}}_{k,i},{\mathrm{\&xi;}}_{k,i}^{*}\&GreaterEqual;0$ for?i＝1，2，…，n

Wherein k is 1,2 ..., m, ξ _KiAnd ξ _Ki ^*Be slack variable.

16. Forecasting Methodology as claimed in claim 15, wherein, the following equation of described optimization problem utilization is expressed as dual problem:

$\max {\mathrm{\&lambda;}}_k,{\mathrm{\&lambda;}}_k^{*}\{-\frac{1}{2}{\mathrm{\&Sigma;}}_{i=1}^n{\mathrm{\&Sigma;}}_{j=1}^n({\mathrm{\&lambda;}}_{k,i}-{\mathrm{\&lambda;}}_{k,j}^{*}){\mathrm{\&theta;}}_i^T{\mathrm{\&theta;}}_j+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}[{\mathrm{\&lambda;}}_{k,i}(y_i^{\left(k\right)}-{\mathrm{\&epsiv;}}_k)-{\mathrm{\&lambda;}}_{k,j}^{*}(y_i^{\left(k\right)}-{\mathrm{\&epsiv;}}_k)\left]\right\}$

$s.t.0\&le;{\mathrm{\&lambda;}}_{k,i},{\mathrm{\&lambda;}}_{k,j}^{*}\&le;C_k$ for?i＝1，2，…，n

$\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}({\mathrm{\&lambda;}}_{k,i}-{\mathrm{\&lambda;}}_{k,j}^{*})=0$

Wherein k is 1,2 ..., m.

17. Forecasting Methodology as claimed in claim 16, wherein, with Lagrange multiplier λ _{K, 1}And λ ^* _KiIn the following equation of substitution to determine ORL about the k time output variable associating support vector regression (AASVR) certainly:

$f_k\left(\mathrm{\&theta;}\right)={w^{*}}_k^T\mathrm{\&theta;}+b_k^{*}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}({\mathrm{\&lambda;}}_{k,i}-{\mathrm{\&lambda;}}_{k,j}^{*}){\mathrm{\&theta;}}_i^T\mathrm{\&theta;}+b_k^{*}.$

18. Forecasting Methodology as claimed in claim 17, wherein, if the Nonlinear Mapping result from the raw data to the higher dimensional space is called vectorial Φ (), the kernel function that the optimization non-linear regression line utilization of following equation is defined as the inner product (being the result of Nonlinear Mapping) of Φ () obtains:

$f_k\left(0\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}({\mathrm{\&lambda;}}_{k,i}-{\mathrm{\&lambda;}}_{k,j}^{*})K({\mathrm{\&theta;}}_i,\mathrm{\&theta;})+b_k^{*}.$

19. Forecasting Methodology as claimed in claim 17 wherein, is utilized θ by following equation _rAnd θ _s(being any support vector) calculates bias term:

$b_k^{*}=-\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}({\mathrm{\&lambda;}}_{k,i}-{\mathrm{\&lambda;}}_{k,j}^{*})[K({\mathrm{\&theta;}}_i,{\mathrm{\&theta;}}_r)+K({\mathrm{\&theta;}}_i,{\mathrm{\&theta;}}_s)].$

Images