Genetic Programming Method for Satellite System Topology and Parameter Optimization

In this paper, a genetic programming method for satellite system design is proposed to simultaneously optimize the topology and parameters of a satellite system. Firstly, the representation of satellite system design is defined according to the tree structure. The genetic programming method is designed based on that representation. Secondly, according to the tree structure of different satellite schemes, different multiscale satellite models are established, in which various physical fields couple together. Then, an evaluation system is also proposed to test the performances of different satellite schemes. Finally, the application to the design of an earth observation satellite demonstrates the effectiveness of the proposed method.


Introduction
The optimization of satellite system design can improve the design quality of satellites and plays an important role in ensuring the satellite's advancement, economic efficiency, and reliability. Since the late 1990s, many researches have been implemented and presented. Fukunaga et al. [1] introduced artificial intelligence and machine learning related technologies and proposed a general heuristic algorithm for spacecraft optimization design. Jilla et al. [2] proposed a multidisciplinary design optimization method that combines the life cycle cost model with a sensitivity analysis tool and applied it to the conceptual design of a distributed satellite system. Dong et al. [3] proposed an optimization method based on double clustering and principal component analysis and applied it to the optimization design of an earth observation satellite. Richie et al. [4] introduced an efficient, structured approach for the optimization of satellite attitude and orbit control as well as power subsystem. Fazeley et al. [5] applied a multiobjective and multidisciplinary optimization method to the design of a bipropellant propulsion system of a spacecraft. In conclusion, the current researches mainly focus on three parts: (1) the optimization of the system parameters such as orbit, payload, mass, and power; (2) the selection of important schemes for the subsystems; and (3) the component level optimization of each subsystem. However, the system composition and structure optimization (i.e., the system topology optimization) are less considered. To solve this problem, two aspects of theoretical and technical basis are essential, namely, the system modeling method and optimization algorithm suitable for satellite system topology optimization.
In the aspect of the satellite system modeling method, text-based system engineering (TBSE) is mainly used in the traditional mode. In that mode, the system design team and the subsystem design team communicate and coordinate through various normative documents. The disadvantage of this mode is that it is difficult to reuse the analysis model of the system design and subsystem design [6,7]. Therefore, simplified models, empirical formulas, and statistical data are commonly used in system design [8,9]. Such a simplified model cannot reflect the changes of the scheme in the component scale, which may have a significant influence on the system design results. Due to the lack of couplings among different disciplines, the simplified model cannot fully reflect the impact of system composition and structure changes on the whole satellite [10]. The feasible solution derived from the simplified model does not guarantee the global optimal solution [9]. Aiming at the problem above, the model-based system engineering (MBSE) research has risen rapidly. In engineering application, Carvalho studied the application of MBSE on IRIS low-cost satellite [11]. NASA applied the MBSE method to the optimization design of the thermal control subsystem of the James Weber space telescope [12] and the MARS 2020 project [13]. As a modeling tool, Digital Satellite (DS) is a new digital achievement where MBSE is applied to satellite [14]. In short, the application of MBSE to the system design optimization of satellites is a crucial research direction to further improve the system design level. One of the key tasks is establishing a multiscale and multiphysical coupling model (hereinafter referred to as the MSMPC model) from the component model to the system model [15,16]. However, how to realize the rapid construction of a MSMPC model in the satellite system design optimization requires further research.
From the aspect of satellite system optimization algorithm, genetic algorithm [17], ant colony algorithm [18], simulated annealing algorithm [19], and other optimization methods are currently widely used in engineering to solve the problem of satellite system parameter optimization. In the algorithms above, the structure of the optimization object is invariant. Therefore, these algorithms are not suitable for the optimization of the composition and structure of the satellite system. On the contrary, by using a tree structure representation method, Genetic Programming (GP) proposed by Professor Koza [20] can express complex structured expressions. This feature makes it suitable for solving optimization problems with variable expression structures [21,22]. Sripramong and Toumazou [23] designed an automatic circuit design system for a CMOS amplifier by GP. Sun et al. [24] realized the automatic evolution of satellite attitude control law expression by GP. The applications of GP in simple engineering practices have shown its ability to solve the optimization problem with variable expression structure. However, as far as we know, there is currently no effective method to express a complex, large, and highly integrated satellite system into a tree structure suitable for GP.
The main contributions of this paper are summarized as follows. (1) By introducing the GP into the satellite system design, a GP expression method of the satellite system design scheme is established, which solves the problem in describing the composition and structure of the satellite system. (2) By combining the GP scheme with the Digital Satellite [14] tools, the problem of fast construction of MSMPC models in system design optimization is solved.
The remainder of this paper is organized as follows. Section 2 introduces the GP optimization framework of satellite system design and proposes the expression method of a satellite system design scheme suitable for GP. In Section 3, the method of establishing MSMPC satellite models is described. The calculation method of the objective function and constraints according to the satellite design scheme is presented in Section 4. The case study and analysis are given in Section 5. Finally, the conclusions are provided in Section 6.

GP Optimization Framework for Satellite
System Design 2.1. Problem Description of Satellite System Design Optimization. The system design scheme set S tree of a satellite is obtained by the dimensional decomposition method. The design variable of the optimization problem is denoted as S, which is an element of the set S tree . For a specified satellite scheme S, a MSMPC model can be uniquely determined, which is denoted as f X . The relationship between S and f ðXÞ can be expressed as S ! D S f ðXÞ, where D S represents the modeling method based on the Digital Satellite. X is the set where all the states of the satellite system described in f ðXÞ are included. The system state equation is expressed as _ X = f ðXÞ. The objective function and constraints of the satellite system design are expressed as JðXÞ and CðXÞ, respectively. The optimization problem of the satellite system design can be described as  Figure 1, which is divided into two modules: GP and the Digital Satellite. The optimization process starts with the GP module. The first step is to initialize the GP parameters. In the second step, the initial population of GP is generated, and the tree structure of the satellite individuals is obtained. In the third step, the tree structure is transformed into the satellite scheme description based on the component library. Then, the optimization process gets into the Digital Satellite module. In the fourth step, for each individual in the population, the MSMPC model is obtained, and the satellite evaluation condition set is designed from the condition library. In the fifth step, the batch processing method is used to complete the calculation in the verification simulation condition set. In the sixth step, the evaluation index of each scheme in the population is obtained by statistics of the simulation results of the condition set. In the seventh step, the fitness value is calculated according to the evaluation index and returned to the GP module. In the eighth step, individuals in the population are ranked according to their fitness. The ninth step is to judge whether the optimal fitness is better than the expected value. If so, stop the optimization to get the optimal result. If not, the process step in the tenth step, the next-generation population of design schemes is generated by the genetic operation of replication, crossover, and mutation. The individual tree structure of the next-generation population is iterated from the third step to the tenth step until the termination conditions are satisfied, and the ideal satellite design results are obtained.

2
International Journal of Aerospace Engineering In the following two subsections: two key steps of the flowchart will be described, namely, the description of satellite scheme S and the basic GP operations.

Description of Satellite System
Design Scheme Based on the Component Library. The satellite system design scheme S consists of environment, function, composition, structure, mode, procedure, and operation. In the system design of satellites, different satellite scheme S can be decomposed to obtain different tree structure representations. The detailed statements of those elements are presented as follows.
(1) System environment parameters contain the orbit elements and task start time (2) For the system function, the task type supported by satellite payload is considered In this paper, taking the system composition as the trunk, the tree structure description of the satellite system scheme is established, as shown in Figure 2. Environment and function parameters can be represented in the subsystem parameter nodes. In addition, the subsystem parameter node can also describe the technical route of the satellite subsystem. The technical route refers to the system structure, mode, program, and operation. The 2nd level nodes indicate subsystems, including attitude determination and control subsystem (ADCS), power subsystem (Power), propulsion subsystem (Propulsion), thermal control subsystem (Thermal), Communication subsystem (TTC), structure subsystem (Structure), and payload subsystem (Payload). In the 3rd level node, h is the orbital altitude, and Tc is the control period. According to the device category nodes in the third level and the device type name node in the fifth level, the components' model can be uniquely determined from the component library.

GP Operations of Satellite System Design Scheme.
According to the tree structure of the satellite system design scheme, the tree structure of individuals in the population of GP is established. As shown in Figure 2, the nodes of the parse tree are divided into two categories. One category is called "operator," which locates in the internal of the tree. The other category is called "terminator" located in the terminal of the tree. The definition of the operator set and terminator set is shown in Table 1.
The basic operations of GP include replication, crossover, and mutation. In this paper, according to the particularity of the satellite scheme tree structure, the crossover and  [14,25]. FitðXÞ is the fitness function.
3 International Journal of Aerospace Engineering mutation operations in the GP operation are adapted. In terms of crossover operation, for the "ParaInfo" nodes, only the parameters with the same name can cross and exchange their values. For other nodes, only two nodes whose superior nodes share the same name can cross. In terms of mutation operation, for the "NUM" and "Name" nodes, the mutation can reinitialize the number and type of components. For the "ParaInfo" node, the mutation can only change the parameter value, not the parameter type. For other nodes, the mutation can only reinitialize its child nodes, but it cannot change the properties of the nodes themselves. Figure 3 takes the "DeviceCategory" node as an example to show the process of the ðN + 1Þ th generation schemes is produced by crossing and mutating two schemes in the N th generation.

Multiscale and Multiphysical Coupling Models of Satellite
The satellite system model provides support for the quantitative evaluation of satellite schemes. In order to reflect the real impact of the system composition and structure on the optimization objectives, it is necessary to establish a MSMPC model [26].

Analysis of the Main Coupling
Relationship. The scale coupling relationship between celestial bodies and satellites is shown in Figure 4. The celestial scale coupling models are represented in the outer frame. The Earth has 4 coupling relationships as follows. (1) It is coupled with the ADCS subsystem through the resultant force/moment formed by gravity, aerodynamic force, and magnetic force. (2) It is coupled with the Thermal subsystem through infrared thermal radiation and reflection thermal radiation of the sun. (3) It is coupled with the Power subsystem through the shading of the sun by the earth shadow area. (4) It is coupled with the TTC subsystem through the effect of the earth radiation belt and ionosphere on the wireless signal. The coupling model of the sun includes the coupling of the resultant force/moment formed by gravity and light pressure with the ADCS subsystem, coupling with the Thermal subsystem by solar thermal radiation, and coupling with the Power subsystem through the photoelectric effect. The moon and other planets are coupled with the ADCS subsystem through gravitation. Ground stations, satellites, and other entities are coupled with the TTC subsystem through wireless signals. In other aspects, the damage from space debris is coupled with the structure subsystem. The cosmic ray and space noise are coupled with the TTC subsystem through the effect of the    Operator set "Satellite" The 1 st level nodes of the tree structure "Subsystem" The subsystem nodes in the 2 nd level "Parameters" The subsystem parameter branch nodes in the 2 nd level "Devices" The devices branch nodes in the 2 nd level "DeviceCategory" The device category nodes in the 3 rd level "DeviceType" The device type nodes in the 4 th level Terminator set "Name" The device type name nodes in the 5 th level "NUM" The device number nodes in the 5 th level "ParaInfo" The parameter information nodes in the 3 rd level 5 International Journal of Aerospace Engineering the torques of wheel, flexible attachment, liquid-filled tank, and flexible manipulator on rigid satellite body.
The main coupling is reflected in the coupling between the working state of the actuator, the power supply, and the thermal control state. For example, when the wheels run in the torque mode, the axle temperature affects the friction coefficient of the wheel bearings, which affects the output torque. The equation is established as follows [27].
where M M is the motor output torque, K M is the torque voltage ratio coefficient, V C is the input voltage, K f is the friction coefficient, C is the specific heat of the lubricating oil, ρ is the lubricating oil density, G is the amount of oil required for circulating lubricating oil, D is the momentum wheel bearing diameter, and ω w is the rotor speed. Considering the friction factor, M f 0 is the static friction moment of the bearing; the installation unit vector of the wheel in the satellite is r. The output torque of the wheels in the body coordinate is [27].
Considering all factors, the orbit and attitude dynamic equations are established as follows [27].    Figure 4: Coupling relationship in celestial-and satellite-scale. TM is short for telemetry. TC is the abbreviation for telecommand. T represents the temperature of subsystems. 6 International Journal of Aerospace Engineering

3.2.2.
Coupling Models in the Thermal Subsystem. The system state variable set of the Thermal subsystem is a subset of satellite system state X. It includes node temperature T i , infrared radiation angle coefficient of the earth ϕ 3i , infrared reflection angle coefficient of the earth ϕ 2i , solar radiation angle coefficient ϕ 1i, and other external heat calculationrelated parameters. The Thermal subsystem is coupled with the Power subsystem by P i which is the sum of the thermal power of all working components belonging to node i. The main coupling of the model is that ϕ 3i , ϕ 2i, and ϕ 1i are coupled with dynamics [28]. k is defined as k = R E /ðR E + hÞ, where R E is the average radius of the earth and h is the orbital altitude of the satellite. β i is the angle between the normal direction of the external surface of node i and the earth satellite line. When β i satisfies 0 ≤ β i ≤ cos −1 k, the infrared radiation angle coefficient is k 2 cos β i . However, when β i satisfies ðπ − cos −1 kÞ ≤ β i ≤ π, it is 0. When β i satisfies cos −1 k ≤ β i ≤ ðπ − cos −1 kÞ, it is shown in equation (5).
The formula of the infrared reflection angle coefficient and the solar radiation angle coefficient is shown in equation (6).
where β s is the angle between the normal direction of the sun and the illuminated surface and Φ is the angle between the sun and the earth satellite line, which can be calculated by the position and attitude information provided by the ADCS subsystem. For different satellite schemes, different geometric parameters of thermal analysis node space are obtained according to different main bearing structures and satellite installation layouts. Considering the diversity of satellite structure material and surface material, different physical parameters of thermal analysis nodes are obtained. 3.3.1. FOG. The state variable set of FOG mainly includes satellite angular velocity value ω, phase difference φ s , optical coherence output matrix of coupler C OUT , the total transmission matrix of the optical path E, and the transmission matrix of the light beam propagating clockwise and counter-clockwise E CW and E CCW . The principle model is as follows [29].

Multiphysical Coupling Model in
where λ is the wavelength of the light source, c is the speed of light in vacuum, L fiber is the total length of the optical fiber, D fiber is the diameter of the optical fiber loop, ΔI is the intensity difference of positive and negative half period modulation, α f is the splitting ratio of the coupler, I 0 is the light source intensity, C IN is the light source matrix, C T is the coupler through arm matrix, C P is the polarizer matrix, C C is the coupling arm matrix, T LP is the rotation matrix of the optical path coordinate system, T CW and T CCW are the clockwise and anticlockwise transmission light wave matrix in fiber ring, and M CW and M CCW are the clockwise and anticlockwise transmission light wave matrix in phase modulator. The error model includes the intensity fluctuation and wavelength fluctuation of the light source. The measured values are, respectively, stated as follows [29].
where I 0real and λ real are the real light intensity and the center wavelength, and I 0d and λ d are the light intensity deviation and the random error of the center wavelength. Bring the values in the above formula into the principle model of FOG; the measured angular velocity with error considered can be obtained.

RLG.
The state variable set of RLG mainly includes the satellite angular velocity ω and the frequency difference Δf between the positive and negative beams. The principle model is as follows [30].
where λ is the wavelength of the light source, n r is the normalized refractive index of the optical path, and A laser and L laser are the area and perimeter surrounded by the closed optical path. The error model considers the zero deviation error B 0 and the random walk error B 1 caused by mechanical vibration used to suppressing the lock area. The measured value of Δf is [30]. 7 International Journal of Aerospace Engineering where a 0 , a 1 , and a 2 are the zero bias compensation coefficients obtained by fitting the measurement data, T laser is the gyro temperature, Ω L is the lock zone threshold, Ω Dm is the peak jitter rate, S K is the laser gyro scale factor, Δf m is the frequency difference measurement value, and Δf real is the true value of the frequency difference. Bring the measured value Δf m in the above formula into the principle model of the laser gyro; the measured value of angular velocity after considering the error can be obtained.

Objective Function and Constraints
According to the previous section, different satellite models have different X. But we can consider a general form of objective function and constraints model. Their input variables only include the common parameters of the satellite and subsystems. It is a proper subset of X and does not vary with the satellite model. In the following subsections, the objective function and constraints are introduced from the whole satellite to each subsystem.

Whole
Satellite. J * ðXÞ is the subobjective function, C * ðXÞ is the subconstraint condition, and ω * is the weight coefficient. In this paper, two levels of index system are established, namely, the system index and subsystem index. The system indices include mass, cost, reliability, costeffectiveness ratio, and performance, which are represented by subscripts Mass, Cost, Relia, CER, and Perf. By further decomposing the performance index, the subsystem index is obtained. The total objective function JðXÞ and total constraint CðXÞ are expressed as follows where the objective function and constraint condition algorithm of mass, cost, reliability, cost-effectiveness ratio, and performance are as follows where the total mass of the satellite is expressed as m, the minimum mass of the satellite in the design requirements is expressed as m min , and the maximum mass in the design requirements is expressed as m max . The total cost of the satellite is expressed as C total , the minimum cost of the satellite in the design requirement is expressed as C min , and the maximum cost of the satellite in the design requirement is expressed as C max . The total reliability of the satellite is expressed as R total [31], and the minimum reliability of the satellite in the design requirements is expressed as R min . The main performance indicators of each subsystem are described in the following subsections.

Payload Subsystem.
The indices and constraints of the Payload subsystem are as follows where k Payload is a design constant, P Payload is the total power consumption of the Payload subsystem, D Rc is data rate [32], GSD is ground sample distance, ψ is ground coverage half center angle [33], I d is total information transmission [34], f cov is the coverage performance of target [35], S w is the coverage bandwidth [36], SNR is signal to noise ratio of image [37], S w min is the minimum value of S w , and SNR min is the minimum value of SNR.

ADCS Subsystem. The indices and constraints of the ADCS subsystem are as follows
International Journal of Aerospace Engineering in which k ADCS is a design constant, P ADCS is the total power consumption of the ADCS subsystem [38], δ measure is the attitude measurement accuracy [39], δ point is the attitude direction accuracy [40], δ stable is the attitude stability accuracy, ω maneuver is the attitude maneuver angular velocity, t stable is the stability time after maneuvering, and δ orbit is the orbit control accuracy. a oc is the orbit control acceleration, T w is the flywheel control torque [37], and a oc min is the minimum value of a oc .

Power Subsystem.
The indices and constraints of the Power subsystem are as follows where k Power is a design constant, P Power is the total power consumption of the Power subsystem [32], d SA is the mass to power ratio [41], δ wing is the solar wing orientation accuracy, λ margin is the solar array power margin [42], P BOL is the output power at the beginning of life, P EOL is the output power at the end of life [37], c s is the battery capacity [38], DOD is the discharge depth [43], k e is the average eclipse factor [36], P BOL min is the minimum value of P BOL , P EOL min is the minimum value of P EOL , c r is the electricity demand during the eclipse period, DOD max is the maximum value of DOD, and k e max is the maximum value of k e .

Propulsion Subsystem.
The indices and constraints of the Propulsion subsystem are as follows where k Propulsion is a design constant, m Propulsion is propellant mass [38], I SP is the specific impulse, EOP is the propellant extrusion efficiency, I total SP is the total impulse [40], EOP min is the minimum value of EOP, and I total SP min is the minimum value of I total SP . 4.6. TTC Subsystem. The indices and constraints of the TTC subsystem are as follows and GT > GT min and EIRP > EIRP min , in witch k TTC is a design constant, P TTC is the total power consumption of the TTC subsystem [32], D Trans is the data transmission rate, f TTC cov is the ground station TTC coverage performance [40], t max ttc is the maximum TTC time [44], GT is the quality of the TTC receiving system, EIRP is the equivalent isotropically radiated power [40], t Trans is the data transmission time, GT min is the minimum value of GT, and ERIP min is the minimum value of EIRP. 4.7. Thermal Subsystem. The indices and constraints of the Thermal subsystem are as follows where P Thermal is the total power consumption of the Thermal subsystem and k Thermal is a design constant. T min i and T max i are the minimum and maximum temperatures of typical components, T bal is the equilibrium temperature of the satellite, T min sat and T max sat are the minimum and maximum temperature of satellite in the orbital period [28], T i min and T i max are the upper and lower boundary of typical components' working temperature, T bal min and T bal max are the upper and lower boundary of the satellite equilibrium temperature, and T sat min and T sat max are the upper and lower boundary of the satellite temperature. where σ eq is the structural strength [36], f axi is axial frequency, f lat is transverse frequency [45], d b is the envelope diameter, h b is the envelope height, l w is the body width, d cylinder is the diameter of bearing cylinder [37], V sat is the satellite volume, and Index str is the structural reliability factor [36]. σ max is the ultimate stress required by strength, σ cr is the critical stress required by stability, f axi min is the axial frequency minimum value, f lat min is the transverse frequency minimum value, d min is the effective space diameter of the fairing, h min is the effective height of the fairing, V min is the satellite volume minimum value, and Index min is the structural reliability factor minimum value.

Application of GP Optimization into an Earth Observation Satellite
Using the earth observation satellite described in [46] as an example, the ADCS, Power, Thermal, TTC, Propulsion, Payload, and Structure subsystems are selected to optimize the system composition, structure, and parameters. The weighted coefficients of the main optimization objectives are set as ω Mass = 0:35, ω Cost = 0:2, ω Relia = 0:1, ω CER = 0:1, and ω Perf = 0:25. For GP optimization, the optimizer is set as follows: the number of individuals in each generation is 10, the maximum number of iterations is 100, the number of individuals performing the replication operation is 2, the probability of crossover is 50%, the probability of mutation is 20%, and the evolution termination condition is that the fitness is greater than 0.9. The tree structure of the satellite scheme is taken as the design variable. It mainly includes the parameters and composition of the satellite system and the parameters of the satellite subsystem. System composition includes component type selection and the number of components of each subsystem, and its screening scope is the component library described in subsection 2.2. System structure information refers to the coupling relationship of subsystems, and its model construction and parameter selection are implemented according to subsection 2.3. The main constraints in the satellite subsystem parameter information are shown in Table 2.
In the selected range of design variables shown in column 4 of Table 3, we designed the initial values as shown in column 5. After the optimization of the method proposed in this paper, the optimization result is given in the sixth column, and the change percentage is given in the last column. It can be seen that most of the parameters have been optimized significantly.
The optimization results of the system parameters and composition are shown in Table 4. By comparing the initial scheme with the optimized scheme, we can see that the system composition has significantly changed. The main changes are the type and number of components. It can also be seen from Table 4 that the parameters of components vary with the component models. To sum up, the results show that the GP method in this paper is feasible to realize the simultaneous optimization of satellite system topology and parameters.
In the optimization iterations of GP, the change of the optimal fitness value of each generation is shown in Figure 5(a). The results show that the algorithm has good  convergence and the optimal fitness value decreased rapidly in the early stage, and then, the trend of change slowed down and tended to be flat in the later stage. The normalized value optimization results of performance, quality, cost, reliability, and cost-effectiveness ratio are shown in Figure 5(b). It can be seen that almost all the indices of the optimized scheme are better than those of the initial scheme. The above results show the effectiveness of the GP method in this paper. Table 3 and Figure 5(a) show that the GP methods can also be used to obtain traditional system parameter optimization results. Table 4 shows that the GP method is used to optimize the composition and structure of the satellite system, which is difficult to be realized by traditional algorithms such as GA in [46].
According to the results of the optimal scheme, the sensitivity analysis studies are carried out on two groups of contrast conditions. One is to analyze the influence of the main system parameters on fitness. The other group is to analyze the influence of the system composition and structure on fitness. In terms of satellite system parameters, Figure 6(a) shows the influence of satellite orbital altitude on fitness, and Figure 6(b) shows the influence of satellite total mass on fitness. In terms of satellite system structure, Figure 6(c) shows the influence of satellite camera composition scheme on fitness, and Figure 6(d) shows the influence of different momentum wheel composition scheme on fitness. To sum up, the sensitivity of system parameters to fitness is between 10% and 20%, and the sensitivity of the satellite system structure to fitness is between 10% and 40%. The above results, to some extent, quantitatively prove the importance of satellite composition and structure optimization to the system design.
To analyze the performance of the optimization tool, three groups of working conditions with different population sizes were designed. Under the same computing resources (3 Beihang High Performance Computing nodes, each has 42 physical cores with a frequency of 3.2 GHz), their respective fitness changes are compared, as shown in Figure 7. Three conclusions can be seen from the figure: (1) the larger the population size, the greater the initial fitness; (2) the larger the population size, the less the number of fitness evolution under the same time and resources; and (3) when the population size is 10, the fitness is the best in the current limited time. The results show that the selection of population fitness has a great influence on the optimization results and convergence speed. In future research, further analysis of this relationship can be done to find out more universal rules.

12
International Journal of Aerospace Engineering

Conclusions
In this paper, a GP method of the satellite system design is proposed given the problem that the general optimization algorithms are unable to optimize the system composition structure. Taking an earth observation satellite as an instance, a GP optimization model for satellite system design is developed. The optimization results show that the GP method of satellite system design achieves the simultaneous optimization of satellite system parameters, parameters, and composition. The tree structure and model of satellite established in this paper may have good expansibility. It can be extended to large space systems such as carrier, launch, and ground, or from single satellite to cluster system. This paper can provide a reference for the system design of the aerospace system or cluster system. Furthermore, the method of dimension decomposition combined with GP can be widely used in other systems, which provides a reference scheme for the optimization of the parameters and composition of other systems.
In this paper, a high-precision satellite system model is directly used for optimization, which requires heavy computation. Thus, the GP population size and number of iterations have to be limited. It is recommended that further research could be undertaken as follows: the idea of granular computing can be introduced into GP to solve the problem of optimization convergence caused by a high-precision model.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.