# Application of the Improved Grey Wolf Algorithm in Spacecraft Maneuvering Path Planning

Authors: Changqing Wu; Xiaodong Han; Weiyu An; Jianglei Gong; Nan Xu
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

In many space missions, spacecraft are required to have the ability to avoid various obstacles and finally reach the target point. In this paper, the path planning of spacecraft attitude maneuver under boundary constraints and pointing constraints is studied. The boundary constraints and orientation constraints are constructed as finite functions of path evaluation. From the point of view of optimal time and shortest path, the constrained attitude maneuver problem is reduced to optimal time and path solving problem. To address this problem, a metaheuristic maneuver path planning method is proposed (cross-mutation grey wolf algorithm (CMGWO)). In the CMGWO method, we use angular velocity and control torque coding to model attitude maneuver, which increases the difficulty of solving the problem. In order to deal with this problem, the grey wolf algorithm is used for mutation and evolution, so as to reduce the difficulty of solving the problem and shorten the convergence time. Finally, simulation analysis is carried out under different conditions, and the feasibility and effectiveness of the method are verified by numerical simulation.


---

## Body

## 1. Introduction

As a low-cost spacecraft, small satellites have the characteristics of light, agility, and large coverage. They have become an ideal spacecraft for studying space, which has put forward new requirements for the movement of satellites [1]. With the development of space exploration technology and space commercial activities, the number of satellites in space has increased dramatically, which will pose a great threat to the safety of existing satellites and the next generation of space missions. There were 19,404 large objects and millions of pieces of debris in Earth orbit, which leads to huge challenges for space resources and the environment [2]. Therefore, owing to ensuring the safety of the spacecraft, it is necessary to avoid the debris when the spacecraft is flying. With the continuous development of space missions, it is extremely urgent to require the ability of rapid and large-angle attitude change, especially for the application of satellite positioning, disaster warning, scientific exploration, and other related tasks [1, 3].When the spacecraft is maneuvering, it is necessary to prevent the light of bright celestial bodies (such as the sun) from entering the field of view of some optical sensors (such as infrared sensors or low light sensitive elements). Otherwise, the optical sensors will be temporarily blind or damaged. In the maneuvering, the direction vector of the solar array must be kept within a certain required range, so as to continuously provide electric energy. These directional constraints greatly limit the feasible region of spacecraft attitude maneuver [4]. In addition, the limitation of angular velocity and control torque will also affect the attitude maneuver path [5, 6]. Due to these complex constraints, time-optimized maneuvering path planning based on spacecraft guidance, navigation, and control systems is very challenging.We transform the satellite orbit planning problem into the optimization problem of finding the best time and shortest path and establish the spacecraft attitude maneuver model by angular velocity-time-control torque coding. Compared with the latest and most advanced research [7, 8], our control model also takes into account more comprehensive constraints. The model proposed in this paper has better performance in solving the subsequent optimal path. In recent years, the metaheuristic algorithm has become more and more popular because of its simplicity advantages. In addition, the rapid development of computing technology makes up for the cost of computing. Many metaheuristic algorithms have been proposed. They can be subdivided into two categories: bionic algorithm and physical algorithm. The first is related to animal biology and behavior, while the second is based on the laws of physics. Among the most commonly used algorithms are the genetic algorithm [9], which is inspired by Darwin’s theory of evolution; ant colony algorithm [10], which is inspired by the behavior of ants; particle swarm optimization algorithm [11], based on human social behavior and flock of birds; and differential evolution algorithm [12], based on the concept of mutation, reorganization, and selection applicable to swarms. Secondly, Takahama and Sakai proposed an improved constrained DE (differential evolution) algorithm [13] to generate the optimal feasible path. Furthermore, an artificial bee colony algorithm combined with evolutionary programming [14], fuzzy quadtree framework [15], improved traveling salesman problem algorithm [16], and improved pulse coupled neural network model [17] is also applicable to the path planning of agents. There are other algorithms for solving similar optimization problems [18–22]. Moreover, researchers have implemented many methods to find the optimal path in different environments. Yang et al. [23] proposed a finite angle A∗ algorithm for path planning on satellite images. The algorithm uses a modified line of sight with a branching factor of 16 to improve the efficiency of the finite angle A∗ algorithm. Duchon [24] and others proposed an improved A∗ algorithm based on theta∗, phi∗, rectangular symmetry reduction, and hop search, which was implemented in symmetric and asymmetric environments. Lou and Tang [25] summarized the motion planning algorithms from the perspective of UGV guidance. The algorithm of path planning is discussed from the aspects of complexity and mapping with existing standards. One problem of metaheuristic algorithm is how to set some internal parameters correctly, such as maximum iteration times, population size, and internal coefficients. These parameters may change according to different problems that need to be optimized, because they will significantly affect the results of the optimization process. Therefore, this limits the use of such algorithms, which are only used by expert users in this field.We have noticed that various improved grey wolf algorithms are being applied in various fields to solve various problems. For example, the novel random walking grey wolf optimizer proposed a new crossover operator-double distribution crossover operator (DDX). The performance of DDX is compared with that of the existing real-coded crossover operator, namely, Laplace crossover operator [26]. The random walk grey wolf optimizer for constrained engineering optimization problems proposes a simple constraint processing technique-constrained version GWO and proposes random walk GWO (RW-GWO) [27] by pointing out some shortcomings of the original GWO leaders in the process of searching prey. A new algorithm OCS-GWO is proposed by the opposition-based chaotic grey wolf optimizer for global optimization task. The algorithm effectively uses the search area by introducing opposition-based learning to approach the search candidate solution closer to the global optimum and chaotic local search, thus improving the performance of the original GWO algorithm [28]. The Cauchy grey wolf optimizer for continuous optimization problems proposes an improved classical GWO algorithm, which is named Cauchy-GWO. In the Cauchy-GWO operator, two new wolves are generated by the Cauchy distribution random number, and then, another new wolf is generated by convex combination of these new wolves [29]. The grey wolf optimizer with enhanced leadership inspiration for global optimization problems proposes an improved leadership-based GWO, namely, GLF-GWO. In GLF-GWO, the leader is updated through the Levy flight search mechanism. The proposed GLF-GWO algorithm improves the search efficiency of leading hunters in GWO and provides better guidance for accelerating the search process of GWO [30]. A hybrid grey wolf optimizer with a mutation operator hybridizes the grey wolf optimizer with differential evolution (DE) mutation and puts forward two versions of DE-GWO and GDE-GWO to avoid the stagnation of solution [31], as well as the reliability-redundancy allocation using the random walking grey wolf optimizer [32], memory-based grey wolf optimizer for global optimization task [33], and over-current relay optimization using the improved leadership-inspired grey wolf optimizer [34].Because the traditional optimization algorithm has the disadvantages of slow convergence, inspired by predecessors [35, 36], we optimized the grey wolf algorithm [37] and applied it to satellite orbit planning. The grey wolf algorithm is a swarm intelligence optimization algorithm which simulates the predation behavior of grey wolves. The grey wolf algorithm has the characteristics of simple operation, few adjustment parameters, and easy programming. Compared with other swarm intelligence optimization algorithms, it has obvious advantages in function optimization. However, at the same time, it also has some shortcomings such as being easy to fall into the local optimum and low solution accuracy. Inspired by previous studies, this paper proposes a cross-mutation grey wolf algorithm in view of the shortcomings of the current algorithm. The algorithm uses the basic GWO algorithm to calculate the fitness value of each grey wolf. With the update of the wolf position, the individual in the wolf group constantly changes, which makes the population constantly update and effectively solves the problem of entering the local optimum. At the same time, the high efficiency of the grey wolf algorithm allows the time for the algorithm to be controlled. In order to improve the search ability of the algorithm, the fitness values of the grey wolf individuals are compared, and a penalty function is imposed on the grey wolf individuals with poor fitness values. By testing five benchmark functions, the results are compared with the traditional GWO algorithm to verify the superiority of the algorithm.


## 2. Spacecraft Constraints

In this section, how to formulate the satellite’s maneuvering path into a behavior schedule optimization problem is presented.Let us imagine a spacecraft working on the space station, some of those parameters are shown in Table1.Table 1
Spacecraft operating parameters.
The angular velocitywControl torqueuAttitude (quaternion)qThe moment of inertiaJSpacecraft to the object attitude cosine matrixRThe angle between some instruments and the sunαThe orientation angle of the instrument itselfβOnly when some constraints [38, 39] are met can the spacecraft work better in space. That is, only under these conditions can the mission be carried out. Therefore, constraints are a prerequisite for planning trajectories.


### 2.1. Boundary Constraint

In the implementation process, in order to prevent the control torque from exceeding the boundary and causing adverse effects on the satellite, we define that the control torque is bounded. Furthermore, the following constraints should be satisfied:u≤γu,i=1,2,3.During attitude maneuver, the angle of some sensors should not rotate too much; otherwise, it will produce excessive force. Therefore, the angular velocity needs to be limited within a certain range and the following constraints should be satisfied:w≤γw,i=1,2,3.Boundary constraints are the prerequisite to ensure that it works properly, so when the model is initialized, boundary constraints should be added, such as setting a range of angular velocity:−w to w, and keeping the angular velocity behind in this range.


### 2.2. Pose Kinematic and Dynamic Constraints

When it comes to the attitude of rigid spacecraft, if the four-dimensional space is considered, the unit quaternion is needed to represent the maneuver value of the parameters, and the purpose of avoiding singular value can be achieved. According to the relationship between moment of inertia, angular velocity, and torque, the constraint is(1)J=diagJ1,J2,J3,(2)w=w1,w2,w3T,(3)w×=0−w3w2w30−w1−w2w10,(4)u=u1,u2,u3T,whereJ is the moment of inertia, defined as the diagonal matrix of J1,J2,J3; w is the angular velocity, w× prime is defined as the cross-product matrix of w, and u is the control torque.After the calculation, the spacecraft should satisfy(5)Jw.=u−w×Jw.Andq=1/2Lw and q̇=1/2Mq.
(6)L=−q1−q2−q3q0−q3q2q3q0−q1−q2q1q0.(7)M=0−w1−w2−w3w10w3−w2w2−w30w1w3w2−w10.In conclusion, the attitude constraint is(8)Jw.=u−w×Jw,q̇=12Mq.


### 2.3. Pointing to the Constraint

The attitude orientation constraint is the constraint that the angle relation of a certain angle relative to a certain attitude must satisfy. It depends not only on the angle of its own devices but also on the external environment, such as the position of the sun. If the limit is exceeded, the spacecraft will be affected. It is generally classified into two types: forbidden constraint and mandatory constraint.


#### 2.3.1. Forbidden Constraint

Some of the instruments in the spacecraft during the mission must avoid bright light; otherwise, it will damage the device. In mathematical terms, the instrument and the sun should be greater than the angle of view of the instrument itself. That is the forbidden constraint.m represents the direction vector of an optical sensitive element in the body coordinate system, and n=n1,n2,n3T represents the direction vector of a bright celestial body in the inertial system. α is the field of view of the sensitive instrument. Constraints should satisfy mTCBn≤cosα, among them:
(9)m=m1,m2,m3T.R is the cosine matrix from the spacecraft inertial coordinate system to the object attitude. CB is the attitude cosine matrix of spacecraft from its own system to the inertial system. After simplification, the above inequality can be expressed as qTHq≤0. After calculation,
(10)H=mTn−cosαm×nTTm×nTnmT+mnT−nTm+cosαI3.Among them,m× is the cross-product of m. Therefore, the forbidden constraint is qTHq≤0.


#### 2.3.2. Mandatory Constraint

On the other hand, some instruments may need to point at a directional angle when operating. For example, taking the solar panel of spacecraft as an example,vB represents the direction vector of the solar panel in the body coordinate system, and the angle between n and vB should be less than β. Namely, vBTCBn≥cosβ.This is called a mandatory constraint. The constraints are as follows:qTKq≥0:
(11)K=nTvB−cosβvB×nTTvB×nTnvBT+vBnT−nTvB+cosβI3.All in all, the spacecraft needs to satisfy motion constraints, boundary constraints, forbidden constraints, and mandatory constraints during the mission. This is also a prerequisite for the spacecraft to complete its mission.To sum up, the constraint to be satisfied is(12)Jw.=u−w×Jw,q̇=12Mq,ui≤γu,i=1,2,3,wi≤γw,i=1,2,3,qTHq≤0,qTKq≥0.


## 2.1. Boundary Constraint

In the implementation process, in order to prevent the control torque from exceeding the boundary and causing adverse effects on the satellite, we define that the control torque is bounded. Furthermore, the following constraints should be satisfied:u≤γu,i=1,2,3.During attitude maneuver, the angle of some sensors should not rotate too much; otherwise, it will produce excessive force. Therefore, the angular velocity needs to be limited within a certain range and the following constraints should be satisfied:w≤γw,i=1,2,3.Boundary constraints are the prerequisite to ensure that it works properly, so when the model is initialized, boundary constraints should be added, such as setting a range of angular velocity:−w to w, and keeping the angular velocity behind in this range.


## 2.2. Pose Kinematic and Dynamic Constraints

When it comes to the attitude of rigid spacecraft, if the four-dimensional space is considered, the unit quaternion is needed to represent the maneuver value of the parameters, and the purpose of avoiding singular value can be achieved. According to the relationship between moment of inertia, angular velocity, and torque, the constraint is(1)J=diagJ1,J2,J3,(2)w=w1,w2,w3T,(3)w×=0−w3w2w30−w1−w2w10,(4)u=u1,u2,u3T,whereJ is the moment of inertia, defined as the diagonal matrix of J1,J2,J3; w is the angular velocity, w× prime is defined as the cross-product matrix of w, and u is the control torque.After the calculation, the spacecraft should satisfy(5)Jw.=u−w×Jw.Andq=1/2Lw and q̇=1/2Mq.
(6)L=−q1−q2−q3q0−q3q2q3q0−q1−q2q1q0.(7)M=0−w1−w2−w3w10w3−w2w2−w30w1w3w2−w10.In conclusion, the attitude constraint is(8)Jw.=u−w×Jw,q̇=12Mq.


## 2.3. Pointing to the Constraint

The attitude orientation constraint is the constraint that the angle relation of a certain angle relative to a certain attitude must satisfy. It depends not only on the angle of its own devices but also on the external environment, such as the position of the sun. If the limit is exceeded, the spacecraft will be affected. It is generally classified into two types: forbidden constraint and mandatory constraint.


### 2.3.1. Forbidden Constraint

Some of the instruments in the spacecraft during the mission must avoid bright light; otherwise, it will damage the device. In mathematical terms, the instrument and the sun should be greater than the angle of view of the instrument itself. That is the forbidden constraint.m represents the direction vector of an optical sensitive element in the body coordinate system, and n=n1,n2,n3T represents the direction vector of a bright celestial body in the inertial system. α is the field of view of the sensitive instrument. Constraints should satisfy mTCBn≤cosα, among them:
(9)m=m1,m2,m3T.R is the cosine matrix from the spacecraft inertial coordinate system to the object attitude. CB is the attitude cosine matrix of spacecraft from its own system to the inertial system. After simplification, the above inequality can be expressed as qTHq≤0. After calculation,
(10)H=mTn−cosαm×nTTm×nTnmT+mnT−nTm+cosαI3.Among them,m× is the cross-product of m. Therefore, the forbidden constraint is qTHq≤0.


### 2.3.2. Mandatory Constraint

On the other hand, some instruments may need to point at a directional angle when operating. For example, taking the solar panel of spacecraft as an example,vB represents the direction vector of the solar panel in the body coordinate system, and the angle between n and vB should be less than β. Namely, vBTCBn≥cosβ.This is called a mandatory constraint. The constraints are as follows:qTKq≥0:
(11)K=nTvB−cosβvB×nTTvB×nTnvBT+vBnT−nTvB+cosβI3.All in all, the spacecraft needs to satisfy motion constraints, boundary constraints, forbidden constraints, and mandatory constraints during the mission. This is also a prerequisite for the spacecraft to complete its mission.To sum up, the constraint to be satisfied is(12)Jw.=u−w×Jw,q̇=12Mq,ui≤γu,i=1,2,3,wi≤γw,i=1,2,3,qTHq≤0,qTKq≥0.


## 2.3.1. Forbidden Constraint

Some of the instruments in the spacecraft during the mission must avoid bright light; otherwise, it will damage the device. In mathematical terms, the instrument and the sun should be greater than the angle of view of the instrument itself. That is the forbidden constraint.m represents the direction vector of an optical sensitive element in the body coordinate system, and n=n1,n2,n3T represents the direction vector of a bright celestial body in the inertial system. α is the field of view of the sensitive instrument. Constraints should satisfy mTCBn≤cosα, among them:
(9)m=m1,m2,m3T.R is the cosine matrix from the spacecraft inertial coordinate system to the object attitude. CB is the attitude cosine matrix of spacecraft from its own system to the inertial system. After simplification, the above inequality can be expressed as qTHq≤0. After calculation,
(10)H=mTn−cosαm×nTTm×nTnmT+mnT−nTm+cosαI3.Among them,m× is the cross-product of m. Therefore, the forbidden constraint is qTHq≤0.


## 2.3.2. Mandatory Constraint

On the other hand, some instruments may need to point at a directional angle when operating. For example, taking the solar panel of spacecraft as an example,vB represents the direction vector of the solar panel in the body coordinate system, and the angle between n and vB should be less than β. Namely, vBTCBn≥cosβ.This is called a mandatory constraint. The constraints are as follows:qTKq≥0:
(11)K=nTvB−cosβvB×nTTvB×nTnvBT+vBnT−nTvB+cosβI3.All in all, the spacecraft needs to satisfy motion constraints, boundary constraints, forbidden constraints, and mandatory constraints during the mission. This is also a prerequisite for the spacecraft to complete its mission.To sum up, the constraint to be satisfied is(12)Jw.=u−w×Jw,q̇=12Mq,ui≤γu,i=1,2,3,wi≤γw,i=1,2,3,qTHq≤0,qTKq≥0.


## 3. Problem Statement

After studying the mathematical expression of the constraint, the improved grey wolf optimizer can be applied to solve the optimal control problem, and the optimal path can be screened through the improved crossover mutation. The evaluation function is set to the minimum value type, and the optimal path is obtained by searching for the minimum value of the evaluation function. The evaluation function is mainly divided into four parts: maneuver time, path direction constraint, energy bounded constraint, and terminal constraint.(13)Y1=t,(14)Y2=exp−qnTHfqn−1incaseqnTHfqn<0,otherwiseitis0,(15)Y=∑n=1nmaxY1+Y2.The total function is simplified toy=Y+Y1, a function to be minimized.Therefore, we transformed the time-constrained attitude maneuver problem into the optimization problem of finding the minimum value vectorm to minimize the total evaluation function and then adopted the improved GWO algorithm for optimization iteration in order to find the optimal solution.In this paper, an improved GWO algorithm is proposed. Fitness function based on motion constraints and boundary constraints is established, and the optimal path is screened through cross-variation, which is conducive to solving the problem of the local optimum in the previous trajectory planning. Last but not least, cross-variation can provide a possibility for the spacecraft to jump out of the local optimum.


## 4. Cross-Mutation Grey Wolf Optimizer (CMGWO) Method

Since the grey wolf algorithm was proposed in 2014, it has been widely used in various engineering fields due to its advantages of simple parameter setting and strong optimization ability. As a young intelligent optimization algorithm, the grey wolf algorithm inevitably has similar defects. However, as far as the current research status is concerned, there is still room for further improvement in both the convergence accuracy and the convergence speed of the grey wolf algorithm. Therefore, it is necessary to improve the grey wolf algorithm in order to make it obtain more ideal optimization results.There is a common problem in swarm intelligence optimization algorithms: the adjustment of local search performance and global search performance. An excellent swarm intelligence algorithm must have an excellent mechanism to adjust the global search performance and the local search performance. As a new member of the swarm intelligence algorithm, the grey wolf algorithm also has similar problems. Therefore, a good balance between global search performance and local search performance is indispensable for the grey wolf algorithm. Inspired by the differential evolution optimization algorithm, we merge the cross-mutation part of the algorithm into the iterative update process of the grey wolf algorithm. With the update of the position of the wolves, the individuals in the wolves are constantly mutating, so that the population is constantly updated, which can effectively solve the problem of adjusting the local search performance and global search performance of the original algorithm.


### 4.1. Improved Cross-Mutation Grey Wolf Algorithm

After comparing the optimization results of different algorithms, how to avoid falling into the local optimization and get the optimization results faster is only the main problem concerned. In this paper, a new grey wolf algorithm called the grey wolf genetic algorithm is proposed. This strategy is enlightened from two aspects of genetic mutation and grey wolf algorithm. With the updating of the wolf position, the individuals in the wolf group continue to mutate, so that the population keeps updating and the problem of getting into the local optimum can be solved effectively. Meanwhile, the high efficiency of the grey wolf algorithm makes the time of algorithm implementation controlled, and the improved grey wolf algorithm flow is shown in Figure1.Figure 1
Algorithm flow diagram.


#### 4.1.1. Grey Wolf Algorithm

According to the previous idea, the common grey wolf is described as follows: at the top of the grey wolf hierarchy is the supreme leader of the grey wolf tribe, known as alphaα wolf, who is responsible for overall planning and decisions on hunting. The second layer of the grey wolf’s social hierarchy is the second leader of the wolf pack, which is called beta β wolf. It mainly plays an assisting role in the overall planning and decision-making of alpha α wolf and is second only to alpha α wolf in the domination of the wolf population. At the third level of the grey wolf social hierarchy is delta δ wolf, also known as scouts, who are responsible for reconnaissance, watchkeeping, and guarding. The fourth level of the grey wolf social rank is the lowest in the grey wolf tribe, known as the omega (ω) wolf, which is essential for balancing the internal relations of the grey wolf tribe. The most suitable solution is considered to be the alpha (α) wolf. Therefore, the second and third best solutions were named beta β and omega ω wolves. The remaining candidate solutions were assumed to be omega wolves. And the three groups instructed the other wolves (omega) to search for the target. During the optimization process, the wolves update the positions of alpha, beta, delta, and omega; the grey wolf social hierarchy in the grey wolf algorithm is shown in Figure 2.Figure 2
Grey wolf social hierarchy.When the grey wolf population surrounds its prey, its predation behavior is defined as follows:(16)D⇀=C⇀·X⇀pt−X⇀t,(17)X⇀t+1=X⇀pt−A⇀·D⇀.Among them,D⇀ represents the distance between the individual wolf and the prey. The second formula is the updated formula of the individual wolf’s position. X⇀pt is the position vector of the prey in generation t. X⇀t is the position vector of the individuals in the t generation wolf tribe. A⇀,C⇀ are coefficient vectors, and their calculation formula is as follows:
(18)A⇀=2a⇀·r⇀1−a⇀.(19)C⇀=2·r⇀2.Among them,a⇀ is the convergence factor. With the number of iterations decreasing linearly from 2 to 0, r⇀1,r⇀2 express the random number between 0 and 1.The location update of grey wolf is shown in Figure3.Figure 3
Schematic diagram of grey wolf location update.When grey wolf catches prey, alpha (α) wolf at the first level leads other wolves to surround the prey. In the grey wolf population, alpha (α) wolf, beta (β) wolf, and delta (δ) wolf are the closest and most perceptive individuals. The location update of other grey wolf individuals is determined by the location of these three wolves. The schematic diagram and mathematical model of individual grey wolf searching for prey are as follows.
(20)D⇀i,αt=D⇀1·X⇀αt−X⇀it,(21)D⇀i,βt=D⇀2·X⇀βt−X⇀it,(22)D⇀i,δt=D⇀3·X⇀δt−X⇀it,(23)X⇀i,αt=X⇀α−A⇀·D⇀i,αt,(24)X⇀i,βt=X⇀β−A⇀·D⇀i,βt,(25)X⇀i,δt=X⇀δ−A⇀·D⇀i,δt,(26)X⇀it+1=X⇀i,αt+X⇀i.βt+X⇀i,δt3.Among them,D⇀i,αt represents the distance between the first generation of wolves and alpha (α) wolves, D⇀i,βt represents the distance between the second generation of grey wolves and beta (β) wolves, D⇀i,δt represents the distance between the second generation of wolves and delta (δ) wolves, the last three formulas define the step size and direction of omega (ω) wolves moving to the optimal solution alpha (α) wolves, beta (β) wolves, and delta (δ) wolves, and the last expression represents the new generation of grey wolf individuals after the location update.Pseudocode 1: Pseudocode of the grey wolf algorithm.
Initialize the grey wolf populationXii=1,2,⋯,nInitialize a, A and CCalculate the fitness of each search agentXα=the best search agentXβ=the second-best search agentXδ=the third best search agentwhile (t < Max number of iterations)for each search agentUpdate the position of the current search agent by equationend forUpdate a, A, and CCalculate the fitness of all search agentsUpdateXα,Xβ,and Xδt=t+1end whilereturnXα


### 4.2. Cross-Mutation Grey Wolf Algorithm

Two different grey wolf individual vectors are randomly selected to subtract and generate the difference vector. After the difference vector is given the weight, it is added to the third randomly selected grey wolf individual vector to generate the variation vector. This operation is called variation. The variation vector and the target vector are mixed to generate the test vector, which is called crossover. With the renewal of wolves, new individuals of grey wolves are produced, which make up for the problem of population limit in the optimization process.


#### 4.2.1. Mutation Operation

The most basic variation component is the difference vector of the parent generation. Each vector pair includes two different individualsxr1t,xr2t in the parent generation (generation t). And the difference vector is defined as
(27)Dr1,2=xr1t−xr2t,wherer1and r2 represent index numbers of two different individuals in the population. The variation vector is generated by adding the difference vector to another randomly selected individual vector. For each target vector xit, the mutation operation is that
(28)vit+1=xr3t+F∗xr1t−xr2t.r1,r2,r3∈2,⋯NP are different integers, and r1,r2,r3 are different from the current target vector index i, so the population size NP≥4. F is the scale factor, and the value range of F is 0,2 to control differential vector scaling.


#### 4.2.2. Cross-Operation

For the target vector individualsxit in the population, cross-operation vit+1 will be carried out with the variation vector to generate the experimental individuals uit+1. In order to ensure the evolution of an individual xit, first by random selection, at least one contributor vit+1 can be selected randomly in uit+1. For others, a cross-probability factor CR can be used to determine which bit is contributed by. The equation of cross-operation is
(29)uit+1=vit+1,randj≤CRorj=randni,xit,randnj>CRorj≠randni,whererandnj∈0,1 is the random number with uniform distribution. j represents the j-th variable (gene).CR is the cross-probability constant, and its value range is 0,1. The size is determined in advance. randni∈1,2,⋯,D is the index of randomly selected dimension variables to ensure that the test vector has at least dimension variables contributed by the variation vector. Otherwise, the test vector may be the same as the target vector and cannot generate new individuals. It can be seen from the formula that the larger the contribution CR, the more contribution vit+1 to uit+1. When CR=1, uit+1=vit+1. This is conducive to local search and accelerating convergence rate. The smaller the CR, the more contribution xit to uit+1. When CR=0, uit+1=xit, which is conducive to maintaining the diversity of the population and global search.


### 4.3. Grey Wolf Genetic Algorithm Flow

(1)
InitializationGm is the maximum number of iterations. nn is the number of nodes, N=5∗nn+1 is the dimension of the problem, NP is the population size, G=1 is the initial algebra, CR is hybrid parameters, u is the control moment, w is the angular velocity, and j is the moment of inertia.
(2)
Initialize the three leadersCalculate the fitness function of each grey wolf, screen out the smallest three fitness values, and use their corresponding individuals as the corresponding alpha wolf, beta wolf, and delta wolf.(3)
Reinitialize particlesUpdate the lead wolf position:(30)D⇀i,αt=D⇀1·X⇀αt−X⇀it,(31)D⇀i,βt=D⇀2·X⇀βt−X⇀it,(32)D⇀i,δt=D⇀3·X⇀δt−X⇀it,(33)X⇀i,αt=X⇀α−A⇀·D⇀i,αt,(34)X⇀i,βt=X⇀β−A⇀·D⇀i,βt,(35)X⇀i,δt=X⇀δ−A⇀·D⇀i,δt,(36)X⇀it+1=X⇀i,αt+X⇀i.βt+X⇀i,δt3.(4)
Cross(37)uit+1=vit+1,randj≤CRorj=randni,xit,randnj>CRorj≠randni.(5)
Mutation(38)Dr1,2=xr1t−xr2t.(39)vit+1=xr3t+F∗xr1t−xr2t.(6)
CompetitionThe particles are brought into the constraint function to calculate their fitness function values, and the minimum three are filtered out, which is the minimum value of this generation.The constraints to be satisfied are(40)Jw.=u−w×Jw,q=12Mq,ui≤u≤uj,wi≤w≤wj,qTHq≤0,qTKq≥0.Iteratively update and calculate the minimum value of each generation. Then, bring the minimum value of each generation into the fitness function to calculate the minimum value, which is the required value.(7)
EndPseudocode 2: Pseudocode of the cross-mutation grey wolf algorithm.
Initialize the grey wolf populationXii=1,2,⋯,nInitialize a, A and CCalculate the fitness of each search agentpopulation[t] = initialize Population (population Size)evaluate Population(population[t]Xα=the best search agentXβ=the second-best search agentXδ=the third best search agentwhile (t < Max number of iterations)for each search agentUpdate the position of the current search agent by equationend forpopulation[t+1] = crossover(population[t])population[t+1] = mutate(population[t+1])evaluate Population(population[t])Update a, A, and CCalculate the fitness of all search agentsUpdateXα,Xβ,and Xδt=t+1end whilereturnXα


## 4.1. Improved Cross-Mutation Grey Wolf Algorithm

After comparing the optimization results of different algorithms, how to avoid falling into the local optimization and get the optimization results faster is only the main problem concerned. In this paper, a new grey wolf algorithm called the grey wolf genetic algorithm is proposed. This strategy is enlightened from two aspects of genetic mutation and grey wolf algorithm. With the updating of the wolf position, the individuals in the wolf group continue to mutate, so that the population keeps updating and the problem of getting into the local optimum can be solved effectively. Meanwhile, the high efficiency of the grey wolf algorithm makes the time of algorithm implementation controlled, and the improved grey wolf algorithm flow is shown in Figure1.Figure 1
Algorithm flow diagram.


### 4.1.1. Grey Wolf Algorithm

According to the previous idea, the common grey wolf is described as follows: at the top of the grey wolf hierarchy is the supreme leader of the grey wolf tribe, known as alphaα wolf, who is responsible for overall planning and decisions on hunting. The second layer of the grey wolf’s social hierarchy is the second leader of the wolf pack, which is called beta β wolf. It mainly plays an assisting role in the overall planning and decision-making of alpha α wolf and is second only to alpha α wolf in the domination of the wolf population. At the third level of the grey wolf social hierarchy is delta δ wolf, also known as scouts, who are responsible for reconnaissance, watchkeeping, and guarding. The fourth level of the grey wolf social rank is the lowest in the grey wolf tribe, known as the omega (ω) wolf, which is essential for balancing the internal relations of the grey wolf tribe. The most suitable solution is considered to be the alpha (α) wolf. Therefore, the second and third best solutions were named beta β and omega ω wolves. The remaining candidate solutions were assumed to be omega wolves. And the three groups instructed the other wolves (omega) to search for the target. During the optimization process, the wolves update the positions of alpha, beta, delta, and omega; the grey wolf social hierarchy in the grey wolf algorithm is shown in Figure 2.Figure 2
Grey wolf social hierarchy.When the grey wolf population surrounds its prey, its predation behavior is defined as follows:(16)D⇀=C⇀·X⇀pt−X⇀t,(17)X⇀t+1=X⇀pt−A⇀·D⇀.Among them,D⇀ represents the distance between the individual wolf and the prey. The second formula is the updated formula of the individual wolf’s position. X⇀pt is the position vector of the prey in generation t. X⇀t is the position vector of the individuals in the t generation wolf tribe. A⇀,C⇀ are coefficient vectors, and their calculation formula is as follows:
(18)A⇀=2a⇀·r⇀1−a⇀.(19)C⇀=2·r⇀2.Among them,a⇀ is the convergence factor. With the number of iterations decreasing linearly from 2 to 0, r⇀1,r⇀2 express the random number between 0 and 1.The location update of grey wolf is shown in Figure3.Figure 3
Schematic diagram of grey wolf location update.When grey wolf catches prey, alpha (α) wolf at the first level leads other wolves to surround the prey. In the grey wolf population, alpha (α) wolf, beta (β) wolf, and delta (δ) wolf are the closest and most perceptive individuals. The location update of other grey wolf individuals is determined by the location of these three wolves. The schematic diagram and mathematical model of individual grey wolf searching for prey are as follows.
(20)D⇀i,αt=D⇀1·X⇀αt−X⇀it,(21)D⇀i,βt=D⇀2·X⇀βt−X⇀it,(22)D⇀i,δt=D⇀3·X⇀δt−X⇀it,(23)X⇀i,αt=X⇀α−A⇀·D⇀i,αt,(24)X⇀i,βt=X⇀β−A⇀·D⇀i,βt,(25)X⇀i,δt=X⇀δ−A⇀·D⇀i,δt,(26)X⇀it+1=X⇀i,αt+X⇀i.βt+X⇀i,δt3.Among them,D⇀i,αt represents the distance between the first generation of wolves and alpha (α) wolves, D⇀i,βt represents the distance between the second generation of grey wolves and beta (β) wolves, D⇀i,δt represents the distance between the second generation of wolves and delta (δ) wolves, the last three formulas define the step size and direction of omega (ω) wolves moving to the optimal solution alpha (α) wolves, beta (β) wolves, and delta (δ) wolves, and the last expression represents the new generation of grey wolf individuals after the location update.Pseudocode 1: Pseudocode of the grey wolf algorithm.
Initialize the grey wolf populationXii=1,2,⋯,nInitialize a, A and CCalculate the fitness of each search agentXα=the best search agentXβ=the second-best search agentXδ=the third best search agentwhile (t < Max number of iterations)for each search agentUpdate the position of the current search agent by equationend forUpdate a, A, and CCalculate the fitness of all search agentsUpdateXα,Xβ,and Xδt=t+1end whilereturnXα


## 4.1.1. Grey Wolf Algorithm

According to the previous idea, the common grey wolf is described as follows: at the top of the grey wolf hierarchy is the supreme leader of the grey wolf tribe, known as alphaα wolf, who is responsible for overall planning and decisions on hunting. The second layer of the grey wolf’s social hierarchy is the second leader of the wolf pack, which is called beta β wolf. It mainly plays an assisting role in the overall planning and decision-making of alpha α wolf and is second only to alpha α wolf in the domination of the wolf population. At the third level of the grey wolf social hierarchy is delta δ wolf, also known as scouts, who are responsible for reconnaissance, watchkeeping, and guarding. The fourth level of the grey wolf social rank is the lowest in the grey wolf tribe, known as the omega (ω) wolf, which is essential for balancing the internal relations of the grey wolf tribe. The most suitable solution is considered to be the alpha (α) wolf. Therefore, the second and third best solutions were named beta β and omega ω wolves. The remaining candidate solutions were assumed to be omega wolves. And the three groups instructed the other wolves (omega) to search for the target. During the optimization process, the wolves update the positions of alpha, beta, delta, and omega; the grey wolf social hierarchy in the grey wolf algorithm is shown in Figure 2.Figure 2
Grey wolf social hierarchy.When the grey wolf population surrounds its prey, its predation behavior is defined as follows:(16)D⇀=C⇀·X⇀pt−X⇀t,(17)X⇀t+1=X⇀pt−A⇀·D⇀.Among them,D⇀ represents the distance between the individual wolf and the prey. The second formula is the updated formula of the individual wolf’s position. X⇀pt is the position vector of the prey in generation t. X⇀t is the position vector of the individuals in the t generation wolf tribe. A⇀,C⇀ are coefficient vectors, and their calculation formula is as follows:
(18)A⇀=2a⇀·r⇀1−a⇀.(19)C⇀=2·r⇀2.Among them,a⇀ is the convergence factor. With the number of iterations decreasing linearly from 2 to 0, r⇀1,r⇀2 express the random number between 0 and 1.The location update of grey wolf is shown in Figure3.Figure 3
Schematic diagram of grey wolf location update.When grey wolf catches prey, alpha (α) wolf at the first level leads other wolves to surround the prey. In the grey wolf population, alpha (α) wolf, beta (β) wolf, and delta (δ) wolf are the closest and most perceptive individuals. The location update of other grey wolf individuals is determined by the location of these three wolves. The schematic diagram and mathematical model of individual grey wolf searching for prey are as follows.
(20)D⇀i,αt=D⇀1·X⇀αt−X⇀it,(21)D⇀i,βt=D⇀2·X⇀βt−X⇀it,(22)D⇀i,δt=D⇀3·X⇀δt−X⇀it,(23)X⇀i,αt=X⇀α−A⇀·D⇀i,αt,(24)X⇀i,βt=X⇀β−A⇀·D⇀i,βt,(25)X⇀i,δt=X⇀δ−A⇀·D⇀i,δt,(26)X⇀it+1=X⇀i,αt+X⇀i.βt+X⇀i,δt3.Among them,D⇀i,αt represents the distance between the first generation of wolves and alpha (α) wolves, D⇀i,βt represents the distance between the second generation of grey wolves and beta (β) wolves, D⇀i,δt represents the distance between the second generation of wolves and delta (δ) wolves, the last three formulas define the step size and direction of omega (ω) wolves moving to the optimal solution alpha (α) wolves, beta (β) wolves, and delta (δ) wolves, and the last expression represents the new generation of grey wolf individuals after the location update.Pseudocode 1: Pseudocode of the grey wolf algorithm.
Initialize the grey wolf populationXii=1,2,⋯,nInitialize a, A and CCalculate the fitness of each search agentXα=the best search agentXβ=the second-best search agentXδ=the third best search agentwhile (t < Max number of iterations)for each search agentUpdate the position of the current search agent by equationend forUpdate a, A, and CCalculate the fitness of all search agentsUpdateXα,Xβ,and Xδt=t+1end whilereturnXα


## 4.2. Cross-Mutation Grey Wolf Algorithm

Two different grey wolf individual vectors are randomly selected to subtract and generate the difference vector. After the difference vector is given the weight, it is added to the third randomly selected grey wolf individual vector to generate the variation vector. This operation is called variation. The variation vector and the target vector are mixed to generate the test vector, which is called crossover. With the renewal of wolves, new individuals of grey wolves are produced, which make up for the problem of population limit in the optimization process.


### 4.2.1. Mutation Operation

The most basic variation component is the difference vector of the parent generation. Each vector pair includes two different individualsxr1t,xr2t in the parent generation (generation t). And the difference vector is defined as
(27)Dr1,2=xr1t−xr2t,wherer1and r2 represent index numbers of two different individuals in the population. The variation vector is generated by adding the difference vector to another randomly selected individual vector. For each target vector xit, the mutation operation is that
(28)vit+1=xr3t+F∗xr1t−xr2t.r1,r2,r3∈2,⋯NP are different integers, and r1,r2,r3 are different from the current target vector index i, so the population size NP≥4. F is the scale factor, and the value range of F is 0,2 to control differential vector scaling.


### 4.2.2. Cross-Operation

For the target vector individualsxit in the population, cross-operation vit+1 will be carried out with the variation vector to generate the experimental individuals uit+1. In order to ensure the evolution of an individual xit, first by random selection, at least one contributor vit+1 can be selected randomly in uit+1. For others, a cross-probability factor CR can be used to determine which bit is contributed by. The equation of cross-operation is
(29)uit+1=vit+1,randj≤CRorj=randni,xit,randnj>CRorj≠randni,whererandnj∈0,1 is the random number with uniform distribution. j represents the j-th variable (gene).CR is the cross-probability constant, and its value range is 0,1. The size is determined in advance. randni∈1,2,⋯,D is the index of randomly selected dimension variables to ensure that the test vector has at least dimension variables contributed by the variation vector. Otherwise, the test vector may be the same as the target vector and cannot generate new individuals. It can be seen from the formula that the larger the contribution CR, the more contribution vit+1 to uit+1. When CR=1, uit+1=vit+1. This is conducive to local search and accelerating convergence rate. The smaller the CR, the more contribution xit to uit+1. When CR=0, uit+1=xit, which is conducive to maintaining the diversity of the population and global search.


## 4.2.1. Mutation Operation

The most basic variation component is the difference vector of the parent generation. Each vector pair includes two different individualsxr1t,xr2t in the parent generation (generation t). And the difference vector is defined as
(27)Dr1,2=xr1t−xr2t,wherer1and r2 represent index numbers of two different individuals in the population. The variation vector is generated by adding the difference vector to another randomly selected individual vector. For each target vector xit, the mutation operation is that
(28)vit+1=xr3t+F∗xr1t−xr2t.r1,r2,r3∈2,⋯NP are different integers, and r1,r2,r3 are different from the current target vector index i, so the population size NP≥4. F is the scale factor, and the value range of F is 0,2 to control differential vector scaling.


## 4.2.2. Cross-Operation

For the target vector individualsxit in the population, cross-operation vit+1 will be carried out with the variation vector to generate the experimental individuals uit+1. In order to ensure the evolution of an individual xit, first by random selection, at least one contributor vit+1 can be selected randomly in uit+1. For others, a cross-probability factor CR can be used to determine which bit is contributed by. The equation of cross-operation is
(29)uit+1=vit+1,randj≤CRorj=randni,xit,randnj>CRorj≠randni,whererandnj∈0,1 is the random number with uniform distribution. j represents the j-th variable (gene).CR is the cross-probability constant, and its value range is 0,1. The size is determined in advance. randni∈1,2,⋯,D is the index of randomly selected dimension variables to ensure that the test vector has at least dimension variables contributed by the variation vector. Otherwise, the test vector may be the same as the target vector and cannot generate new individuals. It can be seen from the formula that the larger the contribution CR, the more contribution vit+1 to uit+1. When CR=1, uit+1=vit+1. This is conducive to local search and accelerating convergence rate. The smaller the CR, the more contribution xit to uit+1. When CR=0, uit+1=xit, which is conducive to maintaining the diversity of the population and global search.


## 4.3. Grey Wolf Genetic Algorithm Flow

(1)
InitializationGm is the maximum number of iterations. nn is the number of nodes, N=5∗nn+1 is the dimension of the problem, NP is the population size, G=1 is the initial algebra, CR is hybrid parameters, u is the control moment, w is the angular velocity, and j is the moment of inertia.
(2)
Initialize the three leadersCalculate the fitness function of each grey wolf, screen out the smallest three fitness values, and use their corresponding individuals as the corresponding alpha wolf, beta wolf, and delta wolf.(3)
Reinitialize particlesUpdate the lead wolf position:(30)D⇀i,αt=D⇀1·X⇀αt−X⇀it,(31)D⇀i,βt=D⇀2·X⇀βt−X⇀it,(32)D⇀i,δt=D⇀3·X⇀δt−X⇀it,(33)X⇀i,αt=X⇀α−A⇀·D⇀i,αt,(34)X⇀i,βt=X⇀β−A⇀·D⇀i,βt,(35)X⇀i,δt=X⇀δ−A⇀·D⇀i,δt,(36)X⇀it+1=X⇀i,αt+X⇀i.βt+X⇀i,δt3.(4)
Cross(37)uit+1=vit+1,randj≤CRorj=randni,xit,randnj>CRorj≠randni.(5)
Mutation(38)Dr1,2=xr1t−xr2t.(39)vit+1=xr3t+F∗xr1t−xr2t.(6)
CompetitionThe particles are brought into the constraint function to calculate their fitness function values, and the minimum three are filtered out, which is the minimum value of this generation.The constraints to be satisfied are(40)Jw.=u−w×Jw,q=12Mq,ui≤u≤uj,wi≤w≤wj,qTHq≤0,qTKq≥0.Iteratively update and calculate the minimum value of each generation. Then, bring the minimum value of each generation into the fitness function to calculate the minimum value, which is the required value.(7)
EndPseudocode 2: Pseudocode of the cross-mutation grey wolf algorithm.
Initialize the grey wolf populationXii=1,2,⋯,nInitialize a, A and CCalculate the fitness of each search agentpopulation[t] = initialize Population (population Size)evaluate Population(population[t]Xα=the best search agentXβ=the second-best search agentXδ=the third best search agentwhile (t < Max number of iterations)for each search agentUpdate the position of the current search agent by equationend forpopulation[t+1] = crossover(population[t])population[t+1] = mutate(population[t+1])evaluate Population(population[t])Update a, A, and CCalculate the fitness of all search agentsUpdateXα,Xβ,and Xδt=t+1end whilereturnXα


## 5. Simulation and Analysis

Basing on the research of previous scientists, we found some questions. In the first segment of this chapter, we will introduce the cross-mutation grey wolf optimizer algorithm and conduct simulation analysis. Later in this section, we will show the simulation results to get the advantages of this algorithm. Finally, we will compare the cross-mutation grey wolf Optimizer with the traditional GWO to verify the validity of the cross-mutation grey wolf optimizer.


### 5.1. Parameter Initialization

In the process of spacecraft, the vector components of four strong light celestial bodies in the inertial system arerA, rB, rC, and rD and position vectors are r1 and r2. The initial time of spacecraft is at position r. The minimum angle between r and r1 is required to be α1, and the minimum angle between r and r2 is required to be α2. The included angle of self-instrument is β. The initial attitude and angular velocity of the spacecraft are, respectively, q and ω, as shown in Table 2.Table 2
Numerical simulation of spacecraft operating parameters.
VariableSimulation valuew[0, 0, 0] rad/sq[-0.383, 0, 0.926, 0]Jdiag(1500,1500,1500)kg·m2r10.56,0.31,−0.18Tr20.65,−0.19,0.54Tα115°α230°β15°rA0,0.2,0.5rB0.8,0.3,−0.4rC0.5,0.3,−0.1rD0.8,0.6,0.2r0,0,1Tγw0.512γu0.114Besides, after data deduction, the condition of forbidden constraint is that the range ofω is from -0.512 to 0.512 and the range of u is from -0.114 to 0.114.


### 5.2. Drawing Results and Result Analysis

The following is the simulation results of the cross-mutation grey wolf optimizer:Using the cross-mutation grey wolf optimizer algorithm to solve, when△T is 1 s, the total time of the maneuvering process is 210 s. The number of path node is 210. Figure 4 shows the attitude maneuver path of spacecraft in the spherical coordinate system of spacecraft, in which two red stars are the starting point and the ending point. We can see that in the process of movement, spacecraft is out of strong light and avoid being directly at the beam of strong light celestial bodies. So, the optimization result is wonderful; the trajectory of the infrared telescope vector in the celestial coordinate system of the cross-mutation grey wolf optimizer is shown in Figure 4, where blue dots indicate obstacles.Figure 4
Trajectory of the infrared telescope vector in the celestial coordinate system of the cross-mutation grey wolf optimizer.Figures5 and 6 are the change of w, q, and u of the satellite in the application of the cross-mutation GWO algorithm.Figure 5
Angular velocity time history based on the cross-mutation grey wolf optimizer.Figure 6
Control torque time history based on the cross-mutation grey wolf optimizer.According to the boundary constraint, when the spacecraft is working,w should be within a constraint. The attitude model can predict the change of w. Experimental data confirm this, and we find that the range of change is small and is in the boundary constraint, as shown in Figure 5.Control torque is also within the constraint range and finally tends to be stable, which is a good thing for the spacecraft and is conducive to stable execution, as shown in Figure6.We found that both angular velocity and control torque are within their range; that is to say, boundary constraints are satisfied. It is a good proof of feasibility. Secondly, the curves of angular velocity and moment are also very good, which conform to the normal working parameters of spacecraft and will not cause failure to spacecraft.By boundary constraint and attitude constraint, the above two results are obtained, which strongly supports the key assumption that the cross-mutation grey wolf optimizer can optimize the flight stability of spacecraft.The change of quaternion is very smooth, and the range is large, so it is easy to face all kinds of situations. And it is helpful to avoid singular value, as shown in Figure7.Figure 7
Attitude quaternion time history based on the cross-mutation grey wolf optimizer.The calculation time curve of the cross-mutation grey wolf optimizer is shown in Figure8.Figure 8
The computation time curve of the cross-mutation grey wolf optimizer.The cross-mutation grey wolf optimizer algorithm solves some shortcomings of the previous algorithm. For example, the previous algorithm is unable to cope with wide range of parameters and the influence of nonlinear factors on the system. Besides, the cross-mutation grey wolf optimizer improved computing and motion speed. The improved curve is very stable, and it can obviously improve the steady and dynamic performance of the system and has better control effect.The control effect of the cross-mutation grey wolf optimizer algorithm is verified by numerical simulation. From the result chart, the effect is very useful, and it can effectively control the spacecraft to work. At the same time, the cross-mutation operation is good for jumping out of the local optimal value. In trajectory planning, the optimal route can be planned faster and better.It is not difficult to verify the optimization stability of the cross-mutation grey wolf optimizer algorithm through pictures. The curves are smooth and the straight-line transition is very smooth, which proves the stability of the cross-mutation grey wolf optimizer algorithm.Besides, we choose to simulate spacecraft flights many times to be samples to draw the pictures, which can predict most situations. By doing that, we can adjust parameters to ensure the comprehensiveness of information collection. This will improve the stability of the spacecraft.


### 5.3. The Comparison of the Cross-Mutation Grey Wolf Optimizer and Traditional GWO

In order to verify the effectiveness of the cross-mutation grey wolf optimizer algorithm, we choose the traditional GWO algorithm to be applied to the model and compare it with the cross-mutation algorithm.Under the same simulation conditions, we also make experiment with the traditional GWO algorithm which takes 265 s. With the increase of the number of experiments, the traditional GWO is easy to fall to the local optimum. And it can only avoid most celestial bodies. Trajectory cannot avoid all objects.Because of the randomness of the algorithm, we have carried out several experiments and recorded the time, as shown in Table3.Table 3
Algorithm simulation experiment data comparison.
Times maneuver time algorithmGWOCross-mutation grey wolf optimizer1258.38212.202261.24209.353266.15214.344256.36213.225260.53214.976255.87216.757259.61208.678261.98209.669264.62207.4810259.41209.11Average time260.415211.575In addition, we conducted 100 more experiments. We used the GWO and cross-mutation grey wolf optimizer algorithm to calculate the satellite’s maneuver time. We can draw the conclusion that the maneuvering time of GWO with cross-variation is about 210 s. The same GWO maneuver time is between 260 s and 270 s, as shown in Table4.Table 4
Evaluating indicator.
AlgorithmBoundary constraintsForbidden constraintComputing timeMovement timeThe change ofw and uTraditional GWO×√105.2 s265 s×Cross-mutation grey wolf optimizer√√95 s210 s√


#### 5.3.1. Interpretations

(1)
According to the image of angular velocity of the traditional GWO, we find that the maximum value ofw exceeds the boundary constraint 5.12; that is, boundary constraints are not satisfied. Spacecraft is prone to failure. And the cross-mutation grey wolf optimizer solves this problem very well, as shown in Figure 4(2)
In terms of time, no matter it is computing time or maneuvering time, the cross-mutation grey wolf optimizer is shorter than the traditional GWO showing that the cross-mutation grey wolf optimizer algorithm is more efficient. It is also very stable in execution and not prone to errors(3)
Another possibility is that the traditional GWO algorithm produces a local optimal value in operation. Therefore, it is not possible to quickly iterate the head wolf and filter out the minimum value, thus increasing the calculation time. However, the cross-mutation grey wolf optimizer solves this problem well(4)
For the change ofω and u, the change of the cross-mutation grey wolf optimizer is smoother than that of common GWO range which is also more reasonable, so it has better application effectIn addition, we also explore the local optimization of the two algorithms, because the cross-mutation grey wolf optimizer can effectively improve the local optimum, so we did 100 experiments and recorded the times of falling into the local optimum, as shown in Figure9.Figure 9
Discussion on the situation of the algorithm falling into the local optimum.It can be seen from the chart that cross and mutation can significantly reduce the number of times falling into the local optimum, which is effective. In practical application, it also performs well.The following figures are the simulation diagram under the grey wolf optimizer algorithm.The time history analysis of angular velocity based on the grey wolf optimization algorithm is shown in Figure10.Figure 10
Angular velocity time history based on the grey wolf optimizer.The control torque time history based on the grey wolf optimizer is shown in Figure11.Figure 11
Control torque time history based on the grey wolf optimizer.The attitude quaternion time history based on the grey wolf optimizer is shown in Figure12.Figure 12
Attitude quaternion time history based on the grey wolf optimizer.The computation time curve of the grey wolf optimizer is shown in Figure13.Figure 13
The computation time curve of the grey wolf optimizer.The trajectory of the infrared telescope vector in the celestial coordinate system of the grey wolf optimizer is shown in Figure14.Figure 14
Trajectory of the infrared telescope vector in the celestial coordinate system of the grey wolf optimizer.


## 5.1. Parameter Initialization

In the process of spacecraft, the vector components of four strong light celestial bodies in the inertial system arerA, rB, rC, and rD and position vectors are r1 and r2. The initial time of spacecraft is at position r. The minimum angle between r and r1 is required to be α1, and the minimum angle between r and r2 is required to be α2. The included angle of self-instrument is β. The initial attitude and angular velocity of the spacecraft are, respectively, q and ω, as shown in Table 2.Table 2
Numerical simulation of spacecraft operating parameters.
VariableSimulation valuew[0, 0, 0] rad/sq[-0.383, 0, 0.926, 0]Jdiag(1500,1500,1500)kg·m2r10.56,0.31,−0.18Tr20.65,−0.19,0.54Tα115°α230°β15°rA0,0.2,0.5rB0.8,0.3,−0.4rC0.5,0.3,−0.1rD0.8,0.6,0.2r0,0,1Tγw0.512γu0.114Besides, after data deduction, the condition of forbidden constraint is that the range ofω is from -0.512 to 0.512 and the range of u is from -0.114 to 0.114.


## 5.2. Drawing Results and Result Analysis

The following is the simulation results of the cross-mutation grey wolf optimizer:Using the cross-mutation grey wolf optimizer algorithm to solve, when△T is 1 s, the total time of the maneuvering process is 210 s. The number of path node is 210. Figure 4 shows the attitude maneuver path of spacecraft in the spherical coordinate system of spacecraft, in which two red stars are the starting point and the ending point. We can see that in the process of movement, spacecraft is out of strong light and avoid being directly at the beam of strong light celestial bodies. So, the optimization result is wonderful; the trajectory of the infrared telescope vector in the celestial coordinate system of the cross-mutation grey wolf optimizer is shown in Figure 4, where blue dots indicate obstacles.Figure 4
Trajectory of the infrared telescope vector in the celestial coordinate system of the cross-mutation grey wolf optimizer.Figures5 and 6 are the change of w, q, and u of the satellite in the application of the cross-mutation GWO algorithm.Figure 5
Angular velocity time history based on the cross-mutation grey wolf optimizer.Figure 6
Control torque time history based on the cross-mutation grey wolf optimizer.According to the boundary constraint, when the spacecraft is working,w should be within a constraint. The attitude model can predict the change of w. Experimental data confirm this, and we find that the range of change is small and is in the boundary constraint, as shown in Figure 5.Control torque is also within the constraint range and finally tends to be stable, which is a good thing for the spacecraft and is conducive to stable execution, as shown in Figure6.We found that both angular velocity and control torque are within their range; that is to say, boundary constraints are satisfied. It is a good proof of feasibility. Secondly, the curves of angular velocity and moment are also very good, which conform to the normal working parameters of spacecraft and will not cause failure to spacecraft.By boundary constraint and attitude constraint, the above two results are obtained, which strongly supports the key assumption that the cross-mutation grey wolf optimizer can optimize the flight stability of spacecraft.The change of quaternion is very smooth, and the range is large, so it is easy to face all kinds of situations. And it is helpful to avoid singular value, as shown in Figure7.Figure 7
Attitude quaternion time history based on the cross-mutation grey wolf optimizer.The calculation time curve of the cross-mutation grey wolf optimizer is shown in Figure8.Figure 8
The computation time curve of the cross-mutation grey wolf optimizer.The cross-mutation grey wolf optimizer algorithm solves some shortcomings of the previous algorithm. For example, the previous algorithm is unable to cope with wide range of parameters and the influence of nonlinear factors on the system. Besides, the cross-mutation grey wolf optimizer improved computing and motion speed. The improved curve is very stable, and it can obviously improve the steady and dynamic performance of the system and has better control effect.The control effect of the cross-mutation grey wolf optimizer algorithm is verified by numerical simulation. From the result chart, the effect is very useful, and it can effectively control the spacecraft to work. At the same time, the cross-mutation operation is good for jumping out of the local optimal value. In trajectory planning, the optimal route can be planned faster and better.It is not difficult to verify the optimization stability of the cross-mutation grey wolf optimizer algorithm through pictures. The curves are smooth and the straight-line transition is very smooth, which proves the stability of the cross-mutation grey wolf optimizer algorithm.Besides, we choose to simulate spacecraft flights many times to be samples to draw the pictures, which can predict most situations. By doing that, we can adjust parameters to ensure the comprehensiveness of information collection. This will improve the stability of the spacecraft.


## 5.3. The Comparison of the Cross-Mutation Grey Wolf Optimizer and Traditional GWO

In order to verify the effectiveness of the cross-mutation grey wolf optimizer algorithm, we choose the traditional GWO algorithm to be applied to the model and compare it with the cross-mutation algorithm.Under the same simulation conditions, we also make experiment with the traditional GWO algorithm which takes 265 s. With the increase of the number of experiments, the traditional GWO is easy to fall to the local optimum. And it can only avoid most celestial bodies. Trajectory cannot avoid all objects.Because of the randomness of the algorithm, we have carried out several experiments and recorded the time, as shown in Table3.Table 3
Algorithm simulation experiment data comparison.
Times maneuver time algorithmGWOCross-mutation grey wolf optimizer1258.38212.202261.24209.353266.15214.344256.36213.225260.53214.976255.87216.757259.61208.678261.98209.669264.62207.4810259.41209.11Average time260.415211.575In addition, we conducted 100 more experiments. We used the GWO and cross-mutation grey wolf optimizer algorithm to calculate the satellite’s maneuver time. We can draw the conclusion that the maneuvering time of GWO with cross-variation is about 210 s. The same GWO maneuver time is between 260 s and 270 s, as shown in Table4.Table 4
Evaluating indicator.
AlgorithmBoundary constraintsForbidden constraintComputing timeMovement timeThe change ofw and uTraditional GWO×√105.2 s265 s×Cross-mutation grey wolf optimizer√√95 s210 s√


### 5.3.1. Interpretations

(1)
According to the image of angular velocity of the traditional GWO, we find that the maximum value ofw exceeds the boundary constraint 5.12; that is, boundary constraints are not satisfied. Spacecraft is prone to failure. And the cross-mutation grey wolf optimizer solves this problem very well, as shown in Figure 4(2)
In terms of time, no matter it is computing time or maneuvering time, the cross-mutation grey wolf optimizer is shorter than the traditional GWO showing that the cross-mutation grey wolf optimizer algorithm is more efficient. It is also very stable in execution and not prone to errors(3)
Another possibility is that the traditional GWO algorithm produces a local optimal value in operation. Therefore, it is not possible to quickly iterate the head wolf and filter out the minimum value, thus increasing the calculation time. However, the cross-mutation grey wolf optimizer solves this problem well(4)
For the change ofω and u, the change of the cross-mutation grey wolf optimizer is smoother than that of common GWO range which is also more reasonable, so it has better application effectIn addition, we also explore the local optimization of the two algorithms, because the cross-mutation grey wolf optimizer can effectively improve the local optimum, so we did 100 experiments and recorded the times of falling into the local optimum, as shown in Figure9.Figure 9
Discussion on the situation of the algorithm falling into the local optimum.It can be seen from the chart that cross and mutation can significantly reduce the number of times falling into the local optimum, which is effective. In practical application, it also performs well.The following figures are the simulation diagram under the grey wolf optimizer algorithm.The time history analysis of angular velocity based on the grey wolf optimization algorithm is shown in Figure10.Figure 10
Angular velocity time history based on the grey wolf optimizer.The control torque time history based on the grey wolf optimizer is shown in Figure11.Figure 11
Control torque time history based on the grey wolf optimizer.The attitude quaternion time history based on the grey wolf optimizer is shown in Figure12.Figure 12
Attitude quaternion time history based on the grey wolf optimizer.The computation time curve of the grey wolf optimizer is shown in Figure13.Figure 13
The computation time curve of the grey wolf optimizer.The trajectory of the infrared telescope vector in the celestial coordinate system of the grey wolf optimizer is shown in Figure14.Figure 14
Trajectory of the infrared telescope vector in the celestial coordinate system of the grey wolf optimizer.


## 5.3.1. Interpretations

(1)
According to the image of angular velocity of the traditional GWO, we find that the maximum value ofw exceeds the boundary constraint 5.12; that is, boundary constraints are not satisfied. Spacecraft is prone to failure. And the cross-mutation grey wolf optimizer solves this problem very well, as shown in Figure 4(2)
In terms of time, no matter it is computing time or maneuvering time, the cross-mutation grey wolf optimizer is shorter than the traditional GWO showing that the cross-mutation grey wolf optimizer algorithm is more efficient. It is also very stable in execution and not prone to errors(3)
Another possibility is that the traditional GWO algorithm produces a local optimal value in operation. Therefore, it is not possible to quickly iterate the head wolf and filter out the minimum value, thus increasing the calculation time. However, the cross-mutation grey wolf optimizer solves this problem well(4)
For the change ofω and u, the change of the cross-mutation grey wolf optimizer is smoother than that of common GWO range which is also more reasonable, so it has better application effectIn addition, we also explore the local optimization of the two algorithms, because the cross-mutation grey wolf optimizer can effectively improve the local optimum, so we did 100 experiments and recorded the times of falling into the local optimum, as shown in Figure9.Figure 9
Discussion on the situation of the algorithm falling into the local optimum.It can be seen from the chart that cross and mutation can significantly reduce the number of times falling into the local optimum, which is effective. In practical application, it also performs well.The following figures are the simulation diagram under the grey wolf optimizer algorithm.The time history analysis of angular velocity based on the grey wolf optimization algorithm is shown in Figure10.Figure 10
Angular velocity time history based on the grey wolf optimizer.The control torque time history based on the grey wolf optimizer is shown in Figure11.Figure 11
Control torque time history based on the grey wolf optimizer.The attitude quaternion time history based on the grey wolf optimizer is shown in Figure12.Figure 12
Attitude quaternion time history based on the grey wolf optimizer.The computation time curve of the grey wolf optimizer is shown in Figure13.Figure 13
The computation time curve of the grey wolf optimizer.The trajectory of the infrared telescope vector in the celestial coordinate system of the grey wolf optimizer is shown in Figure14.Figure 14
Trajectory of the infrared telescope vector in the celestial coordinate system of the grey wolf optimizer.


## 6. Conclusions

This paper describes and analyzes the boundary and pointing constraints of the spacecraft trajectory planning problem based on time and path optimization. Firstly, a mathematical model based on boundary constraints and pointing constraints is established, and then, they are constructed as a finite function of path evaluation. In order to solve the path planning problem under the condition of satisfying the constraints, we propose a new attitude maneuver path planning method—cross-mutation grey wolf algorithm, which improves the diversity of the wolves’ offspring through cross-mutation and obtains the nearest optimal path based on iterative rules. Simulation results show that other forms of the GWO algorithm have good applicability in the process of spacecraft path planning, and the simulation results prove the feasibility of applying grey wolf to the spacecraft control model to solve the shortest path, although there are some problems such as unsatisfying constraints or taking too long to optimize. The CMGWO method performs well in simulation experiments, and it fits perfectly with the control model and not only satisfies the boundary constraints and direction constraints but also greatly shortens the attitude maneuver time and significantly reduces the times of falling into the local optimum. In addition, the angular velocity and control torque curves obtained by this method are smoother and gentler, which has a good application prospect in the future. In the future research, a single satellite may not be able to complete the task accurately and quickly in some cases, which indicates that the difficulty of solving the problem will continue to increase, and it will take more time to plan the attitude maneuver path. The CMGWO method also has the potential to be applied to the trajectory planning task of multiple satellites. Finally, the proposed numerical test is aimed at a specific situation. For some reasons, the simulation cannot be extended, but its universality is worthy of recognition



---

*Source: 8857584-2022-01-10.xml*