CN104571088B - Satellite control system Multipurpose Optimal Method based on fault diagnosability constraint - Google Patents
Satellite control system Multipurpose Optimal Method based on fault diagnosability constraint Download PDFInfo
- Publication number
- CN104571088B CN104571088B CN201410829123.6A CN201410829123A CN104571088B CN 104571088 B CN104571088 B CN 104571088B CN 201410829123 A CN201410829123 A CN 201410829123A CN 104571088 B CN104571088 B CN 104571088B
- Authority
- CN
- China
- Prior art keywords
- particle
- matrix
- control system
- row
- design index
- 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.)
- Active
Links
Classifications
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
- G05B23/00—Testing or monitoring of control systems or parts thereof
- G05B23/02—Electric testing or monitoring
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
- G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
- G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
- G05B13/0265—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric the criterion being a learning criterion
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
- G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
- G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
- G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
Abstract
Satellite control system multiple-objection optimization collocation method of the present invention based on fault diagnosability constraint, step are as follows:Build the even neighbouring matrix of satellite control system;Initialize particle swarm parameter;Generate each subassembly selection vector of satellite control system;Calculate and judge each component costs of satellite control system and measure whether constraint meets to require;Calculate measurement accuracy and diagnosticability Measure Indexes and particle comprehensive Design index;Judge whether comprehensive Design index is more than the optimal synthesis design index of current particle record, if the optimal record of more new particle more than if;Judge whether particle comprehensive Design index is more than the optimal synthesis design index of particle group records, if the optimal record that population is updated more than if;Judge the optimal synthesis design index of population, if meet that regulation requires;Update particle swarm parameter;The present invention provides the part allocation optimum for meeting that diagnosticability requires by the redundancy relationship between analyzing each part output, has filled up both at home and abroad in the blank of the technical field.
Description
Technical field
The present invention relates to a kind of satellite control system Multipurpose Optimal Method based on fault diagnosability constraint, belong to and defend
Star overall design technique field.
Background technology
At present, each part configuration of satellite control system more considers the function and performance of system, and is determined for basic
Allocation plan carry out method for diagnosing faults design, due to be supplied to the information content of fault diagnosis and quality be it is constant, because
The satellite control system trouble diagnosibility that this is realized in this way is very limited.In addition, even if part satellite model exists
Design phase have also contemplated that the demand of fault diagnosis, but simply by rule of thumb, come minor details be changed by analyzing failure one by one
Substantially the configurations of components determined, not yet proceed from the situation as a whole by analyze the redundancy relationship between the output of each part provide satisfaction can
The part allocation optimum of diagnostic requirement.
The content of the invention
The technology of the present invention solves problem:In view of the shortcomings of the prior art, overcome the deficiencies in the prior art, there is provided a kind of
Based on the satellite control system Multipurpose Optimal Method of fault diagnosability constraint, by the measurement essence for considering satellite control system
The indexs such as degree, performance of fault diagnosis and cost, provide each part distributes result rationally, ensure that satellite control system integrates
Best performance.
The present invention technical solution be:
It is as follows based on the satellite control system Multipurpose Optimal Method of fault diagnosability constraint, including step:
(1) the even neighbouring matrix of satellite control system is built;
(2) particle swarm parameter is initialized, the particle swarm parameter includes particle number, the Position And Velocity of particle, particle
With the optimal synthesis design index of population;
(3) each subassembly selection vector of satellite control system is generated according to particle position;
(4) the subassembly selection vector obtained according to step (3), calculates and whether judges each component costs of satellite control system
Less than defined cost and measurement constraint whether meet to require, if be less than and meets measurement constraint requirements if enter step (5),
Otherwise step (8) is entered;The measurement is constrained to each part of satellite control system and disclosure satisfy that normal measurement request;
(5) the even neighbouring matrix that the subassembly selection vector sum step (1) obtained according to step (3) obtains, calculates measurement essence
Degree and diagnosticability Measure Indexes, and obtain particle comprehensive Design index;
(6) whether the particle comprehensive Design index that judgment step (5) obtains is more than the optimal synthesis that current particle records and sets
Index is counted, if being more than, the optimal record of current particle is updated using the comprehensive Design index and its corresponding component selection vector,
The optimal record of described particle includes optimal synthesis design instruction and its corresponding subassembly selection vector, otherwise into step (8);
(7) whether the particle comprehensive Design index that judgment step (5) obtains is more than the optimal synthesis design of particle group records
Index, it is described using the optimal record of the comprehensive Design index and its corresponding component selection vector renewal population if being more than
The optimal record of population include optimal synthesis design instruction and its corresponding subassembly selection is vectorial, otherwise into step (8);
(8) if all particles of population all executed step (3)-(7) are (if bar of the particle due to being unsatisfactory for step (4)
Part and be not carried out step (5)-(7), be also considered as particle executed step (3)-(7)), then judge that the optimal synthesis of population is set
Whether meter index meets regulation requirement, if meeting to enter step (10), step (9) is entered if being unsatisfactory for;If in population
Still suffer from particle and be not carried out step (3)-(7), then therefrom select a particle and perform step (3)-(7);
(9) particle swarm parameter, repeat step (3)-(8) are updated;
(10) terminate.
The concrete mode of even neighbouring matrix structure is as follows in step (1):
(1a) establishes following mould according to sensor measurement model, executing agency's model, satellite motion and kinetics equation
Type:
Wherein, CiNumbered for equation, i=1,2 ..., ns+na+ 6, nsRepresent sensor number, naRepresent executing agency
Number, m is total number of parts m=ns+na,The measurement model of sensor is represented,Table
Show executing agency's model,Represent satellite motion and kinetics equation, GiFor letter
Number expression formula (for example, for different sensor models, function expression is different, but is known in the art);KμTo be each
Sensor exports, μ=1,2 ..., ns, xτFor the known variables of each model, τ=1,2 ..., na+ 6, fεFor the failure letter of each part
Number, ε=1,2 ..., m;
(1b), using known variables as row, builds even neighbouring matrix H ' (i, j) using the equation in step (1a) model as row,
If known variables xjIt is present in equation CiIn, then claim known variables xjWith equation CiBetween side S (C be presenti,xj), juxtaposition element H '
(i, j) is 1, is otherwise 0, wherein j=1 ..., na+6。
Step (2) population parameter initialization includes:
(2a) particle number initializes basic principle:Particle is more, and search space is bigger, and searching process is shorter, can basis
Actual demand is specifically set;
The Position And Velocity of (2b) particle:It is arranged to 0 to 1 random number, dimension m;
The optimal synthesis design index of (2c) particle and population:Initial value is arranged to 0.
The subassembly selection vector L generated in step (3) according to particle positiondFor:
Ld=[Ld,1,Ld,2,...,Ld,z,...,Ld,m]T
Wherein, Ld,zRepresent subassembly selection vector Ld z dimension datas, z=1 ..., m, d=1 ..., Nparticle,
NparticleFor particle number,vd,zRepresent the z dimension datas of particle d speed, rand ()
Represent the number randomly generated.
The component costs of step (4) and measurement are constrained to:
(4a) is according to following formula calculating unit cost Cd,component:
Cd,component=CLd
Wherein, Cd,componentRepresent component costs corresponding to particle d associated components selection vector;C=[R1,...,Rm],
[R1,...,Rm] be each part cost;
(4b) it is described measurement constraints be:Πd,e≠ 0 (e=1,2 ... 6);
Wherein,Πd,e(b) representing matrix Πd,eB rows, ΦeFor the appearance of a rows m row
State determines scheme matrix, and wherein a represents to realize that X-axis, Y-axis or Z axis angle or angular speed determine the number of scheme, ΦeMatrix
Element be 0 or 1, e=1,2 ... 6 correspond to respectively X-axis angle-determining scheme, Y-axis angle-determining scheme, Z axis angle-determining scheme,
X-axis angular speed determines that scheme, Y-axis angular speed determine that scheme, Z axis angular speed determine scheme, b=1,2 ... a.
Measurement accuracy and diagnosticability index of correlation are calculated in step (5), and obtains the concrete mode of comprehensive Design index
It is as follows:
(5a) calculates measurement accuracy according to following formula:
Wherein:Pd,e=[Q1,...,Qa], [Q1,...,Qa] represent that posture determines the measurement accuracy of scheme, describe various
Accuracy corresponding to measurement scheme, it can be obtained according to expertise or mathematical simulation, [Q1,...,Qa]∈(0,1);
(5b) calculates diagnosticability Measure Indexes;Described diagnosticability Measure Indexes include detectable rate and separated
Rate;
(5b1) builds analytical redundancy relation;
(5b1a) is according to subassembly selection vector Ld, antithesis is adjacent to matrix H ' modify to obtain H:If in subassembly selection vector
Ld,z=0, then even neighbouring matrix H ' z rows delete;
(5b1b) obtains the maximum set of matches of even neighbouring matrix H using DM decomposition techniques;Maximum set of matches is known variables
With equation CiThe set on side;If M is the subset of side collection in figure Ξ, if any both sides all do not have same vertices in M, M is claimed to be
A Ξ matching;If in Ξ is schemed be not present another matching N, make in N while number be more than M in while number, then M is referred to as most
Big set of matches.
(5b1c) adds 1 column element into even neighbouring matrix H, judges side S (Ci,xj) whether belong to maximum set of matches, if category
In then element H (i, j) is changed into -1 from 1 corresponding to the side, while put element H (i, na+ 7) it is 1, otherwise by H (i, na+ 7) it is set to 0;
(5b1d) for matrix H the i-th row, if known variables xjCorresponding element H (i, j) is -1, and removes step
Other elements beyond 1 row of (5b1c) addition are 0, then xjExpression formula can be according to corresponding equation CiAnd in step (1)
Mathematical modeling obtainxjExpression formula contain known quantity and failure,Table
Show function expression GiInverse transformation;For the i-th row of matrix H, if known variables xjCorresponding element H (i, j) is -1, is deposited simultaneously
In H (i, w)=1, and xwExpression formula pass through q row equatioies CqObtain
Then by xwExpression formula substitute into equation CiIn obtain xjExpression formula isWherein w=
1,2,...,na+ 6, and w ≠ j, q=1,2 ..., ns+na+ 6 and q ≠ i, aforesaid operations are repeated, until all known variables obtain
Expression formula;
Last column element of the i-th row of (5b1e) matrix H is 0, and the row nonzero element corresponds to known variables and all obtained
Expression formula, then substitute into corresponding expression formula in equation corresponding to the row that (row corresponds to equation can be according to the mathematical modeling of step (1)
Obtain), obtained equation is analytical redundancy relation;The row for being 0 to last all column element performs similar operations, obtains institute
Some analytical redundancy relations;
(5b1f) carried out any one edge contract in maximum set of matches, repetition (5b1b)~(5b1e) with last time
The analytical redundancy relation set that (5b1b)~(5b1e) is obtained merges, and obtains new analytical redundancy relation set, until inciting somebody to action
Each edge in maximum set of matches completes analysis, and above-mentioned steps are not repeated, obtain final analytical redundancy relation set;
(5b2) builds incidence matrix;
Using the analytical redundancy relation obtained in step (5b1) as row, using unit failure as row, structure incidence matrix Ω:
Wherein, Ω represents NARRThe matrix of row m row, Ω elements are 0 or 1;NARRRepresent the number of analytical redundancy relation;U=
1 ..., m, y=1 ..., NARR;fuPart u failure function is represented, l represents y-th of analytical redundancy relation;
(5b3) diagnosability analysis;
In incidence matrix Ω, if unit failure function fuThe all elements of respective column are all 0, then the unit failure can not
Detection, otherwise can detect;In incidence matrix Ω, if the incomplete phase of element of all unit failure respective columns of satellite control system
Together, then each part has separability, if the identical unit failure of the element that respective column in all unit failures be present, claims
The collection of these unit failures composition is combined into ambiguity group;
(5b4) calculates detectable rate FDR:
Wherein, DuDegree, when unit failure can detect, D are can detect for part u failureu=1, otherwise Du=0, λuFor portion
Part u fault rate, m are total number of parts;
(5b5) calculates separable rate FIR:
Wherein, IuFor part u failure degree of isolation, when unit failure separates, Iu=1, otherwise Iu=0, IuFor portion
Part u fault rate, m are total number of parts;
The measurement accuracy of acquisition, detectable rate and separable rate are obtained comprehensive Design index by (5c) by weighted average.
The concrete mode of renewal particle swarm parameter is as follows in step (5):
Particle d Position And Velocity is updated using following formula:
vd(k+1)=ω vd(k)+c1r1(pd(k)-sd(k))+c2r2(pg(k)-sd(k))
sd(k+1)=sd(k)+vi(k+1)
Wherein, vdRepresent particle d speed, sdRepresent particle d position, pdRepresent particle d optimal synthesis design index
Corresponding particle position, pgRepresent particle position corresponding to the optimal synthesis design index of population, ω, c1、r1、r2Represent with
Machine number.
The present invention has the following advantages that compared with prior art:
(1) present invention is using the diagnosable rate and measurement accuracy of satellite control system as optimization aim, using component costs as about
Beam condition, the Optimal Configuration Method of satellite control system is provided, the design phase will be advanced to the consideration of fault diagnosis, be failure
Diagnostic method research provides advantage, the present invention proceed from the situation as a whole by analyze the redundancy relationship between the output of each part and to
Go out to meet the part allocation optimum of diagnosticability requirement, filled up both at home and abroad in the blank of the technical field.
(2) the inventive method carries out the diagnosability analysis of satellite control system, implementation process based on maximum matching technique
It is simple and clear, it is easy to accomplish, versatility is stronger.
(3) present invention optimizes solution using particle cluster algorithm, accelerates search procedure by intelligent search technique, carries
The high efficiency of fault diagnosis and satellite operation, greatlys save manpower and hardware cost.
Brief description of the drawings
Fig. 1 is flow chart of the present invention.
Embodiment
Just the present invention is described further with reference to accompanying drawing below.
As shown in figure 1, the satellite control system Multipurpose Optimal Method based on fault diagnosability constraint, including step is such as
Under:
(1) the even neighbouring matrix of satellite control system is built;
The concrete mode of even neighbouring matrix structure is as follows:
(1a) according to sensor measurement model (present invention used in model refer to sensor model general on satellite, such as
The sensors such as gyro, earth sensor), executing agency's model, satellite motion and kinetics equation, establish such as drag:
(1a) establishes following mould according to sensor measurement model, executing agency's model, satellite motion and kinetics equation
Type:
Wherein, CiNumbered for equation, i=1,2 ..., ns+na+ 6, nsRepresent sensor number, naRepresent executing agency
Number, m is total number of parts m=ns+na,The measurement model of sensor is represented,Table
Show executing agency's model,Represent satellite motion and kinetics equation, GiFor letter
Number expression formula (for example, for different sensor models, function expression is different, but is known in the art);KμTo be each
Sensor exports, μ=1,2 ..., ns, xτFor the known variables of each model, τ=1,2 ..., na+ 6, fεFor the failure letter of each part
Number, ε=1,2 ..., m;
(1b), using known variables as row, builds even neighbouring matrix H ' (i, j) using the equation in step (1a) model as row,
If known variables xjIt is present in equation CiIn, then claim known variables xjWith equation CiBetween side S (C be presenti,xj), juxtaposition element H '
(i, j) is 1, is otherwise 0, wherein j=1 ..., na+6。
Even neighbouring matrix is really a kind of description form of bigraph, if the vertex set of figure is divided into two nonvoid subsets X and Y,
And each edge has a summit in X, and another summit is in Y, then this figure is referred to as bigraph.Two summits difference of bigraph
The corresponding occasionally existing side between the row and column of matrix, non-zero two vertex sets of element representation.For satellite control system mould
For type, equation and variable are the vertex set of bigraph respectively, if variable is present in the equation, then it is assumed that variable and equation it
Between a line be present, therefore even neighbouring matrix relevant position is 1.
Illustrate the method for building up of even neighbouring matrix with a specific embodiment below:For example, it is configured with three orthogonal gyros, rolling
The satellite control system of moving axis infrared earth sensor and three orthogonal momenttum wheels, the measurement model of comprehensive sensor, executing agency
The kinematics and kinetics equation of model and satellite, establish such as drag:
C1:g1=ωx+f1
C2:g2=ωy+f2
C3:g3=ωz+f3
C5:Tx=ux+f5
C6:Ty=uy+f6
C7:Tz=uz+f7
Wherein, gλExported for gyro, λ=1,2,3, ωπFor the axis angular rate of satellite three, π=x, y, z,For the quick sensor in ground
Output,For three shaft angle degree, uπFor three axle desired output torques, ω0For orbit angular velocity, IπRotated for each axle of satellite used
Amount, TπFor momenttum wheel actual output torque,Represent variable ξ differential form.
It is n corresponding to formula (1) each parameters=4, na=3, K1=g1, K2=g2, K3=g3,
The even neighbouring matrix H obtained according to step (1b) ' as shown in table 1.
The even neighbouring matrix H of table 1 '
For example, parameter in table 1With equation C4Matrix element for summit is 1, then explanation has side.
(2) particle swarm parameter is initialized;The particle swarm parameter includes particle number, the Position And Velocity of particle, particle
With the optimal synthesis design index of population;
Population parameter initialization includes:
(2a) particle number initializes basic principle:Particle is more, and search space is bigger, and searching process is shorter, can basis
Actual demand is specifically set;
The Position And Velocity of (2b) particle:It is arranged to 0 to 1 random number, dimension m;
The optimal synthesis design index of (2c) particle and population:Initial value is arranged to 0.
(3) each subassembly selection vector of satellite control system is generated according to particle position;
Subassembly selection vector LdFor:
Ld=[Ld,1,Ld,2,...,Ld,z,...,Ld,m]T
Wherein, Ld,zRepresent subassembly selection vector LdZ dimension datas, z=1 ..., m, d=1 ..., Nparticle,
NparticleFor particle number,vd,zRepresent the z dimension datas of particle d speed, rand ()
Represent the number randomly generated.
For example, the satellite control system that step (1) citing is mentioned relates generally to gyro 1, gyro 2, gyro 3, infrared earth
Sensor, momenttum wheel 1, momenttum wheel 2, momenttum wheel 3, totally 7 parts, if the subassembly selection vector for calculating to obtain according to particle rapidity
For Ld=[1 11011 0]T, then it represents that gyro 1, gyro 2, gyro 3, momenttum wheel 1 and momenttum wheel are only selected during optimization design
2。
(4) the subassembly selection vector obtained according to step (3), calculates and whether judges each component costs of satellite control system
Less than defined cost and measurement constraint whether meet to require, if be less than and meets measurement constraint requirements if enter step (5),
Otherwise step (8) is entered;The measurement is constrained to each part of satellite control system and disclosure satisfy that normal measurement request;
Component costs and measurement are constrained to:
(4a) is according to following formula calculating unit cost Cd,component:
Cd,component=CLd
Wherein, Cd,componentRepresent component costs corresponding to particle d associated components selection vector;C=[R1,...,Rm],
[R1,...,Rm] be each part cost;
(4b) it is described measurement constraints be:Πd,e≠ 0 (e=1,2 ... 6);
Wherein,Πd,e(b) representing matrix Πd,eB rows, ΦeFor the appearance of a rows m row
State determines scheme matrix, and wherein a represents to realize that X-axis, Y-axis or Z axis angle or angular speed determine the number of scheme, ΦeMatrix
Element be 0 or 1, e=1,2 ... 6 correspond to respectively X-axis angle-determining scheme, Y-axis angle-determining scheme, Z axis angle-determining scheme,
X-axis angular speed determines that scheme, Y-axis angular speed determine that scheme, Z axis angular speed determine scheme, b=1,2 ... a.
For example, illustrated exemplified by satellite control system and X-axis angle-determining that step (1) citing is mentioned.Three kinds of schemes
The determination of X-axis angle can be realized, the first scheme is that the output of gyro 1 carries out integration acquisition, and second scheme is infraredly quick
The output of sensor is X-axis angle, and the third scheme is that X-axis angle speed is calculated according to the rotating speed of momenttum wheel 1 and kinetics equation
Degree, X-axis angle, Φ are obtained by integration1It is represented by:For the portion obtained in step (2)
Part selects vectorial Ld=[1 11011 0]T, it is known that Πd,1(1)=1, Πd,1(2)=0, Πd,1(3)=1, it is possible to recognize
To meet that X-axis angular surveying constrains.
(5) the even neighbouring matrix that the subassembly selection vector sum step (1) obtained according to step (3) obtains, calculates measurement essence
Degree and diagnosticability Measure Indexes, and obtain particle comprehensive Design index;
Measurement accuracy and diagnosticability index of correlation are calculated, and the concrete mode for obtaining comprehensive Design index is as follows:
(5a) calculates measurement accuracy according to following formula:
Wherein:Pd,e=[Q1,...,Qa], [Q1,...,Qa] represent that posture determines the measurement accuracy of scheme, describe various
Accuracy corresponding to measurement scheme, it can be obtained according to expertise or mathematical simulation, [Q1,...,Qa]∈(0,1);
For example, the precision for three kinds of X-axis angle-determining schemes that step (4) citing is mentioned is respectively 0.7,0.8 and 0.5, because
This obtains Pd,e=[0.7 0.8 0.5], the Π obtained due to step (4)d,1=[1 0 1]T, P can be obtainedd,1Πd,1=1.2.
(5b) calculates diagnosticability Measure Indexes;Described diagnosticability Measure Indexes include detectable rate and separated
Rate;
(5b) calculates diagnosticability Measure Indexes;Described diagnosticability Measure Indexes include detectable rate and separated
Rate;
(5b1) builds analytical redundancy relation;
(5b1a) is according to subassembly selection vector Ld, antithesis is adjacent to matrix H ' modify to obtain H:If in subassembly selection vector
Ld,z=0, then even neighbouring matrix H ' z rows delete;
(5b1b) obtains the maximum set of matches of even neighbouring matrix H using DM decomposition techniques;Maximum set of matches is known variables
With equation CiThe set on side;If M is the subset of side collection in figure Ξ, if any both sides all do not have same vertices in M, M is claimed to be
A Ξ matching;If in Ξ is schemed be not present another matching N, make in N while number be more than M in while number, then M is referred to as most
Big set of matches.
(5b1c) adds 1 column element into even neighbouring matrix H, judges side S (Ci,xj) whether belong to maximum set of matches, if category
In then element H (i, j) is changed into -1 from 1 corresponding to the side, while put element H (i, na+ 7) it is 1, otherwise by H (i, na+ 7) it is set to 0;
(5b1d) for matrix H the i-th row, if known variables xjCorresponding element H (i, j) is -1, and removes step
Other elements beyond 1 row of (5b1c) addition are 0, then xjExpression formula can be according to corresponding equation CiAnd in step (1)
Mathematical modeling obtainxjExpression formula contain known quantity and failure,Table
Show function expression GiInverse transformation;For the i-th row of matrix H, if known variables xjCorresponding element H (i, j) is -1, is deposited simultaneously
In H (i, w)=1, and xwExpression formula pass through q row equatioies CqObtain
Then by xwExpression formula substitute into equation CiIn obtain xjExpression formula isWherein w=
1,2,...,na+ 6, and w ≠ j, q=1,2 ..., ns+na+ 6 and q ≠ i, aforesaid operations are repeated, until all known variables obtain
Expression formula;
Last column element of the i-th row of (5b1e) matrix H is 0, and the row nonzero element corresponds to known variables and all obtained
Expression formula, then substitute into corresponding expression formula in equation corresponding to the row that (row corresponds to equation can be according to the mathematical modeling of step (1)
Obtain), obtained equation is analytical redundancy relation;The row for being 0 to last all column element performs similar operations, obtains institute
Some analytical redundancy relations;
(5b1f) carried out any one edge contract in maximum set of matches, repetition (5b1b)~(5b1e) with last time
The analytical redundancy relation set that (5b1b)~(5b1e) is obtained merges, and obtains new analytical redundancy relation set, until inciting somebody to action
Each edge in maximum set of matches completes analysis, and above-mentioned steps are not repeated, obtain final analytical redundancy relation set;
For example, the subassembly selection vector obtained for the even neighbouring matrix H ' and step (3) shown in table 1 is Ld=[1 11
0 1 1 0]T, then antithesis adjacent to matrix H ' the 4th row and the 7th row delete after obtain H, carry out maximum obtained after DM decomposition
It is with collection:{C9,θ},{C10,ψ},{C13,Tz, according to step
Suddenly (5b1c) obtains even neighbouring matrix shown in table 2 after being identified.
Table 2 be identified after even neighbouring matrix H
For the 1st row in table 2, in addition to being classified as 1 except last, only ωxCorresponding element is -1, so ωx's
Expression formula can be according to C1:g1=ωx+f1Obtain, as ωx=g1-f1, similarly obtain ωy=g2-f2.And for the 9th in table 2
OK, except last is classified as 1, element corresponding to variable θ is outside -1, and element corresponding to variable θ is 1, while ωyAccording to C2Obtain
Expression formula be ωy=g2-f2, so according toθ=∫ (ω can be derived0+g2-f2)dt.Utilize similar step
Suddenly variable ω is obtainedzAnd TxExpression formula after, by ωx, ωy, ωzAnd TxExpression formula substitute into equation C11In an i.e. available solution
Analyse redundancy relationship.
(5b2) builds incidence matrix;
Using the analytical redundancy relation obtained in step (5b1) as row, using unit failure as row, structure incidence matrix Ω:
Wherein, Ω represents NARRThe matrix of row m row, Ω elements are 0 or 1;NARRRepresent the number of analytical redundancy relation;U=
1 ..., m, y=1 ..., NARR;fuPart u failure function is represented, l represents y-th of analytical redundancy relation;
(5b3) diagnosability analysis;
Diagnosticability refers to the degree that the failure in system can be identified accurately and effectively, wherein accurately referring in event
Barrier can unambiguously be detected every time when occurring and separation failure, and effectively refers to that resource needed for fault reconstruction can be carried out
Optimization, fault diagnosability include detectability and separability.
Failure has detectability, in a period of time after failure that and if only if generation, even if disturbance be present, can also be used
System inputs the detection that failure is realized with output information.During detectability analysis, if the output of fault impact being capable of structure
Analytical redundancy relation is built, then the failure has detectability.The isolabilily refers to for the different faults in given set, is
Ability of the system with different manifestations.During separability analysis, if in the failure collection of analysis, different output is built not
Same analytical redundancy relation, then institute is faulty has separability.
In incidence matrix Ω, if unit failure function fuThe all elements of respective column are all 0, then the unit failure can not
Detection, otherwise can detect;In incidence matrix Ω, if the incomplete phase of element of all unit failure respective columns of satellite control system
Together, then each part has separability, if the identical unit failure of the element that respective column in all unit failures be present, claims
The collection of these unit failures composition is combined into ambiguity group;
(5b4) calculates detectable rate FDR:
Wherein, DuDegree, when unit failure can detect, D are can detect for part u failureu=1, otherwise Du=0, λuFor portion
Part u fault rate, m are total number of parts;
(5b5) calculates separable rate FIR:
Wherein, IuFor part u failure degree of isolation, when unit failure separates, Iu=1, otherwise Iu=0, IuFor portion
Part u fault rate, m are total number of parts;
The measurement accuracy of acquisition, detectable rate and separable rate are obtained comprehensive Design index by (5c) by weighted average.
(6) whether the particle comprehensive Design index that judgment step (5) obtains is more than the optimal synthesis that current particle records and sets
Index is counted, if being more than, the optimal record of current particle is updated using the comprehensive Design index and its corresponding component selection vector,
The optimal record of described particle includes optimal synthesis design instruction and its corresponding subassembly selection vector, otherwise into step (8);
(7) whether the particle comprehensive Design index that judgment step (5) obtains is more than the optimal synthesis design of particle group records
Index, it is described using the optimal record of the comprehensive Design index and its corresponding component selection vector renewal population if being more than
The optimal record of population include optimal synthesis design instruction and its corresponding subassembly selection is vectorial, otherwise into step (8);
(8) if all particles of population all executed step (3)-(7) are (if bar of the particle due to being unsatisfactory for step (4)
Part and be not carried out step (5)-(7), be also considered as particle executed step (3)-(7)), then judge that the optimal synthesis of population is set
Whether meter index meets regulation requirement, if meeting to enter step (10), step (9) is entered if being unsatisfactory for;If in population
Still suffer from particle and be not carried out step (3)-(7), then therefrom select a particle and perform step (3)-(7);
(9) particle swarm parameter, repeat step (3)-(8) are updated;
The concrete mode for updating particle swarm parameter is as follows:
Particle d Position And Velocity is updated using following formula:
vd(k+1)=ω vd(k)+c1r1(pd(k)-sd(k))+c2r2(pg(k)-sd(k))
sd(k+1)=sd(k)+vi(k+1)
Wherein, vdRepresent particle d speed, sdRepresent particle d position, pdRepresent particle d optimal synthesis design index
Corresponding particle position, pgRepresent particle position corresponding to the optimal synthesis design index of population, ω, c1、r1、r2Represent with
Machine number.
(10) terminate.
If the method described in the present invention is used in the case where considering the current highest configuration of each part, the side so reached
Method effect is optimum efficiency.
Unspecified part of the present invention belongs to general knowledge as well known to those skilled in the art.
Claims (7)
1. the satellite control system Multipurpose Optimal Method based on fault diagnosability constraint, it is characterised in that step is as follows:
(1) the even neighbouring matrix of satellite control system is built;
(2) particle swarm parameter is initialized, the particle swarm parameter includes particle number, the Position And Velocity of particle, particle and grain
The optimal synthesis design index of subgroup;
(3) each subassembly selection vector of satellite control system is generated according to particle rapidity;
(4) the subassembly selection vector obtained according to step (3), calculates and judges whether each component costs of satellite control system are less than
Defined cost and measurement constraint whether meets to require, if be less than and meets measurement constraint requirements if enter step (5), otherwise
Into step (8);The measurement is constrained to each part of satellite control system and disclosure satisfy that normal measurement request;
(5) the even neighbouring matrix that the subassembly selection vector sum step (1) obtained according to step (3) obtains, calculate measurement accuracy with
Diagnosticability Measure Indexes, and obtain particle comprehensive Design index;
(6) whether the particle comprehensive Design index that judgment step (5) obtains is more than the optimal synthesis design that current particle records and refers to
Mark, it is described using the comprehensive Design index and its optimal record of corresponding component selection vector renewal current particle if being more than
The optimal record of particle to include optimal synthesis design index and its corresponding subassembly selection vectorial, otherwise into step (8);
(7) whether the particle comprehensive Design index that judgment step (5) obtains is more than the optimal synthesis design index of particle group records,
If being more than, the optimal record of the comprehensive Design index and its corresponding component selection vector renewal population, described grain are utilized
The optimal record in subgroup includes optimal synthesis design index and its corresponding subassembly selection vector, otherwise into step (8);
(8) if all particles of population all executed step (3)-(7), judging the optimal synthesis design index of population is
It is no to meet that regulation requires, if meeting to enter step (10), step (9) is entered if being unsatisfactory for;If grain is still suffered from population
Son is not carried out step (3)-(7), then therefrom selects a particle and perform step (3)-(7);Wherein, when particle is due to being unsatisfactory for
The condition of step (4) and when being not carried out step (5)-(7), it is believed that particle executed step (3)-(7);
(9) particle swarm parameter, repeat step (3)-(8) are updated;
(10) terminate.
2. the satellite control system Multipurpose Optimal Method as claimed in claim 1 based on fault diagnosability constraint, it is special
Sign is:The concrete mode of even neighbouring matrix structure is as follows in the step (1):
(1a) establishes such as drag according to sensor measurement model, executing agency's model, satellite motion and kinetics equation:
Wherein, CiNumbered for equation, i=1,2 ..., ns+na+ 6, nsRepresent sensor number, naExecuting agency's number is represented, m is
Total number of parts m=ns+na,The measurement model of sensor is represented,Represent to perform
Mechanism model,Represent satellite motion and kinetics equation, GiFor function representation
Formula;KμExported for each sensor, μ=1,2 ..., ns, xτFor the known variables of each model, τ=1,2 ..., na+ 6, fεFor each part
Failure function, ε=1,2 ..., m;
(1b), using known variables as row, builds even neighbouring matrix H ' (i, j), if not using the equation in step (1a) model as row
Know variable xjIt is present in equation CiIn, then claim known variables xjWith equation CiBetween side S (C be presenti,xj), juxtaposition H ' (i, j) is 1,
Otherwise it is 0, wherein j=1 ..., na+6。
3. the satellite control system Multipurpose Optimal Method as claimed in claim 1 based on fault diagnosability constraint, it is special
Sign is:Step (2) the population parameter initialization includes:
(2a) particle number initializes basic principle:Particle is more, and search space is bigger, and searching process is shorter, can be according to reality
Border demand is specifically set;
The Position And Velocity of (2b) particle:It is arranged to 0 to 1 random number, dimension m;
The optimal synthesis design index of (2c) particle and population:Initial value is arranged to 0.
4. the satellite control system Multipurpose Optimal Method as claimed in claim 3 based on fault diagnosability constraint, it is special
Sign is:The subassembly selection vector L generated in the step (3) according to particle rapiditydFor:
Ld=[Ld,1,Ld,2,...,Ld,z,...,Ld,m]T
Wherein, Ld,zRepresent subassembly selection vector LdZ dimension datas, z=1 ..., m, d=1 ..., Nparticle, NparticleFor
Particle number,vd,zThe z dimension datas of particle d speed are represented, rand () represents random production
Raw number.
5. the satellite control system Multipurpose Optimal Method as claimed in claim 4 based on fault diagnosability constraint, it is special
Sign is:The component costs of the step (4) and measurement are constrained to:
(4a) is according to following formula calculating unit cost Cd,component:
Cd,component=CLd
Wherein, Cd,componentRepresent component costs corresponding to particle d associated components selection vector;C=[R1,...,Rm],
[R1,...,Rm] be each part cost;
(4b) described measurement is constrained to:Πd,e≠ 0 (e=1,2 ... 6);
Wherein,Πd,e(b) representing matrix Πd,eB rows, ΦePosture for a rows m row determines
Scheme matrix, wherein a represent to realize that X-axis, Y-axis or Z axis angle or angular speed determine the number of scheme, ΦeMatrix element be 0
Or 1, e=1,2 ... 6 correspond to X-axis angle-determining scheme, Y-axis angle-determining scheme, Z axis angle-determining scheme, X-axis angle speed respectively
Degree determination scheme, Y-axis angular speed determine that scheme, Z axis angular speed determine scheme, b=1,2 ... a.
6. the satellite control system Multipurpose Optimal Method as claimed in claim 5 based on fault diagnosability constraint, it is special
Sign is:Measurement accuracy and diagnosticability Measure Indexes are calculated in the step (5), and obtain the specific of comprehensive Design index
Mode is as follows:
(5a) calculates measurement accuracy according to following formula:
Wherein:Pd,e=[Q1,...,Qa], QμRepresent that posture determines the measurement accuracy of scheme, it is corresponding to describe various measurement schemes
Accuracy, can be obtained according to expertise or mathematical simulation, [Q1,...,Qa]∈(0,1);
(5b) calculates diagnosticability Measure Indexes;Described diagnosticability Measure Indexes include detectable rate and separable rate;
(5b1) builds analytical redundancy relation;
(5b1a) is according to subassembly selection vector Ld, antithesis is adjacent to matrix H ' modify to obtain H:If L in subassembly selection vectord,z=
0, then even neighbouring matrix H ' z rows delete;
(5b1b) obtains the maximum set of matches of even neighbouring matrix H using DM decomposition techniques;Maximum set of matches be known variables with etc.
Formula CiThe set on side;If M is the subset of side collection in figure Ξ, if any both sides all do not have same vertices in M, it is Ξ to claim M
One matching;If in Ξ is schemed be not present another matching N, make in N while number be more than M in while number, then M be referred to as maximum
Set of matches.
(5b1c) adds 1 column element into even neighbouring matrix H, judges side S (Ci,xj) whether belong to maximum set of matches, if belonging to
Element H (i, j) is changed into -1 from 1 corresponding to the side, while puts element H (i, na+ 7) it is 1, otherwise by H (i, na+ 7) it is set to 0;
(5b1d) for matrix H the i-th row, if known variables xjCorresponding element H (i, j) is -1, and except step (5b1c) is added
1 row beyond other elements be 0, then xjExpression formula can be according to corresponding equation CiAnd the mathematical modeling in step (1)
ObtainxjExpression formula contain known quantity and failure,Representative function table
Up to formula GiInverse transformation;For the i-th row of matrix H, if known variables xjCorresponding element H (i, j) is -1, while H (i, w) be present
=1, and xwExpression formula pass through q row equatioies CqObtainThen will
xwExpression formula substitute into equation CiIn obtain xjExpression formula isWherein w=1,
2,...,na+ 6, and w ≠ j, q=1,2 ..., ns+na+ 6 and q ≠ i, aforesaid operations are repeated, until all known variables obtain table
Up to formula;
Last column element of the i-th row of (5b1e) matrix H is 0, and the row nonzero element corresponds to known variables and all expressed
Formula, then corresponding expression formula is substituted into equation corresponding to the row, obtained equation is analytical redundancy relation;To it is all last
The row that column element is 0 performs similar operations, obtains all analytical redundancy relations;
(5b1f) carried out any one edge contract in maximum set of matches, repetition (5b1b)~(5b1e) with last time
The analytical redundancy relation set that (5b1b)~(5b1e) is obtained merges, and obtains new analytical redundancy relation set, until inciting somebody to action
Each edge in maximum set of matches completes analysis, and above-mentioned steps are not repeated, obtain final analytical redundancy relation set;
(5b2) builds incidence matrix;
Using the analytical redundancy relation obtained in step (5b1) as row, using unit failure as row, structure incidence matrix Ω:
Wherein, Ω represents NARRThe matrix of row m row, Ω elements are 0 or 1;NARRRepresent the number of analytical redundancy relation;U=
1 ..., m, y=1 ..., NARR;fuPart u failure function is represented,Represent y-th of analytical redundancy relation;
(5b3) diagnosability analysis;
In incidence matrix Ω, if unit failure function fuThe all elements of respective column are all 0, then the unit failure is undetectable,
Otherwise can detect;In incidence matrix Ω, if the element of all unit failure respective columns of satellite control system is incomplete same,
Each part has separability, if the element that respective column in all unit failures be present is identical, claims these unit failures
The collection of composition is combined into ambiguity group;
(5b4) calculates detectable rate FDR:
Wherein, DuDegree, when unit failure can detect, D are can detect for part u failureu=1, otherwise Du=0, λuFor part u's
Fault rate, m are total number of parts;
(5b5) calculates separable rate FIR:
Wherein, IuFor part u failure degree of isolation, when unit failure separates, Iu=1, otherwise Iu=0, m are that part is total
Number;
The measurement accuracy of acquisition, detectable rate and separable rate are obtained comprehensive Design index by (5c) by weighted average.
7. the satellite control system Multipurpose Optimal Method as claimed in claim 1 based on fault diagnosability constraint, it is special
Sign is:The concrete mode of renewal particle swarm parameter is as follows in the step (9):
Particle d Position And Velocity is updated using following formula:
vd(k+1)=ω vd(k)+c1r1(pd(k)-sd(k))+c2r2(pg(k)-sd(k))
sd(k+1)=sd(k)+vd(k+1)
Wherein, vdRepresent particle d speed, sdRepresent particle d position, pdRepresent that particle d optimal synthesis design index is corresponding
Particle position, pgRepresent particle position corresponding to the optimal synthesis design index of population, ω, c1、c2、r1、r2Represent random
Number.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410829123.6A CN104571088B (en) | 2014-12-26 | 2014-12-26 | Satellite control system Multipurpose Optimal Method based on fault diagnosability constraint |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410829123.6A CN104571088B (en) | 2014-12-26 | 2014-12-26 | Satellite control system Multipurpose Optimal Method based on fault diagnosability constraint |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104571088A CN104571088A (en) | 2015-04-29 |
CN104571088B true CN104571088B (en) | 2018-01-05 |
Family
ID=53087414
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201410829123.6A Active CN104571088B (en) | 2014-12-26 | 2014-12-26 | Satellite control system Multipurpose Optimal Method based on fault diagnosability constraint |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104571088B (en) |
Families Citing this family (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105678423B (en) * | 2016-01-28 | 2019-09-24 | 西北工业大学 | Fault diagnosis system Optimum sensor placement method based on D-M (Determiner-Measure) construction model |
CN106647270B (en) * | 2016-12-21 | 2019-07-12 | 北京控制工程研究所 | For the STABLE ADAPTIVE FUZZY Vibration Active Control method of the close frequency structure in space |
CN111324036B (en) * | 2020-01-19 | 2020-11-20 | 北京空间飞行器总体设计部 | Diagnosability quantification method for time-varying system under influence of bounded interference |
CN111290366B (en) * | 2020-02-12 | 2022-05-27 | 北京科技大学顺德研究生院 | Multi-fault diagnosis method for attitude control system of spacecraft |
CN115616919B (en) * | 2022-11-15 | 2023-03-17 | 中国航空工业集团公司金城南京机电液压工程研究中心 | Electromechanical product sensor optimal configuration method |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6341249B1 (en) * | 1999-02-11 | 2002-01-22 | Guang Qian Xing | Autonomous unified on-board orbit and attitude control system for satellites |
CN101859146A (en) * | 2010-07-16 | 2010-10-13 | 哈尔滨工业大学 | Satellite fault prediction method based on predictive filtering and empirical mode decomposition |
CN102735877A (en) * | 2012-06-18 | 2012-10-17 | 北京控制工程研究所 | Optimization method for measuring point of measuring sensor for satellite angular rate based on correlation matrix |
CN102736616A (en) * | 2012-06-18 | 2012-10-17 | 北京控制工程研究所 | Dulmage-Mendelsohn (DM)-decomposition-based measuring point optimal configuration method for closed loop system |
CN103116357A (en) * | 2013-03-14 | 2013-05-22 | 郭雷 | Sliding-mode control method with anti-interference fault-tolerance performance |
CN103676941A (en) * | 2013-12-24 | 2014-03-26 | 北京控制工程研究所 | Satellite control system fault diagnosis method based on kinematics and dynamics model |
CN103697915A (en) * | 2013-12-24 | 2014-04-02 | 北京控制工程研究所 | Diagnostic evaluation method considering disturbing influence for satellite sensor fault |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7458264B2 (en) * | 2004-09-10 | 2008-12-02 | Honeywell International Inc. | Generalized inertial measurement error reduction through multiple axis rotation during flight |
-
2014
- 2014-12-26 CN CN201410829123.6A patent/CN104571088B/en active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6341249B1 (en) * | 1999-02-11 | 2002-01-22 | Guang Qian Xing | Autonomous unified on-board orbit and attitude control system for satellites |
CN101859146A (en) * | 2010-07-16 | 2010-10-13 | 哈尔滨工业大学 | Satellite fault prediction method based on predictive filtering and empirical mode decomposition |
CN102735877A (en) * | 2012-06-18 | 2012-10-17 | 北京控制工程研究所 | Optimization method for measuring point of measuring sensor for satellite angular rate based on correlation matrix |
CN102736616A (en) * | 2012-06-18 | 2012-10-17 | 北京控制工程研究所 | Dulmage-Mendelsohn (DM)-decomposition-based measuring point optimal configuration method for closed loop system |
CN103116357A (en) * | 2013-03-14 | 2013-05-22 | 郭雷 | Sliding-mode control method with anti-interference fault-tolerance performance |
CN103676941A (en) * | 2013-12-24 | 2014-03-26 | 北京控制工程研究所 | Satellite control system fault diagnosis method based on kinematics and dynamics model |
CN103697915A (en) * | 2013-12-24 | 2014-04-02 | 北京控制工程研究所 | Diagnostic evaluation method considering disturbing influence for satellite sensor fault |
Non-Patent Citations (2)
Title |
---|
卫星姿态控制系统的故障可诊断性评价研究;李文博等;《航天控制》;20141031;第40卷(第5期);第8-13页全文 * |
故障可诊断性评价与设计研究进展;刘文静等;《航天控制》;20111231;第29卷(第6期);第72-78页、87页全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN104571088A (en) | 2015-04-29 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104571088B (en) | Satellite control system Multipurpose Optimal Method based on fault diagnosability constraint | |
CN109583092B (en) | Intelligent mechanical system fault diagnosis method based on multi-level and multi-mode feature extraction | |
CN102647292B (en) | Intrusion detecting method based on semi-supervised neural network | |
CN100543775C (en) | The method of following the tracks of based on the 3 d human motion of many orders camera | |
CN104036102B (en) | Calculation method and device for product assembly deviation | |
CN106599580A (en) | Reconfigurable degree-based satellite on-orbit health state assessment method and assessment system | |
CN109255440B (en) | Method for predictive maintenance of power production equipment based on Recurrent Neural Networks (RNN) | |
CN106355011A (en) | Geochemical data element sequence structure analysis method and device | |
CN107103164A (en) | Unmanned plane performs the distribution method and device of multitask | |
CN110216680B (en) | Cloud-ground cooperative fault diagnosis system and method for service robot | |
CN105425583B (en) | The control method of penicillin production process based on coorinated training LWPLS | |
Kaveh et al. | An efficient two‐stage method for optimal sensor placement using graph‐theoretical partitioning and evolutionary algorithms | |
CN113404502B (en) | Shield hob abrasion monitoring device and method based on ballast piece morphology | |
CN104408760A (en) | Binocular-vision-based high-precision virtual assembling system algorithm | |
CN102736616B (en) | Dulmage-Mendelsohn (DM)-decomposition-based measuring point optimal configuration method for closed loop system | |
CN104636486B (en) | A kind of user characteristics abstracting method and draw-out device based on the conversion of non-negative alternating direction | |
CN106533759B (en) | A kind of link prediction method based on path entropy in multitiered network | |
CN111738420A (en) | Multi-scale sampling-based electromechanical equipment state data completion and prediction method | |
Hu et al. | Integrating CART algorithm and multi-source remote sensing data to estimate sub-pixel impervious surface coverage: a case study from Beijing Municipality, China | |
CN109584668A (en) | A kind of rock tunnel(ling) machine training platform based on virtual reality and big data | |
CN101901317B (en) | Growing hierarchical self-organizing maps (GHSOM)-based intrusion detection method for neural network | |
CN109377810A (en) | The monitoring of stope airflow parameter and automatic adjustment experiment porch and system | |
CN110930541B (en) | Method for analyzing working condition state of agricultural machine by using GPS information | |
CN103699121B (en) | Analytical redundancy relationship-based diagnosability determination method for satellite control system sensors | |
CN104176273B (en) | Target asteroid selection method for manned asteroid exploration |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |