Lei Mo

Member, IEEE Angeliki Kritikakou
email: angeliki.kritikakou@irisa.fr

Member, IEEE Olivier Sentieys
email: olivier.sentieys@irisa.fr.

Member, IEEE Xianghui Cao
email: xhcao@seu.edu.cn

† A Kritikakou

Real-time Imprecise Computation Tasks Mapping for DVFS-Enabled Networked Systems

Keywords: Networked Systems, Task Mapping, Imprecise Computation, Quality-of-Service

Networked systems are useful for a wide range of applications, many of which require distributed and collaborative data processing to satisfy real-time requirements. On the one hand, networked systems are usually resource constrained, mainly regarding the energy supply of the nodes and their computation and communication abilities. On the other hand, many realtime applications can be executed in an imprecise way, where an approximate result is acceptable as long as the baseline Quality-of-Service (QoS) is satisfied. Such applications can be modeled through Imprecise Computation (IC) tasks. To achieve a better trade-off between QoS and limited system resources, while meeting application requirements, the IC-tasks must be efficiently mapped to the system nodes. To tackle this problem, we firstly construct an IC-tasks mapping problem that aims to maximize system QoS subject to real-time and energy constraints. Dynamic Voltage and Frequency Scaling (DVFS) and multi-path routing are explored to further enhance real-time performance and reduce energy consumption. Secondly, based on the problem structure, we propose an optimal approach to perform IC-tasks mapping and prove its optimality. Furthermore, to enhance the scalability of the proposed approach, we present a heuristic IC-tasks mapping method with low computation time. Finally, the simulation results demonstrate the effectiveness of the proposed methods in terms of the solution quality and the computation time.

I. INTRODUCTION

W IRELESS Sensor and Actuator Networks are networked systems that cannot only measure the physical environment through the sensor nodes, but can also modify it, through the actions performed by the actuator nodes. This characteristic is one of the key elements of the Internet of Things (IoT) [START_REF] Chiang | Fog and IoT: an overview of research opportunities[END_REF], [START_REF] Mo | Distributed node coordination for real-time energy-constrained control in wireless sensor and actuator networks[END_REF]. Typical application requirements in such networked systems are low energy consumption, low task execution delay, and high system Quality of Service (QoS) [START_REF] Razzaque | Middleware for internet of things: a survey[END_REF]. However, enhancing system QoS often requires more energy consumption and execution time. To balance these contradictory requirements, an efficient mapping of the application tasks on the nodes is required. By properly mapping the tasks on the nodes, we can avoid sending all the data to a central controller that executes the tasks. Part of task execution can be done on-site, on the nodes that have computation and communication capabilities. As a result, only a small part of data is required to be sent, reducing the network traffic and the node workload. This model of computation is also known as "Fog/Edge-computing" [START_REF] Zahaf | Energy-efficient scheduling for moldable real-time tasks on heterogeneous computing platforms[END_REF].

In many application domains [START_REF] Aydin | Optimal rewardbased scheduling for periodic real-time tasks[END_REF], such as multimedia processing, mobile target tracking, real-time heuristic search, information gathering and control engineering, an approximate result, obtained before the deadline, is usually acceptable. For example, in video streaming, frames with a lower quality are better than totally missing frames. In target tracking, an approximate estimation of target's location in time is better than an accurate location, obtained too late. In these domains, the applications are usually modeled as Imprecise Computation (IC) tasks [START_REF] Liu | Imprecise computations[END_REF], where a task is logically decomposed into a mandatory subtask and an optional subtask. The mandatory subtask must be completed before the task deadline to generate the minimum acceptable QoS. Then, the optional subtask is executed, if there are enough free resources in the system. With the IC-tasks model, the longer the optional subtasks are executed, the better the QoS of the result. Dynamic Voltage and Frequency Scaling (DVFS) [START_REF] Chen | Energy optimization for real-time multiprocessor system-on-chip with optimal DVFS and DPM combination[END_REF] is an efficient technique that controls both voltage and frequency, thus, the energy and the time required to execute the tasks. By properly mapping IC-tasks on DVFS-enabled nodes of networked systems, the QoS can be further improved, under the limited system sources and the application requirements.

Task mapping is a well-known problem in embedded systems [START_REF] Chen | Energy optimization for real-time multiprocessor system-on-chip with optimal DVFS and DPM combination[END_REF]- [START_REF] Mendez-Diaz | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF]. However, there are few works that deal with task mapping on the nodes of networked systems [START_REF] Li | Adaptive energy-efficient scheduling for hierarchical wireless sensor networks[END_REF]- [START_REF] Mo | Mapping imprecise computation tasks on cyber-physical systems[END_REF]. The majority of these approaches focuses on precise computation tasks (i.e., the tasks without optional subtasks) and systems without DVFS capabilities, see Table I. When both IC-tasks and DVFS are considered, the way to execute IC-tasks is decided by 1) the task mapping, which refers to both the task allocation (on which node a task is executed) and the task scheduling (when each task starts to execute and how long its optional subtask is executed), and 2) the voltage and frequency of the node, which executes the task. In this context, there are three main differences between our work and the existing works:

1) Compared with the task mapping on embedded systems [START_REF] Chen | Energy optimization for real-time multiprocessor system-on-chip with optimal DVFS and DPM combination[END_REF]- [START_REF] Mendez-Diaz | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF], the task mapping on networked systems is constrained, since some tasks (e.g., sensing and control) have one-to-one allocation on the nodes, whereas the allocation constraints of other tasks (e.g., data processing) may not be restricted. 2) In embedded systems, the task communication cost (time and energy) is usually very small, compared with the task execution cost. Therefore, the impact, that the task mapping decisions have on the task communication cost, is limited. However, in networked systems, the task communication cost must be taken into account, and this cost highly depends on the task mapping decision. When tasks are allocated to different nodes, communication delay is introduced, while the nodes consume energy for transmitting the data. 3) In networked systems, the task mapping approaches usually assume that the tasks are precise [START_REF] Li | Adaptive energy-efficient scheduling for hierarchical wireless sensor networks[END_REF]- [START_REF] Zhang | A trusted real-time scheduling model for wireless sensor networks[END_REF]. However, since our work considers IC-tasks, a set of new variables (relevant to the execution of optional subtasks) has to be introduced into the task mapping problem. The objective is to maximize the QoS by increasing the execution cycles of the optional subtasks. The longer the optional subtasks are executed, the more energy and time are spent. Hence, the adjustment of optional subtasks affects the objective function as well as the energy and the time constraints of the task mapping problem.

A. Related Work

Existing task mapping approaches focus either on energyaware mapping or QoS-aware mapping. The majority of energyaware task mapping approaches usually aim at minimizing energy consumption, under system resource and application constraints [START_REF] Chen | Energy optimization for real-time multiprocessor system-on-chip with optimal DVFS and DPM combination[END_REF]- [START_REF] Leung | Simultaneous task allocation, scheduling and voltage assignment for multiple-processors-core systems using mixed integer nonlinear programming[END_REF]. In order to enhance energy efficiency, DVFS is used. An example considering independent real-time tasks is the work of [START_REF] Li | Minimizing energy consumption for frame-based tasks on heterogeneous multiprocessor platforms[END_REF], where the task allocation problem with DVFS is formulated as an Integer Linear Programming (ILP) and solved by a relaxation method based on Linear Programming (LP). When dependent real-time tasks are considered, Mixed Integer Programming (MIP) is usually used to formulate the task mapping problem [START_REF] Chen | Energy optimization for real-time multiprocessor system-on-chip with optimal DVFS and DPM combination[END_REF], [START_REF] Mahmood | Energy-aware real-time task scheduling in multiprocessor systems using a hybrid genetic algorithm[END_REF]- [START_REF] Leung | Simultaneous task allocation, scheduling and voltage assignment for multiple-processors-core systems using mixed integer nonlinear programming[END_REF]. In order to solve an ILP-based task mapping problem, a hybrid Genetic Algorithm (GA) is designed in [START_REF] Mahmood | Energy-aware real-time task scheduling in multiprocessor systems using a hybrid genetic algorithm[END_REF], and an optimal method based on Benders Decomposition (BD) [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF] is presented in [START_REF] Emeretlis | A logicbased Benders decomposition approach for mapping applications on heterogeneous multicore platforms[END_REF]. Combining DVFS and Dynamic Power Management (DPM), a Mixed-Integer Linear Programming (MILP)-based task mapping problem is formulated in [START_REF] Chen | Energy optimization for real-time multiprocessor system-on-chip with optimal DVFS and DPM combination[END_REF]. The number of variables is further reduced by problem refinement, and then, the refined MILP problem is optimally solved by the CPLEX solver. The Mixed-Integer Non-Linear Programming (MINLP)-based task mapping problem in [START_REF] Leung | Simultaneous task allocation, scheduling and voltage assignment for multiple-processors-core systems using mixed integer nonlinear programming[END_REF] is first relaxed to a MILP by linear approximation, and then, is optimally solved by the Branch and Bound (B&B) method [START_REF] Boyd | Branch and bound methods[END_REF].

Existing QoS-aware task mapping approaches usually consider the IC-task model [START_REF] Cortes | Quasi-static assignment of voltages and optional cycles in imprecise-computation systems with energy considerations[END_REF]- [START_REF] Mendez-Diaz | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF]. Their aim is to maximize system QoS, under real-time and/or energy constraints. In [START_REF] Wei | Cost-constrained QoS optimization for approximate computation real-time tasks in heterogeneous MPSoCs[END_REF], the problem of mapping independent tasks is solved, but the frequency of each processor is decided upfront, whereas task allocation and optional subtasks adjustment are solved one after the other. A similar task mapping problem is studied in [START_REF] Mo | Decomposed task mapping to maximize QoS in energy-constrained real-time multicores[END_REF], where an optimal approach is proposed based on problem decomposition. The work in [START_REF] Yu | Quality-driven dynamic scheduling for real-time adaptive applications on multiprocessor systems[END_REF] focuses on task scheduling and optional subtasks adjustment of dependent tasks, whereas task allocation is fixed and given in advance. Other existing approaches focus on bi-objective optimization, e.g., increasing the QoS while reducing the frequency changes under real-time constraints for independent tasks [START_REF] Mendez-Diaz | Energy-aware scheduling mandatory/optional tasks in multicore real-time systems[END_REF]. Dependent tasks are mapped on a DVFS-enabled uni-processor platform in [START_REF] Cortes | Quasi-static assignment of voltages and optional cycles in imprecise-computation systems with energy considerations[END_REF]. This problem is first formulated as an Integer Non-Linear Programming (INLP), and then, is relaxed to a convex problem. However, the aforementioned approaches mainly focus on embedded systems, thus, the communication cost is not taken into account. In the networked systems, the nodes also consume energy and time for data transmission.

In networked systems, approaches exist that map precise dependent tasks on the nodes of Wireless Sensor Networks (WSNs) [START_REF] Li | Adaptive energy-efficient scheduling for hierarchical wireless sensor networks[END_REF]- [START_REF] Zhang | A trusted real-time scheduling model for wireless sensor networks[END_REF]. For instance, the work of [START_REF] Li | Adaptive energy-efficient scheduling for hierarchical wireless sensor networks[END_REF] minimizes the overall energy consumption and balances the workload of the system while meeting task deadline through a three-phase heuristic. The task allocation problem with metric -balance the energy consumption of the nodes -is considered in [START_REF] Pathak | Energy-efficient task mapping for datadriven sensor network macroprogramming[END_REF]. This problem is first formulated as an INLP, and then, it is transformed to an ILP. Finally, the transformed problem is solved by a greedy algorithm. In [START_REF] Li | Efficient allocation of resources in multiple heterogeneous wireless sensor networks[END_REF], [START_REF] Edalat | An auction-based strategy for distributed task allocation in wireless sensor networks[END_REF], the lifetime of system is maximized (i.e., minimize the energy consumption of the node with the lowest energy level) by allocating dependent tasks on the nodes. On this basis, a heuristic method [START_REF] Li | Efficient allocation of resources in multiple heterogeneous wireless sensor networks[END_REF] and a game theory method [START_REF] Edalat | An auction-based strategy for distributed task allocation in wireless sensor networks[END_REF] are designed to solve task allocation problems, respectively. However, DVFS is not taken into account in the above approaches. The mapping of dependent tasks and DVFS are jointly addressed in [START_REF] Tian | Cross-layer collaborative in-network processing in multihop wireless sensor networks[END_REF], where a twophase heuristic is proposed. In [START_REF] Zhang | A trusted real-time scheduling model for wireless sensor networks[END_REF], the energy consumption for communication and task execution is minimized, by allocating the tasks to DVFS-enabled nodes through the ant colony and the bee colony algorithms. However, since the aforementioned approaches mainly focus on the precise computation tasks, the improvement of system QoS is limited. In [START_REF] Mo | Mapping imprecise computation tasks on cyber-physical systems[END_REF], the dependent IC-tasks mapping problem is formulated as a MILP with the goal of maximizing system QoS and optimally solved by a BD-based approach. Compared with our preliminary results in [START_REF] Mo | Mapping imprecise computation tasks on cyber-physical systems[END_REF], the current work takes DVFS and multi-path routing into account and proposes a novel heuristic method to further reduce the computation time. By using DVFS and multi-path routing, we can achieve a better trade-off between real-time performance and energy efficiency.

B. Contributions

Complementary to the state-of-the-art, this work solves the following mapping problem: given a networked system with DVFS-enabled wireless nodes and a set dependent IC-tasks, we aim to determine: 1) which node should the task be executed on;

2) what voltage/frequency level should be used for each task;

3) what is the execution sequence of the tasks on each node; 4) which paths transmit the data required for the task execution; 5) when should a task start its execution; and 6) how many cycles of the optional subtasks are needed to be executed, such that the system QoS is maximized. while meeting the energy supply and the task deadline constraints. In this context, the task mapping decides task allocation, frequency assignment, task sequence, multi-path routing, task start time and task adjustment. Our main contributions are summarized as follows:

1) We formulate the IC-tasks mapping problem that simultaneously optimizes task allocation, frequency assignment, task sequence, multi-path routing, task start time and task adjustment as a MINLP. The objective is to maximize QoS without violating the real-time and energy constraints. 2) We prove that, by replacing the nonlinear items with some auxiliary variables and adding additional linear constraints into the optimization problem, the MINLP problem can be equivalently transformed to a MILP. 3) Based on the BD framework and the idea of closed-loop control, we propose an Optimal Task Mapping (OTM) algorithm to efficiently solve the transformed problem. This method decomposes the transformed problem into two subproblems, each with fewer constraints and variables. The first subproblem is an ILP, which is responsible for task allocation, frequency assignment, task sequence and multi-path routing. The second subproblem is an LP, which determines the start time and the optional cycles of each task. The proposed OTM algorithm iterates using the solutions of these two subproblems. We prove that OTM converges to the optimal solution of the transformed problem. 4) In order to enhance the scalability of the proposed approach, we present a Heuristic Task Mapping (HTM) algorithm, which reduces the computation time of OTM algorithm by solving the subproblems in a sequence. 5) Finally, we provide extensive experimental results to evaluate the performance of the proposed optimal and heuristic task mapping algorithms. The obtained results show significant performance improvements compared with the state-of-the-art task mapping methods in terms of solution quality and computation time. The rest of this paper is organized as follows. Section II presents the system model and formulates the problem under study. Section III, Section IV and Section V describe the problem linearization method, the optimal and the heuristic task mapping algorithms, respectively. Finally, Section VI shows the simulation results and Section VII concludes this work.

II. SYSTEM MODELS AND PROBLEM FORMULATION

A. Motivational Example

Let's consider a Heating, Ventilation, and Air Conditioning (HVAC) system [START_REF] Vedantham | Hazard avoidance in wireless sensor and actor networks[END_REF] as an example. This system includes eight wireless nodes {θ 1 , . . . , θ 8 }, as shown in Fig. 1. The nodes θ 1 and θ 8 equipped with temperature sensors. They monitor the average temperature within their sensing range and use this information to determine if a fire has occurred. The nodes θ 3 and θ 7 are equipped with the humidity sensors, able to measure the humidity in the air. The node θ 4 is connected to a sprinkler, capable of extinguishing a fire within its acting range. The node θ 5 is connected to a fan, capable of reducing the humidity in the air. The nodes equipped with temperature sensor, humidity sensor, sprinkler and fan are marked with S t , S h , A s and A f , respectively. The HVAC system periodically measures the temperature and the humidity of the environment, compare the measured values with given thresholds, and decide the corresponding actions. When a fire has been detected, the sprinkler must be activated first, and then, the fan. If these actuators are activated in the opposite order, it will cause serious problems. In the HVAC system, a task can be of sensing, processing, or acting type. For example, τ 1 and τ 2 are the temperature measurement tasks, while τ 3 and τ 4 are the humidity measurement tasks. These sensing tasks generate a set of data, with sizes s 15 , s 25 , s 35 , and s 45 , that have to be processed by task τ 5 . After processing, task τ 5 transmits the result to task τ 6 , which determines the control action to be taken. Finally, tasks τ 7 and τ 8 act upon the decision generated by task τ 6 and control the actions of the actuators. The example of Fig. 1 illustrates the task mapping problem addressed in this work.

• Firstly, the allocations of sensing tasks {τ 1 , τ 2 , τ 3 , τ 4 } and acting tasks {τ 7 , τ 8 } are restricted only to the nodes that have the corresponding capabilities. However, the data processing tasks {τ 5 , τ 6 } can be allocated on any node. The dependencies between the tasks define the task graph. Considering that the communication range of each node is limited, the communication between the nodes is modeled by a node graph. Both task graph and node graph influence the task allocation decision. • Secondly, each task τ i is composed of a mandatory subtask and an optional subtask. The optional subtask is executed immediately after its corresponding mandatory subtask. The QoS of the result highly depends on the number of execution cycles of optional subtask. When deciding the task allocation and the task scheduling (i.e., the start time and the end time of each task), the following goals must be achieved at the same time: 1) meet the deadline of the tasks (i.e., each task must be completed within a predefined deadline to generate a result in time), 2) the energy consumption of the nodes (i.e., the energy spent for computation, communication, sensing and acting) should not exceed their energy supply, and 3) increase the execution cycles of optional subtasks to maximize the QoS. These three goals are important, but they contradict with each other, since the real-time and energy constraints may require to sacrifice the system QoS.

B. System Model 1) Task Model: We consider a task set T of N periodic real- time IC-tasks {τ 1 , . . . , τ N }. Each task τ i is described by a tuple {o i , M i , O i , t s i , d i , l i }. o i , M i and O i are measured in Worst Case Execution Cycles (WCEC). o i is the actual WCEC of the optional subtask, M i is the WCEC of the mandatory subtask, and O i is the maximum WCEC of the optional subtask o i (i.e., 0 ≤ o i ≤ O i ). t s
i and d i are the start time and the deadline of task τ i , respectively. l i is the period of task τ i , which is also equal to the scheduling horizon H. We introduce a linear QoS function f i (o i ) for each task τ i , where the generated QoS increases uniformly with the number of optional cycles, decided for actual execution [START_REF] Aydin | Optimal rewardbased scheduling for periodic real-time tasks[END_REF].

We assume that the tasks are dependent and non-preemptive. They are released at the same time 0 and share a common scheduling horizon H. The task set T is modeled by a Directed Acyclic Graph (DAG) G t (V t , E t ), where the vertexes V t and the edges E t represent the set of tasks to be executed and the data dependencies between the tasks, respectively. We consider that the system operates in rounds, which means in each round (i.e., scheduling horizon H) all the tasks in task set T are executed once. Based on the dependency between the tasks, we introduce a task execution order matrix p = [p ij ] N ×N . If p ij = 1, task τ i precedes task τ j and τ j is the closest task of τ i , otherwise, p ij = 0. As the task graph shown in Fig. 1, we have p 15 = p 25 = p 35 = p 45 = p 56 = p 67 = p 68 = 1.

2) Energy Model: We consider a networked system that contains M wireless DVFS-enabled nodes {θ 1 , . . . , θ M }, where M s nodes {θ 1 , . . . , θ Ms } equipped with sensors and M a nodes {θ Ms+1 , . . . , θ Ms+Ma } connected to actuators (M s +M a ≤ M ). The processor of each node has L different Voltage/Frequency (V/F) levels {(v 1 , f 1 ), . . . , (v L , f L )}. The power consumption of the processor under the l th V/F level is computed as

P c l = P s l + P d l , (1)
where P s l = v l K 1 e K2v l e K3v bs + |v bs |I j is the static power, and

P d l = C ef f v 2 l f l is the dynamic power [7]. The constants K 1 ,
K 2 and K 3 are technology dependent. I j is the approximately constant junction leakage current. v bs is the reverse bias voltage used to reduce the leakage power and can be treated as constant. C ef f is the average effective switching capacitance. We assume that the processor of node θ k can operate in two modes: active mode with the power consumption P c k and idle mode with the power consumption P 0 k [START_REF] Li | Minimizing energy consumption for frame-based tasks on heterogeneous multiprocessor platforms[END_REF]. The processor goes into idle mode immediately when it has no task to execute. The transition time and energy is incorporated into the task execution time and energy, since they are very small compared to those required to execute a task [START_REF] Li | Minimizing energy consumption for frame-based tasks on heterogeneous multiprocessor platforms[END_REF]. For each task τ i , τ i starts and ends its execution on the same processor (i.e., no task migration), and the frequency cannot be changed during the execution of the task (i.e., inter-task DVFS).

3) Data Routing Model: Since each node can only communicate directly with the nodes that are within its communication range, we introduce a node graph G n (V n , E n ), where the vertexes V n represent the nodes, while the edges E n represent the communication cost between the nodes. In this paper, we consider multi-path data routing, which means a pair of nodes can communicate with each other through several paths. Specifically, we consider two routing options: the 1 st routing path is energy-oriented, while the 2 nd routing path is time-oriented, since the system under study is energy and time constrained. For the energy (time)-oriented path, an edge represents the energy (time) required for transmitting and receiving a unit of data between the corresponding two nodes. Therefore, the aim of energy (time)-oriented routing is to find the shortest path. The shortest path can be easily found through the existing methods, such as Dijkstra's algorithm [START_REF] Zhang | Distributed workload dissemination for makespan minimization in disruption tolerant networks[END_REF]. As the examples illustrated in Fig. 2, e ij and t ij is the weight of the edge between the nodes θ i and θ j in the energy-oriented and time-oriented graph, respectively. If a data is required to be transmitted from θ 1 to θ 2 , we have θ 1 → θ 3 → θ 2 for the energy-oriented path, and θ 1 → θ 4 → θ 2 for the timeoriented path. To describe the time and the energy related to path selection, we introduce 1) a routing energy matrix r = [r βγkh ] N ×N ×N ×2 , where r βγkh represents the energy consumed per unit of data at node θ k , when routing the messages from θ β to θ γ through the h th path, and 2) a routing time matrix t = [t βγh ] N ×N ×2 , where t βγh denotes the time required to transmit unit of data from θ β to θ γ through the h th path. Note that during the data transmission between θ β to θ γ , if θ k is not included in the h th path, we have r βγkh = +∞, e.g., r 1241 = +∞ in Fig. 2(a).

e

C. Problem Formulation

The IC-tasks mapping problem has the objective of maximizing system QoS subject to real-time and energy constraints, by determining task allocation, frequency assignment, task sequence, multi-path routing, task start time and task adjustment. In order to formulate this problem, we introduce the following variables: 1) binary variable q ik = 1 if task τ i is allocated to node θ k , otherwise, q ik = 0; 2) binary variable c il = 1 if task τ i is executed with frequency f l , otherwise, c il = 0; 3) binary variable u ij = 1 if task τ i proceeds τ j , otherwise, u ij = 0; 4) binary variable x eh = 1 if data over edge e in E t is routed along the h th path, otherwise, x eh = 0; 5) continuous variable t s i denotes the start time of task τ i ; 6) continuous variable o i represents the optional subtask of task τ i . We assume that the following matrices are known: 1) task execution order matrix p, 2) data size matrix s, 3) routing energy matrix r, and 4) routing time matrix t. The parameters and the variables used in the problem formulation are summarized in Table II. For tractability reasons, during the problem formulation, we consider o i as continuous variables. When the problem is solved, we round the result down. Since the tasks are typically executed in hundreds of thousands of cycles, one cycle is a very fine-grained unit [START_REF] Cortes | Quasi-static assignment of voltages and optional cycles in imprecise-computation systems with energy considerations[END_REF].

Let N = {1, . . . , N }, M = {1, . . . , M }, M s = {1, . . . , M s }, M a = {M s + 1, . . . , M s + M a } and H = {1, 2}.
The constraints involved in the task mapping problem can be explained as follows.

1) Task Allocation Constraint: Since each task τ i is assigned to only one node and the task allocation decisions related to the sensors and the actuators are restricted, we have

k∈K q ik = 1, ∀i ∈ N , (2)
q ik = 1, ∀(i, k) ∈ M m , (3)
where M m denote the set of index pairs (tasks and nodes) that have fixed matching. As the example shown in Fig. 1, since q 11 = q 28 = q 33 = q 47 = q 75 = q 84 = 1, we have M m = {(1, 1), (2, 8), [START_REF] Razzaque | Middleware for internet of things: a survey[END_REF][START_REF] Razzaque | Middleware for internet of things: a survey[END_REF], (4, 7), [START_REF] Chen | Energy optimization for real-time multiprocessor system-on-chip with optimal DVFS and DPM combination[END_REF][START_REF] Aydin | Optimal rewardbased scheduling for periodic real-time tasks[END_REF], (8, 4)}.

2) Frequency Assignment Constraint: Since each task τ i is executed with only one V/F level, we get

l∈L c il = 1, ∀i ∈ N . (4)
3) Data Routing Constraint: Since the data over edge e in E t is routed through one path, we obtain h∈H x eh = 1, ∀e = (i, j) ∈ E t .

(5)

4) Task Sequence Constraints: Before executing a task τ j , we need to collect all the data generated from its previous dependent tasks. If p ij = p kj = 1, tasks τ i and τ k precede τ j and they are the tasks closest to τ j . When these tasks are allocated to different nodes, the node that executes τ j should collect the data from the other nodes. In order to avoid communication collision [START_REF] Lee | Towards energy-optimal and reliable data collection via collision-free scheduling in wireless sensor networks[END_REF], the data are received in sequence, since one node cannot receive the data from multiple nodes simultaneously. For the dependent tasks, e.g., τ i and τ j in Fig. 3 (p ij = 1), if we allocate these tasks to different nodes θ β

p ij =     
1 if task τ i proceeds τ j and τ j is the nearest task of τ i 0 else s ij size of data that task τ i produces for task τ j r βγkh energy consumed of node θ k when routing unit of data from θ β to θ γ through the h th path t βγh time required to transmit unit of data from θ β to θ γ through the h th path Binary Variables

q ik = 1 if task τ i is allocated to node θ k 0 else c il = 1 if task τ i is executed with frequency f l 0 else u ij = 1 if task τ i proceeds task τ j 0 else x eh =     
1 if data over edge e in E t is routed along with the h th path 0 else

Continuous Variables

o i optional cycles of task τ i t s i start time of task τ i and θ γ (i.e., q iβ = q jγ = 1) and choose the h th path (i.e., x eh = 1, e = (i, j) ∈ E t ) to transmit the data from θ β to θ γ , the time required for θ γ to receive the data from θ β is s ij p ij q iβ q jγ x eh t βγh . Therefore, the time spent for receiving the data required by the execution of task τ j is t r j = e∈Et β∈M γ∈M h∈H

s ij p ij q iβ q jγ x eh t βγh , ∀j ∈ N . ( 6
)
For the dependent tasks, e.g., τ i and τ j in Fig. 3 (p ij = 1), no matter if they are allocated to the same node or to different nodes, the execution sequence between these tasks is fixed. Therefore, the start and end time of the tasks is bounded by

t s j ≥ p ij t e i + t r j , ∀i, j ∈ N , i = j, (7)
where

t c i = l∈L c il Mi+oi f l
and t e i = t s i + t c i are the execution time and the end time of task τ i , respectively. If p ij = 1, we have t s j ≥ t e i + t r j , else, ( 7) is always satisfied. On the other hand, for the independent tasks, e.g., τ j and τ k in Fig. 3 (p jk = 0), if they are allocated to different nodes, the execution of τ j doesn't affect the execution of τ k . However, if τ j and τ k are allocated to the same node, we need to schedule these tasks so as to make sure that their executions do not overlap with each other, since one processor executes no more than one task at the same time. To this end, we introduce the following constraint:

𝜏 j 𝜏 k 𝜏 l 𝜏 i
t e i + t r j ≤ t s j + (2 -q ik -q jk )H + (1 -u ij )H, ∀i, j ∈ N , i = j, ∀k ∈ M. (8)
If τ i and τ j are executed on the different nodes (i.e., q ik + q jk ≤ 1), (8) is always satisfied, else, ( 8) is meaningful: with the task allocation decision q ik = q jk = 1, if u ij = 1 (i.e., τ i precedes τ j ), we have t e i + t r j ≤ t s j ; if u ij = 0, we get u ji = 1 due to u ji + u ij = 1. Note that (8) can be rewritten as t e j + t r i ≤ t s i +(2-q jk -q ik )H +(1-u ji )H = t s i +(2-q jk -q ik )H +u ij H. Therefore, we have t e j + t r ≤ t s i . 5) Task Deadline Constraint: Since each task τ i must be finished within a given time threshold d i , we have

t e i ≤ d i , ∀i ∈ N . (9)
6) Energy Supply Constraint: With the task allocation decision q ik , the time required to execute all the tasks assigned to node θ k is i∈N q ik t c i , and the idle time of node θ k being in the scheduling horizon H is H -i∈N q ik t c i . Based on the energy model (1), the computation energy of node θ k in each round is

E c k = i∈N q ik e c i + H - i∈N q ik t c i P 0 k , ∀k ∈ M. ( 10
)
where e c i = l∈L c il Mi+oi f l P c l is the energy required to execute task τ i . On the other hand, the energy consumed for node θ k to transmit data in each round is

E t k = e∈Et β∈M γ∈M h∈H s ij q iβ q iγ x eh r βγkh , ∀k ∈ M. ( 11
)
Since the total energy consumed by node θ k during the scheduling horizon H cannot exceed the energy supply E l k , we have

E c k + E t k + E s k ≤ E l k , ∀k ∈ M s , (12)
E c k + E t k + E a k ≤ E l k , ∀k ∈ M a , (13)
E c k + E t k ≤ E l k , ∀k ∈ M, ∀k / ∈ M s , M a , (14)
where E s k is the sensing energy consumption and E a k is the acting energy consumption.

Taking the objective (maximize QoS function i∈N f i (o i )) and all the constraints mentioned above into account, the Primal Problem (PP) can be formulated as

PP : min q,c,u,x, t s ,o - i∈N f i (o i ) (15)
s.t.

     (2) -(14), q ik , c il , u ij , x eh ∈ {0, 1}, 0 ≤ t s i ≤ H, 0 ≤ o i ≤ O i .

III. PROBLEM LINEARIZATION

Since the products q ik c il , q iβ q jγ x eh , c il o i and q ik c il o i are included in the constraints ( 6)-( 14), ( 15) is a MINLP problem, which is hard to solve directly. Note that q ik , c il and x eh are binary variables, while o i is a continuous variable. In order to linearize above nonlinear terms, we introduce the following lemmas.

Lemma 3.1: Let b 1 , b 2 and g denote the binary variables. The nonlinear constraint g = b 1 b 2 can be equivalently replaced by the linear constraints:

g ≤ b 1 , g ≤ b 2 and g ≥ b 1 + b 2 -1.
All the proofs of lemmas and theorems are presented in the Appendices for better readability of the main manuscript.

Based on Lemma 3.1, we first introduce an auxiliary (binary) variable g ikl to replace the nonlinear term q ik c il . Then, we add the following constraints into the PP:

{g ikl ≤ q ik , g ikl ≤ c il , g ikl ≥ q ik + c il -1}, ∀i ∈ N , ∀k ∈ M, ∀l ∈ L. (16)
Similarly, the nonlinear term q iβ q jγ in q iβ q jγ x eh is first replaced by the auxiliary variable b iβjγ and the following constraints:

{b iβjγ ≤ q iβ , b iβjγ ≤ q jγ , b iβjγ ≥ q iβ + q jγ -1}, ∀e ∈ E t , ∀β, γ ∈ M. (17)
Then, the nonlinear term b iβjγ x eh is replaced by the auxiliary variable w iβjγh and the following constraints:

{w iβjγh ≤ b iβjγ , w iβjγh ≤ x eh , w iβjγh ≥ b iβjγ +x eh -1}, ∀e ∈ E t , ∀β, γ ∈ M, ∀h ∈ H. (18)
Note that q ik c il o i = g ikl o i , where g ikl is a binary variable, while o i is a continuous variable bounded by 0 ≤ o i ≤ O i . To deal with the nonlinear terms c il o i and g ikl o i , we introduce the following lemma.

Lemma 3.2: The spaces

S 1 = {[y, b, x]|y = bx, -s 1 ≤ x ≤ s 2 } and S 2 = {[y, b, x]| -bs 1 ≤ y ≤ bs 2 , y + bs 1 -x -s 1 ≤ 0, y -bs 2 -x + s 2 ≥ 0} are equivalent,
where x is a continuous variable, b is a binary variable, and s 1 , s 2 > 0 are constants.

Based on Lemma 3.2, the nonlinear terms c il o i and g ikl o i are replaced by the auxiliary (continuous) variables y il and z ikl , and the following constraints:

{y il ≤ c il O i , y il -o i ≤ 0, y il -c il O i -o i + O i ≥ 0}, ∀i ∈ N , ∀l ∈ L, (19) {z ikl ≤ g ikl O i , z ikl -o i ≤ 0, z ikl -g ikl O i -o i + O i ≥ 0}, ∀i ∈ N , ∀k ∈ M, ∀l ∈ L. (20)
Therefore, with the auxiliary variable w iβjγh , ( 6) and [START_REF] Leung | Simultaneous task allocation, scheduling and voltage assignment for multiple-processors-core systems using mixed integer nonlinear programming[END_REF] are rewritten as follows:

t r j = e∈Et β∈M γ∈M h∈H s ij p ij w iβjγh t βγh , ∀j ∈ N , (21)
E t k = e∈Et β∈M γ∈M h∈H s ij w iβjγh r βγkh , ∀k ∈ M. ( 22
)
On this basis, with the auxiliary variables g ikl , y il and z ikl , ( 7)-( 10) are reformulated as follows:

t s j ≥ p ij t s i + l∈L c il M i + y il f l +t r j , ∀i, j ∈ N , i = j, (23)
t s i + l∈L c il M i + y il f l +t r j ≤ t s j +(2-q ik -q jk )H +(1-u ij )H, ∀i, j ∈ N , i = j, ∀k ∈ M, (24)
t s i + l∈L c il M i + y il f l ≤ d i , ∀i ∈ N , (25)
E c k = HP 0 k + i∈N l∈L g ikl M i + z ikl f l (P c l -P 0 k ), ∀k ∈ M. (26)
Finally, ( 15) is transformed to the following MILP problem:

PP1 : min q,c,u,x,g,b,w, o,t s ,y,z

- i∈N f i (o i ) (27)
s.t.

    

(2) -( 5), ( 12) -( 14), ( 16) -(26), q ik , c il , u ij , x eh , g ikl , b iβjγ , w iβjγh ∈ {0, 1}, 0 ≤ o i , y il , z ikl ≤ O i , 0 ≤ t s i ≤ H. Remark 3.1: Lemma 3.1 and Lemma 3.2 show that the linearization does not change the feasible region and the objective function of the problem. Therefore, solving PP and solving PP1 are equivalent.

IV. OPTIMAL TASK MAPPING ALGORITHM

In this section, we present an OTM algorithm to optimally solve the PP1. Our method decomposes the PP1 into two smaller subproblems with less variables: a Master Problem (MP) for task allocation q, frequency assignment c, task sequence u, and data routing x, and a Slave Problem (SP) for task start time t s and optional cycles adjustment o. Then, the optimal solution of PP1 is found by solving the subproblems iteratively. With this decomposition structure, the computational complexity can be greatly reduced.

For the sake of presentation, ( 27) is reformulated as follows:

PP2 : min x,y Φ = f T y (28) s.t. Ax b 1 , Cx + Dy b 2 ,
where x and y represent the vector of binary and continuous variables, respectively. The vector f represents the coefficients in the objective function. The matrices A, C, D and the vectors b 1 , b 2 represent the coefficients in the constraints.

A. MP and SP formulation

According to the structure of PP2, at the k th iteration, the MP and the SP have the forms:

MP : Φ l (k) = min x, Φ Φ (29) s.t. Ax b 1 , C 1 , C 2 , SP : Φ u (k) = min y 0 f T y (30) s.t. Cx(k) + Dy b 2 ,
where x(k) is the MP solution at the k th iteration,

C 1 : Φ ≥ µ(i) T (Cx -b 2 ), ∀i ∈ A, C 2 : 0 ≥ μ(j) T (Cx -b 2 ), ∀j ∈ B.
Lemma 4.1: Φ l (k) and Φ u (k) are the lower and upper bounds on Φ * , respectively, where Φ * is the optimal value of Φ.

The MP considers all the binary variables x, and the associated part of the constraints (i.e., Ax b 1 ). It also includes the information regarding the SP through a set of constraints (i.e., C 1 and C 2 ) called Benders cuts. Note that the objective function of the PP2 only considers the continuous variables y. To facilitate the iterations between the MP and the SP, an auxiliary (continuous) variable Φ is introduced into the MP as the objective function, where Φ and Φ have the same physical meaning. At each iteration, based on the SP solution, a new constraint is added into C 1 or C 2 to reduce the gap between the upper and the lower bounds. The problem iterations stop when Φ u (k) -Φ l (k) ≤ , where is a small positive tolerance.

B. Iterations Between MP and SP

Initially, we set the iteration counter k = 0 and the MP solution x(0) that satisfies Ax(0) b 1 . In addition, we assume that the constraint sets C 1 and C 2 are null, the lower bound Φ l (0) = -∞ and the upper bound Φ u (0) = +∞.

Since the SP is a LP problem, the strong duality [START_REF] Boyd | Convex optimization[END_REF] exists between the SP and its dual problem (DSP), i.e., solving SP and solving DSP are equivalent. Instead of solving the SP directly, we solve its dual problem:

max µ 0 µ T (Cx(k) -b 2 ) (31)
Φ ≥ µ(k) T (Cx -b 2 ), (32)
is added into C 1 . On the other hand, if (31) has an unbounded solution, i.e., µ(k) T (Cx(k) -b 2 ) = +∞, (30) is infeasible under given x(k). Therefore, B ← {k}∪B and a new constraint:

0 ≥ μ(k) T (Cx -b 2 ), (33)
is added into C 2 , where μ(k) is the optimal solution to (45) at the k th iteration. Based on the solution of (31), at the k th iteration, the upper bound is updated by

Φ u (k) = min{Φ u (k -1), µ(k) T (Cx(k) -b 2 )}. ( 34
)
Finally, with the updated constraint sets C 1 and C 2 , ( 29) is solved again to obtain a solution x(k + 1) for the next round iteration.

Theorem 4.1: Constraints ( 32) and ( 33) exclude the nonoptimal and infeasible solutions of the binary variables x, respectively.

Lemma 4.2: With updating the constraints in C 1 and C 2 , the lower bound sequence {Φ l (0), . . . , Φ l (k)} is increasing while the upper bound sequence {Φ u (0), . . . , Φ u (k)} is decreasing.

Theorem 4.2: The solution found by OTM converges to the global optimal one within a finite number of iterations.

V. HEURISTIC TASK MAPPING ALGORITHM

At each iteration, a new constraint is generated and added into the MP. Therefore, with an increasing number of iterations, the computational complexity and the size of MP both increase. In order to enhance the scalability of the proposed approach, we provide a novel heuristic approach HTM to efficiently solve the PP1. The basic idea of HTM is similar to OTM: both of them are based on problem decomposition. However, the HTM solves the master and the slave problems in sequence (without iteration), thus, the HTM contains two steps. During the first step, we only consider the mandatory subtasks. By balancing the energy consumption of the nodes, while meeting deadline constraints, we obtain the feasible task allocation, frequency assignment, task sequence, and data routing decisions. Under these decisions, in the second step, we determine the start time of the tasks and the cycles of optional subtasks so as to maximize QoS function.

A. Task Allocation Problem

Since the mandatory subtasks must be executed, we initially only consider the allocation, the frequency, the sequence and the routing of the mandatory subtasks (i.e., o i = 0, ∀i ∈ N ). The difficulty during the problem formulation is how to deal with the continuous variables {t s 1 , . . . , t s N } (i.e., the start time of the tasks). If the task start time is determined, the PP1 reduces to an ILP problem, since o i = 0. (9) shows that the task deadline constraint must be satisfied (i.e., t e i ≤ d i ). In the worst case, we have t e i = d i . Thus, the start time of task τ i is

t s i = d i - l∈L c il M i f l , ∀i ∈ N . (35)
On the other hand, under different optional cycles o i , task τ i generates a set of data with size s ij for task τ j . Therefore, the time required to collect the data for the execution of task τ j is given by [START_REF] Edalat | An auction-based strategy for distributed task allocation in wireless sensor networks[END_REF], and the energy consumed for node θ k to transmit the task data is calculated as [START_REF] Zhang | A trusted real-time scheduling model for wireless sensor networks[END_REF]. With o i = 0 and g ikl = q ik c il , ( 7), ( 8) and [START_REF] Emeretlis | A logicbased Benders decomposition approach for mapping applications on heterogeneous multicore platforms[END_REF] have the forms:

t s j ≥ p ij t s i + l∈L c il M i f l + t r j , ∀i, j ∈ N , i = j, (36)
t s i + l∈L c il M i f l + t r j ≤ t s j + (2 -q ik -q jk )H + (1 -u ij )H, ∀i, j ∈ N , i = j, ∀k ∈ M, (37)
E c k = HP 0 k + i∈N l∈L g ikl M i f l (P c l -P 0 k ), ∀k ∈ M. (38)
Hence, the total energy consumption of node θ k is

E all k =      E c k + E t k + E s k , ∀k ∈ M s , E c k + E t k + E a k , ∀k ∈ M a , E c k + E t k , ∀k ∈ M, ∀k / ∈ M s , M a .
In order to make sure that we have enough energy to enhance the QoS in the next step, the aim of the first step is to balance the energy consumption of the nodes. Therefore, the Task Allocation Problem (TAP) can be formulated as the following ILP:

TAP : min q,c,u,x, g,b,w max ∀k∈M E all k E l k (39) s.t.
(2) -( 5), ( 12) -( 14), ( 16) -( 18), ( 21), ( 22), ( 35)

-(38), q ik , c il , u ij , x eh , g ikl , b iβjγ , w iβjγh ∈ {0, 1}.
When solving the TAP, the polynomial-time method such as the Feasibility Pump (FP) method [START_REF] Fischetti | The feasibility pump[END_REF] can be used to reduce the computational complexity. In particular, we introduce a set of auxiliary (continuous) variables {q, ĉ, û, x, ĝ, b, ŵ}, which are assumed to be within the range [0, 1], to replace the original integer variables {q, c, u, x, g, b, w}. Then, we solve the relaxed TAP (LP) and round the solution to the nearest binary matrix that is feasible to the TAP.

B. Task Scheduling Problem

Based on the solution of TAP, we obtain the decisions for the allocation, the frequency, the sequence and the routing of the mandatory subtasks (i.e., q, c, u and x). When the values of binary variables are determined, the PP1 reduces to a LP problem. Note that the continuous variables, i.e., the task start time t s i and the optional subtask cycles o i , affect the constraints ( 7)- [START_REF] Mo | Decomposed task mapping to maximize QoS in energy-constrained real-time multicores[END_REF]. With the aim to maximize QoS under real-time and energy constraints, the Task Scheduling Problem (TSP) can be formulated as

TSP : min o,t s - i∈N g i (o i ) (40) s.t.
                             t s j ≥ p ij t e i + t r j , ∀i, j ∈ N , i = j, t e i + t r j ≤ t s j + (2 -q ik -q jk )H + (1 -u ij )H, ∀i, j ∈ N , i = j, ∀k ∈ M, t e i ≤ d i , ∀i ∈ N , E c k + E t k + E s k ≤ E l k , ∀k ∈ M s , E c k + E t k + E a k ≤ E l k , ∀k ∈ M a , E c k + E t k ≤ E l k , ∀k ∈ M, ∀k / ∈ M s , M a , 0 ≤ o i ≤ O i , 0 ≤ t s i ≤ H.
where

t c i = l∈L c il Mi+oi f l
and t e i = t s i + t c i . E c k and E t k are given by ( 10) and [START_REF] Leung | Simultaneous task allocation, scheduling and voltage assignment for multiple-processors-core systems using mixed integer nonlinear programming[END_REF], respectively.

For the HTM algorithm, the binary variables {q, c, u, x} and the continuous variables {o, t s } are solved in the TAP and the TSP separately. Therefore, we can avoid adding the auxiliary variables {y, z} and the additional constraints [START_REF] Pathak | Energy-efficient task mapping for datadriven sensor network macroprogramming[END_REF] and [START_REF] Li | Efficient allocation of resources in multiple heterogeneous wireless sensor networks[END_REF] into the problem to linearize the nonlinear items c il o i and q ik c il o i . In addition, the problem sizes of TAP and TSP are fixed and these problems are solved only once. Therefore, compared with OTM algorithm, the computation time of HTM algorithm can be greatly reduced.

Remark 5.1: The structure of PP1 shows that the feasibility of this problem is mainly determined by the real-time and the energy related constraints, but these constraints are relaxed in the TAP (i.e., t e i = d i and o i = 0). Therefore, if the PP1 is feasible, the TAP is also feasible, and further, if the TAP is feasible, the TSP is also feasible. This is because in the worst case, we have t s i = d i -t c i and o i = 0, which is consistent with the assumption of the TAP.

VI. SIMULATION RESULTS

For evaluating our approach, we consider a HVAC system with 25 nodes, where 10 nodes equipped with sensors and 5 nodes equipped with actuators. The values and the tuned parameters of the experimental set-up are summarized in Table III. The processor of the node is based on 70 nm technology [START_REF] Chen | Energy optimization for real-time multiprocessor system-on-chip with optimal DVFS and DPM combination[END_REF], and the WCECs of the mandatory subtask M i and the maximum optional subtask O i are calculated from the MiBench and the MediaBench benchmarks [START_REF] Ramkumar | Temperature aware task sequencing and voltage scaling[END_REF]. We consider a linear QoS function g i (o i ) = o i [START_REF] Wei | Cost-constrained QoS optimization for approximate computation real-time tasks in heterogeneous MPSoCs[END_REF] and assume that all data items are unit size (i.e., s ij = 1) [START_REF] Pathak | Energy-efficient task mapping for datadriven sensor network macroprogramming[END_REF].

In order to set a deadline d i of task τ i , we introduce a temporary data receiving time tr i and a relative deadline di . If p ij = 1, τ i precedes τ j and τ j is the closest task of τ i , d j = max ∀i,pij =1 {d i + tr i + di }. If τ j is an entry task, d j = dj . Let min ∀β,γ,h {t βγh } and max ∀β,γ,h {t βγh } denote the minimum and the maximum time required to transmit the data between the nodes. Since the number of tasks precede task τ i is j∈N p ji , the temporary data receiving time tr i is assumed to be within the range [ tr i,min , tr i,max ], where tr i,min and tr i,max are the minimum and the maximum time required to transmit the data related to τ j 's previous dependent tasks, respectively. The relative deadline di is assumed to be within the range [ di,min , di,max ], where di,min and di,max are the minimum and the maximum time required to execute a task with M i + O i cycles, respectively. We assume that the scheduling horizon H = max ∀i {d i } and

E h k = E t k,min + E c k,min
, where E t k,min and E c k,min are the minimum energy required for the node θ k to transmit all the task data and execute all the mandatory subtasks, respectively. The energy supply of node θ k is set to E l k = ηE h k , where η ∈ [0, 1] is an energy efficiency factor. For the sensor and the actuator nodes, we set E l k = 2ηE h k and E l k = 3ηE h k , respectively, since these nodes consume more energy for sensing and acting tasks. Note that different processor and task parameters lead to different values in the parameters {A, C, D, f , b 1 , b 2 } for the PP1. However, the problem structures under different values of parameters are the same. Therefore, the proposed methods are still applicative.

i,max = j∈N pji max ∀β,γ,h {t βγh } di,min = min ∀l { M i +O i f l } di,max = max ∀l { M i +O i f l } tr i ∈ [ tr i,min , tr i,max ] di ∈ [ di,min, di,max] dj = max ∀i,p ij =1 {di + tr i + di} H = max ∀i {di} E t k,min = N min ∀β,γ,h {r βγkh } E c k,min = HP 0 k + i∈N [min ∀l { M i f l (P s l + P d l -P 0 k )}] E l k = ηE h k
Tuned parameters Min/Max/Step N : 25/50/5 η : 0.8/0.9/0.1

Firstly, we compare the system performance (i.e., the system energy consumption and the system QoS) with the proposed task mapping method (i.e., PP) and other task mapping methods [START_REF] Tian | Cross-layer collaborative in-network processing in multihop wireless sensor networks[END_REF], [START_REF] Pathak | Energy-efficient task mapping for datadriven sensor network macroprogramming[END_REF], [START_REF] Mo | Mapping imprecise computation tasks on cyber-physical systems[END_REF]. Secondly, we explore the algorithm performance (i.e., the computation time and the system QoS) of the proposed OTM and HTM methods with: Branch and Bound (B&B) [START_REF] Boyd | Branch and bound methods[END_REF] and Branch and Cut (B&C) [START_REF] Albert | Solving mixed integer linear programs using branch and cut algorithm[END_REF], which are known to provide the optimal solution for the MILP problem. The simulations are performed on a PC with quad-core 2.5 GHz Intel i7 processor and 16 GB RAM, and the algorithms are implemented in Matlab 2016a with CPLEX solver.

Let QoS-WDM, QoS-NDM, EE-ND and EE-WD denote the methods proposed in this paper, [START_REF] Mo | Mapping imprecise computation tasks on cyber-physical systems[END_REF], [START_REF] Pathak | Energy-efficient task mapping for datadriven sensor network macroprogramming[END_REF] and [START_REF] Tian | Cross-layer collaborative in-network processing in multihop wireless sensor networks[END_REF], respectively. Specifically, QoS-WDM and QoS-NDM are QoS-aware task mapping methods, while EE-ND and EE-WD are energy-aware task mapping methods. Compared with QoS-NDM, DVFS and multi-path routing are considered in our QoS-WDM method. In addition, compared with EE-ND, DVFS is considered in EE-WD. Fig. 4 and Fig. 5 show the system performance under these task mapping methods. From Fig. 4, we observe that the QoS increases with the values of N and η in QoS-WDM and QoS-NDM, while the QoS is always equal to 0 in EE-ND and EE-WD. This is because the aims of QoS-WDM and QoS-NDM are to maximize QoS under energy and real-time constraints, while the aims of EE-ND and EE-WD are to minimize energy under real-time constraints. Therefore, with values of N and η increasing, more optional subtasks are executed in QoS-WDM and QoS-NDM, thus, a higher QoS is achieved. The cycles of mandatory subtasks are fixed and they must be always executed. However, the less optional subtask cycles are executed, the less is the energy consumed to execute the tasks. Thus, the execution cycles of the optional subtasks are 0 in EE-ND and EE-WD. Fig. 4 also shows that the QoS achieved by QoS-WDM is higher than QoS-NDM. This is because DVFS and multi-path routing are used in QoS-WDM. Compared with QoS-NDM, where the frequency assignment and the routing path selection are fixed, QoS-WDM is able to find better decisions for frequency assignment and routing path selection, increasing QoS under time and energy constraints. The consumed energy under the different task mapping methods is compared in Fig. 5. Although QoS-WDM and QoS-NDM require more energy than EE-ND and EE-WD, the consumed energy is always smaller than the supplied energy, as the energy related constraints ( 12)-( 14) must be satisfied. From Fig. 4 and Fig. 5, we observe that QoS-aware task mapping method provides a better usage of system resources to enhance QoS. Fig. 5 also shows that the energy efficiency factor η doesn't affect the energy consumption of EE-ND and EE-WD, under the same number of tasks N , since EE-ND and EE-WD aim to minimize the energy consumption. In addition, EE-WD consumed less energy than EE-ND. Typically, for the energyaware task mapping problem, as long as the constraints (e.g., deadline and energy) allow it, methods applying DVFS is able to achieve a better energy efficiency compared to methods without DVFS. In order to further evaluate the behavior of the proposed QoSaware task mapping QoS-WDM, we compare the schedulability of QoS-WDM and QoS-NDM, using the problem feasible ratio as a metric. We set the number of tasks N = 30 and change the value of energy efficiency factor η. Under a given η, we independently construct and solve the task mapping problems 30 times (i.e., N f = 30) for QoS-WDM and QoS-NDM. For each experiment the parameters are randomly generated and same values are used for QoS-WDM and QoS-NDM. Let N w and N n denote the times that the task mapping problems in QoS-WDM and QoS-NDM are feasible, respectively. Therefore, the problem feasible ratios for QoS-WDM and QoS-NDM are defined as N w /N f and N n /N f , respectively. From Fig. 6(a), we observe that with the value of η increasing, the problem feasible ratios of QoS-WDM and QoS-NDM increase as well. This is because with a higher energy supply, a processor can use a higher frequency to execute faster an assigned task. Therefore, the time and the energy related constraints are easier to be satisfied. Fig. 6(a) shows that the problem feasible ratio of QoS-WDM is always higher than that of QoS-NDM. This is due to fact that, by considering DVFS and multi-path routing in QoS-WDM, the explored design space is larger, allowing QoS-WDM to find solutions, even for cases where QoS-NDM is not able to. Similar are the results when we compare the methods with and without DVFS for the energy-aware task mapping problem, as shown in Fig. 6(b). The QoS and the computation time of using OTM, HTM, B&B and B&C to solve task mapping problem PP1 under different number of tasks N and energy efficiency factor η are compared in Fig. 7 and Fig. 8. Fig. 7 shows that the solutions found by OTM, B&B and B&C are the same. This means that OTM is able to find the optimal solution, verifying our analysis about the convergence of OTM. In addition, the achieved QoS increases with the value of η, since more optional subtasks are executed. On the other hand, since HTM just provides a feasible solution, OTM achieves higher QoS than HTM. As shown in Fig. 8(a), with the value of N increasing, more variables and constraints are involved into the problem, thus, the algorithm computation time increases. However, compared with OTM, HTM has a negligible computation time, since HTM only needs to solve two polynomial-time problems in sequence. In addition, OTM has a shorter computation time than that of B&B and B&C. B&C combines the benefits of B&B and Gomory cutting scheme and can better balance optimality, efficiency and stability. Usually, B&C has a faster convergence speed than that of B&B [START_REF] Albert | Solving mixed integer linear programs using branch and cut algorithm[END_REF]. Note that the computational complexity of an optimization problem is highly related to the number of variables and constraints. Solving smaller problems with less variables and constraints (i.e., MP and SP) iteratively is more efficient than solving a single large problem [START_REF] Randazzo | A comparison of optimal methods for local access uncapacitated network design[END_REF]. Fig. 8(b) shows that the influence of the energy efficiency factor η on the computation time of OTM is limited, since η doesn't change the problem size (i.e., the number of variables and constraints).

VII. CONCLUSION

In this paper, we study the IC-tasks mapping problem for the networked system. We consider dependent IC-tasks executed on the wireless DVFS-enabled nodes with real-time and energy constraints. The design objective is to assign the IC-tasks to the nodes and adjust the start and the end time of the tasks so as to maximize QoS without violating the real-time and the energy constraints. By introducing DVFS and multi-path routing, we are able to achieve a better trade-off between realtime performance and energy efficiency. We first develop a MINLP model to describe this task mapping problem. Then, we propose a MILP description of this model without performance degradation. Through the problem transformation, the problem structure can be simplified, thus, the optimal solution is easier to find. This problem is optimally solved by the proposed OTM algorithm. A novel algorithm, HTM, is proposed to reduce the computation time. Our numerical results show that OTM is guaranteed to converge to the optimal solution, while HTM is able to find a feasible solution within a negligible computation time compared with OTM. Moreover, the proposed QoS-aware task mapping strategy outperforms other task mapping strategies in term of QoS-enhancing and energy-utilizing.

Fig. 1 .

 1 Fig. 1. HVAC system.

Fig. 2 .

 2 Fig. 2. (a) Energy-oriented data path. (b) Time-oriented data path.

Fig. 3 .

 3 Fig. 3. Task dependency.

9 Fig. 4 .

 94 Fig. 4. QoS under different task mapping methods with η and N varying.

9 Fig. 5 .

 95 Fig. 5. Energy consumption of the nodes under different task mapping methods with η and N varying.

  Energy-aware task mapping.

Fig. 6 .

 6 Fig. 6. Problem feasible ratio under different task mapping methods with η varying.

9 Fig. 7 .

 97 Fig. 7. QoS under different algorithms with η and N varying.

  Computation time of OTM.

Fig. 8 .

 8 Fig. 8. Computation time under different algorithms with η and N varying.

ACKNOWLEDGMENT

  This work was supported in part by the Open Research Project of the State Key Laboratory of Industrial Control Technology, Zhejiang University, China under Grant ICT20058, in part by the Fundamental Research Funds for the Central Universities, China under Grant 2242020R10059, in part by the Outstanding Young Scholarship of Jiangsu Province under Grant BK20180012, in part by the Southeast University "Zhongying Young Scholars" and "Zhishan Scholars" Projects.

  Task τi characteristics Mi, Oi ∈ [4 × 10 7 , 6 × 10 8 ] Objective function i∈N gi(oi) = i∈N oi Constraints tr i,min = j∈N pji min ∀β,γ,h {t βγh } tr

			TABLE III
			SIMULATION PARAMETERS
		Processor θ k characteristics
	v l (V)	0.65	0.7	0.75	0.8	0.85
	f l (GHz)	1.01	1.26	1.53	1.81	2.10
	P d l (mW) 184.9 P s 246 l (mW) P s 0 (µW)	266.7 290.1	370.4 498.9 340.3 397.6 80	655.5 462.7
	M = 25		Ms = 10	Ma = 5

√ √ √ √

APPENDIX A PROOF OF LEMMA 3.1 Proof: The inequalities g ≤ b 1 and g ≤ b 2 ensure that g = 0 if either b 1 = 0 or b 2 = 0. On the other hand, the inequality g ≥ b 1 +b 2 -1 guarantees that g = 1 if both variables b 1 and b 2 are set to 1. APPENDIX B PROOF OF LEMMA 3.2 Proof: Since y = bx and -s 1 ≤ x ≤ s 2 , we have -bs 1 ≤ y ≤ bs 2 . And further, we obtain (b -1)(x + s 1 ) ≤ 0 and (b -1)(x -s 2 ) ≥ 0 due to -s 1 ≤ x ≤ s 2 and b ∈ {0, 1}. Therefore, we have y+bs 1 -x-s 1 ≤ 0 and y-bs 2 -x+s 2 ≥ 0. S 1 → S 2 holds.

If b = 0, since -bs 1 ≤ y ≤ bs 2 , y + bs 1 -x -s 1 ≤ 0 and y -bs 2 -x + s 2 ≥ 0, we have y = 0 and -s 1 ≤ x ≤ s 2 . On the other hand, if b = 1, we have -s 1 ≤ y = x ≤ s 2 . S 2 → S 1 holds.

APPENDIX C PROOF OF LEMMA 4.1

Proof: Note that the MP can be reformulated as follows

It can be solved by only considering the binary variables x.

Comparing (41) with the following problem

we have Φ(k) ≥ Φ (i). Without loss of generality, we assume that

where x * is the optimal value of x, we get Φ l (k) = Φ(k) is a lower bound of Φ * and Φ has the same physical meaning as Φ.

On the other hand, according to (34), we have

In addition, due to the strong duality between the SP and its dual problem, we get

(43) Based on (42) and (43), Φ u (k) is an upper bound of Φ * . APPENDIX D PROOF OF LEMMA 4.2 Proof: Note that the MP is a minimization problem, and the non-optimal values of binary variables x * have been excluded by the constraints in the set C 1 . In addition, with iteration number k increasing, more constraints are added into the MP. Therefore, the feasible region of the MP will shrink. Accordingly, Φ l (k + 1) = Φ(k + 1) is larger than the previous lower bounds {Φ l (0), . . . , Φ l (k)}. On the other hand, based on (34), Φ u (k + 1) is smaller than the previous upper bounds {Φ u (0), . . . , Φ u (k)}.

APPENDIX E PROOF OF THEOREM 4.1

Proof: If (31) has a bounded solution, the SP is feasible under the given MP solution x(k). However, since x(k) is a non-optimal solution, we have

On the other hand, if [START_REF] Ramkumar | Temperature aware task sequencing and voltage scaling[END_REF] has an unbounded solution, the SP is infeasible under the given MP solution x(k). However, this problem is feasible if the positive variables ξ = [ξ i ] v×1 are introduced to relax the constraints. In order to minimize ξ, we construct the following problem:

Since (44) is a LP problem, we solve its dual problem:

where μ = [μ i ] v×1 are the dual variables. Let ξ(k) and μ(k) denote the solutions of ( 44) and (45), respectively. Since the relaxation variables with respect to infeasible constraints are non-zero, we have 1 T ξ(k) > 0. On this basis, due to the strong duality between (44) and (45), we get

APPENDIX F PROOF OF THEOREM 4.2

Proof: Note that the gap between the lower and upper bounds gradually reduces. In addition, the dimension of binary variables x is finite, and the non-optimal and infeasible values are excluded. Therefore, the solution (x(k), y(k)) converges to optimal one (x * , y * ) within a finite number iterations.


domain: info phys spi

text:
Laurane Charrier

Pierre Godet

Clément Rambour

Flora Weissgerber

Simon Erdmann

Elise Colin Koeniguer

Analysis of dense coregistration methods applied to optical and SAR time-series for ice flow estimations

Keywords: SAR images, optical images, coregistration, optical flow, deep learning, dense, multi-modal, ice flow, glaciers

I. INTRODUCTION

A common critical step in glaciology is the automated computation of velocity maps of glaciers from remote sensing images. Velocity maps are necessary to precisely monitor ice dynamics, to infer sub-glacial processes and/or ocean forcing and to derive other products such as mass-balance or strain rates when the amount of data is sufficient.

Operational methods produce velocity maps with a lower spatial resolution than the spatial resolution of the images [START_REF] Fahnestock | Rapid large-area mapping of ice flow using landsat 8[END_REF], [START_REF] Gardner | Increased west antarctic and unchanged east antarctic ice discharge over the last 7 years[END_REF], [START_REF] Millan | Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach[END_REF]. Indeed, most of the existing coregistration methods in glaciology rely on the correlation estimation between two images. Usually, these methods are time-consuming, so the estimation is limited on a sub-grid of pixels. The sampling, i.e., the spatial frequency at which the algorithm estimates the offset tracking, is higher than the pixel size [START_REF] Millan | Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach[END_REF] [START_REF] Fahnestock | Rapid large-area mapping of ice flow using landsat 8[END_REF]. However, a low resolution can be a limit to analyze spatial discontinuities of ice speed [START_REF] William H Armstrong | Spatial patterns of summer speedup on south central alaska glaciers[END_REF] or to derive other products such as strain rates [START_REF] Altena | Weekly glacier flow estimation from dense satellite time series using adapted optical flow technology[END_REF]. Here, we propose to analyze results obtained with dense algorithms that estimate a displacement at each pixel. Such displacement is expected to increase the spatial resolution of the velocity time-series.

Besides, operational methods only tackle SAR or optical images [START_REF] Fahnestock | Rapid large-area mapping of ice flow using landsat 8[END_REF], [START_REF] Gardner | Increased west antarctic and unchanged east antarctic ice discharge over the last 7 years[END_REF], [START_REF] Millan | Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach[END_REF]. In this paper, we explore multi-modal algorithms that can deal with both optical and SAR images. The complementary between SAR and optical images is expected to increase the temporal resolution of the velocity time-series.

The first considered method is the coregistration algorithm GeFolki [START_REF] Brigot | Adaptation and evaluation of an optical flow method applied to coregistration of forest remote sensing images[END_REF], an optical flow method that has shown high accuracy, high resolution, and high robustness on urban [START_REF] Plyer | A new coregistration algorithm for recent applications on urban sar images[END_REF] and forest [START_REF] Brigot | Adaptation and evaluation of an optical flow method applied to coregistration of forest remote sensing images[END_REF] remote sensing images. Optical flow methods have already shown relevance and limits to compute glacier velocity on optical images [START_REF] Altena | Weekly glacier flow estimation from dense satellite time series using adapted optical flow technology[END_REF], [START_REF] Vogel | Optical flow for glacier motion estimation[END_REF]. Moreover, a study on Antarctica images demonstrated the ability of GeFolki to deal with glaciers SAR images [START_REF] Erdmann | Glaciers velocity maps generation through optical flow method[END_REF]. Therefore, GeFolki appears to be an interesting method to explore.

The second considered method, more recent, is a deeplearning-based dense coregistration algorithm dedicated to remote sensing images [START_REF] Le | A coregistration tool for optical and sar images by deep learning methods[END_REF]. Trained on SAR and optical images obtained both with airborne and spaceborne sensors, this algorithm called PWC-Net-multimodal shown better results than GeFolki on SAR, optical, and SAR/optical image coregistration on continental areas. PWC-Net-multimodal has been trained with synthetic displacement fields computed from a digital elevation model, to correct offsets due to relief. Therefore, this algorithm is initially not intended to operate on terrain or glacier movements. However, to re-train this type of algorithm on glacier flows would be a lengthy process. The generation of realistic synthetic glacier flows would require a full study that is beyond the scope of this paper. Therefore, in the first place, we propose to test this network on glaciers without another training step. This study aims to carry out the analysis of the two different methods, in the area of the Totten glacier, a stable and major outlet, where GeFolki has already shown convincing results on radar images [START_REF] Erdmann | Glaciers velocity maps generation through optical flow method[END_REF]. We consider dense coregistration on optical and radar images, to assess the potential of the different methods available in terms of producing velocity series estimates. Several criteria are evaluated:

• the ability to produce estimates using SAR/SAR pairs; optical/optical, or SAR/optical. • the accuracy on the estimated velocity • the spatial resolution of the velocity fields • the robustness of the estimates.

II. DATASET AND PREPROCESSING

We have created different Sentinel-1 and Sentinel-2 datasets on Totten glacier (67°S 116°E) in Antarctica. Totten glacier is one of the major contributor of mass loss in East-Antarctica. It drains a sector 537.900 km 2 in size, mostly grounded below sea level. [START_REF] Gardner | Increased west antarctic and unchanged east antarctic ice discharge over the last 7 years[END_REF], [START_REF] Li | Ice flow dynamics and mass loss of totten glacier, east antarctica, from 1989 to 2015[END_REF] reveal that Totten glacier increased in velocity between 2001 and 2007 but remained stable since, likely due to a prograde bed [START_REF] Morlighem | Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the antarctic ice sheet[END_REF]. This stability allows us to validate our measured velocities with previous works. Because the methods have never been applied to glaciers, Totten glacier would be a relevant first case to study.

The creation of a dataset was subject to various constraints related to the footprints and dates of acquisitions. Figure 1 shows the areas related to the two footprints used. The Sentinel-1 EW (about 40 m per pixel) data have large footprints that allow us to see the entire glacier in motion. The Sentinel-1 IW and Sentinel-2 data are better resolved (about 10 m per pixel), but the footprints do not cover such a large area of the glacier. Furthermore, the dates of acquisitions are not necessarily overlapping. Sentinel-1 EW images are mainly available for the year 2017 whereas Sentinel-2 images without clouds were mainly acquired in 2018. As for Sentinel-1 IW images, it covers the years 2017 and 2018. Taking all these constraints into account, we created two sets of data:

• The first set is composed of a dozen of Sentinel-1 images covering a large area, in EW mode, for the year 2017 (cf.

Figure 2). • The second set is composed of Sentinel-1 images, which cover a much more restricted area, in IW mode, over 18 months from June 2017 to the end of 2018, and 4 Sentinel-2 images acquired in 2018 (cf. Figure 7).

Images were downloaded using the platform Google Earth Engine, in Level 1 GRD for radar images and Level 1C for optical images. The Terrain Correction step offered by SNAP is already applied to the images.

On the first stack, we make a comparison of the two dense co-registration algorithms GeFolki and PWC-Net-multimodal, while on the second stack, we investigate mainly GeFolki and its capacity to handle optical images.

A colored composition was produced using the REACTIV algorithm [START_REF] Colin-Koeniguer | Colored visualization of multitemporal SAR data for change detection: issues and methods[END_REF] to quickly visualize the second stack of IW images. Figure 3 allows us to view all the images in the stack. A pixel where the signal is not stable in all the images trough the time is saturated. Each color is representative of one specific date.

III. METHODS

In this section, we briefly present the two evaluated algorithms to understand their respective advantages and drawbacks for glacier displacements estimation.

The first considered method, GeFolki, is an optical flow algorithm that has already shown accuracy as good as correlation methods [START_REF] Brigot | Adaptation and evaluation of an optical flow method applied to coregistration of forest remote sensing images[END_REF]. The second considered algorithm, PWC-Netmultimodal, is a more modern CNN method. Both are dense and multi-modal coregistration methods.

A. GeFolki

GeFolki is a dense flow estimation algorithm based on optical flow equations. Its principles are widely described in [START_REF] Brigot | Adaptation and evaluation of an optical flow method applied to coregistration of forest remote sensing images[END_REF], [START_REF] Plyer | A new coregistration algorithm for recent applications on urban sar images[END_REF], but we emphasize here some of them.

GeFolki is based on the Lucas-Kanade paradigm [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF]. It belongs to local or window-based approaches of optical flow: the flow estimated in one pixel is defined as the minimizer of a criterion computed over a local window centered on this pixel. Minimization is done by an iterative Gauss-Newton strategy based on first-order Taylor expansion of the image intensity around a previous guess of optical flow. The first-order Taylor expansion assumes a constant flow over the local window used, therefore spatial discontinuities inside this window are likely to be smoothed. GeFolki is also multi-resolution, i.e., uses a pyramid of images to compute the optical flow at varying scales following a coarse-to-fine strategy. The minimization at each level of resolution, called the level of the pyramid, is initialized with the flow from the previous level. Moreover, at each pyramid level, the window's radius can also iteratively vary. Finally, a Rank transform is applied to the images before minimization. Each pixel is replaced by the number of neighboring pixels with an intensity lower than its intensity. This specificity is of a broad interest in glaciers applications. It allows us to be more robust to radiometric changes. Indeed, the radiometric changes often cause failure in glaciers applications: [START_REF] Altena | Weekly glacier flow estimation from dense satellite time series using adapted optical flow technology[END_REF], [START_REF] Vogel | Optical flow for glacier motion estimation[END_REF] show that, on optical images, irradiance variations between autumn and spring at high latitudes or reflection intensity changes due to different snow "ages" were major sources of failure for optical flow methods.

In practical terms, GeFolki has several parameters to tune: the number of pyramid levels, the radius of the windows, the size of the rank filter, and the number of iterations. The number of pyramid levels should be set according to the expected maximum of flow magnitude. Regarding the radius, the higher it is, the more robust is the algorithm, but the smoother is the result. When increasing the radius, the precision increases, but the resolution, i.e. the spatial frequency at which the precision can be maintained, decreases. We will discuss later a strategy to adjust precisely GeFolki parameters according to each type of image.

In previous studies, GeFolki has been tested on highresolution radar/radar or radar/optical images of urban areas [START_REF] Plyer | A new coregistration algorithm for recent applications on urban sar images[END_REF] and optical/radar images of vegetated areas [START_REF] Brigot | Adaptation and evaluation of an optical flow method applied to coregistration of forest remote sensing images[END_REF]. It has proven to be fast in comparison to other image correlation strategies. [START_REF] Brigot | Adaptation and evaluation of an optical flow method applied to coregistration of forest remote sensing images[END_REF] obtain with GeFolki a flow for a 5,000 x 5,000 pixel image on a Intel(R) Core(TM) i7-2620M @ 2.80GHz with a GPU implementation in hundreds of milleseconds.

B. PWC-Net-multimodal

The PWC-Net-multimodal algorithm is based on the PWC-Net architecture, a state-of-the-art Convolutional Neural Network (CNN) for optical flow estimation in computer vision [START_REF] Sun | PWC-Net: CNNs for optical flow using pyramid, warping, and cost volume[END_REF].

The original PWC-Net uses a siamese encoder for the two input images. On the opposite, PWC-Net-multimodal has two different encoders: one specialized for optical images, the other for radar images. It is then able to deal with heterogeneous image pairs. PWC-Net-multimodal is trained on a remote sensing image dataset composed of radar and optical images from different resolutions and frequencies (low-resolution Sentinel-1 and Sentinel-2 images, and on higher-resolution aerial images from American UAVSAR and NAIP programs), covering continental areas. The training step mixes the coregistration of optical/optical, radar/radar, and radar/optical images. The generalization ability of the algorithm has proved to be very good on X band images with higher spatial resolution. However, the network has been trained with synthetic flow indexed on the SRTM DEM. The goal was to fix errors in the coregistration between images, assumed to be caused by errors in projection over the topography. Therefore, it has never been trained neither on ice-sheets and glaciers nor when the flow is decorrelated from the relief.

IV. ANALYSIS

A. Comparison of GeFolki and PWC-Net-multimodal on Sentinel-1 EW images

First, we study the performances of GeFolki and PWC-Net-multimodal in controlled conditions. We try to find an emulated displacement to assess the ability of both algorithms to retrieve a known flow. Indeed no ground truth measurements are available near Totten glacier. Therefore, we apply a flow from the literature to a Sentinel-1 EW master image to obtain an emulated slave image. Here, we use a master image from the 2 May 2017 and a time-average velocity map called ITS LIVE published by [START_REF] Gardner | Increased west antarctic and unchanged east antarctic ice discharge over the last 7 years[END_REF] with a resolution of 120 m (roughly 3 pixels on Sentinel-1 EW images). The covering area is represented in Figure 9. The velocity map is reprojected and resampled, and we convert its values to obtain a displacement expressed in pixels in the geometry of the master image, for a period equivalent to 6 months.

On the one hand, performances of GeFolki depend on the chosen parameters, especially the radius and number of pyramid levels. Therefore, we pick the parameters which minimize the criterion MSE (Mean Square Error) and MAE (Mean Absolute Error) between the reference flow and the resulted flow. We obtain the highest MSE and MEA when we set the number of pyramid L such as 2 L < D, where D is the maximum of expected flow magnitude. For Sentinel-1 EW images, the best parameters appear to be: radius= [64,32,16,8], levels=4, rank=4, iteration=2. It allows us to reach a MAE of 6 10 -2 pixels and MSE of 9 10 -2 pixels for the flow magnitude. To quantify errors in angles which are circular variables varying modulo 2π , the use of circular statistics is more suitable. Thus the mean of the circular difference between flow directions is 6 10 -4 rad. This gives an idea of the accuracy of GeFolki. The results show that GeFolki can retrieve an emulated flow in most parts of the image (cf. Figure 4). The most difficult areas correspond to the border of the ice stream and the ice front where a break in velocity occurs. This is likely because GeFolki is locally constrained. On the other hand, PWC-Net-multimodal seems to retrieve the known flow with less accuracy, in particular in the Dome area and near the ice front (cf. Figure 5). The MAE is of 1.49 pixels and MSE of 10.41 pixels for the flow magnitude. The mean of the circular difference between flow directions is of 2.6 10 -2 rad. Then, we measure a real displacement between two Sentinel-1 EW images from the 2 May 2017 and the 3 November 2017 (pair 1) using the previous GeFolki parameters. The results are shown in Figure 6. According to [START_REF] Li | Ice flow dynamics and mass loss of totten glacier, east antarctica, from 1989 to 2015[END_REF], we expect the flow magnitude to be around 1800 m/yr near the ice front and 800 m/yr in the ice stream. We convert the estimated displacement in pixel to a displacement in m/yr. For this analysis, we assume that geometric distortion and errors due to ortho-projection are negligible. We find a flow magnitude from 1600 to 2000 m/yr near the ice front, and from 500 to 1000 m/yr in the ice stream. The direction, which is South-West / North-East in most of the image, is also a consistent result: the ice stream flows toward the ocean. However, we notice that the direction found in the Western Southern part of the image is the opposite of the expected one. We try to compute a displacement on other pairs: 1 July 2017 -3 November 2017, 1 July 2017 -4 August 2017 and 2 May 2017 -4 August 2017 and find the same behavior. This point needs further investigation. Furthermore, several criteria can be used to evaluate an estimated displacement. Here, we propose to have a look at the back-forward flow, also called residual. This is the difference between the computed displacement from the image 1 to 2 and its opposite, the displacement from image 2 to 1. For the pair 1, the residual of the x-component of the flow has a median value of 0.14 pixels and a standard deviation of 0.35 pixels. Regarding the y-component of the flow, the median value of the residual is 0.14 pixels, and its standard deviation is 0.25 pixels. By doing the same process for all the four evaluated pairs of images, we find, on average, a residual with a median value of 0.14 pixels for both the x-component and the ycomponent. That gives an idea of the precision of GeFolki for these pairs of images. This criterion can also help to reject the outliers, i.e. pixels where the residue is above a given threshold. Finally, we compare the performances of GeFolki and PWC-Net-multimodal on the same pair of real images (cf. Figure 6). PWC-Net-multimodal finds a high displacement (between 1900 and 2300 m/yr) near the ice front, which is quite consistent with GeFolki's results. However, this high displacement seems to be too concentrated to be realistic. We also notice a huge displacement at the bottom of the image. Since this latest result doesn't have any glaciology meaning, it could be due to an over-learning of relief elements by the algorithm.

To conclude, GeFolki has proven its ability to retrieve a given glacier flow. Even if the reference flow used has a resolution of 120 m and corresponds to a time-average, we can assume that it remains closed to a real glacier displacement. The emulation of such displacement could, therefore, be a way to choose GeFolki parameters. These parameters need to be changed depending on the type of image (optical or SAR) and its resolution. The application of GeFolki on a real pair of images leads to a coherent flow magnitude in the whole image and a coherent direction on most of the ice-sheet.

The results give a first idea of the precision and accuracy of the algorithm for ice flow estimations. However, PWC-Netmultimodal does not seem to be suitable to compute a glacier displacement without another dedicated training. It is likely because the training dataset uses a SRTM digital elevation model to emulate the deformation. This choice was relevant to align two images with no natural source of displacement, but need to be adapted with a glacier flow database for glaciology purpose. In this case, we think that PWC-Net-multimodal could even be more robust than GeFolki, as it was already the case in our previous experiments [START_REF] Le | A coregistration tool for optical and sar images by deep learning methods[END_REF]. Both algorithms still appear to be high-speed methods to retrieve a flow and allow a dense coregistration.

B. Comparaison of GeFolki on optical and SAR images

The GeFolki algorithm can handle heterogeneous images. Therefore, we compare the results of GeFolki on Sentinel-2, Sentinel-1 IW, and Sentinel-1 EW images for the same period. The studied footprint appears in yellow in Figure 1. One of the corresponding Sentinel-2 images is shown in Figure 7.

For Sentinel-2, we use the blue channel (10 m of spatial resolution and 12 bit of radiometric resolution). For Sentinel-1 EW and IW, the channel corresponds to the polarimetry HH. We choose GeFolki parameters using the methodology previously presented in section A. For Sentinel-1 IW, the spatial resolution is higher than the one of Sentinel-1 EW, therefore we should increase the levels of pyramids (cf. Table I). For Sentinel-2, the pair of images corresponds only to a small baseline, so we pick the parameters listed in Table I.

Results on Sentinel-1 IW are similar to results on Sentinel-1 EW: we measure roughly the same flow magnitudes and directions in the same locations (cf. Figure 8). Since the resolution of the Sentinel-1 IW image is higher, we notice more details on the resulted velocity map. Results on Sentinel-2 and Sentinel-1 IW are also very correlated on the ice-sheet parts (Bottom Left part of the image) (cf. Figure 9). We have also tried to coregistrate a pair of radar/optical images. The results are promising but still need to be improved. Indeed, the flow direction is consistent but the flow magnitude is not. This is probably due to a difference in the acquisition geometry. When considering qualitatively the superposition of close images, it appears that the shadows are not projected in the same way. Therefore, GeFolki tries to compensate for the misalignment of the shadow projections instead of a real feature.

V. CONCLUSION

This paper aimed to investigate the capabilities of current coregistration methods designed to handle various images and produce dense displacements for the task of glacier displacements estimation. The goal is to explore the possibility of increasing both the spatial and temporal resolution of velocity maps. We have chosen and evaluated two algorithms: GeFolki and PWC-Net-multimodal. GeFolki has already been successfully applied to several other cases, while PWC-Netmultimodal is a more recent algorithm that gave remarkable results on specific cases but still requires further evaluations.

This study confirmed the ability of GeFolki to retrieve an ice flow using optical and SAR images. Five main advantages appear to us compared to other classical coregistration methods such as image matching methods. First, GeFolki already proved to be a fast method in comparison with other correlation methods [START_REF] Brigot | Adaptation and evaluation of an optical flow method applied to coregistration of forest remote sensing images[END_REF]. It may quickly provide a displacement at each pixel without, for example, the need for previous sparse research as developed in [START_REF] Gardner | Increased west antarctic and unchanged east antarctic ice discharge over the last 7 years[END_REF]. Secondly, the flow computation relies on several window sizes, for the same pixel. It offers the possibility to take into account different levels of details and to be more robust. Thirdly, GeFolki can directly retrieve a subpixel offset, contrary to cross-correlation methods which need to oversample the images and/or the correlation functions. For instance, the method proposed in [START_REF] Michel | Flow of glaciar moreno, argentina, from repeat-pass shuttle imaging radar images: comparison of the phase correlation method with radar interferometry[END_REF] uses an oversampling of the image in the Fourier domain through Shanon interpolation. In [START_REF] Joughin | Ice-sheet velocity mapping: a combined interferometric and speckle-tracking approach[END_REF], the authors propose to oversample the correlation function by a factor of 10. Then, to deal with contrast variations, GeFolki includes a rank filter processing at each level of the pyramid, instead of previously filtering the image with a lowpass filter [START_REF] Fahnestock | Rapid large-area mapping of ice flow using landsat 8[END_REF] or a Wallis operator [START_REF] Gardner | Increased west antarctic and unchanged east antarctic ice discharge over the last 7 years[END_REF]. Finally, GeFolki can even be applied without interferometric conditions. Indeed, it is not sensitive to loss of coherence which is common on SAR images of glaciers and ice-sheets. PWC-Net-multimodal is a much more recent algorithm than GeFolki. Even if it has never been trained on real displacement cases, this CNN can retrieve an ice flow that is consistent with GeFolki results in areas where big displacements occur, for instance on the ice fronts. However, several artifacts are spotted in areas of weak displacements outside of the ice stream. Therefore, the results still have to be improved. For that, a good solution could be to design a new training process of CNN. As PWC-Net-multimodal has never been trained on glaciers or ice-sheets, it appears to us that this algorithm could provide very promising results with a dedicated training dataset. A strategy could be to use glacier flows estimated by other algorithms. However, the training step would be impacted by the accuracy and resolution of these velocity maps. Besides, data augmentation techniques could improve performances at a lower cost.

To go further, it would be interesting to apply GeFolki and PWC-Net-multimodal to glaciers with more complexity such as mountain glaciers like alpine glaciers. Furthermore, GPS data from GLACIOCLIM observatory are available in the Alps (https://glacioclim.osug.fr/). It would allow us to give another validation and a measure of the accuracy of our glacier surface velocities. Last but not least, it will be very relevant to compare the performance of GeFolki to other algorithms such us auto RIFT [ [START_REF] Gardner | Increased west antarctic and unchanged east antarctic ice discharge over the last 7 years[END_REF]], Pycorr [START_REF] Fahnestock | Rapid large-area mapping of ice flow using landsat 8[END_REF] or the modified version of ampcor used by [START_REF] Millan | Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach[END_REF].

Fig. 1 .

 1 Fig. 1. Totten Glacier footprints. The red footprint corresponds to the stack of Sentinel-1 EW images and the yellow footprint to the stack of Sentinel-1 IW and Sentinel-2 images.

Fig. 2 .

 2 Fig. 2. Sentinel-1 EW (40m) image from the 3 November 2017 on Totten Glacier. The level of the product is the Level1 GRD.

Fig. 3 .

 3 Fig.3. Colored composition of the Sentinel-1 IW images using REACTIV algorithm[START_REF] Colin-Koeniguer | Colored visualization of multitemporal SAR data for change detection: issues and methods[END_REF]

Fig. 4 .

 4 Fig. 4. Difference between the flow from ITS LIVE and the flow computed by GeFolki on an emulated pair of Sentinel-1 EW images where the master image is from the 2 May 2017. GeFolki parameters are radius= [64,32,16,8], levels=4, rank=4, iteration=2. The difference in magnitude is displayed on the left and the circular difference in direction on the right. Histograms are displayed below the corresponding difference.

Fig. 5 .

 5 Fig. 5. Difference between the ITS LIVE flow and the flow computed using PWC-Net-multimodal on an emulated pair of images. The master image is from the 2 May 2017. The difference in magnitude is displayed on the left and the circular difference in direction on the right. Histograms are displayed below the corresponding difference.

Fig. 6 .

 6 Fig. 6. Flow computation on a real pair of Sentinel-1 images. The master image is from the 2 May 2017 and the slave image from the 3 November 2017. Results of GeFolki (left), and of PWC-Net-multimodal (right).

Fig. 7 .

 7 Fig. 7. Sentinel 2 image (10 m) from the 31 January 2018 on Totten Glacier. The level of the product is the Level 1C. The values correspond to Top Of the Atmosphere (TOA) reflectance.

Fig. 8 .

 8 Fig. 8. Comparison between a flow estimation with GeFolki on real Sentinel-1 EW images (left) and Sentinel-1 IW images (right) covering roughly the same period. The GeFolki parameters used are listed in TableI. For Sentinel-1 IW, the master image is from the 7 July 2017 and the slave image from the 5 October 2017. For Sentinel-1 EW, the slave image is from the 3 July 2017 and the slave image from the 3 October 2017.

Fig. 9 .

 9 Fig. 9. Comparison between a flow estimation with GeFolki on real Sentinel-2 images (left) and Sentinel-1 IW images (right) overing roughly the same period of time. The GeFolki parameters used are listed in TableI. For Sentinel-1 IW, the master image is from the 27 January 2018 and the slave image from the 8 February 2018. For Sentinel-2, the slave image is from the 31 January 2018 and the slave image from the 10 February 2018.


Aida Todri-Sanial

Toward New Era in Computing: From Devices to Applications

This talk aims to give a bottom-up overview on novel computing paradigms covering aspects from devices, interconnects, architectures to applications. Non Von-Neumann architectures such as neuromorphic and quantum computing progress and challenges are discussed including some of our current research on these topics.


Think Digital

•

  SNCF Voyages : L'IA au service de la Relation Client -Olivier Beaurepaire -Directeur de la Digitalisation de l'Innovation, services et RH chez SNCF Voyages -Didier Bourhis -Directeur Produits Digitaux chez SNCF Voyages -Joffrey Martinez -Leader de la practice Watson IA, IBM France • Géodis : Réussir dans l'e-commerce • Supply Chain : Comment un expert de la logistique souhaite devenir le partenaire privilégié du développement e-Commerce -Ashwani Nath -VP e-Commerce, Geodis -Arnaud Delande -Directeur AI Applications, IBM France • PSA : La planification à long terme, un enjeu stratégique pour répondre à la demande future du marché -Antoine Dewerdt -Directeur des Prévisions de Vente Long Terme Monde chez Groupe PSA -Frédéric Donneaud -Associate Partner, IBM France Journey to Cloud Durant une heure, avec le témoignage de nos clients et la participation des experts IBM, nous aborderons les grands enjeux des infrastructures IT et du Cloud ainsi que les solutions mises en oeuvre. • Optimiser l'IT pour anticiper et s'adapter rapidement à la demande des clients -Ouafaa El Moumouhi -Directeur IBM Public Cloud, IBM France -Stéphane Libman -Directeur des services Hybrid et Multicloud, IBM France -Mimoun Lakri -CTO IBM Public Cloud, IBM France -Lionel Mace -Senior Architecte Cloud, IBM France • Migrer et gérer les workloads critiques dans le cloud de manière intelligente -Stéphane Libman -Directeur des services Hybrid et Multicloud, IBM France -Mimoun Lakri -CTO IBM Public Cloud, IBM France -Arnaud Le Fur -Architecte, expert SAP, IBM France • Automatisation des infrastructures IT assistée par l'intelligence artificielle -Samia Benali -Architecte IBM Services, IBM France -André Deville -IBM / RedHat Hybrid Cloud Synergy Sales, IBM France • Cybersécurité : Sécurité des Clouds -Erwan Vezin -DSI, John Paulo -Aurore Ominetti -Directeur IBM Security, IBM France Think Digital Summit France 6 octobre 2020 Innovation Technologies et innovation : préparez-vous au monde de demain. IBM aide ses clients à faire face à l'émergence de nouveaux défis, chaque jour plus complexes. Dans cette session, nous parlerons technologies innovantes pour ouvrir le champ des possibles et vous préparer au mieux au monde de demain. • Informatique quantique : repousser les limites du calculable -Xavier Vasques -Directeur du Centre Technologique Européen, IBM France -Olivier Hess -Chargé des activités quantiques pour la France & IBM Q Ambassador, IBM France -Aida Todri-Sanial -Directrice de Recherche CNRS, Laboratoire d'Informatique de Robotique et de Microélectronique de Montpellier • Automatisation intelligente : l'IA au coeur des processus pour améliorer la performance opérationnelle des entreprises -Harley Davis -VP IBM France Lab & Automation -Christian de Sainte-Marie -Directeur du Centre des Etudes Avancées, IBM France Lab • IBM Research : l'innovation au service de la planète -Sheila O'Hara -Senior Architect, IBM Global Industry Solutions Center • The Mayflower AS ou la révolution du transport autonome avec l'IA et l'edge computing -Éric Aquaronne -Responsable stratégie mondiale produits matériels, IBM France -Guilhem Molines -ODM Product Architect, IBM France et orateurs des chaînes Smarter Business, Journey to Cloud et Innovation lors de interactives