A WEAKLY COUPLED MODEL OF DIFFERENTIAL EQUATIONS FOR THIEF TRACKING

. In this work we introduce a novel model for the tracking of a thief moving through a road network. The modeling equations are given by a strongly coupled system of scalar conservation laws for the road traﬃc and ordinary diﬀerential equations for the thief evolution. A crucial point is the characterization at intersections, where the thief has to take a routing decision depending on the available local information. We develop a numerical approach to solve the thief tracking problem by combining a time-dependent shortest path algorithm with the numerical solution of the traﬃc ﬂow equa- tions. Various computational experiments are presented to describe diﬀerent behavior patterns.


Introduction.
In the past few years, time-dependent models for transportation problems have been of huge interest for many researchers of various disciplines and application areas. For instance, these models are most commonly used to analyze gas or water flow in connected pipe systems [2,8], evacuation planning [7,16], manufacturing systems [9], or vehicular traffic on road networks [4,19].
The focus of the presented work lies on continuous traffic flow network models that are coupled with discrete decisions at the junctions. In particular, we deal with a thief, who has committed a crime and tries to escape. The time a thief needs to pass a certain street depends on how crowded the street is. Thus, the street density can be interpreted as a measure (or weight) in the network. The thief starts his journey at the initial time when the crime was committed and at the location where the crime took place. From there on, he follows a route which is updated from time to time depending on the current traffic density at his current position.
The simulation of traffic on roads is made using the Lighthill-Whitham-Richards (LWR) model [22,25]. For the simulation of traffic on a network, a model by Coclite, Garavello, and Piccoli [4] is used, which extends the LWR model to networks by incorporating conditions for the coupling of roads. Equipped with the evolution of density in the network, the path of the thief can be modeled. This is a different focus compared to other already existing approaches for crime modeling. Short et al. [27] present an agent-based model for the detection of crime in an urban region. Due to observations it is known that clusters or hotspots of crime evolve in space and time on a two-dimensional domain. Consequently, a hydrodynamic limit of the 448 SIMONE GÖTTLICH AND CAMILL HARTER discrete model leads to a parabolic system of partial differential equations. Theoretical investigations include local existence and uniqueness results of solutions for the PDE system [26] and a bifurcation analysis [28] as well. Recent extensions of the model consider the dynamic deployment of police to changing crime patterns [24,30] or crime probabilities with spatial data information to predict residential burglaries [29]. In particular, in [24] variations of Dijkstras shortest path algorithm are used to calculate the distance between crime locations for a given road network.
However, our intention is to model the probable temporal evolution of a single thief through a road network once the crime was committed. The dynamic of the thief is hereby influenced by the load of streets determined by (1). That means for simulation purposes, that the equations for the traffic and the thief can be solved independently in the sense that the solution of the traffic model is solved first and then used as input to solve the thief tracking model. This type of interplay can be seen as a weakly coupled model, where information flows in one direction only (in contrast to strongly coupled models where solutions affect each other). Weakly coupled models are a preliminary step to a more realistic models, where real data (e.g. from statistics) can be included. Their application is therefore limited and they mainly serve as motivation to study the coupling of already available models. However, in a strongly coupled model the thief would influence the solution of the traffic model. From an application point of view, one may think of carjacking, manipulated traffic flow systems or disregarding traffic rules.
Psychological components as described in [1] are also omitted so far. Similar modeling ideas have been already discussed in another context, see [3,5,6,11,21]. For instance, Delle Monache and Goatin [11] or Lattanzio, Maurizio and Piccoli [21] presented a model, in which a single slow vehicle influences the surrounding traffic. The regular vehicles queue up behind it, since overtaking is not possible without loss of time. The authors use a strongly coupled PDE-ODE model, which is not the same in the present case, because the solution of the LWR model is not affected by the thief. In fact, an ODE is required to capture the dependence of the thief's velocity on the surrounding traffic density in our approach.
Extending the evolution of the thief's path to networks, a criterion for the behavior of the thief at junctions has to be defined. When the thief arrives at the end of a road, he has to decide about which of the succeeding roads to pass next. This decision can be established in numerous ways. For instance, the thief can observe the traffic on the outgoing roads and either choose that way which has the lowest density and therefore promises the highest velocity to reach his target or take the shortest distance to the closest hideout into account. Perturbations of these decisions are also investigated.
The essential aim of the paper is to present a weakly coupled model for thief tracking and its numerical simulation. For the numerical study, the model assumptions are as simple as possible but not realistic. Since the focus is on numerical simulations, there is no well-posedness result.
The article is organized as follows: Section 2 is devoted to the thief tracking model. First, the analysis for a single road strip is made. The model is introduced theoretically and solutions are given for typical Riemann problems. Furthermore, the theory is extended to the use on networks. Therefore, coupling conditions for the traffic as well as for the thief are introduced. Section 3 adds a numerical perspective on the problem. The computational model is set up, analyzed numerically, and tested by some meaningful examples. Moreover, an excursus on time-dependent shortest path (TDSP) problems borrowed from graph theory is given. We present different strategies to choose a new path (or road) at intersections. These decisions rely on the thief's individual behavior and depend on several criteria. Since the dynamic of the thief is influenced by the dynamics of the LWR network model and not vice versa, we use the information on road densities to fully describe the evolution of the thief on roads and junctions as well. Once the network data is fixed, we explain how to compute the thief's time-dependent shortest path to a hideout. Section 4 contains several computational experiments that assess the characteristics of the decision criteria in a network.

2.
A coupled model for a thief on a road strip. In this section, a model which enables us to track the path of a thief moving within a road network is designed. The main traffic on the road behaves according to the LWR model. For this, a solution ρ(t, x) of the LWR model (1) with appropriate initial conditions ρ 0 (x) is assumed to be given: where ρ is the traffic density and f = ρv(ρ) is the flow function. The relation v(ρ) = v max (1 − ρ ρmax ) is the speed of the main traffic depending on ρ, where ρ max is the jam density and v max is the free speed.
The approach for tracking a person moving along a single road strip is then connected to (1). Therefore, the LWR model is coupled to an ordinary differential equation describing the position of the thief. Moreover, the theoretical background of the model is addressed. In particular, theoretical results of the ODE are analyzed for cases in which the conservation law has Riemann initial data.
The velocity of the thief is described by a smooth, decreasing function ω(ρ(t, x)). Let further y : R + → R be a function representing the thief's position on the piece of road at time t. It is commonly known that the derivation of the position with respect to timeẏ(t) represents the velocity. Thus, we can establish the Cauchy problem ẏ(t) = ω(ρ(t, y(t))), t ∈ R + , y(t 0 ) = y 0 , for an appropriate initial position y 0 . The density ρ(t, y(t)) reveals at the point x = y(t), where the thief is located at time t.
2.1. Theoretical background. The main purpose of the next paragraph is to study theoretical properties of the PDE-ODE system (1)- (2). The standard results from the theory of ODEs do not apply here, since existence and uniqueness results require a continuous right-hand side of the ODEẏ = f (t, y). As ω depends on the solution ρ of the conservation law, which exhibits shock waves and therefore discontinuities, w can neither be assumed to be continuous nor Lipschitz continuous. Thus, the theorems of Peano and Picard-Lindelöf are not applicable. The wellposedness of the problem above was analyzed by Colombo and Marson [5,6] from a theoretical point of view. They derived a result by showing forward uniqueness of Filippov solutions and their Hölder continuous dependence on the initial data of the ODE. Bretti and Piccoli [3] further provide a convergence result for wave front tracking approximate solutions. The well-posedness of ODEsẏ(t) = f (t, y) with discontinuous right-hand side in general is analyzed by Filippov [14].
The theoretical embedding of the considered problem (1)-(2) comes with Theorem 2.2 in [6]. According to this theorem, the thief's velocity function must satisfy for an α ∈ (0, 1). The violation of this condition results either in the failure of the well-posedness or in the thief traveling backwards. For the purposes of this work, it is reasonable to assume that the thief has a speed function which is similar or equal to the main traffic speed in the LWR model v(ρ) = v max (1 − ρ ρmax ). It is not very likely that an escaping thief behaves like a regular road user, as he is probably faster being on the run. Thus, his maximum velocity is increased by a parameter k ≥ 1 describing the factor by which the thief can move faster than the average vehicular traffic: As shown above, (4) satisfies (3). Thus, the problem (1)-(2) is well-posed for such velocity functions. In general, the velocity function can be adjusted in many different ways without violating the well-posedness of the problem, e.g. by a fixed velocity difference k between the thief and the average vehicular traffic. The actual choice of the thief's velocity function does not have an increased relevance for our purposes. Consequently, his velocity is set equal to the main traffic velocity from here on by setting k = 1 in (4) for simplicity reasons.
2.2. Riemann problems. The solution of the LWR model with Riemann initial data is well known, e.g. [4]. In the following paragraph, solutions of Riemann problems for the coupled PDE-ODE model (1)-(2) are investigated. Consider standard Riemann initial data Then, two cases are considered for the LWR model (1) are known. If ρ l ≤ ρ r , a shock wave results and if ρ l > ρ r , a rarefaction wave results. Since we intend to analyze the interaction of the thief's trajectory with shock and rarefaction waves, we assume him to be located at point x = −c left to the discontinuity point x = 0, such that he approaches the wave after some time. Starting at t = 0, the following system is obtained: Throughout this work, we assign ρ max = v max = 1 in order to keep the analysis clear. As a result, the velocity function of the thief looks as follows: ω(ρ) = v(ρ) = 1 − ρ and the flow function is f = ρv(ρ) = ρ(1 − ρ).
The solution approach for this setting consists of the piecewise solution of the ODE by considering different initial value problems (IVPs) for the different levels of ρ, which depend on x and t. The same cases as for the LWR model are considered: The thief departs at x = −c, where ρ(0, −c) = ρ l at the initial time. As ρ is constant on the left of the discontinuity, the ODEẏ(t) = (1−ρ(t, y(t))) is independent of t and can be solved by integration: The solution on the right-hand side of the discontinuity is derived analogously. However, the integration constant cannot be determined immediately, since the time when the discontinuity is hit and the corresponding location have to be determined first in order to obtain an initial condition. Therefore, the location of the discontinuity f (ρr)−f (ρ l ) ρr−ρ l t and the location of the thief (1 − ρ l )t − c are set equal and solved for t * .
Thus, the complete solution to the Riemann problem (6) under the condition ρ l < ρ r looks as follows: where Figure 1 illustrates how the thief's trajectory is affected by a shock wave.    Case 2. (ρ l > ρ r ) This case is slightly more complex due to the appearance of the rarefaction wave. There is not a single discontinuity with constant velocity to both sides of the discontinuity, but there is a transition interval, in which the density decreases linearly from ρ l to ρ r . This interval gets larger over time and the whole wave moves as well. Thus, there are three different pieces an IVP has to be solved for, one of them with non-constant ρ. As in case 1, the IVP {ẏ(t) = (1 − ρ(t, y(t))), y(0) = −c} for the leftmost part, where ρ(t, x) = ρ l , can be solved via integration.
Time and location when the left part of the rarefaction wave is reached is determined analogously to case 1 in order to obtain the initial condition for the second IVP. The location of the wave's left end f (ρ l )t is compared with the location of the thief and solved for t * : The corresponding location is then The second IVP is a linear, inhomogeneous IVP. The density ρ(t, x) is not constant in the center of the wave between f (ρ l )t and f (ρ r )t: The resulting IVP has the following solution: The third IVP is analogous to the first one and so is its solution: where c 3 = −c ρ l ρr . Finally, the complete solution to (6) can be formulated: where t * = c ρ l and t * * = c ρ l ρ 2 r . Figure 2 illustrates how the thief's trajectory is affected by a rarefaction wave. Note that figures 1 and 2 will be used for numerical comparisons in Section 4.1.
The upcoming section extends the theory and the results of this section to road networks. Therefore, the situation at junctions is considered in particular.
3. Network extension. Next, the model is applied to road networks by adapting the coupled PDE-ODE model to the theory of the LWR network model [4] and appropriate decision rules at nodes for the thief.
For the simulation of a thief's path in a road network, we first have to adapt the traffic equations (1) for the use in networks. Compared to the analysis on single edges, the traffic behavior at junctions has to be considered additionally. This issue has been studied widely in recent literature. In this paper, the approach of Coclite, Garavello, and Piccoli [4,15] is applied. The road network is modeled by a directed graph G = (V, E), where V describes the set of vertices representing junctions and E describes the set of edges representing roads. The length of road i is denoted by l i ∈ R + . The sets δ + n and δ − n describe the set of incoming and outgoing edges, respectively, for a vertex n ∈ V . Time is captured by the variable t > 0 and the position on a certain road i is described by x ∈ [0, l i ]. The traffic on the single edges is simulated as before using the LWR conservation law (1). In contrast to the single edge analysis, edges are bounded by their adjacent junctions at x = 0 and x = l i . For this, boundary conditions have to be imposed for solving the PDE. These are derived implicitly by the use of coupling conditions.
Consider an arbitrary but fixed vertex n ∈ V . The densities of the incoming road i and the outgoing road j are denoted byρ i (t) := ρ(l i , t) andρ j (t) := ρ(0, t). An unambiguous requirement to the system at each junction is the conservation of mass, i.e. the total incoming flow has to be equal to the total outgoing flow. Otherwise, mass could appear or disappear at a junction. The conservation of flow is described by the Rankine-Hugoniot relation For the coupling flows, the notationγ i (t) := f (ρ i (t)) andγ j (t) := f (ρ j (t)) is used. The Rankine-Hugoniot relation does not sufficiently describe the behavior at a junction. For instance, a coupling flow of zero would always be feasible. As a consequence, further coupling conditions have to be introduced for reasons of wellposedness. On the one hand, flow distribution parameters 0 ≤ α ji ≤ 1 are used to describe the distribution of the traffic at a junction. As α ji denotes the percentage of flow that move from the incoming road i to the outgoing road j, j∈δ − n α ji = 1 andγ can be stated. On the other hand, an entropy condition is imposed in order to maximize the flow at junctions. For each junction n ∈ V , the boundary densities (ρ i ,ρ j ) for the connected roads i ∈ δ + n and j ∈ δ − n are derived from the local fluxes (γ i ,γ j ), which satisfy the coupling conditions (11)-(13), see [4].
The admissible boundary flows corresponding to these constraints lie within the following intervals: andγ In a nutshell, the behavior of the traffic at junctions according to [4] can completely be described by formulating the following linear maximization problem: This linear program contains the maximization of flow (13) as objective function, it contains the conservation of flow and the traffic distribution rules (12) as a linear constraint (12), and it makes the flow satisfy the admissible ranges (14)- (15). A solution of (16) {γ i , i ∈ δ + n ;γ j , j ∈ δ − n } allows for determining the right boundary densities (ρ i ,ρ j ) for (1).

3.1.
The thief tracking network model. In order to simulate the path of a thief escaping from a certain crossing in a road network, the model presented above has to be combined with a thief behavior framework. The framework combines the tracking of the thief's trajectory on the edges with his decisions on which outgoing edge to choose at the vertices. The thief starts at an initial vertex representing the crossing where he committed the crime. Then he chooses one of the possible outgoing edges, passes that edge, and faces another edge choice at the next vertex. He decides again which edge to pass, and so on. His journey ends when the time exceeds the predefined time horizon, or when he arrives at one of the predefined destination vertices D ⊂ V , see Definition 3.1. The required edge passing times are determined using solutions of the thief tracking ODE (2).
Before introducing the thief tracking model for networks, some notations have to be clarified in order to provide a clear framework. As before, we consider a network (V, E). The edge arrival time function a e (t) is the time when a thief arrives at the end of an edge e ∈ E, if he started at time t. The edge arrival time can alternatively be referred to as a ij (t), where i = α(e) is the start vertex of edge e and j = β(e) is its destination vertex. The edge arrival time function can be defined in arbitrarily many ways. Here, the thief tracking ODE (2) is used for determining the travel time on the edges. The path of a thief can be described in two equivalent ways: by a sequence of k + 1 vertices V = (n 0 , . . . , n k ) or by a sequence of k edges E = (e 1 , ..., e k ). Thereby, the sequence A = (a 0 , . . . , a k ) contains the arrival times. a i describes the arrival time at n i or at the end of e i . Definition 3.1. Consider a network (V, E) with a set of hideouts D ⊂ V . Let T max (time horizon) be the maximum allowed escape time. A thief is said to be escaped if the following holds for the path V = (n 1 , ..., n k ), n 1 , ..., n k ∈ V he has completed: (i) The last vertex on the path is a hideout vertex: Otherwise, the thief is said to be caught.
Considering a thief progressing through a network, the path and the corresponding arrival times have to be updated each time he crosses a vertex. For a vertex n i , a i can be written as a composition of multiple edge arrival time functions: The determination of the edge arrival times a e (t) of a thief traveling through a network is explained in the following. The thief tracking ODE (2) has to be adjusted such that it yields the desired output. When the thief leaves a certain vertex n = α(e) for a chosen edge e ∈ δ − n with destination β(e) at time t 0 , the problem of determining the arrival time a e (t 0 ) is as follows: This ODE is solved for values of y within the bounds of the edge, i.e. y ∈ [0, l e ]. The solution y(t) describes the position of the thief on edge e at time t assuming that he started at the initial point x = 0 of the edge at time t 0 . For the arrival time at the end of the edge a e (t) it holds that the thief's position y(a e (t)) corresponds to the length of the edge l e , i.e. y(a e (t)) = l e (19) must be true for the arrival time. The right-hand side of the ODE ω(ρ e (t, y(t))) representing the thief's velocity is non-negative, the function representing his position y(t) is weakly increasing. Moreover, y(t) is a continuous function, which is intuitively clear thinking of a road user moving through a street. The existence of an arrival time a e (t) is guaranteed if the right hand side of the ODE is such that the thief arrives at the end of the road before the time horizon T is exceeded, i.e. the following must hold for ω: If an infinite time horizon is considered, the upper bound of the integral must be set to infinity. Due to the non-strict increase, there might exist a range of times for which y(t) = l e holds. This contradicts the uniqueness of the arrival time. However, there is only one relevant arrival time, which is expressed as follows: After clarifying the trajectory of the thief on edges in a network, we will have a closer look on the thief's behavior at junctions. Therefore, the criteria applied by the thief for selecting one of the outgoing edges at a junction are introduced.

3.2.
Route decision at junctions. Arriving at the end of an edge, the thief has to decide which of the successor edges to take next. This decision process is executed instantaneously when he arrives at the vertex and is assumed not to require time. The criteria for choosing an appropriate edge are the crucial factor for the thief's escaping performance. Five approaches for the decision at a vertex n ∈ V coming from an incoming edge e ∈ δ + n are provided in this work. It is assumed that the thief never returns to the vertex he just came from. Thus, if there is an edge leading back to the previous vertex, the thief is not allowed to take it. However, if this is the only outgoing edge from a vertex, the thief can pass that edge to prevent himself from being stuck at the vertex. Moreover, inflow and outflow edges can never be passed by the thief since the thief's hideouts are nodes that can be anywhere in the system, see Definition 3.1. The set of outgoing edges from a vertex n ∈ V being accessible for a thief coming from an edge e ∈ E with β(e) = n is denoted bỹ Given the arrival time at the current vertex t a , choose that edge k ∈δ − e,n , which has minimal c DC2 k (t a ). (DC3) In addition to c DC2 k (t), the distance from the thief's current location to the closest destination via edge k is taken into account. Let s k ij (t) be the length of the shortest path with respect to length from a vertex i to a vertex j postulating that edge k ∈ δ − i is chosen first and departure time at i is t. The following definition describes the distance from n to the closest destination d ∈ D supposing that edge k is chosen first: c DC2 Note that solutions to (DC2) and (DC3) are unique due to the fact that the edge with the lowest index is chosen in the case of multiple edges with the same ranking.
(DC2) is motivated by the aim for high velocity. The velocity function of the thief is weakly decreasing. Hence, the smaller the density, the faster the thief. By calculating the integral over the edge, the total density on the edge is captured.
The result describes both the average density and the average velocity on the edge. However, the edge which is chosen by (DC2) is not necessarily the one that can be passed with the highest average velocity, since the integral is calculated with respect to the density at a fixed time. Thus, the fluctuation of the density while the thief is passing the edge is not captured and the actual passing time is not perfectly reproduced. In reality, this is a reasonable error, because the thief does not know how traffic will change. He simply looks at the outgoing edges and chooses the one which seems to be the least crowded. The main drawback of (DC2) is that the thief  is not determined towards the hideouts. He makes local decisions at each vertex. As a consequence, it might get hard to find one of the hideouts. This behavior is rational in the case of an unprepared, spontaneous crime with the thief not having any knowledge of the neighborhood. (DC3) adds an additional perspective to the thief's decision criterion at junctions. Thereby, the distance to the closest hideout is taken into account. It is obvious that the thief wants to be at a secure place as fast as possible. Thus, it is reasonable to take the distance to possible hideouts into account in order to get the thief out of the system as fast as possible. However, this criterion assumes the thief to have good knowledge of the area where the crime is committed, and to be able to orientate quickly. A possible drawback of this approach is that the thief can be trapped on an edge where the density is at its maximum. This is possible, when the thief is close to one of the destinations and all other options would result in a huge detour.
One could even go further and assume the thief to have information about the densities on all edges in the network. However, this is not evident, since he cannot observe more than those roads that are incident to his current location. Therefore, the only available information about other edges is their length, which makes the thief able to determine shortest paths independently from the traffic volume. An approach on finding a global shortest path to the hideouts with respect to the fluctuating densities is presented below.
The criteria (DC4) and (DC5) are the stochastic counterparts of (DC2) and (DC3). They include a random number u k ∼ U([0, 1]) in order to simulate partly random decisions.
The differences between the two approaches (DC2) and (DC3) presented above are outlined in the following by an example in a small network.
Example. This network consists of three branches that force the thief to a single decision, see Figure 3. The density on the edges is assumed to be constant over the whole edge and there is no interaction between the edges at the junctions, such that the densities are constant over time. This enables us to calculate the travel times on the edges without using (18), since the thief has constant velocity ω = 1 − ρ e . The set of hideout vertices is D = {4, 6, 8}. The thief starts at the first vertex at t = 0. As the density is 0 on the first edge, he can pass the edge with length 1 in one time unit. Consequently, the thief is at n = 2 after a time of 1 and has to decide, which of the edges leading to the vertices 3, 5, or 7 to pass next. (DC2) selects that edge with the lowest average density. As densities are constant over the edges, the values for c DC2 k (t) are directly given by the constant value, see column 4 in Table 1. Consequently, the thief chooses the edge leading to n = 3. The required time for passing that edge is l e /ω = 1/0.7 ≈ 1.43. At n = 3, there is no further decision, as there is only one predecessor. The thief has to take the connection (3 → 4), which yields a travel time of ≈ 12.5. The total travel time is then 13.93, cf. column 6 in Table 1.  Table 1. In contrast to approach (DC2), the thief chooses the edge leading to n = 5 due to the shorter residual distance to the hideout. The travel time on this edge is 1.67. Adding the travel time 2.5 on the edge (5 → 6), a total travel time of 4.17 results. This is far less than the travel time over n = 3 proposed by (DC2).
The connection (2 → 7) is never chosen. However, this route yields the fastest path to a hideout. Its total travel time for the last two edges is ≈ 3.67.
Note that normally densities fluctuate over time due to the formation of shock and rarefaction waves and there is an interaction between the edges at the junctions.

Time-dependent shortest path (TDSP).
This section concludes with a brief introduction of an approach for optimally determining the thief's route. Even though the thief is not able to apply this method in reality as it would require him to have total information, it serves as a performance indicator for the model above. Therefore, a so-called time-dependent shortest path (TDSP) algorithm is used. The most famous TDSP algorithm by Dreyfuß [12] uses the same techniques as the Dijkstra algorithm, but is applicable to time-dependent edge weights. The major restriction in the method of Dreyfuß is the First-In-First-Out (FIFO) rule. First-in-first-out means that someone who enters an edge before someone else, has to arrive at the end of the edge first, supposed that both use the same dynamic edge weights, i.e. overtaking is not included in the model. Thus, it is never preferable to wait at a vertex aiming for a lower edge passing time.
The relevant networks for this work, whose edge weights correspond to the model (1)-(2) do fulfill the FIFO property. Thus, by doing some adaptions to the TDSP algorithm such that the PDE-ODE model for tracking the thief's position can be integrated, the fastest possible route that the thief can possibly take to reach one of the destinations can be determined. Therefore, the edge arrival times a e (t) have to be derived. a e (t) is the time when one arrives the end β(e) of edge e after starting at its initial point α(e) at time t. Consider a network (V, E) and suppose a solution ρ e (t, x) of the conservation law to be given for each edge e and an appropriate time horizon T . Imagine the thief having arrived at a vertex n at time t 0 and having chosen edge e to pass next. Then the edge arrival time a e (t 0 ) for edge e can be determined by solving the IVP (18). From the solution y(t) of this IVP, the arrival time a e (t 0 ) at the end of the edge can be derived. The actual travel time on the edge is a e (t 0 ) − t 0 .
Given that the dynamic edge weights can be determined by solving the ODE (18), Dreyfuß' algorithm can be applied without further restrictions. When the algorithm calls the edge weight function a e (t) at a time t 0 , (18) has to be solved for t ≥ t 0 with initial condition y(t 0 ) = 0.
Algorithm 1 illustrates the procedure of finding an optimal escape path in pseudocode. Besides the density function ρ e (t, x), the set of possible destinations D ⊂ V \ {n s }, and the description of the thief's speed ω(ρ(t, x)), the algorithm requires a starting vertex n s ∈ V and a departure time t s . An earliest arrival time EA v will be the output. Moreover, the algorithm returns the descriptions E, V, and A of the optimal escape path. solve y(t c ) = l k for t c ; if t c ≤ EA β(k) then EA β(k) = t current ; previous(β(k)) = α(k); Select d ∈ D with minimal EA d ; Determine E, V, and A recursively via backward induction using previous and EA.
It can be shown in a straightforward way that Algorithm 1 satisfies the optimality conditions from Lemma 3.1 in [10] and is therefore correct. Examples of optimal escape paths are given in the next section, as the resulting paths can be compared directly to the paths generated by the different decision criteria (DC1)-(DC5).
Step 2. The thief starts his journey at the initial point of a certain road at time t = 0. The thief's position y(t) on a road is determined by solving the ODE (18) with an appropriate solver. When the thief reaches a junction, Definition (3.2) (Decision Criteria) is applied to determine the road to be passed next. The calculation ends when the thief arrives at a hideout vertex or t ≥ T max . Before assessing the thief tracking model itself, the staggered Lax-Friedrichs method is briefly introduced, see [20]. This method is used for the numerical solution of the LWR network model. The benefit of this method is the direct incorporation of the boundary conditions by boundary flux values (γ i (t),γ j (t)) arising from (16) while the corresponding boundary densities are omitted. This is preferable for the coupling of roads, as the boundary densities (ρ i ,ρ j ) do not have to be determined.
The right-hand side of the thief tracking ODE (2) reveals a special difficulty making the solution more complex. Note that y(t) is not given explicitly, but acts as a parameter for the density function ρ(t, x), which is the solution of the LWR model (1). For this, the ODE and its solution depend not only on its own initial value, but also on the setup of the LWR initial value problem. Here, a distinction has to be made with respect to the way how the LWR model is solved. On the one hand, ρ(t, x) could be assumed to be an analytical solution, which can for example be derived for piecewise constant initial data. Such a solution exhibits rarefaction and shock waves, which results in a discontinuous right-hand side of the ODE. On the other hand, acting on the assumption that the LWR model is solved numerically with the staggered Lax-Friedrichs scheme, one faces a piecewise linear function ρ(t, x), since the solution at the grid points is linearly interpolated. In this case, the right-hand side is continuous. Both cases combine the fact that their respective solutions to the LWR model are not generally expressible by a single equation, which could be inserted in the ODE. Consequently, it is hard to derive general statements about stability, convergence, and consistency of numerical solution schemes, as Lipschitz continuity of the right-hand side cannot be stated for general initial data. The only statement that can be done without any restrictions isẏ(t) ∈ [0, 1], since ρ ∈ [0, 1]. As a result, the ODE is bounded and numerical solutions to the ODE are of bounded variation. Despite the lack of theoretical background for the existing numerical ODE solution schemes, they still apply to the underlying problem. As the solvers evaluate the right-hand side of the ODE at discrete points, there is no problem with piecewise linear functions. Moreover, the schemes do not recognize discontinuities in ρ(t, x) for the same reason.

Performance analysis.
A small study is conducted in this paragraph aiming at evaluating error levels and computation times of different ODE solvers being applied on various initial data settings of the LWR model. Analytical and numerical solutions to the LWR model are both used and their impact on the error is compared. We choose a setting, in which a thief walks along an infinite road. The road has Riemann initial densities, i.e. there is a single discontinuity at x = 0 which produces a wave. The thief starts at x = −5 and his trip is observed until he reaches x = 5. The crucial performance indicator for the respective approaches is the difference in the arrival times at x = 5 compared with the analytical arrival time (cf. Section 2.2), since we are particularly interested in the total travel time on an edge. Nevertheless, truncation errors are considered as well. Two different initial settings for the LWR model are chosen. The first setting generates a shock wave (ρ l = 0.1, ρ r = 0.6) while the second one generates a rarefaction wave (ρ l = 0.9, ρ r = 0.5). The discretization for the staggered Lax-Friedrichs scheme is done with step sizes ∆t = 0.05 and ∆x = 0.1.   The results of the study can be viewed in Table 2. To sum up, it can be stated that the choice of a particular solver is not crucial for the problem. The Figures  4 and 6 show plots of several thief trajectories for a shock and a rarefaction wave, respectively. The Figures 5 and 7 show zooms on critical regions pointing out where the numerical solutions differ the most from the analytical solution.
The major insight from the results is that the largest part of the error seems to come from the numerical solution of the LWR model rather than from the ODE  Table 2. Results of numerical study for the thief tracking ODE.    solution. Thus, the accuracy of the staggered Lax-Friedrichs scheme is the crucial factor for an accurate solution. In order to state that, another experiment is conducted. ode23 is applied multiple times to LWR solutions produced by the staggered Lax-Friedrichs scheme, whose mesh size is halved each time. Both ∆x and ∆t must be halved in order to satisfy the CFL condition. We start with ∆x = 0.1 and ∆t = 0.05 like above, and halve the step sizes twice. The results are shown in Table 3.
It can be clearly seen that the quality of the solution to the ODE depends strongly on the mesh size used by the staggered Lax-Friedrichs scheme. The improvement potential of refining the step sizes in the traffic simulation is so remarkable that it makes the improvement potential of better ODE solvers almost irrelevant. However,

Staggered LF Solution
Mesh ρ l = 0.1, ρ r = 0.6 ρ l = 0.9, ρ r = 0.5  Table 3. Furthermore, as the step sizes are mostly chosen such that the CFL condition is just met, both discretizations have to be refined in order to satisfy the CFL condition. This effect is further amplified when dealing with road networks. Hence, there is a trade-off between solution quality and computation time, which is particularly an issue when large road networks are simulated as is done in this work.

4.2.
Computational experiments. Now, the different decision criteria at junctions are addressed numerically. The main issue of (DC2) is the calculation of the integral of the density on the outgoing edges at the time, when the thief makes his decision at the junction. The solution ρ(t, x) of the LWR model produced by the staggered Lax-Friedrichs scheme is solely provided at the discretization points of the time grid T and the spatial grid X . As explained above, the solution is linearly interpolated in order to be able to access ρ-values for arbitrary parameters t and x. Instead of computing the integral for each spatial cell analytically, the trapezoidal rule is used, which is exact for piecewise linear functions. Thus, the following calculation rule is used: According to (DC2), the edge with minimal c DC2 k (t) is chosen by the thief. If there are multiple edges with the same c k (t), the one with the smallest index is chosen. As a result, the solution of the thief's path is guaranteed to be unique.
(DC3) additionally addresses the distances to the closest destination for the assessment of the edges, i.e. for each edge k ∈δ − e,n , the distances k n of the shortest path including k to one of the destinations is additionally considered for the decision. It is necessary to determine s k nd (t) for all destinations d ∈ D in order to find the distance to the closest d. As the shortest path is defined in terms of length and not in terms of travel time, the Dijkstra algorithm can be applied here without restrictions. A problem equivalent to the all sources-all destinations shortest path problem is solved once at the beginning. The results can be used throughout the whole process of path simulation as the edge lengths are static.
In the next step, a study is conducted, which aims at analyzing the performance of the decision criteria in terms of escape time. We aim at giving an illustrative insight on the dynamics of the thief tracking model in networks. For this, the five decision criteria (DC1) − (DC5) are applied to a test network which is illustrated in Figure 8.
The green vertex 19 is the starting vertex, the three red vertices 24, 29, and 34 are hideout vertices. All interior edges are of length 1. Three different simulations If the thief exceeds the time horizon before arriving a hideout, the information from that iteration is not used, since the density beyond the time horizon is not valid. As a consequence, the results are biased because only information from journeys lasting less than 100 time units counts for the results. This issue cannot be resolved having a finite time horizon. However, a time horizon of 100 proved to be sufficiently large to capture all escape paths that do not end up on a completely jammed edge. An overview on the results of those experiments is given in Table 4. The values of the criteria with random component are average values. The second column contains arrays with four entries, which show the distribution of the arrivals at the different hideouts. The fourth entry is the number of times the thief failed to escape within the time horizon. The third, fourth, and fifth column reveal (average) escape time, speed, and length of the escape path.  Table 4. Results of decision criteria performance analysis.
The time-dependent-shortest-path serves as a performance reference. The results are not surprising. First of all, the hideout vertex 29 is the least popular hideout, as it is the farthest hideout from the starting vertex. The other hideouts have equal distance to the starting vertex. Their different popularity results from the different traffic volume in the area around those hideouts. In simulation (a), the different objectives at junctions between the deterministic approaches (DC2) and (DC3) can be observed. (DC2) has a higher speed, but this comes at the cost of a higher distance as well as a longer escape time. The more targeted approach (DC3) reveals shorter escape distances. The crucial indicator is the escape time. (DC3) is able to escape in less than half the time of (DC2). Here, (DC3) even generates the same solution as the time-dependent-shortest-path.
The corresponding approaches with random component verify the behaviour of their deterministic counterparts. It is clear that the performances are worse due to the partly random choices, but it is remarkable that (DC4)'s results are much closer to (DC5) than those of (DC2) compared to (DC3) in all aspects. The random walks produced by (DC1) exhibit the worst performance by far. In 6 of 100 cases, the thief does not arrive at all within the time horizon. Simulation b's results are obviously better due to the initial densities being halved. In each iteration, even those with (DC1), a hideout is arrived within the time horizon and speed increases for all (DCx). Moreover, the escape time decreases for all (DCx) except from (DC2), which is not surprising since (DC2)'s paths do only have a local focus which can results in huge detours. One would assume that the escape distances decrease as well for all (DCx). Thus, it is surprising that (DC1), (DC2), (DC4), and (DC5) have increased average escape distances. In fact, this is a reasonable result, as (DC1), (DC2), and (DC4) do not take the distance into account for the decision at junctions. Moreover, the influence of the residual distance is small for (DC5), when the thief is at the beginning of his path and all outgoing edges have similar residual distances. Consequently, the effect of the density on the outgoing edges is decisive for the edge decision leading to the same argumentation for the increased escape distances than before. Note that in simulation (c), a hideout is not even reached by the time-dependent shortest path within the time horizon. Hence, there is no possibility for any decision criterion to reach the hideout before the time runs out due to the high initial densities which provoke jams on large parts of the roads. Figure 9 and Figure 10 show the paths of (DC2) (blue), (DC3)(green), and the time-dependent-shortest-path (purple) from simulation (a) and (b), respectively. (DC3) and TDSP generate the same path in simulation (a) (green).  To conclude with this study, two characteristic situations are shown up, which point out the dynamics between the traffic and the thief's path in an illustrative way. First, the susceptibility of decision criterion (DC2) for including circles in its paths is investigated. As we know, (DC2) bases its decisions solely on density. Thus, it can happen that the thief walks along a circle if the average densities at the outgoing edges are such that the thief chooses an edge which leads back to an edge he has already passed. If the fluctuation of the density is very small, it may even happen that he passes a circle multiple times. The occurrence of the first circle (19 → 18 → 1 → 20) in (DC2)'s path in simulation (b) (see Fig. 10) is investigated. Vertex 18 is passed twice. The first time, the thief moves to vertex 1 and returns to 18 via 20 and 19. The two arrival times at vertex 18 are t 1 = 1.10 and t 2 = 5.64. Figure 11 and 12 illustrate, why the thief decides differently the second time. At t 1 the pink line representing c DC2 18→1 (t) is far below the green line representing c DC2 18→17 (t). Thus, the thief chooses the connection (18 → 1). While he is on the circle, the average density on (18 → 1) increases, whereas the average density on (18 → 17) is approximately constant. At t 2 the green line is slightly above the pink line, making the thief choose (18 → 17) this time. However, this is a very tight decision, as the difference is very small. If the thief arrived 0.16 time units later at t = 5.8, he would go on a second loop on the circle.  Second, it is shown, how a thief can be trapped on an edge. This happens when a shock wave with negative velocity enters the edge from the right-hand side and increases the density to 1 while the thief passes the edge. We consider an edge of length 1, which has a density of 0.4 in the interval [0, 0.8] and a density of 1 in the interval [0.8, 1] when the thief arrives at the vertex at the left-hand side of the edge. The average density on the edge is 0.8 · 0.4 + 0.2 · 1 = 0.52. If all other feasible edges have higher average densities, the thief will choose that particular edge. While he passes the edge, the shock wave moves towards him from the end of the edge. Thus, he hits the shock wave after some time and cannot move forward. Oftentimes, such a jam does not dissolve within the time horizon. Figure 13 shows the trajectory of a thief passing the edge. He hits the shock wave after a time of 0.8 and is stuck at 0.48 due to the jam.  Usually, situations in which a thief is trapped on an edge do not occur using the TDSP algorithm, whereas it can well happen for all (DCx). However, the network in simulation (c) is so crowded that entering a completely jammed edge can not even be avoided by the TDSP algorithm.

5.
Conclusion. The present work on the topic of crime modeling adds a new perspective to the research in the field of continuous problems coupled with discrete decisions. Therefore, an approach on how to model the escape path of a thief in an urban road network is provided. On the one hand, the determination of the thief's location is a continuous problem since an ODE is introduced describing the velocity of the thief depending on the current density at his location. The solution of this ODE describes the thief's actual location. As the formulation of the ODE contains the solution of a PDE, the whole approach is a PDE-ODE coupled system. On the other hand, the thief has to make discrete decisions on which outgoing edge to pass next when he arrives at a junction. These decisions are mainly based on the traffic density on those edges.
In conclusion, this work provides a concrete approach on how to simulate the path of single road users in a network having the objective to escape. The thief tracking model is a generally applicable method that could be transferred to other scenarios, e.g. the description of emergency vehicles. Future research includes strategies to catch the thief. Göttlich [18,17] introduced an approach on flow-optimal traffic light setting by expressing the LWR model with linear constraints and steering the traffic through binary variables that are used to change the color of the traffic lights. Expressed by mixed integer constraints, the setting of the binary variables can be optimized such that the flow is maximized. A similar approach could be used to place road barriers such that the thief's path is influenced and he can be caught easily. Moreover, the thief versus police setting could be refined by including game theoretical aspects like in the famous game Scotland Yard, compare for instance [13].