Clément Pira
email: pira@laas.fr

Christian Artigues
email: artigues@laas.fr

Propagation mechanism for a non-preemptive strictly periodic scheduling problem

Keywords: Periodic scheduling, equilibrium, propagation mechanism

We study a non-preemptive strictly periodic scheduling problem. This problem, introduced in [6, 4], arises for example in the avionic field where a set of N periodic tasks (measure of a sensor, etc.) has to be scheduled on P processors distributed on the plane. In the related field of cyclic scheduling [5], some notable success [3] has been recently achieved through the use of specific constraint propagation. In this article, we consider a heuristic, first proposed in [1, 2], which is based on the notion of equilibrium. Following a game theory analogy, each task tries successively to optimize its own schedule and therefore to produce the best response, given the other schedules. We present a propagation mechanism for non-overlapping constraints which significantly improves this heuristic.

• This work was

Periodic scheduling problem and its MILP formulation

We consider a non-preemptive strictly periodic scheduling problem introduced in [START_REF] Korst | Periodic multiprocessors scheduling[END_REF][START_REF] Eisenbrand | Solving an Avionics Real-Time Scheduling Problem by Advanced IP-Methods[END_REF][START_REF] Sheikh | Resource allocation in hard real-time avionic systems -Scheduling and routing problems[END_REF][START_REF] Sheikh | Strictly periodic scheduling in IMA-based architectures[END_REF]. Such a problem arises in the avionic field, where a set of N periodic tasks (measure of a sensor, etc.) has to be scheduled on P processors distributed on the plane. In this problem, each task i has a fixed period T i and a processing time p i . They are subject to non-overlapping constraints : no two tasks assigned to the same processor can overlap during any time period (see Figure 1). A solution is given by an assignment of the tasks to the processors and, for each task by the start time t i (offset) of one of its occurrences. In the following, unlike the papers mentioned above, we adopt a more general model in which processing times p i are generalized by positive latency delays l i,j ≥ 0 (or time lags). The former case is the particular case where l i,j = p i for all other tasks j.

The proposed heuristic adapts easily to this generalization.

In this paper, we only consider the case where the offsets t i are integers. This hypothesis is important to prove convergence of the method described in section 2.1. Since the problem is periodic, all the periods T i also need to be integers 1 .

However the latency delays need not be integer a priori, especially in the case of the optimization problem presented in section 1.1. Non-overlapping constraints In the formulation of a monoprocessor periodic scheduling problem, a latency delay l i,j ≥ 0 has to be respected whenever an occurrence of a task j starts after an occurrence of a task i. Said differently, the smallest positive difference between two such occurrences has to be greater than l i,j . Using Bézout identity, the set (t j + T j Z) -(t i + T i Z) of all the possible differences is equal to (t j -t i ) + g i,j Z where g i,j = gcd(T i , T j ). The smallest positive representative of this set is (t j -t i ) mod g i,j . This is the only representative of (t j -t i ) + g i,j Z which belongs to [0, g i,j -1] (in particular, we consider a classic positive modulo, and not a signed modulo). Therefore, we simply consider the following constraint : (t j -t i ) mod g i,j ≥ l i,j , ∀(i, j), i = j

Classically, processing times are strictly positive, however our model allows zero delays. Since equation ( 1) is trivially satisfied when l i,j = 0, we introduce an additional graph G containing the arcs (i, j) for which l i,j > 0, and we consider equation [START_REF] Sheikh | Resource allocation in hard real-time avionic systems -Scheduling and routing problems[END_REF] only for those couples. As we will see, the constraints associated with (i, j) and (j, i) work naturally together and some complications are introduced when only one of the couples has a strictly positive delay. To avoid special cases, we will suppose in the following that the graph G is symetric2 : two delays l i,j and l j,i are either both zero, or both strictly positive. With this hypothesis, a consequence of proposition 1 (see Appendix) is that equation (1) can be replaced by the following interval constraint for each couple (i, j) ∈ G with i < j :

l i,j ≤ (t j -t i ) mod g i,j ≤ g i,j -l j,i , ∀(i, j) ∈ G, i < j (2)
Note that except for the modulo, these constraints have a clear analogy with classical precedence constraints of the form t j -t i ≥ l i,j . However, they are much harder to handle since they induce an exclusion between tasks3 .

Objective to maximize In an uncertain context, a task may last longer than expected, due to failure of the processor. This increase in the duration is naturally proportional to the original processing time. To withstand the worst possible slowdown, it is natural to make all the delays l i,j proportional to a common factor α ≥ 0 that we try to optimize (see [START_REF] Sheikh | Resource allocation in hard real-time avionic systems -Scheduling and routing problems[END_REF][START_REF] Sheikh | Strictly periodic scheduling in IMA-based architectures[END_REF]).

max α (3)
s.t. l i,j α ≤ (t j -t i ) mod g i,j ≤ g i,j -l j,i α ∀(i, j) ∈ G, i < j (4) t i ∈ Z ∀i (5) α ≥ 0 (6)
The addition of this objective is also a good way to obtain feasible solutions. In the context of precedence constraints, a way to implement the Bellman-Ford algorithm is to cyclically choose a task and to push it forward, up to the first feasible offset. In section 2.2 we present a method to find the next feasible offset in the case of non-overlapping constraints. However there isn't necessarily a solution and a similar propagation mechanism which would cyclically move one task at a time to the first feasible location would have only a small chance to converge towards a feasible solution. However, with the optimization problem, there is always at least a solution having a zero α-value. We can then try to inflate α, while moving a task, what we call the best-response method : instead of choosing the first feasible location, we search for the best location with respect to this additional objective α. Whenever a solution with α ≥ 1 is found, we have a solution for the feasibility problem. This is essentially the principle of the heuristic presented in section 2.

Structure of the solution set

An optimal solution is completely defined by the offsets, since given the offsets (t i ) we can always compute the best compatible α-value :

α = min (i,j)∈G (t j -t i ) mod g i,j l i,j (7)
If we try to draw the function α = (t i -t j ) mod g i,j /l j,i representing the optimal value of α respecting one constraint (t i -t j ) mod g i,j ≥ l j,i α when t i takes rational values, we obtain a curve which is piecewise increasing and discontinuous since it becomes zero on t j + g i,j Z. In the same way, if we draw α = (t j -t i ) mod g i,j /l i,j , we obtain a curve which is piecewise decreasing and becomes zero at the same points. It is therefore more natural to consider the two constraints jointly, i.e. (4), which gives a continuous curve (see Figure 2) 4 .

Since G is symetric, it can be seen as an undirected graph. We will write G(i) for the set of neighbors of i, i.e. the set of tasks with which i is linked through non-overlapping constraints. With this notation, and given some fixed Fig. 2. Possible values for (ti, α) when constrained by a single task j offsets (t j ) j∈G(i) , the set of solutions (t i , α) is the intersection of the curves decribed above for each j ∈ G(i) (see Figure 3). It is composed of several adjacent polyhedra. We can give an upper bound on the number n poly of such polyhedra. A polyhedron starts and ends at zero points. For a given constraint j, there is T i /g i,j zero points in [0, T i -1], hence n poly is bounded by T i j =i 1 gi,j . This upper bound can be reached when the offsets are fractional, since in this case, we can always choose the offsets t j such that the sets of zero points (t j + g i,j Z) j∈G(i) are all disjoint. In the case of integer offsets, there is obvioulsy at most T i zero points in the interval [0, T i -1] and therefore, at most T i polyhedra. We obtain the shape of a two dimensional slice of the solution set, for a given task i, and some fixed (t j ) j∈G(i) . We will see how these two-dimensional slices help to understand the concept of equilibrium and the heuristic : when solving the best-response method, we try to find the new location t i which maximize the α-value. Following the same idea, the complete full dimensional solution set is composed of several adjacent N + 1-dimensional polyhedra.

MILP formulation of the multiprocessor problem MILP formulation

In this section, we present a MILP formulation of the problem which will be used to compare the performance of our heuristic. As explained above, the modulo (t j -t i ) mod g i,j is the only element t j -t i + g i,j q i,j in the interval [0, g i,j ). Here, q i,j is an additional 'quotient variable' which can always be computed by :

q i,j = ⌈(t i -t j )/g i,j ⌉ (8)
Hence for a given couple (i, j) ∈ G, with i < j, the non-overlapping constraint (4) can be rewritten as l i,j α ≤ t j -t i +g i,j q i,j ≤ g i,j -l j,i α. In the multiprocessor case, the scheduling problem is coupled with an assignment problem. Thus, we introduce binary variables a i,k which indicate if a task i is assigned to a processor k, and variables x i,j which indicate if i and j are on different processors. These variables must satisfy (10) and (11). Since the non-overlapping constraints have to be satisfied whenever two tasks are on the same processor, we replace the term l i,j α by l i,j α -l i,j α max x i,j . When i and j are on the same processor (x i,j = 0), we find back the original term l i,j α, otherwise we obtain a negative term l i,j αl i,j α max which makes the constraint trivially satisfiable. This yields constraints (12) and (13). max α (9)

k a i,k = 1 ∀i (10) x i,j ≤ 2 -a i,k -a j,k ∀k, ∀(i, j) ∈ G, i < j (11) t j -t i + g i,j q i,j ≥ l i,j α -α max l i,j x i,j ∀(i, j) ∈ G, i < j (12)
t j -t i + g i,j q i,j ≤ g i,j -l j,i α + α max l j,i x i,j ∀(i, j) ∈ G, i < j (13)

t i ∈ Z ∀i (14) q i,j ∈ Z ∀(i, j) ∈ G, i < j (15) x i,j ∈ [0, 1] ∀(i, j) ∈ G, i < j (16) a i,k ∈ {0, 1} ∀k, ∀i (17)
Bounds on the offsets and quotients We can see that the previous system is invariant under translation t i → t i + n i T i : if (t i ) is feasible and if (n i ) is a vector of integers, then the solution (t ′ i ) defined by t ′ i = t i + n i T i is also feasible. Since t i mod T i = t i + n i T i for some n i ∈ Z, we deduce that (t i mod T i ) is also a solution. We can therefore impose additionnal bounds t i ∈ [0, T i -1]. More generally, this remains true if we replace T i with an updated period :

T * i = lcm((g i,j ) j∈G(i) ) = gcd(T i , lcm((T j ) j∈G(i) )) (18)
The updated period T * i always divide T i (which we denote T * i |T i ) and we can have a strictly smaller value for example if a prime factor p only occurs in one period T i (or more generally if a factor p occurs with multiplicity m in T i but with strictly smaller multiplicities in any other periods of a task connected to i). Indeed, even if the (T i ) are the initial parameters, only their GCD appears in the constraints. Therefore, if a prime factor occurs in only one period, then it completely disappear in the g i,j . By replacing the initial periods (T i ) by the updated ones (T * i ), we simply remove some irrelevant factors. In the following, we will suppose that the periods have been updated to have no proper factors, i.e. T i = T * i . Computing the image of the interval [0, T i -1] by expression (8) immediately gives associated bounds on the quotients :

1 - T j g i,j ≤ q i,j ≤ T i g i,j (19)
In particular, if T j |T i then q i,j ≥ 0 (in fact q i,j ∈ [0, T i /T j ]), which shows that in the harmonic case, we can impose positive variables. Conversely if T i |T j then q i,j ≤ 1. Finally, if T i = T j , then q i,j ∈ {0, 1} which shows that the monoperiodic case correspond to the case of binary variables.

Upper bound on the α-value Let i and j be two tasks on the same processor. Then, constraint (4) gives the following upper bound : α ≤ g i,j /(l i,j +l j,i ). Taking the integrality assumption into account, we get an even tighter upper bound :

α ≤ α i,j max = max 1 l i,j g i,j l i,j + l j,i l i,j , 1 l j,i g i,j l i,j + l j,i l j,i (20)
These bounds are illustrated on Figure 2. In the monoprocessor case, we deduce that α max = min i,j α i,j max is an upper bound on the value of α. In the multiprocessor case, we have at least the trivial upper bound α max = max i,j α i,j max . However, since α max is used as a 'big-M' in the MILP formulation, we would like the lowest possible value in order to improve the efficiency of the model. For this, we solve a preliminary model, dealing only with the assignment (hence without variables (t i ) and (q i,j )), in which non-overlapping constraints (12-13) are replaced by the following weaker constraint :

α ≤ α i,j max + (α max -α i,j max )x i,j , ∀(i, j) ∈ G, i < j (21)
Intuitively, this constraint indicates that α should be less than α i,j max if i and j are assigned to the same processor (x i,j = 0), otherwise it is bounded by α max which is the trivial upper bound. Solving this model is much faster than for the original one (less that 1s for instances with 4 processors and 20 tasks). It gives us a new value α max which can be used in the original model.

An equilibrium-based heuristic

Heuristic and the mutiprocessor best response method

The multiprocessor best response method The main component of the algorithm is called the best response procedure. It takes its name from a game theory analogy. Each task is seen as an agent which tries to optimize its own assignment and offset, while the other assignments and offsets are fixed. Instead of using a binary vector (a i,k ) as in the MILP formulation, we represent an assignment more compactly by a variable a i ∈ [1, P ]. For each agent i, we want to define a method MultiProcBestResponse i which returns the best assignment and offset for the task i, given the current assignment and offsets (t j , a j ) of all the tasks. In order to choose the assignment, an agent simply tries every possible processor. On a given processor p, it tries to find the best offset, taking into account the non-overlapping contraints with tasks currently one this processor.

More formally, the BestOffset p i procedure consists in solving the following program :

(BO p i ) max α (22) s.t. l j,i α ≤ (t i -t j ) mod g i,j ≤ g i,j -l i,j α ∀j∈G(i) aj =p (23)
t i ∈ Z (24) α ≥ 0 (25)
Note that there are only two variables, t i and α, since the other offsets and assignments (t j , a j ) j∈G(i) are parameters. A method to solve this program will be presented in section 2.2. Following the discussion in section 1.2, the previous system is invariant under tranlation by T i , and even by the possibly smaller value T p i = lcm((g i,j ) j∈G(i)|aj =p ). Hence, we can impose t i to belong to [0, T p i -1]. In the same way, we can compute an upper bound on this program : α p max = min j∈G(i)|aj =p α i,j max . If the current best solution found on previous processors is already better than this upper bound, there is no way to improve the current solution with this processor, hence we can skip it. The multiprocessor bestresponse procedure is summarized in Algorithm 1.

Algorithm 1

The multiprocessor best-response The concept of equilibrium and principle of the method to find one Definition 1. A solution (t i , a i ) is an equilibrium iff no task i can improve its assignment or offset using procedure MultiProcBestResponse i .

The heuristic uses a counter N stab to count the number of tasks known to be stable, i.e. which cannot be improved. It starts with an initial solution (for example randomly generated) and tries to improves this solution by a succession of unilateral optimizations. On each round, we choose cyclically a task i and try to optimize its schedule, i.e. we apply MultiProcBestReponse i . If no improvement was found, then one more task is stable, otherwise we update the assignment and offset of task i and reinitialize the counter of stable tasks. We continue until N tasks are stable. This is summarized in Algorithm 2.

Remark 1. The main reason for the use of integers is that it allows to guarantee the convergence of the heuristic. The termination proof relies on the fact that there is only a finite number of possible values for α (we refer to [START_REF] Sheikh | Resource allocation in hard real-time avionic systems -Scheduling and routing problems[END_REF] for the proof of termination and correction). This is not the case with fractional offsets, for which the heuristic is unlikely to converge in a finite number of steps.

Algorithm 2

The heuristic 1: procedure ImproveSolution((tj)j∈I , (aj)j∈I ) 2:

N stab ← 0 ⊲ The number of stabilized tasks 3:

i ← 0 ⊲ The task currently optimized 4:

while N stab < N do ⊲ We run until all the tasks are stable 5:

(newti, newai, αi) ← MultiProcBestResponse(i, (tj)j∈I ) 6:

if newti = ti or newai = ai then 7:

ti ← newti; ai ← newai ⊲ We have a strict improvement for task i 8:

N stab ← 1 ⊲ We restart counting the stabilized tasks 9:

α ← αi 10: else ⊲ We do not have a strict improvement 11:

N stab ← N stab + 1 ⊲ One more task is stable 12:

α ← min(αi, α) 13:

end if 14:

i ← (i + 1) mod N ⊲ We consider the next task 15:

end while 16:

return (α, (tj)j∈I , (ai)j∈I ) 17: end procedure An equilibrium is only an approximate notion of optimum. Hence, in order to find a real optimum, the idea is now to run the previous heuristic several times with different randomly generated initial solutions, and to keep the best result, following a standard multistart scheme.

The best-offset procedure on a given processor

We now need to implement the BestOffset p i method, i.e. solve (BO p i ). Since t i is integer and can be supposed to belong to [0, T p i -1], we can trivially solve this program by computing the α-value for each of these offsets, using expression [START_REF] Megiddo | Linear-Time Algorithms for Linear Programming in R 3 and Related Problems[END_REF], and select the best one. This procedure runs in O(T i N ), hence any method should at least be faster. In [START_REF] Sheikh | Resource allocation in hard real-time avionic systems -Scheduling and routing problems[END_REF], the authors propose a method consisting in computing the α-value only for a set of precomputed intersection points (in the next section we will see that a fractional optimum is at the intersection of an increasing and a decreasing line). This already improves the trivial procedure.

In the following, we present an even more efficient method.

Local (two dimensional) polyhedron We want to compute the local polyhedron which contains a reference offset t * i . Locally, the constraint (t i -t j ) mod g i,j ≥ l j,i α is linear, of the form t i -o j ≥ l j,i α. Here, o j is the largest x ≤ t * i such that (x -t j ) mod g i,j = 0. In the same way, we can compute the decreasing constraint o ′ j -t i ≥ l i,j α. In this case o ′ j is the smallest x ≥ t * i such that (o ′ j -x) mod g i,j = 0. By proposition 2 (see Appendix), we have :

o j = t * i -(t * i -t j ) mod g i,j and o ′ j = t * i + (t j -t * i ) mod g i,j (26)
Therefore, we obtain a local polyhedron (see figure 4). Note that when (t * it j ) mod g i,j > 0, we simply have o ′ j = o j + g i,j by proposition 1 (see Appendix). However, when (t * i -t j ) mod g i,j = 0, we have o j = o ′ j = t * i . In this case, the polyhedron is degenerated since it contains only {t * i }, and the α-value at t * i is zero. Instead of computing this polyhedron, we prefer to choose either the polyhedron on the right, or on the left (see figure 5). If we choose the one on the right, this amounts to defining o ′ j to be the smallest x > t * i which sets the constraint to zero. By proposition 2, we have o ′ j = t * i + g i,j -(t * i -t j ) mod g i,j . Therefore, this simply amounts to enforcing o ′ j = o j + g i,j . Choosing the polyhedron on the left amounts to defining o ′ j using (26) and to enforcing o j = o ′ j -g i,j . Solving the local best offset problem Once the local polyhedron has been defined, the problem is now to solve the following MILP :

(Loc-BO p i ) max α (27) s.t. t i -l j,i α ≥ o j ∀j ∈ G(i), a j = p (28) t i + l i,j α ≤ o ′ j ∀j ∈ G(i), a j = p (29) t i ∈ Z (30)
We can first search for a fractional solution, and for this we can use any available method of linear programming. However, since the problem is a particular two dimensional program, we can give special implementations of these methods. In the following, we present a simple primal simplex approach, which runs in O(N 2 ) in the worst case but has a good behaviour in practice. A local polyhedron is delimited by increasing and decreasing lines, and the fractional optimum is at the intersection of two such lines. Hence a natural way to find the optimum is to try all the possible intersections between an increasing line, of the form x -l j,i α = o j , and a decreasing line, of the form x + l i,k α = o ′ k , and to select the one with the smallest α-value. The coordinates of these intersection points are given by5 :

x = l j,i o ′ k + l i,k o j l j,i + l i,k and α = o ′ k -o j l j,i + l i,k (31)
Since there is at most N -1 lines of each kinds, the algorithm runs in O(N 2 ). In practice, a better approach (with the same worst case complexity) is to start with a couple of increasing and decreasing lines, and alternatively to try to improve the decreasing line (see Figure 6), then the increasing one (see Figure 7), and so on, until no improvement is made. The overall solving process is illustrated on Figure 8. Remark 2. Suppose we just improved the decreasing line (Figure 6) and we get a new intersection point. We know that the whole local polyhedron lies inside the cone oriented below and defined by the increasing and the decreasing line. Then all the decreasing lines with a smaller slope than the new decreasing line, i.e. with a larger delay l j,i , lie completely outside this cone, and therefore cannot be active at the optimum. We can therefore drop these lines in the subsequent rounds. The same is true for the increasing lines. This gives us a method to compute an integral solution in O(N 2 ) since once a fractional solution has been found, we can deduce an integral solution with an additional computation in O(N ). Indeed, if x is integer, then (x, α) is the desired solution. Otherwise we can compute the α-values α -and α + associated with ⌊x⌋ and ⌈x⌉ and take the largest one. Note that since ⌊x⌋ (resp. ⌈x⌉) is on the increasing phase (resp. decreasing phase), only the corresponding constraints are needed to compute α -(resp. α + ) :

α -= min j∈G(i) aj =p (⌊x⌋ -o j )/l j,i and α + = min k∈G(i) a k =p (o ′ k -⌈x⌉)/l i,k (32)
We call this a primal approach since this is essentially the application of the primal simplex algorithm. However the primal simplex is not applied on (Loc-BO p i ) but on its dual. A dual approach is illustrated on Figure 9. From a theoretical perspective, this approach is outperformed by Megiddo algorithm [START_REF] Megiddo | Linear-Time Algorithms for Linear Programming in R 3 and Related Problems[END_REF] which allows to find a solution in O(N ). However Megiddo algorithm is more complex and generally slower in practice. The case of processing times In the case of processing times, no such refinement as Megiddo algorithm is needed in order to obtain a complexity in O(N ).

In fact, in this case the primal simplex method already runs in O(N ). The reason is that in this case, the algorithm stop after three phases. Indeed, we first search for a better decreasing line. However, in the case of processing times, we have l i,j = p i . Since the delays l i,j gives the slope of the decreasing lines, this implies that all of them have the same slope -1/p i (see Figure 10). Thus, even if initially there are N -1 decreasing lines with equations o ′ j -x = αp i , the one with the smallest o ′ j will be selected after this first phase. In a second phase, we search for a better increasing line. In the third phase, we will not find a better decreasing line, therefore we have the optimal value. Finding the next improving offset In section 1.1, we explained how the equilibrium-based heuristic replaces a mechanism to push a task to the next feasible location, by one which search for the best location with respect to an additional objective α. The former has the drawback that it does not always produce a solution. However, such a procedure is interesting anyway to implement the best-response method. Indeed, after a first local optimization, we obtain a current solution x ref with a value α min . In the rest of the procedure, we are only interested by polyhedra which could improve this value. As α min is improved, more and more polyhedra will lie completely below this level and will be skipped (see darker polyhedra in Figure 11). Starting from a current solution x ref with value α min , we want to find a new integer solution x greater than the current solution which strictly improves its value. Hence, this solution should satisfy :

α min l j,i < (x -t j ) mod g i,j < g i,j -α min l i,j ∀j∈G(i) aj =p (33)
However, since the middle expression gives an integer value, we can round the bounds and obtain the equivalent interval constraint (35). Moreover, among all the possible solutions, we would like to find the smallest one :

min x (34)
s.t. ⌊α min l j,i ⌋ + 1 ≤ (x -t j ) mod g i,j ≤ g i,j -⌊α min l i,j ⌋ -1 ∀j∈G(i) aj =p (35) x ∈ [x ref , x end ) (36)
Here, x end (which will be defined in the next section) indicates when to stop the search (since there is possibly no strictly improving solution). In order to solve this program, we will start with x = x ref and propagate until we find a feasible solution. Since we are only interested by contraints associated with tasks j ∈ G(i) currently on processor p, i.e. satisfying a j = p, we define N p i to be the number of such contraints, and we suppose that these contraints are numbered from 0 to N p i -1. On a given round, we choose a constraint j and check if it is satisfied. We could check the two inequalities of the interval constraint (35), however we can also remark that this is equivalent to :

(x -t j -⌊α min l j,i ⌋ -1) mod g i,j ≤ g i,j -⌊α min l i,j ⌋ -⌊α min l j,i ⌋ -2 (37)
Therefore we can compute r = g i,j -(x -t j -⌊α min l j,i ⌋ -1) mod g i,j . If r < ⌊α min l i,j ⌋+⌊α min l j,i ⌋+2, the constraint is violated. In this case, we compute the smallest offset x ′ strictly greater than the current one, and which satisfies the current constraint, i.e. we compute x ′ > x such that (x ′ -t j -⌊α min l j,i ⌋-1) mod g i,j = 0. By proposition 2, we have x ′ = x + r. For this offset x ′ , the current constraint is now verified. We set x = x ′ and continue the process with the next constraint. We cycle along the constraints, until no update is made during N p i successive rounds or the offset x becomes greater or equal than x end . In the former case, all the constraints are satisfied by the current offset, otherwise there is no solution and we return the special value ∅. This procedure is summarized in algorithm 3. It runs in O(N p i n ′ poly ) where n ′ poly is the number of polyhedra on the interval [x start , x end ) 6 . In the worst case, all the tasks are on the same processor and we traverse the whole period, hence a complexity in O(N n poly ) which is bounded by O(N T i ).

Algorithm 3 A propagation procedure to find a feasible integral solution

Solving the best response problem

We are now able to describe procedure which solves (BO p i ). We saw that t i can be supposed to belong to [0, T p i -1]. More generally we can start at any initial offset x start , for example the current value of t i , and we run on the right until we reach the offset x end = x start + T p i . We can compute the local polyhedron (on the right) which contains the current offset. Using the primal simplex method, we solve the associated problem (Loc-BO p i ). We obtain a new local optimum (x ref , α min ). We then use the propagation procedure to reach the next improving solution. We are in a new polyhedron and we restart the local optimization at this point, which gives us a better value. We continue until the propagation mechanism reaches x end (see Figure 11). In the end, we obtain a best offset t i ∈ [x start , x end ). If needed, we can consider t i mod T p i which is an equivalent solution in [0, T p i -1]. If we use Megiddo algorithm to solve the local problems, or if we use the primal simplex approach in the case of processing times, this procedure runs in O(N n poly ). If we use the primal simplex approach, the local optimizations run in O(N 2 ), however most of the polyhedra are skipped by the propagation mechanism which runs in O(N n poly ).

Results

We test the method on non-harmonic instances generated using the procedure described in [START_REF] Eisenbrand | Solving an Avionics Real-Time Scheduling Problem by Advanced IP-Methods[END_REF] : the periods are choosen in the set {2 x 3 y 50 | x ∈ [0, 4], y ∈ [0, 3]} and the processing times are generated following an exponential distribution and averaging at about 20% of the period of the task.

Table 1 presents the results on 15 instances with N = 20 tasks and P = 4 processors. Columns 8-10 contain the results of our new version of the heuristic. The value start 10s represents the number of time the heuristic was launched with different initial solutions during 10s. The value start sol represents the number of starts needed to obtain the best result. Hence, the quantity time sol = 10start sol /start 10s gives approximately the time needed to obtain this best result. Column time sol of the MILP formulation (columns 2-3) represents the time needed by the Gurobi solver to obtain the best solution during a period of 200s 7 . In addition to being much faster, the heuristic sometimes obtains better results (see instances 1 and This table also includes the results of the original heuristic presented in [START_REF] Sheikh | Resource allocation in hard real-time avionic systems -Scheduling and routing problems[END_REF][START_REF] Sheikh | Strictly periodic scheduling in IMA-based architectures[END_REF] (columns 4-7). In fact, we test the same instances that were used by these authors. However, the way they measure performances differs from ours. For their results, the heuristic was started several times with different initial solutions until a bayesian test was satisfied, which gives a value time stop . We abandoned this bayesian test because the results on column 4 show that on a lot of instances, the process stops with a solution which is far from the best solution found by the MILP. However, since the value time single represents the time needed for a single run of their version of the heuristic, we can compute a value starts 10s = 10/time single which measures the average number of starts performed by the original heuristic in 10s. Compared with the equivalent value start 10s (column 10), we see that our version of the heuristic is incomparably faster (about 3200 times on these instances).

MILP (200s)

Original heuristic [START_REF] Sheikh | Strictly periodic scheduling in IMA-based architectures[END_REF] (bayesian test) New heuristic (10s) id α M I L P t i m e s o l α h e u r i s t i c t i m e s t o p t i m e s i n g l e s t a r t s 1 0 s α h e u r i s t i c s t a r t s s o l s t a r t s These good results have encouraged us to perform additional tests on big instances (50 processors, 1000 tasks). Table 2 presents the results for 10 instances, where starts 1000s is the number of runs performed during 1000s, starts sol is the round during which the best solution was found, time sol is the corresponding time, and time single is the average time for one run. This shows that our heuristic can give feasible solutions (α ≥ 1) in about 1min, while these instances cannot even be loaded by the MILP solver. In order to evaluate the contribution of the propagation mechanism to the solving process, we also present results where the propagation has been replaced by a simpler mechanism : once we are at a local optimum, we follow the decreasing line active at this point, until we reach the x-axis; this gives us next reference offset and therefore a next polyhedron. While the impact of the propagation is quite small in the case of small instances, we see on these big instances that the propagation mechanism accelerates the process by a factor of 37 (average ratio of the two time single values).

With propagation

Without propagation id α h e u r i s t i c s t a r t s s o l s t a r t s 1 0 0 0 s t i m e s o l t i m e s i n g l e α h e u r i s t i c t s s o l s t a r t s

Conclusion

In this paper, we have proposed an enhanced version of a heuristic, first presented in [START_REF] Sheikh | Resource allocation in hard real-time avionic systems -Scheduling and routing problems[END_REF][START_REF] Sheikh | Strictly periodic scheduling in IMA-based architectures[END_REF], and allowing to solve a NP-hard strictly periodic scheduling problem. More specifically, we present an efficient way to solve the best-response problem. This solution procedure alternates between local optimizations and an efficient propagation mechanism which allows to skip most of the polyhedra. The results show that the new heuristic greatly improves the original one and compares favorably with MILP solutions. In particular, it can handle instances out of reach of the MILP formulation.

A Appendix Proposition 1. If x mod a = 0, then (-x) mod a = 0, otherwise (-x) mod a = a -x mod a. In particular x mod a > 0 ⇔ (-x) mod a > 0.

Fig. 1 .

 1 Fig. 1. N = 4 non-overlapping periodic tasks on P = 1 processor

Fig. 3 .

 3 Fig. 3. Possible values for (ti, α) constrained by all tasks j ∈ G(i)

Fig. 4 .Fig. 5 .

 45 Fig. 4. Selection of the polyhedron containing a reference offset t * i

Fig. 6 .Fig. 7 .

 67 Fig. 6. The lowest intersection point of a fixed increasing line with decreasing lines

Fig. 8 .Fig. 9 .

 89 Fig. 8. Finding the fractional optimum with a primal simplex approach

Fig. 10 .

 10 Fig. 10. All the decreasing lines are parallel in the case of processing times

Fig. 11 .

 11 Fig. 11. Propagating up to the next strictly improving offset

Proposition 2 . ( 1 )

 21 The smallest y ≥ a such that (y -b) mod c = 0 is given by y = a+(b-a) mod c. (2) The smallest y > a such that (y -b) mod c = 0 is given by y = a + c -(a -b) mod c. (3) The largest y ≤ a such that (y -b) mod c = 0 is given by y = a -(a -b) mod c. (4) The largest y < a such that (y -b) mod c = 0 is given by y = a -c + (b -a) mod c.

Table 1 .

 1 Results of the MILP, the heuristic of[START_REF] Sheikh | Resource allocation in hard real-time avionic systems -Scheduling and routing problems[END_REF], and the new version of the heuristic on instances with P = 4 processors and N = 20 tasks

	1 0 s	t i m	e s o l

Table 2 .

 2 Results of the new version of the heuristic on instances with P = 50 processors and N = 1000 tasks, with timeout = 1000s

	1 0 0 0 s	t i m	e s o l	t i m	e s i n g l e

Note that this hypothesis can always be enforced. If li,j > 0, then the tasks i and j are constrained not to start at the same time. Hence the same is true for the tasks j and i. Therefore, lj,i can also be supposed strictly positive.

In particular, a classic disjunctive constraint tjti ≥ li,j ∨ titj ≥ li,j can be represented by the two constraints (tjti) mod T ≥ li,j and (titj) mod T ≥ lj,i for a large enough period T .

Note that if G was not symetric, we would need to consider some degenerate cases were one of the slope (increasing or decreasing) would be infinite (since li,j = 0 or lj,i = 0).

Note that in the computation only the α-coordinate is needed since the x-coordinate of the selected intersection point can be computed afterward by x = oj + lj,iα

Let us define a phase to be N p i consecutive iterations and let us show that after 2 phases, the current offset x doesn't lie in the same local polyhedron. Note that after a phase, the algorithm either stops or an update has been made. Consider a constraint j updated during the second phase. If j was not updated during the first phase, then this contraint was satisfied during the first phase (hence the current offset was in a given polyhedron of Figure2), but violated during the second. Therefore, at the end of second phase, the offset has been pushed to the next polyhedron of Figure2. If j was already updated during the first phase, then the current offset has been pushed forward to the next feasible solution two times for this constraint. Therefore, the local polyhedron has also changed.

We fix a timeout of 200s because the solver almost never stop even after 1h of CPU time (except for instances 5 and 9 for which optimality has been proved in 7s and 20s respectively).

funded by the French Midi-Pyrenee region (allocation de recherche post-doctorant n • 11050523) and the LAAS-CNRS


M Perrier

J Dumon -

J Dumon

Retours C Canudas

D Simon

A Kibangoo

A Seuret

J Opderbecke

1.

Ifremer -Various operational underwater vehicles : one strategy for simulation

Ifremer operates two types of underwater vehicles for which modeling and simulation is part of the development strategy: Victor (6000m depth Remotely Operated Vehicle) and AsterX or IdefX (3000m depth Autonomous Underwater Vehicles).

Beside the operation of these vehicles, Ifremer also performs development on each, for performance improvement and evolution. Vehicle control, acoustic communication, tele manipulation, are some of the vehicles key functions that clearly needs, for safe and efficient development purpose prior to real implementation and in-situ testing, to rely on software tools that models and simulates, as close to reality as possible, the concerned equipment, behavior or function.

In order to be consistent in its development procedure, and also to benefit from vehicle to another one, vehicle simulation software has been developed with the same strategy, objectives and tools. Portability, reusability are thus key objectives.

Modeling and simulation can concern only one specific equipment of a vehicle (for example: manipulator arm, thruster, propeller, cable or tether, acoustic modem,…) or the entire vehicle itself, and in particular the complete dynamic model of the vehicle, for control design and mission programming studies for instance.

Operating marine vehicles as fleet of individuals aiming and acting to achieve a common objective retains more and more the attention of marine vehicles operators. Several applications can be imagined (pollutant source seeking and tracking, bathymetry mapping of wide areas, animal habitat mapping,…) but common key functions are required in order to operate those vehicles together, safely and efficiently.

Developing a simulation platform allowing to simulate multiple vehicles within a common environment, having the capability to simulate the communication between them in function of the environment characteristics, to navigate in a common area as a group or a fleet, to sense the same environment in order to exchange, share or compare data, is clearly of high interest.

Mimosa mission programming tool

Mimosa is an integrated software tool designed to bring operation solutions to AUV mission management, in particular: mission preparation, planning and validation, real-time mission supervision, post-mission analysis and playback, as well as cruise and mission data management.

An extended version of Mimosa tool has been developed in order to handle multi-vehicles mission management. This version of the software allows to create fleets of marine vehicles, and to define mission profiles for these fleets, using the Mimosa graphical interface.

The Mimosa software has been also modified and adapted to be capable to automatically create mission files understandable by the Multi-Vehicle simulator described below. When dealing with fleets of marine vehicles, once the fleet mission scenario has been defined, Mimosa defines and produces the mission file for each vehicle that belongs to the fleet.

AsterX AUV Simulation platform

The AUV simulation platform consists in a complete set-up reproducing the hardware and software environment of the real AsterX AUV owned and operated by Ifremer (see figure).

Figure. The Aster x AUV © Ifremer

The platform is built around three main computers having different roles :

-the mission preparation and programming computer allowing the definition of the mission to be achieved by the vehicle. This computer runs the Mimosa software as described before. Mimosa generates the mission file to be run on the AUV. -the surface controller computer connected to the vehicle prior to the vehicle launch (vehicle on the deck) or when the vehicle is at the surface (communication using radio communication link) -the vehicle controller computer in charge of executing the downloaded mission in an autonomous way (during the mission execution, no more communication with the surface controller is possible). On the vehicle controller, several software are running: vehicle control (with high-and low-level control of the mission and vehicle), fault management system.

At this stage, it is important to note that for the simulation platform, the surface and vehicle controller used, are strictly identical to the computers used on the real vehicle. This point ensures first the portability of any development made in the simulation platform to the real world. Second, the development in the platform can be made no matter which vehicle is controlled by the vehicle controller, the simulated vehicle or the real one.

In addition to the vehicle control, the AUV is equipped with scientific sensors, the vehicle payload, controlled by a dedicated computer (the payload computer). Information coming from the payload computer (for example concerning the quality of sensed data, the logged data status,…) can be exploited by the vehicle controller in order to react in some situations.

In the simulation platform, the vehicle behavior is produced by a numerical model calculating the vehicle state evolution in response to actuators and propulsion commands. The produced state is then exploited by a set of sensor models producing formatted data under real sensor specification; these sensor data are then used as feedback for the vehicle controller level in the same way than for the real vehicle. Inputs of the model are the same than on the real vehicle: propeller and fins commands. Outputs are acceleration and velocities taken from the vehicle state evolution during time in response to commanded inputs, and are used and formatted by the sensors simulator in charge of producing the same data, with the same format, than on the real vehicle.

4.

Multi-Vehicles Simulator configuration

The Multi-Vehicle Simulator is manly used as a simulator engine. It allows simulating several vehicles having different dynamics, behavior, characteristics, sensing equipment and control modes.

The simulator also implements realistic environment models, aerial and underwater, with realistic physical characteristics. These simulated environments are used for modeling the displacement environment of the vehicles for instance when considering sea currents, for communication simulation (acoustic and radio signal propagation in particular) and in particular communication between vehicles, for the sensors models (sonar beams, salinity sensors, … are modeled and produce sensing data from the environment where the vehicles evolve).

Each marine vehicle, aerial or underwater, is modeled and simulated as an object within the distributed simulator. The vehicles own specific characteristics: dynamic model, low-level control and sensors. Vehicle models can be parameterized differently from one vehicle to another one, giving thus the capability to the user to create and manipulate homogeneous or heterogeneous fleets of marine vehicles.

Each simulated vehicle, aerial or underwater, can be controlled locally using embedded controllers pre-defined in the simulator that can be enriched or not with more complex control laws or behavior management. The simulator also proposes the possibility to control simulated vehicles by external Basically, each vehicle is closed-loop controlled by a controller, internal or external to the Simulator, which sends commands to the vehicles in response to mission controller and to vehicle states evolution.

Simple, or complex behavior models, can be used for each vehicle. Like for the vehicle controllers, simple dynamic models pre-exist in the Simulator. More complex and realistic models can be added within the Simulator, in order to get closer to real vehicle behavior.

As defined before, the mission programming tool, Mimosa, developed by Ifremer, can be used for defining scenarios and mission description, in a format understandable by the Multi-Vehicle Simulator.

A 3D Graphics User Interface is also used for real-time piloting of vehicles (for instance in the case of tele-operated vehicles) as well as for real-time mission and scenario visualization.

5.

Connection between AsterX Simulator and the Multi-Vehicle Simulator

A software module called ConnectAsterX has been developed in order to allow the connection between the two simulators: the AUV Simulator and the Multi-Vehicle Simulator. This software creates the adequate socket connections between the two simulators. In addition, it is in charge of the conversion of the exchanged frames between the simulators, using the Multi-vehicle Simulator communication library protocol.

Once the two simulators are connected, a fake vehicle that needs to be tagged AsterX is created within the Multi-vehicle Simulator; this tag will be used by the ConnectAsterX software interface to associate the fake vehicle in the Multi-Vehicle Simulator and the AsterX AUV simulated beside. The fake vehicle will then get all its behavior from the AsterX simulator. As a consequence, a fleet of vehicles can be modeled in the Multi-Vehicle Simulator with the AsterX AUV as one of the vehicles.

It is then possible to manage a fleet of heterogeneous vehicles having different dynamics and different equipment, around the AsterX AUV.

6.

Hardware in the loop : integration of hardware equipment in simulation

Knowing that the AsterX simulation platform can either control the AsterX AUV model or the real AsterX AUV (it just depends on what Vehicle Controller is connected to the Surface Controller of the platform), it is then possible to manipulate at the same time and in a common scenario, virtual (simulated) vehicles and real vehicles.

In order to go further in fleet management simulation, one could also imagine to define a fleet of simulated marine vehicles, one of them being the real AsterX AUV as defined above, and to manage the communication between the real AsterX AUV and the simulated vehicles thanks to a real acoustic modem fitted on AsterX. The other modem could be installed on the pier and physically connected to the Multi-Vehicle Simulator.

Then, and because communication management is performed within the Multi-Vehicle Simulator such that simulated vehicles can send and receive data through their virtual modems, one can envisage to make communicate for data exchange the simulated vehicles and the real AUV.

1. Context of this document

9.

Connect Project

CONNECT is a project funded by the ANR (French National Research agency). The project deals with the problem of controlling multi-agent systems, i.e. systems composed of several sub-systems interconnected between them by a heterogeneous wireless communication network. In particular the project target the The full architecture of MUSim can be seen in figure 1. We can speed up this 50ms but it increases data loss between applications and especially critical control data between ConnectSim and Simulink. In MASim, there is no limitation of sample time.

5.2.2. Sample time of dynamic model of AsterX

This part of simulation is simulated by Simulink but in variable-step sample time mode (solver "ode15s"). The sampling period is therefore be smaller for a better approximation of a continuous behavior This function has been implemented and validated in and is planed to be implemented in MUSim.

12.

5.3.3. Fleet scalability for MUSim

For MUSim, any number of AUV's, ASV's and AFV can be adding in the scenario, one by one. But, because of the time performance and the time to build the scenario, it is reasonable to have a total number of vehicles around 10 and not more than 100.

5.3.4. Fleet scalability for MASim

Copy-paste can be done to add new vehicle with all characteristics and equipment. Thanks to its lighter architecture, the maximum number of vehicle is most important. Tests are needed to assess the possible limitations.

5.3.5. Using dynamic model of AsterX for AUV's instead of cinematic one

Dynamic model of AsterX is implemented in MASim. This part of the simulation is done with the variable sample time mode.

Controller will be able to use variable sample time from 50 to 300ms, synchronized with sensors if necessary. This functionality is planned to be implemented in MUSim.

12.3.6. 5.3.6. Positioning drift Using real position and calculated position is not implemented.

A positioning system with a drift can be modeled in simulink. We can first model a drift of 2m/h, different for each agent.

For MUSim, displaying 2 different positions in ConnectIHM is not possible.

5.3.7. Bathymetric sensor

In MUSim, a bathymetric sensor is planed to be implemented in each AUV.

In MASim, such a sensor is possible if there is a non flat seabed implemented. MUSim displays a richer graphics and colorful.

5.5.2. 3D display mode

MUSim displays the simulation in 3D with ConnectIHM. It is possible to display 3D in MASim but it will be with Matlab tools so the graphics will be quite poor. Non flat seabed is display based on a real seabed map from Ifremer near Nice (see figure 5).

On the surface, for the landscape around Nice, the shape of the coastline and the relief are real and the texture is generated randomly. Show this kind of graphics in MASim involve the full development of a new GUI or a direct gateway between ConnectIHM and MASim. This solution requires specific developments by Robosoft. If this type of tool is necessary, we should first conduct a small study of feasibility.

12.6. 5.6. Source management 12.6.1. 5.6.1. Static 3D shape of the source In stationary regime, the diffusion laws give an exponentially decreasing spatial distribution of concentration, represented in 2D by the following equation:

Salinity=100*(1-e^(-K*((X-xs)²+(Y-ys)²)))
Where (xs,ys) is the center of the source and K gives the distribution's spread.

We select a steady-state for the diffusion because we postulate that the transient time is over or the time of the mission is much shorter than the diffusion time.

That means that in our study case, the source is like a fresh water source (here since years) or a source of pollution that has already reached its steady-state or with a relatively low flow rate.

To represent the fact that different currents can distort the circle, we replace the circular form by the sum of three concentric ellipses, while retaining the exponential decay (see section

To extend this to a 3D shape and for reasons of simplicity, we chose to change the concentration of the center of 2D figure linearly from 100% at the source level to 0% at the surface (see section 5.6.2, equation (b)).

To represent the diffusion along z, the spread of ellipses is weighted by depth and the 2 firsts ellipses also are exponentially reduced by the depth (see section 5.6.2, equations (a) and (c)).

This description is both used by MASim and MUSim and is summarized in section 5.6.2 Figure 14.

Change this description in MUSim needs new development for Prolexia and Robosoft. Change this in MASim is possible if necessary.

We can also add a noise to the measurement of salinity (white noise with a variance equal to 1%).

12.6.2. 5.6.

Movement of the source

To get a movement of the source due to the wave, we chose to weave the vertical line formed by centers of 2D shapes. The source position is fixed end while the other end to the surface varies according to a segment whose length, direction and frequency are configurable (see equations (d) and (e) below).

Finally, the concentration distribution related to the source is described by the following algorithm: We can display the 2D map of the source to a depth parameterized (see figure 17) It can also display the full plume in 3D: a number of bubbles appear in the volume around the source whose size, color and / or transparency depends on the salinity at this point (see figure 18)

Description
For MASim, display source plume in 3D with matlab tools seems to be difficult but can be studied if necessary. In MUSim, we can set a range of frequency for each acoustic modem (emission and reception).

Each agent can also have several modems in board, so that multichannel communication, for separating Surface-AUV and AUV-AUV communications for example, can be simulated.

In MASim, such a function can be easily implemented.

12.8. f : Frequency of the signal (kHz) k : the energy spreading factor (we take k=1.5)

For a given frequency, the communication range is computed by setting the transmission power and the reception threshold.

Attenuation does not include temperature and depth.

Online adaptation of the transmit power in order to change the communication range during the simulation can be studied if necessary.

Limited communication range is not implemented in MASim but it can be if necessary. 12.8.10.5.8.9. Delay In MUSim, we can use a function that calculates the delay of the communication packet. This delay is applied to the packet and wrote in order to use it later in simulink.

This delay is the sum of propagation delay and emission delay (=packet size/bitrate).

To define packet size, we choose to consider the transmission of an array of numbers and fixe to 32bits the definition of each number. Bit rate is set for each modem.

For example, if we send the array [id,x,y,z,Head,S,|V|,t] (defined in section 5.8.2.), packet size will be 8*32=256bits.

If modem bitrate=480bits/s, emission delay=256/480=0.533s. If the distance between emitter and receiver is 400m, propagation delay = 400/1500=0.266s. Total delay = 0.533+0.266= 0.8s. The packet is received after 7 sample times (after 0.875s) and contains its exact delay (0.8s).

The exact delay could be used for example to test the robustness of communication algorithms to multipath. We can therefore consider that air communications are instantaneous, lossless and with unlimited bandwidth.

12.8.13.5.8.12. Underwater radio frequency modems To communicate between AUV, only acoustic modems are implemented in MASim and MUSim.

Integrates characteristics of underwater radio frequency modems could be interesting.

12.8.14.5.8.13. Limitations of the acoustic propagation model The equation of attenuation of an acoustic underwater signal (section 5.8.3) does not use temperature, salinity and depth. Data detection using baseband acoustic signals will also not be implemented.

12.9. 5.9. Communication display 12.9.1. 5.9.1. Display coverage area MUSim can display communication area related the limited distance (section 5.8.3). When agent noises is used to limit the communication area(section 5.8.7.), the area display is the maximum one (with only its own noise). Display this area of communication in MASim is not implemented and is possible if necessary.

12.9.2. 5.9.2. Display any 3D cone of coverage MUSim can display communication area related the limited opening angle for reception (section 5.8.4). This function is programmed in MUSim but not tested/validated yet. This is not possible to do this in MASim.

12.9.3. 5.9.3. Display connectivity graph In MUSim, if the function is activated, each successfully received packet produce a green line connecting transmitter and receiver and a green spot is traveling along this line from one to the other. Moreover, each unsuccessfully received packet produce a red line connecting transmitter and receiver.

FeedNetBack page 7/ 36

 36 Description of the multi vehicle simulatorIMN/SM/PRAO/11.091 indice : A du 25/03/11

FeedNetBack page 8/ 36

 36 Description of the multi vehicle simulatorIMN/SM/PRAO/11.091 indice : A du 25/03/11

Figure X describes the

  Figure X describes the architecture of the AUV simulator.

Figure .

 . Figure . Architecture of the AUV simulation platform

  Figure X describes the architecture of the Multi-Vehicle Simulator as described below.

Figure .

 . Figure . Distributed Multi-Vehicle Simulator and other software modules

Figure

  Figure . Software architecture of the simulation platform

FeedNetBack page 15/ 36

 36 Description of the multi vehicle simulator IMN/SM/PRAO/11.091 indice : A du 25/03/11 S Si im mu ul li in nk k : : C Co on nt tr ro ol l A Al lg go or ri it th hm ms s R Re ea al l A As st te er rX X o or r R RT T I If fr re em me er r' 's s s si im mu ul la at to or r M Mi im mo os sa a C Co on nn ne ec ct tI IH HM M Graphical interface Environment and Vehicle simulator C Co on nn ne ec ct tS Si im m © Ifremer

Figure 1 :

 1 Figure 1: MUSim architecture

  Both MASim and MUSim can display the simulation in 2D from above.

Figure 3 :Figure 4 :

 34 Figure 3: 2D MUSim display with ConnectIHM Figure 2: 2D MASim display with Matlab tools

Figure 5 :Figure 7 :

 57 Figure 5: 3D MUSim display with ConnectIHM

FigureFigure 12 :

 12 Figure 11: AUV MiniSub in ConnectIHM Figure 10: AUV AsterX in ConnectIHM

7 . 1 . 5 . 7 . 1 .

 71571 of the multi vehicle simulator IMN/SM/PRAO/11.091 indice : A du 25/03/11 This movement description is implemented in MUSim and has to be validated. It is already validated in . 12.7. 5.7. Source display 12.Display source in 2D Both MASim and MUSim can display the source distribution in 2D from above (See figures 15 and 16).

Figure 14 :12. 7 . 2 . 5 . 7 . 2 .

 1472572 Figure 14: Summary of the source

12. 7

 7 .3. 5.7.3. Source movement Source movement is displayed in MUSim and it is possible to do it in MASim with Matlab tools if necessary.

Figure 17 :12. 7 . 4 . 5 . 7 . 4 .

 1774574 Figure 17: 2D source in 3DFigure18: source plume in 3D

Figure 19 : 1 .

 191 Figure 19: source plume in 3D in movement

12. 8

 8 .3. 5.8.3. Limited range of communicationIn MUSim, one of the management modes of communication consists in defining a transmission power, calculating the corresponding attenuation, and then comparing the attenuated signal to a pre-defined threshold associated with the modem receiver. Below the threshold the signal is considered as lost. At a distance X (in Km), the attenuation A(x) of an acoustic underwater signal of frequency f is given by A(x) = x k 10 x.g ( f ) /10 g(f) (in dB/Km) denoting the absorption coefficient for the frequency f. It can be approximated by the Thorp's formula as follows: FeedNetBack page 30/36 Description of the multi vehicle simulator IMN/SM/PRAO/11.091 indice : A du 25/03/11

FeedNetBack page 34/ 36

 36 Description of the multi vehicle simulatorIMN/SM/PRAO/11.091 indice : A du 25/03/11

Figure 21 :

 21 Figure 21: Communication range display

Figure 22 :

 22 Figure 22: Communication range with opening of 110°

Figure 23 :

 23 Figure 23: Green line and spot for each successfully received packets

Description of the multi vehicle simulator IMN

  , linked to the Multi-Vehicle Simulator thanks to a dedicated communication library. For instance, a real vehicle controller can be connected to the simulator for instance for testing a development prior to real testing on a real vehicle. On-going control developments can also be tested by connecting development tool like Matlab to the Simulator.

	FeedNetBack			page 9/36
			Surface Control	Mission
			Computer (SCC)	Programming
			(QNX)	Tool (MIMOSA)
			Download	mission
	Vehicle Control Computer (VCC)	Rpm,	actuators	AUV Model (Matlab/Windows
	(QNX PC)			PC)
			AUV state evolution
		Sensor data	AUV sensors Model
				(Windows PC)
	Ethernet	NMEA

/SM/PRAO/11.091 indice : A du 25/03/11 controllers

  This project started in May 2007 and will finish in February 2011.

		FeedNetBack	page 14/36
	Partners are NeCS INRIA-GIPSA-lab research team (project leader), Ifremer,
	Prolexia and Robosoft-PGES.
	More information in the website http://www.lag.ensieg.inpg.fr/connect/
	10.	3. MUSim and MASim architecture
	MUSim and MASim have been jointly developed and can now be used either
	independently or in a complementary manner.

Description of the multi vehicle simulator IMN/SM/PRAO/11.091 indice : A du 25/03/11 control of a cluster of agents composed of autonomous underwater vehicles and marine surface vessels. MUSim means Multi-vehicle Full Simulator and MASim means Multi-vehicle Fast Simulator. MASim can be seen as a light version of MUSim that offers less functionality. Scenarios and control algorithms can be tested with MASim and then be implanted in MUSim through a reduced work. MASim is a simple Simulink model and can be implemented very easily from a PC to another. As its name suggests, its lightness also allows it to be faster than MUSim. MUSim is mainly based on the Simulator ConnectSim communicating with a Human Machine interface ConnectIHM and Simulink model. It also can be connected to the operational mission planning tool Mimosa. Communication between theses applications is Ethernet UDP based. MASim and the Simulink part of MUSim are developed by Gipsa-Lab, ConnectSim and Mimosa by Prolexia and ConnectIHM by Robosoft-PGES. An interface between MUSim and a real AsterX has been developed by Ifremer. As all Connect partners are French, ConnectSim and Connect IHM are in French.

Table 1 : Summary of functionalities

 1 Movement: the vertical line formed by centers of 2D shapes can oscillate. The source position is fixed end while the other end to the surface varies according to a segment whose length, direction and frequency are configurable

	Functionalities

Description of the multi vehicle simulator IMN/SM/PRAO/11.091 indice : A du 25/03/11 6.2

Detailed functionalities 12.1. 5.1. System requirements 12

	FeedNetBack	page 18/36
	12.1.3. 5.1.3. Windows XP
	In MUSim, ConnectIHM needs windows XP.
	12.1.4. 5.1.4. Hardware
	For MUSim, to run properly ConnectSim, ConnectIHM and Simulink together, hardware
	needs good performances. Cache of CPU and graphical NVIDIA card is particularly
	important.
	For example, MUSim is installed on a tablet PC Dell Precision M6400 with:
	-Intel® Core™2 Extreme Mobile Processor QX9300 (12MB of cache)
	-Graphical card NVIDIA quadro fx 3700M (1GB)
	It is also possible to install ConnectSim, ConnectIHM and Matlab/Simulink in separates
	Computer to improve performances because the communication between these
	applications is Ethernet UDP based.

.1.1. 5.1.1. Matlab-Simulink For both MUSim and MASim, Matlab-Simulink R2007b is needed to run algorithms. Several "S function" are used so to modify and compile this part of the simulator, Visual C is required. Description of the multi vehicle simulator IMN/SM/PRAO/11.091 indice : A du 25/03/11 12.1.2. 5.1.2. Internet access In MUSim, ConnectSim installation needs internet access. 12.1.5. 5.1.5. Ownership and exploitation rights MUSim is co-owned by the consortium of Connect and is restricted to internal and research use. MASim is owned by Gipsa-Lab and has no specific limitation. 12.2. 5.2. Time characteristics 12.2.1. 5.2.1. Sample time In MUSim, typical sample time are 125ms or 100ms (exception of the dynamic model of AsterX).

3. Vehicles and fleet 12.3.1. 5.3.1. Energy vehicles management

  If necessary, a model of energy consumption of a vehicle can be implemented in Simulink and include to the simulator.

	FeedNetBack	page 19/36
	12.3.2. 5.3.2. Computation vehicles management
	If necessary, a model of computation limitation of a vehicle can be implemented in Simulink
	and include to the simulator.
	2.3. 5.2.3. Simulation speed
	MUSim cannot simulates faster than real time with the typical scenario of 5 AUV, using
	cinematic model for AsterX and Ts=125ms.
	MASim can simulate 10 times faster with the same scenario.

Description of the multi vehicle simulator IMN/SM/PRAO/11.091 indice : A du 25/03/11 12.3. 5.

4. 5.4. Environment 12.4.1. 5.4.1. 3D map of the non flat seabed

  3D map of the non flat seabed is implemented in MUSim. This map, owned by Ifremer, describes the seabed of a real area near Nice. This map is used by ConnectSim for bathymetric sensor and to calculate reflection of acoustic signals (see4.8.7.) This map is only used by ConnectIHM to display it. The implementation of this map in MASim has to be studied if necessary. One current and one wind fixed and uniform is implemented in ConnectSim. This can be defined also in MASim if necessary. Current is planed to affect only the dynamic model of AsterX. Plume of the source is defined by an equation which already includes movement. Using the current to affect the plume of the source would change totally the source definition. If necessary, this could be done only for MASim and has to be discussed with Prolexia and Robosoft for MUSim.Wind doesn't affect anything. To use Wind, dynamic model of ASV and AFV has to be implemented in Simulink.

	FeedNetBack	page 20/36
	12.4.2. 5.4.2. Current and wind
	12.3.8. 5.3.8. Collision Management
	No collision management is implemented.
	Such functionalities seem to be not necessary and hard to implement in simulink.

Description of the multi vehicle simulator IMN/SM/PRAO/11.091 indice : A du 25/03/11 12.

  In MUSim, the modem can be fixed on the AUV at any position and orientation. We can also define its opening angle for reception. Combined with the range communication, it yields a 3D cone of coverage. This function is programmed but not yet tested/validated yet. Programming such a function in Simulink for MASim seems to be difficult. Interference if 2 signals are emitted at the same frequency in the same time is planned to be done for MUSim. This will be done in simulink for allowing an implementation in MASim if necessary.For trajectories 3,4 and 5, depth of reflection is taken from the non flat seabed map at the middle of the trajectory. From a geometrical point of view, such a calculation is correct for the third trajectory. It corresponds to a somewhat valid approximation for the remaining cases. For one emission, if this function is activated, another modem can potentially receive 5 identical frames to different power levels (because of the potential attenuation of reflection and the fact that the length of each trip is different).

	12.8.4. 5.8.4. 3D cone of communication coverage
	12.8.6. 5.8.6 Management of acoustic waves reflection
	Wave reflections induce additional signal attenuation.
	Four different trajectories are calculated for the acoustic wave with 1 and 2 reflections:
	1-Direct path
	2-Surface reflection (no attenuation)
	3-Seabed reflection (-3dB attenuation)
	4-Surface + seabed reflection (-3dB attenuation)
	5-Seabed + surface reflection (-3dB attenuation)
	This function will be used to test the robustness of communication algorithms to multipath.
	This function is programmed in MUSim but not yet tested/validated yet.
	Programming such a function in Simulink for MASim can be done if necessary but it needs
	to define a non-flat seabed.

12.8.5. 5.8.5. Interference at the same frequency This function is programmed but not tested/validated yet. Programming such a function in Simulink for MASim can be done (maybe in a simplified version) if necessary.

11.5.8.10. Shadowing zone for communication

  This function is programmed in MUSim but not tested/validated yet. Programming such a function in Simulink for MASim can be done if necessary. Shadowing zone for communication because of a vehicle or a rock in the direct path is not considered. If necessary, this possibility has to be discussed with Prolexia. This function is programmed in MUSim but not tested/validated yet.As simulation of upper layers like WiFi 802.11b/g/n standards can be used, it seems to be useless to develop this part. Moreover, flows to the conventional types of communication are very large compared to those available in underwater communication.

	FeedNetBack	page 33/36
	12.8.12.5.8.11. Aerial communication management
	In MUSim, aerial electromagnetic waves management is implemented with the following
	equation.
	FSL[dB]=C+20 log(D)+20 log (F)
	D= distance
	F=frequency (MHz)

Description of the multi vehicle simulator IMN/SM/PRAO/11.091 indice : A du 25/03/11 12.8.