Optimal Sensor Planning for SSA Using System Identification Concepts

The rapid growth of man-made objects in Earth’s orbit has created a need for planning and prioritizing sensor measurements of these objects in order to keep track of their motion in orbit. In this paper, we adapt a concept from the field of system identification, namely application-oriented input design, to the sensor-planning problem. Previous methods have predominantly focused on the problem of planning the measurements of different objects such that they are as accurate as possible (in a broad sense) from a given set of possible observations. Here we view the problem differently; the objective is to find the least costly number of measurement such that the estimated orbital parameters are good enough for their intended application. To this end, we introduce the concept of an application set. The application set contains all orbital parameters that are considered of sufficient accuracy for the intended use. The objective of the optimization is hence not to find the most accurate set of parameters but find the cheapest combinations of observations such that the estimated parameters are guaranteed to be within the application set with a high probability. We show how to formulate the planning problem as a convex optimization problem that can be solved efficiently with modern algorithms, even for large-scale problems. We then demonstrate the feasibility of the method in a simulation example. Finally, the paper discusses some interesting topics for future research.


Introduction
The number of space objects in Earth's orbit is growing rapidly.This growth drives a need for improved space situational awareness, both in terms of increased sensor capacity and in terms of a more efficient use of sensor resources.The 47 Page 2 of 20 orbits of space objects are predominantly determined using data from networks of radar and optical sensors, but other techniques such as laser and passive radio ranging are used as well.Intergovernmental organizations, governments, and commercial companies continuously expand their sensor networks to cope with the rapid growth.This provides a necessary increase in sensor capacity, but it also adds further complexity to prioritization and optimization of the sensors network for efficient use.
Planning, tasking and scheduling of sensors for space situational awareness has been studied extensively, see e.g.[1][2][3][4].The problem is however not unique for space situational awareness and has been studied in many related areas.A good overview on the subject can be found in [5] and the references therein.A common approach is to try to pick an assignment of sensors that minimizes some cost.The cost is usually a function of information gain, that is, the inverse of the uncertainty in the estimated orbital parameters and includes covariance measures as well as more complex information measures.The tasking is then performed by distributing all possible observations times among the objects of interest that gives the most accurate parameters, by the selected measure.In Kalandros and Pao [6] the sensors are selected for observation, not to find the most accurate estimate but instead to form the parameter uncertainty matrix into a desired shape.
In this paper, we study the problem of optimal sensor planning, with the objective of finding an efficient method for the scheduling of space object measurements using a set of different sensors in a sensor network, in order to retain sufficient accuracy in the orbital parameters of the object.The network could consist of multiple sensors of different types, with different operational modes, and with different accuracy and performance.The main contribution of this paper is to apply the idea of application-oriented input design from the field of system identification to the sensor-planning problem.The field of system identification handles the problem of constructing a model of a dynamical system based on experimental data [7][8][9].
The concept of application-oriented input design for system identification was first introduced in [10].The idea is to find the cheapest possible identification experiment such that the identified model satisfies the performance specification given by the intended use of the model.The optimization is performed over the properties of the input signal such as the signal spectrum, power, or length.To quantify the required model accuracy an application set is introduced.The application set constrains the accuracy needed of the identified parameters by the application in which the model is used.The objective is hence not to find the most accurate parameters, but parameters that are sufficiently accurate for their intended use.
In the sensor-planning problem, instead of designing the input signal, we select at which time a certain sensor should be used to make an observation.Each observation is associated with a cost describing the effort, for example time or money, to perform the measurement.The objective is thus to find the optimal combination of possible observations, with the smallest total cost, such that the estimated orbital parameters are accurate enough for their intended application.The approach we take in this paper hence differs from the more common approach of trying to find the most accurate parameters (by some measure) utilizing all possible observations.The proposed method is sensor-agnostic and can be applied to planning problems in all orbital domains.
The outline of the rest of the paper is as follows: a brief introduction to orbit determination and covariance analysis is given in Sect.2, the sensor-planning problem is formulated in Sect.3, Sect. 4 introduces the concept of an application set and shows how to formulate the sensor-planning problem as a convex optimization problem; the feasibility of the proposed method is shown in simulations in Sect.5; finally, conclusions and directions for future research are given in Sect.6.

Orbit Determination and Covariance Analysis
In this section we give an introduction to the orbit determination and covariance analysis concepts required for the development of the planning algorithm.We will mainly follow the derivations and the notation from [11].Vectors are denoted by bold face lower case letters (e.g., x ) and matrices by bold face upper case letters (e.g., X).
Let denote a time-dependent parameter vector comprising of an object's position ( r ), velocity ( v ), as well as additional parameters related to the force and to the meas- urement models ( p and q , respectively).The parameter vector may equally well be parameterized differently, for example by using the classical orbital elements instead of the position and velocity.

The position and velocity evolve in time as
The objective is to estimate the parameter vector x(t) at a certain time t = t e from N observations at times t 1 , … , t N .Note that the measurements could come from different sensors.
For simplicity of notation, we here assume that the measurements z i (t i ) are sca- lars.However, the case with multiple measurement from a single sensor can easily be handled by adding multiple scalar measurements per sensor.An estimate of the vector x(t e ) can now be found using a least-squares method by minimizing We denote the solution to the least squares problem as x(t e ) .Details on how to find the solution can be found in e.g.[11] or [12].Under some assumptions on the measurements and the modelling of the dynamics, the estimate is Gaussian distributed with expected value This requires for example that there are more measurements than unknowns to estimate and that the measurements are rich enough, i.e., all free parameters in the model are observable in the data and that there exists a representation of the parameter vector, x 0 , describing the true underlying dynamics.If these requirements are sat- isfied the estimate will hence, on average, be correct.Furthermore, it can be shown ( [11]) that the (inverse) covariance matrix of the estimate is given by where Often there is some a priori information about the parameters, usually from a previous orbit determination.This information can be incorporated when solving the least squares problem [11].If the a priori information and its uncertainty is x apr (t e ) and P apr , respectively, the covariance of the estimate can be written as From the statistical properties of the estimated parameters, a confidence ellipsoid can be derived.With probability the estimated parameters x(t e ) lie within the con- fidence ellipsoid [13] E e(t i )e(t j ) = 2 i if i = j 0 otherwise . (1) (2) where 2 (n) is the -percentile of the 2 -distribution with n degrees of freedom.
Here n is the dimension of the parameter vector x.This set can also be seen as all parameter vectors, x , with squared Mahalanobis distance [14] less than 2 (n) from the normal distribution with mean x 0 and covariance P.

Problem Formulation
From Eqs. ( 2) and ( 4) we see that the inverse of the covariance depends on the sum of the contributions from the individual measurements.Hence, by choosing whether to perform the measurement or not, the shape of the covariance matrix, or the confidence ellipsoid, can be altered.Note that since each contribution in the sum H(t i ) T H(t i ) is positive (semi-) definite, every measurement will contribute and improve the estimate, that is, making the covariance smaller.The idea in this paper is to only make the measurements necessary to achieve the accuracy required by the intended use of the estimated parameters.As ( 2) and ( 4) do not depend on the actual measurements, they can be calculated beforehand and used to plan which measurements to make.
First, all possible observations from the different sensors are determined by propagating the position of the object forward in time and checking when the object is observable.This step requires an initial estimate of the object's orbit of sufficient accuracy.For each possible measurement, we calculate its contribution to the covariance matrix using (3) and ( 4).Next, we introduce the selection variables where N is the number of possible measurements.The selection variables determines if a measurement should be performed (1) or not (0).The inverse covariance matrix (4) can now be written as and the confidence ellipsoid as With each measurement we associate a cost.The cost describes the effort, in a broad sense, to perform the measurement.For example, certain sensors might be costlier to utilize both in time and money than others.Different observations using the same sensor may have different costs depending on when the measurement is made.
47 Page 6 of 20 final selection of the measurement set will be based on this cost, were we will aim to minimize the total cost while still achieving the required accuracy.
The planning problem can thus be stated as: find the cheapest combination of measurement such that the estimated parameters, using the measurement data, satisfy the accuracy requirement for their intended use.This can be expressed as an optimization problem: Note that the user is free to lump together multiple contributions to the inverse covariance matrix.For example, if multiple measurements are possible with one sensor during each passage, the user can choose to define a selection parameter either per measurement or per passage, or even split each passage into multiple parts, if desired.
In the next section we will show how to mathematically quantify when the estimated parameters are sufficiently accurate for the intended application, how to define the cost of a measurement, and how to formulate the planning problem as a convex optimization problem.

Application Cost and Application Set
To quantify if the estimated parameters are accurate enough for the intended application, we define an application cost function and an application set.The application function is a scalar function of the parameter vector x , and is denoted V app (x) .The function measures the degradation in performance when using the estimated parameters instead of the true parameters in the application.The application function should satisfy V app (x) ≥ 0 for all x , and V app (x 0 ) = 0 , i.e., the application cost is zero when the true underlying parameters are used.
The application set is then all the parameters that are accurate enough, i.e., all parameters that degrade the performance within an acceptable level.The set is defined as where is a positive scalar accuracy parameter and 1∕ is subsequently the maxi- mum acceptable degradation.The objective is thus to perform an orbit determination min s i total measurement cost s.t.P −1 is accurate for the application such that the estimated orbital parameters, x , lie within the application set for a given accuracy , i.e., x ∈ Θ app ().

Ellipsoidal Approximation of the Application Set
As discussed in Sect.2, the estimated parameters are Gaussian distributed.Due to the infinite support of the distribution, it is not possible to guarantee that the parameters lie within the application set.Instead we require that the application requirement is satisfied with a certain probability , i.e., The probability ( 6) is in general hard to evaluate and the constraint is non-convex in the decision variables s i .To this end, we perform two approximations to convexify the constraint.The first approximation is to require the estimation ellipsoidal ( 5), with confidence level , to lie within the application set instead of the estimate x itself, that is, The next step is to approximate the application set by an ellipsoid.The Taylor expansion of the application cost function at the true parameters x 0 is Since V app (x 0 ) = 0 and V � app (x 0 ) = 0 by definition, the application set can be approxi- mated by In some cases the Taylor ellipsoidal approximation may be an inaccurate representation of the application set, especially for low accuracy requirements.In this case a scenario approach can be used where a finite number of samples from within the application set is used to formulate the constraints, see Annergren [15].Other ellipsoidal approximations are also possible, e.g., by using Khachiyan's algorithm [16].
With the ellipsoidal approximation, the set requirement ( 7) can be approximated as The problem of estimating the orbital parameters accurately enough for the intended application can hence be achieved by shaping the confidence ellipsoid, by including the appropriate observations, to fit within the application ellipsoid.Using the expression for the confidence ellipsoid (5), the constraint (9) can be written as [10] (6) 47 Page 8 of 20 where A ⪰ B means that A − B is positive semi-definite.Equation ( 10) is a linear matrix inequality (LMI) and is a convex constraint on the decision variables s i [17].
To illustrate the use of the application set we present two examples.

Example 1
In the first example we want to guarantee that the estimated position of a satellite, at a certain epoch, is within 100 m of the true position of the satellite.The application cost function can hence be written as the squared distance from the estimated position to the true position where r = r x , r y , r z T is the position vector and r 0 = r 0,x , r 0,y , r 0,z T is the true position of the satellite.The application set is As the application set already represents an ellipsoid, the ellipsoidal approximation ( 8) is exact and can be written as (10) The application ellipsoid is shown in Fig. 1.Assume now that the error in the x-direction is less important and is only required to be known with an accuracy of 400 m.This could for example represent the along-track direction.The application set can in this case be written as The ellipsoid with the relaxation along the x-axis is also shown in Fig. 1.

Example 2
In the second example, we want to estimate the specific angular momentum of a satellite.The angular momentum is given by [12] where a is the semi-major axis, e the eccentricity, and is the gravitational constant.In this application we consider a satellite in a geostationary transfer orbit ( 350 × 35785 km).The objective is to estimate the semi-major axis and the eccen- tricity with such accuracy that the angular momentum is within 10% of the true angular momentum.An application cost function and the set reflecting this are selected as where a 0 and e 0 are the true values of the semi-major axis and the eccentricity, respectively.The Taylor approximation of the application ellipse can, after some calculations, be written as 47 Page 10 of 20 Figure 2 shows the level curves for the accuracy of the specific angular momentum as a function of the deviation in semi-major axis and eccentricity from their true values.The area where the accuracy is better than the required 10% is shown as well as the approximated ellipse.The application ellipse approximates the application set reasonably well.The semi-major axis of the application ellipse is infinite.This is expected since for each a there exists an e, such that the accuracy requirement is satisfied, and vice versa.

Cost Function
We want to select the cheapest combination of observations such that the estimated parameters satisfy the application requirements.First, we need to define the associated with each observation.The cost could be the monetary cost of performing the measurement but it can be used more generally to describe any effort associated with the observation, such as measurement time.It could also be used as a weighting to prioritize certain observations over others.Let c i be the cost of performing observation i.The total cost of the included observations can then be written as

Optimization Problem
We are now ready to formulate the complete optimization problem.Using ( 11) and ( 10) we get Although the objective function and the first constraint of ( 12) are both convex, the constraint that the decision variables are either zero or one is not and makes the optimization problem non-convex.In fact it is a combinatorial problem in the variables s i and the number of combinations to evaluate is in the order of 2 N , where N is the number of possible observations.Hence, the number of combinations is growing fast with the number of possible observations and quickly becomes computational prohibitive.To make the problem tractable, we propose to relax the binary constraint s i ∈ {0, 1} to 0 ≤ s i ≤ 1 .The new optimization problem then becomes and is now convex in s i .The optimization problem, with the LMI constraint is a semi-definite program (SDP) and can be solved efficiently, even for large problems, with modern optimization algorithms [17].At first, the relaxation might seem unintuitive as the solution may involve parts of observations.However, with the positivity constraint on s i , the cost function can be seen as the weighted 1 -norm of the decision variables, weighted with the observation costs.The 1 is a widely used heuristic to impose sparsity to the solution of a convex optimization problem, i.e., that many of the elements of the solution are zero, see [18] and [17].The solution to the optimization problem (13) will in many cases, as we will see in the simulation example in Sect.5, be sparse and many of ( 11) min 47 Page 12 of 20 the elements will be zero.However, all of the decision variables are not zero or one, although usually most of them are.This is expected since there are only a discrete number of values of yielding solutions containing only zeros or one.How to handle this is still an open question.Possible solutions include thresholding the decision variables or solving the combinatorial optimization problem (12) but only for the elements that are not zero or one.

Extension to Multiple Objects
To keep the notation simple, the derivation of the optimization problem in the previous sections was performed for a single object.However, the extension to finding the optimal observations times for multiple objects is straightforward.We here denote the new decision variable s j i where j = 1, … , N obj and N obj is the number of satellites or objects the observations should be planned for.As the number of possible observations could vary between the different objects, we now denote the number of possible times for observations for object j as N j .The cost function (11) can now be written as The cost coefficients c j i now also depend on the specific object and could potentially be used to reflect different priorities among the objects.
The application cost can be different for the different objects and hence the constraint (10) is similarly extended to The case when two or more observations are mutually exclusive, for example if two objects pass over a ground station at the same time and only one of them can be measured, can be handled by adding additional constraints to the optimization problem.If observation i 1 of object j 1 is mutually exclusive with observation i 2 of object j 2 , the constraint can be formulated as If there are no overlapping, mutually exclusive, observations, the optimization problem can be decomposed into a number of smaller optimization problems, one for each object.

The Unknown True Parameter Vector
The astute reader notes that the application set, and therefore the solution to the optimization problem, depends on the true orbital parameter vector x 0 .However, the objective is to select the observations to improve the knowledge of the parameters.If the true values are already known there is no need to do the observations in the first place.This contradiction is well known within the fields of system identification and experimental design but it is fundamentally true in any planning problem: without any initial knowledge, no planning is possible.The solution to this problem is often to use an initial estimate or a priori information of the parameter vector in the observation planning instead of the true parameter vector [15].The idea is that the shape of the application set is similar regardless of whether the initial estimate or the true parameter vector is used.In the orbit determination problems considered in this paper, the objective is to improve an estimate of the parameter vector from an earlier orbit determination.Hence, an initial estimate is often available.
Additionally, the planning can be updated iteratively.An initial estimate of the parameter vector is used to find an initial observation plan.As the measurements from the first actual observations are available, the data is used to estimate a new parameter vector.Future observations can then be re-planned by solving the optimization problem with the most recent estimate.This procedure can then be repeated when needed.

Simulation
The feasibility of the proposed observation-planning method is demonstrated in a simulation example.The example is idealised and simplified to readily illustrate the important concepts of the method.
The objective is to plan sensor measurements to estimate the orbital parameters of a satellite to sufficient accuracy.Two fictitious radar sensors, one in Kiruna in northern Sweden and the other in French Guiana, are available.The two sensors only measure the range to the satellite.Their location and properties are shown in  1.The simulated satellite orbit is given in Table 2.The orbital parameters are given in an inertial frame, aligned with the Earth-centered, Earth-fixed (ECEF) frame at the start of the simulation.
The objective is to estimate the relative state in a Hill-frame [11] of the satellite at 89900 seconds, or 15.5 orbits, after the start of the simulation.During the 15.5 orbits the satellite is visible from Kiruna at 12 passages and from French Guiana at 4 passages.For each passage, we assume that three measurement times are possible; the first at the highest elevation of the passage and the other two are 2 min before and after this moment.Hence, there are potentially 36 measurements from Kiruna and 12 from French Guiana to be used in the planning.The passages and the possible measurements are shown in Fig. 3, where the number indicates the temporal order of the measurements.
We use the Clohessy-Wiltshire equations [19] to model the dynamics relative to the ideal satellite.To be able to easily illustrate the ellipses we only estimate two parameters: the relative along-track and cross-track positions.The radial position and the velocity are assumed known.In real world these, and possibly additional parameters, would also have to be estimated.The application set is where r y and r z are the relative along-track and cross-track positions, respectively.The accuracy parameter will be varied to show the effect on the solution to the optimization problem.The optimization problem (13) is solved using CVXOPT [20], a Python package for convex optimization.We want to satisfy the application constraint with = 95% thus 2 (n) = 5.99 .The values of the selection parameters, s i , are shown in Fig. 4 as is varied.We see that only a few selection parameters at a time increase from zero to one, as the accuracy requirement increases.Furthermore, once a measurement has been selected ( s i has a value of one) it will in general be used also for higher values of .This shows that the convex relaxation of s i works as the resulting selection vector is sparse.
To gain some insight into why certain measurements are prioritized over others, we show the information content of each possible measurement in Fig. 5.As the sensors only measure the range to the satellite, each measurement is one dimensional.This gives that the information matrix I f (i) in (2) has rank 1.The direction, and the amount of information, that each measurement contributes is calculated as the non-zero eigenvalue times the corresponding eigenvector of I f (i) .For example, measurement 31 from Kiruna mostly contributes information in the along-track direction.The amount of information in a measurement depends on a number of factors.First it depends on the geometry during the observation.In Fig. 3 we see that due to the geometry of measurement 31 from Kiruna, it contains mostly information in the along-track direction.As the measurement noise standard deviation scales with the distance squared, the farther away the satellite is from the sensor the less information the measurement contains.Lastly, the amount of information depends on how much information the observable parameters at the measurement time contain about the orbital parameters at the estimation epoch.
In this example, when the application set is a circle, see Fig. 6, measurement 13 from Kiruna and measurement 8 from French Guiana are the first to be used.Even though the measurements by themselves are not the ones with most information, they are almost perpendicular to each other and hence are the optimal combination for finding the lowest-cost ellipse within the set.This can be seen also for higher values of where measurements with perpendicular information content are selected pairwise.For high values all measurements are selected, even though they have a low information content in the directions of interest.
Figure 6 shows the application and estimation ellipses from the solution to the optimization problem for = 100 .To verify the results we run 500 Monte Carlo simulations where actual measurement data is generated for the selected measurements.From the generated measurement data the orbital parameters are estimated by solving the non-linear least squares problem (1) using Powell's method [21].The resulting 500 estimates of the parameters are shown in Fig. 6.Out of the 500 Fig. 6 Application and the optimized estimation ellipses.The dots show the estimated parameters from the 500 Monte Carlo simulations, inside the estimation ellipse, inside the application ellipse, and outside the application ellipse.Here, the application requirements are equal along the cross-track and along-track directions estimates, 483 ≈ 97% are inside the estimation ellipse and 484 ≈ 97% inside the application ellipse, close to the desired 95%.
We now instead require higher accuracy in the along-track direction, i.e., The selection parameters, as a function of the accuracy , are shown in Fig. 7. Using Fig. 5, it can be clearly seen that measurements with a large amount of information in the along-track direction have higher priority.The resulting application and estimation ellipses as well as the results from 500 Monte Carlo simulations are shown in Fig. 8.

Conclusions and Future Work
In this paper, we presented a method to plan sensor measurement based on ideas from the field of system identification and input design.Previous methods have mainly focused on the problem of planning the measurements of different objects such that they are as accurate, in a wide sense, as possible, from a given set of possible observations.Here we view the problem differently; the objective is to find the least costly number of measurement such that the estimated orbital parameters are good enough for their intended application.We are hence not looking for the highest possible accuracy.To quantify if an estimated parameter vector is good enough we defined the application function and the application set.
We then showed how to formulate the planning problem as a convex optimization problem, which can be solved efficiently.Finally, the feasibility of the method was demonstrated in an idealised simulation example.As this is, to the best of the author's understanding, a novel approach to sensor planning, much work remains to be done.Several questions remain to be answered before it can be applied in a real-life context: • As discussed in Sect.4.6, the optimal planning strategy is based on knowledge of the true parameter vector.One solution to this is to base the planning on some initial estimate of the parameters and update the planning iteratively as more and more information is available.Another interesting approach could be to make the planning robust with respect to the uncertainties in the initial estimate, see [22] for examples from input design.These remedies and the sensitivity of the method to uncertainty in the initial estimate requires further investigation • The application functions and application sets shown in this paper are somewhat idealised and used only to illustrate the method.More realistic and practically useful application functions and sets are needed based on real-world applications.• The expression for the covariance matrix (2) is often too simple to represent the true uncertainty in the estimated parameters [5].Extending the method to use more accurate expressions for the covariance and the contribution of each measurement, should be explored further 4 0 0 1 • In this paper, we only applied the algorithm to relatively small problems with a limited number of objects and observation times.The performance of the method when applied to large-scale problems needs to be investigated • As shown in the simulation example, the decision variables are sometimes between 0 and 1.How to handle this is a topic for future research and the solution will probably depend on the application at hand.
Funding Open access funding provided by Swedish Defence Research Agency.This work was funded by the Swedish Armed Forces' Research and Technology Development programme.

Fig. 1
Fig. 1 Application ellipsoids for Example 1. Ellipsoid with equal requirement on all axis (left) and for a relaxed requirements along the x-axis (right)

Fig. 2
Fig. 2 Application ellipse for Example 2. The black lines represent level curves of the application cost function.The shaded red area is the true application set while the blue area is the ellipse approximation (Color figure online)

Fig. 3
Fig. 3 Satellite passages (blue lines) over the two radar sensors (red crosses) located in Kiruna (left) and French Guiana (right).The potential measurements are indicated by red dots and are enumerated by the corresponding number (Color figure online)

Fig. 4 Fig. 5
Fig.4 The value of the selection parameters s i as a function of the accuracy .The numbers in the figure link the selection parameters to the potential measurements in Fig.3.Note that only the first few numbers are shown to aid readability

Fig. 7 Fig. 8
Fig.7 The value of the selection parameters s i as a function of the accuracy .The numbers in the figure link the selection parameters to the potential measurements in Fig.3.Note that only the first few numbers are shown to aid readability

(t) = r(t) v(t) p q T , ̇x = f(t, x).
The function h i describes how the measurement at time t i depends on the state at the epoch t e , and e(t i ) are the measurement errors.

Table 1
Ground sensors and their properties used in the simulation example