Modeling illegal logging in Brazil

Deforestation is a major threat to global environmental wellness, with illegal logging as one of the major causes. Recently, there has been increased effort to model environmental crime, with the goal of assisting law enforcement agencies in deterring these activities. We present a continuous model for illegal logging applicable to arbitrary domains. We model the practice of criminals under influence of law enforcement agencies using tools from multiobjective optimal control theory and consider non-instantaneous logging events and load-dependent travel velocity. We calibrate our model using real deforestation data from the Brazilian rainforest and demonstrate the importance of geographically targeted patrol strategies.


Introduction
Deforestation, and in particular, illegal logging causes some of the most damaging effects to the world's forests. Modeling and quantifying deforestation has become a recent area of study for ecologists, political scientists, and applied mathematicians. Pfaff [29] has validated the correlation between certain parameters and deforestation in tropical regions such as Brazil. The three dominant categories of parameters are identified by Pfaff and other authors [3,21], as accessibility, population and climate. An effort to control deforestation in Brazil, while exploiting timber in a sustainable way, is through allowing legal concessions for industrial timber harvest in public forests [4]. Companies operating in Brazil, under such concessions, fell an average of only one tree per acre instead of clearcutting to allow tree regrowing [37]. However, as is reported in [37], legal timber companies pull out of concession as uncontrolled wildcat loggers invaded the company's land, illegally toppling and stealing trees. The government's failure to detect and punish illegal loggers leads to even more rampant organized crime and more severe deforestation. It is an essential and urgent task to work out effective tactics to combat illegal loggers. Slough and Urpelainen [36] have studied the influence of geographically targeted deterrence on deforestation. Efficient and effective deployment of law enforcement to threatened areas is the best deterrent for these crimes and has been modeled in the continuum setting [1,2,7]. Effective deterrence also can lead to spatial spillovers, as loggers move away from areas with heavy monitoring. Assessing loggers' responses to policies is important when designing an effective system that minimizes deforestation across forested areas. In this work, we build a game theoretic model to predict interactions between the illegal loggers and law enforcement agencies.

Deforestation in the State of Roraima
In this paper, we focus on the PRODES (Amazon Deforestation Monitoring Program) [24] dataset, which is the official dataset the Brazilian Government uses to make annual statistics relating to deforestation. PRODES uses a mixture of computer and human expert analysis to delineate deforestation regions in Brazilian Amazon annually with a minimum patch size of 6.25 hectares (ha) . In particular, we focus on the state of Roraima, the northernmost state in Brazil. We extract annual deforestation events data and tree coverage data from PRODES. In Fig. 1a, we plot the deforestation events from year 2001 to 2015 as well as the transportation system (including highways and waterways) for Roraima, with the observation that many of the deforestation events occur in the vicinity of roads and rivers. The tree coverage data for year 2015 is shown in Fig. 1b where yellow indicates land covered by trees and blue refers to cleared land. There are fifteen municipalities in this state. In this paper, we assume that loggers originate from and return to these fifteen municipalities.

Previous work
The first continuous game theoretic model of deforestation is attributed to Albers [1], who modeled deforestation events in a circular area with radially symmetric benefit and patrol functions. Criminals enter from the boundary of the area and want to maximize their profit where B is the benefit to the attacker, C is the cost of traveling to depth d, Φ is the cumulative patrol function and (1−Φ(d)) represents the probability of not being captured.
Johnson et al. [17] worked on optimal patrol strategies in the framework of Albers' model [1]. Kamra et al. [18] extended the model by removing the assumption that trees are  Albers' [1] model assumes a circular region and radial symmetric functions so that attackers will only move along the radius. b Arnold et al. [2] generalize the model to arbitrary terrain and applied to the Yosemite national park. In both figures, the white area is pristine while the grey area is affected by criminals homogeneously distributed but maintained the circular area. They considered the game between law enforcement and extractors and applied machine learning techniques to find the optimal or near-optimal patrol strategies. All of these works considered a circular region with the assumptions that extractors come from the boundary of the region and move toward the center. The radial symmetry of the region and the functions is a major restriction.
Arnold et al. [2] generalized Albers' model [1] to any closed, simple region in R 2 . The primary tool employed in the model is the level set method [26]. In their model, the cost represents the effort expended by extracting at any point in the protected area, evaluated by the optimal travel time, where the velocity is allowed to depend on terrain data. They model the impact of patrol by including capture probability in the formulation of a heuristic modified velocity. The validity of this model has not been tested against real-world data, but has been modified and improved by Cartee and Vladimirsky [7]. The authors constructed two models based on whether the authorities use ground patrol, where confiscation takes place as perpetrators are detected, or aerial patrol, so that illegal goods are not confiscated until perpetrators exit the protected area.
Meanwhile, models for illegal extraction by discrete methods have been developed by Fang et al. [11,12] and Kar et al. [19,20]. Both Fang et al. [12] and Kar et al. [20] deployed their models in Queen Elizabeth National Park (QENP), Uganda. Fang et al. [12] developed the PAWS algorithm and described the protected region as nodes connected by edges which are natural pathways such as rivers or roads. Kar et al. [20] used machine learning techniques to predict attacks from extractors. One advantage of such discrete methods is that they can easily incorporate realistic concerns such as detailed terrain information, different types of environmental crime (including animal poaching), or different types of patrol teams [22]. However, they have the disadvantage that they do not track the actual movement of the environmental criminals, and results can be difficult to interpret due to the "black box" nature of parameter estimation methods.

Our contribution
Our work builds upon earlier works by Arnold et al. [2] and Cartee et al. [7] where we use optimal-control theory to model and solve the path planning problem faced by illegal loggers as they balance benefit, travel cost and capture risk. We assume the authorities deploy remote patrols, related to model A in [7], where confiscation of illegal goods is delayed. We introduce several significant improvements to arrive at a more realistic model. We consider non-instantaneous logging activities and positive capture risk while logging on site so that loggers can decide the optimal logging time to maximize profit. We also incorporate load-dependent velocity as loggers return from the forest with illegal goods. We work directly with real-world data from Brazil to calibrate the model. As we provide a more transparent interpretation of "pristine area" and evaluation metrics for patrol efficiency, we simulate and conduct a side-by-side comparison of the predicted outcomes of several geographically targeted and data-driven patrol strategies.
The remainder of the paper is organized as follows. Section 2 describes the model formulation, optimization problem solver and numerical method. Section 3 covers experimental results following a detailed description of the experiment setup. We summarize our conclusions in Sect. 4.

Illegal logging model
Our model, based on the work of Albers [1] and Arnold et al. [2], is a more realistic representation of loggers' decision-making process. The model can be applied to an arbitrary domain as in [2], but balances travel time and capture risk in a more judicious manner by appealing to optimal control theory. We also account for logging time, which was ignored by previous models but makes a great difference to loggers' profit as is shown in Sect. 3. In the remaining part of this section, we will first construct the model. The optimal control problem is posed as a static Hamilton-Jacobi equation and solved using a fast-sweeping method [6,27,38,40,42].

Model construction
Given an arbitrary domain Ω ⊂ R 2 , our goal is to construct an expected profit function P(x) : Ω → R from loggers' perspective. We adopt the basic idea from Albers [1], where Here, B is the benefit that describes the value of the timber that loggers will obtain if they are not captured throughout the entire trip. The variable C represents the cost and is measured by the travel cost of both going in and out of the forest. The term 1 − Φ describes the probability of not being captured, which depends on patrollers' detection ability and loggers' trajectories. We present details of these three components in the following paragraphs. We follow the Stackelberg game model and assume loggers have perfect information about patrol.
The benefit B depends on the value and amount of trees that perpetrators decide to log. We assume each location x in the domain Ω has a fixed amount of timber, but with different total value B(x), depending on the category and quality of timber. Departing from previous models where extracting happens instantaneously, we introduce the notion of logging time t log , which is comparable with travel time. We assume that loggers have a constant production rate 1/T , where T is a global constant representing the time to clear all the trees in one location. The actual benefit ignoring existence of patrollers, is then given by t log T B(x). We always assume that t log ≤ T , so the loggers only extract from one spot x in one trip and will return from the forest when the chosen region is cleared.
We assume loggers can be detected when they are logging or on their path back while returning with their illegal goods. We define the capture intensity as ψ : Ω → R, which we assume to be known to loggers and dependent on patrol resources and strategies. In particular, it satisfies a budget constraint modeled as Here, E represents the budget, d(x) is the Euclidean distance from location x to the major highways and μ is an adjustable weight parameter. The term (1 + μd(x)) 2 models a scenario wherein it is more expensive to patrol deeper into the forest. We tested different capture intensity functions in the experiments in Sect. 3, which exhibit interesting patterns. Following the derivation in [7], the probability of not being captured while logging at x after time t log is e −ψ(x)t log . Longer logging time means larger benefit, but also a larger risk of being detected. The probability of not being captured when they are walking back following path X(s) is then given by e − τ 0 ψ(X(s))ds where τ is the travel time. Here we assume the loggers are only detected when they actually have timber in their possession. In this case, they lose all of the benefit but without extra penalty (see discussion at the end of this subsection). The expected benefit by logging at location x for time t log and returning with path X(s) is then The cost-represented by the travel time-is easy to calculate given a path and the associated velocity. We assume that loggers will embark from and return to one of the fifteen municipalities, and use X in and X out to represent the paths to and from the logging location, respectively. We assume that loggers may return to any municipality-not necessarily the one they embarked from-which leads to a different initial value as is discussed further in Sect. 2.2. In our model, we first define the inward velocity field v : Ω → R 2 following the transportation system, so that loggers travel with highest velocity when they are on major highways and more slowly when they are on water or secondary highways. When loggers are off highways or waterways, their velocity is scaled according to terrain slope following Arnold et al. [2]. When they are returning, we assume their velocity is slower because of the loaded cars or boats. In this case, we set the returning velocity v out = v(x)/(1+c(t log /T ) γ ). Here, t log /T measures the amount of trees loggers carry back. The parameters c and γ model the effect of carrying the trees on the speed of motion. The increased travel cost may be another reason that loggers decide to spend less than the maximal logging time.
The previous analysis leads to a more realistic way of calculating profit for logging at position x for time t log and following paths X in , X out . The profit function is where τ in and τ out are the travel times, and α is a dimensional function that converts time to monetary value. As in [7], α may be constant or may change based on location. This can model loggers' preference for certain areas or represent the variance of unit time travel cost. For example, one may expect the unit time cost of traveling on waterways to be smaller than that of traveling on highways. Rational loggers will then try to solve the optimization problem and may go to the spots with positive profit. It is worth mentioning that the optimal path going in vs. out of the forest can be different because of patrols. In reality, a mixture of different patrol methods are deployed in Brazil, including the ground patrol, using boats and motor vehicles, and the remote patrol, using helicopters, planes and drones. As is pointed out in [7], the ground patrol leads to immediate confiscation, and it would be more reasonable for loggers to switch to the minimal-time path thereafter. In this paper, we always assume that the government deploys remote patrols. Since loggers are unaware of being detected, they will choose optimal return paths that balance capture risk and travel time.

Multiobjective approach and Eikonal equations
We adapt the multiobjective optimal control approach of Cartee and Vladimirsky [7] to our model, and describe its use in solving the optimization problem (4). We consider a trajectory X(s) following the dynamicṡ Here, X 0 denotes the set of possible destinations (the fifteen municipalities). The map a is the control plan, taken from the set of valid control functions where R(x) is the minimal cost traveling from X 0 to x. In fact, since our velocity is isotropic, Recall that along the inward path, loggers do not need to worry about the patrol and travel with velocity v(x). Along the outward paths, their velocity will be decreased if they are carrying more timber, and the amount of timber they are carrying is proportional to the logging time. In our model, we set v out = v/(1 + c(t log /T ) γ ). Plugging this into equation (6) yields We can resolve the optimal profit value using a multiobjective control formulation as in [7]. For any λ is the unique viscosity solution [5,9] of the Eikonal equation While we do not expect u λ to be smooth, under mild conditions on v and K , it will be Lipschitz continuous, hence differentiable almost everywhere, and thus the λ-optimal control is uniquely determined for almost every starting point x ∈ Ω \ X 0 . The value functions corresponding to λ-optimal controls are defined by and given u λ , we can resolve u λ 1 and u λ 2 by solving with boundary conditions u 1 (x) = u 2 (x) = 0, for x ∈ X 0 [7,23]. The optimal profit can be calculated by For the inner maximum, is a function of t log , and it is not easy to find the explicit maximum. In practice, we discretize [0, T ] × [0, 1] into finitely many points (t i , λ j ) and simply choose the maximum among these points. We conclude the discussion of our model with a pair of remarks regarding implementation. Remark 1: Choice of X 0 . In our model, we need to solve the Eikonal equation (7) for the minimal travel time, where X 0 is chosen to be the set of municipalities from which the loggers depart. We also need to solve (10) for the λ-optimal value function u λ , where X 0 represents the set of municipalities that loggers transport timber to. As is briefly mentioned in Sect. 2.1, we do not require loggers to return to the same municipality from which they start. From the patrol perspective, it may not be clear which municipality the loggers will choose. Accordingly, we let X 0 be the set of all fifteen municipalities in both (7) and (10). The resulting optimal profit for loggers reflects their freedom to choose their starting and terminal municipalities. In the case that the loggers are required to return to their starting point, one could perform 15 rounds of calculation, each of which takes X 0 to be a singleton corresponding to one municipality in both Eqs. (7) and (10). In all our simulations, we use the former setup, which always gives profit no less than the latter one, and therefore helps the government to prepare for the worst case scenario. Remark 2: Optimal paths. We find the optimal paths from x ∈ Ω to the set X 0 by following the negative gradient directions for the different value functions. That is, to find the minimum cost path (optimal inward path) from x to X 0 in the absence of patrols, we |∇R(x)| . To find the optimal path (optimal outward path) when loggers carry timbers and go outwards, we integrateẋ = −v(x) ∇u λ (x) |∇u λ (x)| , where λ is the optimal value for the extraction point x.
In our model, the optimal logging time is different at different locations. However, we note that the optimal path is the same regardless of the logging time, because the logging time (and the amount of timbers carried) influences the velocity uniformly. Thus, the path will be traversed more slowly when a larger logging time is used, but the spatial location of the path is the same.
Finally, one of the basic assumptions of the model is that velocity is isotropic, meaning that it depends only on position, not on the direction of motion. Because of this, the optimal path between two points can be determined regardless of which is the starting point and which is the ending point. If we chose an anisotropic velocity-for example, if the downstream and upstream velocities on the river were different-this would no longer be true, and we would need to compute in-coming and out-going paths separately, which would require additional PDEs similar to (7) and (10), but formulated with reversed orientation.

Numerical methods
We need to solve two kinds of PDE in our model, namely the standard Eikonal equation (7), (10) and auxiliary PDEs (12). In this paper, the region Ω is the state of Roraima, which is irregularly shaped. We use a uniform Cartesian grid to discretize a rectangular region in R 2 containing Ω. As mentioned in 2.2, we choose X 0 to be the set of all 15 municipalities in the state of Roraima. This applies to all three equations (7), (10) and (12). To mark the boundary of Ω, we set the velocity to zero outside of Ω, which makes it impossible for paths to leave the region.
To approximate the equations, one can apply standard numerical methods for static Hamilton-Jacobi equations [10]. Two of the most popular methods are fast-marching and fast-sweeping schemes. Fast-marching methods are based on the idea of following characteristic flow and updating values at grid nodes monotonically based on the values at neighboring nodes [33][34][35]39]. With the proper choice for the order of node updates, the fast-marching method can approximate the value function at N grid points with the computational cost of O (N log N ). By contrast, the philosophy of fast-sweeping methods is to account for all possible directions of characteristic flow, and sweep through the grid nodes in alternating directions updating values at nodes in a Gauss-Seidel manner. Each sweep captures the correct characteristic flow for some subset of the nodes, and this process is iterated until convergence [6,27,38,40,42].
We opted for the basic fast-sweeping method presented in [27]. While the fast-marching method may be more efficient, the standard fast-sweeping scheme is sufficient for our purposes and is very easy to implement. If efficiency is a concern, and one still prefers fast-sweeping methods, one may be able to parallelize the computation as in [41], though we did not do this. Alternatively, there are some novel fast methods to solve Eikonal equations such as [8], which uses a hybrid fast-marching and fast-sweeping approach.

Implementation and results
In this section, we apply the optimization solver to our model as described in Sect. 2.2. We start with a detailed description of the benefit function, velocity function and evaluation metrics in Sect. 3.1 and follow with an analysis of the numerical results in Sect. 3.2.

Experimental setup
High fidelity inference of the benefit requires domain knowledge of Brazilian forest, and in this paper, we simply construct the benefit function based on the PRODES dataset [24]. We make the assumption that deforestation for agricultural land clearance only takes place within 50 kilometers of the major highways and treat all the other deforestation events as the result of logging, which are marked using red circles in Fig. 3a. We then design the benefit based on a further assumption that high benefit gives rise to high event frequency within the region. Specifically, we use the same technique as in kernel density estimation [28] by assigning a two-dimensional Gaussian to each event. From the PRODES data, we also construct a binary indicator function of tree coverage of the region. We then generate the logging benefit by linearly combining the generated density function and the binary indicator function as is shown in Fig. 3c. Our approach may give a reasonable but not fully accurate evaluation of the true benefit, as features like distance to municipalities and patrols are not incorporated. However, the simpler benefit model allows us to focus on exploring loggers' behavior under the influence of other factors. The inverse problem of recovering the benefit function based on the deforestation event data, travel distance and capture risk may be of interest by itself. Moreover, as is shown in Fig.  3a, many logging events are in the periphery area of the state of Roraima, providing further evidence that these regions have high benefit worthy of long distance travel. In practice, local governments may be able to design a more realistic benefit function by incorporating more granular data involving types of trees and other vegetation in specific areas.
Next, since the logging model is quite sensitive to the transportation system, we design a velocity field, as shown in Fig. 4, to accurately capture the movement of loggers throughout the region. Using the highway and waterway map from OpenStreetMap [25], we assign the velocity 1, 0.7, 0.4 to major highways, secondary highways and waterways, respectively. This reflects the assumption that loggers use trucks and cargo ships to transport timber. Outside of these regions, we use a velocity model based on local slope of the terrain. Specifically, we use elevation data from the Shuttle Radar Topography Mission [13], and set the slope S(x, y) = |∇E(x, y)| where E(x, y) is the elevation map of the region. The velocity is then given as 0.2 times a function of local slope as described by Arnold et al. [2], who based their velocity function on that of Irmischer and Clarke [15]. Note that in reality, the velocity on waterways may be anisotropic, diverging from the isotropic assumption used in our model. One may generalize the model to incorporate the more realistic scenario and arrive at anisotropic Hamilton-Jacobi equations similar to Eqs. (7) and (10), which can be solved similarly.
The model may be helpful to aid in the design or evaluation of geographically targeted patrol strategies. Many governments, including the Brazilian government, endeavor to combat deforestation by designating geographic areas as protected areas or priority areas for monitoring and enforcement. However, identifying where protection should be targeted presents challenges for policymakers. Our model evaluates the efficiency with which targeted patrol strategies can reduce deforestation using three metrics: 1. Pristine area ratio PA: we define the regions with non-positive profit as pristine area.
PA calculates the ratio of the area of pristine region over the area of the state as . 2. Pristine benefit ratio PB: this metric weighs pristine area by benefit as PB = and represents the ratio of benefit within the pristine area over the total benefit. 3. Weighted profit WP: we interpret the positive part of the profit as the probability density for loggers to choose the logging location. We then define WP as the expected is the non-negative profit and is defined as We run the model on a 600 × 600 grid. We discretize λ and t log into 101 levels and set μ = 2 (5 max x∈Ω d(x)) ≈ 7.33 × 10 7 , T = 2,000,000, and set α to be μ on the highways and 0.7μ otherwise.

Results
We test our model with different patrol budgets and patrol strategies. We also explore the influence of logging time and changing the velocity when traveling with goods.

Example 1: No patrol
We first impose no patrol. Recall that the returning velocity is modified to be v(x)/(1 + c(t log /T ) γ ) to account for the influence of carrying timber. When there is no patrol, i.e., ψ(x) = 0, the optimal paths traveling in and traveling out are the same. When additionally the amount of timber has no influence, i.e., c = 0, the loggers will always use the maximal logging time T . These statements are generally not true when patrol is present or c = 0 for the latter case. The resulting nonnegative profit P + is shown in Fig. 5a. We then test c = 0.5, γ = 1 and c = 1, γ = 1.5. As is shown in Fig. 5, larger c gives a harsher penalty to velocity and thus leads to smaller profit. In all of the following experiments, we fix c = 0.5, γ = 1.

Example 2: Comparison of different budgets
Recall that we impose a budget constraint for patrol following Eq. (1). We assume that patrol uses up all of the resources available and hence we impose equality in Eq. (1) for all our simulations. In this example, we set Recall that d(x) is the Euclidean distance to major highways and the above patrol simply puts more effort on regions closer to major highways. The capture intensity is plotted in Fig. 9c, which will be further discussed in Sect. 3.2.4. We then consider values of E set to be 0.001, 0.003 and 0.005. The resulting nonnegative profit is plotted in Fig. 6. As is expected, the higher budget gives lower profit. In all of the following experiments, we fix E to be 0.003.

Example 3: Influence of patrol on logging time
We use the same experimental setup as in the previous example, and we fix E = 0.003. Recall that we discretize the logging time and search for the best logging time by a parameter sweep. In all of the experiments, we use 101 different levels. We plot the optimal logging time in Fig. 7a, where the numerics represent the proportion of the maximal time T used. We then sample four points in this region that achieve optimal logging time when logging for 50%, 60%, 70% and 80% of T , respectively, and are marked as red points in Fig. 7a. Figure 7b shows the profit as a function of logging time at each point and we see that there are different optimal logging times for each point. Note that when c = 0, i.e., when the timber carried has no influence on the travel velocity, the optimal time is only dependent on the capture intensity and is equal to min{1/ψ(x), T}. When c is nonzero, the solution of the optimal logging time is more complicated, as both benefit and travel time now play a role in addition to patrol.

Example 4: Comparison of different patrol strategies
The enforcement strategy in Brazil during 2003-2012 involved combination of satellite and ground patrols. A reduction in deforestation was achieved. Patrols were active in areas with significant deforestation, the so called "priority" municipalities. Patrols were also sent where the satellite system revealed suspiciously high deforestation [16]. In this example, we compare the patrolling efficiency for different capture intensity functions ψ(x), all of which require budget E (fixed as 0.003), calculated based on Eq. (1). We plot the corresponding capture intensity function, profit P + and optimal time for each experiment and summarize the evaluation based on aforementioned metrics in Table 1. First, we consider a patrol only based on distance to roads by setting where r is chosen from the values 1, 5, 15. We focus on regions that are close to roads as logging and patrol costs are low in these regions. Larger r means the patrol is more concentrated near the highways, while smaller r leads to more uniformly distributed patrol. Figures 8, 9 and 10 exhibit the corresponding capture intensity, profit P + and optimal time. When r = 1, the optimal logging time is more uniformly distributed compared with that of larger values of r. In high benefit regions (away from major highways), the optimal logging time and the profit markedly increase with increasing r, as loggers are less likely to be captured away from major roads. As all three profit plots in Figs. 8, 9 and 10 have the pattern that high benefit regions are high profit regions, we are inspired to include benefit in the model. Next, we set where w is chosen to be either 1, 0.5, or 0.2 so that the high benefit regions are targeted. Similar to the previous example, larger w represents more concentrated patrol. The results are shown in Figs. 11, 12 and 13. Clearly, the intense patrol in high profit regions makes Previous experiments inform us that we need to balance benefit and distance when designing a patrol strategy. Here, we set the patrol to be We test this patrol with w = 0.2, r = 1, 5, 15, plotted in Figs. 14, 15 and 16. The statistics in Table 1 confirm that forests are better protected when patrol strategies take both distance and benefit into account (see (7)-(9) compared with (1)-(6) in Table 1). Compare the Fig. 11 a Capture intensity is based on benefit only, following equation (17), where w = 1. b Expected nonnegative profit P + over the entire region. c Optimal logging time plotted on regions with positive benefit Fig. 12 a Capture intensity is based on benefit only, following equation (17), where w = 0.5. b Expected non-negative profit P + over the entire region. c Optimal logging time plotted on regions with positive benefit Fig. 13 a Capture intensity is based on benefit only, following equation (17), where w = 0.2. b Expected non-negative profit P + over the entire region. c Optimal logging time plotted on regions with positive benefit patrol results based on distance only (results (1)-(3)) to the results based on strategy (18) (results (7)-(9)) in Table 1; it is clear that the "optimal" attention we should give to smalldistance regions, reflected by r, may vary based on w, i.e., the attention high benefit regions get. When we ignore the benefit, the weighted profit (WP) metric and the protected benefit (PB) metric indicate that smaller r is better. When benefit is taken into consideration and w = 0.2, the statistics suggest that moderately large concentration along highways is more appropriate. Moreover, the three metrics are not necessarily positively or negatively correlated. For example, both protected area (PA) and weighted profit (WP) of Fig. 16 are larger then those of Figs. 14 and 15 though these two metrics move in opposite directions in many cases. This feature adds to the complexity of finding the "optimal" patrol strategy. Finally, we consider the patrol strategy that puts more effort along one specific waterway compared with all previous strategies, in addition to targeting highways and high benefit regions. The selected waterways are marked in Fig. 17a with blue curves, and are shown to be popular trails included in many optimal paths in the next experiment. We define a new distance functiond(x) which calculates the minimum Euclidean distance to major highways and the selected waterways. We modify the patrol defined in Eq. (18) to obtain  the following strategy: In Fig. 17, we plot the results when w = 0.2, r = 15. As loggers have to reroute to avoid heavy patrol (see Fig. 18d), the profit is less than that of the previous strategies.
All numerical experiments show that both distance and potential benefit are important factors for patrol allocation. For now, we do not have a method to find optimal patrol strategies, but our model can be applied to evaluate and compare different strategies.

Example 5: Optimal paths
Finally, we calculate and compare optimal paths for illegal loggers in the state of Roraima. We randomly sample 500 target locations with probability proportional to the expected profit shown in Fig. 8b and plot the optimal path going to each of these points (Fig.  18a). We also plot the optimal paths returning from those points under different patrol strategies (Fig. 18b-d). As previously discussed, the optimal paths going to the targets are the same regardless of the patrol. To manifest the differences, we plot the paths leaving the target points (blue curves) on top of those going to the target points (red curves), as shown in panels (b-d). Figure 18b shows that most of the optimal paths going in and returning are quite similar under the patrol plotted in Fig. 13. As we increase the variance of capture intensity by focusing more on the northwest corner as is shown in Fig. 11,  Fig. 18 We sample 500 target points in the region and plot the optimal paths going into the forest (in a) and going out of the forest (in (b)-(d)), with the deployed patrol marked in the sub-title. We plot the leaving paths (blue) on top of the entering paths (red) in (b)-(d). Yellow dots mark the 15 municipalities a major change is observed in Fig. 18c where loggers choose a different route to avoid patrollers and return to a different municipality than the one from which they start. Still, we see that many optimal paths that go deeper into the forest in the northwest corner cluster into one trajectory; one reason why this happens is that the capture intensity is much more uniform than the velocity field due to the presence of rivers. The fast travel speed along the river outweighs the risk of being captured. With this in mind, we design the patrol exhibited in Fig. 17, where the waterways that attract loggers are targeted. The corresponding optimal paths plotted in Fig. 18d unsurprisingly demonstrate huge differences, which leads to the decrease of WP and the increase of PB in Table 1 and indicates the importance of a spatially targeted patrol.

Conclusion
We have presented a control theoretic model to predict the practice of illegal loggers, including their travel paths and target locations. We consider logging events that are sufficiently lucrative for the criminals that they are willing to incur some risk of being caught. The criminals balance the risk against the benefit in order to find an optimal logging time. Our model quantifies the intensity of a logging event, as opposed to a binary model where deforestation events clear all the trees from an area. We believe this better approximates reality, since the PRODES dataset [24] indicates the occurrence of past deforestation events in many locations where trees are still present. We detailed the underlying mathematical formalism, including numerical schemes which can be used to simulate the model. Finally, we tested the model with different values for the parameters and made observations comparing different patrol strategies.
We discuss a few directions for future work on this model. First, one of the basic assumptions is that the patrol strategy is known precisely by the loggers. This is a standard game-theoretic simplification, but is likely false in reality. Allowing for imperfect knowledge (perhaps using stochastic effects) may more accurately describe the differential game between the criminals and patrol. Some work of relevance exists on surveillance uncertainty in reach-avoid games [14]. Second, while the model can evaluate a suggested patrol strategy, in its current form it does not resolve the optimal patrol strategy. Designing a model that can resolve the optimal strategy, or even suggest a constructive method for improving a given suboptimal strategy would be a large step forward. Finally, the model described here is static. One could envision a time-series model wherein this is one stage in an on-going game, and the patrol strategy could change at discrete times. Describing this scenario in a realistic manner would likely require some qualitative changes to the model. Studying the long-time behavior could provide additional insight to the expected amount of deforestation over long stretches of time.
Our model is premised on the idea that efficient patrols against deforestation should be spatially targeted, rather than uniformly applied across a territory. This assumption comports with the targeted nature of deforestation enforcement policies used by many countries. However, the most efficient patrols we recover in our experiments suggest more precise spatial targeting of enforcement than those specified by most existing public policies. Such policies typically target administrative units (e.g., municipalities in Brazil) or other large swaths of the forest. There are clear trade-offs between the precise and blunt targeting, including challenges in patrol strategy implementation; communication of control strategies such that logging can be deterred; and political costs of targeting. The tools developed in this article may be used to help researchers and policymakers to study these tradeoffs in order improve the efficacy of deforestation control policy.