CN104268350A - Closed-loop quality control simulation method with simulation prediction and actual production integrated - Google Patents

Abstract

The invention belongs to the field of error control during part machining process, and relates to a closed-loop quality control simulation method with simulation prediction and actual production integrated. The method includes the steps of 1), performing multi-operation multi-error source analysis and description of the part machining process; 2), designing a data integration real-time collecting method and device of on-side multi-operation multi-error sources; 3), performing simulation prediction on comprehensive error on the basis of a dynamic error source coupling model; 4), predicting a quality pre-control method on the basis of dynamic error simulation prediction; 5), adjusting the out-of-control process according to quality monitoring situations and realizing closed-loop quality control. Thinking bandage of an existing method is broken through during control over the multi-operation machining process of the parts, prediction of the errors during the machining process is performed firstly, and then quality pre-control and process adjustment is performed on the basis of error prediction. Conventional afterwards quality control is converted into prior pre-control of quality by the method which is ideal in quality control.

Description

The closed loop quality control emulation mode that a kind of simulation and prediction is mutually integrated with actual production

[technical field]

The present invention relates to part manufacture process closed loop quality control field, the closed loop quality control emulation mode that particularly a kind of simulation and prediction is mutually integrated with actual production.

[background technology]

Along with the requirement of market to product quality is more and more higher, enterprise more and more recognizes the importance of production quality control in actual production.Part process is subject to the impact of multiple error source, exist between upstream and downstream operation and transmit cumulative effect, every procedure is except the introducing error of this sequence node, also to consider the cumulative errors associating operation propagation, these error sources all may cause the fluctuation of part crudy, must analyze these error sources, the generation of Mass Control fluctuation from source, these error sources are constantly changes in process simultaneously, they are also constantly changes on the impact of part quality, the situation of change that Accurate Prediction just needs the practical production status considering these error sources to be carried out to part quality.And current method of quality control be by processing after examination and analysb carry out the prevention of part quality accident below, to quality control and monitoring be all control afterwards and monitor, the prevention of later quality accident can only be played, and not and practical production status immediately connect, make quality monitoring and analyze inaccurate, and being difficult to the adjustment of implementation procedure.

[summary of the invention]

The problem to be solved in the present invention be can be implemented in digitalized artificial environment by this closed-loop simulation method of quality control to predict the quality fluctuation problem of manufacture process, monitor and adjust, avoid the quality accident in actual processing, improve working (machining) efficiency, realize quality control in advance.

To achieve these goals, the technical solution adopted in the present invention is:

The closed loop quality control emulation mode that a kind of simulation and prediction is mutually integrated with actual production, comprise the following steps: be coupled under 1) error source that part process produces being transformed into infinitesimal coordinate system, then according to the error after coupling, consider that practical production status is on the impact of error source, sets up forecast model simultaneously; 2) Acquisition Error information, it can be used as the input of forecast model; 3) forecast model is converted to dynamic SOV model, sets up the OSFE equation corresponding with SOV model, obtain independent identically distributed deviation; 4) according to the control chart of independent identically distributed Deviation Design forecast model; 5) Serial regulation model and MGWMA controller is set up; 6) covariance matrix of controller is solved according to MGWMA controller; 7) make MGWMA control chart introduce a deviation on the basis in initial error source, make control chart out of control, then alignment error source, until this statistic is less than this control limit.

Described step 4) in, the method for the control chart of design forecast model is: 4.1) calculate the mean vector of independent same distribution deviation and the estimated value of covariance matrix according to SOV model; 4.2) according to the mean vector of independent same distribution deviation and the estimated value compute statistics of covariance matrix; 4.3) according to the control limit of the distributed problem solving multivariate statistical procedure of statistic obedience.

Described step 4.1) according to following formulae discovery:

${\stackrel{\‾}{\hat{e}}}_k^{d_k}={\left[\begin{array}{cccc}{\stackrel{\‾}{\hat{e}}}_{k1}^{d_k}& {\stackrel{\‾}{\hat{e}}}_{k2}^{d_k}& ...& {\stackrel{\‾}{\hat{e}}}_{\mathrm{ks}}^{d_k}\end{array}\right]}^T$

$S_k=\frac{1}{m-1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}({\hat{e}}_{\mathrm{ki}}-{\stackrel{\‾}{\hat{e}}}_k){({\hat{e}}_{\mathrm{ki}}-{\stackrel{\‾}{\hat{e}}}_k)}^T$

Wherein, what represent is that the deviation that produces of error source is at the mean vector of the independent same distribution deviation of operation k.

Described step 4.2) concrete grammar be: assuming that the vector of a new independent same distribution deviation, and and S _kindependently vectorial; The then T of kth procedure _k ²the expression formula of statistic is:

${T_k}^2={({\hat{e}}_{\mathrm{ku}}-{\stackrel{\‾}{\hat{e}}}_k)}^T{S}_k^{-1}({\hat{e}}_{\mathrm{ku}}-{\stackrel{\‾}{\hat{e}}}_k)~\frac{s(m+1)(m-1)}{m(m-s)}{F}_{s,m-s}.$

Step 4.3) concrete grammar be: assuming that the probability of error of first kind is α, then the control of multivariate statistical procedure is limited to:

${\mathrm{UCL}}_k=\frac{s(m+1)(m-1)}{m(m-s)}{F}_{s,m-s}\left(\mathrm{\α}\right).$

Described step 1) in forecast model be:

P＝A _JFMU _JFMJ _FM+B _WFMU _WFMW _FM，

Wherein, A _jFMrepresenting the correction factor of geometrical deviation source dynamic change in time, is static coefficient U _jFMjacobian matrix, B _wFMrepresent the correction factor of physical error source dynamic change in time, U _jFMrepresent the static Jacobian matrix coefficient in geometric error source, U _wFMrepresent the static Jacobian matrix coefficient in physical error source, J _fMrepresent geometrical deviation, W _fMrepresent physical deflections.

Step 3) in dynamic SOV model be:

$\left\{\begin{array}{cc}{\mathrm{\μ}}_k^{d_k,t_k,l_k,j_k}=A_{k-1}{\mathrm{\μ}}_{k-1}^{d_{k-1},t_{k-1},l_{k-1},j_{k-1}}{B}_k^{d_k}{u}_k^{d_k}+B_k^{f_k}{u}_k^{f_k}+B_k^{t_k}{u}_k^{t_k}+B_k^{j_k}{u}_k^{j_k}+w_k,w_k~N[0,W_k]& \\ P_k^{d_k}=P=C_E{E}_k^{d_k}{\mathrm{\μ}}_k+v_k,& v_k~N[0,V_k]\\ {\mathrm{\μ}}_0|D_0~N[m_0,C_0]& \end{array}\right.$

Wherein, $\left\{\begin{array}{c}{A}_{k-1}=U_{\mathrm{JFM}}{\left(C_{\mathrm{ws}}^O{U}_{\mathrm{wsi}}\right)}^{+}D\\ B_k^{f_k}=-U_{\mathrm{JFM}}{\left(C_{\mathrm{ws}}^O{U}_{\mathrm{wsi}}\right)}^{+}(N\\ B_k^{t_k}=C_T^O{U}_{\mathrm{JFM}}\\ B_k^{j_k}=C_T^{\mathrm{Ps}}R_T^O\end{array}\right.$

what represent is the coordinate conversion of measuring system, C _erepresent the correction factor that measuring system changes with practical production status.

Described linear model has m error originated from input source, p output characteristics value, and according to SOV model, the adjustment model of multiple-input and multiple-output process is described below:

${P_k}^{d_k}={\mathrm{\α}}_k+E_k^{d_k}*U_k+v_k$

Wherein, the vector of one (p × 1), the mass property that representative exports; α is one (p × 1) vector, represents the compensating parameter of each output characteristics; one (p × m) vector, system of representatives matrix number; $U_k=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\φ}}_{k,i}^{(\·)}\left(\begin{array}{cccc}{B}_i^{d_i}& B_i^{f_i}& B_i^{t_i}& I\end{array}\right)*{\left(\begin{array}{cccc}{u}_i^{d_i}& u_i^{f_i}& u_i^{t_i}& w_i\end{array}\right)}^T+{\mathrm{\φ}}_{k,0}^{(\·)}{\mathrm{\μ}}_0$ Be one (m*1) vector, represent the error source of process system; v _kthe vector of one (p × 1), the noise in expression process.

Assuming that the v in SOV model _kbe a mean vector be 0, covariance matrix is the white noise sequence of ∑, then described step 5) in MGWMA controller be described below:

$g_i=\underset{k=1}{\overset{i}{\mathrm{\Σ}}}[q^{{(i-t)}^b}-q^{{(i-t+1)}^b}]P_k,i\∈I^{+}$

G in above formula _ione (p × 1) vector (g ₀=0), design parameter q is a constant and (0≤q<1), and adjustment parameter b is artificial set-point.

The beneficial effect of the invention: adopt the present invention to carry out quality control, by predicting the predictive simulation of part processing course error, carry out quality Pre-control system on this basis, carry out quality monitoring than existing method and can accomplish quality monitoring in advance and process adjusting, the generation of quality fluctuation and quality accident can be prevented, reduce rejection rate and number of rewelding, improve production and processing efficiency, reduce production cost.

[accompanying drawing explanation]

Fig. 1 closed-loop simulation quality control overall plan.

Fig. 2 part processing course error source is analyzed.

Error modeling method under the multiple error source of Fig. 3 drives.

The formation of Fig. 4 surface of the work error.

The quality monitoring flow scheme design that Fig. 5 multi-source multi-working procedure processing course expansion SoV and MSPC is integrated.

In Fig. 6 two procedures O}, ws}, fs}, ps} with { T} coordinate system defines, and wherein, (a) is first operation, and (b) is second operation work.

The foundation of the control chart of Fig. 7 first operation.

The control chart monitored results of Fig. 8 first operation.

The foundation of the control chart of Fig. 9 second operation work.

Figure 10 second operation work control chart monitored results.

Figure 11 first operation wave process adjustment process, wherein, the variation situation of change that (a) is error source, (b) is for controller is to the result after error source adjustment.

[embodiment]

The present invention breaches existing methods constraint, proposes a kind of method of new quality control.The method does not recycle the part status after processing and carries out quality analysis, but by analyzing the mismachining tolerance of the part that the conversion of different error source obtains under digitalized artificial environment, thus calculate the multiple error source coupling predictive simulation error of part multi-working procedure processing course, and adopt the MES of design integration real-time acquisition device to transmit the practical production status information of multiple error source, revise the error source state in simulated environment, thus make predictive simulation error closer to actual conditions, then by quality Pre-control method, predictive simulation error is monitored and adjusted, thus realize the closed-loop control of multi-working procedure processing course quality.

The present invention is by MES integration real-time acquisition device collection and the error source of record and the real-time production status information of workpiece, technological design, processing planning and measure planning is carried out in digitalized artificial environment, set up the dynamic digitalized machining simulation forecast model of workpiece, compare with theoretical model, obtain quality information, by designing quality pre-control method, according to quality Pre-control simulation, judge that whether process is out of control.If process is out of control, then carry out manufacturing equipment adjustment and technique adjustment, until prediction is up-to-standard, then the operational order of simulation optimization is transferred in MES and carries out actual production, reduce trial cut, rate of reducing the number of rejects and seconds, improve working (machining) efficiency.Can the process of pre-control complex parts by the method, improve the qualification rate of product, get through last link of Digital production line.

The present invention specifically comprises the following steps:

1) carry out part process multiple operation multiple error source to analyze.

2) on-the-spot multiple operation multiple error source data integration real-time collecting method and device is designed.

3) the composition error dynamic simulation based on Dynamic Error Source coupling model is predicted.

4) based on the quality Pre-control method of dynamic error simulation and prediction.

5) according to quality monitoring situation, runaway event is adjusted, realize closed loop quality control.

Described step 1) in multiple operation multiple error source is divided into processing before in the clear and definite error source of cause-effect relationship and processing the indefinite error source of cause-effect relationship state, take operation as technique yardstick, the error modeling method considered under the multiple error source driving of practical production status is proposed, for the process system of workpiece-benchmark-cutter-fixture composition, in conjunction with its change of operation and the change of error source in process, the mode of propagation of research error, to cutter, the error that the local error sources such as fixture produce is analyzed, consider the error analysis method that geometry location error source and statics (cutting force) error source combine simultaneously, set up the error source model of workpiece-benchmark-cutter-fixture under digitalized artificial environment.

Described step 2) in the error source kind according to processing site and environmental information, adopt the vibration signal of ACC sensor Real-time Collection process, calculated the natural frequency of workpiece by model analysis simultaneously, adjust the natural frequency of the speed of mainshaft away from workpiece of lathe, or the damping increased between damping material and workpiece, evade the impact of vibration on work pieces process, AE sensor Cutter wear is adopted to monitor, obtain the abrasion condition of cutter Different periods by analysis, adopt the positioning adjusting block indicated through RFID, the thickness information of locating piece is read by rfid interrogator, and then read the change of positioning error of fixture, adopt force sensor measuring cutter X, Y, the cutting force of Z-direction over time, thus calculate average Milling Force, after these information acquisitions are arrived, change through A/D, data cleansing process, as the source of digitalized artificial environment real-time dynamic data.

Described step 3) according to part feature carry out part rigid body infinitesimal divide, the rigid body translation carrying out infinitesimal for the dynamic geometry positioning error source that cause-effect relationship before processing is clear and definite asks for the error of processing rear geometry location error source and producing, the error caused for the indefinite dynamic cutting force of cause-effect relationship in processing asks for the error after corresponding processing by the conversion of the infinitesimal rigid body of cutting force-induced error, the error that dynamic geometry positioning error source and dynamic cutting force error source produce is coupled, the coupling error after performance prediction processing.

Described step 4) in the pre-control of carrying out process quality according to error coupler predictive simulation result, on the basis of mismachining tolerance digitalized artificial prediction, multiple mass property is paid close attention to for part process, but not there is funtcional relationship between these mass propertys, error prediction value and error source is adopted to build statistic, calculate and control limit, set up multivariate process quality control figure, realize the in advance pre-monitoring of process at digitized environment.

Described step 5) in the pre-control situation according to quality, monitored results is fed back to the process optimization stage, adjust according to the error source of multivariate statistical procedure method to runaway event at operation stage, calculate statistics adjustment amount, runaway event is returned stable, realize the closed loop quality control of manufacture process.

Below in conjunction with accompanying drawing, the present invention is described in further detail:

One, Analysis of error source

Referring to shown in Fig. 2, is the error source situation that the present invention considers:

The impact of error source on part quality is divided into two kinds of situations by the present invention, one: in multiple error source variablees of performance process variation, there is cause-effect relationship, final part quality situation is obtained, as not having the rigid body geometry location error source under stressing conditions before processing by cause-effect relationship (rigid body translation); Two: the cause-effect relationship of multiple error source variable is indefinite, cannot set up Causal model, cannot by the quality condition of Causality Analysis part, as stress deformation error source in thin-walled parts process affects situation to part quality.The error source existed in the process that present patent application is considered comprises: the error that (1) positioning datum precision produces; (2) error (fixture) error produced is installed; (3) error of location attitude of the cutter change generation; (4) distortion inaccuracy of cutting force generation; (5) noise error, wherein, first three plants error for there is cause-effect relationship, and latter two is that cause-effect relationship is indefinite.

Two, the forecast model under digitalized artificial environment is set up according to error source

Refer to shown in Fig. 3: be the error modeling method of the present invention under multiple error source drives.

The present invention take operation as technique yardstick, the error modeling method considered under the multiple error source driving of practical production status is proposed, for the process system of workpiece-benchmark-cutter-fixture composition, in conjunction with its change of operation and the change of error source in process, the mode of propagation of research error, to cutter, the error that the local error sources such as fixture produce is analyzed, consider the error analysis method that geometry location error source and statics (cutting force) error source combine simultaneously, set up the error coupler forecast model of workpiece-benchmark-cutter-fixture under digitalized artificial environment.

At operation k, define 5 coordinate systems, be respectively global coordinate system { O}, workpiece coordinate system { ws}, property coordinate system { fs}, and infinitesimal coordinate system { ps} and tool coordinate system { T}.

Process step design model SM set expression: <S _mn, S _mn-1..., S _mk..., S _m1>.

Procedure technology model GM set expression: <G _mn, G _m( _n-1) ..., G _mk..., G _m1>.

Operation realistic model FM set expression: <F _mn, F _m( _n-1) ..., FM _k..., FM ₁>, the dynamic effects in geometry location error source and physical error source considered by this model.

Geometric error source operation emulation set expression: <J _fMn, J _fM( _n-1) ..., J _fMk..., J _fM1>.

Physical error source operation emulation set expression: <W _fMn, W _fM( _n-1) ..., W _fMk..., W _fM1>.

Operation total error set expression: <P _n, P _n-1..., P ₁>.

In order to carry out the prediction of error, under error source being transformed into infinitesimal coordinate system, carrying out the calculating of coupling error, considering that practical production status is on the impact of error source simultaneously, thus the forecast model of following formula (1) can be obtained:

P＝A _JFMU _JFMJ _FM+B _WFMU _WFMW _FM (1)

Wherein A _jFMrepresenting the correction factor of geometrical deviation source dynamic change in time, is static coefficient U _jFMjacobian matrix, B _wFMrepresent the correction factor of physical error source dynamic change in time, U _jFMrepresent the static Jacobian matrix coefficient in geometric error source, U _wFMrepresent the static Jacobian matrix coefficient in physical error source, J _fMrepresent geometrical deviation, W _fMrepresent physical deflections.

Wherein, geometric error and physical error are transformed into the method for infinitesimal coordinate system and are:

By the deviation situation in various dynamic geometry positioning error source under rigid body translation is transformed into infinitesimal coordinate system, the geometrical deviation obtained under infinitesimal coordinate system is:

$\begin{array}{c}{P}^I=A_{\mathrm{JFM}}{U}_{\mathrm{JFM}}{J}_{\mathrm{FM}}=A_{\mathrm{JFM}}\mathrm{\&Delta;}{X}_T^{\mathrm{ps}}={\left(H_{\mathrm{ps}}^O\right)}^{-1}{\left(H_{\mathrm{fs}}^O\right)}^{-1}{\left(H_{\mathrm{Ws}}^O\right)}^{-1}{[C_T^O{\mathrm{\&Delta;}}_T^O-C_{\mathrm{ws}}^O{\mathrm{\&Delta;}}_{\mathrm{ws}}^O]}_{6*1}\left(H_{\mathrm{ws}}^O\right)\left(H_{\mathrm{fs}}^O\right)\left(H_{\mathrm{ps}}^O\right)\\ =C_T^O{U}_{\mathrm{JFM}}{\mathrm{\&Delta;}}_T^O-U_{\mathrm{JFM}}\left[{\left(C_{\mathrm{ws}}^O{U}_{\mathrm{wsi}}\right)}^{+}(R_0^{\mathrm{li}}{}_O{\mathrm{\&Delta;}}_{\mathrm{li}}-R_0^{\mathrm{ws}}{}_O{\mathrm{\&Delta;r}}_{\mathrm{wsi}})\right]=C_T^O{U}_{\mathrm{JFM}}{\mathrm{\&Delta;}}_T^O-U_{\mathrm{JFM}}\left[{\left(C_{\mathrm{ws}}^O{U}_{\mathrm{wsi}}\right)}^{+}(\mathrm{N\&Delta;}{r}_l-\mathrm{D\&Delta;}{r}_{\mathrm{wsd}})\right]\end{array}---\left(2\right)$

In formula:

$U_{\mathrm{JFM}}=\left[\begin{array}{cc}\left(R_0^{\mathrm{ps}}{}_O\right)\left(R_0^{\mathrm{fs}}{}_O\right)\left(R_0^{\mathrm{ws}}{}_O\right)& -\left(R_0^{\mathrm{Ps}}{}_O\right)\left(R_0^{\mathrm{fs}}{}_O\right)(\left(R_0^{\mathrm{ws}}{}_O\right)\left(r\&times;_{\mathrm{ws}}^O\right)+r\&times;_{\mathrm{fs}}^O)\left(R_0^{\mathrm{ws}}{}_O\right)-\left(R_0^{\mathrm{ps}}{}_O\right)\left(r\&times;_{\mathrm{ps}}^O\right)\left(R_0^{\mathrm{fs}}{}_O\right)\left(R_0^{\mathrm{ws}}{}_O\right)\\ 0_{3*3}& \left(R_0^{\mathrm{ps}}{}_O\right)\left(R_0^{\mathrm{fs}}{}_O\right)\left(R_0^{\mathrm{ws}}{}_O\right)\end{array}\right],$

N＝diag(n ₁...n _m)∈R ^3m*m，

$D=\mathrm{diag}(R_{\mathrm{ws}1}^{\mathrm{ws}}{}_O...R_{\mathrm{wsm}}^{\mathrm{ws}}{}_O)\&Element;R^{3m*3m}.$

Wherein, represent the nominal value of infinitesimal coordinate system relative to the rotation matrix of global coordinate system, representation feature coordinate system relative to the nominal value of the rotation matrix of global coordinate system, represent the nominal value of workpiece coordinate system relative to the rotation matrix of global coordinate system, with make difficulties title matrix, for cutter is relative to the departure of global coordinate system, for workpiece is relative to the departure of global coordinate system, for feature is relative to the departure of global coordinate system, represent the departure of infinitesimal relative to global coordinate system.Δ r _lirepresent that workpiece contacts with i-th register pin the deviation of register pin, property coordinate system relative to the nominal value of the rotation matrix of global coordinate system, Δ r _wsirepresent the deviation of the workpiece that workpiece contacts with i-th register pin, Δ r _wlfor register pin deviation, Δ r _wsdfor positioning datum deviation. represent the time dependent correction factor of cutter deviation, represent the time dependent correction factor of workpiece deviation, can be obtained by actual measurement, for the position auto―control of the corresponding global coordinate system O of infinitesimal coordinate system ps, for the position auto―control of the corresponding global coordinate system O of property coordinate system fs, for the position auto―control of the corresponding global coordinate system O of workpiece coordinate system ws, represent the nominal value of i-th register pin coordinate system relative to the rotation matrix of global coordinate system.

As shown in Figure 4, in process, cutting force makes cutter and thin-wall part produce elastic deformation, and after feed, elastic deformation recovers, and causes portion of material not cut, causes the mismachining tolerance of piece surface.The component of machined surface normal direction is the principal element determining surface of the work error, so herein will based on normal component of force when calculating thin-wall part distortion.

Under under workpiece coordinate system, dynamic cutting force is transformed into infinitesimal coordinate system, being then deformed into of the infinitesimal that under infinitesimal coordinate system, cutting force causes:

$P^{\mathrm{II}}=B_{\mathrm{WFM}}{U}_{\mathrm{MFW}}{W}_{\mathrm{FM}}=B_{\mathrm{WFM}}R_T^{\mathrm{Ps}}\frac{F}{K}=C_T^{\mathrm{Ps}}R_T^{\mathrm{PS}}\frac{F}{K}---\left(3\right)$

In formula:

$U_{\mathrm{WFM}}=R_T^{\mathrm{Ps}}={\left[\begin{array}{cc}{\left(R_0^{\mathrm{ps}}{}_O\right)}^T& -{\left(R_0^{\mathrm{ps}}{}_O\right)}^{T}r_{\mathrm{ps}}^O\&times;\\ 0_{3*3}& {\left(R_0^{\mathrm{ps}}{}_O\right)}^T\end{array}\right]}_{6*6}$

Wherein, represent cutter deviation time dependent correction factor in infinitesimal coordinate system, represent the rotation matrix of tool coordinate system relative to infinitesimal coordinate system.

Three, signals collecting

Comprise with lower part and function:

1) all vibration signals of ACC sensor Real-time Collection process are adopted, calculated the natural frequency of workpiece by model analysis simultaneously, the speed of mainshaft of adjustment lathe ensures the natural frequency of natural frequency away from workpiece of machine tool chief axis, or the damping increased between damping material (damping material is arranged between workpiece and fixture) and workpiece, evades the impact of vibration on work pieces process.

2) AE sensor Cutter wear is adopted to monitor in real time.

3) positioning adjusting block (positioning adjusting block is arranged between fixture and worktable) indicated through RFID is adopted, read the thickness information of locating piece by rfid interrogator, and then read the change (namely variation in thickness) of positioning error of fixture.

4) adopt the cutting force of force sensor measuring cutter X, Y, Z-direction over time, thus calculate average Milling Force.

Above-mentioned 2), 3), 4) in three data, 2) and 3) be geometrical deviation, 4) be physical deflections.

After these information acquisitions are arrived, after A/D conversion, denoising, dimensionality reduction, just can as the source of digitalized artificial environment real-time dynamic data, as the input of forecast model geometrical deviation and physical deflections.

Four, the closed loop quality control method of forecast model

Forecast model, before carrying out closed loop quality control, first needs forecast model to be converted to dynamic SoV model, and then carries out closed-loop control, introduce respectively below.

1, forecast model is converted to dynamic SoV model

Formula (1) is write as the form of SoV, thus dynamic SoV model tormulation formula can be obtained:

$\left\{\begin{array}{cc}{\mathrm{\&mu;}}_k^{d_k,t_k,l_k,j_k}=A_{k-1}{\mathrm{\&mu;}}_{k-1}^{d_{k-1},t_{k-1},l_{k-1},j_{k-1}}{B}_k^{d_k}{u}_k^{d_k}+B_k^{f_k}{u}_k^{f_k}+B_k^{t_k}{u}_k^{t_k}+B_k^{j_k}{u}_k^{j_k}+w_k,w_k~N[0,W_k]& \\ P_k^{d_k}=P=C_E{E}_k^{d_k}{\mathrm{\&mu;}}_k+v_k,& v_k~N[0,V_k]\\ {\mathrm{\&mu;}}_0|D_0~N[m_0,C_0]& \end{array}\right.---\left(4\right)$

Wherein:

$\left\{\begin{array}{c}{A}_{k-1}=U_{\mathrm{JFM}}{\left(C_{\mathrm{ws}}^O{U}_{\mathrm{wsi}}\right)}^{+}D\\ B_k^{f_k}=-U_{\mathrm{JFM}}{\left(C_{\mathrm{ws}}^O{U}_{\mathrm{wsi}}\right)}^{+}(N\\ B_k^{t_k}=C_T^O{U}_{\mathrm{JFM}}\\ B_k^{j_k}=C_T^{\mathrm{Ps}}R_T^O\end{array}\right.$

indicate by benchmark d _k, cutter t _k, location l _k, and cutting force j _kthe state value of the error caused, represent observed reading, represent the error transmitted by kth-1 procedure, represent the error of the kth procedure caused by benchmark, represent by the error of locating the kth procedure caused, represent the error of the kth procedure caused by cutter, represent the error of the kth procedure caused by cutting force, w _krepresent process noise, v _krepresent measurement noises, what represent is the coordinate conversion of measuring system, C _erepresent the correction factor that measuring system changes with practical production status.

2, closed loop quality control

As shown in Figure 5, be the integrated monitoring flow process of quality Pre-control, namely, the integrated flow figure of SoV model and quality monitoring.

The present invention, on the basis that mismachining tolerance digitalized artificial is predicted, only pays close attention to a mass property parameter for part process, can regard the special case of multivariate control chart as, devises unit multiple operation control chart and evaluates its performance; Multiple mass property is paid close attention to for part process, but not there is funtcional relationship between these mass propertys, error prediction value and error source is adopted to build statistic, calculate and control limit, set up multivariate control chart, and carry out control chart performance to evaluate the analysis realizing performance in different deviation situation, provide the foundation for carrying out error diagnostics by the method for Statistics decomposition; For the situation between part process mass property with funtcional relationship, the profile as free curve free form surface controls, and the basis of data compression devises nonparametric profile control chart to free curve and evaluates its performance.

2.1 based on dynamic SoV model, set up the OSFE corresponding with it (One-Step ahead Forecast Error, in advance one-step prediction error equation, as follows:

$\left\{\begin{array}{c}{e}_k^{d_k}=(P_k^{d_k}-C_E{E}_k^{d_k}{\mathrm{\&mu;}}_k^{d_k,t_k,f_k,j_k})R_k^{-1/2}\\ {\mathrm{\&mu;}}_k^{d_k,t_k,f_k,j_k}=A_{k-1}^{d_{k-1}}{\mathrm{\&mu;}}_{k-1}^{d_{k-1},t_{k-1},f_{k-1},j_{k-1}}+B_k^{f_k}{u}_k^{f_k}+B_k^{t_k}{u}_k^{t_k}+B_k^{j_k}{u}_k^{j_k}+G_{k-1}{e}_{k-1}\\ R_k={\left(C_E{E}_k^{d_k}\right)}^T\left(C_E{E}_k^{d_k}\right)H_k+{\mathrm{\&Sigma;}}_{v_k}^T{\mathrm{\&Sigma;}}_{v_k}\\ H_k={\left(A_{k-1}^{d_{k-1}}\right)}^T{A}_{k-1}^{d_{k-1}}{H}_{k-1}-A_{k-1}^{d_{k-1}}{E}_k^{d_k}{H}_{k-1}^T{G}_{k-1}+{\mathrm{\&Sigma;}}_{w_k}^T{\mathrm{\&Sigma;}}_{w_k}\\ G_k=A_k^{d_k}{C}_E{D}_k^{d_k}{H}_k{R}_k^{-1}\end{array}\right.---\left(5\right)$

Wherein, initial value meets ξ is initial vector.

Independent same distribution deviation can be solved according to above formula, be designated as this independent same distribution deviation obeys standardized normal distribution.

2.2 according to independent same distribution deviation the control chart of sequences Design forecast model.

2.2.1 the mean vector of independent same distribution deviation and the estimated value of covariance matrix of operation k is calculated according to SOV model, as follows:

${\stackrel{\&OverBar;}{\hat{e}}}_k^{d_k}={\left[\begin{array}{cccc}{\stackrel{\&OverBar;}{\hat{e}}}_{k1}^{d_k}& {\stackrel{\&OverBar;}{\hat{e}}}_{k2}^{d_k}& ...& {\stackrel{\&OverBar;}{\hat{e}}}_{\mathrm{ks}}^{d_k}\end{array}\right]}^T---\left(6\right)$

$S_k=\frac{1}{m-1}\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}({\hat{e}}_{\mathrm{ki}}-{\stackrel{\&OverBar;}{\hat{e}}}_k){({\hat{e}}_{\mathrm{ki}}-{\stackrel{\&OverBar;}{\hat{e}}}_k)}^T---\left(7\right)$

Wherein what represent is that the deviation that produces of error source is at the mean vector of the independent same distribution deviation of operation k.

2.2.2 assuming that the vector of a new independent same distribution deviation, and and S _kindependently vectorial.

The then T of kth procedure _k ²the expression formula of statistic is:

${T_k}^2={({\hat{e}}_{\mathrm{ku}}-{\stackrel{\&OverBar;}{\hat{e}}}_k)}^T{S}_k^{-1}({\hat{e}}_{\mathrm{ku}}-{\stackrel{\&OverBar;}{\hat{e}}}_k)~\frac{s(m+1)(m-1)}{m(m-s)}{F}_{s,m-s}---\left(8\right)$

2.2.3 the probability supposing error of first kind (error of first kind refers to: null hypothesis is true, and is judged as vacation) is α, then the control of multivariate statistical procedure is limited to:

${\mathrm{UCL}}_k=\frac{s(m+1)(m-1)}{m(m-s)}{F}_{s,m-s}\left(\mathrm{\&alpha;}\right)---\left(9\right)$

2.3 process adjusting.

2.3.1 set up Serial regulation model

For the thin-walled parts process of multiple-input and multiple-output, the present invention devises a Serial regulation model, and this linear model has m error originated from input source, p output characteristics value, according to formula (4), the adjustment model of multiple-input and multiple-output process is described below:

$\begin{array}{c}{P}_k^{d_k}=\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{E}_k^{d_k}{\mathrm{\&phi;}}_{k,i}^{(\&CenterDot;)}{B}_i^{d_i}{u}_i^{d_i}+\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{E}_k^{d_k}{\mathrm{\&phi;}}_{k,i}^{(\&CenterDot;)}{B}_i^{f_i}{u}_i^{f_i}+\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{E}_k^{d_k}{\mathrm{\&phi;}}_{k,i}^{(\&CenterDot;)}{B}_i^{t_i}{u}_i^{t_i}+E_k^{d_k}{\mathrm{\&phi;}}_{k,0}^{(\&CenterDot;)}{\mathrm{\&mu;}}_0+\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{E}_k^{d_k}{\mathrm{\&phi;}}_{k,i}^{(\&CenterDot;)}{w}_i+v_k\\ =E_k^{d_k}(\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&phi;}}_{k,i}^{(\&CenterDot;)}\left(\begin{array}{cccc}{B}_i^{d_i}& B_i^{f_i}& B_i^{t_i}& I\end{array}\right)*{\left(\begin{array}{cccc}{u}_i^{d_i}& u_i^{f_i}& u_i^{t_i}& w_i\end{array}\right)}^T+{\mathrm{\&phi;}}_{k,0}^{(\&CenterDot;)}{\mathrm{\&mu;}}_0)+v_k\\ =E_k^{d_k}*U_k+v_k\end{array}---\left(10\right)$

Can be write as:

$P_k^{d_k}={\mathrm{\&alpha;}}_k+E_k^{d_k}*U_k+v_k---\left(11\right)$

Wherein the vector of one (p × 1), the mass property that representative exports; α _kbe one (p × 1) vector, represent the compensating parameter of each output characteristics; one (p × m) vector, system of representatives matrix number; $U_k=\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&phi;}}_{k,i}^{(\&CenterDot;)}\left(\begin{array}{cccc}{B}_i^{d_i}& B_i^{f_i}& B_i^{t_i}& I\end{array}\right)*{\left(\begin{array}{cccc}{u}_i^{d_i}& u_i^{f_i}& u_i^{t_i}& w_i\end{array}\right)}^T+{\mathrm{\&phi;}}_{k,0}^{(\&CenterDot;)}{\mathrm{\&mu;}}_0$ Be one (m*1) vector, represent the error source of process system; v _kthe vector of one (p × 1), the noise in expression process.

2.3.2 the v in (4) formula is supposed _kbe a mean vector be 0, covariance matrix is the white noise sequence of ∑, then the expression formula of MGWMA controller (this controller is used for adjusting adjustment model) is as follows:

$g_i=\underset{k=1}{\overset{i}{\mathrm{\&Sigma;}}}[q^{{(i-t)}^b}-q^{{(i-t+1)}^b}]P_k,i\&Element;I^{+}---\left(12\right)$

G in above formula _ione (p × 1) vector (g ₀=0), design parameter q is a constant and (0≤q<1), and adjustment parameter b is decided by technician.

2.3.3 according to MGWMA controller, g can be obtained _icovariance matrix be:

$\begin{array}{c}{\mathrm{\&Sigma;}}_{g_i}=\mathrm{var}\{\underset{t=1}{\overset{i}{\mathrm{\&Sigma;}}}\left[[q^{(i-t)b}-q^{{(i-t+1)}^b\left]P_t\right\}}\right]\\ =(q^{0^b}-q^{1^b}){\mathrm{\&Sigma;}}_{P^i}(q^{0^b}-q^{1^b})+(q^{1^b}-q^{2^b}){\mathrm{\&Sigma;}}_{P_i}(q^{1^b}-q^{2^b})\\ +...+(q^{{(i-1)}^b}-q^{i^b}){\mathrm{\&Sigma;}}_{P_i}(q^{{(i-1)}^b}-q^{i^b})\\ =Q_i{\mathrm{\&Sigma;}}_{P_i}\\ =Q_i\left[(1+\frac{\mathrm{\&omega;}}{2-\mathrm{\&omega;}}(1-{(1-\mathrm{\&omega;})}^{2(i-1)}))\mathrm{\&Sigma;}\right]\end{array}---\left(13\right)$

Wherein, $\begin{array}{c}{Q}_i=(q^{0^b}-q^{1^b}){\mathrm{\&Sigma;}}_{P_i}(q^{0^b}-q^{1^b})+(q^{1^b}-q^{2^b}){\mathrm{\&Sigma;}}_{P_i}(q^{1^b}-q^{2^b})\\ +...+(q^{{(i-1)}^b}-q^{i^b}){\mathrm{\&Sigma;}}_{P_i}(q^{{(i-1)}^b}-q^{i^b})\end{array}.$

2.3.4MGWMA controller introduces a deviation on the basis in initial error source, makes control chart (namely statistic exceeds control limit) out of control, and then, alignment error source, until this statistic is less than this control limit.

Described statistic is:

$\begin{array}{c}{T_i}^2={g_i}^{\&prime;}{\left[Q_i{\mathrm{\&Sigma;}}_{P_i}\right]}^{-1}{g}_i\\ ={g_i}^{\&prime;}{\left[Q_i(1+\frac{\mathrm{\&omega;}}{2-\mathrm{\&omega;}}(1-{(1-\mathrm{\&omega;})}^{2(i-1)}))\mathrm{\&Sigma;}\right]}^{-1}{g}_i>h_3\end{array}---\left(14\right)$

Wherein upper control limit h ₃(being greater than 0) is the specific ARL according to reaching ₀and selected.

Embodiment:

The validity of extracting method in error modeling in order to prove, the sheet part as shown in Figure 6 of consideration, wherein ps ₁to ps ₄for infinitesimal.Consider several situation that may produce error source, assuming that process is made up of two procedures, first operation working position 1, second operation work working position 2, the machining feature of first operation is the machining benchmark of second operation work.

(1) first operation error prediction

Arranging noise is 0, then can obtain infinitesimal ps by dynamic SoV forecast model formula ₁prediction error value be:

[-0.0035，-0.0334，0.00126，0.0006，0.0000，0.0003] ^T

(2) second operation work error prediction

Second operation work benchmark ps can be known by first operation ₁the infinitesimal deviation produced is under infinitesimal coordinate system:

[-0.0035 ,-0.0334,0.00126,0.0006,0.0000,0.0003] ^t, benchmark ps ₂and ps ₃deviation be approximately [0,0,0,0,0,0].

Second operation work infinitesimal ps ₄total departure is:

${\mathrm{\&Delta;X}}_p^{\mathrm{wsI}}={\mathrm{\&Delta;X}}_p^{\mathrm{wsI}}+\mathrm{\&Delta;}{X}_p^{\mathrm{wsLI}}={[-0.0032,-\mathrm{0.0109,0.0119,0,0,0},]}^T$

The error arranging fixture is 0, and location attitude of the cutter error is 0, and noise is 0, then can obtain final infinitesimal ps by dynamic SoV prediction ₄error amount under workpiece coordinate system is:

[-0.0032,-0.0109,0.0119,0,0,0,] ^T

Can find out the change of the cutting force caused due to geometry location error and cutting-in, thus cause y direction and z direction and create deviation, the error meeting set error source produces rule, and therefore prediction is rational.

(3) quality Pre-control system

The error analysis of first operation:

A error that () benchmark causes

Benchmark is the error that a upper procedure causes, and first operation does not go up a procedure, then its fiducial error is set to 0.

B workpiece machining error that () fixture geometrical deviation causes

Simulate 30 groups of samples, under the analyzing samples supposition steady state (SS) of first 20 groups, the mismachining tolerance Normal Distribution N (0 that fixture geometrical deviation causes, 0.0075*0.0075), latter 10 groups are used for control sample, by filling in the patch of different size, introducing fixture mean bias delt=0.5, carrying out the simulation of fluctuating.

The error of c workpiece that () cutting force causes

Simulate 30 groups of samples, first 20 groups is under analysis sample supposition steady state (SS), the workpiece error Normal Distribution N (0,0.0086*0.0086) that cutting force causes, introduce mean bias delt=-0.6972, simulation obtains the control sample of latter 10 groups.

D error that () cutter path causes

Cutter path error is set to 0.

E () surveying instrument noise error obeys N (0,0.00025*0.00025)

Comprehensive above-mentioned error calculation goes out total predicated error, then adopts control chart to monitor to it.Estimate its covariance by the predicted data obtained, can know that obtained data are relevant, therefore adopt multivariate control chart to monitor it.Because this step simulations is the roughing of blade, so be great fluctuation process, therefore adopt T-square control chart.Analysis control chart is set up, as shown in Figure 7 by front 20 groups of data.Then set up control control chart according to the rear 10 groups of data introducing the data after fluctuating to monitor, obtain the monitored results of first operation as shown in Figure 8.

Using the input of first operation as second operation work, carry out digital simulation.

The error analysis of second operation work:

A error that () benchmark causes

The error calculating first operation is above incorporated in second operation work.

B workpiece machining error that () fixture geometrical deviation causes

Simulate 30 groups of samples, assuming that under steady state (SS), the mismachining tolerance Normal Distribution N (0 that fixture geometrical deviation causes, 0.0075*0.0075), simulating the analysis sample of first 20 groups, by filling in the patch of different size, introducing deviation delt=0.5, carry out the calculating of error, the control sample that simulation is latter 10 groups.

The error of c workpiece that () cutting force causes

Simulate 30 groups of samples, assuming that under steady state (SS), the workpiece error Normal Distribution N (0 that cutting force causes, 0.0086*0.0086), simulating first 20 groups is analysis sample, then introduce deviation delt=-0.6972, calculate the error of workpiece, simulation obtains the control sample of latter 10 groups.

D error that () cutter path causes

Cutter path error is set to 0.

E () surveying instrument noise error obeys N (0,0.00025*0.00025)

Comprehensive above-mentioned error calculation goes out total error, then control chart is adopted to monitor to it, the same with first operation, estimate its covariance by the predicted data obtained, can know that obtained data are relevant, therefore adopt multivariate control chart to monitor it.Because this step simulations is the roughing of blade, so be great fluctuation process, therefore adopt T-square control chart.Set up analysis control chart according to the data under steady state (SS), as shown in Figure 9, set up control control chart according to the data after introducing fluctuation, as shown in Figure 10.

(4) process adjusting

The fluctuation of the process adjusting method adjustment first operation adopting the present invention to propose, as shown in figure 11, can find out that process has returned and stablize, the adjustment process of second operation work is similar, does not repeat at this for result.

As can be seen from above-mentioned demonstration, the simulation result of error simulated prediction method herein and quality Pre-control method well can both react practical condition.

The above; be only the present invention's preferably embodiment, but protection scope of the present invention is not limited thereto, is anyly familiar with those skilled in the art in the technical scope that the present invention discloses; the change that can expect easily or replacement, all should be encompassed within protection scope of the present invention.Therefore protection scope of the present invention should be as the criterion with the protection domain of claim.

Claims (9)

1. the closed loop quality control emulation mode that simulation and prediction is mutually integrated with actual production, is characterized in that, comprise the following steps:

1) be coupled under the error source that part process produces being transformed into infinitesimal coordinate system, then according to the error after coupling, consider that practical production status is on the impact of error source, sets up forecast model simultaneously;

2) Acquisition Error information, it can be used as the input of forecast model;

3) forecast model is converted to dynamic SOV model, sets up the OSFE equation corresponding with SOV model, obtain independent identically distributed deviation;

4) according to the control chart of independent identically distributed Deviation Design forecast model;

5) Serial regulation model and MGWMA controller is set up;

6) covariance matrix of controller is solved according to MGWMA controller;

7) make MGWMA control chart introduce a deviation on the basis in initial error source, make control chart out of control, then alignment error source, until this statistic is less than this control limit.

2. the closed loop quality control emulation mode that a kind of simulation and prediction according to claim 1 is mutually integrated with practical production status, is characterized in that, described step 4) in, the method for the control chart of design forecast model is:

4.1) mean vector of independent same distribution deviation and the estimated value of covariance matrix is calculated according to SOV model;

4.2) according to the mean vector of independent same distribution deviation and the estimated value compute statistics of covariance matrix;

4.3) limit according to the control of the distributed problem solving multivariate statistical procedure of statistic obedience.

3. the closed loop quality control emulation mode that a kind of simulation and prediction according to claim 2 is mutually integrated with practical production status, is characterized in that, described step 4.1) according to following formulae discovery:

${\stackrel{\&OverBar;}{\hat{e}}}_k^{d_k}={\left[\begin{array}{cccc}{\stackrel{\&OverBar;}{\hat{e}}}_{k1}^{d_k}& {\stackrel{\&OverBar;}{\hat{e}}}_{k2}^{d_k}& ...& {\stackrel{\&OverBar;}{\hat{e}}}_{\mathrm{ks}}^{d_k}\end{array}\right]}^T$

$s_k=\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}({\hat{e}}_{\mathrm{ki}}-{\stackrel{\stackrel{\&OverBar;}{^}}{e}}_k){({\hat{e}}_{\mathrm{ki}}-{\stackrel{\stackrel{\&OverBar;}{^}}{e}}_k)}^T$

Wherein, what represent is that the deviation that produces of error source is at the mean vector of the independent same distribution deviation of operation k.

4. the closed loop quality control emulation mode that a kind of simulation and prediction according to claim 3 is mutually integrated with practical production status, is characterized in that, described step 4.2) concrete grammar be:

Assuming that the vector of a new independent same distribution deviation, and and S _kindependently vectorial;

The then T of kth procedure _k ²the expression formula of statistic is:

${T_k}^2={({\hat{e}}_{\mathrm{ku}}-{\stackrel{\stackrel{\&OverBar;}{^}}{e}}_k)}^T{S}_k^{-1}({\hat{e}}_{\mathrm{ku}}-{\stackrel{\stackrel{\&OverBar;}{^}}{e}}_k)~\frac{s(m+1)(m-1)}{m(m-s)}{F}_{s,m-s}.$

5. the closed loop quality control emulation mode that a kind of simulation and prediction according to claim 4 is mutually integrated with practical production status, is characterized in that, step 4.3) concrete grammar be:

Assuming that the probability of error of first kind is α, then the control of multivariate statistical procedure is limited to:

${\mathrm{UCL}}_k=\frac{s(m+1)(m-1)}{m(m-s)}{F}_{s,m-s}\left(\mathrm{\&alpha;}\right).$

6. the closed loop quality control emulation mode that a kind of simulation and prediction according to claim 1 is mutually integrated with practical production status, is characterized in that, described step 1) in forecast model be:

P＝A _JFMU _JFMJ _FM+B _WFMU _WFMW _FM，

Wherein, A _jFMrepresenting the correction factor of geometrical deviation source dynamic change in time, is static coefficient U _jFMjacobian matrix, B _wFMrepresent the correction factor of physical error source dynamic change in time, U _jFMrepresent the static Jacobian matrix coefficient in geometric error source, U _wFMrepresent the static Jacobian matrix coefficient in physical error source, J _fMrepresent geometrical deviation, W _fMrepresent physical deflections.

7. the closed loop quality control emulation mode that a kind of simulation and prediction according to claim 1 is mutually integrated with practical production status, is characterized in that, step 3) in dynamic SOV model be:

$\left\{\begin{array}{cc}{\mathrm{\&mu;}}_k^{d_k,t_k,l_k,j_k}=A_{k-1}{\mathrm{\&mu;}}_{k-1}^{d_{k-1},t_{k-1},l_{k-1},j_{k-1}}{B}_k^{d_k}{u}_k^{d_k}+B_k^{f_k}{u}_k^{f_k}+B_k^{t_k}{u}_k^{t_k}+B_k^{j_k}{u}_k^{j_k}+w_k,w_k~N[0,W_k]& \\ P_k^{d_k}=P=C_E{E}_k^{d_k}{\mathrm{\&mu;}}_k+v_k,& v_k~N[0,V_k]\\ {\mathrm{\&mu;}}_0|D_0~N[m_0,C_0]& \end{array}\right.$

Wherein, $\left\{\begin{array}{c}{A}_{k-1}=U_{\mathrm{JFM}}{\left(C_{\mathrm{ws}}^O{U}_{\mathrm{wsi}}\right)}^{+}D\\ B_k^{f_k}={-U}_{\mathrm{JFM}}{\left(C_{\mathrm{ws}}^O{U}_{\mathrm{wsi}}\right)}^{+}(N\\ B_k^{t_k}=C_T^O{U}_{\mathrm{JFM}}\\ B_k^{j_k}=C_T^{\mathrm{Ps}}R_T^O\end{array}\right.$

what represent is the coordinate conversion of measuring system, and CE represents the correction factor that measuring system changes with practical production status.

8. the closed loop quality control emulation mode that a kind of simulation and prediction according to claim 1 is mutually integrated with practical production status, it is characterized in that, described linear model has m error originated from input source, p output characteristics value, according to SOV model, the adjustment model of multiple-input and multiple-output process is described below:

$P_k^{d_k}={\mathrm{\&alpha;}}_k+E_k^{d_k}*U_k+v_k$

Wherein, the vector of one (p × 1), the mass property that representative exports; α is one (p × 1) vector, represents the compensating parameter of each output characteristics; one (p × m) vector, system of representatives matrix number; $U_k=\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&phi;}}_{k,i}^{(\&CenterDot;)}\left(\begin{array}{cccc}{B}_i^{d_i}& B_i^{f_i}& B_i^{t_i}& I\end{array}\right)*{\left(\begin{array}{cccc}{u}_i^{d_i}& u_i^{f_i}& u_i^{t_i}& w_i\end{array}\right)}^T+{\mathrm{\&phi;}}_{k,0}^{(\&CenterDot;)}{\mathrm{\&mu;}}_0$ Be one (m*1) vector, represent the error source of process system; v _kthe vector of one (p × 1), the noise in expression process.

9. the closed loop quality control emulation mode that a kind of simulation and prediction according to claim 1 is mutually integrated with practical production status, is characterized in that, assuming that the v in SOV model _kbe a mean vector be 0, covariance matrix is the white noise sequence of ∑, then described step 5) in MGWMA controller be described below:

$\begin{array}{cc}{g}_i=\underset{k=1}{\overset{i}{\mathrm{\&Sigma;}}}[q^{{(i-t)}^b}-q^{{(i-t+1)}^b}]P_k& i\&Element;I^{+}\end{array}$

G in above formula _ione (p × 1) vector (g ₀=0), design parameter q is a constant and (0≤q<1), and adjustment parameter b is artificial set-point.