domain: spi phys math info

text:
Clément Rappasse
email: clement.rappasse@gmail.com

Nicolas Merlinge
email: nicolas.merlinge@onera.fr

Baptiste Agez
email: baptiste.agez@onera.fr

Leonard Felicetti
email: leonard.felicetti@cranfield.ac.uk

Multi-Disciplinary Design Optimization for Relative Navigation in Non-cooperative Rendezvous

On-orbit applications, such as active debris removal, satellite refueling, maintenance and satellite health diagnosis require the ability for low-cost spacecraft to closely inspect other orbital objects in a non-cooperative manner. During this kind of mission, the relative navigation process becomes of critical importance for guaranteeing safe and collision-free proximity operations and maneuvers. The selection of relative navigation sensors appears being a critical task as it might drive the design choices in other subsystems (e.g., Attitude and orbit control system (AOCS), payload, power supply). This paper aims at studying the application of a Multidisciplinary Design Optimization (MDO) method to the design of AOCS and Navigation subsystems under small satellite constraints. The porposed MDO process is based on the optimization of navigation performance and mass reduction while respecting volume and power constraints for orbital close rendezvous. The navigation chain, including sensor simulation and navigation filter is simulated and integrated into the design cost function. The proposed MDO process based on a Genetic Algorithm (GA) and on an Extended Kalman Filter (EKF) aims at simplifying satellite design by determining a set of optimal admissible sensor combinations despite contradictory objectives on navigation and payload accuracy, mass reduction, power consumption and volume. Possible key advantages of the inclusion of the relative navigation subsystem within MDO process are the reduction of the design process time, the automation and optimization of the navigation architecture while respecting volume and power constraints of small satellites. A demonstration of the effectiveness of the proposed MDO method is provided on the benchmark of AVANTI mission.

INTRODUCTION

Non-cooperative rendezvous represents one of the key assets that can enable active debris removal operations and onorbit servicing missions. The close approach and docking to a non-cooperative target is indeed the first and necessary step of a sequence of possible operations that future mission can perform for de-orbiting dangerous pieces of debris, for carrying out health diagnoses, repair and refuel satellites in orbit [START_REF] Felicetti | Design of robotic manipulators for orbit removal of spent launchers' stages[END_REF] [START_REF] Shan | Review and comparison of active space debris capturing and removal methods[END_REF].

Compared to cooperative rendezvous, where both the chasing and target satellite cooperatively share information about their mutual positions, the relative navigation task during non-cooperative rendezvous appears to be more complex and critical, because it needs to be carried out only by the chaser satellite only. Only a few non-cooperative rendezvous were performed to date, mainly because of the risk of collisions and the technical performance of space-qualified relative navigation sensors [START_REF] Fehse | Rendezvous with and Capture / Removal of Non-Cooperative Bodies in Orbit: The Technical Challenges[END_REF] [START_REF] Opromolla | A review of cooperative and uncooperative spacecraft pose determination techniques for close-proximity operations[END_REF]. Such kind of issues become even more challenging when the non-cooperative rendezvous tasks need to be accomplished by small-sized satellites, where the selection of the relative navigation hardware becomes the key driver for the overall spacecraft and mission design. Indeed, the selection of a fairly accurate combination of sensors (e.g. cameras, LIDARs etc.) might contribute to the mission success or even enhance the performance of the mission, but in some cases the range of possible options can be limited by volume, power or even cost constraints. Such a selection might also have important consequences on the definition of the different architectures and components for other subsystems (e.g. Attitude and Orbit Control System (AOCS), payload, power subsystem etc.) and in some cases an optimal global design is not possible, due to concurrent needs and requirements imposed by the different subsystems and, consequently, only a compromise can be achieved.

The process of finding a good initial guess for the baseline architecture of the mission is difficult and somehow based on the experience of the system engineers who, in some cases, might not be aware of the latest state-of-art technology available on the market. It appears then evident that a tool addressing these initial challenges, by taking into account the multidisciplinary nature of the design process, might have tangible benefits in terms of reduction of the design process time, automation and optimization of the selection of specific key hardware, especially when the project is carried out in an agile design framework.

The utilization of Multidisciplinary Design Optimization 1 (MDO) tools during the mission and spacecraft design has gained an increased interest from both industrial and scientific applications in the past decades. Providing a trade-off support while accounting for decision-maker preferences has made MDO a critical tool when dealing with high-complexity design problems. It was first introduced and defined in late 1980s by J. Sobieski [START_REF] Sobieszczanski-Sobieski | Multidisciplinary Design Optimization: An Emerging New Engineering Discipline[END_REF] as "an emerging new engineering discipline" allowing the optimization of a system whose subsystems and physical phenomena are interrelated. In literature, many studies mention the use of MDO for designing aerospace vehicles. N. Bérend and S. Bertrand [START_REF] Bérend | MDO approach for early design of aerobraking orbital transfer vehicles[END_REF] developed a MDO methodology to optimize the mass of each discipline for aerobraking orbital transfer vehicles. W. Wu et al. [START_REF] Wu | Satellite multidisciplinary design optimization with a high-fidelity model[END_REF] proposed to include a parametric finite element model in the structure discipline of an MDO model to improve the fidelity of the analysis model and therefore optimize the design of the satellite. X. H. Wang et al. [START_REF] Wang | Comparison of MDO methods for an earth observation satellite[END_REF] compared three MDO methods (Collaborative Optimisation, All-at-Once and Multidisciplinary Feasible Method) to optimize the performance of an Earth observation satellite. J.T. Hwang et al. [START_REF] Hwang | Large-scale MDO of a small satellite using a novel framework for the solution of coupled systems and their derivatives[END_REF] developed a large-scale multidisciplinary design optimization involving seven disciplines and over 25,000 design variables for a small satellite. Recently, R. Shi et al. [START_REF] Shi | Surrogate assisted multidisciplinary design optimization for an allelectric GEO satellite[END_REF], whose objective was to reduce the total mass of an all-electric GEO satellite, proposed a surrogate assisted MDO framework involving most of the sub-systems. S. Berrezzoug et al. [START_REF] Berrezzoug | Interactive design and multidisciplinary optimization of geostationary communication satellite[END_REF] exposed the use of MDO for the design of a geostationary communication satellite. The application of the MDO to the design of small satellites performing non-cooperative rendezvous seems, to the knowledge of the authors, unexplored but worth further investigations.

The goal of this paper is to present a MDO process that can satisfy the specific objectives and multidisciplinary needs of the design of several satellite sub-systems for noncooperative rendezvous. Specifically, a particular attention will be given to the selection process addressing the relative navigation tasks, selecting the optimal combination of navigation sensors (IMU, LIDAR, Star tracker etc.) which feed a standard relative navigation algorithm (e.g. Extended Kalman Filter) with the aim of having a complete and accurate estimation of the relative position, attitude and velocity of the chaser satellite with respect to the target in a standard rendezvous scenario. Key performance index, such as the accuracy of the estimation (i.e. covariance matrix) can be indeed paired with classical system design parameters and included into a cost function in order to enable a more comprehensive and navigation-focused trade-off among the different options. A genetic algorithm is used to select and optimize the combination of sensors from a database of the state-of-art navigation technology alongside with a standardized definition of the key error models for each typology of sensors. Constraints on power, cost and volume, as well as inter-relations with other subsystems (such as AOCS, power etc) are taken into account during the process.

The remaining structure of the article is as follows. In Section 2, a formulation of a suitable MDO problem is proposed, for sizing and selecting components of AOCS, GNC and payloads for performing rendezvous. The section also defines the main coordinate frames and parameters used during the formulation as well as outlines the Kalman Filter that has been selected for the navigation task. Section 3 presents the MDO process that solves the problem by taking into account requirements and constraints imposed by other subsystems and their interactions with the relative navigation components. A database of the main characteristics of commercial-off-the-shelf components is then presented in Section 4, alongside with their error models. Section 5 presents the resulting optimal GNC configurations obtained by the application of proposed MDO algorithm to the test case scenario of the AVANTI mission. The performance analysis and numerical simulations with this section will show the consistency of the selected configuration as well as demonstrate the validity of the proposed MDO methodology. Finally, Section 6 summarizes the main findings, conclusions and future works of the authors in this area.

PROBLEM STATEMENT

The present work proposes a methodology to automatize and optimize the GNC design for a satellite that needs to perform a non-cooperative rendezvous within its mission. The considered model accounts for:

• AOCS actuators (e.g. thrusters, reaction wheels, etc. );

• AOCS and relative navigation sensors (e.g. IMUs, LI-DARs, Star-trackers );

• A relative navigation algorithm (e.g. Extended Kalman Filter) providing a complete navigation solution (relative position and velocity with respect to the target and absolute orientation of the chaser); • A payload whose performance is dependent from the AOCS and/or navigation chain. Without loss of generality, we consider a passive sensor pointing towards the target (e.g. a camera) and aiming to measure distances on the surface of the target.

The following subsections provide a baseline framework in order to define the considered rendezvous scenario (Section 2-1) and the adopted relative navigation algorithm (Section 2-2) that then concur to the formulation and specification the MDO problem for the design of a spacecraft performing noncooperative rendezvous (Section 2-3).

Non-cooperative rendezvous scenario

Figure 1 shows a typical non-cooperative rendezvous scenario where the chaser satellite (denoted with the letter B) is supposed to approach and rendezvous to the target satellite (denoted with the letter A). The chaser is therefore required to estimate its relative position through the use of its sensors (e.g. IMU, star sensor, LIDAR) and consequently perform a rendezvous maneuver by using the on-board actuators (e.g. thrusters and reaction wheels). Three different reference frames are defined: namely, i x, i ŷ, i ẑ shapes the Earth Central Inertial frame (ECI) frame, [ r x, r ŷ, r ẑ ] describes the Local Vertical Local Horizontal (LVLH) frame associated to the target and b x, b ŷ, b ẑ is the body frame associated to the chaser B. The positions of the two satellites are defined by r A and r B in the ECI frame, respectively and, δ δ δ r represents the relative position of the chaser with respect to the target. The latter is usually defined with the components taken with respect to the LVLH frame.

The variables characterizing the chaser satellite with respect to the target satellite are gathered in the state vector x described in Eq. ( 1): δ δ δ r and δ δ δ ṙ define the relative position and velocity in the LVLH frame of the target satellite, the quaternion q and the angular velocity ω ω ω describes the attitude of the chaser satellite in the body frame.

x = [δ δ δ r δ δ δ ṙ q ω ω ω] T

(1) To control the relative position, the relative velocity and the attitude of the chaser satellite, a control command vector u defined in Eq. ( 2) where u a provides the accelerations by the thrusters in the LVLH frame and u t , the torques in the body frame.

u = u a T u t T T (2)
The navigation problem consists of estimating the system's state x ∈ R d from a series of measurements x ∈ R d m The estimation scheme can be defined as an iterative stochastic process made of the state x and the measurements y. The state evolution is described by a dynamical model f :

x k = f (x k-1 , u k ) + w k (3)
where w k represents the process noise, and the measurements are linked to the state by an observation model :

y k = h(x k ) + v k (4)
where w k represents the observation noise.

Extended Kalman Filter

Relative and absolute navigation can be performed on-board by a Kalman Filter [START_REF] Ribeiro | Kalman and extended kalman filters : concept, derivation and properties[END_REF]. The Kalman Filter iteratively fuses the information provided by an a priori defined evolution model of the system with onboard measurements to obtain an estimation of the relative kinematic status of the chaser with respect to the target. The first step of the process consists in predicting the state vector (see sub-section 2-1) propagating the physical model. The second step updates the predicted state using the measurements. Figure 2 summarizes the steps and variables introduced in the Kalman Filter.

Time update-This step consists in predicting each variable of the state vector from time-step k using physical models to obtain a predicted state x k|k-1 . For example, in the case of close orbital rendezvous, the prediction of relative position and velocity is based on the Clohessy-Wiltshire equations [START_REF] Curtis | Orbital Mechanics for Engineering Students[END_REF]. The predicted covariance matrix P k|k-1 is obtained from the dynamical Jacobian matrix, the prior covariance matrix P k-1 and the process noise covariance Q k that reflects approximations made to solve the equations as linearization.

F is obtained by All variables of the state vector are not necessarily measured.

F = ∂ f ∂ x x=x k .
In this simulation, the relative position, the quaternion and the angular rate are measured by a LIDAR and a star sensor. The Kalman gain K k is also introduced in this step. It weights the importance of the measurements and of the model. R k is the measurement covariance matrix. In practice, R k is defined from the sensors performance characteristics and must be diagonal and positive definite.

To propagate the Kalman Filter, the initial state vector x 0 and the initial estimate covariance matrix P 0 are defined as initial conditions, according to the considered scenario.

Application to relative navigation-Generally, a navigation system relies on IMU iterative integration. However, the IMU navigation solution drifts because of gyrometric and accelerometric sensing errors. Therefore, IMU needs to be hybridized with additional sensors to correct the navigation solution. The whole navigation performance depends both on IMU accuracy and hybridization sensors accuracy. In the case of a relative navigation filtering scheme, IMU performance has an impact on covariance prediction (Time update) via the covariance matrix Q k . The hybridization sensors have an impact on the correction step via the covariance matrix R k used to compute the Kalman gain K k (Measurement update).

Multidisciplinary Design Optimization

Usually, the search for an optimal solution in engineering problems is performed thanks to gradient based methods, especially when the number of variables or functions is easy to deal with. Most complex problems, however, such as satellite design, requires suitable tools, capable of operating in the presence of a considerable number of variables and couplings [START_REF] Ravanbakhsh | Multidisciplinary design Optimization application to conceptual design of university-class microsatellite projects[END_REF].

Those tools, mostly used in the field of design optimization, are called stochastic methods. As opposed to deterministic methods, stochastic methods randomly search in the design space, free from specific rule or initial guesses. In the MDO process, this stochastic aspect lies upon the choice of the algorithm.

Considering their ease of use, Evolutionary Algorithms (EA) such as Genetic Algorithms (GA) are most of the time chosen [START_REF] Jafarsalehi | Collaborative optimization of remote sensing small satellite mission using genetic algorithms[END_REF]. Inspired by biological evolution, EA mathematically reproduces the natural mechanisms of mutation, recombination and selection over a set of individuals forming a population [START_REF] Hassan | A comparison of particle swarm optimization and the genetic algorithm[END_REF]. Selection operator gives the fittest solutions the highest probability to be in the next generation; recombination (or crossover) creates a new solution by mixing two solutions from the previous generation; and eventually mutation can randomly modify some solutions. The flow chart of such algorithm [START_REF] Leal-Naranjo | Multiobjective optimization of a parallel manipulator for the design of a prosthetic arm using genetic algorithms[END_REF] is given in Figure 3. In MDO terms, GA can deal with large design space, being able to handle optimization problems with many local optima. This capacity, along with the ease of use and adequation to discrete design variables [START_REF] Jafarsalehi | Collaborative optimization of remote sensing small satellite mission using genetic algorithms[END_REF], makes GA the chosen algorithm for the present optimization problem.

Generation of initial population

Inside the MDO process, the algorithm deals with the design variables, which consitute the population of individuals. At each iteration, the fitness of every candidate is evaluated thanks to mathematical models. The fitness is related to an objective function, in which the user puts different criteria of his choice, and to constraints, if any, that the solution must respect.

A MDO problem for a single objective can thus be mathematically written as follow [START_REF] Balesdent | A survey of multidisciplinary design optimization methods in launch vehicle design[END_REF]:

X * = argmin G(X) ≤ 0 X ∈ Ω ≤ X F(X) (5)
where F(X) is the objective function, G(X) is the constraint function and X the vector of design variables which can move freely inside the design space Ω.

In practice, the process can be divided into several sequential steps :

• Definition and selection of the design variables, which have the major impact on design;

• Evaluation of design for a set a variables thanks to the mathematical model of the system;

• Computation of the objective function for this specific design;

• Iterations on variables via an optimization algorithm to get the best design configuration.

Several approaches exist and use different architectures to tackle the problem, including MDF (Multi-Disciplinary Feasible), CO (Collaborative Optimization), AAO (All-At-Once), BLISS (Bi-Level Integrated Synthesis), CSSO (Concurrent Subspace Optimization) among others. Their major difference lies in their way of proceeding and organizing the different parts (variables, objectives, constraints) of the process. For instance, they can be single or multi-level, meaning that they would perform subsystem's optimization before proceeding to system-level optimization. Interested reader can refer to [START_REF] Martins | Multidisciplinary design optimization: A survey of architectures[END_REF] for more information.

In this study, two main objectives are included in the objective function: minimizing the mass and maximizing the performances (i.e. minimizing their errors) of the satellite embedded sensors. Both are the results of Multidisciplinary Analysis (MDA) of several satellite's subsystems, described in section 3-2. The constraints associated to the problem are the volume occupied by the embedded devices (sensors and actuators) and their power consumption. A single level Allat-Once (AAO) framework is used to formulate our problem :

X * = argmin V (X) -V ≤ 0 P(X) -P ≤ 0 X ∈ Ω F [M(X), σ p (X), σ v (X), σ q (X), ε ε ε(X)] (6)
where X ∈ N N is a vector of discrete variables pointing towards a combination of N embedded devices (e.g. navigation sensors, thrusters, payload sensors), V (X) is the total volume occupied by the selected devices, V is the maximum admissible volume, P(X) is the total power consumption for the selected devices, P is the maximum admissible power, Ω ⊂ N N is the design space (i.e. the list of all possible device combinations), M(X) is the total mass of the selected devices, σ p (X), σ v (X), and σ q (X) are respectively performance indicators for navigation (relative position, relative velocity, and absolute attitude) and ε ε ε(X) contains some performance indicator about the payload in terms of error (e.g. position and velocity error estimation of a point located on the surface of the target).

In the specific case of a device combination involving: one IMU, one star tracker (ST), one LIDAR (L), a payload (P), one type of thrusters (TH) and one type of reaction wheels (RW), the design vector is defined as:

X = [XIMU X ST X L X P X T H X RW ] T (7)

METHODOLOGY

This section presents the proposed MDO methodology to solve the Problem (Eq. ( 6)). An overview of the MDO process is first given. Then, the modeling of each subsystem is detailed in terms of cost function and constraint and represented in Figure 4.

MDO methodology

MDO relies on an optimization algorithm that determines the optimal device configuration (i.e., navigation sensors, payload and actuators) given design constraints (e.g., volume, power). Since the considered design variables are discrete choices of COTS sensors, the optimization algorithm must handle discrete optimization space. As mentioned previously, GA is selected: this type of optimization algorithm iterates on a population of candidate configurations and aims to achieve a consensus. A convergence is met when a consensus is achieved among the population of candidate solutions. A given candidate solution is evaluated by cost functions that are developed in the following sections.

For a given device configuration, the sensor measurements y k are simulated from a reference trajectory and processed by the navigation filter. Scalar performance criteria are defined to quantify the navigation performance at a given time-step k e of the scenario. Actuators and payload performance requirements are also modeled, depending on sensors and navigation performance. Finally, structural and design characteristics are evaluated (e.g. sum of all devices mass, volume and power consumption). This scheme is characterized by strongly nonlinear and interdependent couplings between involved subsystems, which motivates the use of an optimization algorithm.

Figure 4: Overview of the proposed MDO methodology: a discrete optimization algorithm evaluates the cost and constraints of Problem (Eq. ( 6)) by running a temporal loop simulating the whole navigation process. The MDO process aims to determine the best sensor configuration X * .

Subsystems modeling

This section details the modeling of all considered subsystems of Problem [START_REF] Bérend | MDO approach for early design of aerobraking orbital transfer vehicles[END_REF]. The following assumptions are made:

• The nominal trajectory of the chaser is a constraint given by a prior mission preparation process and is not a design variable;

• The flight avionic calculator is not considered as a design variable;

• Each device is dedicated to one specific task;

• The total available volume for all the devices is a constraint;

• The power supply system is fixed and constrains the total energy available on-board;

• The relative positioning of the devices in the body frame is not considered.

AOCS and GNC performance-The attitude and orbit control subsystem is directly related to the guidance and navigation subsystem. The choice of sensors is the key part of the MDO process, and is thus subject to higher fidelity than the other models.

A database of COTS sensors (IMU, LIDAR, Star Tracker), described in Section 4, has been built and provides figures for the sensors and actuators: expected performance, mass, volume, and power consumption. The relative navigation algorithm (e.g., the Extended Kalman Filter presented in Section 2-2) provides a complete navigation solution (relative position and velocity with respect to the target and absolute orientation of the chaser). The navigation performances indicators are defined by the spectral radius ρ (•) of the filter's covariance sub-matrix for each sub-state vector (e.g. position, velocity, attitude) at an a priori defined time of evaluation k e (e.g. at final time):

σ p (X) ρ( P p k e ) σ v (X) ρ( P v k e ) σ q (X) ρ( P q k e ) (8)
where P

p k e , P v k e , and P q k e are the position, velocity and attitude sub-matrices of the filter's covariance P k e . Note that this definition is not dependent of the specific navigation algorithm and can be applied to other recursive filters.

To complete the design of the AOCS subsystem, actuators must also be considered. Reaction wheels, combined with small thrusters, are selected. This set up is often used for 3axis control spacecraft [START_REF] Wertz | Space Mission Engineering: The New SMAD[END_REF]. Thrusters are necessary provide small orbit corrections, while wheels are used for orientation, attitude perturbations and relative maneuvers.

As for sensors, a database of reaction wheels and thrusters for small satellites applications was built, that outputs mass, power consumption, and performance for each device, namely torque capacity for the wheels and thrust and Isp for thrusters. In the frame of the MDO process, actuators come into play when sensors have been selected at each evaluation i. Once the sensors' performance and their respective accuracies are known, it is then possible to compute the requirement in terms of ∆v and torque.

In this paper, the function of thrusters is dual. They desaturate the wheels and provide motion for orbit corrections. These corrections, considered for ∆v calculations, are based on the relative position and velocity uncertainties obtained from the sensors and the navigation filter performance. Specific Impulse requirements are then addressed thanks to Tsiolkovsky's equation :

Isp = ∆v ln( m init m f inal )g 0 ( 9
)
where g 0 = 9.81 m/s 2 is the gravitational constant, m init the initial mass (before thrust) and m f inal the final mass. ∆v must account for both velocity and position errors: considering a distance correction of σ p to achieve within a period t man , two maneuvers (one for catching back the missing distance, one for slowing down back to previous velocity) are required with total absolute value 2σ p , which yields ∆v = σ v + σ p t man .

On the other hand, momentum dumping requires a thrust force F that can be determined by the following formula [START_REF] Wertz | Space Mission Engineering: The New SMAD[END_REF]: [START_REF] Shi | Surrogate assisted multidisciplinary design optimization for an allelectric GEO satellite[END_REF] with h w the stored momentum, t burn the burn time and L the moment arm.

F = h w Lt burn
The same process can be applied for the wheels. The values of the uncertainties on pitch/yaw and roll orientations resulting from sensors' performance are taken as inputs to determine the torque capacity of the wheel necessary to account for those perturbations. The necessary slewing torque is [START_REF] Wertz | Space Mission Engineering: The New SMAD[END_REF]:

T = 4θ I t 2 spin ( 11
)
with I the inertia of the satellite, t spin the spinning period, and θ the slew angle required. In our case, θ constitutes the maximum uncertainty angle among the three Euler angles obtained from the star tracker.

The second characteristic of wheels is their angular momentum capacity, which allows a wheel to store the acquired momentum. This parameter is obtained with the torque T and the acquisition time t acq [START_REF] Wertz | Space Mission Engineering: The New SMAD[END_REF]:

h w = 2T t acq √ 2 (12)
Payload-In addition to the AOCS and GNC sensor chain, a payload device is considered. The payload performance may be coupled with AOCS and/or navigation sensors. In this work, we consider a relative passive sensor pointing towards the target (e.g. a camera) and aiming to determine:

• Distances on the surface of the target;

• Velocity of points on the surface of the target (e.g. if the target is spinning).

The payload is assumed to measure an angle α ∈ [-π, π] between the Lines Of Sights (LOS) associated to two points X 1 and X 2 of the target expressed in the sensor frame, as shown in Figure 4. This corresponds to a distance ∆X ≥ 0 on the surface of the target. The relative range between the chaser and the target is assumed to be significantly larger than ∆X. Therefore, the measurement can be approached by α = arctan ∆X r where r is the relative range. The involved uncertainties are the payload sensor resolution ε α , the evaluation error ε ∆X on the ∆X distance, and the navigation range error ε r . Such uncertainties affect the measurements of the angle α as follows:

α = α ± ε α = ∆X ± ε ∆X r ± ε r (13)
which yields:

ε ∆X ≤ r|ε α | + |ε r | ∆X r + |ε r ||ε α | ( 14
)
where α is assumed to be small. Assuming that the size of the target is known, the range error r can be related to the relative navigation uncertainty in position. In practice, this can be achieved by considering ε r as an error bound at n σ p on the position components (with, e.g., n = 3 and σ p the spectral radius of the filter covariance P k on position components). ε ∆X can be taken as a payload performance evaluation.

A similar model can be designed to quantify the velocity evaluation error for a point X of the target. Consider the traveled distance ∆X of X during a time ∆t > 0. The estimated velocity of X in LVLH frame can be expressed as:

v x ≤ ∆X + |ε ∆x | ∆t + (v s x + |ε v |) + (r + |ε r |)(Ω s z + |ε Ω |) (15)
where v s x is the chaser relative velocity projected to the component x of sensor frame and Ω s z is the chaser angular rate projected to the component z of sensor frame. ε V and ε Ω can be respectively be linked to relative velocity navigation uncertainty and angular rate noise (gyrometers). An upper bound of the velocity payload performance evaluation criterion ε V x can then be derived from Eq. ( 15), for example by taking the worst case for involved uncertainties. The payload cost in Problem ( 6) can be defined by:

ε ε ε = ε ∆X ε v x (16)
Power budget-Usual constraints on small satellites include the power budget dedicated to each subsystem. Power subsystem is considered to be already designed, providing as input the amount of power available for combined AOCS and payload subsystems. The power budget equation is defined in Eq. [START_REF] Leal-Naranjo | Multiobjective optimization of a parallel manipulator for the design of a prosthetic arm using genetic algorithms[END_REF]. Power consumption of devices is taken from the different databases.

P (X) = N ∑ i=1 P X[i] (17)
where P X[i] is the power consumption of device X[i].

Mass and Volume budgets-Mass is a critical characteristic of satellite conception. Directly related to the cost, it most of the times drives (or at least comes into play for) the key choices of the design process. As for power, the masses of the different devices are extracted from their respective databases. The mass and the volume budget equation are also based on a total sum.

M (X) = N ∑ i=1 M X[i] (18)
V (X) = N ∑ i=1 V X[i] (19)
where M X[i] and V X[i] are the mass and the volume of device X[i].

DATA DESIGN OPTIONS

The MDO process determines the most adequate combination of sensors and actuators from databases to perform onorbit rendezvous in LEO. In this section, performances and models of COTS sensors and actuators are presented. In what follows, a unique number is assigned to each sensor to anonymize the brand to which it belongs. Obtained values are taken from datasheets provided by manufacturers.

Navigation sensors database

Inertial Measurements Unit (IMU)-Inertial measurements units are widely used devices in the aerospace field to estimate position, velocity and attitude by numerically integrating specific acceleration and angular rate. IMUs are made up of three accelerometers and three gyroscopes. Due to the iterative integration of measurements, the resulting navigation solution is affected by a drifting error. In the present work, twenty-six COTS IMUs are used.

Gyroscope and accelerometer were modeled by Eq. ( 20).

The acceleration γ γ γ and the angular rate ω ω ω are expressed in the body frame. γ γ γ true and ω ω ω true correspond to the actual accelerations and angular rates. ξ ξ ξ , ν ν ν and τ represent the bias, the noise and the scale factor, respectively.

γ γ γ meas = γ γ γ true (1 + τ acc ) + ξ ξ ξ acc + ν ν ν acc ω ω ω meas = ω ω ω true (1 + τ gyro ) + ξ ξ ξ gyro + ν ν ν gyro ( 20
)
The following figures illustrate the distribution of IMU performance features within the whole database. Figure 6 and7 describe the gyroscope bias and the accelerometers bias with respect to the mass of the IMU itself, respectively. These figures highlight the complexity of selecting a sensor from the full device set in an optimal way. Ideally, the selected sensor should have a low bias and a low mass. However, better accuracies are generally obtained with heavier sensors, while low mass sensors can even have two order of magnitude worse performance (e.g. bias). A similar behavior appears with IMU random walks and scale factors. This motivates the use of a global optimization algorithm to find an appropriate trade-off. Star sensors-Star sensors aim to determine the attitude of the satellite by matching between optical measurements and an embedded celestial map. In recent years, the use of this type of sensor has increased continuously thanks to its high accuracy range (from 1 to 10 arcsec) [START_REF] Wertz | Space Mission Engineering: The New SMAD[END_REF]. However, the high accuracy of these sensors yields high cost, weight and size, which requires a challenging design trade-off. Forty-three COTS star sensors are considered in the database.

Eq. ( 21) models the observation equation of the star sensor. The true quaternion q true represents the actual attitude quaternion of the chaser. The star sensor bias ξ ST and noise ν ST are provided by manufacturers. q ξ and q ν are the quaternions representing the pointing errors due to the bias and random noise of the instrument, respectively.

q meas = q ξ q ν q true (21) respectively correspond to the bias and the noise of the LIDAR on the angular measurements (azimuth and elevation) and on the relative distance (range). Eq. ( 23) describes the conversion from the relative position coordinates to the spherical coordinates.

         θ meas l = arctan δ z √ δ x 2 +δ y 2 ϕ meas l = arctan δ y δ x δ r meas l = δ x 2 + δ y 2 + δ z 2 (23)
Based on open literature, three LIDARs are considered in the database. As the accuracy of LIDARs depends on the distance between the chaser and the target, it was chosen to compare the bias with the mass for a reference range of 100 meters. Nineteen COTS cameras were listed (see e.g. Figure 11, accuracy as function of the sensor mass).

Actuators database

Reaction wheels-Reaction wheels are torque motors with integrated rotors that can spin in either direction, providing one axis control per wheel [START_REF] Wertz | Space Mission Engineering: The New SMAD[END_REF]. They are used to slew the spacecraft, i.e. reorienting the vehicle when needed. Wheels most often deal with disturbances and pointing perturbations, to which they respond by spinning at higher speed, creating a torque. In the present work, twenty-six COTS reaction wheels are listed with their mass, power, and performances. Figures 12 and13 show respectively the torque and the angular momentum of the listed reaction wheels with respect to their mass. Thrusters-In the scope of AOCS subsystem, thrusters are mostly used for orbit maintenance, desaturation of the wheels if required, and for low-∆v maneuvers. The working principle lies in the expulsion of mass which produces a force or a torque. In the database, eighteen thrusters are listed with their mass, power, Isp and thrust. Two types of thruster were selected, cold gas and monopropellant. Cold gas systems consist of pressurized propellant tank, valves, feed pipes and a nozzle, and provide thrust thanks to the stored enthalpy of the gas. On the opposite, liquid chemical propulsion systems include, as its main parts, a combustion chamber and a convergent-divergent nozzle. These two types of chemical propulsion have the advantage of having excellent flight heritage and are reliable and simple [START_REF] Mani | Combined Chemical-Electric Propulsion Design and Hybrid Trajectories for Stand-Alone Deep-Space CubeSats[END_REF]. Figure 14 and15 show the Isp and the thrust of the listed thrusters with respect to their mass.

TEST CASE

This section presents an application of the proposed methodology to a specific test case scenario. The AVANTI mission [START_REF] Ardaens | Flight demonstration of spaceborne real-time angles-only navigation to a noncooperative target in low earth orbit[END_REF] is first introduced and a test case trajectory is defined, based on the open literature. MDO results are then provided to illustrate the behavior of the design process on this testcase.

AVANTI mission: navigation for non collaborative rendezvous

To assess the method and evaluate the performance of the MDO method, the scenario of the AVANTI mission was selected. In 2016, the German Aerospace Center (DLR) conducted the AVANTI mission that aimed at performing two scenarii between non-cooperative satellites. The first experiment allowed an approach from 13 km to 3 km. The second, an approach from 3 km to 50 m. Both of them were successful [START_REF] Ardaens | Flight demonstration of spaceborne real-time angles-only navigation to a noncooperative target in low earth orbit[END_REF].

Physics modeling-Two sets of equations are used to model the absolute and relative dynamics of the chaser satellite. Clohessy-Wiltshire equations [START_REF] Curtis | Orbital Mechanics for Engineering Students[END_REF] are used to model the relative translational motion of the chaser with respect to the LVLH reference frame, defined in Section 2-1.

   δ ẍ -3n 2 δ x -2nδ ẏ = u ax δ ÿ + 2nδ ẋ = u ay δ z -n 2 δ z = u az ( 24
)
where n is the mean motion of the chaser satellite and u ax , u ay and u az are the components of the acceleration provided by the thrusters, defined in the LVLH coordinate frame. Angular velocity and quaternion derivatives (respectively ω ω ω and q) are modeled by the attitude dynamic and kinematic equations (Eq. ( 25)). T ext = T dist + u t is the sum of the external torques that are applied on the satellite. This corresponds to the sum of the environmental torques T dist and of the commanded torques u t defined with respect to the body reference frame of the chaser spacecraft. I is the inertia matrix of the satellite in the body frame. ω ω ω is the angular rate vector of the satellite in the inertial frame and q the current quaternion.

ω ω ω = I -1 [T ext -ω ω ω × (Iω ω ω)] q = 1 2 Ω Ω Ωq (25)
with:

Ω Ω Ω =    0 ω z -ω y ω x -ω z 0 ω x ω y ω y -ω x 0 ω z -ω x -ω y -ω z 0    (26)
The two main disturbances acting on the system are the gravity gradient torque T grav and the aerodynamic drag torque T drag [START_REF] Sunde | Sensor modelling and attitude determination for micro-satellite[END_REF]. Their models is described in Eq. ( 27). r B corresponds to its position in the ECI frame. C d , A and V are respectively the drag coefficient, the cross sectional area and the velocity of the chaser satellite, respectively. C p and C m are the position of the center of pressure and the position of the center of mass in the body frame. ρ is the atmospheric density.

     T dist = T grav + T drag T grav = 3n 2 r B ×I r B r 2 B T drag = ρC d AV 2 2 × (C p -C m ) (27)
The complete set of equations of motion includes Eq. ( 24) and Eq. ( 25) and can be rewritten by adopting a state space representation, where x is defined as in Eq. ( 1) and the control actions u as in Eq. ( 2). The resulting system can be then rewritten as in Eq. ( 28). In the following equations, the inertia matrix I is considered as a diagonal matrix in the body frame.

ẋ = F x + G u + G T (28)
where:

F = F 1 0 6×7 0 7×6 F 2 G = G 1 0 6×7 0 7×6 G 2 T = 0 1×6 T T dist T (29)
F 1 =       0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 3n 2 0 0 0 2n 0 0 0 0 -2n 0 0 0 0 -n 2 1 0 0       G 1 = 0 3×3 0 3×3 0 3×3 I 3×3 (30)
F 2 = F kin 0 4×3 0 3×4 F dyn G 2 = 0 3×3 0 3×3 0 3×3 I -1 (31)
F kin = 1 2 Ω Ω Ω F dyn = 1 2     0 ω z (Iyy-Izz)
I zz 0     (32)
Initial conditions and modeled trajectory-The close range rendezvous experiment of the AVANTI mission was chosen as a baseline test case for the simulations [START_REF] Gaias | Flight demonstration of autonomous noncooperative rendezvous in low earth orbit[END_REF]. Specifically, the chaser named BIROS is initially placed in a sunsynchronous orbit. The orbital elements of the orbit and relevant physical characteristics of BIROS are described in Table 1 and 2 respectively. The target satellite, BEESAT-4, is launched from BIROS using a picolauncher [START_REF] Wermuth | Operational Concept of a Picosatellite Release From a Leo[END_REF]. To avoid a collision situation, a spiral-like approach was performed. In this work, we derived a scenario representing a portion of this phase. Initial position and velocity of the chaser, as well as the accelerations, provided during a maneuver were derived from the analysis of [START_REF] Gaias | Flight demonstration of autonomous noncooperative rendezvous in low earth orbit[END_REF] and are reported in Figure 16, which shows the results of the Clohessy-Wiltshire equations propagation.

Item
q 0 =     q 0 x q 0 y q 0 z q 0 scalar     =        1 4q 0 scalar R 0 2,3 -R 0 3,2 1 4q 0 scalar R 0 3,1 -R 0 1,3 1 4q 0 scalar R 0 1,2 -R 0 2,1 1 2 1 + R 0 1,1 + R 0 2,2 + R 0 3,3        (33)
This DCM is calculated thanks to the Eq. ( 34) where b R 0 cw corresponds to the initial rotation from the body frame to the LVLH frame. It was arbitrarily considered that the x, y and z axes of body frame are respectively aligned with the z, x and y axes of the LVLH frame, respectively. cw R 0 i is computed thanks to the initial position and velocity of the chaser in the inertial frame.

R 0 = b R 0 i = b R 0 cw cw R 0 i ( 34
)
The initial angular rate was set to [0.05 0 0] • .s -1 .

Numerical results

The MDO methodology presented in Section 3-1 is tested on the AVANTI scenario defined in Section 5-1 and the database introduced in Section 4. The design space corresponds to a total of 29,823,768 device configurations (i.e., 18

× 26 × 19 × 3 × 43 × 26).
Optimization's tools definition-To perform the optimization, a Python library called Platypus is used [START_REF] Hadka | Platypus-multiobjective optimization in python[END_REF]. This library provides optimization algorithms and analysis tools for multiobjective optimization. It includes NSGAII [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: Nsgaii[END_REF], a Nondominated Sorting Genetic Algorithm II, which is selected to run the simulation.

The visualization of the results is made throughout the utilization of the Pareto front. To explain the reading of this criterion, two usual terms are defined: feasible solution and dominated solution.

A feasible solution is a vector included in the design space, which respects the constraints. A feasible solution a dominates a feasible solution b if and only if [START_REF] Balesdent | Multidisciplinary Design Optimization, application to aerospace vehicle design[END_REF]:

∀i ∈ [1, .., m], F i (a) ≤ F i (b) and ∃ j ∈ [1, .., m], F j (a) < F j (b) (35)
with F the objective function and m the length of a solution vector.

The Pareto front (also called optimal Pareto set) is then the set of the feasible solutions that are non dominated, and can be represented, in the case of a two-objective optimization, with a single two-dimensional plot.

The convergence of the simulation is obtained via the hypervolume indicator. This tool allows to perform quantitative assessment of the Pareto front, by measuring the size of the portion of the objective space that is dominated by the solutions [START_REF] While | A faster algorithm for calculating hypervolume[END_REF]. This measure is made relative to a reference point, ideally far from the Pareto front, from which this portion can be computed. Represented on Figure 17, the hypervolume I H (A) of the solution set A can mathematically be defined as the Lebesgue measure λ of the set H(A, R) : of the plots should be like this : a value of 2 on the Y-axis for an individual on, for instance, the torque plot (bottom middle graph on Figure 21) means that this particluar devices' configuration has in capacity twice the necessary torque. More than just a unique optimum, the Pareto front outputs a set of optimal solutions. These solutions are best compromises between the two antagonists objectives of mass reduction and performance.

I H (A) = λ (H(A, R)) = R 2 1 H(A,R) (z)dz (36)
In the present work, equal weight has been given to both inside the objective function: however, it is possible to modifiy this weight, and value more one over another.

But even in such case, a whole set of optimal solutions is presented to the decision-maker. This leaves the possibility of choice among the solutions, and especially the ability to visually find which one suits the best the needs.

To this extent, two solutions are proposed in this work : one minimizing mass and the other focusing on performance (errors' minimization). They are shown in Figure 18, as mass driven solution and performance driven solution, respectively. Tables 3 and4 summarizes these two possible combinations of components obtained with MDO.

Within the design space (i.e. all available COTS components), the cost to gain mass (or performance) while still respecting the constraints expands fast when reaching the limits of the design area (Figure 18). Moreover, and as it could be expected, some components have greater influence on the results than others. Choices of actuators stay limited to fit in the different constraints. On the other hand, sensors and payload offer broader ranges of possibilities. It is worth noting that some components can sometimes have barely the same performances or characteristics, which can cause final values to be shared among those equivalent components. Performance indicators are navigation filter accuracy in relative position σ p , relative velocity σ v and attitude σ q (in terms of navigation filter standard deviation), the payload performance in position detection ε ∆X and velocity evaluation ε v x .

A significant impact of the choice of the sensors' configuration can be noticed on all criteria, especially the navigation velocity, which is a critical point for collision avoidance. Since the actual value of relative velocity is about 0.1 m/s, it needs to be estimated with an accuracy of one or two orders of magnitude below.

Figure 22 illustrates the filter's relative position estimation error and covariance during the whole trajectory for two sensor configurations.

Finally, one of the key advantage of MDO is the gain of computational time. Whereas it would have taken the computation of all the 29,823,768 possible devices configurations, only 320,000 (80 individuals for 4000 evaluations) were necessary to reach acceptable convergence and results for the MDO process. This represents a decrease of computational time of 2 orders of magnitude (98%), which is a significant time reduction.

Low quality configuration (position error):

MDO solution output 2 (position error):

CONCLUSION

In this paper, a MDO methodology including navigation sensors choices has been introduced for satellite design in the context of non cooperative rendezvous. The optimization problem is formulated in terms of mass reduction and navigation and payload performance maximization, subject to volume and power consumption constraints. The optimization process is based on a Genetic Algorithm that aims to find the optimal sensors and actuators configuration among ). The MDO solutions 1 and 2 respectively correspond to the configurations described in Tables 3 and4 The navigation performance evaluation involves the simulation of the whole navigation chain, including sensor emulation and navigation filter evaluation (e.g. an Extended Kalman Filter). Accounting for the whole navigation chain in the design process is critical in the context of non cooperative rendezvous, since sensor uncertainties are coupled in a nonlinear way. Furthermore, the navigation performance is highly related to the mission success and to the spacecraft safety, as well as the other sub-systems.

The proposed MDO scheme offers the interest of providing not only an optimal solution among a large combinatory of devices, but also a subset of admissible solutions using the Pareto front criterion. This subset of design solutions can be used to orient the work of design engineers as a primary hardware solution selection. The identified design solutions can then be further studied or included in more extensive MDO schemes.

The methodology was applied to the AVANTI mission scenario. Obtained results demonstrate that the proposed MDO process allows for an optimization of the relative navigation equipment while satisfying the spacecraft design constraints.

Future works may involve more extensive design variables, such as other subsystems, trajectory and maneuver optimization, or financial cost minimization. There could be some interest in adapting the algorithm if more subsystems were to be added, such as switching from NSGAII (single-level) to Collaborative Optimization (multi-level) to reduce simulation time.

The proposed approach can be considered as a linchpin to automatize design processes that include recursive estimation algorithms and can be extended to other fields of study such as aeronautics, automotive, medical or mobile robotics.

Figure 1 :

 1 Figure 1: Position of chaser B relative to target A

Figure 2 :

 2 Figure 2: Extended Kalman Filter flowchart

Figure 3 :

 3 Figure 3: General flow chart of a genetic algorithm.

Figure 5 :

 5 Figure 5: Payload relative sensor scheme. The range r is assumed to be significantly larger than the distance ∆X between two points of the target.

Figure 6 :

 6 Figure 6: IMU gyroscope: Bias vs. mass

Figure 7 :

 7 Figure 7: IMU accelerometer: Bias vs. mass

Figure 8

 8 Figure 8 and 9 describe the pitch, yaw and roll bias as a function of the mass of the star sensor. As for the IMUs, the relation between accuracy and mass reduction describes the complexity of realizing trade-offs. The large number of star sensors increases the complexity of the selection.

Figure 8 :Figure 9 :

 89 Figure 8: Star sensor: Pitch and yaw bias vs. mass

Figure 10 :

 10 Figure 10: LIDAR: Bias vs. mass at 100 m

Figure 11 :

 11 Figure 11: Payload: Accuracy vs. mass

Figure 12 :Figure 13 :

 1213 Figure 12: Reaction wheel: Torque vs. mass

Figure 14 :Figure 15 :

 1415 Figure 14: Thruster: Isp vs. mass

Figure 16 :

 16 Figure 16: Modelled trajectory

Figure 20 :

 20 Figure 20: Convergence of simulation represented by Hypervolume.

Figure 21 :

 21 Figure 21: Configurations final values for constraints margins. Along with Power and Volume, we also show the margins left for the actuators capacities.

3 .

 3 Performance of the MDO algorithm Table 5 presents the navigation performance for two specific sensor configurations. The best configuration corresponds to the output of an unconstrained optimization (IMU 36, Star Tracker 87, LIDAR 186, Payload 171). The low quality configuration illustrates the impact of low accuracy sensors on navigation performance (IMU 18, Star Tracker 92, LIDAR 185, Payload 155).

Figure 22 :

 22 Figure 22: Illustration of the relative position performance for one of the MDO solutions and the low quality configuration. Blue curves represents the filter's covariance (±3σ ) and red curves plot the filter's estimation error.

  . The low quality configuration corresponds to low accuracy sensors (IMU 18, Star Tracker 92, LIDAR 185, Payload 155). a COTS device database (e.g., IMU, Star Tracker, LIDAR, Optical camera, reaction wheels, thrusters).

Table 1 :

 1 Initial orbital elements

			Value	Unit
	Semi-major axis		6884	km
	Eccentricity		0.0012	-
	Inclination		97.5	deg
	Right ascension of the ascending node	233.6	deg
	Argument of perigee		246.3	deg
	Item	Value	Unit
	Volume	58 × 88 × 68	cm
	Mass	140	kg
	Ballistic coefficient	7.9.10 -3	m 2 .kg -1
	Drag coefficient	1.3	-
	Inertia matrix	diag(9,6,9)	kg.m 2

Table 2 :

 2 Relevant chaser satellite characteristics

Table 3 :

 3 First solution of the MDO process : focus on performance

		Component	ID or value
		Sensors' error	11.5
		Mass	4167g
		IMU	33
	Solution 2	Star Tracker	75
		LIDAR	185
		Payload	155
		Reaction Wheel	10
		Thruster	3

Table 4 :

 4 Second solution of the MDO process : focus on mass

Table 5 :

 5 Navigation and payload sub-systems performance at final time for four specific sensor configurations. The best configuration corresponds to the output of unconstrained optimization (IMU 36, Star Tracker 87, LIDAR 186, Payload 171

ACKNOWLEDGMENTS

The authors would like to thank the Master in Astronautics and Space Engineering of Cranfield University (UK) for having provided a grant for the participation to this conference.

with :

being the set of objective vectors that are enclosed by the Pareto front p(A); 1 H(A,R) is the characteristic function of H(A,R); and eventually R is the reference set. This definition is given in [START_REF] Auger | Theory of the hypervolume indicator : optimal µdistributions and the choice of the reference point[END_REF], where more information can be found on hypervolume issues. Numerical settings-The values of the different settings and parameters that come into play for this simulation (equations 9 to 12) are as follow:

𝐹(X)

• The moment arm is evaluated thanks to the satellite dimensions provided in Table 2. Taking the optimum spot to mount the thruster (diagonally) yields a moment arm of L = 0.5m.

• The spinning period t spin in Eq. ( 11), the burn time t burn in Eq. ( 10) and the maneuvering time t man are set to t spin = t burn = t man = 1s.

• The fuel consumption for the orbit correction is 1% of the satellite wet mass, meaning a mass ratio of m f inal /m init = 1.01 in Eq. ( 9).

• The trajectory lasts t acq = 6000s (Eq. ( 12)).

Power and volume have been attributed limits, represented by constraints in the process. The power constraint includes all the devices, whereas the volume constraint only takes into account the sensors. Considering all the previous assumptions on the size and mass of the satellite, we can roughly estimate the typical power that could be available for such a mission, based on data from [START_REF] Wertz | Space Mission Engineering: The New SMAD[END_REF]: for a small satellite of 200W total power consumption, payload and AOCS would account for approximately half of the power consumption. A power constraint of 100W was thus chosen for the simulation. The volume limit is set up to 10000 cm 3 .

Eventually, the population is set to 80 individuals developed over 50 generations (yielding 4000 evaluations).

Results-The Pareto front of the simulation is shown in Figure 18. The two objectives were the minimization of the total mass (payload + AOCS) and the minimization of the AOCS sensors' error. The error of the sensors is represented in Figure 18 as a normalized scalar criterion for the sake of readibility, calculated with equation 38.

where σ p 0 , σ v 0 and σ q 0 are constants corresponding to the smallest values that can be taken by σ p (X), σ v (X) and σ q (X), respectively. As a consequence, E(X) ≥ 1.

The resulting design variables configurations of every final solution are also given in Figure 19. LIDAR, reaction wheel and thruster appear to stay the same (i.e only one component for all the solutions), while IMU, star tracker and payload propose more options.