Sanjoy K Baruah
email: baruah@cs.unc.edu

Vincenzo Bonifaci
email: vincenzo.bonifaci@iasi.cnr.it

Renato Bruni
email: bruni@diag.uniroma1.it

Alberto Marchetti-Spaccamela
email: alberto@diag.uniroma1.it

ILP models for the allocation of recurrent workloads upon heterogeneous multiprocessors

Keywords: task partitioning, sporadic tasks, unrelated machines, speedup bound, ILP rounding 1 Introduction

come

common case. The presence of such heterogeneity requires choices to be made when mapping software components onto processing elements. The need to resolve such choices adds considerable complexity to resource allocation, and inhibits the adoption of such platforms by the embedded computing industry despite significant potential benefits in terms of balancing performance and energy.

We consider here real-time systems that are modeled as collections of independent sporadic tasks (the model is described in detail in Section 2). We seek to devise algorithms for implementing such systems upon heterogeneous multiprocessor platforms under the partitioned paradigm. To our knowledge, this topic has not been studied much previously:

-On the one hand, most prior real-time scheduling research that considers heterogeneous platforms (see, e.g., [START_REF] Chwa | Optimal real-time scheduling on two-type heterogeneous multicore platforms[END_REF][START_REF] Raravi | Task assignment algorithms for two-type heterogeneous multiprocessors[END_REF][START_REF] Raravi | Assigning real-time tasks on heterogeneous multiprocessors with two unrelated types of processors[END_REF][START_REF] Raravi | A PTAS for assigning sporadic tasks on two-type heterogeneous multiprocessors[END_REF][START_REF] Wiese | Partitioned EDF scheduling on a few types of unrelated multiprocessors[END_REF]; [START_REF] Raravi | Real-Time Scheduling on Heterogeneous Multiprocessors[END_REF] has a nice survey) has restricted attention to implicit-deadline sporadic tasks. -On the other hand, prior research that does address the partitioned scheduling of task systems represented using models that are more general than the implicit-deadline model considers identical multiprocessor platforms only (see, e.g., [START_REF] Baruah | The partitioned multiprocessor scheduling of deadline-constrained sporadic task systems[END_REF][START_REF] Chen | Resource augmentation for uniprocessor and multiprocessor partitioned scheduling of sporadic real-time tasks[END_REF]), or restricts the amount of heterogeneity, in the sense that the models allow only a very small number of types of processors (see, e.g., [START_REF] Andersson | Scheduling constraineddeadline parallel tasks on two-type heterogeneous multiprocessors[END_REF]).

In this paper, we initiate a methodical study of the problem of partitioning, upon arbitrarily heterogeneous multiprocessor platforms, task systems that are represented using the constrained-deadline sporadic task model 1 . We assume that once the partitioning has been performed and tasks assigned to the processors, run-time scheduling is done on each processor using the earliest deadline first (EDF) scheduling algorithm, which is known to be optimal for this purpose [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard real-time environment[END_REF][START_REF] Dertouzos | Control robotics : the procedural control of physical processors[END_REF]. On the other hand, we remark that the task partitioning problem subsumes the NP-hard unrelated machines scheduling problem [START_REF] Lenstra | Approximation algorithms for scheduling unrelated parallel machines[END_REF][START_REF] Shmoys | An approximation algorithm for the generalized assignment problem[END_REF].

Our approach. We will derive various approaches to task partitioning. These algorithms share the commonality that they are all based upon formulating the task partitioning problem as an integer linear program (ILP).

For implicit-deadline task systems, this is not particularly difficult to do; indeed most of the research on partitioning implicit-deadline sporadic task systems (including the works [START_REF] Raravi | Task assignment algorithms for two-type heterogeneous multiprocessors[END_REF][START_REF] Raravi | Assigning real-time tasks on heterogeneous multiprocessors with two unrelated types of processors[END_REF][START_REF] Raravi | A PTAS for assigning sporadic tasks on two-type heterogeneous multiprocessors[END_REF][START_REF] Wiese | Partitioned EDF scheduling on a few types of unrelated multiprocessors[END_REF] cited above) has been based upon first formulating such an ILP, and then seeking polynomial-time algorithms for obtaining approximate solutions to these ILPs (solving an ILP is known to be NP-hard [START_REF] Karp | Reducibility among combinatorial problems[END_REF], and hence unlikely to be solvable exactly in polynomial time). Despite this inherent intractability of solving ILPs, however, the optimization community has recently been devoting immense effort to devise extremely efficient implementations of ILP solvers, and highly-optimized libraries with such efficient implementations are widely available today. Modern ILP solvers, particularly when running upon powerful computing clusters, are often capable of solving ILPs with tens of thousands of variables and constraints. We therefore believe that it is reasonable to attempt to solve ILPs exactly rather than only approximately, and seek to obtain ILP formulations that we will seek to solve exactly to solve the partitioning problem for constrained-deadline sporadic task systems. Since the running time of ILP solvers tends to increase with the number of variables and constraints in the ILP to be solved, we seek to develop ILPs for task partitioning in which the number of variables and constraints are restricted to be low-order polynomials of the representation of the task system. While the number of constraints may not always be a good indicator of the complexity of an ILP formulation, we use it as a first approximation: indeed, the best known complexity bounds for solving ILPs do increase with the number of linear constraints [START_REF] Frank | An application of simultaneous diophantine approximation in combinatorial optimization[END_REF]Theorem 5.3]. Possibly more refined metrics, such as the constrained induced-width [START_REF] Freuder | A sufficient condition for backtrackbounded search[END_REF] or the constraint density [START_REF] Achlioptas | On the solution-space geometry of random constraint satisfaction problems[END_REF], have also been suggested in other settings, but not in the context of the problem of partitioning tasks onto heterogeneous processorsnot even implicit-deadline tasks [START_REF] Kamath | Unrelated parallel machine schedulingperspectives and progress[END_REF]; these refinements fall outside the scope of this work.

of presentation we do not explore this issue any further in this paper, but leave it for future work.

Our results. In partitioning implicit-deadline sporadic task systems, an ILP represents an exact solution to the partitioning problem -solving an ILP exactly therefore constitutes an optimal algorithm for performing such partitioning. For partitioning constrained-deadline systems, however, we do not know how to obtain such an exact ILP representation of this problem with only polynomially many constraints -this difficulty was previously identified even for partitioning upon identical multiprocessors in [START_REF] Baruah | Partitioned scheduling of sporadic task systems: an ILP-based approach[END_REF]. Instead, our goal here is to obtain polynomial-sized ILP representations of the problem of partitioning constrained-deadline sporadic task systems upon heterogeneous multiprocessor platforms with the property that exact solutions of the ILP constitute approximate solutions to the partitioning problem. Our metric of effectiveness of such an approximate solution is the speedup factor -an ILP formulation has speedup factor f , f ≥ 1, if any constrained-deadline sporadic task system that can be partitioned upon a particular heterogeneous platform by a hypothetical optimal algorithm can be partitioned using this ILP formulation upon a platform in which each processor is at least f times as fast.

We have derived several different ILP representations for the problem of partitioning constrained-deadline sporadic task systems upon heterogeneous multiprocessor platforms, all of which have number of variables and constraints polynomial in the representation of the task system. All these ILP formulations have nm integer variables, where n is the number of tasks and m the number of processors, but specify different numbers of constraints and offer different speedup guaranteesthey are summarized in Table 1.

Outline of the paper. The rest of the paper is structured as follows. In Section 2, we formally define our model and introduce some notation. In Section 3, we present a first ILP formulation for the partitioning problem, with a guaranteed speedup bound of 4, and show one possible generalization. In Section 4, we consider a strengthened ILP formulation that trades off the number of constraints with the speedup guarantee. Section 5 discusses a variant of the ILP of Section 3 that is used as the basis of a polynomial-time partitioning algorithm with a guaranteed speedup bound of 7.83, which improves the current bound of 12.9 [START_REF] Marchetti-Spaccamela | Assigning sporadic tasks to unrelated machines[END_REF]. In Section 6, we report the results of a large number of experimental results that test our ILP formulations on synthetically generated instances. We give our conclusions in Section 7.

Model Number of constraints Speedup factor Comments

1. n + m + m log 2 d max 3 In Section 3.2 1'. n + m + m log ρ d max 1 + ρ In Section 3.3. ρ is a constant > 1. A generalization of Model 1 (which is obtained if ρ = 2) 2. n + m + n m k 1 + 1 k
In Section 4. k is any integer ≥ 1.

3.

n

+ m + m log ρ d max 2 + ρ + ρ 2 ρ-1
In Section 5. A "poorer" version of Model 1' -same number of constraints, larger speedup factor. But more amenable to polynomial-time approximation -see Section 5 for details

System model, background, and notation

We seek to partition a sporadic task system τ comprising the n independent sporadic tasks τ 1 , τ 2 , . . . , τ n upon an unrelated multiprocessor platform π comprising the m processors π 1 , π 2 , . . . , π m . The i'th sporadic task τ i is characterized by a period p i and a relative deadline d i , and m worst-case execution time (WCET) parameters c i,1 , c i,2 , . . . c i,m , with c i,j denoting the WCET of τ i if it executes upon processor π j . In this paper, we restrict attention to task systems in which d i ≤ p i for each task τ i ∈ τ -such sporadic task systems are called constrained-deadline sporadic task systems. During run-time, τ i releases a sequence of jobs at timeinstants that are not known beforehand, but with the constraint that successive jobs are released at least p i time units apart. Each job of τ i has a deadline d i time units after its release time; the amount of execution needed by this job depends upon the identity of the processor on which it executes. More specifically, since we are studying partitioned scheduling in this paper, given task system τ and multiprocessor platform π, our objective is to obtain a partitioning of the tasks upon the processors. Let f (i) ∈ {1, 2, . . . , m} denote the index of the processor to which each τ i is assigned under such a partitioning; each job of τ i needs to execute for up to c i,f (i) time units upon processor π f (i) .

Since the preemptive Earliest Deadline First scheduling algorithm (EDF) is known to be optimal for scheduling a single preemptive processor, we will use EDF as the scheduling algorithm upon each individual processor during run-time. The demand bound function (dbf) [START_REF] Baruah | Preemptively scheduling hard-real-time sporadic tasks on one processor[END_REF] of a sporadic task is widely used to quantify the computational demand of such a task, where the dbf(τ i , t) of sporadic task τ i with period p i , relative deadline d i , and WCET c i for an interval of duration t is defined as follows:

dbf(τ i , t) := t + p i -d i p i c i .
It is known that a collection of sporadic tasks can be scheduled to always meet all deadlines upon a preemptive uniprocessor by EDF if and only if for all t ≥ 0, the sum of the dbf's of all the tasks in the collection for an interval of duration t does not exceed t. Some additional notation:

-Let N := {1, 2, . . . , n} denote the task index set.

-Let M := {1, 2, . . . , m} denote the processor index set. -Let u i,j := c i,j /p i denote the utilization of task τ i on processor π j . -Let d max := max 1≤i≤n {d i } denote the largest relative deadline parameter of any task in τ .

As above, let f : N → M denote a partitioning of the task system τ upon multiprocessor platform π. We use the notation dbf f,j (t) to denote the sum of the dbf's of all the tasks in τ that have been assigned to processor π j under the partitioning f , for interval duration t:

dbf f,j (t) := i∈N :f (i)=j t + p i -d i p i c i,j .

A simple ILP model for task partitioning

Marchetti-Spaccamela et al. have previously [START_REF] Marchetti-Spaccamela | Assigning sporadic tasks to unrelated machines[END_REF] derived a polynomial-time algorithm for assigning sporadic tasks to heterogeneous processors, and shown that this algorithm has a speedup bound of (8 + 2 √ 6) or ≈ 12.9. An intermediate step in [START_REF] Marchetti-Spaccamela | Assigning sporadic tasks to unrelated machines[END_REF] is the derivation of a 0/1 ILP representation of the partitioning problem with nm variables and (n + m + m log d max ) linear constraints, and a proof that this ILP representation has a speedup factor of 6. In this section, we present two results:

1. We derive, in Theorem 1 below, an improved ILP with the same number of variables and constraints and show (Corollary 2) that it has a superior (i.e., smaller) speedup factor of 3. 2. In Theorem 2, we generalize the derivation of this improved ILP in the following manner. For any constant ρ > 1, we can derive an ILP with nm variables and (n + m + m log ρ d max ) linear constraints and speedup bound 2ρ; by choosing ρ to be smaller than two, a smaller speedup bound than 3 is thus obtained at the cost of needing to solve an ILP with a larger number of constraints.

3.1 A review of some results from [START_REF] Marchetti-Spaccamela | Assigning sporadic tasks to unrelated machines[END_REF] First a preliminary definition. Let D denote the set of values 0, 1, 2, 4, . . . , 2 log 2 dmax .

We call D the deadline checkpoint set for task system τ .

The starting point of the reasoning in [START_REF] Marchetti-Spaccamela | Assigning sporadic tasks to unrelated machines[END_REF] is the following lemma.

Lemma 1 (from [START_REF] Marchetti-Spaccamela | Assigning sporadic tasks to unrelated machines[END_REF]) Let f : N → M denote an assignment of the tasks in task system τ to the processors of unrelated multiprocessor platform π such that, for each j ∈ M i:f (i)=j u i,j ≤ β, and for each j ∈ M and k, 1

≤ k ≤ log 2 d max , i:(f (i)=j)∧(2 k-1 <di≤2 k ) c i,j ≤ β • 2 k .
Then for each j ∈ M , dbf f,j (t) ≤ 6βt for all t ≥ 0.

Let us try and understand what this lemma means.

The first set of inequalities requires that the cumulative utilization assigned to each processor not exceed β; the second, that the sum of the WCETs of all tasks assigned to a processor that have relative deadlines between two successive powers of two not exceed β times the larger power of two. (For example, considering k = 6, the second constraint mandates that the sum of the WCETs of all tasks with relative deadline in the range (32, 64] not exceed 64β.) The lemma asserts that if these conditions are satisfied by an assignment, then this assignment constitutes a valid partitioning upon processors of speed 6 × β. Now, suppose that τ is feasible upon π under partitioned scheduling, i.e., there is an assignment f : N → M of τ upon π such that all jobs of all tasks will always complete by their deadlines if tasks are assigned according to this partitioning, and each processor scheduled during run-time by EDF. For this assignment f , it is evident that the utilization constraints of Lemma 1 are satisfied for β = 1. The second set of constraints in Lemma 1 will also be satisfied for β = 1, by the following reasoning:

-Since the partitioning f is feasible, the sum of the dbf's of all the tasks assigned to the jth processor for interval duration

2 k is no more than 2 k . -Each task with relative deadline in (2 k-1 , 2 k ] must have dbf(τ i , 2 k ) ≥ c i,j .
-Hence the sum of the c i,j 's for all tasks τ i that have relative deadline in (2 k-1 , 2 k ] and are assigned to the jth processor must be ≤ 2 k .

Hence if τ is feasible upon π, there exists an f : N → M for which the conditions of Lemma 1 are satisfied with β = 1; therefore dbf f,j (t) ≤ 6t for all t ≥ 0 for each j.

Thus [START_REF] Marchetti-Spaccamela | Assigning sporadic tasks to unrelated machines[END_REF] derives the following speedup result:

Corollary 1 Let f : N → M denote an assignment of tasks to processors satisfying the conditions of Lemma 1.

Then f constitutes a feasible partitioning of τ upon π under a speedup factor of 6β. In particular, if the conditions of the lemma are satisfied with β ≤ 1/6, then assignment f constitutes a feasible partitioning of τ upon the given platform π.

In summary, the conditions of Lemma 1 are necessary when β = 1 and sufficient when β = 1/6, i.e., they yield a speedup factor of 6 for partitioning constraineddeadline sporadic tasks upon heterogenous multiprocessor platforms.

An ILP with speedup factor 3

We now prove an improved version of Lemma 1 that yields a superior speedup bound (of three, rather than six). The conditions specified in our theorem below differs from those in Lemma 1 only in that the summation of c i,j 's in the second inequality is over all tasks with relative deadline ≤ 2 k (rather than only those with relative deadline > 2 k-1 and ≤ 2 k ).

Theorem 1 Let f : N → M denote an assignment of the tasks in task system τ to the processors of unrelated multiprocessor platform π such that, for each j ∈ M

i:f (i)=j u i,j ≤ β (1)
and for each j ∈ M and each k,

1 ≤ k ≤ log 2 d max , i:(f (i)=j)∧(di≤2 k ) c i,j ≤ β • 2 k . ( 2
)
Then for each j, 1 ≤ j ≤ m, dbf f,j (t) ≤ 3βt for all t ≥ 0.

Proof Consider any t ≥ 0 and let s := 2 k denote the smallest integer power of 2 that is not smaller than t (i.e. s = 2 k ≥ t > s/2).

Consider the assignment f : N → M of the tasks to the processors as in the hypothesis of the theorem. For any j ∈ M , we have

dbf f,j (t) = i:f (i)=j∧di≤t t + p i -d i p i c ij ≤ i:f (i)=j∧di≤t t c ij p i + c ij ≤ t i:f (i)=j c ij p j + i:f (i)=j∧di≤s c ij = t i:f (i)=j u ij + β • 2 k ≤ βt + βs ≤ 3βt,
and the theorem is proved.

This implies that any f satisfying Inequalities 1 and 2 of Theorem 1 constitutes a feasible partitioning of the tasks in τ upon the set of heterogeneous processors π, when the processors receive a speedup factor of 3β.

Corollary 2 Let f : N → M denote an assignment of tasks to processors satisfying Inequalities 1 and 2 of Theorem 1. Then f constitutes a feasible partitioning of τ upon π under a speedup factor of 3β. In particular, if Inequalities 1 and 2 are satisfied with β ≤ 1/3, then assignment f constitutes a feasible partitioning of τ upon the given platform π.

We now apply the result of Theorem 1 above to construct a 0/1 integer linear program (ILP) for specifying a feasible partitioning of sporadic task system τ upon the platform π. That is, we will construct a 0/1 ILP, a solution to which will yield a partitioning f : N → M that satisfies Inequalities 1 and 2. This ILP is constructed as follows:

-For each i ∈ N, j ∈ M , we have a 0/1 integer variable (i.e., an integer variable that takes on either the value zero or the value one) x i,j , denoting whether τ i is to be assigned to processor π j .

-We specify that each task gets assigned to exactly one processor; this is done by the following n constraints:

∀i ∈ N j∈M x i,j = 1 (3)
-We next specify that Inequality 1 of Theorem 1 should be satisfied; this is achieved by the following m constraints:

∀j ∈ M i∈N x i,j u i,j ≤ β (4)
-Finally, we specify Inequality 2 of Theorem 1 by the following (log d max × m) constraints:

∀k ∈ D, ∀j ∈ M (i∈N )∧(di≤2 k ) c i,j x i,j ≤ β • 2 k (5)
By Theorem 1, any assignment of integer values to the {x i,j } variables satisfying the Constraints 3-5 above bears witness to the schedulability of τ , with a speedup of 3β, upon the platform π. Moreover, for a τ that is feasible upon π the model always admits a solution with β = 1, since all inequalities are clearly valid; thus, this ILP guarantees a speedup factor of at most 3. When the ILP model is feasible with β ≤ 1/3, Theorem 1 guarantees schedulability on the original platform; hence, a reasonable objective function for the ILP with Constraints (3)-( 5) would be to minimize β.

Throughout the rest of this paper, we will still refer to the problem of the minimization of β subject to (3)-( 5) as ILP, even if β is not strictly constrained to be integer. However, we adopt this terminology for extension, and also because β could equivalently be fixed, for example in a binary search fashion.

The ILP model consisting of Contraints (3)-( 5) will be referred to as Model 1 in the remainder of the paper.

A generalization

Above, we derived an ILP model for the problem of partitioned scheduling of constrained-deadline sporadic task systems upon unrelated multiprocessors, such that any solution to the ILP immediately yields a partitioning algorithm at a speedup bound of 4. We also saw that this ILP has (n + m + m × log 2 d max ) constraints; we now briefly describe how to reduce the speedup bound by increasing the number of constraints.

Recall that we had defined the deadline checkpoint set D as powers of two: D = {0, 2 0 , 2, 2 2 , . . . , 2 log 2 dmax }. For any given constant ρ > 1, we could instead have chosen to define it as

D ρ = {0, 1, ρ, ρ 2 , ρ 3 , . . . , ρ log ρ dmax },
The following generalization of Theorem 1 is easily proved via a proof analogous to the proof of Theorem 1:

Theorem 2 Let f : N → M denote an assignment of the tasks in task system τ to the processors of unrelated multiprocessor platform π such that, for each j ∈ M i:f (i)=j u i,j ≤ β, and for each j ∈ M and k,

1 ≤ k ≤ log ρ d max , i∈N :(f (i)=j)∧(di≤ρ k ) c i,j ≤ β • ρ k Then for each j ∈ M , dbf f,j (t) ≤ (1 + ρ)βt for all t ≥ 0.
It follows, from arguments virtually identical to those of Corollary 2, that the speedup bound for the ILP constructed based on Theorem 2 above is (1 + ρ); hence by choosing ρ to be smaller than 2 a speedup bound smaller than 3 is obtained. The tradeoff is that the number of constraints increases to (n+m+m× log ρ d max ); this is > log 2 d max for ρ < 2, becoming larger as ρ → 1.

Theorem 2 thus suggests one approach for obtaining speedup bounds arbitrarily close to 2, by simply selecting a denser deadline checkpoint set. In Section 4 below, we explore another approach, that allows for speedup bounds arbitrarily close to one (once again at the cost of having additional constraints).

A strengthened ILP formulation

We now explore a different idea that also trades off an increase in the number of constraints in the ILP for a superior speedup bound. Specifically, for any positive integer constant k we will derive an ILP model with (n + m + m n k) constraints, finding a feasible solution to which corresponds a partitioning at a speedup bound of 1 + 1 k . Approximation schemes have been defined for computing the value of dbf to any desired degree of accuracy (see, e.g. [START_REF] Albers | An event stream driven approximation for the analysis of real-time systems[END_REF]). Equation 6below gives such an approximation scheme; for any fixed positive integer value of k, dbf (k) (τ i , t) defines an approximation of dbf(τ i , t) that is exact for the first k "steps" of dbf(τ i , t), and an upper bound for larger values of t:

dbf (k) (τ i , t) = c i,j × t+pi-di pi if t ≤ (k -1)p i + d i c i + (t -d i )u i otherwise (6)
The following lemma (see, e.g., [START_REF] Baruah | Partitioned scheduling of sporadic task systems: an ILP-based approach[END_REF]) provides a quantitative bound on the degree by which dbf (k) may deviate from dbf:

Lemma 2 For all t ≥ 0 dbf(τ i , t) ≤ dbf (k) (τ i , t) < 1 + 1 k dbf(τ i , t).
That is, dbf (k) (τ i , t) provides an upper bound on dbf (k) (τ i , t) that is no more than a fraction 1/k greater than the actual value of dbf (k) (τ i , t).

As previously stated, it is known that a collection of sporadic tasks can be scheduled to always meet all deadlines upon a preemptive uniprocessor by EDF if and only if for all t ≥ 0, the sum of the dbf's of all the tasks in the collection over an interval of duration t does not exceed t. For schedulability, it is clearly necessary that the utilizations of all the tasks in the collection sum to no more than 1. Since dbf (k) is an upper bound on dbf, it follows that a sufficient uniprocessor EDFschedulability test for a collection of tasks is that the sum of the dbf (k) functions of all the tasks in the collection over an interval of duration t not exceed t. Albers and Slomka showed [START_REF] Albers | An event stream driven approximation for the analysis of real-time systems[END_REF]Lemma 4] that it suffices to validate this fact only for those values of t at which one or more of the dbf (k) functions has a step discontinuity. More precisely, let

S i,k = {t : t = d i + hp i , h = 0, 1, . . . , k} and S k = all i S i,k
It suffices to test that the sum of the dbf (k) functions of all the tasks in the collection over an interval of duration t not exceed t, only for values of t ∈ S k .

We can use this result to define a revised ILP formulation for modeling the partitioned scheduling of sporadic task systems upon heterogeneous multiprocessors. The first part of this revised ILP is identical to the one constructed in Section 3:

-As in Section 3, we will have a zero-one integer variable x i,j , denoting whether τ i is to be assigned to processor π j , for each i ∈ N, j ∈ M . -Again as in Section 3, the following n constraints specify that each task gets assigned to exactly one processor:

∀i ∈ N j∈M x i,j = 1 (7)
-The following m constraints bound the total utilization on each processor:

∀j ∈ M i∈N x i,j u i,j ≤ β (8)
The final set of constraints replace the Inequalities 5 of the ILP in Section 3 with constraints based upon the dbf (k) approximation, that express the requirement that the sum of the dbf (k) functions of all tasks assigned each processor over an interval duration not exceed the duration. As we had stated above, this condition only needs to be validated for interval durations in S k ; this motivates the following set of constraints:

∀t ∈ S k , ∀j ∈ M i∈N x i,j × dbf (k) j (τ i , t) ≤ β • t (9)
(where for each j ∈ M , dbf

j (τ i , t) denotes the function dbf (k) (τ i , t) when the WCET of τ i is set equal to c i,j .) Since each τ i contributes at most k distinct points to S k , it follows that |S k | ≤ nk; hence there are at most mnk such constraints.

Note that the Inequalities 9 are constructed for specified values of t; i.e., for each t ∈ S k . For each such specified t, it is straightforward to observe that the inequality is indeed a linear one, since inspection of Expression 6 reveals that for a given value of t for each τ i the expression dbf

(k) j (τ i , t) is a constant.
Theorem 3 A feasible solution to the ILP on the 0/1 variables {x i,j }, i ∈ N, j ∈ M , with Constraints 7-9 (defined above) yields a feasible partitioning of the tasks in τ to a platform in which each processor in π is speeded up by a multiplicative factor of (1 + 1/k)β. In particular, if the inequalities are satisfied with β ≤ (1 + 1/k) -1 , then a feasible solution to the ILP yields a feasible partitioning on the given platform.

Proof It is evident from the result of Albers and Slomka that satisfying Constraints 7-9 is sufficient for feasibility upon a speed-β(1 + 1/k) platform. Additionally we conclude from the lower bound in Lemma 2 that failure to satisfy these conditions implies infeasibility upon a speed-β platform.

The ILP model consisting of Constraints ( 7)-( 9) will be referred to as Model 2 in the remainder of the paper. The quality of the solution that is obtained by solving Model 2 depends on the value of k: the larger this value, the better is the quality (i.e., the lower the speedup factor) of the obtained solution. However we observe that the number of constraints increases with k. It follows that large values of k lead to an ILP that is not solvable with state of the art packages.

An ILP that is amenable to polynomial time approximation

In the previous sections, we discussed assignment algorithms based on solving reasonably-sized ILPs. However, in some scenarios the solution of an ILP may be a computational bottleneck -solving an ILP is, after all, a prototypical NP-hard problem. Therefore, the design of efficient (polynomial time) assignment and schedulability algorithms retains interest.

A standard technique that has been developed for efficiently obtaining an approximate solution to an ILP is to first consider the linear program (LP) obtained by "relaxing" (i.e., ignoring) the integrality requirement, solve this LP (this can be done in polynomial time), and then "rounding" the solution so obtained to obtained an integral solution as desired. The main challenge in designing the rounding procedure is to ensure that such rounding does not degrade the feasibility or the quality of the solution (i.e., the value of the objective function that was optimized) too much.

Recently, a new approach to rounding, known as iterative rounding [START_REF] Lau | Iterative Methods in Combinatorial Optimization[END_REF], has been shown to provide improved rounding guarantees. Analogously to prior rounding approaches, the first step requires that an LP relaxation be solved and a non-integer solution (say, X) be obtained. However, instead of rounding all non integral values of X at the same time, only one variable is rounded; assume, for example, that the value of variable x 1 is set to x1 . Then the method iterates and solves a new LP-relaxation that is obtained from the original LP by fixing the value of x 1 to x1 . In this way, a new solution X is obtained; as in the previous case, the method now rounds one variable of X ; the method is iterated until all variables satisfy the given integrality constraints.

In this section, we seek to construct an ILP formulation of the problem of partitioning sporadic tasks upon heterogeneous multiprocessors that is more amenable to iterative rounding than the ILP formulations we have seen above. It will turn out (Theorem 4 below) that this ILP has the same number of variables and constraints, but a poorer (larger) speedup factor than the one described in Section 3.3 (Theorem 2). Hence from the perspective of just developing an ILP, this is not a particularly useful result. However, we will see that this ILP can in fact be rounded iteratively in a manner that we were unable to pull off with the earlier ILP formulations, resulting in a polynomial-time algorithm for partitioning sporadic tasks upon heterogeneous multiprocessors that has speedup bound of ≈ 7.83, thereby improving the ≈ 12.9 speedup bound of Marchetti-Spaccamela et al. [START_REF] Marchetti-Spaccamela | Assigning sporadic tasks to unrelated machines[END_REF].

LP-rounding based approach

As before, we use variables x ij for each pair (i, j) ∈ N ×M , modeling the assignment of τ i to π j . Apart from the usual assignment constraints, the first constraints we consider are the utilization bounds on the tasks assigned to the same processor. That is, we require that

i∈N u ij x ij ≤ 1 ∀j ∈ M. (10)
Now, let ρ denote any constant, ρ > 1. We define the function r(x) := ρ log ρ x , and the set D ρ := {ρ 0 , ρ 1 , . . . , r(d max )}. We want to express the requirement that for all tasks assigned to the same processor with deadline at most ρ k , the sum of their WCETs is at most ρ k . Note that this is the set of tasks with r(d i ) ≤ ρ k . For technical reasons that will become apparent later (in Lemma 3), we adopt the slightly weaker constraint

i∈N : r(di)≤d c ij 1 - d i p i x ij ≤ d ∀d ∈ D ρ , ∀j ∈ M. (11)
We call these constraints [START_REF] Chwa | Optimal real-time scheduling on two-type heterogeneous multicore platforms[END_REF] the relaxed dbf constraints.

It is clear that these constraints have to be fulfilled by any feasible task assignment. (In particular, if

d i ≤ ρ k and x ij = 1, then dbf f,j (ρ k ) ≥ c ij > c ij (1 -d i /p i ),
where f is the assignment represented by x). We therefore arrive at the following ILP, denoted poly-ILP.

j∈M x ij = 1 ∀i ∈ N (12a) i∈N u ij x ij ≤ 1 ∀j ∈ M (12b) i∈N : r(di)≤d c ij 1 - d i p i x ij ≤ d ∀d ∈ D ρ , ∀j ∈ M (12c) x ij ∈ {0, 1} ∀i ∈ N, ∀j ∈ M s.t. c ij ≤ d i . (12d)
If poly-ILP is infeasible, then there can be no feasible task assignment. Now assume that it is feasible and consider its relaxation, which is obtained by replacing each constraint (12d) by

x ij ≥ 0 ∀i ∈ N, ∀j ∈ M s.t. c ij ≤ d i . ( 13
)
Since it is an LP and not an ILP, the relaxation can be solved in polynomial time. Let x * denote its solution.

For each j ∈ M and deadline d ∈ D ρ , we compute the value

b j,d := i∈N : r(di)=d c ij (1 -d i /p i )x * ij .
Note that, by (12c),

d ≤d b j,d ≤ d ∀d ∈ D ρ , ∀j ∈ M. ( 14
)
Based on these computed values, we define a variation of poly-ILP, denoted by sparse-ILP in the sequel. We obtain the latter by replacing the constraints (12c) with the following set of constraints:

i∈N : r(di)=d c ij 1 - d i p i x ij ≤ b j,d ∀d ∈ D ρ , ∀j ∈ M.
By dividing both sides of the inequality by d, these constraints can also be written as

i∈N : r(di)=d cijd 1 - d i p i x ij ≤ bj,d ∀d ∈ D ρ , ∀j ∈ M. ( 15
)
where bj,d = b j,d /d ≤ 1 (since x * is feasible for the relaxation of poly-ILP) and

cijd = c ij /d ≤ 1 (since if r(d i ) = d then c ij ≤ d i ≤ d).
Let A be the set of vectors x satisfying (12a). We can now express sparse-ILP in matrix notation as

{x ∈ A ∩ {0, 1} N ×M : Ax ≤ b, x ≥ 0},
where A and b are, respectively, the matrix of coefficients and the vector of right hand sides of the constraints (12b) and [START_REF] Kamath | Unrelated parallel machine schedulingperspectives and progress[END_REF]. Note that all entries of A and b take values between 0 and 1. By construction, if x * is a feasible solution for the LP relaxation of poly-ILP it is also feasible for the LP relaxation of sparse-ILP, and if the LP relaxation of sparse-ILP is infeasible, then no feasible task assignment exists. Our goal will be to round x * to an integral vector which approximately satisfies the constraints of sparse-ILP.

A reason for preferring sparse-ILP to poly-ILP is that the former is an ILP formulation in which the constraint matrix is sparse: each variable appears in only a small number of constraints. This sparsity gives the potential to derive efficient rounding schemes which result in integral solutions, violating the relaxed dbfconstraints only by constant factors. We present such a rounding scheme below; to this end, the following theorem shows that-even when violated up to constant factors-the relaxed dbf constraints (15), together with the utilization constraints (12b), are approximately sufficient.

Theorem 4 Let β ≥ 1 and let f : N → M be an assignment encoded by a vector x ∈ A ∩ {0, 1} N ×M such that, for each j ∈ M and

d ∈ D ρ , i∈N u ij xij ≤ β ( 16
)
and

i∈N : r(di)=d c ij 1 - d i p i xij ≤ b j,d + (β -1) • d. ( 17
)
Then dbf f,j (s) ≤ (β + ρ + (β -1)ρ 2 /(ρ -1))s for all s ≥ 0. In particular, if β ≤ 2, f is a feasible assignment under a speedup factor of (2 + ρ + ρ 2 /(ρ -1)).

Proof For any s ≥ 0 and j ∈ M , we bound dbf f,j (s) as follows:

dbf f,j (s) = i∈N : di≤s s + p i -d i p i c ij xij ≤ i∈N : di≤s s c ij p i + c ij 1 - d i p i xij ≤ i∈N : r(di)≤r(s) s c ij p j + c ij 1 - d i p i xij ≤ s i∈N c ij p i xij + i∈N : r(di)≤r(s) c ij 1 - d i p i xij = s i∈N u ij xij + log ρ (r(s)) k=0 i∈N : r(di)=ρ k c ij (1 - d i p i ) xij (16)
≤ βs + log ρ (r(s))

k=0 i∈N : r(di)=ρ k c ij (1 - d i p i ) xij (17)
≤ βs

+ log ρ (r(s)) k=0 b j,ρ k + (β -1)ρ k (14)
≤ βs + ρ log ρ (r(s)) + log ρ (r(s))

k=0 (β -1)ρ k = βs + r(s) + (β -1) log ρ (r(s)) k=0 ρ k = βs + r(s) + (β -1) • ρ log ρ (r(s))+1 -1 ρ -1 = βs + r(s) + (β -1) ρ • r(s) -1 ρ -1 ≤ (β + ρ + (β -1) ρ 2 ρ -1 )s.
The last inequality follows from r(s) ≤ ρ • s.

To construct x, we adopt an iterative rounding procedure that is similar to the procedure presented in [START_REF] Karp | Global wire routing in two-dimensional arrays[END_REF][START_REF] Marchetti-Spaccamela | Assigning sporadic tasks to unrelated machines[END_REF]. The idea of the iterative rounding procedure is the following. In each iteration k, we first compute an extreme point solution x k of a linear program LP k , where LP 0 is the relaxation of sparse-ILP, and each LP k is obtained by fixing the value for some variables or removing some constraints of LP k-1 .

Given a feasible fractional solution x k , to define LP k+1 we first fix all variables which are integral in x k , i.e., those variables are not allowed to be changed anymore in the remainder of the procedure. Let v be the number of variables in LP k and let w a , w b and w b be the number of constraints of types (12a), (12b) and [START_REF] Kamath | Unrelated parallel machine schedulingperspectives and progress[END_REF], respectively. Let w = w a + w b + w c . To obtain LP k+1 we either delete one or more variables, in case v > w, or delete a constraint while ensuring that in the final solution that constraint will not be violated too much. Along the way we ensure that the constraints of type (12a) are always satisfied exactly, so that x k ∈ A at all times.

Note that if there is some variable x ij that is fixed at value 1 and removed from the program, then for all j ∈ M \ {j}, x ij will be set to 0 and also be removed from the program. The constraint of type (12a) corresponding to this i is then superfluous and will also be removed.

To derive the bounds ( 16)-( 17), we need to study the coefficient matrix A in more detail. Let γ be the maximum, over all x ij , of the sum of the values of the coefficients of variable x ij in constraints (12b) and [START_REF] Kamath | Unrelated parallel machine schedulingperspectives and progress[END_REF]. We first derive a bound on γ.

Lemma 3 For any task set τ , γ ≤ 1.

Proof Observe that γ is just the maximum value of u ij + cijd (1 -d i /p i ) across all variables x ij in the program. Recall that for all such pairs (i, j), c ij ≤ d i , i.e., cijd ≤ 1, otherwise the variable x ij is forced to 0 and removed from the LP. We can now bound

u ij + cijd 1 - d i p i ≤ c ij p i + 1 - d i p i ≤ c ij p i + 1 - c ij p i = 1.
The following technical lemma is instrumental to our rounding procedure. It is a specialization to our setting of a more general rounding result for assignment LPs [START_REF] Bonifaci | Algorithms for hierarchical and semi-partitioned machine 0[END_REF]. Lemma 4 Let LP k be the linear program that is solved in iteration k of the rounding procedure, with s variables and r constraints. Let x k be an extreme point solution to this LP. Then either, (i) x k has at least one integer component; or (ii) there is j ∈ M and a corresponding constraint of type (12b) such that i∈N u ij z ij -i∈N u ij x k ij ≤ γ for any integer solution z; or (iii) there are j ∈ M , d ∈ D ρ and a corresponding constraint of type [START_REF] Kamath | Unrelated parallel machine schedulingperspectives and progress[END_REF]

such that i∈N : r(di)=d cijd 1 - d i p i (z ij -x k ij ) ≤ γ
for any integer solution z.

Proof Let A be the coefficient matrix of LP k . If v > w, the null space of A is nontrivial, so let x 0 be a nonzero vector in the null space. Since x k is an extreme-point solution to LP k , it cannot be expressed as the convex combination of two (or more) solutions to LP k . If

x k does not have any integral entry, then we can find a value δ > 0 such that x k + δx 0 and x kδx 0 are both solutions to LP k (since A(x ± δx 0 ) = Ax) and, in particular, x k is a convex combination of these two solutions. Therefore x k must have at least one integral entry.

If v ≤ w, we show that there always exists a constraint of type (12b) or type [START_REF] Kamath | Unrelated parallel machine schedulingperspectives and progress[END_REF] such that max z∈S {(Az) l -(Ax) l } ≤ γ, where γ is the maximum sum of coefficients in a column of constraints (12b) and ( 15) and where S is the integer solution space for all remaining variables, i.e., S = {0, 1} v .

We show the statement by contradiction. Assume that the statement is not true, that is, for each constraint l of type (12b) or ( 15) it holds that there exists a vector z ∈ S such that

(Az) l -(Ax) l > γ. (18)
Note that all variables still present in the linear program correspond to a processor j ∈ M and a task i ∈ N that is not yet assigned fully to one processor, but fractionally to multiple processors. Hence, the constraint of type (12a) corresponding to each τ i is still present in the linear program. It follows that j∈M i∈N : xij ∈(0,1)

x ij = w a , (19)
where w a is the number of constraints of type (12a) remaining in LP k . Define L as the set of constraints of types (12b) and ( 15) present in the current linear program, and let w b and w c be their number, respectively (so w a + w b + w c = w). For any q = (i, j), let L q denote the set of these constraints containing variable x q ; by definition of γ, we have

max q l∈L q a lq ≤ γ. (20)
Then,

γ(w -w a ) = γ(w b + w c ) (18)
< l∈L max z∈S ((Az) l -(Ax) l ) as all a lq ≥0 = l∈L ((A1) l -(Ax) l ) = l∈L q a lq (1 -x q ) = q l∈L q a lq (1 -x q ) ≤ q γ(1 -x q ) = γv - q γx q = γ(v -r a ). (21)
The second inequality follows from [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard real-time environment[END_REF].

The chain of inequalities implies that γ(ww a ) < γ(ww a ) ⇒ r < s which is a contradiction to being in the case that v ≤ w. Hence we conclude that if v ≤ w, there must be a constraint l of type (12b) or [START_REF] Kamath | Unrelated parallel machine schedulingperspectives and progress[END_REF] for which

max z∈S {(Az) l -(Ax) l } ≤ γ.
Lemma 4 is used to guide the rounding process. If Case (i) applies, the variables that have an integer value are fixed at that value and removed from the LP. If a variable x ij is fixed at value 1, then for all j ∈ M \{j}, the variables x ij are fixed at value 0 and the constraint of type (12a) corresponding to i is removed. If we are in Case (ii) or (iii), the constraint for which the claim holds can be found in polynomial time by checking, for each constraint l ∈ L of type (12b) or [START_REF] Kamath | Unrelated parallel machine schedulingperspectives and progress[END_REF], whether q a lq (1 -x q ) ≤ γ. This is sufficient since all a lq ≥ 0 and the maximum value any variable x q can take is 1. If such a constraint is of type (12b) (Case (ii)) or [START_REF] Kamath | Unrelated parallel machine schedulingperspectives and progress[END_REF] (Case (iii)), the final task assignment will satisfy [START_REF] Karp | Reducibility among combinatorial problems[END_REF] or [START_REF] Karp | Global wire routing in two-dimensional arrays[END_REF] for that constraint, respectively, even if the constraint is dropped ; thus, we drop the constraint, obtaining the next (smaller) LP.

After either all constraints have been removed or the values of all variables have been fixed at an integer value, we obtain an integral vector x which satisfies

i∈N u ij xij ≤ 1 + γ for each j ∈ M and i∈N : r(di)=d cijd 1 - d i p i x ij ≤ bj,d + γ
for all j ∈ M and all deadlines d ∈ D ρ . Hence, the vector x satisfies constraints ( 16), [START_REF] Karp | Global wire routing in two-dimensional arrays[END_REF] with β := 1 + γ.

We are now in the position to invoke Theorem 4 to obtain our final guarantee.

Theorem 5 There is a polynomial-time partitioning algorithm with a speedup bound of (5 + 2 √ 2) ≈ 7.83 for the problem of assigning constrained-deadline tasks to heterogeneous processors.

Proof All steps required to construct x can be carried out in polynomial time. The assignment induced by x satisfies ( 16)-( 17) with β = 1 + γ ≤ 2 (Lemma 3). Thus, by Theorem 4 with ρ = 1+1/ √ 2, the assignment induced by x is feasible with speedup

2 + ρ + ρ 2 ρ -1 = 5 + 2 √ 2 ≈ 7.83.

Practical variant of the LP-rounding model

The above LP-rounding approach was designed with the aim of minimizing the worst-case speedup bound in Theorem 5. However, we observed that certain steps required for theoretical soundness -namely, sparsification of the constraints-can be avoided in practice. Thus, in this subsection we briefly discuss the practical variant of the LP-rounding model that is the one we adopted in the schedulability experiments, and that we will refer to as Model 3.

The starting point is the following minimization version of equations (12a)-(12d):

min β (22a) j∈M x ij = 1 ∀i ∈ N (22b) i∈N u ij x ij ≤ β ∀j ∈ M (22c) i∈N : r(di)≤d c ij d 1 - d i p i x ij ≤ β ∀d ∈ D ρ , ∀j ∈ M (22d) x ij ∈ {0, 1} ∀i ∈ N, ∀j ∈ M s.t. c ij ≤ d i , (22e)
where β is the speedup parameter to be minimized. We apply to this model the iterative rounding procedure as described in the previous subsection, except that when we need to drop a constraint, we drop the constraint l for which the potential violation (i.e., q a lq (1 -x q ) where A is the coefficient matrix of Constraints (22c)-(22d)) is the smallest. We also keep track of the largest potential violation value during the whole process (we call it γ, in analogy with the previous section). To determine whether the instance is schedulable, at the end of the process it is sufficient to compare the final value of β +γ against (1+ρ) -1 , due to the following Theorem.

c ij 1 - d i p i xij ≤ (β + γ)d. (24)
Then dbf f,j (s) ≤ (1 + ρ)(β + γ)s for all s ≥ 0. In particular, if β + γ ≤ (1 + ρ) -1 , then f is a feasible assignment on unit speed processors.

Proof The proof is similar to that of Theorem 4. For any s ≥ 0 and j ∈ M , we bound dbf f,j (s) as follows:

dbf f,j (s) = i∈N : di≤s s + p i -d i p i c ij xij ≤ i∈N : di≤s s c ij p i + c ij 1 - d i p i xij ≤ i∈N : r(di)≤r(s) s c ij p j + c ij 1 - d i p i xij ≤ s i∈N c ij p i xij + i∈N : r(di)≤r(s) c ij 1 - d i p i xij = s i∈N u ij xij + i∈N : r(di)≤r(s) c ij 1 - d i p i xij ≤ s(β + γ) + (β + γ)r(s) ≤ s(β + γ) + ρs(β + γ) = (1 + ρ)(β + γ)s.
The last inequality follows from r(s) ≤ ρ • s.

Schedulability experiments

In the sections above, we saw how the problem of partitioned scheduling of sporadic task systems upon unrelated multiprocessors could be modeled by ILPs. Our motivation for doing so is that the optimization community has devoted immense effort to coming up with extremely efficient (although still exponential-time, since solving ILPs is NP-hard) algorithms for solving ILPs, and highly-optimized libraries implementing these efficient algorithms are widely available. This is particularly true for ILPs like the ones we have constructed above, in which each variable is further constrained to take in only the values zero or one. In this section, we validate the performance of our ILP-based schedulability tests against synthetic workloads in terms of the percentage of schedulable task sets.

Our results show that Model 2 (the ILP model discussed in Section 4) provides solutions having better quality, while Model 1 (the ILP model discussed in Section 3) generates solutions of lower quality but in less time. Model 3, the polynomial time LP-rounding approach discussed in Section 5 has a solution quality comparable to that of the Model 1, but runs even faster.

Generation of the task sets and solutions

We developed a parametric framework with several parameters (m, κ, p, Ū , α -they are detailed below) to randomly generate our workloads; this framework is general enough to support our entire range of experiments.

We consider m-processor platforms and n = κm tasks, with κ ≥ 1 an integer-valued parameter. We randomly generate an affinity mask R i,j for each i ∈ N and j ∈ M : R i,j ← 1 with probability p and 0 with probability (1 -p). These affinity masks help define the C i,j values: C i,j has a value < ∞ if and only if R i,j = 1. If the generated mask does not allow a particular task to be processed on any processor, we then allow that task to be processed upon a randomly chosen processor.

We then generate utilization values for every allowed pair (i, j) for which R i,j = 1. Tasks are grouped into m groups of size κ each: tasks τ 1+(k-1)κ to τ kκ form the kth group. For each group we use the UUniSort algorithm [START_REF] Bini | Measuring the performance of schedulability tests[END_REF] to generate randomly distributed utilizations with total value Ū for the allowed task-processor pairs in the group. Note that since there are m groups, each of total utilization Ū , the value Ū represents the average load that a processor can expect if tasks are randomly assigned. Note that Ū ≤ 1 is not a necessary condition for schedulability (indeed some of our ILP formulations are able to schedule task sets with Ū > 1).

We generate the periods by setting p i = 2 ∆i , with each ∆ i uniformly distributed in the range 3...10. The worst-case execution times are computed directly from the periods and utilizations. Finally, the relative deadline of each task is sampled uniformly in the range

[(1 -α) • (max j c ij ) + αp i , p i ], where α ∈ [0, 1].
In the experiments we consider the three models: Model 1 (that is, the ILP model discussed in Section 3.2), Model 2 (that is, the ILP model discussed in Section 4) with k = 3, and Model 3 (that is, the iterative LP-rounding approach discussed in Section 5). All optimization models have been solved by using a branch-and-cut approach implemented in the mathematical programming solver Gurobi 6.50 [START_REF]Gurobi Optimizer[END_REF] on a PC with Intel i7-4770 CPU at 3.4 GHz and 16Gb RAM. Collectively, the experiments consist of 5200 instances of the task partitioning problem.

6.2 Discussion of the results

Schedulability

Experiment 1: Variation of Ū (Figure 1) In the first type of experiment, we use the average load Ū as the independent variable. We vary Ū from 0.2 to 1.5. We fix m = 10, κ = 10, p = 0.5, α = 0.2. We generate 30 task sets for each value of Ū . When the percentage of schedulable task sets was strictly between 0 and 1, we generated 20 additional task sets to achieve a higher precision. Figure 1 shows the dependency on Ū of the percentage of task sets that are guaranteed to be schedulable by Corollary 2 (for Model 1), Theorem 3 (for Model 2) and Theorem 6 (for Model 3). Indeed, as could be expected, schedulability decreases when the load factor Ū increases. It is interesting to note that for Model 2, more than 90% of the generated task sets are schedulable as long as Ū ≤ 1.

Experiment 2: Variation of m (Figure 2) In the second type of experiment, we look at the impact of the number of processors on schedulability, so we vary m from 2 to 12. We fix κ = m, so n = m 2 . In this set of experiments, we fix Ū = 1, p = 0.8, α = 0.2. For each experiment, we generated 30 task sets for each value of m. In all cases, schedulability increases with the number of processors even though Ū is fixed. This is due to the fact that, for a fixed total utilization value of each group, a higher number of tasks in a group implies lower utilization values for the individual tasks, and therefore a more efficient partitioning. The best quality is again obtained by Model 2, while Model 1 and Model 3 achieve a similar solution quality.

Experiment 3: Variation of p (Figure 3) In the third type of experiment, we control the sparsity of the processor affinity matrix R, by varying p from 0.2 to 0.9.

We fix m = 10, κ = 10, Ū = 1, α = 0.2. We generate 30 task sets for each value of p. We find out that the sparsity has a high impact on schedulability: there are clear schedulability thresholds around p = 0.7 (for Models 1 and 3) and p = 0.45 (for Model 2). This is not entirely unexpected, as when the affinity matrix is sparser, it may happen that several tasks of large combined utilization can only be assigned to a small set of processors.

Experiment 4: Variation of the number of processor types (Figure 4) In the fourth type of experiment, we control the effect of similarity among processors. We group the processors into types; we ensure that c i,τ = c i ,τ whenever i and i are of the same type. We then vary the number of types of processors. We fix m = 8 and for simplicity we consider equal-sized groups, so we vary the number of types in the set {1, 2, 4, 8}. We also fix κ = 8, Ū = 0.6, p = 0.5, α = 0.2, and we generate 30 task sets for each value of the independent variable. In this case, it seems more significant to plot the average speedup that is required to ensure schedulability (thus, lower values are preferable, and values below 1 ensure schedulability without speedup). As expected, similarity is helpful, i.e., an increase in the number of types is associated with an increase in the speedup required, particularly for Models 1 and 3.

Running times

In this subsection we discuss the solution time required to solve the instances for the experiments described in the previous section. These times are reported in the Tables 2, 3, 4, and 5. As a general observation, running times for solving Models 1 and 3 are considerably shorter than those needed for Model 2, but as we have seen the latter one is clearly more effective in terms of schedulability. We also note that, for Models 1 and 2, each set of instances generated with the same parameters may have rather variable solution times: the maximum time is often more than 10 times the average time. Model 3, on the other hand, is extremely fast and its maximum times are quite close to the averages. For Experiment 1, where we vary the average load Ū (Table 2), we observe a slight decrease of the running times as the instances become less schedulable. This is likely due to the solver being able to rule out quickly those instances that happen to be markedly overloaded.

For Experiment 2, where we vary the number m of processors (Table 3), the running times increase rapidly with m, but this is in large part due to the fact that in this scenario we generate n = m 2 tasks per instance, so even the size of the input becomes of the order of m 3 .

Theorem 6

 6 Let β, γ > 0 and let f : N → M be an assignment encoded by a vector x ∈ A ∩ {0, 1} N ×M such that, for each j ∈ M and d ∈ D ρ , i∈N u ij xij ≤ β + γ (23) and i∈N : r(di)≤d

Fig. 1 Fig. 2

 12 Fig.1Experimental results 1: percentage of task sets found to be schedulable by the three models. Variation of Ū (m = 10, κ = 10, p = 0.5, α = 0.2).

Fig. 3 Fig. 4

 34 Fig.3Experimental results 3: percentage of task sets found to be schedulable by the three models. Variation of p (m = 10, κ = 10, Ū = 1, α = 0.2).

Table 1

 1 Summary of ILP models and results

Although we expect that most of our results will also extend to the arbitrary-deadline sporadic task model, for ease

After discounting this fact, the running times are still positively correlated with m, in particular for Model 2.

For Experiment 3, where we vary the affinity probability p (Table 4), we observe a sharp transition from unschedulable to schedulable instances, and running times progressively increase. Again, this is likely due to the solver being able to rule out quickly those instances where the affinity relation is too sparse to allow schedulability.

For Experiment 4, where we vary the number of types of processors (Table 5), all running time are limited to a few seconds at most and we do not find any significant correlation between the running time and the number of types.

Summary and conclusions

In this work, we proposed a partitioning approach for constrained-deadline tasks on heterogenous (unrelated) processors. The approach is based on integer linear programming formulations and allows the derivation of guaranteed speedup bounds and consequently, sufficient schedulability tests.

Experiments using randomly generated task workloads clearly show that the approach based on the solution of Model 2 is viable in terms of computation time, and quite effective in its scheduling capacity, especially when the task-processor affinity relation is dense. On the other hand, the approach based on the solution of Model 1, despite its reduced scheduling capacity, may be of use in the case of instances so large that they cannot be solved with the previous approach in reasonable times. The LP-rounding variant, Model 3, achieves a comparable quality to Model 1 and appears to be extremely fast in practice.

Work supported by NSF grants CNS 1115284, CNS 1218693, CNS 1409175, and CPS 1446631, AFOSR grant FA9550-14-1-0161, ARO grant W911NF-14-1-0499, and a grant from General Motors Corp. This work has also been partially supported by the research project "Designing Human-Agent Collectives for Sustainable Future Societies" (C26A15TXCF) of Sapienza University of Rome, by MIUR PRIN grant 2012JXB3YF_003, and by the National Group of Computing Science (GNCS-INDAM).


domain: spi phys

text:
Xavier Dupuis

Valentine Wasik

Alexandre Alakian

Dominique Dubucq

MULTI-BAND SUPERVISED CLASSIFICATION FOR POLARIMETRIC SAR

Keywords: SAR, classification, multi-band

This work addresses the potential of multi-band polarimetric SAR imaging for terrains and vegetation classification. A classic supervised Wishart classifier is adapted to polarimetric multi-band datasets, and is applied on the X-, L-and UHF-band acquisitions done during the NAOMI campaign (ONERA-Total) in Gabon (Africa) in 2015. The contributions of the different frequencies are shown and discussed. It is shown that the use of the multi-band dataset improves significantly the classification result.

INTRODUCTION

The NAOMI project (New Advanced Observation Method Integration) is a Total and ONERA partnership initiated in 2014. Among the objectives of this project, one is to analyze the capacities of remote sensing techniques to characterize surfaces (biophysical and geological) and to detect hydrocarbons (directly or indirectly through environmental descriptors such as vegetation health and state). In this context, the NAOMI campaign that took place in Gabon in 2015 allowed ONERA and Total to construct a particularly diversified dataset composed by polarimetric SAR and VNIR-SWIR hyperspectral acquisitions along with in-situ measurements over the Mandji peninsula. This site presents an alternation between dense tropical forests, sparse and short vegetation, backwaters and man-made structures, with oil deposits that have been exploited for over 50 years in this area. On such equatorial area, optical acquisitions are difficult because of clouds and heavy rains season. Radar is a good alternative to overcome the weather problem. Moreover the swath of a SAR sensor is wider than optical one, which makes SAR imaging more suitable to cover large areas in such a place. Landscape classification from SAR imaging seems to be an interesting technique to characterize terrains and vegetation for the purpose of the NAOMI project. Moreover, the Mandji dataset gives us the opportunity to explore the contribution of several bands in the results of the classification products, since acquisitions have been done in X-, L-and UHF-band. In this work, we propose to adapt and apply a classic supervised Wishart classifier [START_REF] Baqué | SETHI / RAMSES-NG New performances of the flexible multi-spectral airborne remote sensing research platform[END_REF] to our multi-band SAR dataset. After a brief presentation of the data, the classifier methodology is explained. The reference regions from in-situ measurements are then detailed, and the classification results are finally discussed.

DATASET

Data were acquired with SETHI [START_REF] Baqué | SETHI / RAMSES-NG New performances of the flexible multi-spectral airborne remote sensing research platform[END_REF], an ONERA sensor, in Gabon over Mandji area in 2015, July 8 th . Two flights were performed, one to acquire simultaneously X and L bands and a second to acquire simultaneously UHF and L bands. X and UHF bands were acquired in full polarimetric mode. Lband was acquired in dual-polarimetric mode. The SAR dataset is coregistered and projected into the Xband native slc geometry for classification purpose. Flights dedicated to optical acquisitions were operated the same day. The characteristics of the hyperspectral cameras are presented in SAR and hyperspectral (VNIR, SWIR) data are groundcoregistered [START_REF] Plyer | A New Coregistration Algorithm for Recent Applications on Urban SAR Images[END_REF] in a common WGS84 referential.

MULTIBANDS CLASSIFIER

The Wishart classifier is an interesting and largely used technique for the classification of the landscapes with polarimetric SAR images [START_REF] Lee | Classification of multi-look polarimetric SAR imagery based on complex Wishart distribution[END_REF]. This iterative method focuses on the statistical properties of the response. For fullypolarimetric SAR data, the complex response for each pixel can be written in the Pauli vector 𝒌 = 1

√2

[𝑆 𝐻ℎ + 𝑆 𝑉𝑣 , 𝑆 𝐻ℎ -𝑆 𝑉𝑣 , 2𝑆 𝐻𝑣 ] 𝑇 , where 𝑆 𝑋𝑦 is the complex backscattered signal for an emitted polarization 𝑋 and a direction of reception 𝑦 (𝑋 = 𝐻, 𝑉 ; 𝑦 = ℎ, 𝑣). In case of fully-developed speckle, 𝒌 is assumed to be a random vector following a complex circular multi-variate Gaussian law with zero-mean for which the statistical properties are entirely described by its covariance matrix 𝐓 = 〈𝒌𝒌 † 〉. For a set of 𝑁 measurements {𝒌 1 , … , 𝒌 𝑁 }, the maximum-likelihood estimator for 𝐓 (ie 𝐓 ̂= 1 𝑁 ∑ 𝑘 𝑛 𝑘 𝑛 † 𝑁 𝑛=1

) follows a complex Wishart law [START_REF] Baqué | SETHI / RAMSES-NG New performances of the flexible multi-spectral airborne remote sensing research platform[END_REF]. The principle of the Wishart classifier is to use this property to construct a measure for classification. For each pixel of the image and its associated empirical covariance matrix 𝐓 ̂, the following measure is calculated for each class 𝑤 𝑖 [START_REF] Baqué | SETHI / RAMSES-NG New performances of the flexible multi-spectral airborne remote sensing research platform[END_REF]: In the Mandji data presented, the acquisitions have been done with three different frequency bands: X, L and UHF. The enhanced behaviors of the scene will be different from one frequency band to the other since the interaction between the media and the illuminating wave is very frequency dependent. The contribution of this frequency diversity in the classification process is an interesting point to focus on. For this purpose, we have modified the Wishart classifier so that multi-band polarimetric SAR data could be processed. For each pixel, the idea is to estimate one covariance matrix for each considered band (𝐓 ̂𝑋, 𝐓 ̂𝐿, 𝐓 ̂𝑈𝐻𝐹 ) and to calculate independently the associated Wishart measures for a common ensemble of classes. The pixel is attributed to the class 𝑤 𝑗 that minimizes the sum of the three measures: 𝑤 𝑗 = min 𝑖 (𝑑 𝑤 𝑋 (𝐓 ̂𝑋, 𝑤 𝑖 ) + 𝑑 𝑤 𝐿 (𝐓 ̂𝐿, 𝑤 𝑖 ) + 𝑑 𝑤 𝑈𝐻𝐹 (𝐓 ̂𝑈𝐻𝐹 , 𝑤 𝑖 )).

𝑑 𝑊 (𝐓 ̂, 𝑤 𝑖 ) = ln|𝐕 𝑖 | + 𝑡𝑟(𝐕 𝑖 -1 𝐓 ̂),
The algorithm is then iterated as previously.

One difficulty comes from the fact that the L-band acquisitions are dual-polarized (emission in 𝐻 only). To overcome this limitation, the coefficient 𝑆 𝑉𝑣 𝐿 in the Pauli vector for the L-band is replaced by the ratio 𝑆 𝐻𝑣 𝐿 𝑆 𝐻ℎ 𝐿 ⁄ .

GROUND TRUTH

Ground truth is based on in-situ photography and ASD reflectance ground measurements. These measurements are used to identify the classes, but they do not cover a sufficient area to get a reference ground truth map for the classification (see Figure 5 (a)). To extend the ground truth map, a classification on both VNIR and SWIR data taking into account the shape similarities of spectral reflectance (spectral angle), the reflectance amplitude similarities and the specific absorption is performed. The spectral angle criterion is the most sensitive parameter, a strict constraint is 1° threshold and less strict is 3°. The results are presented in Figure 5.

SAR CLASSIFICATION RESULTS

Prior to the classification, a bilateral filtering is processed on each SAR images to attenuate the speckle noise [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF]. The classification algorithm studied in this paper is supervised, so we first define the reference classes based on the extended ground truth map (Figure 5 (c)). The "old backwater" class was not considered for SAR classification since it was not sufficiently different from other classes. A class of "shadow" (colored in black in the following) has been added due to the particularity of the SAR acquisition geometry as well as a class for "man-made structures" in red A classification is performed for each band independently (X, L, and UHF) to analyze the impact of each radar band on the results. Then multi-band classifications are processed for "X plus L" and "X plus L plus UHF" combinations. Note that all the results presented below are based on the same class definition. Figure 6 presents the classification result for X-band dataset.

It appears noisy because X-band is very sensitive to roughness of small objects, accordingly to its wavelength, and cannot distinguish vegetation types as trees, dense and short vegetation. At the opposite, the water (sea and backwater in blue), and the sand (in grey) are well detected. Figure 7 presents the L-band classification result. Compared to X-band, there are confusions for water (sea and backwater) but vegetation is less noisy except for the trees (dark green). Figure 8 presents the UHF classification result. Except for the strong echoes from buildings and trees (red class) and for the sand on the beach, the classification is very noisy due to non-sensitivity of such a wavelength to short and medium vegetation. As it can be seen on Figure 10, combining the three bands, X, L and UHF, improves significantly the classification of tall vegetation (in the white ellipse) and building (in red). This result takes advantage of the different bands: X-band for water and bare surfaces, L-band for short and medium vegetation and UHF-band for tall vegetation and buildings.

Figure 10 : X+L+UHF band classification.

CONCLUSION

For the Mandji dataset acquired by ONERA during the NAOMI campaign, multi-band SAR classification improves significantly the results for vegetated areas with different height and state of vegetation in comparison with using only mono-band data. It takes advantages form the different wavelength.

To complete this first qualitative analysis, a quantitative study with confusion matrices is going to be processed.

Figure 1 to

 1 Figure 1 to Figure 3 present the SAR images of the area at X-band, L-band and UHF-band respectively.

Figure 1 :

 1 Figure 1: X-band data, Hh polarization is coded on the Red channel, Vh on the Green channel and Vv on the Blue one.

Figure 2 :

 2 Figure 2: L-band data, Hh polarization is coded on the Red channel, Hv on the Green channel and Hh on the Blue one.

Figure 3 :

 3 Figure 3: UHF-band data, Hh polarization is coded on the Red channel, Vh on the Green channel and Vv on the Blue one.

Figure 4 :

 4 Figure 4 : Optical image of the Mandji site.

  where 𝑉 𝑖 is the covariance matrix of the pixels in the class (i.e. class center), |. | is the determinant, and 𝑡𝑟(. ) is the trace operator. The pixel is assigned to the class 𝑤 𝑗 for which the Wishart measure is the smallest. This algorithm is iterative. In our study the initialization is based on the supervised definition of the classes and for each step the class centers are updated. The algorithm is stopped after 5 iterations.

Figure 5 :

 5 Figure 5 : Ground truth map : (a) in-situ ground truth (from ASD measurement), (b) classes definition, (c) computed map with 1° spectral angle threshold and (d) computed map with 3° spectral angle threshold.

Figure 6 :

 6 Figure 6 : X-band classification.

Figure 7 :

 7 Figure 7 : L-band classification.

Figure 8 :

 8 Figure 8 : UHF-band classification.A first result of multi-band combination is presented Figure9. Qualitatively, the "Xplus L-"band classification seems to provide a good result as it takes advantages of both bands. The sea and backwaters appears qualitatively well classified; the different classes of vegetation seem close to the extended ground truth map of Figure5 (d). But the buildings (in red) and the trees (in the white ellipse) are not very well delimited.

Figure 9 :

 9 Figure 9 : X+L band classification.

Table 1

 1

	summarizes the characteristics of the data.
	SAR data	UHF/VHF band	L band	X band
	Bandwidth	200 MHz	150 MHz	300 MHz
	Range resolution	0.66m	0.88m	0.44m
	Azimuth resolution	1.33m	0.88m	0.44m
	Windowing rectangular	rectangular	rectangular
	Polarization	full	Dual Hh-Hv	full

Table 1 :

 1 SAR dataset main characteristics.

Table 2 .

 2

		Hyspex VNIR FOV expander	Hyspex SWIR
	Imaging	Pushbroom	Pushbroom
	Spatial pixels	1600	320
	Field of view	34°	14°
	Spectral range	0.4 -1 µm	1 -2.5 µm
	Spectral bands	160	256
	Spectral sampling distance	3.7 nm	6 nm
	Ground sampling distance	1.3m	2.6m

Table 2 :

 2 Hyperspectral cameras characteristics.


text:
M-J Cantalloube

IMAGING OF SNOW/ICE SUBSURFACE FEATURES FROM AIRBORNE SAR AT UHF, L AND X BAND. THE ONERA SAR CAMPAIGN IN SOUTH GREENLAND

Keywords: Ground penetration radar, snow, ice, stereo radar imaging

An acquisition campaign of the ONERA airborne synthetic aperture radar system RAMSES have been set up on the Greenland ice-shelf to assess the detectability of buried objects and subsurface feature through dry snow/ice cover. Three band were investigated UHF band, L band and X band, forecasting a better penetration at lower band, but an higher radar cross section for a given object/feature at higher band, thus a trade off for the optimal frequency band to use that was not clear prior to experiment. Three sites were imaged, one in the accumulation zone North of the polar circle (East of Kangerlussuaq airport), one in the accumulation zone in a warmer southern area (East of Narsassuaq airport) and one in the intermediate zone between accumulation and ablation zones (plus a calibration zone and a front glacier close to Kangerlussuaq). Processing of the data revealed that both UHF and X band are indeed providing complementing results (deeper bigger feature versus shallow smaller features), but intermediate L band combining both drawbacks of higher attenuation and lower features RCS did not seem to yield any further information.

INTRODUCTION

Satellite SAR imaging of icecaps and ice-shelves revealed an high radar return at C or X band, though the surface appears visually very smooth and homogeneous. SAR imaging of icy celestial bodies by space probes also revealed unexpectedly high surface reflection [START_REF] Campbell | Galilean satellites: 1976 radar results[END_REF]. The supposed reason is that dry snow/ice is not very conductive, and large radio wave penetration allows buildup of radar cross section from several ice/snow layers. This was corroborated by comparison between the ice shelf altitude derived from LIDAR measurements (such as ICEsat arctic terrain mapping) and the altitude derived from interferometric SAR (such as RadarSAT at C-band or TerraSAR at Thanks to DGAC (French Air Directorate) for funding, to the Danish GEUS and Arctic Command for the valuable ground truth, target setup & operational support. ICEsat DEM are provided by the Polar Geospatial Center under NSF OPP awards 1043681, 1559691 & 1542736 X-band) showing that the radar "surface" measured is some (tens of) meters below the optical surface [START_REF] Dall | Insar elevation bias caused by penetration into uniform volumes[END_REF]. The campaign presented here intends to use the higher spatial resolution and viewing angle flexibility of an airborne SAR sensor to investigate the though snow/ice imaging capability of synthetic aperture radar.

EXPERIMENTAL SETUP

The RAMSES airborne SAR system of ONERA was operated either at UHF+L bands in full polarisation mode with a range resolution of 63 and 100 cm resp, or at X-band alone in full polarisation with a range resolution of 20 cm). For calibration, a corner reflector was set close to the Kangerlussuaq airport (Fig. 1). Its position was measured with a few cm accuracy by differential GPS, allowing for fine estimation of the internal radar delays. A similar corner reflector was buried under packed snow at one of the work area (area number 2, the "warm" accumulation zone). In this zone, a Luneburg lens with a partially emerging surface was set up (during the experiment, it ended up buried under 50 cm of fresh snow), as well as a few "artifacts" (hole dug and then filled back with packed snow and an ice lens made of a plastic sled filled with water buried at about 2 m depth).

RESULTS FROM NADIR LOOKING (SOUNDING) SAR IMAGE MODE

Though the antennae were pointed to the side (45 • down wing at UHF & L band and 30 or 50 • at X band), significant side lobe (especially for V polarisation at X band) and the huge specular reflection at nadir allows nadir imaging (Fig. 2).

The sounding image at X band (with a depth resolution in vacuum of 20 cm) shows a strong discontinuity at an apparent depth of about 2 m. This is most probably the transition between the fresh snow fallen during the winter (acquisition is made early April) and the firn (percolated refrozen snow anterior to the previous summer). No further details appear at X band below this discontinuity. Note that the vertical scale inside the snow/ice cover should be downscaled by the refraction index, which is 1.25 to 1.35 in the upper layer and increases up to 1.76 for the solid ice at large depth.

RESULTS FROM SIDE LOOKING (USUAL) SAR IMAGE MODE

Usual (side look focused) SAR images such as fig. 3 yield similar conclusions: on X band image, the edges of shallow crevices are visible, but deeper feature such as refrozen lakes are invisible. The ice cascade on the right edge (darker vertical zone) appears at X band, but the snow cover above ice is very shallow or even absent on this feature. The RCS of snow/ice cover at X band is indeed high: it ranges (for co-polar) approximatively from -3 dBm 2 /m 2 to +1 dBm 2 /m 2 (for down wing angle from 25 • up to 65 • ) which is about 7 dB above that of a typical temperate forest. At L band, the deeper crevices parallel to acquisition path becomes more visible and the buried frozen lakes appear blurred by the diffusion through the snow/ice upper layers. Comparatively, the UHF image is considerably richer in deeper features: the tip of the air cavity in the crevices on the middle left of the image are buried under 40 m of snow/ice (compared with 2 m for the more recent one on the right) and the dense refrozen "floor" is around 65 m deep. The UHF resolution, however, is less than a third of that of Fig. 2. Nadir (sounding) SAR mode images. X band (top), L band (middle), UHF band (bottom). Vertical span is 300 m, length is 10 km the X-band, and RCS of small targets such as out setup is very low at this band. On fig. 3 the white area at the top and bottom edges of L and X band images corresponds to high noise to signal ratio where the aircraft roll motions shift the antenna pattern main lobe away from this edge of the swath. The hue vertical gradient on the UHF image is physical: the cross-polar to co-polar ratio increases with incident angle, which is expected from a "semi infinite random media" model ( [START_REF] Tishkovets | Coherent opposition effects for semi-infinite discrete random medium in the double-scattering approximation[END_REF] for example). The antenna orientation for the L-band was not precisely determined (due to its low interest), hence the antenna pattern compensation is poor for this image [This should be corrected for, in the final article].

OBSERVATION OF BURIED TARGET

Observation of the target buried on work area number 2 by GEUS yield unexpected results:

1. For the buried corner reflector, the measured RCS is from 25 to 35 dB below the one expected from an uncovered corner reflector of similar orientation at H-pol and more that 27 dB below at V-pol.

2. For the Luneburg lens, that was set up at the surface, but covered under 50 cm of fresh snow after the last radar acquisitions, there is a significant drop in RCS (32 to 38 dB) between the first and the later acquisitions (when the lens was probably covered with snow). There is also, in the later images, a duplication of the echo in the range direction (fig. 4).

The reason is that a Luneburg Sphere is designed to act as a strong reflector when in a media of index close 1 (as air or vacuum), when buried in dense snow, the index of the surrounding media might be higher than the outer layer of the LSp (of index very close to 1). This has two consequences:

1. the snow/sphere interface behave as a big air bubble within the snow (hence a specular reflection at this surface).

2. the incoming waves are no more focused at the opposite metallized surface (the SPh becomes sort of farsighted) thus attenuating strongly the RCS of the target echo (which is localized somewhere beyond the bottom pole of the sphere opposite to the incoming radiation).

On fig. 4 the offset in position is due mis-compensation of stereo effect to calibrating target apparent depth: unlike the buried targets, the Luneburg sphere, intended to act as a reference point for ice drift, was emerging at the surface (at the beginning of the experiment).

MASSIVE MULTILOOKING FOR SHALLOW TARGET DETECTION

Work area number 2 was imaged from 4 headings with 2 antenna pointing and 3 horizontal separations. Acquisitions from two headings were repeated with the aim of performing repeat-pass interferometry. However, the strong volume scattering made the repeated acquisition incoherent, with the inconvenient that interferometry (e.g. for ice drift or for target depth measure) was ruined, but the unexpected consequence is that the multilooking provides an higher equivalent number of looks (fig. 5-right) thus allowing a considerable smoothing of the clutter speckle without the resolution degradation of the conventional "spatial" speckle filters. This was favourable for detection of buried target with low RCS (consequence of the strong RCS damping described above).

Figure 5 shows on the left a detail of the native resolution (20×20 cm) image. Detail area is 200×265 m on the ground. The area is selected with shallow crevices for illustration purpose.

The image on the right of fig. 5 is the ENL for all acquisitions on work zone number 2 at a 100 times larger scale (full area is 20×20 km). The detail image is situated in the upper right of the white area.

Thermal & system noise is about 20 dB below the clutter level for a SAR image computed in the middle of the main lobe in the antenna pattern, but when we compute 6 looks per channel and acquisition (a total amount of 432 looks for example on work zone number 2), the extreme looks angle correspond to image illumination of about 15% of main lobe level, and the Doppler ambiguities are some 7 times stronger. We thus mix images with marginal (0.5%) noise + ambiguities together with images with about 50% noise + ambiguities. The combination law is derived from [START_REF] Cantalloube | How to compute a multilook sar image?[END_REF], with the extra twist that the clutter RCS is known (and not recursively estimated) and thus is directly injected into the algorithm. Combination law is also slightly altered for compositing the Combination of images from several headings (illumination directions) however, requires compensation of the stereo effect (consequence of the depth of the imaged feature). Because of refraction index of the snow/ice cover, the apparent stereo depth is different from the one deduced in the absence of refraction. Simplistic modelling by an homogeneous volume of constant index representative of the upper layers (1.25 to 1.35) showed that the transverse mis-positioning is below the 20 cm resolution cell, but the apparent/actual depth varies strongly with incident angle (from the refractive index in sounding mode, increasing with grazing incidence). Without knowledge of actual refractive index variation in the upper layers during acquisitions (there was a strong snow fall of an estimated 50 cm during the 5 days separating the North/South and the East/West acquisitions at X-band) the apparent depth used for image registration was adjusted using a calibrating known target (yielding an apparent depth of 3.4 m for the earlier acquisitions at 60 • incidence and 1.9 m for the later acquisitions at 40 • incidence).

CONCLUSION

Airborne SAR imaging at X and UHF band are a valuable tool for ice-shelf sub-surface observation: Remote observation of large features at UHF band to a depth of 50 m and of small objects at X band for depth of 3 m are possible. This requires, however, a large number of acquisitions and a significant computation load to yield results comparable to surface measurement with ground penetrating radars. Due to the snow accumulation on the ice-shelf, this means that airborne SAR exploration has visibility horizon of roughly a couple of years for small features at X band and around 30 years for larger features at UHF bands. For example, the WW2 "lost squadron" present depth is estimated at around 100 m making it most probably beyond the range of this technique [START_REF] Krogager | Results from the daloex 2015 campaign with f-sar in greenland[END_REF].

Fig. 1 .

 1 Fig. 1. Target setup: Calibrating corner reflector setup in Kangerlussuaq (top). Corner reflector during burial (bottom left), Luneburg sphere (bottom centre) and Ice lens during fabrication (bottom right)

Fig. 3 .

 3 Fig. 3. Usual (side looking) SAR mode images. X band, Hpolarisation (top), L band H-polarisation (middle), UHF band R/G/B=Cross-pol/V-pol/H-pol (bottom).

Fig. 4 .

 4 Fig. 4. Single acquisition images (6 looks × 4 polar) of the Luneburg lens probably while not completely covered with snow (left) and later probably covered with about 50 cm fresh snow (right).

Fig. 5 .

 5 Fig. 5. Sample of composite image at X band showing very low clutter speckle (left). Equivalent number of looks for the full work zone number 2 (right, false colour: yellow=100 looks, white= above 200 looks).

Fig. 6 .

 6 Fig. 6. Detection of the target setup ice lens (plastic sled filled with ice at the top of the window) and an empty hole packed with snow (bottom) on a composite image of all acquired X band images.