Na Geng

Xiaolan Xie

Zhibin Jiang

Prof Xiaolan Xie

Implementation Strategies of A Contract-based MRI Examination Reservation Process for Stroke Patients

Keywords: healthcare, hospital, capacity allocation, capacity reservation, implementation strategies

Timely imaging examinations are critical for stroke patients due to the potential life threat. We have proposed a contract-based Magnetic Resonance Imaging (MRI) reservation process [1] in order to reduce their waiting time for MRI examinations. Contracted time slots (CTS) are especially reserved for Neural Vascular Department (NVD) treating stroke patients. Patients either wait in a CTS queue for such time slots or are directed

to Regular Time Slot (RTS) reservation. This strategy creates "unlucky" patients having to wait for lengthy RTS reservation. This paper proposes and analyzes other contract implementation strategies called RTS reservation strategies. These strategies reserve RTS for NVD but do not direct patients to regular reservations. Patients all wait in the same queue and are served by either CTS or RTS on a FIFO (First In First Out) basis. We prove that RTS reservation strategies are able to reduce the unused time slots and patient waiting time.

Extensive numerical results are presented to show the benefits of RTS reservation and to compare various RTS reservation strategies.

I. INTRODUCTION

This paper is motivated by our collaboration with a large French university teaching hospital in order to reduce the waiting time of stroke patients treated in Neural Vascular Department (NVD). A stroke is the sudden loss of brain function, which is usually caused by interruption of the blood supply (ischemic stroke) and rupture of a blood vessel (hemorrhagic stroke). The brain is deprived of oxygen and nerve cell death will occur in the area of the brain with no blood flow. The patients will suffer from severe physical and cognitive deficits. It is crucial for stroke patients to have timely imaging examinations in order to have appropriate diagnosis and treatment. However, significant delays are observed because many key examinations rely on expensive and heavily used imaging facilities such as Magnetic Resonance Imaging (MRI). In France, the average waiting time for MRI examinations is about 30~40 days. Stroke patients, just like the other routine patients, have to reserve the time slots via fax or via telephone for the emergency patients. Long waiting time has a negative impact on quality of care and patient service [START_REF] Van Sambeek | Reducing MRI access times by tackling the appointment-scheduling strategy[END_REF].

In order to reduce the waiting time of stroke patients, we proposed in [START_REF] Geng | A Monte Carlo Optimization and Dynamic Programming Approach for Managing MRI Examinations of Stroke Patients[END_REF] a contract-based MRI examination reservation process. The imaging department managing the MRI facilities reserve each week some Contracted Time Slots (CTS) for NVD treating stroke patients.

Stroke patients can be served by either CTS or regular time slots (RTS) in case of arrival surges of stroke patients. CTS decisions and RTS assignment policy determine the efficiency of the reservation process. The former is the number of CTS and its distribution over the time, whereas the latter refers to the policy for assigning incoming patients to RTS. We proposed a method combining stochastic programming model and Markov Decision Process (MDP) to simultaneously determine the two decisions. Structural properties of the optimal RTS assignment policy were established by an average cost MDP model. It is proved that there exists a threshold L i for each day i and the optimal RTS assignment control consists in sending patients to RTS by keeping CTS queue length below L i . Then Monte Carlo approximation and local improvement were used to determine CTS decisions and RTS assignment policy.

Computational results showed that the reservation process can greatly reduce the average waiting time of stroke patients. However, there exist some "unlucky" patients who are directed to RTS and have to wait much longer time than those who wait for CTS.

In order to solve this problem, this paper proposes three new MRI reservation processes. In these processes, RTS is reserved for NVD and added to the list of CTS. CTS and RTS time slots are grouped according to their day of availability and filled by patients. No patient is directed to wait for lengthy RTS. All patients wait in the same queue for MRI time slots and are scheduled to either CTS or RTS in a FIFO (First In First Out) order and according to the release dates of CTS and RTS.

Three RTS reservation policies are proposed in this paper based on patient queue length information and some service ratio criterion. We first prove that RTS reservation policies outperform RTS assignment policy with shorter patient waiting time and less unused MRI time slots. We then perform extensive numerical experiments to compare the performances of different contract implementation strategies.

The rest of this paper is organized as follows. Literature review is described in Section II.

Section III defines and compares different MRI-contract implementation strategies. Formal proofs that RTS reservation strategies improve RTS assignment strategy are given in the Appendix. Section IV presents computational results to show the efficiency of RTS reservation strategies and the impact of different problem parameters. Conclusions and perspectives are given in Section V.

II. LITERATURE REVIEW

Capacity reservation and scheduling of imaging facilities, such as computer tomography (CT) and MRI scanners, have received limited coverage. The two earliest contributions are [START_REF] Walter | A comparison of appointment schedules in a hospital radiology department[END_REF] and [START_REF] Lev | Patient flow analysis and the delivery of radiology service[END_REF]. Simulation was used in [START_REF] Walter | A comparison of appointment schedules in a hospital radiology department[END_REF] to investigate the effect of scheduling rules on patient waiting time and physician idle time for X-rays in a chest radiology department. It was demonstrated in [START_REF] Lev | Patient flow analysis and the delivery of radiology service[END_REF] that more technicians and orderlies do not imply better service quality. To improve the radiology services, the emphasis was the design of the management systems and scheduling techniques.

A new scheduling method was proposed in [START_REF] Van Sambeek | Reducing MRI access times by tackling the appointment-scheduling strategy[END_REF] to identify improvement potentials in order to reduce MRI access time. Finite-horizon dynamic program was used in [START_REF] Green | Managing patient service in a diagnostic medical facility[END_REF] to effectively allocate the expensive imaging diagnosis capacity among several classes of patients during a day. Properties of the optimal policies are identified in order to design the outpatient appointment schedule, and establish dynamic priority rules for admitting patients into services.

[6] proposed a simple approach for dividing the available diagnosis capacity between emergency and inpatients on the one hand and outpatients on the other. MDP was used in [START_REF] Patrick | Dynamic multi-priority patient scheduling for a diagnostic resource[END_REF] for the admission of multi-priority patients on a waiting queue to a diagnostic resource. An approximate dynamic programming approach was proposed to overcome the state space explosion problem. The same problem was solved in [START_REF] Erdelyi | Computing protection level policies for dynamic capacity allocation problems by using stochastic approximation methods[END_REF] with protection level policies by protecting a part of the capacity from the lower priority jobs in order to make it available for the future higher priority jobs. A perturbation analysis technique was proposed to evaluate sample-path gradients with respect to the protection levels. A stochastic approximation approach was used to determine the optimal protection levels.

MDP was used in [START_REF] Kolisch | Providing radiology health care services to stochastic demand of different customer classes[END_REF] to allocate two CT-scanners to three patient groups with different arrival patterns and different cost-structures in order to maximize the total expected reward. [START_REF] Sickinger | A generalized bailey-welch rule and simple tabu search procedure for outpatient appointment scheduling[END_REF] determined the optimal number of outpatients to schedule and the assignment of outpatients to a variable-block/fixed-interval appointment schedule. An MDP approach was proposed in [START_REF] Schutz | Capacity allocation for magnetic resonance imaging scanners[END_REF] to decide whether to accept requests for MRI examinations from patients with different priorities such as inpatients and outpatients. Different examination types, cancellations, no-shows and over-booking, and same-day demand were considered. A continuous-time Markov decision process was used in [START_REF] Schütz | Approximate dynamic programming for capacity allocation in the service industry[END_REF] to model the problem of accepting or rejecting the reservation of different services by different classes of customers. The solution strategy was proposed by using simulation-based approximate dynamic programming (ADP) combined with a discrete event simulation of the service period. Numerical experiments show that the heuristic ADP algorithm performs very well in terms of objective function value, solution time, and memory requirements.

Queueing theory was used in [START_REF] Vasanawala | Accommodation of requests for emergency US and CT: Applications of queueing theory to scheduling of urgent studies[END_REF] to determine the number of schedule slots to open in a radiology department for urgent CT and ultrasonography in order to keep the rescheduling rate of routine studies to accommodate emergencies below a certain level. [START_REF] Shumsky | Dynamic Capacity Management with Substitution[END_REF] examined a multi-period capacity allocation model with upgrading by considering multiple product types and multiple classes of demand. The optimal allocation policy was shown to be a simple two-step algorithm: meet demands with available capacity of the same-class, and then upgrade customers until a protection limit of a class-dependent capacity. Bounds of optimal protection limits were proposed.

Most existing studies focus on medical service capacity allocation and scheduling from the perspective of service provider side except for [START_REF] Geng | A Monte Carlo Optimization and Dynamic Programming Approach for Managing MRI Examinations of Stroke Patients[END_REF][START_REF] Geng | Capacity Reservation and Cancellation of Critical Resources[END_REF][START_REF] Geng | Optimizing contracted resource capacity with two advance cancelation modes[END_REF]. [START_REF] Geng | A Monte Carlo Optimization and Dynamic Programming Approach for Managing MRI Examinations of Stroke Patients[END_REF] proposed contracted based MRI reservation process, and determined CTS decisions and RTS assignment policy by combining stochastic programming model and average-cost MDP. It differs from the previous studies on critical resource scheduling by investigating the problem from a totally different perspective and explores solutions from the perspective of a given class of customers, i.e., stroke patients. CTS are pre-reserved for stroke patients by taking their average waiting time and unused CTS into account. On this basis, [START_REF] Geng | Capacity Reservation and Cancellation of Critical Resources[END_REF][START_REF] Geng | Optimizing contracted resource capacity with two advance cancelation modes[END_REF] considered the possibility of canceling some CTS and proposed average-cost MDP approaches. [START_REF] Geng | Capacity Reservation and Cancellation of Critical Resources[END_REF] determined simultaneously RTS assignment policy and one-day advance cancellation policy for CTS; whereas [START_REF] Geng | Optimizing contracted resource capacity with two advance cancelation modes[END_REF] determined the RTS assignment policy, one-day and two-day advance cancellation policy for CTS. This paper investigates different strategies for implementing an MRI-contract in order to reduce waiting time variation. Advance cancellation of CTS is not considered in this paper.

III. PROBLEM AND STRATEGIES

A. MRI time slot contracting and reservation

In order to reduce the waiting times of stroke patients, we proposed in [START_REF] Geng | A Monte Carlo Optimization and Dynamic Programming Approach for Managing MRI Examinations of Stroke Patients[END_REF] a contract-based MRI examination reservation process. In this reservation process, NVD treating stroke patients is allocated a certain number of time slots every day which is called Contracted Time Slots (CTS). The distribution of CTS over time will be called the contract. When the queue for CTS is too long, stroke patients can also use additional time slots via regular reservation which is called Regular Time Slots (RTS). Patients either wait for MRI time slots in a CTS queue or are directed to RTS via lengthy regular reservation process (30-40 days). This strategy will be called RTS assignment policy. Although the combination of contract and RTS assignment policy greatly reduces the average waiting time of stroke patients, the variation of patients' waiting time is quite large. The "unlucky" patients assigned to RTS have to wait much longer than those waiting for CTS.

In order to improve waiting time distribution, this paper proposes three new implementation strategies. These new strategies still make use of the same contract. With these new strategies, no patient is directed to wait for RTS reservation, and all patients wait in the same queue for MRI time slots. Additional MRI time slots are reserved by NVD through the regular reservation process. The reservation depends on the queue length of patients, and the corresponding time slots will be available in T time periods later. These additional time slots and CTS are pooled together and given to strokes patients on a FIFO basis. These new reservation processes are called RTS reservation policies.

B. System dynamics under contract implementation strategies

In this paper, the following assumptions are made: A1. Only MRI examination is considered and each patient can use either one CTS or one RTS.

A2. Only one class of patients is considered and all patients have equal priority.

A3. The arrival of patients varies during a week but is stationary from one week to another.

Further, the number of arrivals in one day is independent of the arrivals of other days.

A4. The contract is known and remains the same from one week to another. A5. The waiting time for a time slot via regular reservation is a constant T.

Remark 1:

The average waiting time of the stroke patients served by RTS is approximated by a constant T in this paper. In practice, the use of contract-based solution is expected to modify T and reserving more CTS for stroke patients will lead to longer T. However, assumption A5 is reasonable for the following reasons. First, in practice MRI time slots are shared among multiple departments. The proportion of MRI usage by NVD is usually not too high. For example, the hospital we collaborated with has 5 MRI machines shared by 61 medical services of the hospital, many external services and the imaging department itself for various diagnosis, researches and teaching activities. The change of reservation process of NVD alone is expected to have limited impact on the delay of regular reservation. Second, our simulation studies show that change of waiting time for regular reservation due to the use of well-designed contract is quite limited when the overall workload of MRI facilities is reasonable and is for example below 90%. Third, the impact of CTS on T can be further reduced by appropriate advance CTS cancellation strategies not addressed in this paper but investigated in [START_REF] Geng | Capacity Reservation and Cancellation of Critical Resources[END_REF] and [START_REF] Geng | Optimizing contracted resource capacity with two advance cancelation modes[END_REF]. 

The MRI reservation problem is characterized by the following notation:

t index of days with t=1,2, …;
T number of days for obtaining a time slot through regular reservation with T > 1;

a t number of patients arrived in day t. By assumption 3, daily arrivals a t are mutual independent random variables and weekly arrivals (a 7w+1 , a 7w+2 , …, a 7w+7 ) are identically distributed for all w = 0, 1, …. As a result, the arrival process is characterized by probability matrix [r ij ] for i = 1, …, 7 and for all j ≥ 0 with r ij denoting the probability of j arrivals in day i; n t number of CTS of day t. By assumption 4), n t = n t+7 and the contract can be denoted by Four implementation strategies P i with i =0, 1, 2, 3 will be considered in this paper. Each policy P i is associated with the following notation: The capital letter of each notation denotes the cumulative total from 0 to t.

y i,
Notation A t , D i,t
and U i,t will be used. The following notation is also used:

  , , , 1, 1
i t i t i i t q x Y t T t y       (1)
where Y i (t',t) = y i,t' +… + y i,t , Y i (t-T+1,t-1) denotes the outstanding RTS assignment or reservation in day t. If P i is a RTS assignment policy, q i,t corresponds to the CTS queue length and , , , 1

t i t i t i tT x q y        . If P i is a RTS reservation policy, x i,t
is the queue length of waiting patients and q i,t equals with x i,t minus the total number of outstanding RTS reservations.

With a given contract n and a given implementation strategy P i , the sequence of events in each day is as follows. At the beginning of the day, the queue lengths q i,t-1 and x i,t-1 as well as outstanding RTS assignments / reservations y i,t' for t' = t -T to t-

Fig. 1. RTS assignment and RTS reservation

Without loss of generality, we assume that the system starts empty with

  , , , , , , 0 i
t i t i t x q y t    0 (2)
The system dynamics can be characterized by the following:

  , , 1 , i t i t t t i t q q a n y       , for RTS assignment (3)   , , 1 , i t i t t t i t T x x a n y       , for RTS reservation (4)   , , 1 i t t i t t u n q a      , for RTS assignment (5)   , , , 1 i t t i t T i t t u n y x a       , for RTS reservation (6) ,, i t t i t D A x  (7)
The performance measures considered in this paper include: long-run average cost i C , average unused time slots i U , average patient waiting time i W and average queue length i X .

 

 

is the average daily arrival rate. This implies that higher global queue length x i,t implies higher average waiting time.

C. Contract implementation strategies

One RTS assignment policy denoted P 0 and three RTS reservation policies P 1 , P 2 and P 3 are considered in this paper.

P 0 is the optimal RTS assignment policy which is shown in [START_REF] Geng | A Monte Carlo Optimization and Dynamic Programming Approach for Managing MRI Examinations of Stroke Patients[END_REF] to be a policy characterized by a threshold L t associated with each day with L t = L t+7 . This policy keeps the CTS queue length q 0t at the end of each day t below or equal to L t . As a result,

  0, 0, 1 t t t t t y q a n L       (10)     0, 0, 1 min , t t t t t q q a n L      (11)
With a RTS reservation policy P i (i =1, 2, 3), no one is directed to regular reservation and all patients wait in the same patient queue and are served in FIFO order. At the beginning of period t, the length of the patient queue is x i,t-1 , the total number of available time slots is n t +y i,t-T . The number of RTS to reserve y i,t is determined via different methods. P 1 is a RTS reservation policy such that y 1,t = y 0,t . More precisely, P 1 keeps track of an artificial CTS queue length q 0,t as if P 0 were used and it determines y 0,t and hence y 1,t with the artificial queue and relation [START_REF] Sickinger | A generalized bailey-welch rule and simple tabu search procedure for outpatient appointment scheduling[END_REF]. That is, P 1 is similar to P 0 but with additional RTS reserved instead of sending patients directly to RTS. In the following, P 1 is called RTS reservation with artificial queue. P 2 is a RTS reservation policy defined as follows:

    2, 2, 2 1, 1 t t t y x Y t T t L        (12)
and is called RTS reservation with real queue.

      ' 1 1, 1 tT t t t T s s st x Y t T t y L n E a             With     ' 1        tT t t T s s st L L
n E a and the nonnegativity of y t , the above relation leads to [START_REF] Schütz | Approximate dynamic programming for capacity allocation in the service industry[END_REF].

Combining with ( 12) with (1),

2,tt qL  ( 13
)
P 3 is also a RTS reservation policy but is service ratio-oriented and will be called RTS reservation with service ratio . Here the service ratio at the end of a day is defined as the probability that all existing patients are served. Each day t, the number of RTS to reserve is determined such that the service ratio at the end of day t+T is at least , i.e.

 

3,

1 3 t T t T t T t P x a n y         
This policy requires the determination of the probability distribution of 3, 1

t T t T xa     which
depends on the current queue length and all outstanding RTS reservations. In this paper, this probability is determined by Monte Carlo simulation from day t+1 to t+T with a fixed but large number of random samples of patient arrivals.

Remark 2:

From the description of the problem, the dynamics of an RTS reservation policy can be modeled as a Markov chain but its state space is huge and beyond practical limit as the state at the beginning of each day includes not only the length of the patient queue but also the number of RTS reserved in each of the T-1 previous days, i.e. each state corresponds to a T-dimension vector. Direct analysis of the underlying Markov chain is impossible. For this reason, we compare the implementation strategies by sample path comparisons and simulation studies. 

D. Performance comparison of contract implementation strategies

This sub-section compares policies P 0 , P 1 and P 2 under the condition of common arrivals and with the same thresholds L t for all policies. Sample-path comparison is used in the Appendix to establish the results in Table 1.

From The control parameters of the four policies are determined as follows. For any contract n, the optimal thresholds of P 0 can be determined by linear programming as shown in [START_REF] Geng | A Monte Carlo Optimization and Dynamic Programming Approach for Managing MRI Examinations of Stroke Patients[END_REF] and are denoted L 0 . For P 3 , a line search with golden point insertion in interval [0, 1] is used to determine the optimal service ratio . Two parameter settings are considered for P 1 and P 2 .

P 1 L 0 and P 2 L 0 are policies P 1 and P 2 with thresholds L 0 . Another threshold vector is also derived with local search (increasing or decreasing one CTS in one day) to improve P 1 and P 2 by starting from L 0 . Let L 1 and L 2 be the corresponding threshold vectors, and P 1 L 1 and P 2 L 2 be the corresponding policies.

A common randomly generated sample path a of patient arrivals over a time horizon of 10000 days is used in control parameter optimization of policies P 3 , P 1 L 1 and P 2 L 2 . All policies are then evaluated by ten replications of the simulation of the MRI reservation process over 10000 days. Further, with policy 3, at the end of each day, 1000 replications of the on-line simulation over the next T days are performed in order to estimate the patient queue length distribution at the end of day t+T.

B. RTS assignment vs. RTS reservation

This section considers the base case and compares the RTS assignment policy P 0 and RTS reservation policies P 1 L 0 and P 2 L 0 in order to confirm theoretical results of Section IV. The contract under consideration is n 0 = [1,1,1,1,2,0,0] which is the optimal policy with 6 CTS.

Remind that the average weekly demand is 5.74. The optimal thresholds of P 0 are L 0 = Both range between 0 to 24 days with 89.25% served in 10 days and 99.41% in 15 days for P 1 L 0 and 89% in 10 days and 99.35% in 15 days for P 2 L 0 .

Table 2 compares the average daily cost, the percentage of time slots that are not used that will be called unused ratio, the average, standard deviation and maximum of waiting times of the three policies. These results as well as the remaining experimental results confirm the results of Section IV: P 1 and P 2 both improve P 0 , both P 1 and P 2 have smaller waiting time variation and smaller maximal waiting time than P 0 , P 2 has smaller unused ratio but P 1 has shorter waiting time. These results also show that P 1 and P 2 are actually very close with respect to all performance measures as well as waiting time distribution. Remark 3: We have discussed the results with the managers of the collaborating hospital.

They are very interested in the proposed methodology and intend to change the current way.

However, there are several issues they are concerned about: 1) the impacts of this method on the waiting times of other patients; 2) how to implement the strategies if advance CTS cancellation is considered and 3) whether it benefits from using contracts for multiple medical departments. Decisions of real implementation will made based on the results of the above issues. The implementation has also been delayed due to serious financial problems of the hospital, a common situation for most French public hospitals. Cutting cost has become top priority of this hospital while the work of this paper mainly aims at improving service quality to stroke patients. We expect the implementation of the contract-based solutions resumed after the current financial pressure.

C. Sensitivity analysis of unused time slot cost c

This section considers the base case with contract n 0 = [1,1,1,1,2,0,0] and performs sensitivity analysis of the six policies under consideration by changing unused time slot cost c in {1, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, 100}.

Table 3 shows the thresholds L 0 , L 1 and L 2 and the service ratio  of the policies. When c increases, L 0 , L 1 and L 2 increase and  decreases leading to less RTS assignments / reservations and hence smaller unused ratio. Local search reduces the thresholds and hence more RTS reservations and more unused time slots. Figures 3 to 7 compare the average daily cost, the unused ratio, the average value, the standard deviation and the maximum of patients waiting times of the six policies P 0 , P 1 L 0 , P 1 L 1 , P 2 L 0 , P 2 L 2 and P 3 . The following observations can be made.

For average daily cost, (i) all daily cost increases when c increases; (ii) P 3 has the lowest daily cost; (iii) P 0 has the highest daily cost; (iv) P 1 L 0 and P 2 L 0 have almost the same daily cost;

(v) local search improves P 1 and P 2 ; (vi) daily cost are almost the same for all policies when c > 20.

For unused time slots, (i) the unused ratio decreases when c increases, (ii) P 2 L 0 has the lowest unused ratio, (iii) P 0 and P 1 L 0 have the same unused ratio that is slightly higher than P 2 L 0 , (iv) local search leads to significantly higher unused ratio, (v) P 3 often has the highest unused ratio.

For the average waiting times, (i) the average delay increases from half day to over 6 days when c increases, (ii) P 0 has the highest delay, (ii) P 2 L 0 has shorter delay than P 0 but slightly longer delay than P 1 L 0 , (iii) local search significantly reduces the delay, (iv) P 3 often has the shortest waiting time.

For the standard deviation of waiting times, (i) P 0 always has the largest waiting time variation but it decreases when c increases, (ii) the waiting time variations of other policies increase when c increases, (iii) P 1 L 0 and P 2 L 0 have almost the same variation, (iv) local search reduces this variation, (v) P 3 has often the smaller waiting time variation.

For the maximal waiting time, the sensitivity with respect to c is almost the same as that of the standard deviation except that the maximal delay of P 0 remains the same.

D. Sensitivity analysis of RTS delay

This section considers the base case with contract n 0 = [1,1,1,1,2,0,0] and performs sensitivity analysis of the six policies under consideration by changing RTS delay T in {2, 10, 20, 30, 35, 40, 50, 60, 70, 80, 90, 100}.

Table 4 shows the thresholds L 0 , L 1 and L 2 and the service ratio  of the policies. When T increases, L 0 , L 1 and L 2 increase while  is insensitive to the change of T. Local search reduces the thresholds and hence more RTS reservations and more unused time slots. For average daily cost, all daily cost increases when T increases. Except for the case T = 2, P 0 has the highest daily cost, P 3 has the lowest cost, P 1 L 0 and P 2 L 0 has almost the same cost.

Local search does improve P 1 and P 2 and P 1 L 1 and P 2 L 2 have almost the same cost. Cost reduction of RTS reservation over RTS assignment become larger when T increases.

For unused time slots, the unused ratios of P 0 , P 1 L 0 and P 2 L 0 are almost the same and decreases as T increases. P 3 has the highest unused ratio while the unused ratios of P 1 L 1 and P 2 L 2 are smaller than that of P 3 but significantly higher than that of P 1 L 0 and P 2 L 0 especially for large T. There is no clear trend of unused ratios of P 3 , P 1 L 1 , P 2 L 2 .

For the average waiting times, (i) the average delay increases when c increases, (ii) P 0 has the highest delay, (ii) P 1 L 0 and P 2 L 0 has almost the same waiting times, (iii) local search significantly reduces the delay especially for large T, (iv) P 3 often has the shortest waiting time, (v) P 1 L 1 , P 2 L 2 and P 3 are able to limit the average waiting time below 3.5 days even for very large T.

For the standard deviation of waiting times, (i) the waiting time variations increase when T increases, (ii) P 0 always has the largest waiting time variation, (iii) the improvement of RTS reservation policies becomes larger for large T, (iv) P 1 L 0 and P 2 L 0 have almost the same variation, (v) local search significantly reduces this variation, (vi) P 3 has the smaller waiting time variation, (vii) P 1 L 1 , P 2 L 2 and P 3 are able to limit the standard deviation below 3.5 days even for very large T.

For the maximal waiting time, the maximal delay increases as T increases, P 0 has the worse maximal delay equal to T, P 1 L 0 and P 2 L 0 have almost the same maximal delay that is significantly lower than T, P 1 L 1 , P 2 L 2 and P 3 have even smaller maximal delay and are able to limit it below 22 days even for very large T.

E. Sensitivity analysis of patient arrival rate

This section performs the sensitivity analysis of patient arrival rate by multiplying the patient arrival rates of the base case by 1, 5 and 10. This leads to three cases that are called

T P0 P1L0 P1L1 P2L0 P2L2 P3
respectively Low, Medium, and High demand cases. For all cases, the contract n 0 is the optimal contract with RTS assignment determined by the Monte Carlo approach in [START_REF] Geng | A Monte Carlo Optimization and Dynamic Programming Approach for Managing MRI Examinations of Stroke Patients[END_REF] in order to minimize the average daily cost but without limiting the number of CTS. Local search of P 1 and P 2 is not considered here and the impact is similar to policy P 3 .

Table 5 shows the contract n 0 , the thresholds L 0 and the service ratio  of the policies. As the demand increases, n 0 and L 0 increase while  remains stable.

Table 6 gives the performances of different strategies where each performance for each strategy is given in the order of Low, Medium and High demand cases. The followings are observed. When the arrival rate increases, the daily cost increases but all other performance measures decreases. While the reduction of unused ratio and average waiting time is expected.

The reduction of the waiting time variation and maximal waiting time is not straightforward.

RTS reservation strategies do improve the waiting time variations and the maximal waiting time for all cases. This confirms the results of Section III.

Finally, contracts of this section use more CTS than the contract used the previous section which limits the number of CTS to the average weekly arrival rate. With more CTS, the average daily cost, unused ratio and average waiting time are almost the same for all policies.

This suggests that the benefits of RTS reservation strategies with respect to these criteria are significant only when the contracted number of CTS is near or below the weekly arrival rate.

The optimal contract meets this condition if the unused time slot cost c is high.

                  
From definition, we have

Combining the above with (4) for i=1, [START_REF] Shumsky | Dynamic Capacity Management with Substitution[END_REF], the induction assumption and x 0,t-1 ≥q 0,t-1 + y 0,t-T , we have x 0,t ≥x 1,t which completes the proof.  Remark 6: From the proof of Property 1, with the same number of time slots, more patients can be served with RTS reservation P 1 than with RTS assignment P 0 . The deviation, i.e. x 0,t >

x 1,t , happens when n t > q 0,t-1 +a t . In this case, in P 0 , there are not enough patients to fill CTS and patients directed to RTS cannot be redirected. In P x y q x y q t    .

Property 5. The proofs of the above properties need the following basic relations of policy P 2 .

   

* 0, 0 , 0, 0, 2, 2, 2, 2, , , , , , , , , t t
    2, 2, 1 2
1, 1

t t t t u n q a Y t T t          (15)   2, 2, 1 t t t t t y q a n L       (16)       2, 2. 1 2 2 max , , 1, t t t t q x a n Y t T t Y t T t          (17)     2, 2, 1 2, 2 max , 1, t t t t t q q a n y Y t T t         (18)
where ( 15) is from ( 1) and ( 6), ( 16) from ( 4), ( 12) and ( 1), ( 17) from ( 1) and ( 4) , ( 18) from ( 17) and ( 1). The proof of property 7 is similar with that of Property 3 and is neglected.

Proof of Property 4 by induction:

The property clearly holds for t = 0. Assume that it holds up to day t-1. We now prove that it holds for day t.

(A) 0, 2, tt yy  as, by induction and relations [START_REF] Sickinger | A generalized bailey-welch rule and simple tabu search procedure for outpatient appointment scheduling[END_REF] where the first inequality is from (18) and the fact y 2,t ≥0, the second from the induction assumption, the third from the second inequality, ( 13) and [START_REF] Schutz | Capacity allocation for magnetic resonance imaging scanners[END_REF].

Finally x 0,t ≥x x y q u x y q u t t    by induction. It clearly holds for t = 0. Assume that it holds up to day t-1. We now prove that it holds for day t.

(i) Proof of u 0,t =u 2,t . By definition of t*, if u 0,t > 0, then Y 0 (t-T+1,t-1)=0 and, from ( 5) and ( 15), u 0,t =u 2,t . Otherwise, from (15), 0≤u 2,t ≤(n t -q 2,t-1 -a t ) + =u 0,t . Since u 0,t = 0 in this case, u 0,t = u 2,t .

(ii) From ( 16) and ( 10), the induction assumption implies y 0,t = y 2,t .

(iii) Proof of q 0,t = q 2,t . From relations (10), ( 11) and (18)

  0, 0, 1 0, t t t t t q q a n y

      (19)     2, 2, 1 2 2, max , 1, 1
t t t t t q q a n Y t T t y          (20)
If   0 1, 1 0 Y t T t     , by induction and (ii), q 0,t =q 2,t . Otherwise, u 0,t =0, from (5), q 0,t-1 +a t -n t ≥0. By induction assumption, (19) and (20), q 0,t =q 2,t .

(iv) From (i), induction assumption and (ii) and (iii), x 0,t =x

a 7 -

 7 integer vector n = [n 1 , …, n 7 ]; c penalty factor of an unused MRI time slot. It serves as a weighting factor in order to balance waiting time and unused MRI time slots. The cost structure is as follows. Each patient incurs a cost of 1 for each waiting day for either CTS or RTS. Each unused CTS or RTS time slot incurs a cost of c. The goal is to minimize the long run average cost incurred by patient waiting and unused time slots of both CTS and reserved RTS.

  numerical experiments performed to evaluate the performance of several MRI reservation processes. We then perform sensitivity analysis to show how the reservation processes depend on different factors such as T and c, and patient arrival rates.A. Numerical experimentsWe first describe the base case corresponding to our real case study. From the data collected from the NVD of our study, the average numbers of patient arrivals during the week are as follows: {1, 0.89, 0.95, 1.16, 1.53, 0.16, 0.05} which sum up to 5.74 patients per week. The number of patient's arrival each day is assumed to follow a Poisson distribution. The average waiting time for RTS is in the range of 30~40 days with an average of T = 35 days. The weight, c, is set to 15. The base case is then modified to investigate the impact of parameters c, T, and patient arrival rates.

Fig. 2 .

 2 Fig. 2. Waiting time distribution

  of P 1 L 0 and P 2 L 0 are much smoother without unlocking patients.

Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .

 34567 Fig. 3. Daily cost vs. unused time slot cost c

Fig. 8 .Fig. 9 .Fig. 10 .Fig. 11 .Fig. 12 .

 89101112 Fig. 8. Daily cost vs. RTS delay T

  t number of patients directed to RTS or number of RTS reserved at the end of day t;

u i,t number of unused time slots in day t;

x i,t total number of patients waiting for a time slot at the end of day t including those directed to RTS but not yet served. It will be called global queue length; d i,t number of patients having received their time slots and hence left in day t.

Table 1 ,

 1 P 1 and P 2 both improve the RTS assignment policy P 0 . More specifically, the global queue length under P 0 is the longest one, whereas that under P 1 is the shortest, among P 2 uses the least number of additional time slots, whereas P 0 and P 1 use the same number of additional time slots for every period t. Therefore, P 2 has the least number of unused time slots while P 1 and P 0 have the same average unused time slots. With the longest waiting time and highest number of unused time slots, P 0 has the highest cost.

	policies P 0 , P 1 and P 2 . This together with (9) implies 0 W W W 2 1 . In terms of additional
	time slots reserved,

Table 1 :

 1 Performance comparison of policies P 0 , P 1 and P 2

	Performance indicator	Policy ordering
	Global queue length	x 0,t ≥x 2,t ≥x 1,t
	Reservation of additional time slots	y 2,t ≤y 1,t =y 0,t
	Ratio of unused time slots

Table 2 :

 2 Performances of P 0 , P 1 and P 2 for the base case

		Daily cost % unused	Average	Waiting time (days) Standard deviation	Maximum
	P 0	5.03	8.84	4.67	7.37	35
	P 1 L 0	4.74	8.84	4.31	3.73	24
	P 2 L 0	4.77	8.66	4.39	3.77	24

Table 3 :

 3 Thresholds and service ratio v.s. penalty cost c

	c	L 0	L 1	L 2	
	1	{5,5,5,5,4,5,5}	{4,6,2,0,0,5,5}	{5,5,5,5,-18,5,5}	0.69
	10	{6,6,6,6,5,5,6}	{6,6,3,2,2,1,6}	{6,6,6,0,-1,-2,6}	0.34
	15	{6,6,6,6,5,6,6}	{3,4,4,4,2,6,6}	{6,6,2,2,0,6,6}	0.31
	20	{7,6,6,7,6,6,7}	{6,6,5,4,3,4,7}	{5,5,5,4,2,3,7}	0.25
	30	{7,7,7,7,6,7,7}	{7,7,6,5,5,7,7}	{6,8,5,5,3,7,7}	0.22
	40	{8,8,8,8,7,8,8}	{7,7,6,7,6,7,8}	{6,7,7,7,5,8,8}	0.13
	50	{9,9,9,9,8,8,9}	{9,10,9,8,6,8,9}	{9,10,9,8,6,8,9}	0.09
	60	{9,9,9,9,8,9,9}	{8,8,9,9,7,8,9}	{8,8,9,9,7,8,9}	0.09
	70	{10,10,10,10,9,9,10}	{10,11,11,10,8,8,10}	{9,9,10,10,8,9,10}	0.07
	80	{10,10,10,11,10,10,10}	{10,9,11,11,9,10,10}	{10,9,11,11,9,10,10}	0.05
	90	{11,11,11,11,10,11,11}	{11,11,10,12,10,10,11}	{11,11,10,12,10,10,11}	0.05
	100	{12,12,12,12,11,11,12}	{12,12,11,13,11,11,12}	{12,12,11,13,11,11,12}	0.02

Table 4 :

 4 Thresholds and service ratio vs. RTS delays T

	T	L 0	L 1	L 2	
	2	{3,3,3,3,2,3,3}	{3,3,3,3,2,4,4}	{3,3,3,3,2,3,3}	0.31
	10	{4,4,4,4,3,4,4}	{4,3,5,4,3,3,4}	{4,3,5,4,3,3,4}	0.24
	20	{5,5,5,5,4,4,5}	{6,3,5,4,3,3,5}	{5,3,5,4,3,2,5}	0.25
	30	{6,6,6,6,5,5,6}	{4,2,8,3,2,2,6}	{6,6,7,3,2,1,6}	0.27
	35	{6,6,6,6,5,6,6}	{3,4,4,4,2,6,6}	{6,6,2,2,0,6,6}	0.31
	40	{7,6,6,7,6,6,7}	{7,7,6,5,4,2,7}	{7,7,6,2,1,-1,7}	0.31
	50	{7,7,7,7,6,7,7}	{7,8,4,3,1,7,7}	{7,8,7,-1,-1,7,7}	0.31
	60	{8,8,8,8,7,8,8}	{9,9,8,2,1,8,8}	{9,9,8,8,-3,8,8}	0.29
	70	{9,9,9,9,8,8,9}	{9,10,10,3,2,1,9}	{9,10,10,9,8,-5,9}	0.29
	80	{9,9,9,9,8,9,9}	{10,9,9,3,1,9,9}	{10,9,9,9,-5,9,9}	0.31
	90	{10,10,10,10,9,9,10}	{10,10,10,10,2,1,10}	{10,10,10,10,9,-8,10}	0.31
	100	{10,10,10,11,10,10,10}	{11,10,1,11,10,11,10}	{10,10,-9,11,10,11,10}	0.31
	Figures 8 to 12 compare different performance measures of the six policies. The following
	observations can be made.

Table 5 :

 5 Contracts, thresholds and service ratio vs. arrival rate

	Arrival Rate	n 0	L 0	
	Low	{1,1,1,1,3,0,0}	{11,11,11,11,9,10,10}	0.25
	Medium	{5,5,5,6,9,1,0}	{21,21,21,21,19,19,21}	0.27
	High	{10,9,10,12,17,2,1} {31,31,31,31,28,29,30}	0.24

Table 6 :

 6 Performances vs. patient arrival rate The first strategy is built on an artificial queue, the second one use actual patient queue information and the third one is service ratio oriented. Apart from the theoretical analysis, an extensive numerical experiment is performed to show the effectiveness of the RTS reservation strategies.Future research can be pursued in several directions. First, it would be interesting to combine the RTS reservation strategies with advance CTS cancelation. Queueing analysis of RTS reservation strategies is another interesting direction. Other directions are related to the impacts of contracting MRI time slots for the neural vascular department on the waiting times of other patients. Contract design for multiple medical departments and for patients requiring multiple diagnostic facilities is an important challenging issue. The maximal waiting time under P 1 is smaller or equal to T if the thresholds t L The conditions of Property 3 basically assume that patients in CTS queue under P 0 wait no longer than T and a patient in the CTS queue will not be directed to RTS. The last condition holds for any optimal RTS assignment policy. By relaxing the first condition, it is still possible to show that the maximal waiting time of P 1 does not exceed that of P 0 .  Proof of Property 1 by induction: The property clearly holds for t = 0. Assume that it holds for day t-1. We now prove that it holds for day t. First, for P 0 and for day t, q 0,t-1 patients wait for CTS and x 0,t-1 patients wait for the CTS and RTS. The number of new patients is a t . During the day, MIN{n t , q 0,t-1 +a t } patients are served by CTS and y 0,t-T patients by RTS. From the

				Waiting time (days)
	Daily cost	% unused	Average	Standard Dev	Maximum
	P0 4.48 9.78 13.93 18.16 7.71 6.06 2.13 1.14 0.73 2.74 2.12 1.55 35 35 35
	P1L0 4.48 9.75 13.88 18.16 7.71 6.06 2.12 1.13 0.72 2.45 1.46 1.06 15.4 9.9 6.9
	P2L0 4.48 9.75 13.88 18.16 7.71 6.06 2.13 1.14 0.73 2.45 1.47 1.07 15.4 9.9 6.9
	P3 4.47 9.80 13.91 18.32 8.49 6.27 2.09 1.01 0.69 2.40 1.35 1.03 14.9 8.6 6.8

  L t , y 0,t ≤ q 0,t-1 +a t ≤ L t-1 +a t , and The finiteness of x 1,t follows from Property 1.  Proof of Property 3: Consider the last patient arriving in day t under P 1 , i.e. the x 1,t -th patient in the patient queue. Assume by contradiction that its waiting time exceeds T. As patients are

															tt
															0, x	t		0, q	t		* y TL s   where L* is the 1 s s t T s t T a      
	maximum t L for all t. Remark 7: From the proofs of Properties 1 and 2,	tt 0, 11 ss 1, ss uu   	0	. Policy P 1 has less
	unused time slots even though the average underutilization cost is the same.
	served in FIFO order,
															tT 
															1, xa  1 tT t     
															tT 
															0, xa  1 tT t     
	which contradicts Property 1 and concludes the proof.	
	B. RTS assignment vs. RTS reservation with real queue
	Property 4. 	0,	t	0, , , t	0,	t	 	2, t	2, , , t		2,	t		,	0
									0, x	t	0 x  	t  1   t T  1 s s s a  	0, y	s		t  1  s		0, n u s 	s	
									1, x	t	0 x  	tt 11 s ss a   		1, n y s 	s T 		1, u	s	
	As 0, y	1, tt y 	0,    , subtracting the two relations leads to: 0 t
															tt
												0, x	t		1, x	t		0, 11 s ss u   	1, u	s

[START_REF] Geng | A Monte Carlo Optimization and Dynamic Programming Approach for Managing MRI Examinations of Stroke Patients[END_REF]

, as no patients are directed RTS, these extra CTS time slots can be filled by patients that were directed to RTS in P 0 .  Proof of Property 2: By definition of global queue, which proves the property if both x 0,t and x 1,t are finite numbers. The finiteness of x 0,t is true as q 0,t ≤ Consider P 0 . From relation (1) and q 0,t ≤L t , the conditions of the property ensure that all patients in x 0,t have been served by time t + T. As a result,

  and (16),

		0, y	t			    2, 1 t t t t q a n L  0, 1 t t t t q a n L       		2, y	t
	(B) 0, qq  tt 2,	as
								2. t q			2, 1 t q 	t a n t  		
								2, q	t			0, 1 t q 	t a n t  		
						2. t q		min			0, 1 t q 	t a n t  			,	t L			0, q	t

  where the first equality is from (1), the second inequality from (A), (B) and induction, the third equality from[START_REF] Geng | A Monte Carlo Optimization and Dynamic Programming Approach for Managing MRI Examinations of Stroke Patients[END_REF]. 

	2, x	t		2, q	t		 Y t T 2  	1,	 t q 0, 	t		 Y t T 0  	1,	 t x 0, 	t
	Proof of Property 5: We first prove 	0,	t	0, , , , 0, t t	0,	t	 	2,	t	2, , , , 2, t t	2,	t		,	*
	2,t as

  2,t .

	Proof of Property 6: 2, DD  0, tt	is a direct consequence of (7) and 2, xx  0, tt	. Since the
	system is initially empty,			 ,, 1,  i t t i i t D N Y t T U    	which together with Property 4 implies
	2, UU  tt 0,	and	20 UU 		.		20 WW 	flows from (9) and Property 4.	20 CC 	is a
	consequence of Property 4, 2, UU  tt 0,	and
						i C		 11  , , , 11 lim i t i t i tt cu x cU     lim     , i t x      	
	C. RTS reservation with artificial queue vs. with real queue
	From the definition of P 1 and Property 4,
												2, yy  1, tt	(21)
	By induction and by using (4) and (21), we have
												2, xx  1, tt	(22)
	From Properties 2 and 6, 21
	Consider now day t*. 0, * t u		* n q 0, * 1 t t   	*   , 0 t a	 Y t T 0 *  	 1, * 1 0 t   . From (15),
							u	2, * t			* n q 2, * 1 t t  	 a Y t T * 2 * t    	 1, * 1 t 		
												* n q 0, * 1 t t  	 a Y t T * 0 * t    	 1, * 1 t 		
												 tt   0, * 0 0, * * 1, * 1 u Y t T t u      
	Since	* n q 0, * 1 t t  	*   t a	0		, from (10), (11) and (1),	0, * t q 	0	,	t 0, * y 	0	,
	0, * t x		 Y t T 0 *  	 1, * 1 t  . Similar to the proof of (ii), 2, * tt 0, * yy  . From (18), 0
						2, * t q			max		2, * 1 t q 	 a n Y t T * * 2 , * t t     	 1, * 1 t 	
								max		0, * 1 t q 	 a n Y t T * * 0 , * t t     	 1, * 1 t 	
								max		u 	0, * t	,	 Y t T 0 *   	 1, * 1 t   	0
	Combining with (1), 0, * tt 2, * xx 	.	

UU  . From relation (

9

) and (

22

), 21 WW  .

Acknowledgement

The authors acknowledge the support of the National Natural Science Foundation of China (NSFC) under Grant 71131005, 61104173, Specialized Research Fund for the Doctoral Program of Higher Education (20110073120048), and Shanghai Pujiang Program.


Robust TDOA Passive Location Using Interval Analysis and Contractor Programming

Keywords: Passive Location, TDOA, Quimper, Interval Analysis, Constraint Propagation, Contractor Programming, Bounded Error I

This paper presents a new approach for solving non-linear passive location problems. It is based on a high level interval modeling language named Quimper. Whereas classical passive location resolutions do not provide any guarantee of convergence to a solution, interval analysis, constraint propagation and contractor programming allow us to avoid any approximations and any linearization. Besides, Quimper naturally provides guarantees on location and bounded error. TDOA passive location configurations are discussed to prove Quimper's efficiency.

INTRODUCTION

Passive location has been intensively studied in the past years. Numerous devices may actually use passive location including wireless or cellular mobiles [START_REF]An Enhanced Two-Step Least Squared Approach for TDOA/AOA Wireless Location[END_REF][START_REF] Cong | Nonline-Of-Sight Error Mitigation In Mobile Location[END_REF], sonar and radar systems [START_REF] Jiabing | Study on Moving Target Detection to Passive Radar Based on FM Broadcast Transmitter[END_REF], vehicles localization systems. Passive location can take advantage of the growingly number of emitters that are present in environment. It offers a complementary alternative to GPS when indoor, guaranteed or accurate location is required. Moreover, quickly locating an emitter or finding its own location is an important challenge for electronic warfare and unmanned autonomous vehicles (UAV) [START_REF] Drake | Intelligent Sensors, Sensor Networks and Information Processing Conference[END_REF].

Sometimes, decisions have to be taken from results of passive location. When human people are involved, many extra-parameters can be taken into account to choose right actions. For example, symmetrical image location ambiguity in in-line three receivers configuration may be removed from geographical considerations or past observations. In the case of UAVs, there is no human decision maker process. Decisions must be taken from imprecise detectors measurements. That's why exact error estimation is as important as exact location.

This paper intends to show that Quimper software is able to produce an exact bounded error from measurements and, in the same time, to give guarantees on location. Passive location classical approaches give certainly precise results and error estimation. Nevertheless, this error estimation often comes from a probability model which does not care about real-time measured uncertainty. Besides, optimization processes are used to reach a precise result. Theses processes guarantee neither to converge nor to avoid local minimum. That's why an UAV may believe in its location with high level of precision and confidence, while actually being elsewhere. This paper demonstrates the possibility to solve passive location nonlinear equations, without optimization, approximations or linearization.

Second section of this paper gives necessary notions about interval analysis, constraint propagation and contractor to understand Quimper software. Third section is dedicated to a simple TDOA passive location configuration which is solved by Quimper's natural contractors. Last section deals with contractor programming concept and is illustrated by an in-line three receivers configuration.

II. INTERVAL ANALYSIS AND CONSTRAINT PROPAGATION

A. Short Introduction to Interval Analysis

Assume x is a random variable x of IR. An interval approach introduces [x] to represent an interval [x -,x + ] which encloses the support of the probability function of x. Standard Operators +, -, ., / and elementary functions like exp, log, sin or cos can be defined and applied on intervals.

To understand basic idea of interval analysis, consider two variables x and y. Assume that they belong to some prior intervals [x -,x + ] and [y -,y + ]. Addition and difference of intervals can be defined as follow:

] , [ ] , [ ] , [ + + - - + - + - + + = + y x y x y y x x (1) and ] , [ ] , [ ] , [ - + + - + - + - + + = - y x y x y y x x (2)
Multiplication and inversion of intervals are:

)] , , ,

] , [ * ] , [ + + - + + - - - + + - + + - - - + - + - = y x y x y x y x y x y x y x y x y y x x (3) and ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∈ = - + + - + - . 1 , 1 ] , [ 0 ] , [ 1 otherwise x x x x if IR x x (4) max( ), , , , [min(
Let's give some numerical examples of interval computations:

Addition: ] 13 , 2 [ ] 8 , 3 [ ] 5 , 1 [ - = - + Difference: ] 8 , 7 [ ] 8 , 3 [ ] 5 , 1 [ - = - - Multiplication: ] 40 , 15 [ ] 8 , 3 [ * ] 5 , 1 [ - = - Inversion: [ ] ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 3 1 , 8 1 8 , 3 1
Interval arithmetic tools have been developed 1 and it is now possible and simple to build robust programs directly handling uncertainties [START_REF] Jaulin | Applied Interval Analysis with Examples in Parameter and State Estimation[END_REF].

B. Interval Constraint Propagation

Interval constraint propagation combines interval computation [START_REF] Moore | Interval Analysis[END_REF] and constraint propagation [START_REF] Mackworth | Consistency in Networks of Relations[END_REF]. The combination of these two tools has been first presented independently by Clearly [START_REF] Cleary | Logical arithmetic[END_REF] and Davis [START_REF] Davis | Constraint propagation with interval labels[END_REF]. Constraints can be used to contract the prior feasible domains by removing inconsistent values in the domains of the variables.

For example, consider the three variables x,y, and z and assume they belong to prior feasible domains as follow:

[] [ ] [ ] ⎪ ⎩ ⎪ ⎨ ⎧ - ∈ - ∈ ∈ 2 , 5 8 , 3 5 , 1 z y x ( 5
)
Assume also that these three variables are linked by the ADD constraint:

y x z + = (6)
Then, feasible domains for the variables can easily be contracted by taking into account [START_REF] Moore | Interval Analysis[END_REF]. As this constraint is a primitive constraint, a simple projection algorithm exists:

Our numerical example leads to:

1 See for example INRIA COPRIN project. An in-line interval analysis solver is available at http://www-sop.inria.fr/coprin/ael/form.html. Projection is a simple way to contract initial domains. In this case, it is possible because the constraint addition is a trivial constraint. Fortunately, more complex propagation processes have been found. If constraints like equations or inequalities exist between the random variables, a constraint propagation process can produce efficient polynomial algorithms to solve constraints and compute precise solution intervals [START_REF] Jaulin | Applied Interval Analysis with Examples in Parameter and State Estimation[END_REF]. Besides, this process guarantees that these computed intervals enclose all solutions for given initial intervals.

[] [] [ ] [] [ ]

There are many algorithms to generate constraint propagation. Nevertheless, one of the most efficient way is forward and backward algorithm [START_REF] Jaulin | Scientific Computing, Validated Numerics, Interval Methods, Proceedings of SCAN 2000[END_REF]. Quimper Software implements forward and backward algorithm. But, in addition, it focuses on another concept: contractors [START_REF] Chabert | QUIMPER, A Language For Quick Interval Modeling And Programming In A Bounded-Error Context[END_REF].

C. Quimper Software and Natural Contractor

Quimper Software is a recent high-level language for QUick Interval Modeling and Programming in a bounded-ERror context (QUIMPER). It is based on Profil/BIAS 2 and IBEX libraries. It has been created to democratize contractor programming, constraint propagation process and interval analysis. Even if Quimper syntax is simple, manual and examples are available on Quimper's web site 3 .

Quimper's paver algorithm is a generic solver. It takes a list of contractor, an initial box and follows a classical recursion: the contractors are successively called on the current box until either it gets empty or no more contraction could be done. In the latter case, the box is bisected and contractors are called back again.

Contractor's definition is given in [START_REF] Chabert | QUIMPER, A Language For Quick Interval Modeling And Programming In A Bounded-Error Context[END_REF]. Thinking of contractor programming as an extension of constraint programming is valid to the extent that contractors help in modelling the output of a problem. But, fundamentally, there is not such an extension since constraints basically tell the "what" whereas contractors tell the "how".

If no special contractor is given, then Quimper uses natural contractors which result from constraints broken into several trivial constraints. Forward-backward propagation is then applied by contractor to reduce boxes. Nevertheless, Quimper presents specific operators to program ad-hoc contractors. A solver can then be programmed, rather than configured, by combining different contractors. These new contractors faster reduce initial domains. But let first try to solve TDOA hyperbolic equations using Quimper language.

Algorithm ADD_PROJECTION(inout: [z],[x],[y]) [] [] [] [] ( ) [] [] [][] ( ) [] [] [][] ( ) x z y y y z x x y x z z - ∩ = - ∩ = + ∩ = : : : III. TDOA PASSIVE LOCATION WITH QUIMPER SOFTWARE
TDOA passive location requires two stages: the first stage requires TDOA estimation for each pair receiver-receiver. The second stage uses TDOA estimation to build a set of non-linear hyperbolic equations. In this paper, we suppose that the first stage, TDOA estimation, is already done, for example by correlation techniques [START_REF] Carter | Time delay estimation for passive sonar signal processing[END_REF][START_REF] Knapp | The Generalized Correlation Method for Estimation of Time Delay[END_REF].

A. TDOA Hyperbolic Equations

Consider the classical TDOA passive location configuration composed with one emitter and three receivers as sketched on figure 1. Let (x,y) be the location of the emitter, and (x i ,y i ) the location of the receivers. Distance from emitter to receiver i is:

( ) ( ) 2
2 i i i y y x x D - + - = (7) If j i ij D D d - =
, we have:

( ) ( ) ( ) ( ) 2 2 2 2 j j i i ij y y x x y y x x d - + - - - + - = (8)
Let t ij be the measured Time Difference Of Arrival (TDOA) of the signal between receiver i and j. As

ij ij ct d = , hyperbolic
TDOA equations are:

( ) ( ) ( ) ( ) 2 2 2 2 j j i i ij y y x x y y x x ct - + - - - + - = ( 9
)
where c is the speed of the signal and ()(

) ( ) () { } 0 , 2 , 2 , 1 , 1 , 0 , = j i .
Solving these non-linear equations is not a trivial problem. Closed solutions of these equations are known [START_REF] Fang | Simple Solutions for Hyperbolic and Related Position Fixes[END_REF][START_REF] Chan | A Simple and Efficient Es tim ator for Hyperbolic Location[END_REF]. However, these closed forms do not provide any error estimation. Others accurate approaches linearize this set of equations through second order Taylor-series expansion [START_REF] Foy | Postion-Location Solutions by Taylor-Series Estimation[END_REF][START_REF] Torrieri | Statistical Theory of Passive Location Systems[END_REF][START_REF] Abel | A Divide and Conquer Approach to Least-Squares Estimation[END_REF], but need a good initial guess and intensive computation. Our approach based on interval analysis, constraint propagation and contractor programming allow us to avoid any approximations and give true error estimation. So, let us try to model TDOA passive location with Quimper.

B. Quimper Example of TDOA Passive Location

Unknown variables of this problem are x and y, or, more exactly, [x] and [y]. Indeed, more than an accurate result which might be obtained by others approaches, we are looking for domains which might contain solutions and domains which do not contain solutions. Suppose that the bounded receiver sensibility we are using allows us to define some prior feasible domains for these two variables. Let [x] and [y] belong to [-10000,10000] m, which corresponds to a 400 km 2 initial searching area, about four times Paris urban area. Suppose that transmission medium is made of air and that electromagnetic waves are used. The signal speed c is then 3e8 km.s -1 . The three receivers are located at R 0 (-9000, 0) m, R 1 (9000,-9000) m and R 2 (9000,9000) m. Emitter is located at (1540,345) m.

In this simulation, measurement errors are introduced by specifying intervals for t ij , instead of a simple real value. We first extract the exact values of time of arrival from geometrical knowledge. Then, fixed or random values may be added or subtracted to exact t ij value to generate errors. In this first example, we supposed that t ij is known with an uncertainty of plus or minus τ = 500 ns. Initial domain of t ij is then 1 µs wide. This uncertainty corresponds to an analog to digital converter with bad precision and a basic signal correlation. Therefore, ct ij belong to domain [c(t ij -τ), c(t ij + τ)]. We can now write the following Quimper input file (see Script I). A Quimper file is divided in several parts, including constants and variables declarations, and contractors list. It should be noticed that constants may be intervals but that variables must be an intervals.

At the end of this Quimper's file are listed all contractors. Here, we simply put hyperbolic equations: Quimper builds natural contractors corresponding to these constraints. Each contractor removes boxes that are not feasible, one after another. Then, if no more domain contraction is observed, domains are bisected and contractors applied on these new domains. Last contractor is a special one called "Thickness contractor" which uses built-in operator maxdiamGT: isThick removes all boxes that are not thick enough. It is useful if we want solutions to be wrapped into intervals with a size lower than a limit. In this example, this lower size limit is 10m. It means that boxes which contain solutions are contracted by isThick and have a maximum dimension of 10m.

C. TDOA Passive Location Results

Figure 2 shows results extracted from Quimper. Time calculation is about 0.088 s on a Intel Core Duo at 2GHz. Boxes removed by hyperbolas contractors are painted in black and grey. Indiscernible boxes removed by isThick which may contain solutions appear in white, in the center of the figure near exact emitter location (1540,345).

It can be seen that wide domains are removed with only one operation: if a domain do not belong to any hyperbola, then it is no use to search in. The corresponding contractor removes it. All in all, 1745 boxes have been created to solve this problem. Initial searching area is 400 km². From this area, hyperbolas contractors remove 99.990 % and guarantee that there is no solutions in. isThick removes 0.001% of this area which may contain solutions and represents about 37834 m². This result tells that the emitter is located in a 195x195 m² area centred in (1540,342) m. Concretely, it means that we have located somebody in a big stadium, whereas we were looking for him on an area four times greater than Paris.

D. Discussion

These results are fairly remarkable for several reasons: first, no approximations are made, though non-linear equations are quickly solved. These computations may be further accelerated, since parallel-computing may be used by interval analysis. Secondly, removed boxes are guaranteed without any solution: no global minimum has been missed because of a local one and no initial guess is needed. Thirdly, error estimation naturally results from computation and may come from real-time error measurements. If an isThick maxdiamGT parameter of 1 m is chosen instead of 10 m, isThick area is computed at 35097 m² instead of 37834 m². It means that no time consuming small bisections are needed to reach precise results. Time uncertainty can be tuned. If τ equals 50 ns instead of 500 ns and maxdiamGT parameter remains 10 m, then isThick area is 671 m², which represents 0.0002% of the initial are. Location area is now reduced to a big tennis court centred at (1538,344). isThick area can hence be seen as a true characterization of time measurement error.

IV. TDOA AND CONTRACTOR PROGRAMMING

A. Constraint List and Built-In Contractors

Quimper's key idea is that a solver can be programmed rather than configured by combining different contractors. The previous example shows that natural contractors can work concurrently. Let's now design new contractors via Quimper language.

Let's try to compute again the three receivers and one emitter configuration of section III.B, but in a more efficient way (see Script II). Constants and variables are the same as in first script. However, instead of using natural contractors built by Quimper, three stages are proposed to build better ad-hoc contractors. First stage gives a constraint list which takes into account all links between variables. Here, only three equations are used, but others equations may easily be added in this list. For example, equations from others redundant receivers, phase interferometers or radiogoniometers can be introduced in terms of (x,y). No instruction change is needed to handle inconsistent equations system: adding equation to the list is sufficient. Moreover, the greater the number of different constraints is, the faster boxes are removed from initial domain.

Second stage specifies that boxes which do not satisfy intersection of the three constraints can not be solution of our problem. This stage is Quimper's translation of geometrical intersection of the three TDOA hyperbolas. It is done thanks to inter built-in contractor. In the script, hyperbolas(i) is implicitly interpreted as the natural contractor associated to the i th equation. Hence, this stage creates a new contractor clhyp which is the intersection of the three natural contractors. First script's strategy applies natural contractor only one time and one after another. However, repeating several times the same contractor can contract domains in a more efficient way. Therefore, third stage uses the built-in contractor propag to repeat the intersection contractor until no more contraction is observed. Then, isThick contractor is applied. As a last resort, indiscernible domains are bisected and entire process restart on these new domains.

B. Contractor Programming Results

Figure 3 shows computation results of Quimper for the prior example, but with contractor programming as described in Script II. propInter contractor removes boxes and leads to the same result, but has created fewer boxes (1207 instead of 1217). Computation time is about the same than with natural contractors (0.088 s). It means that formula (9) provides a way to build an efficient natural contractor to solve TDOA equations.

C. In-Line TDOA Passive Location Configuration

Consider now the following in-line configuration: emitter is at (8000,7000) m and receivers are located at R 0 (-9000,0) m, R 1 (0,0) m and R 2 (9000,0) m. This configuration is a bad one, because in-line receiver x-location generates y-ambiguity: two symmetrical solutions are indiscernible. Quimper's results are shown figure 4. isThick maxdiamGT parameter is 10 m. Time measurement uncertainty τ is 500 ns. This example demonstrates that interval analysis and contractor programming do not remove any solution. In a pedagogical way, it also shows that this configuration generates bigger error location: initial area is 400 km² and isThick area is about 260557 m² (i.e. 0.06 %). Reasons are time uncertainty, receivers' proximity and in-line topology.

V. CONCLUSION

Quimper is a new powerful tool. To our mind, it is the first time that passive location non-linear equations as hyperbolic TDOA equations are solved without any approximations thanks to interval analysis and contractor programming. Bounded-error estimation naturally results from Quimper's computation. Contractor programming is a way to create adhoc and more efficient contractors. Finally, this high level interval modeling language provides means to mix different passive location approaches. Hyperbolic TDOA, goniometry or phase interferometry may easily be introduced in the Quimper's simple formalism.

  and [z] have been reduced following ADD constraint. No solution has been removed from these intervals.

2 See

 2 Profil/BIAS web site at http://www.ti3.tuharburg.de/Software/PROFILEnglisch.html 3 See IBEX QUIMPER site at http://ibex-lib.org/

Figure 1 .

 1 Figure 1. TDOA classical configuration : three receivers and one emitter. Emitter is located at the intersection of the three hyperbolas for which ti-tj is a constant.

Figure 2 .

 2 Figure 2. Solution of TDOA equations with natural contractors: removed boxes from hyperbola1 (dark grey), hyperbola2 (black), hyperbolas3 (light grey) and isThick (white). Receivers are sketched with crosses: R0 (-9000, 0) m, R1 (9000,-9000) m and R2 (9000,9000) m. Emitter is located at (1540,345) m.

Figure 3 .

 3 Figure 3. Solution of TDOA equations with programmed contractors: removed boxes propInter (light grey) and isThick (white). Receivers are sketched with crosses: R0 (-9000, 0) m, R1 (9000,-9000) m and R2 (9000,9000) m. Emitter is located at (1540,345) m

Figure 4 .

 4 Figure 4. Solution of TDOA equations with programmed contractors: removed boxes from propInter (dark grey) and isThick (white). Receivers are sketched with crosses: R0 (-9000, 0) m, R1 (0,0) m and R2 (9000,0) m. Emitter is located at (8000,7000) m


Olivier Reynet
email: olivier.reynet@no-spam@ensieta.fr

Jean-Christophe Le Lann

Benoît Clément

JOG : une approche haut niveau des systèmes embarqués via Armadeus et Java

Keywords: ARM, FPGA, Java, Linux, capteurs, automate, conception

JOG est un robot mobile indoor dédié à l'enseignement des systèmes embarqués à base de processeur ARM, du système d'exploitation Linux et d'une machine virtuelle Java. JOG permet d'aborder le développement d'un système embarqué complet avec un langage de haut niveau, Java, sans pour autant faire abstraction des aspects matériels du système. Une API Java Jarmadeus a été développée afin de rendre possible ce lien matériel-machine virtuelle. Elle a été testée lors d'un module d'enseignement présenté sous forme de compétition de robots aux étudiants.

Introduction

JOG est un robot mobile indoor dédié à l'enseignement des systèmes embarqués à base de processeur ARM, d'un FPGA, du système d'exploitation Linux et d'une machine virtuelle Java. Ces quatre éléments sont représentatifs des technologies phares des systèmes embarqués actuels et sont disponibles simultanément sur la plateforme open source Armadeus.

Vocations de JOG

JOG a deux vocations initiales. La première vocation de JOG est pédagogique : il s'agit de présenter les concepts des systèmes embarqués à des étudiants de première année d'école d'ingénieur. Ceux-ci ont des notions dans les domaines de l'automatique, de la programmation Java et de l'électromécanique. Cependant, ils ne connaissent pas le langage C et n'ont quasiment aucune notion d'électronique numérique.

La seconde vocation de JOG est liée à la recherche : JOG est une plateforme peu onéreuse et facile à construire. Plusieurs exemplaires de JOG permettront d'étudier des comportements multi-agents, dans les domaines de la modélisation et de la localisation.

Origine de JOG et approche

L'approche des systèmes embarqués adoptée ici est résumée sur la figure 1. À l'origine de cette approche et de la conception de ce robot, plusieurs constats d'enseignants-chercheurs :

1. Pour la majorité des étudiants, il est important de concrétiser les apprentissages au travers de projets motivants, transverses et cohérents avec leurs enseignements précédents. La robotique fournit un terrain idéal tant au niveau de la concrétisation qu'au niveau de la transversalité. C'est pourquoi nous avons présenté le module aux étudiants comme un challenge de robotique.

Figure 1: Approche JOG des systèmes embarqués : les étudiants ne programment qu'en Java mais accèdent néanmoins aux registres des périphériques à la manière des langages de bas niveau.

2. Les étudiants doivent appréhender les concepts les plus importants, sans pour autant se noyer dans l'aspect technique. Il est donc essentiel de leur masquer les difficultés techniques, sans pour autant leur ôter le contrôle et la vision matérielle du système. Concrètement, il nous apparaît important que les étudiants consultent les datasheets des composants électroniques qu'ils ont à programmer, et ne s'appuient pas sur une quelconque interface graphique qui masque totalement le composant. Aussi avons-nous décidé de créer et de fournir aux étudiants une interface logicielle souple constituée d'une API Java utilisant des bibliothèques C et Java Native Interface (JNI). Cette API mandataire préserve la vision bas niveau registre des matériels électroniques tout en simplifiant la lecture et l'écrire sur ces registres.

3. L'importance croissante de la conception conjointe dans les systèmes embarqués rend nécessaire une approche des systèmes embarqués dans toute leur verticalité, et non pas uniquement par le bas ou par le haut niveau comme cela est souvent présenté 1 . Cela signifie que la modélisation des comportements tout comme la programmation des éléments matériels de bas niveau doivent trouver leur place dans l'enseignement des systèmes embarqués.

La section suivante présente en détail les éléments constitutifs matériels et logiciels de JOG. La troisième section décrit l'API Java Jarmadeus dont les étudiants de première année se sont servis pour réaliser le JOG Challenge. La quatrième partie décrit les différents contextes pédagogiques qui mettent en oeuvre le JOG. Grâce à la carte Armadeus APF27 et la carte de développement APF27-devfull, le JOG dispose de nombreuses interfaces, dont deux bus I2C, deux bus SPI, et des liens série. Des capteurs faibles coûts ont été connectés sur le bus I2C (cf. figure 3) : 5 télémètres sonar I2C SRF02 (20-200 cm) et une boussole I2C CMP03. Par ailleurs, la carte de développement comporte un convertisseur analogique numérique Max1027. Quatre détecteurs d'obstacles IR analogiques GP2D120 (0-40cm) ont donc été connectés sur les voies 2,4,5 et 6 du Max1027. Une carte d'adaptation a été developpée, afin de réhausser les tensions des bus pour les capteurs.

Ces capteurs sont tous faibles coût et leurs performances sont donc très moyennes. Pour plus de précision dans l'exploration de l'environnement, un laser Hokuyo peut facilement être interfacé avec la carte Armadeus. Il est également possible de connecter une centrale inertielle, par exemple via le port USB de l'APF27 ou un GPS.

Linux, modules Armadeus et bibliothèque as devices

Le système d'exploitation de la carte Armadeus est Linux (noyau 2.6.29.6). Le système de fichier est au format Ubifs et le Bios est u-boot. Le temps de démarrage est d'environ 10 secondes. Busybox permet ensuite de commander simplement le système d'exploitation. Des modules spécifiques à Armadeus peuvent être chargés pour faciliter l'accès aux périphériques, notamment pour piloter le convertisseur analogique numérique ou pour gérer les GPIO. À noter qu'il existe une bibliothèque as devices au sein du projet Armadeus, bibliothèque qui a inspiré notre approche : il s'agit d'une bibliothèque en C et C++ qui permet d'interfacer aisément tous les périphériques de la carte de développement.

FPGA pour des périphériques à la demande

Le FPGA est utilisé comme un fournisseur de périphérique à la demande5 . En l'occurence, nous avons choisi de créer 4 PWM (2 pour les moteurs, 2 pour des servomoteurs pour des capteurs/actuateurs auxiliaires), un timer hardware, un registre qui pilote 4 diodes et un système de comptage d'impulsions provenant des odomètres des moteurs du JOG. Pour concevoir ce système, nous avons utilisé GHDL, GTKWave, make et les outils de synthèse Xilinx.

L'interaction avec le FPGA a été très appréciée par les étudiants : ceux-ci ont réalisé l'intérêt d'avoir un système capable de fonctionnalités de haut niveau et indépendant du processeur. Le système de comptage des impulsions de l'odomètre a permis aux étudiants de comprendre le fonctionnement interne du FPGA. Le registre pilotant des diodes a renforcé l'acquisition des concepts de l'électronique numérique de base (base 2 et 16, masques, opérations logiques) et l'usage des opérateurs bits à bits.

3 Jarmadeus, une API Java qui facilite l'accès au matériel 3.1 Talker, un intermédiaire à qui parler D'une manière générale, nous avons choisi d'utiliser le concept de mandataire pour créer une interface Java-matériel. Ce mandataire a pour mission de réaliser le dialogue avec des pseudofichiers Linux créés par des modules Linux (par exemple pour le convertisseur analogique numérique), avec un bus I2C ou avec le FPGA. Un package Java formé de classes de type Talker (cf. figure 4) a été développé, dans lequel chaque classe propose une interface simple d'écriture et de lecture sur des systèmes particuliers. Les classes Talker ne suffisent pas pour programmer un composant, car il est nécessaire de transmettre les bonnes commandes, mais elles facilitent grandement la communication avec le composant.

Dialogue ARM-FPGA : memory mapping

Le dialogue avec entre le processeur et le FPGA s'effectue via un memory mapping en langage C. La zone mémoire représentant les registres du FPGA est synchronisée avec le contenu réél des registres. L'accès au FPGA est donc très rapide, même en Java grâce à JNI.

Pseudofichiers Linux et Java

Les pseudo-fichiers Linux sont accessibles comme n'importe quel fichier depuis Java. Cette classe aurait donc pu être écrite par les étudiants. Elle renforce le concept tout fichier Linux. D'une manière plus générale, de nombreux fichiers ont été ouverts et fermés périodiquement au travers de ces mandataires. L'appel aux fonctions C étant momentané, on ne peut pas sauvegarder de descripteur de fichier d'une manière simple tout au long de l'utilisation des fichiers, sans que ce soit apparent pour les étudiants et sans modifier profondément les classes Talker. Néanmoins, l'accès Java au système de fichier étant très rapide, JOG ne souffre pas de cette limitation.

Conclusion

JOG est un outil pédagogique simple de conception, peu onéreux et multifonction. Il a le mérite de présenter aux étudiants un système complet, sans masquer la réalité technologique des composants électroniques et en proposant une API simple d'utilisation fondée sur le principe du mandataire. Il peut être utilisé en traitement du signal, en automatique, en programmation temps réél ou en électronique numérique.

Figure 2 :Figure 3 :

 23 Figure 2: Photographie du robot JOG

Figure 4 :

 4 Figure 4: Package Jarmadeus : intefaces des classes Talker

  Le robot JOG (cf. figure2) est très simple sur le plan mécanique. Il s'agit principalement d'une plateforme Stinger, disponible sur internet au prix de 100 euros environ. Il comporte un châssis en aluminium, 2 moteurs DC, 2 odomètres, 1 roue folle et 1 interupteur. Les batteries 12V permettent d'obtenir une autonomie d'environ 1 heure. La valeur totale d'un JOG est d'environ 800 euros, ce qui reste très raisonable par rapport aux capacités de ce robot.2.2 Carte électronique, capteurs et busL'importance croissante de Linux et des logiciels libres pour l'embarqué rend la carte open source Armadeus 2 relativement incontournable pour des projets innovants dans le domaine. Des concur-rents existent, notamment la carte Fox 3 , ou la Beagle Board 4 . Les avantages de l'Armadeus, en plus d'être un projet open source, sont notamment d'utiliser Buildroot, Busybox, Binutils, Linux et GCC.

2 Anatomie d'un JOG 2.1 Mécanique

De nombreuses sociétés proposent des solutions clés en mains, avec interfaces graphiques exceptionnelles : Lego MINDSTORMS NXT, Robotino de Festo, POB Robotics Suite. . .

http://www.armadeus.com/francais/ et http://www.armadeus.com/wiki/

http://foxlx.acmesystems.it/

http://beagleboard.org/

Cette vision est également portée par le projet POD sur le wiki Armadeus.