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 PDF

Info

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
Application number
CN201410829123.6A
Other languages
Chinese (zh)
Other versions
CN104571088A (en
Inventor
王大轶
刘文静
邢琰
何英姿
刘成瑞
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Institute of Control Engineering
Original Assignee
Beijing Institute of Control Engineering
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Beijing Institute of Control Engineering filed Critical Beijing Institute of Control Engineering
Priority to CN201410829123.6A priority Critical patent/CN104571088B/en
Publication of CN104571088A publication Critical patent/CN104571088A/en
Application granted granted Critical
Publication of CN104571088B publication Critical patent/CN104571088B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G05CONTROLLING; REGULATING
    • G05BCONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B23/00Testing or monitoring of control systems or parts thereof
    • G05B23/02Electric testing or monitoring
    • GPHYSICS
    • G05CONTROLLING; REGULATING
    • G05BCONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/0265Adaptive 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
    • GPHYSICS
    • G05CONTROLLING; REGULATING
    • G05BCONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04Adaptive 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

Satellite control system Multipurpose Optimal Method based on fault diagnosability constraint
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:g1x+f1
C2:g2y+f2
C3:g3z+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:g1x+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.
CN201410829123.6A 2014-12-26 2014-12-26 Satellite control system Multipurpose Optimal Method based on fault diagnosability constraint Active CN104571088B (en)

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)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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

Patent Citations (7)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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