Margaux Nattaf
email: nattaf@laas.fr

Christian Artigues
email: artigues@laas.fr

Pierre Lopez
email: lopez@laas.fr

Cumulative scheduling with variable task profiles and concave piecewise linear processing rate functions

Keywords: continuous scheduling, continuous resources, concave piecewise linear functions, energy constraints, energetic reasoning

Cumulative scheduling with variable task profiles and concave piecewise linear processing rate functions

Margaux Nattaf, Christian Artigues, Pierre Lopez

This paper deals with a scheduling problem involving a set of tasks and a continuously-divisible renewable resource of limited capacity shared by the tasks. We consider the case where tasks resource requirement is not fixed but is part of the problem and has to be determined. Thus, the duration of a task is not fixed either but is determined by the resource requirement function as the task is finished once it has received a necessary amount of energy. Furthermore, we consider that the total energy received by a task is not equal to the total amount of the resource used by it. Instead we have processing rate (or efficiency) functions, which translates the required resource amounts into energy.

We perform an analysis of the structural properties of the problem for concave piecewise linear processing rate functions. To our knowledge, this variant of the cumulative constraint has never been considered in the literature. We show first that the resource demand profile of a task can be restricted to a piecewise constant function with break points at the starts and ends of tasks. From these theoretical properties, we are able to compute the minimal resource consumption of a task inside an interval in O(1) and we prove that the set of the relevant intervals of polynomial size that was shown sufficient for energetic reasoning with linear functions is also sufficient in our case. Furthermore, a full characterization of relevant intervals for time-window adjustments is provided, completing the work for linear function. We also define a new propagation algorithm together with a satisfiability test, which relies on the time-table disjunctive reasoning for the cumulative constraint [START_REF] Gay | Time-Table Disjunctive Reasoning for the Cumulative Constraint[END_REF] for the first one and on a flow-based linear program for the latter one. Finally, complementarity of the algorithms is assessed via their integration in a hybrid branch-and-bound and computational experiments on small-size instances.

Problem statement, properties and context

In this section, we formally define the Continuous Energy-Constrained Scheduling Problem (CECSP). Then, we present a foundry application in details and finally, we exhibit some properties of the CECSP, which we will use throughout the paper.

Problem definition

In the CECSP, a set of tasks A = {1, . . . , n} has to be scheduled using a continuous, cumulative and renewable resource of capacity B. The resource amount that a task requires during its processing time is not fixed but instead the resource usage of a task i ∈ A is a function of time, b i (t), that must be determined in a continuous time setting (t ∈ R + ). Once the task is started and until its finishing time, its resource usage is constrained to lie between a maximum and a minimum requirement, b max i and b min i = 0, respectively 1 . In addition, when a task uses a part of the resource, it receives a certain amount of energy and we say that a task is finished when it has received an energy W i . This energy is computed using a continuous, non-decreasing, concave and piecewise linear power processing rate function [START_REF] Machowiak | Scheduling malleable tasks on parallel processors to minimize the makespan[END_REF]:

f i : [b min i , b max i ] -→ R b i (t) -→ f i (b i (t))
Furthermore, each task has a release date r i and a deadline d i , and has to be fully executed in [r i , d i ].

The CECSP consists of finding, for each task, a start time

st i ∈ [est i , lst i ], an end time et i ∈ [eet i , let i ], and its resource usage b i (t) ∈ R + , ∀t ∈ R + such that: ri ≤ sti < eti ≤ di ∀i ∈ A (1)
1 Although all results presented in this paper can be adapted when b min i = 0, for ease of notation, we assume b min i = 0.

b min i ≤ bi(t) ≤ b max i ∀i ∈ A, ∀t ∈ [sti, eti] (2)
et i st i fi(bi(t))dt = Wi ∀i ∈ A (3) i∈A bi(t) ≤ B ∀t (4)
Example 1. In the example of Fig. a: Instance data. Throughout this paper, we will use the following expression for function f i (see Figure 2):

1 3 B = 5 2 b: A solution. f3(b) =    2 • b + 1 if b ∈ [1, 2[ b + 3 if b ∈ [2, 3[ 1/4 • b + 21/4 if b ∈ [3, 5]
f i (b i (t)) =          a i1 • b i (t) + c i1 if b i (t) ∈ [b min i , x i 1 [ a i2 • b i (t) + c i2 if b i (t) ∈ [x i 1 , x i 2 [ . . . a iPi • b i (t) + c iPi if b i (t) ∈ [x i Pi-1 , b max i ]
with P i being the number of pieces of f i and x i p , p ∈ {1, . . . , P i -1} its breakpoints.

Context

The CECSP comes from an industrial problem. This problem, presented in [START_REF] Artigues | The energy scheduling problem: Industrial case study and constraint propagation techniques[END_REF], arises in a pipe-manufacturing plant and more precisely, the foundry where metal

b fi(b) b min i x i 1 x i 2 x i 3 x i 4 b max i a i1 • b i (t) + c i1 a i2 • b i (t) + c i2 a i3 • b i (t) + c i3 a i4 • b i (t) + c i4 a i5 • b i (t) + c i5
Fig. 2: Illustration of the notations for function f i with P i = 5.

is melted in induction furnaces. In this department, melting and heating use a huge amount of energy especially electricity. The expenses of the plant for electricity represent more than half of the annual energy costs. The cost of electricity depends on the total energy consumed and on penalties for power overrun in reference to a subscribed maximal power. The goal is to minimize the energy bill.

The foundry has several production lines (furnaces) and each metal operation has to be assigned to a furnace and scheduled within a time window. Furthermore, an operation has a variable duration that depends on the power given to the furnace. Thus, the electrical power of the furnaces, which can be adjusted at any time to avoid exceeding a maximum limit, can be seen as a continuous function of time to be determined. However, the function must lie within a limit (due to physical and operational considerations); thus, a minimum and a maximum power level must be satisfied for the melting operation. A melting job is composed of three sequential parts: loading, heating and unloading. The duration of loading and unloading are known but heating duration has to be determined. The heating operation can be stopped once the necessary energy has been received.

To solve this problem, the authors suppose that jobs have a piecewise constant consumption profile and define a two-step method. In the first step, scheduling of jobs on the furnaces is performed, using a constraint programming model, with fixed job durations. The task sequence and assignment are then used as data in the second step and a MILP model is used to determine their consumption profile and duration. The algorithm then returns to the first step with new jobs durations and it stops when the objective function is no more improved.

Unfortunately, the initially considered energy model was not sufficient to achieve a good energy consumption. Therefore, a new problem, which can be used to determine task starting and finishing times as well as their consumption profile was introduced in [START_REF] Artigues | The energy scheduling problem: Industrial case study and constraint propagation techniques[END_REF], the Energy Scheduling Problem (EnSP). Due to the complexity of the problem, the proposed model considers a time discretization, which can lead to suboptimal or infeasible solutions by over-constraining the problem. Furthermore, efficiency functions were not considered. In a continuous time setting but still without considering the efficiency functions, constraint propagation algorithms based on the energetic reasoning concept were proposed in [START_REF] Artigues | Energetic reasoning for energy-constrained scheduling with a continuous resource[END_REF] for the CECSP, which is the continuous version of the EnSP.

Despite the fact that processing rate functions were not considered in these papers, actual processing rate functions in scheduling problems that involve energy-consuming tasks are intrinsically continuous and non-linear. As a typical example, the non-linear function given in Fig. 3 gives the shape of the processing rate function of a fuel cell [START_REF] Ngueveu | Scheduling under a non-reversible energy source: An application of piecewise linear bounding of non-linear demand/cost functions[END_REF]. An extension of the work of [START_REF] Artigues | Energetic reasoning for energy-constrained scheduling with a continuous resource[END_REF] to linear processing rate functions as well as several solution methods was proposed in [START_REF] Nattaf | Energetic reasoning and mixed-integer linear programming for scheduling with a continuous resource and linear efficiency functions[END_REF][START_REF] Nattaf | A hybrid exact method for a scheduling problem with a continuous resource and energy constraints[END_REF].

Properties and remarks

In [START_REF] Nattaf | Energetic reasoning and mixed-integer linear programming for scheduling with a continuous resource and linear efficiency functions[END_REF], the authors use linear functions to approximate real-world ones. Despite the fact that there exist many real-world processing functions, which are concave [START_REF] Hung | Non-convex aggregative technology and optimal economic growth[END_REF][START_REF] Lewis | Algebra symposium: Optimizing fuel consumption[END_REF], approximating a function by a concave piecewise linear one is always at least as good as an approximation by a linear function (see Example 2 below). Furthermore, the authors of [START_REF] Nattaf | Energetic reasoning and mixed-integer linear programming for scheduling with a continuous resource and linear efficiency functions[END_REF] prove the NP-hardness of their problem using a reduction from the cumulative problem. This proof is still valid in our case and then the CECSP with concave functions is NP-hard. Clearly, the approximation is tighter with the concave piecewise linear function. Note also that this result is also valid for approximation from above or at the median point.

In order to define an efficient solution method for our problem, we prove that if there exists a solution for the CECSP, then a solution where all resource usage functions, b i (t), are piecewise constant exists. This is the statement of the following theorem: Theorem 1. Let I be a feasible instance of CECSP, with non-decreasing, concave piecewise linear functions f i . A solution such that, for all i ∈ A, b i (t) is piecewise constant, exists. Furthermore, ∀i ∈ A the breakpoints of b i (t) can be restricted to the start and end times of the tasks.

Let us assume a solution S with non-piecewise constant b i (t) exists. Then, S can be transform into a new solution S with piecewise constant b i (t), ∀i ∈ A.

Before proving Theorem 1, we prove that, if ∃[t 1 , t 2 ] where b i (t) is not constant, then b i (t) can be set to its mean value over [t 1 , t 2 ]. Doing so, i will consume the same resource quantity in [t 1 , t 2 ] and will received more energy.

Lemma 1. Let b iq = t 2 t 1 bi(t)dt t2-t1
. Then, we have:

t2 t1 b iq dt = t2 t1 b i (t)dt (5) t2 t1 f i (b iq )dt ≥ t2 t1 f i (b i (t))dt (6)
Proof. Equation ( 5) is trivially verified by replacing b iq by its value. To prove that equation ( 6) is satisfied, we can use the following theorem, due to Jensen [START_REF] Jensen | Sur les fonctions convexes et les inégalités entre les valeurs moyennes[END_REF].

Theorem 2 (Jensen). Let α(t) and g(t) be two integrable functions on

[t 1 , t 2 ] such that α(t) ≥ 0, ∀t ∈ [t 1 , t 2 ]. We have: φ t2 t1 α(t)g(t)dt t2 t1 α(t)dt ≥ t2 t1 α(t)φ(g(t))dt t2 t1 α(t)dt (7)
where φ is a continuous concave function in

[ min t∈[t1,t2] g(t), max t∈[t1,t2] g(t)].
Replacing φ(t) by f i (t), g(t) by b i (t) and α(t) by the constant function equal to 1 gives the desired result. Consider an interval [t 1 , t 2 = t 1 + 6] and functions: Proof (Theorem 1). Let S be a feasible solution of I and let (t q ) {q=1..Q} be the increasing series of distinct start and end time values (Q ≤ 2n). S can be defined as follows:

b i (t) = 3 if t ∈ [t 1 , t 1 + 3[ 1 if t ∈ [t 1 + 3, t 2 ] which yields f i (b i (t)) = 6 if t ∈ [t 1 , t 1 + 3[ 3 if t ∈ [t 1 + 3, t 2 ] t bi(t) t1 i t2 1 2 3 -→ t bi(t) t1 i t2 1
-b i (t) =      b i0 if t ∈ [t 0 , t 1 ] . . . b i(Q-1) if t ∈ [t Q-1 , t Q ] with b iq = t q+1 tq bi(t)dt tq+1-tq -st i = st i -et i = min(τ | τ sti f i (b i (t))dt = W i ) and then b i (t) = 0, ∀t > et i t1 t2 t3 t4 t bi(t) 1 2 4 3 a: bi(t) is not piecewise constant.
-→ t bi(t)

1 2 4 3 t1 t2 t3 t4 b: bi(t) is set to its mean value over [tq, tq+1].
-→ t bi(t)

1 2 4 3 t1 t2 t3 t4
c: eti is adjusted, so each task receives exactly Wi units of energy. S clearly verifies the energy constraints (3) since it is defined in this way. S also satisfies the time window constraints (1) since st i ≤ st i and et i ≥ et i . In addition, as S is a feasible solution, we have ∀q ∈ {1, . . . , Q} and ∀t ∈ [t q , t q+1 ]:

i∈A b i (t) ≤ B ⇒ i∈A tq+1 tq b i (t)dt ≤ B(t q+1 -t q ) ⇒ i∈A b i (t) ≤ i∈A b iq = i∈A t q+1 tq bi(t)dt tq+1-tq
≤ B and S also verifies the resource capacity constraints (4). Finally, we can show that S satisfies the resource requirement constraints (2) in a similar way.

An interesting remark can be made about Theorem 1. Actually, in order to find a solution to CECSP, we only have to find, for each task, its start time st i , its end time et i and the quantity of resource allocated to i between two consecutive start/end times b iq .

Time-Table based reasoning

In this section, we first describe an algorithm which allows us to detect some infeasibilities and then, another algorithm performing time-bound adjustments is presented. Both are based on the well-known Time-Table reasoning [START_REF] Lahrichi | Ordonnancements. La notion de parties obligatoires et son application aux problèmes cumulatifs[END_REF] but the first one combines it with a flow-based algorithm, whereas the second one is an adaptation of the Time-Table Disjunctive reasoning for the cumulative constraints [START_REF] Gay | Time-Table Disjunctive Reasoning for the Cumulative Constraint[END_REF].

The Time-Table reasoning is based on the following observation: if the earliest end time eet i of a task i is higher than its latest start time lst i then i must be in process during interval [lst i , eet i ]; this interval is called the compulsory part of i. For the CECSP, as the resource usage of i is not fixed, we can only deduce that the task will consume at least b min i units of resource (see Figure 6). Aggregating all compulsory parts, we can compute in O(n 2 ) the minimum profile of the resource, denoted by T T (t), ∀t and use it to detect infeasibility and to perform time-bound adjustments.

Time-Table Flow based reasoning

The first algorithm described in this paper embeds the concept of compulsory part into a flow-based linear program. This linear program allows us to detect some infeasibility. Indeed, if the program has no solution, then the considered CECSP instance is infeasible. To describe this program, let (t q ) q∈Q be the increasing series of distinct variable domain bounds, i.e. t q gather all distinct latest/earliest start/end times of all tasks. Clearly, |Q| ≤ 4 • n. Then, ∀[t q , t q+1 ], let β iq (resp. w iq ) represent the quantity of consumed resource (resp. received energy) in this interval. Thus, the following linear program can be used as a checker:

i∈A β iq ≤ B(t q+1 -t q ) ∀q ∈ Q (8)
β iq ≥ b min i (t q+1 -t q ) ∀i ∈ A ; ∀q ∈ Q | s max i ≤ t q ≤ e min i (9)
β iq ≤ b max i (t q+1 -t q ) ∀i ∈ A ; ∀q ∈ Q (10)
β iq = 0 ∀i ∈ A ; ∀q ∈ Q | t q ∈ [r i , d i ] (11)
w iq ≤ a ip β iq + c ip (t q+1 -t q ) ∀i ∈ A ; ∀q ∈ Q ; ∀p ∈ {1, . . . , P i } (12)
w iq ≤ M β iq ∀i ∈ A ; ∀q ∈ Q (13) q∈Q w iq = W i ∀i ∈ A ( 14
)
β iq ≥ 0, w iq ≥ 0 ∀i ∈ A ; ∀q ∈ Q ( 15
)
where M is some large enough constant.

The compulsory part constraints are expressed by constraints [START_REF] Jensen | Sur les fonctions convexes et les inégalités entre les valeurs moyennes[END_REF]. Constraints (8) model the resource capacity limitations. Constraints [START_REF] Lahrichi | Ordonnancements. La notion de parties obligatoires et son application aux problèmes cumulatifs[END_REF] impose that the maximum resource requirements are satisfied. Constraints [START_REF] Lewis | Algebra symposium: Optimizing fuel consumption[END_REF] set the resource consumption of task i in [t q , t q+1 ] to be equal to 0, if [t q , t q+1 ] ⊆ [r i , d i ]. Constraints [START_REF] Nattaf | Energetic reasoning and mixed-integer linear programming for scheduling with a continuous resource and linear efficiency functions[END_REF] together with constraints (13) ensure a correct resource conversion. Indeed, as f i is concave and piecewise linear, the first constraints ensure that w iq can take the value f i ( βiq tq+1-tq ) • (t q+1 -t q ) whereas the second ones set w iq to zero if β iq = 0. Finally, constraints [START_REF] Nattaf | Flow and energy based satisfiability tests for the Continuous Energy-Constrained Scheduling Problem with concave piecewise linear functions[END_REF] state that the tasks received the required energy. Note that this test ensures the same level of consistency than the Time-Table Overload Check [START_REF] Wolf | O(n log n) Overload Checking for the Cumulative Constraint and its application[END_REF]. Indeed, it does the Overload Check at the same time as the Time-Table. This is the same consistency check that performs Time-Table Edge-Finding [START_REF] Vilím | Timetable Edge Finding Filtering Algorithm for Discrete Cumulative Resources. Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems[END_REF]. We now describe the adaptation of the Time-Table Disjunctive reasoning for the cumulative constraint [START_REF] Gay | Time-Table Disjunctive Reasoning for the Cumulative Constraint[END_REF] to our problem. This reasoning will be used to perform time-bound adjustments on the start and end time variables.

Time-Table Disjunctive reasoning

The idea of the Time-Table Disjunctive reasoning is to take advantage of the minimum resource profile to detect disjunctive pairs of tasks, i.e. tasks that cannot be processed in parallel, dynamically. Indeed, the disjunctive reasoning only considers pairs of tasks that exceed the capacity of the resource if they are scheduled in parallel, i.e. b min i + b min j > B. However, tasks in A \ {i, j} may not leave B units of resource available during the overlap of i and j. Furthermore, if task i has no compulsory part then an additional filtering can be done when starting a task j at est j would make it overlap i in every schedule. This is due to the fact that j cannot contain a time interval that i must overlap. Therefore, the authors in [START_REF] Gay | Time-Table Disjunctive Reasoning for the Cumulative Constraint[END_REF] defined the smallest interval overlapping i in every solution, called the minimum overlapping interval and denoted by moi i . For the CECSP, this interval is exactly [est i , lst i ] \ [est i , eet i [, i.e the smallest interval containing [eet i , lst i ] (cf. Fig. 7). Note that, if a task has a compulsory part, then moi i = ∅. ). Then, using the minimum resource profile to compute the available resource and to search for disjunctive pairs of tasks, the time-table disjunctive reasoning can be stated as follows:

Proposition 1. Let i = j be two tasks having no compulsory consumption and such that b min i

+ b min j + min t∈moii T T (t) > B. If moi i ⊆ [est j , et max j ]
we must have eet i ≤ st j and so est j can be adjusted to eet i . Example 4. Consider the example of Figure 8 adapted from [START_REF] Gay | Time-Table Disjunctive Reasoning for the Cumulative Constraint[END_REF]. Due to the minimum profile, i and j cannot overlap and, since j cannot be scheduled before i (moi i ⊆ [est j , et max j ]), j has to be scheduled after and then we can set est j to eet i . A symmetrical reasoning allows to make adjustment of let j . Furthermore, Proposition 1 can easily be extended to the case where tasks i and j has no compulsory consumption following the method presented in [START_REF] Gay | Time-Table Disjunctive Reasoning for the Cumulative Constraint[END_REF] for the cumulative constraint.

Energetic reasoning

In this section, we present the extension of the well-known energetic reasoning to the CECSP with concave and piecewise linear processing rate functions. This work extends the work of [START_REF] Nattaf | A hybrid exact method for a scheduling problem with a continuous resource and energy constraints[END_REF] for linear processing rate functions. The full characterization of relevant intervals for the time-bound adjustments is also provided, extending the result for the cumulative constraint [START_REF] Derrien | A new characterization of relevant intervals for energetic reasoning[END_REF]. This characterization was not presented in [START_REF] Nattaf | A hybrid exact method for a scheduling problem with a continuous resource and energy constraints[END_REF]. First, we present the central idea of the reasoning and we use it to provide a satisfiability test as well as time-bound adjustments. The second part of this section will be dedicated to the characterization of relevant intervals.

Satisfiability test and time-bound adjustments

The idea of the energetic reasoning is to test whether the available resource within an interval [t 1 , t 2 ] is sufficient to provide the minimum resource quantity needed by each task i in this interval, denoted by β(i, t 1 , t 2 ). If not, then the corresponding instance of the CECSP is infeasible.

Theorem 3 ([5]

). Let I be an instance of CECSP. If there exists (t 1 , t 2 ) such that B(t 2 -t 1 ) -i∈A β(i, t 1 , t 2 ) < 0 then I is infeasible.

In order to compute the minimum resource consumption of task i in [t 1 , t 2 ], we first have to compute the minimum energy requirement of i in this interval, denoted by w(i, t 1 , t 2 ). This minimum requirement is obtained by scheduling as much energy as possible outside the interval [t 1 , t 2 ] while satisfying constraints ( 1)-( 4). This always corresponds to a case where the activity is:

left-shifted: the task starts at est i and is scheduled at its maximum requirement between est i and t 1 , right-shifted: the task ends at let i and is scheduled at its maximum requirement between t 2 and let i , or both: when scheduling at minimum requirement inside [t

expression condition Figure 0 t 1 ≥ eet i ∨ t 2 ≤ est i - W i t 1 ≤ est i ∧ let i ≤ t 2 - W i -f i (b max i )(t 1 -est i ) est i ≤ t 1 ≤ eet i ∧ let i ≤ t 2 8a min   W i -f i (b max i )(t 1 -est i ), max W i -f i (b max i )(t 1 -est i + let i -t 2 )), f i (b min i )(t 2 -t 1 )   est i ≤ t 1 < t 2 ≤ let i 8b max W i -f i (b max i )(let i -t 2 )), f i (b min i )(t 2 -est i ) t 1 ≤ est i ∧ t 2 < let i 8c
The three last cases are illustrated in Figure 9. Thus, the minimum energy requirement in [t 1 , t 2 ] is:

w(i, t 1 , t 2 ) = min ω LS (i, t 1 , t 2 ), ω RS (i, t 1 , t 2 ) (16)
I f i (b i (t))dt ≥ w(i, t 1 , t 2 ) (18) b i (t) ≥ b min i (19)
Indeed, the goal is to find the minimum resource quantity (equation ( 17)) i has to consume in I to receive an energy w(i, t 1 , t 2 ) (equation ( 18)). In addition, we have to make sure that the minimum requirement constraints are satisfied by the solution of the program (equation ( 19)).

Then, Lemma 1 can be used to simplify the program. Actually, the lemma applied to the problem involving only task i implies a solution with constant b i (t), say b i (t) = bi exists and then the program can be rewritten as follows:

minimize bi |I| subject to f i ( bi )|I| ≥ w(i, t 1 , t 2 ) bi ≥ b min i
Thus, we have two cases to consider: either the task can be scheduled at b min i during whole interval I or not. In the first case, we can remove the second constraint and we obtain:

bi = min(b ∈ [0, b max i ] | f i (b) ≥ w(i, t 1 , t 2 )/|I|) ⇒ bi = f -1 i (w(i, t 1 , t 2 )/|I|)
Since f i is concave and piecewise linear, f -1 i can easily be computed 2 and bi

|I| = f -1 i (w(i, t 1 , t 2 )/|I|) • |I| = β(i, t 1 , t 2 ) is equal to: f -1 i w(i, t 1 , t 2 ) |I| =              w(i,t1,t2)-ci0|I| ai0|I| if w(i,t1,t2) |I| ∈ [f i (b min i ), f i (γ i 1 )] w(i,t1,t2)-ci1|I| ai1|I| if w(i,t1,t2) |I| ∈ [f i (γ i 1 ), f i (γ i 2 )] . . . w(i,t1,t2)-c iP i |I| a iP i |I| if w(i,t1,t2) |I| ∈ [f i (γ i Pi ), f i (b max i )]
The second case corresponds to the case where executing the task at b min i during the whole interval I give too much energy to the task, i.e. w(i, t

1 , t 2 ) < |I| • f i (b min i ). In this case, β(i, t 1 , t 2 ) is equal to b min i • w(i, t 1 , t 2 )/f i (b min i ).
2 Indeed, suppose that fi is constant in some pieces the first one being p ∈ {1, . . . , Pi}, which would cause problem to compute f -1 i . In that case, fi must be constant on all pieces {p, . . . , Pi} due to concavity. It follows that b max i can be reduced to the start of piece p.

Finally, the expression of the minimum resource consumption is:

β(i, t1, t2) = max b min i w(i, t1, t2) fi(b min i ) , max p∈{1,...,P i } 1 aip (w(i, t1, t2) -|I|cip) (20)
Example 5. Consider the instance described in Example 1. In particular, we compute w(i, t 1 , t 2 ) and β(i, t 1 , t 2 ) for task

β(3, 0, 6) = w(3, 0, 6) -1 • 6 2 • 6 • 6 = 21.5 -6 12 • 6 = 31 4
For interval [0, 4], we have w(3, 0, 4

) = W 3 -2 • f 3 (b max 3 ) = 8.5.
Here, the computation of β(3, 0, 4) falls into the case where w(3, 0, 4

) < f 3 (b min 3 ) • • • |I| = 12.
Thus, we have (see Figure 9b): We now describe how this quantity is used to perform time-bound adjustments. These adjustments are similar to the ones for linear functions described in [START_REF] Nattaf | Energetic reasoning and mixed-integer linear programming for scheduling with a continuous resource and linear efficiency functions[END_REF], so we just briefly present them. We start by defining some notation. We denote by β LS (i, t 1 , t 2 ) (respectively

β(3, 0, 4) = b min 3 • w(3, 0, 4) f 3 (b min 3 ) = 1 • 8.5 3 = 17
β RS (i, t 1 , t 2 )) the minimal resource consumption corresponding to ω LS (i, t 1 , t 2 ) (resp. ω RS (i, t 1 , t 2 )).
We have:

β LS (i, t 1 , t 2 ) = max b min i ω LS (i, t 1 , t 2 ) f i (b min i ) , max p∈{1,...,Pi} 1 a ip (ω LS (i, t 1 , t 2 ) -|I|c ip )
and a similar expression for β RS (i, t 1 , t 2 ).

Lemma 2. If t 1 > r i and there exists [t 1 , t 2 ] such that:

j∈A j =i β(j, t 1 , t 2 ) + β RS (i, t 1 , t 2 )) > B(t 2 -t 1 )
then,

lst i ≤ t 1 -1 b max i ( j∈A j =i β(j, t 1 , t 2 ) + β RS (i, t 1 , t 2 ) -B(t 2 -t 1 ))) and, if b min i 0, let i ≤ t 1 + 1 b min i (B(t 2 -t 1 ) -j∈A j =i β(j, t 1 , t 2 ))
Indeed, the only configuration for which task i starts after t 1 and leading to the minimum resource consumption inside [t 1 , t 2 ] is if the task is right-shifted. Therefore, j∈A; j =i β(j, t 1 , t 2 ) + β RS (i, t 1 , t 2 ) is the total minimum resource consumption in [t 1 , t 2 ] when task i starts after t 1 . Hence, if this quantity is greater than the quantity of available resource in [t 1 , t 2 ], i has to start before t 1 . Otherwise j∈A; j =i β(j, t

1 , t 2 ) + t2 t1 b i (t) ≥ j∈A; j =i β(j, t 1 , t 2 ) + β RS (i, t 1 , t 2 ) ≥ B(t 2 -t 1 ). Furthermore, j∈A; j =i β(t 1 , t 2 , j) + β RS (i, t 1 , t 2 ) -B(t 2 -t 1 )
is the minimum amount of resource that has to be allocated to i before t 1 . Hence, we can divide this number by b max i to obtain a valid upper bound of the start time of i. Similar arguments lead to the adjustment on let i . An example of such adjustments can be found in [START_REF] Nattaf | Energetic reasoning and mixed-integer linear programming for scheduling with a continuous resource and linear efficiency functions[END_REF].

Relevant intervals for the energetic reasoning

In [START_REF] Nattaf | Energetic reasoning and mixed-integer linear programming for scheduling with a continuous resource and linear efficiency functions[END_REF], the authors present a full characterization of the relevant intervals for the satisfiability test in the case where processing rate functions are linear. These intervals are exactly the ones that are relevant for the case of concave piecewise linear functions. Furthermore, as these intervals are included in the set of relevant intervals for the time-bound adjustments, we only present the second results. Recall that an adjustment can be performed on task i if B(t 2 -t 1 )i =j β(j, t 1 , t 2 ) -β RS (i, t 1 , t 2 ) < 0. Thus, we are looking for all intervals [t 1 , t 2 ] such that the function B(t 2 -t 1 )i =j β(j, t 1 , t 2 ) -β RS (i, t 1 , t 2 ) is negative. Theorem 4 ([4]). B(t 2 -t 1 )i =j β(j, t 1 , t 2 ) -β RS (i, t 1 , t 2 ) is locally minimum in interval [t 1 , t 2 ] only if one of the four following conditions is satisfied:

∃(k, ), δ + β(k, t 1 , t 2 ) δt 1 < δ -β(k, t 1 , t 2 ) δt 1 ∧ δ + β( , t 1 , t 2 ) δt 2 < δ -β( , t 1 , t 2 ) δt 2 (21) ∃k, δ + β(k, t 1 , t 2 ) δt 1 < δ -β(k, t 1 , t 2 ) δt 1 ∧ δ + β RS (i, t 1 , t 2 ) δt 2 < δ -β RS (i, t 1 , t 2 ) δt 2 (22) ∃ , δ + β RS (i, t 1 , t 2 ) δt 1 < δ -β RS (i, t 1 , t 2 ) δt 1 ∧ δ + β( , t 1 , t 2 ) δt 2 < δ -β( , t 1 , t 2 ) δt 2 (23) δ + β RS (i, t 1 , t 2 ) δt 1 < δ -β RS (i, t 1 , t 2 ) δt 1 ∧ δ + β RS (i, t 1 , t 2 ) δt 2 < δ -β RS (i, t 1 , t 2 ) δt 2 (24)
with δ + f δt2 (resp. (see [START_REF] Derrien | A new characterization of relevant intervals for energetic reasoning[END_REF]) and

W i ≤ f i (b min i )(let i -est i ).
First, let H (resp. I ) be the intersection point of line

f i (b min i )(t 2 -t 1 ) = W i - (let i -t 2 )f i (b max i ) and f i (b min i )(t 2 -t 1 ) = W i -(t 1 -est i )f i (b max i ) (resp. W i - (t 1 -est i + let i -t 2 )f i (b max i ) = f i (b min i )(t 2 -t 1 ) and t 2 = d i ). Also, let U (t 1 ) (resp. D(t 1 )) be the point t 2 such that f i (b min i )(t 2 -t 1 ) = W i - (t 1 -est i )f i (b max i ) (resp. f i (b min i )(t 2 -t 1 ) = W i -(let i -t 2 )f i (b max i ) Lemma 3. Let i be a task s.t. W i > f i (b min i )(let i -est i ).
Then, for any fixed t 1 , at most two intervals [t 1 , t 2 ] satisfying the second condition of (21) exist and at most two intervals [t 1 , t 2 ] satisfying the second condition of (22) exist. These intervals are described in Table 1.

function relevant intervals condition β(i, t1, t2) [t1, leti] if t1 ≤ esti [t1, leti] and [t1, D(t1)] else if t1 ≤ I t 1 [t1, U (t1)] and [t1, D(t1)] else if t1 < lsti ∨ t1 < Ht 1 [t1, leti + esti -t1] else if eeti < lsti ∧ t1 ≥ Ht 1 [t1, U (t1)] else if lsti ≤ eeti ∧ t1 ≥ lsti none else if t1 ≥ eeti βRS(i, t1, t2) [t1, leti] if t1 ≤ esti ∧ leti > t1 ≥ lsti [t1, leti] and [t1, D(t1)] if t1 > esti ∧ t1 < lsti none otherwise Table 1: Relevant t 2 for case W i > f i (b min i )(let i -est i ).
with H t1 (resp. I t1 ) being the projection on the x-axis of point H (resp. I ).

Proof. We only present how to obtain relevant t 2 for the second line of Table 1.

The other cases can be obtained in a similar way. In order to prove the lemma, we analyze the variation of t 2 → β(i, t 1 , t 2 ). Fig. 11 represents these variations.

where

expr ip = 1 aip (W i -f i (b max i )(let i -t 2 + t 1 -est i ) -c ip (t 2 -t 1
)) The intervals for which the left derivative is smaller than the right are [t 1 , let i ] and [t 1 , D(t 1 )]. Indeed, since f i is a concave piecewise linear function, we have a ip > a ip+1 and c ip < c ip+1 , and we have

δ -exprip δt2 < δ + exprip+1 δt2
resource used by task i and for the energy received by i between events t e and t e+1 . Since the event-based MILP used for our case is almost the same than the one used for linear functions, we are not describing the model here. The main difference lies in the constraints converting resource into energy. These constraints were as follows in the previous model: To generate instances with concave piecewise linear processing rate functions, we use instances of [START_REF] Nattaf | Energetic reasoning and mixed-integer linear programming for scheduling with a continuous resource and linear efficiency functions[END_REF] with identical functions. First, the instances were solved using the time-indexed mixed integer linear program described in [START_REF] Nattaf | Energetic reasoning and mixed-integer linear programming for scheduling with a continuous resource and linear efficiency functions[END_REF]. In this formulation, the planning horizon is discretized in T time periods of size 1 and a variable b it is used to represent the resource consumption of task i in period t. ] is then divided into P i parts. For each piece p, a random coefficient a ip such that a ip < a ip-1 is generated and c ip is computed to ensure the continuity of the function. Finally, W i is set to W i = T t=1 f i (b it ). We repeat this process until we obtain 80 instances with 10 tasks and 140 instances with 20 tasks. Table 2 presents the results of the hybrid branch-and-bound with a time limit of 7200 seconds. The first row describes the results of the hybrid branch-and-bound with the time-table flow based reasoning (TTFlow), the second row with the energetic reasoning (ER), the third row with the time-table disjunctive reasoning (TTDR) and the last row with both ER and TTFlow tests. For each row, the first column presents the time needed to solve the instances, the second column shows the time spent in the tree, the third column the percentage of solved instances, the fourth column the number of solved MILPs (if the instance is solved), the fifth column the number of explored nodes (if the instance is solved), and the last column describes the percentage of "out of memory". We can see that the time-table flow and the time-table disjunctive reasoning provide the best results. For the 10-task instances, the time-table flow solves the instances faster but the time spent in the tree is higher than for the time- for both algorithms. For the 20-task instances, both algorithms achieved similar performances both in terms of solving time and solved instances.

W ie ≤ a i B ie + c i (t e+1 -t e ) ∀i
The energetic reasoning has the worst performances and the combination of the time-table flow and the energetic reasoning improves significantly the performances of the solution method in comparison with the use of the energetic reasoning alone. However, the time-table flow and time-table disjunctive reasoning is better than these two algorithms. Another remark we can do concerns the huge difference in terms of number of MILPs and nodes between the resolution of the 10-task and 20-task instances with the energetic reasoning. This can be explained by the small number of solved instances for the 20-task instances. Indeed, the number of nodes/MILPs in the table is associated only with solved instances.

In addition, we can see that most of the time is spent in the MILP resolution. Therefore, improving the formulation is one of the main perspectives of this work.

Conclusions

This paper presents the extension of the energetic reasoning and the time-table disjunctive reasoning to the CECSP with concave and piecewise linear processing rate functions. A full characterization of the relevant intervals on which the timebound adjustments of the energetic reasoning has to be applied is provided. Both methods are then embedded in a hybrid branch-and-bound and tested on smallsize instances. The interest of considering concave piecewise processing rate functions is shown through examples and experiments. Furthermore, all the new results described in this paper are still valid for linear functions.

For the perspectives, time has to be spent on the instance generation and on the resolution of larger instances. One way of doing this will be by improving the MILP. Also, the method can be improved by the use of dedicated branching heuristics. Finally, the consideration of objective functions is an important perspective of this work.

  Function f3(b).

Fig. 1 :

 1 Fig. 1: An example of instance and corresponding solution of CECSP.

Example 2 .

 2 Consider the neither concave nor convex function described in Figure 3. We approximate the function from below (to obtain a relaxation) with a linear function (left side of the figure) or with a concave piecewise linear one (right side of the figure).

Fig. 3 :

 3 Fig. 3: Approximation from below of a processing rate function by a linear one (left) and by a concave piecewise linear one (right).

Example 3 .

 3 Consider a task i and the concave piecewise linear efficiency function described in Example 1, f 3 (b).

2 3Fig. 4 :

 24 Fig. 4: Illustration of Lemma 1.

Fig. 5 :

 5 Fig. 5: Construction of S from S.

  i has to be in process within[lsti, eeti]

Fig. 6 :

 6 Fig. 6: Compulsory part of a task i.

Fig. 7 :

 7 Fig. 7: Minimum overlapping interval.

Fig. 8 :

 8 Fig. 8: Example of time-table disjunctive adjustments.

  as much energy as possible is scheduled after t2.

Fig

  Fig.9: Computation of ω LS (i, t 1 , t 2 ).

(b min 3 )

 3 3 and interval [0, 4] and [0, 6]. For interval [0, 6], we have w(3, 0, 6) = W = 21.5 since [est 3 , let 3 ] ⊆ [t 1 , t 2 ]. Then, the computation of β(3, 0, 6) falls into the case where w(3, 0, 6) ≥ f 3 • • • |I| = 18. Thus, w(3, 0, 6)/|I| = 21.5/6 3.583 belonging to interval [f 3 (1), f 3 (2)[= [2, 5[, we have (see Figure 9a):

  of β(3, 0, 6).

3 b:

 3 Computation of β(3, 0, 4).

Fig. 10 :

 10 Fig. 10: Example of computation of β(i, t 1 , t 2 ).

  The efficiency functions f i are generated by randomly selecting a number of pieces P i . The interval [b min i , b max i

  . 9: Computation of ω LS (i, t 1 , t 2 ).The minimum energy requirement is then used to compute the minimum resource requirement. First, let I = [est i , let i ] ∩ [t 1 , t 2 ], then β(i, t 1 , t 2 ) can be computed by solving the following program:

	minimize	b i (t)dt	(17)
	I
	subject to

  δ -f δt2 ) the right (resp. left) derivative of f w.r.t. t 2 .This theorem is then used to characterize, for a task i and a fixed t 1 , the value of function t 2 → β(i, t 1 , t 2 ) or t 2 → β RS (i, t 1 , t 2 ) for which its left derivative is greater than its right, and reciprocally for fixed t 2 . Then, the two results are combined to obtain list of relevant intervals [t 1 , t 2 ]. Since there are many cases to consider, we only present how we obtain relevant t 2 for the case where W i > f i (b min

	b min i	= b max i	i	)(let i -est i ). The other cases to consider are

  ∈ A ; ∀e ∈ E and are replaced by the following ones in the new model:W ie ≤ a ip B ie + c ip (t e+1 -t e ) ∀i ∈ A ; ∀p ∈ P i ; ∀e ∈ EExperiments The experiments are conducted on an Intel Core i7-4770 processor with 4 cores and 8 gigabytes of RAM under the 64-bit Ubuntu 12.04 operating system. The hybrid branch-and-bound algorithm is coded in C++ and uses CPLEX 12.6 with 2 threads at each leaf. The heuristic tested for choosing the variable on which the algorithm will branch is the following one: we choose the variable corresponding to the smallest size interval among all [est i , lst i ] and [eet i , let i ]. The parameter used in the experiments is equal to 5 since this parameter value provides the best results for the case of linear efficiency functions[START_REF] Nattaf | Energetic reasoning and mixed-integer linear programming for scheduling with a continuous resource and linear efficiency functions[END_REF].

Table 2 :

 2 table disjunctive reasoning while the number of nodes and MILP solved is equivalent Results of the hybrid branch-and-bound with = 5.

	#tasks	reasoning	total time(s) tree time(s) %solved #MILP #nodes %OOM
		TTFlow	124.4	28.3	96.6	1	14.9	3.4
	10	ER TTDR	1176 244.3	8.26 4.80	84.7 96.6	3346.6 6732.4 1 14.9	1.7 3.4
		ER + TTFlow	542.1	25.4	86.4	1	16.8	3.4
		TTFlow	83.8	27.3	99.2	1.01	37.5	0.8
	20	ER TTDR	4961.9 79.1	9.71 22.6	31.1 99.2	1.03 1	36.9 37.0	3.8 0.8
		ER + TTFlow	834.4	38.5	70.1	1.01	100.8	8.4

Acknowledgments

The authors thank José Verschae for enlightening discussions.

Fig. 11: Relevant intervals for case est i < t 1 ≤ I t1 .

Computational Results

In this section, we start by presenting how we have conducted our experiments on the algorithms stated in this paper. Then, we describe the results of these experiments.

Our propagation algorithms and satisfiability tests were embedded in a hybrid branch-and-bound combining branching scheme and mixed-integer linear programming (MILP). This procedure is an adaptation of the one in [START_REF] Nattaf | Energetic reasoning and mixed-integer linear programming for scheduling with a continuous resource and linear efficiency functions[END_REF] for linear functions.

Hybrid branch-and-bound At first, a branch-and-bound algorithm is used to reduce the size of possible start and end intervals (until their size is less than a given > 0) and, then, an event-based MILP is used in order to find exact task start and end times and to determine the quantity of resource allocated to i between two consecutive events. The branching procedure is as follows. At the beginning, a task can start (resp. end) at any time

). The idea is, at each node, to reduce the size of one of these intervals by creating two nodes splitting the interval into two parts of equal size.

At each node, we apply one or both of the satisfiability tests described above and, if the test does not fail, we perform the corresponding time-window adjustments. We continue this procedure using a depth-first strategy until all intervals are smaller than an . When it happens, the remaining solution space is searched via the event-based MILP.

The MILP used in our algorithm is based on the on/off event-based formulation for the CECSP with linear processing rate functions [START_REF] Nattaf | Energetic reasoning and mixed-integer linear programming for scheduling with a continuous resource and linear efficiency functions[END_REF]. In this formulation, an event corresponds either to a task start or a task end time. These events are represented by a set of continuous variables t e and E = {1, . . . , 2n} represents the index set of these events. We use a binary variable z ie to assign the different event dates to the start and end times of the tasks. Indeed, z ie is equal to 1 if and only if task i is in process during interval [t e , t e+1 ]. Finally, two continuous variables b ie and w ie are also defined. These variables stand for the quantity of


est destinée au dépôt

Boosting probabilistic coastal flood hazard assessment by combining extreme value analysis, full-process based models & metamodels

Jeremy Rohmer, Jessie Louisor, Déborah Idier, Thomas Bulteau, François

Paris, Rodrigo Pedreros

To cite this version:

Third international conference on Advances in Extreme Value Analysis and Application to Natural Hazard (EVAN), 5 to 7 September 2017, Southampton, UK

Boosting probabilistic coastal flood hazard assessment by combining extreme value analysis, full-process based models & metamodels

Jeremy Rohmer 1,* , Jessie Louisor The estimation of marine flooding probability relies on the necessary integration of the probability laws, including extremes, of the sea condition variables (wave characteristics, offshore water level, etc., with potential correlation among them) over the function (denoted f) that defines the inland consequences, typically in terms of water height. To do so, the ideal strategy should combine two key ingredients: (1) high fidelity full process based numerical models to provide accurate simulations of storm responses and a proper estimation of f;

(2) Monte Carlo-based sampling methods to enable a practical integration. The bottleneck of this strategy is the large number (>10,000s) of combinations of offshore hydrodynamic conditions whose responses in terms of inundation should be evaluated using the full-process-based model. A direct approach is hardly achievable due to the computation time cost of a single full process based model run (typically ranging from a few minutes to several hours depending on the flooding processes). A possible alternative can rely on the use of a metamodel (aka surrogate, response surface, etc.), that replaces the longrunning numerical code by a statistical emulator that learn f using a limited number of computer experiments (typically 100-200).

To study the feasibility of this strategy, we focus on the Boucholeurs site along the French Atlantic coast (close to La Rochelle), which is particularly exposed to overflow processes; the most illustrative event being the 2010 winter storm Xynthia. The following steps are followed: (1) Coastal flooding processes are simulated using the finite-difference MARS model including some in-house developments to account for specificities of local flooding (hydraulic processes around culverts and weirs; breaching effect). A single model run takes ~30-60 minutes of computation; (2) Two hundreds of synthetic scenarios are simulated by varying the main characteristics of time-varying offshore conditions (simplified but realistic): surge peak, tide maximum level, phase difference between surge and tide, time characteristics of the surge increase/decrease denoted TIDec; (3) A meta-model of type "random forest" is then built to learn the relationship between these offshore conditions and the water height at four given key locations (camping, holiday camp, marshland, road); (4) A stochastic generator is set up to jointly sample the extreme surge values with the TIDec features using the probabilistic law and dependency model learned from the analysis of the La Rochelle-La Pallice tide gauge measurements; (5) Due to the low computation time cost of the metamodel, the storm response of 10,000s of random samples is evaluated to finally estimate the flooding probability at the four locations.

In addition to the technical practicability, we further discuss the pros and cons of such a strategy compared to the most common approach, consisting in selecting the scenarios to be run from the extreme statistics of the offshore conditions.

1,** , Déborah Idier 1 , Thomas Bulteau 2 , François Paris 2 , Rodrigo Pedreros 1

	1 BRGM, 3, av. Claude Guillemin, BP 36009, 45060 Orleans Cedex 2, France
	*E-mail: j.rohmer@brgm.fr; **Presenting author: j.louisor@brgm.fr
	2 BRGM, Parc Technologique Europarc, 24, avenue Leonard de Vinci, 33600 Pessac, France

Jeremy Rohmer, Jessie Louisor, Déborah Idier, Thomas Bulteau, François Paris, et al.. Boosting probabilistic coastal flood hazard assessment by combining extreme value analysis, full-process based models & metamodels. Third international conference on Advances in Extreme Value Analysis and Application to Natural Hazard (EVAN) " Sep 2017, Southampton, United Kingdom. 2017. hal-


V Boitier
email: vboitier@laas.fr

M Bafleur

P Durand Estebe

J-M Dilhac

Thermogeneration harvesting in an aeronautical application, from specification to realization

Keywords: thermogenerator, energy harvesting, structural health monitoring, aeronautical. I

In this paper, we present results related to a thermoelectric harvesting system developed for powering a SHM datalogger installed in an Airbus A380 closed to the engine.

INTRODUCTION

Aircraft manufacturers and more precisely their R&D department are very interested in the development of wireless solutions for Structural Health Monitoring (SHM) applications. For safety reasons and due to extreme temperatures (high or low), batteries are prohibited in favor of supercapacitors for the storage stage. An energy harvester must be added to balance the energy budget and compensate the SHM datalogger energy consumption. In this work, a solution using a thermogenerator is presented. The SHM sensors are located in the Aft Pylon Fairing (APF) of the A380 aircraft; the pylon holds the engine to the wing (fig. 1). The TEG is attached to the pylon primary structure; this structure is heated by the engine at temperature T H . However, the ambient air temperature in the pylon T c is lower than those of the primary structure giving a thermal gradient estimated between 100 and 200C during flight (cruise) In a previous paper [START_REF]Thermoelectric harvester to power SHM sensors in aircraft pylon[END_REF], specifications have been presented. We summarize hereafter the main requirements: temperature (up to 300°C), strong vibrations (12g), lightning protection and no destruction in case of engine fire (1100°C) according to DO160-G Standard. More details can be found in [START_REF]Energy Autonomy of Batteryless and Wireless Embedded Systems -Aeronautical Applications[END_REF].

II. TEG

Very few commercial TEGs are able to handle 300°C. We chose a 50x50x4.4mm 3 Kelk thermogenerator. The TEG assembly with its heatsink has been previously presented in [START_REF]Thermoelectric harvester to power SHM sensors in aircraft pylon[END_REF]. The TEG is installed (fig. 2) in the pylon primary structure. Ambient air and aircraft skin temperature sensors are implemented close to the TEG. The TEG (standalone and with its heatsink) has been tested, characterized and modeled before flight in laboratory under expected thermal conditions (T = T H -T C between 100 and 200°C). Results are summarized in Table 1. We must notice, the gradient (T TEG ) across the TEG is smaller than T and in laboratory conditions we get T TEG  1/10.T. As shown in fig. 3, to maximize energy transfer, we simply chose a fixed working point (see LT3105 working principle) at 0.35V leading to harvest always more than 30 mW if T > 80°C at ambient pressure.

III. HMEC BOARD

Distant from the TEG, the management circuit called HMEC (Health Monitoring Energy Converter), previously described in [START_REF]Thermoelectric harvester to power SHM sensors in aircraft pylon[END_REF], is installed in a less hot zone. Its main components are listed hereafter. It includes an input common mode filter to serve as an EMI protection against lightning risks and a LTC3105 circuit (harvesting function, storage voltage regulation). A LDO allows powering the output enabling comparator. The storage tank is composed of two inseries supercapacitors Maxwell BCAP0003 3.3F, 2.3V (equivalent capacitance 1.65F with 4.6V maximum voltage). A low power LT6700 comparator circuit is used to enable the boost stage when storage voltage is high enough. A LTC3539 boost performs output 5V regulation. We added an output protection circuit LTC4362 in order to increase the HMEC reliability.

The HMEC operates successfully between -50°C and 125°C.

IV. RESULTS

Both components (TEG and HMEC) have been correctly operating in flight for 1 year. The SHM datalogger is powered by a local power network or by the harvester as soon as the supercapacitor voltage reaches 4.4V and until it drops down to 1.5V. Unfortunately, the datalogger consumption initially supposed to be around 2.4 mW is in practice 10 times higher at 25 mW. So the harvester cannot power permanently the datalogger as seen on fig. 4 (supercapacitor voltage decreases when the harvester supplies the datalogger as the harvested power cannot compensate the consumption). On this figure, the TEG output power is around 12mW and the harvested power in supercapacitor is around 6 mW: LT3105 efficiency is around 44%. The datalogger is powered by the harvester 15% of the time.

The flight analysis allowed highlighting the main influence of the altitude in the convection phenomenon (fig. 5). In fact despite a high temperature gradient (around 100°C), the harvested power in supercapacitor drops to 0 mW if P < 200 hPa. As a consequence, the supercapacitor voltage decreases because of the self-consumption of the enable electronic circuit ( 90µW as Vsupercap < 2.4V). Pressure effect is well seen too on fig. 6. On this figure, T is always between 65 and 80°C. The designed energy harvester module is robust and operates properly although temperatures and vibrations are significant which is a good result. Unfortunately, due to imprecise starting specifications and false environmental data it could not power the load all the time long despite the implemented safety margin in the design. At high altitude, convection reduction must be taken into account when ventilation is limited.

FIGURE 1 .

 1 FIGURE 1. Localization of the thermogenerator (TEG) and the associated electronic device (HMEC). Electronics is deported in a cooler zone next to SHM dataloggers.

FIGURE 2 .

 2 FIGURE 2. TEG and its copper heatsink installed into the A380.

FIGURE 3 .

 3 FIGURE 3. Power vs Voltage of the TEG. Experimental characterization (dotted lines) and simulations results (colored lines).

FIGURE 4 .

 4 FIGURE 4. Flight 1134 Altitude, pressure, temperature gradient vs time during a low altitude flight.

FIGURE 6 .

 6 FIGURE 6. Harvested power stored vs flight pressure under different gradients of temperature. Il manque la valeur des ∆T V. CONCLUSION

TABLE 1 .

 1 TEG PARAMETERS

VI. ACKNOWLEDGEMENTS

This work has been funded by the French program "Investissements d'Avenir" within the framework of CORALIE project.

VII. REFERENCE


Adel Belkadi
email: adel.belkadi@univ-lorraine.fr

Hernan Abaunza
email: habaunza@hds.utc.fr

Laurent Ciarletta
email: laurent.ciarletta@loria.fr

Pedro Castillo
email: castillo@hds.utc.fr

D Theilliol
email: didier.theilliol@univ-lorraine.fr

Distributed path planning for controlling a fleet of UAVs : application to a team of quadrotors

Keywords: Fleet control, Flight Training, Optimization, PSO algorithm, Path generation

come

INTRODUCTION

In recent years, the control of multiple vehicles has attracted increasing attention from scholars around the world, owing to its potential applications whether in the military or civilian fields and by the fact that many benefits can be obtained when a single complicated vehicle is equivalently replaced by multiple yet simpler vehicles. To tackle this issue, two approaches are commonly adopted in literature, the first denoted centralized approach [START_REF]Flight formation control[END_REF] [START_REF] Chiaramonti | Formation control laws for autonomous flight vehicles[END_REF] which is a direct extension of the traditional single-vehicle control is based on the assumption that a central station controls a whole group of vehicles. The second approach, called a distributed approach, which is the one considered in this paper does not require a central station for control and allows for a better adaptation to many physical constraints such as limited resources, short chains of wireless communication and autonomy problems [START_REF] Olfati-Saber | Flocking for multi-agent dynamic systems: Algorithms and theory[END_REF] [START_REF] Park | Cucker-Smale flocking with inter-particle bonding forces[END_REF]. The recent results in this area are categorized into several directions, such as consensus (synchronization [START_REF] Li | Consensus of multiagent systems and synchronization of complex networks: a unified viewpoint[END_REF] [START_REF] Schenato | Average TimeSynch: A consensus-based protocol for clock synchronization in wireless sensor networks[END_REF], rendezvous [START_REF] Corts | Robust rendezvous for mobile autonomous agents via proximity graphs in arbitrary dimensions[END_REF][8] [START_REF] Hui | Finite-time rendezvous algorithms for mobile autonomous agents[END_REF]), formation control and flocking [START_REF] Reynolds | Flocks, herds and schools: A distributed behavioral model[END_REF]- [START_REF] Belkadi | Robust flocking con-trol design for a fleet of autonomous agents[END_REF], optimization [START_REF] Nedic | Distributed subgradient methods for multi-agent optimization[END_REF] [START_REF] Rabbat | Distributed optimization in sensor networks[END_REF], and task assignment [START_REF] Michael | Distributed multi-robot task assignment and formation control[END_REF] [START_REF] Berman | Optimized stochastic policies for task allocation in swarms of robots[END_REF].

⋆ Research supported by the Conseil Regional de Lorraine, the Ministére de L'Education Nationale, de l'Enseignement Supérieur et de La Recherche, and the National Network of Robotics Platforms (ROBOTEX), from France

In the operation of Unmanned Aerial Vehicles (UAVs), the ability to coordinate their trajectories in a multiple vehicle system is vital to carry out different type of missions. In this paper, the issue of trajectory planning is addressed for UAVs operating in a hostile environment. Specifically, the objective is to get a set of UAVs to fly on a 2D plane while converging to a predefined spatial configuration around a target point while ensuring collision and obstacle avoidance between team members. The approach is formulated as a distributed optimization problem. Each UAV is equipped with an independent optimization algorithm which, by minimizing a cost function for each UAV allows the fleet to reach the defined objectives. This distributed planning path was introduced in previous works [START_REF] Belkadi | UAVs team flight training based on a virtual leader: Application to a fleet of Quadrirotors[END_REF] [START_REF] Belkadi | UAVs fleet control design using distributed particle swarm optimization: A leaderless approach[END_REF]. In order to assume an on-line approach, it is necessary to consider an efficient optimization algorithm in terms of resource or time consuming for an embedded system. One of the best performing algorithms in this area is the Particle Swarm Optimization algorithm (PSO) [START_REF] Kennedy | Swarm intelligence[END_REF]. The PSO has demonstrated optimal performance by minimum search term with a minimum calculation time and a relatively easy implementation. In [START_REF] Eberhart | A new optimizer using particle swarm theory[END_REF] authors reported that PSO algorithm could yield significant improvement in solving Trajectory Planning Problem (TPP) and in [START_REF] Guo | Trajectory planning of redundant robot manipulators using QPSO algorithm[END_REF] authors introduced Quantum-PSO (QPSO) to design a collision-free trajectory for planar redundant manipulator. In [START_REF] Chyan | Obstacle avoidance control of redundant robots using variants of particle swarm optimization[END_REF] four variants of Particle Swarm Optimization are proposed to solve the obstacle avoidance control problem of redundant robots.

The rest of paper is organized as follows. The details of the problem formulation are discussed in Section 2. In Section 3, we present the proposed distributed path planning solution. Then, Section 4 illustrates the performance of the method based on various experiments. Finally, a conclusion and perspectives are presented in the last section.

PRELIMINARIES AND PROBLEM STATEMENT

In general, an agent refers to a dynamic system. In this paper, the term 'agent' is interchangeable with "quadrotor", which corresponding non-linear dynamics. Each agent is described in the state space as follows:

Ẋi (t) = f (X i , U i ) (1) where X i (t) ∈ R m and U i (t) ∈ R k : m, k ∈ N A fleet
consists of a set of agents N = 1, ..., n where n is the number of quadrotors. Information exchange among agents is modeled by an undirected graphs. Let describe system by a weighted graph G(N, ϵ(t), A(t)) where N = 1, ..., n and ϵ(t) ∈ N × N are the set of nodes and the set of arcs respectively of the graph G and A(t) a matrix of M n (R) of elements in R as a ij (t) > 0 if : (j, i) ∈ ϵ(t). At each sample, i is an incoming neighbor of j and j is an outgoing neighbor of i. G is assumed to be an undirected graph such as :

∀(i, j) ∈ ϵ(t) ⇒ (j, i) ∈ ϵ(t). In this case, matrix A(t) is symmetrical with a ij (t) = a ji (t).
First of all, let us define for each i node a set of neighbors Ξ i (t) such that:

Ξ i (t) = {j ∈ N : ∥x j (t) -x i (t)∥ ≤ l} (2) where l defines the scope of the neighborhood and Ξ i (t) is called the metrical neighborhood of the agent i.

Influence between agents is represented by a ij (t) elements of the matrix A(t) such that:

a ij (t) = 1 + exp ( c -d ij (t) σ ) (3)
The value of the a ij (t) function depends on the difference between d ij (t) the distance between two agents i and j with the safety distance c (where c > l), and it is all the greater when a ij (t) < c and converges to the value 1 when

d ij (t) >> c
∂(t) = ∂ d (4)
(2) Collisions between interaction agents are avoided ∀i, j ∈ N, i ̸ = j :

∥x i (t) -x j (t)∥ > c (5) (3

) Ensuring the obstacle avoidance

The issue is formulated as an online optimization problem where by minimizing a certain cost function it can control a fleet of UAVs. A distributed planning path strategy based on Particle Swarm Optimization (PSO) algorithm is implemented to achieve this goal and was introduced in our previous work [18][19]. In this paper, experimental tests are carried out to consolidate the theoretical part and at the same time adjustments are added to the optimization block. This includes the main contribution of this paper with obstacle avoidance part included in the control of the fleet.

CONTROLLERS DESIGN

Definition of the cost function

In order to define in deep the developed method, let us consider a fleet of n UAVs (of the quadrotor type) operating in the R 2 space. This multiple rotary wing UAV has four motors arranged in the form of a cross of equal sizes. Its structure offers it the advantage of vertical take off and landing, and the possibility of stationary flights.

At first, let us define a cost function Λ i (t) for each U AV i such that:

Λ i (t) = ρ ( ∥P d -[x i (t) + h]∥ -d d ip ) + q ∑ j=1 a ij (t) ( ∥x j (t) -[x i (t) + h]∥ -d d ij )
)

+ m ∑ k=1 Obs ik (t) (6)
with i ̸ = j, q = card(Ξ(t)), h ∈ R 2 , and ρ >> 1, d d ip representing the desired distance between U AV i and target point with m the number of obstacles detected by the ith UAV.

The main goal is to find the best vector h for each agent i minimizing the cost function Λ i (t) such that:

∀i ∈ N : lim t→∞ Λ i (t) = 0 ⇒ lim t→∞ ∂(t) = ∂ d (7)
The reference trajectory at time t + τ for the agent i will then be:

Y id (t + τ ) = Y i (t) + h (8)
The ρ coefficient in the cost function should be larger than 1. This choice is due to the fact that we favor each agent i position that tends to be in the direction of the target point P d . During the evolution of the agents in R 2 space, noncollision of the agents should be ensured. This constraint is introduced in the cost function Λ i (t) through function a ij (t) which is all the greater when d ij (t) < c and converges to the value 1 when a ij (t) >> c. So, each agent i is more likely to favor the values of h which avoids the need for the distances d ij (t) < c and minimizes the cost function Λ i (t), thanks to the PSO algorithm.

To take into account obstacle avoidance in the fleet control strategy, an additional term is introduced into the cost function such that:

Obs ik (t) = exp   c - ( ∥x kobs (t) -[x i (t) + h]∥) ) σ   (9)
The value of the Obs ik (t) function also depends on the difference between the distance between agent i and obstacle k with the safety distance c, and it converges to zero when the distance is much greater than the safety distance. It is assumed that the obstacle positions are unknown and are detected at a distance d obs .

The algorithm 1 shows, in more detail, the most important steps considered in the proposed strategy to manage the fleet of UAVs.

x id (t + τ ) = x i (t) + h x y id (t + τ ) = y i (t) + h y end if until End of experimentation

Definition of the control law

In order to validate the proposed PSO algorithm in real UAVs, an appropriate position control law needs to be selected such that the Y id (8) is correctly followed.

The selected scheme is a quaternion-based control, which separates the rotation and translation dynamics. Consider the vehicle's dynamic model as

d dt    p ṗ q ω    =       ṗ q ⊗ bF th m ⊗ q * + g 1 2 q ⊗ ω J -1 (τ -ω × J ω)       (10)
where p ∈ R 3 and ṗ ∈ R 3 are the vehicle's position and velocity vectors, F th defines the thrust generated by the quadrotor motors and b = [ 0 0 1 ] T its direction, m and g represent the vehicle's mass and gravity vector, respectively, q describes the vehicle's orientation as a quaternion, ω denotes is the angular velocity, J introduces the inertia matrix with respect to the body-fixed frame and τ is the torque vector caused by the action of the motors.

A position law which stabilizes [START_REF] Reynolds | Flocks, herds and schools: A distributed behavioral model[END_REF] to a desired reference is then defined as u pos = -K pos (x pos -x pos d ), [START_REF] Ahn | Oncollisionavoiding initial configurations to cucker-smale type flocking models[END_REF] where x pos = [ p ṗ ] T is the vehicle's translational state vector, K pos ∈ R 3 denotes a diagonal gain matrix with all positive entries, and

x pos d = [ Y d (t) Ẏd (t) ]
T is the desired reference vector, in which Y d (t) is the PSO trajectory computed from equation [START_REF] Su | Rendezvous of multiple mobile agents with preserved network connectivity[END_REF], and its derivate Ẏd (t) = 0 is neglected.

The translational control law is then incorporated using a quaternion attitude trajectory q d as

q ′ = (b • u pos + ||u pos ||) + b × u pos q d = q ′ ||q ′ || F th = ||u p || , ( 12
)
Then, an attitude control law is defined as:

u = -K att (x att -x att d ) (13
) where K att ∈ R 3 denotes the control gains, x att = [ (2 ln q) T ω T ] T , and the attitude reference is given by

x att d = [ (2 ln q d ) T ( d dt 2 ln q d ) T ]T
For more details, please refer [START_REF] Carino | Quadrotor quaternion control[END_REF].

PRACTICAL RESULTS

Our proposed algorithm was first introduced on a single quadrotor to validate its feasibility and performance using real time experiments, afterwards, three quadrotor were considered in a simulated environment to corroborate the algorithm's behavior on a fleet.

Firstly, the PSO technique is used by one quadrotor to plan its path required to converge to a desired position. Next, unknown disturbances are added. The PSO algorithm re-computes the required trajectory to recover from the perturbations. Then, obstacles are added such that the quadrotor changes its path to avoid collisions, the location of the obstacles is considered to be known.

finally, three quadrotor are introduced, each one computes its own trajectory such that a uniform fleet is formed around a given target while avoiding collisions amongst each other. This last scenario was validated in a simulated environment.

The simulation parameters considered in this paper are illustrated in Table 1.

ρ σ c d obs 5 0.4 0.8m 3m h xi , h xi τ d d ij d d ∈ [-0.4m, 0.4m] 0.1s 1.732 0 Table 1. PSO considered parameters
The quadrotor used in all the experiments is a Parrot ARDrone2, which firmware was modified to be capable of running our custom code. The framework used in this drone is part of the platforms created by the team of the Heudiasyc Laboratory at the Universit de Technologie de Compigne.

One agent, no obstacles, no disturbances

In this first experiment, three consecutive targets are given to the quadrotor using a ground station. The optimization algorithm computes the desired position taking into account the dynamic constraints of the UAV, since the position is calculated in each iteration, it forms a trajectory which makes the UAV reach the desired point in a smooth manner.

Figure 1 illustrates the desired target points, the computed trajectory, and the real position of the UAV in a 2D space. The trajectory generation and tracking is represented in Figure 2, note that, even if when the target's position changes abruptly, the PSO trajectory converges smoothly to the desired reference. The simulations and experiments can be watched at the following link:

https://www.youtube.com/watch?v=QLRl0Js72L0

One agent with disturbances

To assess the robustness of the trajectory generator, the quadrotor was subjected to external disturbances. Instead of only compensating the increased error, the PSO algorithm re-computes the trajectory to recover from the perturbation.

a) The first experiment considers the same scenario as in section 4.1, but now, a disturbance is added while the UAV is following the trajectory, see Figure 3. Figure 6 illustrates a member of our team displacing the drone to a position far from the desired reference (which was set to [0,0]) with one hand.

Observe on figure 7 that the PSO algorithm computes a trajectory to return to the reference target, but since it stays close to the UAV's position, the recovery from the disturbance is smooth.

One agent in the presence of obstacles

In the next experiment, a barrier is considered to be across the path of the quadrotor. When the vehicle approaches near the obstacle, the PSO algorithm adapts the trajectory to avoid a collision, see Figure 8.

Fleet of three quadrotors

Finally, the PSO algorithm was applied in a fleet of UAVs consisting of three quadrotors in a simulated environment, see Figure 10.

In this simulation, the target point's position P d and the desired configuration for the UAVs from the matrix ∂ d are fixed. The elements d ij (t) are defines such that the desired distance between UAVs is homogeneous, thus arriving to a triangle configuration around the target point. A second simulation was then performed to illustrate the trajectory followed by each agent to arrive to the formation of the fleet. Figure 12 illustrates the path followed by each UAV.

Note the triangular formation of the fleet when the target is reached, and the agents reaction when they come close to each other. An important issue that needs to be addressed in the fleet formation problem is convergence of the distance between agents to a desired value. In this case, an uniform triangular formation is expected with a desired distance of 1m towards the target, thus implying using simple geometry that the distance between agents is d agents = 2(1m) cos(π/6) = 1.732m Fig. 13. Desired distance between agents in a simulated experiment Figure 13 illustrates the convergence of all the distances between agents to the desired value, thus ensuring the fleet formation around the target. It should be noticed that the safety distance of the neighborhood is respected throughout the simulation.

Position in X axis [m]

CONCLUSION

This paper has presented a distributed path planning approach for the control of a fleet of Unmanned Aerial Vehicles. The control objective has been formulated as an optimization problem based on a Particle Swarm Optimization algorithm (PSO). The minimization of a cost function allowed to generate in a distributed way the trajectories of the UAVs which guarantee the control of the fleet. The main advantage of the proposed approach is that it is only dependent on the trajectory generation part and it is independent from the model or from the control law of the UAVs, which makes it easier to use.

Simulation and experiment results have clearly shown that the proposed method is effective for the surface coverage and have also shown a robust approach against the collisions and obstacle avoidance. From those first results, we intend to develop this algorithm more in terms of formalization and implement it soon on control of fleet of real quadrotors in a 3D space. Another perspective would be to adapt the algorithm parameters to consider the case of a moving target.

Fig. 1 .

 1 Fig. 1. 2D targets, trajectory and drone's position for the first experiment.

Fig. 2 .

 2 Fig. 2. Target reference, trajectory, and UAV's position for each axis.

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. One of the members of our team pushing the UAV to generate a disturbance in flight experiments.

Fig. 6 .Fig. 7 .

 67 Fig.6. A member of our team inducing large perturbations on the quadrotor.

Fig. 8 .

 8 Fig. 8. Vehicle's trajectory correction to avoid collision with an obstacle.

Figure 9 Fig. 9 .

 99 Figure 9 illustrates the generation of the trajectory for both x and y axes in the previous scenario, note the corrections made to the trajectory mainly in the y axis to avoid the obstacle.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Simulated environment for fleet of three dronesThe simulated experiment consisted on a set of four reference points, given to each quadrotor using a ground station. Each quadrotor then computes independently its trajectory taking into account the position of the other agents.A first simulation consisted in initializing the quadrotors in random position, then a sequence of four target points are given at different times. The UAVs then form a fleet around the reference ensuring a safety distance between each other. Figure11represents this behavior.

Fig. 12 .

 12 Fig. 12. A second simulation shows the trajectory followed by each UAV to converge to the desired formation

  .

	The sets of elements d ij (t) representing the distances be-
	tween agents forms the matrix ∂(t) at time t, and d d ij the
	desired distances between agents forms the matrix ∂ d . It
	is also assumed that each agent can detect an obstacle at
	a distance l ′ .
	The control objective is then:
	(1) To bring the fleet from an initial geometrical configu-
	ration ∂(t

0 ) to a desired geometrical configuration ∂ d around a target point lim t→+∞

  Minimization of the cost function Λ i (t) by PSO Algorithm . Chose the vector h[h x , h y ] that minimizes the cost function Λ i (t) . Calculate the next position x id and y id at time t + τ :

Algorithm 1 Trajectory generation . Initialize parameters . Define the desired configuration for the UAVs with the matrix ∂ d repeat {at each step time} . Calculate the matrix ∂(t) if ( (∃d ij (t) < c)or(t%τ == 0) ) then . Find all topological neighbors for quadrotor i by set Ξ i (t) . Detect obstacles for quadrotor i .