CN104570742B - Feedforward PID (proportion, integration and differentiation) control based rapid high-precision relative pointing control method of noncoplanar rendezvous orbit - Google Patents

Abstract

The invention relates to a rapid high-precision relative pointing control method of a noncoplanar rendezvous orbit, in particular to a feedforward PID (proportion, integration and differentiation) control based rapid high-precision relative pointing control method of the noncoplanar rendezvous orbit, and aims to solve the problem that rapid high-precision tracking-pointing control methods related to satellite attitudes under the noncoplanar rendezvous orbit do not exist in the prior art. The feedforward PID control based rapid high-precision relative pointing control method including adopting an Euler angle for attitude description of a spacecraft, creating dynamical and kinematical equations of the spacecraft, and acquiring an precise expected angle theta according to a noise-containing expected angle z by a Kalman filtering algorithm of a spaceborne computer; designing an attitude control law of each axis according to a formula of ; selecting two parallelly-placed single-frame control moment gyroscopes to control a yaw axis, and selecting two flywheels to control a rolling axis and a pitch axis respectively; calculating gyroscopic moment T and flywheel actual output moment uw to complete rapid high-precision relative pointing control of the noncoplanar rendezvous orbit. The feedforward PID control based rapid high-precision relative pointing control method is applicable to rapid high-precision relative pointing control of the noncoplanar rendezvous orbit.

Description

Fast and high-precision relative pointing of fast-changing tracks crossing different planes based on feed-forward PID control Control Method

technical field

The invention relates to a fast and high-precision relative pointing control method for intersecting fast-changing tracks with different planes.

Background technique

With the development of aerospace science and technology, space technology has been widely used in various fields. It can be said that space has existed as a resource that human beings rely on for survival and development. In order to make the satellite payload point to a specific target with a certain accuracy, it is first necessary to calculate the desired attitude when the satellite load line-of-sight axis points to the target according to the target position, and then track the desired attitude so that the satellite is in line with the commanded attitude under the influence of orbital motion and disturbance torque The deviation and relative rotation speed are kept within the allowable range.

"Robust Backstepping Adaptive Attitude Maneuvering and Active Vibration Suppression of Flexible Spacecraft" studies the problem of attitude maneuvering and active vibration suppression of flexible spacecraft, and designs a dual-loop robust control method, which is effective for parameter uncertainties Strong robustness, increasing the damping of the flexible structure, can realize the high-precision attitude control and vibration suppression of the flexible spacecraft. However, the algorithm enters the steady state slowly, which affects the rapidity of the system.

"Application of Variable Structure Control in Rapid Attitude Maneuvering Control of Flexible Spacecraft" focuses on the variable structure control law of flexible spacecraft. The parameter setting method of the fastest maneuvering down-variable structure is analyzed. This algorithm has strong robustness, simple algorithm, does not depend on model parameters, has better dynamic performance than PID control, and improves the rapidity and accuracy of satellite attitude maneuvering. . However, the algorithm relies too much on the accuracy of satellite attitude dynamics modeling, which limits its application in practical engineering.

"Research on Attitude Control Method for Satellites with Solar Panels" studies the problems of complex attitude dynamics of flexible satellites, modeling uncertainty and limited torque output of actual satellite actuators, and designs an inaccurate model based on support vector machines. Inverse Model Control Law. The control law has good steady-state effect and dynamic quality, effectively reduces the influence of flexible mode vibration on attitude control, and has a certain ability to suppress disturbance. However, the stability analysis and calculation process of this algorithm are relatively complicated, and the application is still relatively few, and there is no literature on the combination of support vector machine and flexible satellite attitude.

At present, there is almost no research on the problem of fast attitude tracking and high-precision tracking of satellites in different plane crossing orbits.

Contents of the invention

The present invention aims to solve the problem that there is no control method for the attitude fast and high-precision tracking and pointing of satellites under different plane crossing orbits in the existing control methods.

A fast and high-precision relative pointing control method for fast-changing tracks across different planes based on feed-forward PID control, including the following steps:

Step 1: Determine the desired pose:

The tracking star and the target star are located on different orbits, and the tracking star needs to detect the position of the target autonomously; in order to make the laser emitter or observation equipment of the tracking star always point to the target star, the desired attitude must be determined first, so that the problem of pointing to the target can be easily solved. It becomes a problem of attitude tracking; due to the limited detection range, laser attack or observation and monitoring can only be carried out when the distance between the tracking star and the target star is relatively close. , at this time, the expected attitude often changes quickly (the embodiment of "fast change"), so it is necessary to have a faster response speed for the pointing tracking of the target; and because the pointing task of the target star often requires a certain accuracy, the The accuracy puts forward higher requirements;

Assuming that the line of sight axis of the tracking star coincides with the x-axis of the body, the x-axis of the desired attitude coordinate system points to the target star, the y-axis is perpendicular to the plane formed by the x-axis and the tracking star’s geocentric vector, and is opposite to the orbital angular velocity, and the z-axis and x , y-axis to form a right-handed coordinate system, as shown in Figure 1; the desired attitude coordinate system designed in this way is almost obtained by rotating the center-of-mass orbital coordinate system around the z-axis by a certain angle, and the orbital coordinate system of the mass center is used as the attitude reference coordinate system, then the expected attitude The change is approximately a rotation around the z-axis;

In the earth-centered inertial coordinate system, the unit vector of each coordinate axis of the desired attitude coordinate system is expressed as:

Where r _t and r _c are the position vectors of the tracking star and the target star in the earth-centered inertial coordinate system respectively; i, j, k are the unit vectors of each coordinate axis of the desired attitude coordinate system;

Let x ₁ , y ₁ , z ₁ be the unit vectors of each coordinate axis of the earth-centered inertial coordinate system, x ₂ , y ₂ , z ₂ be the unit vectors of each coordinate axis of the center-of-mass orbital coordinate system, and each coordinate axis of the desired attitude coordinate system The unit vectors i, j, k of the earth-centered inertial coordinate system and the unit vectors x ₁ , y ₁ , z ₁ of each coordinate axis of the earth-centered inertial coordinate system are obtained to obtain the cosine transformation matrix R _di of the desired attitude coordinate system relative to the earth-centered inertial coordinate system, and the mass center orbit The unit vector x ₂ , y ₂ , z ₂ of each coordinate axis of the coordinate system and the unit vector x ₁ , y ₁ , z ₁ of each coordinate axis of the earth-centered inertial coordinate system are used to obtain the cosine of the center-of-mass orbital coordinate system relative to the earth-centered inertial coordinate system Transformation matrix R _oi ,

Then the transformation matrix of the desired attitude coordinate system relative to the centroid orbital coordinate system is

Step 2: The Kalman filter estimation uses the second order derivative of the desired attitude angle in the control law, which cannot be directly measured; the measured desired attitude angle is affected by noise, and the direct second difference will amplify the noise, so the design of the standard continuous type Kalman The Mann filter determines the expected attitude angle and its derivative; the expected attitude angle is the attitude angle obtained under the coordinate axis where the i, j, and k unit vectors are located as the expected attitude coordinate system, and the expected attitude angle is θ, is the first derivative of the desired attitude angle θ, is the second derivative of the desired attitude angle θ; take the state vector as

The expected angle z containing noise information is obtained by measurement:

z=H _z x+v (3)

In the formula, H _z =[1,0,0], v is the unestimated high-frequency noise;

Input the desired angle z with noise into the Kalman filter; in the Kalman filter,

Let the estimators for the three axes take the same form; the state equation for each axis is:

in, w is a virtual white noise with non-zero variance to reflect the change characteristics of the actual angular acceleration information;

The estimator is

K=PH _z ^T R ^-1 (6)

PA ^T +AP-PH _z ^T R ^-1 _Hz P+BQB ^T ＝0 (7)

is the state estimate for the state variable x; The first derivative of the state estimation value of the state variable x; Q is the semi-positive definite system noise variance matrix; R is the positive definite measurement noise matrix; P is the state The variance matrix of

According to the expected angle z containing noise, the precise expected angle θ is obtained through the Kalman filter algorithm of the on-board computer;

Step 3: Design the attitude control law of each axis. The specific expression of the attitude control law is the control torque u; it has the following form:

Where u is the control torque, I is the principal inertia of the axis, K _p is the proportional amplification factor, K _i is the integral control coefficient, K _d is the differential control coefficient; a _d is the attitude angle in the center of mass coordinate system of the body, is the second derivative of a _d ; e is the Euler attitude angle error;

Step 4: Select the actuator, use two parallel single-frame control moment gyroscopes to control the yaw axis, that is, the z-axis; use two flywheels to control the roll axis and pitch axis, namely the x-axis and y-axis.

The present invention has the following beneficial effects:

1. In actual engineering, the thruster has an upper limit of amplitude, which is not continuous. Compared with some existing orbit transfer schemes under continuous control, the present invention considers the amplitude limit of the thruster, which is more in line with actual engineering application.

2. Compared with some schemes that only consider the optimal time or fuel consumption, the present invention considers both the transfer time and the fuel consumption, and can find the time-fuel maximum by adjusting the proportion of the two. Excellent control scheme.

3. Different from some studies on orbit transfer in inertial system, the present invention studies orbit transfer in relative orbit coordinate system, and the orbit transfer start point and end point in relative orbit coordinate system are not fixed points in inertial system, making this The invention has greater practical application value in the problem of relative orbit transfer between spacecraft.

Description of drawings

Figure 1 is a schematic diagram of the desired attitude coordinate system;

Fig. 2 is a schematic diagram of body and track coordinate system;

Figure 3 is a schematic diagram of the installation configuration of the control moment gyroscope;

Figure 4 is a schematic diagram of the PD controller of the flywheel in torque mode;

Figure 5 is a curve diagram of desired posture and body posture variation;

Fig. 6 attitude angle error change curve;

Figure 7 satellite angular velocity change curve;

Figure 8 CMG and flywheel output torque curve;

Fig. 9 is a curve diagram of total control torque change.

detailed description

Specific implementation mode one:

A fast and high-precision relative pointing control method for fast-changing tracks across different planes based on feed-forward PID control, including the following steps:

Step 1: Determine the desired pose:

The tracking star and the target star are located on different orbits, and the tracking star needs to detect the position of the target autonomously; in order to make the laser emitter or observation equipment of the tracking star always point to the target star, the desired attitude must be determined first, so that the problem of pointing to the target can be easily solved. It becomes a problem of attitude tracking; due to the limited detection range, laser attack or observation and monitoring can only be carried out when the distance between the tracking star and the target star is relatively close. , at this time, the expected attitude often changes quickly (the embodiment of "fast change"), so it is necessary to have a faster response speed for the pointing tracking of the target; and because the pointing task of the target star often requires a certain accuracy, the The accuracy puts forward higher requirements;

Assuming that the line of sight axis of the tracking star coincides with the x-axis of the body, the x-axis of the desired attitude coordinate system points to the target star, the y-axis is perpendicular to the plane formed by the x-axis and the tracking star’s geocentric vector, and is opposite to the orbital angular velocity, and the z-axis and x , y-axis to form a right-handed coordinate system, as shown in Figure 1; the desired attitude coordinate system designed in this way is almost obtained by rotating the center-of-mass orbital coordinate system around the z-axis by a certain angle, and the orbital coordinate system of the mass center is used as the attitude reference coordinate system, then the expected attitude The change is approximately a rotation around the z-axis;

In the earth-centered inertial coordinate system, the unit vector of each coordinate axis of the desired attitude coordinate system is expressed as:

Where r _t and r _c are the position vectors of the tracking star and the target star in the earth-centered inertial coordinate system respectively; i, j, k are the unit vectors of each coordinate axis of the desired attitude coordinate system;

Let x ₁ , y ₁ , z ₁ be the unit vectors of each coordinate axis of the earth-centered inertial coordinate system, x ₂ , y ₂ , z ₂ be the unit vectors of each coordinate axis of the center-of-mass orbital coordinate system, and each coordinate axis of the desired attitude coordinate system The unit vectors i, j, k of the earth-centered inertial coordinate system and the unit vectors x ₁ , y ₁ , z ₁ of each coordinate axis of the earth-centered inertial coordinate system are obtained to obtain the cosine transformation matrix R _di of the desired attitude coordinate system relative to the earth-centered inertial coordinate system, and the mass center orbit The unit vector x ₂ , y ₂ , z ₂ of each coordinate axis of the coordinate system and the unit vector x ₁ , y ₁ , z ₁ of each coordinate axis of the earth-centered inertial coordinate system are used to obtain the cosine of the center-of-mass orbital coordinate system relative to the earth-centered inertial coordinate system Transformation matrix R _oi ,

Then the transformation matrix of the desired attitude coordinate system relative to the centroid orbital coordinate system is

Step 2: The Kalman filter estimation uses the second-order derivative of the desired attitude angle in the control law, which cannot be directly measured; the measured desired attitude angle is affected by noise, and if the direct secondary difference will amplify the noise, so the design of the standard continuous Kalman The Mann filter determines the expected attitude angle and its derivative; the expected attitude angle is the attitude angle obtained under the coordinate axis where the i, j, and k unit vectors are located as the expected attitude coordinate system, and the expected attitude angle is θ, is the first derivative of the desired attitude angle θ, is the second derivative of the desired attitude angle θ; take the state vector as

The expected angle z containing noise information is obtained by measurement:

z=H _z x+v (3)

In the formula, H _z =[1,0,0], v is the unestimated high-frequency noise;

Input the desired angle z with noise into the Kalman filter; in the Kalman filter,

Let the estimators for the three axes take the same form; the state equation for each axis is:

in, w is a virtual white noise with non-zero variance to reflect the change characteristics of the actual angular acceleration information;

The estimator is

K=PH _z ^T R ^-1 (6)

PA ^T +AP-PH _z ^T R ^-1 _Hz P+BQB ^T ＝0 (7)

is the state estimate for the state variable x; The first derivative of the state estimation value of the state variable x; Q is the semi-positive definite system noise variance matrix; R is the positive definite measurement noise matrix; P is the state The variance matrix of

According to the expected angle z containing noise, the precise expected angle θ is obtained through the Kalman filter algorithm of the on-board computer;

Step 3: Design the attitude control law of each axis. The specific expression of the attitude control law is the control torque u; it has the following form:

Where u is the control torque, I is the principal inertia of the axis, K _p is the proportional amplification factor, K _i is the integral control coefficient, K _d is the differential control coefficient; a _d is the attitude angle in the center of mass coordinate system of the body, is the second derivative of a _d ; e is the Euler attitude angle error;

Step 4: Select the actuator, use two parallel single-frame control moment gyroscopes to control the yaw axis, that is, the z-axis; use two flywheels to control the roll axis and pitch axis, namely the x-axis and y-axis.

Specific implementation mode two: the implementation process of step 3 of this embodiment mode is:

The angular velocity ω of the tracking star in inertial space is equal to the sum of the angular velocity ω _br of the tracking star body coordinate system relative to the attitude reference coordinate system and the implicated angular velocity ω _ri of the attitude reference coordinate system relative to the center-of-mass inertial coordinate system, namely

ω＝ω _br +ω _ri (8)

In inertial orientation flight mode

At this time, the reference coordinate system is the center-of-mass inertial coordinate system, ω _ri =0, we have

C _x ( ), C _y ( ), and C _z ( ) are the rotation matrices in the coordinate system conversion process;

ψ is the yaw angle; is the pitch angle; is the roll angle; ω _x , ω _y , and ω _z are the components of the tracking star angular velocity ω in the x-axis, y-axis, and z-axis in the center-of-mass inertial coordinate system;

As shown in Figure 2, the Euler angle ψ, The geometric meaning of is:

ψ is the yaw angle—the angle between the projection of the tracking star roll axis O _b X _b (pointing to the tracking star velocity direction) on the local horizontal plane and the orbital coordinate system O _o X _o axis;

is the pitch angle—the angle between the tracking star rolling axis O _b X _b and its projection on the local horizontal plane;

Roll angle—the angle between the pitch axis O _b Y _b of the tracking star and its projection on the local horizontal plane;

The three angles known in the art

Get the Euler angle attitude kinematic equation:

Since the attitude control task is similar to single-axis rotation, the attitude tracking control law is designed with three-axis decoupling; when the tracking star approximately rotates around the z-axis, the x-axis and y-axis attitude angles are approximately 0, and the orbital angular velocity is compared with the tracking star angular velocity is also a small amount, it can be known from formula (10) The angular velocity measured by the gyro can be directly used for feedback control, and the form of 3-1-2 sequence Euler angles is selected to describe the attitude of the tracking star;

Convert the expected attitude angle θ to the attitude angle a _d in the body center of mass coordinate system according to the transformation matrix a _d = θR _do ; The attitude angle of the body is recorded as a _b , and the Euler attitude angle error e=a _d -a _b ; when the expected attitude changes rapidly around the z-axis, the traditional PID control law cannot meet the tracking accuracy requirements, so the second order of the expected Euler angle is added The derivative (obtained by filtering and estimation in step 2) is a feed-forward term, which is used to predict changes in angular acceleration and improve pointing and tracking capabilities;

The attitude control law for each axis has the following form:

Among them, u is the control torque, I is the main inertia of the shaft, K _p is the proportional amplification factor, K _i is the integral control coefficient, and K _d is the differential control coefficient;

In order to prevent integral saturation, the idea of "integral separation" is adopted, that is, the integral action is only introduced when the error is small, and when |e| is less than a certain threshold value, the integration starts, and when |e| The threshold is determined according to the actual Euler attitude angle error requirements.

Other steps and parameters are the same as in the first embodiment.

Specific implementation mode three: the implementation process of step 4 of this embodiment mode is:

Due to the large change in the expected attitude near the intersection point of the orbits of the two stars, the requirements for the control torque are relatively high; although the thruster can provide a large torque, the output torque is not continuous, and the command torque cannot be accurately tracked, which destroys the control System performance; In the feedforward PID attitude tracking control, the accuracy of the control quantity is required to be higher, and it is more likely to cause chattering when the thruster is used as the actuator, and the ideal control effect cannot be achieved;

In the present invention, only the torque requirement of the yaw axis is large, and the torque requirement of the other two axes is very small. The flywheel and the CMG can be used in combination, that is, the yaw axis is controlled by two parallel single-frame control moment gyros, that is, the z-axis; The two flywheels control the roll axis and the pitch axis respectively, that is, the x-axis and the y-axis;

Control Momentum Gyroscope, referred to as CMG; the installation configuration of Control Momentum Gyroscopes (CMGs) in the tracking star body coordinate system is shown in Figure 3, h ₁ and h ₂ are the angular momentum of the two gyroscope rotors; in In the installation configuration, the relative coordinate system of the control moment gyro installation is established, and the three axes point to Fx, Fy, and Fz respectively; this configuration can only output the torque in the Fx and Fz directions, and the two-dimensional vectors below represent the directions of Fx and Fz;

In order to reduce the influence of the coupling moment, try to make the angular momentum of the entire tracking star zero, that is, when the angular velocity of the tracking star body is zero, make h ₁ and h ₂ in opposite directions, and define the direction frame angle shown in the figure as zero;

The yaw axis is controlled by two parallel single-frame control moment gyros, the two gyro rotor angular momentums h ₁ and h ₂ are both h, and the frame angles of the CMG are δ ₁ and δ ₂ respectively, according to the yaw axis control Calculate the gyro frame angles δ ₁ and δ ₂ inversely from the control torque u designed by the law, and then calculate the gyro torque T according to δ ₁ and δ ₂

The total angular momentum of the gyroscope group in the tracking star system is

The synthetic gyro torque T generated by the rotation of each gyro frame can be expressed as

is the first derivative of H;

In addition to the yaw axis, the roll axis and the pitch axis use the flywheel as the actuator;

The roll axis and the pitch axis are respectively controlled by a flywheel. In the working mode of the control torque u designed by the control law of the roll axis and the pitch axis, the flywheel is used as a simple first-order inertial link. The analysis and controller design of the flywheel are as follows (14)

In the formula, t _s is the time constant of the first-order inertial system;

The PD controller is used for design, and its form is shown in Figure 4;

Flywheel actual output torque u _w

t _s is the time constant of the first-order inertial system, k _p is the integral time constant, k _d is the differential time constant;

According to the synthesized gyro torque T generated by the rotation of the gyro frame and the actual output torque u _w of the flywheel, the fast and high-precision relative pointing control of crossing fast-changing tracks in different planes is completed.

Other steps and parameters are the same as those in Embodiment 1 or Embodiment 2.

Example

Assuming that the target star is in a 600km sun-synchronous orbit, the initial orbit elements are: the semi-major axis of the target star orbit a ₁ =6978140m, the target star orbit eccentricity e ₁ =0, the target star orbit inclination angle i ₁ =97.7597°, the target star orbit elevation The longitude of intersection Ω ₁ =0, the angular distance of the pericentric point of the target star ω ₁ =0, the true periapsis f ₁ of the target star =-6.2°; The numbers are: tracking star orbit semi-major axis a ₂ =6978140m, tracking star orbit eccentricity e ₂ =0, tracking star orbit inclination i ₂ =i ₁ +5°, tracking star orbit ascending node longitude Ω ₂ =0, tracking star orbit Angular distance of periapsis ω ₂ =0, true periapsis of tracking star f ₂ =f ₁ -0.0821°.

tracking star inertia matrix

The spatial disturbance moments added in the simulation include: sunlight pressure moment, gravity gradient moment, geomagnetic moment, and aerodynamic moment.

Controller Design Parameters

PID parameter K _p ＝1.8I, K _i ＝1, K _d ＝1.8I, I is the main inertia of the corresponding axis, the integral threshold is 0.01°, the filter parameter R＝(10/3600/57.3)^2, Q＝(500 /5000/100)^2; set the angular momentum of a single gyro rotor h=30Nms, the upper limit of the angular velocity of the gyro frame is 30°/s, and the steering law weight coefficient α=e ^-10D .

simulation analysis

The simulation images obtained by the feedforward PID attitude pointing tracking control algorithm are shown in Figure 5-Figure 9:

It can be seen from Figure 5 that the desired attitude yaw angle changes rapidly from -80° to 80°, and the roll and pitch angles are almost zero, which is in line with the idea of approximate single-axis rotation, and the body attitude can track the change of expected attitude well. It can be seen from Figure 6 that except for the initial large deviation adjustment process, the tracking accuracy is relatively good, and the error will increase near the track intersection. From Figure 8, it can be seen that the control torque is provided by the control moment gyro CMG and the flywheel during the maneuvering process, and the yaw axis torque requirement near the track intersection is relatively large, which is provided by the control moment gyro CMG.

Claims (3)

1. The fast and high-precision relative pointing control method based on feed-forward PID control of different plane crossing fast-changing tracks, characterized in that: comprising the following steps: Step 1: Determine the desired pose: Assuming that the line of sight axis of the tracking star coincides with the x-axis of the body, the x-axis of the desired attitude coordinate system points to the target star, the y-axis is perpendicular to the plane formed by the x-axis and the tracking star’s geocentric vector, and is opposite to the orbital angular velocity, and the z-axis and x , y-axis form a right-handed coordinate system, and the center-of-mass orbital coordinate system is used as the attitude reference coordinate system, then the change of the expected attitude is approximately a rotation around the z-axis; In the earth-centered inertial coordinate system, the unit vector of each coordinate axis of the desired attitude coordinate system is expressed as: $\begin{array}{c}\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}}}{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}}<span\; class="notranslate">\; |</span>}\\ \mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; \&times;</span>{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}}}{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; \&times;</span>{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}}<span\; class="notranslate">\; |</span>}\\ \mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; \&times;</span>\mathrm{<span\; class="notranslate">\; j</span>}\end{array}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$ Where r _t and r _c are the position vectors of the tracking star and the target star in the earth-centered inertial coordinate system respectively; i, j, k are the unit vectors of each coordinate axis of the desired attitude coordinate system; Let x ₁ , y ₁ , z ₁ be the unit vectors of each coordinate axis of the earth-centered inertial coordinate system, x ₂ , y ₂ , z ₂ be the unit vectors of each coordinate axis of the center-of-mass orbital coordinate system, and each coordinate axis of the desired attitude coordinate system The unit vectors i, j, k of the earth-centered inertial coordinate system and the unit vectors x ₁ , y ₁ , z ₁ of each coordinate axis of the earth-centered inertial coordinate system are obtained to obtain the cosine transformation matrix R _di of the desired attitude coordinate system relative to the earth-centered inertial coordinate system, and the mass center orbit The unit vector x ₂ , y ₂ , z ₂ of each coordinate axis of the coordinate system and the unit vector x ₁ , y ₁ , z ₁ of each coordinate axis of the earth-centered inertial coordinate system are used to obtain the cosine of the center-of-mass orbital coordinate system relative to the earth-centered inertial coordinate system Transformation matrix R _oi , Then the transformation matrix of the desired attitude coordinate system relative to the centroid orbital coordinate system is ${\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; o</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; i</span>}}{\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; i</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$ Step 2: Design a standard continuous Kalman filter to determine the expected attitude angle and its derivative; the expected attitude angle is the attitude angle obtained under the coordinate axis where the i, j, and k unit vectors are located as the expected attitude coordinate system, and the expected attitude angle is θ, is the first derivative of the desired attitude angle θ, is the second derivative of the desired attitude angle θ; take the state vector as The expected angle z containing noise information is obtained by measurement: z=H _z x+v (3) In the formula, H _z =[1,0,0], v is the unestimated high-frequency noise; Input the desired angle z with noise into the Kalman filter; in the Kalman filter, Let the estimators for the three axes take the same form; the state equation for each axis is: $\stackrel{<span\; class="notranslate">\; \&Center\; Dot;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; A</span>}\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; B</span>}\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 4</span>}<span\; class="notranslate">\; )</span>$ in, w is a virtual white noise with non-zero variance to reflect the change characteristics of the actual angular acceleration information; The estimator is $\stackrel{<span\; class="notranslate">\; ^</span>}{\stackrel{<span\; class="notranslate">\; \&Center\; Dot;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; A</span>}\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; K</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; z</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; h</span>}\mathrm{<span\; class="notranslate">\; z</span>}\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 5</span>}<span\; class="notranslate">\; )</span>$ K=PH _z ^T R ^-1 (6) PA ^T +AP-PH _z ^T R ^-1 _Hz P+BQB ^T ＝0 (7) is the state estimate for the state variable x; The first derivative of the state estimation value of the state variable x; Q is the semi-positive definite system noise variance matrix; R is the positive definite measurement noise matrix; P is the state The variance matrix of According to the expected angle z containing noise, the precise expected angle θ is obtained through the Kalman filter algorithm of the on-board computer; Step 3: Design the attitude control law of each axis. The specific expression of the attitude control law is the control torque u; it has the following form: $\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; =</span>{\stackrel{<span\; class="notranslate">\; \&CenterDot;\&CenterDot;</span>}{\mathrm{<span\; class="notranslate">\; a</span>}}}_{\mathrm{<span\; class="notranslate">\; d</span>}}\mathrm{<span\; class="notranslate">\; I</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}\mathrm{<span\; class="notranslate">\; e</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; \&Integral;</span>\mathrm{<span\; class="notranslate">\; e</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; d</span>}}\stackrel{<span\; class="notranslate">\; \&Center\; Dot;</span>}{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 11</span>}<span\; class="notranslate">\; )</span>$ Where u is the control torque, I is the principal inertia of the axis, K _p is the proportional amplification factor, K _i is the integral control coefficient, K _d is the differential control coefficient; a _d is the attitude angle in the center of mass coordinate system of the body, is the second derivative of a _d ; e is the Euler attitude angle error; Step 4: Select the actuator, use two parallel single-frame control moment gyroscopes to control the yaw axis, that is, the z-axis; use two flywheels to control the roll axis and pitch axis, namely the x-axis and y-axis. 2. the fast and high-precision relative pointing control method based on feed-forward PID control of different planes crossing fast-changing tracks according to claim 1, is characterized in that, the realization process of step 3 is: The angular velocity ω of the tracking star in inertial space is equal to the sum of the angular velocity ω _br of the tracking star body coordinate system relative to the attitude reference coordinate system and the implicated angular velocity ω _ri of the attitude reference coordinate system relative to the center-of-mass inertial coordinate system, namely ω＝ω _br +ω _ri (8) In inertial orientation flight mode At this time, the reference coordinate system is the center-of-mass inertial coordinate system, ω _ri =0, we have C _x ( ), C _y ( ), and C _z ( ) are the rotation matrices in the coordinate system conversion process; ψ is the yaw angle; is the pitch angle; is the roll angle; ω _x , ω _y , and ω _z are the components of the tracking star angular velocity ω in the x-axis, y-axis, and z-axis in the center-of-mass inertial coordinate system; Get the Euler angle attitude kinematic equation: Since the attitude control task is similar to single-axis rotation, the attitude tracking control law is designed with three-axis decoupling; when the tracking star approximately rotates around the z-axis, the x-axis and y-axis attitude angles are approximately 0, and the orbital angular velocity is compared with the tracking star angular velocity is also a small amount, it can be known from formula (10) The angular velocity measured by the gyro can be directly used for feedback control, and the form of 3-1-2 sequence Euler angles is selected to describe the attitude of the tracking star; Convert the expected attitude angle θ to the attitude angle a _d in the body center of mass coordinate system according to the transformation matrix a _d = θR _do ; Body attitude angle is denoted as a _b , Euler attitude angle error e=a _d -a _b ; The attitude control law for each axis has the following form: $\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; =</span>{\stackrel{<span\; class="notranslate">\; \&CenterDot;\&CenterDot;</span>}{\mathrm{<span\; class="notranslate">\; a</span>}}}_{\mathrm{<span\; class="notranslate">\; d</span>}}\mathrm{<span\; class="notranslate">\; I</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}\mathrm{<span\; class="notranslate">\; e</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; \&Integral;</span>\mathrm{<span\; class="notranslate">\; e</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; d</span>}}\stackrel{<span\; class="notranslate">\; \&CenterDot;</span>}{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 11</span>}<span\; class="notranslate">\; )</span>$ Among them, u is the control torque, I is the main inertia of the shaft, K _p is the proportional amplification factor, K _i is the integral control coefficient, and K _d is the differential control coefficient; When |e| is less than a certain threshold, the integration starts, and when |e| is greater than this value, there is no integration. 3. according to claim 1 or 2 described based on feed-forward PID control different plane cross fast changing track fast high-precision relative pointing control method, it is characterized in that, the realization process of step 4 is: Use two parallel single-frame control moment gyroscopes to control the yaw axis, that is, the z axis; use two flywheels to control the roll axis and pitch axis, namely the x axis and y axis; Control Moment Gyroscope, referred to as CMG; The yaw axis is controlled by two parallel single-frame control moment gyros, assuming that the angular momentums h ₁ and h ₂ of the two gyro rotors are both h, and the frame angles of the CMG are δ ₁ and δ ₂ respectively, according to the yaw axis control Calculate the gyro frame angles δ ₁ and δ ₂ inversely from the control torque u designed by the law, and then calculate the gyro torque T according to δ ₁ and δ ₂ The total angular momentum of the gyroscope group in the tracking star system is $\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; h</span>}\left[\begin{array}{c}\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span>\\ <span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; the\; s</span>}\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; the\; s</span>}\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span>\end{array}\right]<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 12</span>}<span\; class="notranslate">\; )</span>$ The synthetic gyro torque T generated by the rotation of each gyro frame can be expressed as $\mathrm{<span\; class="notranslate">\; T</span>}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>\stackrel{<span\; class="notranslate">\; \&CenterDot;</span>}{\mathrm{<span\; class="notranslate">\; h</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; h</span>}\left[\begin{array}{cc}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; sin\&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}& {\mathrm{<span\; class="notranslate">\; sin\&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}\\ <span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; cos\&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}& {\mathrm{<span\; class="notranslate">\; cos\&delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}\end{array}\right]\left[\begin{array}{c}{\stackrel{<span\; class="notranslate">\; \&CenterDot;</span>}{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}\\ {\stackrel{<span\; class="notranslate">\; \&CenterDot;</span>}{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}\end{array}\right]<span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; h</span>}\mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&delta;</span>}<span\; class="notranslate">\; )</span>\stackrel{<span\; class="notranslate">\; \&CenterDot;</span>}{\mathrm{<span\; class="notranslate">\; \&delta;</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 13</span>}<span\; class="notranslate">\; )</span>$ is the first derivative of H; In addition to the yaw axis, the roll axis and the pitch axis use the flywheel as the actuator; The roll axis and the pitch axis are respectively controlled by a flywheel. In the working mode of the control torque u designed by the control law of the roll axis and the pitch axis, the flywheel is used as a simple first-order inertial link. The analysis and controller design of the flywheel are as follows (14) $\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 14</span>}<span\; class="notranslate">\; )</span>$ In the formula, t _s is the time constant of the first-order inertial system; Design with PD controller; Flywheel actual output torque u _w ${\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; w</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; d</span>}}\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; )</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}\mathrm{<span\; class="notranslate">\; e</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 15</span>}<span\; class="notranslate">\; )</span>$ t _s is the time constant of the first-order inertial system, k _p is the integral time constant, k _d is the differential time constant; According to the synthesized gyro torque T generated by the rotation of the gyro frame and the actual output torque u _w of the flywheel, the fast and high-precision relative pointing control of crossing fast-changing tracks in different planes is completed.