CN109522519A - A kind of dependence evaluation method between multiple performance parameters of ammunition parts - Google Patents
A kind of dependence evaluation method between multiple performance parameters of ammunition parts Download PDFInfo
- Publication number
- CN109522519A CN109522519A CN201811313259.6A CN201811313259A CN109522519A CN 109522519 A CN109522519 A CN 109522519A CN 201811313259 A CN201811313259 A CN 201811313259A CN 109522519 A CN109522519 A CN 109522519A
- Authority
- CN
- China
- Prior art keywords
- parameter
- function
- performance
- formula
- model
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000011156 evaluation Methods 0.000 title claims abstract description 6
- 238000009826 distribution Methods 0.000 claims abstract description 63
- 241000039077 Copula Species 0.000 claims abstract description 60
- 238000000034 method Methods 0.000 claims abstract description 52
- 230000007850 degeneration Effects 0.000 claims abstract description 30
- 238000007476 Maximum Likelihood Methods 0.000 claims abstract description 17
- 230000003412 degenerative effect Effects 0.000 claims abstract description 14
- 238000006731 degradation reaction Methods 0.000 claims description 33
- 230000015556 catabolic process Effects 0.000 claims description 31
- 238000005315 distribution function Methods 0.000 claims description 21
- 230000001186 cumulative effect Effects 0.000 claims description 17
- 238000001514 detection method Methods 0.000 claims description 11
- 238000006243 chemical reaction Methods 0.000 claims description 7
- 238000005457 optimization Methods 0.000 claims description 7
- 101150052147 ALLC gene Proteins 0.000 claims description 6
- 238000013461 design Methods 0.000 claims description 5
- 241000208340 Araliaceae Species 0.000 claims description 4
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 claims description 4
- 235000003140 Panax quinquefolius Nutrition 0.000 claims description 4
- 235000008434 ginseng Nutrition 0.000 claims description 4
- 238000005259 measurement Methods 0.000 description 9
- 238000012360 testing method Methods 0.000 description 8
- 239000011159 matrix material Substances 0.000 description 6
- 230000002068 genetic effect Effects 0.000 description 4
- 241000255789 Bombyx mori Species 0.000 description 3
- 238000013459 approach Methods 0.000 description 3
- 239000010902 straw Substances 0.000 description 3
- 238000006467 substitution reaction Methods 0.000 description 3
- 238000013519 translation Methods 0.000 description 3
- 230000001133 acceleration Effects 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 239000004744 fabric Substances 0.000 description 2
- 241001269238 Data Species 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000000611 regression analysis Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Life Sciences & Earth Sciences (AREA)
- Operations Research (AREA)
- Probability & Statistics with Applications (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Algebra (AREA)
- Evolutionary Biology (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Bioinformatics & Computational Biology (AREA)
- Complex Calculations (AREA)
Abstract
The invention proposes the dependence evaluation methods between a kind of multiple performance parameters of ammunition parts, and the reliability model of two-parameter degeneration is given on the basis of Sklar determines power.Further directed to the model, a kind of two-step method with higher efficiency is proposed to estimate the unknown parameter in model, modeling analysis is carried out to two performance parameter degenerative processes respectively first, estimates the unknown parameter in its edge distribution；Then, edge distribution estimates of parameters is substituted into whole maximum likelihood function.By being compared and analyzed with the measured value that storehouse is stored, reach 85.92% using the model prediction accuracy that 1 this parameter of parameter is calculated in the method for the present invention, the model prediction accuracy of 2 this parameter of parameter reaches 88.60%, and the model prediction accuracy of final two-parameter Copula functional value reaches 93.75%.
Description
Technical field
The present invention relates to a kind of ammunition component degenerative process models.
Background technique
For the product of failure mechanism complexity, there may be multiple performance parameters, and may be deposited between each performance parameter
In dependence.These performance parameters affect the reliability of product simultaneously, any one performance parameter value reaches or beyond design
Prescribed limit can all lead to failure.Now it is necessary to consider shadow of the degeneration to product reliability of each performance parameter simultaneously
It rings, carries out reasonable Comprehensive Assessment from multiple angles, a kind of common method is exactly to be joined using Copula function to more performances
Dependence relation between number is portrayed.
It is necessary to the related and interdependent differences that first gets across before introducing Copula function.There are two words in English
For describing the incidence relation of stochastic variable: Correlation and Dependence.The former corresponds to correlation, and the latter corresponds to
It is interdependent.
Correlation is generally only to work as X for linear relationship, when Y is uncorrelated, might not represent X, Y is mutually indepedent.System
Stochastic variable X is measured frequently with Pearson product-moment correlation coefficient in meter, the linear relationship of Y:
Wherein, ρXYFor Pearson product-moment correlation coefficient, E () and D () are respectively the expectation and variance of stochastic variable.
Pearson product-moment correlation coefficient is only used for metric linear relationship, for example, stochastic variable X obeys standard normal point
Cloth, stochastic variable Y=X2, the scatter plot of (X, Y) is as shown in Figure 1, it is clear that X, Y have strong dependence, and the value of Y can be completely by X
It determines, but correlation coefficient ρXYIt is 0, thus Pearson product-moment correlation coefficient is deposited in terms of measuring stochastic variable incidence relation
In defect.
Interdependent (Dependence) is the antonym of independent (Independence), and dependence refers to that stochastic variable X, Y exist
Do not have independence in probability performance, it had not only contained linear relationship but also had contained non-linear relation, was a kind of more generalized
Description.In real life, non-linear relation is a large amount of existing, therefore, here using " interdependent " description stochastic variable it
Between incidence relation.
Tradition depends on related coefficient based on the modeling method of multivariate joint probability distribution (such as multiple normal distribution), for two dimension
For normal random variable (X, Y), the related coefficient of X, Y are equivalent to X for 0, and Y is mutually indepedent, this feature causes it not applicable
In non-linear interdependent occasion.And the method based on Copula function can handle non-linear interdependent data, ask to solve this
Topic provides a strong tool.
The definition and Sklar theorem of dimensional Co pula function
Copula is intended that " connection ", and the modeling approach using Copula function is to decompose the Joint Distribution of random vector
For edge distribution and Copula function two parts, edge distribution is connected by Copula function, as Joint Distribution.This
So that the modeling of Joint Distribution is more flexible, it is no longer limited to existing classical Joint Distribution.The previously described random vector (X, Y)
Joint Distribution, generally using the classical Joint Distribution such as Two dimension normal distribution, this kind of Joint Distribution is single often through one
Family of distributions constructs to obtain, and the structure of the family of distributions is also relatively easy, and otherwise the expression formula of Joint Distribution can be excessively complicated.This
Outside, define that Joint Distribution is equivalent to implicitly define edge distribution, this is being many times unreasonable.It is based on
The Joint Distribution modeling method of Copula function can respectively specify that the form of edge distribution and co-ordinative construction, so that solution to model
It is stronger to release power.
The Joint Distribution that can be used for constructing Multivariate Random Vector on Copula function theory, due to Copula function itself
In some limitations of higher-dimension situation, dimensional Co pula function is only introduced.The definition of dimensional Co pula function and correlation theorem are general
It includes as follows:
It defines 1 dimensional Co pula function and is defined in [0,1]2On function C (u, v), it meets following two condition:
1 couple of any u, v ∈ [0,1], C (u, 0)=C (0, v)=0, C (u, 1)=u, C (1, v)=v is set up.
2 couples of arbitrary u1≤u2,v1≤v2, u1,u2,v1,v2∈ [0,1], has
C(u2,v2)-C(u2,v1)-C(u1,v2)+C(u1,v1)≥0
If there is 1 joint distribution function H () of theorem continuous boundary to be distributed F1(·),F2() certainly exists one only
One Copula function C (), to arbitrary x1,x2∈ R, has
H(x1,x2)=C (F1(x1),F2(x2))\*MERGEFORMAT(2)
It is the Sklar theorem that Sklar is proposed in nineteen fifty-nine above.The theorem proving existence of Copula function and only
One property, gives sufficient and necessary condition existing for Copula function, illustrates to solve the problems, such as that dependence is with Copula function
It is feasible.Sklar theorem shows that any one two-dimentional Joint Distribution can be decomposed into a Copula function and corresponding two
A edge distribution, this makes Two dimensional Distribution modeling more flexible.The theorem is Copula function in multi-variate statistical analysis field
Using having established theoretical basis.
Common Copula function
Common Copula function has oval a small bundle of straw, etc. for silkworms to spin cocoons on Copula function, FGM (Farlie-Gumbel-Morgenstern) a small bundle of straw, etc. for silkworms to spin cocoons on
Copula function and Archimedes's a small bundle of straw, etc. for silkworms to spin cocoons on Copula function.Gaussian distribution is one kind of family of ellipses distribution, although this distribution
There is no closed expression, but relatively simple when construction multivariate joint probability distribution, so being widely used in engineering.It is oval
Race is distributed maximum the disadvantage is that requiring its edge distribution all is same family of ellipses, in order to overcome this difficulty, proposes ellipse
Race Copula.It does not have strict requirements to edge distribution, may belong to different family of ellipses.
Gaussian Copula function is defined as:
Wherein, Σ indicates that the element on diagonal line is 1, and other elements are the symmetric positive definite matrix of ρ；ΦρIndicate phase relation
Matrix number is the standard multiple normal distribution cumulative distribution function of Σ；Φ-1() indicates standardized normal distribution cumulative distribution function
Inverse function；I indicates unit matrix.
Summary of the invention
The object of the present invention is to provide the dependence evaluation methods between a kind of multiple performance parameters of ammunition parts, to realize
This purpose, scheme are as follows:
It is characterized by:
It is assumed that establishing the two-parameter degradation model based on Copula under the precondition of Copula parameter constant, and provide phase
The method for parameter estimation answered:
Assuming that certain product has two performance characteristic parameters, the available degeneration rail of the degenerative process of each performance characteristic parameter
Mark model or random process model description, are denoted as { X1(t),X2(t)}；
{X1(t),X2(t) } joint distribution function is
H(X1(t),X2(t))=C (F1(X1(t)),F2(X2(t))\*MERGEFORMAT(5)
Wherein, Fp(XpIt (t)) is Xp(t) cumulative distribution function, p=1,2；
For the product with two performance characteristic parameters, Reliability Function is
R (t)=C (R1(t),R2(t))\*MERGEFORMAT(6)
Wherein, Rp(t) be p-th of performance parameter Reliability Function, p=1,2；
When any one performance parameter value reaches or when threshold value beyond design code, will lead to product failure；Fail threshold
Value is denoted as { D1,D2}；Therefore, for the product with two performance characteristic parameters, reliability is
R (t)=P ((X1(t),X2(t)) ∈ S)=P (X1(t) < D1,X2(t) < D2)\*MERGEFORMAT(7)
S represents security domain in formula, and when all properties parameter is respectively positioned in security domain S, product is worked normally；
Due to P (X1(t) < D1)=F1(D1| t), P (X2(t) < D2)=F2(D2|t)；Formula (7) can be write as
R (t)=P (X1(t) < D1,X2(t) < D2)
=C (F1(D1|t),F2(D2|t))\*MERGEFORMAT(8)
=C (R1(t),R2(t))
F in formulap(XpIt (t)) is Xp(t) cumulative distribution function, p=1,2；Rp(t) be p-th of performance parameter reliability
Function, p=1,2；
Particularly, if performance characteristic parameter is mutually indepedent, the Reliability Function in formula (8) can simplify for
R (t)=R1(t)R2(t)\*MERGEFORMAT(9)
Wherein RpIt (t) is the marginal Reliability Function of p-th of performance characteristic parameter, p=1,2；
For two performance degenerative process { X1(t),X2(t) }, enabling it in the Reliability Function of t moment is respectively u=R1
(t) and v=R2(t)；Then
(P1) for allC (u, 0)=C (0, v)=0；
(P2) for allC (u, 1)=u, C (1, v)=v；
(P3) for allMax (u+v-1,0) < C (u, v) < min (u, v)；
Property (P1) illustrates that any one performance parameter is more than that failure threshold can all lead to product failure, property (P2) explanation
If one of performance parameter is not degenerated, the reliability of product will be determined by another performance parameter completely；Property (P3) description
Fr é chet-Hoeffding section, by R (t)=C (R1(t),R2(t)) it substitutes into property (P3) and obtains Fr é chet-
The section Hoeffding are as follows:
max(R1(t)+R2(t)-1,0)≤R(t)≤min(R1(t),R2(t))\*MERGEFORMAT(10)
The property illustrates that the reliability of two-parameter performance degradation failure type product is centainly in above formula section；
Estimate model parameter method:
Step 1: solving edge distribution parameter
1) for the covariant model in nonlinear drift Wiener-Hopf equation, parameter to be estimated and contiguous function μzIt is related,
When contiguous function is μz=ae-b/zWhen, parameter to be asked is βp=[ap,bp,σp], p=1,2；Degeneration increment xpkijProbability density
Function is
Log-likelihood function is
In formulaIt is the degeneration increment observed, Δ (tkij)=tkij-tki,j-1It is twice
The interval time of detection, τkij=τ (tkij)-τ(tki,j-1) be conversion after interval time,For the probability of standardized normal distribution
Density function；
2) for the covariant model in inverse Gaussian process, parameter to be estimated is βp=[ap,bp,σp], p=1,2,
Middle ap,bpIt is contiguous function μzIn parameter；Degeneration increment xpkijProbability density function be
Log-likelihood function is
In formulaIt is the degeneration increment observed, τkij=τ (tkij)-τ(tki,j-1) it is to turn
Detection interval time after changing,It is about ap,bpFunction；
Step 2: solving Copula function parameter；
The first step has found out the parameter of two edge distributionsSubstituted into formula (32), available pass
In the likelihood function of parameter alpha:
In formulaWithIt is the cumulative distribution function and probability density function of p-th of edge distribution respectively,Be C (；Density function α), xpkijIt is that the degeneration observed increases
Amount, p=1,2；
After saving constant term, formula (37) can simplify are as follows:
The estimation parameter value of parameter alpha can indicate are as follows:
Θ is parameter space in formula；
Since, only comprising one to two parameters, maximum likelihood estimation can pass through Newton- in likelihood function
The optimization algorithms such as Raphson acquire.
Figure of description
Fig. 1 is that stochastic variable X obeys standardized normal distribution, stochastic variable Y=X2, the scatter plot of (X, Y), wherein phase relation
Number is 0, but X and Y be not independent；
Fig. 2 is equivalent serial model；
Fig. 3 is 1 Degradation path of parameter at 70 DEG C；
Fig. 4 is 2 Degradation path of parameter at 70 DEG C；
Fig. 5 is 1 Degradation path of parameter at 85 DEG C；
Fig. 6 is 2 Degradation path of parameter at 85 DEG C；
Fig. 7 is 1 Degradation path of parameter at 100 DEG C；
Fig. 8 is 2 Degradation path of parameter at 100 DEG C
Fig. 9 is the 1 degraded data average value of parameter (70 DEG C) under power function fitting；
Figure 10 is the 1 degraded data average value of parameter (85 DEG C) under power function fitting；
Figure 11 is the 1 degraded data average value of parameter (100 DEG C) under power function fitting.
Specific embodiment
With reference to the accompanying drawing and specific embodiment the present invention will be further described.
1.1 reliability model
The present invention describes the dependence between two performance parameters using Copula function.It initially sets up and is increased based on non-stationary
The one-parameter degradation model of amount process, then it is assumed that being established under the precondition of Copula parameter constant double based on Copula
Parameter degradation model, and provide corresponding method for parameter estimation.
For a tool there are two the product of critical performance parameters, the present invention has done following hypothesis:
The degeneration behavior of two performance parameters can be modeled with the non-stationary incremental process hereafter established.
There are dependence, Dependence Structures can be described with Copula function for two performance parameters.
The range that any one performance parameter value reaches or allows beyond design all will lead to product failure.
The parameter of Copula function is not at any time and stress variation.
Assuming that certain product has two performance characteristic parameters, the available degeneration rail of the degenerative process of each performance characteristic parameter
Mark model or random process model description, are denoted as { X1(t),X2(t)}.According to theorem 1 it is found that { X1(t),X2(t) } joint point
Cloth function is
H(X1(t),X2(t))=C (F1(X1(t)),F2(X2(t))\*MERGEFORMAT(5)
Wherein, Fp(XpIt (t)) is Xp(t) cumulative distribution function, p=1,2.
For the product with two performance characteristic parameters, Reliability Function is for proposition 1
R (t)=C (R1(t),R2(t))\*MERGEFORMAT(6)
Wherein, Rp(t) be p-th of performance parameter Reliability Function, p=1,2.
For the product with two performance parameters, reliability block diagram can be equivalent to series connection mould shown in Fig. 2
Type.
When any one performance parameter value reaches or when threshold value beyond design code, will lead to product failure.Fail threshold
Value is denoted as { D1,D2}.Therefore, for the product with two performance characteristic parameters, reliability is
R (t)=P ((X1(t),X2(t)) ∈ S)=P (X1(t) < D1,X2(t) < D2)\*MERGEFORMAT(7)
S represents security domain in formula, and when all properties parameter is respectively positioned in security domain S, product is worked normally.
Due to P (X1(t) < D1)=F1(D1| t), P (X2(t) < D2)=F2(D2|t).Formula (7) can be write as
R (t)=P (X1(t) < D1,X2(t) < D2)
=C (F1(D1|t),F2(D2|t))
=C (R1(t),R2(t))
\*MERGEFORMAT(8)
F in formulap(XpIt (t)) is Xp(t) cumulative distribution function, p=1,2.Rp(t) be p-th of performance parameter reliability
Function, p=1,2.
Particularly, if performance characteristic parameter is mutually indepedent, the Reliability Function in formula (8) can simplify for
R (t)=R1(t)R2(t)\*MERGEFORMAT(9)
Wherein RpIt (t) is the marginal Reliability Function of p-th of performance characteristic parameter, p=1,2.
Property 1Copula function has good statistical property, for two performance degenerative process { X1(t),X2(t) } it, enables
It is respectively u=R in the Reliability Function of t moment1(t) and v=R2(t).Then (P1) is for allC (u, 0)=C
(0, v)=0；
(P2) for allC (u, 1)=u, C (1, v)=v；
(P3) for allMax (u+v-1,0) < C (u, v) < min (u, v).
Property (P1) illustrates that any one performance parameter is more than that failure threshold can all lead to product failure, property (P2) explanation
If one of performance parameter is not degenerated, the reliability of product will be determined by another performance parameter completely.Property (P3) description
Fr é chet-Hoeffding section, by R (t)=C (R1(t),R2(t)) it substitutes into property (P3) and obtains Fr é chet-
The section Hoeffding are as follows:
max(R1(t)+R2(t)-1,0)≤R(t)≤min(R1(t),R2(t))\*MERGEFORMAT(10)
The property illustrates that the reliability of two-parameter performance degradation failure type product is centainly in above formula section.
Two-parameter degradation model parameter Estimation
Method for parameter estimation under 1.1 specified stress
Assuming that putting into n sample in total, each sample existsMoment is detected, the pth of i-th of sample
Performance parameter is in tijThe performance detection value at moment is denoted asI=1,2 ..., n, j=0,1,2 ..., mi, p=1,2；
Wherein n is sample size, miIt is the detection number of i-th of sample.Degraded data can be expressed as form:
Remember that α is the parameter in Copula function, βpFor the model parameter of p-th of performance degradation process,
For i-th of sample p-th of performance parameter in time interval [ti,j-1,tij] in degeneration increment, whereinIt degenerates
The joint probability density function of increment is
F in formulap(·；βp) and fp(·；βp) be respectively p-th of edge distribution cumulative distribution function and probability density letter
Number,Be C (；Density function α).
According to formula (12) and independent increment property, the likelihood function that test data can be obtained is
Log-likelihood function is accordingly
It include multiple unknown parameters in formula (14), by taking simplest model as an example, β1,β2It respectively include three parameters, α includes
One to two parameters.It is by no means easy to directly adopt Maximum Likelihood Estimation solving model parameter, it is even infeasible.Observation
It was found that parameter to be estimated consists of two parts: first part is parameter beta=[β in edge distribution1,β2], second part is
Parameter alpha in Copula function.The method of two steps estimation of the invention, is successively fitted each edge distribution first, solves edge distribution
In parameter beta=[β1,β2]；Then Copula structure appropriate is selected, by the obtained edge distribution estimates of parameters of the first stepLikelihood function shown in substitution formula (14) obtains the estimation of Copula parameter alpha using Maximum Likelihood Estimation
Value.Since, all only comprising a small amount of parameter, difficulty in computation is remarkably decreased, and solution efficiency will also get a promotion in each step.Xu, Joe
The validity and asymptotic normality of this two steps estimation are had studied with Wang Fang et al., and are carried out with whole Maximum Likelihood Estimation
Comparison, it is indicated that this method efficiency with higher, edge distribution estimates of parameters and Copula estimates of parameters all ten tap
Nearly actual parameter.
Step 1: solving edge distribution parameter.
The parameter of two edge distributions is estimated respectively
1) for the naive model X (t) in nonlinear drift Wiener-Hopf equation=σpB(t)+μpτp(t) for, βp=[μp,
σp,γp], p=1,2, wherein γpFor τp(t) unknown parameter in, such asDegeneration increment xpijProbability density function
For
Log-likelihood function is
In formulaIt is the degeneration increment observed, τij=τ (tij)-τ(ti,j-1) be conversion after
Detection interval time twice, Δ tij=tij-ti,j-1It is the interval time detected twice,For the probability of standardized normal distribution
Density function.
2) for the random-effect model in nonlinear drift Wiener-Hopf equation
For, βp=[ωp,δp,σp,γp], p=1,2, γpFor τp(t)
In unknown parameter, such asDegeneration increment xpijProbability density function be
Log-likelihood function is
In formulaIt is the degeneration increment observed, τij=τ (tij)-τ(ti,j-1) be conversion after
The interval time detected twice, Δ tij=tij-ti,j-1It is the interval time detected twice.
3) for naive model IG (μ τ (t), the σ τ in inverse Gaussian process2(t)) for, βp=[μp,σp,γp], p=1,
2, γpFor τp(t) unknown parameter in, such asDegeneration increment xpijProbability density function be
Log-likelihood function is
In formulaIt is the degeneration increment observed, τij=τ (tij)-τ(ti,j-1) be conversion after
Detection interval time.
4) random-effect model IG (μ τ (t), the σ τ in inverse Gaussian process2(t)),μ-1~TN (ω, δ-2) for, βp=
(σp,γp,ωp,δp), p=1,2, γpFor τp(t) unknown parameter in, such asDegeneration increment xpiJoint probability
Density function is
Log-likelihood function is
In formulaIt is the degeneration increment observed, τij=τ (tij)-τ(ti,j-1),
Step 2: solving Copula function parameter.
The first step has found out the parameter of two edge distributionsSubstituted into formula (14) available pass
In the likelihood function of parameter alpha:
In formulaWithIt is the cumulative distribution function and probability density function of p-th of performance parameter respectively,Be C (；Density function α), xpijIt is that the degeneration observed increases
Amount, p=1,2.
Latter half is unrelated with parameter alpha in formula (23), does not influence final result, can cast out.New likelihood function are as follows:
The maximum likelihood estimation of parameter alpha can indicate are as follows:
Θ is parameter space in formula.
It only include a parameter alpha in likelihood function, maximum likelihood estimation can for one-parameter Copula function
To be easy to find out.By taking Frank Copula as an example, Copula density function are as follows:
Copula density function substitution formula (25) is obtained:
In other one-parameter Copula functionsIt can also be solved by similar approach.
It include two parameter alpha=[α in likelihood function for two-parameter Copula function1,α2].With Joe-
For Clayton Copula, Copula density function are as follows:
Copula density function in formula (28) is substituted into formula (25) to obtain:
Parameter in other two-parameter Copula functionsIt can also be solved by similar approach.
It only include two parameters in formula (29), optimal solution can be acquired using common optimization algorithm, such as Newton-
Raphson algorithm.
Method for parameter estimation under 2.2 Constant Acceleration stress
The method for parameter estimation under specified stress has been discussed above, on this basis, has continued discussing under Constant Acceleration stress
Method for parameter estimation.Assuming that there is a product that there is two performance characteristic parameters, the degenerative process of each performance characteristic parameter
Available degradation path model or random process model description.Assuming that stress level zkLower investment nkA sample, total number of samples amount areL is stress total number.Sample is in tkijMoment is detected, and p-th of performance parameter of i-th of sample exists under k-th of stress
tkijThe performance detection value at moment is denoted asK=1,2 ..., l, i=1,2 ..., nk, j=0,1,2 ..., mki, p=
1,2。nkIt is the sample size under k-th of stress level, mkiIt is the detection number of i-th of sample under k-th of stress level.It degenerates
Data can be expressed as form as shown in table 1:
The two-parameter product accelerated degradation test data of table 1
Remember that α is the parameter in Copula function, βpFor the model parameter of p-th of performance degradation process,
It is performance parameter in time interval [tki,j-1,tkij] in the degeneration increment that observes, whereinThe connection of degeneration increment
Closing probability density function is
F in formulap(·；βp) and fp(·；βp) be respectively p-th of edge distribution cumulative distribution function and probability density function,Be C (；Density function α).
According to independent increment property, the likelihood function that test data can be obtained is
Log-likelihood function is accordingly
Equally model parameter is estimated using two step estimation methods.
Step 1: solving edge distribution parameter.
3) for the covariant model in nonlinear drift Wiener-Hopf equation, parameter to be estimated and contiguous function μzIt is related,
When contiguous function is μz=ae-b/zWhen, parameter to be asked is βp=[ap,bp,σp], p=1,2.Degeneration increment xpkijProbability density
Function is
Log-likelihood function is
In formulaIt is the degeneration increment observed, Δ (tkij)=tkij-tki,j-1It is twice
The interval time of detection, τkij=τ (tkij)-τ(tki,j-1) be conversion after interval time,For the probability of standardized normal distribution
Density function.
4) for the covariant model in inverse Gaussian process, parameter to be estimated is βp=[ap,bp,σp], p=1,2,
Middle ap,bpIt is contiguous function μzIn parameter.Degeneration increment xpkijProbability density function be
Log-likelihood function is
In formulaIt is the degeneration increment observed, τkij=τ (tkij)-τ(tki,j-1) it is to turn
Detection interval time after changing,It is about ap,bpFunction.
Step 2: solving Copula function parameter.
The first step has found out the parameter of two edge distributionsSubstituted into formula (32), available pass
In the likelihood function of parameter alpha:
In formulaWithIt is the cumulative distribution function and probability density function of p-th of edge distribution respectively,Be C (；Density function α), xpkijIt is the degeneration observed
Increment, p=1,2.
After saving constant term, formula (37) can simplify are as follows:
The estimation parameter value of parameter alpha can indicate are as follows:
Θ is parameter space in formula.
Since, only comprising one to two parameters, maximum likelihood estimation can relatively easily pass through in likelihood function
The optimization algorithms such as Newton-Raphson acquire.
The assessment of 2 ammunition assembly reliabilities
2.1 accelerated degradation test data
For ammunition component, constant stress temperature accelerated degradation test is carried out, has chosen 70 DEG C, 85 DEG C, 100 DEG C and be used as 3
A accelerated stress is horizontal, and 5 samples are respectively arranged at each temperature.For each sample, its parameter 2 of periodic measurement and parameter 1 this
Two performance degradation parameters measure 9,12,9 times respectively under three stress.Complete performance degradation parameter measurement data such as table
Shown in 2-4 and Fig. 3-8.
Performance parameter degraded data at 2 70 DEG C of table
Performance parameter degraded data at 3 85 DEG C of table
Performance parameter degraded data at 4 100 DEG C of table
According to two-parameter model set forth above and two-step method method for parameter estimation, next first to parameter 2 and parameter 1
The two performance parameters carry out translation processing, i.e., are 0 by the translation of initial drift amount, next carry out respectively to the parameter after translation
Modeling is then based on Copula function and models to the dependence relation of two parameter, and carries out parameter Estimation and reliability assessment.
2.2 parameters 1
For this single performance parameter of parameter 1, modeling analysis is carried out using common degradation path model:
Y (t)=λ1γ1Λ1(t)+ε1\*MERGEFORMAT(40)
Wherein, Λ1It (t) is nonlinear function, γ1It is covariant parameter,For error term.
Since the reliability test of ammunition component is completed under temperature accelerated stress, common Allan Buddhist nun is used
This black model models covariant parameter, i.e.,Then formula (40) becomes
ByIt is found that t at any one time,If its probability
Density function is fY(y):
Wherein,For the probability density function of standardized normal distribution,For Λ1(t) unknown parameter in.
It has furthermore been found that under same temperature stress, E [Y (t)] ∝ Λ1(t), therefore by the institute under simple stress
There is the average value of 5 sample degraded datas to carry out regression analysis, can tentatively judge Λ1(t) citation form.Here with linear
Alternately model is analyzed for three kinds of function, exponential function and power function forms, as a result as shown in table 5 and Fig. 9 to Figure 11, hair
The fitting effect of existing power function is preferable, therefore selects here
The coefficient R of 1 degraded data average value of parameter under the different fitting functions of table 52
Λ is being determined1(t) after citation form, so that it may using Maximum Likelihood Estimation in degradation model not
Know that parameter is estimated, maximum likelihood function are as follows:
Wherein, (ti,j,k,yi,j,k) it is time point and parameter 1 of i-th of sample in jth time measurement under k-th of stress
Amount, K are total stress number, nkFor the sample number under k-th of stress, miFor the pendulous frequency of i-th of sample.
Parameters are obtained to above-mentioned maximum likelihood function maximizing using intelligent optimization methods such as genetic algorithms
Estimated value is as shown in the table:
6 parameter of table, 1 degradation model parameter Estimation
It follows that ammunition parameterWherein ε1~N (0,
0.22072)。
Therefore, 1 average value of parameter stored under 25.7 ° by model prediction 10 years is
Wherein, -1.2333 be parameter 1 average value of all samples in on-test.
And mechanical 1 amount of ammunition component parameter of storehouse storage (25.7 °) is sent out by actual measurement two, record is stored 10 years
Measured value is respectively -0.7 and -1.1, i.e., the actual measurement of parameter 1 average value is
Therefore, for this parameter of parameter 1, the precision of prediction of model
2.3 parameters 2
For this performance degradation parameter of parameter 2, handled using similar degradation model, i.e. X (t)=λ2γ2Λ2
(t)+ε2, whereinFor nonlinear function,It is covariant related with accelerated stress,For error term.
T at any one time,If its probability density function is fX(x):
The unknown parameter in degradation model is estimated also with Maximum Likelihood Estimation, maximum likelihood function
Are as follows:
Wherein, (ti,j,k,xi,j,k) it is time point and parameter 2 of i-th of sample in jth time measurement under k-th of stress
Amount, K are total stress number, nkFor the sample number under k-th of stress, miFor the pendulous frequency of i-th of sample.
Parameters are obtained to above-mentioned maximum likelihood function maximizing using intelligent optimization methods such as genetic algorithms
Estimated value is as shown in the table:
7 parameter of table, 2 degradation model parameter Estimation
| Parameter | λ2 | σ2 | η2 | β2 |
| Estimated value | 7.6911 | 0.1808 | 0.2259 | 0.0075 |
It follows that parameterWherein ε2~N (0,0.18082).Cause
This, 2 average value of parameter stored under 25.7 ° 10 years by model prediction is
Wherein, -0.4 is parameter 2 average value of all samples in on-test.
And 2 amount of parameter of storehouse storage (25.7 °) is sent out by actual measurement two, record stores 10 years measured value difference
For -0.5 and -0.4, i.e. the actual measurement of parameter 2 average value is
Therefore, for this parameter of parameter 1, the precision of prediction of model
2.4 two-parameter modelings
When considering the correlation of parameter 2 and parameter 1, two-parameter modeling is carried out using Gaussian Copula function,
Citation form is
Wherein, Σ indicates that the element on diagonal line is 1, and other elements are the symmetric positive definite matrix of ρ；ΦρIndicate phase relation
Matrix number is the standard multiple normal distribution cumulative distribution function of Σ；Φ-1() indicates standardized normal distribution cumulative distribution function
Inverse function；I indicates unit matrix.
For the two-dimentional Gaussian Copula in two-parameter situation, functional form can be write as
Wherein,Indicate that related coefficient is the standard dyadic normpdf of ρ.
Based on this Gaussian Copula function and parameter 1 and 2 degradation model of parameter, available total likelihood function is
Wherein FY(xijk),FX(yijk) be respectively parameter 1 and parameter 2 cumulative distribution function.The likelihood function is taken pair
Number, obtains
By two-step method set forth above it is found that can first estimate the unknown ginseng in 2 two edge distributions of parameter 1 and parameter
Number, as a result as shown in front 6 and table 7；It is as follows by new likelihood function is obtained after these parameter estimation results substitution (56).
The maximum value of above-mentioned likelihood function is solved by intelligent optimization methods such as genetic algorithms, can further be estimated
Unknown parameter ρ=0.0089 in Gaussian Copula function.
Therefore, model prediction in 10 years is stored for the two-parameter degenerative process there are dependence relation, under 25.7 °
Copula mean value functions are
The ammunition component parameter 1 and 2 measured value of parameter under storehouse storage be respectively (- 0.5, -0.7) and (- 0.4, -
1.1), then Copula function actual measurement average value is
Therefore, for this two-parameter degenerative process, the precision of prediction of model
For there are the two-parameter degenerative processes of dependence relation, present invention firstly provides a kind of multiple property of ammunition parts
Dependence evaluation method between energy parameter, gives the reliability model of two-parameter degeneration on the basis of Sklar determines power.Into one
Step is directed to the model, proposes a kind of two-step method with higher efficiency to estimate the unknown parameter in model, right respectively first
Two performance parameter degenerative processes carry out modeling analysis, estimate the unknown parameter in its edge distribution；Then, by edge distribution
Estimates of parameters substitutes into whole maximum likelihood function, can estimate to obtain the unknown parameter in Copula function.
For actual ammunition component degenerative process, its crucial parameter 1 and parameter 2 two interdependent degeneration ginsengs are considered
Number is carried out modeling analysis using the two-parameter degradation model of proposition and is asked respectively using genetic algorithm Sum Maximum Likelihood Estimate method
Obtained the unknown parameter in the edge distribution parameter and Gaussian Copula function of two performance amount of degradations.By with library
The measured value of room storage compares and analyzes, and the model prediction accuracy that 1 this parameter of parameter is calculated reaches 85.92%, ginseng
The model prediction accuracy of 2 this parameter of number reaches 88.60%, and the model prediction accuracy of final two-parameter Copula functional value reaches
To 93.75%.
Claims (1)
1. the dependence evaluation method between a kind of multiple performance parameters of ammunition parts, it is characterised in that:
It is assumed that establishing the two-parameter degradation model based on Copula under the precondition of Copula parameter constant, and provide corresponding
Method for parameter estimation:
Assuming that certain product has two performance characteristic parameters, the available Degradation path mould of the degenerative process of each performance characteristic parameter
Type or random process model description, are denoted as { X1(t),X2(t)}；{X1(t),X2(t) } joint distribution function is
H(X1(t),X2(t))=C (F1(X1(t)),F2(X2(t))
Wherein, Fp(XpIt (t)) is Xp(t) cumulative distribution function, p=1,2；
For the product with two performance characteristic parameters, Reliability Function is
R (t)=C (R1(t),R2(t))
Wherein, Rp(t) be p-th of performance parameter Reliability Function, p=1,2；
When any one performance parameter value reaches or when threshold value beyond design code, will lead to product failure；Failure threshold note
For { D1,D2}；Therefore, for the product with two performance characteristic parameters, reliability is
R (t)=P ((X1(t),X2(t)) ∈ S)=P (X1(t) < D1,X2(t) < D2)
S represents security domain in formula, and when all properties parameter is respectively positioned in security domain S, product is worked normally；
Due to P (X1(t) < D1)=F1(D1| t), P (X2(t) < D2)=F2(D2|t)；Formula (7) can be write as
R (t)=P (X1(t) < D1,X2(t) < D2)
=C (F1(D1|t),F2(D2|t))
=C (R1(t),R2(t))
F in formulap(XpIt (t)) is Xp(t) cumulative distribution function, p=1,2；Rp(t) be p-th of performance parameter Reliability Function,
P=1,2；
Particularly, if performance characteristic parameter is mutually indepedent, the Reliability Function in formula (8) can simplify for
R (t)=R1(t)R2(t)
Wherein RpIt (t) is the marginal Reliability Function of p-th of performance characteristic parameter, p=1,2；
For two performance degenerative process { X1(t),X2(t) }, enabling it in the Reliability Function of t moment is respectively u=R1(t) and v
=R2(t)；Then
(P1) for allC (u, 0)=C (0, v)=0；
(P2) for allC (u, 1)=u, C (1, v)=v；
(P3) for allMax (u+v-1,0) < C (u, v) < min (u, v)；
Property (P1) illustrates that any one performance parameter is more than that failure threshold can all lead to product failure, and property (P2) is if illustrate it
In a performance parameter do not degenerate, the reliability of product will be determined by another performance parameter completely；Property (P3) describes Fr é
The section chet-Hoeffding, by R (t)=C (R1(t),R2(t)) it substitutes into property (P3) and obtains the section Fr é chet-Hoeffding
Are as follows: max (R1(t)+R2(t)-1,0)≤R(t)≤min(R1(t),R2(t))
The property illustrates that the reliability of two-parameter performance degradation failure type product is centainly in above formula section；
Estimate model parameter method:
Step 1: solving edge distribution parameter
1) for the covariant model in nonlinear drift Wiener-Hopf equation, parameter to be estimated and contiguous function μzIt is related, work as connection
Function is μz=ae-b/zWhen, parameter to be asked is βp=[ap,bp,σp], p=1,2；Degeneration increment xpkijProbability density function be
Log-likelihood function is
In formulaIt is the degeneration increment observed, Δ (tkij)=tkij-tki,j-1It is to detect twice
Interval time, τkij=τ (tkij)-τ(tki,j-1) be conversion after interval time,For the probability density of standardized normal distribution
Function；
2) for the covariant model in inverse Gaussian process, parameter to be estimated is βp=[ap,bp,σp], p=1,2, wherein ap,
bpIt is contiguous function μzIn parameter；Degeneration increment xpkijProbability density function be
Log-likelihood function is
In formulaIt is the degeneration increment observed, τkij=τ (tkij)-τ(tki,j-1) be conversion after
Detection interval time,It is about ap,bpFunction；
Step 2: solving Copula function parameter；
The first step has found out the parameter of two edge distributionsIt is substituted into formula (32), it is available about ginseng
The likelihood function of number α:
In formulaWithIt is the cumulative distribution function and probability density function of p-th of edge distribution respectively,Be C (；Density function α), xpkijIt is that the degeneration observed increases
Amount, p=1,2；
After saving constant term, formula (37) can simplify are as follows:
The estimation parameter value of parameter alpha can indicate are as follows:
Θ is parameter space in formula；
Since, only comprising one to two parameters, maximum likelihood estimation can pass through Newton-Raphson etc. in likelihood function
Optimization algorithm acquires.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811313259.6A CN109522519B (en) | 2018-11-06 | 2018-11-06 | Dependency evaluation method among multiple performance parameters of ammunition component |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811313259.6A CN109522519B (en) | 2018-11-06 | 2018-11-06 | Dependency evaluation method among multiple performance parameters of ammunition component |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109522519A true CN109522519A (en) | 2019-03-26 |
| CN109522519B CN109522519B (en) | 2023-02-03 |
Family
ID=65774292
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811313259.6A Active CN109522519B (en) | 2018-11-06 | 2018-11-06 | Dependency evaluation method among multiple performance parameters of ammunition component |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109522519B (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111160713A (en) * | 2019-12-06 | 2020-05-15 | 中国南方电网有限责任公司超高压输电公司广州局 | Composite insulator reliability assessment method based on multidimensional joint distribution theory |
| CN113361142A (en) * | 2021-07-12 | 2021-09-07 | 西安工程大学 | Spinning machine reliability evaluation method with multiple parameter degradation dependence |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103020423A (en) * | 2012-11-21 | 2013-04-03 | 华中科技大学 | Copula-function-based method for acquiring relevant characteristic of wind power plant capacity |
| CN103746370A (en) * | 2013-12-20 | 2014-04-23 | 河海大学 | Wind-power-plant reliability modeling method |
| CN106227906A (en) * | 2016-05-20 | 2016-12-14 | 广州韵脉质量技术服务有限公司 | A kind of appraisal procedure of the intelligent manufacturing equipment reliability analyzed based on performance degradation |
| CN106548418A (en) * | 2016-12-09 | 2017-03-29 | 华北电力大学（保定） | Power system small interference stability appraisal procedure |
| US20170300911A1 (en) * | 2016-04-13 | 2017-10-19 | Abdullah Abdulaziz I. Alnajem | Risk-link authentication for optimizing decisions of multi-factor authentications |
| CN107436963A (en) * | 2016-12-19 | 2017-12-05 | 鲍兆伟 | A kind of O-shaped rubber seal life-span prediction method based on the polynary degeneration of Copula functions |
| CN107506337A (en) * | 2017-10-12 | 2017-12-22 | 中国人民解放军海军航空工程学院 | Reliability statistics estimating method based on polynary acceleration degraded data |
-
2018
- 2018-11-06 CN CN201811313259.6A patent/CN109522519B/en active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103020423A (en) * | 2012-11-21 | 2013-04-03 | 华中科技大学 | Copula-function-based method for acquiring relevant characteristic of wind power plant capacity |
| CN103746370A (en) * | 2013-12-20 | 2014-04-23 | 河海大学 | Wind-power-plant reliability modeling method |
| US20170300911A1 (en) * | 2016-04-13 | 2017-10-19 | Abdullah Abdulaziz I. Alnajem | Risk-link authentication for optimizing decisions of multi-factor authentications |
| CN106227906A (en) * | 2016-05-20 | 2016-12-14 | 广州韵脉质量技术服务有限公司 | A kind of appraisal procedure of the intelligent manufacturing equipment reliability analyzed based on performance degradation |
| CN106548418A (en) * | 2016-12-09 | 2017-03-29 | 华北电力大学（保定） | Power system small interference stability appraisal procedure |
| CN107436963A (en) * | 2016-12-19 | 2017-12-05 | 鲍兆伟 | A kind of O-shaped rubber seal life-span prediction method based on the polynary degeneration of Copula functions |
| CN107506337A (en) * | 2017-10-12 | 2017-12-22 | 中国人民解放军海军航空工程学院 | Reliability statistics estimating method based on polynary acceleration degraded data |
Non-Patent Citations (7)
| Title |
|---|
| MARIUS HOFERT等: "Likelihood inference for Archimedean copulas in high dimensions under known margins", 《JOURNAL OF MULTIVARIATE ANALYSIS》 * |
| 夏毅锐等: "基于Copula函数的机械结构可靠性评估研究", 《自动化与仪器仪表》 * |
| 张世艳等: "惯性器件加速试验环境谱编制与数据处理方法研究", 《装备环境工程》 * |
| 王孟: "基于Copula函数的多元加速退化试验方法研究", 《中国优秀硕士学位论文全文数据库工程科技Ⅱ辑》 * |
| 王浩伟等: "随机环境应力冲击下基于多参数相关退化的导弹部件寿命预测", 《航空学报》 * |
| 马小兵等: "基于伪寿命分布的退化数据可靠性评估方法", 《系统工程与电子技术》 * |
| 鲍兆伟等: "基于Copula函数选择的多元退化产品可靠性模型", 《科学技术与工程》 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111160713A (en) * | 2019-12-06 | 2020-05-15 | 中国南方电网有限责任公司超高压输电公司广州局 | Composite insulator reliability assessment method based on multidimensional joint distribution theory |
| CN113361142A (en) * | 2021-07-12 | 2021-09-07 | 西安工程大学 | Spinning machine reliability evaluation method with multiple parameter degradation dependence |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109522519B (en) | 2023-02-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Le Son et al. | Remaining useful lifetime estimation and noisy gamma deterioration process | |
| Chwialkowski et al. | A kernel test of goodness of fit | |
| Liu et al. | Application of Bayesian networks for diagnostics in the assembly process by considering small measurement data sets | |
| Fang et al. | Inverse Gaussian processes with correlated random effects for multivariate degradation modeling | |
| Zhang et al. | Phase I analysis of multivariate profiles based on regression adjustment | |
| Beasley | Seemingly unrelated regression (SUR) models as a solution to path analytic models with correlated errors | |
| Han et al. | Control of the false discovery rate under arbitrary covariance dependence | |
| Bücher et al. | Multiplier bootstrap of tail copulas with applications | |
| CN109522519A (en) | A kind of dependence evaluation method between multiple performance parameters of ammunition parts | |
| Huang et al. | Two-mode network autoregressive model for large-scale networks | |
| Song et al. | Analysis of structural equation model with ignorable missing continuous and polytomous data | |
| ul Haq et al. | McDonald power function distribution with theory and applications | |
| Pang et al. | RUL prediction for bivariate degradation process considering individual differences | |
| Bibi et al. | Yule–Walker type estimators in periodic bilinear models: strong consistency and asymptotic normality | |
| Tripathi et al. | An approximate Bayes analysis of ARMA model for Indian GDP growth rate data | |
| Jiao et al. | A corrected and more efficient suite of MCMC samplers for the multinomal probit model | |
| Wang et al. | Generalized extreme value regression for ordinal response data | |
| Xiao et al. | Analysis of reliability experiments with blocking | |
| Xu et al. | A Stein Goodness of Fit Test for Exponential Random Graph Models | |
| Ma et al. | A variational Bayesian approach to modelling with random time-varying time delays | |
| Iliopoulos et al. | Bayesian estimation in Kibble's bivariate gamma distribution | |
| Zhang et al. | Estimation method for extremely small sample accelerated degradation test data | |
| Zhang et al. | Bayesian local influence analysis of general estimating equations with nonignorable missing data | |
| Franke et al. | Network modularity in the presence of covariates | |
| Romero | Relationships between ASME VV 10 & 20 AIAA CFD and Real Space Model Validation Frameworks. |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CP03 | Change of name, title or address |
Address after: 400039 Chongqing Jiulongpo Yuzhou Road No. 33 Patentee after: NO.59 Institute of China Ordnance Industry Country or region after: China Address before: 400039 Chongqing Jiulongpo Shiqiaopu Yuzhou Road No. 33 Patentee before: NO 59 Research Institute OF CHINA ORDNACE INDUSTRY Country or region before: China |