Ran Bao

Christian Attiogbe

Benoît Delahaye

Paulin Fournier

Didier Lime

Parametric Statistical Model Checking of UAV Flight Plan

Keywords: UAV, Formal Model, Markov Chain, Parametric statistical model checking

Unmanned Aerial Vehicles (UAV) are now widespread in our society and are often used in a context where they can put people at risk. Studying their reliability, in particular in the context of flight above a crowd, thus becomes a necessity. In this paper, we study the modeling and analysis of UAV in the context of their flight plan. To this purpose, we build a parametric probabilistic model of the UAV and use it, as well as a given flight plan, in order to model its trajectory. This model takes into account parameters such as potential filter or sensor (like GPS) failure as well as wind force and direction. Because of the nature and complexity of the successive obtained models, their exact verification using tools such as PRISM or PARAM is impossible. We therefore develop a new approximation method, called Parametric Statistical Model Checking, in order to compute failure probabilities. This method has been implemented in a prototype tool, which we use to resolve complex issues in a practical case study.

Introduction

Unmanned Aerial Vehicles (UAV) are more and more present in our lives through entertainment or industrial activities. They can be dangerous for their environment, for instance in case of a failure when an UAV (aka a drone) is flying above a crowd. Unfortunately until today, there does not exist any kind of UAV regulation around the world. Only some recommendations are used; for instance in order to avoid accidents in case of malfunctioning, a drone should never fly above a crowd.

In this context, we are working with PIXIEL group to build a reliable UAV control system. PIXIEL group is a company expert in safety drones and public performances including UAVs. For example, PIXIEL is in particular known for developing a public performance in the French entertainment park called "Puy Supported by PIXIEL and Association Nationale Recherche Technologie (ANRT) du Fou" that includes both human actors and drones. The company is strongly attached to the safety of the public. Therefore, ensuring that its UAV systems are secure for humans during the performances is a priority. As for the current practices, the performances including UAVs are only allowed to occur when the weather is sunny and when the area above which the UAVs fly is unauthorized for actors and public. However, there is no certification proving that the UAVs always follow their intended flight plan.

The management of performances indeed requires to pay close attention to the drone trajectory computation as well as to the accuracy of the measurements concerning its immediate position in space and its movements. However, a rigorous study is necessary to ensure reliability of the drone control system, for instance by decreasing the risks of failure using the appropriate tuning of the drone flying parameters which impact the computation of its trajectory. Accordingly, the questions are how to prove that the UAV failure probability is low and which parameters have to be taken into account to ensure human safety during performances including UAVs.

High-quality aircrafts such as Hexarotors can easily avoid the majority of minor failures related to hardware because they can fly with only five motors and the probability of concurrent failure of more than two motors is in general negligible. In the same way, in case of battery failure, the UAV is able to land down on a specified area without any safety issue for the environment as long as it is situated in a safe zone where humans are not endangered. However, software failure may be a lot more problematic and complex to study. In this case, the UAV behavior might become unpredictable. One critical issue in this context is the potential inaccuracy of position estimation in drone systems, either as a result of inaccurate sensor measurements or of misinterpretation of data coming from those sensors. Besides aircraft system failure consideration, there is also a far more critical aspect to take into account: the weather environment. Therefore, a general approach to improve UAV safety is to study the impact of inaccuracy in position measurements on the resulting flight path compared to a given, fixed, flight plan while taking into account weather conditions.

There are many works dedicated to the UAV domain. In [20] Koppány Máthé and Lucian Buşoniu basically explain the functioning of a drone. UAV movement recognition is studied in [10]. Automatic landing on target is described in [START_REF] Kyristsis | Towards Autonomous Modular UAV Missions: The Detection, Geo-Location and Landing Paradigm[END_REF] and monitoring and conservation are dealt with in [START_REF] Gonzalez | Unmanned aerial vehicles (uavs) and artificial intelligence revolutionizing wildlife monitoring and conservation[END_REF]. Some works also try to detect breakdowns and malfunctions that can impact drones. We can mention inter alia, the detection of communication errors in a multi-drone framework studied in [START_REF] Heredia | Multi-Unmanned Aerial Vehicle (UAV) Cooperative Fault Detection Employing Differential Global Positioning (DGPS), Inertial and Vision Sensors[END_REF] or the development of a basic diagnosis model for solving system issues in [START_REF] Freddi | A model-based fault diagnosis system for unmanned aerial vehicles[END_REF]. Our work is closer to this second category of topics. However, to the best of our knowledge, there are no existing works on the parametric study of the impact of component inaccuracy on UAV trajectory. In [START_REF] Chang | Secure estimation based Kalman Filter for cyber-physical systems against sensor attacks[END_REF][START_REF] Zhou | Efficient path planning algorithms in reach-avoid problems[END_REF] the authors study through the secure estimation problems how to estimate the true states of an UAV system when the measurements from sensors are corrupted, for instance by attackers. In their work, these authors reformulate the estimation problem into the error correction problem and then they use the successively observed measurement anomalies to reconstruct the correct states of the system. While the used techniques are completely different, the objectives of avoiding bad states is similar to ours, in avoiding the states reaching bad security zones.

The purpose of our work is therefore to provide means to study the reliability of UAVs in the context of a given flight plan. In order to do that, we have to build a formal (mathematical) model which will allow us [START_REF] Alur | Parametric real-time reasoning[END_REF] to analyze the drone system and detect the most important parameters, and (2) to tune those parameters in order to reduce the system failure probability. To this intent, we thoroughly study the UAV system, formalize it and analyze it with using parametric probabilistic methods. Among the components of a drone system, we particularly focus on the Flight Control System (FCS), which is responsible for computing estimations of the UAV position during its flight in order to adapt its trajectory to a given predefined flight plan. We therefore build a formal model of the flight controler in terms of parametric probabilistic models that takes into account the potential inaccuracy of the position estimation. Since UAVs are particularly sensitive to the weather environment (and in particular to wind conditions), we also enhance our model in order to take into account potential wind perturbations. Since wind force can drastically vary from one point of a given flight plan to another, we also use parameters to encode the wind force and allow our model to adapt to particular weather conditions.

The contributions of this paper are:

a method to build a parametric model of UAV systems; the parameters can then be finely tuned until reaching values that ensure defined safety thresholds; a parametric statistical model checking technique; this enables us to formally analyze the parametric models build for the drones. Indeed because of the complexity of the built models, tools such as PRISM [START_REF] Kwiatkowska | PRISM: Probabilistic Model Checking for Performance and Reliability Analysis[END_REF][START_REF] Kwiatkowska | PRISM 4.0: Verification of probabilistic real-time systems[END_REF] and PARAM [START_REF] Hahn | PARAM: A Model Checker for Parametric Markov Models[END_REF] were limited for their analysis. an illustration of the use of our method on a complex industrial case study.

The paper is organized as follows. In Section 2 we provide the essential background to understand UAV functioning and then we build a formal model that support their behaviours. Section 3 is an introduction to parametric Markov chains and Statistical Model Checking. Implementations of the models and experimentations are presented in Section 4; finally Section 5 draws conclusions and further work.

Building a Formal Model of UAV

In this section, we present our method to build the UAV model. Recall that we are interested in studying UAV safety, i.e. studying the probability that a UAV encounters dangerous situations. These situations are of two kinds: either the UAV can stop flying and fall, or it can enter a "forbidden" zone were it endangers humans. As explained earlier, professional UAV can handle the falling risk through material redundancy. Moreover, as long as a UAV stays in a "safe" zone, it will not endanger human even in case of falling. The aim of our model is therefore to evaluate the probability for a UAV to enter a "forbidden" zone.

We start by explaining how the zones are computed with respect to the given flight plan. We then show how the UAV software can be decomposed into components and focus on the most important ones. Finally, we detail how the formal (mathematical) models for the important components are built and present the resulting global model.

Safety zones

Fig. 1: Safety zones

In the context of software, considerations in airborne systems and equipment certification (named DO-178C) defined five levels of safety zones, the most secure being Zone 1 and the most dangerous being Zone 5. These zones are characterized by their distance from the intended flight plan, as shown in Figure 1.

The size of each safety zone is not definitely fixed; it can be defined for a specific requirement or for a given application. In practice the safety zones are specifically defined for a flight environment and for a given flight plan. The main principle is that no human should be present in Zones 1 to 3, while a few people can be present in Zone 4 and most people can be present in Zone 5. As a consequence, the probability that the UAV endangers humans is directly proportional to the probability that it enters Zones 4 or 5. In the following of the paper, our target will therefore be to compute this probability.

Drone components

We now move to the decomposition of the UAV hardware and UAV software into components and introduce the most important component in the UAV system: the flight controller (FC). The FC is responsible for collecting data from various sensors, using this data to compute the precise position and attitude of the drone and adjust the attitude in order to follow the given flight plan to the best of its ability.

Fig. 2: Attitude coordinates

Notice the difference between position and attitude: while the position of the UAV is defined by 3-dimensional coordinates x, y and z, its attitude is the collection of yaw , pitch and roll measurements for the UAV compared to the vertical (see Figure 2). The attitude allows to control the movement Fig. 3: Flight Control Overview of the UAV: by controlling the speed of each motor, one can control which motor will be the highest, and hence control the direction the UAV will fly to.

Flight Controller As explained above, the FC is the central component in any UAV as it is responsible for collecting data from sensors and translating them to the UAV attitude. An overview of the FC of an UAV is given in Figure 3. Remark that the FC can be linked to components responsible for communicating with a remote control. While these components are necessary in order to allow a pilot to take over when the automatic flight mode of the UAV fails, we will consider in the following that this is not the case and that the UAV we study are always in automatic flight mode.

As one can see from Figure 3, the intuitive behavior of the FC is as follows. The filter uses sensors measurements in order to compute the current drone position and attitude. Since the data can be noisy and inaccurate, the filter uses complex algorithms in order to clean the noises in the measurements and compute a realistic position and attitude. Remark that in some cases, the filter can itself introduce inaccuracy in the computed position and attitude, which can be problematic. Once the estimated current position and attitude are computed, the Proportional Integral Derivative (PID) uses this information to compute the local trajectory that the drone has to follow in order to be as close as possible to its intended flight plan. This local trajectory is then transformed into a new value for the attitude of the drone. Finally, Modulation transforms this attitude into signal to the Electronic Speed Controller (ESC) which is responsible for controlling each motor's speed.

Recall that we are interested in computing the probability that a UAV enters a forbidden zone while following its flight plan. By construction, as long as the position and attitude measurements are perfect, there is no reason why the UAV should deviate from its intended trajectory, and therefore the probability that it enters a forbidden zone is null. However, as explained above, the data gathered from sensors can be noisy and inaccuracy can sometimes be introduced through filtering. In this case, the estimated position and attitude of the UAV can be faulty, resulting in a deviation from the intended flight plan and potentially leading to a forbidden zone. It is therefore of paramount importance to study how the filters work and to take into account in our formal model the potential inaccuracy of position and attitude measurement.

Filter The role of the filter is to use sensors measurements in order to compute the UAV position and attitude with the highest possible precision. However, the high precision comes with a cost in terms of complexity: in order to gain precision, filters have to run complex algorithms which takes time. As a consequence, the most precise filters are also the slowest, which implies that the position can be estimated less often, which itself results in inaccuracy.

There exists a large amount of filters in UAV industry, among which one can find Extended Kalman Filter (EKF) [START_REF] Sabatelli | A Double-Stage Kalman Filter for Orientation Tracking With an Integrated Processor in 9-D IMU[END_REF], Explicit Complement Filter [START_REF] Euston | A complementary filter for attitude estimation of a fixed-wing UAV[END_REF], Gradient Descent [START_REF] Madgwick | Vision and Control for UAVs: A Survey of General Methods and of Inexpensive Platforms for Infrastructure Inspection[END_REF], Conjugate Gradient, and a more accurate but slower filter: Unscented Kalman Filter (UKF) [START_REF] De Marina | UAV Attitude Estimation Using Unscented Kalman Filter and TRIAD[END_REF], etc. Usually, researchers use EKF as a fundamental to compare to other kinds of filters and explain precision and speed differences. All filters improve their accuracy during the flight through training, in particular by recording recurrent noises and correcting them. However, this training is only valid through a single flight and is lost as soon as the UAV lands.

Since the accuracy of the estimated position and attitude is of paramount importance for computing the probability of entering a forbidden zone, and since the choice of filter has a direct impact on this measurement, we chose to implement this accuracy as a parameter of our model. This will be explained in more details in Section 2.4. We use a flow diagram to present our global approach for formalizing the UAV functioning (See Figure 4). After a step where the filter computation reflects the precision of position and attitude estimation, we consider the computations of the probabilities to reach the given safety zones in the next time-step; accordingly, the idea is to adapt the next attitude according to the original flight plan in order to be more secure. The last step allows to incorporate wind perturbations and compute the next UAV position.

Formal Model of the UAV in its Environment

As explained above, the filter is one of the most crucial components and its ability to estimate the UAV position precisely has a huge impact on the probability of reaching a forbidden zone. For this reason, we choose to represent the accuracy of the estimated position of the UAV (therefore including both sensor measurements and filter correction) as a parameter of our model. In the following, we show how the next position of the UAV is computed according to the current estimated position, and how errors in the estimation can lead to the drone entering forbidden zones.

Computation of the next position We now explain how the next position is computed according to the estimated current position. In particular, we show that inaccuracy in the estimation can lead the UAV to entering a forbidden zone.

Fig. 5: Issue on drone location and misleading positions

For the sake of simplicity, we assume here that the UAV moves in 2 dimensions only and that inaccuracy only occurs on one of them. Figure 5 illustrates the situation. Assume that the intended flight plan consists in going from point A to point B. Assume also that the current position of the UAV is exactly on A but that the estimated position (taking into account sensors and filter inaccuracy) is on A . As a consequence, the PID will try to correct the current deviation by changing the angle of the UAV in order to lead it back to B. However, since the UAV is really on A, the correction will instead lead the UAV to a position B , in the forbidden zone. Fortunately, the position estimation takes place several times between A and B, according to the filter frequency f . Therefore, a new position will be estimated before reaching B , hopefully with a better accuracy, which will allow the PID to again correct the trajectory. We should also take into account that the speed of the UAV is also computed according to the flight plan, which precises the remaining time and distance before the next checkpoint. We now show how we can compute the safety zone where the UAV ends before the position is estimated again. In Figure 5, this zone is represented by the distance S n .

Let S answer be the distance that the UAV covers before a new estimation of the position. Let V be the velocity of the UAV, which is computed by the PID in order to reach B on time, i.e. in precisely T time units. We therefore have V = A B/T , and S answer = V /f = (A B)/(T * f ). Finally, AA /A B = S n /S answer , and therefore

S n = AA T * f .
Remark that the resulting deviation is directly proportional to AA /f , hence the necessity to take into account the trade-off between accuracy and filter speed in order to optimize the probability of never entering any dangerous zone.

Taking into account wind perturbations follows a similar computation than the one presented above. This allows us to incorporate wind parameters as well in our model. The global model of the UAV flight control system is depicted in Figure 6. The purpose of this model is to represent the computations taking place in the FCS in order to adapt the UAV trajectory to the intended flight plan according to inaccurate position and attitude estimations as well as wind perturbations. In this model, the exact position of the UAV is encoded using 3-d coordinates. These coordinates are then compared to the intended flight plan in order to decide to which safety zone they belong. As soon as the UAV reaches one of the forbidden zones (4 or 5), the computation stops.

Resulting Global Model

The model uses several probabilistic parameters. Parameters Filter-Proba1, FilterProba2, FilterProba3, FilterProba4 and FilterProba5 represent the accuracy of the position and attitude estimation by both the filter and the sensors. The resulting probabilistic choice depicted in the box labelled Filter Computation therefore dictates the distance between the exact and estimated position of the UAV. This choice is followed by a computation in the box labelled Safety Zone Computation that computes the exact coordinates of the next position of the drone and allows to decide the safety zone to which this position belongs. When the wind is not taken into account, the result of this computation is enough to decide whether the model should pursue its execution. When the wind is taken into account, another step follows, depicted in the box labelled Wind Computation, where other probabilistic parameters are used in order to decide the wind strength (we assume that the direction is constant) and a new position taking into account these perturbations is computed. Finally, the zone to which this last position belongs is computed and, depending on whether this zone is safe, the model goes on to another position estimation.

Remark that the filter frequency and the position and distance of checkpoints in the flight plan are given as inputs to the model. The position of checkpoints in the flight plan allows to compute the required UAV speed, while the frequency of the filter allows to fix the number of position estimations that will happen in a given flight plan (i.e. the number of loops the model goes through, at most).

Parametric Statistical Model Checking

As explained above, we have developed a parametric probabilistic model in order to represent the behaviour of our UAV according to a given flight plan. We now introduce the necessary theory to formally compute the probabilities of a given UAV entering a forbidden zone in the context of its flight plan. We start by recalling a classical verification technique called Statistical Model Checking (SMC), then introduce the modeling formalism we use: parametric Markov Chains (pMCs) and finally show how SMC can be adapted to this formalism.

Standard Statistical Model Checking

Recall that a Markov Chain (MC) is a purely probabilistic model M = (S, s 0 , P ), where S is a set of states, s 0 ∈ S is the initial state, and P : S × S → [0, 1] is a probabilistic transition function that, given a pair of states (s 1 , s 2 ), yields the probability of moving from s 1 to s 2 .

Given a MC M, one can define a probability measure on the infinite executions of M using a standard construction based on the σ-algebra of cylinders.

A run of a MC is a sequence of states s 0 , s 1 . . . such that for all i, P (s i , s i+1 ) > 0. Given a finite run ρ = s 0 s 1 . . . s l , its length, written |ρ| represents the number of transitions it goes through (including repetitions). Here |ρ| = l. We write Γ M (l) (or simply Γ (l) when M is clear from the context) for the set of all finite runs of length l, and Γ M for all finite runs i.e. Γ M = ∪ l∈N Γ M (l). As usual we define the probability measure, written P M on runs based on the sigma-algebra of cylinders (see e.g. [START_REF] Baier | Principles of model checking[END_REF]). This gives us that for any finite run ρ = s 0 s 1 . . . s n , P M (ρ) = n i=1 P (s i-1 , s i ). In the rest of the section, we only consider f inite runs. Given a reward function r : Γ (l) → R, we write E l M (r) for the expected value of r on the runs of length l of a given MC M.

Statistical Model Checking [START_REF] Sen | On statistical model checking of stochastic systems[END_REF] is an approximation technique that allows to compute an estimation of the probability that a purely probabilistic systems satisfies a given property 4 . In particular, the Monte Carlo technique uses samples of the runs of length l, Γ (l), of a given Markoc chain M in order to estimate the probability that M satisfies a given bounded linear property. It can also be used for approximating the expected value of a given reward function r on the runs Γ (l) of M. In order to provide some intuition, we briefly recall how standard Monte Carlo analysis works in the context of statistical model checking of MC. In this context, a set of n samples of the runs of the MC. These runs are generated at random using the probability distribution define through the Markov chain. Each of these samples is evaluated, yielding a reward value according to the reward function r. According to the law of large numbers (see e.g. [START_REF] Rubinstein | Simulation and the Monte Carlo method[END_REF]), the mean value of the samples provides a good estimator for the expected value of the reward function r on the runs of the given MC. Moreover, the central limit theorem provides a confidence interval that only depends on the number of samples (provided this number is large enough).

Parametric Markov Chains(pMC)

Markov Chains are inadequate in the context of drone flight plan analysis. Indeed, the models we develop in this context are subject to uncertainties that we model using parameters, such as precision of the position and attitude estimations and wind strength. The resulting models are therefore not purely probabilistic since they contain parameters. As a consequence, we need to use a more expressive type of model that allows to take into account probabilistic parameters, such as Parametric Markov Chains (see e.g. [START_REF] Alur | Parametric real-time reasoning[END_REF]).

A pMC is a tuple M = (S, s 0 , P, X) such that S is a finite set of states, s 0 ∈ S is the initial state, X is a finite set of parameters, and P : S × S → P oly(X) is a parametric transition probability function, expressed as a polynomial on X. A parameter valuation is a function v : X → [0, 1] that assigns values to parameters. A parameter valuation v is valid w.r.t. a given pMC M if, when replacing parameters with their assigned values, the resulting object is a MC (i.e. the outgoing probabilities of all states sum up to 1). If v is a valid parameter valuation with respect to M, the resulting Markov chain is written M .

Given a pMC M, a run ρ of M is a sequence of states s 0 s 1 . . . such that for all i ≥ 0, P (s i , s i+1 ) = 0 (i.e. the probability is either a strictly positive real constant or a function of the parameters). As for MCs, we write Γ M (l) for the set of all finite runs of length l and Γ M for the set of all finite runs.

Observe that for any valid parameter valuation v, Γ M v (l) ⊆ Γ M (l) since v may assign 0 to some transition probabilities.

Parametric SMC

As it is, standard SMC cannot be used in the context of pMC because of their parametric nature. Indeed, we cannot produce samples according to the parametric transition probabilities. Luckily, the underlying theory used in SMC can be extended in order to take into account parameters. The method we propose in the following is in line with a technique called importance sampling (see [START_REF] Rubinstein | Simulation and the Monte Carlo method[END_REF] for a description). The purpose of this technique is to sample a stochastic system using a chosen probability distribution (which is not the original distribution present in this system) and "compensate" the results using a likelihood ratio in order to estimate a measure according to the original distribution. In the context of SMC, importance sampling has mainly been used in order to estimate the probability of rare events [START_REF] Barbot | Coupling and importance sampling for statistical model checking[END_REF] and/or to reduce the number of required samples in order to obtain a given level of guarantee [14]. It has also been used in the context of parametric continuous-time Markov chains in order to estimate the value of a given objective function on the whole parameter space while using a reduced number of samples [START_REF] Bortolussi | Smoothed model checking for uncertain continuous-time markov chains[END_REF]. However, to the best of our knowledge, importance sampling has never been used in order to produce symbolic functions of the parameters as we do here.

The intuition of the method we propose here is to fix the transition probabilities to an arbitrary function f , which we call normalization function, and to use these transition probabilities in order to produce samples of the pMC M. However, instead of evaluating the obtained runs by directly using the desired reward function r, we define a new (parametric) reward function r that takes into account the parametric transition probabilities. We show that, under any parameter valuation v, the evaluation of the mean value of r on the set of samples is a good estimator for the expected value of the reward r on M v . The central limit theorem (see e.g. [START_REF] Rubinstein | Simulation and the Monte Carlo method[END_REF]) also allows to produce parametric confidence intervals, but we do not go into details here (see [START_REF] Delahaye | Statistical model checking for parameterized models[END_REF] for more details on this topic).

Remark. The choice of the normalization function is crucial. In particular, the results presented below require that the graph structure of the MC obtained with this normalization function is identical to the graph structure of the MC obtained using the chosen parameter valuation. This is discussed in more details in [START_REF] Delahaye | Statistical model checking for parameterized models[END_REF]. In the following, we only consider parameter valuations that assign non-zero probability to parameterized transitions. Since we use the uniform normalization function, the graph structures of the obtained MCs are indeed identical, which ensures that the results presented below hold as expected.

Let P a : Γ M → P oly(X) be a parametric reward function. For any valid valuation v and any run ρ ∈ Γ M v we have P M v (ρ) = P a(ρ)(v).

Given any valid normalization function f and any run ρ ∈ Γ M , let parametric reward function r be r (ρ) = P a(ρ) P M f (ρ) r(ρ).

We now prove that the expected values are equal. Let ρ ∈ Γ M f (l) be a random sample of M f and let Y be the random variable defined as follows Y = r (ρ). The following computation shows that, under any valid parameter valuation v such that M f and M v have the same structure, we have

E(Y )(v) = E l M v (r). E(Y )(v) = ρ∈Γ M f (l) P M f (ρ)y(ρ) (v) = ρ∈Γ M f (l) P M f (ρ) P a(ρ) P M f (ρ) r(ρ) (v) = ρ∈Γ M f (l) P a(ρ)(v)r(ρ) = ρ∈Γ M f (l) P M f (ρ)r(ρ) = ρ∈Γ M v (l) P M f (ρ)r(ρ) = E l M v (r)
Our adaptation of the Monte Carlo technique for pMC is thus to estimate the expected value of Y in order to obtain a good estimator for the expectation of r. Let ρ 1 , . . . , ρ n be a set of n runs of length l of M f . Let Y i be the random variable with values in P oly(X) such that Y i = r (ρ i ). Notice that the Y i are independent copies of the random variable Y . Y i are therefore independent and identically distributed. Let γ be the parametric function giving their mean value.By the results above, for all valid parameter valuation v such that M v and M f have the same structure,

E l M v (r) = E(Y )(v) = E( n i=1 Y i /n)(v) = γ(v).
Our parametric approximation of the expected value is therefore:

γ = n i=1 Y i /n.
In the sequel we will this use Parametric Statistical Model Checking (PSMC) to check the formal model we will implement for the UAV.

Implementation, Experimentations and Results

While our complete formal model has been introduced in Section 2 in the form of an automata, we now explain how we successively implemented and improved the model by considering different formalisms and model checking tools. At each step, we show the limitations of the related model which leads to the next step of the implementation. The different steps of the model implementations are depicted in Figure 7.

To start, a first partial version of the formal model of Section 2.3 was implemented as a PRISM model using the PRISM tool [START_REF] Kwiatkowska | PRISM 4.0: Verification of probabilistic real-time systems[END_REF], without parameters. This first version, as depicted in 7a, corresponds to a very simple UAV flight plan, going in a straight line from point A to point B in T time units. In this context, the intermediate positions are estimated T * f times, where f is the frequency of the filter. The sizes used for the five security zones are respectively 20m, 40m, 60m, 80m and 100m. As explained in Section 2, the filter removes the noise corrupting data coming from sensors. In this first version, we only consider potential deviations along the y-axis. At each computation step, the inaccurate position given by the filter is computed using the accuracy of the filter and sensors (as a single real-valued variable), and compared to the intended position as given by the flight plan. The safety zone is deduced from the distance between the estimated position and the intended position. If the UAV enters Zones 4 or 5, the computation stops.

In this first model, the accuracy of the filter is probabilistic but not parametric, i.e. probability values have been encoded directly in the model. These values are the results of a set of experiments performed by using a flight controller plugged on a production line with a predefined path with a loop. We launched several runs of the device on the production line path and measured the outputs of the EKF filter. These measures then allowed us to compute the estimated position, which can then be compared to the exact position on the production line. We consequently obtained probabilities for the accuracy of the position estimation using an EKF filter and sensors coming from an industrial UAV. However, the major drawback of these experimentations is that they did not reflect a realistic UAV environment. In particular, since the experiment was conducted indoor using a fixed production line, the precision of some of the sensors (GPS for instance) is not representative of the precision one could obtain in a realistic flight environment. Although we were able to verify this model using PRISM, the results are not representative and can only be considered as a proof-of-concept. Since our aim is to study the same problem for different accuracy probabilities, we changed the exact probability values to parameters and submitted this new model to the PRISM Model Checker. However, because of the real-valued variables used in the model and of the numerous intermediate computations, PRISM was not able to handle this model and timed-out after 2 hours of unsuccessful computations.

Facing these shortcomings with the PRISM tool, we considered the implementation of our model with the PARAM tool [START_REF] Hahn | PARAM: A Model Checker for Parametric Markov Models[END_REF] which is a model checker for parametric discrete-time Markov chains. PARAM is efficient and allows to compute the probability of satisfying given properties as polynomials or rational functions of the parameters. As PRISM, PARAM also failed to model check our current version of the model. At this stage, since both PARAM and PRISM failed to verify our simplest model because of its complexity, we considered using a different approach based on Parametric Statistical Model Checking. For this purpose, we developed a prototype tool 5 . In this context, our model was expressed as a python program using real-valued variables both for the position of the UAV and for the probabilistic parameters. It appears that PSMC is particularly efficient in this context, and was able to verify our model (by performing more than 20k simulations) in less than 1 minute. We therefore chose to pursue our experimentation using this prototype tool and refined versions of our model.

In the second version of the model, depicted in Figure 7b, we allowed deviations to also occur along the x-axis. This is not problematic when considering a straight line flight plan, but could become important as soon as the flight plan is curved (as the one in Figure 1). Indeed, in this context, deviations along the x-axis (for example if the drone is "late") could result in the PID deciding to cut the trajectory, i.e. going straight to point C before reaching point B, therefore promoting a trajectory that might colide with the forbidden safety zones. Again, our tool managed to verify this model in a very short time.

For the third version of the model as depicted in Figure 7c, we add a third target point to the flight plan, which is not aligned with the first to points, i.e. like in Figure 1.In this third version, the inacurracy of the position estimation along the x-axis also allows the UAV to be "late" and decide to cut the flight plan as explained above.

Finally, the last version, as depicted in Figure 7d, takes into account wind perturbations. We assumed here that the wind direction was constant but that the wind force was again parametric. This will allow us to study the right tradeoff between filter capacity and frequency depending on the weather conditions. This last version is the most complex we studied, and therefore took more time to verify than the previous ones. With our prototype tool, it took 190 seconds to perform the verification using 10k simulations in this context while the same amount of simulations only took 28 seconds for the previous model (without wind parameters). 6The outputs of our prototype tool are multivariate polynomials on the parameters of our model. Given the number of parameters, the size of the model and the length of the considered simulations, these polynomials are quite complex and therefore difficult to report in this paper. As an example, below is the output polynomial representing the probability that a UAV enters Zones 4 or 5 using our last version of the model: 0.43 * P robaF ilter 3 * P robaW ind 1 + 0.16 * P robaF ilter 3 * P robaW ind 2 + 0.17 * P robaF ilter 3 * P robaW ind 3 + 0.28 * P robaF ilter 3 * P robaW ind 4 + 0.85 * P robaF ilter 4 * P robaW ind 1 + . . . Instead of showing the resulting polynomials, we will only present the evaluation of these polynomials using realistic values for the parameters. We defined two scenarios (Scenario 1, Scenario 2) with one set of values of parameters for each scenario. For these two scenarios, ProbaF0 (resp F1, F2, F3, F4) models the probability that the estimated position is from 0 to 2m (resp. 2-4m, 4-6m, 6 -8m, 8 -10m) from the real position. In the first (resp. second) scenario, we have set these values to 0.15/0.3/0.4/0.1/0.05 (resp. 0.1/0.25/0.35/0.2/0.1). According to experiments done at PIXIEL, the first scenario is more realistic than the second one. Similarly, the wind parameters correspond to the probability of having a wind force of 0 -20km/h,20 -30km/h, 30 -50km/h and 50 -70km/h respectively and have been set to 0.55/0.43/0.01/0.01 (which corresponds to typical weather conditions in Nantes, France) for the numerical evaluation. In both scenarios, Zone 4 (resp. 5) is situated 8m (resp. 50m) from the flight plan.

In Table 1, we gather the results for running the simulation for the two considered scenarios; the simulation with PSMC is performed with 10k, 20k and 50k samples. Each time, a polynomial is computed and then evaluated using the parameter values given above. In order to illustrate the stability of our results despite their statistical nature, each complete scenario was performed two times (labelled V1 and V2 in the table). The value reported in the table represents the probability of the UAV eventually reaching Zones 4 or 5 during its flight. Experiments were performed using the formal models presented in Figure 7c (without wind) and Figure 7d (including wind perturbations) on a flight plan resembling the one shown in Figure 1, with a total flight duration of 5s and a filter frequency of 1Hz. We considered two versions of the model from Figure 7d: 7d(np) where wind strength is directly input as a constant probability in the model (resulting in a polynomial where the only variables represent the precision of position estimation), and 7d(p) where wind strength is input as parameter variables in the model (allowing to evaluate/optimize the resulting polynomial according to any wind strength). Remark that the results in the first case are more precise because there are less variables in the polynomial, and obtained in a more efficient manner. Depending on whether we are interested in specific or generic information concerning the weather environment, we can chose to use the first of the second version. Remark that the probabilities of entering the forbidden zones are quite high. This is not surprising as Zone 4 is situated 8m from the intended trajectory and the precision of position estimation can be up to 10m. These values have been made deliberately high for the purpose of this study but can be chosen more realistically when verifying the real model.

Conclusion and Future Work

In this paper, we have presented a formal model to study the safety of a UAV in automatic flight following a predefined flight plan. This formal model consists in a parametric Markov Chain and takes that takes into account the precision of position and attitude estimation using sensors and filters as well as potential wind perturbations. We have also proposed a new verification technique for parametric probabilistic model: parametric Statistical Model Checking. This new technique has been implemented in a prototype tool. While state of the art tools such as PRISM and PARAM have timed out on the verification of the simplest version of our formal model, our prototype tool has been able to successfully verify the most complex version in less than 12 minutes.

In the future, we plan to keep enhancing our model in order to include filter frequency to be used as a parameter in the model. Using these parameters will allow us to obtain the parametric probability to enter dangerous zones depending on both the filter frequency and the precision probabilities. Studying/optimizing this parametric probability will allow PIXIEL to work on the trade-off between frequency and precision in order to choose their components wisely depending on their intended flight plan.

Fig. 4 :

 4 Fig. 4: A flow diagram of the formalization steps

Fig. 6 :

 6 Fig. 6: Global behaviour of the FCS

Fig. 7 :

 7 Fig. 7: Incremental development of the SMC model

Table 1 :

 1 Results of the experiments

		Model	10k			20k	50k
			V1	V2	V1	V2	V1	V2
	Running time 7c	28s		51-54s	142-143s
	Scenario 1	7c	4.99% 5.09% 4.74% 5.10% 4.91% 4.98%
	Scenario 2	7c 10.38% 10.04% 9.82% 10.05% 9.95% 9.81%
	Running time 7d(np)	28s		53-54s	149-155s
	Scenario 1	7d(np) 5.43% 5.31% 5.61% 5.21% 5.59% 5.47%
	Scenario 2	7d(np) 10.8% 10.9% 10.8% 10.8% 10.9% 10.7%
	Running time 7d(p)	185-190s	311-314s	612-621s
	Scenario 1	7d(p) 4.95% 5.97% 5.28% 6.62% 4.16% 5.61%
	Scenario 2	7d(p) 9.55% 9.87% 10.3% 11.3% 9.57% 10.7%

Particular SMC techniques also allow to estimate the satisfaction of qualitative properties[START_REF] Legay | Statistical model checking: An overview[END_REF].

available at https://github.com/paulinfournier/MCpMC

We do not share the exact models used in our prototype tool for confidentiality reasons, but the models used in PRISM and PARAM can be found here: https://github.com/br4444/modelPrism/tree/master


Thomas T Hildebrandt

Declarative Choreographies and Liveness

Keywords: Choreographies, Liveness, Declarative Models

We provide the first formal model for declarative choreographies, which is able to express general omega-regular liveness properties. We use the Dynamic Condition Response (DCR) graphs notation for both choreographies and end-points. We define end-point projection as a restriction of DCR graphs and derive the condition for end-point projectability from the causal relationships of the graph. We illustrate the results with a running example of a Buyer-Seller-Shipper protocol.

Introduction

Choreographies are an important tool for the development of highly distributed applications. Using an "Alice-talks-to-Bob" notation, they permit to abstract away details of a distributed implementation and focus on how the different components interact. This has been a fundamental driver for the adoption of choreographies in industry standards such as Message Sequence Charts (MSC) [START_REF]Message Sequence Chart (MSC)[END_REF], UML Sequence Diagrams [START_REF]Object Management Group UML Technical Committee: Unified Modeling Language[END_REF], WS-CDL, and BPMN Choreography Notation [START_REF]Object Management Group BPMN Technical Committee: Business Process Model and Notation, version 2.0[END_REF]. Moreover, choreography notations have been used in a range of application areas, including web-service development [START_REF] Carbone | Structured communication-centered programming for web services[END_REF][START_REF] Carbone | Deadlock-freedom-by-design: Multiparty asynchronous global programming[END_REF], synthesis of protocol behaviour [START_REF] Lange | From communicating machines to graphical choreographies[END_REF], monitoring [START_REF] Bocchi | Monitoring networks through multiparty session types[END_REF], parallel programming [START_REF] López | Protocol-based verification of message-passing parallel programs[END_REF] and cyber-physical systems [START_REF] López | Enforcing availability in failure-aware communicating systems[END_REF]. Paired with static analysis techniques (e.g., behavioural type systems), they are capable of deriving distributed (endpoint) implementations where endpoints generated from a choreography ascribe to all and only the behaviors defined by it. In practice, choreographies have "graduated" from the academic world, and several industrial programming languages implementing choreographies exist, e.g., [START_REF] Kouzapas | Typechecking protocols with mungo and stmungo: A session type toolchain for java[END_REF][START_REF] Carbone | Deadlock-freedom-by-design: Multiparty asynchronous global programming[END_REF][START_REF] Honda | Scribbling interactions with a formal foundation[END_REF].

A central aspect of choreography languages is the notion of interaction. An interaction is a first class citizen in any choreography language and, as a minimum, it collects information regarding the sender, the receivers, and the action used to synchronize participants. In Fig. 1, we show an exemplary BPMN choreography, based on a variant of the Buyer-Seller protocol [START_REF] Carbone | Structured communication-centered programming for web services[END_REF]. The choreography involves three participants, a Buyer, a Seller and a Shipper. After asking the Seller for a quote and getting the reply, the Buyer may either Accept, Reject or Ask again. If the Buyer accepts, the Seller sends an Order to the Shipper, which subsequently sends the detailed confirmation directly to the Buyer.

The second aspect considered in the design of choreography languages is the ordering of interactions. Typically, choreography languages are seen as imperative programs as the BPMN choreography above, that describe how interactions should occur. Any other flow not explicitly written in the language is considered forbidden. It has been observed that imperative notations are often insufficiently flexible for modelling business processes [START_REF] Van Der Aalst | DecSerFlow: Towards a truly declarative service flow language[END_REF][START_REF] Reijers | Declarative modeling-an academic dream or the future for bpm?[END_REF]. An imperative notation focuses on describing a small number of ideal flows through a process. Adding more flows to represent edge cases and less common solutions to the model tends to increase its complexity significantly. While this approach works for processes where the ideal case is all we are interested in, this does not suffice for knowledge and case work: knowledge workers tend to deal with highly variant scenarios for which they need to determine unique solutions. For instance, in the choreography above, we may in practice really want the liveness property, that a Quote is eventually followed by a decision, but that the Sellers can provide new quotes or the Buyer can ask again, any finite number of times, before accepting or rejecting a quote. The present imperative choreography languages do not allow to specify such general liveness property.

In the present paper, we propose using the declarative Dynamic Condition Response (DCR) graph notation [START_REF] Mukkamala | A Formal Model For Declarative Workflows: Dynamic Condition Response Graphs[END_REF][START_REF] Slaats | Flexible Process Notations for Cross-organizational Case Management Systems[END_REF][START_REF] Hildebrandt | Declarative event-based workflow as distributed dynamic condition response graphs[END_REF][START_REF] Debois | Safety, liveness and run-time refinement for modular process-aware information systems with dynamic sub processes[END_REF][START_REF] Mukkamala | Towards trustworthy adaptive case management with dynamic condition response graphs[END_REF] as a formal declarative notation for both choreographies and end-point specifications, allowing the specification of both safety and general liveness properties. The DCR graphs notation has been developed for the formalisation and digitalisation of collaborative, adaptive case management processes. The notation is both supported by a range of formal techniques, and serves as the formal base for the industrial (dcrgraphs.net) design and simulation tool. During the recent years, the DCR graphs technology has been employed in major industrial case management systems used in the public sector in Denmark. DCR graphs have been extended to include both data [START_REF] Strømsted | Dynamic evaluation forms using declarative modeling[END_REF], time [START_REF] Hildebrandt | Contracts for crossorganizational workflows as timed dynamic condition response graphs[END_REF] and sub-processes [START_REF] Debois | Safety, liveness and run-time refinement for modular process-aware information systems with dynamic sub processes[END_REF]. In the present paper we consider only the core notation, which is expressive enough to represent both regular and omega-regular languages [START_REF] Debois | Safety, liveness and run-time refinement for modular process-aware information systems with dynamic sub processes[END_REF] as well as so-called true concurrency [START_REF] Debois | Concurrency and asynchrony in declarative workflows[END_REF]. This means that we provide the first choreography model supporting end-point projection and general liveness properties. Definition and simulation of DCR graphs is supported by the on-line DCR Workbench [START_REF] Debois | The DCR Workbench: Declarative Choreographies for Collaborative Processes[END_REF] available at http://dcr.tools/ forte19. DCR diagrams in this paper were all produced using the workbench.

One of the important reasons for using choreography languages is their correctness-by-design guarantees. Message-passing distributed systems consist of communicating endpoints whose behaviours are defined in terms of input/output actions. So if the choreography is to be implemented by a message-passing distributed system, it is necessary to translate choreographies into code that can be executed by these endpoints. Such a translation is referred to as an endpoint projection. The endpoint projection, paired with the global properties of the choreography, warrants the safety of the distributed execution of the endpoints (e.g. deadlock-freedom). A catch of this approach is, however, that the choreography language often allows specifications that are not well-formed, meaning that it is not possible to realise the choreography as the composition of end-point processes. A key result for any choice of choreography language and end-point language is therefore to provide criteria for the choreography to be well-formed.

A core property is that of local causality. Intuitively, local causality means that if a participant initiates an interaction, it must not have direct dependencies (or causal relationships) to interactions in which this participant is not involved. The criteria for end-point projectability is, however, highly dependent of the chosen languages. In BPMN 2.0.2 it is formulated as a constraint on sequencing "The Initiator of a Choreography Activity MUST have been involved (as Initiator or Receiver) in the previous Choreography Activity." as well as a number of more complex constraints on the use of so-called branching gateways in BPMN for choices. Proving the correctness of such criteria requires a formal semantics, which is not yet provided for BPMN Choreographies, but for similar notations [START_REF] Carbone | Structured communication-centred programming for web services[END_REF]. In the present paper, we will build upon the formal semantics and theory of safe projections [START_REF] Hildebrandt | Safe distribution of declarative processes[END_REF][START_REF] Hildebrandt | Contracts for crossorganizational workflows as timed dynamic condition response graphs[END_REF] for DCR graphs to provide endpoint projections for DCR choreographies.

Summary of contributions: We provide a general end-point projection result for: 1) a declarative choreography model, 2) that can represent general omega-regular liveness properties [START_REF] Debois | Safety, liveness and run-time refinement for modular process-aware information systems with dynamic sub processes[END_REF], 3) supports a broad range of extensions such as dynamic process spanning and refinement [START_REF] Debois | Safety, liveness and run-time refinement for modular process-aware information systems with dynamic sub processes[END_REF], true-concurrency semantics [START_REF] Debois | Concurrency and asynchrony in declarative workflows[END_REF] and time [START_REF] Hildebrandt | Contracts for crossorganizational workflows as timed dynamic condition response graphs[END_REF], and 4) is supported by both academic and industrial design and simulation tools.

Interactions and Dynamic Condition Response Graphs

In this section we first define the general concept of interactions, which are common to previous work on choreographies. We then recall the model of Dynamic Condition Response (DCR) graphs.

Interactions

Assume a fixed set of actions A, ranged over by a, b, c and a fixed set of roles R, ranged over by r, r , r 1 , r 2 , ... (referred to as participants in [START_REF] Carbone | Structured communication-centred programming for web services[END_REF]). Definition 1. An interaction is a triple (a, r → R), in which the action a ∈ A is initiated by the role r and received by the roles R ⊂ fin R\{r}, i.e a finite set of roles distinct from r. Define Initiator((a, r → R)) = r. We use the shorthand (a, r → r ) for interactions between two participants (a, r → {r }). We denote by IA the set of all interactions.

We proceed to define projections of interactions to actions for end-point processes. End-point processes describe the view of the process from a single participant r synchronising with the other participants via messages on channels: For each interaction (a, r → R) in the choreography, there will be channels (a, r → r ) from r to r for each r ∈ R. To ease the definition of projections and avoid introducing new notation, we describe actions for an end-point also as interactions. That is, for the end-point process at role r, we use the interaction (a, r → R ) to represent the action !(a, r → R ) for sending a message on the channels (a, r → r ) for all participants r ∈ R . The interaction (a, r → r), represents the action ?(a, r → r) for receiving a message on the channel (a, r → r). We apologize to the reader for the inconvenience this reuse of notation may cause.

Definition 2. For an interaction α = (a, r → R ), define the end-point projection of α at r by:

α |r =    (a, r → r) when r ∈ R (a, r → R ) when r = r τ otherwise (1)
We extend end-point projections to sets and sequences of interactions by pointwise projection and removing τ actions, and finally to sets of sequences of interactions in the obvious way.

Definition 3. A choreographic language is a triple (C, A, R) where A ⊆ A, R ⊆ f in R, and C ⊆ IA ∞ = IA * ∪ IA ω .
That is, a set C of finite and infinite sequences of interactions for a given set of actions A and roles R.

When A and R are obvious from the context, we shall take C as defining a choreographic language.

Definition 4. The end-point projection of a choreographic language (C, A, R) is the family of languages (C |r ) r∈R .

DCR Graphs

In this section we recall Dynamic Condition Response (DCR) graphs [START_REF] Debois | Replication, refinement & reachability: complexity in dynamic condition-response graphs[END_REF][START_REF] Hildebrandt | Declarative event-based workflow as distributed dynamic condition response graphs[END_REF][START_REF] Mukkamala | A Formal Model For Declarative Workflows: Dynamic Condition Response Graphs[END_REF][START_REF] Debois | Concurrency and Asynchrony in Declarative Workflows[END_REF][START_REF] Debois | The DCR Workbench: Declarative Choreographies for Collaborative Processes[END_REF]. This paper follows the set-based formulation of [START_REF] Hildebrandt | Declarative event-based workflow as distributed dynamic condition response graphs[END_REF][START_REF] Mukkamala | A Formal Model For Declarative Workflows: Dynamic Condition Response Graphs[END_REF][START_REF] Debois | Concurrency and Asynchrony in Declarative Workflows[END_REF].

As formally defined below, a DCR graph consists of a directed graph and a marking. The nodes of the graph are labelled events and the edges are relations of five kinds: conditions (→•), responses (•→), inclusions (→+), exclusions (→%) and milestones (→ ).

Definition 5. A DCR graph is a tuple (E, M, L, , →•, •→, → , →+, →%), where -E is a set of events -M ⊆ E × E × E is a marking -L is a set of labels -: E → L is a labelling function -φ ⊆ E × E for φ ∈ {→•, •→, → , →+, →%} are relations between events.
A DCR graph defines a process whose executions are finite and infinite sequences of (labelled) events. Note that an event may be executed several times. The three sets of events in the marking M = (Ex, Re, In) defines the state of the DCR graph process, and are referred to as the executed events (Ex), the pending response4 events (Re) and the included events (In). The relations define effects of the execution of events and constrain the executions of the process defined by the DCR graph as defined formally below. Briefly:

-An inclusion (respectively exclusion) relation e →+ e (respectively e →% e ) means that if e is executed, then e is included (respectively excluded). -A condition relation e →• e means that e is a condition for e , i.e. if e is included, then e must have been executed for e to be enabled for execution. -A response relation e •→ e means that whenever e is executed, e becomes a pending response. During a process execution, a pending event must eventually be executed (which makes it no longer pending, unless it has a response relation to itself) or be excluded. We refer to e as a response to e. -A milestone relation e → e means that if e is included it must not be pending for e to be enabled for execution. We refer to e as a milestone for e . Milestones are typically used in cyclic behaviour, when some earlier executed event e may be required to be executed again, i.e. it becomes pending, before the process can proceed executing event e .

For DCR graph G with events E and marking M = (Ex, Re, In) and event e ∈ E we write (→• e) for the set {e ∈ E | e →• e}, write (e •→) for the set {e ∈ E | e •→ e } and similarly for (e →+), (e →%) and (→ e). We can now define when the events of a DCR graph are enabled. That is, enabled events (a) are included, (b) their included conditions have already been executed, and (c) have no pending included milestones.

Example 7. We give an example of a DCR graph in Fig. 2 as visualised by the online-tool dcr.tools/forte19. Events are indicated by boxes with solid borders and collections of events are shown with dashed boxes. Relations are shown as arrows between the boxes. As formalised in [START_REF] Hildebrandt | Nested Dynamic Condition Response Graphs[END_REF], such collections are referred to as "nestings" and are just a visual shorthand, understanding arrows to (from) nestings to represent arrows to (from) every event inside the nesting.

Traditionally, the labels of DCR graphs only consist of an action and possibly a set of roles that may perform the action. In the present paper, however, the labels of the DCR graphs are interactions rather than just actions. Instead of labelling the boxes representing the events simply with the label, e.g., (Ask, Buyer → {Seller1, Seller1}), we have split the label in three fields in the visualisation similarly to the notation for BPMN choreographies: The initiator (Buyer) is written in the field at the top of the box; the action (Ask) is written in the middle field, and the receiver(s) (Seller1, Seller2) in the bottom field with a grey background. When no confusion is possible, we refer to events by the action shown in middle field, speaking of, e.g., "the event Ask" rather than the more precise "the event labelled Buyer, Ask, Seller1,Seller2."

The marking of the graph and whether events are enabled or not is indicated visually: If the background of all fields is grey for a box, the event is included, but not enabled. For instance, the Ask event is enabled, but the two Quote events are not enabled (because the Ask event is a condition for the events and not yet executed). A box which is made opaque/dimmed out, such as the two Order events and the Details event, represents an event which is not included. When explaining Fig. 3 we describe the visualisation of executed and pending events.

The response relation (•→) from the event Ask to the nesting box around the Quote events means that the Quote events become pending when Ask is executed. This expresses the liveness constraint, that if the buyer asks, a quote must eventually be given by both sellers. The milestone relation (→ ) from the box around the Quote events to the nesting box labelled Decide means that the events (Accept1, Accept2 and Reject) inside the Decide box cannot be executed if any of the Quote events are pending, not even if Quote happened in the past. This expresses the safety condition, that the buyer can not accept or reject if one of the sellers has not responded after the last time the buyer asked for a quote. The inclusion relation from e.g. Accept1 to Order for Seller 1 means that order event will be included if the buyer accepts the quote from Seller 1. The circular exclusion arrow in the Decide box means that any event inside the box is related by an exclude relation to any event inside the box, i.e. they are mutually and self-exclusive. That is, whenever one of the events inside the box happens, all three events inside the box are excluded. Moreover, due to the exclude relation from the Decide box to the Negotiate box, also the Ask and Quote events are excluded when a decision is made.

Below we formalise how the marking changes when an enabled event e is executed : (a) the event e is added to the set of executed events, (b) e is removed from the set of pending response events, and the responses to e are added to the set of pending response events, (c) the events excluded by e are removed from the set of included events, and the events included by e are added to the set of included events. Example 9. In the graph in Fig. 2, we may execute the event Ask. Following the relations in the graph, this will make the two Quote events enabled and pending: they were previously not enabled due to their condition relation from Ask and once Ask becomes "executed" in the marking, that condition is fulfilled, and Quote becomes enabled. Altogether, executing Ask yields the graph shown in Fig. 3. Red text and an exclamation mark after the action in the middle field represents an event which is pending, as is the case for the two Quote events. A box with a check mark after the label represents an event which is executed, as can be seen for the event Ask. Note that Ask may be executed again immediately, leaving the process in the same state, but the two Quote events must eventually be executed (without any intermediate Ask) in order for the run to be accepting. From the definition of execution we can define a transition semantics for DCR graphs using labelled event transition system with responses. We say that two LETSR T and T are isomorphic, written T ≡ T , if there is an isomorphism between the sets of states preserving and respecting transitions and the response function. 5We define the language of a DCR graph as all finite and infinite sequences of such executions, where we demand that all pending responses are either eventually executed or excluded.

Definition 11 (Language of a DCR graph). Let G = (E, M, L, , →•, •→ , → , →+, →%) be a DCR graph. A run of G is a finite or infinite sequence of events e 0 , e 1 , . . . such that e i ∈ enabled(G i ), execute(G i , e i ) = G i+1 , and G 0 = G. We call a run accepting iff for each G i with marking M i = (Ex i , Re i , In i ) and e ∈ Re i ∩ In i there exists a j ≥ i such that e j = e or e ∈ Re j ∩ In j .

The language lang(G) ⊆ L ∞ of G is the set of finite and infinite sequences of labels l 0 l 1 • • • such that there is an accepting run e 0 , e 1 , . . . where (e i ) = l i .

It has been proven in [START_REF] Debois | Safety, liveness and run-time refinement for modular process-aware information systems with dynamic sub processes[END_REF] that DCR graphs can express exactly the languages that are the union of a regular and an ω-regular language. This means that one can express regular safety and liveness properties in DCR graphs.

Since the definition of accepting runs only depends on the included pending responses in the markings of the graphs and the events being executed during a run, it is easy to see that if two DCR graphs have isomorphic transition systems with responses then they also have the same languages.

DCR Choreographies

Below we first account for how DCR Choreographies and DCR End-points can be defined using DCR graphs. We then derive the criteria for end-point projectability and provide the operational correspondence between an end-point projectable DCR graph and the synchronous composition of its end-points. Definition 13. Let A ⊆ A and R ⊂ f in R be sets of roles and actions, respectively. A triple (G, A, R) is then a DCR choreography when G is a deadlock-free DCR graph such that the labels L ⊆ IA of G are interactions with actions in A and participants in R. For a role r ∈ R, a tuple (G, A, R, r) is a DCR Endpoint when the labels L of G are interactions either of the form (a, r → R ) or (a, r → r) with a ∈ A.

Example 14. The DCR graph in Fig. 2 is the DCR graph G of a DCR choreography (G, A, R) where the actions A and roles R are given by: A = {Ask, Quote, Accept1, Accept2, Reject, Order, Details} R = {Buyer, Seller, Shipper} Note that by virtue of being a general declarative notation, one may specify DCR graphs with deadlocks, e.g. by having a cycle of condition relations. It is easy to prove, that if (→• ∪ → ), i.e. the union of the condition and milestone relations, is acyclic, then the DCR graph is free of deadlocks. Moreover, such graphs can express all languages expressed by general DCR graphs, and in particular the complex behaviour in our running example.

We now turn to the key question for any choreography language: How do we project a global choreography description onto the intended behaviour of individual participants? And in particular, is this operation always possible, or are some global descriptions in fact not realisable by individual end-points?

Projections and distributed execution have been studied for DCR graphs in previous work [START_REF] Hildebrandt | Safe distribution of declarative processes[END_REF][START_REF] Hildebrandt | Declarative modelling and safe distribution of healthcare workflows[END_REF][START_REF] Hildebrandt | Contracts for crossorganizational workflows as timed dynamic condition response graphs[END_REF], but in a rather different setting where events are only labeled by actions and initiating roles, not recieving roles. For this reason, it is safe in [START_REF] Hildebrandt | Safe distribution of declarative processes[END_REF][START_REF] Hildebrandt | Declarative modelling and safe distribution of healthcare workflows[END_REF][START_REF] Hildebrandt | Contracts for crossorganizational workflows as timed dynamic condition response graphs[END_REF] to leave out an event in the projection to an end-point, if the execution of this event does not impact the state or enabledness of any event initiated by the participant responsible for that end-point. An example of such an event in our running example is the Details event for the Buyer end-point.

In the present paper we have explicit receivers and need to preserve all receiving events for a participant. Consequently, we can not directly use the notion of projection given in [START_REF] Hildebrandt | Safe distribution of declarative processes[END_REF][START_REF] Hildebrandt | Declarative modelling and safe distribution of healthcare workflows[END_REF][START_REF] Hildebrandt | Contracts for crossorganizational workflows as timed dynamic condition response graphs[END_REF]. However, we may yet build on these projections to obtain one useful for DCR choreographies. The core idea in those papers was to project a graph G with events E to a network of local graphs for any division (not necessarily disjoint) of the events δ 1 ∪ δ 2 ∪ . . . ∪ δ n = E, then define synchronous composition of such networks of DCR graphs. Intuitively, shared events, i.e. events occurring in more than one graph, are executed synchronously in the network, representing communication. The projection then ensured that execution of the network formed by the local graphs would have a transition system isomorphic to that of the global graph, and thus in particular exhibit the same language as the global graph.

In the following we reconcile the notion of projection from [START_REF] Hildebrandt | Safe distribution of declarative processes[END_REF][START_REF] Hildebrandt | Declarative modelling and safe distribution of healthcare workflows[END_REF][START_REF] Hildebrandt | Contracts for crossorganizational workflows as timed dynamic condition response graphs[END_REF] and then subsequently define end-point projections.

First, we characterise when the execution of an event may change the marking or enabledness of another. To this end, we define the notion of direct dependency. That is, e e iff either (1) they are the same, (2) there is a relation from e to e, (3) e includes or excludes an event which is itself a condition or milestone for e, or (4) e has a response to a milestone for e.

The following proposition states that an event e must be directly dependent on any event e whose execution may change the marking or enabledness of e. Proof. (Sketch) For lack of space we just show why condition 2 above implies e e, the other conditions follow from a similar inspection of the definitions. First note that the set Ex of executed events always grows, i.e. once executed an event can never become not executed. So we only need to consider the case e ∈ Ex and e ∈ Ex . From the Def. 8 it is clear that the only event that can be included in the set Ex during execution is the event being executed, so it follows that e = e and thus e e.

We note that this implication is not a bi-implication, e.g., in the DCR graph comprising just the two events e, f and the single relation e →• f in a marking where e is already executed, clearly have e f (by Definition 15(2) because there is a relation from e to f ), yet executing e in fact cause no changes to marking or enabledness of f .

Intuitively, we will obtain the end-point projection for a participant r by keeping (a) events labelled with interactions involving r, as well as (b) the direct dependencies of the events for which r is the initiator. We then interpret the interactions as end-point actions as described in Sec. 2.1. In order for the interactions to make sense as actions for the end-point process at r, the role r must be involved in its direct dependencies. We formalise this as follows.

Definition 17. Let (G, A, R) be a DCR choreography and the labelling function of G; and let r ∈ R be a role. This choreography is end-point projectable for r iff for all e, if Initiator(e) = r and e e, then (e ) |r = τ , Example 18. Referring again to the example DCR choreography in Fig. 2, we find that this choreography is in fact not end-point projectable for the participants Seller1 and Seller2. We see in Fig. 2 that the Accept1 event causes Accept2 to be excluded, however, Seller2 is initiator of Accept2, but not participating in Accept1. To be precise, because of the exclusion we have Accept1 Accept2 and (Accept2) = (Accept2, Buyer → Seller2) yet (Accept1) = (Accept1, Buyer → Seller1), so (Accept1) |Shipper = τ .

We fix this by redefining the choreography such that Seller2 is included in the Accept1 interaction, that is, so that Seller2 is notified that he lost the contract; and vice versa including Seller1 in the Accept2 interaction. We show the projectable process in Fig. 4. This choreography is end-point projectable.

We proceed by defining the end-point projection. We start by recalling the definition of projection in [START_REF] Hildebrandt | Safe distribution of declarative processes[END_REF][START_REF] Hildebrandt | Declarative modelling and safe distribution of healthcare workflows[END_REF][START_REF] Hildebrandt | Contracts for crossorganizational workflows as timed dynamic condition response graphs[END_REF], adapted to keep all labels.

Definition 19 (Adapted DCR δ-Projection cf. [START_REF] Hildebrandt | Safe distribution of declarative processes[END_REF][START_REF] Hildebrandt | Declarative modelling and safe distribution of healthcare workflows[END_REF]). Given a DCR graph

G = (E, M, L, , →•, •→, → , →+, →%)
and a set of events δ ⊆ E, define the projection of G to the events δ as the graph

G| δ = (E| δ , M | δ , L| δ , | δ , →•| δ , •→| δ , → | δ , →+| δ , →%| δ )
given by:

= Ex ∩ E| δ (b) Re| δ = Re ∩ E| δ (c) In| δ = In ∩ ((→• δ) ∪ (→ δ) ∪ δ) ∪ E| δ \ ((→•δ) ∪ (→ δ) ∪ δ) . 3. | δ (e) = (e), 4. L| δ = img( ) 5. →•| δ = →• ∩ (→• δ) × δ 6. → | δ = → ∩ (→ δ) × δ 7. •→| δ = •→ ∩ ((•→→ δ) × (→ δ)) ∪ ((•→ δ) × δ) 8. →+| δ =→+ ∩(((→+→• δ)×(→• δ))∪((→+→ δ)×(→ δ))∪((→+ δ)×δ)) 9. →%| δ =→% ∩(((→%→• δ)×(→• δ))∪((→%→ δ)×(→ δ))∪((→% δ)×δ))
The complexity in these rules arises mostly from the necessity of including events that may affect milestones or conditions for the events in δ. We see this particularly in the right-most half of 2(c), in the second clause in 7-9, and in the fourth clause in 8-9.

We now define the end-point projection for a DCR choreography with respect to a role r. The projection comes in two steps: first we compute the δ-projection, taking δ to be the set of events for which r is the initiator. Second, we simply add all events where r is a receiver. The latter step does not really change the behaviour in terms of sequences of actions, but it ensures that all receiving roles will be present for an interaction, even when it has no effect on other events in the end-point. As also described in the beginning of the section, this was not essential for the previous work, since receivers were not explicit. The end-point projection of (G, A, R) for r is then defined as the DCR end-point (G |r , A |r , R, r) where

G |r = (E |δ ∪ E , M ∪ M , img( ), , →•, •→, → , →+, →%).
Example 21. The result of end-point projecting the corrected choreography in Fig. 4 can be seen in Fig. 5, Fig. 6 and Fig. 7.

Lemma 22. Let (G, A, R) be a DCR choreography, let r ∈ R be a role of R, and let (G |r , A |r , R, r) be the projection of that choreography to r. If (G, A, R) is projectable for r, then every label in G |r is an interaction which has r as a participant.

Proof. The set of events of G |r consists of the events E |δ of the δ-projection for δ = {e ∈ E 0 | Initiator(e) = r} and the events E = {e ∈ E | ∃a r . (e) |r = (a, r , r)}. Clearly, the events in E by definition all have the role r among the receivers and thus as participant. According to Def. 19, we have E| δ = {e ∈ E | ∃e ∈ δ. e e }. Now, since Initiator(e ) = r for all e ∈ δ it follows from the definition of end-point projectability in Def. 17 that (e) |r = τ when e e and thus r is also a participant in the interaction for all e ∈ E| δ .

As shown below, it follows easily, that if an event is shared between two end-points, it has the same initiator. Proof. According to Def. 20, the label of an event e in an end-point projection for role r is given by the restriction (e) |r . According to Def. 2 the initiator role is preserved by the restriction or the label is τ . However, by Def. 17 the label cannot be τ . Thus, being end-point projections of the same end-point projectable choreography C the event in the two end-point projections will have the same initiator role.

We now define the synchronous composition of a finite set of DCR end-points, for which the labels of shared events agree on the Initiator role. Intuitively, an event e is enabled in the synchronous composition, if it is enabled in all of the end-points in which it occurs. The execution of an event is then defined simply by executing the event in all of the components it occurs. Finally, the label is the interaction obtained by taking the union of receivers.

Definition 24 (Synchronous composition of DCR end-points). For R = {r 1 , r 2 , . . . , r n } and DCR end-points P i = (G i , A i , R, r i ) for i ∈ {1, .., n} we write the synchronous parallel composition as P = Π i∈{1,..,n} P i . Define -E = i ∈ {1, .., n}E i , where E i is the events of G i .

e ∈ enabled(P ) iff e ∈ E i implies e ∈ enabled(G i ) for all i ∈ {1, .., n}.

execute(P, e) = Π i∈{1,..,n} P i , if e ∈ enabled(P ) and

P i = (G i , A i , R, r i ) and G i = execute(G i , e), if e ∈ E i and P i = P i otherwise. -P (e) = (a, r → R ) if e ∈ E i implies i (e) = (a, r → R i ) and R = i ∈ I, where I = {i ∈ {1, . . . , n} | e ∈ E i }. -ρ P (P ) = i∈{1,...,n} Re i ∩ In i if P = Π i∈{1,..,n} P i , P i = (G i , A i , R, r i ) and the marking of G i is (Ex i , Re i , In i ).
We now define the LETSR for P by T (P ) = (P, P, E, L, P , → P , ρ P ), where (P , e, P ) ∈→ P if e ∈ enabled(P ) and P = execute(P , e), and P = {P | P → * P }.

The following theorem establishes the key property, that the synchronous composition of the end-points yields a transition system with responses isomorphic to the transition system for the choreography, and thus is in particular deadlock free.

Theorem 25. Let C = (G, A, R) be an end-point projectable DCR choreography, R = {r 1 , . . . , r n } and P i = (G i , A i , R, r i ) for i ∈ {1, . . . , n} the DCR endpoints resulting from end-point projection of C. Then T (C) ≡ T (Π i∈{1,..,n} P i ) and thus Π i∈{1,..,n} P i is deadlock free.

Proof. (Sketch) The proof follows the same approach as the proof of Thm.5.1 in [START_REF] Hildebrandt | Contracts for crossorganizational workflows as timed dynamic condition response graphs[END_REF] where a bisimulation is constructed between the original graph (in this case the choreography) and the network of synchronous parallel composition of projections. The reason why the same approach can be used is that the main difference between the present work and the work in [START_REF] Hildebrandt | Contracts for crossorganizational workflows as timed dynamic condition response graphs[END_REF] is that we have included also receiving events in the end-points that have no effect on the synchronous product, such as the Details event in Fig. 5.

We note that the isomorphism by Prop. 12 implies that the language of the choreography is the same as the language of the composition of the end-points.

Conclusion and Related Work

Based on the formal process notation of DCR graphs, we have provided the first declarative model for choreographies able to describe general liveness properties. We identified the local causality criteria for end-point projectability and defined end-point projections, using previous work on distributions of DCR graphs. We showed that the synchronous product of the end-point projections had the same behaviour as the original choreography. As future work we intend to extend the results to declarative timed choreographies, benefiting from projections already being defined for timed DCR graphs in [START_REF] Hildebrandt | Contracts for crossorganizational workflows as timed dynamic condition response graphs[END_REF]. Related Work. Properties for guaranteeing projectability are proposed in various settings and depend on the chosen choreography language. The results in [START_REF] Carbone | Structured communication-centered programming for web services[END_REF] require three main properties: connectedness, well-threadedness, and coherence. While well-threadedness and coherence concern the behaviour of replicated servers, connectedness is the same as the projectability criterium of BPMN 2.0.2, which we also adopt in the present paper. The connectedness property occurs also in other works on choreography, e.g., [START_REF] Lanese | Bridging the gap between interaction-and process-oriented choreographies[END_REF][START_REF] Cruz-Filipe | Communications in choreographies, revisited[END_REF]. In the theory of multiparty session types [START_REF] Honda | Multiparty asynchronous session types[END_REF] and in Chor [START_REF] Carbone | Deadlock-freedom-by-design: Multiparty asynchronous global programming[END_REF], such property is omitted at a price of a more flexible interpretation of sequencing.

To the authors' knowledge, this is the first work considering general liveness properties at the choreography-level. Other works in the literature have studied liveness for multiparty interactions from session types and contract development: Padovani et al. [START_REF] Padovani | Typing liveness in multiparty communicating systems[END_REF] propose a type system for session types to control liveness properties. However, the model considered is roleless since types describe interactions but without specifying which roles implement them. The work in [START_REF] Debois | Type-checking liveness for collaborative processes with bounded and unbounded recursion[END_REF] extends binary session types to specify response properties, that is applied to a variant of a collaborative BPMN process language to verify whether liveness for dead-lock free processes can be achieved. A recent paper by Lange et al. [START_REF] Lange | Fencing off go: liveness and safety for channel-based programming[END_REF] investigates a bounded liveness property for GO programs, where protocols are specified as global types. Such property resembles a progress property and is not as general as our liveness. For instance, requiring that after a Quote we have eventually an Accept or a Reject cannot be expressed as bounded liveness. Honesty is a variant of liveness used in contract-oriented programming [START_REF] Bartoletti | Honesty by typing[END_REF]. In short, an endpoint is honest if it abides the sequence of actions stipulated in its contract. Honesty will fail if the contract promises the execution of an action and the endpoint does not execute it. Contracts in this sense correspond to DCR responses. We differ from [START_REF] Bartoletti | Honesty by typing[END_REF] in the sense that we do not require a session-type to verify liveness. It is specified in the model as a behavioural constraint only in the places that is required.

Our previous work in [START_REF] Carbone | A logic for choreographies[END_REF], presented a proof system for choreographies where properties such as liveness and connectedness can be expressed in terms of modal (may/must) operators. Apart from the difference on the languages explored (the global calculus in [START_REF] Carbone | A logic for choreographies[END_REF] and DCR in the present work), we differ from being able to express one-to-many communications. For DCR graphs, projections were studied in a different context in [START_REF] Hildebrandt | Safe distribution of declarative processes[END_REF][START_REF] Hildebrandt | Declarative modelling and safe distribution of healthcare workflows[END_REF][START_REF] Hildebrandt | Contracts for crossorganizational workflows as timed dynamic condition response graphs[END_REF], where all participants were implicit receivers of actions and projections thus always defined.

Acknowledgments:

  Work supported by the Innovation Fund Denmark project EcoKnow (7050-00034A); the Danish Council for Independent Research project Hybrid Business Process Management Technologies (DFF-6111-00337), and the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement BehAPI No.778233.

Fig. 1 .

 1 Fig. 1. BPMN Choreography for Buyer-Seller-Shipper example.

Definition 6 (Fig. 2 .

 62 Fig. 2. Example DCR choreography

Definition 8 (

 8 Execution). Let G = (E, M, L, , →•, •→, → , →+, →%) be a DCR graph, with marking M = (Ex, Re, In). When e ∈ enabled(G), the result of executing e, written execute(G, e) is a new DCR graph G with the same events, labels, labelling function and relations, but a new marking M = (Ex , Re , In ), where (a) Ex = Ex∪{e} (b) Re = (Re\{e})∪(e •→), and (c) In = (In\(e →%))∪ (e →+)

Fig. 3 .

 3 Fig. 3. DCR choreography after execution of Ask.

Definition 10 (

 10 Transition semantics). Let G = (E, M, L, , →•, •→, → , →+ , →%) be a DCR graph. The Labelled Event Transition System with Responses (LETSR) for G is defined as T (G) = (G, G, E, L, , →, ρ), where the DCR graph G is the initial state, E is the set of events, L is the set of labels, is the labelling function, →⊆ G × E × G is the transition relation, defined by (G, e, G ) ∈→ iff e ∈ enabled(G) and G = execute(G, e), and G = {G | G → * G }, the set of states, is the graphs reachable from the initial graph G by execution of events, and finally ρ is the response function defined on DCR graphs by ρ(G ) = Re ∩ In, if (Ex, Re, In) is the marking of G .

Proposition 12 .

 12 Let G and G be DCR graphs. If T (G) ≡ T (G ) then lang(G) = lang(G ).

Definition 15 .

 15 Let G = (E, M, L, , →•, •→, → , →+, →%) be a DCR graph and let e, e ∈ E be events of E. Then there is a direct dependency e e from e to e iff either of the following conditions are true 1. e = e, 2. e (→• ∪ •→ ∪ →+ ∪ →% ∪ → )e, 3. ∃e . e (→+ ∪ →%)e (→• ∪ → )e, 4. ∃e . e •→ e → e.

Proposition 16 .

 16 Let G be a DCR graph with marking M = (Ex, Re, In). Suppose e ∈ enabled(G), and let G = execute(G, e ) and M = (Ex , Re , In ) be the marking of G . If either of the following hold, then e e. 1. e ∈ enabled(G) ⇔ e ∈ enabled(G ), 2. e ∈ Ex ⇔ e ∈ Ex , 3. e ∈ Re ⇔ e ∈ Re , 4. e ∈ In ⇔ e ∈ In .

Fig. 4 .

 4 Fig. 4. End-point projectable DCR choreography

Definition 20 (

 20 DCR end-point projection). Let (G, A, R) be a DCR choreography with events E and labelling function . For any r ∈ R, define δ = {e ∈ E | Initiator(e) = r} , and let G |δ = (E |δ , M, L, |δ , →•, •→, → , →+, →%) be the δ-projection of G for r. Suppose M = (Ex, Re, In) and define E = {e ∈ E | ∃a r . (e) |r = (a, r , r)} M = (∅, ∅, E \(E\In)) (e) = (e) |r if e ∈ E |δ ∪ E undefined otherwise

Lemma 23 .

 23 Let C = (G, A, R) be an end-point projectable DCR choreography, r, r ∈ R be roles of R, and (G |r , A |r , R, r) and (G |r , A |r , R, r ) be the projections of C to r and r . If e ∈ E ∩ E , where E and E are the events of G |r and G |r respectively, then Initiator( (e)) = Initiator( (e )).

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. End-point projection of Fig. 4 for Buyer.

We often simply say "pending" instead of "pending response".

Isomorphism could be defined more generally by also having an isomorphism on the set of events, but the given definition is sufficient for the present paper.