CN108536014B - A Model Predictive Control Method for Spacecraft Attitude Avoidance Considering Flywheel Dynamics - Google Patents

Abstract

A model predictive control method for spacecraft attitude avoidance considering the dynamic characteristics of the flywheel, comprising the following steps: establishing a prediction model based on the attitude dynamics of the spacecraft and the dynamics of the flywheel; secondly, establishing a constrained control method according to the dynamic characteristics of the flywheel and the sight angle of the instrument Mathematical model; then design performance index functions for different task requirements, and convert the control problem into the problem of finding the extreme value of the objective function under the constraints of equality and inequality; finally, through the optimization method based on real-time iteration, quickly solve the above problems , this method can well deal with the attitude avoidance problem of the spacecraft under the constraints of the actuator, achieve the comprehensive optimization of energy and time through the design of the objective function, and through the processing of real-time iteration and hot start method, and can reduce Small amount of computation to solve optimization problems.

Description

A Model Predictive Control for Spacecraft Attitude Avoidance Considering Flywheel Dynamics method

technical field

The invention belongs to the technical field of spacecraft control, is mainly applied to the avoidance control of the attitude maneuver of the spacecraft controlled by the reaction flywheel, and particularly relates to a model prediction control method for the attitude avoidance of the spacecraft considering the dynamic characteristics of the flywheel.

Background technique

In recent years, with the rapid development of space technology, the mission requirements of on-orbit spacecraft are also increasing. Therefore, various optical instruments, such as CCD cameras and infrared interferometers, are equipped on the spacecraft. During the work of these instruments It is necessary to avoid direct facing strong light to protect the components that are sensitive to light and temperature in the instrument. Therefore, it is required that during the attitude maneuver of the spacecraft, the pointing of these instruments should avoid the direction of the strong light. At the same time, the actuator in most spacecraft attitude control systems is a combination of reaction flywheels. The reaction flywheel has the advantages of high control accuracy and only consumption of electric energy for the output torque. However, the reaction flywheel also has problems such as relatively small reaction torque and easy saturation during use, which affects the attitude control of the spacecraft. Therefore, it is very important to study the attitude avoidance problem of spacecraft under the constraints of actuator performance.

Regarding the method of attitude avoidance, there is the potential function method mentioned in the patent 201710521561.X. The potential function method can effectively avoid the constraint area, but using the model predictive control method, not only the attitude avoidance constraints can be considered, but also the actuator's In addition, the optimization of energy and time can be realized according to the performance function, which is more suitable for attitude maneuver control under the performance constraints of the actuator.

SUMMARY OF THE INVENTION

The purpose of the present invention is to overcome the deficiencies of the prior art, provide a model predictive control method for spacecraft attitude avoidance considering the dynamic characteristics of the flywheel, and the problems of actuator performance constraints and attitude constraints existing in the attitude maneuvering process of the on-orbit spacecraft , the present invention provides a model prediction control method for spacecraft attitude avoidance considering the dynamic characteristics of the flywheel. This method is a model predictive control method that can deal with the performance constraints of actuators and attitude avoidance at the same time. It uses real-time iteration and hot start methods in the optimization process to improve the efficiency of solving optimization problems in MPC, and design performance The function achieves a comprehensive optimization of energy consumption and time during attitude maneuvering.

The invention provides a model predictive control method for spacecraft attitude avoidance considering the dynamic characteristics of the flywheel, comprising the following steps:

(1) According to the attitude kinematics and dynamic characteristics of the spacecraft and the rotational dynamic characteristics of the flywheel, the attitude model of the spacecraft including the actuator is established as the prediction model of the MPC;

(2) Establish a mathematical model of attitude avoidance constraints according to the sight angle and avoidance vector of the spacecraft onboard instruments, and establish a mathematical model of actuator constraints according to the angular momentum saturation and torque saturation characteristics of the flywheel;

(3) Design the corresponding optimization objective function according to the task requirements, including the quadratic type of the actuator input and the quadratic type of the system state error, so as to comprehensively consider the execution time and energy consumption;

(4) Convert the control problem into the problem of finding the extreme value of the objective function under the constraints of the equations and states of the system dynamics equations and the input-limited inequality constraints, and then use the real-time iterative optimization processing method to quickly solve the optimization to obtain The solution is output as the control quantity of the system.

The spacecraft attitude dynamics model described in step (1) is obtained by combining the dynamic equation of the flywheel into the spacecraft attitude dynamics equation and discretizing it, and its representation is as follows:

表示ω对时间的导数；J是航天器总的惯量矩阵，简化表示为对角阵J＝diag(J₁,J₂,J₃)，J₁,J₂,J₃为绕惯量主轴的转动惯量；S(ω)是斜对称矩阵，其形式为τ表示执行机构的输出力矩；q＝[q₀,q₁,q₂,q₃]^T表示航天器的姿态单位四元数，

表示航天器的姿态单位四元数中的标量部分，θ表示绕着欧拉轴转过的一个角度，

e_x,e_y,e_z代表欧拉轴三个方向上的旋转轴，且满足

表示q对时间的导数；Ω(ω)是斜对称矩阵，其形式如下：Among them, ω=[ω ₁ , ω ₂ , ω ₃ ] ^T represents the attitude angular velocity vector of the spacecraft relative to the inertial coordinate system in the body coordinate system, ω ₁ , ω ₂ , ω ₃ are the transverse direction of the spacecraft in the body system, respectively Angular velocity components on the roll, yaw and pitch axes;

Represents the derivative of ω with respect to time; J is the total inertia matrix of the spacecraft, simplified as a diagonal matrix J=diag(J ₁ , J ₂ , J ₃ ), J ₁ , J ₂ , J ₃ are the rotations around the main axis of inertia Inertia; S(ω) is an obliquely symmetric matrix of the form τ represents the output torque of the actuator; q=[q ₀ , q ₁ , q ₂ , q ₃ ] ^T represents the attitude unit quaternion of the spacecraft,

represents the scalar part of the attitude unit quaternion of the spacecraft, θ represents an angle rotated around the Euler axis,

e _x , e _y , e _z represent the rotation axes in the three directions of the Euler axis, and satisfy

Represents the derivative of q with respect to time; Ω(ω) is an obliquely symmetric matrix with the following form:

The model of the reaction flywheel combination is as follows:

为飞轮组合的角动量相对时间的导数，在由四个飞轮组成的飞轮组合中其角动量与转速的关系如下：where H _rw is the angular momentum of the flywheel combination,

is the derivative of the angular momentum of the flywheel combination with respect to time, and the relationship between its angular momentum and the rotational speed in a flywheel combination consisting of four flywheels is as follows:

H _rw =CJ _rw N

where C is the 3×4 flywheel installation matrix, N=[n ₁ ,n ₂ ,n ₃ ,n ₄ ] ^T is the angular velocity vector of the flywheel, n ₁ ,n ₂ ,n ₃ ,n ₄ represent the angular velocity of each flywheel respectively ; J _rw represents the moment of inertia matrix of the flywheel combination, and its form is J _rw =J _α I _4×4 , J _α represents the moment of inertia of a single flywheel, and I _4×4 is the fourth-order unit matrix;

Since the actuator controls the attitude by exchanging angular momentum with the spacecraft, the total angular momentum of the system is conserved:

H=H _rw +Jω

H is the total angular momentum of the system, which is a constant when there is no external torque.

Integrating the actuator model and the spacecraft attitude dynamics model yields:

The above model is discretized, and the sampling interval is set to Δt to obtain, at the kth time:

ω _k+1 =J ^-1 ΔtS(ω _k )H+J ^-1 J _α CΔN _k +ω _k

The subscript k represents the value of the corresponding variable at the kth time, ΔN _k =N _k −N _k-1 , and I _3×3 is a third-order unit matrix.

In step (2):

(a) Reaction flywheel maximum output torque constraint:

The torque output of the reaction flywheel is realized by changing the angular momentum of the flywheel, so it is expressed in the following form:

where T _max is the maximum output torque vector. After discretizing the above formula, we get:

After sorting, we get:

(b) Reaction flywheel maximum angular momentum constraint:

The angular momentum saturation of the flywheel is reflected in the fact that the rotational speed of the flywheel rotor reaches the upper limit, so the angular momentum saturation constraint can be expressed by the flywheel angular velocity constraint:

-N _max ≤ΔN _k +N _k-1 ≤N _max

where N _max is the maximum angular velocity vector of the flywheel.

After sorting, we get:

(c) Line-of-sight angle constraints for spacecraft attitude pointing:

Considering that the orientation of the spacecraft should avoid some conical sight areas, the designed attitude constraints are as follows:

Among them, α represents the unit vector that the spacecraft points to in the body coordinate system, β represents the unit vector of the direction the spacecraft needs to avoid in the body coordinate system, and θ represents the size of the sight angle of the avoidance area.

In step (3), the optimization objective function function V(x _k , u _k ) is expressed as:

where x _k =[q _k ,ω _k ] ^T represents the state quantity of the system, and the physical meaning is the attitude and angular velocity of the spacecraft; the change in the combined speed of the flywheel is used as the input, that is, u _k =ΔN _k ; _NP is the prediction of the MPC Scope; Q and P are the weight matrices of state variables and input variables. If Q is larger than P, it means that the optimization goal is more focused on stable time, and if Q is larger than P, it means that the optimization goal is more focused on energy consumption ;

Organized into concise form:

in

In step (4), the control problem is transformed into the problem of finding the extreme value of the objective function under constraints, and the control variable is solved by real-time iterative method. Transform the problem of solving the control input into the following mathematical problem:

x _i,k+1 =F( _xi,k ,ui _,k )

L(x _i,k ,u _i,k )≤0 _17×1 k=0,1,…,N _P -1

表示在当前时刻i系统的状态反馈，x_i,k表示由当前时刻状态x_i,0推算的第k时刻的系统状态；u_i,k表示第k时刻的预计输入，0_17×1表示17×1的零矩阵；F(x_i,k,u_i,k)为系统状态方程，具体表示为：in represents the optimization problem at time i,

Represents the state feedback of the system at the current moment i, x _i,k represents the system state at the kth moment calculated from the current moment state x _i,0 ; u _i,k represents the expected input at the kth moment, 0 _17×1 means 17 ×1 zero matrix; F(x _i,k ,ui _,k ) is the state equation of the system, specifically expressed as:

L(x _i,k ,u _i,k ) is the system constraint, specifically expressed as:

求解出来的求出来的

和

进而求解优化问题

具体分为两个阶段：准备阶段和响应阶段：The processing process of the method using real-time iteration is: the optimization problem is known at the previous moment

solved, solved

and

to solve the optimization problem

It is divided into two phases: the preparation phase and the response phase:

根据

和

计算敏感矩阵A_i,_k,B_i,k,C_i,k,D_i,k和误差l_i,k,r_i,k，具体表达式为：The process of the preparation stage is: shift x _i-1 and u _i-1 by a sampling moment, keep the last element unchanged, and get and

according to

and

Calculate the sensitivity matrices A _i , _k , B _i,k ,C _i,k ,D _i,k and errors l _i,k ,r _i,k , the specific expressions are:

By linearization, a nonlinearly constrained optimization problem is transformed into a linearly constrained quadratic programming problem:

Δu _i,k =u _i,k -u _i-1,k

Δx _i,k+1 =A _i,k Δx _i,k +B _i,k Δu _i,k +r _i,k

C _i,k Δx _i,k +D _i,k Δu _i,k +l _i,k ≤0 k=0,1,…,N _P -1

带入求解问题得到Δx_i,k和Δu_i,k通过式子：The process of the response phase is: to obtain the status feedback of the system

Bring in to solve the problem Obtain Δx _i,k and Δu _i,k by formula:

的解(x_i,k,u_i,k)，最后再将u_i＝u_i,0作为系统的控制输入。get question

(x _i,k ,ui _,k ), and finally u _i =ui _,0 is used as the control input of the system.

The model prediction control method for the attitude avoidance of the spacecraft considering the dynamic characteristics of the flywheel of the present invention can make the spacecraft perform attitude avoidance in the process of attitude maneuvering and the performance of the actuator is limited, and can satisfy the comprehensive optimization of energy and time at the same time. . The given optimization method can improve the solution efficiency of MPC. The idea of hot start is used to accelerate the speed of solving optimization problems. The idea of real-time iteration is divided into two stages can output the state feedback at the latest moment, which improves the system execution. Timeliness of control. In summary, the present invention has the following advantages:

(1) Integrate the model of the spacecraft and the model of the actuator, and use the change in the rotational speed of the flywheel as the input, which is different from the traditional control method of first obtaining the control torque and then assigning the control torque to each flywheel. Acting directly on the actuator, the process of distributing the control quantity is eliminated, and the entire system process is directly optimized.

(2) Compared with the existing potential function methods for attitude avoidance of spacecraft, the model predictive control method can simultaneously consider the performance constraints of the actuator and the attitude pointing constraints during the satellite attitude maneuvering process. The corresponding performance function is set so that the process of attitude maneuvering can achieve the comprehensive optimization of energy and time.

(3) Compared with the general attitude avoidance optimal control, the method using model predictive control is closed-loop control, which can re-optimize the trajectory according to the state feedback after each execution action, making the control robust sex.

(4) The optimization idea of real-time iteration and hot start used in the model predictive control of the present invention can reduce the amount of computation and optimization efficiency compared with the ordinary nonlinear model predictive control method, and can make a decision based on the system state at the latest moment. In response, the output of the control is more time-sensitive.

Description of drawings

Fig. 1 is a flow chart of a model predictive control method for spacecraft attitude avoidance considering flywheel dynamic characteristics;

Fig. 2 is a process block diagram for the control system to solve the control input;

Detailed ways

The specific implementation of the present invention will be described in detail below. It is necessary to point out that the following implementation is only used for the further description of the present invention, and should not be construed as a limitation on the protection scope of the present invention. Some non-essential improvements and adjustments made still belong to the protection scope of the present invention.

The present invention provides a model predictive control method for spacecraft attitude avoidance considering the dynamic characteristics of the flywheel. The attitude model of the spacecraft including the actuator is established as the prediction model of the MPC; secondly, the mathematical model of the attitude avoidance constraint is established according to the sight angle of the spacecraft mounted instruments and the avoidance vector, and the actuator constraint is established according to the angular momentum saturation and torque saturation characteristics of the flywheel. mathematical model. Then the corresponding optimization objective function is designed according to the task requirements, including the quadratic form of the actuator input and the quadratic form of the system state quantity error, so as to comprehensively consider the execution time and energy consumption. Finally, the control problem is transformed into a mathematical problem of finding the extreme value of the objective function under the constraints of the equations and states of the system dynamics equations and the input-limited inequality constraints; the real-time iterative optimization processing method is used to solve the problem quickly, and the optimized The solution is output as the control quantity of the system. The block diagram of the control system solving the control input process is shown in Figure 2.

The specific implementation steps are as follows:

The first step is to establish the mathematical model of the attitude control system. Combine the flywheel dynamics equations into the spacecraft attitude dynamics equations:

表示ω相对时间的导数；J是航天器总的惯量矩阵，简化表示为对角阵J＝diag(J₁,J₂,J₃)，J₁＝30kg/m²,J₂＝50kg/m²,J₃＝40kg/m²为绕惯量主轴的转动惯量；S(ω)是斜对称矩阵，其形式为τ表示执行机构的输出力矩；q＝[q₀,q₁,q₂,q₃]^T表示航天器的姿态单位四元数，

表示航天器的姿态单位四元数中的标量部分，与绕欧拉轴旋转的角度有关，θ表示绕着欧拉轴转过的一个角度，

e_x,e_y,e_z代表欧拉轴三个方向上的旋转轴，且满足

表示q对时间的导数；Ω(ω)是斜对称矩阵，其形式为

Among them, ω=[ω ₁ , ω ₂ , ω ₃ ] ^T represents the attitude angular velocity vector of the spacecraft relative to the inertial coordinate system in the body coordinate system, ω ₁ , ω ₂ , ω ₃ are the transverse direction of the spacecraft in the body system, respectively Angular velocity components on the roll, yaw and pitch axes;

Represents the derivative of ω relative to time; J is the total inertia matrix of the spacecraft, which is simplified as a diagonal matrix J=diag(J ₁ , J ₂ , J ₃ ), J ₁ =30kg/m ² , J ₂ =50kg/m ² , J ₃ =40kg/m ² is the moment of inertia around the main axis of inertia; S(ω) is an obliquely symmetric matrix whose form is τ represents the output torque of the actuator; q=[q ₀ , q ₁ , q ₂ , q ₃ ] ^T represents the attitude unit quaternion of the spacecraft,

represents the scalar part of the attitude unit quaternion of the spacecraft, which is related to the angle of rotation around the Euler axis, θ represents an angle rotated around the Euler axis,

e _x , e _y , e _z represent the rotation axes in the three directions of the Euler axis, and satisfy

represents the derivative of q with respect to time; Ω(ω) is an obliquely symmetric matrix of the form

The model of the reaction flywheel combination is as follows:

为H_rw对时间的导数，在由四个飞轮组成的飞轮组合中其角动量与转速的关系如下：where H _rw is the angular momentum of the flywheel combination,

is the derivative of H _rw with respect to time, and the relationship between its angular momentum and rotational speed in a flywheel combination consisting of four flywheels is as follows:

H _rw =CJ _rw N

in Install the matrix for the flywheel, N=[n ₁ , n ₂ , n ₃ , n ₄ ] ^T is the angular velocity vector of the flywheel, n ₁ , n ₂ , n ₃ , n ₄ represent the angular velocity of each flywheel respectively; J _rw represents the flywheel The combined moment of inertia matrix is in the form of J _rw =J _α I _4×4 , J _α =0.01608kg/m ² represents the moment of inertia of a single flywheel, and I _4×4 is a fourth-order unit matrix.

Since the actuator controls the attitude by exchanging angular momentum with the spacecraft, the total angular momentum of the system is conserved:

H=H _rw +Jω

H is the total angular momentum of the system, which is a constant when there is no external torque disturbance.

Integrating the actuator model and the spacecraft attitude dynamics model yields:

The above model is discretized, and the sampling interval is set to Δt=0.2s to obtain, at the kth moment:

ω _k+1 =J ^-1 ΔtS(ω _k )H+J ^-1 J _α CΔN _k +ω _k

The subscript k represents the value of the corresponding variable at the kth time, ΔN _k =N _k -N _k-1 , and I _3×3 is a third-order identity matrix.

The initial attitude of the spacecraft, the initial angular velocity and the initial angular velocity of the flywheel combination are q _intial =[-0.9524,-0.3048,0,0] ^T ω _intial =[0,0,0] ^T and N _intial =[0,0 ,0,0] ^T , the desired attitude and angular velocity of the control are respectively q _d =[1,0,0,0] ^T ω _d =[0,0,0] ^T . Therefore, the total initial angular momentum of the system is H=0.

The second step is to establish the attitude avoidance constraints of the spacecraft during the attitude maneuvering process and the performance constraints of the actuators:

(1) Reaction flywheel maximum output torque constraint:

The torque output of the reaction flywheel is realized by changing the angular momentum of the flywheel, so it is expressed in the following form:

where T _max =[1 1 1 1] ^T Nm is the maximum output torque vector. After discretizing the above formula, we get:

After sorting, we get:

(2) Reaction flywheel maximum angular momentum constraint:

The angular momentum saturation of the flywheel is reflected in the fact that the rotational speed of the flywheel rotor reaches the upper limit, so the angular momentum saturation constraint can be expressed by the flywheel angular velocity constraint:

-N _max ≤ΔN _k +N _k-1 ≤N _max

Wherein N _max =[200π 200π 200π 200π] ^T rad/s is the maximum angular velocity vector.

After sorting, we get:

(3) The line-of-sight angle constraint of the spacecraft attitude pointing:

Considering that the orientation of the spacecraft should avoid some conical sight areas, the designed attitude constraints are as follows:

表示航天器在本体坐标系下需规避方向的单位向量，

表示规避区域的视线角的大小。where α=[0 0 1] ^T indicates that the spacecraft points to the unit vector in the body coordinate system,

is the unit vector representing the direction the spacecraft needs to avoid in the body coordinate system,

Indicates the size of the sight angle of the avoidance area.

The third step is to design the optimization objective function based on the model predictive controller:

Where x _k =[q _k ,ω _k ] ^T represents the state quantity of the system, and the physical meaning is the attitude and angular velocity of the spacecraft; the change in the combined speed of the flywheel is used as the input, that is, u _k =ΔN _k ; _NP =5 is the MPC The prediction range of ; Q=I _7×7 and P=I _4×4 are the weight matrices of state variables and input variables, and I _7×7 is the 7th-order unit matrix. If Q is larger than P, it means that the optimization goal is more important In the stable time, if Q is larger than P, it means that the optimization goal is more focused on energy consumption. Organized into concise form:

in

The fourth step is to convert the control problem into the problem of finding the extreme value of the objective function under constraints, and use the real-time iterative method to solve the control variables. The optimization objective function of , the problem of solving the control input is transformed into the following mathematical problem:

x _i,k+1 =F( _xi,k ,ui _,k )

L(x _i,k ,u _i,k )≤0 _17×1 k=0,1,…,N _P -1

表示在时刻i的最优化问题，

表示在当前时刻i系统的状态反馈，x_i,k表示由当前时刻状态x_i,0推测在第k时刻的系统状态；u_i,k表示在第k时刻的预计输入；F(x_i,k,u_i,k)为第一步中的系统状态方程具体表示为：in

represents the optimization problem at time i,

Represents the state feedback of the system at the current moment i, _xi,k represents the system state at the kth moment estimated from the current moment state _xi,0 ; u _i,k represents the expected input at the kth moment; F( _{xi, k} ,ui _,k ) is the system state equation in the first step, which is specifically expressed as:

L(x _i,k ,ui _,k ) is the system constraint in the second step, specifically expressed as:

求解出来的求出来的

和

进而求解优化问题

的过程，具体分为两个阶段：准备阶段和响应阶段The processing process of the method using real-time iteration is: the previous optimization problem is known

solved, solved

and

to solve the optimization problem

The process is divided into two phases: the preparation phase and the response phase

和

根据

和

计算敏感矩阵A_i,k,B_i,k,C_i,k,D_i,k和误差l_i,k,r_i,k，具体表达式为：The process of the preparation stage is: shift x _i-1 and u _i-1 by a sampling moment, and the last position remains unchanged to obtain

and

according to

and

Calculate the sensitivity matrices A _i,k ,B _i,k ,C _i,k ,D _i,k and errors l _i,k ,r _i,k , the specific expressions are:

Through the above linearization, a nonlinear constrained optimization problem is transformed into a linear constrained quadratic programming problem:

Δu _i,k =u _i,k -u _i-1,k

Δx _i,k+1 =A _i,k Δx _i,k +B _i,k Δu _i,k +r _i,k

C _i,k Δx _i,k +D _i,k Δu _i,k +l _i,k ≤0 k=0,1,…,N _P -1

带入求解问题

得到Δx_i,k和Δu_i,k通过式子：The process of the response phase is: to obtain the status feedback of the system

Bring in to solve the problem

Obtain Δx _i,k and Δu _i,k by formula:

get question the solution of (x _i,k ,ui _,k ). Finally, u _i =u _i,0 is used as the control input of the system.

Using the above control method, the spacecraft can perform attitude avoidance in the process of attitude maneuvering and the performance of the actuator is limited, and can meet the comprehensive optimization of energy and time at the same time. The given optimization method can improve the solution efficiency of MPC. The idea of hot start is used to accelerate the speed of solving optimization problems. The idea of real-time iteration is divided into two stages can make control output according to the state feedback at the latest moment, which improves the system performance. The timeliness of execution control, the optimized process is shown in Figure 2.

Although exemplary embodiments of the present invention have been described for purposes of illustration, workers skilled in the art will recognize that changes may be made in form and detail without departing from the scope and spirit of the invention as disclosed in the accompanying claims. Carry out various modifications, additions and substitutions, etc., and all these changes should belong to the protection scope of the appended claims of the present invention, and the various steps in the various departments and methods of the products claimed in the present invention can be combined arbitrarily. form together. Accordingly, the description of the embodiments disclosed in the present invention is not intended to limit the scope of the present invention, but to describe the present invention. Accordingly, the scope of the present invention is not limited by the above embodiments, but is defined by the claims or their equivalents.

Claims (1)

1. A model prediction control method for spacecraft attitude evasion considering dynamic characteristics of a flywheel is characterized by comprising the following steps:

(1) establishing a spacecraft attitude model containing an actuating mechanism as a prediction model of the MPC according to the attitude kinematics and the dynamic characteristics of the spacecraft and the rotation dynamic characteristics of the flywheel;

(2) establishing a mathematical model of attitude avoidance constraint according to a line-of-sight angle and an avoidance vector of a spacecraft carrying instrument, and establishing a mathematical model of actuating mechanism constraint according to angular momentum saturation and moment saturation characteristics of a flywheel;

(3) designing a corresponding optimization objective function according to task requirements, wherein the optimization objective function comprises a quadratic form input by an execution mechanism and a quadratic form of system state errors so as to comprehensively consider execution time and energy consumption;

(4) converting the control problem into an extreme value of an objective function under the constraint conditions of an equality and a state of a system dynamic equation and an inequality with limited input, rapidly solving by using a real-time iterative optimization processing method, and outputting a solution obtained by optimization as a control quantity of the system;

combining a kinetic equation of the flywheel into a spacecraft attitude kinetic equation, performing discretization processing, and establishing a prediction model of the MPC;

the spacecraft attitude dynamics model in the step (1) is expressed as follows:

wherein ω is [ ω ═ ω [ [ ω ]₁,ω₂,ω₃]^TRepresenting the attitude angular velocity vector omega of the spacecraft relative to an inertial coordinate system under a body coordinate system₁,ω₂,ω₃Angular velocity components of the spacecraft about the roll axis, yaw axis and pitch axis in the system, respectively;

represents the derivative of ω with time; j is the total inertia matrix of the spacecraft, simplified as the diagonal matrix J ═ diag (J)₁,J₂,J₃)，J₁,J₂,J₃Is the moment of inertia around the inertia main shaft; s (omega) is a skew symmetric matrix of the formτ represents the output torque of the actuator; q ═ q₀,q₁,q₂,q₃]^TRepresents a unit quaternion of the attitude of the spacecraft,

represents the scalar part of the quaternion of the attitude unit of the spacecraft, theta represents an angle of rotation about the euler axis,

e_x,e_y,e_zrepresents the rotation axis in three directions of the Euler axis and satisfies

Represents the derivative of q with respect to time; Ω (ω) is a diagonally symmetric matrix of the form:

the model of the reaction flywheel combination is as follows:

wherein H_rwIs the angular momentum of the flywheel combination,

for its derivative with respect to time, the angular momentum is related to the rotational speed in a flywheel combination consisting of four flywheels as follows:

H_rw＝CJ_rwN

where C is a 3 × 4 flywheel mounting matrix and N ═ N₁,n₂,n₃,n₄]^TIs the angular velocity vector of the flywheel, n₁,n₂,n₃,n₄Respectively representing the angular velocity of each flywheel; j. the design is a square_rwA moment of inertia matrix, of the form J, representing a flywheel assembly_rw＝J_αI_4×4，J_αRepresenting the moment of inertia of a single flywheel, I_4×4Is an identity matrix of 4 orders;

since the actuator controls attitude by exchanging angular momentum with the spacecraft, the total angular momentum of the system is conserved:

H＝H_rw+Jω

h is the total angular momentum of the system and is a constant when no external moment interference exists;

integrating the model of the actuating mechanism and the model of the spacecraft attitude dynamics to obtain:

discretizing the model, setting a sampling interval to be delta t, and obtaining the model at the kth moment by:

ω_k+1＝J^-1ΔtS(ω_k)H+J^-1J_αCΔN_k+ω_k

where the index k denotes the value of the corresponding variable at the k-th instant, Δ N_k＝N_k-N_k-1，I_3×3Is a third-order identity matrix;

in the step (2):

(1) constraint of maximum output torque of reaction flywheel:

the torque output of the reaction flywheel is achieved by changing the angular momentum of the flywheel, so it is expressed in the form:

wherein T is_maxIs the maximum output torque vector; discretizing the formula to obtain:

after finishing, obtaining:

(2) constraint of maximum angular momentum of reaction flywheel:

the angular momentum saturation of the flywheel is reflected in that the rotational speed of the flywheel rotor reaches an upper limit, so the angular momentum saturation constraint can be expressed by the flywheel angular speed constraint:

-N_max≤ΔN_k+N_k-1≤N_max

wherein N is_maxIs the maximum angular velocity vector of the flywheel;

after finishing, obtaining:

(3) and (3) restricting the view angle of the spacecraft attitude pointing:

considering the pointing direction of the spacecraft to avoid certain conical sight areas, the designed attitude is constrained in the following form:

α represents a unit vector pointed by the spacecraft in a body coordinate system, β represents a unit vector of a direction to be avoided of the spacecraft in the body coordinate system, and theta represents the size of a line-of-sight angle of an avoidance area;

optimizing the objective function V (x) in the step (3)_k,u_k) Expressed as:

wherein x_k＝[q_k,ω_k]^TRepresenting the state quantity of the system, wherein the physical meaning is the attitude and the angular speed of the spacecraft; using flywheel combined speed variation as input, i.e. u_k＝ΔN_k；N_PA prediction range for MPC; q and P are weight matrixes of the state variables and the input variables, if Q is larger than P, the optimization target is more focused on stable time, and if Q is larger than P, the optimization target is more focused on energy consumption;

arranging into a concise form:

wherein

Converting the control problem into an extreme value problem of solving an objective function under a constraint condition, and solving a control quantity by using a real-time iteration method;

converting the problem of solving the control input into the following mathematical problem:

x_i,k+1＝F(x_i,k,u_i,k)

L(x_i,k,u_i,k)≤0_17×1k＝0,1,…,N_P-1

wherein

Is shown at the timeThe optimization problem of the moment i is solved,

indicating the state feedback of the system at the current time i, x_i,kIndicating the state x from the current moment_i,0The estimated system state at the k-th time; u. of_i,kIndicating the expected input at time k, 0_17×1Represents a zero matrix of 17 × 1; f (x)_i,k,u_i,k) The system state equation is specifically expressed as:

L(x_i,k,u_i,k) For system constraints, it is specifically expressed as:

the process of using the real-time iterative method is as follows: knowing the optimization problem at the previous moment

Solved for

And

and then solve the optimization problem

The method is specifically divided into two stages: preparation phase and response phase:

the preparation stage comprises the following processes: x is to be_i-1And u_i-1Shifting and keeping the last element unchanged to obtain

And

according toAnd

calculating the sensitivity matrix A_i,k,B_i,k,C_i,k,D_i,kAnd error l_i,k,r_i,kThe specific expression is as follows:

converting a nonlinear constrained optimization problem into a linear constrained quadratic programming problem through linearization:

Δu_i,k＝u_i,k-u_i-1,k

Δx_i,k+1＝A_i,kΔx_i,k+B_i,kΔu_i,k+r_i,k

C_i,kΔx_i,k+D_i,kΔu_i,k+l_i,k≤0k＝0,1,…,N_P-1

the response phase process is as follows: of the acquisition systemState feedback

Carry-over solving problem

To obtain Δ x_i,kAnd Δ u_i,kBy the formula:

get a problemSolution of (x)_i,k,u_i,k) Finally, u is again added_i＝u_i,0As a control input to the system.

Images