# Genetic Programming Method for Satellite System Topology and Parameter Optimization

Authors: Jinhui Li; Yunfeng Dong; Ming Xu; Hongjue Li
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

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.


---

## Body

## 1. 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. Section2 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.


## 2. GP Optimization Framework for Satellite System Design

### 2.1. Problem Description of Satellite System Design Optimization

The system design scheme setStree 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 Stree. For a specified satellite scheme S, a MSMPC model can be uniquely determined, which is denoted as fX. The relationship between S and fX can be expressed as S→DSfX, where DS represents the modeling method based on the Digital Satellite. X is the set where all the states of the satellite system described infX are included. The system state equation is expressed as Ẋ=fX. The objective function and constraints of the satellite system design are expressed as JX and CX, respectively. The optimization problem of the satellite system design can be described as
(1)find:S∈Stree,min:JX,s.t.:S→DSfX,Ẋ=fX,CX≤0.


### 2.2. Flowchart of Satellite System Design Optimization Based on GP and the Digital Satellite

The optimization framework is shown in Figure1, which is divided into two modules: GP and the Digital Satellite.Figure 1
Flowchart of two-level optimization of satellite system topology and parameter. The component library consists of all component models that can support the satellite system design [14, 25]. FitX is the fitness function.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.In the following two subsections: two key steps of the flowchart will be described, namely, the description of satellite schemeS and the basic GP operations.


### 2.3. Description of Satellite System Design Scheme Based on the Component Library

The satellite system design schemeS 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)
Systemenvironment parameters contain the orbit elements and task start time(2)
For the systemfunction, the task type supported by satellite payload is considered(3)
The systemcomposition mainly describes the subsystems contained in the satellite, the component types under each subsystem, several component models under each component type, and the number of each model(4)
The systemstructure includes the controller structure of ADCS, power supply and distribution relationship of the power subsystem, pipeline and valve definition of the propulsion subsystem, TTC link definition, main bearing structure, and data flow definition(5)
Systemmode refers to the control mode of ADCS(6)
Systemprogram information indicates the ADCS mode switching program and the daily operation process of the satellite(7)
Systemoperation means remote command informationIn this paper, taking the systemcomposition 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.Figure 2
Tree structure topology of satellite system composition schemes.


### 2.4. 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 Figure2, 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.Table 1
Definition of operator set and terminator set.
Node categoryNode symbolCorresponding nodes in the tree structureOperator set“Satellite”The 1st level nodes of the tree structure“Subsystem”The subsystem nodes in the 2nd level“Parameters”The subsystem parameter branch nodes in the 2nd level“Devices”The devices branch nodes in the 2nd level“DeviceCategory”The device category nodes in the 3rd level“DeviceType”The device type nodes in the 4th levelTerminator set“Name”The device type name nodes in the 5th level“NUM”The device number nodes in the 5th level“ParaInfo”The parameter information nodes in the 3rd levelThe 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 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. Figure3 takes the “DeviceCategory” node as an example to show the process of the N+1th generation schemes is produced by crossing and mutating two schemes in the Nth generation.Figure 3
Schematic diagram of satellite scheme crossover and mutation operation.


## 2.1. Problem Description of Satellite System Design Optimization

The system design scheme setStree 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 Stree. For a specified satellite scheme S, a MSMPC model can be uniquely determined, which is denoted as fX. The relationship between S and fX can be expressed as S→DSfX, where DS represents the modeling method based on the Digital Satellite. X is the set where all the states of the satellite system described infX are included. The system state equation is expressed as Ẋ=fX. The objective function and constraints of the satellite system design are expressed as JX and CX, respectively. The optimization problem of the satellite system design can be described as
(1)find:S∈Stree,min:JX,s.t.:S→DSfX,Ẋ=fX,CX≤0.


## 2.2. Flowchart of Satellite System Design Optimization Based on GP and the Digital Satellite

The optimization framework is shown in Figure1, which is divided into two modules: GP and the Digital Satellite.Figure 1
Flowchart of two-level optimization of satellite system topology and parameter. The component library consists of all component models that can support the satellite system design [14, 25]. FitX is the fitness function.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.In the following two subsections: two key steps of the flowchart will be described, namely, the description of satellite schemeS and the basic GP operations.


## 2.3. Description of Satellite System Design Scheme Based on the Component Library

The satellite system design schemeS 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)
Systemenvironment parameters contain the orbit elements and task start time(2)
For the systemfunction, the task type supported by satellite payload is considered(3)
The systemcomposition mainly describes the subsystems contained in the satellite, the component types under each subsystem, several component models under each component type, and the number of each model(4)
The systemstructure includes the controller structure of ADCS, power supply and distribution relationship of the power subsystem, pipeline and valve definition of the propulsion subsystem, TTC link definition, main bearing structure, and data flow definition(5)
Systemmode refers to the control mode of ADCS(6)
Systemprogram information indicates the ADCS mode switching program and the daily operation process of the satellite(7)
Systemoperation means remote command informationIn this paper, taking the systemcomposition 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.Figure 2
Tree structure topology of satellite system composition schemes.


## 2.4. 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 Figure2, 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.Table 1
Definition of operator set and terminator set.
Node categoryNode symbolCorresponding nodes in the tree structureOperator set“Satellite”The 1st level nodes of the tree structure“Subsystem”The subsystem nodes in the 2nd level“Parameters”The subsystem parameter branch nodes in the 2nd level“Devices”The devices branch nodes in the 2nd level“DeviceCategory”The device category nodes in the 3rd level“DeviceType”The device type nodes in the 4th levelTerminator set“Name”The device type name nodes in the 5th level“NUM”The device number nodes in the 5th level“ParaInfo”The parameter information nodes in the 3rd levelThe 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 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. Figure3 takes the “DeviceCategory” node as an example to show the process of the N+1th generation schemes is produced by crossing and mutating two schemes in the Nth generation.Figure 3
Schematic diagram of satellite scheme crossover and mutation operation.


## 3. 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]. For different satellite schemes, the MSMPC models are different. The MSMPC model of satellite mentioned can be decomposed into multiphysical field coupling model of mechanics, electricity, thermodynamics, optics, and electromagnetism in three space scales of a celestial body (coupling relation model on kilometer-scale), satellite (coupling relation model on meter-scale), and component (coupling relation model on centimeter-scale), i.e., celestial-scale, satellite-scale, and component-scale.The analysis of the main coupling relationship is given in subsection3.1. In subsection 3.2, the coupling model in celestial-scale and satellite-scale are described by taking the ADCS and thermal subsystem model as examples. In subsection 3.3, the principle and error model of fiber optic gyro (FOG) and ring laser gyro (RLG) are taken as examples to introduce the multiphysical field coupling model in component-scale.


### 3.1. Analysis of the Main Coupling Relationship

The scale coupling relationship between celestial bodies and satellites is shown in Figure4. 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 wireless signal. The main satellite scale coupling models are represented in the inner frame. The ADCS subsystem outputs the space position R and attitude Q of the satellite which provides the input parameters for the Power, Thermal, and TTC. The thruster switch command from the ADCS subsystem provides the input parameters for the Propulsion subsystem. The Thermal subsystem provides temperature control services for ADCS, Power, Propulsion, and TTC, including monitoring the temperature of each subsystem component through sensors and adjusting the temperature using heaters and cooling windows. The TTC subsystem receives the telemetry data from other subsystems and relays the telecommand data to other subsystems. The Power subsystem distributes power for other subsystems and estimates the power consumption of each subsystem to ensure the power supply balance. The Propulsion subsystem generates force and torque through the thrusters to support the ADCS subsystem. The Structure subsystem influences the ADCS, Propulsion, TTC, and Thermal subsystem through the change of shape.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.


### 3.2. Multiphysical Coupling Model in Celestial- and Satellite-Scale

The celestial-scale and multiphysical coupling model mainly includes the electromechanical thermo-optic-magnetic coupling relationship model between the earth, moon, sun, other vehicles, and satellites. The satellite scale model mainly refers to the coupling model among subsystems, such as ADCS and Thermal in the following paragraphs.


#### 3.2.1. Coupling Models in the ADCS Subsystem

The system state variable set of the ADCS subsystem is a subset of satellite system stateX. It consists of variables in the orbit dynamics and attitude dynamics. In the orbit dynamics, m is the satellite total mass, V is the satellite’s velocity vector, Fout is the combined external force of the environment, and Fflexible, Ftank, and Farm are the forces of flexible attachment, liquid-filled tank, and flexible manipulator on rigid satellite body, respectively. In the attitude dynamics, I is the total inertia, ω is the attitude angular velocity, Mout is the resultant torque of the environment, and Mwheel, Mflexible, Mtank, and Marm are 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].
(2)MM=KMVC,Kf=ΔT×C×ρ×G1.047×10−1×D×ωw,whereMM is the motor output torque, KM is the torque voltage ratio coefficient, VC is the input voltage, Kf 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, Mf0 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].
(3)Mwheel=MM−Mf0−KfωWr.Considering all factors, the orbit and attitude dynamic equations are established as follows [27].
(4)dVdt=1mFout+Fflexible+Ftank+Farm,dωdt=I−1Mout+Mwheel+Mflexible+Mtank+Marm−ω×Iω.According to the different composition of satellite components, the dynamical equations have different expressions. There existsMwheel in the model only when there are momentum exchange devices in the satellite. There exist Fflexible and Mflexible in the model only when there are large flexible solar panels or large antennas in the satellite. The characteristics of the panels and antennas determine the specific expressions. There exist Ftank and Mtank in the model only when there are liquid sloshing tanks in the satellite. The specific expressions are determined by the characteristics of the tank. There exist Farm and Marm in the model only when there are manipulators in the satellite determined. The configurations of the manipulators determine the specific expressions.


#### 3.2.2. Coupling Models in the Thermal Subsystem

The system state variable set of the Thermal subsystem is a subset of satellite system stateX. It includes node temperature Ti, 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 calculation-related parameters. The Thermal subsystem is coupled with the Power subsystem by Pi 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=RE/RE+h, where RE 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−1k, the infrared radiation angle coefficient is k2cosβi. However, when βi satisfies π−cos−1k≤βi≤π, it is 0. When βi satisfies cos−1k≤βi≤π−cos−1k, it is shown in equation (5).
(5)ϕ3i=k2cosβi+1ππ2−1−k2k2−cos2βi−sin−11−k2sinβi−k2cosβicos−11−k2ktanβi.The formula of the infrared reflection angle coefficient and the solar radiation angle coefficient is shown in equation (6).
(6)ϕ2i=ϕ3icosΦ,ϕ1i=cosβs,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. Multiphysical Coupling Model in Component-Scale

The component model includes the principle model, error model, failure model, and telecommunication interface model. Each of the 4 models may be multiphysical coupled. Taking the principle model and error model of FOG and RLG as comparison examples, the multiphysics coupling model of component scale is introduced.


#### 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 COUT, the total transmission matrix of the optical path E, and the transmission matrix of the light beam propagating clockwise and counter-clockwise ECW and ECCW. The principle model is as follows [29].
(7)ω=λc2πLfiberDfiberφs,φs=arcsin−ΔI2αf1‐αf2×traceCOUTI0,COUT=ECINET,E=ECW+ECCW,ECW=CTCPCTTLPTCWMCWCCCPCC,ECCW=CTCPCCTLPTCCWMCCWCCCPCC,where λ is the wavelength of the light source, c is the speed of light in vacuum, Lfiber is the total length of the optical fiber, Dfiber 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, I0 is the light source intensity, CIN is the light source matrix, CT is the coupler through arm matrix, CP is the polarizer matrix, CC is the coupling arm matrix, TLP is the rotation matrix of the optical path coordinate system, TCW and TCCW are the clockwise and anticlockwise transmission light wave matrix in fiber ring, and MCW and MCCW 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].
(8)I0m=I0real+I0d,λm=λreal+λd,where I0real and λreal are the real light intensity and the center wavelength, and I0d 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.


#### 3.3.2. 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].
(9)ω=λLlaserΔf4Alasernr,where λ is the wavelength of the light source, nr is the normalized refractive index of the optical path, and Alaser and Llaser are the area and perimeter surrounded by the closed optical path.The error model considers the zero deviation errorB0 and the random walk error B1 caused by mechanical vibration used to suppressing the lock area. The measured value of Δf is [30].
(10)B0=a0+a1Tlaser+a2Tlaser2,B1=ΩLSK/2πΩDm,Δfm=Δfreal+B0+B1, where a0, a1, and a2 are the zero bias compensation coefficients obtained by fitting the measurement data, Tlaser is the gyro temperature, ΩL is the lock zone threshold, ΩDm is the peak jitter rate, SK is the laser gyro scale factor, Δfm is the frequency difference measurement value, and Δfreal is the true value of the frequency difference. Bring the measured value Δfm 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.


## 3.1. Analysis of the Main Coupling Relationship

The scale coupling relationship between celestial bodies and satellites is shown in Figure4. 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 wireless signal. The main satellite scale coupling models are represented in the inner frame. The ADCS subsystem outputs the space position R and attitude Q of the satellite which provides the input parameters for the Power, Thermal, and TTC. The thruster switch command from the ADCS subsystem provides the input parameters for the Propulsion subsystem. The Thermal subsystem provides temperature control services for ADCS, Power, Propulsion, and TTC, including monitoring the temperature of each subsystem component through sensors and adjusting the temperature using heaters and cooling windows. The TTC subsystem receives the telemetry data from other subsystems and relays the telecommand data to other subsystems. The Power subsystem distributes power for other subsystems and estimates the power consumption of each subsystem to ensure the power supply balance. The Propulsion subsystem generates force and torque through the thrusters to support the ADCS subsystem. The Structure subsystem influences the ADCS, Propulsion, TTC, and Thermal subsystem through the change of shape.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.


## 3.2. Multiphysical Coupling Model in Celestial- and Satellite-Scale

The celestial-scale and multiphysical coupling model mainly includes the electromechanical thermo-optic-magnetic coupling relationship model between the earth, moon, sun, other vehicles, and satellites. The satellite scale model mainly refers to the coupling model among subsystems, such as ADCS and Thermal in the following paragraphs.


### 3.2.1. Coupling Models in the ADCS Subsystem

The system state variable set of the ADCS subsystem is a subset of satellite system stateX. It consists of variables in the orbit dynamics and attitude dynamics. In the orbit dynamics, m is the satellite total mass, V is the satellite’s velocity vector, Fout is the combined external force of the environment, and Fflexible, Ftank, and Farm are the forces of flexible attachment, liquid-filled tank, and flexible manipulator on rigid satellite body, respectively. In the attitude dynamics, I is the total inertia, ω is the attitude angular velocity, Mout is the resultant torque of the environment, and Mwheel, Mflexible, Mtank, and Marm are 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].
(2)MM=KMVC,Kf=ΔT×C×ρ×G1.047×10−1×D×ωw,whereMM is the motor output torque, KM is the torque voltage ratio coefficient, VC is the input voltage, Kf 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, Mf0 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].
(3)Mwheel=MM−Mf0−KfωWr.Considering all factors, the orbit and attitude dynamic equations are established as follows [27].
(4)dVdt=1mFout+Fflexible+Ftank+Farm,dωdt=I−1Mout+Mwheel+Mflexible+Mtank+Marm−ω×Iω.According to the different composition of satellite components, the dynamical equations have different expressions. There existsMwheel in the model only when there are momentum exchange devices in the satellite. There exist Fflexible and Mflexible in the model only when there are large flexible solar panels or large antennas in the satellite. The characteristics of the panels and antennas determine the specific expressions. There exist Ftank and Mtank in the model only when there are liquid sloshing tanks in the satellite. The specific expressions are determined by the characteristics of the tank. There exist Farm and Marm in the model only when there are manipulators in the satellite determined. The configurations of the manipulators determine the specific expressions.


### 3.2.2. Coupling Models in the Thermal Subsystem

The system state variable set of the Thermal subsystem is a subset of satellite system stateX. It includes node temperature Ti, 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 calculation-related parameters. The Thermal subsystem is coupled with the Power subsystem by Pi 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=RE/RE+h, where RE 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−1k, the infrared radiation angle coefficient is k2cosβi. However, when βi satisfies π−cos−1k≤βi≤π, it is 0. When βi satisfies cos−1k≤βi≤π−cos−1k, it is shown in equation (5).
(5)ϕ3i=k2cosβi+1ππ2−1−k2k2−cos2βi−sin−11−k2sinβi−k2cosβicos−11−k2ktanβi.The formula of the infrared reflection angle coefficient and the solar radiation angle coefficient is shown in equation (6).
(6)ϕ2i=ϕ3icosΦ,ϕ1i=cosβs,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.2.1. Coupling Models in the ADCS Subsystem

The system state variable set of the ADCS subsystem is a subset of satellite system stateX. It consists of variables in the orbit dynamics and attitude dynamics. In the orbit dynamics, m is the satellite total mass, V is the satellite’s velocity vector, Fout is the combined external force of the environment, and Fflexible, Ftank, and Farm are the forces of flexible attachment, liquid-filled tank, and flexible manipulator on rigid satellite body, respectively. In the attitude dynamics, I is the total inertia, ω is the attitude angular velocity, Mout is the resultant torque of the environment, and Mwheel, Mflexible, Mtank, and Marm are 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].
(2)MM=KMVC,Kf=ΔT×C×ρ×G1.047×10−1×D×ωw,whereMM is the motor output torque, KM is the torque voltage ratio coefficient, VC is the input voltage, Kf 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, Mf0 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].
(3)Mwheel=MM−Mf0−KfωWr.Considering all factors, the orbit and attitude dynamic equations are established as follows [27].
(4)dVdt=1mFout+Fflexible+Ftank+Farm,dωdt=I−1Mout+Mwheel+Mflexible+Mtank+Marm−ω×Iω.According to the different composition of satellite components, the dynamical equations have different expressions. There existsMwheel in the model only when there are momentum exchange devices in the satellite. There exist Fflexible and Mflexible in the model only when there are large flexible solar panels or large antennas in the satellite. The characteristics of the panels and antennas determine the specific expressions. There exist Ftank and Mtank in the model only when there are liquid sloshing tanks in the satellite. The specific expressions are determined by the characteristics of the tank. There exist Farm and Marm in the model only when there are manipulators in the satellite determined. The configurations of the manipulators determine the specific expressions.


## 3.2.2. Coupling Models in the Thermal Subsystem

The system state variable set of the Thermal subsystem is a subset of satellite system stateX. It includes node temperature Ti, 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 calculation-related parameters. The Thermal subsystem is coupled with the Power subsystem by Pi 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=RE/RE+h, where RE 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−1k, the infrared radiation angle coefficient is k2cosβi. However, when βi satisfies π−cos−1k≤βi≤π, it is 0. When βi satisfies cos−1k≤βi≤π−cos−1k, it is shown in equation (5).
(5)ϕ3i=k2cosβi+1ππ2−1−k2k2−cos2βi−sin−11−k2sinβi−k2cosβicos−11−k2ktanβi.The formula of the infrared reflection angle coefficient and the solar radiation angle coefficient is shown in equation (6).
(6)ϕ2i=ϕ3icosΦ,ϕ1i=cosβs,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. Multiphysical Coupling Model in Component-Scale

The component model includes the principle model, error model, failure model, and telecommunication interface model. Each of the 4 models may be multiphysical coupled. Taking the principle model and error model of FOG and RLG as comparison examples, the multiphysics coupling model of component scale is introduced.


### 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 COUT, the total transmission matrix of the optical path E, and the transmission matrix of the light beam propagating clockwise and counter-clockwise ECW and ECCW. The principle model is as follows [29].
(7)ω=λc2πLfiberDfiberφs,φs=arcsin−ΔI2αf1‐αf2×traceCOUTI0,COUT=ECINET,E=ECW+ECCW,ECW=CTCPCTTLPTCWMCWCCCPCC,ECCW=CTCPCCTLPTCCWMCCWCCCPCC,where λ is the wavelength of the light source, c is the speed of light in vacuum, Lfiber is the total length of the optical fiber, Dfiber 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, I0 is the light source intensity, CIN is the light source matrix, CT is the coupler through arm matrix, CP is the polarizer matrix, CC is the coupling arm matrix, TLP is the rotation matrix of the optical path coordinate system, TCW and TCCW are the clockwise and anticlockwise transmission light wave matrix in fiber ring, and MCW and MCCW 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].
(8)I0m=I0real+I0d,λm=λreal+λd,where I0real and λreal are the real light intensity and the center wavelength, and I0d 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.


### 3.3.2. 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].
(9)ω=λLlaserΔf4Alasernr,where λ is the wavelength of the light source, nr is the normalized refractive index of the optical path, and Alaser and Llaser are the area and perimeter surrounded by the closed optical path.The error model considers the zero deviation errorB0 and the random walk error B1 caused by mechanical vibration used to suppressing the lock area. The measured value of Δf is [30].
(10)B0=a0+a1Tlaser+a2Tlaser2,B1=ΩLSK/2πΩDm,Δfm=Δfreal+B0+B1, where a0, a1, and a2 are the zero bias compensation coefficients obtained by fitting the measurement data, Tlaser is the gyro temperature, ΩL is the lock zone threshold, ΩDm is the peak jitter rate, SK is the laser gyro scale factor, Δfm is the frequency difference measurement value, and Δfreal is the true value of the frequency difference. Bring the measured value Δfm 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.


## 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 COUT, the total transmission matrix of the optical path E, and the transmission matrix of the light beam propagating clockwise and counter-clockwise ECW and ECCW. The principle model is as follows [29].
(7)ω=λc2πLfiberDfiberφs,φs=arcsin−ΔI2αf1‐αf2×traceCOUTI0,COUT=ECINET,E=ECW+ECCW,ECW=CTCPCTTLPTCWMCWCCCPCC,ECCW=CTCPCCTLPTCCWMCCWCCCPCC,where λ is the wavelength of the light source, c is the speed of light in vacuum, Lfiber is the total length of the optical fiber, Dfiber 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, I0 is the light source intensity, CIN is the light source matrix, CT is the coupler through arm matrix, CP is the polarizer matrix, CC is the coupling arm matrix, TLP is the rotation matrix of the optical path coordinate system, TCW and TCCW are the clockwise and anticlockwise transmission light wave matrix in fiber ring, and MCW and MCCW 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].
(8)I0m=I0real+I0d,λm=λreal+λd,where I0real and λreal are the real light intensity and the center wavelength, and I0d 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.


## 3.3.2. 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].
(9)ω=λLlaserΔf4Alasernr,where λ is the wavelength of the light source, nr is the normalized refractive index of the optical path, and Alaser and Llaser are the area and perimeter surrounded by the closed optical path.The error model considers the zero deviation errorB0 and the random walk error B1 caused by mechanical vibration used to suppressing the lock area. The measured value of Δf is [30].
(10)B0=a0+a1Tlaser+a2Tlaser2,B1=ΩLSK/2πΩDm,Δfm=Δfreal+B0+B1, where a0, a1, and a2 are the zero bias compensation coefficients obtained by fitting the measurement data, Tlaser is the gyro temperature, ΩL is the lock zone threshold, ΩDm is the peak jitter rate, SK is the laser gyro scale factor, Δfm is the frequency difference measurement value, and Δfreal is the true value of the frequency difference. Bring the measured value Δfm 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.


## 4. Objective Function and Constraints

According to the previous section, different satellite models have differentX. 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.


### 4.1. 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, cost-effectiveness 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 JX and total constraint CX are expressed as follows
(11)JX=ωMassJMassX+ωCostJCostX+ωReliaJReliaX+ωCERJCERX+ωPerfJPerfX,CX≤0⇔CADCSX≤0,CPowerX≤0,CThermalX≤0,CTTCX≤0CProlusionX≤0,CPayloadX≤0,CStructureX≤0CMassX≤0,CCostX≤0,CReliaX≤0,where the objective function and constraint condition algorithm of mass, cost, reliability, cost-effectiveness ratio, and performance are as follows
(12)JMassX=m−mmin/mmax−mmin,CMassX≤01CMassX>0,CMassX≤0⇔mmin≤m≤mmax,JCostX=Ctotal−Cmin/Cmax−Cmin,CCostX≤01CCostX>0,CCostX≤0⇔Cmin≤Ctotal≤Cmax,JReliaX=Rtotal−Rmin/1−Rmin,CReliaX≤01CReliaX>0,CReliaX≤0⇔Rmin≤Rtotal≤1,JCERX=2πarctanJCostXJPerfX,JPerfX=ωADCSJADCSX+ωPowerJPowerX+ωThermalJThermalX+ωTTCJTTCX+ωProlusionJProlusionX+ωPayloadJPayloadX+ωStructureJStructureX, where the total mass of the satellite is expressed as m, the minimum mass of the satellite in the design requirements is expressed as mmin, and the maximum mass in the design requirements is expressed as mmax. The total cost of the satellite is expressed as Ctotal, the minimum cost of the satellite in the design requirement is expressed as Cmin, and the maximum cost of the satellite in the design requirement is expressed as Cmax. The total reliability of the satellite is expressed as Rtotal [31], and the minimum reliability of the satellite in the design requirements is expressed as Rmin. The main performance indicators of each subsystem are described in the following subsections.


### 4.2. Payload Subsystem

The indices and constraints of the Payload subsystem are as follows(13)JPayloadX=2arctankPayloadPPayloadGSD/ψDRcIdfcov/π,CPayloadX≤01CPayloadX>0,CPayloadX≤0⇔Sw>Sw_min,SNR≥SNRmin,where kPayload is a design constant, PPayload is the total power consumption of the Payload subsystem, DRc is data rate [32], GSD is ground sample distance, ψ is ground coverage half center angle [33], Id is total information transmission [34], fcov is the coverage performance of target [35], Sw is the coverage bandwidth [36], SNR is signal to noise ratio of image [37], Sw_min is the minimum value of Sw, and SNRmin is the minimum value of SNR.


### 4.3. ADCS Subsystem

The indices and constraints of the ADCS subsystem are as follows(14)JADCSX=2arctankADCSPADCSδmeasureδpointδstableδorbittstable/ωmaneuver/π,CADCSX≤01CADCSX>0,CADCSX≤0⇔Tw>Tturb,aoc>aoc_min,in which kADCS is a design constant, PADCS 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, tstable is the stability time after maneuvering, and δorbit is the orbit control accuracy. aoc is the orbit control acceleration, Tw is the flywheel control torque [37], and aoc_min is the minimum value of aoc.


### 4.4. Power Subsystem

The indices and constraints of the Power subsystem are as follows(15)JPowerX=2arctankPowerPPowerdSAδwing/λmargin/π,CPowerX≤01CPowerX>0,CPowerX≤0⇔PBOL≥PBOL_minandPEOL≥PEOL_minandcs≥crandDOD≤DODmaxandke≤ke_max,where kPower is a design constant, PPower is the total power consumption of the Power subsystem [32], dSA is the mass to power ratio [41], δwing is the solar wing orientation accuracy, λmargin is the solar array power margin [42], PBOL is the output power at the beginning of life, PEOL is the output power at the end of life [37], cs is the battery capacity [38], DOD is the discharge depth [43], ke is the average eclipse factor [36], PBOL_min is the minimum value of PBOL, PEOL_min is the minimum value of PEOL, cr is the electricity demand during the eclipse period, DODmax is the maximum value of DOD, and ke_max is the maximum value of ke.


### 4.5. Propulsion Subsystem

The indices and constraints of the Propulsion subsystem are as follows(16)JPropulsionX=2arctankPropulsion/mPropulsionISP/π,CPropulsionX≤01CPropulsionX>0,CPropulsionX≤0⇔EOP>EOPmin,ISPtotal>ISP_mintotal,where kPropulsion is a design constant, mPropulsion is propellant mass [38], ISP is the specific impulse, EOP is the propellant extrusion efficiency, ISPtotal is the total impulse [40], EOPmin is the minimum value of EOP, and ISP_mintotal is the minimum value of ISPtotal.


### 4.6. TTC Subsystem

The indices and constraints of the TTC subsystem are as follows(17)JTTCX=2arctankTTCPTTC/DTransfcovTTC/π,CTTCX≤01CTTCX>0,CTTCX≤0⇔tTrans≤tttcmaxandGT>GTminandEIRP>EIRPmin,in witch kTTC is a design constant, PTTC is the total power consumption of the TTC subsystem [32], DTrans is the data transmission rate, fcovTTC is the ground station TTC coverage performance [40], tttcmax is the maximum TTC time [44], GT is the quality of the TTC receiving system, EIRP is the equivalent isotropically radiated power [40], tTrans is the data transmission time, GTmin is the minimum value of GT, and ERIPmin is the minimum value of EIRP.


### 4.7. Thermal Subsystem

The indices and constraints of the Thermal subsystem are as follows(18)JThermalX=2arctanPThermalkThermal/π,CThermalX≤01CThermalX>0,CThermalX≤0⇔Ti_min≤Timin,Timax≤Ti_maxandTbal_min≤Tbal≤Tbal_maxandTsat_min≤TsatminandTsatmax≤Tsat_max,where PThermal is the total power consumption of the Thermal subsystem and kThermal is a design constant. Timin and Timax are the minimum and maximum temperatures of typical components, Tbal is the equilibrium temperature of the satellite, Tsatmin and Tsatmax are the minimum and maximum temperature of satellite in the orbital period [28], Ti_min and Ti_max are the upper and lower boundary of typical components’ working temperature, Tbal_min and Tbal_max are the upper and lower boundary of the satellite equilibrium temperature, and Tsat_min and Tsat_max are the upper and lower boundary of the satellite temperature.


### 4.8. Structure Subsystem

The indices and constraints of the Structure subsystem are as follows(19)JStructureX=0,X≤01CStructureX>0,JStructureX≤0⇔σeq<σmaxandσeq<σcrandflat>flat_minandfaxi>faxi_minanddb≤dminandhb≤hminanddcylinder≥1/5lwanddcylinder≤1/3lwandVsat≥VminandIndexstr>Indexmin, where σeq is the structural strength [36], faxi is axial frequency, flat is transverse frequency [45], db is the envelope diameter, hb is the envelope height, lw is the body width, dcylinder is the diameter of bearing cylinder [37], Vsat is the satellite volume, and Indexstr is the structural reliability factor [36]. σmax is the ultimate stress required by strength, σcr is the critical stress required by stability, faxi_min is the axial frequency minimum value, flat_min is the transverse frequency minimum value, dmin is the effective space diameter of the fairing, hmin is the effective height of the fairing, Vmin is the satellite volume minimum value, and Indexmin is the structural reliability factor minimum value.


## 4.1. 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, cost-effectiveness 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 JX and total constraint CX are expressed as follows
(11)JX=ωMassJMassX+ωCostJCostX+ωReliaJReliaX+ωCERJCERX+ωPerfJPerfX,CX≤0⇔CADCSX≤0,CPowerX≤0,CThermalX≤0,CTTCX≤0CProlusionX≤0,CPayloadX≤0,CStructureX≤0CMassX≤0,CCostX≤0,CReliaX≤0,where the objective function and constraint condition algorithm of mass, cost, reliability, cost-effectiveness ratio, and performance are as follows
(12)JMassX=m−mmin/mmax−mmin,CMassX≤01CMassX>0,CMassX≤0⇔mmin≤m≤mmax,JCostX=Ctotal−Cmin/Cmax−Cmin,CCostX≤01CCostX>0,CCostX≤0⇔Cmin≤Ctotal≤Cmax,JReliaX=Rtotal−Rmin/1−Rmin,CReliaX≤01CReliaX>0,CReliaX≤0⇔Rmin≤Rtotal≤1,JCERX=2πarctanJCostXJPerfX,JPerfX=ωADCSJADCSX+ωPowerJPowerX+ωThermalJThermalX+ωTTCJTTCX+ωProlusionJProlusionX+ωPayloadJPayloadX+ωStructureJStructureX, where the total mass of the satellite is expressed as m, the minimum mass of the satellite in the design requirements is expressed as mmin, and the maximum mass in the design requirements is expressed as mmax. The total cost of the satellite is expressed as Ctotal, the minimum cost of the satellite in the design requirement is expressed as Cmin, and the maximum cost of the satellite in the design requirement is expressed as Cmax. The total reliability of the satellite is expressed as Rtotal [31], and the minimum reliability of the satellite in the design requirements is expressed as Rmin. The main performance indicators of each subsystem are described in the following subsections.


## 4.2. Payload Subsystem

The indices and constraints of the Payload subsystem are as follows(13)JPayloadX=2arctankPayloadPPayloadGSD/ψDRcIdfcov/π,CPayloadX≤01CPayloadX>0,CPayloadX≤0⇔Sw>Sw_min,SNR≥SNRmin,where kPayload is a design constant, PPayload is the total power consumption of the Payload subsystem, DRc is data rate [32], GSD is ground sample distance, ψ is ground coverage half center angle [33], Id is total information transmission [34], fcov is the coverage performance of target [35], Sw is the coverage bandwidth [36], SNR is signal to noise ratio of image [37], Sw_min is the minimum value of Sw, and SNRmin is the minimum value of SNR.


## 4.3. ADCS Subsystem

The indices and constraints of the ADCS subsystem are as follows(14)JADCSX=2arctankADCSPADCSδmeasureδpointδstableδorbittstable/ωmaneuver/π,CADCSX≤01CADCSX>0,CADCSX≤0⇔Tw>Tturb,aoc>aoc_min,in which kADCS is a design constant, PADCS 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, tstable is the stability time after maneuvering, and δorbit is the orbit control accuracy. aoc is the orbit control acceleration, Tw is the flywheel control torque [37], and aoc_min is the minimum value of aoc.


## 4.4. Power Subsystem

The indices and constraints of the Power subsystem are as follows(15)JPowerX=2arctankPowerPPowerdSAδwing/λmargin/π,CPowerX≤01CPowerX>0,CPowerX≤0⇔PBOL≥PBOL_minandPEOL≥PEOL_minandcs≥crandDOD≤DODmaxandke≤ke_max,where kPower is a design constant, PPower is the total power consumption of the Power subsystem [32], dSA is the mass to power ratio [41], δwing is the solar wing orientation accuracy, λmargin is the solar array power margin [42], PBOL is the output power at the beginning of life, PEOL is the output power at the end of life [37], cs is the battery capacity [38], DOD is the discharge depth [43], ke is the average eclipse factor [36], PBOL_min is the minimum value of PBOL, PEOL_min is the minimum value of PEOL, cr is the electricity demand during the eclipse period, DODmax is the maximum value of DOD, and ke_max is the maximum value of ke.


## 4.5. Propulsion Subsystem

The indices and constraints of the Propulsion subsystem are as follows(16)JPropulsionX=2arctankPropulsion/mPropulsionISP/π,CPropulsionX≤01CPropulsionX>0,CPropulsionX≤0⇔EOP>EOPmin,ISPtotal>ISP_mintotal,where kPropulsion is a design constant, mPropulsion is propellant mass [38], ISP is the specific impulse, EOP is the propellant extrusion efficiency, ISPtotal is the total impulse [40], EOPmin is the minimum value of EOP, and ISP_mintotal is the minimum value of ISPtotal.


## 4.6. TTC Subsystem

The indices and constraints of the TTC subsystem are as follows(17)JTTCX=2arctankTTCPTTC/DTransfcovTTC/π,CTTCX≤01CTTCX>0,CTTCX≤0⇔tTrans≤tttcmaxandGT>GTminandEIRP>EIRPmin,in witch kTTC is a design constant, PTTC is the total power consumption of the TTC subsystem [32], DTrans is the data transmission rate, fcovTTC is the ground station TTC coverage performance [40], tttcmax is the maximum TTC time [44], GT is the quality of the TTC receiving system, EIRP is the equivalent isotropically radiated power [40], tTrans is the data transmission time, GTmin is the minimum value of GT, and ERIPmin is the minimum value of EIRP.


## 4.7. Thermal Subsystem

The indices and constraints of the Thermal subsystem are as follows(18)JThermalX=2arctanPThermalkThermal/π,CThermalX≤01CThermalX>0,CThermalX≤0⇔Ti_min≤Timin,Timax≤Ti_maxandTbal_min≤Tbal≤Tbal_maxandTsat_min≤TsatminandTsatmax≤Tsat_max,where PThermal is the total power consumption of the Thermal subsystem and kThermal is a design constant. Timin and Timax are the minimum and maximum temperatures of typical components, Tbal is the equilibrium temperature of the satellite, Tsatmin and Tsatmax are the minimum and maximum temperature of satellite in the orbital period [28], Ti_min and Ti_max are the upper and lower boundary of typical components’ working temperature, Tbal_min and Tbal_max are the upper and lower boundary of the satellite equilibrium temperature, and Tsat_min and Tsat_max are the upper and lower boundary of the satellite temperature.


## 4.8. Structure Subsystem

The indices and constraints of the Structure subsystem are as follows(19)JStructureX=0,X≤01CStructureX>0,JStructureX≤0⇔σeq<σmaxandσeq<σcrandflat>flat_minandfaxi>faxi_minanddb≤dminandhb≤hminanddcylinder≥1/5lwanddcylinder≤1/3lwandVsat≥VminandIndexstr>Indexmin, where σeq is the structural strength [36], faxi is axial frequency, flat is transverse frequency [45], db is the envelope diameter, hb is the envelope height, lw is the body width, dcylinder is the diameter of bearing cylinder [37], Vsat is the satellite volume, and Indexstr is the structural reliability factor [36]. σmax is the ultimate stress required by strength, σcr is the critical stress required by stability, faxi_min is the axial frequency minimum value, flat_min is the transverse frequency minimum value, dmin is the effective space diameter of the fairing, hmin is the effective height of the fairing, Vmin is the satellite volume minimum value, and Indexmin is the structural reliability factor minimum value.


## 5. 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 subsection2.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.Table 2
Subsystem level constraints.
SubsystemParameter nameSymbolBoundsADCSAttitude measurement accuracy requirements/°δmeasure≤0.001Attitude pointing accuracy requirements/°δpoint≤0.001Attitude stability accuracy requirements/°/sδstable≤0.001Attitude maneuver angular velocity requirements/°/sωmaneuver≥0.01Orbit control accuracy requirements/mδorbit≤1Minimum orbit control acceleration/m/s2aoc≥0.05PayloadMinimum coverage bandwidth/mSw≥50000Minimum image SNR/dbSNR≥50PowerMaximum discharge depth/%DOD≤35Average eclipse factorke≤0.9PropulsionPropellant extrusion efficiencyEOP≥0.9TTCMinimum quality requirement/dbGT≥0.1Minimum requirement ERIP/dbwERIP≥1.0e-6ThermalSatellite equilibrium temperature/°CTbal0~50Satellite temperature range/°CTi0~50StructureAxial frequency/Hzfaxi≥24Lateral frequency/Hzflat≥10Effective space diameter of fairing/mdb≤3.5Effective space height of fairing/mhb≤6Satellite volume/m3Vsat≥1Structural reliability factorIndexstr≥1In the selected range of design variables shown in column 4 of Table3, 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.Table 3
Bounds, initial value, and optimization results of system-level design variable.
SubsystemParameter nameSymbolRangesInitialResultRelative changes (%)ADCSOrbital altitude/kmh[600, 1000]979.3736.1-24.83Orbit inclination/°i[95, 105]104.8100.4-4.20PowerLocal time of descending node/hdnt[8, 12]9.167109.09Life/aLife[3, 10]7.3836.614-10.42StructureTotal mass of satellite/kgm[500, 2000]2000726.2-63.69Envelope diameter/mdb[0.4, 4]21.414-29.30Shell thickness/mts[0.005, 0.5]0.020.01-50.00Envelope height/mhb[0.5, 6]1.4141-29.28Material typeMaterial{1, 2, 3, 4}12\Safety factorIndexmin[1, 1.1]1.011.0786.73Bearing cylinder diameter/mdcylinder[0.1, 2]0.25180.2367-6.00The optimization results of the system parameters and composition are shown in Table4. 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.Table 4
Optimization results of key components in the tree structure.
Component categoryInitial schemeOptimized schemeTypeNumberParametersTypeNumberParametersCameraPL-SC-5m-37.51Focal length 900 mm, view angle 28°, pixel size 20μmPL-SC-5m-501Focal length 1000 mm, view angle 30°, pixel size 10μmWheelMW-1Nm3Angular momentum 0.957 N·m·sMW-2Nm5Angular momentum 1.88 N·m·sSolar wingGaAs-152Area 13.55 m2GaAs-182Area 10.39 m2ThrusterTH-1N6Force 1 NTH-2N8Force 0.1 NAntennaAT-S20002Area 1 m2AT-L15001Area 1.5 m2Hot windowHW-R0902Area 0.305 m2HW-R0882Area 0.3229 m2In the optimization iterations of GP, the change of the optimal fitness value of each generation is shown in Figure5(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.Figure 5
Optimization results of fitness and system indexes. (a) Relationship between GP’s optimal fitness value and generation number of population. (b) Optimization results of normalization system indexes.
(a)(b)The above results show the effectiveness of the GP method in this paper. Table3 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, Figure6(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.Figure 6
Sensitivity analysis of typical system parameters and structure. (a) Sensitivity analysis of orbital altitude. (b) Sensitivity analysis of total mass. (c) Sensitivity analysis of camera composition scheme. (d) Sensitivity analysis of wheel composition scheme.
(a)(b)(c)(d)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 Figure7. 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.Figure 7
Effect of population size on fitness convergence.


## 6. 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.



---

*Source: 6673848-2020-12-14.xml*