Agile satellite attitude control method suitable for river imaging

Abstract

The invention discloses an agile satellite attitude control method suitable for river imaging, and the method comprises the steps: firstly planning the maneuvering time of each attitude and the imaging initial time and terminal time of each imaging strip according to an imaging task, generating a fold line strip, and calculating the attitude, angular velocity and angular acceleration information in the maneuvering process; then defining a system state variable, and further establishing an attitude system discretization model; solving an LQR problem to obtain a feedback gain matrix, further calculating a minimum robust invariant set according to the feedback gain matrix, finally solving an optimization problem to generate an optimal nominal control sequence, and combining with a feed-forward torque compensation amount calculation module to output instruction torque; and finally, calculating a frame angle instruction angular velocity of the control moment gyroscope group according to the instruction moment, and realizing moment output through a singular robust inverse instruction operation law. According to the method, the winding track imaging process of the river basin is improved to a great extent, the whole satellite can be rapidly maneuvered and rapidly stabilized, and a better imaging effect can be obtained.

Description

Technical Field

[0001] The present invention relates to the technical field of remote sensing satellite attitude control, and specifically to an agile satellite attitude control method suitable for river imaging. Background Art

[0002] A single remote sensing satellite can cover a river basin of thousands of kilometers, realizing synchronous acquisition of hydrological parameters in the whole basin and large-scale real-time monitoring, which is of great significance for environmental protection, disaster warning and water resource management. Since the swath width of high-resolution satellites is relatively small, generally between 10 km and 14 km, in order to image a large area, agile satellites are used to quickly adjust the attitude using control moment gyros. However, due to the fact that large actuators have greater torque noise when outputting large torques, and in addition, during large angular velocity maneuvers, the measurement errors of satellite attitude measurement sensors such as star sensors and gyros will become worse. During the imaging of meandering rivers, continuous maneuvering imaging is required, and it is difficult to achieve good attitude control effects, which affects remote sensing imaging. Using the traditional proportional-integral-derivative control (hereinafter referred to as PID) for attitude maneuvering to the stable imaging transition time is long. Using the method of standard model predictive control is based on the working conditions of accurate models, and the attitude control effect will also be affected in the case of large external disturbances, and the robustness is very challenging. Summary of the Invention

[0003] The purpose of the present invention is to overcome the deficiencies of the prior art and propose an agile satellite attitude control method suitable for river imaging, so as to solve the technical problems of fast and stable imaging of meandering river channels and maintaining robustness. This method greatly improves the imaging process of the meandering trajectory in the river basin, and the entire satellite can quickly maneuver and quickly stabilize, and can obtain better imaging effects.

[0004] In order to achieve the above purpose, the specific technical solution adopted by the present invention is as follows:

[0005] An agile satellite attitude control method suitable for river imaging, comprising the following steps:

[0006] Step 1, imaging trajectory planning, planning the attitude maneuver time for each time and the initial and terminal times of each imaging strip according to the imaging task, generating a polyline strip, and calculating the attitude, angular velocity and angular acceleration information during the maneuver;

[0007] Step 2, designing a tube model predictive controller. First, establish the spacecraft attitude dynamics and kinematics models; then define the system state variables, and then establish the discrete model of the attitude system; solve the LQR problem to obtain the feedback gain matrix, and then calculate the minimum robust invariant set according to the feedback gain matrix. Finally, solve the optimization problem to generate the optimal nominal control sequence, and calculate the command torque by combining with the feedforward torque compensation amount calculation module;

[0008] Step 3: Calculate the frame angle command angular velocity of the control moment gyro group according to the command torque, and realize the torque output through the singular robust inverse command operation law.

[0009] Preferably, the imaging trajectory planning in the above Step 1 specifically includes:

[0010] Step 1.1: Obtain the mission requirement information, satellite orbit information and ground station parameter information, determine the maximum available time window of the observation area, generate multiple broken line strips, and record the starting and ending longitudes and latitudes, altitudes and lengths of each strip;

[0011] Step 1.2: Based on the satellite inertia and the maneuverability of the control moment gyro, calculate the maneuver time and the corresponding attitude, angular velocity, and angular acceleration information for each maneuver;

[0012] Step 1.3: Adopt the seventh-order polynomial trajectory planning method to judge whether the angular momentum exceeds the angular momentum envelope threshold of the control moment gyro during the maneuver. If it exceeds, terminate the maneuver; otherwise, perform subsequent attitude maneuvers.

[0013] Preferably, the specific method of the above Step 1.2 is as follows: According to the satellite inertia and the maneuverability of the control moment gyro, the relationship table between the maneuver time and the maneuver angle can be pre-calculated. Interpolate the maneuver time for each maneuver according to the required attitude angle; divide the length of the imaging strip by the moving speed of the sub-satellite point to obtain the imaging time of each strip, and then the initial and terminal moments of imaging for each strip can be confirmed in turn, so as to calculate the maneuver time for each maneuver and the attitude, angular velocity, and angular acceleration information at the initial and terminal moments of imaging for each imaging strip; at the starting moment and the terminal moment of the strip, take two sets of orbital elements near these two points, and obtain the numerical values of the attitude, angular velocity, and angular acceleration at the initial and terminal moments through iterative calculation.

[0014] Preferably, the above Step 2 further includes discretizing the attitude dynamics rewritten model as follows:

[0015] ;

[0016] ;

[0017] ;

[0018] In the formula, k represents the current moment, and k + 1 represents the next moment. and respectively refer to the system state variables at the next moment and the current moment. , represents the influence of environmental disturbances on the system dynamics at each time step k, that is, the disturbance. is the upper limit of the perturbation, is the lower limit of the perturbation.

[0019] Preferably, in the rewritten discretized model, both the system state variables and the system control variables are constrained. Among them, the constraint of the system control variable is subject to the angular momentum sphere envelope and the maximum torque, and its constraint formula is:

[0020] ;

[0021] ;

[0022] In the formula and are constraint matrices, is formed by stacking the identity matrix to form a 12 6 matrix, is formed by stacking the identity matrix to form matrix and are respectively the upper and lower limits of the state variables. and are respectively the upper and lower limits of the control variables.

[0023] Preferably, in step 2, the method for solving the feedback gain matrix by solving the LQR problem is as follows: within the prediction horizon N, the system state can be decomposed into a nominal part and an error , the error represents the actual state relative to the nominal state after the time step , and we get , where is the nominal control input, and the feedback gain matrix is designed offline to make the system quadratically stable. The solution of the feedback gain matrix is a classical LQR problem that satisfies , where T is the device. Given , and and , we can directly solve , corresponds to 6 state variables, namely attitude and angular velocity; is the control output weight matrix , which reflects the torque limitation.

[0024] Preferably, in step 2, the method for solving the minimum robust invariant set is:

[0025] According to the calculated feedback gain matrix , the minimum robust invariant set is calculated by the Minkowski sum calculation method:

[0026] ; ;

[0027] is the Minkowski sum calculation,

[0028] stop when, and finally .

[0029] Preferably, when calculating the minimum robust invariant set, a constrained tightening calculation is introduced:

[0030] ;

[0031] wherein, is the Pontryagin set difference calculation.

[0032] Preferably, the objective function of the tube model predictive controller optimization problem in step 2 is:

[0033] ;

[0034] wherein, and are the nominal state and input, is the nominal state at the end of the prediction horizon; the obtained by solving the optimization problem is the optimal nominal control sequence the control action at this moment.

[0035] Preferably, the optimal nominal control sequence satisfies: ; .

[0036] Preferably, the feedforward torque compensation amount is calculated by the following formula:

[0037] ;

[0038] wherein, is the expected angular velocity of the previous beat, is the control period, is the adjustment coefficient.

[0039] Preferably, the control moment gyro group is installed in a pyramid configuration and includes five control moment gyros, and its combined torque is calculated by the following formula:

[0040] ;

[0041] Among them, Each column of the matrix is the angular momentum direction of the high-speed rotation of the control moment gyro when the gimbal angle of the control moment gyro is at 90°; Each column of the matrix is the angular momentum direction of the high-speed rotation of the control moment gyro when the gimbal angle of the control moment gyro is at the zero position; is the Jacobi matrix of the torque output of the control moment gyro group; is the gimbal angle vector.

[0042] Preferably, the singular robust inverse command operation law is:

[0043] ;

[0044] Among them, is the adjustment coefficient.

[0045] Preferably, the method updates the nominal state and control input in real time through a rolling horizon control process, specifically including: measuring the current system state, initializing the nominal state, solving the optimization problem, executing the control, and updating the system state.

[0046] The present invention has the following characteristics and beneficial effects:

[0047] The present invention proposes a tube model predictive control agile satellite attitude control method, in which the tube model predictive control designs the feedback gain matrix offline and only needs to optimize the nominal trajectory online. The influence of disturbances is compensated in real time through feedback, and the computational complexity is significantly reduced; when dealing with bounded disturbances, by designing a robust positive invariant set, it is ensured that the closed-loop system can still be asymptotically stable under disturbances. Thus, the technical problem of fast, stable, and robust imaging of meandering rivers is effectively solved. This method greatly improves the imaging process of meandering trajectories in river basins. The entire satellite can quickly maneuver and quickly stabilize, and better imaging effects can be obtained. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 is the block diagram of the tube satellite attitude controller of the present invention;

[0049] Figure 2 is the strip generation diagram of the meandering waters of the Yangtze River in this embodiment;

[0050] Figure 3 is the simulation result of the attitude control process of the meandering waters of the Yangtze River of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0051] The present invention will be described in detail below with reference to specific embodiments. The following embodiments will help those skilled in the art to further understand the present invention, but do not limit the present invention in any form. It should be noted that, without conflict, the embodiments in the present invention and the features in the embodiments can be combined with each other.

[0052] An agile satellite attitude control method applicable to river imaging, as Figure 1 shown, includes the following steps:

[0053] Step 1, imaging trajectory planning. According to the imaging task, plan the attitude maneuver time for each time, the initial time and the terminal time of each imaging strip, generate a broken-line strip, and calculate the attitude, angular velocity and angular acceleration information during the maneuvering process.

[0054] In this embodiment, for the area between Figure 2 latitude 26 degrees and 28 degrees where the Yangtze River Basin meanders as shown, generate a broken-line strip.

[0055] Specifically, it includes the following sub-steps:

[0056] Step 1.1. First, obtain the task requirement information, satellite orbit information and ground station parameter information to find the maximum available time window of the observation area. This example is based on the confirmed existence of the observation time and the observation imaging angle. According to the main standard points of the given Yangtze River imaging area, calculate the azimuth angle between the first two points as the reference azimuth angle base_azimuth.

[0057]

[0058] In the formula , are the latitudes of the starting point and the ending point, is the azimuth angle of the starting point and the ending point, is the longitude difference between the starting point and the ending point.

[0059] Traverse the subsequent trajectory segments and perform the following operations: calculate the azimuth angle of the two points in the current segment and calculate the absolute deviation from the reference azimuth angle: , as is greater than the threshold. In this embodiment, the threshold is set to 5 degrees, then stop merging. If it is less than the threshold, merge until the last point is processed. In this way, multiple broken-line strips can be obtained, and record the longitude and latitude, altitude, and strip length information of the starting and ending points of each strip.

[0060] Step 1.2. Based on the satellite inertia and the maneuvering ability of the control moment gyro, calculate the maneuvering time for each time and the corresponding attitude, angular velocity, and angular acceleration information.

[0061] Specifically, according to the inertia of the satellite and the maneuvering ability of the control moment gyroscopes, a relationship table of maneuvering time and maneuvering angle (as shown in Table 1) can be pre-simulated and calculated. The duration of each maneuver is interpolated based on the required maneuvering attitude angle. The imaging time of each strip can be obtained by dividing the length of the imaging strip by the moving speed of the sub-satellite point. In this way, the attitude, angular velocity, and angular acceleration information at the initial and terminal times of each maneuver and each imaging strip can be calculated. At the starting time of the strip and the terminal time Two sets of orbital elements are taken near these two points for interpolation in the entire calculation process.

[0062]

[0063] Table 1 Table of attitude maneuvering time correspondence

[0064] It can be understood that the moving speed of the sub-satellite point is obtained according to the orbital parameter values in this embodiment. In this embodiment, the moving speed of the sub-satellite point is 7.2 km / s.

[0065] The numerical values of the attitude, angular velocity, and angular acceleration at the initial time are obtained through iterative calculation. Similarly, the target point is changed to the sub-satellite target point at the terminal time. Define as the temporary iteration time, rounded to the nearest second; is the Greenwich sidereal time angle at time

[0066] In this embodiment, the iteration condition is The attitude, angular velocity, and angular acceleration values at the terminal time of the strip can also be obtained by subtracting the start time of the iterative calculation from the terminal time. The calculation method is as follows:

[0067] (Attitude maneuvering time - start time of iterative calculation));

[0068] ;

[0069] ;

[0070] ;

[0071] and are the angle and angular velocity values, and are the temporary storage values of the angle and angular velocity, saving the data of the previous step for differential calculation; ; , is the imaging duration ; In the formula is The distance between the geocenter and point P on the strip at the moment, is The distance between the geocenter and the starting point of the strip at the moment, is The distance between the geocenter and the ending point of the strip at the moment.

[0072] ; is The geocentric subtended angle at the moment, is the unit vector pointing from the geocenter to the starting point of the strip, perpendicular to , , is the unit vector pointing from the geocenter to the ending point of the strip:

[0073] The unit vector pointing from the geocenter to point P on the strip ;

[0074] The vector pointing from the geocenter to point P on the strip ;

[0075] According to the orbit extrapolation result, the orbital elements are obtained according to and interpolated to obtain which are respectively The semi-major axis, eccentricity, orbital inclination, argument of latitude, right ascension of the ascending node, true anomaly, geocentric distance, and the components of the unit vector pointing from the geocenter to the satellite in the Earth-fixed coordinate system at the moment.

[0076] The orbital angular velocity ; is the Earth's gravitational constant 398600.5.

[0077] The coordinate transformation matrix from the Earth-fixed system to the orbital system ; is the coordinate transformation matrix from the Earth-fixed system to the inertial system is the coordinate transformation matrix from the inertial system to the orbital system;

[0078] The components of the vector pointing from the geocenter to the satellite in the Earth-fixed coordinate system at the moment ;

[0079] The components of the satellite velocity in the orbital system at the moment ;

[0080] The components of the velocity of the pointing point in the orbital coordinate system ;

[0081] The components of the vector pointing from the satellite to point P on the strip at the moment in the orbital coordinate system ;

[0082] Modulo operation ;

[0083] Roll angle ;

[0084] Pitch angle ;

[0085] Coordinate transformation matrix from orbit system to reference attitude system is generated by matrix through 123 transformation sequence;

[0086] Differential calculation of roll angular velocity and pitch angular velocity ;

[0087] Satellite reference angular velocity ; where is the orbital angular velocity

[0088] ;

[0089] ;

[0090]

[0091] Components of the linear velocity of the ground pointing point caused by the satellite reference angular velocity in the reference attitude coordinate system ;

[0092] Components of the linear velocity of the pointing point relative to the satellite in the reference attitude coordinate system .

[0093] ;

[0094]

[0095]

[0096] Angular acceleration vector .

[0097] Step 1.3: Perform trajectory planning based on the attitude values at the imaging start time and end time, and then determine whether the maneuverability will exceed the limit during the trajectory. In this embodiment, the seventh-order polynomial trajectory planning method is used to determine whether the angular momentum exceeds the angular momentum envelope threshold of the control moment gyro during the maneuver. If it exceeds, the maneuver is terminated; otherwise, the subsequent attitude maneuver is carried out.

[0098] First, define the trajectory planning coefficients , and the expression of the trajectory planning coefficients is as follows:

[0099]

[0100] ;

[0101] Among them, is the duration of attitude maneuver; are polynomial coefficients, where , can be set as the three-axis attitude angles respectively. The calculation formulas are the same. The subscript 0 represents the initial attitude angle of attitude maneuver, and the subscript represents the terminal attitude angle of attitude maneuver; is the three-axis attitude angle at the initial moment of attitude maneuver, is the three-axis attitude angle at the end moment of attitude maneuver.

[0102] According to the polynomial coefficients, calculate the target attitude at time , angular velocity , and angular acceleration :

[0103] ;

[0104] where t is relative to the initial moment of attitude maneuver, 0 ≤ t ≤ tm, can be set as the three-axis attitude angles roll angle , pitch angle , and yaw angle respectively, that is , , .

[0105] Iteratively calculate the three-axis attitude ( is the pitch angle, is the roll angle, is the yaw angle) and Euler angular velocity from the set maneuver process time:

[0106] (3)

[0107] Then calculate the angular momentum of the stereo imaging satellite:

[0108]

[0109] Among them is the value of the principal axis of the overall satellite inertia, is the value of the angular momentum of the satellite body during the maneuver process. Calculate the angular momentum during the maneuver process and judge whether it exceeds the threshold set by the control moment gyro angular momentum envelope. If it exceeds, it reflects that the angular velocity maneuver ability cannot meet the requirements of the maneuver time, and stop the subsequent attitude maneuver actions. If it does not exceed, the subsequent attitude maneuver can be carried out.

[0110] Step 2: Design a tube model predictive controller. First, establish the spacecraft attitude dynamics and kinematics models; then define the system state variables, and further establish the discretized model of the attitude system; solve the LQR problem to obtain the feedback gain matrix, and then calculate the minimum robust invariant set according to the feedback gain matrix. Finally, solve the optimization problem to generate the optimal nominal control sequence, and calculate the command torque by combining the feedforward torque compensation amount calculation module.

[0111] It can be understood that as Figure 1 shown, the satellite attitude control includes a maneuvering process target attitude determination module, a feedforward torque compensation amount calculation module, and a tube model predictive control module. The maneuvering process target attitude determination module gives the real-time target attitude angle, angular velocity, and angular acceleration according to the imaging requirements; the feedforward torque compensation amount calculation module calculates the feedforward torque; the tube model predictive control module calculates the optimal control sequence according to the tube model predictive control principle, and uses the control quantity and the feedforward torque to superimpose and output the control torque to the satellite body through the control moment gyro group for attitude control; the attitude measurement and determination module calculates the satellite attitude and angular velocity in real time through the on-board gyro and star sensor.

[0112] The determination of the target attitude reference trajectory in the maneuvering process uses the polynomial trajectory planning calculation in Step 1.3. According to the three-axis attitude angle , angular velocity , and angular acceleration at the initial imaging moment, the target attitude angle , angular velocity , and angular acceleration at the end of imaging, and the time deviation relative to the starting moment, calculate the expected three-axis attitude angle , angular velocity , and angular acceleration at each moment through the polynomial formula. After the target attitude is generated, first convert the attitude through the 123 rotation sequence to the attitude transformation matrix from the target pointing to the orbital system, calculate the attitude transformation matrix from the orbital system to the inertial system through the orbital data, convert the attitude transformation matrix and the attitude transformation matrix from Euler angles to the reference trajectory quaternion, and the reference trajectory angular velocity is:

[0113]

[0114] Specifically, it includes the following sub-steps:

[0115] Step 2.1: Establish the spacecraft attitude dynamics and kinematics models as follows:

[0116] ;

[0117]

[0118] In the formula, is the inertial tensor of the entire satellite; is the angular velocity of the spacecraft about three axes; represents the control torque about three axes of the satellite, is the angular momentum of the actuator; A torque is generated by changing the angular momentum of the actuator. Here, the actuator is a group of control moment gyros; is the disturbance torque. In engineering, the Euler angles are often used to describe the motion model of the satellite attitude, , which are the roll angle, pitch angle, and yaw angle respectively.

[0119] Step 2.2. Define the system state variables (6 1), and establish the discretized model of the attitude system as follows:

[0120] The tube model predictive control rewrites the attitude dynamics and kinematics models in 2.1 into the discretized model in the following form:

[0121]

[0122] In the formula , represents the influence on the system dynamics at each time step k due to environmental disturbances. It is defined as an independent and identically distributed random variable, which is polyhedral and convex, and its value is also bounded.

[0123]

[0124] Both the state variables and control variables of this system are subject to constraints. The constraints on the control variables are mainly restricted by the angular momentum sphere envelope and the maximum torque, and their constraint formulas are:

[0125] ;

[0126] ;

[0127] In the formula , is the constraint matrix, is formed by stacking the identity matrix to form a 12 6 matrix, is formed by stacking the identity matrix to form matrix and are the upper and lower limits of the state variables respectively. and They are the upper and lower limits of the control variables, respectively.

[0128] Step 2.3: Solve the LQR problem to obtain and :

[0129]

[0130] where , and are the weight matrices; is the reference trajectory.

[0131] Within the prediction horizon N, the system state can be decomposed into the nominal part and the error , and the error represents the deviation of the actual state relative to the nominal state at the time step after, resulting in , where is the nominal control input, and the feedback gain matrix is designed offline to make the system quadratically stable and obtained by solving the LQR.

[0132] The solution of the feedback gain matrix is a classical LQR problem satisfying , given , as well as and can be directly solved after , corresponds to 6 state variables, namely the attitude quaternion and three angular velocities; the control output weight matrix reflects the torque limit.

[0133] Step 2.4: Solve the minimum robust invariant set:

[0134] According to the calculated feedback gain matrix , calculate the minimum robust invariant set through the Minkowski sum calculation method:

[0135] ; ;

[0136] is the Minkowski sum calculation, represents the i-th power of the closed-loop state matrix. For example, for simplified calculation, it can be considered that = is an invariant;

[0137] Adopt an iterative calculation method: ;

[0138] Stop when .

[0139] Calculation of constraint tightening

[0140] ,

[0141] wherein, is the Pontryagin set difference calculation .

[0142] The nominal state Z and input V need to satisfy the tightened constraints to ensure that the actual state X and input U are always within the original constraints.

[0143] Step 2.5, Finally, solve the optimization problem. The objective function of the tube model predictive controller optimization problem is:

[0144] ;

[0145] wherein, and are the nominal state and input, is the nominal state at the end of the prediction horizon; obtained by solving the optimization problem is the optimal nominal control sequence the control action at this moment. This sequence satisfies: The optimal nominal control sequence satisfies:

[0146] ;

[0147] .

[0148] It should be noted that the algorithm flow for each cycle is as follows:

[0149] (1) Measure the state: Obtain the real system state;

[0150] (2) Initialize the nominal state: Set (usually aligned with );

[0151] (3) Solve the optimization problem: Minimize the objective function ;

[0152] (4) Execute the control: ;

[0153] (5) Update the state: ;

[0154] (6) Roll the time domain: Move to the next moment , repeat steps 1 - 5.

[0155] Step 2.6, calculate the commanded torque by combining the feed - forward torque compensation amount calculation module output:

[0156] The feed - forward torque compensation amount is calculated by the following formula:

[0157] ;

[0158] where, is the expected angular velocity value of the previous beat, is the control period, is the adjustment coefficient.

[0159] Based on the above, the feed - forward model predictive control law is:

[0160] ;

[0161] is the commanded torque of the satellite, used to output to the control moment gyro group; is the feed - forward control amount at the k - th moment; is the tube model predictive control amount at the k - th moment.

[0162] Step 3, calculate the gimbal angle commanded angular velocity of the control moment gyro group according to the commanded torque, and realize the torque output through the singular robust inverse command operation law.

[0163] It can be understood that, according to the satellite's body three - axis commanded torque obtained in step 2, calculate the gimbal angle commanded angular velocity of the control moment gyro group.

[0164] Specifically, in this embodiment, the control moment gyro group includes five control moment gyros (CMGs) and is installed in a pyramid configuration.

[0165] The control moment gyro group includes five control moment gyros (CMGs) and is installed in a pyramid configuration. The nominal angular momentum of each CMG is , and the synthetic angular momentum of the CMGs :

[0166] ;

[0167] where, is the CMGs gimbal angle vector array, is the gimbal angle of the i - th control moment gyro; and are matrices related to the CMGs installation orientation. Each column of is the angular momentum direction of the CMG rotating at high speed when the gimbal angle of each CMG is 90°, Each column represents the angular momentum direction of the high-speed rotation of the CMG when the frame angle of each CMG is at zero position. is an n×1 unit vector.

[0168] The resultant torque generated by the control moment gyro group :

[0169]

[0170] where Each column of the matrix represents the angular momentum direction of the high-speed rotation of the control moment gyro when the frame angle of the control moment gyro is at 90°; Each column of the matrix represents the angular momentum direction of the high-speed rotation of the control moment gyro when the frame angle of the control moment gyro is at zero position; is the Jacobi matrix of the torque output of the control moment gyro group;

[0171] The operation rate of the control moment gyro group is:

[0172]

[0173] where is the angular velocity command of the frame angle; is the singular robust inverse command operation law; is the zero motion singularity avoidance operation law; is the nominal frame angle regression operation law;

[0174] It should be noted that there are various implementation methods for the zero motion singularity avoidance operation law and the frame angle regression operation law, which are well-known content. In this embodiment, the singular robust inverse command operation law is mainly described as follows:

[0175] Singular robust inverse command operation law:

[0176] ;

[0177] is the adjustment coefficient; according to the extreme value condition, it can be known that: .

[0178] According to the implementation process provided in this embodiment, the tube model predictive control effect diagram in the embodiment of the present invention is obtained. In the imaging of the meandering Yangtze River Basin in this example, a complex maneuvering imaging process can be completed.

[0179] In summary, by using the satellite attitude control method based on broken line strip planning and tube model predictive control of the present invention, the entire satellite can quickly maneuver and quickly stabilize, giving full play to the satellite's effectiveness.

[0180] The above has shown and described the basic principles, main features and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited by the above embodiments. The above embodiments and the descriptions in the specification are only preferred examples of the present invention and are not used to limit the present invention. Without departing from the spirit and scope of the present invention, the present invention will have various changes and improvements, and these changes and improvements all fall within the scope of the present invention claimed. The scope of protection claimed by the present invention is defined by the appended claims and their equivalents.

Claims

1. An agile satellite attitude control method applicable to river imaging, characterized in that, It includes the following steps: Step 1, imaging trajectory planning: According to the imaging task, plan the attitude maneuver time for each time, the initial and terminal times of each imaging strip, generate a polyline strip, and calculate the attitude, angular velocity, and angular acceleration information during the maneuver, and then judge whether to perform subsequent attitude maneuvers; Step 2, design a tube model predictive controller: First, establish a spacecraft attitude dynamics and kinematics model; then define system state variables, and then establish a discretized model of the attitude system; solve the LQR problem to obtain the feedback gain matrix, and then calculate the minimum robust invariant set according to the feedback gain matrix. Finally, solve the optimization problem to generate the optimal nominal control sequence, and calculate the command torque by combining the feedforward torque compensation amount calculation module; Step 3, calculate the frame angle command angular velocity of the control moment gyro group according to the command torque, and realize the torque output through the singular robust inverse command operation law.

2. The method according to claim 1, wherein The imaging trajectory planning in Step 1 specifically includes: Step 1.1, obtain the task requirement information, satellite orbit information, and ground station parameter information, determine the maximum available time window of the observation area, generate multiple polyline strips, and record the starting and ending longitudes, latitudes, altitudes, and lengths of each strip; Step 1.2, based on the satellite inertia and the maneuverability of the control moment gyro, calculate the maneuver time for each time and the corresponding attitude, angular velocity, and angular acceleration information; Step 1.3, adopt the seventh-order polynomial trajectory planning method, calculate the angular momentum according to the maneuver time for each time and the corresponding attitude, angular velocity, and angular acceleration information, and then judge whether the angular momentum exceeds the angular momentum envelope threshold of the control moment gyro during the maneuver. If it exceeds, terminate the maneuver; otherwise, perform subsequent attitude maneuvers.

3. The method according to claim 2, wherein The specific method of Step 1.2 is: According to the inertia of the satellite and the maneuverability of the control moment gyro, pre-calculate the relationship table between the maneuver time and the maneuver angle, and interpolate the corresponding relationship table according to the attitude angle to be maneuvered to calculate the duration of each maneuver; The imaging time of each strip is obtained by dividing the length of the imaging strip by the moving speed of the sub-satellite point. The initial and terminal moments of imaging for each strip are confirmed in sequence, that is, the maneuvering time for each time, as well as the attitude, angular velocity, and angular acceleration information at the initial and terminal moments of imaging for each imaging strip are calculated; at the starting moment of the strip and the terminal moment Two sets of orbital elements are taken near the two points. The numerical values of the attitude, angular velocity, and angular acceleration at the initial and terminal moments are obtained through iterative calculation.

4. The method according to claim 1, wherein In Step 2, it also includes rewriting the attitude dynamics to obtain a discretized model: ； ； ； where k represents the current time step, and k + 1 represents the next time step, and refer to the system state variables at the next time step and the current time step respectively, represents the system control variable, , represents the impact of environmental disturbances on the system dynamics at each time step k, i.e., the perturbation; is the upper bound of the perturbation, .

5. The method according to claim 4, characterized in that In the discretized model, both the system state variables and the system control variables are constrained, The constraint formula for the system state variables is: ； The constraint formula for the system control variables is: ; where and are constraint matrices, is formed by stacking the identity matrix to form a 12 × 6 matrix, is formed by stacking the identity matrix to form the matrix , and are the upper and lower limits of the state variables, respectively. and are the upper and lower limits of the control variables, respectively.

6. The method according to claim 5, characterized in that, In the said step 2, the method for solving the LQR problem to obtain the feedback gain matrix is as follows: within the prediction time domain N, the system state can be decomposed into a nominal part and an error . The error represents the deviation of the actual state relative to the nominal state at the time step later, so . Among them, is the nominal control input, and the feedback gain matrix is designed offline to make the system quadratically stable. The solution of the feedback gain matrix satisfies , where T is the device. Given , and and , can be directly solved , corresponding to the state variables; is the control output weight matrix , reflecting the torque limitation.

7. The method according to claim 6, wherein In Step 2, the solution method for the minimum robust invariant set is: According to the calculated feedback gain matrix , the minimum robust invariant set is calculated by the Minkowski sum calculation method: ； ； For Minkowski sum calculation, representing the i-th power of the closed-loop state matrix; stop at that time, and finally .

8. The method according to claim 7, wherein When calculating the minimum robust invariant set, introduce the calculation of constraint tightening: ； Among them, is the Pontryagin set difference calculation.

9. The method according to claim 7, wherein The objective function of the optimization problem of the tube model predictive controller in Step 2 is: ； wherein, and are the nominal state and the input, is the nominal state at the end of the prediction horizon; obtained by solving the optimization problem is the optimal nominal control sequence which is the control action at this moment.

10. The method according to claim 9, wherein The optimal nominal control sequence satisfies: ; .

11. The method according to claim 9, characterized in that, The feedforward torque compensation amount is calculated by the following formula: ； Among them, is the expected angular velocity value of the previous beat, is the control period, is the adjustment coefficient.

12. The method according to claim 1, wherein The control moment gyro group is installed in a pyramid configuration, including five control moment gyros, and its combined torque is calculated by the following formula: ； Among them, Each column of the matrix is the angular momentum direction of the high-speed rotation of the control moment gyro when the gimbal angle of the control moment gyro is at 90°; Each column of the matrix is the angular momentum direction of the high-speed rotation of the control moment gyro when the gimbal angle of the control moment gyro is at the zero position; is the Jacobi matrix of the torque output of the control moment gyro group; is the gimbal angle vector.

13. The method according to claim 12, characterized in that, The singular robust inverse command operation law is: ； Among them, is an adjustment coefficient.

14. The method according to claim 1, characterized in that, This method updates the nominal state and control input in real time through a rolling horizon control process, specifically including: measuring the current system state, initializing the nominal state, solving the optimization problem, executing the control, and updating the system state.

Images