Florent Altché

Xiangjun Qian

Arnaud De

La Fortelle

Time-optimal Coordination of Mobile Robots along Specified Paths

In this paper, we address the problem of timeoptimal coordination of mobile robots under kinodynamic constraints along specified paths. We propose a novel approach based on time discretization that leads to a mixed-integer linear programming (MILP) formulation. This problem can be solved using general-purpose MILP solvers in a reasonable time, resulting in a resolution-optimal solution. Moreover, unlike previous work found in the literature, our formulation allows an exact linear modeling (up to the discretization resolution) of second-order dynamic constraints. Extensive simulations are performed to demonstrate the effectiveness of our approach.

I. INTRODUCTION

The deployment of autonomous mobile robots is expected to bring major benefits in many applications, and their number is likely to grow dramatically in the next decades. Therefore, the need to coordinate these robots, which means finding a way for each of them to reach its target without colliding with another robot, will become increasingly important. A vast literature on motion planning, which generalizes the problem of coordination, already exists (see e.g. [START_REF] Latombe | Robot motion planning[END_REF]).

We consider the problem of optimally coordinating multiple robots along specified paths with variable speed under kinodynamic constraints. The simplifying fixed path assumption is notably suited for structured environments, for instance traffic intersections or warehouses, where robots are generally bound to navigate inside lanes or aisles. Various methods exist to quickly find a feasible solution to this problem, and many of these (see e.g. [START_REF] Siméon | Path coordination for multiple mobile robots: A resolution-complete algorithm[END_REF]) make use of the so-called coordination (or configuration) space introduced in [START_REF] Lozano-Perez | Spatial planning: A configuration space approach[END_REF]. However, the optimal coordination problem is known to be NP-hard [START_REF] Dasler | On the complexity of an unregulated traffic crossing[END_REF], even in the absence of dynamic constraints, and few methods exist to provide a good solution with a high number of robots. As shown in [START_REF] Gregoire | Priority-based coordination of mobile robots[END_REF], this algorithmic complexity stems from the implicit decision of choosing which of any two potentially colliding robots should be the first to pass, making the problem inherently combinatorial with complexity scaling as high as 2 N (N -1)/2 for N robots. Exhaustive enumeration could be used for small instances (see e.g. [START_REF] Chitsaz | Exact Pareto-optimal coordination of two translating polygonal robots on an acyclic roadmap[END_REF], [START_REF] Murgovski | Convex modeling of conflict resolution at traffic intersections[END_REF]), but would be difficult to scale up with a larger number of robots.

A possible way of handling this complexity is to prune the corresponding decision tree by removing provably nonoptimal branches. Mixed-integer programming (MIP) is a 60 Bd St Michel 75006 Paris, France [florent.altche, xiangjun.qian, arnaud.de_la_fortelle] @mines-paristech.fr widely-used framework that allows efficient handling of such combinatorial problems. General MIP problems involving arbitrary functions are very hard, but good techniques exist for a subclass of these problems, called mixed-integer second-order cone programming, or MISOCP [START_REF] Alizadeh | Second-Order Cone Programming[END_REF]. In these problems, a convex quadratic objective function is minimized with quadratic positive semi-definite or affine constraints. A better-known subclass of MISOCP is mixed-integer linear programming (MILP), where the objective function and constraints are linear. These techniques have already been applied to trajectory planning in general, and to the coordination problem in particular.

In [START_REF] Richards | Spacecraft Trajectory Planning with Avoidance Constraints Using Mixed-Integer Linear Programming[END_REF], a MILP formulation is used to compute fueloptimal trajectories for multiple spacecrafts, while avoiding collisions and exhaust gases from other crafts. However, this model is inherently different from ours as paths are not specified in advance in this formulation, leading to a much higher computational complexity. On the opposite end of the spectrum, Wang et al. [START_REF] Wang | Optimal trajectory planning for trains using mixed integer linear programming[END_REF] use MILP to find an optimal velocity profile under piecewise-constant dynamic constraints for a train on a fixed path, but do not consider conflicts with other trains.

In [START_REF] Cafieri | Aircraft deconfliction with speed regulation: New models from mixed-integer optimization[END_REF], a mixed-integer nonlinear program is used to ensure safe separation of aircrafts evolving along fixed paths, with minimal deviation from an original flight plan. Even though dynamic constraints are not considered, nonlinearities stemming from the euclidean distance constraints render the problem very hard to solve in reasonable time, with a practical limit of 4 aircrafts.

Peng and Akella [START_REF] Peng | Coordinating Multiple Robots with Kinodynamic Constraints Along Specified Paths[END_REF] consider a problem almost identical to ours. First, they propose to discretize robot paths into "collision-free" and "conflicting" segments. In a second step, by identifying the time instants when robots can enter and exit collision segments, they formulate collision avoidance constraints for every pair of robots into a nonlinear mixedinteger problem. However, the dynamic constraints are nonconvex in this formulation, rendering the problem unsuitable for general-purpose solvers. To handle this difficulty, additional constraints are introduced to transform the problem into two linear subproblems that are solved as MILPs and give an upper and lower bound for the solution, but the actual optimal value is only approximated.

The problem of non-convexity of the dynamic constraints is also encountered in [START_REF] Murgovski | Convex modeling of conflict resolution at traffic intersections[END_REF], where authors propose a convexification method using linearization along a reference speed. However, this approximation underestimates the maximum and overestimates the minimum acceleration, leading to a sub-optimal solution.

Our main contribution in this paper is the introduction of a new formulation of the optimal coordination problem as a MILP by discretizing over time instead of space, and using bounding polygons to transcribe safety requirements as linear constraints. This formulation allows an exact linear modeling of second-order dynamics up to the discretization resolution; solving this problem gives optimal velocity profiles for the robots under these constraints. For illustration purposes, we use the mean sojourn time as the minimization objective, but other linear optimization criteria could also be chosen.

The paper is articulated as follows. In Section II, we describe the modeling of the system of robots and the general coordination problem, and we give some theoretical insights on the problem from previous work in Section III. In Section IV, we present our formulation of the time-optimal coordination problem as a MILP problem. This formulation is then validated using computer simulation on the example of automated vehicles in an intersection, as presented in Section V. Finally, Section VI concludes the study.

II. PROBLEM STATEMENT

We consider a set N of N robots evolving on predetermined paths, and we assume that coordination between robots is only needed inside a bounded region, which we call the coordination region. In the case of automated driving, for instance, coordination is mostly needed in the middle of road intersections, as illustrated in Fig. 1, while vehicles only need to keep a safe distance with the vehicle in front of them when they are not in the intersection. In this example, the coordination region would be chosen as the center of the intersection, including a portion of the roads leading to, and exiting from this center. Note that the predefined path assumption allows to only consider longitudinal movement of the robots, thus reducing the computational complexity, and is classical in robots coordination problems in a constrained environment [START_REF] Peng | Coordinating Multiple Robots with Kinodynamic Constraints Along Specified Paths[END_REF], [START_REF] Pecora | On mission-dependent coordination of multiple vehicles under spatial and temporal constraints[END_REF].

Each robot i ∈ N is supposed to follow a predetermined path γ i inside the coordination region, so that we need only consider the longitudinal comportment of the robots. The dynamics of robot i are described as a double integrator:

( ṡi , vi ) = (v i , a i ) (1)
where s i is the curvilinear position of robot i along its path, v i its longitudinal velocity and a i its acceleration. The origin of s i is chosen so that s i = 0 when the front of the robot enters the coordination region, and s i = s out i > 0 when it fully exits the coordination region. s i can therefore be interpreted as the distance traveled by robot i inside the coordination region. The velocity is assumed to be nonnegative and bounded, such that v i ∈ [0, v i ]. The state of robot i is noted x i = (s i , v i ); we call trajectory of robot i the (continuous) function mapping a given time t to the state of i at time t, noted x i (t) = (s i (t), v i (t)). Boldface x denotes the vector (x 1 , x 2 , ...), representing the state of the multi-robot system and x(t) is the system trajectory. To account for the dynamic constraints on the robots, we assume that the longitudinal acceleration a i of robot i is bounded to an interval a i , a i with a i < 0 < a i .

We assume that robot i enters the coordination region at time t in i with speed

v in i ∈ [0, v i ],
and is required to leave the coordination region with speed v out i ; we let t out i denote the corresponding exit time. Note that v out i should be properly chosen to avoid collisions outside of the coordination region, for instance when a fast vehicle exits after a slower one.

For a pair of distinct robots i and j, we call collision set between i and j, noted C ij ⊂ [0, s out i ] × 0, s out j the set of positions (s i , s j ) inside the coordination region where i and j would collide; we say that i and j are conflicting if

C ij = ∅.
To simplify the rest of the presentation, we assume that the collision set between two robots is connected, i.e. we exclude the case of non-conflicting segments between two conflict segments, for instance when two paths intersect multiple times. If this is not the case, the presented results still hold provided every connected component C p ij of C ij is considered individually. Moreover, we approximate the exact collision set by a (minimal) bounding polygon with edges either parallel to the horizontal or vertical axis, or to the s i = s j line. Under our hypotheses, four possible types of conflicts can exist: robots can follow one another throughout the coordination region (a), or have crossing (b), merging (c) or diverging (d) paths, as illustrated in Fig. 2. Note that cases (c) and (d) are not mutually exclusive, and paths can merge and then diverge.

s i s j C ij (a) Following paths s i s j C ij (b) Crossing paths s i s j C ij (c) Merging paths s i s j C ij (d) Diverging paths Fig. 2.
Possible cases for the collision set between two conflicting rectangular robots i and j (assuming paths intersect orthogonally). The edges of C ij are straight lines with equations

s i = A, s j = B or |s i -s j | = C where A, B, C are constants.
The above bounding polygon approximation is well suited when robots have nearly rectangular shapes, and induces a negligible loss of optimality in this case, as illustrated in Fig. 3. Note that arbitrary convex polygons can also be used, at the cost of introducing additional binary variables. We now define the time-optimal coordination problem as: Definition 1: The time-optimal coordination problem is that of finding the optimal system trajectory x * (t) minimizing the mean sojourn time

1 N i∈N (t out i -t in i )
, under the following constraints for every robot i:

• initial conditions: x i (t in i ) = (0, v in i ) • kinematics: ṡi = v i ∈ [0, v i ] • dynamics: vi = a i ∈ [a i , a i ] • safety: for all t ≥ t in i and for every robot j = i, (s i (t), s j (t)) / ∈ C ij • liveness: there exists t out i < +∞ with s i (t out i ) = s out i • exit speed: v i (t out i ) = v out i .
Note that in the above definition, the initial condition x i (t in i ) = (0, v in i ) corresponding to a fixed entry time can be replaced by a fixed initial state constraint x i (0) = (s 0 i , v 0 i ). III. THEORETICAL ANALYSIS For an arbitrary pair of robots with non-empty collision set, one necessarily passes before the other to avoid collisions. For two conflicting robots i and j and a given collision-free system trajectory, we say that i has priority over j if i goes before j, and we note i j in this case. The collection of all priorities for all pairs of robots can be encoded as a priority graph where robots are the nodes and priorities are directed edges. It has been shown in [START_REF] Gregoire | Priority-based coordination of mobile robots[END_REF] that such priority graphs can be bijectively mapped to homotopy classes of trajectories for the multi-robot system. For a given continuous optimization criterion, there exists at least one optimal trajectory in any non-empty homotopy class. The global optimal trajectory x * (t) can then be found by enumerating all optimal trajectories for all priority graphs.

Therefore, the time-optimal coordination problem has both a discrete (enumerating all priority graphs) and a continuous part (optimizing the trajectories of every robots under assigned priorities, which is a continuous optimal control problem). The discrete part of the problem is combinatorial since, for N robots inducing p pairs of conflicts, there are up to 2 p possible priority graphs. In general, p can be as high as 1 2 N (N -1); however, many of these graphs can be discarded for poor performance or for being incompatible with the robots entry times. For this reason, branch-and-bound algorithms seem particularly well-suited to our problem as they are designed to find global optima without needing to explore the whole decision tree.

IV. TIME-OPTIMAL COORDINATION

Using the above theoretical results, we formulate the time-optimal coordination problem as a mixed-integer linear program (MILP) which can be solved by widely-available solvers using branch-and-bound techniques.

A. Completed Collision Set

For conflicting robots i and j, we call completed collision set C i j the set of configurations (s i , s j ) leading to a collision or a violation of priority i j; mathematically,

C i j = C ij + (R -× R + ).
Since we have approximated each collision set C ij by a minimal bounding polygon with edges parallel to the coordinate axes or to the s i = s j lines, the completed set C i j is also a polygon with the same properties. Note that if C ij is not connected, a completed collision set has to be defined for each connected component.

To formulate safety requirements as linear constraints, we first define a partition of C i j as C i j ∪ C ⊥ i j . C i j is the subset of C i j with boundary parallel to s i = s j , and C ⊥ i j the subset with boundary parallel to the s i axis as illustrated for the merging case in Fig. 4. Subset C i j is the set of priority violations (or collisions) that could happen when j should follow i, and C ⊥ i j is that of violations (or collisions) that could happen when j should wait for i to pass. In what follows, we note S ij (•) and S ij (•) the lower and upper bounds of the projection of C i j on the s • axis (for • = i or j), and we define S ⊥ ij (•) and S ⊥ ij (•) similarly. However, to ensure those subsets have empty intersection, we exclude

S ⊥ ij (i), s j : s j ∈ [0, s out j ] from C ⊥ i j . If either C i j or C ⊥
i j is empty, we let the corresponding lower and upper bounds be equal to 0; note that in this particular case, those subsets formally do not form a partition.

An important remark is that C i j and C ⊥ i j are both connected and left invariant by translation of a vector from {0} × R + . Therefore, with this decomposition, if S ij (i) ≤

s i ≤ S ij (i), condition (s i , s j ) / ∈ C i j is equivalent to s j ≤ s i -a ij
where a ij is a constant corresponding to a following distance and a potential offset of curvilinear abscissa.

If S ⊥ ij (i) ≤ s i < S ⊥ ij (i), condition (s i , s j ) / ∈ C i j is equivalent to s j ≤ S ⊥ ij (j)
. If s i is not in either of these intervals, collisions or priority violations cannot occur. Note that, if needed, a finer decomposition could be used to better approximate the exact shape of the collision set, by using bounding convex polygons with edges having slopes different than 0 or 1, and defining one subset of C i j per edge.

B. Discretization

If priorities have been chosen, finding the best trajectories in continuous time and space is equivalent to a nonlinear time-optimal control problem in a non-convex space, which ⊥ ij (i), S ij (i) and S ij (i). s j = S ⊥ ij (j) when j reaches B in Fig. 4a, and s i = S ij (i) when i reaches C. Similarly, s i = S ⊥ ji (i) when i reaches B . Also note that the equation for the lower-right boundary of the collision set is

i j B B C D S ⊥ ij = S ij (a) Merging paths s i s j C ⊥ i j C i j C i j S ⊥ ij (i) S ⊥ ij (i) = S ij (i) S ij (i) S ⊥ ij (j) (
s i ≥ s j + S ij (i) -S ⊥ ij (j).
is usually difficult to solve. Several authors, among which [START_REF] Peng | Coordinating Multiple Robots with Kinodynamic Constraints Along Specified Paths[END_REF] and [START_REF] Murgovski | Convex modeling of conflict resolution at traffic intersections[END_REF], have used spatial discretization to overcome this difficulty and formulate a simpler problem. In this setting, paths are divided into segments and the average speed v avg and occupancy time t occ for each robot in each segment are chosen as variables. One of the major issues arising from this discretization is that the second-order dynamics of the robot imply nonlinear (and non-convex) relations between the length of a segment and these two variables. Indeed, these constraints are written as = v avg t occ which translates to ≤ v avg t occ ∧ ≥ v avg t occ ; one of those inequalities defines a non-convex set in the (v avg , t occ ) plane that has to be approximated in order to be treated in most optimization frameworks, thus degrading solution quality. Instead of using spatial segments, we propose a temporal discretization which, to the best of our knowledge, has not yet been used to solve the coordination problem with a fixedpaths assumption. In what follows, we note τ > 0 the fixed duration of a time step; for an integer k ≥ 0 and a robot i, we note s k i = s i (kτ ), and use similar notations for v i . We note K the maximum number of time steps in the computation.

C. Variables

Besides variables s k i and v k i , we need to introduce a few supplementary variables, all of them binary. For two conflicting robots i and j, we let π ij = 1 if and only if (iff) i j, and for all time steps k we define several variables ε ij that indicate if, at time step k, robots i and j have entered and/or exited C ⊥ i j and C i j . Specifically, we let ε

•,in ij (i, k) = 1 iff s k i ≥ S • ij (i) and ε •,out ij (i, k) = 1 iff s k i ≥ S • ij (i)
, where • is either or ⊥, and we define similarly ε

•,in ij (j, k) = 1 iff s k j ≥ S • ij (j) and ε •,out ij (j, k) = 1 iff s k j ≥ S • ij (
j). We also introduce, for every robot, the variables µ k i = 1 iff s k i ≥ 0 and σ k i = 1 iff s k i ≥ s out i , respectively indicating if the robot has entered and exited the coordination region at time step k.

D. Objective function

As stated earlier, we wish to minimize the average sojourn time, which is equivalent (since entry times are prescribed) to minimizing the average exit times. Since our formulation does not use time as a variable, we use

O = 1 N i=1...N k=0...K σ k i (2)
for our objective function. O is the average number of time steps spent after exiting the coordination region, which approximates (to the duration of a time step τ ) the total amount of time spent after exiting the coordination region, divided by τ . More precisely,

O = 1 N i=1...N (K -k out i + 1)
where k out i is the time step at which robot i exits the coordination region, and thus maximizing O is equivalent to minimizing the average value of k out i .

E. Constraints

Many of the constraints needed for the coordination problem are conjunctions (noted ∧), or logical implications (⇒); we use the binary variables introduced in the previous section as indicators, and use a "big-M " formulation [START_REF] Richards | Mixed-integer programming for control[END_REF] to enforce those constraints. In what follows, left-hand side terms are variables, and right-hand side terms are problem parameters.

a) Binary variables: To ensure the additional binary variables are set to the correct value, the following helper conditions are enforced for all 0 ≤ k ≤ K, every robot i ∈ N and every robot j ∈ N conflicting with i:

µ k i = 0 ⇒ s k i ≤ 0 (h1)
µ k i = 1 ⇒ s k i ≥ 0 (h2)
σ k i = 0 ⇒ s k i ≤ s out i (h3)
σ k i = 1 ⇒ s k i ≥ s out i (h4) π ij +π ji = 1 (h5)
Constraints similar to (h1)-(h2) are used for ε

•,in ij (i, k), ε •,in ij (j, k), ε •,out ij (i, k) and ε •,out ij (j, k)
, where • is either ⊥ or . Constraint (h5) ensures that exactly one of the variables π ij and π ji is set to 1. Note that strict inequalities cannot be enforced in a MILP framework; therefore, the above constraints do not specify the value of each indicator variables at its point of discontinuity. This limitation, however, has little effect on the solutions.

b) Initial and terminal values, bounds: To account for the initial values of the variables and the different bounds on the problem, we use the following boundary constraints for all 0 ≤ k ≤ K -1 and each robot i:

µ k i = 0 ⇒ v k+1 i = v in i ( b1
)
σ k+1 i = 1 ⇒ v k i = v out i ( b2
)
s 0 i = -v in i t in i (b3) s K i ≥ s out i (b4) v 0 i = v in i (b5) v k i ∈ [0; v i ] (b6)
Conditions (b1), (b3) and (b5) enforce the initial constraint x i (t in i ) = (0, v in i ), (b2) the final speed constraint, (b6) enforces the bounds on speed and (b4) ensures the liveness constraint, as all robots are required to exit in finite time. c) Kinodynamic constraints: We assume that robots use a constant acceleration during each time step. Under this assumption, we enforce the kinodynamic constraints using the conditions, for all 0 ≤ k ≤ K -1 and every robot i:

σ k i = 0 ⇒ s k+1 i -s k i - 1 2 τ v k+1 i + v k i = 0 (k1) σ k i = 0 ⇒ v k+1 i -v k i ≤ a i τ (k2) σ k i = 0 ⇒ v k+1 i -v k i ≥ a i τ (k3)
Condition (k1) enforces the kinematic constraints (under the constant acceleration assumption); (k2) and (k3) account for the dynamic constraints. Note that the constant acceleration hypothesis could be relaxed using third-order dynamics, or more general second-order cone programming, as

x k i -v k i τ + a i τ a i -a i y k i - 1 2(a i -a i ) y k i 2 ≤ a i a i τ 2 2(a i -a i ) x k i -v k i τ + a i τ a i -a i y k i - 1 2(a i -a i ) y k i 2 ≥ a i a i τ 2 2(a i -a i )
where

x k i = s k+1 i -s k i and y k i = v k+1 i -v k i .
A justification for this extension can be found in [START_REF] Johnson | Optimal acceleration-bounded trajectory planning in dynamic environments along a specified path[END_REF] which shows that these constraints exactly describe the set of reachable positions and speeds under second-order integrator dynamics with bounded acceleration in a given time τ .

d) Safety constraints: Using the previously-defined indicator variables ε and the results from Section IV-A, we translate the safety constraints as follows: for every pair of conflicting robots (i, j) and for all 0 ≤ k ≤ K -1:

π ij = 1 ∧ ε ⊥,out ij (i, k) = 0 ⇒ ε ⊥,in ij (j, k + 1) = 0 (s1) π ij = 1 ∧ ε ,in ij (j, k) = 1 ∧ ε ,out ij (i, k) = 0 ⇒ s k+1 i -s k+1 j ≥ a ij (s2) π ij = 1 ∧ ε ,in ij (j, k) = 1 ∧ ε ,out ij (i, k) = 0 ⇒ s k+1 i -s k+1 j + τ 2 v k+1 i -v k+1 j ≥ a ij (s3) with a ij = d + S ij (i) -S ⊥ ij (j)
, in which d is a following distance (from front of the follower to rear of the leader) between two robots. The term S ij (i)-S ⊥ ij (j) in a ij accounts for the potential offset in curvilinear position between two robots in the case of merging paths (and vanishes in other cases), as illustrated in Fig. 4. Condition (s1) can be phrased as "if i has priority and has not yet passed the collision set, then j cannot go in". The additional condition (s3) is used to ensure that no collision occurs between two time steps; for the same reason, (s1) involves ε ⊥,out ij (i, k) and ε ⊥,in ij (j, k + 1). This approach can be seen as a systematization of the concept of spatio-temporal trajectory envelopes presented in [START_REF] Pecora | On mission-dependent coordination of multiple vehicles under spatial and temporal constraints[END_REF], [START_REF] Cirillo | Integrated Motion Planning and Coordination for Industrial Vehicles[END_REF], where the temporal extent of these envelopes is automatically adjusted during resolution.

Note that the above formulation can be used in the case of a collision set with multiple connected components, by introducing a set of variables π and ε, and of parameters S and S for each of these components. This method allows simultaneous resolution over multiple conflict areas, thus ensuring the solution is a global optimum and not comprised of several local optima, which would be the case if each connected component was treated separately.

F. Optimization problem

To simplify notations, we note X the tuple of all the variables described above. The optimization problem (in which indices have been omitted for readability) of finding

max X O(X) (3)
s.t. (h1) -(h5), (b1) -(b6), (k1) -(k3), (s1) -(s3)
either gives a solution to the discretized time-optimal coordination problem, or is infeasible. Specifically, infeasibilities can either be caused by the choice of a too small value for the number of time steps K, or because the initial states of the robots do not allow a safe passage through the intersection.

The sub-optimality caused by the time discretization vanishes as the time step goes to zero. Since the limit on computation time effectively sets a lower bound on the time step duration, it is desirable to choose a value providing good quality solutions in reasonable time. This issue will be discussed in Section V. However, an additional artifact arises from time discretization with a finite resolution, as objective function O does not distinguish solutions with sojourn times differing by less than the duration of a time step for at least one robot. To correct this issue, we make use of the fact that maximizing the speed of a robot allows to minimize its (nondiscretized) sojourn time. Therefore, adding an "averaged normalized speed" term 1

N K i,k v k i
vi to function O allows to choose, among all solutions of (3), the one with highest average speed and thus smallest average sojourn time. Note that the weighting of the added term ensures this solution is still optimal for problem (3), and the modified objective function has been used in our simulations.

Moreover, our formulation could also be used for continuous arrivals of robots in real-time. Let T be an upper bound for the computation time for N T robots: at a given time t we consider the set N t of robots entering the coordination region between times t and t+T , and we assume |N t | ≤ N T . If optimal trajectories have been assigned for the robots of N t-T and taking those as constraints for the robots of N t , we can compute optimal trajectories for those robots before they reach the entry of the coordination region. Constraint (b1) ensures those trajectories remain feasible by that time and can therefore be assigned to the robots of N t . Note that this time-receding method may, however, cause sub-optimalities compared to considering all robots at once. Also note that optimizing robots speed profiles for minimum sojourn time comes in pair with reducing safety margins to a minimum. For actual applications, this may cause problems in the case of unexpected events such as failure of a robot, which could make other robots unable to avoid a collision. Future work will study how to balance efficiency with the ability to cope with contingencies, to improve the robustness of the coordination.

V. SIMULATION RESULTS

The use of the above optimization problem to find a time-optimal coordination has been validated by computer simulation on the example of autonomous vehicles in the intersection of Fig. 1. The simulation is based on the free traffic modeling tool SUMO [START_REF] Krajzewicz | Recent development and applications of SUMO -Simulation of Urban MObility[END_REF] and uses its path generation algorithm to compute collisions sets.

Vehicles are generated either deterministically (first simulation), or using random Poisson arrival times with normallydistributed entry speeds truncated to a minimum and a maximum speed (second and third simulations). Optimization problem (3) (using the modified objective function) is then run into the commercial MILP solver Gurobi [START_REF] Gurobi Optimization | Gurobi optimizer reference manual[END_REF], using its Python interface to generate the constraints. Lastly, if the problem is feasible, the solution trajectories are simulated in SUMO using the TraCI interface to verify that they do not generate collisions.

In all simulations, vehicles are modeled as rectangles of 5 m length by 2 m width, with [a i , a i ] = [-3, +4] m s -2 . The exit speed for all i ∈ N is set as v out i = v = 15 m s -1 , which ensures the absence of collisions outside of the coordination region. The entry speed v in i is deterministically chosen in the first simulation and is normally distributed with average 12 m s -1 and standard deviation 3 m s -1 , truncated to [START_REF] Wang | Optimal trajectory planning for trains using mixed integer linear programming[END_REF][START_REF] Johnson | Optimal acceleration-bounded trajectory planning in dynamic environments along a specified path[END_REF] m s -1 in the second and third simulations.

Simulations were performed on a personal computer running on a 3.60 GHz Intel Core i7-4790 CPU with 16 GB of RAM. A replay of some of our simulations is available in the accompanying video submission 1 .

A. Microscopic simulation

We first demonstrate the ability of our method to find the global optimum on a simple example with three vehicles in the case of the intersection displayed in Fig. 1: vehicle 1 goes from south to west, vehicle 2 from west to east and vehicle 3 from north to south. Vehicles initially start at (s 0 1 , s 0 2 , s 0 3 ) = (0, 0, 25) m with speeds (v in 1 , v in 2 , v in 3 ) = (5, 15, 10) m s -1 . Fig. 5a shows the globally optimal trajectories for each vehicle, which lies in the homotopy class represented by priorities 3 2, 2 1, 3 1. For comparison purposes, Fig. 5b shows the (locally) optimal trajectories when suboptimal priorities 1 3, 3 2, 1 2 are enforced. The optimum average sojourn time found by our algorithm is 6.5 s, and the example sub-optimal one is 8.4 s.

B. Influence of time step duration

The discretization time step has a double effect on the solution: first, we assume constant acceleration during one time step. Second, the safety constraints require that one robot of each conflicting pair leaves the conflict area one time step before the other can enter.

In Fig. 6 we show the average optimality loss caused by choosing larger time step durations for a random set of 85 initial configurations of 15 vehicles. For each instance, the resolution-optimal average time t τ opt is computed for time 1 Also available at https://youtu.be/RiW2OFsdHOY step durations τ ranging from 0.125 s to 5 s. The relative loss of optimality is computed as t τ opt -t 0.125 opt t 0.125 opt . Interestingly, the averaged values fit closely to an affine function with slope 7.2% per second for the above set of parameters. The loss of optimality remains less than 6% when the time step is smaller than 1 s; moreover, the solution of (3) converges as the time step duration vanishes.

The kinodynamic parameters, i.e. the maximum speed v and the acceleration bounds a i and a i influence the magnitude of the optimality loss. Fig. 7 shows a comparison of the losses of optimality for three different scenarios, namely "reference", "lower speed" and "higher acceleration". Reference parameters are, as before, (v, a i , a i ) = (15 m s -1 , -3 m s -2 , 4 m s -2 ). The same set of 85 initial configurations is used for those three scenarios, and we only vary the kinodynamic parameters; in the lower speed scenario, initial speeds are also reduced by 33%. We find that an increase of time step duration causes higher losses of optimality in instances with more dynamic robots (higher speeds or higher absolute values of acceleration bounds) than those with less dynamic ones. As a result, the time step duration can be adapted to the dynamic characteristics of robots; for instance using longer time steps for slower robots.

C. Computation time

To measure the computation time required to solve problem (3) for the intersection displayed in Fig. 1, we run the simulator for different numbers of vehicles with a fixed time step of 1 s (Fig. 8). Computation time remains below 1 s for up to 8 vehicles, which would make our approach suitable for real-world applications. Note that the same 30 s time horizon is considered across all simulations to provide fair comparison, while a shorter horizon could be used for problems involving fewer robots, thus reducing computation time; moreover, the presented MILP problem has been formulated for clarity rather than execution speed, notably by introducing redundant variables. As a result, performance can likely be further improved to process more robots in the same time frame.

VI. CONCLUSION

This article presents a new approach to compute timeoptimal trajectories for the coordination of mobile robots in a structured environment, taking into account kinodynamic constraints. Assuming each robot follows a predetermined path and using well-chosen bounding polygons, it is possible to formulate this problem as a mixed-integer linear program. By using temporal instead of spatial discretization, our approach allows an exact modeling of second-order integrator dynamics with piecewise-constant acceleration. The formulated problem mixes the discrete choice of relative priorities between conflicting robots with the continuous optimization of speed profiles respecting these priorities.

Extensive computer simulations have been run, on the example of a road intersection with autonomous vehicles, using the open-source traffic simulator SUMO and commercial MILP solver GUROBI. These simulations demonstrate the ability of our approach to compute optimal, collision-free trajectories in reasonable time.

Using consumer-grade hardware, the proposed method can treat up to 8 robots in less than a second. It is therefore suitable for real-time use, for instance for robots in an automated warehouse or autonomous vehicles at an intersection. In situations where computation time is more important than optimality, our formulation could also be used to design more efficient heuristics, with provable performance bounds, or to assess potentially dangerous situations where no safe trajectories exist. More general quadratic programming techniques could allow to use more complex objective functions, and relax the constant acceleration assumption while still providing an exact modeling of the dynamic constraints, opening many perspectives for future research.

Fig. 1 .

 1 Fig. 1. Example of paths inside and outside the coordination region for a two-lane road intersection.

Fig. 3 .

 3 Fig. 3. Exact shape of the collision region C ij for (polygonal) car-like robots. The blue rectangles represent the corresponding bounding polygon.

Fig. 4 .

 4 Fig. 4. Illustration of the completed collision set C i j in the merging case, and the corresponding values of S ⊥ ij , S

Fig. 5 .Fig. 6 .

 56 Fig. 5. Optimal trajectories within given priority classes, corresponding to a global (left) and local (right) optimum

Fig. 7 .

 7 Fig.7. Average relative optimality loss (compared to a 0.25 s time step), depending on time step duration, for 85 instances of 15 vehicles with reference parameters (red), lower maximum speed of 10 m s -1 (blue) and higher absolute values of accelerations [a i , a i ] = [-6, +8] m s -2 (green).

Fig. 8 .

 8 Fig. 8. Averaged computation time over 10 instances for a time step duration of 1 s and a varying number of vehicles.

MINES ParisTech, PSL Research University, Centre for robotics,

This work has received support from the European FP7 project Au-toNet2030 (Grant Agreement NO. 610542).


J J Liu
email: jjliu@seu.edu.cn

G Nakamura

M Sini

RECONSTRUCTION OF THE SHAPE AND SURFACE IMPEDANCE FROM ACOUSTIC SCATTERING DATA FOR ARBITRARY CYLINDER

Keywords: Inverse scattering, far-field, impedance boundary, singularity analysis, numerics AMS subject classifications. 35P25, 35R30, 78A45

The inverse scattering for an obstacle D ⊂ R 2 with mixed boundary condition can be considered as a prototype for radar detection of complex obstacles with coated and non-coated parts of the boundary. We construct some indicator functions for this inverse problem using the far-field pattern directly, without the necessity to transform the far-field to the near field. Based on the careful singularity analysis, these indicator functions enable us to reconstruct the shape of the obstacle and distinguish the coated from the non-coated part of the boundary. Moreover, an explicit representation formula for the surface impedance in the coated part of the boundary is also given. Our reconstruction scheme reveals that the coated part of the obstacle is less visible than the non-coated one, which corresponds to the physical fact that the coated boundary absorbs some part of the scattered wave. Numerics are presented for the reconstruction formulas, which show that both the boundary shape and the surface impedance in the coated part of boundary can be reconstructed accurately. The theoretical reconstruction techniques proposed in this work can be applied in the practical 3-dimensional electromagnetic inverse scattering problems with hopeful numerical performances, which are of great importance in the design of non-detectable obstacles.

Introduction and examples.

The inverse scattering problems aim to identify some properties of an obstacle such as the boundary shape and type from the information contained in the scattered wave for given incident waves. The optimization techniques are well-known for reconstructing the obstacle, up to some accuracy, by minimizing the objective functional for unknown obstacle from given inversion input data by iterations procedures. However, it seems that a good initial guess is needed.

In the recent years, some new inversion methods for the reconstruction of obstacle boundary have been proposed. The general and common idea of these methods is the construction of some indicator functions from given inversion input data, which depend on some detecting point (a parameter) varying inside or outside the obstacle. When this point approaches the obstacle, these indicator functions blowup. The linear sampling method [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF], the factorization method [START_REF] Kirsch | Characterization of the shape of a scattering obstacle using the spectral data of the far field operator[END_REF] and the singular sources method [START_REF] Potthast | Point Sources and Multipoles in Inverse Scattering Theory[END_REF] construct the indicator functions from the far-field pattern directly, while the probe method [START_REF] Ikehata | Reconstruction of the shape of the inclusion by boundary measurements[END_REF][START_REF] Ikehata | Reconstruction of obstacles from boundary measurements[END_REF] constructs the indicator in terms of the near field. The near field can be obtained from the far field by some regularization procedures [START_REF] Potthast | Point Sources and Multipoles in Inverse Scattering Theory[END_REF]. However, we can also state the natural version of the probe method directly from the far field data, without reducing the far field to the near field, see [START_REF] Honda | The no-response approach and its relation to non-iterative methods for the inverse scattering[END_REF]. For a review of these methods, the readers are referred to [START_REF] Potthast | Sampling and probe methods -an algorithmical view[END_REF][START_REF] Potthast | A survey on sampling and probe methods for inverse problems[END_REF] and for some relations between them to [START_REF] Honda | The no-response approach and its relation to non-iterative methods for the inverse scattering[END_REF][START_REF] Nakamura | Unification of the probe and singular sources methods for the inverse boundary value problem by the non-responce test[END_REF].

If the scattering is caused by multiple obstacles with different types of boundary or with mixed boundary condition, one should identify both the boundary shape, boundary type and the surface impedance. These kinds of problems come from some industry designs such as radar detection by electromagnetic wave scattering, see [START_REF] Hoppe | Impedance boundary condition in electromagnetism[END_REF]. The obstacle is illuminated by an electromagnetic wave coming from an antenna. The wave is scattered by the obstacle and received by an antenna located in a different place. One of the objectives is to design the shape of the obstacle such that the reflected wave can be avoided or minimized in some directions. One possible way for this goal is to introduce a coating on the surface of the obstacle or on some of its parts. This is motivated by the fact that reflections are minimized by applying such a surface coating. The surface coating is modelled by introducing an impedance boundary condition on a part or on the whole surface of the scatterer, which gives a relation between the electric and the magnetic field through a coefficient called surface impedance.

Due to this practical importance, the reconstruction of boundary impedance has been addressed by many authors. In [START_REF] Akduman | Direct and inverse scattering problems for inhomogeneous impedance cylinders of arbitrary shape[END_REF], the authors construct the inhomogeneous boundary impedance for a cylinder obstacle with known shape using only one incident wave, assuming that the surface impedance is distributed along the whole boundary of the obstacle. In this case, the scattering of electromagnetic waves can be described by the 2-dimensional Helmholtz equation. We also refer to [START_REF] Kress | Inverse scattering for shape and impedance[END_REF], where an optimization method is applied. After reducing the far field to the near field, a moment method is suggested in [START_REF] Cheng | Recovery of the shape of an obstacle and the boundary impedance from the far-field pattern[END_REF] to reconstruct the surface impedance approximately in the case of completely coated obstacle, and the identification of different types of multiple obstacles is given in [START_REF] Cheng | Recovery of boundaries and types for multiple obstacles from the far-field pattern[END_REF] in the case where on each obstacle we have one type of boundary condition.

The problem whether a part of the surface of the obstacle is coated or not is of relevant importance. Answering to this question and reconstructing the surface impedance, in case of coating, from far field measurements is our main object. In this work, we restrict ourself to the acoustic wave scattering governed by 2-dimensional Helmholtz equation, noticing that the 3-dimensional electromagnetic wave scattering in the cylinder case can be modelled by the 2-dimensional Helmholtz equation [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. These issues have been firstly considered in [START_REF] Cakoni | The determination of the surface impedance of a partially coated obstacle from far field data[END_REF][START_REF] Cakoni | The determination of the surface conductivity of a partially coated dielectric[END_REF], by the linear sampling method, where the authors simultaneously reconstruct the obstacle and compute the L ∞ -norm of the surface impedance. This can be used to answer the question of existence or absence of coating and give the value of the surface impedance in case it is known to be constant.

Motivated by these last works, our aim is to give another way to consider these issues and give further information on the obstacle. We proceed by constructing some indicator functions giving a direct link between the far-field pattern and the physical parameters of the obstacle. More precisely, we establish pointwise formulas which enable us to detect the boundary of the scatterer and distinguish and recognize the coated and the non-coated parts of the obstacle surface. In addition, on the coated part of the obstacle, the indicator functions give explicitly the pointwise values of this surface impedance as a functional of the far fields. These types of results have been initiated in [START_REF] Nakamura | Pointwise determination of the surface impedance from scattering data[END_REF], where the theoretical justification of these formulas in 3-dimensional acoustic scattering is given. Since we need more singularity analysis in the 2-dimensional case than in 3-dimensional case, which is due to the use of a more singular point source, we give the theoretical justification of the steps where it is necessary and refer to [START_REF] Nakamura | Pointwise determination of the surface impedance from scattering data[END_REF] for the rest of proof. We would like to emphasize that we are reconstructing the obstacle, localizing the eventual coated part and reconstruct the surface impedance in on step, i.e. simultaneously, compare with [[1], [START_REF] Cakoni | The determination of the surface impedance of a partially coated obstacle from far field data[END_REF], [START_REF] Cakoni | The determination of the surface conductivity of a partially coated dielectric[END_REF], [START_REF] Cheng | Recovery of boundaries and types for multiple obstacles from the far-field pattern[END_REF], [START_REF] Cheng | Recovery of the shape of an obstacle and the boundary impedance from the far-field pattern[END_REF], [START_REF] Kress | Inverse scattering for shape and impedance[END_REF]]. Also, since the analysis is done pointwise, then we can also consider the multiple obstacles and give similar results.

The validity of the theoretical reconstruction formula presented in this paper is also checked by numerical tests with satisfactory performances. We would like to mention the following observations from the numerics. The coated part of an obstacle with larger impedance is less visible than the other part in terms of the value of the indicator. This explains the practical motivation for introducing the coating, i.e, to avoid or perturb the detection of an obstacle by applying the absorbing boundary layer. On the other hand, for non-convex obstacles with mixed boundary conditions, the inversion formulas proposed in this paper also generate a satisfactory reconstruction by combining different blowing-up criterion together. These reconstruction performances are supported by our numerical implementations given in the last section of this paper.

The rest of the paper is organized as follows. In section 1, we state the problem mathematically. In section 2, we present the results which we prove in section 3. The section 4 is devoted to the numerical tests.

Statement of the problem.

Let D be a bounded domain of R 2 such that R 2 \ D is connected. We assume that its boundary ∂D is of class C 2 and has the following form

∂D = ∂D I ∪ ∂D D , ∂D I ∩ ∂D D = ∅,
where ∂D D and ∂D I are open surfaces in ∂D.

The propagation of time-harmonic acoustic fields in homogeneous cylinder media can be modeled by the Helmholtz equation

∆u + κ 2 u = 0 in R 2 \ D, (2.1)
where κ > 0 is the wave number. At the part ∂D I of the obstacle boundary, we assume the total field u to satisfy the impedance boundary condition while on the part ∂D D to satisfy the Dirichlet boundary condition. That is,

∂u ∂ν + iκσu = 0 on ∂D I (2.2)
with some impedance function σ and

u = 0 on ∂D D , (2.3)
where ν is the outward unit normal of ∂D. We assume that σ is a real valued Holder continuous function of order β ∈ (0, 1] and has a uniform lower bound σ -> 0 on ∂D I . The part ∂D I is referred to by the coated part of ∂D and ∂D D is the non-coated part. For a given incident plane wave u i (x, d) = e iκd•x with incident direction d ∈ S 1 , where S 1 is a unit circle in R 2 , we look for a solution u := u i + u s of (2.1), (2.2) and (2.3), where the scattered field u s satisfies the Sommerfeld radiation condition

lim r→∞ √ r( ∂u s ∂r -iκu s ) = 0 (2.4)
with r = |x| and the limit is uniform for all directions x ∈ S

       (∆ + κ 2 )u = 0, in R 2 \ D, u = f on ∂D D , ∂u ∂ν + iκσu = h, on ∂D I , lim r→∞ √ r( ∂u ∂r -iκu) = 0, (2.5)
and the solution satisfies

u H 1 (ΩR∩(R 2 \D)) ≤ C R ( f H 1/2 (∂DD) + h H -1 2 (∂D I ) ) (2.6)
where Ω R is a disk of radius R and C R is positive constant depending on R, see [START_REF] Cakoni | Qualitative Methods in Inverse scattering Theory, Interaction of Mechanics and Mathematics[END_REF] for more details.

It is well known, see [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF], that the scattered wave has the asymptotic behavior:

u s (x, d) = e iκr √ r u ∞ (x, d) + O(r -3/2 ), r := |x| → ∞, (2.7)
where the function u ∞ (•, d) defined on S 1 is called the far-field of the scattered wave u s corresponding to incident direction d. We introduce a constant γ 2 := e iπ/4 √ 8πκ and Φ(x, y)

:= i 4 H (1) 0 (κ|x -y|), x = y, x, y ∈ R 2 ,
the fundamental solution to the Helmholtz equation in R 2 , where H

(1) 0

is the Hankel function of the first kind of order zero. In this paper, we will consider the following Inverse scattering problem for an obstacle with mixed boundary type: Given u ∞ (•, •) on S 1 × S 1 for the scattering problem (2.1) -(2.4), reconstruct the shape of obstacle D, identify ∂D I and ∂D D and reconstruct the surface impedance σ(x) on ∂D I .

Presentation of the results.

It is well known also, see [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF], that the scattered field associated with the Herglotz incident field v i g := v g defined by

v g (x) := S 1 e iκx•d g(d) ds(d), x ∈ R 2 (3.1)
with g ∈ L 2 (S 1 ) is given by

v s g (x) := S 1 u s (x, d)g(d) ds(d), x ∈ R 2 \ D, (3.2)
and its far field is

v ∞ g (x) := S 1 u ∞ (x, d)g(d) ds(d), x ∈ S 1 . (3.3)
We will need the following identity, see [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF],

u ∞ (x, d) = -γ 2 ∂D ∂u s (y, d) ∂ν e -iκx•y - ∂e -iκx•y ∂ν u s (y, d) ds(y) (3.4)
Φ s (x, z) = - ∂D ∂Φ s (y, z) ∂ν(y) Φ(x, y) -Φ s (y, z) ∂Φ(x, y) ∂ν(y) ds(y), x, z ∈ R 2 \ D.
(3.5) Assume that D ⊂⊂ Ω for some known Ω with smooth boundary. For a ∈ Ω \ D, denote by {z p } ⊂ Ω \ D a sequence tending to a. For any z p , set D p a a C 2regular domain such that D ⊂ D p a with z q ∈ Ω \ D p a for every q = 1, 2, • • • , p and that the Dirichlet interior problem on D p a for the Helmholtz equation is uniquely solvable, see Figure 3.1 for the configuration. In this case, the Herglotz wave operator H defined from L 2 (S 1 ) to L 2 (∂D p a ) by

H[g](x) := v g (x) = S 1 e iκx•d g(d) ds(d) (3.6)
is injective, compact with dense range, see [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. Let z * p be a point on ∂D p a near z p such that z * p → a as z p → a, as chosen in Figure 3.1. Denote by ν(z * p ) the outward normal of ∂D p a at z * p . Now we consider the sequence of point sources Φ(•, z p ). For every p fixed, we construct two density sequences {g p n } and {f p m } in L 2 (S 1 ) by the Tikhonov regularization such that

v g p n -Φ(•, z p ) L 2 (∂D p a ) → 0, n → ∞ (3.7) v f p m - ∂ ∂ν(z * p ) Φ(•, z p ) L 2 (∂D p a ) → 0, m → ∞, (3.8)
where

∂ ν(z * p ) Φ(•, z p ) := ∇ x Φ(x, z p )•ν(z * p ). Since both v g p n and Φ(•, z p ) satisfy the same Helmholtz equation in D p a , (3.7) implies that v g p n -Φ(•, z p ) H 1 2 (∂D) → 0, n → ∞ (3.9)
and

∂ ∂ν v g p n - ∂ ∂ν Φ(•, z p ) H -1 2 (∂D) → 0, n → ∞ (3.10)
Similarly, it follows from (3.8) that

v f p m - ∂ ∂ν(z * p ) Φ(•, z p ) H 1 2 (∂D) → 0, m → ∞ (3.11)
and

∂ ∂ν v f p m - ∂ ∂ν ( ∂ ∂ν(z * p ) Φ(•, z p )) H -1 2 (∂D) → 0, m → ∞ (3.12)
Multiplying (3.4) by f p m (d)g p n (x) and integrating over S 1 × S 1 , we have

- S 1 S 1 u ∞ (-x, d)f p m (d)g p n (x) ds(x)ds(d) = γ 2 ∂D S 1 ∂u s (y, d) ∂ν f p m (d) ds(d) • S 1 e iκx•y g p n (x) ds(x) - S 1 ∂e iκx•y ∂ν g p n (x) ds(x) • S 1 u s (y, d)f p m (d) ds(d) ds(y) = γ 2 ∂D ∂v s f p m ∂ν (y)v i g p n (y) - ∂v i g p n ∂ν (y)v s f p m (y) ds(y). (3.13)
From (3.9), (3.10) and (3.13), we have

lim n→∞ S 1 S 1 u ∞ (-x, d) f p m (d) g p n (x) ds(x)ds(d) = γ 2 ∂D v s f p m ∂Φ(y, z p ) ∂ν(y) - ∂v s f p m ∂ν(y) Φ(y, z p ) ds(y) = γ 2 v s f p m (z p ) (3.14)
from the Green formula, where v s

f p m (•) is the scattered wave corresponding to incident wave v i f p m (x) = H[f p m ](x)
. Denote by E s (x, z p ) the scattered wave corresponding to the incident wave ∂Φ(x,zp) ∂ν(z * p ) , which is well defined for every x ∈ R 2 \ D. Then it follows from (3.11), (3.12), the well posedness of the direct scattering problem and the use of interior estimate that

E s (x, z p ) = lim m→∞ v s f p m (x), x ∈ R 2 \ D. (3.15)
Finally, it follows from (3.14) that lim

m→∞ lim n→∞ S 1 S 1 u ∞ (-x, d) f p m (d) g p n (x) ds(x)ds(d) = γ 2 E s (z p , z p ). (3.16)
The reconstruction of ∂D as well as its surface impedance in the coating part can be established based on (3.16). For this purpose, an analysis of E s (x, z) near ∂D is the key point. We need the natural C 2 smoothness assumption on the regularity of ∂D. Precisely, for every point a ∈ ∂D, there exists a rigid transformation of coordinates under which the image of a is 0 and a function f ∈ C 2 (-r, r) such that

f (0) = df dx (0) = 0, D ∩ B(0, r) = {(x, y) ∈ B(0, r); y > f (x)} (3.17)
in terms of the new coordinates where B(0, r) is the 2-dimensional ball of center 0 with radius r.

For the points a ∈ ∂D, we choose the sequence {z p } p∈N included in C a,θ , where C a,θ is a cone with center a, angle θ ∈ [0, π 2 ) and axis ν(a). The main theoretical result of this paper is as follows.

Theorem 3.1. Assume that the boundary ∂D is of class C 2 and σ is a real valued Holder continuous function with positive lower bound. Then the boundary properties of the obstacle D can be identified by the following indicator functions.

1. The obstacle boundary ∂D can be constructed from the following property

lim p→∞ lim m,n→∞ Re γ -1 2 S 1 S 1 u ∞ (-x, d)f p m (d)g p n (x) ds(x)ds(d) = +∞, a ∈ ∂D, < +∞, a ∈ Ω \ D. (3.18)
Precisely, we have the blowup rate:

lim m,n→∞ Re γ -1 2 S 1 S 1 u ∞ (-x, d)f p m (d)g p n (x) ds(x)ds(d) = ±1 4π|(z p -a) • ν(a)| + (3.19) +O(| ln |z p -a|| 2 ),
where z p := (z p,1 , z p,2 ) and a = (a 1 , a 2 ). The sign (+) is for a ∈ ∂D D while the sign (-) is for a ∈ ∂D I .

2. The coating and the non-coating parts of ∂D can also be distinguished from the following properties

lim p→∞ lim m,n→∞ Im γ -1 2 S 1 S 1 u ∞ (-x, d)f p m (d)g p n (x) ds(x)ds(d) | ln |(z p -a) • ν(a)|| s = +∞, a ∈ ∂D I , 0, a ∈ ∂D D (3.20
) by choosing any fixed s ∈ (0, 1).

3. The impedance coefficient on ∂D I can be detected by the following formula:

lim p→∞ lim m,n→∞ Im γ -1 2 S 1 S 1 u ∞ (-x, d)f p m (d)g p n (x) ds(x)ds(d) | ln |(z p -a) • ν(a)|| = κ π σ(a), a ∈ ∂D I . (3.21) Remark 3.2. The formula (3.18) is also true if f p m is replaced by g p n .
That is, the singularity of Φ(x, z p ) is theoretically enough for identifying ∂D. However, as the blowup rate in this 2-dimensional case is of logarithmic order, it is not suitable to localize the obstacle clearly in the numerical experiments. Due to this reason, we introduced the density f p m which is related to a stronger singularity ∂ ∂ν(z * p ) Φ(•, z p ) to get a blowup rate of order |z p -a| -1 . For the formulas (3.20) and (3.21), the stronger singularity of ∂ ∂ν(z * p ) Φ(•, z p ) is necessary. Moreover, we can in fact use ν(a) instead of ν(z * p ), since ∂D has been determined in terms of (3.18). The formula (3.19) can also be used to distinguish the coated part ∂D I from the non-coated ∂D I .

Remark 3.3. Theorem 3.1 is stated for the case of single obstacle. But these results are still true for the multiple obstacle case with coated and non-coated parts.

Remark 3.4. If a ∈ D, then the limit in (3.18) is conjectured to be ∞, see [START_REF] Honda | The no-response approach and its relation to non-iterative methods for the inverse scattering[END_REF]. However up to now, we do not have the full answer. The approach in [START_REF] Ikehata | A new formulation of the probe method and related problems[END_REF] can be used to justify it in the case where the frequency κ is small enough.

4. Proof of Theorem 3.1.. For any given point a ∈ ∂D, we firstly take the rotation R a and the translation M a such that R a (ν(a)) = (0, 1), R a (a) + M a = 0 in the new coordinate system x. Under the transform x := T(x) := R a (x) + M a , it follows that T(ν(a)) = (0, 1), T(a) = 0.

Define σ(x) := σ(x) and consider the following two problems in the coordinate x = (x 1 , x2 ) for any given z = (z 1 , z2 ) ∈ R 2 + . We set w+ σ(0) (x, z) and w+ D (x, z) to be two functions satisfying

∆ w+ σ(0) = 0, x ∈ R 2 + ( ∂ ∂ x2 w+ σ(0) + iκσ(0) w+ σ(0) )(x, z)| x2=0 = -( ∂ ∂ x2 + iκσ(0)) ∂ ∂x2 Γ(x, z)| x2=0 , (4.1)
∆ w+

D = 0, x ∈ R 2 + , w+ D (x, z)| x2=0 = -∂ ∂x2 Γ(x, z)| x2=0 (4.2)
respectively, where Γ(x, z) = 1 2π ln 1 |x-z| and the subscript D in w+ D (x, z) refers to the Dirichlet boundary condition in (4.2).

We give explicit solutions to these two problems in the following proposition. Proposition 4.1. We have the explicit form of w + σ(0) (x, z)

w+ σ(0) (x, z) = 1 4π R e i(x1-z1)ξ1 e -(x2+z2)|ξ1| |ξ 1 | + iκσ(0) |ξ 1 | -iκσ(0) dξ 1 , (4.3)
while w+ D (x, z) has the form

w+ D (x, z) = - 1 4π R e i(x1-z1)ξ1 e -(x2+z2)|ξ1| dξ 1 . (4.4)
This proposition can be proven by expressing z)dξ 1 and computing the density functions φ ± from the boundary value problems (4.1), (4.2), where φ is the 1-dimensional Fourier transform of φ, see [START_REF] Nakamura | Pointwise determination of the surface impedance from scattering data[END_REF] for explicit computations. Define

w+ σ(0) (x, z) = (U + [x 2 ]φ + )(x 1 ), w+ D (x, z) = (U + [x 2 ]φ -)(x 1 ) in R 2 + with (U ± [x 2 ]φ)(x 1 ) := 1 2π R e ix1ξ1∓x2|ξ1| φ(ξ 1 ,
w + σ(a) (x, z) = w+ σ(0) (Tx, Tz), w + D (x, z) = w+ D (Tx, Tz) for x, z ∈ R 2 \ D near a,
|Im E s (x, z) -Im w + σ(a) (x, z)| ≤ C, for (x, z) ∈ B + (a, δ(a)) ∩ C a,θ , (4.5
) Proof of Theorem 3.1.

|Re E s (x, z)-Re w + σ(a) (x, z)| ≤ C| ln |x-a||•| ln |z p -a||, for (x, z) ∈ B + (a, δ(a))∩C a,θ , (4.6
Step A: It follows from Proposition 4.2 that

|Re E s (x, z p ) -Re w + (x, z p )| ≤ C ln 1 |z p -a| (4.7)
uniformly for all x, z p near any fixed point a ∈ ∂D, where w + (z p , z p ) may be w + σ(a) (z p , z p ) or w + D (z p , z p ), depending on the position of a. For w

+ (z p , z p ) = w + D (z p , z p ), it follows from (4.4) that Re w + (z p , z p ) = - 1 4π R e -2|zp,2-a2||ξ1| dξ 1 = 1 4π|z p,2 -a 2 | , while for w + (z p , z p ) = w + σ(a) (z p , z p ) it holds from (4.3) that Re w + (z p , z p ) = 1 4π R e -2|zp,2-a2||ξ1| |ξ 1 | 2 -κ 2 σ 2 (a) |ξ 1 | 2 + κ 2 σ 2 (a) dξ 1 = - 1 4π|z p,2 -a 2 | + O(1)
, Suppose that a is outside D. We can construct z * p , z p tending to a as we did for a ∈ ∂D. Recall that E s (x, z p ) satisfies

where z p = (z p,1 , z p,2 ) → a = (a 1 , a 2 ) ∈ ∂D as p → ∞.
         (∆ + κ 2 )E s (x, z p ) = 0 in R 2 \ D, E s (•, z p ) = -∂Φ ∂ν(z * p ) (x, z p ) on ∂D D ( ∂ ∂ν + iκσ(x))E s (x, z p ) = -(∂ ν + iκσ(x)) ∂Φ ∂ν(z * p ) (x, z p ) on ∂D I E s (•, z p ) satisfies the Sommerfeld radiation conditions,
where ν(z * p ) is the unit outward normal on ∂D p a at the point z * p . Hence the boundary condition is bounded with respect to x in H 1/2 (∂D D ) and H -1 2 (∂D I ), respectively, for z * p , z p near a. It follows from the well-posedness of the direct problem and interior estimates near a (i.e. away from ∂D) that E s (x, z p ), and then E s (z p , z p ), is bounded.

Step B. Let a ∈ ∂D I . From (4.3) we have

w+ σ(0) (z, z) = 1 4π R e -2z2|ξ1| |ξ 1 | + iκσ(0) |ξ 1 | -iκσ(0) dξ 1 . (4.8)
By taking the imaginary part and setting z = (z

1 , z2 ) = R a (z) + M a for z ∈ C(a, θ),
we get The rest of this section is devoted to the proof of Proposition 4.2. As we said in the introduction, in the 2D case, we need more singularity analysis than in [START_REF] Nakamura | Pointwise determination of the surface impedance from scattering data[END_REF]. This is due to the use of the more singular point source ∂ ∂ν(z * p ) Φ(•, z p ). We give the detailed analysis and refer to [START_REF] Nakamura | Pointwise determination of the surface impedance from scattering data[END_REF] for the steps which do not need important changes.

Im(4πw + σ(a) (z, z)) = 4κσ(a) +∞ 0 e -2(z-a)•ν(a)r r r 2 + κ 2 σ(a) 2 dr = 4κσ(a) [-ln(κσ(a)) -ln((z -a) • ν(a))+ 2 +∞ 0 ln(r 2 + κ 2 |(z -a) • ν(a)| 2 σ 2
4.1. Proof of Proposition 4.2.. We give the proof for a ∈ ∂D I . The proof for a ∈ ∂D D is similar.

Let Ẽs (x, z p ) be the solution of

     (∆ + κ 2 ) Ẽs (x, z p ) = 0 in R 2 \ D, ( ∂ ∂ν + iκσ(x)) Ẽs (x, z p ) = -(∂ ν + iσ(x)) ∂ ∂ν(z * p ) Φ(x, z p ) on ∂D Ẽs (•, z) satisfies the Sommerfeld radiation condition. (4.10) Hence (E s -Ẽs )(x, z p ) satisfies          (∆ + κ 2 )(E s -Ẽs )(x, z p ) = 0 in R 2 \ D, ( ∂ ∂ν + iκσ(x))(E s -Ẽs )(x, z p ) = 0 on ∂D I (E s -Ẽs )(•, z p ) = -∂ ∂ν(z * p ) Φ(x, z p ) -Ẽs on ∂D D (E s -Ẽs )(•, z) satisfies the Sommerfeld radiation condition. (4.11) We state H σ (x, z) := Ẽ(x, z) + ∂ ν(z * p ) Φ(x, z). Hence H satisfies    (∆ + κ 2 )H σ (x, z) = -∇δ(x, z) • ν(z * p ) in R 2 \ D, ( ∂ ∂ν + iκσ(x))H σ (x, z) = 0 on ∂D H σ (•, z) satisfies the Sommerfeld radiation condition.
(4.12)

We have the following estimates:

       |G σ (x, z)| ≤ c| ln |x -z|| |∇G σ (x, z)| ≤ c|x -z| -1 |H σ (x, z)| ≤ c|x -z| -1 |∇H σ (x, z)| ≤ c|x -z| -2 in R 2 \ D where c is a positive constant. (4.13)
The justification of these properties can be derived following, for instance, the approach of [START_REF] Solonnikov | On Green's matrices for elliptic boundary value problems(I)[END_REF] and [START_REF] Solonnikov | The Green's matrices for elliptic boundary value problems(II), Boundary Value Problems of Mathematical Physics[END_REF] since an explicit form of a local fundamental solution for the half space case can be derived as we did in Proposition4.1. See also [START_REF] Littman | Regular points for elliptic equations with discontinuous coefficients[END_REF] and [START_REF] Gruter | The Green function for uniformly elliptic equations[END_REF] for the case of elliptic problems with rough coefficients. From these estimates, we deduce that Ẽ(•, z p ) and its derivatives are bounded for x ∈ ∂D D and z p near a ∈ ∂D I . The well posed-ness of (4.11) implies that (E s -Ẽs )(•, z p ) is bounded in H 1 loc (R 2 \D) for z p near a. Introducing a cutoff function near the point a and using (4.11), we deduce that (E s -Ẽs )(•, z p ) is bounded for x near ∂D I and z p near a.

This means that we can replace in Proposition 4.2 E s by Ẽs . We introduce w s σ (•, z p ) as the solution of

   (∆ + κ 2 )w s σ (x, z p ) = 0 in R 2 \ D, ( ∂ ∂ν + iκσ(x))w s σ (x, z p ) = -(∂ ν + iσ(x)) ∂ ∂ν(a) Φ(•, z p ) on ∂D w s σ (•,
|( Ẽs -w s σ )(x, z p )| ≤ C(R)| ln d(x, ∂D)| • | ln d(z p , ∂D)|, |Im( Ẽs -w s σ )(x, z p )| ≤ C(R), |(w s σ -w s σ(a) )(x, z p )| ≤ C(R) for z p ∈ B(a, δ(a)) ∩ C a,θ and x ∈ (R 2 \ D) ∩ B(0, R), for any R > 0 fixed. Let w s σ(a),Φ (•, z) be the solution of      (∆ + κ 2 )w s σ(a),Φ (x, z) = 0 in Ω \ D, ( ∂ ∂ν + iκσ(a))w s σ(a),Φ (x, z) = -( ∂ ∂ν + iκσ(a)) ∂ ∂ν(a) Φ(x, z p ) on ∂D w s σ(a),Φ (•, z) = -∂ ∂ν(a) Φ(x, z p ) on ∂Ω
+ σ(a) )(x, z)| ≤ C for (x, z) ∈ B(a, δ(a)) ∩ C a,θ . 2. |(Rew s,0 σ(a) -Rew + σ(a) )(x, z)| ≤ C| ln d(z p , ∂D)| for (x, z) ∈ B(a, δ(a)) ∩ C a,θ
. By combining all the lemmas stated above, we end the proof of Proposition 4.2.

In the proofs of these lemmas we do not, in general, specify the interdependency of the constants appearing in the estimates. However we distinguish the constant depending on the angle θ and the ones which do not depend.

Proof of Lemma 4.4. The function Ẽs (x,

z p ) -w s σ (x, z p ) satisfies    (∆ + κ 2 )( Ẽs -w s σ )(x, z p ) = 0 in R 2 \ D, ( ∂ ∂ν + iκσ)( Ẽs -w s σ )(x, z p ) = -( ∂ ∂ν + iκσ)∇Φ • [ν(z * p ) -ν(a)] on ∂D, ( Ẽs -w s σ )(•, z p )
(i) a ∈ ∂D l(0) a , D ⊂ D l(0) a ⊂ Ω. (ii) ∂D l(1) a = ∂Ω, D l(s) a ⊂ D l(s ′ ) a (0 ≤ s < s ′ ≤ 1). (iii) l intersects ∂D l(s) a at l(s). (iv) ∂D l(s) a depends C 2 smoothly on s ∈ [0, 1].
For every s ∈ [0, 1] we denote ν(l(s)) to be the unit normal of ∂D l(s) a at l(s). From (iv), the map s ∈ [0, 1] → ν(l(s)) is C 1 . Choosing l to be a one to one curve near a, we deduce that the map l(s) ∈ l([0, 1]) → ν(l(s)) is also C 1 near l(0). Now let {s p } ⊂ [0, 1] be such that s p → 0 (p → ∞) and l be a C 3 smooth curve such that z * p := l(s p ) and {z p } p∈N ⊂ l([0, 1]). Since for every p, z p and z * p are in Ω \ D, then we can always choose l such that l(s) ∈ Ω \ D, (s ∈ [0, 1]). We set D p a := D l(sp) a

. Hence the sequence ν(z * p ) satisfies Claim 4.8.

The function Im( Ẽs

-w s σ ) satisfies    (∆ + κ 2 )Im( Ẽs -w s σ )(x, z) = 0 in R 2 \ D, ∂ ∂ν Im( Ẽ-w s σ )(x, z) = [κσRe( Ẽs -w s σ ) -[ ∂ ∂ν Im∇Φ + κσRe∇Φ] • [ν(z * p ) -ν(a)]](x, z) on ∂D.
(4.17

-( Ẽs -w s σ )(x, z p ) = ∂D G σ (x, y)( ∂ ∂ν + iκσ(y))[∇Φ • (ν(z * p ) -ν(a))](y, z p )ds(y).
Hence

|( Ẽs -w s σ )(x, z p )| ≤ c ∂D (ln |x -y|)|y -z p | -
C α > 0 such that |( Ẽs -w s σ )(x, z p )| ≤ C α |x -z p | -α .
From the explicit form of Φ, we have Similarly to the proof of Lemma 4.4, the proofs are based on the use of integral representations and the pointwise estimates of the Green's functions of the corresponding problems. So, we omit the details.

| ∂ ∂ν Im∇Φ(x, z p )|•|ν(z * p )-ν(a)| ≤ c|x-z p | -1 |z * p - a| and Re∇Φ(x, z p )|•|ν(z * p )-ν(a)| ≤ c|x-z p | -
:= w s σ (x, z)-w s σ(a) (x, z). Then it satisfies    (∆ + κ 2 )R(x, z) = 0 in R 2 \ D, ∂R(x,z) ∂ν + iκσ(a)R(x, z) = -iκ(σ(x) -σ(a))(w s λ (x, z) + ∂ ∂ν(a) Φ(x, z)) on ∂D, R(•,
Proof of Lemma 4.7. Since w s,0 σ(a) satisfies

     ∆w s,0 σ(a) (x, z) = 0 in Ω \ D, ( ∂ ∂ν + iκσ(a))(w s,0 σ(a) (•, z)) = -( ∂ ∂ν + iκσ(a)) ∂ ∂ν(a) Γ on ∂D, w s,0 σ(a) (•, z) = -∂ ∂ν(a) (Γ) on ∂Ω, (4.22)
then it is clear that G 0 σ(a) := w s,0 σ(a) (x, y)

+ ∂ ν(a) Γ(x, y) satisfies      ∆(G 0 σ(a) )(x, z) = -∂ ∂ν(a) δ(x -y) in Ω \ D, ( ∂ ∂ν + iκσ(a))(G 0 σ(a) )(•, z) = 0 on ∂D, (G 0 σ(a) )(•, z) = 0 on ∂Ω. (4.23)
We can assume without loss of generality that a = (0, 0) and ν(a) = (0, 1) by using the rigid transformation of coordinates [R a (ν(a)) + M a ]. Let ξ = F (x) be the local change of variables

ξ 1 = x 1 , ξ 2 = x 2 -f (x 1 ), (4.24)
where f is the function defined in the introduction. We have the following properties:

   c 1 |x -z| ≤ |F (x) -F (z)| ≤ c 2 |x -z|, |F (x) -x| ≤ c 3 |x| 2 , |DF (x) -I| ≤ c 4 |x| (4.25)
for x, z near the point a, where c i (i = 1, 2, 3, 4) are positive constants, which is due to hypothesis on the regularity of ∂D. Let x, z be points near a. From (4.23), we deduce that G0 σ(a) (ξ, η) = G 0 σ(a) (x, z) satisfies:

∇ ξ • B(ξ)∇ ξ G0 σ(a) = -J T (ξ)∇ ξ δ(ξ -η) • (0, 1) near F (a), |J -T ν|B(ξ)∇ ξ G0 σ(a) • ν + iκσ(a) G0 σ(a) = 0 on ∂R 2 + near F (a), (4.26)
where ξ := F (x) and η := F (z), B := JJ T and J := ∂ξ ∂x (F -1 (ξ)) and ν := (0, 1) is the unit normal to ∂R 2 + . We denoted by J -T the adjoint of J -1 . We have from (4.25) that

|J T (ξ) -J T (0)| ≤ c|ξ|, |B(ξ) -B(0)| ≤ c|ξ| and J(0) = B(0) = I. We set Γ σ(a) (x, z) := (w + σ(a) + ∂ ∂x2 Γ)(x, z) and write R(ξ, η) := G0 σ(a) (ξ, η) - Γ σ(a) (ξ, η). Then the function R(•, η) satisfies ∇ ξ • B(ξ)∇ ξ R = ∇ ξ • (I -B)∇ ξ Γ σ(a) + (I -J T (ξ))∇ ξ δ(ξ -η) • (0, 1), B(ξ)∇ ξ R • ν + iκσ(a) R = (I -B)∇ ξ Γ σ(a) • ν + iκσ(a)(1 -|J -T ν| -1 ) R, (4.27)
where the first relation holds in R 2 + near F (a), while the second one is satisfied on ∂R 2 + near F (a). We remark that (I -J T (ξ)) is equal to the matrix given by the two line vectors (0, -f (ξ)) and (0, 0). Hence we have

(I -J T (ξ))∇ ξ ImΓ σ(a) (ξ, η) • (0, 1) = 0.
With this remark, the problem (4.27) is exactly the one studied in [START_REF] Nakamura | Pointwise determination of the surface impedance from scattering data[END_REF]. We write ∂B + r = S r ∪ S c r with S r := ∂B + r ∩ ∂F (D). Arguing as in [START_REF] Nakamura | Pointwise determination of the surface impedance from scattering data[END_REF], we obtain the following estimate We go back to R(x, z)

:= G 0 σ(a) (x, z) -Γ σ(a) (x, z). We have R(x, z) = G 0 σ(a) (x, z) -Γ σ(a) (F (x), F (z)) + Γ σ(a) (F (x), F (z)) -Γ σ(a) (x, z),
which can be rewritten as

R(x, z) = R(F (x), F (z)) + [Γ σ(a) (F (x), F (z)) -Γ σ(a) (F (x), z)]+ + [Γ σ(a) (F (x), z) -Γ σ(a) (x, z)]. (4.28)
By the same argument as in [START_REF] Nakamura | Pointwise determination of the surface impedance from scattering data[END_REF], we end up with the estimate:

|Im R(x, z)| ≤ c(θ) (4.29) for x ∈ B(a, δ(a)) such that F (x) ∈ S r and z ∈ C a,θ ∩ B(a, δ(a)). For z ∈ C a,θ ∩ B(a, δ(a) 2 ) and x ∈ [∂B(a, δ(a))] ∩ R 2 \ D, we have |Im R(x, z)| ≤ c (4.30)
with some positive constant c, because C a,θ ∩ B(a, δ(a)

2 ) and [∂B(a, δ(a))] ∩ R 2 \ D are separated sets. Since in B(a, δ(a)) ∩ (R 2 \ D), we have ∆ x Im R(x, z) = 0, then using (4.29) and (4.30

) we have |Im R(x, z)| ≤ c(θ) for x ∈ [R 2 \ D] ∩ B(a, δ(a)) and z ∈ C a,θ ∩ B(a, δ(a)
2 ), by the maximum principle. Similarly we have |ReR(x, z)| ≤ C| ln |x -z|| for x ∈ B(a, δ(a)) such that F (x) ∈ S r and z ∈ C a,θ ∩ B(a, δ(a)). Hence |ReR(x, z)| ≤ C| ln d(z, ∂D)| for x ∈ B(a, δ(a)) such that F (x) ∈ S r and z ∈ C a,θ ∩ B(a, δ(a)) which is the counterpart of (4.30) for Re R. The rest of the proof is the same as the one for the imaginary part.

Numerical tests.

In this section, we consider two reconstruction model problems for numerical tests based on Theorem 3.1. In the first model we take the obstacle to be a disc and on its whole boundary we impose the impedance boundary condition. The purpose of considering such a model is to show the influence of wave number κ and singularity strength on the inversion scheme. In the second model, we consider a non-convex obstacle with mixed boundary conditions. We check our inversion formulas fully, especially for the identification of different type of boundary. Also we show the effect of the non-convex part on the inversion performance.

In the reconstruction scheme, the approximation of Φ(x, z p ) and Φ ν(a) (x, z p ) by the Herglotz wave function plays a key role. To do this, we need to construct the approximate domain D p a in the way mentioned in Section 2. Then their density functions can be determined by the standard argument of minimum norm solution of the integral equation of the first kind. In this way the numerators in (3.18), (3.20) and (3.21) can be computed for every sequence of points (z p ) p∈N approaching to ∂D in terms of the far-field pattern.

In testing our inversion scheme, we simulate the inversion input data (far-field pattern) by solving the direct problem using the combined single-and double-layer potential method, see [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF] and [START_REF] Colton | Integral Equation Methods in Scattering theory[END_REF].

In the subsections 4.1 and 4.2, we consider the first model, while subsection 4.3 is dedicated to a non-convex obstacle with mixed boundary conditions, from which we can test our theoretical results.

Reconstruction of ∂D. For the numerical test we take

∂D = {x = (x 1 (t), x 2 (t)) = 1.2(cos t, sin t) : t ∈ [0, 2π]}. We can construct D p
a in a special way by

∂D p a = {(x 1 (t), x2 (t)) = (1.2 + δ 0 (p))(cos t, sin t) : t ∈ [0, 2π]}
and take z p = (x 1 (t 0 ), x2 (t 0 ))+δ 1 (p)(cos t 0 , sin t 0 ) outside of D p a for a = 1.2(cos t 0 , sin t 0 ) in ∂D. The smallness of δ 1 (p) determines the singularity of Φ(x, z p ), Φ ν(a) (x, z p ) on ∂D p a near z p . Example 1. Construction of ∂D with unknown impedance σ(x) for given far-field data.

It follows from the first point of Theorem 3.1 that the boundary ∂D can be reconstructed from the blow up behavior of the approximate values of the indicator function:

I m,n (z p ) :≈ π N 2 2N -1 i,j=1 Re γ -1 2 U ∞ (i, j)F p m (j)G p m (i) (5.1)
for large m, n, where

G p n (i) := g p n (x i ), F p m (j) := f p m ( dj ) and U ∞ (i, j) represents u ∞ (-x i , d j ).
Therefore, we can test the numerical performance of this formula by taking z p approaching to ∂D. If I m,n (z p ) is greater than some a-priori given large value, for fixed m, n large enough, we consider z p to be almost in ∂D. Notice that only the singularity in Φ ν(zp) (x, z p ) is needed in this procedure, we do not need the boundary ∂D.

In our numerical implementations, we take κ = 0.6, N = 32 and keep the singularity of Φ(•, z p ) unchanged near ∂D p a by fixing δ 1 = 0.015. When z p approaches to ∂D by decreasing δ 0 from different directions t 0 , we get different values of I m,n (z p ). We choose the same constant CB as the criterion for the blowing-up of the indicator I m,n (z p (t 0 )) at different direction t 0 . That is, we choose z p (t 0 ) as the approximation of the point a ∈ ∂D in the direction t 0 when I m,n (z p (t 0 )) is larger than CB. Then the constructed approximate position of ∂D is given by interpolating those points. In this way, we can draw the approximate shape of ∂D by choosing all directions t 0 using the same given CB. As obtained in the theoretical result, the larger the value CB is, the better the approximation of ∂D should be.

Consider the boundary reconstruction problem with non-constant impedance

σ(x) = 2 + x 1 x 2 (3 + x 2 ) 2 , x ∈ ∂D.
Under this configuration, the reconstructions for different blowing-up criteria are shown in Figure 5.1 (left), while the distribution of indicator value is also given in the right picture.

-

-0.5 0 0.5 There is an interesting phenomena in the numerics. For a reasonable indicator value CB = 2, we can see the whole rough shape of ∂D. However, for larger values, i.e. CB = 3, most of the part of ∂D can be seen with more satisfactory accuracy, but some part, i.e. ( AB), is not visible. In this part the indicator value is less than 3 (but, of course, bigger than 2). This numerical performance is closely related to the impedance distribution in ∂D. It can be seen from the right part of Figure 5.1 that the indicator value is obviously smaller in the narrow domain at each radius layer. Therefore as r decreases, the part of ∂D related to these angles can not be detected with the same accuracy. Considering the distribution of boundary impedance, this part corresponds to σ(x) with large value, so it can not be seen by using the same criterion values CB, as clearly as the other part. This may be explained by the fact that the scattered wave along these directions will be much absorbed. Another, but related, reason is the property (3.19)

1 1.5 CB=0.8 CB=2.0 CB=3.0 Exact O A B 1. 2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2
, i.e. lim m,n→∞ Re γ -1 2 S 1 S 1 u ∞ (-x, d)f p m (d)g p n (x) ds(x)ds(d) = ±1 4π|(z p -a) • ν(a)| + +O(| ln |z p -a||),
where O(| ln |z p -a||) may be large if σ is large near the point a. Hence numerically the second term can weaken the blowup of the first term. Physically, this is the reason why we introduce the coated part of an obstacle, i.e. to avoid or perturb the detection of the obstacle.

The other special property of this example is that the whole boundary shape can be well-reconstructed using only one blowing-up criterion CB. This comes from the special geometric shape and the fact that we have a complete impedance boundary condition. However, in the case of general problems, as for non-convex obstacle, with mixed boundary condition, we need to use multiple blowing-up values to get a satisfactory reconstruction, see subsection 4.3.

κ π σ(a) ≈ 1 | ln((z p -a) • ν(a))| π n 2 2n-1 i,j=1 Im γ -1 2 U ∞ (i, j)F p m (j)G p n (i) (5.2)
for large m, n and small |z p -a|. Example 2. Consider the variable impedance distribution given in Example 1.

We take n = 32. First, let us keep the singularity unchanged by fixing δ 1 and shrink the radius of ∂D p a such that z p → a ∈ ∂D. The computed values of σ(x(t)) with fixed δ 1 = 0.015 and different δ 0 = 1.22, 0.82, 0.12, 0.02 for two wave number κ = 0.6, κ = 0.7 are shown in Figure 5.2. It can be seen from this figure that when δ 0 → 0 (δ 0 = 0.02), the reconstruction results are satisfactory for both wave numbers.

It is interesting to see that when δ 0 is large enough (δ 0 = 1.22), the reconstruction is invalid, even if here we use a strong singularity δ 1 = 0.015. This is reasonable since the approximation to the strong singularity of the fundamental solution contain much error. (3+x 2 ) 2 for δ 0 = 1.22, 0.82, 0.12, 0.02 with two different wave numbers κ = 0.6(left) and κ = 0.7(right), where we fix δ 1 = 0.015.

Finally, we fix ∂D p a near ∂D and take z p tending to ∂D p a . In our configuration, this means δ 1 → 0 for fixed small δ 0 . The tests with δ 0 = 0.002 for two wave number κ = 0.6, 1.0 at different δ 1 = 0.025, 0.020, 0.015, 0.010 are given in Figure 5.3, which show the influence of the singularity. It can be seen that the approximation is sensitive to the wave number κ. The reason is that we ignore the remained term C/| ln |(z p -a)•ν(a)||, where the constant C comes from (4.5). Theoretically, this terms tend to 0 as z p → a. However, this procedure causes a difficulty in approximating the fundamental solution. We expect that the constant C becomes large as κ becomes small. This is naturally related to the following relation, in the 2-dimensional, case between the fundamental solution to the Helmholtz and Laplace equations:

H (1) 0 (κ|x -y|) = - 1 2π ln(|x -y|) -ln κ + O(1)
locally for x, y ∈ R 2 and κ small enough, see [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. This implies that the difference between the fundamental solution behaves as ln κ. So the constant C should have a similar behavior with respect to κ. These remarks on the dependency on wave number κ are also observed in the tests for detecting ∂D. However, we think that this is just a 2D phenomenon. (3+x 2 ) 2 for δ 1 = 0.010, 0.015, 0.025, 0.035 with two different wave numbers κ = 0.6(left) and κ = 1.0(right), where we fix δ 0 = 0.002.

5.3.

Reconstruction of an obstacle with mixed type boundary. Since the main advantage of our inversion method is its ability to identify the full complex obstacle simultaneously, here we consider the numerical behavior of our inversion method acting on a non-convex obstacle with mixed boundary condition.

Example 3. Consider a non-convex obstacle D with the boundary

∂D = {x : x(t) = (x 1 (t), x 2 (t)) = (cos t + 0.65 cos 2t -0.65, 1.5 sin t), t ∈ [0, 2π]}.
Let ∂D be composed of sound-soft part ∂D D for t ∈ [0, 1.42π] and impedance part ∂D I for t ∈ [1.42π, 2π]. In ∂D I , we assume the impedance coefficient σ(x(t)) ≡ 3. We also choose ∂D p a and z p in a special way. For two small constants δ 0 , δ 1 > 0, we take ∂D p a = {y(t) := x(t) + δ 0 × (cos t, sin t) : x(t) ∈ ∂D, t ∈ [0, 2π]}, z p (t) = y(t) + δ 1 × (cos t, sin t) for t ∈ [0, 2π].

(5.3)

In terms of Theorem 3.1, the inversion schemes contain the following three steps:

Step 1. Identify the location of ∂D using (3.18);

Step 2. Distinguish the different parts of ∂D in terms of (3.20);

Step 3. Reconstruction of σ(x) in ∂D I from (3.21).

We present the numerical results with fixed wave number κ = 0.9 and δ 1 = 0.01.

Step 1. we take n = 16 and decrease δ 0 = l × 0.05 by taking l from 20 to 2. The indicator values in (3.18) for ∂D are computed for z p (t) and ∂D p a specified here for every direction t j and different δ 0 . Then we draw the contour line to obtain an approximation of ∂D. As for example 1 and example 2, we choose some appropriate value CB for the stopping rule of l. In the case if the indicator is always less than CB in some direction, we take z p for the initial guess (the largest l) as an approximate location of points in ∂D.

In this case, the situation is different from the examples given in the previous subsections, the reconstruction of the boundary using only one blowing-up value CB is not sufficient, see the left picture in Figure 5.4 for the reconstruction of CB = 3.0, where the kite-shape in red color is the exact obstacle. Enlarging CB can improve the reconstruction along some directions, but the approximation to the whole boundary is still not satisfactory. The reason for this phenomena is that our theoretical results do not guarantee the uniform blowing-up property for all directions. To overcome this difficulty and apply our reconstruction formula, we propose to combine the reconstruction results for different CB together and take the concave hull. Then ∂D can be approximated very well, see Figure 5.4 where we compare the shape given by using one CB only and the one obtained by using and 4 CB values. In Figure 5.5, we show the reconstruction results by using eight and twelve CB values. Precisely, the reconstructions in Figure 5 It can be seen that the reconstruction is satisfactory. This means that using the technic of combining several CB's, the theoretical formulas provide good reconstructions.

We give the indicator value distribution in Figure 5.6 for all directions t with different l at each direction. We can see how the indicator near t = 0.58π, π, 1.42π has some special property, which explains the difficulty of reconstructing these parts shown in Figure 5. Different from the formula (3.18), these two formulas need the boundary shape ∂D which is theoretically obtained from (3.18). Since we can only get some approximation to ∂D numerically in terms of (3.18), it is necessary to check the approximate versions of (3.20) and (3.21) for distinguishing the boundary type and recovering σ(x) on the coating part.

Step 2. We express the quantitative behavior for distinguishing ∂D D and ∂D I by given the indicator distribution. As explained in Step 1, by using different blowing-up criterion in the shape reconstruction, we can get a good approximation to ∂D. In this step, we specify ∂ D ≈ ∂D with an explicit expression given by ∂ D = cos δ *sin δ * sin δ * cos δ * [∂D + {δ * × (cos t, sin t) : t ∈ [0, 2π]}] (

with small constant δ * > 0, for simplicity. In this way, D is no more symmetric with respect to x 1 and the location for corner part also differs from that of ∂D. To check the effect of this domain approximation on (3.20) and (3.21), we also evaluate them using the exact ∂D.

We generate (∂ Dp a , ∂D p a ) from (∂ D, ∂D) and therefore the sequences ({z p }, {z p }) as in (5.3). In this way, we have zp → ã ∈ ∂ D, z p → a ∈ ∂D and δ 0 , δ 1 → 0. In the computation, we take n = 32 and decrease δ 0 = l × 0.05 → 0 by taking l from 10 to 2. The indicator behavior using the same far-field pattern for ∂D and ∂ D with δ * = 0.05 are given in Figure 5.7. Noticing the fact that the sound-soft part corresponds to the parameter t ∈ [0, 1.42π], while the impedance part is related to t ∈ [1.42π, 2π], the above numerical behavior supports (3.21) strongly with a large difference in [0, 1.42π] and [1.42π, 2π], that is, we can distinguish the boundary type in terms of the obvious difference of indicator values when z p approaches to the boundary (for small l), even if the boundary shape is known with a relative error.

Step 3. We compute the impedance coefficient on ∂D I by applying the formulas for ∂D and also for its approximation ∂ D, respectively. The reconstruction behavior for exact ∂D as well as its approximation ∂ D with δ * = 0.05 is shown in Figure 5.8. It can be seen that, for given exact boundary shape ∂D, the theoretical result (3.21) for the impedance is valid (left figure), except near the end points of ∂D I . The rough approximation in this part is reasonable, since this part is near to the sound-soft boundary. Using the approximate domain ∂ D, the impedance can still be captured (right figure). Of course, the reconstruction is less accurate. Due to the nonconvex property of the obstacle and the mixed boundary condition, we think these results are satisfactory. However, if the perturbation for ∂D is very large, then the reconstruction of σ(x(t)) is much contaminated. The reason is due to the sensitivity of the approximate reconstruction to the boundary shape, especially for the corner part of the non-convex domain, noticing that in the formula (3.21), the normal direction ν(a) appears. The perturbation scheme (5.4) moves the position of the corner part by rotation. An inversion result for the impedance with δ * = 0.1 is shown in Figure 5.9. From the left picture of this figure, we see how the corner part of ∂ D with δ * = 0.1 has been much moved from that of ∂D with a relative error almost 10%. We conclude that the theoretical results (3.20) and (3.21) are well supported by our numerical tests even for non-convex domains. If the approximate domain ∂ D is used in these two formulas, then we can still distinguish the boundary type in terms of the obvious blowing-up property of the indicator. However, the quantitative identification of impedance depends on the error level of ∂D.

Fig. 3 . 1 .

 31 Fig. 3.1. Geometric configuration

  which is well-defined by the definition of T. The next proposition gives the relation between E s (x, z) and w + σ(a) (x, z), w + D (x, z) near the point a. Proposition 4.2. If a ∈ ∂D I , then there exist δ(a) > 0 and C > 0 such that

  ) where B + (a, δ(a)) := B(a, δ(a)) ∩ (R 2 \ D) and B(a, δ(a)) is the ball of center a and radius δ(a). Similarly, if a ∈ ∂D D , we obtain (4.5) and (4.6) by replacing w + σ(a) by w + D . Remark 4.3. The estimate of |Re E s (x, z) -Re w + σ(a) (x, z)| is not optimal. We do not need the term | ln |x -a||. The upper bound in (4.6) can be replaced by C α |z -a| -α , for any α > 0 where C α depends on α. But to prove Theorem 3.1 the estimate given in (4.6) is enough. Now we can prove Theorem 3.1 based on these propositions.

  (a))e -2r dr], (4.9) which leads to the first relation in (3.20) by dividing by | ln((za) • ν(a))| s for 0 < (za) • ν(a) < 1 with 0 < s < 1 using (4.5) and (3.16). The representation (3.21) for σ(a) can be gotten from the above relation by dividing by | ln((za) • ν(a))| for 0 < (za) • ν(a) < 1. Step C. Let a ∈ ∂D D . Proposition 4.2 for w + D (x, z) and (3.14) imply the second relation in (3.20), noticing the fact that Imw + D (z, z) = 0.

  z) satisfies the Sommerfeld radiation condition, (4.14) and denote by w s σ(a) (•, z p ) the solution of (4.14) replacing σ(x) by σ(a). Then we have Lemma 4.4. There exist δ(a) > 0 and C(R) > 0 such that

(4. 15 )

 15 and w s σ(a),Γ (•, z) be the solution of (4.15) replacing Φ by Γ. Then we have Lemma 4.5. There exists C > 0 such that |(w s σ(a)w s σ(a),Φ )(x, z)| ≤ C, |(w s σ(a),Φw s σ(a),Γ )(x, z)| ≤ C for z ∈ Ω \ D near D and x ∈ Ω \ D. We define w s,0 σ(a) to be the solution of (4.15) replacing Φ by Γ and the Helmholtz equation by the Laplace equation. Then we have Lemma 4.6. There exists C > 0 such that |(w s σ(a),Γw s,0 σ(a) )(x, z)| ≤ C, for z ∈ Ω \ D near D and x ∈ Ω \ D. Finally, we have Lemma 4.7. There exist C > 0, δ(a) > 0 such that 1. |(Im w s,0 σ(a) -Im w

  satisfies the Sommerfeld radiation condition.(4.16) We need the following lemma.Lemma 4.8. We have the estimate:|ν(z * p )ν(a)| ≤ C|z * p -a|, for z *p near a where C is a positive constant. Proof of Claim 4.8. Take a point b ∈ ∂Ω and connect it to a by a C 3 smooth curve l such that l(0) = a, l(1) = b, l(s) ∈ Ω \ D (s ∈ (0, 1)). By Theorem 7.1 of[START_REF] Honda | The no-response approach and its relation to non-iterative methods for the inverse scattering[END_REF], there is a C 2 strict deformation family {D l(s) a } of a and ∂Ω. That is, each ∂D l(s) a is C 2 diffeomorphic to the unit circle and {D l(s) a } satisfies the following properties:

  ) Hence we have -Im( Ẽsw s σ )(x, z p ) = ∂D F (y, z p )G N (y, x)ds(y)++ ∂B R Im( Ẽsw s σ )(y, z p ) ∂ ∂ν G N (y, x)ds(y)(4.18)whereF (x, z p ) := κσ(x)Re( Ẽsw s σ )(x, z) + [-∂ ∂ν Im∇Φκσ(x)Re∇Φ](x, z) • [ν(z * p )ν(a)]and G N (x, y) is the Green's function of the problem given by the Helmholtz equation in B R \D with Neumann boundary condition on ∂D and Dirichlet boundary condition on ∂B R . The normal ν is oriented outside B \ D. From (4.16), we have

  z) satisfies the Sommerfeld radiation condition. (4.20) From (4.20), we have the representation: R(x, z) = -∂D iκ(σ(y)-σ(a))G σ(a) (y, x)(w s σ + ∂ ∂ν(a) Φ)(y, z)ds(y), for (x, z) ∈ R 2 \D.

(4. 21 )

 21 We define K(x, z) := -( ∂ ∂ν + iκσ(x)) ∂ ∂ν(a) Φ(x, z). From (4.14) we have the following representation w s σ (x, z) = ∂D G σ(a) (y, x)K(y, z)ds(y), hence due to the estimates of the Green's function G σ(a) (x, y) and Φ(x, y), we have |w s σ (x, z)| ≤ ∂D | ln(|y -x|)||z -y| -2 ds(y) ≤ c |x -z| . From (4.21) and the Holder regularity of σ(x), we deduce that |R(x, z)| ≤ c ∂D |y -a| β ln(|y -x|)||z -y| -1 ds(y). From the inequality |y -a| ≤ c(θ)|y -z| for y ∈ ∂D and z ∈ C a,θ ∩ B(a, δ(a)) we have |y -a| β |y -z| ≤ c(θ) β C |y -z| 1-β , which implies |R(x, z)| ≤ ∂D c(θ) β C| ln |y -x|| |y -z| 1-β dy and therefore |R(x, z)| = O(1) for x ∈ R 2 \ D and z ∈ C a,θ ∩ B(0, R). Proofs of Lemma 4.5 and Lemma 4.6.

  |Im R(ξ, η)| < c for ξ ∈ S r and η ∈ C F (a),θ . and |Re R(ξ, η)| < c| ln |ξ -η|| for ξ ∈ S r and η ∈ C F (a),θ .

Fig. 5 . 1 .

 51 Fig. 5.1. Reconstruction of ∂D for variable impedance with κ = 0.6 with different blowing-up constants CB, the values of indicator are small in the directions corresponding to large impedance.

5. 2 .

 2 Reconstruction of σ(x) for known ∂D. For a given a ∈ ∂D, we takez p ∈ R 2 \ Dp a with D ⊂⊂ D p a . By the theoretical result given (3.21), the approximate formula for σ(a) is

Fig. 5 . 2 .

 52 Fig. 5.2. Reconstruction of σ(x) = 2+x 1 x 2(3+x 2 ) 2 for δ 0 = 1.22, 0.82, 0.12, 0.02 with two different wave numbers κ = 0.6(left) and κ = 0.7(right), where we fix δ 1 = 0.015.

Fig. 5 . 3 .

 53 Fig. 5.3. Reconstruction of σ(x) = 2+x 1 x 2(3+x 2 ) 2 for δ 1 = 0.010, 0.015, 0.025, 0.035 with two different wave numbers κ = 0.6(left) and κ = 1.0(right), where we fix δ 0 = 0.002.

Fig. 5 . 4 .

 54 Fig. 5.4. Reconstruction of ∂D with CB = 3.0 (left) and 4 CB ′ s (right) by concave-hull.

Fig. 5 . 5 .

 55 Fig. 5.5. Reconstruction of ∂D by concave-hull using 8 CB ′ s (left) and 12 CB ′ s (right).

  .4 and Figure 5.5 correspond the following CB ′ s values: Figure 5.4 (Right): 4 values-CB = 0.8, 2.5, 3.0, 3.5; Figure 5.5 (Left): 8 values-CB = 0.8, 2.5, 3.0, 3.5, 1.2, 1.5, 5.5, 6.5; Figure 5.5 (Right): 12 values-CB = 0.8, 2.5, 3.0, 3.5, 1.2, 1.5, 5.5, , 6.5, 2.0, 3.2, 4.5, 6.0.

Fig. 5 . 6 .

 56 Fig. 5.6. The indicator values for boundary shape.

Fig. 5 . 7 .

 57 Fig. 5.7. Indicator values for different part of the boundary using exact ∂D (left) and ∂ D with δ * = 0.1 (right).

Fig. 5 . 8 .

 58 Fig. 5.8. Recovery of σ using exact ∂D (left) and approximate ∂ D with δ * = 0.05 (right).

Fig. 5 . 9 .

 59 Fig. 5.9. Approximate ∂ D with δ * = 0.05, 0.1 (left) and reconstruction of σ for different l from approximate ∂ D with δ * = 0.1 (right).

  For y ∈ ∂D and z p ∈ C a,θ we have the inequality |z p -a| ≤ C(θ)|z p -y|. Applying this inequality to z *

	2 |z * p -a|ds(y)

p , enlarging θ if necessary, we have

|( Ẽsw s σ )(x, z p )| ≤ c ∂D (ln |x -y|)|yz p | -1 ds(y) (4.19)

since |z * p -y| ≤ |z p -y|, for y ∈ ∂D and z p near a. Hence for every α > 0 there exists

  1 |z * p -a|. Hence F (y, z p ) ≤ C α |y-z p | -α for y ∈ ∂D. Using the estimate |G N (x, y)| ≤ C α ln |x -y| and choosing α < 1 we deduce from (4.18) that |Im( Ẽsw s σ )(x, z p )| ≤ C α .

We have the first estimate of Lemma 4.4 from (4.19), i.e

|( Ẽsw s σ )(x, z p )| ≤ C ln d(x, ∂D)| • | ln d(z p , ∂D)|.

Now consider the third estimate of Lemma 4.4. We set R(x, z)

Acknowledgments. The authors thank the anonymous referees for their valuable comments and suggestions, which led to a modified version of this paper. The first author is supported by NSFC (No. 10371018). The second author is partially supported by Grant-in-Aid for Scientific research (B)(2) (N0.14340038) of Japan Society for Promotion of Science. The third author is supported by IIRC of Kyung Hee University, Korea, via the SRC/ERC program of MOST/KOSEF (R11-2002-103).