CN104698976B - A kind of depth diagnostic method of predictive control model hydraulic performance decline - Google Patents

Abstract

The invention discloses a kind of depth diagnostic method of predictive control model hydraulic performance decline, using production run data calculation process disturbing signal, then creation data and the step-response coefficients computation model predicated error in Design of Predictive stage are utilized, the model performance index judgment models overall performance constructed by the two is good and bad.For the model of penalty, new model performance index further is calculated with the method for removal input variable one by one, the corresponding submodel performance of input of removal is judged by the situation of change of performance indications, so as to realize being monitored each submodel performance.The present invention is merely with production process data and design data, assessment can not only be provided to umlti-variable finite elements model overall performance, the corresponding submodel performance of all inputs can be more estimated, suggestion is given for engineer is controlled system maintenance, predictive controller maintenance cost can be greatly reduced.

Description

A kind of depth diagnostic method of predictive control model hydraulic performance decline

Technical field

The invention belongs to manufacturing forecast control system performance monitoring field, it is related to based on model performance index and removes input The forecast model monitoring technology of variable depth inspection.

Background technology

Because production process dynamic characteristic changes with the time, industrial model predictive control (MPC) system is generally being thrown Fortune was accomplished by being safeguarded that maintenance process needs object to model again, cycle high cost long less than 1 year, and generally needed life Product process is stopped, and causes very big economic loss.Moreover, judge whether to be necessary before maintenance to model whole object again, Or model again and whether determine that whole control performance can be improved, effective analysis method is lacked at present.In order to avoid building again The unnecessary input that mould and controller are safeguarded, it is necessary to study a kind of method to judge where department pattern performance is occurred in that down object Drop, and any degree have decreased to.

Whether control performance monitoring can need maintenance to provide suggestion MPC to enterprise engineering teachers.Control performance is monitored Conventional method be minimum variance performance indications that Harris is proposed, actual output variance and theoretical minimum variance are compared Compared with so as to judge actual controller performance.Subsequently multi-input multi-output system, time-varying system, Constrained system are occurred in that in succession The minimum variance Performance Evaluation index of system.Because the minimum variance assessment technology of multi-input multi-output system needs to calculate association square Battle array, and the calculating of incidence matrix generally requires to obtain relatively difficult time lag and impulse response coefficient.To avoid incidence matrix meter Calculate, occur in that some more practical evaluation indexes, such as Linear-Quadratic-Gauss benchmark, user-defined benchmark, design objective and The contrast benchmark of actual index, history top performance benchmark, etc..

Control performance monitoring technology gives the method for assessing whole control system performance, but, control performance declines can Can decline due to controller performance, forecast model is mismatched, new interference appearance and many factors such as actuator is non-linear.Closely Many research work concentrate on the reason for declining to controller performance and diagnose over year.Such as document Xuemin Tian, GongquanChen,Sheng Chen.A data-based approach for multivariate model Predictive control performance monitoring.Neurocomputing 74 (2011) 588-597, propose The species that control performance declines is screened using the sorting technique of data-driven, document Keck Voon Ling, Weng Khuen Ho,Yong Feng,and Bingfang Wu.Integral-Square-Error Performance of Multiplexed Model Predictive Control.IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS,2011,7(2):196-203, it is proposed that the ISE indexs of controller hydraulic performance decline in itself are judged, so as to recognize Go out whether control system hydraulic performance decline is caused by controller.Also there are some research works for the diagnosis of forecast model hydraulic performance decline Make, such as by design control baseline error judgment models performance (such as Abhijit S.Badwe, Rohit S.Patwardhan, Sirish L.Shah,Sachin C.Patwardhan,Ravindra D.Gudi.Quantifying the impact of model-plant mismatch on controller performance.Journal of Process Control 20 (2010) 408-425), the process frequency by testing responds benchmark decision model performance (G.Ji, K.Zhang, Y.Zhu.A method of MPC model error detection.Journal of Process Control,22(2012)635– 642) method such as contrast decision model performance, or by forecast model surplus and actual interference amount (Zhijie Sun, S.Joe Qin,Ashish Singhal,Larry Megan.Performance monitoring of model-predictive controllers via model residual assessment.Journal of Process Control,23(2013) 473-482)。

For the assessment of forecast model performance, there is problems with the studies above work：(1) actual industrial process is mostly The multi-variable system of multiple-input and multiple-output (MIMO), the model of corresponding each input variable of each output variable is that a list is defeated Enter single output (SISO) model, current research work only gives the assessment result whether whole model performance declines, but work Cheng Shi it is interested be which model performance of specific which variable declines and result in control performance decline, so as in control Workload is reduced during system maintenance, it is cost-effective；(2) definition of performance indications generally requires the process mechanism knowledge for being difficult to obtain Or to carry out plant characteristic test, it is difficult to promoted to practical application.Accordingly, it would be desirable to a kind of only using production process information without shadow Ring production process, can on-line analysis go out model overall performance whether decline and diagnose specifically which submodel hydraulic performance decline Method, for engineers provide maintenance suggestion before system maintenance, reduce unnecessary maintenance cost.

The content of the invention

It is an object of the invention to provide a kind of method of predictive control model performance monitoring, the method was only using producing Journey service data is normally run without influence production process, can provide model overall performance evaluation index, and in globality Specific which submodel hydraulic performance decline and decline degree are further diagnosed to be when can decline, are that the maintenance of Predictive Control System subtracts Few workload, it is cost-effective.

In order to solve the above technical problems, the present invention provides a kind of depth diagnostic method of predictive control model hydraulic performance decline, Comprise the following steps：

Step 1：For certain multiple-input and multiple-output process, k moment output variable y (k) corresponds to N_uIndividual input variable u_j (k) (j=1,2 ..., N_u), e^oK () is process disturbance signal, p is time window length, and m is output variable order, and n becomes for input Amount order, defines output vector y in time window p_p(k), input vector u_jp(k) (j=1,2 ..., N_u), process disturbance signal to AmountAnd the joint vector constructed by last time input/output variable

Step 2：By output vector y_pK vector is combined in () and input and outputCalculating process disturbing signal sequence vector

Step 3：Using input variable u_jThe conventional operation data and forecast model coefficient a of (k), output variable y (k)_ij, point Not Ji Suan k to the k-p moment model predictive error, obtain model predictive error sequence e_p(k)=[e (k) e (k-1) ... e (k-p)]；

Step 4：Process disturbance signal sequence [e is obtained by step 2^o(k) e^o(k-1) ... e^o(k-p)] obtained with step 3 The model predictive error sequence [e (k) e (k-1) ... e (k-p)] for arriving, computation model performance indications η_R；

Step 5：Comparison model performance indications η_RWith predetermined model performance indicator threshold value η_R0If, model performance Index η_RHigher than the threshold value, then it is assumed that the model performance of the output variable is good, otherwise makes l=1, goes to step 6；

Step 6：Remove l-th input variable u_l, calculate corresponding second process disturbance signal

Step 7：Remove l-th input variable u_l, second model predictive error at k to k-p moment is calculated respectively, obtain the Two model predictive error sequence e_pl=[e_l(k) e_l(k-1) ... e_l(k-p)]；

Step 8：The the second process disturbance signal obtained by step 6 The the second model predictive error sequence [e obtained with step 7_l(k) e_l(k-1) ... e_l(k-p) the second model performance], is calculated Index η_RlAnd performance indications compare κ_l；

Step 9：Compare κ according to performance indications_lJudge input variable u_l(k) corresponding submodel performance, if κ_l>1, represent Input variable u_l(k) corresponding submodel penalty, it is necessary to safeguard, if instead κ_l<1, represent that corresponding submodel performance is good It is good；

Step 10：Make l=l+1, return to step 6, until l=N_u, all of input processing finishes；

Step 11：Return to step 1 processes next output variable y (k), until the assessment of all of output variable is finished.

Further, in the step 1, vector process in k moment time windows p is constructed as follows：

y_p(k)=[y (k) y (k-1) ... y (k-p)]

Wherein j=1 ... N_u

u_jp(k)=[u_j(k) u_j(k-1) ... u_j(k-p)]

Wherein, Y_m(k-1) by the y at k-1 moment to k-m moment_pK () constructs, U_n(k-1) by the k-1 moment to k-n moment u_jpK () constructs, the order of m, n difference output variable and input variable.

Further, in the step 2, process disturbance signal vectorCalculating process is as follows：

DefinitionThe orthocomplement, orthogonal complement projection in row space：

Wherein I is (m+n*N_u) dimension unit matrix.

When p is sufficiently large, is projected by orthocomplement, orthogonal complement and calculate process disturbance signal vector

Calculated to simplify, it is rightImplement LQ to decompose

Wherein vector Q₁, Q₂It is orthogonal,The inferior triangular flap obtained after being decomposed for LQ.

Then calculating process disturbing signal is vectorial

Wherein symbol is whereinRepresent broad sense inverse operation.

Further, in the step 3, the process for calculating k moment model predictive error e (k) is as follows：

Computation model prediction output y_m(k)：

Wherein N₀It is forecast model length, N_uIt is the input variable number corresponding to output variable y, a_ji(i=1 ... N_u, j= 1,...,N₀) it is to determine in the Design of Predictive stage for j-th i-th step-response coefficients of input variable, Δ u_j (k-i) it is variation delta u between two neighbouring sample moment of j-th input variable_j(k-i)=u_j(k-i)-u_j(k-i-1)；

According to output variable y (k), model prediction output y_m(k) computation model predicated error

E (k)=(1-q^-1)(y(k)-y_m(k)) (7)

Wherein q^-1It is One-step delay operator.

Further, in the step 4, model performance index η_RCalculating process is as follows：

Further, calculated in step 6 and remove input variable u_lCorresponding second process disturbance signal afterwardsProcess is as follows：

Remove input variable u_l, similar formula (1) reconfigures vectorial in time window：

y_p(k)=[y (k) y (k-1) ... y (k-p)]

Wherein j=1 ..., l-1, l+1 ... N_u

u_jp(k)=[u_j(k) u_j(k-1) ... u_j(k-p)]

Wherein T symbols represent transposition.It is then carried out LQ decomposition

It is wherein vectorialWithIt is orthogonal,The inferior triangular flap obtained after being decomposed for LQ.

Calculate the process disturbance signal after removing input variable ul

Further, in the step 7, input variable u is removed_lAfterwards, the model predictive error e at k moment is calculated_l(k) process It is as follows：

Remove input variable u_l, a is removed in formula (6) Section 2_liΔu_l(k-i) item calculates y_ml(k)

The second model predictive error e is calculated by formula (13)_l(k)：

e_l(k)=(1-q^-1)(y(k)-y_ml(k)) (13)

Further, in the step 8, the second model performance index η_RlAnd performance indications compare κ_lCalculating process is as follows：

The beneficial effect that the present invention is reached：The present invention provides a kind of predictive control model deep monitored for engineers Method, the method need not interrupt production process, it is not required that technique priori, only using production process data to pre- observing and controlling The performance of block mold processed and each submodel is estimated, so that judge where department pattern performance occurs in that decline to object, and Any degree is have decreased to, makes whether engineers pre- perception model before Predictive Control System maintenance needs to safeguard, and which portion Dividing needs to safeguard, saves unnecessary modeling cost again.

Brief description of the drawings

Fig. 1 is model performance index method schematic diagram；

Fig. 2 is depth diagnosis submodel performance methodology schematic diagram；

Fig. 3 is Wood-Berry destilling tower experimental result pictures.

Specific embodiment

The present invention is described in detail with reference to the accompanying drawings and detailed description.

The technical solution adopted in the present invention principle is as shown in Figure 1.

Consider the Predictive Control System shown in Fig. 1, wherein G^o(q) and H^oQ () represents process model and process interference mould respectively Type, G_mQ () and H (q) represent predictive control model and predicted interference model, G respectively_cQ () is predictive controller, e^o(k) and d^o K () is process interference signal and disturbance quantity, e (k) and d (k) is predicated error and Controlling model error, and r (k) is reference locus, u K () is control input (performance variable), y (k) is output variable (controlled variable),For prediction is exported.

It is assumed that Controlling model and Disturbance Model are matched, it is clear that predicated error e (k) should be equal to process disturbance signal e^o (k)., whereas if there is model mismatch, because predicated error e (k) includes model mismatch information, process disturbance signal is should be greater than e^o(k), therefore Definition Model performance indications

Wherein N is assessment time domain data length, e^oK () can be by output variable y (k) and the conventional fortune of input variable u (k) Row data are obtained, and e (k) can be obtained by model step-response coefficients and service data.η_RScope for (0,1], if η_RIt is close Show that model performance is good in 1, η_RShow that model performance deteriorates close to zero.

Performance indications η_RAssessment can be given to block mold performance, but for multiple-input and multiple-output process, forecast model is It is made up of the submodel of multiple single-input single-outputs, the decline of model overall performance does not necessarily mean that all of submodel Can all decline.Therefore, for the model of hydraulic performance decline, the present invention proposes the method for further diagnosing each submodel performance, such as Shown in Fig. 2.

It is assumed that judging for l-th input u_lSubmodel performance, due to being difficult to that correspondence is obtained from output variable y (k) In u_lComponent, therefore for u_lSubmodel performance cannot be by η_RCalculating formula be directly calculated.Side proposed by the present invention Method is to remove input variable u one by one_lNew model performance index is calculated afterwards, if model performance indications improve after removal, explanation The corresponding submodel poor-performing of input variable of removal is, it is necessary to be safeguarded.

Due to input variable u_lEffect to exporting cannot be removed from output variable y (k), therefore u_lEffect from defeated Enter passage and be transferred to disturbance passage, so all input channels and disturbance passage are constant to the resultant action for exporting, y (k) still can be with For calculating e^o(k), the e being so calculated^oK () includes removed u_lEffect.During calculating e (k) also from forecast model Removal u_lContinuous item, the also effect comprising l-th removed submodel in result of calculation e (k).

Implement step as follows：

Step 1：For certain multiple-input and multiple-output process, k moment output variable y (k) corresponds to N_uIndividual input variable u_j (k) (j=1,2 ..., N_u), e^oK () is process disturbance signal, p is time window length, and m is output variable order, and n becomes for input Amount order, defines output vector y in time window p_p(k), input vector u_jp(k) (j=1,2 ..., N_u), process disturbance signal to AmountAnd the joint vector constructed by last time input/output variable

Step 2：By output vector y_pK vector is combined in () and input and outputCalculating process disturbing signal sequence vector

Step 3：Using input variable u_jThe conventional operation data and forecast model coefficient a of (k), output variable y (k)_ij, point Not Ji Suan k to the k-p moment model predictive error, obtain model predictive error sequence e_p(k)=[e (k) e (k-1) ... e (k-p)]；

Step 4：Process disturbance signal sequence [e is obtained by step 2^o(k) e^o(k-1) ... e^o(k-p)] obtained with step 3 The model predictive error sequence [e (k) e (k-1) ... e (k-p)] for arriving, computation model performance indications η_R；

Step 5：Comparison model performance indications η_RWith predetermined model performance indicator threshold value η_R0If, model performance Index η_RHigher than the threshold value, then it is assumed that the model performance of the output variable is good, otherwise makes l=1, goes to step 6；

Step 6：Remove l-th input variable u_l, calculate corresponding second process disturbance signal

Step 7：Remove l-th input variable u_l, second model predictive error at k to k-p moment is calculated respectively, obtain the Two model predictive error sequence e_pl=[e_l(k) e_l(k-1) ... e_l(k-p)]；

Step 8：The the second process disturbance signal obtained by step 6 The the second model predictive error sequence [e obtained with step 7_l(k) e_l(k-1) ... e_l(k-p) the second model performance], is calculated Index η_RlAnd performance indications compare κ_l；

Step 9：Compare κ according to performance indications_lJudge input variable u_l(k) corresponding submodel performance, if κ_l>1, represent Input variable u_l(k) corresponding submodel penalty, it is necessary to safeguard, if instead κ_l<1, represent that corresponding submodel performance is good It is good；

Step 10：Make l=l+1, return to step 6, until l=N_u, all of input processing finishes；

Step 11：Return to step 1 processes next output variable y (k), until the assessment of all of output variable is finished.

Further, in the step 1, vector process in k moment time windows p is constructed as follows：

y_p(k)=[y (k) y (k-1) ... y (k-p)]

Wherein j=1 ... N_u

u_jp(k)=[u_j(k) u_j(k-1) ... u_j(k-p)]

Wherein, Y_m(k-1) by the y at k-1 moment to k-m moment_pK () constructs, U_n(k-1) by the k-1 moment to k-n moment u_jpK () constructs.

Further, in the step 2, process disturbance signal vectorCalculating process is as follows：

DefinitionThe orthocomplement, orthogonal complement projection in row space：

Wherein I is (m+n*N_u) dimension unit matrix.

When p is sufficiently large, is projected by orthocomplement, orthogonal complement and calculate process disturbance signal vector

Calculated to simplify, it is rightImplement LQ to decompose

Wherein vector Q₁, Q₂It is orthogonal,The inferior triangular flap obtained after being decomposed for LQ.

Then calculating process disturbing signal is vectorial

Wherein symbol is whereinRepresent broad sense inverse operation.

Further, in the step 3, the process for calculating k moment model predictive error e (k) is as follows：

Computation model prediction output y_m(k)：

Wherein N₀It is forecast model length, N_uIt is the input variable number corresponding to output variable y, a_ji(i=1 ... N_u, j= 1,...,N₀) it is to determine in the Design of Predictive stage for j-th i-th step-response coefficients of input variable, Δ u_j (k-i) it is variation delta u between two neighbouring sample moment of j-th input variable_j(k-i)=u_j(k-i)-u_j(k-i-1)；

According to output variable y (k), model prediction output y_m(k) computation model predicated error

E (k)=(1-q^-1)(y(k)-y_m(k)) (7)

Wherein q^-1It is One-step delay operator.

Further, in the step 4, model performance index η_RCalculating process is as follows：

Further, calculated in step 6 and remove input variable u_lCorresponding second process disturbance signal afterwardsProcess is as follows：

Remove input variable u_l, similar formula (1) reconfigures vectorial in time window：

y_p(k)=[y (k) y (k-1) ... y (k-p)]

Wherein j=1 ..., l-1, l+1 ... N_u

u_jp(k)=[u_j(k) u_j(k-1) ... u_j(k-p)]

Wherein T symbols represent transposition.It is then carried out LQ decomposition

It is wherein vectorialWithIt is orthogonal,The inferior triangular flap obtained after being decomposed for LQ.

Calculate and remove input variable u_lProcess disturbance signal afterwards

Further, in the step 7, input variable u is removed_lAfterwards, the model predictive error e at k moment is calculated_l(k) process It is as follows：

Remove input variable u_l, a is removed in formula (6) Section 2_liΔu_l(k-i) item calculates y_ml(k)

The second model predictive error e is calculated by formula (13)_l(k)：

e_l(k)=(1-q^-1)(y(k)-y_ml(k)) (13)

Further, in the step 8, the second model performance index η_RlAnd performance indications compare κ_lCalculating process is as follows：

Experiment simulation and analysis：

Forecast model depth diagnostic method proposed by the present invention has carried out emulation experiment in Wood-Berry destilling towers.The mistake Transfer function matrix G (s) of journey comes from bibliography R.K.Wood, M.W.Berry, Terminal composition control of a binary distillation column,Chemical Engineering Science28(1973) In 1707-1717.

Capacity of returns and the input variable (performance variable) that steam flow is the process, are designated as u respectively₁And u₂, unit lb/ Min, tower top and bottom product component are two controlled variables, and y is designated as respectively₁And y₂, unit mol%.It is assumed that the sampling period is 1 Minute, process transfer function matrix G after discretization^oQ () is

It is assumed that actual interference process model H^oK () is

The disturbing signal e of applying^oK () is that variance is diag { 0.744²,0.13²Independent white noise.

MPC predicts that time domain and control time domain elect 100 and 10 as respectively.Weight matrix elects Q=diag { 1,10 }, S=diag as {1,10}.Two output variable setting values are respectively

The validity of research method in order to verify, using different object models and Disturbance Model, walks according to preceding method Suddenly tested, model output order m and input order n take 30, and time window length p takes 3940.

Situation 1：Accurate object model and Disturbance Model；

Situation 2：Accurate object model, Disturbance Model is the model form that DMC is commonly used

Situation 3：Object model and the equal mismatch of Disturbance Model, process model are elected as

Disturbance Model is (19) formula；

Situation 4：Object model and the equal mismatch of Disturbance Model, process model are elected as

Disturbance Model is (19) formula；

Situation 5：First input variable u is removed in situation 4₁, diagnose y₂Model；

Situation 6：Second input variable u is removed in situation 4₂, diagnose y₂Model；

Corresponding MQI values are shown in Table 1, and corresponding histogram compares sees Fig. 2.

The Wood-Berry destilling tower model evaluation results of table 1

Be can be seen that from table 1 and Fig. 2：

(1) in situation 1 during Model Matching, two MQI indexs of controlled variable show that MQI indexs being capable of essence all close to 1 Really reflect model performance.

If hydraulic performance decline (such as situation 2, situation 3, situation 4) occur in process model or Disturbance Model, MQI indexs are occur Decline, and letdown procedure is corresponding with unmatched models degree.Such as the MQI2 of situation 3 and 4, the gain of situation 3 is 1.2 times of realities Border target gain, the gain of situation 4 is then 3 times of practical object gains, and corresponding MQI2 situations 3 are 0.6090, and situation 4 is 0.1188, therefore the numerical values recited of MQI can reflect the mismatch degree of model.

(2) the 4th sub- model gain is 3 times of actual gains in situation 4, therefore model is seriously mismatched.It is calculated MQI2 be 0.1188, much smaller than 1, therefore can determine whether second output variable model performance degradation, this and actual conditions It is consistent.This belongs to the first stage of the inventive method, that is, judge block mold performance, can not still judge y₂Under model performance Drop is by u₁Correspondence submodel or u₂Correspondence submodel causes.

(3) situation 5 and situation 6 belong to second stage of the invention, i.e., removing input using shifting division calculation one by one becomes MQI indexs after amount, the corresponding submodel performance of variable that degree is lifted according to its index to judge to remove.Due in situation 4 Judge y₂Corresponding model performance deteriorates, and situation 5 and situation 6 calculate removal u respectively₁With removal u₂MQI indexs afterwards.

Compared with situation 4, MQI2 drops to 0.0664, corresponding κ values 0.5589 from 0.1188 in situation 5<1, show to move Except a submodel of good performance.In fact, being input into u with first in situation 4₁Related submodel is accurate submodule Type, therefore, judged result matches with actual conditions.

Conversely, compared with situation 4, MQI2 is substantially increased in situation 6,0.3727, corresponding κ values are risen to from 0.118 3.1372>1, show to remove a submodel for penalty.In fact, being input into u with second in situation 4₂Related son Model is implicitly present in serious mismatch, therefore, judged result matches with actual conditions.

The above is only the preferred embodiments of the present invention, and any formal limitation is not made to the present invention, It is every according to technical spirit of the invention to any simple modification made for any of the above embodiments, equivalent variations and modification are belonged to In the range of technical solution of the present invention.

Claims (8)

1. a kind of depth diagnostic method of predictive control model hydraulic performance decline, it is characterised in that comprise the following steps：

Step 1：For certain multiple-input and multiple-output process, k moment output variable y (k) corresponds to N_uIndividual input variable u_j(k)(j =1,2 ..., N_u), e^oK () is process disturbance signal, p is time window length, and m is output variable order, and n is input variable rank It is secondary, define output vector y in time window p_p(k), input vector u_jp(k) (j=1,2 ..., N_u), process disturbance signal vectorAnd the joint vector constructed by last time input/output variable

Step 2：By output vector y_pK vector is combined in () and input and outputCalculating process disturbing signal sequence vector

Step 3：Using input variable u_jThe conventional operation data and forecast model coefficient a of (k), output variable y (k)_ij, count respectively The model predictive error at k to k-p moment is calculated, model predictive error sequence e is obtained_p(k)=[e (k) e (k-1) ... e (k- p)]；

Step 4：Process disturbance signal sequence [e is obtained by step 2^o(k) e^o(k-1) ... e^o(k-p)] obtained with step 3 Model predictive error sequence [e (k) e (k-1) ... e (k-p)], computation model performance indications η_R；

Model performance index η_RComputing formula is

${\mathrm{\&eta;}}_R=\frac{\underset{i=0}{\overset{p}{\&Sigma;}}{\&lsqb;e^o(k-i)\&rsqb;}^2}{\underset{i=0}{\overset{p}{\&Sigma;}}{\&lsqb;e(k-i)\&rsqb;}^2}---\left(8\right);$

Step 5：Comparison model performance indications η_RWith predetermined model performance indicator threshold value η_R0If, model performance index η_RHigher than the threshold value, then it is assumed that the model performance of the output variable is good, otherwise makes l=1, goes to step 6；

Step 6：Remove l-th input variable u_l, calculate corresponding second process disturbance signal

Step 7：Remove l-th input variable u_l, second model predictive error at k to k-p moment is calculated respectively, obtain the second mould Type predicated error sequence e_pl=[e_l(k) e_l(k-1) ... e_l(k-p)]；

Step 8：The the second process disturbance signal obtained by step 6And step Rapid 7 the second model predictive error sequence [e for obtaining_l(k) e_l(k-1) ... e_l(k-p) the second model performance index], is calculated η_RlAnd performance indications compare κ_l；

Step 9：Compare κ according to performance indications_lJudge input variable u_l(k) corresponding submodel performance, if κ_l>1, represent input Variable u_l(k) corresponding submodel penalty, it is necessary to safeguard, if instead κ_l<1, represent that corresponding submodel is functional；

Step 10：Make l=l+1, return to step 6, until l=N_u, all of input processing finishes；

Step 11：Return to step 1 processes next output variable y (k), until the assessment of all of output variable is finished.

2. the depth diagnostic method of predictive control model hydraulic performance decline according to claim 1, it is characterised in that described In step 1, vector process in k moment time windows p is constructed as follows：

Wherein j=1 ... N_u

$Y_m(k-1)=\left[\begin{array}{c}{y}_p(k-1)\\ y_p(k-2)\\ .\\ .\\ .\\ y_p(k-m)\end{array}\right],U_n(k-1)=\left[\begin{array}{c}{u}_{1p}(k-1)\\ .\\ .\\ .\\ u_{1p}(k-n)\\ .\\ .\\ .\\ u_{N_{u}p}(k-1)\\ .\\ .\\ .\\ u_{N_{u}p}(k-n)\end{array}\right],{\stackrel{\&OverBar;}{Z}}_p\left(k\right)=\left[\begin{array}{c}{Y}_m(k-1)\\ U_n(k-1)\end{array}\right]---\left(1\right)$

Wherein, Y_m(k-1) by the y at k-1 moment to k-m moment_pK () constructs, U_n(k-1) by the u at k-1 moment to k-n moment_jp(k) Construction.

3. the depth diagnostic method of predictive control model hydraulic performance decline according to claim 1, it is characterised in that described In step 2, process disturbance signal vectorCalculating process is as follows：

DefinitionThe orthocomplement, orthogonal complement projection in row space：

${\mathrm{\&Pi;}}_{{\stackrel{\&OverBar;}{Z}}_p\left(k\right)}^{\&perp;}=I-{\stackrel{\&OverBar;}{Z}}_p{\left(k\right)}^T{\&lsqb;{\stackrel{\&OverBar;}{Z}}_p\left(k\right){\stackrel{\&OverBar;}{Z}}_p{\left(k\right)}^T\&rsqb;}^{-1}{\stackrel{\&OverBar;}{Z}}_p\left(k\right)---\left(2\right)$

Wherein I is (m+n*N_u) dimension unit matrix；

Projected by orthocomplement, orthogonal complement and calculate process disturbance signal vector

$e_p^o\left(k\right)=y_p\left(k\right){\mathrm{\&Pi;}}_{{\stackrel{\&OverBar;}{z}}_p\left(k\right)}^{\&perp;}---\left(3\right).$

4. the depth diagnostic method of predictive control model hydraulic performance decline according to claim 3, it is characterised in that described In step 2, process disturbance signal vectorSimplified calculating process it is as follows：

It is rightImplement LQ to decompose

$\left[\begin{array}{c}{\stackrel{\&OverBar;}{Z}}_p\left(k\right)\\ y_p\left(k\right)\end{array}\right]=\left[\begin{array}{cc}{L}_{11}& \\ L_{21}& L_{22}\end{array}\right]\left[\begin{array}{c}{Q}_1\\ Q_2\end{array}\right]---\left(4\right)$

Wherein vector Q₁, Q₂It is orthogonal,The inferior triangular flap obtained after being decomposed for LQ,

Then calculating process disturbing signal is vectorial

Wherein symbol is whereinRepresent broad sense inverse operation.

5. the depth diagnostic method of predictive control model hydraulic performance decline according to claim 1, it is characterised in that described In step 3, the process for calculating k moment model predictive error e (k) is as follows：

Computation model prediction output y_m(k)：

$y_m\left(k\right)=\underset{j=1}{\overset{N_u}{\&Sigma;}}\underset{i=1}{\overset{N_0-1}{\&Sigma;}}{a}_{ji}{\mathrm{\&Delta;u}}_j(k-i)+\underset{j=1}{\overset{N_u}{\&Sigma;}}{a}_{{\mathrm{jN}}_0}{u}_j(k-N_0)---\left(6\right)$

Wherein N₀It is forecast model length, N_uIt is the input variable number corresponding to output variable y, a_ji(i=1 ... N_u, j= 1,...,N₀) it is to determine in the Design of Predictive stage for j-th i-th step-response coefficients of input variable, △ u_j (k-i) it is variable quantity △ u between two neighbouring sample moment of j-th input variable_j(k-i)=u_j(k-i)-u_j(k-i-1)；

According to output variable y (k), model prediction output y_m(k) computation model predicated error

E (k)=(1-q^-1)(y(k)-y_m(k)) (7)

Wherein q^-1It is One-step delay operator.

6. the depth diagnostic method of predictive control model hydraulic performance decline according to claim 1, it is characterised in that step 6 Middle calculating removes input variable u_lCorresponding second process disturbance signal afterwardsProcess is as follows：

Remove input variable u_l, similar formula (1) reconfigures vectorial in time window p：

Wherein j=1 ..., l-1, l+1 ... N_u

$Y_m(k-1)=\left[\begin{array}{c}{y}_p(k-1)\\ y_p(k-2)\\ .\\ .\\ .\\ y_p(k-m)\end{array}\right],$

$U_{\mathrm{nl}}(k-1)={\left[\begin{array}{cccccccccccccc}{u}_{1p}(k-1)& ...& u_{1p}(k-n)& ...& u_{(l-1)p}(k-1)& ...& u_{(l-1)p}(k-n)& u_{(l+1)p}(k-1)& ...& u_{(l+1)p}(k-1)& ...& u_{N_{u}p}(k-1)& ...& u_{N_{u}p}(k-n)\end{array}\right]}^T,$

${\stackrel{\&OverBar;}{Z}}_{pl}\left(k\right)=\left[\begin{array}{c}{Y}_m(k-1)\\ U_{nl}(k-1)\end{array}\right]---\left(9\right)$

Wherein symbol T represents transposition, is U with step 1 difference_nl(k-1) it is removal u_lVector afterwards, is then carried out LQ points Solution

$\left[\begin{array}{c}{\stackrel{\&OverBar;}{Z}}_{pl}\left(k\right)\\ y_p\left(k\right)\end{array}\right]=\left[\begin{array}{cc}{L}_{11}^l& \\ L_{21}^l& L_{22}^l\end{array}\right]\left[\begin{array}{c}{Q}_1^l\\ Q_2^l\end{array}\right]---\left(10\right)$

It is wherein vectorialWithIt is orthogonal,The inferior triangular flap obtained after being decomposed for LQ, calculates and removes input variable u_lAfterwards Process disturbance signal

7. the depth diagnostic method of predictive control model hydraulic performance decline according to claim 1, it is characterised in that described In step 7, input variable u is removed_lAfterwards, the model predictive error e at k moment is calculated_lK () process is as follows：

Remove input variable u_l, a is removed in formula (6) Section 2_li△u_l(k-i) item calculates y_ml(k), i.e.,

$y_{ml}\left(k\right)=\underset{j=1}{\overset{l-1}{\&Sigma;}}\underset{i=1}{\overset{N_0-1}{\&Sigma;}}{a}_{ji}{\mathrm{\&Delta;u}}_j(k-i)+\underset{j=l+1}{\overset{N_u}{\&Sigma;}}\underset{i=1}{\overset{N_0-1}{\&Sigma;}}{a}_{ji}{\mathrm{\&Delta;u}}_j(k-i)+\underset{j=1}{\overset{l-1}{\&Sigma;}}{a}_{{\mathrm{jN}}_0}{u}_j(k-N_0)+\underset{j=l+1}{\overset{N_0}{\&Sigma;}}{a}_{{\mathrm{jN}}_0}{u}_j(k-N_0)---\left(12\right)$

The second model predictive error e is calculated by formula (13)_l(k)：

e_l(k)=(1-q^-1)(y(k)-y_ml(k)) (13)。

8. the depth diagnostic method of predictive control model hydraulic performance decline according to claim 1, it is characterised in that described In step 8, the second model performance index η_RlAnd performance indications compare κ_lForm is

${\mathrm{\&eta;}}_{Rl}=\frac{\underset{i=0}{\overset{p}{\&Sigma;}}{\left(e_l^o\right(k-i\left)\right)}^2}{\underset{i=0}{\overset{p}{\&Sigma;}}{\left(e_l\right(k-i\left)\right)}^2}---\left(14\right)$

${\mathrm{\&kappa;}}_l=\frac{{\mathrm{\&eta;}}_{Rl}}{{\mathrm{\&eta;}}_R}---\left(15\right).$