Satellite attitude control system fault diagnosis and fault-tolerant reconstruction method and system

Abstract

The invention provides a satellite attitude control system fault diagnosis and fault-tolerant reconstruction method and system, and relates to the technical field of satellites, and the method comprises the steps: obtaining attitude data through a Kalman filtering state observer, building a fault diagnosis model, generating a residual sequence, and carrying out the adaptive dynamic threshold calculation to obtain a fault feature vector; and inputting the feature vector into a pre-trained deep belief network for fault diagnosis, and determining the fault type and degree. According to a diagnosis result, a mixed integer programming model considering attitude control precision and actuator power constraint is established, an improved genetic algorithm is adopted for solving, an actuator reconstruction scheme and controller parameters are obtained, and fault-tolerant control is achieved. According to the invention, faults can be diagnosed rapidly and accurately, fault-tolerant reconstruction is carried out according to fault types and degrees, the reliability and stability of a satellite attitude control system are effectively improved, and the satellite attitude control precision is ensured.

Description

Technical Field

[0001] This invention relates to satellite technology, and more particularly to a method and system for fault diagnosis and fault-tolerant reconfiguration of satellite attitude control systems. Background Technology

[0002] The satellite attitude control system is one of the key systems for satellite operation in orbit, and its performance directly affects whether the satellite can successfully complete its planned mission. With the increasing complexity of satellite missions and the extension of their on-orbit lifespan, the probability of failure in the satellite attitude control system also increases. Therefore, fault diagnosis and fault-tolerant reconfiguration of the satellite attitude control system are of great significance.

[0003] Traditional fault diagnosis methods for satellite attitude control systems primarily rely on hardware redundancy, such as using multiple identical sensors or actuators. While simple and reliable, this approach increases the satellite's weight and cost. Furthermore, traditional fault-tolerant reconfiguration methods typically employ pre-defined reconfiguration schemes, lacking the flexibility to adjust to actual fault conditions, resulting in limited reconfiguration accuracy and efficiency.

[0004] The existing technology has the following main defects and shortcomings:

[0005] 1. Fault diagnosis methods based on hardware redundancy will increase the weight and cost of the satellite, and are not suitable for microsatellites that are sensitive to weight and cost.

[0006] 2. Traditional fault-tolerant refactoring methods lack flexibility and cannot be adjusted according to actual fault conditions, resulting in limited refactoring accuracy and efficiency.

[0007] 3. Traditional fault diagnosis methods are difficult to effectively handle complex faults, such as simultaneous failures of multiple components or situations where the degree of fault is uncertain. Summary of the Invention

[0008] The present invention provides a method and system for fault diagnosis and fault-tolerant reconfiguration of satellite attitude control systems, which can at least solve some of the problems existing in the prior art.

[0009] A first aspect of the present invention provides a method for fault diagnosis and fault-tolerant reconfiguration of a satellite attitude control system, comprising:

[0010] Attitude angle data, angular velocity data, reaction wheel speed data, and torque data of the satellite attitude control system are collected to establish a Kalman filter state observer. The Kalman filter state observer is embedded into a fault diagnosis model, and a residual sequence is generated using the fault diagnosis model. An adaptive dynamic threshold is then calculated on the residual sequence to obtain a fault feature vector. The fault feature vector is input into a pre-trained deep belief network in the fault diagnosis model to obtain a fault diagnosis result, which includes the fault type and fault severity.

[0011] Based on the fault diagnosis results, the faulty component is identified, and a mixed integer programming model is established based on the fault type and fault degree of the faulty component. Attitude control accuracy constraints and actuator power constraints are set in the mixed integer programming model. An improved genetic algorithm is used to solve the mixed integer programming model to obtain a fault-tolerant control strategy, which includes an actuator reconfiguration scheme and controller parameters.

[0012] The actuator operating mode in the satellite attitude control system is reconfigured according to the actuator reconfiguration scheme in the fault-tolerant control strategy; the controller of the satellite attitude control system is updated using the controller parameters in the fault-tolerant control strategy; the satellite attitude is controlled based on the updated controller to determine the satellite attitude error; when the satellite attitude error exceeds a preset threshold, the fault diagnosis model is used to re-diagnose the fault.

[0013] Establish a Kalman filter state observer, including:

[0014] The attitude angle data is used as the attitude quaternion observation, and the angular velocity data is used as the angular velocity observation. A state vector containing the attitude quaternion observation and the angular velocity observation is constructed. A satellite attitude dynamics model is established. According to the Euler motion equation, the reaction wheel speed data and the torque data are substituted into the satellite attitude dynamics model to obtain a state space expression. The state space expression is discretized within the sampling time interval to obtain a discrete state equation.

[0015] An extended Kalman filter state observer is constructed. Based on the discrete state equations, the Jacobian matrix is ​​calculated. The Jacobian matrix is ​​linearized at the current state estimate to obtain the state transition matrix. The observation matrix is ​​calculated based on the observation equations. The system noise covariance matrix and observation noise covariance matrix of the extended Kalman filter state observer are initialized. Based on the initialized system noise covariance matrix and observation noise covariance matrix, a one-step state prediction value and an error covariance prediction value are calculated. The Kalman gain matrix is ​​calculated using the one-step state prediction value and the error covariance prediction value. The Kalman gain matrix is ​​multiplied by the observation residual to obtain the state estimate correction. The state estimate correction is added to the one-step state prediction value to obtain the state estimate.

[0016] A new error covariance matrix is ​​calculated based on the state estimate. The state estimation accuracy is calculated based on the error covariance matrix. When the state estimation accuracy meets the preset accuracy requirement, the state estimate is output. The observation residual sequence is calculated. The actual covariance of the observation residual sequence is compared with the theoretical covariance to obtain an adaptive factor. The system noise covariance matrix and the observation noise covariance matrix are adjusted using the adaptive factor. The adjusted system noise covariance matrix and the observation noise covariance matrix are substituted into the extended Kalman filter state observer to obtain the Kalman filter state observer.

[0017] Establish a satellite attitude dynamics model, including:

[0018] Quaternions are used to describe the satellite attitude. A quaternion vector containing scalar and vector parts is constructed. The attitude kinematic equation is established by associating the quaternion vector with the angular velocity vector. The derivative of the quaternion vector is represented by the product of the quaternion correlation matrix and the angular velocity vector. The attitude dynamic equation is established based on the Euler rigid body dynamics principle. The satellite rotational inertia matrix, the cross product of the angular velocity vector, the control torque, and the disturbance torque are substituted into the attitude dynamic equation. The control torque is generated by the reaction wheel system.

[0019] A dynamic equation for the reaction wheel is established by substituting the angular momentum vector, driving torque, and cross product of the angular velocity vector and angular momentum vector of the reaction wheel into the dynamic equation. The attitude kinematics equation, the attitude dynamic equation, and the reaction wheel dynamic equation are combined to form a complete satellite attitude dynamic model, which includes quaternion derivative equations, angular velocity derivative equations, and reaction wheel angular momentum derivative equations. The satellite attitude dynamic model is discretized using the fourth-order Runge-Kutta numerical integration method to obtain discrete state equations containing quaternion state update functions, angular velocity state update functions, and reaction wheel angular momentum state update functions.

[0020] A process noise term is added to the discrete state equation, and the quaternion vector, angular velocity vector, and reaction wheel angular momentum vector are combined to form the state vector. The reaction wheel driving torque is used as the control input to establish a state-space expression containing process noise. An observation equation is established based on the attitude quaternion measurement value of the star sensor, the angular velocity measurement value of the gyroscope, and the angular momentum measurement value of the reaction wheel speed sensor. Measurement noise is added to the observation equation, and the state-space expression and the observation equation are combined to form a satellite attitude dynamics model.

[0021] The fault feature vector is input into a pre-trained deep belief network in the fault diagnosis model to obtain fault diagnosis results, including:

[0022] The fault feature vector is mapped to the zero-to-one interval using the minimum-maximum normalization method to obtain a normalized feature vector. The covariance matrix is ​​calculated based on the normalized feature vector, and the covariance matrix is ​​decomposed into principal components. Principal components with a cumulative contribution rate greater than a preset threshold are selected to construct a dimensionality reduction mapping matrix. The normalized feature vector is then subjected to dimensionality reduction processing through the dimensionality reduction mapping matrix to obtain a dimensionality reduction feature vector. The dimensionality reduction feature vector is input into the visible layer of a deep belief network, wherein the deep belief network contains multiple stacked restricted Boltzmann machines, and the number of neurons in each layer of the restricted Boltzmann machine decreases.

[0023] The dimensionality-reduced feature vector is processed layer by layer through the Restricted Boltzmann Machine (RBM) of the deep belief network. The product of the RBM input vector and the weight matrix, and the sum of the bias vectors are calculated. The sum of the bias vectors is activated by a normalization function to obtain the hidden layer output. The neurons in the hidden layer output are set to zero according to a preset random deactivation probability. For the hidden layer outputs that are not set to zero, numerical compensation is performed according to a scaling factor. The hidden layer output after random deactivation is used as the visual layer input of the next layer of the RBM. The feature transformation, normalization activation, random deactivation and numerical compensation process is repeated until the hidden layer output of the last layer of the RBM is obtained. The random deactivation probability and the scaling factor are determined according to the number of neurons in each layer.

[0024] The hidden layer output of the last layer of the restricted Boltzmann machine is connected to a normalized classifier. The original score of each fault type is calculated based on the product of the hidden layer output and the classifier weight matrix and the bias vector. The original score is then exponentially calculated to obtain a non-negative value. This non-negative value is divided by the sum of the exponents of all category scores to obtain the conditional probability distribution of different fault types. The category with the highest conditional probability that is greater than a preset probability threshold is selected as the fault type. The activation values ​​of the output layer of the normalized classifier are weighted and combined to obtain the fault severity. The fault type and the fault severity are combined to generate a fault diagnosis result.

[0025] Based on the fault diagnosis results, the faulty component is identified, and a mixed-integer programming model is established based on the fault type and fault severity of the faulty component, including:

[0026] Based on the fault diagnosis results, fault type, fault severity, and fault location information are obtained. A fault impact assessment matrix is ​​established based on the fault type and fault severity. The fault impact assessment matrix includes the impact coefficient of the faulty component on attitude control accuracy and the impact coefficient of the faulty component on actuator power. The fault impact assessment matrix is ​​input into a mixed integer programming model, and attitude control variables are established in the mixed integer programming model. The attitude control variables include pitch angle, yaw angle, and roll angle. The product of the attitude control variables and the fault impact assessment matrix is ​​used as the predicted value of attitude control accuracy under fault conditions.

[0027] An actuator power demand model is established based on the dynamic characteristics of the faulty component. The product of the actuator power demand model and the fault impact assessment matrix is ​​used as the actuator power prediction value under the fault state. Actuator power constraints are set in the mixed integer programming model. An attitude control accuracy threshold is determined according to the equipment attitude control accuracy requirements. An attitude control accuracy constraint is established by the difference between the attitude control accuracy prediction value and the attitude control accuracy threshold. The attitude control accuracy constraint is set in the mixed integer programming model.

[0028] An objective function is established based on the fault impact assessment matrix. The objective function includes a term for minimizing attitude control deviation, a term for minimizing actuator power consumption, and a term for fault impact compensation. A mixed integer programming model is constructed by combining the objective function, the attitude control variables, the actuator power constraints, and the attitude control accuracy constraints. The attitude control deviation minimization term is calculated based on the deviation between the attitude control variables and the desired attitude. The actuator power consumption minimization term is calculated based on the predicted actuator power value. The fault impact compensation term is calculated based on the impact coefficients in the fault impact assessment matrix.

[0029] An improved genetic algorithm is used to solve the mixed-integer programming model to obtain a fault-tolerant control strategy. This strategy includes an actuator reconfiguration scheme and controller parameters, including:

[0030] In a mixed-integer programming model, actuator reconfiguration variables and controller parameter variables are encoded as chromosomes. The actuator reconfiguration variables include actuator operating status flags, actuator output allocation coefficients, and actuator switching timing. The controller parameter variables include proportional coefficient matrices, integral coefficient matrices, and differential coefficient matrices. A multi-objective fitness function is established based on the actuator reconfiguration variables and controller parameter variables. The multi-objective fitness function includes actuator resource utilization evaluation items, control performance evaluation items, fault tolerance evaluation items, and energy consumption evaluation items.

[0031] The chromosome population is initialized according to the multi-objective fitness function; the chromosome population is stratified to construct a population stratification structure, an elite population library is established in the population stratification structure, the fitness value of the chromosomes is calculated according to the multi-objective fitness function, and non-dominated solution chromosomes are stored in the elite population library; an adaptive crossover operator is used to perform a crossover operation on the chromosome population, the crossover probability of the adaptive crossover operator is determined based on the chromosome distribution in the elite population library and the multi-objective fitness function value, and the crossover chromosomes are generated and added to the chromosome population;

[0032] A hybrid mutation operator is used to mutate the crossover chromosome population to obtain an elite population pool. Cluster analysis is performed on the chromosomes in the elite population pool, and a clustering algorithm is used to group the chromosomes. From each group, the chromosome with the optimal multi-objective fitness function value is selected to form an optimal solution set. The chromosomes in the optimal solution set are decoded to obtain a fault-tolerant control strategy. The fault-tolerant control strategy includes an actuator operating mode switching timing table, an actuator output allocation matrix, an actuator cooperative control scheme, a proportional coefficient matrix, an integral coefficient matrix, and a differential coefficient matrix.

[0033] A hybrid mutation operator is used to mutate the chromosome population after crossover to obtain an elite population pool, including:

[0034] A hybrid mutation operator is used to mutate the crossover chromosome population. This hybrid mutation operator includes discrete and continuous mutation operations. The discrete mutation operation changes the values ​​of the actuator's operating status flags by bit reversal, adjusts the distribution of the actuator's output allocation coefficients by random permutation, and modifies the timing of the actuator switching sequence by binary encoding conversion. The continuous mutation operation adjusts the parameters of the proportional coefficient matrix by Gaussian distribution perturbation, the parameters of the integral coefficient matrix by uniform distribution perturbation, and the parameters of the differential coefficient matrix by non-uniform distribution perturbation, thus obtaining the mutated chromosome population.

[0035] The mutated chromosome population is evaluated based on the multi-objective fitness function, and chromosomes with better evaluation results are selected into the elite population pool. The simulated annealing algorithm is used to perform local search optimization on the chromosomes in the elite population pool. The initial temperature of the simulated annealing algorithm is determined according to the value range of the multi-objective fitness function, the cooling coefficient of the simulated annealing algorithm is calculated in an adaptive manner, and the termination temperature of the simulated annealing algorithm is dynamically adjusted based on the number of iterations.

[0036] During the local search process of the simulated annealing algorithm, a discrete neighborhood structure is adopted for the actuator reconstruction variables, which includes bit-flipped neighborhoods, swapped neighborhoods, and insertion neighborhoods; a continuous neighborhood structure is adopted for the controller parameter variables, which includes Gaussian neighborhoods, uniform neighborhoods, and non-uniform neighborhoods; the acceptance probability of a new solution is calculated based on the multi-objective fitness function, and the acceptance probability decreases as the temperature decreases during the iteration process; the chromosome optimized by the local search is then updated to the elite population library.

[0037] A second aspect of the present invention provides a fault diagnosis and fault-tolerant reconfiguration system for a satellite attitude control system, comprising:

[0038] The first unit is used to collect attitude angle data, angular velocity data, reaction wheel speed data, and torque data of the satellite attitude control system, and establish a Kalman filter state observer; embed the Kalman filter state observer into a fault diagnosis model, use the fault diagnosis model to generate a residual sequence, and perform adaptive dynamic threshold calculation on the residual sequence to obtain a fault feature vector; input the fault feature vector into a pre-trained deep belief network in the fault diagnosis model to obtain a fault diagnosis result, which includes the fault type and fault degree;

[0039] The second unit is used to determine the faulty component based on the fault diagnosis result, establish a mixed integer programming model based on the fault type and fault degree of the faulty component, set attitude control accuracy constraints and actuator power constraints in the mixed integer programming model, and solve the mixed integer programming model using an improved genetic algorithm to obtain a fault-tolerant control strategy, which includes an actuator reconfiguration scheme and controller parameters.

[0040] The third unit is used to reconfigure the actuator operating mode in the satellite attitude control system according to the actuator reconfiguration scheme in the fault-tolerant control strategy; update the controller of the satellite attitude control system using the controller parameters in the fault-tolerant control strategy; control the satellite attitude based on the updated controller to determine the satellite attitude error; and when the satellite attitude error exceeds a preset threshold, re-diagnose the fault using the fault diagnosis model.

[0041] Third aspect of the present invention

[0042] An electronic device is provided, comprising:

[0043] processor;

[0044] Memory used to store processor-executable instructions;

[0045] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0046] Fourth aspect of the embodiments of the present invention,

[0047] A computer-readable storage medium is provided, having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0048] The beneficial effects of this application are as follows:

[0049] 1. Improve fault diagnosis accuracy: By combining Kalman filter state observer and deep belief network, the fault type and degree can be identified more accurately, avoiding misdiagnosis and missed diagnosis.

[0050] 2. Achieve fault-tolerant reconfiguration of attitude control system: Based on fault diagnosis results, through mixed integer programming model and improved genetic algorithm, actuator reconfiguration scheme and controller parameters can be optimized to achieve fault-tolerant control of attitude control system and ensure satellite attitude stability.

[0051] 3. Enhance the robustness of the attitude control system: Through adaptive dynamic threshold calculation and closed-loop fault diagnosis process, the satellite attitude error can be monitored in real time, and the control strategy can be adjusted according to the actual situation, thereby improving the robustness and adaptability of the attitude control system. Attached Figure Description

[0052] Figure 1 This is a flowchart illustrating the fault diagnosis and fault-tolerant reconfiguration method for a satellite attitude control system according to an embodiment of the present invention.

[0053] Figure 2 This is a schematic diagram of the fault diagnosis and fault-tolerant reconfiguration system for the satellite attitude control system according to an embodiment of the present invention. Detailed Implementation

[0054] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0055] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0056] Figure 1This is a flowchart illustrating the fault diagnosis and fault-tolerant reconfiguration method for a satellite attitude control system according to an embodiment of the present invention. Figure 1 As shown, the method includes:

[0057] S101. Collect attitude angle data, angular velocity data, reaction wheel speed data, and torque data of the satellite attitude control system, and establish a Kalman filter state observer; embed the Kalman filter state observer into the fault diagnosis model, use the fault diagnosis model to generate a residual sequence, and perform adaptive dynamic threshold calculation on the residual sequence to obtain a fault feature vector; input the fault feature vector into a pre-trained deep belief network in the fault diagnosis model to obtain a fault diagnosis result, which includes the fault type and fault degree;

[0058] S102. Based on the fault diagnosis results, determine the faulty component, and establish a mixed integer programming model based on the fault type and fault degree of the faulty component; set attitude control accuracy constraints and actuator power constraints in the mixed integer programming model; use an improved genetic algorithm to solve the mixed integer programming model to obtain a fault-tolerant control strategy, which includes an actuator reconfiguration scheme and controller parameters;

[0059] S103. Reconfigure the actuator operating mode in the satellite attitude control system according to the actuator reconfiguration scheme in the fault-tolerant control strategy; update the controller of the satellite attitude control system using the controller parameters in the fault-tolerant control strategy; control the satellite attitude based on the updated controller to determine the satellite attitude error; when the satellite attitude error exceeds a preset threshold, re-diagnose the fault using the fault diagnosis model.

[0060] In one alternative implementation, a Kalman filter state observer is established, including:

[0061] The attitude angle data is used as the attitude quaternion observation, and the angular velocity data is used as the angular velocity observation. A state vector containing the attitude quaternion observation and the angular velocity observation is constructed. A satellite attitude dynamics model is established. According to the Euler motion equation, the reaction wheel speed data and the torque data are substituted into the satellite attitude dynamics model to obtain a state space expression. The state space expression is discretized within the sampling time interval to obtain a discrete state equation.

[0062] An extended Kalman filter state observer is constructed. Based on the discrete state equations, the Jacobian matrix is ​​calculated. The Jacobian matrix is ​​linearized at the current state estimate to obtain the state transition matrix. The observation matrix is ​​calculated based on the observation equations. The system noise covariance matrix and observation noise covariance matrix of the extended Kalman filter state observer are initialized. Based on the initialized system noise covariance matrix and observation noise covariance matrix, a one-step state prediction value and an error covariance prediction value are calculated. The Kalman gain matrix is ​​calculated using the one-step state prediction value and the error covariance prediction value. The Kalman gain matrix is ​​multiplied by the observation residual to obtain the state estimate correction. The state estimate correction is added to the one-step state prediction value to obtain the state estimate.

[0063] A new error covariance matrix is ​​calculated based on the state estimate. The state estimation accuracy is calculated based on the error covariance matrix. When the state estimation accuracy meets the preset accuracy requirement, the state estimate is output. The observation residual sequence is calculated. The actual covariance of the observation residual sequence is compared with the theoretical covariance to obtain an adaptive factor. The system noise covariance matrix and the observation noise covariance matrix are adjusted using the adaptive factor. The adjusted system noise covariance matrix and the observation noise covariance matrix are substituted into the extended Kalman filter state observer to obtain the Kalman filter state observer.

[0064] First, acquire the satellite's attitude quaternion observations, angular velocity observations, reaction wheel rotation speed data, and torque data. The attitude quaternion observations describe the satellite's attitude direction, the angular velocity observations describe the satellite's rotational rate, the reaction wheel rotation speed data describes the rotational speed of the reaction wheel, and the torque data describes the external torque acting on the satellite. For example, the attitude quaternion observations are (0.1, 0.2, 0.3, 0.9), the angular velocity observations are (0.01, 0.02, 0.03) radians per second, the reaction wheel rotation speed data are (1000, 1200, 1500) revolutions per minute, and the torque data are (0.001, 0.002, 0.003) Newton-meters.

[0065] Next, the state vector is constructed. The attitude quaternion observations and angular velocity observations are combined into a state vector, which fully describes the satellite's attitude and motion state. For example, the state vector can be represented as [0.1,0.2,0.3,0.9,0.01,0.02,0.03].

[0066] Next, a satellite attitude dynamics model is established. Based on Euler's equations of motion, and combined with the reaction wheel speed and torque data, the change in satellite attitude over time is described. This can be understood as predicting the satellite's attitude at the next moment based on the current attitude, the applied torque, and the reaction wheel speed.

[0067] The established satellite attitude dynamics model is discretized. Since computer data processing is discrete, the continuous dynamics model needs to be transformed into discrete state equations to facilitate computation. This process is similar to converting a smooth curve into a series of discrete points. For example, assuming a sampling time interval of 0.01 seconds, the dynamics model can be discretized into a series of state transition relationships with 0.01-second intervals.

[0068] An Extended Kalman Filter (EKF) state observer is constructed. The EKF is used to estimate the satellite's attitude. It continuously refines the attitude estimate through two steps: prediction and update.

[0069] The state transition matrix is ​​obtained by calculating and linearizing the Jacobian matrix. The Jacobian matrix describes the interaction between each element in the state vector. By linearizing the Jacobian matrix at the current state estimate, the state transition matrix is ​​obtained, which describes the transition relationship of the state vector from the current time step to the next time step.

[0070] Calculate the observation matrix. The observation matrix describes the relationship between observations and state vectors. For example, attitude quaternion observations and angular velocity observations are both components of the state vector, so the observation matrix can be used to transform the state vector into the corresponding observations.

[0071] Initialize the system noise covariance matrix and the observation noise covariance matrix. These two matrices describe the errors of the system model and the observation data, respectively. For example, the system noise covariance matrix can be initialized as an identity matrix multiplied by a small coefficient, and the observation noise covariance matrix can be initialized as a diagonal matrix, with the diagonal elements being the squares of the measurement errors of each observation.

[0072] Calculate the one-step state prediction and the error covariance prediction. Based on the current state estimate and the system model, predict the state value at the next time step and calculate the covariance matrix of the prediction error.

[0073] Calculate the Kalman gain matrix. The Kalman gain matrix is ​​used to balance the influence of predicted and observed values ​​on the final state estimate.

[0074] The state estimate correction is calculated and added to the one-step state prediction to obtain the state estimate. The state estimate correction is calculated using the Kalman gain matrix and the observation residual (the difference between the observed and predicted observations), and then added to the one-step state prediction to obtain the final state estimate.

[0075] Calculate the new error covariance matrix and the state estimation accuracy. Update the error covariance matrix based on the state estimates, and calculate the state estimation accuracy based on the error covariance matrix.

[0076] Determine whether the state estimation accuracy meets the preset accuracy requirements. If it does, output the state estimate. For example, if the state estimation accuracy is higher than a preset threshold, the attitude estimation result is considered reliable, and the state estimate is output.

[0077] The actual covariance and theoretical covariance of the observed residual sequence are calculated to obtain the adaptive factor. The observed residual sequence reflects the difference between the observed values ​​and the predicted observed values. By comparing the actual covariance and the theoretical covariance, the adaptive factor can be obtained, which is used to adjust the system noise covariance matrix and the observed noise covariance matrix.

[0078] The system noise covariance matrix and the observation noise covariance matrix are adjusted using adaptive factors, and the adjusted matrices are substituted into the extended Kalman filter state observer to obtain the final Kalman filter state observer.

[0079] Beneficial effects:

[0080] The accuracy of attitude estimation is improved. By fusing attitude quaternion and angular velocity information and using an extended Kalman filter for optimal estimation, the influence of measurement noise can be effectively reduced, thus improving the accuracy of attitude estimation.

[0081] The adaptive capability of the attitude estimation system is enhanced. By introducing an adaptive factor, the system noise covariance matrix and the observation noise covariance matrix can be adjusted in real time according to the observation residual sequence, enabling the filter to adapt to different environments and operating conditions, thereby improving the robustness and reliability of the system.

[0082] This method reduces computational complexity. Compared to other attitude estimation methods, it requires less computation, is easy to implement on spaceborne computers, and can meet real-time requirements.

[0083] In one alternative implementation, a satellite attitude dynamics model is established, including:

[0084] Quaternions are used to describe the satellite attitude. A quaternion vector containing scalar and vector parts is constructed. The attitude kinematic equation is established by associating the quaternion vector with the angular velocity vector. The derivative of the quaternion vector is represented by the product of the quaternion correlation matrix and the angular velocity vector. The attitude dynamic equation is established based on the Euler rigid body dynamics principle. The satellite rotational inertia matrix, the cross product of the angular velocity vector, the control torque, and the disturbance torque are substituted into the attitude dynamic equation. The control torque is generated by the reaction wheel system.

[0085] A dynamic equation for the reaction wheel is established by substituting the angular momentum vector, driving torque, and cross product of the angular velocity vector and angular momentum vector of the reaction wheel into the dynamic equation. The attitude kinematics equation, the attitude dynamic equation, and the reaction wheel dynamic equation are combined to form a complete satellite attitude dynamic model, which includes quaternion derivative equations, angular velocity derivative equations, and reaction wheel angular momentum derivative equations. The satellite attitude dynamic model is discretized using the fourth-order Runge-Kutta numerical integration method to obtain discrete state equations containing quaternion state update functions, angular velocity state update functions, and reaction wheel angular momentum state update functions.

[0086] A process noise term is added to the discrete state equation, and the quaternion vector, angular velocity vector, and reaction wheel angular momentum vector are combined to form the state vector. The reaction wheel driving torque is used as the control input to establish a state-space expression containing process noise. An observation equation is established based on the attitude quaternion measurement value of the star sensor, the angular velocity measurement value of the gyroscope, and the angular momentum measurement value of the reaction wheel speed sensor. Measurement noise is added to the observation equation, and the state-space expression and the observation equation are combined to form a satellite attitude dynamics model.

[0087] First, we define the representation of satellite attitude. Quaternions are used to describe satellite attitude; a quaternion contains a scalar part and a three-dimensional vector part. To describe attitude changes, we need to relate the quaternion to the satellite's angular velocity. Here, the rate of change (derivative) of the quaternion can be represented by the product of a quaternion correlation matrix and the angular velocity vector, establishing the attitude kinematic equations.

[0088] Next, the attitude dynamics equations are established. According to Euler's rigid body dynamics principles, the satellite's attitude changes are caused by torques acting on the satellite. These torques include control torques and disturbance torques. The control torque is generated by the satellite's reaction wheel system and is used to control the satellite's attitude. The disturbance torque is caused by environmental factors, such as atmospheric drag and solar radiation pressure. The attitude dynamics equations relate the satellite's moment of inertia matrix, the cross product of angular velocities, the control torque, and the disturbance torque.

[0089] Then, the dynamic equations of the reaction wheel are established. The change in angular momentum of the reaction wheel is determined by the cross product of the driving torque and angular velocity of the reaction wheel with its angular momentum. The driving torque of the reaction wheel is the control input, used to change the rotational speed of the reaction wheel, thereby changing the attitude of the satellite.

[0090] Combining the attitude kinematics equations, attitude dynamics equations, and reaction wheel dynamics equations forms a complete satellite attitude dynamics model. This model contains three differential equations: the quaternion derivative equation, the angular velocity derivative equation, and the reaction wheel angular momentum derivative equation.

[0091] To perform simulation and control on a computer, the continuous dynamic model needs to be discretized. Here, a fourth-order Runge-Kutta numerical integration method is used to discretize the dynamic model. This method yields discretized state update functions, including quaternion state update functions, angular velocity state update functions, and reaction wheel angular momentum state update functions. The discrete state equations describe how the satellite's state changes at each time step.

[0092] To more realistically simulate satellite motion, process noise needs to be incorporated into the discrete state equations. Process noise represents factors not considered in the model, such as model errors and unmodeled disturbance forces. By combining quaternions, angular velocities, and reaction wheel angular momentum into a state vector, and using the reaction wheel driving torque as the control input, a state-space expression incorporating process noise is established.

[0093] To estimate a satellite's attitude and angular velocity, sensors are needed for measurement. Star sensors measure the satellite's attitude quaternions, gyroscopes measure its angular velocity, and reaction wheel speed sensors measure the angular momentum of the reaction wheel. Because sensor measurements have inherent errors, measurement noise needs to be added to the observation equations. Combining the state-space expression and the observation equations yields a complete satellite attitude dynamics model.

[0094] Example: Assume the satellite's moment of inertia matrix is ​​a diagonal matrix, with each diagonal element being 10 kg·m. 2 12kg·m 2 and 15kg·m 2 The initial quaternion is [1,0,0,0], the initial angular velocity is [0.1 rad / s, 0.2 rad / s, 0.3 rad / s], and the initial angular momentum of the reaction wheel is [0.01 N·m·s, 0.02 N·m·s, 0.03 N·m·s]. The driving torque of the reaction wheel is assumed to be [0.001 N·m, 0.002 N·m, 0.003 N·m]. Numerical integration is performed using the fourth-order Runge-Kutta method with a time step of 0.01 seconds. The quaternion, angular velocity, and angular momentum of the reaction wheel at the next moment can be calculated.

[0095] Beneficial effects:

[0096] I. High model accuracy: The model uses quaternions to describe attitude, avoiding singularity problems, and combines Euler dynamics equations and reaction wheel dynamics equations to accurately describe the satellite's attitude motion.

[0097] Second, the model has strong integrity: the model includes attitude kinematic equations, attitude dynamic equations and reaction wheel dynamic equations, forming a complete attitude dynamic model that can be used for the design and simulation of satellite attitude control systems.

[0098] Third, it is highly practical: The model is discretized using the fourth-order Runge-Kutta method, which is easy to implement on a computer. It also incorporates process noise and measurement noise, making it closer to the actual situation. It can be used to verify satellite attitude determination and control algorithms.

[0099] In one optional implementation, the fault feature vector is input into a pre-trained deep belief network in the fault diagnosis model to obtain a fault diagnosis result, including:

[0100] The fault feature vector is mapped to the zero-to-one interval using the minimum-maximum normalization method to obtain a normalized feature vector. The covariance matrix is ​​calculated based on the normalized feature vector, and the covariance matrix is ​​decomposed into principal components. Principal components with a cumulative contribution rate greater than a preset threshold are selected to construct a dimensionality reduction mapping matrix. The normalized feature vector is then subjected to dimensionality reduction processing through the dimensionality reduction mapping matrix to obtain a dimensionality reduction feature vector. The dimensionality reduction feature vector is input into the visible layer of a deep belief network, wherein the deep belief network contains multiple stacked restricted Boltzmann machines, and the number of neurons in each layer of the restricted Boltzmann machine decreases.

[0101] The dimensionality-reduced feature vector is processed layer by layer through the Restricted Boltzmann Machine (RBM) of the deep belief network. The product of the RBM input vector and the weight matrix, and the sum of the bias vectors are calculated. The sum of the bias vectors is activated by a normalization function to obtain the hidden layer output. The neurons in the hidden layer output are set to zero according to a preset random deactivation probability. For the hidden layer outputs that are not set to zero, numerical compensation is performed according to a scaling factor. The hidden layer output after random deactivation is used as the visual layer input of the next layer of the RBM. The feature transformation, normalization activation, random deactivation and numerical compensation process is repeated until the hidden layer output of the last layer of the RBM is obtained. The random deactivation probability and the scaling factor are determined according to the number of neurons in each layer.

[0102] The hidden layer output of the last layer of the restricted Boltzmann machine is connected to a normalized classifier. The original score of each fault type is calculated based on the product of the hidden layer output and the classifier weight matrix and the bias vector. The original score is then exponentially calculated to obtain a non-negative value. This non-negative value is divided by the sum of the exponents of all category scores to obtain the conditional probability distribution of different fault types. The category with the highest conditional probability that is greater than a preset probability threshold is selected as the fault type. The activation values ​​of the output layer of the normalized classifier are weighted and combined to obtain the fault severity. The fault type and the fault severity are combined to generate a fault diagnosis result.

[0103] First, collect the equipment's operational data, which may include sensor readings, system logs, and so on. For example, the operational data of a rotating machine may include vibration signals, temperature, current, and so on.

[0104] Then, the collected operational data undergoes preprocessing, such as noise removal and outlier handling, to improve data quality. For example, wavelet transform can be used to remove high-frequency noise from vibration signals. The preprocessed data is then used to extract fault features.

[0105] Next, fault features are extracted from the preprocessed data. These features reflect the health status of the equipment. For example, frequency domain features, such as peak frequency and root mean square value, can be extracted from vibration signals. The extracted features form a feature vector. Suppose five features are extracted: feature A, feature B, feature C, feature D, and feature E, with values ​​of 10, 20, 30, 40, and 50, respectively. Then the fault feature vector is [10, 20, 30, 40, 50].

[0106] To eliminate the influence of different dimensions between features, the extracted fault feature vectors are normalized. For example, the minimum-maximum normalization method is used to map the feature values ​​to the interval between zero and one. Assuming that the minimum value of feature A is 5 and the maximum value is 15, the normalized value of feature A is (10-5) / (15-5) = 0.5. Other features are also normalized in the same way to obtain normalized feature vectors, such as [0.5, 0.6, 0.7, 0.8, 0.9].

[0107] The covariance matrix is ​​calculated based on normalized eigenvectors, reflecting the correlation between different features. Then, eigenvalue decomposition is performed on the covariance matrix to obtain principal components. Principal components with a cumulative contribution rate greater than a preset threshold are selected, for example, those with a cumulative contribution rate greater than 85%. These principal components constitute a dimensionality reduction mapping matrix.

[0108] The normalized eigenvectors are reduced in dimensionality using a dimensionality reduction mapping matrix to obtain dimensionality-reduced eigenvectors. For example, assuming the dimensionality reduction mapping matrix reduces a five-dimensional eigenvector to three dimensions, the dimensionality-reduced eigenvectors might be [0.3, 0.5, 0.7].

[0109] The dimensionality-reduced feature vectors are input into a pre-trained deep belief network for fault diagnosis. The deep belief network consists of multiple stacked Restricted Boltzmann Machines (RBMs), with the number of neurons decreasing in each layer. For example, the first layer has 100 neurons, the second layer has 50 neurons, and the third layer has 25 neurons.

[0110] The dimensionality-reduced feature vectors are first input into the first layer of the deep belief network, which is the visible layer. The output of the visible layer is then processed by calculating the weight matrix and bias vector, and then activated by a normalization function to obtain the hidden layer output.

[0111] To prevent overfitting, the hidden layer outputs are randomly deactivated. A preset random deactivation probability, such as 20%, is used to determine whether to zero out neurons in the hidden layer outputs. For hidden layer outputs that are not zeroed out, numerical compensation is performed using a scaling factor, such as 1.25.

[0112] The hidden layer output after random deactivation is used as the visual layer input of the next restricted Boltzmann machine. The process of feature transformation, normalized activation, random deactivation and numerical compensation is repeated until the hidden layer output of the last restricted Boltzmann machine is obtained.

[0113] The hidden output of the last layer of the Restricted Boltzmann Machine is connected to the normalized classifier. The raw score for each fault type is calculated based on the product of the hidden output and the classifier weight matrix, and the bias vector.

[0114] The original scores are exponentially calculated to obtain non-negative values. These non-negative values ​​are then divided by the sum of the exponents of all category scores to obtain the conditional probability distributions for different fault types. The category with the highest conditional probability that is greater than a preset probability threshold is selected as the fault type. For example, assuming there are three fault types, A, B, and C, with conditional probabilities of 0.7, 0.2, and 0.1 respectively, and a preset probability threshold of 0.6, then the fault type is A.

[0115] The activation values ​​of the normalized classifier output layer are weighted and combined to obtain the fault severity. For example, a fault severity of 0.8 indicates a relatively severe fault. Finally, the fault type and fault severity are combined to generate a fault diagnosis result, such as "Fault type A, fault severity 0.8".

[0116] Beneficial effects:

[0117] Improved accuracy of fault diagnosis: Through the deep feature learning capability of deep belief networks, fault features can be effectively extracted, thereby improving the accuracy of fault diagnosis.

[0118] Enhanced robustness of fault diagnosis: By employing methods such as random deactivation and principal component analysis, the impact of noise and outliers can be effectively reduced, thereby enhancing the robustness of fault diagnosis.

[0119] Improved efficiency of fault diagnosis: By using dimensionality reduction and pre-trained models, the amount of computation can be reduced, thus improving the efficiency of fault diagnosis.

[0120] In one optional implementation, the faulty component is determined based on the fault diagnosis result, and a mixed-integer programming model is established based on the fault type and fault severity of the faulty component, including:

[0121] Based on the fault diagnosis results, fault type, fault severity, and fault location information are obtained. A fault impact assessment matrix is ​​established based on the fault type and fault severity. The fault impact assessment matrix includes the impact coefficient of the faulty component on attitude control accuracy and the impact coefficient of the faulty component on actuator power. The fault impact assessment matrix is ​​input into a mixed integer programming model, and attitude control variables are established in the mixed integer programming model. The attitude control variables include pitch angle, yaw angle, and roll angle. The product of the attitude control variables and the fault impact assessment matrix is ​​used as the predicted value of attitude control accuracy under fault conditions.

[0122] An actuator power demand model is established based on the dynamic characteristics of the faulty component. The product of the actuator power demand model and the fault impact assessment matrix is ​​used as the actuator power prediction value under the fault state. Actuator power constraints are set in the mixed integer programming model. An attitude control accuracy threshold is determined according to the equipment attitude control accuracy requirements. An attitude control accuracy constraint is established by the difference between the attitude control accuracy prediction value and the attitude control accuracy threshold. The attitude control accuracy constraint is set in the mixed integer programming model.

[0123] An objective function is established based on the fault impact assessment matrix. The objective function includes a term for minimizing attitude control deviation, a term for minimizing actuator power consumption, and a term for fault impact compensation. A mixed integer programming model is constructed by combining the objective function, the attitude control variables, the actuator power constraints, and the attitude control accuracy constraints. The attitude control deviation minimization term is calculated based on the deviation between the attitude control variables and the desired attitude. The actuator power consumption minimization term is calculated based on the predicted actuator power value. The fault impact compensation term is calculated based on the impact coefficients in the fault impact assessment matrix.

[0124] This method utilizes real-time fault information provided by a fault diagnosis system to dynamically adjust the attitude control strategy, thereby achieving attitude stability and mission execution of the spacecraft under fault conditions.

[0125] First, fault diagnosis information is acquired. The fault diagnosis system continuously monitors the status of various components of the spacecraft and outputs diagnostic results. These results include the specific location information of the faulty component, the type of fault (e.g., sensor drift, actuator jamming), and the severity of the fault (e.g., the magnitude of the drift, the range of the jamming angle). For example, the diagnostic system might output, "The attitude control system's X-axis reaction wheel is partially jammed, with a jamming angle range of ±10 degrees."

[0126] Next, a fault impact assessment matrix is ​​constructed. Based on the acquired fault type and severity, a matrix is ​​established to quantify the impact of the fault on the spacecraft's attitude control performance. This matrix contains two types of impact coefficients: impact coefficients on attitude control accuracy, such as the impact on pitch, yaw, and roll angles; and impact coefficients on actuator power, such as the impact on reaction wheel and thruster power. Taking a partial jamming of the X-axis reaction wheel as an example, the row corresponding to the X-axis reaction wheel in the fault impact assessment matrix might show a 20% decrease in roll angle control accuracy in the attitude control accuracy impact coefficient, while the actuator power impact coefficient might show a 15% increase in X-axis reaction wheel power consumption.

[0127] Then, an attitude control accuracy prediction model is established. By multiplying the attitude control variables (pitch angle, yaw angle, roll angle) with the fault impact assessment matrix, the attitude control accuracy under fault conditions can be predicted. For example, if the desired roll angle is 0 degrees, but the X-axis reaction wheel partially jams, causing a 20% decrease in roll angle control accuracy, the predicted roll angle may fluctuate between ±0.2 degrees.

[0128] Simultaneously, an actuator power demand prediction model is established. Based on the dynamic characteristics of the faulty component, an actuator power demand model is created. Multiplying this model by the fault impact assessment matrix allows for the prediction of the actuator power demand under fault conditions. For example, if the X-axis reaction wheel partially jams, causing its power consumption to increase by 15%, then the predicted X-axis reaction wheel power demand will increase by 15% compared to the normal demand.

[0129] Next, a mixed-integer programming model is constructed. The objective function of this model comprises three parts: a term for minimizing attitude control deviation, i.e., minimizing the deviation between the actual and desired attitude; a term for minimizing actuator power consumption, i.e., minimizing actuator power consumption; and a term for compensating for fault effects, i.e., compensating for the negative impacts of faults. The model's constraints include actuator power constraints, such as the maximum output torque of the reaction wheel, and attitude control accuracy constraints, such as the roll angle accuracy must be controlled within ±0.5 degrees.

[0130] Finally, the mixed-integer programming model is solved. An optimization algorithm is used to solve this model to obtain the optimal attitude control strategy, which minimizes the objective function while satisfying all constraints. For example, the solution might indicate increasing the output of the Y-axis and Z-axis reaction wheels to compensate for the effects of a faulty X-axis reaction wheel, while adjusting the thruster output to maintain spacecraft attitude stability.

[0131] Beneficial effects:

[0132] 1. Improved robustness of spacecraft attitude control system. This method can dynamically adjust attitude control strategy based on fault diagnosis results, effectively responding to various fault conditions and preventing spacecraft instability.

[0133] 2. Optimized spacecraft energy management. This method minimizes actuator power consumption and extends the spacecraft's on-orbit lifespan while ensuring attitude control accuracy.

[0134] 3. Enhanced spacecraft mission execution capabilities. Even in the event of a malfunction, this method can ensure that the spacecraft completes its intended mission as much as possible.

[0135] In one optional implementation, an improved genetic algorithm is used to solve the mixed-integer programming model to obtain a fault-tolerant control strategy. The fault-tolerant control strategy includes an actuator reconfiguration scheme and controller parameters, including:

[0136] In a mixed-integer programming model, actuator reconfiguration variables and controller parameter variables are encoded as chromosomes. The actuator reconfiguration variables include actuator operating status flags, actuator output allocation coefficients, and actuator switching timing. The controller parameter variables include proportional coefficient matrices, integral coefficient matrices, and differential coefficient matrices. A multi-objective fitness function is established based on the actuator reconfiguration variables and controller parameter variables. The multi-objective fitness function includes actuator resource utilization evaluation items, control performance evaluation items, fault tolerance evaluation items, and energy consumption evaluation items.

[0137] The chromosome population is initialized according to the multi-objective fitness function; the chromosome population is stratified to construct a population stratification structure, an elite population library is established in the population stratification structure, the fitness value of the chromosomes is calculated according to the multi-objective fitness function, and non-dominated solution chromosomes are stored in the elite population library; an adaptive crossover operator is used to perform a crossover operation on the chromosome population, the crossover probability of the adaptive crossover operator is determined based on the chromosome distribution in the elite population library and the multi-objective fitness function value, and the crossover chromosomes are generated and added to the chromosome population;

[0138] A hybrid mutation operator is used to mutate the crossover chromosome population to obtain an elite population pool. Cluster analysis is performed on the chromosomes in the elite population pool, and a clustering algorithm is used to group the chromosomes. From each group, the chromosome with the optimal multi-objective fitness function value is selected to form an optimal solution set. The chromosomes in the optimal solution set are decoded to obtain a fault-tolerant control strategy. The fault-tolerant control strategy includes an actuator operating mode switching timing table, an actuator output allocation matrix, an actuator cooperative control scheme, a proportional coefficient matrix, an integral coefficient matrix, and a differential coefficient matrix.

[0139] Fault-tolerant control strategies aim to enable the system to maintain stable operation even when a fault occurs. This implementation proposes a fault-tolerant control strategy generation method based on an improved genetic algorithm. This method uses an improved genetic algorithm to solve a mixed integer programming model that includes actuator reconfiguration and controller parameter optimization, and finally obtains an optimized fault-tolerant control strategy.

[0140] First, the variables in the mixed-integer programming model are encoded to form chromosomes in the genetic algorithm. A chromosome consists of two parts: executor reconfiguration variables and controller parameter variables. Executor reconfiguration variables include executor status flags (indicating whether the executor is working), executor output allocation coefficients (determining the output percentage of each executor), and executor switching timing (specifying the timing of executor switching). Controller parameter variables include proportional coefficient matrices, integral coefficient matrices, and differential coefficient matrices; these parameters determine the controller's performance. For example, in a system with three executors, the executor status flags can be represented by three binary digits, such as "101" indicating that the first and third executors are working, and the second executor is not working. The output allocation coefficients can be represented as a vector, such as "[0.6,0,0.4]" indicating that the first executor has a 60% output percentage, the second executor does not output, and the third executor has a 40% output percentage. The executor switching timing can be represented as a series of time points, such as "[10s,20s]" indicating that the executor switches at 10 seconds and 20 seconds. The controller parameter matrix is ​​determined based on the system characteristics. For example, the proportional coefficient matrix can be "[2,0,1; 0,3,0; 1,0,4]".

[0141] Next, a multi-objective fitness function is established to evaluate the quality of chromosomes. This function includes four evaluation terms: actuator resource utilization (e.g., fewer working actuators mean higher utilization), control performance (e.g., smaller overshoot and settling time mean better performance), fault tolerance (e.g., shorter system recovery time after a fault means better fault tolerance), and energy consumption (e.g., lower actuator power consumption means lower energy consumption). Each evaluation term is assigned a weight according to specific system requirements. For example, for systems with high resource utilization requirements, the resource utilization evaluation term can be given a higher weight. Assuming a chromosome has an actuator resource utilization score of 0.8, a control performance score of 0.9, a fault tolerance score of 0.7, and an energy consumption score of 0.6, and the weights of the four evaluation terms are 0.2, 0.3, 0.4, and 0.1 respectively, then the fitness value of this chromosome is 0.2*0.8 + 0.3*0.9 + 0.4*0.7 + 0.1*0.6 = 0.76.

[0142] Then, the chromosome population is initialized and stratified to construct a population hierarchy and establish an elite population pool. The fitness value of each chromosome is calculated, and non-dominated solution chromosomes (i.e., solutions that are not inferior to other chromosomes on any objective) are stored in the elite population pool. For example, if the initial population contains 100 chromosomes, after fitness value calculation and non-dominated solution screening, 10 chromosomes may enter the elite population pool.

[0143] Subsequently, an adaptive crossover operator is used to perform a crossover operation on the chromosome population. The crossover probability is determined based on the chromosome distribution in the elite population pool and the multi-objective fitness function value. For example, if chromosomes with high fitness values ​​are concentrated in a certain region in the elite population pool, the crossover probability in that region will be higher. The crossover chromosomes are then added to the population. For example, if the two parent chromosomes are "[1,0,1]" and "[0,1,1]", and the crossover position is the second position, then the two generated offspring chromosomes will be "[1,1,1]" and "[0,0,1]".

[0144] Next, a hybrid mutation operator is used to mutate the crossover population, generating a new population and updating the elite population pool. For example, if a chromosome is "[1,0,1]" and the mutation bit is the first position, then the mutated chromosome will be "[0,0,1]".

[0145] Cluster analysis is performed on the chromosomes in the elite population database to group them into groups. For example, the k-means algorithm is used to divide the chromosomes in the elite population database into 5 groups. The chromosome with the best multi-objective fitness function value is selected from each group to form the optimal solution set.

[0146] Finally, the chromosomes in the optimal solution set are decoded to obtain the final fault-tolerant control strategy, including the actuator operating mode switching timing table, actuator output allocation matrix, actuator cooperative control scheme, proportional coefficient matrix, integral coefficient matrix, and differential coefficient matrix. For example, if the actuator output allocation matrix obtained after decoding a chromosome is "[0.6,0,0.4]" and the proportional coefficient matrix is ​​"[2,0,1;0,3,0;1,0,4]", it means that the first actuator's output accounts for 60%, the third actuator's output accounts for 40%, and the controller uses the corresponding proportional coefficient matrix for control.

[0147] Beneficial effects:

[0148] 1. Improve the reliability of the control system: By reconfiguring actuators and optimizing controller parameters, the system can maintain stable operation even if some actuators fail, thus improving the system's reliability and fault tolerance.

[0149] 2. Optimize resource utilization: This method can dynamically adjust the working status and output allocation of the actuator according to the system's operating status and fault conditions, thereby avoiding resource waste and improving the actuator utilization rate.

[0150] 3. Improve control performance: By optimizing controller parameters, the dynamic and steady-state performance of the system can be improved, such as reducing overshoot and shortening settling time, thereby improving the overall performance of the control system.

[0151] In one optional implementation, a hybrid mutation operator is used to mutate the crossover chromosome population to obtain an elite population pool, including:

[0152] A hybrid mutation operator is used to mutate the crossover chromosome population. This hybrid mutation operator includes discrete and continuous mutation operations. The discrete mutation operation changes the values ​​of the actuator's operating status flags by bit reversal, adjusts the distribution of the actuator's output allocation coefficients by random permutation, and modifies the timing of the actuator switching sequence by binary encoding conversion. The continuous mutation operation adjusts the parameters of the proportional coefficient matrix by Gaussian distribution perturbation, the parameters of the integral coefficient matrix by uniform distribution perturbation, and the parameters of the differential coefficient matrix by non-uniform distribution perturbation, thus obtaining the mutated chromosome population.

[0153] The mutated chromosome population is evaluated based on the multi-objective fitness function, and chromosomes with better evaluation results are selected into the elite population pool. The simulated annealing algorithm is used to perform local search optimization on the chromosomes in the elite population pool. The initial temperature of the simulated annealing algorithm is determined according to the value range of the multi-objective fitness function, the cooling coefficient of the simulated annealing algorithm is calculated in an adaptive manner, and the termination temperature of the simulated annealing algorithm is dynamically adjusted based on the number of iterations.

[0154] During the local search process of the simulated annealing algorithm, a discrete neighborhood structure is adopted for the actuator reconstruction variables, which includes bit-flipped neighborhoods, swapped neighborhoods, and insertion neighborhoods; a continuous neighborhood structure is adopted for the controller parameter variables, which includes Gaussian neighborhoods, uniform neighborhoods, and non-uniform neighborhoods; the acceptance probability of a new solution is calculated based on the multi-objective fitness function, and the acceptance probability decreases as the temperature decreases during the iteration process; the chromosome optimized by the local search is then updated to the elite population library.

[0155] A hybrid method for optimizing control system parameters and actuator configuration, aimed at improving the performance of the control system, is proposed, which is mainly achieved through the global search capability of genetic algorithm and the local search capability of simulated annealing algorithm.

[0156] First, the chromosome population is initialized. Each chromosome represents a potential configuration scheme for the control system, encoding the actuator's operating state, output allocation coefficients, switching timing, and controller parameters (proportional, integral, and derivative coefficient matrices). For example, assuming the system has three actuators, the actuator operating state flag can be represented by a 3-bit binary number, where 1 indicates operation and 0 indicates inactivity. "101" represents that the first and third actuators are operating, and the second actuator is not operating. The output allocation coefficients can be represented by an array, such as [0.4, 0, 0.6], indicating that the first actuator allocates 40% of the output, the second actuator is not operating, and the third actuator allocates 60% of the output. The switching timing can be represented by an array, such as [10, 20, 30], indicating that the three actuators switch states at the 10th, 20th, and 30th time units, respectively. The controller parameter matrix can be represented by a multi-dimensional array, for example, the proportional coefficient matrix is ​​[[0.1, 0.2], [0.3, 0.4]]. The initial population size is set to 100 chromosomes, and the parameter values ​​of each chromosome are randomly generated within a preset range.

[0157] Next, a crossover operation is performed. A two-point crossover is used, randomly selecting two chromosomes from the population as parents. Two crossover points are randomly selected on the chromosomes, and gene segments between the two crossover points are exchanged to generate two new offspring chromosomes. For example, if the actuator switching sequence of parent chromosome 1 is [10,20,30], and the actuator switching sequence of parent chromosome 2 is [15,25,35], and the crossover points are the second and third genes, then the actuator switching sequence of the generated offspring chromosome 1 will be [10,25,35], and the actuator switching sequence of offspring chromosome 2 will be [15,20,30].

[0158] Then, mutation operations are performed on the crossover chromosome population. A hybrid mutation operator is used, including discrete and continuous mutations. For discrete variables, such as actuator operating status flags, bit reversal is used for mutation, i.e., changing 0 to 1 or 1 to 0. For output allocation coefficients, random permutation is used, i.e., randomly swapping the positions of elements in the allocation coefficient array. For actuator switching timing, binary encoding conversion is used for mutation, i.e., converting time points to binary numbers, randomly selecting one bit for reversal, and then converting back to decimal numbers. For continuous variables, such as controller parameters, different types of perturbation methods are used for mutation. The proportional coefficient matrix uses Gaussian distribution perturbation, adding a random number following a Gaussian distribution to the original values. The integral coefficient matrix uses uniform distribution perturbation, adding a random number following a uniform distribution to the original values. The differential coefficient matrix uses non-uniform distribution perturbation, with the perturbation range decreasing as the number of iterations increases.

[0159] The mutated chromosome population is evaluated based on a multi-objective fitness function. The fitness function measures the performance of a control system, such as stability, response speed, and anti-interference capability. Chromosomes with better evaluation results are selected for an elite population pool.

[0160] Simulated annealing is used to perform local search optimization on chromosomes in the elite population pool. The initial temperature of the simulated annealing algorithm is determined based on the range of the multi-objective fitness function. The cooling coefficient is calculated adaptively, for example, adjusted according to the current temperature and the number of iterations. The termination temperature is dynamically adjusted based on the number of iterations, for example, decreasing as the number of iterations increases. During the local search process, discrete neighborhood structures are used for actuator reconstruction variables, such as bit-flip neighborhoods, swap neighborhoods, and insertion neighborhoods. Continuous neighborhood structures are used for controller parameter variables, such as Gaussian neighborhoods, uniform neighborhoods, and non-uniform neighborhoods. The acceptance probability of a new solution is calculated based on the multi-objective fitness function, and the acceptance probability decreases as the temperature decreases during iteration. The locally search-optimized chromosomes are then updated in the elite population pool.

[0161] The crossover, mutation, evaluation, and local search operations are repeated until a termination condition is met, such as reaching the maximum number of iterations. Finally, the chromosome with the highest fitness value is selected from the elite population pool as the optimal solution, i.e., the optimal configuration of the control system.

[0162] Beneficial effects:

[0163] 1. Improve control system performance: This method can simultaneously optimize the parameters and actuator configuration of the control system, thereby effectively improving the overall performance of the control system, such as stability, response speed and anti-interference ability.

[0164] 2. Enhance the search capability of the algorithm: By combining the global search capability of the genetic algorithm and the local search capability of the simulated annealing algorithm, it is possible to effectively avoid getting trapped in local optima, thereby finding a better control system configuration scheme.

[0165] 3. Adaptable to complex control systems: This method can handle both discrete and continuous variables and is applicable to various types of control systems, including complex nonlinear systems.

[0166] Figure 2 This is a schematic diagram of the fault diagnosis and fault-tolerant reconfiguration system for the satellite attitude control system according to an embodiment of the present invention, as shown below. Figure 2 As shown, the system includes:

[0167] The first unit is used to collect attitude angle data, angular velocity data, reaction wheel speed data, and torque data of the satellite attitude control system, and establish a Kalman filter state observer; embed the Kalman filter state observer into a fault diagnosis model, use the fault diagnosis model to generate a residual sequence, and perform adaptive dynamic threshold calculation on the residual sequence to obtain a fault feature vector; input the fault feature vector into a pre-trained deep belief network in the fault diagnosis model to obtain a fault diagnosis result, which includes the fault type and fault degree;

[0168] The second unit is used to determine the faulty component based on the fault diagnosis result, establish a mixed integer programming model based on the fault type and fault degree of the faulty component, set attitude control accuracy constraints and actuator power constraints in the mixed integer programming model, and solve the mixed integer programming model using an improved genetic algorithm to obtain a fault-tolerant control strategy, which includes an actuator reconfiguration scheme and controller parameters.

[0169] The third unit is used to reconfigure the actuator operating mode in the satellite attitude control system according to the actuator reconfiguration scheme in the fault-tolerant control strategy; update the controller of the satellite attitude control system using the controller parameters in the fault-tolerant control strategy; control the satellite attitude based on the updated controller to determine the satellite attitude error; and when the satellite attitude error exceeds a preset threshold, re-diagnose the fault using the fault diagnosis model.

[0170] A third aspect of the embodiments of the present invention,

[0171] An electronic device is provided, comprising:

[0172] processor;

[0173] Memory used to store processor-executable instructions;

[0174] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0175] Fourth aspect of the embodiments of the present invention,

[0176] A computer-readable storage medium is provided, having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0177] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0178] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the above embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for fault diagnosis and fault-tolerant reconstruction of a satellite attitude control system, characterized in that: include: Collect attitude angle data, angular velocity data, reaction wheel speed data, and torque data of the satellite attitude control system to establish a Kalman filter state observer; The Kalman filter state observer is embedded in a fault diagnosis model, a residual sequence is generated using the fault diagnosis model, and an adaptive dynamic threshold value is calculated on the residual sequence to obtain a fault feature vector; Inputting the fault feature vector into the pre-trained deep belief network in the fault diagnosis model to obtain a fault diagnosis result, wherein the fault diagnosis result includes a fault type and a fault degree; Determine the fault component according to the fault diagnosis result, establish a mixed integer programming model based on the fault type and the fault degree of the fault component; set attitude control accuracy constraint conditions and actuator power constraint conditions in the mixed integer programming model; use an improved genetic algorithm to solve the mixed integer programming model to obtain a fault-tolerant control strategy, wherein the fault-tolerant control strategy includes an actuator reconstruction scheme and controller parameters; The actuator working mode in the satellite attitude control system is reconfigured according to the actuator reconstruction scheme in the fault-tolerant control strategy; the controller of the satellite attitude control system is updated using the controller parameters in the fault-tolerant control strategy; the satellite attitude is controlled based on the updated controller to determine the satellite attitude error, and when the satellite attitude error exceeds a preset threshold, the fault diagnosis model is used to re-perform fault diagnosis.

2. The method according to claim 1, characterized in that Establish a Kalman filter state observer, including: The attitude angle data is used as the attitude quaternion observation quantity, and the angular velocity data is used as the angular velocity observation quantity, to construct a state vector including the attitude quaternion observation quantity and the angular velocity observation quantity; a satellite attitude dynamics model is established, and the reaction wheel speed data and the torque data are substituted into the satellite attitude dynamics model according to the Euler motion equation to obtain a state space expression, and the state space expression is discretized within a sampling time interval to obtain a discrete state equation; Construct an extended Kalman filter state observer, calculate the Jacobian matrix based on the discrete state equation, linearize the Jacobian matrix at the current state estimate to obtain a state transfer matrix, and calculate the observation matrix based on the observation equation; initialize the system noise covariance matrix and the observation noise covariance matrix of the extended Kalman filter state observer, and calculate the state one-step prediction value and the error covariance prediction value based on the initialized system noise covariance matrix and the observation noise covariance matrix; calculate the Kalman gain matrix using the state one-step prediction value and the error covariance prediction value, multiply the Kalman gain matrix by the observation residual to obtain a state estimation correction, and add the state estimation correction value to the state one-step prediction value to obtain a state estimation value; A new error covariance matrix is ​​calculated based on the state estimation value, and the state estimation accuracy is calculated according to the error covariance matrix. When the state estimation accuracy meets the preset accuracy requirement, the state estimation value is output; an observation residual sequence is calculated, and the actual covariance of the observation residual sequence is compared with the theoretical covariance to obtain an adaptive factor, and the system noise covariance matrix and the observation noise covariance matrix are adjusted using the adaptive factor, and the adjusted system noise covariance matrix and the observation noise covariance matrix are substituted into the extended Kalman filter state observer to obtain a Kalman filter state observer.

3. The method according to claim 2, characterized in that Establish satellite attitude dynamics model, including: Quaternions are used to describe the satellite attitude, a quaternion vector including a scalar part and a vector part is constructed, and the quaternion vector is associated with an angular velocity vector to establish an attitude kinematic equation, wherein the derivative of the quaternion vector is represented by the product of the quaternion correlation matrix and the angular velocity vector; an attitude dynamics equation is established based on the Euler rigid body dynamics principle, and the satellite moment of inertia matrix, the cross product term of the angular velocity vector, the control torque and the interference torque are substituted into the attitude dynamics equation, wherein the control torque is generated by a reaction wheel system; Establish a reaction wheel dynamics equation, substitute the reaction wheel angular momentum vector, the reaction wheel driving torque, and the cross product of the angular velocity vector and the reaction wheel angular momentum vector into the reaction wheel dynamics equation; combine the attitude kinematics equation, the attitude dynamics equation, and the reaction wheel dynamics equation to form a complete satellite attitude dynamics model, wherein the satellite attitude dynamics model includes a quaternion derivative equation, an angular velocity derivative equation, and a reaction wheel angular momentum derivative equation; discretize the satellite attitude dynamics model using a fourth-order Runge-Kutta numerical integration method to obtain a discrete state equation including a quaternion state update function, an angular velocity state update function, and a reaction wheel angular momentum state update function; A process noise term is added to the discrete state equation, a quaternion vector, an angular velocity vector and a reaction wheel angular momentum vector are used to form a state vector, a reaction wheel driving torque is used as a control input, and a state space expression containing process noise is established; an observation equation is established based on the attitude quaternion measurement value of the star sensor, the angular velocity measurement value of the gyroscope and the angular momentum measurement value of the reaction wheel speed sensor, measurement noise is added to the observation equation, and the state space expression and the observation equation are combined to form a satellite attitude dynamics model.

4. The method according to claim 1, characterized in that: Inputting the fault feature vector into the pre-trained deep belief network in the fault diagnosis model to obtain a fault diagnosis result includes: The fault feature vector is mapped to the interval between zero and one by the minimum and maximum value normalization method to obtain a normalized feature vector, a covariance matrix is ​​calculated according to the normalized feature vector, the covariance matrix is ​​subjected to eigenvalue decomposition to obtain principal components, and principal components whose cumulative contribution rate is greater than a preset threshold are selected to construct a dimensionality reduction mapping matrix; the normalized feature vector is subjected to dimensionality reduction processing by the dimensionality reduction mapping matrix to obtain a reduced dimensionality feature vector, and the reduced dimensionality feature vector is input into a visual layer of a deep belief network, wherein the deep belief network comprises a plurality of stacked restricted Boltzmann machines, and the number of neurons of each layer of restricted Boltzmann machines decreases; The reduced-dimensional feature vector is subjected to layer-by-layer feature conversion by a restricted Boltzmann machine of the deep belief network, the product of the restricted Boltzmann machine input vector and the weight matrix, and the sum of the bias vectors are calculated, the sum of the bias vectors is activated by a normalization function to obtain a hidden layer output, whether to set the neurons in the hidden layer output to zero is determined according to a preset random inactivation probability, and the hidden layer output that is not set to zero is numerically compensated according to a proportional coefficient, and the hidden layer output after random inactivation is used as the visual layer input of the next layer of restricted Boltzmann machine, and the feature transformation, normalization activation, random inactivation and numerical compensation processes are repeatedly performed until the hidden layer output of the last layer of restricted Boltzmann machine is obtained, and the random inactivation probability and the proportional coefficient are determined according to the number of neurons in each layer; The hidden layer output of the last layer of restricted Boltzmann machine is connected to the normalized classifier, the original score of each fault type is calculated according to the product of the hidden layer output and the classifier weight matrix and the bias vector, the original score is subjected to exponential operation to obtain a non-negative value, the non-negative value is divided by the exponential sum of all category scores to obtain the conditional probability distribution of different fault types, and the category with the largest conditional probability and greater than a preset probability threshold is selected as the fault type; the activation value of the normalized classifier output layer is weightedly combined to obtain the fault degree; the fault type and the fault degree are combined to generate a fault diagnosis result.

5. The method according to claim 1, characterized in that Determining a faulty component according to the fault diagnosis result, and establishing a mixed integer programming model based on the fault type and the fault degree of the faulty component, including: According to the fault diagnosis result, the fault type, fault degree and fault location information are obtained, and a fault impact assessment matrix is ​​established based on the fault type and the fault degree, wherein the fault impact assessment matrix includes the influence coefficient of the fault component on the attitude control accuracy and the influence coefficient of the fault component on the actuator power; the fault impact assessment matrix is ​​input into a mixed integer programming model, and attitude control variables are established in the mixed integer programming model, wherein the attitude control variables include a pitch angle, a yaw angle and a roll angle, and the product of the attitude control variable and the fault impact assessment matrix is ​​used as a predicted value of the attitude control accuracy under the fault state; An actuator power demand model is established based on the dynamic characteristics of the faulty component, the product of the actuator power demand model and the fault impact assessment matrix is ​​used as the actuator power prediction value under the fault state, and an actuator power constraint condition is set in the mixed integer programming model; a posture control accuracy threshold is determined according to the equipment posture control accuracy requirement, a posture control accuracy constraint condition is established based on the difference between the posture control accuracy prediction value and the posture control accuracy threshold, and the posture control accuracy constraint condition is set in the mixed integer programming model; An objective function is established according to the fault impact assessment matrix, wherein the objective function includes an attitude control deviation minimization item, an actuator power consumption minimization item and a fault impact compensation item, and the objective function, the attitude control variable, the actuator power constraint condition and the attitude control accuracy constraint condition are combined to construct a mixed integer programming model, wherein the attitude control deviation minimization item is calculated based on the deviation between the attitude control variable and the desired attitude, the actuator power consumption minimization item is calculated based on the actuator power prediction value, and the fault impact compensation item is calculated based on the influence coefficient in the fault impact assessment matrix.

6. The method according to claim 1, characterized in that The mixed integer programming model is solved by using an improved genetic algorithm to obtain a fault-tolerant control strategy, which includes an actuator reconstruction scheme and controller parameters, including: The actuator reconstruction variables and controller parameter variables in the mixed integer programming model are encoded as chromosomes, wherein the actuator reconstruction variables include actuator working state flags, actuator output allocation coefficients, and actuator switching timings, and the controller parameter variables include a proportional coefficient matrix, an integral coefficient matrix, and a differential coefficient matrix; a multi-objective fitness function is established based on the actuator reconstruction variables and the controller parameter variables, and the multi-objective fitness function includes an actuator resource utilization evaluation item, a control performance evaluation item, a fault tolerance capability evaluation item, and an energy consumption evaluation item; Initializing the chromosome population according to the multi-objective fitness function; stratifying the chromosome population to construct a population stratification structure, establishing an elite population library in the population stratification structure, calculating the fitness value of the chromosome according to the multi-objective fitness function and storing the non-dominated solution chromosome in the elite population library; performing a crossover operation on the chromosome population using an adaptive crossover operator, wherein the crossover probability of the adaptive crossover operator is determined based on the chromosome distribution in the elite population library and the multi-objective fitness function value, generating a chromosome after crossover and adding it to the chromosome population; A hybrid mutation operator is used to perform mutation operation on the chromosome population after crossover to obtain an elite population library; a clustering analysis is performed on the chromosomes in the elite population library, and the chromosomes are grouped using a clustering algorithm, and the chromosomes with the best multi-objective fitness function value are selected from each group to form an optimal solution set; the chromosomes in the optimal solution set are decoded to obtain a fault-tolerant control strategy, which includes an actuator working mode switching timing table, an actuator output allocation matrix, an actuator collaborative control scheme, a proportional coefficient matrix, an integral coefficient matrix, and a differential coefficient matrix.

7. The method according to claim 6, characterized in that A hybrid mutation operator is used to perform mutation operation on the chromosome population after crossover to obtain an elite population library, including: A hybrid mutation operator is used to perform mutation operation on the chromosome population after crossover, and the hybrid mutation operator includes a discrete mutation operation and a continuous mutation operation; wherein the discrete mutation operation changes the value of the actuator working state flag bit by bit reversal, adjusts the distribution of the actuator output allocation coefficient by random perturbation, and modifies the time point of the actuator switching timing by binary code conversion; the continuous mutation operation adjusts the parameter value of the proportional coefficient matrix by Gaussian distribution perturbation, adjusts the parameter value of the integral coefficient matrix by uniform distribution perturbation, and adjusts the parameter value of the differential coefficient matrix by non-uniform distribution perturbation to obtain the mutated chromosome population; Based on the multi-objective fitness function, the mutated chromosome population is evaluated, and the chromosomes with better evaluation results are selected into the elite population library; the chromosomes in the elite population library are locally searched and optimized using a simulated annealing algorithm, wherein the initial temperature of the simulated annealing algorithm is determined according to the value range of the multi-objective fitness function, the cooling coefficient of the simulated annealing algorithm is calculated in an adaptive manner, and the termination temperature of the simulated annealing algorithm is dynamically adjusted based on the number of iterations; In the local search process of the simulated annealing algorithm, a discrete domain structure is used for actuator reconstruction variables, and the discrete domain structure includes a bit flip neighborhood, a swap neighborhood, and an insertion neighborhood; a continuous domain structure is used for controller parameter variables, and the continuous domain structure includes a Gaussian neighborhood, a uniform neighborhood, and a non-uniform neighborhood; the acceptance probability of the new solution is calculated based on the multi-objective fitness function, and the acceptance probability decreases as the temperature decreases during the iteration process, and the chromosome optimized by the local search is updated to the elite population library.

8. A satellite attitude control system fault diagnosis and fault-tolerant reconstruction system, used to implement the method described in any one of claims 1 to 7, characterized in that: include: The first unit is used to collect attitude angle data, angular velocity data, reaction wheel speed data, and torque data of the satellite attitude control system to establish a Kalman filter state observer; The Kalman filter state observer is embedded in a fault diagnosis model, a residual sequence is generated using the fault diagnosis model, and an adaptive dynamic threshold value is calculated on the residual sequence to obtain a fault feature vector; Inputting the fault feature vector into the pre-trained deep belief network in the fault diagnosis model to obtain a fault diagnosis result, wherein the fault diagnosis result includes a fault type and a fault degree; The second unit is used to determine the faulty component according to the fault diagnosis result, establish a mixed integer programming model based on the fault type and the fault degree of the faulty component; set attitude control accuracy constraint conditions and actuator power constraint conditions in the mixed integer programming model; use an improved genetic algorithm to solve the mixed integer programming model to obtain a fault-tolerant control strategy, wherein the fault-tolerant control strategy includes an actuator reconstruction scheme and controller parameters; The third unit is used to reconfigure the actuator working mode in the satellite attitude control system according to the actuator reconstruction scheme in the fault-tolerant control strategy; update the controller of the satellite attitude control system using the controller parameters in the fault-tolerant control strategy; control the satellite attitude based on the updated controller to determine the satellite attitude error, and when the satellite attitude error exceeds a preset threshold, re-diagnose the fault using the fault diagnosis model.

9. An electronic device, characterized in that: include: processor; a memory for storing processor-executable instructions; The processor is configured to call the instructions stored in the memory to execute the method described in any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that: When the computer program instructions are executed by a processor, the method according to any one of claims 1 to 7 is implemented.

Images