Thermally induced deformation prediction and optimization control method for spaceborne antenna

Abstract

The invention discloses a thermally induced deformation prediction and optimization control method for a spaceborne antenna, and belongs to the technical field of space antenna control. The method solves the problem that an existing control method is low in control precision due to the fact that an existing method cannot effectively cope with antenna thermal deformation caused by a non-constant thermal field. According to the method, the orbit prediction model of the space heat flow is established, the heat flow borne by the whole antenna at different orbit positions and time points can be accurately calculated, the thermal deformation rule and the key stress area of the antenna structure are comprehensively analyzed according to the thermally induced deformation characteristic caused by a time-varying thermal field, the position distribution layout of the actuator is optimally designed, and the space heat flow prediction accuracy is improved. According to the invention, the thermal deformation peak value area and the structure easily affected part are preferentially covered, so that more efficient shape control and vibration suppression are realized, the antenna is ensured to maintain a stable structural form and high-precision working performance in a complex thermal environment, and the control precision is effectively improved. The method can be applied to thermally induced deformation prediction and optimization control of the antenna.

Description

Technical Field

[0001] The present invention belongs to the technical field of space antenna control, and in particular relates to a method for predicting and optimizing the thermal deformation of a space-borne antenna. Background Art

[0002] Spaceborne antennas are a core component of modern aerospace technology, widely used in deep space exploration, global communications, Earth observation, and other fields. They are of great strategic significance for the efficient utilization of space resources and the exploration of the mysteries of the universe. As critical aerospace infrastructure, spaceborne antennas play an irreplaceable role, and their performance directly affects the success or failure of space missions. In recent years, with the continuous growth of space mission requirements, spaceborne antennas are moving towards larger sizes, more complex structures, and distributed configurations. This trend not only improves the communication capabilities and mission adaptability of antennas, but also introduces many new technical challenges.

[0003] At the same time, as the size of satellite-borne antenna structures increases, the impact of the space environment on antenna deformation becomes increasingly significant. Ultra-large satellite-borne antennas, due to their large surface area, make the phenomenon of thermally induced vibration more pronounced. Therefore, how to ensure the high stability and high-precision control of ultra-large satellite-borne antennas under thermal space environment conditions has become a key issue that needs to be urgently addressed in the current field of aerospace technology. For small satellite-borne antennas with an aperture of no more than 10 meters, some advanced thermal-structural coupling analysis methods have been developed. For the analysis of these thermal-structural coupling deformations, some passive and active suppression methods have also been developed. However, these studies ignore the impact of solar radiation evolving over time. For thermal deformation caused by non-constant thermal fields, this will bring additional challenges to antenna shape control and vibration suppression.

[0004] In summary, since the existing control method still cannot effectively deal with the thermal deformation of the antenna caused by the non-constant thermal field, the control accuracy of the existing control method is still low. It is very necessary to propose a new method to solve the above problem. Summary of the Invention

[0005] The purpose of the present invention is to solve the problem that existing methods cannot effectively cope with antenna thermal deformation caused by non-constant thermal fields, resulting in low control accuracy of existing control methods, and propose a method for predicting and optimizing the thermal deformation of satellite-borne antennas.

[0006] The technical solution adopted by the present invention to solve the above technical problems is: a method for predicting and optimizing the thermal deformation of a space-borne antenna, the method specifically comprising the following steps:

[0007] Step 1: Establish the geocentric equatorial inertial coordinate system and the orbital coordinate system to obtain the position vector of the sun in the geocentric equatorial inertial coordinate system at each moment;

[0008] Step 2: Calculate the heat flux per unit area of ​​solar radiation on the satellite orbit based on the sun position vector;

[0009] Step 3: Use the finite element method to build a satellite antenna model. Then, perform a thermal deformation analysis on the satellite antenna based on the heat flux per unit area of ​​solar radiation. Based on the analysis results, determine the optimal placement of the distributed control actuators on the satellite antenna.

[0010] Step 4: Establish an objective function according to the position of the actuator arrangement, and perform shape optimal distributed model predictive control based on the objective function.

[0011] Furthermore, the geocentric equatorial inertial coordinate system is specifically:

[0012] The origin E of the geocentric equatorial inertial coordinate system is located at the center of the earth, z i The axis is perpendicular to the equatorial plane and z i The positive direction of the axis points to the North Pole; i axis and y i The axes are all located in the equatorial plane, where x i The positive direction of the axis points to the vernal equinox, x i Axis, y i axis and z i The axes form a right-handed coordinate system.

[0013] Furthermore, the origin O of the orbital coordinate system is located at the center of mass of the spacecraft, z o The positive direction of the axis points to the center of the earth, x o Axis and z o The axes are vertical and x o The x axis is in the orbital plane of the spacecraft. o The positive direction of the y axis points to the direction of satellite flight, o The axis is perpendicular to the orbital plane of the spacecraft, x o Axis, y o axis and z o The axes form a right-handed coordinate system.

[0014] Furthermore, the position vector of the sun in the geocentric equatorial inertial coordinate system is:

[0015] The transformation matrix M from the orbital coordinate system to the geocentric equatorial inertial coordinate system is:

[0016]

[0017] Where Ω is the right ascension of the ascending node, i is the orbital inclination, μ is the argument of latitude, μ = ω + θ, ω is the argument of perigee, and θ is the true anomaly.

[0018] α s=arctan(sinεsinλ s / cosλ s ), δ s =arcsin(sinεsinλ s ) (2)

[0019] where ε is the inclination of the ecliptic, λ s is the ecliptic longitude, α s is the right ascension, δ s is declination;

[0020] The position component of the sun in the geocentric equatorial inertial coordinate system (x i ,y i ,z i )for:

[0021] x i =r s cosα s cosδ s , z i =r s sinα s cosδ s ,y i =r s sinδ s (3)

[0022] Among them, r s It represents the distance from the center of the sun to the center of the earth.

[0023] Furthermore, the specific process of step 2 is as follows:

[0024] Step 2.1 Calculate the angle β between the satellite and the earth e :

[0025]

[0026] Among them, R e is the radius of the Earth, R sat is the distance from the center of gravity of the satellite antenna to the center of the Earth;

[0027] The angle between the line connecting the center of gravity of the satellite antenna and the center of the earth and the line connecting the center of gravity of the satellite antenna and the center of the sun is denoted as β. Compare β with β e Size:

[0028] (1) When β is greater than β e When the satellite antenna is illuminated, the solar irradiance S c for:

[0029]

[0030] in, represents the average solar irradiance; R s is the distance from the center of the sun to the center of the earth; is the average distance between the center of the sun and the center of the earth;

[0031] (2) When β is less than or equal to β e When the satellite antenna is in the shadow area, the solar irradiance S c =0;

[0032] Step 2.2. Solar radiation heat flux per unit area q s for:

[0033]

[0034] where n is the normal of the antenna surface mesh in the orbital coordinate system.

[0035] Furthermore, the specific process of step three is:

[0036] The parabolic model of the satellite antenna is established using the commercial software Abaqus. The heat flux per unit area calculated in step 2 is then loaded using the Load module in Abaqus to obtain the deformation of each area of ​​the satellite antenna. The optimal position of the distributed control actuator is then determined based on the deformation of each area of ​​the satellite antenna.

[0037] Furthermore, the specific process of step 4 is as follows:

[0038] Step 4.1. The structural dynamics equation of the satellite-borne antenna is:

[0039]

[0040] Where M, C′ and K are the inertia matrix, damping matrix and stiffness matrix of the satellite antenna respectively, B′ is the actuator arrangement matrix, and F c is the cable control force vector in physical space, u d is the force generated by thermal radiation in physical space, x represents the vector composed of the coordinates of the nodes where each actuator is arranged, represents the first-order derivative of x, represents the second derivative of x;

[0041] Perform modal coordinate transformation on the structural dynamics equations of the spaceborne antenna:

[0042] x=Φη (8)

[0043] Where η is the N-dimensional modal coordinate vector of the flexible structure, Φ is the first N-order regular modal matrix;

[0044] Then the form of the distributed cable dynamic equation is:

[0045]

[0046] in, is the second-order derivative of η, is the first-order derivative of η, Λ a is the frequency matrix, ξ is the damping coefficient matrix, B c is the distributed cable control matrix, Φ T is the transpose of Φ;

[0047] Frequency matrix Λ a The damping coefficient matrix ξ is in the form of:

[0048]

[0049] Where, ω a1 、ω a2 ,…,ω aN Determined by formula (11):

[0050]

[0051] Where r = 1, 2, ..., N, Φ r is the first r-order regular mode matrix, is the rth order main frequency corresponding to the structural vibration, and:

[0052]

[0053] Then the damping coefficient ξ r Determined by formula (13):

[0054]

[0055] Among them, a0 and a1 are constants;

[0056] Step 42: Convert equation (9) into the system state equation:

[0057]

[0058] y=CX (15)

[0059] Where, is the first-order derivative of X, u represents the system control quantity, y represents the system output, and I is the unit matrix;

[0060]

[0061] Distributed cable control matrix B cThe solution is related to the force on the cable, and the control force at both ends of the cable is decomposed into three axes:

[0062]

[0063] Where, (x l ,y l ,z l ) represents the position of the l-th actuator node, l = 1, 2, …, L′, L′ represents the total number of actuators arranged;

[0064] Then the distributed cable control matrix B of the lth cable is cl for:

[0065]

[0066] Among them, φ xi 、φ yi and φ zi They are the three degrees of freedom corresponding to the i-th order main mode, i=1,2,…,N,B c =[B c1 ,B c2 ,…,B cL′ ];

[0067] Step 4.3. Convert the system state equations of Equation (14) and Equation (15) into discrete state space equations:

[0068]

[0069] Among them, X(k|k) represents the measured value of the state variable at time k, X(k+1|k+1) represents the measured value of the state variable at time k+1, u(k) represents the control quantity at time k, η(k+1|k+1) represents the output quantity at time k+1, and u d (k) represents the force generated by thermal radiation at time k in physical space, and the discrete state matrices H(k), F(k) and G(k) are respectively:

[0070] F(k)=I+AT s ,G(k)=BT s ,H(k)=DT s (20)

[0071] Where, T s is the sampling time, I is the identity matrix;

[0072] The state variables in the prediction domain are:

[0073]

[0074] Where N pis the prediction time domain; X(k+N p |k) is the system at time k, the prediction time domain N p Internal system state quantity; η(k+N p |k) is the system at time k, the prediction time domain N p The internal system predicts the output; u(k+j|k) is the prediction of the control quantity at time k+j at time k; u d (k+j|k) is the prediction of the force generated by thermal radiation at time k+j at time k;

[0075] Rewrite Equation (21) into matrix form:

[0076] E(k+1)=Ψ(k)X(k|k)+Θ(k)U(k)+Ξ(k)U d (k) (23)

[0077] Γ(k+1)=L(k)E(k+1) (24)

[0078] In the formula, E(k+1), Ψ(k), Θ(k), Ξ(k), Γ(k+1), L(k), U(k) and U d The expression of (k) is as follows:

[0079]

[0080]

[0081] Construct the objective function J of the lth subsystem l (k) in the following form:

[0082] J l (k)=U l (k) T R l U l (k)+(Γ l (k+1)-Γ rl (k+1)) T Q l (Γ l (k+1)-Γ rl (k+1)) (31)

[0083] Among them, U l (k) represents the control quantity corresponding to the lth subsystem in U(k), R l represents the control constraint weight matrix corresponding to the lth subsystem; Γ l (k+1) represents the output modal coordinates of the lth subsystem in Γ(k+1); Γ rl (k+1) represents the target modal coordinate corresponding to the lth subsystem; Ql Represents the curve tracking weight matrix corresponding to the lth subsystem;

[0084]

[0085] Where, represents the optimal control quantity sequence of the lth subsystem, U max Indicates the maximum control amount;

[0086] The optimal output sequence u of the lth cable l (k) is expressed as:

[0087]

[0088] Among them, u(k)=[u1(k),u2(k),…,u L′ (k)].

[0089] Furthermore, the force generated by thermal radiation in the physical space is:

[0090]

[0091] Among them, q s is the solar radiation heat flux per unit area, S flux is the sun-receiving area of ​​the antenna, and c is the speed of light.

[0092] The beneficial effects of the present invention are:

[0093] 1. This invention addresses the significant impact of spatial heat flux variations on very large satellite antennas due to their enormous size. By developing an orbital prediction model for spatial heat flux, this model accurately calculates the heat flux experienced by the entire antenna at different orbital positions and time points. Compared to traditional static or empirical models, this method more accurately reflects the dynamic changes in heat flux, improving the accuracy of thermal environment prediction and the reliability of the model, thereby providing stronger support for antenna structural optimization and thermal control.

[0094] 2. This invention is applicable to low-orbit, medium-orbit, high-orbit, and even deep-space exploration missions, possessing excellent adaptability and broad application value. Its model parameters can be flexibly adjusted to match different antenna configurations, orbital environments, and mission requirements, while being unrestricted by antenna shape and size. Simply providing the antenna's initial state and time information allows for precise calculation of the heat flux at any moment, providing reliable thermal environment prediction support for various space missions.

[0095] 3. This invention addresses the thermally induced deformation characteristics caused by a time-varying thermal field, comprehensively analyzes the thermal deformation patterns and critical stress-bearing areas of the antenna structure, and optimizes the distribution and design of the actuators. By rationally configuring the actuator positions, they prioritize coverage of peak thermal deformation areas and vulnerable structural parts, thereby achieving more efficient shape control and vibration suppression. This design not only improves the actuator's control efficiency but also reduces unnecessary energy consumption, ensuring that the antenna maintains a stable structural form and high-precision operating performance in complex thermal environments, effectively improving control accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0096] Figure 1 This is a flow chart of a method for predicting and optimizing the thermal deformation of a spaceborne antenna according to the present invention;

[0097] Figure 2 It is a schematic diagram of two coordinate systems;

[0098] Figure 3 It is a schematic diagram of cable force analysis;

[0099] In the figure, T ax2 =-T ax1 ,T ay2 =-T ay1 ,T az2 =-T az1 ;

[0100] Figure 4 It is the curve of solar radiation heat flux changing in two days;

[0101] Figure 5(a) shows the amplitude of the total node displacement in the result of thermal deformation analysis of the antenna model using Abaqus;

[0102] -9.800e-03 means -9.800×10 -3 ;

[0103] Figure 5(b) shows the amplitude of displacement along the X-axis in the result of thermal deformation analysis of the antenna model using Abaqus;

[0104] Figure 5(c) shows the amplitude of displacement along the Y-axis in the result of thermal deformation analysis of the antenna model using Abaqus;

[0105] Figure 5(d) shows the amplitude of displacement along the Z axis in the result of thermal deformation analysis of the antenna model using Abaqus;

[0106] Figure 6 It is a schematic diagram of the actuator arrangement position;

[0107] Figure 7 is the control result of the first six modal coordinates;

[0108] Among them, (a) is the coordinate control result of the first-order mode, (b) is the coordinate control result of the second-order mode, (c) is the coordinate control result of the third-order mode, (d) is the coordinate control result of the fourth-order mode, (e) is the coordinate control result of the fifth-order mode, and (f) is the coordinate control result of the sixth-order mode.

[0109] Figure 8 is the comparison of deformation suppression results of the three selected nodes;

[0110] Among them, (a) is the deformation displacement at node 103, (a) is the deformation displacement at node 150, and (c) is the deformation displacement at node 166. DETAILED DESCRIPTION

[0111] Specific implementation method 1: Combination Figure 1 This embodiment describes a method for predicting and optimizing the thermal deformation of a space-borne antenna, the method specifically comprising the following steps:

[0112] Step 1: Establish the geocentric equatorial inertial coordinate system and the orbital coordinate system to obtain the position vector of the sun in the geocentric equatorial inertial coordinate system at each moment;

[0113] Step 2: Calculate the heat flux per unit area of ​​solar radiation on the satellite orbit based on the sun position vector;

[0114] Step 3: Use the finite element method to build a satellite antenna model. Then, perform a thermal deformation analysis on the satellite antenna based on the heat flux per unit area of ​​solar radiation. Based on the analysis results, determine the optimal placement of the distributed control actuators on the satellite antenna.

[0115] Step 4: According to the positions of the actuators, an objective function is established by minimizing the error between the actual surface shape and the ideal surface shape, and shape optimal distributed model predictive control (SOD-MPC) is performed according to the objective function.

[0116] Specific implementation method 2: Combination Figure 2 This embodiment is different from the first embodiment in that the geocentric equatorial inertial coordinate system is:

[0117] The origin E of the geocentric equatorial inertial coordinate system is located at the center of the earth, z i The axis is perpendicular to the equatorial plane and z i The positive direction of the axis points to the North Pole; i axis and y i The axes are all located in the equatorial plane, where x i The positive direction of the x-axis points to the vernal equinox (i.e., the direction of the line connecting the earth and the sun at the vernal equinox (around March 21)). iAxis, y i axis and z i The axes form a right-handed coordinate system.

[0118] Other steps and parameters are the same as those in the first embodiment.

[0119] Specific implementation method three: Combination Figure 2 The difference between this embodiment and the first or second embodiment is that the origin O of the orbital coordinate system is located at the center of mass of the spacecraft. o The positive direction of the axis points to the center of the earth, x o Axis and z o The axes are vertical and x o The x axis is in the orbital plane of the spacecraft. o The positive direction of the y axis points to the direction of satellite flight, o The axis is perpendicular to the orbital plane of the spacecraft, x o Axis, y o axis and z o The axes form a right-handed coordinate system.

[0120] Other steps and parameters are the same as those in the first or second embodiment.

[0121] In the geocentric equatorial inertial coordinate system, the position of the satellite can be expressed in rectangular coordinates x i ,y i and z i It can also be expressed in spherical coordinates r, α, and δ. In spherical coordinates, r is the distance from the satellite to the center of the Earth; α is the right ascension, which is the angle from the vernal equinox eastward to the projection of the radial vector r on the equatorial plane; δ is the declination, which is the angle from the equatorial plane to the radial vector r.

[0122] Specific embodiment 4: This embodiment differs from any one of specific embodiments 1 to 3 in that the position vector of the sun in the geocentric equatorial inertial coordinate system is:

[0123] The orbital coordinate system and the geocentric equatorial inertial coordinate system are related by the ascending node right ascension Ω, orbital inclination i, argument of perigee ω and true anomaly θ. The transformation matrix M from the orbital coordinate system to the geocentric equatorial inertial coordinate system is:

[0124]

[0125] Where Ω is the right ascension of the ascending node, i is the orbital inclination, μ is the argument of latitude, μ = ω + θ, ω is the argument of perigee, and θ is the true anomaly.

[0126] Due to the rotation of the Earth, the position of the Sun changes over time. Assuming the Sun's right ascension α s and declination δ s Remaining constant throughout the day, they can be expressed as follows:

[0127] α s =arctan(sinεsinλ s / cosλ s ), δ s =arcsin(sinεsinλ s ) (2)

[0128] where ε is the inclination of the ecliptic, λ s is the ecliptic longitude, α s is the right ascension, δ s is declination;

[0129] The position component of the sun in the geocentric equatorial inertial coordinate system (x i ,y i ,z i )for:

[0130] x i =r s cosα s cosδ s , z i =r s sinα s cosδ s ,y i =r s sinδ s (3)

[0131] Among them, r s It represents the distance from the center of the sun to the center of the earth.

[0132] The other steps and parameters are the same as those in the first to third embodiments.

[0133] Specific embodiment 5: This embodiment differs from specific embodiments 1 to 4 in that the specific process of step 2 is as follows:

[0134] Step 2.1 Calculate the angle β between the satellite and the earth e :

[0135]

[0136] Among them, R e is the radius of the Earth, R sat is the distance from the center of gravity of the satellite antenna to the center of the Earth;

[0137] The angle between the line connecting the center of gravity of the satellite antenna and the center of the earth and the line connecting the center of gravity of the satellite antenna and the center of the sun is denoted as β. Compare β with β e Size:

[0138] (2) When β is greater than βe When the satellite antenna is illuminated, the solar irradiance S c for:

[0139]

[0140] in, Indicates the average solar irradiance (equal to 1367W / m 2 );R s is the distance from the center of the sun to the center of the earth; is the average distance between the center of the sun and the center of the earth;

[0141] (2) When β is less than or equal to β e When the satellite antenna is in the shadow area, the solar irradiance S c =0;

[0142] Step 2.2. Solar radiation heat flux per unit area q s for:

[0143]

[0144] where n is the normal of the antenna surface mesh in the orbital coordinate system.

[0145] The other steps and parameters are the same as those in the first to fourth embodiments.

[0146] Specific embodiment 6: This embodiment differs from any one of specific embodiments 1 to 5 in that the specific process of step 3 is as follows:

[0147] The parabolic model of the satellite antenna is established using the commercial software Abaqus. The heat flux per unit area calculated in step 2 is then loaded using the Load module in Abaqus to obtain the deformation of each area of ​​the satellite antenna. The optimal position of the distributed control actuator is then determined based on the deformation of each area of ​​the satellite antenna.

[0148] The other steps and parameters are the same as those in the first to fifth embodiments.

[0149] The Abaqus simulation duration was set to 86,400 seconds, and the thermal-displacement coupling analysis step was selected. The deformation distribution of the satellite antenna surface was analyzed to determine the extent of thermally induced deformation at different locations on the antenna surface over a day. Based on the degree of deformation, actuators were prioritized for locations with the most significant deformation and the greatest impact on antenna performance, enabling more precise shape control.

[0150] Specific embodiment 7: This embodiment differs from any one of specific embodiments 1 to 6 in that the specific process of step 4 is as follows:

[0151] Step 4.1. The standard form of the structural dynamics equation of the spaceborne antenna is:

[0152]

[0153] Where M, C′ and K are the inertia matrix, damping matrix and stiffness matrix of the satellite antenna respectively, B′ is the actuator arrangement matrix, and F c is the cable control force vector in physical space, u d is the force generated by thermal radiation in physical space, x represents the vector composed of the coordinates of the nodes where each actuator is arranged, represents the first-order derivative of x, represents the second derivative of x;

[0154] Perform modal coordinate transformation on the structural dynamics equations of the spaceborne antenna:

[0155] x=Φη (8)

[0156] Where η is the N-dimensional modal coordinate vector of the flexible structure, Φ is the first N-order regular modal matrix;

[0157] Then the form of the distributed cable dynamic equation is:

[0158]

[0159] in, is the second-order derivative of η, is the first-order derivative of η, Λ a is the frequency matrix, ξ is the damping coefficient matrix, B c is the distributed cable control matrix, Φ T is the transpose of Φ;

[0160] Frequency matrix Λ a The damping coefficient matrix ξ is in the form of:

[0161]

[0162] Where, ω a1 、ω a2 ,…,ω aN Determined by formula (11):

[0163]

[0164] Where r = 1, 2, ..., N, Φ r is the first r-order regular mode matrix, is the rth order main frequency corresponding to the structural vibration, and:

[0165]

[0166] The damping coefficient is calculated using the Rayleigh damping method, so the damping coefficient ξ r Determined by formula (13):

[0167]

[0168] Among them, a0 and a1 are constants;

[0169] Step 4.2: Convert equation (9) into the system state equation:

[0170]

[0171] y=CX (15)

[0172] Where, is the first-order derivative of X, u represents the system control quantity, y represents the system output, and I is the unit matrix;

[0173]

[0174] Distributed cable control matrix B c The solution is related to the cable stress. The stress analysis of a single cable is as follows: Figure 3 As shown, the control forces at both ends of the cable are decomposed into three axes:

[0175] L=[T axl T ayl T azl -T axl -T ayl -T azl ] (17)

[0176] Where, (x l ,y l ,z l ) represents the position of the l-th actuator node, l = 1, 2, …, L′, L′ represents the total number of actuators arranged;

[0177] Then the distributed cable control matrix B of the lth cable is cl for:

[0178]

[0179] Among them, φ xi 、φ yi and φ zi They are the three degrees of freedom corresponding to the i-th order main mode, i=1,2,…,N,B c =[B c1 ,B c2 ,…,B cL′ ];

[0180] Step 4.3. Convert the system state equations of Equation (14) and Equation (15) into discrete state space equations:

[0181]

[0182] Among them, X(k|k) represents the measured value of the state variable at time k, X(k+1|k+1) represents the measured value of the state variable at time k+1, u(k) represents the control quantity at time k, η(k+1|k+1) represents the output quantity at time k+1, and u d (k) represents the force generated by thermal radiation at time k in physical space, and the discrete state matrices H(k), F(k) and G(k) are respectively:

[0183] F(k)=I+AT s ,G(k)=BT s ,H(k)=DT s (20)

[0184] Where, T s is the sampling time, I is the identity matrix;

[0185] The state variables in the prediction domain are:

[0186]

[0187] Where N p is the prediction time domain; X(k+N p |k) is the system at time k, the prediction time domain N p Internal system state quantity; η(k+N p |k) is the system at time k, the prediction time domain N p The internal system predicts the output; u(k+j|k) is the prediction of the control quantity at time k+j at time k; u d (k+j|k) is the prediction of the force generated by thermal radiation at time k+j at time k;

[0188] Rewrite Equation (21) into matrix form:

[0189] E(k+1)=Ψ(k)X(k|k)+Θ(k)U(k)+Ξ(k)U d (k) (23)

[0190] Γ(k+1)=L(k)E(k+1) (24)

[0191] In the formula, E(k+1), Ψ(k), Θ(k), Ξ(k), Γ(k+1), L(k), U(k) and U d The expression of (k) is as follows:

[0192]

[0193] It should be noted that the matrix C in formula (29) does not change with time;

[0194] In order to achieve the optimal control effect of the system in the prediction time domain, the objective function J of the lth subsystem (each subsystem corresponds to an actuator) is constructed. l (k) in the following form:

[0195] J l (k)=U l (k) T R l U l (k)+(Γ l (k+1)-Γ rl (k+1)) T Q l (Γ l (k+1)-Γ rl (k+1)) (31)

[0196] Among them, U l (k) represents the control quantity (predicted value) corresponding to the lth subsystem in U(k), R l represents the control constraint weight matrix corresponding to the lth subsystem; Γ l (k+1) represents the output modal coordinates of the lth subsystem in Γ(k+1); Γ rl (k+1) represents the target modal coordinate corresponding to the lth subsystem; Q l Represents the curve tracking weight matrix corresponding to the lth subsystem;

[0197] Considering the output characteristic constraints of the cable, the above optimization function can be converted into a quadratic programming problem, and the expression can be written as:

[0198]

[0199] Where, represents the optimal control quantity sequence of the lth subsystem, U max Indicates the maximum control amount;

[0200] The optimal output sequence u of the lth cable l (k) is expressed as:

[0201]

[0202] Among them, u(k)=[u1(k),u2(k),…,u L′ (k)].

[0203] The other steps and parameters are the same as those in the first to sixth embodiments.

[0204] Specific embodiment eight: This embodiment differs from any one of specific embodiments one to seven in that the force generated by thermal radiation in the physical space is:

[0205]

[0206] Among them, q s is the solar radiation heat flux per unit area, S flux is the sun-receiving area of ​​the antenna, c is the speed of light, and 3×10 8 m / s.

[0207] The other steps and parameters are the same as those in the first to seventh embodiments.

[0208] Experimental part

[0209] The algorithm of the present invention is tested in a simulation environment. Assume that the initial time and orbit parameters are as shown in Table 1:

[0210] Table 1 Initial time and orbit parameters

[0211]

[0212]

[0213] Through simulation calculation, the change results of the solar radiation heat flux received by the unit surface in the positive direction of the z-axis and the negative direction of the z-axis of the orbital coordinate system within two days are obtained, as shown in the following figure: Figure 4 Then, the antenna model was established in Abaqus, the calculated heat flux was imported, the thermal-displacement coupling analysis step was set, and the simulation time was set to 86400s. The results of the antenna thermal deformation within 1 day were obtained, as shown in the figure below. Figure 5(a) to Figure 5(d) As shown, in order to actively maintain the shape accuracy, the axisymmetric configuration of the cables is performed in the solved antenna according to the local peak value of thermal displacement, as shown in Figure 6 shown.

[0214] The controller proposed by the method of the present invention is used to control the antenna model and compared with the PD control. The control results of the first six modal coordinates are as follows: Figure 7 As shown, take Figure 6 The deformation suppression results at the three nodes 103, 150 and 166 marked in FIG are as follows: Figure 8 As shown in the results, it can be seen that the method of the present invention has more advantages.

[0215] The above examples are merely illustrative of the calculation model and process of the present invention and are not intended to limit the embodiments of the present invention. Persons skilled in the art will readily appreciate that other variations or modifications based on the above description are possible. This list of embodiments is not exhaustive; however, any obvious variations or modifications derived from the technical solution of the present invention remain within the scope of protection of the present invention.

Claims

1. A method for predicting and optimizing the thermal deformation of a spaceborne antenna, characterized in that: The method specifically comprises the following steps: Step 1: Establish the geocentric equatorial inertial coordinate system and the orbital coordinate system to obtain the position vector of the sun in the geocentric equatorial inertial coordinate system at each moment; Step 2: Calculate the heat flux per unit area of ​​solar radiation on the satellite orbit based on the sun position vector; Step 3: Use the finite element method to build a satellite antenna model. Then, perform a thermal deformation analysis on the satellite antenna based on the heat flux per unit area of ​​solar radiation. Based on the analysis results, determine the optimal placement of the distributed control actuators on the satellite antenna. Step 4: Establish an objective function according to the position of the actuator arrangement, and perform shape optimal distributed model predictive control based on the objective function.

2. The method for predicting and optimizing the thermal deformation of a space-borne antenna according to claim 1, wherein: The geocentric equatorial inertial coordinate system is specifically: The origin E of the geocentric equatorial inertial coordinate system is located at the center of the earth, z i The axis is perpendicular to the equatorial plane and z i The positive direction of the axis points to the North Pole; i axis and y i The axes are all located in the equatorial plane, where x i The positive direction of the axis points to the vernal equinox, x i Axis, y i axis and z i The axes form a right-handed coordinate system.

3. The method for predicting and optimizing the thermal deformation of a space-borne antenna according to claim 2, wherein: The origin O of the orbital coordinate system is located at the center of mass of the spacecraft, z o The positive direction of the axis points to the center of the earth, x o Axis and z o The axes are vertical and x o The x-axis is in the orbital plane of the spacecraft. o The positive direction of the y axis points to the direction of satellite flight, o The axis is perpendicular to the orbital plane of the spacecraft, x o Axis, y o axis and z o The axes form a right-handed coordinate system.

4. The method for predicting and optimizing thermal deformation of a space-borne antenna according to claim 3, wherein: The position vector of the sun in the geocentric equatorial inertial coordinate system is: The transformation matrix M from the orbital coordinate system to the geocentric equatorial inertial coordinate system is: Where Ω is the right ascension of the ascending node, i is the orbital inclination, μ is the argument of latitude, μ = ω + θ, ω is the argument of perigee, and θ is the true anomaly. a s =arctan(sinεsinλ s / cosλ s ),d s =arcsin(sinεsinλ s ) (2) where ε is the inclination of the ecliptic, λ s is the ecliptic longitude, α s is the right ascension, δ s is declination; The position component of the sun in the geocentric equatorial inertial coordinate system (x i ,y i ,z i )for: x i =r s cosα s cosδ s ,z i =r s sinα s cosδ s ,y i =r s sinδ s (3) Among them, r s It represents the distance from the center of the sun to the center of the earth.

5. The method for predicting and optimizing the thermal deformation of a space-borne antenna according to claim 4, wherein: The specific process of step 2 is: Step 2.1 Calculate the angle β between the satellite and the earth e : Among them, R e is the radius of the Earth, R sat is the distance from the center of gravity of the satellite antenna to the center of the Earth; The angle between the line connecting the center of gravity of the satellite antenna and the center of the earth and the line connecting the center of gravity of the satellite antenna and the center of the sun is denoted as β. Compare β with β e Size: (1) When β is greater than β e When the satellite antenna is illuminated, the solar irradiance S c for: in, represents the average solar irradiance; R s is the distance from the center of the sun to the center of the earth; is the average distance between the center of the sun and the center of the earth; (2) When β is less than or equal to β e When the satellite antenna is in the shadow area, the solar irradiance S c =0; Step 2.

2. Solar radiation heat flux per unit area q s for: where n is the normal of the antenna surface mesh in the orbital coordinate system.

6. The method for predicting and optimizing the thermal deformation of a space-borne antenna according to claim 5, wherein: The specific process of step three is: The parabolic model of the satellite antenna is established using the commercial software Abaqus. The heat flux per unit area calculated in step 2 is then loaded using the Load module in Abaqus to obtain the deformation of each area of ​​the satellite antenna. The optimal position of the distributed control actuator is then determined based on the deformation of each area of ​​the satellite antenna.

7. The method for predicting and optimizing the thermal deformation of a space-borne antenna according to claim 6, wherein: The specific process of step 4 is as follows: Step 4.

1. The structural dynamics equation of the satellite-borne antenna is: Where M, C′ and K are the inertia matrix, damping matrix and stiffness matrix of the satellite antenna respectively, B′ is the actuator arrangement matrix, and F c is the cable control force vector in physical space, u d is the force generated by thermal radiation in physical space, x represents the vector composed of the coordinates of the nodes where each actuator is arranged, represents the first-order derivative of x, represents the second derivative of x; Perform modal coordinate transformation on the structural dynamics equations of the spaceborne antenna: x=Φη (8) Where η is the N-dimensional modal coordinate vector of the flexible structure, Φ is the first N-order regular modal matrix; Then the form of the distributed cable dynamic equation is: in, is the second-order derivative of η, is the first-order derivative of η, Λ a is the frequency matrix, ξ is the damping coefficient matrix, B c is the distributed cable control matrix, Φ T is the transpose of Φ; Frequency matrix Λ a The damping coefficient matrix ξ is in the form of: Where, ω a1 、ω a2 ,…,ω aN Determined by formula (11): Where r = 1, 2, ..., N, Φ r is the first r-order regular mode matrix, is the rth order main frequency corresponding to the structural vibration, and: Then the damping coefficient ξ r Determined by formula (13): Among them, a0 and a1 are constants; Step 4.2: Convert equation (9) into the system state equation: y=CX (15) Where, is the first-order derivative of X, u represents the system control quantity, y represents the system output, and I is the unit matrix; Distributed cable control matrix B c The solution is related to the force on the cable, and the control force at both ends of the cable is decomposed into three axes: L=[T axl T ayl T azl -T axl -T ayl -T azl ] (17) Where, (x l ,y l ,z l ) represents the position of the l-th actuator node, l = 1, 2, …, L′, L′ represents the total number of actuators arranged; Then the distributed cable control matrix B of the lth cable is cl for: Among them, φ xi 、φ yi and φ zi They are the three degrees of freedom corresponding to the i-th order main mode, i=1,2,…,N,B c =[B c1 ,B c2 ,…,B cL′ ]; Step 4.

3. Convert the system state equations of Equation (14) and Equation (15) into discrete state space equations: Among them, X(k|k) represents the measured value of the state variable at time k, X(k+1|k+1) represents the measured value of the state variable at time k+1, u(k) represents the control quantity at time k, η(k+1|k+1) represents the output quantity at time k+1, and u d (k) represents the force generated by thermal radiation at time k in physical space, and the discrete state matrices H(k), F(k) and G(k) are respectively: F(k)=I+AT s ,G(k)=BT s ,H(k)=DT s (20) Where, T s is the sampling time, I is the identity matrix; The state variables in the prediction domain are: Where N p is the prediction time domain; X(k+N p |k) is the system at time k, the prediction time domain N p Internal system state quantity; η(k+N p |k) is the system at time k, the prediction time domain N p The internal system predicts the output; u(k+j|k) is the prediction of the control quantity at time k+j at time k; u d (k+j|k) is the prediction of the force generated by thermal radiation at time k+j at time k; Rewrite Equation (21) into matrix form: E(k+1)=Ψ(k)X(k|k)+Θ(k)U(k)+Ξ(k)U d (k) (23) Γ(k+1)=L(k)E(k+1) (24) In the formula, E(k+1), Ψ(k), Θ(k), Ξ(k), Γ(k+1), L(k), U(k) and U d The expression of (k) is as follows: Construct the objective function J of the lth subsystem l (k) in the following form: J l (k)=U l (k) T R l U l (k)+(Γ l (k+1)-C rl (k+1)) T Q l (C l (k+1)-C rl (k+1)) (31) Among them, U l (k) represents the control quantity corresponding to the lth subsystem in U(k), R l represents the control constraint weight matrix corresponding to the lth subsystem; Γ l (k+1) represents the output modal coordinates of the lth subsystem in Γ(k+1); Γ rl (k+1) represents the target modal coordinate corresponding to the lth subsystem; Q l Represents the curve tracking weight matrix corresponding to the lth subsystem; Where, represents the optimal control quantity sequence of the lth subsystem, U max Indicates the maximum control amount; The optimal output sequence u of the lth cable l (k) is expressed as: Where, u(k) = [u1(k),u2(k),...,u L′ (k)].

8. The method for predicting and optimizing the thermal deformation of a space-borne antenna according to claim 7, wherein: The force generated by thermal radiation in the physical space is: Among them, q s is the solar radiation heat flux per unit area, S flux is the sun-receiving area of ​​the antenna, and c is the speed of light.

Images