CN101520651A - Analysis method for reliability of numerical control equipment based on hidden Markov chain - Google Patents

Abstract

The invention provides an analysis method for reliability of numerical control equipment based on a hidden Markov chain. The method particularly comprises the following steps: 1, monitoring dynamic performance signals of the numerical control equipment, and abstracting the performance characteristic parameter values showing the change of the reliability of the numerical control equipment; 2, constructing a predictive model of the performance characteristic parameter values; 3, using the predictive model to predict the performance characteristic parameter values within the time needed for vectorization, and adopting the source coding method to vectorize the predictive values of the performance characteristic parameter values; 4, adopting the discrete hidden Markov chain model to identify the state transition probability of the numerical control equipment; and 5, utilizing the Chapman-Kolmogorov differential equation to establish a relational expression of the operational state and the state transition probability so as to deduce the probability of the numerical control equipment in different operating states, namely obtaining the reliability of the numerical control equipment. The method can accurately analyze, evaluate and predict the change of the reliability of the numerical control equipment before the numerical control equipment goes wrong, thereby avoiding the fault of the numerical control equipment and improving the operational reliability of the numerical control equipment.

Description

A kind of analysis method for reliability of numerical control equipment based on hidden Markov chain

Technical field

The present invention belongs to the numerical control equipment technical field, specifically is a kind of reliability of numerical control equipment variation analysis method based on discrete Hidden Markov chain model.

Background technology

Reliability is a kind of inherent attribute of numerical control equipment, is meant that numerical control equipment carries out it in assignment of mission and condition thereof in the next stipulated time and require the ability of the various key functions of task.It has reflected the mass property in the numerical control equipment use, time factor and environmental factor when having considered the numerical control equipment operation.

At present, the general thinking of reliability of numerical control equipment analysis is by gathering and analyzing numerically controlled equipment composition parts and the design of machine system and the historical failure and the lifetime data of manufacturing information, field failure data and like product, start with from the Analysis on Fault Diagnosis of numerical control equipment, according to numerical control equipment work and functional characteristics, trouble-shooting pattern and reason, assessment reliability of numerical control equipment level, and the reliability innovative approach is proposed.Wherein, mostly the approach of field failure data acquisition is to follow the tracks of for a long time tens numerical control equipments, and fault data derives from the field notes and the maintenance report of numerical control equipment.This shows, the same model numerical control equipment that existing analysis method for reliability of numerical control equipment needs and works long hours, quantity is abundant, and a large amount of faults and lifetime data will be arranged; By the analysis of failure data, and based on certain statistics judgment criterion, the failure probability distributions of identification equipment, and then adopt mathematical statistics method to infer reliability of numerical control equipment, find out the weak link in the numerical control equipment use.

Yet this subject matter that faces based on the reliability of numerical control equipment research of fault data has: on the one hand, it needs a large amount of fault and lifetime data, to find the life-span distribution curve of numerical control equipment.May obtain the data of a large amount of numerical control equipment fault and lifetime data, particularly machine system in the actual engineering hardly.In addition, numerical control equipment breaks down through long-time running, and this moment, equipment may be in end of lifetime, caused the reliability result of gained to lack actual directive significance.On the other hand, it has only considered the normal and fault two states of numerical control equipment, does not consider performance state variation feature in the numerical control equipment use.For example, numerical control equipment breaks down and goes to toward having experienced a series of deterioration state, is embodied in the output performance parameter and progressively departs from its setting, until final generation disabler.This causes the user can't understand the real-time change conditions of reliability of numerical control equipment.

Therefore, press for a kind of new method that is used for the reliability of numerical control equipment field of design and development, it does not need to wait for that numerical control equipment breaks down through long-play, just can analyze, assess and predict the reliability of numerical control equipment change conditions exactly, thereby guides user is taked rational preventive measure in advance, prevents that fault from taking place.

Summary of the invention

The objective of the invention is to be at the deficiencies in the prior art, a kind of reliability of numerical control equipment variation analysis method based on discrete Hidden Markov chain model is provided, before breaking down, numerical control equipment analyzes, assess and predict the reliability of numerical control equipment change conditions exactly, find out weak link in the numerical control equipment use etc., avoid the generation of numerical control equipment fault, improve the operational reliability of numerical control equipment.

A kind of analysis method for reliability of numerical control equipment based on hidden Markov chain may further comprise the steps:

Dynamic property signal when (1) the monitoring numerical control equipment moves extracts the performance characteristic parameter value that more than one reflection reliability of numerical control equipment changes at more than one time point;

(2), make up the forecast model of performance characteristic parameter value according to the performance characteristic parameter value that extracts;

(3) use forecast model prediction T ₀Each performance characteristic parameter value of individual time point carries out vector quantization to these predicted values, obtains the output observation sequence of performance characteristic parameter;

(4), adopt discrete Hidden Markov chain model to determine the running status transition probability matrix of numerical control equipment according to the output observation sequence;

(5) utilize Qie Puman-Andrei Kolmogorov differential equation to set up the relational expression of running status and status change probability matrix, infer that numerical control equipment is in the probability under the different running statuses.

Described step (2) adopts non-stationary autoregression integration running mean algorithm to set up the forecast model of j performance characteristic ginseng value $y_{\mathrm{tj}}=\underset{z=1}{\overset{p}{\mathrm{\Σ}}}{\mathrm{\φ}}_z{y}_{(t-z)j}+{\mathrm{\ϵ}}_{\mathrm{tj}}-\underset{z=1}{\overset{q}{\mathrm{\Σ}}}{\mathrm{\φ}}_z{\mathrm{\ϵ}}_{(t-z)j},$ Wherein, t represents the predicted time point, and p is the autoregression order, and q is the running mean order, φ _zBe running mean coefficient, measuring error ${\mathrm{\ϵ}}_{\mathrm{tj}}~N(0,{\mathrm{\σ}}_{\mathrm{\ϵ}}^2).$

Described step (3) is specially: the output observation sequence of performance characteristic parameter is expressed as $O^j=\{o_1^j,\&CenterDot;\&CenterDot;\&CenterDot;,o_{T_0}^j\},$ $o_{\tilde{t}}^j\&Element;(v_1^j,,v_N^j),$ $\tilde{t}=1,\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{<figure-callout\; id="0"\; label="T"\; filenames="A200910060949E00181.png"\; state="\{\{state\}\}">T</figure-callout>}}_0,$ J=1 ..., m, m are the performance characteristic number of parameters, t represents the predicted time point, wherein,

$v_1^j=1,v_N^j=N,$ N is the zone of dispersion number;

$v_{i_0}^j=\mathrm{Index}\left(y_{t_{1}j}\right)=\left\{\begin{array}{c}1,y_{t_{1}j}\&le;\mathrm{partition}(1,j)\\ i_0,\mathrm{partition}(i_0-1,j)\&le;y_{t_{1}j}\&le;\mathrm{partition}(i_0,j);\\ N,\mathrm{partition}(N-1,j)\&le;y_{t_{1}j}\end{array}\right.$

i ₀=2 ..., N-1, t ₁∈ (1 ..., T ₀),

Represent that j performance characteristic parameter is at t ₁The predicted value of time point, $\mathrm{partition}({\mathrm{<figure-callout\; id="1"\; label="i"\; filenames="A200910060949E00181.png,A200910060949E00182.png"\; state="\{\{state\}\}">i</figure-callout>}}_1,j)=h_{i_{1}j},$ i ₁=1 ..., N-1,

Represent that j performance characteristic parameter is in i ₁The value of defining of zone of dispersion.

Described step (4) is specially: make numerical control equipment comprise 0～k+1 kind running status, the output observation sequence and the predefined original state transition probability matrix A that obtain according to step (3) ₀, initial observation value probability matrix B ₀With initial probability distribution vector π ₀, utilize Bao Mu-Wei Erqi (Baum-Welch) algorithm iteration calculate numerical control equipment constantly status change model λ of t=(π), π represents the probability distribution vector for A, B,

The status change probability matrix $A={\left[a_{\mathrm{cd}}\right]}_{(k+1)\&times;(k+1)}=\left[\begin{array}{ccccc}{a}_{00}& a_{01}& \&CenterDot;\&CenterDot;\&CenterDot;& a_{0k}& a_{0(k+1)}\\ 0& a_{11}& \&CenterDot;\&CenterDot;\&CenterDot;& a_{1k}& a_{1(k+1)}\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ 0& 0& \&CenterDot;\&CenterDot;\&CenterDot;& a_{\mathrm{kk}}& a_{k(k+1)}\\ 0& 0& \&CenterDot;\&CenterDot;\&CenterDot;& 0& a_{(k+1)(k+1)}\end{array}\right],$ a _Cd{ c} is in state { probability of d}, 0≤a constantly at t+1 for numerical control equipment is in state constantly at t _Cd≤ 1, and $\underset{d=0}{\overset{k+1}{\mathrm{\&Sigma;}}}{a}_{\mathrm{cd}}=1;$

The observed reading probability matrix $B={\left[b_{\mathrm{cj}}\right]}_{(k+1)\&times;m}=\left[\begin{array}{cccc}{b}_{01}& b_{02}& \&CenterDot;\&CenterDot;\&CenterDot;& b_{0m}\\ b_{11}& {\mathrm{<figure-callout\; id="12"\; label="b"\; filenames="A200910060949E00192.png"\; state="\{\{state\}\}">b</figure-callout>}}_{12}& \&CenterDot;\&CenterDot;\&CenterDot;& b_{1m}\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ b_{k1}& b_{k2}& \&CenterDot;\&CenterDot;\&CenterDot;& b_{\mathrm{km}}\\ b_{(k+1)1}& b_{(k+1)2}& \&CenterDot;\&CenterDot;\&CenterDot;& b_{(k+1)m}\end{array}\right],$ b _CjFor numerical control equipment is in state constantly at t { observed reading appears during c} $\left\{v_c^j\right\},c=0,\&CenterDot;\&CenterDot;\&CenterDot;,k+1$ Probability, 0≤b _Cj≤ 1, and $\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{b}_{\mathrm{cj}}=1.$

The beneficial effect that the present invention has is embodied in:

Numerical control equipment is similar to human body, and its function is decline gradually along with time lengthening, and the function that is embodied in numerical control equipment and performance be deterioration progressively, causes the no plan downtime rate increase of numerical control equipment generation, has reduced the numerical control equipment operational reliability.Patent of the present invention is monitored, is analyzed and predict by performance degradation and reliability change to numerical control equipment, can help enterprise to find the potential failure mode of numerical control equipment, formulate rational maintenance schedule and sudden inefficacy takes place to avoid numerical control equipment, thereby realize the reliability growth of numerical control equipment, this all has great importance for improving numerical control equipment utilization factor, minimizing numerical control equipment maintenance cost, prolongation numerical control equipment serviceable life and structure numerical control equipment health status monitoring system etc.

Description of drawings

Fig. 1 general principles figure of the present invention.

Fig. 2 vector quantization figure of the present invention.

Fig. 3 status change topology diagram of the present invention.

Fig. 4 reliability change of the present invention result of calculation figure.

Embodiment

Below in conjunction with drawings and Examples technical solutions according to the invention are described further.

With reference to Fig. 1, a kind of reliability of numerical control equipment change computing method based on discrete Hidden Markov chain model mainly comprise the performance characteristic parameter monitoring and the modules such as extraction, performance degradation model, the quantification of performance characteristic parameter vector, status change model and Calculation of Reliability of numerical control equipment.

1, performance characteristic parameter monitoring and extraction.

The performance characteristic parameter monitoring is by being installed on the numerical control equipment or on every side testing tools such as the sensor various dynamic property signals (as vibration signal, motor message, hydraulic pressure signal etc.) when coming the testing equipment operation, extracting characteristic parameter through signal Processing and analysis with extracting.The performance characteristic parameter of being extracted is relevant with the type and the disposal route thereof of signal, but it should follow following principle generally: (1) performance characteristic parameter should have clear physical meaning, (2) the performance characteristic parameter is easy to measure and is relatively stable, the key parameter that (3) performance characteristic parameter changes for the reflection reliability of numerical control equipment etc.For example,, and adopt time domain method to analyze the time history of vibration signal if gather the vibration signal of numerical control equipment with acceleration transducer, then the performance characteristic parameter that can extract have that peak-peak, root mean square refer to, peak factor and the nargin factor etc.

According to the testing tool of prior layout, put and measure respectively the l in the numerical control equipment use preset time, obtains the dynamic property signal.The signal of gathering is handled and is analyzed, extract the performance characteristic parameter that can reflect timely and accurately that the numerical control equipment performance state changes, promptly D={ (T, Y) | (t _i, y _Ij); I=1,2 ..., l; J=1,2 ..., m}, wherein t _iBe i Measuring Time point, y _IjFor at t _iJ the performance characteristic parameter that the moment measures.

2, set up the performance degradation model of numerical control equipment.

The discrete-time series { (t of check feature characteristic parameter _i, y _Ij); I=1,2 ..., l; J=1,2 ..., the stationarity of m}.If non-stationary then needs to adopt { the y of d jump point-score with non-stationary _IjTime series is converted into { w stably _IjTime series, promptly $w_{\mathrm{ij}}={\&dtri;}^d{y}_{\mathrm{ij}}.$

Adopt the nARIMA algorithm to { w _IjTime series carries out modeling.Because j performance characteristic parameter w _IjEach value w in the not only preceding p of value step with it _{(i-1) j}, w _{(i-2) j}..., w _{(i-p) j}Relevant, and the random perturbation value ε in q preceding with it step _{(i-1) j}, ε _{(i-2) j}..., ε _{(i-q) j}, then obtain p rank autoregression q rank moving average model:

$w_{\mathrm{ij}}=\underset{z=1}{\overset{p}{\mathrm{\&Sigma;}}}{\mathrm{\&theta;}}_z{w}_{(i-z)j}-\underset{z=1}{\overset{q}{\mathrm{\&Sigma;}}}{\mathrm{\&phi;}}_z{\mathrm{\&epsiv;}}_{(i-z)j}+{\mathrm{\&epsiv;}}_{\mathrm{ij}}---\left(1\right)$

Wherein, θ _zBe autoregressive coefficient, φ _zBe the running mean coefficient, ${\mathrm{\&epsiv;}}_{\mathrm{ij}}~N(0,{\mathrm{\&sigma;}}_{\mathrm{\&epsiv;}}^2)$ Be measuring error,

Be variance.Introduce backward shift operator B, then formula (1) is converted to formula (2):

$\left\{\begin{array}{c}\mathrm{\&theta;}\left(B\right)w_{\mathrm{ij}}=\mathrm{\&phi;}\left(B\right){\mathrm{\&epsiv;}}_{\mathrm{ij}}\\ \mathrm{\&theta;}\left(B\right)=1-\underset{z=1}{\overset{p}{\mathrm{\&Sigma;}}}{\mathrm{\&theta;}}_z{B}^z,\mathrm{\&phi;}\left(B\right)=1-\underset{z=1}{\overset{q}{\mathrm{\&Sigma;}}}{\mathrm{\&phi;}}_z{B}^z\end{array}\right.---\left(2\right)$

Wherein, the order p of model and q and unknown parameter θ ₁..., θ _p, φ ₁..., φ _qAnd σ _εDirectly have influence on model accuracy.Adopt minimal information criterion (AIC criterion) to determine model order, adopt nonlinear least square method to come the estimation model unknown parameter.In case determine model order and model parameter, then obtained non-stationary nARIMA (p, d, q) performance degradation model:

$y_{\mathrm{tj}}=\underset{z=1}{\overset{p}{\mathrm{\&Sigma;}}}{\mathrm{\&phi;}}_z{y}_{(t-z)j}+{\mathrm{\&epsiv;}}_{\mathrm{tj}}-\underset{z=1}{\overset{q}{\mathrm{\&Sigma;}}}{\mathrm{\&phi;}}_z{\mathrm{\&epsiv;}}_{(t-z)j}---\left(3\right)$

Wherein, φ _zBe running mean coefficient, measuring error ${\mathrm{\&epsiv;}}_{\mathrm{tj}}~N(0,{\mathrm{\&sigma;}}_{\mathrm{\&epsiv;}}^2)$ P is the autoregression order, and the representation model predicted value is relevant with preceding p rank performance parameter value, and q is the running mean order, and expression is carried out the correction on q rank to the model predication value error.

3, the performance characteristic parameter vector quantizes

In case set up the performance degradation model, then can predict the performance characteristic parameter situation of change of numerical control equipment in following a period of time.For example, { y _Tj, j=1 ..., m} represents one group of t moment performance characteristic parameter that do as one likes energy degradation model forecasting institute gets.These performance characteristic parameters have reflected the implicit state variation feature of numerical control equipment from the side.

Yet because numerical control equipment itself is very complicated, the performance characteristic parameter that is observed in its use can not be corresponding simply one by one with state.Therefore, the present invention uses discrete Hidden Markov chain model, performance state variation characteristic and reliability change conditions thereof when setting up numerical control equipment status change model and coming analyzing numerically controlled equipment operation.

Because the input quantity of status change model is necessary for limited discrete value, and the performance characteristic parameter that is obtained is a continuous real number value, adopts Lloyd (Laue moral) algorithm that the performance characteristic parameter is carried out vector quantization for this reason and handle to form discrete code set.

With reference to Fig. 2, the vector quantization process of performance characteristic parameter is: according to subregion vector (partition) and codebook vectors (codebook), vector quantization is divided into N zone in proper order with the characteristic parameter span, discrete observation value of each zone mapping.The required predicted time of the vector quantization T that counts ₀Generally get three to six unit interval points, at this T ₀In the individual predicted time point section first is to T ₀Time point, j characteristic ginseng value is at each regional index value

Be defined as:

$v_1^j=1,v_N^j=N,$

$v_{i_0}^j=\mathrm{Index}\left(y_{t_{1}j}\right)=\left\{\begin{array}{c}1,y_{t_{1}j}\&le;\mathrm{partition}(1,j)\\ i_0,\mathrm{partition}(i_0-1,j)\&le;y_{t_{1}j}\&le;\mathrm{partition}(i_0,j)\\ N,\mathrm{partition}(N-1,j)\&le;y_{t_{1}j}\end{array}\right.$

i ₀＝2，…，N-1，t ₁∈(1，…，T ₀)；

Wherein, subregion vector (partition) is defined by the performance failure canonical matrix H of equipment, that is:

$H={\left[h_{i_{1}j}\right]}_{(N-1)\&times;m}=\left[\begin{array}{ccc}{h}_{11}& \&CenterDot;\&CenterDot;\&CenterDot;& h_{1m}\\ h_{21}& \&CenterDot;\&CenterDot;\&CenterDot;& h_{2m}\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ h_{(N-1)1}& \&CenterDot;\&CenterDot;\&CenterDot;& h_{(N-1)m}\end{array}\right],$ $\mathrm{partition}({\mathrm{<figure-callout\; id="1"\; label="i"\; filenames="A200910060949E00181.png,A200910060949E00182.png"\; state="\{\{state\}\}">i</figure-callout>}}_1,j)=h_{i_{1}j},$

Being used to of expression original definition weighed j characteristic parameter and is in i ₁The value of defining in zone, i ₁=1 ..., N-1

The numerical control equipment that forecasting institute is got is at 1～T ₀Performance characteristic parameter set { y constantly _TjCarry out vector quantization, obtain the output observation sequence of performance characteristic parameter:

$O^j=\{o_1^j,\&CenterDot;\&CenterDot;\&CenterDot;,o_{T_0}^j\},$ $o_{\tilde{t}}^j\&Element;(v_1^j,,v_N^j),$ $\tilde{t}=1,\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{<figure-callout\; id="0"\; label="T"\; filenames="A200910060949E00181.png"\; state="\{\{state\}\}">T</figure-callout>}}_0,j=1,\&CenterDot;\&CenterDot;\&CenterDot;,m$

Wherein, the required predicted time point T of vector quantization ₀Generally get the 3rd to the 6th unit interval point.

4, status change model

With reference to Fig. 3, set up the status change model of equipment.Running status transition probability is meant that the equipment running status is constant or turns to the probability of poorer running status.

In the status change topology diagram, state space S={0,1,2 ..., k, k+1} are numerical control equipment performance running status collection.Running status comprises normal condition, deterioration state and failure state.Circle is represented latent state among Fig. 3,0} represents normal condition, and 1,2 ..., k} represents k deterioration state, { k+1} represents failure state; Directed arc is called the transition arc, the transfer of expression state.Under the situation of not considering to keep in repair, numerical control equipment performance degradation process is irreversible; Each state can only be to the state transitions of its higher numbering in right side, and each state also can shift to self simultaneously.Therefore, make numerical control equipment be expressed as at t status change probability matrix A constantly:

$A={\left[a_{\mathrm{cd}}\right]}_{(k+1)\&times;(k+1)}=\left[\begin{array}{ccccc}{a}_{00}& a_{01}& \&CenterDot;\&CenterDot;\&CenterDot;& a_{0k}& a_{0(k+1)}\\ 0& a_{11}& \&CenterDot;\&CenterDot;\&CenterDot;& a_{1k}& a_{1(k+1)}\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ 0& 0& \&CenterDot;\&CenterDot;\&CenterDot;& a_{\mathrm{kk}}& a_{k(k+1)}\\ 0& 0& \&CenterDot;\&CenterDot;\&CenterDot;& 0& a_{(k+1)(k+1)}\end{array}\right]$

Wherein, a _Cd{ c} is in state { probability of d}, 0≤a constantly at t+1 for numerical control equipment is in state constantly at t _Cd≤ 1, and $\underset{d=0}{\overset{k+1}{\mathrm{\&Sigma;}}}{a}_{\mathrm{cd}}=1.$

The state do as one likes of numerical control equipment can the characteristic parameter observed reading be come perception.According to numerical control equipment at T0 output observation sequence constantly $O^j=\{o_1^j,\&CenterDot;\&CenterDot;\&CenterDot;,o_{T_0}^j\},$ Make numerical control equipment be expressed as at t observed reading probability matrix B constantly:

$B={\left[b_{\mathrm{cj}}\right]}_{(k+1)\&times;m}=\left[\begin{array}{cccc}{b}_{01}& b_{02}& \&CenterDot;\&CenterDot;\&CenterDot;& b_{0m}\\ b_{11}& {\mathrm{<figure-callout\; id="12"\; label="b"\; filenames="A200910060949E00192.png"\; state="\{\{state\}\}">b</figure-callout>}}_{12}& \&CenterDot;\&CenterDot;\&CenterDot;& b_{1m}\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ b_{k1}& b_{k2}& \&CenterDot;\&CenterDot;\&CenterDot;& b_{\mathrm{km}}\\ b_{(k+1)1}& b_{(k+1)2}& \&CenterDot;\&CenterDot;\&CenterDot;& b_{(k+1)m}\end{array}\right]$

Wherein, b _CjFor numerical control equipment is in state constantly at t { observed reading appears during c} $\left\{v_c^j\right\},c=0,\&CenterDot;\&CenterDot;\&CenterDot;,k+1$ Probability, 0≤b _Cj≤ 1, and $\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{b}_{\mathrm{cj}}=1.$

Under the starting condition, numerical control equipment is in normal operating condition, and then the probability distribution of original state vector is π ₀=(1,0 ..., 0), status change probability matrix A ₀With observed reading probability matrix B ₀Take random device to choose, thereby obtain the status change model λ of numerical control equipment under the starting condition ₀=(A ₀, B ₀, π.)。

After running status collection, status change probability matrix A and observed reading probability matrix B definition, adopt discrete hidden Markov chain model solution status change probability matrix A and observed reading probability matrix B, be specially: with the performance characteristic parameter at 1～T ₀Output observation sequence constantly $O^j=\{o_1^j,\&CenterDot;\&CenterDot;\&CenterDot;,o_{T_0}^j\}$ Be input to original state transition model λ ₀=(A ₀, B ₀, π.) in, utilize Baum-Welch (Bao Mu-Wei Erqi) algorithm carries out iterative computation, makes model parameter progressively tend to more reasonably than the figure of merit, thus obtain numerical control equipment constantly status change model λ of t=(A, B, π).In order to verify the rationality of this status change model, adopt Forward-Backward (front and back to) probability of algorithm computation observation sequence O under given λ, i.e. P (O| λ).If P (O| λ) surpasses expectation value 0.8, then think status change model λ=(A, B are feasible π).

5, Calculation of Reliability

In case determine numerical control equipment t constantly status change model λ=(A, B π), then calculate the reliability change conditions of numerical control equipment.Detailed process is: make P _c(t)=P (q _t=c) the expression numerical control equipment is in the probability of c state constantly at t.According to Qie Puman-Andrei Kolmogorov differential equation, have:

P′(t)＝P(t)·A (4)

Wherein, P (t)=(P ₀(t), P ₁(t) ..., P _k(t), P _(k+1)(t)) be state vector, P ' is the single order differential state vector of P (t) (t), and A is a t status change matrix constantly.Formula (4) is carried out Laplace (Laplce) conversion:

$s\&CenterDot;\left[\begin{array}{c}{P}_0\left(s\right)\\ P_1\left(s\right)\\ \text{\&CenterDot;}\\ \&CenterDot;\\ \&CenterDot;\\ P_k\left(s\right)\\ P_{(k+1)}\left(s\right)\end{array}\right]-\left[\begin{array}{c}{P}_0\left(0\right)\\ P_1\left(0\right)\\ \text{\&CenterDot;}\\ \&CenterDot;\\ \&CenterDot;\\ P_k\left(0\right)\\ P_{(k+1)}\left(0\right)\end{array}\right]\left[\begin{array}{ccccc}{a}_{00}& 0& \&CenterDot;\&CenterDot;\&CenterDot;& 0& 0\\ a_{01}& a_{11}& \&CenterDot;\&CenterDot;\&CenterDot;& 0& 0\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ a_{0k}& a_{1k}& \&CenterDot;\&CenterDot;\&CenterDot;& a_{\mathrm{kk}}& 0\\ a_{0(k+1)}& a_{1(k+1)}& \&CenterDot;\&CenterDot;\&CenterDot;& a_{k(k+1)}& a_{(k+1)(k+1)}\end{array}\right]\left[\begin{array}{c}{P}_0\left(s\right)\\ P_1\left(s\right)\\ \text{\&CenterDot;}\\ \&CenterDot;\\ \&CenterDot;\\ P_k\left(s\right)\\ P_{(k+1)}\left(s\right)\end{array}\right]---\left(5\right)$

If numerical control equipment is in normal condition usually in starting condition, then have:

P(0)＝(P ₀(0)，P ₁(0)，…，P _k(0)，P _(k+1)(0))＝(1，0，…，0，0)

Simultaneously, in t moment status change matrix A substitution formula (5), then calculate:

P(s)＝(P ₀(s)，P ₁(s)，…，P _k(s)，P _(k+1)(s))

Then, P (s) is carried out the Laplace inverse transformation, then obtain numerical control equipment is in different conditions constantly at t probability P (t)=(P ₀(t), P ₁(t) ..., P _k(t), P _(k+1)(t)), calculate numerical control equipment t reliability R (t)=1-P constantly _(k+1)Reliability index such as (t).And then, the reliability change conditions in the analyzing numerically controlled equipment use, the weak link of discovery numerical control equipment.

Embodiment:

This embodiment has provided the specific implementation process of the present invention in engineering practice, simultaneous verification should the invention validity.

With reference to Fig. 4, be research object with Ou Tai OTM-650 CNC milling machine, by the circular motion precision deterioration that simulating assembly causes because of backlass, infer reliability of numerical control equipment change rule.Test unit comprises Ou Tai OTM-650 CNC milling machine, Central China " century star " digital control system HMC-21M, Heidenhain KGM182 plane grating measurement mechanism and ordinary PC etc.Wherein, plane grating KGM182 is a kind of high-accuracy testing tool that is used for non-contact measurement equipment two dimensional motion contour accuracy.Making worktable by NC programming is that 50mm and speed of feed are that 2000mm/min does counterclockwise/circular motion clockwise with the radius in XOY plane, adopts the circular motion track of plane grating on-machine measurement equipment.Allow numerical control equipment move 480 hours, every interval sampling in about 20 hours 1 time obtains 24 groups of circular traces.According to the Circular Test of NC Machine Tools test stone that ISO230-4:1996 provides, 24 groups of circular traces of check and analysis to measure gained calculate corresponding circular motion accuracy characteristic amounts such as circle hysteresis, circle deviation and radius deviation, and are as shown in table 1.

Table 1 circular motion accuracy characteristic amount

Because every kind of characteristic quantity only has 24 measured values, in order to improve match and precision of prediction, adopt interpolation 96 data of sampling generation again according to measured value, 120 circular motion accuracy characteristic amounts are arranged altogether.For every kind of characteristic quantity, choose preceding 100 data as the model training collection, back 20 data are as the model checking collection.Find circular motion accuracy characteristic amount right and wrong stably through check, adopt first difference method that they are converted into time series { w stably for this reason _Ij.Lagging behind with circle is example, determines the order of performance degradation model according to AIC criterion, when (p, when q) being (5,2), the information function of model reaches minimum value 1.7939.And then, adopt nonlinear least square method to determine the model parameter of nARIMA (5,1,2), that is:

$\left\{\begin{array}{c}\mathrm{\&theta;}\left(B\right)=1-{0.6847\mathrm{<figure-callout\; id="1"\; label="B"\; filenames="A200910060949E00181.png,A200910060949E00182.png"\; state="\{\{state\}\}">B</figure-callout>}}^1+{0.8273B}^2-{0.7267B}^3+{0.0877B}^4+{0.1747B}^2\\ \mathrm{\&phi;}\left(B\right)=1-{0.0498B}^1-{0.9502B}^2\end{array}\right.$

Thereby foundation circle hysteresis degradation model is:

y _ij＝0.6847y _(i-1)j-0.8273y _(i-2)j+0.7267y _(i-3)j-0.0877y _(i-4)j-0.1747y _(i-5)j+…

+ε _ij-0.0498ε _(i-1)j-0.9502ε _(i-2)j；i＝5，…，100，j＝4

Utilize back 20 data to verify this model, the FPE of model (final prediction deviation) is 1.5459e-5.Reach a conclusion thus: this model satisfies requirement of actual application.

The process of hysteresis degradation model is justified in similar foundation, can set up the performance degradation model of following radius deviation and circle deviation respectively:

(1) radius deviation (F _Max) degradation model is:

y _ij＝y _(i-1)j+0.9165y _(i-2)j-1.43y _(i-3)j+0.5338y _(i-4)j+0.0033y _(i-5)j+…

+ε _ij-0.0048ε _(i-1)j-0.9952ε _(i-2)j；i＝5，…，100，j＝1

(2) radius deviation (F _Min) degradation model is:

y _ij＝1.267y _(i-1)j-0.1327y _(i-2)j-0.723y _(i-3)j+0.4271y _(i-4)j-0.1187y _(i-5)j+0.2803y _(i-6)j+…

+0.2942y _(i-7)j-0.1943y _(i-8)j+0.1601y _(i-9)j-0.1513y _(i-10)j+0.4735y _(i-11)j-0.166y _(i-12)j+…

-0.4402y _(i-13)j+ε _ij-0.7509ε _(i-1)j+0.2458ε _(i-2)j；i＝13，…，100，j＝2

(3) justifying deviation (G) degradation model is:

y _ij＝0.5102y _(i-1)j+0.5102y _(i-2)j+ε _ij-0.0043ε _(i-1)j+0.207ε _(i-2)j；i＝2，…，100，j＝3

According to above-mentioned performance degradation model, the processing technology of combining with digital control equipment requires and uses rules, with the numerical control equipment failure state be divided into normally, slight degradation, seriously deterioration and inefficacy one of four states, and definite circular motion precision failure criteria, as shown in table 2.Then, adopt the Lloyd algorithm that the performance characteristic parameter of numerical control equipment is carried out vector quantization, obtain numerical control equipment t circular motion precision output observation sequence 0 constantly.

Table 2 circular motion precision failure criteria

Because numerical control equipment is in normal operating condition under starting condition, the probability distribution of original state vector π then ₀=[1,0,0,0] ^T, status change probability matrix A ₀With observed reading probability matrix B ₀Take random device to choose respectively, have:

$A_0=\left[\begin{array}{cccc}0.3853& 0.2959& 0.2063& 0.1125\\ 0& 0.2198& 0.4064& 0.3738\\ 0& 0& 0.4989& 0.5011\\ 0& 0& 0& 1\end{array}\right],$ ${\mathrm{<figure-callout\; id="0"\; label="B"\; filenames="A200910060949E00181.png"\; state="\{\{state\}\}">B</figure-callout>}}_0=\left[\begin{array}{cccc}0.2343& 0.1668& 0.1296& 0.4693\\ 0.2575& 0.3889& 0.1219& 0.2317\\ 0.1961& 0.2752& 0.0932& 0.4355\\ 0.2710& 0.2810& 0.2115& 0.2365\end{array}\right]$

Thereby set up the status change model λ under the original state ₀=(A ₀, B ₀, π ₀).With circular motion precision output observation sequence O substitution λ ₀=(A ₀, B ₀, π ₀) in, utilize the Baum-Welch algorithm, come the physical training condition transition model by iterative computation.It is 100 that total training iterations is set, and along with the increase of iterations, model output logarithm probable value progressively increases, enter the convergence state of saturation through model after 10 iteration, the output probability value almost no longer increases, and obtains t status change model λ=(A constantly, B, π), that is:

$\mathrm{\&pi;}=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right],$ $A=\left[\begin{array}{cccc}0.9443& 0.0455& 0.0102& 0\\ 0& 0.6142& 0.3857& 0.0001\\ 0& 0& 0.9665& 0.0335\\ 0& 0& 0& 1\end{array}\right],$ $B=\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 0& 0& 0\\ 1& 0& 0& 0\\ 0.1021& 0.4290& 0.1546& 0.3143\end{array}\right]$

Simultaneously, adopt the Forward-Backward algorithm to calculate t circular motion constantly precision output observation sequence O in that λ=(π) the probable value P under (O| λ) is 80.84% for A, B.Reach a conclusion thus: the status change model λ of foundation=(A, B π) can satisfy application request.

At last, with λ=(π) π in and A are updated in the formula 5 for A, B, find the solution this formula, and obtain t state vector P (t)=(P of numerical control equipment constantly through the Laplace inverse transformation ₀(t), P ₁(t), P ₂(t), P ₃(t)), wherein:

P ₀(t)＝exp(0.9443t)

$P_1\left(t\right)=0.2757\exp \left(0.7793t\right)\frac{\exp \left(0.1651t\right)-\exp (-0.1651t)}{2}$

P ₂(t)＝13.77exp(0.9443t)-13.93exp(0.9397t)+0.1633exp(0.6142t)

P ₃(t)＝-8.28exp(0.9443t)+7.74exp(0.9397t)-0.0141exp(0.6142t)+0.5556exp(t)

Then numerical control equipment is R (t)=1-P in t fiduciary level constantly ₃(t), with reference to Fig. 4.As shown in Figure 4, fiduciary level had been lower than 0.8 after numerical control equipment moved about 300 hours, and numerical control equipment is in serious deterioration state, and the probability that breaks down is very high.

Claims (4)

1, a kind of analysis method for reliability of numerical control equipment based on hidden Markov chain may further comprise the steps:

Dynamic property signal when (1) the monitoring numerical control equipment moves extracts the performance characteristic parameter value that more than one reflection reliability of numerical control equipment changes at more than one time point;

(2), make up the forecast model of performance characteristic parameter value according to the performance characteristic parameter value that extracts;

(3) use forecast model prediction T ₀Each performance characteristic parameter value of individual time point carries out vector quantization to these predicted values, obtains the output observation sequence of performance characteristic parameter;

(4), adopt discrete Hidden Markov chain model to determine the running status transition probability matrix of numerical control equipment according to the output observation sequence;

(5) utilize Qie Puman-Andrei Kolmogorov differential equation to set up the relational expression of running status and status change probability matrix, infer that numerical control equipment is in the probability under the different running statuses.

2, a kind of analysis method for reliability of numerical control equipment based on hidden Markov chain according to claim 1 is characterized in that, described step (2) adopts non-stationary autoregression integration running mean algorithm to set up the forecast model of j performance characteristic ginseng value $y_{\mathrm{tj}}=\underset{z=1}{\overset{p}{\mathrm{\&Sigma;}}}{\mathrm{\&phi;}}_z{y}_{(t-z)j}+{\mathrm{\&epsiv;}}_{\mathrm{tj}}-\underset{z=1}{\overset{q}{\mathrm{\&Sigma;}}}{\mathrm{\&phi;}}_z{\mathrm{\&epsiv;}}_{(t-z)j},$ Wherein, t represents the predicted time point, and p is the autoregression order, and q is the running mean order, and φ z is the running mean coefficient, measuring error ${\mathrm{\&epsiv;}}_{\mathrm{tj}}~N(0,{\mathrm{\&sigma;}}_{\mathrm{\&epsiv;}}^2).$

3, a kind of analysis method for reliability of numerical control equipment based on hidden Markov chain according to claim 1 is characterized in that, described step (3) is specially: the output observation sequence of performance characteristic parameter is expressed as $O^j=\{o_1^j,\&CenterDot;\&CenterDot;\&CenterDot;,o_{T_0}^j\},$ $o_{\tilde{t}}^j\&Element;(v_1^j,,v_N^j),$ $\tilde{t}=1,\&CenterDot;\&CenterDot;\&CenterDot;,T_0,j=1,\&CenterDot;\&CenterDot;\&CenterDot;,m,$ M is the performance characteristic number of parameters, and t represents the predicted time point, wherein,

$v_1^j=1,$ $v_N^j=N,$ N is the zone of dispersion number;

$v_{i_0}^j=\mathrm{Index}\left(y_{t_{1}j}\right)=\left\{\begin{array}{c}1,y_{t_{1}j}\&le;\mathrm{partition}(1,j)\\ i_0,\mathrm{partition}(i_0-1,j)\&le;y_{t_{1}j}\&le;\mathrm{partition}(i_0,j);\\ N,\mathrm{partition}(N-1,j)\&le;y_{t_{1}j}\end{array}\right.$

i ₀=2 ..., N-1, t ₁∈ (1 ..., T ₀), Represent that j performance characteristic parameter is at t ₁The predicted value of time point, $\mathrm{partition}(i_1,j)=h_{i_{1}j},$ i ₁=1 ..., N-1,

Represent that j performance characteristic parameter is in i ₁The value of defining of zone of dispersion.

4, a kind of analysis method for reliability of numerical control equipment according to claim 1 based on hidden Markov chain, it is characterized in that, described step (4) is specially: make numerical control equipment comprise 0～k+1 kind running status, the output observation sequence and the predefined original state transition probability matrix A that obtain according to step (3) ₀, initial observation value probability matrix B ₀With initial probability distribution vector π ₀, utilize Bao Mu-Wei Erqi (Baum-Welch) algorithm iteration calculate numerical control equipment constantly status change model λ of t=(π), π represents the probability distribution vector for A, B,

The status change probability matrix $A={\left[a_{\mathrm{cd}}\right]}_{(k+1)\&times;(k+1)}=\left[\begin{array}{ccccc}{a}_{00}& a_{01}& \&CenterDot;\&CenterDot;\&CenterDot;& a_{0k}& a_{0(k+1)}\\ 0& a_{11}& \&CenterDot;\&CenterDot;\&CenterDot;& a_{1k}& a_{1(k+1)}\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ 0& 0& \&CenterDot;\&CenterDot;\&CenterDot;& a_{\mathrm{kk}}& a_{k(k+1)}\\ 0& 0& \&CenterDot;\&CenterDot;\&CenterDot;& 0& a_{(k+1)(k+1)}\end{array}\right],$ a _Cd{ c} is in state { probability of d}, 0≤a constantly at t+1 for numerical control equipment is in state constantly at t _Cd≤ 1, and $\underset{d=0}{\overset{k+1}{\mathrm{\&Sigma;}}}{a}_{\mathrm{cd}}=1;$

The observed reading probability matrix $B={\left[b_{\mathrm{cj}}\right]}_{(k+1)\&times;m}=\left[\begin{array}{cccc}{b}_{01}& b_{02}& \&CenterDot;\&CenterDot;\&CenterDot;& b_{0m}\\ b_{11}& b_{12}& \&CenterDot;\&CenterDot;\&CenterDot;& b_{1m}\\ \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;\\ b_{k1}& b_{k2}& \&CenterDot;\&CenterDot;\&CenterDot;& b_{\mathrm{km}}\\ b_{(k+1)1}& b_{(k+1)2}& \&CenterDot;\&CenterDot;\&CenterDot;& b_{(k+1)m}\end{array}\right],$ b _CjFor numerical control equipment is in state constantly at t { observed reading appears during c}

C=0 ..., the probability of k+1,0≤b _Cj≤ 1, and $\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{b}_{\mathrm{cj}}=1.$

Images