domain: spi info math phys

text:
Samuel Squillaci
email: samuel.squillaci@onera.fr

Stéphanie Roussel
email: stephanie.roussel@onera.fr

Cédric Pralet
email: cedric.pralet@onera.fr

Managing complex requests for a constellation of Earth observing satellites

Nowadays, the development of Earth observation systems in volving multiple satellites and multiple ground stations al lows Lo satisfy in parallel many observation requests formu lated by end-users. One difficulty though is that in this con text, the mission planning system must often manage various kinds of observation requests (one-shot requests, periodic re quests, etc.), and various kinds of requirements over each re quest (collect data as quickly as possible, download data only towards specific stations, transmit stereoscopic observations during the same downlink window, etc.). This paper intro duces a generic approach to deal with such complex requests, based on so-called modes, where each mode explicitly defines a set of observation and download activities that allows are quest to be fulfilled. The goals arc then to select a subset of requests among the candidate ones, Lo select a mode for each of these requests, and to schedule ail activities covered by the selected modes. For this, the paper first introduces a MILP mode! and then a fast incomplete search procedure based on iterative request insertion and plan repair. The experiments performed show that the approach scales well on realistic in stances involving 16 or 32 satellites.

Introduction

In this paper, we consider an Earth observation system under development that will be composed of about 30 satellites, one objective being to decrease the duration between two possible successive visits of a given ground area. The system will receive observation requests coming from several end users, some ofwhich might have their own reception stations in addition to the main ground station facilities. For this new system, a key point concerns the development of a mission planning tool taking as an input a set of candidate imaging requests and returning as an output a mission plan for each satellite over the next time period.

Severa! techniques were proposed in the last decade to de fine such an observation selection and scheduling tool for a satellite constellation. Sorne authors proposed to use a unique global optimization problem to define the activities of the satellites. Such approaches are based on Branch & Price (Han et al. 2018-12-01;Hu 2019) or MILP (Kim et al. 2019;Xiao et al. 2019). Severa! local search based meth ods were also tested such as tabu search [START_REF] Bianchessi | A heuristic for the multi-satellite, multi-orbit and multi-user management ofEarth observation satellites[END_REF], evolutionary algorithms (Wang et al. 2015; Eddy and Kochenderfer 2020), or simulated annealing (Vongsan tivanich et al. 2018). Pareto Search is also considered when multiple objective functions must be optimized (Sun 2012;Li and Li 2019). Other approaches explicitly identify on one hand an assignment problem that consists in determining the satellite that should be used to perform a given request, and on the other hand a scheduling problem that consists in or dering ail activities assigned to a given satellite. For the as signment part, a first approach selects at each step an ob servation task that maximizes a given scoring fonction and assigns this task to a satellite whose workload is the low est [START_REF] Dishan | A Dy namic Scheduling Method of Earth-Observing Satellites by Employing Rolling Horizon Strategy[END_REF]. Another approach assigns tasks to satellites based on randomized heuristics that take into account various features (size of the observation windows, Nadir position, etc.), and couples this process with an Adap tive Large Neighborhood Search (ALNS) that uses destroy and repair operations to try and decrease the number of un scheduled tasks (He et al. 2018). A last approach assigns each task to the satellite that offers the highest scheduling success probability according to a neural network success prediction mode!, and then performs tabu search to build a plan f r om the task assignment strategy (Du et al. 2020). Fi nally, strategies for on-line planning have also been studied. An approach consists in using heuristics for measuring the difficulty of inserting an observation in the plan of a given satellite (He et al. 2019).

With regards to these previous works, one issue is that for the system under development that we consider, the cen tral mission planner must be able to deal with observation requests that are more complex than those encountered in the literature. More specifically, this planner must be able to deal with several kinds of requests that are relevant for the end-users, including standard one-shot simple observation requests, one-shot stereoscopic requests (two observations of a given ground area, using a forward and a backward pointing respectively), periodic imaging requests (observa tion every X hours over a given a r ea), periodic and stereo scopic imaging requests, etc.

The mission planning system to define must also deal with various constraints and preferences on the way observation and data downlink should be performed. For instance, there might be constraints on the choice of reception stations for downloading collected data. Also, for stereoscopic requests requiring two observations, there might be a need to down load the two associated data files during the same downlink window. As a result, in addition to standard constraints hold ing on each individual activity, there can be constraints hold ing over several activities composing a given request. The same remark holds for preferences in the sense that there might be preferences holding on individual activities (e.g., preference related to cloud cover conditions) or preferences holding on multiple activities (preference for a minimum temporal distance between observation and data transmis sion, preference for the smallest deviation from the periocl for a periodic observation request, etc.).

One contribution of this paper in the frame of Earth Ob servation Satellites (EOSs) consists in the use of so-called request modes to manage such a variety of requirements. Basically, each mode corresponds to a set of observation and download activities that allow to cover a given request white satisfying the request observation and download con straints. To cover complex preferences, a reward is also as sociated with each mode to measure the quality of service obtained for the request with that mode. This reward can aggregate features related to cloud cover forecast or pref erences over the temporal distance between activities cov ered by the mode. One key advantage of this approach is that instead of having to explicitly manage user constraints and preferences, the generic global planning system only needs to reason about the modes and the constraints expressing the physical capabilities of each satellite (kinematic constraints and resource constraints). Doing so, it is easy to integrate new kinds of requests if needed during the lifetime of the system, without changing the core mission planner.

The basic idea of using modes for describing alternatives to carry out an activity is not new in the scheduling domain. For instance, in Multi-mode Resource Constraint Scheduling Problems (MRCPSP (Talbot 1982;[START_REF] Coelho | The Multi-Mode Resource-Constrained Project Scheduling Problem[END_REF]), the goal is to schedule a set of activities over a set of machines that have limited capacities, given that each ac tivity has a set of possible modes. Each mode requires using a subset of the machines during a certain duration and with a certain consumption level for each machine. In MRCPSP, the activities to perform are however not as complex as ob servation requests encountered for EOSs. In another direc tion, the concept of switch groups recently introduced for the Mars 2020 mission also allows to defined several ways of performing an activity [START_REF] Agrawal | En abling Limited Resource-Bounded Disjunction in Schedul ing[END_REF]. Each switch group consists of a set of execution alternatives that con sume more or less time and resources and that are more or Jess preferred, and switch groups can be very quickly used onboard to adapt the plan to the real context.

The rest of the paper is organized as follows. Section 2 describes the problem considered and formalizes the request modes. Section 3 defines a MILP mode) for this problem. Section 4 introduces a fast incomplete search strategy based on iterative request insertion and plan repair. Section 5 pro vides experimental results on scenarios involving heteroge neous requests to be performed by 16 or 32 satellites. Sec tion 6 gives some perspectives for this work.

Problem description

We consider a constellation of Low Earth Orbit satellites that can perform observations and download collected data to ground reception stations. The optical instrument and the data emission antenna are body-mounted on each satellite, meaning that during each observation, a satellite must be pointed to the associated ground target, and during a data download activity, it must be pointed to a reception sta tion. As a result, acquisition and data transfer cannot be per fonned in parallel. In the following, we denote by S the set of satellites, and the objective is to plan the activities of the satellites over a temporal horizon [O, H].

Observation requests

At a given planning phase, the mission centre considers a set of observation requests R. Each request r E R is as sociated with a specific ground area. Satisfying a request means ( I) performing observations over the corresponding ground area based on the available observation opportuni ties, (2) downloading observation data towards ground re ception stations based on the available download opportuni ties. As explained Iater, we consider a set of alternatives to decompose each request into basic observation and down load activities, each alternative being referred to as a mode (more details below).

Observation opportunities

For each satellite s E S, there is a set of observation oppor tunities O ., over relevant ground targets. Each observation oppmtunity o E O s is defined by: • a window [Start 0 , End 0 ] during which the ground target is visible f r om satellite s;

• an observation duration l:,,.

0 E]O, End o -Start o ];
• a download duration 0 0 E ]R+, representing the time re quired for downloading data collected during opportunity o if the latter is used;

• a memory usage Mem 0 E N, representing the memory space used by the data collected if an observation is actu ally performed during opportunity o.

Download opportunities

For each satellites ES, there is a set of download opportu nities V

Past observations

The activities of satellites are regularly planned over a rolling horizon. At a given step, for each satellite s E S, there exists a set of past observations P ., that are stored on board but not downloaded to a ground station yet. These past observations must be taken into account when planning the activities of the satellites over the next scheduling period.

Each past observation o E P s is defined by:

• a download duration 0 0 E ]R + , representing the time re quired for downloading o;

• a memory usage lvlem 0 E N, representing the memo.ry space used for storing o onboard the satellite.

Requests modes

We now describe the link between observation requests, ob servation opp011unities, download opp011unities, and past obst;rvations. Basical ly, for each observation request r E R, there exists a set of candidate modes M r . Each mode rep resents an alternative for completing request r. A mode is a set of pairs ( o, d) where:

• o is an observation opportunity or a past observation con tributing to request r; • d is a download opportunity to be used for downloading data associated with o (with the assumption that dis asso ciated with the same satellite as o). More formally, each mode m E M r is a set of pairs (o, d) E U s E S ( (O s U P s ) x V .., ), with the assumption that an obser vation opportunity or a past observation o does not appear more than once in a single mode (contra.rily to a download opp011unity, that can be used for several observations). Note that an observation or download opportunity can be part of several modes for a given request.

A request r is completed as soon as ail observation and download opportunities associated with one mode m E M. are used for collecting and transmitting data for r. Figure l illustrates modes for a toy example involving two satellites s1 and s2. Satellite s1 bas two observation opportunities (o 1 and 02) and one download opportunity (d 1 ). Satellite s 2 has three observation opportunities (o 3 , o 4 and o 5 ), two down load opportunities (d2 and d3 ) and one past observation (p 1 ). We consider two requests r 1 and r2 that have two modes each:

This means that to satisfy request r 1 , it is possible (1) ei ther to perform an observation within observation opportu nity o 1 and download the associated data within download opportunity d 1 , (2) or to perform an observation within op portunity 04 and download the associated data within op portunity d 3 . Similarly, to satisfy request r 2 , (1) past ob servation Pl must be downloaded during opp011unity d 2 , (2) an observation must be performed during opportunity o 2 and downloaded during opportunity d 1 , and (3) one of the observation-download pairs (o 3 , d 2 ), (o 5 , d 3 ) must be used. As mentioned in the introduction, the concept of mode al lows to caver one-shot and periodic requests, monoscopic and stereoscopic observations, and more generally any con straint over the observation-download pattern .

Last, a reward Wrn E JR + is associated with each mode m E M,.. Such a reward can caver any preferences on fea tures like the weather conditions for the observations, the time at which the observation and download activities take place, the ground station to which data is downloaded, etc. Compared to many models used in the literature, attaching the rewards to modes and not to the basic observation and download activities is more expressive, since the reward of a request might not be the sum of the rewards associated with involving 2 satellites, 2 requests, 5 observation opportunities, 3 download oppOI1unities, and1 past observation its basic activities. One typical example is given by periodic observations that should be separated by a temporal distance that is neither to low nor to high.

A potential drawback of the approach is however that the number of possible modes for a given request rnight be qui te large. For this, we can consider only the most relevant com binations of observation-download pairs for each request, and new modes might be generated on-the-fly if needed.

Maneuver model

Each satellite s performs a sequence of observation and download activities. The temporal distance between succes sive activities must be consistent with the transition time re quired by s to maneuver between the successive target point ings. To model this point, let us denote by A s the set of ac tivity opportunities available for satellite s (A s = C\ U 'D 8 ).

Then, for each pair of activities ( a, b) E A;, we define a minimum transition duration Ta .b E R + that represents the minimum maneuver time from the satellite attitude required for performing a to the satellite attitude required for per forming b.

The previous definition can be extended to the case where a is equal to a fictitious activity o: 8 placed at the beginning of the plan of s, and to the case where b is equal to a fictitious activity f3 s placed at the end of the plan of s.

7 Memory resource model

At any time, the total mass memory used by a satellite must not exceed a limit capacity Mem Cap 8 E N. More formally, for each satellite s and each time t over the scheduling hmi zon, the total memory used by the observations that have been fully performed before t but not fully downloaded yet at t must not exceed Mem Cap 8 • For the system considered, this capacity is not the main bottleneck of the problem.

Solution plan

A solution plan a is defined by : • a subset R(a) Ç R of requests that are selected in a;

• for each request r E R( a), a mode m r (a) chosen for r in a;

• for each satellite s E S, a sequence of opportunities used IIs (a) = [o:s , Ils , l , ... ,IIs ,k ,f3s ] , where each opportu nity II s ,i belongs to A s , and where the set of activities involved in sequences Il s ( ( 7) is consistent with the set of opportunities covered by the modes selected for the re quests in R((7).

From this, the earliest start time of the first activity a = II s ,l of a given satellite sis given by e arlie st 0 _ ( (7) = St art a , and the earliest start time of an activity a = II s ,i that follows an activity b = II s ,i-l is recursively computed by: e arlie st a (OE) = max(St art a , e arlie st b ((7) + dur b ((7) + T&, a )

where durb( (7) stands for the duration of activity b (either the observation duration, when b cor r esponds to an obser vation opportunity, or the sum of the download duration for data downloaded within b, when b corresponds to a down load opportunity).

A plan is said to be feasible if and only if, for each satel lite,

• data observation precedes data download for every observation-download pair ( o, d) covered by a mode used;

• the mass memory capacity constraint is satisfied;

• each activity a can be finished within its allocated time window ( e arlie st a ( OE) + dur a ( (7) :S End a ).

A feasible plan OE is said to be optimal if it maximizes the reward w ( (7) defined as the sum of the rewards provided by the selected request modes (w((7

) = � rER ( rr ) W m ,-(rr ) )
The problem obtained is referred to as an EOSCPM problem (EOS Constellation Planning with Modes).

MILPmodel

In this section, we provide a MILP mode! that can be used to find optimal solutions to small EOSCPM instances.

Additional notations To make the modeling easier, we in troduce a few additional notations:

• 0 = U sES O s denotes the set of observation opportuni- ties;

• D = U sES D s denotes the set of download opportunities; • A = U sES A s denotes the set of activity opportunities; • P = U s ES P s denotes the set of past observations; • M = U,-E n M r denotes the set of request modes;

• for each activity opportunity a E A s for a satellite s, we denote by At (resp. A;;-) the set of activity opportunities that can be performed just after (resp. just before) a in se quence Il s . Variables. We introduce the following variables:

• Vr E R, Vm E M r , Y m E {O, 1} specifies whether mode m is selected for performing request r;

• Va E A:

-Za E {O, 1} takes value 1 if and only if activity oppor tunity a is used; -ta E [ St art 0 , End a] represents the start time of the activity pe1formed during opportunity a; -dur a E [O, End a -St art a ] represents the duration of the activity performed within opportunity a;

• Vs ES, Va E A s U {a 8 }, Vb E At U {/3 s }, x a ,b E {O, 1}
indicates whether b is the successor of a in the sequence of activity opportunities used for s.

MlLP model

The MILP model is given in Equations to 15. The objective fonction (Eq. 1) consists in maximiz ing the global reward provided by the selected modes. Con straint 2 imposes that at most one mode can be selected for each request. Constraint 3 imposes that each activity must be finished before the end of its associated time window.

The duration of an observation is equal to the duration given in the input data (Eq. 4). The download duration within a download opportunity d corresponds to the sum of down load durations for observation opportunities and past ob servations that are associated with d in the selected modes (Eq. 5). Note that contrarily to observation activities, the duration of the download activities is not fixed beforehand. Constraints 6-7 enforce the consistency between the pres ence of activities in the plan (z a variables) and the sequences of activities used for each satellite (x a ,b variables). Con straints 8-9 impose that the initial and final activities for each satellite s are a 8 and /Js. Constraint 10 expresses that an observation opportunity is present if and only one mode containing it is selected. Constraints 11-12 ensure the con sistency between the presence of download opportunities in the plan (z a variables) and the selection of modes (Y m variables): if a mode is selected, then every download op pmtunity of this mode is present, and a download oppor tunity is present in the plan only if at least one mode con taining it is selected. Constraint 13 expresses that the start and end dates of two successive activities must be consis tent with the minimum maneuver times (big-M formulation using constant r a ,b = End a + T a ,b -St art b to ensure that the constraint is inactive when X a ,b takes values 0). Con straint 14 imposes that observation must precede data down load (big-M formulation using constants r a, b again to ensure that the constraint is inactive when Y m takes values 0). Con straint 15 models the memory capacity limitation for each satellite: for each download opportunity d, the memory re quired to store ail observations possibly pe1formed but pos sibly not downloaded yet must not exceed the mass mem ory capacity. Checking the memory capacity between two successive download slots is useless since the memory con sumption peeks are necessarily located at the beginning or during download windows. Note that this constraint is con servative as it uses dates of visibility windows of activities instead of start and end dates decided in the plan.

maximize L L WmYm rER mEM,.

(

Iterative insertion and repair

The MILP mode! defined in the previous section can be used as a reference approach on instances containing a few re quests. For realistic instances however, there !s a nee� for algorithms th at can quickly deliver good-quahty solut10ns. This is why we introduce a fast iterative search procedure. The latter starts from an empty plan and tries to insert at each step one more request r into the plan, using the best candi date mode left for r. Requests are inserted until the so1ution becomes infeasible. In this case, the solution is repaired ei ther by changing the ordering of activities over some satel lites, or by removing requests from the cu1Tent plan. Thi � re pair procedure is app1ied until the solution becomes feas1ble again. To choose a request to inse11 and remove at each step, the alg01ithm respectively uses an insertion heuristic and a repair heuristic. On top of that, a perturbation phase is intro duced to try and find better solutions. The search scheme ob tained can be seen as a Large Neighborhood Search (LNS) that alternates between constructive and destructive search phases. We provide below a more detailed description of each component of the algorithm. The approach directly searches for a plan in the form proposed in Section 2.8 (no manipulation of the MILP variables defined before).

Request insertion procedure

Algorithm 1 gives the basic procedure used to insert one re quest r into a solution plan a using one mode from a set of allowed modes µ r Ç M, .. The algorithm first gets the best mode available in µ r based on function getBestMode (line 2). The latter returns a best mode according to features like the mode reward or the cu1Tent load of each satellite (more details later on these points). The mode chosen for r is stored in a global variable refen-ed to as m ,. (a).

Algorithm 1 then inserts one by one ail activities cov ered by mode m ,, (a). More precisely, for each observat1on download pair ( o, d)

Request removal procedure

At some steps, the search strategy needs to remove requests from the current plan. This is done based on Algorithm 2. The latter ensures that when a request r is removed, ail ob servation and download activities associated with the cur rent mode m ,, (a) of r are removed from the plan. For an observation-download pair ( o, d) in m ,, (a ) , removing the download activity means first decreasing the duration asso ciated with d, and then removing d if its duration becomes null (case in which no more download activity is using d).

Function removeActivity simply removes an activity from the sequence of activities of a satellite, without changing the ordering among the other activities. if durd(a) = 0 then s f-removeActivity(a, d) 6: return a

Plan repair procedure

As mentioned before, after inserting a given request r, the current solution plan CT might become infeasible. In this case, we iteratively use two strategies to repair the solution fol lowing the insertion of r (Algorithm 3).

• The fi r st one performs a local search on the sequences of opportunities used by the satellites whose plan is infeasi ble (line 4 ). This strategy does not change the content of the solution as no activity is removed or inserted.

• The second strategy removes one request from the plan following a repair heuristic denoted H re p (lines 6-8). The algorithm returns the new solution plan CJ" and the set of requests removed f r om the plan to repair the inconsistency. With th is procedure, requests inserted into the plan at a given step can be removed for the sake of the requests that are not satisfied yet.

Iterative insertion and repair

Algorithm 4 gives an iterative algorithm built fro m the pre vious basic components. The algorithm takes as an input a cmTent plan CJ" (possibly an empty plan), an insertion heuris tic H ins , and a repair heuristic H re p. Initially, there is a set of requests R(a) that are already integrated into the plan, and ail requests that are not in R( CJ") are put in the set of requests C that are candidate for inser tion. For every candidate request r in C, ail modes are avail able initially (line 3).

Then, at each step, the algorithm selects one candidate request r in C and inserts this request into the plan, follow

ing insertion heuristic H in s (line 5-6). Request r is removed from the set of candidate requests, and if the solution ob tained is not feasible, it is repaired based on the Repair fonc tion presented before. The latter removes a set of requests E from the cuITent plan to restore consistency (line 9). The current mode m e ( CJ") of each request e E E is then removed f r om the set of modes available for e, and if there is a mode left for e, the latter becomes a candidate for insertion again.

The algorithm can be stopped if a maximum CPU time is reached, and it terminales since at least one more mode is accepted or removed after each repair (and there is a finite set of modes).

LNS algorithm

Finally, Algorithm 5 corresponds to the high-level algorithm that we use for EOS constellation planning with modes. This Algorithm 4 IterativelnsertRepair with H ins = insertion heuristic and H 1 •ep = repair heuristic ] : function ITERATIVE[NS ERTREPAIR(cr, H ins , !-l'" ep ) 2: C +---R\R(cr) 3:

for ail r E C do /h +---Mr 4:

while C /= 0 do 5 :

r +---setectlnsertion(cr, C, H ins ) 6:

s +---insertRequest(cr, r, H ins ) 7:

C+---C\{r} 8:

if -, f easible (cr) then 9:

(cr, E) f---Repair(cr, r, Hre p ) 10:

for all e E [, do 11 : µ,. (c---µ,. \ { me (cr)} 12: if µ,. /= 0 then C +---C U { e} 13: return cr algorithm is a kind of Large Neighborhood Search (LNS), that alternates between constructive and destructive phases. More precisely, the algorithm starts from an empty plan and quickly builds a first solution by using the Itemtivelnserti onRepair procedure seen before together with a simple insertion heuristics H � ns that systematical ly selects a request whose best mode has the highest reward, and a simple repair heuristic H� e p that systematically re moves the last request inserted to repair the plan (line 2).

After that, while there is some CPU time left, the algo rithm tries to apply the iterative insertion and repair pro cedure again, based on inse11ion and repair heuristics H ins and H re p that are typically less greedy than H�n s and H� e p (line 6). In our experiments, H ins is the lowest satellite load heuristic detailed in Section 4.6 and H re p is the highest conjf.ict-index heuristic detailed in Section 4.7.

If the new solution a obtained is strictly better that the best solution found since the last destructive phase, the number of iterations without improvement (it) is reset to 0 (line 8) and the best solution found is updated if needed (line 9). Otherwise, the number of iterations without im provement is updated and if it is greater than parameter maxStablelt given as an input, function Perturb is called to randomly remove some requests from current solution a (line 13). This allows to diversify sea.rch. Finally, the best solution found along ail iterations is returned.

Insertion heuristics

Insertion heuristic H ins mentioned before is used to select a candidate request at each step. We consider two versions of this heuristic.

• Highest reward heuristic: this heuristic chooses a request r whose best allowed mode m offers the highest reward Wm -This criteria does not depend on the cu1Tent state of the solution. CT t-Pert urb(CT , H pert urb ) 14:

(w Siab le , it) f-(w( r7 ) , 0)

15:

return O"• the future. Basically, the average Joad of satellite s at time t in solution cris defined by : '°' dur•a ( Cl )

Loa ds (t, Cl ) = L-End a -Starta ] (Starta ,Enda l (t)

aEA, (a)

In other words, over its allowed time window, each activ ity a consumes a resource ratio that depends on the dura tion of a in a. Without going into further implementation details, the satellite Joad heuristic selects at each step a request r* and a mode m * so as to minimize the height of the maximum peek that is obtained in the satellite Joad profiles after inserting the mode. This criteria depends on the content of the cuITent solution plan O".

7 Repair heuristics

Repair heuristic H re p is used to choose a mode to remove when the current solution becomes infeasible. We consider three versions of this heuristic.

• Highest tardiness heuristic: when the solution becomes infeasible due to temporal constraints, it is possible to de fine , for each activity a of the plan, a tardiness measure T a ( a ) = max ( 0 , earliest a ( a ) + dur a ( a ) -End a )

This tardiness quantifies by how much activity a is late with regards to the end time of its allowed time window. Then, the tardiness of a mode m used in the plan can be defined as:

The highest tardiness heuristic removes a request whose mode rn ,. ( a ) in the current solution has the highest tardi ness.

• Lowest re ward heuristic: among the requests whose mode is involved in the infeasibility explanation ( T m ( a ) > 0), this hemistic removes a request r whose mode m r ( u ) has the lowest reward.

• Highest conj(.ict-index heuristic: instead of measuring the precise tardiness associated with a mode, it is possible to count the number of times each mode occurs in an in feasible part of the solution. To do this, for each activity a= Ils , i of the plan such that T a ( a ) > 0, the conflict-set explaining the tardiness of a can be defined as the short est subsequence exp a ( u ) = [ Ils , i-k, . . . , Il s , i ] such that earliestn . i-k ( Cl ) = Startn . i-k. From this, the conflict index of an' activity a in solution ais defi ned as the num ber of times it belongs to a conflict-set:

Cl a( u ) = card {b E A ( a ) ! a E expb ( a )}

where A ( a ) stands for the set of activity opportunities in cluded in u. Then, the conflict-index of each mode used in the plan is defined by:

The highest conflict-index heuristic removes a request whose mode has the highest conflict-index.

Experiments 5.1 Instances description

The instances generated correspond to constellations com posed of 16 and 32 satellites. The scheduling horizon covers 48 hours, the last 24 hours being dedicated only to down Ioad opportunities. Areas of interest for the requests are ail located in Europe and paititioned into 4 types: ( 1 ) monotonie requests requiring a single observation; (2) stereoscopic re quests requiring two observations performed using two dif ferent painting angle ranges; ( 3 ) systematic re quests requir ing the scheduling of ail possible observations of the target; ( 4 ) periodic requests requüing the scheduling of one obser vation within each period ( 4 periods are considered in the instances ) . We make the number of requests vary from 20 to 500. Each instance is named Rl-R2-R3-R4 where Ri rep resents the number of requests of type i as described pre viously. The number of modes and the size of the modes directly depend on the type of the associated request. In or der to limit the number of modes for each request, we only consider the 10 best modes for each request. Moreover , we only consider one download opportunity for each observa tion. For each request, we randomly generate two bounds a and b in [ D , 1 ] , and a random reward in [ a , b ] is assigned to each mode of this request.

Results analysis

We have implemented the LNS approach in Python 3.6.9 and encoded the MILP approach in CPLEX 20. 1. Both approaches have been run on 20-core Intel ( R ) Xeon ( R ) CPU ES-2660 v3 @ 2.60GHz , 62GB RAM. For the LNS approach, we respectively use the Lowest satellite load heuristic and the Highest conflict-index heuristic as inser tion and repair heuristics ( Hi ns an d H re p ) _ Hperturb consists in removing 5% of the modes. We arbitrarily choose MaxStablelt = 1 and the maximum time for the LNS algo rithm is set to MaxTime = 5 minutes. Maximum time for the MILP approach is set to lO minutes. Table 1: Global reward obtained for each approach along with the number of satisfied requests in parentheses.

Representative results are presented in Ta ble 1. The reward value of the first solution obtained by IterativelnsertRepair (line 2 in Algorithm 5 ) is given in the column named IIR 0 .

A bold font is used for optimal solutions ( MILP approach) .

Results show that the MILP approach is effective on small instances but does not scale up. lndeed, for most instances involving more than 200 requests, the approach returns an empty plan or does even finish pre-processing. ln compari son, LNS produces results for ail generated instances, what ever the types of the requests considered. The global reward varies with respect to the types of requests in the instance since satis fy ing requests having small modes ( monotonie and stereoscopic requests ) is easier than satis fy ing requests having large modes (systematic or periodic requests). For instance, the LNS approach satisfies ail the requests for in stance 250-250-0-0 but not for instance 125-/25-/25-125. Moreover, the global reward increases with the number of observation opportunities ( more satellites ) or with the num ber of non-systematic requests. Solutions returned by the LNS approach are quite close f r om the optimal solutions for small instances and have a better global reward for instances with 200 requests or more. Note that the global reward ob tained by the first iterative insert repair (IIR o ) is improved during LNS. As illustrated in Figure 2a, most requests that are satisfied are fulfilled through one of their best modes but some requests are fulfilled through their 30th best mode. Re quests that are not fulfilled with their best mode have their associated reward decreased. As illustrated in Figure 2b, the reward Joss is smaller than 40% w. 1: t. the best mode.

Figure 3 illustrates the scheduling of 500 requests for a constellation of 16 satellites. The observation, download and

Conclusion and future works

This paper presented a mode-based approach for managing complex requests for a constellation of satellites, together with two algorithmic approaches (MILP and LNS) . Exper iments pe,formed on generated instances show that LNS scales well, which is promising with respect to large cata logs of requests expected for future satellite constellations.

A first perspective of this work is to improve the LNS approach. To do so, it would be possible to limit the num ber of modes and allow on-the-fly creation of new modes for a given request when current modes for this request cannot be inserted. Parallelization techniques could also be used by splitting the scheduling horizon in order to obtain smaller independent problems. A second perspective con sists in making the problem more realistic by considering time-dependent maneuver duration models or by using re ward models that take into account weather forecasts deliv ery time for the observation, requests pri ority, etc. '

Figure 1 :

 1 Figure 1: Toy example involving 2 satellites, 2 requests, 5 observation opportunities, 3 download oppOI1unities, and 1 past observation

  Formally, At = {b E A.!St art a + T a ,b '.S End1,} and A;;-= {b E AslSt art b + T&,a '.S End a }; on the example provided in Figure 1, if we suppose that the transition time between activities is null, A;;. = { d 2 , 00, d 3 } and A;;-. = { 03, d 2, 05 }; • for each mode m, O m and D m denote the sets of obser vation and download opportunities involved in mode m respectively. Formally, O m = {o E Ol:l(o,d) E m} and D m = {d E Dj:l(o,d) E m}; • for a modem E M and a satellites ES, m s is the set of observation-download pairs (o, d) in m that are associated with satellites. Formally, m s = {(o,d) E mlo E 0 5 U P s ,d ED ., }.

  -125-125 239. 16 ( 320 ) 254.95 ( 348 ) 0 (0)

  Loss w.r.t. the best mode

Figure 2 :Figure 3 :

 23 Figure 2: Selected modes for a 500-request instance sat id " l4 l3 i, ll lO -• 1 . -t . . t • . l:H D-I l . ! t -

  8 • Each download opportunity d E V ., corresponds to a time window [Start d , End d ] during which data can be transmitted to a ground reception station.

  covered by the mode, • function greedylnsert is used to insert o into the sequence of opportunities used by the satellite associated with o; this function follows an insertion heuristic H"' 8 (line 4) that selects a position for o in sequence Il s so as to mini rnize the sum of maneuver times for s; • if download opportunity d is already used in the plan, the procedure simply updates the current duration of d in a (line 5); otherwise, download opportunity d is added to the sequence of oppmtunities used by satellite s, based on function greedylnsert ag ain (line 8).

	Algorithm 1 Insert a request
	1: function INSERTREQUEST(u, r, µ r , H ins ) 2: m r (a) f-getBestMode (r, µ, , u)
	3:	for ail (o, d) E m,.(a) do
	4:	u f-greedylnsert(u, o, H' n")
	5: 6:	if d Eu then durd(u) f-durd(a) + 80 else
	7: 8:	durd(u) f-8 0 u f-greedylnsert(a, d, H ins )
	9:	return u

  H�•ep , Hins , H re v, maxSt ablelt , It erativelnsertRepair( O", H in •, H'"eP)

	Algorithm S Large Neighborhood Search
	2:	H pertu ,•b, Max Tim e) O" f-llerat:ivelnsertRepair(©, Ht'", H; e p)
	3: 4: 5:	(T* f-(7 (w Stable , it) f-(w ( CT ), 0) while cpu Time() < MaxTime do
	6: cr f-7: if w(O") > w st<ible then 8: (w Stable , it) f-(w( r7), 0) 9: if w( cr) > w( cr• ) then cr* f-O"
	10: 11 : 12:	else it t-it + 1 if it > MaxStablelt then
	13:
		• Lowest satellite load heuristic: this heuristic takes into ac
		count the Joad of the satellites during the time windows
		associated with the candidate modes of each candidate re
		quest. The main idea is that a low satellite load increases
		the chance to get a feasible solution after inse1tion and de
		creases the chance to get conflicts with modes inserted in

1: fonction LN S(Hr" ,

Acknowledgments

This work has been perfonned with the support of the French government in the context of the "Programme d'lnvestissements d' Avenir", namely by the BPI PSPC pro j ect LiChIE led by AIRBUS Defence and Space and in volving several p artners, including ONERA.


Sophie Thery
email: sophie.thery@univ-grenoble-alpes.fr

On the links between observed and theoretical convergence rates for Schwarz waveform relaxation algorithm for the time-dependent problems

In the framework of linear problems, a usual approach to study the convergence of Schwarz algorithms is to calculate the convergence rate in the frequency domain. However, for time-dependent problems, this tool provides results not fully representative of the observed behaviour of the algorithm. In this article we highlight differences between the theoretical convergence rate and the convergence observed in the physical space-time domain. We also explain how the theoretical convergence rate can be used to provide bounds to the observed convergence rate obs . For problems defined on time windows of finite size, we recall that the bounds usually considered to study the convergence are empirical estimates albeit robust. In conclusion of this paper, numerical experiments are carried out to illustrate the relevance of the theoretical analysis.

       ( , ) -A ( , ) = ( , ) on Ω ×]0, [ B ( , ) = ( , ) on Ω ext ×]0, [ ( , 0) = ,0 ( ) in Ω (1a) C 1,1 1 | Γ ( ) = C 1,2 2 | Γ ( ) on [0, [ C 2,2 2 | Γ ( ) = C 2,1 1 | Γ ( ) on [0, [ (1b)
where can be a finite or infinite time. The Schwarz waveform relaxation algorithm is applied on problem (1a) with interface conditions (1b). For given first guess 0 | Γ ( ) on the interface Γ, the state of the algorithm is given at each iteration ∈ N by [START_REF] Martin | An optimized schwarz waveform relaxation method for the unsteady convection diffusion equation in two dimensions[END_REF].

We suppose here the well-posedness of the initial problem (1) and of the algorithm [START_REF] Martin | An optimized schwarz waveform relaxation method for the unsteady convection diffusion equation in two dimensions[END_REF]. This means there exist a unique solution to (1) in L 2 (0, ; L (Ω )) noted and there exist a unique ∈ L 2 (0, ; L 2 (Ω )) for all iterations . Some results on the well-posedness of such kind of problems can be found in [START_REF] Gander | Optimized schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF][START_REF] Martin | An optimized schwarz waveform relaxation method for the unsteady convection diffusion equation in two dimensions[END_REF] (for problems on finite time window) and in a more general framework in [START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology[END_REF] (for problems on finite or infinite time window).

           1 ( , ) -A 1 1 ( , ) = 1 ( , ) on Ω 1 × [0, [ B 1 1 ( , ) = 1 ( , ) on Ω ext 1 × [0, [ 1 ( , 0) = 1,0 ( ) in Ω 1 C 1,1 1 ( , ) = C 1,2 -1 2 ( , ) on Γ × [0, [ (2a)
           2 ( , ) -A 2 2 ( , ) = 2 ( , ) on Ω 2 × [0, [ B 2 2 ( , ) = 2 ( , ) on Ω ext 2 × [0, [ 2 ( , 0) = 2,0 ( ) in Ω 2 C 2,2 2 ( , ) = C 2,1 1 ( , ) on Γ × [0, [ (2b)
From now on we also suppose ( ) ∈ L 2 (]0, [) for all ∈ Ω . To quantify and possibly optimize the convergence of algorithm [START_REF] Martin | An optimized schwarz waveform relaxation method for the unsteady convection diffusion equation in two dimensions[END_REF], it is relevant to calculate a convergence rate as obs A {1,2} , B {1,2} , C {{1,2},{1,2}} , , = / -1 , where = -| Ω is the error at each iteration . In the rest of the paper, indicies A, B, C are neglected to simplify the notation.

Remark 1

We consider from now on that Ω are one-dimensional domains. Since all convergence factors are calculated in Fourier space, all results explained here can be extended to higher space dimensions parallel to the interface . Also, we consider here Schwarz algorithms applied to multiphysics problems (for nonoverlapping domains) but the following results are also valid in the presence of an overlap.

For example, for parabolic problems we need to have ∈ L 2 (0, ; L 2 (Ω )) and ,0 ∈ L 2 (Ω ) For example, for parabolic problems we need to have ∈ L 2 (0, ; H 1 (Ω )), that it satisfied if and first guess are regular enouth (see [START_REF] Lions | Non-Homogeneous Boundary Value Problems and Applications[END_REF]) This involves applying Fourier transforms in all directions parallel to the interface. Fourier transforms on spatial dimensions do not give rise to the problem that we expose here which is specific to the temporal dimension

Convergence for problems on an infinite time window

We first consider that the simulation is made on an infinite time window, i.e. = +∞.

Convergence factor in Fourier space: For time-dependent problems, the observed convergence factor cannot be calculated analytically. Thus a usual approach consists in applying a time Fourier transform to the error system. In the case where = +∞ and considering that the error is equal to zero for negative times, the convergence is determined in the Fourier space by solving the following system:

       1 ( , ) -A 1 1 ( , ) = 0 on Ω 1 × R B 1 1 ( , ) = 0 on Ω 1 × R C 1,1 1 ( , ) = C 1,2 -1 2 ( , ) on Γ × R (3a)        2 ( , ) -A 2 2 ( , ) = 0 on Ω 2 × R B 2 2 ( , ) = 0 on Ω 2 × R C 2,2 2 ( , ) = C 2,1 1 ( , ) on Γ × R (3b)
Suppose that (3) can be solved for any ∈ R, then the convergence factor and convergence rate in the Fourier space can be calculated as:

( ) := | Γ ( ) -1 | Γ ( ) ( ) := | ( )| (4)
It can be shown that is independent of the space variable and of the domain . Without lack of generality, from now on suppose is calculate from the errors at the interface Γ. General methods to study the convergence of Schwarz algorithms can be found in [START_REF] Carlenzoli | Adaptive domain decomposition methods for advectiondiffusion problems[END_REF][START_REF] Gander | Optimized schwarz methods[END_REF][START_REF] Gander | Schwarz methods over the course of time[END_REF].

Observed convergence factor: From the well-posedness properties of the algorithm, we assume that ( , •) ∈ L 2 (R) for all ∈ Ω . Then ( , •) ∈ L 2 (R) and the following inequality is obviously satisfied:

inf ∈R ( ) -1 ( , •) 2 ≤ ( , •) 2 ≤ sup ∈R ( ) -1 ( , •) 2 (5)
Thanks to Parseval's theorem, we can thus provide the following bounds to the observed convergence factor:

inf ∈R ( ) ≤ ( , •) 2 -1 ( , •) 2 =: obs , ( ) ≤ sup ∈R ( ) . (6)
Thus, if sup

∈R ( ) < 1 then algorithm (2) converges in L 2 (]0, +∞[) norm : ( , •) 2 -→ →∞ 0.
Moreover, because given by ( 4) is the same for all ∈ Ω , previ-ous bounds ( 5) and ( 6) are also valid for obs , :=

L 2 ( ]0,∞[, L 2 (Ω)) / -1 L 2 ( ]0,∞[, L 2 (Ω))
.

Finally the theoretical convergence rate provides bounds for the observed convergence in the L 2 (]0, ∞[, L 2 (Ω )) norms.

Discrete in time problems : Let first assume that we simulate the semi-discrete problem over an infinite time window with a similar time step in both subdomains. The observed numerical error is denoted , ( ) with ∈ N, and can be seen as the result of a Dirac comb on the continuous error ( , •):

, ( ) = , ( ) -( , ) and , • ( ) = Δ ( , •)
with the solution of the discrete problem, the Dirac comb of period and ( , ) the error on the continuous problem. Frequencies higher than / are not generated by the temporal grid [START_REF] Dubois | The optimized Schwarz method with a coarse grid correction[END_REF]. Applying Shannon theorem leads to restrict the study of the errors in Fourier space , • on an interval := -; and

, • ( , ) = ( , ) for all ∈ . Details of the process can be found in [START_REF] Gasquet | Fourier Analysis and Applications[END_REF].

Thus L 2 norm of , • can be calculate using result of the continuous case. Parseval's theorem can be used to obtain bounds on observed convergence rate in L 2 norm:

min | | ≤ / ( ) ≤ obs , ( ) := , • ( ) 2 -1, • ( ) 2 ≤ max | | ≤ / ( ) (7)
Consequently, the same bounds apply on convergence rate obs , in L 2 (]0, ∞[, L 2 (Ω )) norm.

Convergence for problems on a finite time window

Bold notation is used to describe the solution u and the error e over a finite window of time [0, ] with 0 < < +∞. We will consider ∈ with = R if we consider a continuous simulation, and defined in section 2 if we consider a discrete problem.

Difficulties expressing error over a finite time window:

Applying the Fourier transform to the windowed signal would lead to search for the solution of an equation of the type:

e ( , ) + A ( e ( , )) = -e ( , ) exp(- ) (8)
Here the discrete solution (resp. the discrete error ) is obtained by discretizing the continuous solution (resp. the continous error ). The discrete signal optained with a numerical simulation is an approximation of depending on the numerical scheme.

Without more knowledge about the error at time , one cannot solve the differential equation [START_REF] Dubois | The optimized Schwarz method with a coarse grid correction[END_REF], therefore cannot express the error e ( , ) according only to the parameters of the equation. Nevertheless, the error can be expressed by e = * [0, ] ( ) where is the rectangular function on [0, ] and [0, ] = exp(-/2) sinc( /2) The convergence rate for the error e for a given frequency thus reads:

( ) = +1 | Γ * [0, ] ( ) | Γ * [0, ] ( ) = ∫ ( ) | Γ ( ) [0, ] ( -) ∫ | Γ ( ) [0, ] ( -) (9)
which clearly shows that and are different functions, except in the exceptional case where ( ) is a constant. Also, definition [START_REF] Gasquet | Fourier Analysis and Applications[END_REF] supports that function cannot be seen as the convergence factor at a given frequency:

( ) ≠ | e ( )|/| e -1 ( )|.
Bound on observed convergence : The bound on the convergence factor given by ( 9) is complicated to determined. However it is possible to directly bound the error :

Theorem 1 (Bound on the L 2 norm) e ( , •) 2 ≤ sup ∈ ( ) e 0 ( , •) 2 ( 10
)
wich implies a bound for the -product on the observed convergence rate:

Π =1 , ( ) ≤ sup ∈ ( ) ∀ ∈ Ω (11)
This ensures the convergence of the error for the windowed algorithm as long as sup ∈ ( ) < 1. This bounds also works for e ( ]0, [, L 2 (Ω ))

.

Proof It is possible to link the convergence of the windowed problem to a corresponding infinite-in-time problem. We can bound the L 2 norm of error on the windowed problem by the corresponding error of a inifite in time problem :

e ( , •) 2 ≤ ( , •) 2 ≤ sup ∈ ( ) 0 ( , •) 2
where 0 is any first guess extended to infinite time. Using the particular extension 0 | [0, ] = e 0 and 0 | ] ,∞[ = 0 leads to [START_REF] Thery | Analysis of Schwarz Waveform Relaxation for the Coupled Ekman Boundary Layer Problem with Continuously Variable Coefficients[END_REF]. Then combining e ( , •) 2 / e 0 ( , •) 2 = Π =1 , ( ) and ( 10) leads to [START_REF] Gander | Dirichlet-neumann waveform relaxation methods for parabolic and hyperbolic problems in multiple subdomains[END_REF].

Remark 2 A bound on the convergence factor given by ( 9) was already calculated in [START_REF] Thery | Analysis of Schwarz Waveform Relaxation for the Coupled Ekman Boundary Layer Problem with Continuously Variable Coefficients[END_REF]. This bound is complicated to calculate and then hardly usable. In this paper, a global remark on the possible influencing of the time windowing is done. It is explained why the method used in [START_REF] Gander | Optimized schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF][START_REF] Gander | Dirichlet-neumann waveform relaxation methods for parabolic and hyperbolic problems in multiple subdomains[END_REF] needs special conditions and cannot be applied in a general context .

Range of influencing frequencies:

For a given problem discretized in time over a time window [0. ], we estimate that in the general framework: min

/ ≤ | | ≤ / ( ) ≤ obs , ( ) ≤ max / ≤ | | ≤ / ( ) (12)
and the interval | | ∈ , is called influencing frequencies. This interval of frequencies is usually considered for optimizing the convergence rate. As discussed in section 2, it is justified to consider that / is the maximum frequency. However the choice of the minimum frequency is justified for time-independent problems but is an empirical estimate for time-dependent problems. We can still find justifications for this choice by considering the definition [START_REF] Gasquet | Fourier Analysis and Applications[END_REF]. First, it shows that convergence is influenced by ∫ ( ) more than by its value at a given frequency . Moreover, thanks to the property of ( -), frequencies such that | | / have a low impact on the convergence . That said, relevance of the minimum influncing frequency / is to be proved.

Remark 3 (minimal frequency for time-independent problem)

The frequency min is justified in some cases of time-independent problem. If the conditions on border of the dimension parallel to the interface (the one where the Fourier transform is made) are determined, then the corresponding error system is periodic and a Discrete Fourier Transform (DFT) can be applied (for example see [START_REF] Gander | Optimized schwarz methods without overlap for the helmholtz equation[END_REF]). In our case, we can apply a Fourier transform on our discretised signal but, for the reasons evoked section 3, we cannot guarantee that | (e +1 ) ( )|/| (e ) ( )| is equal to ( ).

Remark 4 (optimization)

Usually the optimisation of the convergence speed is made by choosing interface conditions C { {1,2}, {1,2} } under such conditions , such that inf

C {{1,2},{1,2}} ∈ max / ≤ | | ≤ /
C {{1,2},{1,2}} ( ). From [START_REF] Gander | Optimized schwarz methods without overlap for the helmholtz equation[END_REF] this guarantee an minimal upper bound to the observed convergence rate and consequently a fast convergence of

L 2 (0, , L 2 (Ω))
to zero.

Numerical illustration

We propose to illustrate the previous properties on the coupling of two diffusion equations, with Dirichlet-Neumann interface conditions with a non-overlapping interface in = 0. it requires the determination of the inverse Fourier transform F -1 1 which may not exist or can be hard to calculate for a given , frequencies such that | -| / are drown in the integration in ( 9)

       ( , ) - 2 ( , ) = 0 on Ω ×]0, [ ( , ) = 0 on Ω ext ×]0, [ ( , = 0) = 0 on Ω (13a) 1 (0, ) = 2 (0, ) on [0, [ 2 2 (0, ) = 1 1 (0, ) on [0, [ (13b)
with

Ω 1 = [ℎ 1 , 0] and Ω 2 = [0, ℎ 2 ].
We simulate the problem via an implicit finite difference scheme. Schwarz's algorithm on this problem is performed on 20 iterations, with a time step = 1000s and parameters ℎ 1 = -50m, ℎ 2 = 300m, 1 = 0.12m 2 s -1 and 2 = 0.6m 2 s -1 . In figures 1 and 2 we compare the theoretical convergence rate given in the Fourier domain ( ) with the observed convergence rate

, and the convergence rate measured on the DFT of the error at the interface | (e +1 ) (0, )|/| (e ) (0, )| which can be seen as an approximation of ( ). Fisrt guesses are initialised by a random signal which generates a large frequency spectrum. In these two figures, we find that bounds of the observed convergence verify the estimate (12) thus also verify theorem 1. The evolution of the L 2 norm of the error is not explicitly given here but it can be deduce from , (middle panel in 1 and 2). As expected the convergence observed on a given frequency does not correspond to the theoretical convergence ( ) and conversely we tend towards equality for a window of assumed size infinite. Other examples on such problems were made in [START_REF] Thery | Etude numérique des algorithmes de couplage océan-atmosphère avec prise en compte des paramétrisations physiques des couches limites[END_REF] and corroborate the estimate.

Conclusion

In the context of a time dependent problem, the convergence rate calculated in the Fourier space can only be taken as such on problem considering an infinite time window. Thanks to Parceval theorem, informations on the algorithm in the physical space can be obtain on the L 2 norm of the error. It is therefore possible to bound the observed convergence rate obs with the bounds of the theoretical convergence . For a finite time window, we can no longer consider as a convergence rate for a given frequency. Yet, bounds on the observed convergence rate are still relevant and we can precise these bounds by estimating an interval of influencing frequencies. In a futur work, it may be relevant to determine how to choose optimized interface conditions using the results on the observed convergence rate.

Fig. 1

 1 Fig. 1 For a finite time window with = 200 . Left panel: theoretical convergence rate ( ), influencing frequencies are given by vertical lines and grey zones give the reached values of obs , . The observed convergence factor obs , is given in the middle panel as a function of the iteration number for the two domains. Right panel: the observed rate | (e 1 ) (0, ) |/| (e -1 1 ) (0, ) | is compared to the theoretical convergence rate for the first four iterations.

Fig. 2

 2 Fig. 2 Same as Figure 1 with = 10 5 . It is considered to be close to an infinite time window.


domain: spi phys

text:
Keywords: DC/AC converters, Power Electronics, Conducted electromagnetic emission, Thevenin generalized model

This article deals with circuit modeling of power systems and more particularly of DC/AC converters used for controlling electric motors. The objective of the model is to be able to assess conducted emission currents at both inputs and outputs in order to apply appropriate mitigation design rules such as filtering. First, we propose a linear Thevenin block-model of the converter considering it as a 5-port linear black-box allowing possible interaction between the DC and AC ports. To this extent, the block-model is made of an impedance matrix and a Thevenin voltage generators vector. We also propose theoretical relationships allowing the processing of data in order to be able to generate this model from measurements. Then, we test the linear simplifications on a real DC/AC converter and compare the results to measurements. For this, S-parameter measurements in offmode and currents at the 5 ports in on-mode are used for the characterization of the model on a first test setup. The model is then applied and checked by comparing circuit simulated and measured current responses on a modified test setup for which the lengths of the connection cables are modified. Finally, the capability of this block-model to extrapolate by calculation the response of other installation configurations is showed considering the design of a filter and its insertion at input/output ports of the converter with much longer cables.

I. INTRODUCTION

ew technologies for automobile electric traction and aeronautics electric propulsion are increasingly integrating electronic power systems in on-board equipment for the conditioning of electric motors. Most of electronic power equipment uses switching technologies, which makes it possible the conversion of a direct current (DC) at the input into a three-phase alternative current (AC) at the output for driving and controlling the electric motors. Unfortunately, DC/AC High Voltage (HV) power converters are sources of significant EMI (Electromagnetic Interference) signals propagating in the Electrical Wiring Interconnection System (EWIS). The resulting currents are known in the Aeronautics EMC community as "conducted-emissions" currents. Characterizing those currents in the frequency domain is the specific subject of section 21 (Emission of Radio Frequency Energy) of the L. Guibert, J-P. Parmantier, I. Junqua and M. Ridel are with ONERA/DEMR, Université de Toulouse, F-31055 Toulouse -France (email: laurent.guibert@onera.fr) aeronautic equipment standard RTCA D0-160 [START_REF]Environmental Conditions and Test Procedures for Airborne Equipment[END_REF] for qualification of on-board equipment. It is therefore imperative to be able to control the level of unwanted currents, for example by inserting appropriate filtering devices [START_REF]IEEE Standard for Methods of Measurement of Radio-Frequency Power-Line Interference Filter in the Range of 100 Hz to 10 GHz[END_REF]. In aeronautics as in the whole transportation industry, the main challenge of high-power filtering generally results in a trade-off between performance and additional mechanical weight.

For an optimal EMC-design accounting for all vehicle integration constraints, being able to model DC/AC devices in the frequency domain and to include this model in an overall circuit model is thereby essential in order to assess the impact of the equipment insertion in various installation configurations. Modelling the full power electronic chain from DC generation to AC activation together with the installation conditions allows for example optimization of filter designs, being able to account for the entire wiring topology. From this perspective, sophisticated filtering techniques can be derived ranging from lumped filters combined with specific cables [START_REF] Jullien | Developement De Harnais De Cable En Vue De L'Optimisation Cem D'une Chaine De Conversion De Puissance Pour L'Aeronautique[END_REF] to active filters switching at specific frequencies [START_REF] Baisden | High frequency modeling of a converter with an RF-EMI filter[END_REF]. A significant research activity is carried out on this topic. However, because the functioning of HV DC/AC converters relies on nonlinear circuitry, full modelling of such devices including at the same time the HVDC inputs and HVAC outputs is a real challenge due to their nonlinear nature. Most converter EMC characterization references investigate "black-box" models [START_REF] Foissac | Black box' EMC model for power electronics converter[END_REF]. The models generally focus on EMI sources and lead to explicit equivalent circuit models of converters, which makes this black-box model concept not "fully black". Besides, the circuit models generally address the problem from the two HVDC input ports for given conditions of the HVAC output only ( [START_REF] Wang | Model conducted EMI emission of switching modules for converter system EMI characterization and prediction[END_REF], [START_REF] Liu | Frequency-domain EM1 Noise Emission Characterization of Switching Power Modules in Converter Systems[END_REF]) or eventually the three HVAC output ports for given conditions of the HVDC ports only [START_REF] Costa | Modeling of conducted common mode perturbations in variable-speed drive systems[END_REF] . Nevertheless, the full modelling of both HVDC and HVAC ports is not addressed at the same time, especially because of the difficulty to establish an explicit and complete equivalent circuit model accounting for the interaction between the HVDC and HVAC ports.

In this article, we thereby want to investigate the possibility to establish a 5-port "fully-black" model for which no knowledge on the DC/AC internal circuit components is required. If applicable, such a decomposition would not have any limitations on the number of ports and would be general enough to include the efficient Differential Mode (DM) and

N

Common Mode (CM) decomposition models as done for HVDC inputs in usual references ( [START_REF] Foissac | Black box' EMC model for power electronics converter[END_REF], [START_REF] Wang | Model conducted EMI emission of switching modules for converter system EMI characterization and prediction[END_REF], [START_REF] Liu | Frequency-domain EM1 Noise Emission Characterization of Switching Power Modules in Converter Systems[END_REF], [START_REF] Amara | Control of conducted emissions of power electronics[END_REF]). Particularly, the full 5-port characterization of a DC/AC would include implicitly the possible reaction between the HVDC input and HVAC output ports, which is not considered in usual practice but which is required for optimized mitigation of EMI once the power chain is installed in a vehicle with its specific geometry. For this, such block models could thereby easily be introduced in larger circuit models for which load conditions as well as connection cables and possible filters could be modified as desired, either at the HVDC input or at the HVAC output.

Such models of course exist for real linear circuits. In this article, despite the highly nonlinear nature of the DC/AC converters, we want to challenge the relevance of deriving such a linear block model from basic measurements only, based on very simple hypothesis. For this purpose, we consider a generalized Thevenin block-model, classically made in our case of a 5x5-impedance matrix and a 5-component voltagegenerator vector [START_REF]Chapitre 2 : Couplage aux structures filaires[END_REF]. As far as the input impedance matrix is concerned, we do not want to investigate the topology of the circuit model inside the black-box: the impedance matrix terms are measured at the 5-port level with a VNA (Vectorial Network Analyzer). Moreover, we want to make the crude a-priori assumption that the impedance can be measured when the system is "off" in order to allow very simple measurement processes. As far as Thevenin generator components are concerned, we want to measure them directly in the frequency domain with a power spectrum analyzer in order to avoid time domain measurements and the subsequent Fourier transform, as in previously mentioned references ( [START_REF] Foissac | Black box' EMC model for power electronics converter[END_REF] to [START_REF]Chapitre 2 : Couplage aux structures filaires[END_REF]). All our a-priori assumptions will be validated a-posteriori for specific power operation conditions by a simulation in which the DC/AC linear block-model will be introduced. The comparison with frequency domain measurements results will allow us to determine the application domain of the converter model and the required modelling practices for this applicability.

From this perspective, this article is organized as follows. In section II, we present the generalized Thevenin model theoretical concepts and relationships that will be used to postprocess specific measurements for the establishment of the block-model. In section III, we apply the methodology to characterize the block-model of a laboratory test-bench consisting of a real 5kW DC/AC converter and its associated loads. We apply the relationships introduced in section II and we evaluate the interactions between the HVDC input and the HVAC output ports. We then assess the robustness of the derived block-model by comparing calculated and measured conducted-emission currents on a specific test-bed configuration. In section IV, we illustrate how this model can be used to define an appropriate filter in order to be compliant with a given conducted emission standard. For this, we introduce the DC/AC block-model into a Multiconductor Transmission Line Network (MTLN) model allowing the management of complex cable bundle models and linear blockmodels [START_REF] Parmantier | CRIPTE : Code de réseaux de lignes de transmission multiconducteur[END_REF].

II. EQUIVALENT BLOCK-MODEL AND MATHEMATICAL PROPERTIES

A. Presentation of the equivalent linear block-model

In order to be representative of our DC/AC multiport circuit problem, we present here our Thevenin linear block-model problem made of 5 ports: 2 ports for the HVDC input and 3 ports for the three-phase HVAC output (Fig. 1). We consider a model for which all voltages at all ports are defined with respect to the same electrical ground reference [START_REF]Chapitre 2 : Couplage aux structures filaires[END_REF]. components Ii, as usually done in linear circuit theory [START_REF] Wing | Classical Circuit Theory[END_REF], at each port i. Note that the Thevenin model is based on the hypothesis of linearity and independence of sources, which implies the theorem of superposition. This means that the [Z], [Y] and [S] matrices are independent from the Thevenin generators. At this stage of our analysis, we may ask ourselves if those matrices must be determined setting the power of the multiport either on "on" or "off". However, there is no doubt that the determination of the Thevenin generator requires the power of the multiport to be "on".

We can then group together all the voltages, currents and voltage generators at the 5 ports in vectors:

[𝑉] 5 = [ 𝑉 1 𝑉 2 𝑉 3 𝑉 4 𝑉 5 ] , [𝐼] 5 = [ 𝐼 1 𝐼 2 𝐼 3 𝐼 4 𝐼 5 ] , [𝐸 𝑡ℎ ] 5 = [ 𝐸 𝑡ℎ , 1 𝐸 𝑡ℎ , 2 𝐸 𝑡ℎ , 3 𝐸 𝑡ℎ , 4 𝐸 𝑡ℎ , 5 ] (1)

B. Relationships between impedance, admittance and Sparameter matrices

In the following of the article, we will determine the [Z] matrix from S-parameter measurements. We thereby remind hereafter some interesting transformation properties between [Z], [Y] and [S] matrices, together with their application domains [START_REF]Chapitre 2 : Couplage aux structures filaires[END_REF]. Whatever the number of ports of the multiport system, [Z], [Y] and [S]-matrices are defined with respect to the voltage and current vectors, [V] and [I], as defined in (1):

{ [𝑉] = [𝑍]. [𝐼] [𝐼] = [𝑌]. [𝑉] [𝑉] -[𝑍𝑟𝑒𝑓]. [𝐼] = [𝑆 𝑍𝑟𝑒𝑓 ]. ([𝑉] + [𝑍𝑟𝑒𝑓]. [𝐼]) (2)
[SZref] is the S-parameter matrix for which [Zref] is the impedance matrix on which the S-parameters are determined (or "referenced"). By elimination of the [V] and [I] vectors in the two first equations of (2), the third equation then makes it possible to determine various relationships between those matrices. Among them, we will focus on two particular ones:

[𝑌] = [𝑍] -1 (3)
and

[𝑍] = ([𝐼𝑑] -[𝑆 𝑍𝑟𝑒𝑓 ]) -1 . ([𝐼𝑑] + [𝑆 𝑍𝑟𝑒𝑓 ]). [𝑍𝑟𝑒𝑓] (4)
where [Id] is the identity matrix.

Note that, whereas the [Z] and [Y] matrices are not always defined (for respectively open and short circuit conditions), the [SZref] matrix has the property to be always defined, provided the [Zref] matrix is itself defined. On the other hand, the [SZref] matrix depends on the reference [Zref] matrix whereas [Z] and [Y] matrices do not depend on any reference matrices. This is the reason why, in the following of this paper, we will write down all [S] matrices as [SZref] to avoid any ambiguity on their reference impedance. Finally, note that, unlike [Z] and [Y] matrices, [S] matrices are not generally symmetric. They are symmetric only for very specific conditions of [Zref], for example when [Zref] is diagonal with all diagonal terms equal to each other. This is for instance the case for S-parameters referenced on 50 for which [Zref]=50. [Id]. This is of particular interest for measurements since the VNA provides the 50 loading on its measurement ports and the possibility of very large-frequency band measurements.

C. Determination of the S-parameter matrix

From the third equation of (2), we can show that the SZref,i,j matrix terms of the j h column of the [SZref] matrix are mathematically defined as follows, for the off-diagonal terms:

𝑆 𝑍 𝑟𝑒𝑓 ,𝑖,𝑗 𝑗≠𝑖 = 2.𝑉 𝑖 𝐸 𝑗 ) [𝑉]+[𝑍 𝑟𝑒𝑓 ].[𝐼]=[0],𝑒𝑥𝑐𝑒𝑝𝑡 𝑓𝑜𝑟 𝑗 (5)
and for the diagonal terms:

𝑆 𝑍 𝑟𝑒𝑓 ,𝑗,𝑗 = 2.𝑉 𝑗 𝐸 𝑗 -1) [𝑉]+[𝑍 𝑟𝑒𝑓 ].[𝐼]=[0],𝑒𝑥𝑐𝑒𝑝𝑡 𝑓𝑜𝑟 𝑗 (6)
where Ej is the voltage generator exciting port j, all other ports ij being with no voltage generator and, Vi and Vj the resulting voltages at ports i and j.

From an electrical circuit point of view, the condition, "[V]+[Zref].

[I]=[0] except for j", in ( 5) and ( 6), can be obtained by loading all ports of the multiport system by [Zref]. Such a method can be applied for measuring any [SZref] matrix terms applying sequentially a generator Ej to each port j of the blockmodel to be characterized. For instance, for determining the 3 rd column of the [SZref]5x5 matrix in Fig. 1, the block-diagram described in Fig. 2 can be applied. In this figure, the voltage generator E3 is applied on port 3 of the equivalent block-model. Same types of equations as ( 5) and ( 6) would apply for the determination of the j h column terms of [Z] and [Y] matrices. Nevertheless, for the Zi,j terms of the column j, all ports except port j have to be left open. For the Yj,j terms of the column j, all ports except port j have to be set in short-circuit. Consequently, the drawback of such methods is that the conditions of loads on the ports depend on the matrix column to be determined. On the contrary, the interest of S-parameters is that the load condition remains the same on all the ports whatever the column to be determined. This is a significant advantage from a practical point of view when the determination of the [S Zref ] matrix is obtained from laboratory measurements on a large number of ports. Note also that the same advantage also happens in numerical techniques for embedding sub-networks as blackboxes with large numbers of ports (see for example [START_REF] Parmantier | Various Ways to Think of the Resolution of the BLT Equation with an LU Technique[END_REF] for MTLN subnetwork compaction applications).

D. Determination of the Thevenin generator vector

The Thevenin generator vector of the block-model is independent of the load conditions of the multiport system. Indeed, its terms can also be defined as the induced open-circuit voltage generators. The vector can therefore be deduced from measurements of voltages on any [Zo] matrix. As said before, unlike for S-parameter determination, the power of the multiport system must this time be turned "on" which generally implies to load the ports by a [Z o ] matrix defining the operating conditions of the system. In Fig. 3

III. TEST-BENCH STUDY-CASE

A. Approach and objective

The objective is now to establish the block-model of a real DC/AC HV converter from measurements using the relationships and load conditions defined in section II. Then we will measure the currents of a specific DC/AC converter setup and we will build the block-model-based circuit of this setup. The calculated currents will be compared to direct measurements. This comparison will have the intent to validate our model as well as the hypothesis made on the derivation of the linear block-model.

B. Presentation of the laboratory test-bench

In the context of EMC in power electronics, we have developed a 5kW laboratory test-bench from the Commercial Off-The-Shelf (COTS) elements of a real high-performance electric mountain biking system. Fig. 4 presents a picture of this test-bench consisting of two brushless types of motors mounted on wheels coupled by a transmission belt system (Fig. 4).

The operating principle is as follows: when one of the two motors receives its power through the DC/AC converter called "BMSWC" (Brushless Motor Sinus Wave Controller), the other motor serves as a variable load operating in alternator mode. Fig. 4 shows in the left-hand side foreground the converter with its two LISNs (Line Impedance Stabilization Networks) in the right-hand side connected to the two converter HVDC inputs. In the background a 72 VDC power supply is connected to the inputs of the LISN for the complete HVDC conditioning of the test-bench.

C. Presentation of the test-bench electrical network model

Fig. 5 presents the electrical scheme of the test-setup. The DC/AC is the 5-port system to be characterized. It is made of 2 HVDC input ports and 3 HVAC output ports. The electrical connection between the HVDC input and the LISN can be ensured by two types of 2-wire cables: cable n°1 is a short-length cable (1.0 m), and cable n°2 is a long-length cable (5.0 m). For the connection between the HVAC output and the motor, two types of 3-wire cables (3 phases) are also possible: cable n°3 is a short-length cable (1.0 m), and cable n°4 is a long-length cable (5.0 m).

The electrical network block-model equivalent to Fig. 5's electrical circuit is represented in Fig. 6. The circuit blockmodel equation writes as follows:

[ 𝐸 1 𝐸 2 𝐸 3 𝐸 4 𝐸 5 ] = ( [ [ 𝑍_ 𝑆ℎ𝑜𝑟𝑡_ 𝐶𝑎𝑏𝑙𝑒_ 𝐿𝐼𝑆𝑁 ] 2𝑥2 [0] 2𝑥3 [0] 3𝑥2 [0] 3𝑥3 ] + [ 𝑍_ 𝐵𝑀𝑆𝑊𝐶_ 𝐶𝑜𝑛𝑣𝑒𝑟𝑡𝑒𝑟_ 𝐻𝑉𝐷𝐶_ 𝐻𝑉𝐴𝐶 ] 5𝑥5 + [ [0] 2𝑥2 [0] 2𝑥3 [0] 3𝑥2 [ 𝑍_ 𝑆ℎ𝑜𝑟𝑡_ 𝐶𝑎𝑏𝑙𝑒_ 𝑀𝑜𝑡𝑜𝑟 ] 3𝑥3 ] ) * [ 𝑖 1 𝑖 2 𝑖 3 𝑖 4 𝑖 5 ] (8)
[Z_Short _Cable_LISN]2x2 is the 2x2 impedance matrix of the 2-port LISN system seen from the HVDC input of the converter, thereby taking into account the 2x2 impedance matrix of cable n°1 connected to the 2x2 impedance matrix of the two LISNs. The 3x3 [Z_Short _Cable_Motor]3x3 impedance matrix is seen from the HVAC output of the converter, thereby taking into account the 3x3 impedance of cable n°3 together with the 3x3 impedance matrix of the motor. Finally, the [Z_BMSWC_Converter_HVDC_HVAC]5x5 impedance matrix is seen from the 5 input and output terminals of the converter.

In [START_REF] Costa | Modeling of conducted common mode perturbations in variable-speed drive systems[END_REF], the building of the matrices is such that the impedance matrices of the LISN and the motor are added in series to the two HVDC input ports and to the three HVAC output ports of the DC/AC converter respectively. To this extent, the 2x2 and 3x3 matrices are inserted as upper and lower diagonal blocks in a null 5x5 matrix ([0]nxm is a null block with n lines and m columns). Then, we make our simplification hypothesis that the impedance of the DC/AC can be determined in "off" mode. This hypothesis implies that the impedance matrix can be considered as unchanged when high power is delivered in real operating conditions with [Z_Short _Cable_LISN] 2x2 and [Z_Short _Cable_Motor]3x3 loading connections. Thereby, the [Z_BMSWC_Converter_HVDC_HVAC]5x5 matrix is also obtained from the derivation of S-parameters measured on the 5 ports of the DC/AC converter loaded on 50 loads, Therefore, [Zref]2x2=50.[Id]5x5 is applied in transformation (4) on the 5x5 [S] matrix. Specific adaptation devices had to be developed on purpose for applying the required 50 loading conditions and mechanical connection of the two VNA measurement ports to the two DC/AC HVDC input-ports connected to cable n°1 on the one hand and the 3 DC/AC HVAC output-ports connected to cable n°3 on the other hand. Electric cable connections have been made as short as possible with respect to the mechanical constraints of the test-setup in order to be as close as possible to the ports under measurement (Fig. 7). Note also that since we used a two-port VNA, particular attention has been brought to load the ports not connected to the VNA to a 50 resistance load. This measure is absolutely required in order to ensure the theoretical coherence of the [S] matrix. Forgetting the 50 loads on "VNA-unused" ports is a common mistake in such measurements. The interest of the Y-parameters is that they can be easily mapped to the admittance values Y e ij of the parallel circuit elements of the equivalent circuit network seen from the 5 ports. Such circuit elements can be obtained by deriving the [Y]5x5 matrix from which they can be easily extracted with the following equations [START_REF] Wang | Model conducted EMI emission of switching modules for converter system EMI characterization and prediction[END_REF]:

𝑌 𝑖𝑗,𝑖≠𝑗 𝑒 = -𝑌 𝑖𝑗 ( 9
)
𝑌 𝑖𝑖 𝑒 = ∑ 𝑌 𝑖𝑗 5 𝑗=1,5 (10)
In Fig. 8, typical admittance values Y e ij are plotted as samples of:

-common mode admittance on HVDC ports (Y All the plots display high resonance peaks that must be exclusively attributed to the electrical behavior of the converter (switching commutations). 𝑌 𝑖∈ [START_REF]Environmental Conditions and Test Procedures for Airborne Equipment[END_REF][START_REF]IEEE Standard for Methods of Measurement of Radio-Frequency Power-Line Interference Filter in the Range of 100 Hz to 10 GHz[END_REF],𝑗∈ [START_REF] Jullien | Developement De Harnais De Cable En Vue De L'Optimisation Cem D'une Chaine De Conversion De Puissance Pour L'Aeronautique[END_REF][START_REF] Foissac | Black box' EMC model for power electronics converter[END_REF] 𝑒 admittance values between input and output ports have significant amplitude, which proves the existence of this interaction.

E. Derivation of the DC/AC block-model source vector

The derivation of the Thevenin voltage generators at the level of the 5 ports of the converter is obtained from the measurement of the 5 currents i1 to i5 on the 5 ports of the converter, with respect to the block-model network represented in Fig. 6. The setup used applies short-cable connections for both HVDC (cable n°1) and HVAC (cable n°3) ports. The derivation is directly obtained from relation [START_REF] Costa | Modeling of conducted common mode perturbations in variable-speed drive systems[END_REF].

Currents in Fig. 5's test-setup are measured with a current clamp, compatible with our wide band analysis and connected to a spectrum analyzer. Unfortunately, the spectrum analyzer only provides amplitudes whereas our block-model-based electrical circuit modelling technique explicitly requires phase management. To this extent, we had to make another straightforward assumption: we post-processed all measured current amplitudes enforcing a-priori and logically expected phase shifts, i.e.:

-for HVDC, Arg(i1)=0 and Arg(i2)=π -for HVAC, Arg(i3)=0, Arg(i4)=2π /3 and Arg(i5)=4π/ 3The results of the Thevenin generators resulting from the derivation of ( 8) are presented in Fig. 9 and Fig. 10. It is important to observe that equivalent voltage sources of similar amplitudes are derived on both HVDC input ports and HVAC output ports.

F. Verifications of the DC/AC block-model

As a matter of fact, the [i]5x1 vector figured in [START_REF] Costa | Modeling of conducted common mode perturbations in variable-speed drive systems[END_REF], with the components i1 to i5, represents the conducted emission currents obtained in a specific configuration when cable n°1 and cable n°3 are used. In order to verify the robustness of our DC/AC block-model, we propose a simulation consisting in calculating those currents in a configuration different from the configuration used for building the converter block model.

The simulation consists in computing the conducted emission currents when we substitute short cables by long cables, i.e. cable n°1 by cable n°2 and cable n°3 by cable n°4. For the derivation, we will use the following relationship, directly updated from (8) with [Z_Long _Cable_LISN]2x2 and [Z_Long _Cable_Motor]3x3 matrices and solve the circuit model in Fig. 11. We obtain: _Cable_Motor]3x3 for cable n°4 loaded by the motor are obtained from VNA S-parameter measurements on 50 in the same way as before, for the short-cables configuration. The results of Fig. 12 and Fig. 13 present the comparison between calculated and measured current amplitudes. Despite some deviations in the high frequency regime, above some few MHz, that we can attribute to the limited precision of the calibration due to our finite size measurement loading device, the agreement is quite satisfactory with respect to the various approximations made in the DC/AC converter block model derivation. Consequently, those results may be considered as fully acceptable for the evaluation of the conducted emission currents. Even if they do not answer the question of the independence of the Thevenin model with the delivered power, they seem to confirm the fact that the determination of the impedance of the converter in "off" mode is acceptable for the derivation of the impedance block-model. In this section, we want to illustrate the interest of such a block-model strategy in terms of modelling the entire system circuit chain. For this, we have chosen the preliminary design of a filter to be introduced at the HVDC input and HVAC output of the converter in order to comply with a given standard level. We have arbitrarily chosen a conducted-emission level (Power Lines -Cat H) provided in section 21 of the aeronautics EMC DO-160G standard [START_REF]Environmental Conditions and Test Procedures for Airborne Equipment[END_REF]. This standard defines the maximum total current (complex sum of the 2 DC port currents on the one side and the 3-phase currents on the other side respectively) not to be overcome on the cables connecting the HVDC input to the RSIL and the HVAC output to its load, between 100kHz and 200MHz. The usual practice to mitigate those currents consists in the insertion of multiport filters in the input and output circuits. In aeronautics, oversizing the filtering performance must be absolutely avoided because the use of filters generally comes together with a significant increase of weight. This is why the design of the filter must account for the converter performance together with all other parameters that also contribute to the current response and the filtering effect. Here, we mainly think of the installation configuration (type of cables, length of cables and cable running paths in to the 3D hosting structure) and end-loads. To this extent, we have used a MTLN modelling tool based on the resolution of the Baum-Liù-Tesche (BLT) equation [START_REF] Parmantier | CRIPTE : Code de réseaux de lignes de transmission multiconducteur[END_REF] that allows us to use S, Z and Yparameter block matrices in large circuit topologies involving large varieties of cable-bundle models and EM driven current and voltage generator models.

[ 𝑖 1 𝑖 2 𝑖 3 𝑖 4 𝑖 5 ] = ( [ [ 𝑍 𝐿𝑜𝑛𝑔 𝐶𝑎𝑏𝑙𝑒 𝐿𝐼𝑆𝑁 ] 2𝑥2 [0] 2𝑥3 [0] 3𝑥2 [0] 3𝑥3 ] + [ 𝑍_ 𝐵𝑀𝑆𝑊𝐶_ 𝐶𝑜𝑛𝑣𝑒𝑟𝑡𝑒𝑟_ 𝐻𝑉𝐷𝐶_ 𝐻𝑉𝐴𝐶 ] 5𝑥5 + [ [0] 2𝑥2 [0] 2𝑥3 [0] 3𝑥2 𝑍 𝐿𝑜𝑛𝑔 𝐶𝑎𝑏𝑙𝑒 𝑀𝑜𝑡𝑜𝑟 3𝑥3 ] ) -1 . [ 𝐸 1 𝐸 2 𝐸 3 𝐸 4 𝐸 5 ] (11)
For the demonstration, we have considered a "Very long cable" configuration in which a 5m-long cable and a 30m-long cable have been connected to the HVDC input and HVAC output ports respectively. At both extremities, those cables have been connected to the previously defined short-cable loads (based on cable n°1 and cable n°3).

First, in Fig. 14 and in Fig. 15 we show the total common mode currents at the HVDC input and HVAC output of the converter, computed by the MTLN tool, including the blackbox model and the different cables, and compared to the DO-160G standard threshold in 3 different configurations:

• With no additional cables,

• With a 30m-long cable at the HVAC output only • With additional cables at both the HVDC input (5mlong cable) and HVAC output (30m-long cable) The figures show the impact of the additional cables together with the DO-160G threshold level. This threshold level is clearly overcome below 10MHz and at specific frequencies over 10 MHz for all configurations. The influence of the long cables combined with the existence of significant interaction between the HVDC input and HVAC output of the transformer signs on the current responses with a larger average level below 1 MHz at the HVDC input and additional frequency peaks (at about 300 kHz, 2 MHz for example) at both the HVDC input and the HVAC output. Second, we have considered the generic filter circuit template model represented in Fig. 16 and we have studied its impact when it is inserted either on the HVAC output or both at the HVDC input and the HVAC output. The values of the circuit element have been set in such a way to lower the total current level under the DO-160G standard, over the whole frequency range of interest, at both the HVDC input and the HVAC output. In Fig. 17 and Fig. 18, we show the common mode currents obtained when 2 additional cables are inserted at input and output of the converter with and without filters and compare them to the DO-160G standard threshold level. When the filter is at the HVAC output only, the attenuation is obviously obtained at the output but there is a significant impact on the HVDC input. The only way to attenuate the total current at the HVDC input is to use another filter at this input as well (we made the hypothesis of the same filter design). However, we note in Fig. 18 that the total current, even if still acceptable, is modified again at the HVAC output. Such results give an idea of the complexity for optimizing filters while accounting for all installation parameters as well as the characteristic of the DC/AC converter. Of course, the optimized design of a real filter should lead to more sophisticated circuit designs, both at the converter HVDC input and the HVAC output. However, our first-level approach shows that any design must account for a precise knowledge of the impedance matrix of the converter and, especially, of the coupling between the HVDC input and HVAC output ports of the converter.

V. CONCLUSIONS

We have presented a frequency-domain generalized Thevenin block-model of a DC/AC converter. This model has been tested and validated on a 5kW-72V DC driven DC/AC converter, using electric COTS components of an electric bicycle. The model has been established considering several a-priori simplified hypotheses that have all proved to be successful in our a-posteriori verifications, which makes our experimental model characterization really efficient by its simplicity, i.e:

-

The system is supposed to be linear which makes it possible the characterization in the frequency domain of its input impedance circuit block model as Z-Y-Sparameter matrices, -The [Z] matrix has been derived from matrix transformation of VNA-measured S-parameters on 50, the converter being in "off" mode, making the assumption that our model is independent from the delivered power.

-

The Thevenin generator vector model has been derived from electrical-current measurements when the DC/AC converter is in "on" mode and when its HVDC input ports and HVAC output ports are loaded on two LISNs and a 3-phase electric motor respectively, both of them being connected to short-length cables (1 m). The amplitudes of those currents have been measured with a spectrum analyzer and with current probes suitable in the 1kHz -400MHz frequency band. The phases of the signal have been appropriately and ideally assigned to the current amplitudes in order to build complex number currents.

At the end, those simplifications allowed us to establish a simple block-model, as we would have done it for a regular linear circuit. The robustness of this model has been verified considering a setup with longer cables. The comparison between calculated and measured currents proved to be very satisfactory. In particular, we put to the fore significantly large input-output interaction terms in the [Y] matrix. The confidence in the block-model finally allowed us to show the interest of block-models for design purposes in order to define the characteristics of a filter while accounting for installation cables. Our design had the simplified objective to maintain the conducted emission current level under a given value that came in our case from the RTCA DO-160G aeronautical standard.

In the future, the conclusions of this paper will have to be further confronted with different types of DC/AC converters and loads, especially for confirming the independence of the DC/AC block-model with operating power. Indeed, the evaluation and the validation showed in this paper have been performed for a specific converter and specific power configuration imposed by the same two LISNs and the same electrical motor as for the characterization of the converter model. Indeed, even if the model has been tested considering electrical cable connections of various types, those configurations do not significantly modify the absorbed power due to their limited losses. Measurement with real filters as the one defined in simulation might give a first answer to this question.

To this extent, the property that our method is based on the measurement of S-parameters referenced to any [Zref] matrix appears as a real advantage. Indeed, for a multiport system requiring continuous polarization or nominal operating power, [Zref] can be chosen as the load ensuring the suitable operating power condition. The application at ports of the voltage generator required for the S-parameter injection can be performed by a transformer (typically a current clamp) [START_REF] Parmantier | S-parameter determination with a pair of current injection and measurement probes[END_REF], which allows the user to keep as such the loading by the [Zref] reference impedance that remains connected during the whole measurement process while imposing the required power operating conditions. By this way, the [S] matrix determination also becomes a fully non-intrusive method which is not the case with a VNA direct-port connection. Besides, [Zref] can be also used to determine the Thevenin generator vector components. Nevertheless, compared to a direct voltage generator application as done in this paper with a VNA (with the requirement of 50 port loading conditions), the use of such transformation devices following procedures as indicated in [START_REF] Parmantier | S-parameter determination with a pair of current injection and measurement probes[END_REF] may limit the capability of wide band analysis. For instance, current clamps have generally a high-pass behavior and the size of those current clamps limits precise calibration at high frequencies.

Fig. 1 .

 1 Fig. 1. Generalized Thevenin block-model, in the case of a 5port multiport system. The equivalent block-model is composed of the following components: -The impedance matrix, [Z]5x5 (or its equivalents Yparameter admittance [Y]5x5 matrix or S-parameter [S]5x5) matrix. -The 5 Thevenin voltage generators, Eth,i, at each port i of the multiport system. -The 5 voltages, Vi, and 5 incoming current vectorcomponents Ii, as usually done in linear circuit theory[START_REF] Wing | Classical Circuit Theory[END_REF], at each port i. Note that the Thevenin model is based on the hypothesis of linearity and independence of sources, which implies the theorem of superposition. This means that the [Z], [Y] and [S] matrices are independent from the Thevenin generators. At this stage of our analysis, we may ask ourselves if those matrices must be determined setting the power of the multiport either on "on" or "off". However, there is no doubt that the determination of the Thevenin generator requires the power of the multiport to be "on".We can then group together all the voltages, currents and voltage generators at the 5 ports in vectors:

Fig. 2 .

 2 Fig. 2. Principle of measurement of the [SZref] matrix referenced on [Zref]: determination of the 3 rd column of Fig. 1's [S] matrix.

  , we show the block-modelbased equivalent circuit with Fig. 1's multiport system loaded by the [Zo] matrix block-model. Under these loading conditions, we can write the simple matrix relation of the circuit in a compact matrix form: [𝐸 𝑡ℎ ] = ([𝑍] + [𝑍 𝑜 ]). [𝐼] (7) Thereby, the determination of the components of the [I] current vector delivered on the [Zo] matrix impedance load provides, by means of the knowledge of the [Z] matrix of the multiport and the [Zo] matrix impedance, the components of the [Eth] Thevenin vector.

Fig. 3 .

 3 Fig. 3. Principle of measurement of the Thevenin generator vector on a [Zo] impedance load block-model in Fig. 1's setup example.

Fig. 4 .

 4 Fig. 4. The 5kW laboratory test-bench setup.

Fig. 5 .

 5 Fig. 5. Test-setup electrical scheme of the 5-port system DC/AC converter

Fig. 6 .

 6 Fig. 6. Block-model-based electrical network with the converter in "short-cable" configuration.

Fig. 7 .

 7 Fig. 7. S-parameter measurement setup. On the right handside, the adaptation device allowing connection of the VNA measurement-ports and the 50 loading of "VNA-unused" ports. Y-parameters of the DC/AC converter are derived from Zparameters as follows: -Calculation of the Z-parameter [Z]5x5 matrix from the [SZref]5x5 matrix as in (2) -Calculation of the Y-parameter [Y] 5x5 matrix by inverting the [Z]5x5 matrix as in (3) Note that the numerical processing for the derivation of the [Y] matrix preserves the symmetry of the [Y] matrix and that the 2 input DC ports and the 3 outputs AC ports have similar behaviors.

e 11 )

 11 -common mode admittance on HVAC ports (Y e 33) -admittance between the 2 input HVDC ports (Y e 12 ) -admittance between the 3 output HVAC ports (Y e 34) -admittance between an input and an output ports (Y e 15)

Fig. 8 .

 8 Fig. 8. DC/AC-converter admittance values Y e ij associated to electrical circuit derived from S-parameter measurements.

Fig. 9 .

 9 Fig. 9. Thevenin generator source terms of the HVDC input of the DC/AC converter.

Fig. 10 .

 10 Fig. 10. Thevenin generator source terms of the HVAC output of the DC/AC converter.

Fig. 11 .

 11 Fig. 11. Block-model-based electrical network with the converter in "long-cable" configuration The new impedance matrices [Z_Long _Cable_LISN]2x2 for cable n°2 loaded by the two LISNs and [Z_Long

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Long cable configuration: comparison between measurements and calculations of the conducted emission currents at the HVDC input.

Fig. 14 .

 14 Fig. 14. Very long cable configuration: RTCA Section 21 DO-160G (Power Linescat H), conducted emission threshold compared with the HVDC total currents.

Fig. 16 .

 16 Fig. 16. Generic filter circuit template.

Fig. 17 .

 17 Fig. 17. Filtering effect on very long cable configuration: RTCA Section 21 DO-160G (Power Linescat H), conducted emission threshold compared with the HVDC total currents.

ACKNOWLEDGMENT

This work has been supported and sponsored, partly with internal resources from ONERA and partly by the French DGAC MECEP SDC 2016 93 0800 CORAC project.