Multi-agent Orbit Design for Visual Perception Enhancement Purpose

Abstract

This paper develops a robust optimization-based method to design orbits on which the sensory perception of the desired physical quantities are maximized. It also demonstrates how to incorporate various constraints imposed by many spacecraft missions, such as collision avoidance, co-orbital configuration, altitude and frozen orbit constraints along with Sun-synchronous orbit constraints. The paper specifically investigates designing orbits for constrained visual sensor planning applications as its case study. For this purpose, the key elements to form an image in such vision systems are considered and effective factors are taken into account to define a metric for perception quality. The method employs a max-min model to ensure robustness against possible perturbations and model uncertainties. While fulfilling the mission requirements, the algorithm devises orbits on which a higher level collective observation quality for the desired sides of the targets is available. The simulation results confirm the effectiveness of the proposed method for several scenarios involving low and medium Earth orbits as well as a challenging space-based space surveillance program application.

1. Introduction

In many practical applications, there is an increasing tendency to use a network of inexpensive agents with more sophisticated data-gathering equipment. This involves exploiting the coordinated operation of the agents to achieve more accurate measurements, and also makes the mission more robust to possible errors and faults. Challenging problems in multi-sensor management systems arise when the sensors operate in a dynamic environment, and they should be repeatedly moved to perceive the desired physical quantities of non-stationary targets. In such a multi-agent system, in order to leverage resources more efficiently, a meticulous plan of action is required to distribute tasks among the agents [1].

In this paper, a special class of multi-agent systems to acquire sensory measurements from physical phenomena of interest is considered. Many orbit design problems belong to this class of multi-agent systems. In fact, one scenario that frequently happens in various space missions is the design of satellite trajectories on which the perception qualities of certain physical quantities of some targets are maximized. In this constrained path planning problem, the orbit should also fulfill particular requirements imposed by the mission. Depending upon the system configuration and requirements, the design problem could be quite elaborate, even for a small number of agents and targets.

There is an extensive literature on the subject of orbit design for different science missions specifically oriented towards special types of orbits [2–5]. In these studies, the main aim is to define the attributes of the mission, usually without providing information about multi-agent path planning. In another interesting field of study, researchers have addressed a coordinated trajectory design for multi-agent, multi-target systems in the context of formation flying. The formation flying research field seeks to develop strategies to control the formation of a satellite fleet in order to meet a global performance objective, and the focus of interest is to employ an adequate control strategy to maintain the desired formation [6–10]. In many applications, the ultimate goal is to effectively allocate each sensor to a target at any point in time [11, 12]. However, before performing the resource allocation process, there are fundamental questions that need to be addressed early in the planning stage. The first and most important concerns how the agents (i.e., the leaders in a cluster formation) should be placed or moved with respect to the targets during the operation to optimize that particular mission [13]. In other words, prior to the allocation problem, the agent trajectories that maximize the sensory perception of the desired physical quantities should be determined by the planning scheme.

In this article, the multi-agent orbit design procedure is formalized as a constrained parametric optimization problem by taking into account different constraints that might be involved in a science mission. In particular, the Sun-synchronous orbit (SSO) and frozen orbit constraints that are widely employed in practical application are defined in terms of optimization parameters. The developed method is not limited to any particular sensor allocation scheme, and it increases the overall system performance since a higher level of perception quality of the targets is achievable on the designed orbits. Nourzadeh and McInroy [11] specifically investigate the sensor management problem and propose a non-myopic (multi-step) scheduling algorithm to address the multi-agent, multi-target sensor allocation problem.

To analyse the performance of the proposed method, a special spaced-based space situational awareness application is chosen. In this application, the main objective is to determine agents’ trajectories on which the observation quality of some resident space objects (RSOs) is maximized. Although the orbit design procedure is performed in the context of inspection quality, the proposed formulation is quite general, and any type of physical quantity can be incorporated and optimized. To deal with large-scale cases, some guidelines are proposed to choose an adequate numerical optimization algorithm.

The contributions of this paper are: (1) The proposed formulation allows us to solve more general multi-agent orbit design problems seeking to optimize the sensory perception of multiple attributes of targets. This extends the existing single scenario optimization schemes. (2) It also provides a way of incorporating many prevailing configurations in spacecraft missions, such as collision avoidance, co-orbital configuration, altitude and frozen orbit constraints as well as the SSO. The latter extends Boain's scheme [2] of SSO mission design to the multi-agent case. (3) To handle possible perturbations and uncertainties in the system, a max-min model is employed to bring robustness into the planning. (4) A perception model for the vision system is derived by considering key contributing factors affecting the quality of the images which encompasses a wide range of circumstances occurring in space situational awareness applications. The model is particularly useful for spaced-based SSA (SBSSA) applications, e.g., a spaced-based space surveillance program whose main goal is to build a network of satellites equipped with powerful vision systems in order to detect and track objects in Earth orbits.

This paper is organized as follows. Section 2 formulates the coordinated trajectory design problem that determine the paths on which the quality of the received information is maximized (in a sense). Section 3 discusses contributing elements in image formation and defines a quality metric for a space situational awareness (SSA) application by considering the dominating factors. Section 4 describes the motion equations of the objects and explains two sets of parameters commonly used to uniquely identify an orbit. In Section 5, the constrained multi-agent, multi-target trajectory optimization problem is constructed, and a description is given of how to incorporate the constraints of different configurations that might be employed in a science mission. Brief remarks on choosing suitable numerical algorithms to solve the optimization problem are given in Section 6. Simulation studies are included in Section 7, and conclusions are drawn in the final section.

2. Problem Formulation

This section focuses on the problem of designing the best agent trajectories on which sensory perception is maximized subject to some constraints imposed by the system's nature. Consider a general path planning problem as illustrated in Figure 1. As depicted, two cube-shaped objects are moving on two different curved paths (AB and CD). The objects passively or else actively send out certain types of signals that can be perceived at any point of the space with different quality levels. At a particular point in time, depending on the system configuration, specific quality fields are formed for any perceptive physical aspects of the objects. In fact, when objects move on their paths, a time-dependent vector or scalar field is formed for any of these physical quantities.

Figure 1.

Path-planning to maximize sensory reception

In the case of a scalar field, given the state of a target i(ⁱX_R) and of a sensor j(^jX_S), and effective physical factors, the total achievable quality along the sensor trajectories (e.g., EF and GH) in the space is the sum of values of the perceptions at all points on the curve. It can be expressed by Equation 1:

$$\begin{array}{rrr}\hfill & \hfill & \hfill i=\mathrm{1,}\cdots \mathrm{,}n\\ \hfill {\displaystyle \underset{t_1}{\overset{t_2}{\int }}}{\mathrm{\Phi }}_l\left({}^{i}X{}_R(t)\mathrm{,}{}^{j}X{}_S(t)\right)dt\mathrm{,}& \hfill & \hfill j=\mathrm{1,}\cdots \mathrm{,}m\\ \hfill & \hfill & \hfill l=\mathrm{1,}\cdots \mathrm{,}L\end{array}$$

(1)

where Φ _l is a function that describes the perception value of quantity l sensed by sensor j from object i at any point in time. The symbols m,n and L denote the number of agents, targets and physical quantities, respectively.

In the discrete case, the perception quality of target i obtained by all the sensors can be written as Equation (2). For i = 1, …, n and l = 1, …, L, (2) gives a perception quality matrix J_sum ∊ 葷^{n × L}. The (i,l) element of J_sum gives the l^th quality attribute of the i_th object along the trajectories:

$$J_{sum}={\displaystyle \sum _{j=1}^m}{\displaystyle \sum _{k=1}^N}{\mathrm{\Phi }}_l\left({}^i{X}_R[k{]\mathrm{,}}^j{X}_S[k]\right)\mathrm{.}$$

(2)

In this paper, it is presumed that each trajectory can be parametrized with a set of parameters. Concretely, ^j X_S moves on the trajectory j that is parametrized by p_j ∊ 葷 ^np_j where np_j is the number of parameters required to represent trajectory j. Let $p\in {葷}^{{\displaystyle \sum _{j=1}^m}n{p}_j}$ be a vector consisting of all the requisite parameters to describe the sensors’ trajectories (Equation (3)):

$$p=\underset{j=1}{\stackrel{m}{\oplus }}{p}_j\mathrm{.}$$

(3)

Here, ⊕ is the vertical concatenation symbol.

It is worthwhile noting that the relation between the parameters p_j and the trajectory j does not have to be explicit, and even a heuristic path planning algorithm that generates a trajectory j based on a tuning parameter p_j can be utilized here.

Depending upon the application and how to maximize the J_sum matrix, various possible optimization strategies can be considered to acquire different sets of sensors’ trajectories. In this article, to deal with worst cases, a max-min optimization model is employed to design trajectories which enforce planning redundancy, and make the result robust to the probable uncertainty sources in the system. Therefore, the optimization formulation determines the sensors’ trajectories that maximize the minimum achievable perception for all objects (the minimum element of J_sum) subject to some constraints imposed by the system (Equation (4)). Figure 1 portrays the case when each agent's trajectory is confined to planes z = 0 and x = 0:

$$\{\begin{array}{ll}\underset{p}{\text{Maximize}}\hfill & min(J_{sum})\hfill \\ \text{subject to}\hfill & :\text{Constraints}\hfill \end{array}$$

(4)

As mentioned above, two prime requirements to construct the optimization problem 4 are the perception function, Φ, and the trajectories parameter vector, p. In the next sections, these elements are defined for a visual measurement SSA application. In this scenario, it is assumed that all the agents are equipped with passive visual sensors and that the physical quantities that are maximized along the trajectories are inspection qualities of the different sides of the objects. Furthermore, each trajectory is parametrized by orbital elements parameters. While the following section elaborates the calculation process of the observation quality given the states of the objects and the cameras in this scenario, Section 4 describes the motion dynamics of agents and targets as well as the trajectories parametrization process.

3. Optical Imaging Quality Calculation

There are many vision systems with different sensors and means of data acquisition, and the process of sensor planning is highly dependent on the way in which information is gathered in that specific application. In some vision systems, only passive sensors are used while both passive and active sensors are used in some other machine vision systems. For simplicity, suppose that all the vision systems in our application are equipped with passive sensors (although the methods developed here also apply to active sensing). To compute the qualities of the resulting images, all the effective factors for forming an image should be considered, i.e., the intrinsic parameters of the sensors and environmental factors. There are a lot of techniques for automatically determining image quality [14–18].

In this section, the preliminary requirements to compute the observation quality in an optical imaging system will be demonstrated in detail by investigating all the effective aspects in forming images and defining a quality metric based upon them.

3.1 Effective Factors on Image Quality

The output of an optical imaging system depends on several physical factors, i.e., lighting, atmospheric attenuation, light diffraction, occlusion, imaging sensor resolution and sensitivity, electronics parts and output devices (Figure 2). The modulation transfer function (MTF) and the contrast transfer function (CTF) are two prevalent techniques for assessing the performance of a vision system. In fact, the MTF enables the quantifying of the resolving ability of each component of the system in different spatial frequencies. The cut-off frequency of MTFs can be used to determine the performance of each component [19, 20].

Figure 2.

Contributing elements on the image formation and quality

In SSA applications, only some physical aspects have a significant influence on the output image. As such, in subsequent sections, only the illumination, light fraction, occlusion and side observation quality will be taken into account.

3.1.1 Light diffraction

In optics, diffraction is categorized into two different classes, i.e., Fraunhofer diffraction (or far-field diffraction) and Fresnel diffraction (or near-field diffraction). Fraunhofer diffraction occurs when waves of a far-field distance are passed through an aperture or slit. This happens when the Fresnel number F « 1 and the parallel rays approximation is applicable.

Because of the Fraunhofer diffraction phenomenon, even an optical lens with perfect quality has limited performance and the MTF of each lens expresses this limitation in a spatial frequency. In fact, light coming from a point light source diffracts and forms an airy disk pattern. Using this pattern, Rayleigh proposed a criterion to find an optical resolution of the lens. The criterion states that two points with a angular separation equal to the angular radius of the airy disk can be resolved. Equation 5 relates the angular resolution θ to a wavelength of light (Λ) and the diameter of a lens aperture D [21]:

$$\sin \theta =1.22\frac{\lambda }{D}$$

(5)

The image-resolving quality for a system with a fixed aperture diameter D and a fixed wavelength Λ can be expressed by Equation (6), where δ is the distance between the object and the lens.

$$q_{res}=\frac{D}{1.22\lambda \delta }$$

(6)

Thus, by increasing the object distance, the resolving ability decreases because of the diffraction.

3.1.2 Illumination

The relationship between the reflected light (radiance) and the incoming illumination depends on the direction from which the light arrives as well as the shape and type of surface. Absorption, transmission, scattering or a combination of these effects applies to the incoming light when it strikes the surface. In other words, the intensity and the colour reflected depends on the illumination and reflection angle and the surface material. Surfaces can be categorized into three different groups, i.e., diffuse, Lambertian and specular surfaces. Specular surfaces behave like a mirror and reflect light into a lobe of the specular direction, and so the reflected light is highly dependent on the illumination direction, while for Lambertian and diffuse surfaces such as cotton cloth, matte paper and matte paint, the radiance leaving the surface does not have any meaningful correlation with the illumination direction, and the bidirectional reflectance distribution function (BRDF) of these surfaces is constant [22]. Since the outer surfaces of man-made satellites are highly reflective, the angle of the light with respect to the viewing axis is extremely critical. Let α _lit be the angle between the viewing axis and incoming sunlight (Figure 4). When α _lit ≃ 0, the image contains much specular glare. On the other hand, if α _lit > 90°, the satellite will be back-lit by the Sun, meaning that it will produce a silhouette with fewer details. α = 45° will provide the best illumination level. A suitable function that describes the illumination effect and which lies between zero and one can be written as Equation (7), where q_lit is the minimum value that the illumination quality takes. Figure 3 depicts the values which q_lit takes on as α _lit changes from zero to 360 degrees in polar coordinates when q_lit = 0.2:

$$q_{lit}=\{\begin{array}{ll}\underset{\_}{q_{lit}}+\left(1-\underset{\_}{q_{lit}}\right){\sin }^2\left(2{\alpha }_{lit}\right)\hfill & {\alpha }_{lit}<\frac{\pi }{2}\hfill \\ \underset{\_}{q_{lit}}\hfill & {\alpha }_{lit}\ge \frac{\pi }{2}\hfill \end{array}$$

(7)

Figure 3.

Illumination quality of a specular surface for different values of the observation angle

Figure 4.

Determining illumination quality in a SSA application

Figure 4 illustrates the illumination quality of an observation for three RSOs which can be observed by two observer satellites (O₁, O₂). Observer 2 inspects RSO 1, and the angle ${\alpha }_{lit}^{12k}$ , the angle between sunlight and the observation direction, is around 45°. Therefore, the illumination quality will be around 1. Observer 1 has two options to inspect: either RSO 1 or RSO 3. The image taken of RSO 3 is of better quality in terms of illumination ( ${\alpha }_{lit}^{21k}\approx 45°$ ), while the RSO 1 image will have a silhouette problem, since ${\alpha }_{lit}^{11k}>90°$ .

3.1.3 Occlusion and Sunlight Quality

Since an occluded object cannot be inspected by the observer, the resulting image quality for that object should be considered to be zero. On the other hand, if there is no obstacle between the object and the observer, the possibility of having a perfect image of that object exists (Equation 8):

$$\{\begin{array}{ll}{q}_{los}=0\hfill & \text{\hspace{1em}occlusion}\hfill \\ q_{los}=1\hfill & \text{\hspace{1em}no occlusion}\hfill \end{array}$$

(8)

In order to observe the object, an acceptable level of illumination is required. For instance, in satellite imaging, sunlit images are more desirable while dark images obtained from shadow areas behind the Earth are not as useful as sunlit images. To incorporate the sunlight quality factor, let q_lum = 1 when the object is in sunlight. Otherwise, q_lum is a small positive number(Equation 9):

$$\{\begin{array}{ll}{q}_{lum}=1\hfill & \text{\hspace{1em}lit image}\hfill \\ q_{lum}\ll 1\hfill & \text{\hspace{1em}dark image}\hfill \end{array}$$

(9)

3.1.4 Sides Observation Quality

In visual sensor planning, the vast majority of studies are mainly concerned with finding the vision system parameters to inspect every side of an object with a minimum number of observations [23]. One way to do this is to assign an outward facing unit normal for each side of interest, i.e., let n_ilk be the l^th outward facing unit normal on an object i at a sample time k. To determine whether a specific side of an object is viewed, let ν _ijk be the unit vector pointing from object i to observer j at sample time k. Thus, the observation quality of side l of object i can be computed as Equation (10):

$$q_{view}^{ilkl}=\{\begin{array}{ll}{n}_{ilk}^T{v}_{ijk}\hfill & n_{ilk}^T{v}_{ijk}>0\hfill \\ 0\hfill & n_{ilk}^T{v}_{ijk}\le 0\hfill \end{array}$$

(10)

Figure 5 depicts a RSO with some assigned sides of interest (D_l). According to Equation (10), (D₁) and (D₂) have a positive quality and the inner product of υ and the other side normals are less and equal to zero. This means that an observer in this direction can partially inspect two sides.

Figure 5.

Assigned normals to specify the sides of interest and the computation of the sides observation quality (inner product)

3.2 Single Quality Metric

In order to form a single number reflecting the overall observation quality of an object i inspected by an observer j during a sample k along a face l, all the calculated qualities should be combined. Since the overall quality varies from application to application, this can be done by multiplying weighted qualities. This method allows us to control the participation of all the constituent parts in a single quality metric. Let q_ijkl be the product of all the mentioned qualities (Equation (11)):

$${\mathrm{\Phi }}_l\left({}^i{X}_R[k{]\mathrm{,}}^j{X}_j[k]\right)=q_{ijkl}=q_{res}^{ijk}{q}_{los}^{ijk}{q}_{lum}^{ijk}{q}_{angle}^{ijk}{q}_{view}^{ijkl}$$

(11)

For different values of i, j, k and l, q_ijkl forms a 4D observation quality array $Q\in {葷}_{+}{}^{n\times m\times N\times L}$ .

4. Orbital Elements and Motion Equations

4.1 Motion Dynamics

Assuming a uniform, spherical Earth, Equation (12) represents the motion equation of the target/object on its orbit:

$$\ddot{\overrightarrow{r}}=-\frac{\mu }{r^3}\overrightarrow{r}+\overrightarrow{d}$$

(12)

where $\overrightarrow{r}=[xyz{]}^T$ is the relative position of the target with respect to the centre of the Earth. The vector $\overrightarrow{d}=[d_x{d}_y{d}_z{]}^T$ denotes exogenous forces, including the attitude control system (ACS) actions and modelled or unmodelled disturbing forces affecting the target, such as solar pressure, atmospheric drag, etc. The Earth gravitational constant $\mu \triangleq 398600.4418{\text{Km}}^3{\mathrm{s}}^{-2}$ , and $r=\parallel \overrightarrow{r}{\parallel }_2$ , together represent the Euclidean distance of vector $\overrightarrow{r}$ .

Letting the state vector $X=[\overrightarrow{r}\dot{\overrightarrow{r}}{]}^T$ and the velocity vector $\overrightarrow{v}=\dot{\overrightarrow{r}}=[\dot{x}\dot{y}\dot{z}{]}^T$ , Equation (12) can be rewritten as a set of first-order ordinary differential equations (ODEs) consisting of only the orbital state vector (Equation (13)):

$$\dot{X}=f(X\mathrm{,}\overrightarrow{d})=\left[\begin{array}{c}\overrightarrow{r}\\ -\mu \frac{\overrightarrow{r}}{r^3}+\overrightarrow{d}\end{array}\right]$$

(13)

The discretized motion equation can be numerically solved by choosing an adequate ODE solver (Equation (14)):

$$X\left[k+1\right]=F\left(X[k]\mathrm{,}\overrightarrow{d}[k]\right)\mathrm{,}X[k]\approxeq X\left(k{T}_s\right)\mathrm{.}$$

(14)

where T_s denotes the sampling rate.

4.2 Orbital Elements

When the ACS actions and perturbing forces cancel each other (i.e., perfect control), the $\overrightarrow{d}$ will be a zero vector and the target trajectory forms an elliptic orbit. This orbit can be uniquely determined by knowing the orbital elements at an epoch (time t). There are several sets of orbital parameters that result in the same motion dynamics. In classical two-body systems, Keplerian orbital elements and the orbital state vector X comprise two commonly used sets of orbital elements. The orbital state vector contains two physical quantities, the position and velocity of the moving target at an epoch written in the reference frame.

On the other hand, the Keplerian orbital elements are directly related to the shape of the resultant orbit and the position of the target on the path. They consist of six parameters: the semi-major axis (a) is the sum of the periapsis (perigee) and apoapsis (apogee) distances divided by two; eccentricity (e) controls the shape of the elliptic orbit and specifies by how much the orbit deviates from a circular orbit; the inclination ι is the angle between the equatorial plane and the plane containing the elliptical orbit (elliptical plane); the longitude of the ascending node (Ω) is the angle between the ascending node vector and the direction of the vernal equinox; the argument of perigee (ω) orients the elliptical orbit in the elliptical plane, and it is the angle between the right ascending node direction and the Perigee's direction; and finally the true anomaly (ν) at epoch t represents the geometric angle of the target in the orbital plane (Figure 6).

Figure 6.

Keplerian and orbital state elements

Since an orbital state vector can be converted into Keplerian orbital parameters by a nonlinear transformation at any epoch [24], throughout this paper, any of these representations might be employed wherever it is more convenient to do so.

5. Optimization Formulation

The problem of maximizing the minimum observation quality for all sides of the targets with respect to the sensor's trajectories parameters can be written as Equation (15):

$$\{\begin{array}{cc}\underset{p}{\mathrm{maximize}}& f(p)\\ \mathrm{subject}\mathrm{to}& g_i(p)\le 0,i=1,\cdots ,n_g\\ & h_j(p)=0,j=1,\cdots ,n_h\end{array}$$

(15)

where n_g and n_h refer to the number of inequality and equality constraints, respectively:

$$p=\underset{j=1}{\overset{m}{\oplus }}{p}_j=\underset{j=1}{\overset{m}{\oplus }}{[a_j,e_j,{\Omega }_j,{\omega }_j,{\iota }_j,v_j]}^T.$$

(16)

The subscript j denotes the number assigned to the agent that moves in orbit with the orbital parameter p_j, and f(p) = min(J_sum(p)).

5.1 Co-orbital Configuration

In many surveillance missions, some agents will be moving on the same orbit. In this case, all the orbital elements of these agents are the same except for ν, which can be freely chosen. To enforce this constraint for agents j₁ and j₂ sharing the same orbit, the following linear equality constraint should be added to the optimization problem:

$$A^{T}p=0$$

(17)

where A^T is a 5 × 6m matrix consisting of m horizontally concatenated 5 × 6 block matrices. All of these sub-matrices are a zero matrix except for j^th₁ and j^th₂, which are replaced by A_s and – A_s, respectively. Mathematically, matrix A can be constructed by the following equations:

$$A=\underset{j=1}{\overset{m}{\oplus }}\{\begin{array}{ccc}{A}_s^T& j=j_1& \\ -A_s^T& j=j_2& ,A_s=\left[I_{5\times 5},0_{5\times 1}\right]\\ 0_{6\times 5}& \mathrm{otherwise}& \end{array}.$$

5.2 Uniform Distribution

Another configuration that might be considered in orbit design would be to uniformly distribute those agents which have a common orbit (trailing formation). Let O be the number of distinct observers’ orbits and let m_o denote the number of satellites in the orbit o. In this case, the true anomaly difference of the initial states of two adjacent satellites in orbit o should be ${\frac{360}{m_o}}^{\circ }$ . It is worthwhile to mention that only in an unperturbed circular orbit will the difference be maintained during motion. Figure 7 illustrates one possible configuration when three observer satellites are evenly distributed with a constant true anomaly distance.

Figure 7.

Spreading satellites in an orbit

To incorporate this scheme in the optimization problem, consider a subset of agents’ indices that share the same trajectory (orbit o) {j₁, j₂, …, j_mo}. Assuming different types of satellites with different characteristics, there will be (m_o – 1)! ways to arrange these m_o distinct objects in orbit o. To find the optimal arrangement, one should solve the optimization problem for all possible orders and choose the one that has the best resultant cost function. This method is impractical when dealing with a large-scale system, and to avoid this cumbersome process, an approximate approach should be taken into account. One heuristic way that results in satisfactory suboptimal solutions can be obtained by assuming identical agents and solving the corresponding optimization problem. Next, check the performance of all possible arrangements with the resulting trajectories and choose the best configuration.

Assuming identical agents, the true anomaly of the first agent (ν) in orbit o(agent j₁) should vary between 0 and ${\frac{360}{m_o}}^{\circ }$ . In addition, the following constraints should be added to the optimization problem for adjacent agents in orbit o (Figure 7).

$$\{\begin{array}{l}{\nu }_{j_2}-{\nu }_{j_1}=\frac{360°}{m_o}\hfill \\ {\nu }_{j_3}-{\nu }_{j_2}=\frac{360°}{m_o}\hfill \\ ⋮\hfill \\ {\nu }_{j_{m_o}}-{\nu }_{j_{m_o-1}}=\frac{360°}{m_o}\hfill \end{array}$$

(18)

5.3 Collision Avoidance Constraint

To avoid collisions between agents and moving objects around the Earth, the optimization problem should guarantee that all the sensors retain a minimum safe distance from other moving objects during the mission. This can be incorporated in the optimization problem (Equation (15)) by heavily penalizing the cost function whenever the minimum distance (d) between agents and any orbiting objects in all time instances is less than an acceptable threshold tr (Equation (19)):

$$f(p)=\{\begin{array}{ll}{J}_{sum}(p)\hfill & d\ge tr\hfill \\ 0\hfill & d<tr\hfill \end{array}$$

(19)

5.4 Altitude Constraints

Satellite altitude is an important parameter that relates to many mission requirements. In the mission design process, many factors, such as launch constraints, type of satellite, desired perception resolution, mission lifetime, ground track repeatability and mission expenses directly affect the altitude of the satellite.

For instance, regardless of the fuel required to retain orbital speeds, depending upon the material technology, a satellite can withstand a limited amount of heating caused by atmospheric resistance. Therefore, for any SSA mission, there is a minimum operational altitude associated with a satellite. This restriction can be included in the trajectories optimization problem by adding the following nonlinear constraint:

$$a_j\left(1-e_j\right)>R_E+{\underset{\_}{c}}_j\mathrm{,}j=\mathrm{1,}\cdots \mathrm{,}m$$

(20)

where R_E = 6371 Km is the Earth's radius and ${\underset{\_}{c}}_j$ denotes the minimum orbital altitude for an agent j.

Depending upon the altitude, an orbiting satellite receives different kinds of radiation (such as cosmic rays, van Allen radiation, solar flares, etc.) which can cause damage to the satellite equipment. To avoid such hazardous exposure, the mission designer can choose different strategies based on any constraints and key requirements. For a satellite in low Earth orbit (LEO), one way to stay away from a region of intense radiation is by keeping the altitude of the satellite lower than the Allen radiation inner belt altitude, which is around 1,000 Km. Another way is to use expensive radiation-tolerant components in the satellite building process to make sure that the equipment operates properly throughout the mission.

Therefore, for a space mission, an upper bound associated with the altitude might be defined. The upper bound constraint can be expressed by Equation (21):

$$a_j\left(1-e_j\right)<R_E+{\overline{c}}_j\mathrm{,}j=\mathrm{1,}\cdots \mathrm{,}m$$

(21)

Here, ${\overline{c}}_j$ denotes the agent j's maximum orbital altitude.

5.5 Specify Orbit Type

A satellite is classified into various categories according to its orbital altitude, inclination, eccentricity and period, etc. Practically every constraint regarding the shape of the orbit expressed by orbital parameters can be included in the proposed trajectory optimization problem. For instance, to design an orbit for agent j within LEO, its altitude should vary within the 0–2,000 km range, i.e., (0 < a_j < R_E + 2000). Similarly, the motion of an agent j moving in a polar orbit can be planned by adding the equality constraint ι _j = 90° to the problem.

5.6 Sun-synchronous Orbit

In science missions, one of the most widely-used types of orbits is the SSO. The SSO is a near-polar and almost-circular geocentric orbit whose nodal precession rate $(\dot{\mathrm{\Omega }})$ is equal to the Earth's mean rotation rate around the Sun. Geometrically, a SSO approximately orients in such a way that the angle between the orbital plane and the vector from the Sun to the Earth remains the same during the mission. Therefore, the illumination angle of the ground track will be constantly maintained throughout the mission. Besides this interesting characteristic, a SSO has other orbital properties that make it highly desirable for various applications [2]. In the subsequent subsections, the SSO orbital parameters’ constraints imposed by a common scientific mission are briefly discussed, and it is shown how the mission requirements greatly restrict the feasible space of the optimization parameters.

5.6.1 Precession Rate Constraint

Due to the out-of-plane gravitational force caused by the Earth's equatorial bulge, an orbit plane gyroscopically precesses. The corresponding nodal precession rate $\dot{\mathrm{\Omega }}$ can be operatively computed by Equation (22). To have a SSO, $\dot{\mathrm{\Omega }}$ should be equal to the Earth's mean orbital rate around the Sun, i.e., $\dot{\mathrm{\Omega }}=\frac{365}{365.242199}degree/day=1.991063802746144\times {10}^{-7}rad/sec$ :

$$\dot{\mathrm{\Omega }}=-\frac{3}{2}{\tilde{J}}_2{\left(\frac{R_E}{a\left(1-e^2\right)}\right)}^2\sqrt{\frac{\mu }{a^3}}\cos (\iota )$$

(22)

where ${\tilde{J}}_2=1.08263\times {10}^{-3}$ denotes the Earth's dimensionless zonal harmonic coefficient.

Figure 8 illustrates the surface on which the precession rate condition is satisfied for a SSO within the operational region of a LEO where R_E + 300Km ≤ a ≤ R_E + 2000Km. Although every point on this surface meets the precession rate condition, a relatively small portion of the surface contains the admissible orbits’ parameters - for instance, most points on the surface do not meet the minimum altitude requirement. In Figure 8, the area surrounded by black lines contains the points satisfying the minimum altitude constraint a(1 – e) > R_E + 300Km.

Figure 8.

Surface of all the SSO within a LEO that satisfy the precession rate condition

5.6.2 Frozen Orbit Constraint

Owing to the perturbing forces caused by the oblateness of the Earth, the satellite nominal trajectory based on the selected orbital elements does not remain fixed, instead varying as a function of time. In a SSO orbit, this fact is exploited to achieve the gyroscopic precession of the orbital plane. On the other hand, perturbing forces also adversely affect the eccentricity, e, and argument of the perigee, ω, of the orbit. Perturbation theory states that systematic choices of the orbital parameters can minimize the drift from the selected initial values. In a frozen SSO orbit, the parameters e, ι and ω are picked in such a way that the secular perturbations of J₂ and J₃ cancel out, and parameters only undergo a relatively small periodic perturbation with a period equal to the orbit period (Equation (23)). Theoretically, a satellite in such an orbit has a minimum propellant usage requirement during a long-term mission:

$$\{\begin{array}{l}\tilde{\omega }=90\text{ or }270\text{ degrees}\hfill \\ \tilde{e}=-\frac{J_2}{J_3}\frac{\sin (\iota )}{2a\left(1-e^2\right)}\hfill \end{array}$$

(23)

where $\tilde{\omega }$ and $\tilde{e}$ are the mean value of the desired frozen orbit parameters.

Figure 9 depicts the locus of the all the SSO parameters that satisfy the frozen orbit constraint. As highlighted, a very small part of the precession rate surface results in the minimum perturbed orbital parameters set. This has a direct practical implication for the optimization algorithm performance by drastically decreasing the size of the feasible sets.

Figure 9.

The locus of the frozen SSO parameters within the LEO

6. Numerical Optimization Algorithm

This section considers the problem of choosing an adequate optimization strategy to solve the trajectory design problem governed by Equation (15). The optimization algorithm selection is the most important part of solving a practical optimization problem. This process is performed by considering different factors that can affect the overall performance. The performance of an algorithm is commonly assessed according to its convergence speed, solution optimality, robustness to perturbation and uncertainty, required resources, computational complexity and implementation difficulty, etc. Quantifying all the effective factors is extremely difficult and occasionally impossible for an optimization method without any knowledge about the problem itself.

Equation (15) is a deterministic constrained nonlinear optimization problem with continuous variables. For such a large-scale system with several targets and agents, the idea of utilizing global optimization methods seems contrived. In fact, current known global approaches offer successful performance in small- and medium-sized problems, as the solution space size exponentially increases with respect to the number of decision variables. Therefore, in this paper, a local optimization algorithm is employed to tackle (15).

In practice, variations of sequential quadratic programming (SQP) and the interior point method (IPM) are commonly-used algorithms to solve large-scale general constrained optimization problems. In each iteration, a SQP method generates admissible steps towards a local minimum by solving a quadratic model of the objective function subject to a linearized version of the problem constraints. SQP-type algorithms have been quite successful, especially in dealing with nonlinear constraints, and has shown more robustness in relation to the badly-scaled optimization problem than IPM algorithms.

Like SQPs, IPM-type algorithms are also well-known for their superior performance in solving nonlinear constrained problems. IPM algorithm's promising performance in solving linear problems has motivated experts to utilize its key ideas, i.e., primal-dual steps, to devise powerful nonlinear optimization algorithms. IPMs often outperform SQPs in dealing with large-scale applications, especially when the system has a block-structure and sparse representation. On the other hand, IPMs have shown a lack of robustness in relation to initial point selection, problem scaling and barrier parameters. Recently, successful attempts have been made to develop more robust IPM-type algorithms [25].

In this paper, since the trajectory design problem formulation is block-structured in terms of optimization parameters, to treat large-scale optimization problems, IPM-type solvers are employed, and for small-scale problems, SQP methods are used.

7. Simulations

To overcome obstacles in ground-based SSA in detecting and tracking space objects, SBSSA has begun to build a network of satellites equipped with different sensors. To fulfil this aim, on 26 September 2010, Minotaur IV performed its first orbital launch of the SBSS system. The main objective of the SBSS program is to search, detect and track objects in Earth orbits, especially geosynchronous orbit (GEO) objects. The SBSS program is about to build a larger constellation of observer satellites to cover wider areas of space. With the development of technology, it is expected that there will be more SBSSA missions in the near future. The simulation part of this paper is intended to exploit the proposed multi-agent trajectory design method to plan agents’ movement in the context of SBSSA.

This section consists of two parts. In the first subsection, an illustrative simulation study for a fairly small-scale system with five RSOs and two agents is presented. In this system, all the targets and agents are orbiting within a LEO and in the equatorial plane. In section 7.2, multi-agent orbit design is performed for more general scenarios involving unperturbed a SSO.

To perform the numerical integration of the motion dynamics of all the objects (Equation (13)), the fourth-order Runge-Kutta (RK4) method is employed with an adequate step size. Since the chosen method to solve the trajectory optimization problem results in the largest function value in some feasible neighbourhood (i.e., the local optimum solution), to achieve a decent sub-optimal solution for each problem the optimization algorithm will be run from different initial parameters.

7.1 Illustrative Case Study

In this case study, five RSOs are moving in two elliptic orbits in the equatorial plane with the following Keplerian elements at time t = 0. The objective is to plan two agents’ motion so that the visual information gathered from two faces of the targets is maximized. The simulation duration is equal to five-times the period of the slowest orbit $(t_2-t_1=5\times 2\pi \sqrt{a^3/\mu }=35830.7Sec)$ , and the sunlight direction is assumed to be from the positive sides of the x-axis to its negative side, i.e., [−1 0 0] ^T . The threshold for the collision avoidance constraint is considered as tr = 3Km.

	a	e	ι	Ω	ω	ν
Orbit 1	8033	0.126	0	0	68.02	247, 62, 237
Orbit 2	7898	0.057	0	0	225.2	281, 271

Figure 10 depicts how the observation quality field changes in five time instances as the system evolves. As shown, the quality field varies in a highly nonlinear manner. When the Sun shines on one side of the Earth, it casts a shadow on the opposite side. To provide a basis for comparison, two different scenarios are considered. In the first situation, it is assumed that the agents are co-orbital and uniformly distributed on an equatorial orbit, while in the second case, two different equatorial orbits are designed. For both cases, 12 unknown variables should be determined by solving the optimization problem, i.e., $p=[a_1\mathrm{,}{e}_1\mathrm{,}{\mathrm{\Omega }}_1\mathrm{,}{\omega }_1\mathrm{,}{\iota }_1\mathrm{,}{\nu }_1\mathrm{,}{a}_2\mathrm{,}{e}_2\mathrm{,}{\mathrm{\Omega }}_2\mathrm{,}{\omega }_2\mathrm{,}{\iota }_2\mathrm{,}{\nu }_2{]}^T$ . The p vector will be limited to the upper and lower $\overline{p}$ and $\underset{\_}{p}$ vectors in Equation 24. Since the systems are small-scale in these scenarios, a SQP-type algorithm is employed to numerically solve the corresponding optimization problem.

Figure 10.

Observation quality field evolution in time

$$\begin{array}{l}\overline{p}={\left[\mathrm{8033,0.15,360,360,10,360,8033,0.15,360,360,10,360}\right]}^T\hfill \\ \underset{\_}{p}={\left[\mathrm{7898.3,0,0,0,0,0,7898.3,0,0,0,0,0}\right]}^T\hfill \end{array}$$

(24)

7.1.1 Co-orbital Case

Equation (25) is the nonlinear optimization problem corresponding to the first scenario.

$$\{\begin{array}{ll}\underset{p}{\text{maximize}}\hfill & min\left(f(p)\right)\hfill \\ \text{subject to:}\hfill & \hfill \\ \hfill & \underset{\_}{p}\le p\le \overline{p}\hfill \\ \hfill & a_j\left(1-e_j\right)>R_E+\mathrm{300,}j=\mathrm{1,2}\hfill \\ \hfill & \begin{array}{l}{\mathrm{\Omega }}_1=0\to \left[0_{1\times 2}10_{1\times 11}\right]p=0\hfill \\ {\mathrm{\Omega }}_2=0\to \left[0_{1\times 8}10_{1\times 3}\right]p=0\hfill \end{array}\hfill \\ \hfill & \begin{array}{l}{\iota }_1=0\to \left[0_{1\times 4}10_{1\times 7}\right]p=0\hfill \\ {\iota }_2=0\to \left[0_{1\times 10}10\right]p=0\hfill \end{array}\hfill \\ \hfill & A^{T}p=\left[I_{5\times 5}{0}_{5\times 1}-I_{5\times 5}{0}_{5\times 1}\right]p{=0}_{5\times 1}\hfill \\ \hfill & \begin{array}{l}{\nu }_1\le 180\to \left[0_{5\times 1}10_{5\times 1}0\right]p\le 180\hfill \\ {\nu }_2-{\nu }_1=180\to \left[0_{5\times 1}-10_{5\times 1}1\right]p=180\hfill \end{array}\hfill \end{array}$$

(25)

Assuming the minimum operational altitude for both agents to be 300 Km, the altitude constraints are a_j(1 – e_j) > R_E + 300 for j = 1, 2. Equality constraints Ω₁ = Ω₂ = ι₁ = ι₂ = 0 force the orbits to be in the equatorial plane. Constraint A^Tp = 0 makes the two agents co-orbital, and constraints ν₁ ≤ 180 and ν₂ – ν₁ = 180 result in equally-spaced agents in the orbit:

7.1.2 Agents with Different Orbits

The second scenario formulation is like Equation (25) with only the six first constraints to make sure that the resultant orbits satisfy the boundary conditions, minimum altitude constraints and position in the equatorial plane. Both optimization problems are run from several feasible initial parameters to obtain satisfactory sub-optimal solutions. Figure 11 describes how the best trajectory achieved evolved from a feasible initial solution in the co-orbital case. As can be seen, the agents are uniformly distributed along the trajectories in all iterations.

Figure 11.

The best obtained suboptimal solution of Equation (25)

Figures 13 and 14 illustrate how the scaled optimization variables and the minimum observation quality of all the targets’ sides converge for each scenario. As highlighted, after convergence, the optimization cost function in the co-orbital case is 90.449, which is clearly less than the corresponding obtained observation quality along two other trajectories (146.826). This indicates that, in systems in which agents have more degrees of freedom, a higher perception level can be achieved.

Figure 12 depicts the sum of the quality pattern for each side of the targets, as well as the initial and optimized trajectories for both design scenarios. The graph shows how the initial trajectories are altered and developed to be well-positioned in areas with a higher perception level on both sides.

Figure 12.

Optimal trajectory and sides quality pattern of all RSOs for both design scenarios. (Co-orbital case (Top) and two distinct agents’ trajectories design case (Bottom)).

Figure 13.

Optimization progress in the co-orbital scenario

Figure 14.

Optimization progress in the scenario with two distinct orbits

7.2 Frozen SSO Design

In this part, the multi-agent orbit design procedure is performed for larger SBSSA systems. In the first part, because of the desirable characteristics of the SSO reported for the SBSS program [26], Sun-synchronous frozen orbits are found for a SBSSA system with five agents which must characterize 60 RSOs. In this scenario, 51 RSOs are travelling on a GEO and the rest are moving within a LEO (Figure 15).

Figure 15.

RSOs arrangement in a typical SBSSA application

In this design problem, the Keplerian parameters of the desired orbits should satisfy the SSO and frozen orbit constraints (Equation (22) and (23)). Furthermore, the lower and upper band for altitude are assumed to be ${\underset{\_}{c}}_j=300Km$ and ${\overline{c}}_j=1000Km$ , the orbital inclination lies in the [96.5, 102.5] interval and the argument of the perigee ω _j = 90. Thus, the optimization problem can be written as Equation (26):

$$\mathrm{.}\{\begin{array}{ll}\underset{p}{\text{maximize}}\hfill & min\left(f(p)\right)\hfill \\ \text{subject to:}\hfill & \hfill \\ \hfill & \underset{\_}{p}\le p\le \overline{p}\hfill \\ \hfill & a_j\left(1-e_j\right)>R_E+\mathrm{300,}j=\mathrm{1,}\cdots \mathrm{,5}\hfill \\ \hfill & a_j\left(1-e_j\right)<R_E+\mathrm{1000,}j=\mathrm{1,}\cdots \mathrm{,5}\hfill \\ \hfill & -\frac{3}{2}{\tilde{J}}_2{\left(\frac{R_E}{a_j(1-e_j^2)}\right)}^2\sqrt{\frac{\mu }{a_j^3}}\cos ({\iota }_j)=1.99\times {10}^{-7}\hfill \\ \hfill & j=\mathrm{1,}\cdots \mathrm{,5}\hfill \\ \hfill & e_j=-\frac{J_2}{J_3}\frac{\sin ({\iota }_j)}{2a_j\left(1-e_j^2\right)}\mathrm{,}j=\mathrm{1,}\cdots \mathrm{,5}\hfill \\ \hfill & {\omega }_j=\mathrm{90,}j=\mathrm{1,}\cdots \mathrm{,5}\hfill \end{array}$$

(26)

Figure 16 and 17 report the procedure to achieve a sup-optimal solution for the optimization problem in Equation (26). Figure 16 depicts 200 runs of the optimization from different initial parameters. The best minimum inspection quality for the desired sides of RSOs is 0.633, and Figure 17 indicates how the initial objective function is evolved in each optimization iteration. Figure 18 illustrates the designed trajectories on which the minimum perception quality is maximized, and it corresponds to the best sub-optimal solution acquired.

Figure 16.

The best obtained sub-optimal solution from 200 optimization runs for the SSO design case

Figure 17.

The cost function evolution of the best obtained solution in the SSO design scenario

Figure 18.

The optimized frozen SSOs for a multi-agent system

8. Conclusion

In this paper, a systematic means of optimizing trajectories for a sensor network system for perception enhancement purposes is presented. This has a major impact on the performance of the resource allocation stage in any sensor management system. Specifically, the proposed method is of interest for large-scale systems characterized by highly dynamic behaviour because of its robustness in relation to possible perturbations and model uncertainties. The case study results for SBSSA applications indicate the effectiveness of the method in improving the received perception qualities along the designed trajectories where the agents cooperatively meet the mission requirements, such that less effort is needed to serve the main purpose.