A state inversion method for electromechanical actuators based on digital twins

Abstract

The present invention provides a state inversion method for electromechanical actuators based on digital twins, including obtaining multi-source data information under EMA multi-health working states, and dividing the multi-source data information into different test data sets in the time dimension according to health indicators; preprocessing the multi-source data information in the data set, merging similar data to form a verified and adjusted multi-source information data set; injecting each type of multi-source data information in the multi-source information data set into the established EMA digital twin model; optimizing the digital twin model with a genetic algorithm, identifying, estimating and inverting the parameters of the actual operation to the twin space. Based on the consistency judgment of the virtual and real states, the state parameters that are unmeasurable or difficult to measure in the physical space are reliably approximated in the virtual space. The present invention can improve the consistency between the state inversion of complex electromechanical equipment and the actual operating state and the accuracy of parameter inversion.

Description

Technical Field

[0001] The present invention relates to the technical field of aviation electromechanical actuator control, and in particular to an electromechanical actuator state inversion method based on digital twins. Background Art

[0002] The development of more-electric / all-electric aircraft (MEA / AEA) has attracted more and more attention, and the related power-by-wire (PBW) technology has also been widely used in the aerospace field. In this technology, power-by-wire actuators are widely used. Electromechanical actuator (EMA) is a new type of power-by-wire actuator, which has many advantages over traditional hydraulic and electro-hydraulic actuators. It not only provides efficient energy conversion, but also reduces energy loss and environmental pollution, improves system reliability and maintainability, and reduces weight and space. Due to its superior performance, electromechanical actuators are increasingly widely used in more-electric / all-electric aircraft. As a potential new type of electromechanical equipment, electromechanical actuators have broad application prospects and provide more opportunities and challenges for the future aerospace industry.

[0003] The high power density design of the electromechanical actuation (EMA) system is integrated with a closed transmission structure and usually operates in limited space and harsh environments. Due to the limited installation space, EMA cannot be fully configured with sensors, and some state parameters related to health status are difficult to obtain directly through sensors and need to be obtained with the help of state inversion and parameter identification technology. Therefore, state inversion and health status assessment have a strong correlation, which in turn determines whether it can be widely used on aircraft.

[0004] Digital twin technology can assist in the state inversion of electromechanical actuation (EMA) systems. This is based on the basic requirements for the operation and optimization of digital twins throughout the life cycle. Digital twins are information complexes that exist relatively independently of physical entities in a virtual scope that is different from the physical space and can operate throughout the entire life cycle. The operation of the entire life cycle can be understood as when there is a state of a physical entity, there must be a digital twin corresponding to it, and there is a virtual state corresponding to it. This evolution helps to provide effective assistance and grip in the inversion, monitoring and predictive maintenance of equipment operation status. At the same time, based on the injection, integration and fusion of data, DT will map the actual operating physical state information inversion to the virtual space, which will also be of great help in optimizing decisions, monitoring health status, and predicting the remaining useful life (RUL) to achieve condition-based maintenance (CBM) and other fields. Summary of the invention

[0005] In order to solve the above technical problems, the present invention proposes a state inversion method of an electromechanical actuator based on digital twins, which obtains multi-source data information under multiple health working states of EMA, and divides the multi-source data information into different test data sets in the time dimension according to health indicators; pre-processes the multi-source data information in the data set, merges similar data to form a verified and adjusted multi-source information data set; injects each type of multi-source data information in the multi-source information data set into the established EMA digital twin model; uses a genetic algorithm to optimize the digital twin model, and identifies, estimates and inverts the actual operating parameters into the twin space.

[0006] The object of the present invention is to provide a state inversion method for an electromechanical actuator based on digital twins, comprising obtaining multi-source data information under multiple health working states of EMA, and dividing the multi-source data information into different test data sets according to health indicators in the time dimension, and also comprising the following steps:

[0007] Step 1: Preprocess the multi-source data information, merge similar data from the time-frequency domain feature dimension to form a verified and adjusted multi-source information data set;

[0008] Step 2: Inject each type of multi-source data information in the multi-source information data set into the established EMA digital twin model, and use a genetic algorithm to optimize the digital twin model to identify, estimate, and inversely map the actual operating parameters to the twin space;

[0009] Step 3: Determine the consistency between the actual running state and the twin mapping state.

[0010] Preferably, the multi-source data information includes a motor current signal (I 1 ,I 2 ,I 3 )、Motor voltage signal (U 1 ,U 2 ,U 3 ), speed signal (V), torque signal (T), vibration signal (A) and displacement sensor signal (X), a total of 10 dimensions of raw data information.

[0011] In any of the above solutions, preferably, the test data set includes a wear four-state data set D 1～4 , expressed as: D 1 ={I1 1,I1 2,I1 3,U1 1,U1 2,U1 3,V 1 ,T 1 ,A 1 ,X 1}、D 2 ={I2 1,I2 2,I2 3,U21,U22,U2 3,V2 ,T 2 ,A 2 ,X 2}、D 3 ={I3 1,I3 2,I3 3,U3 1,U3 2,U3 3,V 3 ,T 3 ,A 3 ,X 3} and D 4 ={I41,I4 2,I4 3,U4 1,U4 2,U4 3,V 4 ,T 4 ,A 4 ,X 4}, where D 1 is the normal data set of the system, D 2 is a light wear data set, D 3 is the moderate wear data set, D 4 This is a severe wear data set.

[0012] In any of the above schemes, preferably, the four wear states include normal system, light wear, moderate wear and heavy wear.

[0013] In any of the above schemes, preferably, step 1 includes the following sub-steps:

[0014] Step 11: Extract time domain features;

[0015] Step 12: Extract frequency domain features;

[0016] Step 13: Extract time-frequency domain features;

[0017] Step 14: Calculate the cosine distance dis(x i ,x j ), which measures the difference between two vectors.

[0018] In any of the above solutions, preferably, the time domain characteristics include a waveform factor S, a crest factor C, an impulse factor I, a kurtosis factor K and a margin factor L,

[0019] The calculation formula of the waveform factor S is:

[0020]

[0021] The calculation formula of the crest factor C is:

[0022]

[0023] The calculation formula of the pulse factor I is:

[0024]

[0025] The calculation formula of the kurtosis factor K is:

[0026]

[0027] The calculation formula of the margin factor L is:

[0028]

[0029] Among them, x(t) is the data set A 1 t is the number of detection signals, and N is the total number of collected signals.

[0030] In any of the above solutions, preferably, the frequency domain features include centroid frequency FC, mean square frequency MSF and frequency variance VF,

[0031] The calculation formula of the center of gravity frequency FC is:

[0032]

[0033] The calculation formula of the mean square frequency MSF is:

[0034]

[0035] The calculation formula of the frequency variance VF is:

[0036]

[0037] Where s(f) is the power spectrum function, f is the frequency, and df is the frequency differential symbol.

[0038] In any of the above solutions, preferably, the calculation formula of the power spectrum function s(f) is:

[0039]

[0040] Wherein, F[] represents Fourier transform, and m=1,2,…,N / 2.

[0041] In any of the above solutions, preferably, step 13 includes calculating the energy value E(j,i) of the i-th node on the j-th layer, and the formula is:

[0042]

[0043] The energy percentage characteristics D of each frequency band in the signal energy spectrum after wavelet packet decomposition i As the fault diagnosis characteristic value, the calculation formula is:

[0044]

[0045] Among them, pv is the wavelet transform coefficient, is the square of the norm, D i It is the energy percentage characteristic of each frequency band in the signal energy spectrum after wavelet packet decomposition.

[0046] In any of the above solutions, preferably, the cosine distance dist(x i ,x j ) is calculated as

[0047]

[0048] Among them, x i is the feature vector of the i-th sample, x j is the feature vector of the jth sample, x ik is the k-th eigenvalue of the eigenvector of the i-th sample, x jk is the k-th eigenvalue of the eigenvector of the j-th sample, k is the dimension of the eigenvector, θ is the angle between the eigenvectors, and n is the total number of dimensions of the eigenvector.

[0049] In any of the above solutions, preferably, the EMA digital twin model is

[0050]

[0051]

[0052] Among them, u a is the phase a voltage of the motor three-phase voltage, u b is the b-phase voltage of the motor three-phase voltage, u c is the c-phase voltage of the motor three-phase voltage, R s (T) is the resistance of the motor stator winding, i a is the a-phase current of the motor three-phase current, i b is the b-phase current of the motor three-phase current, i c is the c-phase current of the motor three-phase current, d is the d-phase current of the motor three-phase current, dt is the differential sign of time, ψ a is the flux linkage of motor phase a, ψ b is the flux linkage of motor phase b, ψ c is the flux linkage of phase c of the motor, L a is the self-inductance of phase a winding, L b is the self-inductance of phase b winding, L c is the self-inductance of phase c winding, M ab is the mutual inductance of the ab phase winding, M ac is the mutual inductance of the ac phase winding, M ba is the mutual inductance of phase b winding, M bc is the mutual inductance of bc phase winding, M cais the ca phase winding mutual inductance, M cb is the mutual inductance of cb phase winding, ψ f is the permanent magnet flux, θ e The electrical angle.

[0053] In any of the above solutions, preferably, step 2 further includes applying Park transformation to the voltage equation and flux equation in the three-phase stationary coordinate system, and taking into account the influence of temperature on the winding resistance, to obtain two voltages and flux in the dq axis.

[0054] The voltage equation is

[0055]

[0056] The magnetic flux equation is

[0057]

[0058] Among them, u d is the d-axis voltage after park transformation, u q is the q-axis voltage after park transformation, is the correction value of the stator winding resistance, i d is the d-axis current, i q is the q-axis current, ψ d is the d-axis magnetic flux, ψ q is the q-axis magnetic flux, ω m is the mechanical angular velocity of the motor rotor, L d is the d-axis self-inductance, L q is the q-axis self-inductance, ψ f (T) is the magnetic flux, is the correction value of permanent magnet flux, ψ f (T) is the magnetic flux.

[0059] In any of the above schemes, preferably, the electromagnetic torque equation of the motor is:

[0060]

[0061] Where p is the number of pole pairs of the motor.

[0062] In any of the above solutions, preferably, the kinematic equation of the motor is:

[0063]

[0064] Among them, T L is the load torque, J is the moment of inertia of the motor load converted to the motor output shaft end, B is the damping coefficient of the resistor, dω is the differential sign of the angular velocity, and ω is the angular velocity.

[0065] In any of the above solutions, preferably, in the EMA dynamic model, the torque equation at the input end of the ball screw is:

[0066]

[0067] The dynamic equation of the screw is:

[0068]

[0069] Among them, T L2 is the ball screw input torque, F G is the axial force between the screw and the nut, c bs is the axial equivalent stiffness between the ball screw and the nut, T f3 is the equivalent friction torque between the ball screw and the nut, θ n is the rotation angle, ε 2 is the clearance between the ball screw and the nut, x ema is the displacement of the screw end, d bs is the axial equivalent viscous friction coefficient between the ball screw and the nut, ω n is the angular velocity, v ema is the moving speed of the screw end, and l is the lead of the ball screw.

[0070] In any of the above schemes, preferably, step 2 further includes generating data in the virtual space after the physical acquisition data is injected into the EMA digital twin, and the generated data is generated by D DT ={IDT 1,IDT2,IDT 3,UDT 1,UDT 2,UDT3,V DT ,T DT ,A DT ,X DT} to indicate that

[0071] Among them, IDT 1, IDT 2, IDT 3, UDT 1, UDT 2, UDT 3 are used as input quantities and are consistent with physical information. DT ,T DT ,A DT ,X DT The quantity to be inverted needs to be inverted through an optimized model and kept consistent with the physical space.

[0072] In any of the above solutions, preferably, the genetic algorithm comprises the following sub-steps:

[0073] Step 21: Binary encode the parameters to be inverted to form genes and chromosomes in the genetic algorithm;

[0074] Step 22: Initialize and select the initial population for optimization;

[0075] Step 23: Generate offspring population using crossover and mutation operators;

[0076] Step 24: Conduct fitness assessment and retain excellent individuals to ensure good genes and chromosomes;

[0077] Step 25: After iteration, select the best individual as the current optimal solution;

[0078] Step 26: Parameter decoding to obtain the optimal inversion parameters.

[0079] In any of the above schemes, preferably, the consistency determination method is to determine dist(D PT ,D DT )<δ holds true,

[0080] Among them, D PT It is the information collected physically. 1～4 , D DT Map the data to the data generated in the digital twin.

[0081] In any of the above solutions, preferably, step 3 further includes quantifying the accuracy of the inversion using the leave-one-out method. When based on an N sample set X={X (i) |X (i) =(X(i)1,X(i)2,…,X(i)3),i=1,2,3,…,N}, the inversion model is represented by r X (X), the formula is

[0082] r x (X) = β T ψ(X)

[0083] Among them, β T is the weight vector of each polynomial of the inversion model, and ψ is each polynomial of the inversion model.

[0084] At the same time, set the leave-one-out sample set to have N-1 samples, denoted as X ～j ={X (i) |X (i) =(X(i)1,X(i)2,…,X(i)3),i=1,2,3,…,j-1,j+1,…,N} and It is defined as the inversion model based on the leave-one-sample set, and the inversion error is defined as r x (X) and The difference

[0085]

[0086] The accumulation of all inversion errors on the sample set is defined as the leave-one-out error, and the formula is:

[0087]

[0088] The inversion accuracy is:

[0089]

[0090] Where Y = {Y 1 ,Y 2 ,…,Y N} is the parameter after inversion. When the inversion accuracy is higher than the preset value, the inversion result is considered to be credible and consistent with the actual operating status.

[0091] The present invention proposes a state inversion method for an electromechanical actuator based on digital twins, which improves the consistency between the state inversion of complex electromechanical equipment and the actual operating state and the accuracy of parameter inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0092] Figure 1 It is a flowchart of a preferred embodiment of the electromechanical actuator state inversion method based on digital twin according to the present invention.

[0093] Figure 2 It is a flow chart of another preferred embodiment of the electromechanical actuator state inversion method based on digital twin according to the present invention.

[0094] Figure 3 A schematic diagram of a merging and partitioning architecture of a feature-based multi-source information data set according to an embodiment of the electromechanical actuator state inversion method based on digital twins of the present invention.

[0095] Figure 4 It is a schematic diagram of the architecture of an embodiment of digital twin parameter inversion based on genetic algorithm according to the digital twin-based electromechanical actuator state inversion method of the present invention.

[0096] Figure 5 It is a schematic diagram for comparing the virtual-real states of an embodiment of the electromechanical actuator state inversion method based on digital twin according to the present invention.

[0097] Figure 6 It is a schematic diagram of an embodiment of the state inversion accuracy convergence based on digital twin of the electromechanical actuator state inversion method based on digital twin according to the present invention. DETAILED DESCRIPTION

[0098] The present invention is further described below in conjunction with the accompanying drawings and specific embodiments.

[0099] Embodiment 1

[0100] like Figure 1As shown, step 100 is executed to obtain multi-source data information under the multi-health working state of EMA, and the multi-source data information is divided into different test data sets in the time dimension according to the health index. The multi-source data information includes the motor current signal (I 1 ,I 2 ,I 3 )、Motor voltage signal (U 1 ,U 2 ,U 3 ), speed signal (V), torque signal (T), vibration signal (A) and displacement sensor signal (X), a total of 10 dimensions of raw data information;

[0101] The test data set includes a wear four-state data set D 1～4 , expressed as: D 1 ={I11,I12,I13,U11,U12,U13,V 1 ,T 1 ,A 1 ,X 1}、D 2 ={I21,I22,I23,U21,U22,U23,V 2 ,T 2 ,A 2 ,X 2}、D 3 ={I31,I32,I33,U31,U32,U33,V 3 ,T 3 ,A 3 ,X 3} and D 4 ={I41,I42,I43,U41,U42,U43,V 4 ,T 4 ,A 4 ,X 4}, where D 1 is the normal data set of the system, D 2 is a light wear data set, D 3 is the moderate wear data set, D 4 This is a severe wear data set;

[0102] The four wear states include normal system, light wear, medium wear and heavy wear.

[0103] Execute step 110 to pre-process the multi-source data information, merge similar data from the time-frequency domain feature dimension to form a verified and adjusted multi-source information data set, including the following sub-steps:

[0104] Execute step 111 to extract time domain features, where the time domain features include a waveform factor S, a crest factor C, an impulse factor I, a kurtosis factor K, and a margin factor L.

[0105] The calculation formula of the waveform factor S is:

[0106]

[0107] The calculation formula of the crest factor C is:

[0108]

[0109] The calculation formula of the pulse factor I is:

[0110]

[0111] The calculation formula of the kurtosis factor K is:

[0112]

[0113] The calculation formula of the margin factor L is:

[0114]

[0115] Among them, x(t) is the data set A 1 t is the number of detection signals, and N is the total number of collected signals.

[0116] Execute step 112 to extract frequency domain features, where the frequency domain features include centroid frequency FC, mean square frequency MSF, and frequency variance VF.

[0117] The calculation formula of the center of gravity frequency FC is:

[0118]

[0119] The calculation formula of the mean square frequency MSF is:

[0120]

[0121] The calculation formula of the frequency variance VF is:

[0122]

[0123] Where s(f) is the power spectrum function, f is the frequency, and df is the frequency differential symbol.

[0124] In any of the above solutions, preferably, the calculation formula of the power spectrum function s(f) is:

[0125]

[0126] Wherein, F[] represents Fourier transform, and m=1,2,…,N / 2.

[0127] Execute step 113 to extract time-frequency domain features, including calculating the energy value E(j,i) of the i-th node on the j-th layer. The formula is:

[0128]

[0129] The energy percentage characteristics D of each frequency band in the signal energy spectrum after wavelet packet decomposition i As the fault diagnosis characteristic value, the calculation formula is:

[0130]

[0131] Among them, p v is the wavelet transform coefficient, is the square of the norm, D i It is the energy percentage characteristic of each frequency band in the signal energy spectrum after wavelet packet decomposition.

[0132] Execute step 114 to calculate the cosine distance dis(x i ,x j ), used to measure the difference between two vectors, the cosine distance dist(x i ,x j ) is calculated as

[0133]

[0134] Among them, x i is the feature vector of the i-th sample, x j is the feature vector of the jth sample, x ik is the k-th eigenvalue of the eigenvector of the i-th sample, x jk is the k-th eigenvalue of the eigenvector of the j-th sample, k is the dimension of the eigenvector, θ is the angle between the eigenvectors, and n is the total number of dimensions of the eigenvector.

[0135] Execute step 120, inject each type of multi-source data information in the multi-source information data set into the established EMA digital twin model, and use a genetic algorithm to optimize the digital twin model, identify, estimate and inversely map the actual operating parameters to the twin space, and the EMA digital twin model is

[0136]

[0137]

[0138] Among them, u a is the phase a voltage of the motor three-phase voltage, u bis the b-phase voltage of the motor three-phase voltage, u c is the c-phase voltage of the motor three-phase voltage, R s (T) is the resistance of the motor stator winding, i a is the a-phase current of the motor three-phase current, i b is the b-phase current of the motor three-phase current, i c is the c-phase current of the motor three-phase current, d is the d-phase current of the motor three-phase current, dt is the differential sign of time, ψ a is the flux linkage of phase a of the motor, ψ b is the flux linkage of motor phase b, ψ c is the motor c phase flux, L a is the self-inductance of phase a winding, L b is the self-inductance of phase b winding, L c is the self-inductance of phase c winding, M ab is the mutual inductance of the ab phase winding, M ac is the mutual inductance of the ac phase winding, M ba is the mutual inductance of phase b winding, M bc is the mutual inductance of bc phase winding, M ca is the ca phase winding mutual inductance, M cb is the mutual inductance of cb phase winding, ψ f is the permanent magnet flux, θ e The electrical angle.

[0139] Apply Park transformation to the voltage equation and flux equation in the three-phase stationary coordinate system, and consider the influence of temperature on the winding resistance, and obtain the two voltages and flux equations in the dq axis.

[0140] The voltage equation is

[0141]

[0142] The magnetic flux equation is

[0143]

[0144] Among them, u d is the d-axis voltage after park transformation, u q is the q-axis voltage after park transformation, is the correction value of the stator winding resistance, i d is the d-axis current, i q is the q-axis current, ψ d is the d-axis magnetic flux, ψ q is the q-axis magnetic flux, ω m is the mechanical angular velocity of the motor rotor, L d is the d-axis self-inductance, L q is the q-axis self-inductance, ψ f (T) is the magnetic flux, is the correction value of permanent magnet flux, ψ f (T) is the magnetic flux.

[0145] The electromagnetic torque equation of the motor is

[0146]

[0147] Where p is the number of pole pairs of the motor.

[0148] The kinematic equation of the motor is

[0149]

[0150] Among them, T L is the load torque, J is the moment of inertia of the motor load converted to the motor output shaft end, B is the damping coefficient of the resistor, dω is the differential sign of the angular velocity, and ω is the angular velocity.

[0151] In the EMA dynamic model, the torque equation at the input end of the ball screw is:

[0152]

[0153] The dynamic equation of the screw is:

[0154]

[0155] Among them, T L2 is the ball screw input torque, F G is the axial force between the screw and the nut, c bs is the axial equivalent stiffness between the ball screw and the nut, T f3 is the equivalent friction torque between the ball screw and the nut, θ n is the rotation angle, ε 2 is the clearance between the ball screw and the nut, x ema is the displacement of the screw end, d bs is the axial equivalent viscous friction coefficient between the ball screw and the nut, ω n is the angular velocity, v ema is the moving speed of the screw end, and l is the lead of the ball screw.

[0156] When the physical data is injected into the EMA digital twin, data is generated in the virtual space. DT ={IDT 1,IDT 2,IDT 3,UDT 1,UDT 2,UDT 3,V DT ,T DT ,A DT ,X DT} to indicate that

[0157] Among them, IDT 1, IDT 2, IDT 3, UDT 1, UDT 2, UDT 3 are used as input quantities and are consistent with physical information. DT ,T DT ,A DT ,X DT The quantity to be inverted needs to be inverted through an optimized model and kept consistent with the physical space.

[0158] The consistency determination method is to determine dist(D PT ,D DT )<δ holds true,

[0159] Among them, D PT It is the information collected physically. 1～4 , D DT Map the data to the data generated in the digital twin.

[0160] The genetic algorithm comprises the following sub-steps:

[0161] Execute step 122 to binary encode the parameters to be inverted to form genes and chromosomes in the genetic algorithm;

[0162] Execute step 123 to initialize and select an initial population for optimization;

[0163] Execute step 123 to generate a progeny population using crossover and mutation operators;

[0164] Execute step 124 to perform adaptability determination, retain excellent individuals to ensure good genes and chromosomes;

[0165] Execute step 125, after iteration, select the best individual as the current optimal solution;

[0166] Execute step 126, parameter decoding, and obtain the optimal inversion parameters.

[0167] Execute step 130 to determine the consistency of the actual operating state and the twin mapping state, and quantify the accuracy of the inversion using the leave-one-out method. The consistency determination method is to determine dist(D PT ,D DT )<δ holds true,

[0168] Among them, D PT It is the information collected physically. 1～4 , D DT Map the data to the data generated in the digital twin.

[0169] When based on an N sample set X={X (i) |X (i)=(X(i)1,X(i)2,…,X(i)3),i=1,2,3,…,N}, the inversion model is represented by r X (X), the formula is

[0170] r x (X) = β T ψ(X)

[0171] Among them, β T is the weight vector of each polynomial of the inversion model, and ψ is each polynomial of the inversion model.

[0172] At the same time, set the leave-one-out sample set to have N-1 samples, denoted as X ～j ={X (i) |X (i) =(X(i)1,X(i)2,…,X(i)3),i=1,2,3,…,j-1,j+1,…,N} and It is defined as the inversion model based on the leave-one-sample set, and the inversion error is defined as r x (X) and The difference

[0173]

[0174] The accumulation of all inversion errors on the sample set is defined as the leave-one-out error, and the formula is:

[0175]

[0176] The inversion accuracy is:

[0177]

[0178] Where Y = {Y 1 ,Y 2 ,…,Y N} is the parameter after inversion. When the inversion accuracy is higher than the preset value, the inversion result is considered to be credible and consistent with the actual operating status.

[0179] Embodiment 2

[0180] The present invention proposes a state inversion method of an electromechanical actuator (EMA) based on digital twins, and the specific steps are as follows:

[0181] Step 1: First, obtain multi-source data information under EMA multi-health working status, and divide the multi-source data information into different test data sets according to health indicators in the time dimension.

[0182] Step 2: Preprocess the multi-source data information in the collected data set, merge similar data from the time-frequency domain feature dimensions to form a verified and adjusted multi-source information data set.

[0183] Step three: Inject each type of multi-source data information in the multi-source information data set into the established EMA digital twin model, and use a genetic algorithm to optimize the digital twin model to identify, estimate and inversely map the actual operating parameters to the twin space.

[0184] Step 4: Verify the consistency between the actual operating state and the twin mapping state to ensure that the state parameters that are unmeasurable or difficult to measure in the physical space are reliably approximated in the virtual space.

[0185] The advantages and positive effects of the present invention are:

[0186] (1) Aiming at the problem of incomplete sensor layout of electromechanical actuator (EMA), a state inversion method of electromechanical actuator (EMA) based on digital twin is proposed. This method uses the digital twin method to reliably approximate the state parameters that are unmeasurable or difficult to measure in the physical space in the virtual space.

[0187] (2) Use multi-sensor information features to merge data sets to ensure the relevance of data and health status. Merge similar data through time-frequency domain feature dimensions to form a verified and adjusted multi-source information data set, and then inject data under different health states into the digital twin model for inversion.

[0188] Embodiment 3

[0189] like Figure 2 As shown in the figure, the process of the electromechanical actuator (EMA) state inversion method based on digital twin is as follows:

[0190] Step 1: Obtain multi-source data information under multiple healthy working states of EMA. Taking the EMA wear state as an example, collect signals representing normal system, light wear, moderate wear and heavy wear respectively. When EMA operates under different wear states, it will affect the change of its own sensor signal, and at the same time, it will also affect the sensor signals of other parameters in its coupling loop, causing the operating state of the entire device to change. Therefore, it is possible to obtain sensor signals with high sensitivity to the health state in the EMA actuator system, including the motor current signal (I 1 ,I 2 ,I 3 )、Motor voltage signal (U 1 ,U 2 ,U 3 ), speed signal (V), torque signal (T), vibration signal (A) and displacement sensor signal (X), a total of 10 dimensions of raw data information. It is divided into different data sets according to the preset health state. Taking the four wear states as an example, the data sets are represented as D 1～4 , D 1={I1 1,I1 2,I1 3,U1 1,U1 2,U1 3,V 1 ,T 1 ,A 1 ,X 1}、D 2 ={I2 1,I2 2,I2 3,U2 1,U2 2,U2 3,V 2 ,T 2 ,A 2 ,X 2}、D 3 ={I3 1,I3 2,I3 3,U3 1,U3 2,U3 3,V 3 ,T 3 ,A 3 ,X 3} and D 4 ={I4 1,I4 2,I4 3,U4 1,U4 2,U4 3,V 4 ,T 4 ,A 4 ,X 4}.

[0191] Step 2: Time-sequence the system multi-source data training set obtained in step 1 to form a time series. Since the data is collected and segmented according to the subjective health status, it is necessary to use algorithms to merge and adjust similar data in the feature dimension. The specific implementation process is as follows: Figure 3 As shown. It can be seen that the adjusted dataset avoids data confusion at different stages and retains data with unique characteristics, ensuring the usability of the data.

[0192] In order to determine the required multi-source information data set, each type of information is preprocessed and feature extracted. The present invention focuses on extracting signal features that are sensitive to health status from three aspects: time domain, frequency domain and time-frequency domain. There are five time domain features, including waveform factor S, crest factor C, pulse factor I, kurtosis factor K and margin factor L. Frequency domain feature extraction is based on the frequency domain feature parameters of common power spectrum analysis, including centroid frequency, mean square frequency and frequency variance. These three dimensionless parameters are sensitive to changes in health status. Time-frequency domain feature extraction extracts signal energy features as time-frequency domain features based on the wavelet packet decomposition method, solves the signal energy at different decomposition scales, and adaptively selects the corresponding frequency band to match the signal spectrum according to the signal characteristics and analysis requirements, and arranges these energy values ​​into feature vectors in scale order for classification.

[0193] (1) Time domain feature extraction

[0194] The calculation formula is as follows, where x(t) (t = 1, 2, ..., N) refers to the data set A1 The detection signal in, N is the total number of signals collected.

[0195] ① Waveform factor S:

[0196]

[0197] ②Crest Factor C:

[0198]

[0199] ③ Pulse Factor I:

[0200]

[0201] ④ Kurtosis factor K:

[0202]

[0203] ⑤Margin factor L:

[0204]

[0205] (2) Frequency domain feature extraction:

[0206] Among them, s(f) is the power spectrum function. It can be expressed as:

[0207]

[0208] Where F[] represents Fourier transform, where t = 1, 2, ..., N / 2. The frequency domain feature calculation formula is as follows:

[0209] ⑥Center of gravity frequency FC

[0210]

[0211] ⑦ Mean square frequency MSF:

[0212]

[0213] ⑧Frequency variance VF:

[0214]

[0215] (3) Time-frequency domain feature extraction

[0216] The specific calculation formula is:

[0217]

[0218] Where E(j,i) represents the energy value of the i-th node on the j-th layer, p v are the wavelet transform coefficients, The present invention adopts db4 wavelet to perform three-layer decomposition, and adaptively divides it into 8 frequency bands from low frequency to high frequency, that is, j=3, i=1,2,…,8.

[0219] ⑨ After decomposing the wavelet packet, the energy percentage characteristics of each frequency band in the signal energy spectrum are obtained. i As the fault diagnosis characteristic value. The calculation formula is:

[0220]

[0221] In summary, we can get the multi-source information dataset D 1～4 16 features of each detection signal in the . By analyzing the collected sensor signals in multiple dimensions of time domain, frequency domain and time-frequency domain, in order to provide a basis for merging signal data sets, it is necessary to calculate based on the characteristic distance of the data. The cosine distance can be used to measure the difference between two vectors. The formula is as follows:

[0222]

[0223] The value range of the angle cosine is [-1,1]. The larger the cosine, the smaller the angle between the two vectors, and the smaller the cosine, the larger the angle between the two vectors. When the directions of the two vectors coincide, the cosine takes the maximum value of 1, and when the directions of the two vectors are completely opposite, the cosine takes the minimum value of -1. Therefore, the original data set is rebuilt and adjusted according to the different feature distances to ensure the consistency of the data with the actual operating status.

[0224] Step 3: For the adjusted multi-source information data set obtained in step 2, inject the data into the EMA digital twin, where the digital twin model is as follows:

[0225] EMA motor electromagnetic model:

[0226]

[0227]

[0228] By applying Park transformation to the voltage equation and flux equation in the three-phase stationary coordinate system and taking into account the influence of temperature on the winding resistance, the two voltages and flux equations in the dq axis can be obtained.

[0229] The stator voltage equation is:

[0230]

[0231] In the formula, i d 、i q are the d-axis and q-axis currents [A], ψ, ψ q are the d-axis and q-axis magnetic flux [Wb] respectively; ω mis the mechanical angular velocity of the motor rotor [rad / s]; R s (T) is the resistance of the motor stator winding [Ω], which is affected by temperature; is the correction value of the stator winding resistance [Ω].

[0232] The stator flux equation is:

[0233]

[0234] Where, L d , L q are the d-axis and q-axis self-inductance [H] respectively; ψ f (T) is the permanent magnet flux [Wb], and this value is affected by temperature; is the correction value of permanent magnet flux [Wb].

[0235] The electromagnetic torque equation of the motor is:

[0236]

[0237] Where p is the number of pole pairs of the motor.

[0238] The kinematic equation of the motor is:

[0239]

[0240] Where, T L is the load torque [N·m]; J is the moment of inertia of the motor load converted to the motor output shaft end [kg·m 2 ]; B is the damping coefficient of the resistor [N·s / m].

[0241] EMA dynamics model:

[0242] The torque equation at the input end of the ball screw is:

[0243]

[0244] In the formula, F G is the axial force between the screw and the nut [N]; T f3 is the equivalent friction torque between the ball screw and the nut [N·m].

[0245] The dynamic equation of the screw is:

[0246]

[0247] In the formula, ε 2 is the clearance between the ball screw and the nut [m]; c bs is the axial equivalent stiffness between the ball screw and the nut [N / m]; d bsis the axial equivalent viscous friction coefficient between the ball screw and the nut [Pa·s]; x ema is the displacement of the screw end [m]; v ema is the moving speed of the screw end [m / s]; l is the ball screw lead [m].

[0248] When the physical data is injected into the EMA digital twin, data will also be generated in the virtual space. DT ={IDT 1,IDT 2,IDT 3,UDT 1,UDT 2,UDT3,V DT ,T DT ,A DT ,X DT}. Among them, IDT1, IDT 2, IDT 3, UDT 1, UDT 2, UDT 3 are consistent with the input quantity and physical information, V DT ,T DT ,A DT ,X DT The inversion quantity needs to be inverted by optimizing the model and keep it consistent with the physical space. Figure 4 As shown in the figure, the genetic algorithm is a heuristic optimization algorithm that can adaptively obtain the optimal solution by simulating the process of genetic variation in nature. The specific steps are as follows:

[0249] ① Binary encode the parameters to be inverted to form genes and chromosomes in the genetic algorithm;

[0250] ② Initialize and select the initial population for optimization;

[0251] ③ Use crossover and mutation operators to generate offspring populations;

[0252] ④ Adaptability determination, retaining excellent individuals to ensure that good genes and chromosomes can be retained;

[0253] ⑤After iteration, select the best individual as the current optimal solution;

[0254] ⑥ Parameter decoding to obtain the optimal inversion parameters.

[0255] Step 4: For the parameters or signals inverted in step 3, it is necessary to keep the virtual and real consistent. Therefore, it is necessary to make a consistency judgment, that is, dist(D PT ,D DT )<δ. Where D PT It is the information collected physically. 1～4 , D DT Mapping data to the data generated in the digital twin. Figure 5As shown in the figure, the vibration state mapped to the digital virtual space is compared with the actual vibration state. The comparison shows that the results of state inversion can be consistent with the actual state from the aspects of feature tracking and statistical information.

[0256] In addition, the accuracy of the inversion needs to be quantified. Here, the leave-one-out method is used. (i) |X (i) =(X(i)1,X(i)2,…,X(i)3),i=1,2,3,…,N}. The inversion model can be expressed as X (X).

[0257] r x (X) = β T ψ(X) (21)

[0258] At the same time, set the leave-one-out sample set to have N-1 samples, denoted as X ～j ={X (i) |X (i) =(X(i)1,X(i)2,…,X(i)3), i=1,2,3,…,j-1,j+1,…,N} and is defined as the inversion model based on a leave-one-sample set. Therefore, the inversion error is defined as r x (x) and The difference.

[0259]

[0260] The inversion error accumulation of all the sample sets is defined as the leave-one-out error, as shown below:

[0261]

[0262] The inversion accuracy is:

[0263]

[0264] Where Y = {Y 1 ,Y 2 ,…,Y N} is the parameter after inversion. When the inversion accuracy is higher than the preset value, the inversion result is considered to be credible and consistent with the actual operating status.

[0265] Embodiment 4

[0266] Figure 6The figure shows the convergence diagram of the state inversion accuracy based on digital twins. After iteration, the inversion accuracy can exceed the preset lower limit of accuracy. It can be seen that the EMA state inversion method based on digital twins can reliably approximate the state parameters that are unmeasurable or difficult to measure in the physical space in the virtual space after the consistency judgment of the virtual and real states.

[0267] In order to better understand the present invention, the above is described in detail in conjunction with the specific embodiments of the present invention, but it is not intended to limit the present invention. Any simple modifications made to the above embodiments based on the technical essence of the present invention still fall within the scope of the technical solution of the present invention. Each embodiment in this specification focuses on the differences from other embodiments, and the same or similar parts between the embodiments can be referenced to each other. For the system embodiment, since it basically corresponds to the method embodiment, the description is relatively simple, and the relevant parts can be referred to the partial description of the method embodiment.

Claims

1. A state inversion method for an electromechanical actuator based on digital twins, comprising obtaining multi-source data information under multiple healthy working states of EMA, and dividing the multi-source data information into different test data sets according to health indicators in the time dimension, characterized in that: The following steps are also included: Step 1: Preprocess the multi-source data information, merge similar data from the time-frequency domain feature dimension to form a verified and adjusted multi-source information data set; Step 2: Inject each type of multi-source data in the multi-source information data set into the established EMA digital twin model, and use a genetic algorithm to optimize the digital twin model, identify, estimate and inversely map the actual operating parameters to the twin space. The EMA digital twin model is Among them, u a is the phase a voltage of the motor three-phase voltage, u b is the b-phase voltage of the motor three-phase voltage, u c is the c-phase voltage of the motor three-phase voltage, R s (T) is the resistance of the motor stator winding, i a is the a-phase current of the motor three-phase current, i b is the b-phase current of the motor three-phase current, i c is the c-phase current of the motor three-phase current, d is the d-phase current of the motor three-phase current, dt is the differential sign of time, ψ a is the flux linkage of motor phase a, ψ b is the flux linkage of motor phase b, ψ c is the flux linkage of phase c of the motor, L a is the self-inductance of phase a winding, L b is the self-inductance of phase b winding, L c is the self-inductance of phase c winding, M ab is the mutual inductance of the ab phase winding, M ac is the mutual inductance of the ac phase winding, M ba is the mutual inductance of phase b winding, M bc is the mutual inductance of bc phase winding, M ca is the ca phase winding mutual inductance, M cb is the mutual inductance of cb phase winding, ψ f is the permanent magnet flux, θ e is the electrical angle; Apply Park transformation to the voltage equation and flux equation in the three-phase stationary coordinate system, and consider the influence of temperature on the winding resistance, and obtain the two voltages and flux equations in the dq axis. The voltage equation is The magnetic flux equation is Among them, u d is the d-axis voltage after park transformation, u q is the q-axis voltage after park transformation, is the correction value of the stator winding resistance, i d is the d-axis current, i q is the q-axis current, ψ d is the d-axis magnetic flux, ψ q is the q-axis magnetic flux, ω m is the mechanical angular velocity of the motor rotor, L d is the d-axis self-inductance, L q is the q-axis self-inductance, ψ f (T) is the magnetic flux, is the correction value of permanent magnet flux, ψ f (T) is the magnetic flux; When the physical data is injected into the EMA digital twin, data is generated in the virtual space. DT ={IDT 1,IDT 2,IDT 3,UDT 1,UDT 2,UDT 3,V DT ,T DT ,A DT ,X DT } to indicate that Among them, IDT 1, IDT 2, IDT 3, UDT 1, UDT 2, UDT 3 are used as input quantities and are consistent with physical information. DT ,T DT ,A DT ,X DT The inverted quantity needs to be inverted through the optimization model and kept consistent with the physical space; Step 3: Determine the consistency between the actual running state and the twin mapping state.

2. The electromechanical actuator state inversion method based on digital twin according to claim 1, characterized in that: The multi-source data information includes motor current signals (I1, I 2, I3), motor voltage signal (U1, U 2, U3), speed signal (V), torque signal (T), vibration signal (A) and displacement sensor signal ( ), a total of 10 dimensions of original data information.

3. The electromechanical actuator state inversion method based on digital twin according to claim 2, characterized in that: The test data set includes a wear four-state data set D 1~4 , expressed as: D1={I1 1,I1 2,I1 3,U1 1,U1 2,U1 3,V 1 ,T 1 ,A 1 ,X 1 }, D2={I2 1,I2 2,I2 3,U2 1,U2 2,U2 3,V 2 ,T 2 ,A 2 ,X 2 }, D3={I3 1,I3 2,I3 3,U31,U3 2,U3 3,V 3 ,T 3 ,A 3 ,X 3 } and D4={I4 1,I4 2,I4 3,U4 1,U4 2,U4 3,V 4 ,T 4 ,A 4 ,X 4 }, where D1 is the normal system data set, D2 is the light wear data set, D3 is the moderate wear data set, and D4 is the heavy wear data set.

4. The electromechanical actuator state inversion method based on digital twin according to claim 3, characterized in that: The step 1 includes the following sub-steps: Step 11: Extract time domain features; Step 12: Extract frequency domain features; Step 13: Extract time-frequency domain features; Step 14: Calculate the cosine distance dis(x i , x j ), which measures the difference between two vectors.

5. The electromechanical actuator state inversion method based on digital twin according to claim 4, characterized in that: The time domain characteristics include waveform factor S, crest factor C, pulse factor I, kurtosis factor K and margin factor L. The calculation formula of the waveform factor S is: The calculation formula of the crest factor C is: The calculation formula of the pulse factor I is: The calculation formula of the kurtosis factor K is: The calculation formula of the margin factor L is: Wherein, x(t) is the detection signal in the data set A1, t is the number of detection signals, and N is the total number of signals collected.

6. The electromechanical actuator state inversion method based on digital twin according to claim 5, characterized in that: The frequency domain features include the centroid frequency FC, the mean square frequency MSF and the frequency variance VF. The calculation formula of the center of gravity frequency FC is: The calculation formula of the mean square frequency MSF is: The calculation formula of the frequency variance VF is: Where s(f) is the power spectrum function, f is the frequency, and df is the frequency differential symbol.

7. The electromechanical actuator state inversion method based on digital twin according to claim 6, characterized in that: The calculation formula of the power spectrum function s(f) is: Where F[] represents Fourier transform, and m=1,2,…,N / 2.

8. The electromechanical actuator state inversion method based on digital twin according to claim 7, characterized in that: The step 13 includes calculating the energy value E(j,i) of the i-th node on the j-th layer, and the formula is: The energy percentage characteristics D of each frequency band in the signal energy spectrum after wavelet packet decomposition i As the fault diagnosis characteristic value, the calculation formula is: Among them, p v is the wavelet transform coefficient, is the square of the norm.

9. The electromechanical actuator state inversion method based on digital twin according to claim 8, characterized in that: The cosine distance dist(x i , x j ) is calculated as Among them, x i is the feature vector of the i-th sample, x j is the feature vector of the jth sample, x ik is the k-th eigenvalue of the eigenvector of the i-th sample, x jk is the k-th eigenvalue of the feature vector of the j-th sample, k is the dimension of the feature vector, θ is the angle between the feature vectors, and n is the total number of dimensions of the feature vector.

Images