- Email: [email protected]
Onboard mission planning for agile satellite using modified mixed-integer linear programming
JID:AESCTE AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.1 (1-13)
Aerospace Science and Technology ••• (••••) •••–•••
1
Contents lists available at ScienceDirect
2
67 68
3
Aerospace Science and Technology
4 5
69 70 71
6
72
www.elsevier.com/locate/aescte
7
73
8
74
9
75
10
76
11 12 13
Onboard mission planning for agile satellite using modified mixed-integer linear programming
14 15 16 17 18
Yuchen She
, Shuang Li
a,∗,1
, Yanbin Zhao
b
23 24 25 26 27 28 29 30 31
79 81 82
a
College of Astronautics, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China b Shanghai Institute of Satellite Engineering, Shanghai, 201109, China
83 84 85
20 22
78 80
a, 1
19 21
77
86
a r t i c l e
i n f o
Article history: Received 2 July 2017 Received in revised form 22 October 2017 Accepted 5 November 2017 Available online xxxx Keywords: Agile satellite Onboard mission planning Earth observation Linear programming
32 33 34
a b s t r a c t
87 88
This paper investigates a new Agile Earth Observation Satellite (AEOS) mission planning algorithm. This new developed algorithm is based on Modified Mixed-Integer Linear Programming (MILP) approach, which takes the minimum slew angle and highest priority criterion. The key point of the paper lies in that the planning process is treated as a Dynamical Combinatorial Optimization (DCO) problem due to the time-varying constraint and the requirement of the problem. The design of the new algorithm is also to meet the real-time requirement of onboard mission planning. The mathematical formulation of the problem is modified to satisfy the demand of the LP solver, and a newly developed dynamic database is used to simulate the constantly changing satellite-target relative position with limited observation windows. The new mission planning model is solved by the LP method to test its reliability. Another solution with Generic Algorithm (GA) is also given to check the result. The calculation time of these two methods is compared to show that the LP approach has a better compatibility with the onboard mission planning requirement. © 2017 Elsevier Masson SAS. All rights reserved.
89 90 91 92 93 94 95 96 97 98 99 100
35
101
36
102
37
103
38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
1. Introduction The Agile Earth Observation Satellite (AEOS) is considered as the most promising representative of the next-generation earth observation platform [1]. The Attitude and Orbit Control System (AOCS) of AEOS is significantly improved recent years, and the observation mode is totally different from those of the traditional earth observation spacecraft. To this reason, the new mission analysis, planning and operation system must be developed to meet the more complex mission requirements [2–4]. The idea of mission management for aircraft and spacecraft is a long discussed topic, and many important works in this region have been reported. Lin addressed the issue of space station overall mission planning by use of decomposition approach [5]. Wang presented the robust gain scheduled control and optimal trajectory planning for an autonomous space rendezvous system subject to input saturation [6]. However, the onboard mission planning and self-management technology is not widely adopted yet. The most scientific observation missions still rely on the ground space cycle operation mode [7]. The newly demonstrated CASPER system is able to realize the mission planning process for EO-1 satel-
59 60 61 62 63 64 65 66
* 1
Corresponding author. E-mail address: [email protected] (S. Li). Department of Astronautics Engineering.
https://doi.org/10.1016/j.ast.2017.11.009 1270-9638/© 2017 Elsevier Masson SAS. All rights reserved.
lite considering the time-varying condition, which is one of the most convincing examples of the onboard planning technology being used in real satellite mission [8]. Whereas, CASPER’s planning service only stops at the general activity level, the detailed satellite dynamic model and optimal control criterion are not considered. In the field of earth observation satellites, the development of the onboard mission planning system is hindered by the limited onboard computing capability. The traditional earth observation satellites only have the capability of slewing along the certain fixed directions. The onboard instruments can take images only when the attitude of satellite is stable. To these reasons, the researches in the field of satellite mission planning are mainly focused on the constellation configuration and multi-satellite cooperation. Examples of such works can be found as the framework of operation strategy of space constellation proposed by Chen [9], and the software architecture design presented by Xu to improve the constellation autonomous formation flying capability [10]. With the continuous development of AOCS technology, more and more attention has been paid to the attitude control of complex structures with higher accuracy and higher angular velocity. An example of such works can be given as the flexible structure control strategy addressed by Cao [11,12]. Computational optimization is one of the major research topics in the artificial intelligence and computer science. The GA based
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
AID:4284 /FLA
2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
[m5G; v1.224; Prn:10/11/2017; 11:12] P.2 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
combinatorial optimization algorithm was presented by Yeniay to solve the optimization problem with various constraints [13]. The advanced optimizer is also adopted in the field of aerospace control. For example, Komfeld proposed a GA based algorithm for the optimal attitude planning [14]. However, the computing capability and the strict real-time constraint of the onboard mission planning are rarely considered in the studies mentioned above. Melton presented a Maximum-Likelihood method in order to achieve the optimal attitude control, which proved that the calculation time can be reduced dramatically with the modified formulation and new developed computing algorithms from the attitude control perspective [15]. The results from the work of Melton prove the feasibility of reducing the calculation time cost to meet the realtime requirement. At the same time, this paper also presents the similar idea with the Linear Programming (LP) solver. The LP technology was firstly introduced by Orourke to solve the general optimization problem (not specified to be in aerospace area) [16]. The concept of LP continues to develop and becomes more efficient, and recent research shows that this approach can be combined with neural network and reinforcement training approaches to realize more complex and effective optimization [17]. The use of LP approach in satellite trajectory and attitude planning has also long been carried out. Coverstone-Carroll adopted the traditional LP method to realize an activity planning for a space robot considering the robot’s orbit dynamics and its orientation cost criterion [18]. Richards proposed a mixed integer optimizer for satellite’s trajectory planning with avoidance constraints [19]. Zhang developed a modified LP algorithm for space rendezvous mission planning, in which the problem of rendezvous and docking maneuver planning was solved by dynamical activity distribution approach [20]. The mission planning for agile satellite is a relatively new research topic [21]. Lemaitre proposed a mission selecting and scheduling algorithm using image processing and area target regrouping approach, which mainly focused on the ground area target analysis and the Field of View (FOV) scanning direction planning of onboard instrument [22]. Based on the discussion above, it is clear that more improvements can be made in this research area. The aim of this paper is to propose a new onboard mission planning method for agile satellite using modified mixed-integer linear programming. The rest of this paper is organized as follows. Section 2 describes the detailed formulation of the onboard mission planning problem of agile satellites. In Section 3, the new model is solved by the LP method to test its reliability, and the solution obtained by Generic Algorithm (GA) is also given to conduct comparative analysis. And finally a conclusion is drawn in Section 4. 2. Formulation of the problem
51 52
2.1. Coordinate system definition
55 56 57 58 59 60 61 62 63 64 65 66
The earth observation satellites usually operate on the sunsynchronized orbit with the inclination of 97 degree and the altitude of 500 km. This kind of orbit enables the satellite to have global earth coverage and relative low altitude to obtain clear images [23,24]. The detailed mission scenario and related parameters in this paper are set in Table 1. The standard observation interval defined in Table 1 is obtained according to the fact that the average access time of a target point on ground to a satellite flying in the orbit of 500 km altitude is about 700 s. Since the aim of the paper is more focused on the mission planning level, the simplified dynamic model of satellite is adopted. It is also assumed that the satellite is able to slew with a constant angular velocity in any direction. The acceleration and
68 69
Related parameters
Value
Standard observation sequence interval Onboard mission planning interval
700 s 5 min before the beginning of observation sequence Anti-sun pointing Target area during day time 6◦ /s
Observation pointing constraint Observation interval constraint Standard slew angular velocity
70 71 72 73 74 75 76
deceleration processes at the beginning and end of slew maneuver are neglected. In order to conduct the earth observation mission, the AEOS satellite needs to slew frequently in engineering practice. To properly describe and define the movement and the attitude of the satellite, the following coordinate systems are defined:
77 78 79 80 81 82 83
1) The earth-centered earth-fixed coordinate system o f – X f Y f Z f : This coordinate system is to define the orientation of satellite and the positions of all the ground targets. The satellite’s orientation is presented in the quaternion form, while the ground targets are presented in a two-dimensional coordinate system: pti = [λi ϕi ] with λi and ϕi denoting the target’s longitude and latitude, respectively. For simplicity, the coordinate system is named as the earth fixed coordinate thereafter. 2) The satellite body coordinate system ob − X b Y b Z b : This coordinate system is defined to calculate the pointing direction of the instrument in the inertial space. The center of the coordinate locates at the center of mass of the satellite: The z direction is defined by the bore sight of the instrument. The y direction is defined by the rotational axis of spacecraft’s solar panel and the x direction completes the right hand rule. 3) The earth J2000 inertial coordination system o j − X j Y j Z j : The definition of this coordination system is traditional. It can be used to calculate the orbit of satellite. From the orbit information, one can determine the observation window of all targets. The coordinate system is named as the earth inertial coordinate for simplicity. 4) The orbit coordinate system oo − X o Y o Z o : This coordinate system is used to describe the position of satellite in orbit. The z axis is collinear to the satellite position vector in the earth J2000 coordinate system. The x axis is perpendicular to z axis and points to the satellite’s velocity direction and the y axis completes the right-hand rule.
84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
Fig. 1 shows the relative position between the earth, the satellite and the ground targets. It also demonstrates all the coordinate systems and the concept of scanning path of the satellite’s instrument during the observation. With the definition of the coordinate systems, the following relative positions can be easily obtained:
⎡
53 54
67
Table 1 Related parameters of mission scenario.
dt f
⎤
⎢ =⎣
f xsat f y sat f zsat
⎤
113 114 115 116 117 118 119
(1)
120 121
⎤
122
i xsat i ⎦ = ⎣ y sat i zsat
⎡
f rsat
⎤
xt f cos(ϕi ) · cos(λi ) = ⎣ yt f ⎦ = re · ⎣ sin(ϕi ) · cos(λi ) ⎦ zt f sin(λi )
⎡
i rsat
⎡
112
123
(2)
124 125
⎡
⎥ ⎦ = (t ) · ⎣
i xsat i y sat i zsat
126
⎤ ⎦
127
(3)
128 129 130 131
f
f
rrel = dt f − rsat
(4)
132
JID:AESCTE AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.3 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
3
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
Fig. 1. Demonstration of coordinate systems and the mission scenario.
35
101
36
102
⎡
37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
i rrel
−1
=
f
xrel
⎤
⎢ f ⎥ ⎥ ·⎢ ⎣ y rel ⎦
(5)
f
zrel sat i rrel = (q¯ ) · rrel
(6)
where dt f denotes the position of a ground target in the earth i fixed coordinate, and re represents the radius of earth. rsat stands f
for the position of satellite in the earth inertial coordinate. rsat denotes the position of satellite in earth fixed frame and represents the transformation matrix from the earth inertial coordinate f i to the earth fixed coordinate. rrel , rrel are the relative position between the target and the satellite in the earth fixed and inertial coordinate, respectively. Eq. (6) can convert the relative position vector into the satellite’s body frame with q¯ the attitude of the satellite, which defines the slew requirement for pointing the instrument to this target.
56 57
2.2. Target list definition and observation window generation
58 59 60 61 62 63 64 65 66
The instrument onboard AEOS satellite usually has only a small FOV due to the demand of high spatial resolution in order to take a clear image of a point target from space [25]. To simulate the mission scenario with this constraint, several point targets are generated within the target area. The target list is defined by seven parameters: target number, latitude, longitude, required observation duration, priority, start time of interest, end time of interest. An example of the target queue is given in Table 2. The duration
is set to be a random value around 30 s to imitate the fact that the instrument needs a certain exposure time to obtain a clear image. The priority is set to be an integer value scaled from 1 to 10 to simulate the increase of importance of the target. The start and end values imitate the duration of interest of the target, which means that the observation needs to be finished within this time period [Start, End] if this target is selected to be observed. In order to ensure an acceptable observation quality, the angle between the bore sight of the instrument and the local normal direction of a ground target should be limited within a certain range. The value of this angle depends on the instrument imaging capability and the requirement of image quality. During the simulation in this paper, this angle value is set to be 45 degree. Given the detailed orbit parameter, one can easily calculate the evolution of relative position between the satellite and the target, and then determine the observation window by calculating the angle using the following formula:
αsat−i (t ) = ag
f −rrel_i , dtf _i
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
(7)
where the operator ag(a, b) denotes the calculation of finding the angle between two vectors. Fig. 2 depicts the definition of an observation window by illustrating the geometry relation between the satellite and a ground target. The validated observation interval can be determined from Eq. (7) and the time of interest of all targets. Fig. 3 shows an example of observation window, time of interest and validated interval for 10 targets. From this figure, one can notice that the observation window obtained by Eq. (7) does not necessarily imply the validated interval. Taking the target 6 as an example, the actual validate interval is much shorter
121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
AID:4284 /FLA
1 2 3 4 5 6 7 8
[m5G; v1.224; Prn:10/11/2017; 11:12] P.4 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
4
67
Table 2 Examples of target list.
68
Number
Latitude (rad)
Longitude (rad)
Duration (s)
Priority
TE Start (s)
TE End (s)
1 2 3 4 5
0.9151 0.3747 −0.7730 −0.0552 0.2552
3.4986 2.1092 3.4898 4.7560 4.1297
32 36 38 39 37
2 6 4 8 5
85.40 160.7 184.61 173.34 99.93
645.94 416.7 454.01 408.34 638.64
69 70 71 72 73 74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
Fig. 2. Illustration of observation window.
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49 50
115
Fig. 3. Example of validated observation interval for 10 random targets.
116
51 52 53 54 55 56
117
than one expected because the time of interest ends earlier than the observation window. This kind of effect must be taken into account during the mission planning algorithm design, which inevitably makes the planning process more complex.
57 58
2.3. Formulation of the problem and dynamic database
118
Table 3 Satellite’s basic parameter.
119 120
System parameters
Value
121
Satellite mass (kg) Moment of inertia tensor (kg·m) Slew mode
700 [2500, 0, 0; 0, 3500, 0; 0, 0, 4000] Constant angular velocity
122
59 60 61 62 63 64 65 66
123 124 125
The mission scenario is assumed to be a task of observing a limited target area containing a certain number of point targets whose parameters are defined following the format of Table 1. As mentioned above, the AEOS satellites are usually placed into the orbit with relatively low altitude, so the mission scenario time is set to be at 720 s. The initial position of satellite in the earth inertial coordinate is set to be at [r 0 0] where r is the orbit altitude
plus the earth radius. The satellite position is calculated in time with traditional orbit dynamics. All orbit disturbance terms are neglected due to the small propagation time [26]. The initial attitude of satellite is chosen to be at Nadir pointing to the earth, which means that the quaternion of satellite in the orbit coordinate system is [0 0 0 1] [27]. Other basic parameters of the AEOS satellite are defined in Table 3.
126 127 128 129 130 131 132
JID:AESCTE AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.5 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
5
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46 47
112
Fig. 4. Algorithm flowchart of mission planning using dynamical database approach.
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
113 114
To conduct the simulation and analyze the mission planning algorithm, a traditional database containing all targets’ information must be provided first. The observed targets and the targets with closed observation windows are removed from the traditional database as the observation goes on, which can reduce the total computation time. In order to facilitate solving the mission planning problem, a dynamic database is introduced in this paper. The aim of the dynamic database is to update the observable target queue at the current time point. In this case, the planning solver will be called multiple times, but the number of feasible solutions per call will be greatly reduced. In contrast, the planning solver is usually called only once in the traditional planning algorithm, which results in a great challenge to find an optimal solution in a set with massive candidate solutions. The detailed mission planning algorithm is shown in Fig. 4, which depicts the relationship between different types of database. In this figure, the notation [tag] includes all the target information given in Table 1. For the sake of simplicity, the Matlab notation is directly
adopted to explain the algorithm in this paper. The notations min(tag_s) and max(tag_e) stand for calculating the min(tag(:, 6)) and max(tag(:, 7)), respectively. The tag(:, 6)) stands for the sixth column of the ‘tag’ matrix. The value of max(tag_e)–min(tag_s) denotes the maximum duration of the observation sequence that can be obtained with the current dynamic database. The observation sequence is then generated for each dynamic database, and the last target from the previous sequence is taken as the initial target for the next one.
115 116 117 118 119 120 121 122 123 124
2.4. LP optimization
125 126
In this section, the LP algorithm is adopted and modified to determine the optimal observation sequence using the dynamical database as input. The optimization criterion is based on the maximum priority, maximum observation number and minimum slew angle. The formulation of the optimization problem can be given as follows:
127 128 129 130 131 132
JID:AESCTE
AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.6 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
6
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
⎧ ⎪ ⎪ Seq = [N ⎛1 N 2 . . . Nn ] ⎪ ⎪ n ⎪ ⎪ ⎪ ⎪ ⎝ k · n + k · Tag ( N i , 5) Max ⎪ Seq 0 1 ⎪ ⎪ ⎪ i = 1 ⎪ ⎪ ⎛ ⎡ ⎤⎞⎞ ⎪ ⎪ n −1 0 ⎪ ⎪ f ⎪ ⎪ ⎝−1 (t ) · r ⎣ 0 ⎦⎠⎠ ¯ k ag q ( t ) · − ⎪ 2 i − 1 rel_i ⎪ ⎪ ⎪ 1 i =1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ST : N i = N j , ∀i = j ⎪ ⎪ ⎪ ⎡ ⎤ ⎪ ⎪ ⎪ 0 ⎨ f q¯ i (t ) · ⎣ 0 ⎦ = −1 (t ) · rrel_i , ∀i = 1, 2, 3, 4 . . . n, ⎪ ⎪ 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ n∀t ∈ [t i t i + t obs_i ] ⎪ ⎪ ⎪ ⎪ ⎪ (t obs_i + T slew_i ) ≤ max _e(Tag) − min _s(Tag) , ⎪ ⎪ ⎪ ⎪ i =1 ⎪ ⎪ ⎪ ⎪ ∀i = 1, 2, 3⎛, 4 . . . ⎪ ⎡ ⎤⎞ ⎪ ⎪ ⎪ 0 ⎪ ⎪ f ⎪ ⎪ ag ⎝−1 (t ) · rrel_i (q¯ i −1 (t )) · ⎣ 0 ⎦⎠ ⎪ ⎪ ⎪ 1 ⎪ ⎪ ⎩ T slew_i =
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
68 69 70 71 72 73 74 75 76
65 66
78 79 80 81 82 83 84 85 86 87 88 89
where k0 , k1 , k2 are positive constants as the criterion coefficient. The notation (q¯ (t )) represents the quaternion orientation transform from the satellite body coordinate to the earth inertial coordinate. The notation ag stands for the minimum slew angle criterion and ω denotes the standard angular velocity of satellite defined in Table 1. The pointing requirement is described in the fourth line of Eq. (8). During the observation time t obs_i of target i, the pointing direction of the instrument in the earth inertial coordinate (left part) should be equal to the relative position vector between the satellite and the target i (right part). The fifth line of Eq. (8) presents the time constraint that all selected targets should be observed before the observation window ends. And the definition of slew time is given in the last line of Eq. (8). It can be noticed that Eq. (8) is a standard optimization problem which can be solved by many existing solvers such as neutral network and GA. The Matlab function ‘ga’ can be directly adopted to solve this problem by considering the fourth line as the nonlinear equality constraint and the fifth line as the linear non-equality constraint. However, the LP approach is adopted in order to solve the problem more efficiently. Then, the formula of this optimization problem needs to be re-written as the following form [27]:
minx f · x
⎧ ⎨ A·x≤b st Aeq · x = beq ⎩ lb ≤ x ≤ ub
(9)
where lb and ub indicate the lower and upper boundary of variable x respectively. In order to convert the problem defined in Eq. (8) into the LP form, the following two steps are needed. First, find the most valuable target in the dynamical database while taking the maximum time interval as the only constraint. Second, determinate the optimal slew sequence based on the minimum slew angle criterion, recalculate the total sequence time considering the slew maneuver and remove the unobservable targets. With this rule, the following transformations are made based on the method proposed in the reference [28]:
63 64
77
(8)
ω
23 24
67
The first step: 1) Set the traditional database as Tag, the format of the database is in the matrix form defined in Table 1.
90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
Fig. 5. Creating process of dynamic database.
106 107
Calculate all the validate interval using Eq. (7) and Fig. 3. The vector of validate interval for each target is noted as
I=
Is1 Ie1
Is2 Ie2
· · · Is N · · · Ie N
109 110 111 112
with Isi and Iei denote the start and end time of validate interval, respectively. 2) Set current time = tm. Arrange the matrix using the following equation:
[ t index ] = sort Tag I(1, :) flag1 = sgn I(2, :) − I(1, :) Tagnew = t , flag1, flag2, Tag(index, :)
108
113 114 115 116 117 118
(10)
where the operator sort is able to put the matrix Tag in ascending order of the first line of the matrix I. The flag1 is used to trigger the removal operation of a target from the traditional database, if the flag1 of a target is negative, the target is not observable in the scenario and it will not be picked into the dynamic database in any situation. The flag2 is the second flag which checks whether the target has been observed already. These modifications can make the matrix be easily adopted by the algorithm defined in Fig. 4. 3) Create the dynamic database Tagdyn using the algorithm shown in Fig. 5.
119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.7 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
4) With the creation of the new matrix Tagdyn , the integer linear form of Eq. (8) can be obtained by:
⎧ f = −k1 · Tagdyn (:, 5) ⎪ ⎪ ⎪ ⎪ ⎪ A = Tagdyn (:, 6) ⎪ ⎪ ⎪ ⎪ ⎨ b = maxTag (:, 7) dyn ⎪ Intcon = 1 : size (Tagdyn ) ⎪ ⎪ ⎪ ⎪ ⎪ lb = zeros 1, size(Tagdyn ) ⎪ ⎪ ⎪ ⎩ ub = ones 1, size(Tagdyn )
(11)
where ones(m, n) and zeros(m, n) create the 1-matrix or 0-matrix with m lines and n columns. Here, the Matlab notations are adopted once again to simplify the expression of the algorithm: size(Tagdyn ) returns the dimension of the matrix. The vector Intcon describes the mixed-integer requirement of x. It is clear by comparing Eq. (11) with Eq. (9) that the LP problem is an exhaustion method which goes through all the possible combination of choices. The solution vector x indicates whether the target is chosen to be observed by given a Boolean number 1 or 0. The equality constraint matrixes in Eq. (9) (Aeq and beq) are set to be empty matrix.
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
The second step: The LP solution of the first step noted as xoutput . The first step provides the best combination of most valuable targets. The sum of their observation time requirements fills the maximum duration of current the dynamic database. The slew sequence is not considered yet in the previous part so that the further optimization needs to be done. To realize the operation, several new matrixes are constructed as:
Tag_list = Tagdyn (i ) xoutput (i ) = 1, ∀i Dir = v Tag_list(2 : 3, :)
(12)
In Eq. (12), Tag_list is the sub matrix of dynamic database containing only the targets chosen by the first step. The operation v (Tag_list (2 : 3, :)) stands for converting the two dimension location (longitude and latitude) into a 3D direction vector using Eq. (1). The matrix Rm is defined to denote the rotation matrix from the earth fixed coordinate to the earth inertial coordinate according to Eq. (5):
Rm = −1 (tm) −1 (tm + 1) · · ·
−1 (tm + N )
Dire_list = P − Rm · Dir
⎧ f = k2 [ agslew1 agslew2 · · · agslew(n−1) ] ⎪ ⎪ ⎪ A eq ⎪ ⎪ ⎡ ⎤ ⎪ ⎪ [zeros (1, size(direction_databse, 1) − 2) 1 1] ⎪ ⎪ ⎪ ⎪ ⎢ ⎥ ⎪ 1 0 ⎥ ⎪ ⎢ zeros (2, size(direction_databse, 1) − 3) 1 ⎪ ⎪ ⎢ ⎥ ⎪ 1 0 1 ⎥ ⎪ ⎢ ⎪ ⎪ ⎢ ⎥ ⎪ =⎢ ⎪ [zeros (3, size(direction_databse, 1) − 4) ones(3, 1) E(3)] ⎥ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎢ ⎥ .. ⎪ ⎪ ⎣ ⎦ ⎪ . ⎪ ⎪ ⎪ ⎪ ( m , 1 ) ones ( size ( direction_databse , 1 ) − 1 , 1 ) E m ( ) [zeros ] ⎪ ⎪ ⎪ ⎪ ⎨ beq = 2 · ones(size(direction_databse, 1) − 1, 1) intcon = 1 : size(Dire_list, 2) 3 ⎪ ⎡⎡ ⎤⎤ ⎪ ⎪ ⎪ A1 ⎪ ⎪ ⎪ ⎢ ⎢ A2 ⎥ ⎥ ⎪ ⎪ ⎢⎢ ⎥⎥ ⎪ ⎪ ⎪ A = ⎢ ⎢ .. ⎥ ⎥ ⎪ ⎢ ⎣ ⎪ . ⎦⎥ ⎪ ⎥ ⎢ ⎪ ⎪ ⎦ ⎣ ⎪ A ⎪ M ⎪ ⎪ ⎪ B ⎪ ⎪ ⎪ ⎪ ⎪ − [Tag_list (:, 4)] ⎪ ⎪ b = ⎪ ⎪ N ⎪ ⎪ ⎪ ⎪ ⎪ lb = zeros (size(Dire_list, 2)/3) ⎩ ub = ones (size(Dire_list, 2)/3) (16)
(14)
where P is the matrix containing the position of satellite obtained by Eq. (3). The matrix Dire_list is a (3N , M) matrix here. N is the number of time step defined in Eq. (13), and M is the number of target in Tag_list. The matrix is converted into a linear form, and then all possible combination of the slew sequence can be constructed as:
⎧ = Dire_list(1, :) Dire_list(2, :) · · · Dire_list(end, :) ⎪ ⎨ Dire_list Minx = 1 2 · · · size(Dire_list, 2) 3 ⎪ ⎩ direction_database = C2Minx
(15) The operation C2Minx denotes the combination operation of choosing two elements from the group of Minx . This operation can be easily realized using Matlab function ‘nchoosek’. Then, the linear model of the second step of the optimization problem can be described as follows:
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91
The operation ag has been defined in Eq. (7). The equality constraints A eq and beq imitate the requirement that all the selected targets are observed only once. The inequality constraints A and b are given to consider the following two factors: All selected target should be observed for enough period of time, and the total slew time should be less than the maximum duration of the current dynamic database. The matrixes Ai and B can be constructed as:
Ai = [zeros (1, I (1, i ) − 1) ones (1, I (2, i ) − I (1, i )) zeros (1, end − I (2, i ))] B = ω · ones (1, Dire_list, 2)
(17)
92 93 94 95 96 97 98 99 100 101 102 103 104
This model can be then fitted to a LP solver. The Matlab function ‘linprog’ can be used to solve the problem. The detailed simulation results are represented in the section 3.
105 106 107 108
3. Simulation and comparison analysis
(13)
With this matrix, the relative direction between the satellite and each target can be easily obtained by: T
7
109 110
In order to confirm the effectiveness of the planning algorithm developed in this paper, the simulation tests are carried out in two parts. In the first part, the algorithms presented in Figs. 4 and 5 are tested with the formulation of the problem defined in Eq. (8) through (17). Several study cases were designed to test the reliability of the algorithm. In the second part, the performance of new developed LP based algorithm is compared with that of the GA solver taking Eq. (8) as input. The calculation results are reported, and the efficiency of the algorithm based on the LP approach is discussed.
111 112 113 114 115 116 117 118 119 120 121
3.1. Test of the algorithm and the formulation
122 123
The first study case is to test the rationality of the basic mission scenario design, which includes the observable windows and orbit calculation. The observation window calculation is based on the satellite’s orbit propagator. Since the total mission scenario is set to be a small duration around 700 s, the orbit propagator for the satellite is chosen to be a two-body model. The orbit inclination is set to be 96 degree, and the altitude is set to be 500 km. The target area is assumed as a rectangular area on the earth surface, which can be defined by the longitude and latitude pairs of
124 125 126 127 128 129 130 131 132
JID:AESCTE
8
AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.8 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31 32
97
Fig. 6. Test 1 mission scenario and observable window for each target.
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48 49
114
Fig. 7. Observation sequence of test 2. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
115 116
the four corner points of the area. Several target points are randomly generated within the area with different level of priority. According to the discussion in Section 2.2, the observable window for each target point can then be obtained using the orbit element. Fig. 6 shows the mission scenario calculation results. The relative positions between the satellite and the target area are depicted. The second simulation case is to test the LP problem formulations defined in Eqs. (10) through (15). 10 point targets are generated within a defined area. The ground track of satellite orbit is set to pass through the target area for simplicity. The observation window, relative positions between the satellite and the targets are all considered. The detailed information of simulation input area given in Table 4, and the corresponding simulation results are shown in Figs. 7 and 8 respectively.
117
Table 4 Simulation input.
118
System parameters
Value
Longitude range Latitude range Total target number Time of interest Minimum period of interest
±1◦ −0.5◦ ∼ 5.4◦ 10 Randomly between 0 to 700 s 100 s
119 120 121 122 123 124 125 126
The red rectangles in Fig. 7 present the observation interval while the green rectangles present the slew maneuver. One can notice that the observation time periods are relatively short and are always contained within the validated interval of the corresponding target (blue rectangles). The slew trajectory is represented in the Fig. 8 along with the target area, non-observed targets and
127 128 129 130 131 132
JID:AESCTE AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.9 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
1
9
67
Table 5 Simulation input.
2 3 4 5 6 7 8
68
System parameters
Value
69
Longitude range Latitude range Total target number Time of interest Minimum period of interest
−1◦
70
∼ 2◦
−0.5◦ ∼ 5.4◦
71
60 Randomly between 0 to 700 s 100 s
72
9 11 12 13 14 15
Fig. 8. Slew trajectory of test 2.
17 18 19 20 21 22 23 24
74 75
10
16
73
satellite’s ground track. In the figure, the satellite moves from the bottom to the top. Since only a small part of the ground track is shown, the simulation considers it as a straight-line form. It can be seen that nine targets out of ten are observed, which proves the efficiency of the mission planning result. At the same time, the scanning trajectory moves forward as the satellite flies over
the target area. The simulation tests confirm the validity of the dynamic database and the LP optimizer for solving the mission planning problem of AEOS satellite. In order to further test the performance of the algorithm developed in this paper, one more study case is proposed with more targets within a slightly rotated area. The basic inputs of the simulation are listed in Table 5, and the corresponding simulation results are depicted in Figs. 9 and 10 respectively. It can be observed from Figs. 9 and 10 that the planning algorithm proposed in this paper has a good capability to handle the complex mission planning problem with a huge amount of input targets. The number of total observed targets is going to reach a saturation status due to the limited total mission scenario time (720 s). However, if the FOV of the instrument is taken into con-
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44 45
110
Fig. 9. Observation sequence of test 3.
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65 66
131
Fig. 10. Slew path of test 3.
132
JID:AESCTE
10
AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.10 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27 28
Fig. 11. Target area coverage with the FOV of instrument.
93 94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
Fig. 12. Results of performance test. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
117 118
sideration, and the assumption is made that all the point targets scanned by the FOV of the instrument is considered to be observed, one can still notice that the LP based mission planning algorithm has a good coverage of the target area. This effect of the FOV of the instrument is illustrated in Fig. 11. It can be seen from Fig. 11 that only less than 10 targets out of sixty are completely non-observed by the instrument FOV, which means that the obtained slew path has a good coverage over the target area. In order to obtain the complete quantitative information about the performance of the new developed algorithm, a series of test is carried out with the different number of candidate targets. The simulation results are given in Fig. 12. It can be observed from Fig. 12 that the number of observed targets (blue curve) increases
with the total number of candidate targets. But this number will reach a saturation status when the number of total candidate targets is very large. This saturation number seems to be less than 30. As the number of observed targets will arrive at the saturation point when the number of total candidate targets is around 60, the input of candidate number larger than sixty is of little significance. The total priority value (red curve) increases with the number of candidate targets, and the observed priority (green curve) is fixed around 50% of the total priority value. It is also noted that the value of the observed priority continues to slowly increase even after the observed number reaches its saturation. The simulation result implies a good and stable performance of the algorithm.
119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.11 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
1 2
67
Table 6 Solver comparison.
68
3
Case 1
4
System parameters
Value
Longitude range Latitude range Total target number Time of interest Minimum period of interest
−1◦
5 6 7 8 9
11
69
Case 2
∼ 2◦
−0.5◦
∼ 5.4◦
40 Randomly between 0 to 400 s 50 s
System parameters
Value
Longitude range Latitude range Total target number Time of interest Minimum period of interest
−1◦
70
∼ 2◦
−0.5◦
∼ 5.4◦
60 Randomly between 0 to 700 s 100 s
10 12 13 14 15 16 17 18 19 20 21 22 23 24 25 27
Fig. 13. Simulation result of Case 1 (target number = 40).
28 29 30
3.2. Comparison of the LP approach with GA solutions
31 32 33 34
72 73 74 75 76
11
26
71
To demonstrate the time efficiency of the planning algorithm developed in this paper, the comparison analysis is conducted between the new method and the traditional GA solver. Since the
existing GA solver can already handle the non-linear constraints and non-linear math model, Eq. (8) is directly taken as input. This test is to demonstrate the performance and the calculation time of the new developed algorithm. To start with, two study cases are proposed and the observation sequences are obtained by both approaches. The input parameters are given in Table 6, and the corresponding simulation results are shown in Figs. 13 and 14 respectively. From the simulation results shown in Figs. 13 and 14, it can be seen that both the solvers can handle the optimization problem properly. The GA solver provides a solution which has a slightly higher observed target number and higher observed target priority. However, the advantage is not quite obvious. In fact, the optimal solution of this problem is not unique due to the saturation effect of the number of observed targets. To this reason, the advantage of the GA solver can become even negligible when the number of total candidate target is large, and the planning results offered by different algorithms do not have a huge impact on the mission. To prove this hypothesis, more simulation tests are conducted, and the results are summarized in Figs. 15 and 16. From the simulation results, it can be noticed that the GA algorithm has a slightly better calculation capability to find the general optimal solution in a complex non-linear optimization problem. However, the advantage is not overwhelming, but limited. The observation sequence
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50 51
116
Fig. 14. Simulation result of Case 2 (target number = 60).
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65 66
131
Fig. 15. Comparison of observed target number and priority of two solvers.
132
JID:AESCTE
AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.12 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
12
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
Fig. 16. Comparison of calculation time of two solvers.
15
81
16 17 18 19 20 21 22 23 24 25 26
82
obtained by GA can only observe 2–3 more targets with limited more priority. However, it can be clearly seen from Fig. 16 that the calculation time required of GA solver to find the solution grows exponentially as the number of total candidate target increases. The complexity of the solver makes it need a huge amount of time and a suitable physical environment to run, which may be unacceptable for space engineering applications. Therefore, the new algorithm based on the LP approach proposed in this paper is a more realistic solution for the onboard mission management and planning problem.
27 28
4. Conclusion
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
A new algorithm is presented in this paper to handle the agile satellite earth observation mission planning problem. The motivation for developing this planning algorithm is based on the dramatic improvement of the capability of the attitude control system for agile satellite. The new developed planning algorithm is based on the modified LP approach, and can make the full use of agile satellite’s capability to increases the observation efficiency. The simulation tests and comparative analysis are carried out to confirm the effectiveness of the planning algorithm proposed in this paper. The simulation results show that the calculation time of GA solver is significantly longer than the new method developed in this paper though the GA solver has a slightly better solution without obvious impact on the observation results. Therefore, conclusions can be drawn that the new proposed algorithm has a good performance to meet the mission planning requirements, and it is much more suitable for onboard applications in engineering. Conflict of interest statement
50 51
None declared.
52 53 54 55 56 57 58 59 60 61
Acknowledgements This work is supported by the National Natural Science Foundation of China (Grant No. 11672126), the National Key Research and Development Program of China (Grant No. 2016YFB0500801), and the Opening Grant from the Key Laboratory of Space Utilization, Chinese Academy of Sciences (Grant No. LSU-2016-04-01). The authors fully appreciate their financial supports.
62 63
References
64 65 66
80
[1] F. Bunkheila, E. Ortore, C. Circi, A new algorithm for agile satellite-based acquisition operations, Acta Astronaut. 12 (3) (2016) 121–128.
[2] Y. Zhang, M. Li, Z. Song, J. Shan, Design and analysis of a moment control unit for agile satellite with high attitude stability requirement, Acta Astronaut. 122 (5) (2016) 90–105. [3] R. Xu, H. Chen, X. Liang, Priority-based constructive algorithms for scheduling agile earth observation satellites with total priority maximization, Expert Syst. Appl. 51 (6) (2016) 195–206. [4] M. Lemaître, G. Verfaillie, F. Jouhaud, Selecting and scheduling observations of agile satellites, Aerosp. Sci. Technol. 6 (5) (2002) 367–381. [5] K. Lin, Y. Luo, J. Zhang, G. Tang, Space station overall mission planning using decomposition approach, Aerosp. Sci. Technol. 33 (1) (2014) 26–39. [6] Q. Wang, B. Zhou, G. Duan, Robust gain scheduled control of spacecraft rendezvous system subject to input saturation, Aerosp. Sci. Technol. 42 (5) (2015) 442–450. [7] J.C. Ong, M. Rackley, Autonomous Operation for the SWIFT Mission, paper for SpaceOps 2002 in Houston, TX, https://doi.org/10.2514/6.2002-T3-47. [8] B. Cichy, S. Chien, Validating the EO-1 autonomous science agent, in: International Workshop on Planning and Scheduling for Space, Germany, Darmstadt, 2004. [9] Y. Chen, Y. Tan, R. He, Multi-satellite mission planning for environmental and disaster monitoring satellite system, in: SpaceOps 2008 Conference, vol. 10(6), 2008, pp. 3469–3488. [10] M. Xu, Y. He, K. Yu, A software architecture design for autonomous formation flying control, IEEE Trans. Aerosp. Electron. Syst. (2017), https://doi.org/10.1109/TAES.2017.2721658. [11] H. Kim, Y. Chang, Mission scheduling optimization of SAR satellite constellation for minimizing system response time, Aerosp. Sci. Technol. 40 (1) (2015) 17–32. [12] X. Cao, C. Yue, M. Liu, Flexible satellite attitude maneuver via constrained torque distribution and active vibration suppression, Aerosp. Sci. Technol. 67 (8) (2017) 387–397. [13] Ö. Yeniay, Penalty function methods for constrained optimization with genetic algorithm, Math. Comput. Appl. 10 (1) (2005) 45–56. [14] R.P. Komfeld, On-board autonomous attitude maneuver planning for planetary spacecraft using genetic algorithms, in: AIAA Guidance, Navigation, and Control Conference and Exhibit, vol. 8, 2003, pp. 2003–5784. [15] R.G. Melton, Maximum-likelihood estimation optimizer for constrained, timeoptimal satellite reorientation, Acta Astronaut. 103 (11) (2014) 185–192. [16] N.W. Orourke, Linear programming for life support optimization, AIAA J. 19 (1) (2002) 2852–2853. [17] Z. Arjmandzadeh, M. Safi, A. Nazemi, A new neural network model for solving random interval linear programming problems, Neural Netw. 89 (3) (2017) 11–18. [18] V.L. Coverstone-Carroll, N.M. Wilkey, Optimal control of a satellite-robot system using direct collocation with non-linear programming, Acta Astronaut. 36 (3) (1995) 149–162. [19] A. Richards, T. Schouwenaars, P. Jonathan, Spacecraft trajectory planning with avoidance constraints using mixed-integer linear programming, J. Guid. Control Dyn. 25 (3) (2002) 755–764. [20] J. Zhang, G. Tang, Y. Luo, H. Li, Orbital rendezvous mission planning using mixed integer nonlinear programming, Acta Astronaut. 68 (7–8) (2011) 1070–1078. [21] F. Yao, J. Li, B. Bai, R. He, Earth observation satellites scheduling based on decomposition optimization algorithm, Int. J. Image Graph. Signal Process. 1 (11) (2010) 10–18. [22] M. Lemaître, G. Verfaillie, F. Jouhaud, Selecting and scheduling observations of agile satellites, Aerosp. Sci. Technol. 6 (3) (2002) 367–381. [23] Y.N. Razoumny, Fundamentals of the route theory for satellite constellation design for Earth discontinuous coverage. Part 4: compound satellite structures on orbits with synchronized nodal regression, Acta Astronaut. 129 (12) (2016) 459–465.
83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.13 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
1 2 3 4 5 6
[24] G. Denis, H. de Boissezon, S. Hosford, The evolution of Earth observation satellites in Europe and its impact on the performance of emergency response services, Acta Astronaut. 127 (10) (2016) 619–633. [25] H. Curtis, Orbit Mechanic for Engineering Students, Elsevier Aerospace Engineering Series, 2013. [26] F.L. Markley, J.L. Grassidis, Fundamentals of Spacecraft Attitude Determination and Control, Springer, ISBN 978-1-4939-0801-1, 2014.
13
[27] R.J. Vanderbei, Linear Programming Foundations and Extensions, International Series in Operations Research & Management Science, 2008. [28] M. Diaby, A linear programming formulation of the traveling salesman problem, in: Proceedings of the 11th WSEAS International Conference on Applied Mathematics, Dallas, Texas, USA, March 22–24, 2007.
67 68 69 70 71 72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65
131
66
132
[m5G; v1.224; Prn:10/11/2017; 11:12] P.1 (1-13)
Aerospace Science and Technology ••• (••••) •••–•••
1
Contents lists available at ScienceDirect
2
67 68
3
Aerospace Science and Technology
4 5
69 70 71
6
72
www.elsevier.com/locate/aescte
7
73
8
74
9
75
10
76
11 12 13
Onboard mission planning for agile satellite using modified mixed-integer linear programming
14 15 16 17 18
Yuchen She
, Shuang Li
a,∗,1
, Yanbin Zhao
b
23 24 25 26 27 28 29 30 31
79 81 82
a
College of Astronautics, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China b Shanghai Institute of Satellite Engineering, Shanghai, 201109, China
83 84 85
20 22
78 80
a, 1
19 21
77
86
a r t i c l e
i n f o
Article history: Received 2 July 2017 Received in revised form 22 October 2017 Accepted 5 November 2017 Available online xxxx Keywords: Agile satellite Onboard mission planning Earth observation Linear programming
32 33 34
a b s t r a c t
87 88
This paper investigates a new Agile Earth Observation Satellite (AEOS) mission planning algorithm. This new developed algorithm is based on Modified Mixed-Integer Linear Programming (MILP) approach, which takes the minimum slew angle and highest priority criterion. The key point of the paper lies in that the planning process is treated as a Dynamical Combinatorial Optimization (DCO) problem due to the time-varying constraint and the requirement of the problem. The design of the new algorithm is also to meet the real-time requirement of onboard mission planning. The mathematical formulation of the problem is modified to satisfy the demand of the LP solver, and a newly developed dynamic database is used to simulate the constantly changing satellite-target relative position with limited observation windows. The new mission planning model is solved by the LP method to test its reliability. Another solution with Generic Algorithm (GA) is also given to check the result. The calculation time of these two methods is compared to show that the LP approach has a better compatibility with the onboard mission planning requirement. © 2017 Elsevier Masson SAS. All rights reserved.
89 90 91 92 93 94 95 96 97 98 99 100
35
101
36
102
37
103
38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
1. Introduction The Agile Earth Observation Satellite (AEOS) is considered as the most promising representative of the next-generation earth observation platform [1]. The Attitude and Orbit Control System (AOCS) of AEOS is significantly improved recent years, and the observation mode is totally different from those of the traditional earth observation spacecraft. To this reason, the new mission analysis, planning and operation system must be developed to meet the more complex mission requirements [2–4]. The idea of mission management for aircraft and spacecraft is a long discussed topic, and many important works in this region have been reported. Lin addressed the issue of space station overall mission planning by use of decomposition approach [5]. Wang presented the robust gain scheduled control and optimal trajectory planning for an autonomous space rendezvous system subject to input saturation [6]. However, the onboard mission planning and self-management technology is not widely adopted yet. The most scientific observation missions still rely on the ground space cycle operation mode [7]. The newly demonstrated CASPER system is able to realize the mission planning process for EO-1 satel-
59 60 61 62 63 64 65 66
* 1
Corresponding author. E-mail address: [email protected] (S. Li). Department of Astronautics Engineering.
https://doi.org/10.1016/j.ast.2017.11.009 1270-9638/© 2017 Elsevier Masson SAS. All rights reserved.
lite considering the time-varying condition, which is one of the most convincing examples of the onboard planning technology being used in real satellite mission [8]. Whereas, CASPER’s planning service only stops at the general activity level, the detailed satellite dynamic model and optimal control criterion are not considered. In the field of earth observation satellites, the development of the onboard mission planning system is hindered by the limited onboard computing capability. The traditional earth observation satellites only have the capability of slewing along the certain fixed directions. The onboard instruments can take images only when the attitude of satellite is stable. To these reasons, the researches in the field of satellite mission planning are mainly focused on the constellation configuration and multi-satellite cooperation. Examples of such works can be found as the framework of operation strategy of space constellation proposed by Chen [9], and the software architecture design presented by Xu to improve the constellation autonomous formation flying capability [10]. With the continuous development of AOCS technology, more and more attention has been paid to the attitude control of complex structures with higher accuracy and higher angular velocity. An example of such works can be given as the flexible structure control strategy addressed by Cao [11,12]. Computational optimization is one of the major research topics in the artificial intelligence and computer science. The GA based
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
AID:4284 /FLA
2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
[m5G; v1.224; Prn:10/11/2017; 11:12] P.2 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
combinatorial optimization algorithm was presented by Yeniay to solve the optimization problem with various constraints [13]. The advanced optimizer is also adopted in the field of aerospace control. For example, Komfeld proposed a GA based algorithm for the optimal attitude planning [14]. However, the computing capability and the strict real-time constraint of the onboard mission planning are rarely considered in the studies mentioned above. Melton presented a Maximum-Likelihood method in order to achieve the optimal attitude control, which proved that the calculation time can be reduced dramatically with the modified formulation and new developed computing algorithms from the attitude control perspective [15]. The results from the work of Melton prove the feasibility of reducing the calculation time cost to meet the realtime requirement. At the same time, this paper also presents the similar idea with the Linear Programming (LP) solver. The LP technology was firstly introduced by Orourke to solve the general optimization problem (not specified to be in aerospace area) [16]. The concept of LP continues to develop and becomes more efficient, and recent research shows that this approach can be combined with neural network and reinforcement training approaches to realize more complex and effective optimization [17]. The use of LP approach in satellite trajectory and attitude planning has also long been carried out. Coverstone-Carroll adopted the traditional LP method to realize an activity planning for a space robot considering the robot’s orbit dynamics and its orientation cost criterion [18]. Richards proposed a mixed integer optimizer for satellite’s trajectory planning with avoidance constraints [19]. Zhang developed a modified LP algorithm for space rendezvous mission planning, in which the problem of rendezvous and docking maneuver planning was solved by dynamical activity distribution approach [20]. The mission planning for agile satellite is a relatively new research topic [21]. Lemaitre proposed a mission selecting and scheduling algorithm using image processing and area target regrouping approach, which mainly focused on the ground area target analysis and the Field of View (FOV) scanning direction planning of onboard instrument [22]. Based on the discussion above, it is clear that more improvements can be made in this research area. The aim of this paper is to propose a new onboard mission planning method for agile satellite using modified mixed-integer linear programming. The rest of this paper is organized as follows. Section 2 describes the detailed formulation of the onboard mission planning problem of agile satellites. In Section 3, the new model is solved by the LP method to test its reliability, and the solution obtained by Generic Algorithm (GA) is also given to conduct comparative analysis. And finally a conclusion is drawn in Section 4. 2. Formulation of the problem
51 52
2.1. Coordinate system definition
55 56 57 58 59 60 61 62 63 64 65 66
The earth observation satellites usually operate on the sunsynchronized orbit with the inclination of 97 degree and the altitude of 500 km. This kind of orbit enables the satellite to have global earth coverage and relative low altitude to obtain clear images [23,24]. The detailed mission scenario and related parameters in this paper are set in Table 1. The standard observation interval defined in Table 1 is obtained according to the fact that the average access time of a target point on ground to a satellite flying in the orbit of 500 km altitude is about 700 s. Since the aim of the paper is more focused on the mission planning level, the simplified dynamic model of satellite is adopted. It is also assumed that the satellite is able to slew with a constant angular velocity in any direction. The acceleration and
68 69
Related parameters
Value
Standard observation sequence interval Onboard mission planning interval
700 s 5 min before the beginning of observation sequence Anti-sun pointing Target area during day time 6◦ /s
Observation pointing constraint Observation interval constraint Standard slew angular velocity
70 71 72 73 74 75 76
deceleration processes at the beginning and end of slew maneuver are neglected. In order to conduct the earth observation mission, the AEOS satellite needs to slew frequently in engineering practice. To properly describe and define the movement and the attitude of the satellite, the following coordinate systems are defined:
77 78 79 80 81 82 83
1) The earth-centered earth-fixed coordinate system o f – X f Y f Z f : This coordinate system is to define the orientation of satellite and the positions of all the ground targets. The satellite’s orientation is presented in the quaternion form, while the ground targets are presented in a two-dimensional coordinate system: pti = [λi ϕi ] with λi and ϕi denoting the target’s longitude and latitude, respectively. For simplicity, the coordinate system is named as the earth fixed coordinate thereafter. 2) The satellite body coordinate system ob − X b Y b Z b : This coordinate system is defined to calculate the pointing direction of the instrument in the inertial space. The center of the coordinate locates at the center of mass of the satellite: The z direction is defined by the bore sight of the instrument. The y direction is defined by the rotational axis of spacecraft’s solar panel and the x direction completes the right hand rule. 3) The earth J2000 inertial coordination system o j − X j Y j Z j : The definition of this coordination system is traditional. It can be used to calculate the orbit of satellite. From the orbit information, one can determine the observation window of all targets. The coordinate system is named as the earth inertial coordinate for simplicity. 4) The orbit coordinate system oo − X o Y o Z o : This coordinate system is used to describe the position of satellite in orbit. The z axis is collinear to the satellite position vector in the earth J2000 coordinate system. The x axis is perpendicular to z axis and points to the satellite’s velocity direction and the y axis completes the right-hand rule.
84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
Fig. 1 shows the relative position between the earth, the satellite and the ground targets. It also demonstrates all the coordinate systems and the concept of scanning path of the satellite’s instrument during the observation. With the definition of the coordinate systems, the following relative positions can be easily obtained:
⎡
53 54
67
Table 1 Related parameters of mission scenario.
dt f
⎤
⎢ =⎣
f xsat f y sat f zsat
⎤
113 114 115 116 117 118 119
(1)
120 121
⎤
122
i xsat i ⎦ = ⎣ y sat i zsat
⎡
f rsat
⎤
xt f cos(ϕi ) · cos(λi ) = ⎣ yt f ⎦ = re · ⎣ sin(ϕi ) · cos(λi ) ⎦ zt f sin(λi )
⎡
i rsat
⎡
112
123
(2)
124 125
⎡
⎥ ⎦ = (t ) · ⎣
i xsat i y sat i zsat
126
⎤ ⎦
127
(3)
128 129 130 131
f
f
rrel = dt f − rsat
(4)
132
JID:AESCTE AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.3 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
3
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
Fig. 1. Demonstration of coordinate systems and the mission scenario.
35
101
36
102
⎡
37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
i rrel
−1
=
f
xrel
⎤
⎢ f ⎥ ⎥ ·⎢ ⎣ y rel ⎦
(5)
f
zrel sat i rrel = (q¯ ) · rrel
(6)
where dt f denotes the position of a ground target in the earth i fixed coordinate, and re represents the radius of earth. rsat stands f
for the position of satellite in the earth inertial coordinate. rsat denotes the position of satellite in earth fixed frame and represents the transformation matrix from the earth inertial coordinate f i to the earth fixed coordinate. rrel , rrel are the relative position between the target and the satellite in the earth fixed and inertial coordinate, respectively. Eq. (6) can convert the relative position vector into the satellite’s body frame with q¯ the attitude of the satellite, which defines the slew requirement for pointing the instrument to this target.
56 57
2.2. Target list definition and observation window generation
58 59 60 61 62 63 64 65 66
The instrument onboard AEOS satellite usually has only a small FOV due to the demand of high spatial resolution in order to take a clear image of a point target from space [25]. To simulate the mission scenario with this constraint, several point targets are generated within the target area. The target list is defined by seven parameters: target number, latitude, longitude, required observation duration, priority, start time of interest, end time of interest. An example of the target queue is given in Table 2. The duration
is set to be a random value around 30 s to imitate the fact that the instrument needs a certain exposure time to obtain a clear image. The priority is set to be an integer value scaled from 1 to 10 to simulate the increase of importance of the target. The start and end values imitate the duration of interest of the target, which means that the observation needs to be finished within this time period [Start, End] if this target is selected to be observed. In order to ensure an acceptable observation quality, the angle between the bore sight of the instrument and the local normal direction of a ground target should be limited within a certain range. The value of this angle depends on the instrument imaging capability and the requirement of image quality. During the simulation in this paper, this angle value is set to be 45 degree. Given the detailed orbit parameter, one can easily calculate the evolution of relative position between the satellite and the target, and then determine the observation window by calculating the angle using the following formula:
αsat−i (t ) = ag
f −rrel_i , dtf _i
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
(7)
where the operator ag(a, b) denotes the calculation of finding the angle between two vectors. Fig. 2 depicts the definition of an observation window by illustrating the geometry relation between the satellite and a ground target. The validated observation interval can be determined from Eq. (7) and the time of interest of all targets. Fig. 3 shows an example of observation window, time of interest and validated interval for 10 targets. From this figure, one can notice that the observation window obtained by Eq. (7) does not necessarily imply the validated interval. Taking the target 6 as an example, the actual validate interval is much shorter
121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
AID:4284 /FLA
1 2 3 4 5 6 7 8
[m5G; v1.224; Prn:10/11/2017; 11:12] P.4 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
4
67
Table 2 Examples of target list.
68
Number
Latitude (rad)
Longitude (rad)
Duration (s)
Priority
TE Start (s)
TE End (s)
1 2 3 4 5
0.9151 0.3747 −0.7730 −0.0552 0.2552
3.4986 2.1092 3.4898 4.7560 4.1297
32 36 38 39 37
2 6 4 8 5
85.40 160.7 184.61 173.34 99.93
645.94 416.7 454.01 408.34 638.64
69 70 71 72 73 74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
Fig. 2. Illustration of observation window.
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49 50
115
Fig. 3. Example of validated observation interval for 10 random targets.
116
51 52 53 54 55 56
117
than one expected because the time of interest ends earlier than the observation window. This kind of effect must be taken into account during the mission planning algorithm design, which inevitably makes the planning process more complex.
57 58
2.3. Formulation of the problem and dynamic database
118
Table 3 Satellite’s basic parameter.
119 120
System parameters
Value
121
Satellite mass (kg) Moment of inertia tensor (kg·m) Slew mode
700 [2500, 0, 0; 0, 3500, 0; 0, 0, 4000] Constant angular velocity
122
59 60 61 62 63 64 65 66
123 124 125
The mission scenario is assumed to be a task of observing a limited target area containing a certain number of point targets whose parameters are defined following the format of Table 1. As mentioned above, the AEOS satellites are usually placed into the orbit with relatively low altitude, so the mission scenario time is set to be at 720 s. The initial position of satellite in the earth inertial coordinate is set to be at [r 0 0] where r is the orbit altitude
plus the earth radius. The satellite position is calculated in time with traditional orbit dynamics. All orbit disturbance terms are neglected due to the small propagation time [26]. The initial attitude of satellite is chosen to be at Nadir pointing to the earth, which means that the quaternion of satellite in the orbit coordinate system is [0 0 0 1] [27]. Other basic parameters of the AEOS satellite are defined in Table 3.
126 127 128 129 130 131 132
JID:AESCTE AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.5 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
5
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46 47
112
Fig. 4. Algorithm flowchart of mission planning using dynamical database approach.
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
113 114
To conduct the simulation and analyze the mission planning algorithm, a traditional database containing all targets’ information must be provided first. The observed targets and the targets with closed observation windows are removed from the traditional database as the observation goes on, which can reduce the total computation time. In order to facilitate solving the mission planning problem, a dynamic database is introduced in this paper. The aim of the dynamic database is to update the observable target queue at the current time point. In this case, the planning solver will be called multiple times, but the number of feasible solutions per call will be greatly reduced. In contrast, the planning solver is usually called only once in the traditional planning algorithm, which results in a great challenge to find an optimal solution in a set with massive candidate solutions. The detailed mission planning algorithm is shown in Fig. 4, which depicts the relationship between different types of database. In this figure, the notation [tag] includes all the target information given in Table 1. For the sake of simplicity, the Matlab notation is directly
adopted to explain the algorithm in this paper. The notations min(tag_s) and max(tag_e) stand for calculating the min(tag(:, 6)) and max(tag(:, 7)), respectively. The tag(:, 6)) stands for the sixth column of the ‘tag’ matrix. The value of max(tag_e)–min(tag_s) denotes the maximum duration of the observation sequence that can be obtained with the current dynamic database. The observation sequence is then generated for each dynamic database, and the last target from the previous sequence is taken as the initial target for the next one.
115 116 117 118 119 120 121 122 123 124
2.4. LP optimization
125 126
In this section, the LP algorithm is adopted and modified to determine the optimal observation sequence using the dynamical database as input. The optimization criterion is based on the maximum priority, maximum observation number and minimum slew angle. The formulation of the optimization problem can be given as follows:
127 128 129 130 131 132
JID:AESCTE
AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.6 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
6
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
⎧ ⎪ ⎪ Seq = [N ⎛1 N 2 . . . Nn ] ⎪ ⎪ n ⎪ ⎪ ⎪ ⎪ ⎝ k · n + k · Tag ( N i , 5) Max ⎪ Seq 0 1 ⎪ ⎪ ⎪ i = 1 ⎪ ⎪ ⎛ ⎡ ⎤⎞⎞ ⎪ ⎪ n −1 0 ⎪ ⎪ f ⎪ ⎪ ⎝−1 (t ) · r ⎣ 0 ⎦⎠⎠ ¯ k ag q ( t ) · − ⎪ 2 i − 1 rel_i ⎪ ⎪ ⎪ 1 i =1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ST : N i = N j , ∀i = j ⎪ ⎪ ⎪ ⎡ ⎤ ⎪ ⎪ ⎪ 0 ⎨ f q¯ i (t ) · ⎣ 0 ⎦ = −1 (t ) · rrel_i , ∀i = 1, 2, 3, 4 . . . n, ⎪ ⎪ 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ n∀t ∈ [t i t i + t obs_i ] ⎪ ⎪ ⎪ ⎪ ⎪ (t obs_i + T slew_i ) ≤ max _e(Tag) − min _s(Tag) , ⎪ ⎪ ⎪ ⎪ i =1 ⎪ ⎪ ⎪ ⎪ ∀i = 1, 2, 3⎛, 4 . . . ⎪ ⎡ ⎤⎞ ⎪ ⎪ ⎪ 0 ⎪ ⎪ f ⎪ ⎪ ag ⎝−1 (t ) · rrel_i (q¯ i −1 (t )) · ⎣ 0 ⎦⎠ ⎪ ⎪ ⎪ 1 ⎪ ⎪ ⎩ T slew_i =
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
68 69 70 71 72 73 74 75 76
65 66
78 79 80 81 82 83 84 85 86 87 88 89
where k0 , k1 , k2 are positive constants as the criterion coefficient. The notation (q¯ (t )) represents the quaternion orientation transform from the satellite body coordinate to the earth inertial coordinate. The notation ag stands for the minimum slew angle criterion and ω denotes the standard angular velocity of satellite defined in Table 1. The pointing requirement is described in the fourth line of Eq. (8). During the observation time t obs_i of target i, the pointing direction of the instrument in the earth inertial coordinate (left part) should be equal to the relative position vector between the satellite and the target i (right part). The fifth line of Eq. (8) presents the time constraint that all selected targets should be observed before the observation window ends. And the definition of slew time is given in the last line of Eq. (8). It can be noticed that Eq. (8) is a standard optimization problem which can be solved by many existing solvers such as neutral network and GA. The Matlab function ‘ga’ can be directly adopted to solve this problem by considering the fourth line as the nonlinear equality constraint and the fifth line as the linear non-equality constraint. However, the LP approach is adopted in order to solve the problem more efficiently. Then, the formula of this optimization problem needs to be re-written as the following form [27]:
minx f · x
⎧ ⎨ A·x≤b st Aeq · x = beq ⎩ lb ≤ x ≤ ub
(9)
where lb and ub indicate the lower and upper boundary of variable x respectively. In order to convert the problem defined in Eq. (8) into the LP form, the following two steps are needed. First, find the most valuable target in the dynamical database while taking the maximum time interval as the only constraint. Second, determinate the optimal slew sequence based on the minimum slew angle criterion, recalculate the total sequence time considering the slew maneuver and remove the unobservable targets. With this rule, the following transformations are made based on the method proposed in the reference [28]:
63 64
77
(8)
ω
23 24
67
The first step: 1) Set the traditional database as Tag, the format of the database is in the matrix form defined in Table 1.
90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
Fig. 5. Creating process of dynamic database.
106 107
Calculate all the validate interval using Eq. (7) and Fig. 3. The vector of validate interval for each target is noted as
I=
Is1 Ie1
Is2 Ie2
· · · Is N · · · Ie N
109 110 111 112
with Isi and Iei denote the start and end time of validate interval, respectively. 2) Set current time = tm. Arrange the matrix using the following equation:
[ t index ] = sort Tag I(1, :) flag1 = sgn I(2, :) − I(1, :) Tagnew = t , flag1, flag2, Tag(index, :)
108
113 114 115 116 117 118
(10)
where the operator sort is able to put the matrix Tag in ascending order of the first line of the matrix I. The flag1 is used to trigger the removal operation of a target from the traditional database, if the flag1 of a target is negative, the target is not observable in the scenario and it will not be picked into the dynamic database in any situation. The flag2 is the second flag which checks whether the target has been observed already. These modifications can make the matrix be easily adopted by the algorithm defined in Fig. 4. 3) Create the dynamic database Tagdyn using the algorithm shown in Fig. 5.
119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.7 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
4) With the creation of the new matrix Tagdyn , the integer linear form of Eq. (8) can be obtained by:
⎧ f = −k1 · Tagdyn (:, 5) ⎪ ⎪ ⎪ ⎪ ⎪ A = Tagdyn (:, 6) ⎪ ⎪ ⎪ ⎪ ⎨ b = maxTag (:, 7) dyn ⎪ Intcon = 1 : size (Tagdyn ) ⎪ ⎪ ⎪ ⎪ ⎪ lb = zeros 1, size(Tagdyn ) ⎪ ⎪ ⎪ ⎩ ub = ones 1, size(Tagdyn )
(11)
where ones(m, n) and zeros(m, n) create the 1-matrix or 0-matrix with m lines and n columns. Here, the Matlab notations are adopted once again to simplify the expression of the algorithm: size(Tagdyn ) returns the dimension of the matrix. The vector Intcon describes the mixed-integer requirement of x. It is clear by comparing Eq. (11) with Eq. (9) that the LP problem is an exhaustion method which goes through all the possible combination of choices. The solution vector x indicates whether the target is chosen to be observed by given a Boolean number 1 or 0. The equality constraint matrixes in Eq. (9) (Aeq and beq) are set to be empty matrix.
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
The second step: The LP solution of the first step noted as xoutput . The first step provides the best combination of most valuable targets. The sum of their observation time requirements fills the maximum duration of current the dynamic database. The slew sequence is not considered yet in the previous part so that the further optimization needs to be done. To realize the operation, several new matrixes are constructed as:
Tag_list = Tagdyn (i ) xoutput (i ) = 1, ∀i Dir = v Tag_list(2 : 3, :)
(12)
In Eq. (12), Tag_list is the sub matrix of dynamic database containing only the targets chosen by the first step. The operation v (Tag_list (2 : 3, :)) stands for converting the two dimension location (longitude and latitude) into a 3D direction vector using Eq. (1). The matrix Rm is defined to denote the rotation matrix from the earth fixed coordinate to the earth inertial coordinate according to Eq. (5):
Rm = −1 (tm) −1 (tm + 1) · · ·
−1 (tm + N )
Dire_list = P − Rm · Dir
⎧ f = k2 [ agslew1 agslew2 · · · agslew(n−1) ] ⎪ ⎪ ⎪ A eq ⎪ ⎪ ⎡ ⎤ ⎪ ⎪ [zeros (1, size(direction_databse, 1) − 2) 1 1] ⎪ ⎪ ⎪ ⎪ ⎢ ⎥ ⎪ 1 0 ⎥ ⎪ ⎢ zeros (2, size(direction_databse, 1) − 3) 1 ⎪ ⎪ ⎢ ⎥ ⎪ 1 0 1 ⎥ ⎪ ⎢ ⎪ ⎪ ⎢ ⎥ ⎪ =⎢ ⎪ [zeros (3, size(direction_databse, 1) − 4) ones(3, 1) E(3)] ⎥ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎢ ⎥ .. ⎪ ⎪ ⎣ ⎦ ⎪ . ⎪ ⎪ ⎪ ⎪ ( m , 1 ) ones ( size ( direction_databse , 1 ) − 1 , 1 ) E m ( ) [zeros ] ⎪ ⎪ ⎪ ⎪ ⎨ beq = 2 · ones(size(direction_databse, 1) − 1, 1) intcon = 1 : size(Dire_list, 2) 3 ⎪ ⎡⎡ ⎤⎤ ⎪ ⎪ ⎪ A1 ⎪ ⎪ ⎪ ⎢ ⎢ A2 ⎥ ⎥ ⎪ ⎪ ⎢⎢ ⎥⎥ ⎪ ⎪ ⎪ A = ⎢ ⎢ .. ⎥ ⎥ ⎪ ⎢ ⎣ ⎪ . ⎦⎥ ⎪ ⎥ ⎢ ⎪ ⎪ ⎦ ⎣ ⎪ A ⎪ M ⎪ ⎪ ⎪ B ⎪ ⎪ ⎪ ⎪ ⎪ − [Tag_list (:, 4)] ⎪ ⎪ b = ⎪ ⎪ N ⎪ ⎪ ⎪ ⎪ ⎪ lb = zeros (size(Dire_list, 2)/3) ⎩ ub = ones (size(Dire_list, 2)/3) (16)
(14)
where P is the matrix containing the position of satellite obtained by Eq. (3). The matrix Dire_list is a (3N , M) matrix here. N is the number of time step defined in Eq. (13), and M is the number of target in Tag_list. The matrix is converted into a linear form, and then all possible combination of the slew sequence can be constructed as:
⎧ = Dire_list(1, :) Dire_list(2, :) · · · Dire_list(end, :) ⎪ ⎨ Dire_list Minx = 1 2 · · · size(Dire_list, 2) 3 ⎪ ⎩ direction_database = C2Minx
(15) The operation C2Minx denotes the combination operation of choosing two elements from the group of Minx . This operation can be easily realized using Matlab function ‘nchoosek’. Then, the linear model of the second step of the optimization problem can be described as follows:
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91
The operation ag has been defined in Eq. (7). The equality constraints A eq and beq imitate the requirement that all the selected targets are observed only once. The inequality constraints A and b are given to consider the following two factors: All selected target should be observed for enough period of time, and the total slew time should be less than the maximum duration of the current dynamic database. The matrixes Ai and B can be constructed as:
Ai = [zeros (1, I (1, i ) − 1) ones (1, I (2, i ) − I (1, i )) zeros (1, end − I (2, i ))] B = ω · ones (1, Dire_list, 2)
(17)
92 93 94 95 96 97 98 99 100 101 102 103 104
This model can be then fitted to a LP solver. The Matlab function ‘linprog’ can be used to solve the problem. The detailed simulation results are represented in the section 3.
105 106 107 108
3. Simulation and comparison analysis
(13)
With this matrix, the relative direction between the satellite and each target can be easily obtained by: T
7
109 110
In order to confirm the effectiveness of the planning algorithm developed in this paper, the simulation tests are carried out in two parts. In the first part, the algorithms presented in Figs. 4 and 5 are tested with the formulation of the problem defined in Eq. (8) through (17). Several study cases were designed to test the reliability of the algorithm. In the second part, the performance of new developed LP based algorithm is compared with that of the GA solver taking Eq. (8) as input. The calculation results are reported, and the efficiency of the algorithm based on the LP approach is discussed.
111 112 113 114 115 116 117 118 119 120 121
3.1. Test of the algorithm and the formulation
122 123
The first study case is to test the rationality of the basic mission scenario design, which includes the observable windows and orbit calculation. The observation window calculation is based on the satellite’s orbit propagator. Since the total mission scenario is set to be a small duration around 700 s, the orbit propagator for the satellite is chosen to be a two-body model. The orbit inclination is set to be 96 degree, and the altitude is set to be 500 km. The target area is assumed as a rectangular area on the earth surface, which can be defined by the longitude and latitude pairs of
124 125 126 127 128 129 130 131 132
JID:AESCTE
8
AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.8 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31 32
97
Fig. 6. Test 1 mission scenario and observable window for each target.
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48 49
114
Fig. 7. Observation sequence of test 2. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
115 116
the four corner points of the area. Several target points are randomly generated within the area with different level of priority. According to the discussion in Section 2.2, the observable window for each target point can then be obtained using the orbit element. Fig. 6 shows the mission scenario calculation results. The relative positions between the satellite and the target area are depicted. The second simulation case is to test the LP problem formulations defined in Eqs. (10) through (15). 10 point targets are generated within a defined area. The ground track of satellite orbit is set to pass through the target area for simplicity. The observation window, relative positions between the satellite and the targets are all considered. The detailed information of simulation input area given in Table 4, and the corresponding simulation results are shown in Figs. 7 and 8 respectively.
117
Table 4 Simulation input.
118
System parameters
Value
Longitude range Latitude range Total target number Time of interest Minimum period of interest
±1◦ −0.5◦ ∼ 5.4◦ 10 Randomly between 0 to 700 s 100 s
119 120 121 122 123 124 125 126
The red rectangles in Fig. 7 present the observation interval while the green rectangles present the slew maneuver. One can notice that the observation time periods are relatively short and are always contained within the validated interval of the corresponding target (blue rectangles). The slew trajectory is represented in the Fig. 8 along with the target area, non-observed targets and
127 128 129 130 131 132
JID:AESCTE AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.9 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
1
9
67
Table 5 Simulation input.
2 3 4 5 6 7 8
68
System parameters
Value
69
Longitude range Latitude range Total target number Time of interest Minimum period of interest
−1◦
70
∼ 2◦
−0.5◦ ∼ 5.4◦
71
60 Randomly between 0 to 700 s 100 s
72
9 11 12 13 14 15
Fig. 8. Slew trajectory of test 2.
17 18 19 20 21 22 23 24
74 75
10
16
73
satellite’s ground track. In the figure, the satellite moves from the bottom to the top. Since only a small part of the ground track is shown, the simulation considers it as a straight-line form. It can be seen that nine targets out of ten are observed, which proves the efficiency of the mission planning result. At the same time, the scanning trajectory moves forward as the satellite flies over
the target area. The simulation tests confirm the validity of the dynamic database and the LP optimizer for solving the mission planning problem of AEOS satellite. In order to further test the performance of the algorithm developed in this paper, one more study case is proposed with more targets within a slightly rotated area. The basic inputs of the simulation are listed in Table 5, and the corresponding simulation results are depicted in Figs. 9 and 10 respectively. It can be observed from Figs. 9 and 10 that the planning algorithm proposed in this paper has a good capability to handle the complex mission planning problem with a huge amount of input targets. The number of total observed targets is going to reach a saturation status due to the limited total mission scenario time (720 s). However, if the FOV of the instrument is taken into con-
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44 45
110
Fig. 9. Observation sequence of test 3.
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65 66
131
Fig. 10. Slew path of test 3.
132
JID:AESCTE
10
AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.10 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27 28
Fig. 11. Target area coverage with the FOV of instrument.
93 94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
Fig. 12. Results of performance test. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
117 118
sideration, and the assumption is made that all the point targets scanned by the FOV of the instrument is considered to be observed, one can still notice that the LP based mission planning algorithm has a good coverage of the target area. This effect of the FOV of the instrument is illustrated in Fig. 11. It can be seen from Fig. 11 that only less than 10 targets out of sixty are completely non-observed by the instrument FOV, which means that the obtained slew path has a good coverage over the target area. In order to obtain the complete quantitative information about the performance of the new developed algorithm, a series of test is carried out with the different number of candidate targets. The simulation results are given in Fig. 12. It can be observed from Fig. 12 that the number of observed targets (blue curve) increases
with the total number of candidate targets. But this number will reach a saturation status when the number of total candidate targets is very large. This saturation number seems to be less than 30. As the number of observed targets will arrive at the saturation point when the number of total candidate targets is around 60, the input of candidate number larger than sixty is of little significance. The total priority value (red curve) increases with the number of candidate targets, and the observed priority (green curve) is fixed around 50% of the total priority value. It is also noted that the value of the observed priority continues to slowly increase even after the observed number reaches its saturation. The simulation result implies a good and stable performance of the algorithm.
119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.11 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
1 2
67
Table 6 Solver comparison.
68
3
Case 1
4
System parameters
Value
Longitude range Latitude range Total target number Time of interest Minimum period of interest
−1◦
5 6 7 8 9
11
69
Case 2
∼ 2◦
−0.5◦
∼ 5.4◦
40 Randomly between 0 to 400 s 50 s
System parameters
Value
Longitude range Latitude range Total target number Time of interest Minimum period of interest
−1◦
70
∼ 2◦
−0.5◦
∼ 5.4◦
60 Randomly between 0 to 700 s 100 s
10 12 13 14 15 16 17 18 19 20 21 22 23 24 25 27
Fig. 13. Simulation result of Case 1 (target number = 40).
28 29 30
3.2. Comparison of the LP approach with GA solutions
31 32 33 34
72 73 74 75 76
11
26
71
To demonstrate the time efficiency of the planning algorithm developed in this paper, the comparison analysis is conducted between the new method and the traditional GA solver. Since the
existing GA solver can already handle the non-linear constraints and non-linear math model, Eq. (8) is directly taken as input. This test is to demonstrate the performance and the calculation time of the new developed algorithm. To start with, two study cases are proposed and the observation sequences are obtained by both approaches. The input parameters are given in Table 6, and the corresponding simulation results are shown in Figs. 13 and 14 respectively. From the simulation results shown in Figs. 13 and 14, it can be seen that both the solvers can handle the optimization problem properly. The GA solver provides a solution which has a slightly higher observed target number and higher observed target priority. However, the advantage is not quite obvious. In fact, the optimal solution of this problem is not unique due to the saturation effect of the number of observed targets. To this reason, the advantage of the GA solver can become even negligible when the number of total candidate target is large, and the planning results offered by different algorithms do not have a huge impact on the mission. To prove this hypothesis, more simulation tests are conducted, and the results are summarized in Figs. 15 and 16. From the simulation results, it can be noticed that the GA algorithm has a slightly better calculation capability to find the general optimal solution in a complex non-linear optimization problem. However, the advantage is not overwhelming, but limited. The observation sequence
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50 51
116
Fig. 14. Simulation result of Case 2 (target number = 60).
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65 66
131
Fig. 15. Comparison of observed target number and priority of two solvers.
132
JID:AESCTE
AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.12 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
12
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
Fig. 16. Comparison of calculation time of two solvers.
15
81
16 17 18 19 20 21 22 23 24 25 26
82
obtained by GA can only observe 2–3 more targets with limited more priority. However, it can be clearly seen from Fig. 16 that the calculation time required of GA solver to find the solution grows exponentially as the number of total candidate target increases. The complexity of the solver makes it need a huge amount of time and a suitable physical environment to run, which may be unacceptable for space engineering applications. Therefore, the new algorithm based on the LP approach proposed in this paper is a more realistic solution for the onboard mission management and planning problem.
27 28
4. Conclusion
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
A new algorithm is presented in this paper to handle the agile satellite earth observation mission planning problem. The motivation for developing this planning algorithm is based on the dramatic improvement of the capability of the attitude control system for agile satellite. The new developed planning algorithm is based on the modified LP approach, and can make the full use of agile satellite’s capability to increases the observation efficiency. The simulation tests and comparative analysis are carried out to confirm the effectiveness of the planning algorithm proposed in this paper. The simulation results show that the calculation time of GA solver is significantly longer than the new method developed in this paper though the GA solver has a slightly better solution without obvious impact on the observation results. Therefore, conclusions can be drawn that the new proposed algorithm has a good performance to meet the mission planning requirements, and it is much more suitable for onboard applications in engineering. Conflict of interest statement
50 51
None declared.
52 53 54 55 56 57 58 59 60 61
Acknowledgements This work is supported by the National Natural Science Foundation of China (Grant No. 11672126), the National Key Research and Development Program of China (Grant No. 2016YFB0500801), and the Opening Grant from the Key Laboratory of Space Utilization, Chinese Academy of Sciences (Grant No. LSU-2016-04-01). The authors fully appreciate their financial supports.
62 63
References
64 65 66
80
[1] F. Bunkheila, E. Ortore, C. Circi, A new algorithm for agile satellite-based acquisition operations, Acta Astronaut. 12 (3) (2016) 121–128.
[2] Y. Zhang, M. Li, Z. Song, J. Shan, Design and analysis of a moment control unit for agile satellite with high attitude stability requirement, Acta Astronaut. 122 (5) (2016) 90–105. [3] R. Xu, H. Chen, X. Liang, Priority-based constructive algorithms for scheduling agile earth observation satellites with total priority maximization, Expert Syst. Appl. 51 (6) (2016) 195–206. [4] M. Lemaître, G. Verfaillie, F. Jouhaud, Selecting and scheduling observations of agile satellites, Aerosp. Sci. Technol. 6 (5) (2002) 367–381. [5] K. Lin, Y. Luo, J. Zhang, G. Tang, Space station overall mission planning using decomposition approach, Aerosp. Sci. Technol. 33 (1) (2014) 26–39. [6] Q. Wang, B. Zhou, G. Duan, Robust gain scheduled control of spacecraft rendezvous system subject to input saturation, Aerosp. Sci. Technol. 42 (5) (2015) 442–450. [7] J.C. Ong, M. Rackley, Autonomous Operation for the SWIFT Mission, paper for SpaceOps 2002 in Houston, TX, https://doi.org/10.2514/6.2002-T3-47. [8] B. Cichy, S. Chien, Validating the EO-1 autonomous science agent, in: International Workshop on Planning and Scheduling for Space, Germany, Darmstadt, 2004. [9] Y. Chen, Y. Tan, R. He, Multi-satellite mission planning for environmental and disaster monitoring satellite system, in: SpaceOps 2008 Conference, vol. 10(6), 2008, pp. 3469–3488. [10] M. Xu, Y. He, K. Yu, A software architecture design for autonomous formation flying control, IEEE Trans. Aerosp. Electron. Syst. (2017), https://doi.org/10.1109/TAES.2017.2721658. [11] H. Kim, Y. Chang, Mission scheduling optimization of SAR satellite constellation for minimizing system response time, Aerosp. Sci. Technol. 40 (1) (2015) 17–32. [12] X. Cao, C. Yue, M. Liu, Flexible satellite attitude maneuver via constrained torque distribution and active vibration suppression, Aerosp. Sci. Technol. 67 (8) (2017) 387–397. [13] Ö. Yeniay, Penalty function methods for constrained optimization with genetic algorithm, Math. Comput. Appl. 10 (1) (2005) 45–56. [14] R.P. Komfeld, On-board autonomous attitude maneuver planning for planetary spacecraft using genetic algorithms, in: AIAA Guidance, Navigation, and Control Conference and Exhibit, vol. 8, 2003, pp. 2003–5784. [15] R.G. Melton, Maximum-likelihood estimation optimizer for constrained, timeoptimal satellite reorientation, Acta Astronaut. 103 (11) (2014) 185–192. [16] N.W. Orourke, Linear programming for life support optimization, AIAA J. 19 (1) (2002) 2852–2853. [17] Z. Arjmandzadeh, M. Safi, A. Nazemi, A new neural network model for solving random interval linear programming problems, Neural Netw. 89 (3) (2017) 11–18. [18] V.L. Coverstone-Carroll, N.M. Wilkey, Optimal control of a satellite-robot system using direct collocation with non-linear programming, Acta Astronaut. 36 (3) (1995) 149–162. [19] A. Richards, T. Schouwenaars, P. Jonathan, Spacecraft trajectory planning with avoidance constraints using mixed-integer linear programming, J. Guid. Control Dyn. 25 (3) (2002) 755–764. [20] J. Zhang, G. Tang, Y. Luo, H. Li, Orbital rendezvous mission planning using mixed integer nonlinear programming, Acta Astronaut. 68 (7–8) (2011) 1070–1078. [21] F. Yao, J. Li, B. Bai, R. He, Earth observation satellites scheduling based on decomposition optimization algorithm, Int. J. Image Graph. Signal Process. 1 (11) (2010) 10–18. [22] M. Lemaître, G. Verfaillie, F. Jouhaud, Selecting and scheduling observations of agile satellites, Aerosp. Sci. Technol. 6 (3) (2002) 367–381. [23] Y.N. Razoumny, Fundamentals of the route theory for satellite constellation design for Earth discontinuous coverage. Part 4: compound satellite structures on orbits with synchronized nodal regression, Acta Astronaut. 129 (12) (2016) 459–465.
83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:4284 /FLA
[m5G; v1.224; Prn:10/11/2017; 11:12] P.13 (1-13)
Y. She et al. / Aerospace Science and Technology ••• (••••) •••–•••
1 2 3 4 5 6
[24] G. Denis, H. de Boissezon, S. Hosford, The evolution of Earth observation satellites in Europe and its impact on the performance of emergency response services, Acta Astronaut. 127 (10) (2016) 619–633. [25] H. Curtis, Orbit Mechanic for Engineering Students, Elsevier Aerospace Engineering Series, 2013. [26] F.L. Markley, J.L. Grassidis, Fundamentals of Spacecraft Attitude Determination and Control, Springer, ISBN 978-1-4939-0801-1, 2014.
13
[27] R.J. Vanderbei, Linear Programming Foundations and Extensions, International Series in Operations Research & Management Science, 2008. [28] M. Diaby, A linear programming formulation of the traveling salesman problem, in: Proceedings of the 11th WSEAS International Conference on Applied Mathematics, Dallas, Texas, USA, March 22–24, 2007.
67 68 69 70 71 72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65
131
66
132