An agile satellite attitude control method for river imaging

Abstract

This invention discloses an agile satellite attitude control method suitable for river imaging. First, the time for each attitude maneuver and the initial and terminal imaging times for each imaging strip are planned according to the imaging mission. A broken line strip is generated, and the attitude, angular velocity, and angular acceleration information during the maneuver are calculated. System state variables are then defined, and a discretized model of the attitude system is established. The LQR problem is solved to obtain a feedback gain matrix, and the minimum robust invariant set is calculated based on the feedback gain matrix. Finally, an optimization problem is solved to generate an optimal nominal control sequence, and a command torque is output in conjunction with a feedforward torque compensation calculation module. Finally, the frame angle command angular velocity of the control torque gyro group is calculated based on the command torque, and torque output is achieved through a singular robust inverse command operation law. This method significantly improves the imaging process of meandering trajectories in river basins, enabling rapid maneuvering and stabilization of the entire satellite, resulting in better imaging results.

Description

Technical Field

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

[0002] A single remote sensing satellite can cover thousands of kilometers of river basins, enabling simultaneous collection of hydrological parameters across the entire basin and large-scale real-time monitoring. This is of great significance for environmental protection, disaster warning, and water resource management. Because high-resolution satellites have a relatively small swath width (approximately 10-14 km), agile satellites are used to rapidly adjust attitude using control torque gyros for large-scale imaging. However, large actuators produce greater torque noise when outputting large torques. Furthermore, during high-angular-rate maneuvers, the measurement errors of satellite attitude measurement sensors, such as star sensors and gyros, deteriorate. Continuous maneuvers are required for imaging of winding rivers, making it difficult to achieve effective attitude control and impacting remote sensing imaging. Traditional proportional-integral-derivative (PID) control takes a long time to transition from attitude maneuvers to stable imaging. Standard model predictive control methods, which rely on precise model conditions, can also compromise attitude control effectiveness in the presence of large external disturbances, making robustness a challenge. Summary of the Invention

[0003] The purpose of the present invention is to overcome the shortcomings of the existing technology and propose an agile satellite attitude control method suitable for river imaging, so as to solve the technical problem of fast, stable and robust imaging of winding river channels. This method greatly improves the imaging process of winding trajectories in river basins. The entire satellite can maneuver quickly and stabilize quickly, and better imaging effects can be obtained.

[0004] In order to achieve the above object, the technical solutions specifically adopted by the present invention are as follows:

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

[0006] Step 1: Imaging trajectory planning: According to the imaging task, the time of each attitude maneuver and the imaging initial and terminal moments of each imaging strip are planned, a broken line strip is generated, and the attitude, angular velocity, and angular acceleration information during the maneuver are calculated;

[0007] Step 2: Design a model predictive controller. First, establish the spacecraft attitude dynamics and kinematics model. Then, define the system state variables and establish a discretized model of the attitude system. Solve the LQR problem to obtain the feedback gain matrix. Then, calculate the minimum robust invariant set based on the feedback gain matrix. Finally, solve the optimization problem to generate the optimal nominal control sequence and combine it with the feedforward torque compensation calculation module to output the command torque.

[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 step 1 specifically includes:

[0010] Step 1.1: Obtain 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 latitude and longitude, altitude, and length 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;

[0012] Step 1.3: Use the seventh-order polynomial trajectory planning method to determine 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 step 1.2 is as follows: according to the inertia of the satellite and the maneuverability of the control moment gyro, a relationship table of maneuvering time and maneuvering angle can be calculated in advance, and the duration of each maneuver can be calculated by interpolation of the attitude angle required for maneuvering; 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, and the initial moment and terminal moment of the imaging of each strip can be confirmed in turn, so that the attitude, angular velocity and angular acceleration information of each maneuvering time and the initial moment and terminal moment of the imaging of each imaging strip can be calculated; at the starting moment of the strip and terminal moment Two sets of orbital elements are taken near the two points, and the values ​​of the attitude, angular velocity and angular acceleration at the initial and terminal moments are obtained through iterative calculation.

[0014] Preferably, the step 2 further includes rewriting the discretization model of the posture dynamics as follows:

[0015] ;

[0016] ;

[0017] ;

[0018] In the formula, k represents the current moment, k+1 represents the next moment, and Refers to the system state variables at the next moment and the current moment respectively, , represents the impact of environmental interference on the system dynamics at each time step k, i.e., disturbance; is the upper limit of disturbance, is the lower limit of disturbance.

[0019] As a preference, in the rewritten discretization model, both the system state variables and the system control variables are constrained, wherein the constraints on the system control variables are subject to the angular momentum sphere envelope and the maximum torque, and the constraint formula is:

[0020] ;

[0021] ;

[0022] In the formula 、 is the constraint matrix, By the identity matrix Stacking to form 12 6 matrix, By the identity matrix Stack formation Matrix and are the upper and lower limits of the state variable, respectively. and are the upper and lower limits of the control variable, respectively.

[0023] As a preference, in step 2, the method for solving the LQR problem to obtain the feedback gain matrix is: within the prediction time domain N, the system state can be decomposed into the nominal part and error ,error Indicates the actual status Relative to nominal state At time step The subsequent deviation is ,in is the nominal control input, the feedback gain matrix Offline design makes the system quadratically stable and the feedback gain matrix The solution is the classic LQR problem that satisfies , where T is the device, given 、 as well as and Then we can directly solve , Corresponding to 6 state variables, namely attitude and angular velocity; is the control output weight matrix , reflecting 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 , calculate the minimum robust invariant set by the Minkowski and calculation method:

[0026] ; ;

[0027] Calculate the Minkowski sum,

[0028] When stopped, finally .

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

[0030] ;

[0031] in, Calculate the Pontryagin set difference.

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

[0033] ;

[0034] in, and is the nominal state and input, is the nominal state at the end of the prediction time domain; 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 is calculated by the following formula:

[0037] ;

[0038] in, is the expected value of the angular velocity of the previous beat, To control the cycle, is the adjustment factor.

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

[0040] ;

[0041] in, Each column of the matrix is ​​the direction of the angular momentum of the control moment gyro when the frame angle of the control moment gyro is at 90°; Each column of the matrix is ​​the direction of the angular momentum 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; is the frame angle vector.

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

[0043] ;

[0044] in, is the adjustment factor.

[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 control and updating the system state.

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

[0047] This paper proposes a method for agile satellite attitude control using a tube model predictive control. This method uses offline design of the feedback gain matrix, while online optimization only requires optimizing the nominal trajectory. Disturbance effects are compensated in real time through feedback, significantly reducing computational complexity. When processing bounded disturbances, robust positive invariant sets are designed to ensure that the closed-loop system remains asymptotically stable under disturbances. This effectively addresses the technical challenges of fast, stable, and robust imaging of winding river channels. This method significantly improves the imaging process for winding trajectories in river basins, enabling rapid maneuverability and stability of the entire satellite, resulting in superior imaging results. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

[0050] Figure 3 This is the simulation result of the posture control process in the meandering waters of the Yangtze River according to the present invention. DETAILED DESCRIPTION

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

[0052] An agile satellite attitude control method suitable for river imaging, such as Figure 1 As shown, the following steps are included:

[0053] Step 1: Imaging trajectory planning: According to the imaging task, the time of each attitude maneuver and the initial and terminal imaging moments of each imaging strip are planned, a broken line strip is generated, and the attitude, angular velocity, and angular acceleration information during the maneuver are calculated.

[0054] In this embodiment, Figure 2 The broken line belt is generated in the area of ​​26 degrees to 28 degrees north latitude where the Yangtze River basin meanders.

[0055] The specific steps include the following:

[0056] Step 1.1: First, obtain mission requirements, satellite orbit information, and ground station parameters to determine the maximum available time window for the observation area. This example assumes that the observation time and imaging angle are already confirmed. Based on the key reference points in the Yangtze River imaging area, calculate the azimuth between the first two points as the base azimuth.

[0057]

[0058] In the formula 、 are the starting and ending latitudes, are the azimuths of the start and end points, 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 of the two points in the current segment and calculate the absolute deviation from the reference azimuth: ,like If the value is greater than the threshold, which is set to 5 degrees in this embodiment, the merging stops. If the value is less than the threshold, the merging continues until the last point is processed. In this way, multiple broken line strips can be obtained, recording the latitude and longitude, altitude, and strip length of the starting and ending points of each strip.

[0060] 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.

[0061] Specifically, based on the satellite's inertia and the maneuverability of the control moment gyro, a table of maneuver time and maneuver angle can be pre-simulated (as shown in Table 1), and the duration of each maneuver can be calculated by interpolating the attitude angle required for the maneuver. 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, and the initial and terminal moments of each strip's imaging can be confirmed in turn. In this way, the attitude, angular velocity, and angular acceleration information of each maneuver time and the initial and terminal moments of each imaging strip can be calculated. At the starting moment of the strip and terminal moment Two sets of orbital elements are taken near two points and used for interpolation in the entire calculation process.

[0062]

[0063] Table 1 Attitude maneuver time correspondence table

[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 values ​​of the initial attitude, angular velocity and angular acceleration are obtained through iterative calculation. Similarly, the target point is changed to the sub-satellite target point at the terminal time. Definition Temporary iteration time, in whole seconds; for The Greenwich sidereal hour angle at that moment.

[0066] In this embodiment, the iteration condition is Terminal moment - iterative calculation start time can also get the values ​​of attitude, angular velocity and angular acceleration at the terminal moment of the strip. The calculation method is as follows:

[0067] Attitude motivation moment - iterative calculation start time);

[0068] ;

[0069] ;

[0070] ;

[0071] and are the angle and angular velocity values, and Temporarily save the angle and angular velocity values ​​for the previous beat, which will be used for differential calculation later. ; , is the imaging duration Where for The distance between the center of the earth and point P on the strip at any moment, for The distance between the center of the Earth and the starting point of the strip at that moment, for The distance between the center of the Earth and the end point of the strip at that moment.

[0072] ; for The center of the earth is always expanding, is the unit vector from the center of the Earth to the starting point of the strip, Perpendicular to , , Here is the unit vector pointing from the center of the Earth to the end point of the strip:

[0073] Unit vector of P on the strip pointing to the center of the Earth ;

[0074] The vector of P on the strip pointing to the center of the earth ;

[0075] According to the orbital extrapolation results Get the number of orbital elements and interpolate to get They are The semi-major axis, eccentricity, orbit inclination, latitude argument, right ascension of the ascending node, true anomaly, distance from the center of the Earth, and components of the unit vector pointing from the center of the Earth to the satellite in the Earth-fixed coordinate system at this moment.

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

[0077] Coordinate transformation matrix from ground-fixed system to orbital system ; is the coordinate transformation matrix from the earth-fixed system to the inertial system is the coordinate transformation matrix from inertial system to orbital system;

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

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

[0080] Component of the linear velocity of the pointing point in the orbital coordinate system ;

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

[0082] Modulus ;

[0083] Roll angle ;

[0084] Pitch angle ;

[0085] Orbital system to reference attitude system coordinate transformation matrix yes The matrix is ​​generated by 123 transposition;

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

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

[0088] ;

[0089] ;

[0090]

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

[0092] The component 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: Plan a trajectory based on the attitude values ​​at the initial and final imaging moments, and then determine whether the maneuverability will exceed the tolerance during the trajectory. This embodiment uses a seventh-order polynomial trajectory planning method to determine whether the angular momentum exceeds the angular momentum envelope threshold of the control moment gyro during the maneuver. If so, the maneuver is terminated; otherwise, subsequent attitude maneuvers are performed.

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

[0099]

[0100] ;

[0101] in, is the duration of the attitude maneuver; are polynomial coefficients, where , Can be set as three-axis attitude angle The calculation formula is consistent, the subscript 0 represents the starting attitude angle of the attitude maneuver, and the subscript is represents the attitude angle of the attitude maneuver terminal; is the three-axis attitude angle at the initial moment of attitude maneuver, is the three-axis attitude angle at the end of attitude maneuver.

[0102] According to the polynomial coefficients, calculate Target posture at all times , angular velocity , angular acceleration :

[0103] ;

[0104] Where t is relative to the start time of attitude maneuver, 0≤t≤tm, Can be set as three-axis attitude angle and roll angle respectively , pitch angle , yaw angle ,Right now , , .

[0105] The three-axis attitude is calculated iteratively by the set maneuvering process time ( is the pitch angle, is the roll angle, is the yaw angle) and Euler angular velocity:

[0106] (3)

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

[0108]

[0109] in is the value of the principal axis of inertia of the entire star, The angular momentum of the satellite during the maneuver is calculated to determine whether it exceeds the threshold set by the control moment gyro angular momentum envelope. If so, it indicates that the angular velocity maneuverability cannot meet the maneuver time requirements and subsequent attitude maneuvers are terminated. If not, subsequent attitude maneuvers can be performed.

[0110] Step 2: Design a tube model predictive controller. First, establish the spacecraft attitude dynamics and kinematics model; then define the 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 based on the feedback gain matrix. Finally, solve the optimization problem to generate the optimal nominal control sequence, and combine it with the feedforward torque compensation calculation module to output the command torque.

[0111] Understandably, Figure 1 As shown, the satellite attitude control system includes a maneuvering target attitude determination module, a feedforward torque compensation calculation module, and a tube model predictive control module. The maneuvering target attitude determination module provides real-time target attitude angle, angular velocity, and angular acceleration based on imaging requirements. The feedforward torque compensation calculation module calculates the feedforward torque. The tube model predictive control module calculates the optimal control sequence based on the tube model predictive control principle. The control torque is output to the satellite through the control torque gyro group after superimposing the control variable with the feedforward torque to control the satellite's attitude. The attitude measurement and determination module calculates the satellite's attitude and angular velocity in real time using onboard gyros and star sensors.

[0112] The target attitude reference trajectory of the maneuver process is determined using the polynomial trajectory planning calculation in step 1.3, based on the three-axis attitude angle at the initial moment of imaging. , angular velocity and angular acceleration , target attitude angle at the end of imaging , angular velocity and angular acceleration , and the time deviation relative to the starting time is used to calculate the expected three-axis attitude angle at each moment through the polynomial formula , angular velocity and angular acceleration After the target attitude is generated, the attitude is first converted into the attitude conversion matrix from the target pointing to the orbit system through the 123 conversion sequence. , calculate the attitude transformation matrix from the orbital system to the inertial system through the orbital data , transform the posture matrix and the attitude transformation matrix Converting from Euler angles to reference trajectory quaternion, the reference trajectory angular velocity is:

[0113]

[0114] The specific steps include the following:

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

[0116] ;

[0117]

[0118] Where, is the inertia tensor of the entire star; is the spacecraft's three-axis angular velocity; represents the satellite's three-axis control torque, is the actuator angular momentum; The torque is generated by changing the angular momentum of the actuator, where the actuator is a group of control torque gyroscopes; is the interference torque. Euler angle is often used in engineering A motion model describing the satellite's attitude, , are the roll angle, pitch angle, and yaw angle respectively.

[0119] Step 2.2: Define system state variables (6 1), establish the discrete model of the posture system as follows:

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

[0121]

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

[0123]

[0124] The state variables and control variables of the system are subject to constraints. The constraints of the control variables are mainly subject to the angular momentum ball envelope and the maximum torque. The constraint formula is:

[0125] ;

[0126] ;

[0127] In the formula 、 is the constraint matrix, By the identity matrix Stacking to form 12 6 matrix, By the identity matrix Stack formation Matrix and are the upper and lower limits of the state variable, respectively. and are the upper and lower limits of the control variable, respectively.

[0128] Step 2.3, solve the LQR problem and get and :

[0129]

[0130] in, 、 and is the weight matrix; is the reference trajectory.

[0131] In the prediction time domain N, the system state can be decomposed into the nominal part and error ,error Indicates the actual status Relative to nominal state At time step The subsequent deviation is ,in is the nominal control input, the feedback gain matrix The offline design makes the system quadratically stable, which is obtained by solving the LQR.

[0132] Feedback gain matrix The solution is the classic LQR problem that satisfies , given 、 as well as and Then we can directly solve , Corresponding to 6 state variables, namely attitude quaternion and three angular velocities; controlling the output weight matrix , reflecting the torque limitation.

[0133] Step 2.4, solve the minimum robust invariant set:

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

[0135] ; ;

[0136] Calculate the Minkowski sum, Represents the i-th power of the closed-loop state matrix. If the calculation is simplified, it can be considered = is an invariant;

[0137] Using iterative calculation method: ;

[0138] When stopped, finally .

[0139] Constrained Calculation

[0140] ,

[0141] in, Calculate the difference between the Pontryagin sets .

[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 optimization problem of the model predictive controller is:

[0144] ;

[0145] in, and is the nominal state and input, is the nominal state at the end of the prediction time domain; 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) Measurement state: obtaining the real system state;

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

[0151] (3) Solving optimization problems: Minimizing the objective function ;

[0152] (4) Execution control: ;

[0153] (5) Update status: ;

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

[0155] Step 2.6: Combine the feedforward torque compensation calculation module to output the command torque:

[0156] The feedforward torque compensation is calculated using the following formula:

[0157] ;

[0158] in, is the expected value of the angular velocity of the previous beat, To control the cycle, is the adjustment factor.

[0159] Based on the above, the feedforward model predictive control law is:

[0160] ;

[0161] The command torque of the satellite is used to output to the control torque gyro group; is the feedforward control quantity at the kth moment; is the control quantity predicted by the control model at the kth moment.

[0162] 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.

[0163] It can be understood that the frame angle command angular velocity of the control moment gyro group is calculated based on the three-axis command torque of the satellite body obtained in step 2.

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

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

[0166] ;

[0167] in, is the CMGs frame angle vector matrix, is the frame angle of the i-th control moment gyro; and is the matrix related to the installation orientation of CMGs, Each column is the angular momentum direction of each CMG when the frame angle is at 90°. Each column is the direction of the angular momentum of each CMG when its frame angle is at zero. is an n×1 unit vector.

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

[0169]

[0170] in, Each column of the matrix is ​​the direction of the angular momentum of the control moment gyro when the frame angle of the control moment gyro is at 90°; Each column of the matrix is ​​the direction of the angular momentum 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 operating rate of the control moment gyro group is:

[0172]

[0173] in, is the frame angle command angular velocity; It is a singular robust inverse instruction operation law; It is the zero motion singularity avoidance operation law; is the nominal frame angle regression operation law;

[0174] It should be noted that there are multiple implementations of the zero-motion singularity avoidance operation law and the frame angle regression operation law, which are well-known. This embodiment mainly describes the singularity robust inverse instruction operation law, which is as follows:

[0175] Singular robust inverse instruction operation law:

[0176] ;

[0177] is the adjustment coefficient; according to the extreme value conditions, we can know that: .

[0178] According to the implementation process provided in this embodiment, the pipe model prediction control effect diagram in the embodiment of the present invention is obtained. In the imaging of the winding Yangtze River Basin in this example, a complex maneuver 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 maneuver quickly and stabilize quickly, and give full play to the satellite's performance.

[0180] The above shows and describes 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 to the above embodiments. The above embodiments and descriptions are merely preferred examples of the present invention and are not intended to limit the present invention. Various changes and improvements may be made to the present invention without departing from the spirit and scope of the present invention. Such changes and improvements fall within the scope of the present invention. The scope of protection claimed in the present invention is defined by the appended claims and their equivalents.

Claims

1. An agile satellite attitude control method suitable for river imaging, characterized in that: The steps include: Step 1: Imaging trajectory planning: According to the imaging task, the time of each attitude maneuver and the imaging initial and terminal moments of each imaging strip are planned, a broken line strip is generated, and the attitude, angular velocity and angular acceleration information during the maneuver are calculated to determine the subsequent attitude maneuver; Step 2: Design a model predictive controller. First, establish the spacecraft attitude dynamics and kinematics model. Then, define the system state variables and establish a discretized model of the attitude system. Solve the LQR problem to obtain the feedback gain matrix. Then, calculate the minimum robust invariant set based on the feedback gain matrix. Finally, solve the optimization problem to generate the optimal nominal control sequence and combine it with the feedforward torque compensation calculation module to output the command torque. The objective function of the tube model predictive controller optimization problem is: ； in, and is the nominal state and control input, is the nominal state at the end of the prediction time domain; The optimal nominal control sequence The control action at this moment, is the control output weight matrix, Indicates the nominal state, Indicates the kth moment, A refers to the system state transfer matrix; B refers to the control input matrix, Indicates the length of the prediction time domain, which is an integer; represents the state weight matrix; represents the solution matrix of the Lyapunov equation; represents the state constraint set after tightening; represents the set of control constraints after tightening; is transposed; represents the objective function of the optimization problem of the tube model predictive controller; The step 2 also includes rewriting the posture dynamics to obtain a discretized model: ； ； ； In the formula, k+1 represents the next moment, and Refers to the system state variables at the next moment and the current moment respectively, represents the system control variable, , It represents the influence of environmental interference on the system dynamics at each k-th moment, that is, the disturbance; is the upper limit of disturbance, ; Refers to the rewritten discretized model; It represents the space composed of 6-dimensional real vectors; Represents the space composed of 7-dimensional real vectors; Represents a block matrix formed by stacking two 6-order unit matrices on top of each other; The method for solving the LQR problem and obtaining the feedback gain matrix is: error Indicates the actual status Relative to nominal state The deviation after the kth moment is ,in is the nominal control input, the feedback gain matrix Offline design makes the system quadratically stable and the feedback gain matrix Solution , where T is the transpose, given 、 as well as and Then you can directly solve , Corresponding state variables; is the control output weight matrix , reflecting the limitation of the moment; Indicates nominal state At time step Control variables at time; The solution method of the minimum robust invariant set is: According to the calculated feedback gain matrix , calculate the minimum robust invariant set by the Minkowski and calculation method: ； ； Calculate the Minkowski sum, represents the i-th power of the closed-loop state matrix; represents the closed-loop state matrix; When stopped, finally , in, represents the number of iterations; i represents the power number; represents the T+1 step robust invariant set; represents the T-step robust invariant set; It means that after multiple iterations, it will eventually converge to the value, and the number of iterations tends to infinity. In actual use, it will converge quickly after a finite number of iterations. The optimal nominal control sequence satisfies: ; ; in, Indicates nominal state System state variables at time k; represents the optimal nominal control input obtained through optimization at the kth moment; It means the first k-th moment Interference value; It represents the optimal nominal disturbance compensation at the kth moment; Represents the system state variables at the current moment; Indicates nominal state; 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 agile satellite attitude control method for river imaging according to claim 1, characterized in that: The imaging trajectory planning in step 1 specifically includes: Step 1.1: Obtain 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 latitude and longitude, altitude, and length of each strip; 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; Step 1.3: Use the seventh-order polynomial trajectory planning method to calculate the angular momentum based on the maneuver time and the corresponding attitude, angular velocity, and angular acceleration information. Then, determine whether the angular momentum exceeds the angular momentum envelope threshold of the control moment gyro during the maneuver. If so, terminate the maneuver; otherwise, perform subsequent attitude maneuvers.

3. The agile satellite attitude control method for river imaging according to claim 2, characterized in that: The specific method of step 1.2 is as follows: based on the inertia of the satellite and the maneuverability of the control moment gyro, a relationship table between maneuvering time and maneuvering angle is pre-calculated, and the duration of each maneuver is calculated by interpolating the attitude angle required for the maneuver from the corresponding relationship table; The imaging time of each strip is obtained by dividing the length of the imaging strip by the moving speed of the subsatellite point. The initial and terminal moments of each strip imaging are confirmed in turn, that is, the maneuvering time and the attitude, angular velocity and angular acceleration information of the initial and terminal moments of each imaging strip imaging are calculated. At the starting moment of the strip and terminal moment Two sets of orbital elements are taken near the two points, and the values ​​of the attitude, angular velocity and angular acceleration at the initial and terminal moments are obtained through iterative calculation.

4. The agile satellite attitude control method for river imaging according to claim 1, characterized in that: In the discretization model, both the system state variables and the system control variables are constrained. The system state variable constraint formula is: ； The system control variable constraint formula is: ; In the formula 、 is the constraint matrix, By the identity matrix Stacking to form 12 6 matrix, By the identity matrix Stack formation Matrix and are the upper and lower limits of the state variables respectively; and are the upper and lower limits of the control variable, Represents the system state variables at the current moment, represents the system control variable, is a parameter related to the upper and lower limits of the system state variables, that is, a threshold vector , Is a parameter related to the upper and lower limits of the control variable, that is, a threshold vector , Represents the system state variable constraints, Represents the system control variable constraints.

5. The agile satellite attitude control method for river imaging according to claim 4, characterized in that: When calculating the minimum robust invariant set, a constraint tightening calculation is introduced: ； in, For the Pontryagin set difference calculation, represents the system state variable constraints, represents the set of control constraints after tightening, represents the system control variable constraints, represents the feedback gain matrix, Represents the set of state constraints after tightening.

6. The agile satellite attitude control method for river imaging according to claim 4, characterized in that: The feedforward torque compensation is calculated by the following formula: ； in, is the expected value of the angular velocity of the previous beat, To control the cycle, is the adjustment coefficient, Indicates the value of the principal axis of inertia of the entire satellite, Indicates the expected value of the current angular velocity.

7. The agile satellite attitude control method for river imaging according to claim 1, characterized in that: The control moment gyro group is installed in a pyramid configuration and includes five control moment gyros. The resultant torque is calculated by the following formula: ； in, Each column of the matrix is ​​the direction of the angular momentum of the control moment gyro when the frame angle of the control moment gyro is at 90°; Each column of the matrix is ​​the direction of the angular momentum 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; is the frame angle vector; represents the control moment gyro; represents the angular momentum amplitude of the control moment gyro; It represents the change of the composite angular momentum of the control moment gyro group, which is also the composite torque; represents the first derivative of the frame angle vector.

8. The agile satellite attitude control method for river imaging according to claim 7, characterized in that: The singular robust inverse instruction operation law is: ； in, is the adjustment coefficient, is the Jacobi matrix of the torque output of the control moment gyro group; Transpose operation of the Jacobi matrix for the torque output of the control torque gyro group; Represents the singular robust inverse instruction operation law.

9. The agile satellite attitude control method for river imaging according to claim 1, characterized in that: 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 control and updating the system state.

Images