Priority Maps for Surveillance and Intervention of Wildfires and other Spreading Processes

TL;DR

This paper introduces a scalable optimization framework for creating priority maps to monitor and intervene in wildfire spread using UAVs, integrating landscape, wind, and fire dynamics for effective path planning.

Contribution

It presents a novel method leveraging positive system properties to generate priority maps that adapt to changing wildfire dynamics for UAV surveillance and intervention.

Findings

Successfully applied to 16 and 1000 node landscapes

Demonstrates integration of fire dynamics, landscape, and wind conditions

Shows how priority maps adapt to system changes

Abstract

Unmanned Aerial Vehicle (UAV) path planning algorithms often assume a knowledge reward function or priority map, indicating the most important areas to visit. In this paper we propose a method to create priority maps for monitoring or intervention of dynamic spreading processes such as wildfires. The presented optimization framework utilizes the properties of positive systems, in particular the separable structure of value (cost-to-go) functions, to provide scalable algorithms for surveillance and intervention. We present results obtained for a 16 and 1000 node example and convey how the priority map responds to changes in the dynamics of the system. The larger example of 1000 nodes, representing a fictional landscape, shows how the method can integrate bushfire spreading dynamics, landscape and wind conditions. Finally, we give an example of combining the proposed method with a travelling salesman problem for UAV path planning for wildfire intervention.

Equations

$$\dot{x}(t)=Ax(t)$$

$$a_{ij}=\{\begin{array}{ll}-{\delta }_i\le 0& \text{if}i=j,\\ {\beta }_{ji}\ge 0& \text{if}i\ne j,(i,j)\in E,\\ 0& \text{otherwise.}\end{array}$$

$$J={\int }_0^{\infty }{e}^{-rt}Cx(t)dt$$

$${\int }_0^{\infty }{e}^{-rt}Cx(t)dt=\sum _{i=1}^n{p}_i{x}_i(0)$$

$$p^T=C(rI-A{)}^{-1}.$$

$$\frac{d}{dt}\left(e^{-rt}{p}^{T}x(t)\right)=-e^{-rt}Cx(t),$$

$$\begin{array}{rl}{e}^{-rt}{p}^T\dot{x}(t)-r{e}^{-rt}{p}^{T}x(t)& =-e^{-rt}Cx(t),\\ p^T(rI-A)x(t)& =Cx(t).\end{array}$$

$$\begin{array}{rl}\text{minimize}& |p{|}_1\\ \text{such that}& p\ge 0,p^{T}A-r{p}^T\le -C\end{array}$$

$$\dot{x}(t)=(A-K)x(t)$$

$$\begin{array}{rl}\text{minimize}& |p{|}_1\\ \text{such that}& p\ge 0,p^{T}A-r{p}^T-p^{T}K\le -C\\ & \sum _{ij}K\le \Gamma ,K\ge 0,K\le A\\ & k_{ij}=0\forall (i,j)\notin E\end{array}$$

$$\begin{array}{rl}\text{minimize}& \sum _i^n{P}_{ii}\\ \text{such that}& P\ge 0,Q\ge 0\\ & \vec{1}(PA-rP-Q)\le -C\\ & \sum _{ij}{P}_0^{-1}Q\le \Gamma ,P_0^{-1}Q\le A\\ & q_{ij}=0\forall (i,j)\notin E\end{array}$$

Full text

I INTRODUCTION

The application of Unmanned Aerial Vehicles (UAVs) for disaster response is an active area of research [1, 2, 3]. In the case of wildfires, a.k.a. bushfires, there is an interest in exploring how UAVs can offer quicker and more efficient situational awareness as well as intervention via waterbombing. Their mobility, aerial view and ability to quickly gather and transmit data offer the potential of making them a central tool in future fire assessment and response [4, 5]. Although many fire propagation models and simulators exist (e.g. [6, 7]) and path planning for spatial monitoring is well advanced (e.g. [8, 9, 10]), the link connecting them is still missing. Other spreading processes such as flood and the spread of disease pose similar challenges [11]. In this paper, we aim to bridge this gap by generating priority maps for robotic surveillance or intervention that depend on the dynamics of the spreading process.

Pioneering work in firefront propagation modelling was done by Rothermel in [12, 13] and almost all available fire models build upon this approach. Fire spread models can be split out in two categories regarding their landscape representation: continuous and discrete, of which discrete models are better in dealing with heterogeneous data. In the discrete landscape category, cellular automata (CA) models, which are based on fire transition probabilities from one cell to another, are most commonly used due to their ease of implementation, low computational cost and suitability for heterogenous conditions [14]. When based on accurate information about vegetation, terrain, and weather conditions, CA models have been shown to give quantitatively accurate predictions when compared to real wildfire spread [15, 16]. However, a downside of CA models is that they are not directly amenable to optimization-based analysis and decision-making.

A closely related process is the outbreak and spread of epidemics, for which CA-type models are usually simplified to ordinary differential equation (ODE) models such as the spreading Susceptible-Infected-Susceptible (SIS) model [17, 18]. A linear approximation of the nonlinear model is given in [19, 20, 21]. As proven by [20] the linear model can be considered an upperbound for the non-linear model and can therefore be seen as a worst-case scenario.

A priority map defines critical areas, or nodes on a graph, from which the the spread of a fire (or flood, or disease) would be particularly catastrophic in terms of risk of life and property. In epidemics, many different methods for ranking nodes exist, where the most common ones are based on graph topology, such as degree [22], reproduction number [21] or eigenvector centrality[23]. However, in [24] is shown that considering topological features alone is in most cases not enough to indicate influential nodes. In wildfire models, it is particularly important to consider spatially-varying spreading rates, due to the effects of vegetation, terrain, and weather. These are usually not considered in epidemic models.

Once the critical nodes are defined, an interesting problem is to look at strategies to influence the system and reduce the impact of the spreading process within a budgeted constraint, i.e. develop a cost-to-go map indicating which nodes to ‘attack’. However the problem of how many and which nodes or links to remove turns out to be respectively NP complete and NP hard [25]. Therefore a more suitable approach is to reduce the impact of a node by tuning and controlling the spreading rate, as demonstrated by [26]. Solutions to the optimal control problem for vaccination and patching for regulating infection levels are presented in [27, 28, 29].

In [30] parameters for controllability of a network are identified, where the main result is based on topological features. Similar in [31], leader selection and its influence on stability and optimal performance in directed graphs is explored. Furthermore, combination drug therapy is studied by the authors and solved by influencing the main diagonal of the dynamic matrix. Another method for sparse resource allocation for linear network spread dynamics by modifying the diagonal is given in [32]. The dominant eigenvalue of the dynamics is minimized to achieve this.

The main idea of this paper is to determine priority maps for both surveillance and intervention scenarios by considering a bushfire as a spreading process on a graph with similar dynamics as other wide modeled network spreading processes such as epidemics and computer viruses, e.g. as given in [28, 32]. Because the spreading rate will always be nonnegative, it falls within the category of positive systems [33]. Such systems allow for control [34, 35] and identification [36] algorithms based on linear programming, with significantly improved scalability compared to general linear systems.

In this paper we present a method to prioritize nodes in a network for different scenarios for positive systems with an application to UAV path planning. The presented method distinguishes itself on two main points from previous research: 1) it is not topological dependent and therefore, can deal with network structures, such as grids, where all nodes have similar degrees and 2) it also considers nondiagonal entries for control resources and is therefore more flexible and suitable for different types of spreading dynamics. We explain this further in Section II, where we give a framework for modeling spreading processes of positive systems on a network and develop a method to obtain priority maps for two different scenarios: the surveillance and intervention problem. The first is explained further in Section III, whereas the latter is evaluated in Section IV. With examples we convey that the proposed method is easy applicable and provides quick, good and solid results.

II PROBLEM AND MODEL FORMULATION

Cellular automata bushfire simulations like [6, 7, 15] use an underlying stochastic model with Markov transition probabilities. Let us consider a graph with $n$ nodes, where each node $i\in\{1,2,...,n\}$ has a state $x_{i}(t)$ associated with it. This state now transitions following a set of rules. If a node $i$ is on fire, $x_{i}(t)=1$, it recovers with probablity $\delta_{i}\Delta t$ to $x_{i}(t+\Delta t)=0$ and the fire spreads to a neighbouring node with probability $\beta_{i}\Delta t$.

Using a mean-field estimation and Kolmogorov forward equations the deterministic model can be approximated from this stochastic model [11, 37]. Here, the main assumption taken is that all the random variables have zero covariance. It is proven in [38, 39] that the obtained approximated probabilities upper bound the actual ones, which has a positive effect on controlling the underlying stochastic process.

The linear model is finally obtained by linearizing this deterministic model around the fire-free equilibrium point ($\mathbf{x}=0$) [37]. This model again upper bounds the nonlinear model [20, 26]. Therefore the probabilities determined can be seen as the worst case scenario as shown in Fig. 1. Here it can be seen that the linear model realistically upperbounds the non linear deterministic one, which upperbounds the stochastic model. Closer to the equilibrium point the linear model more closely approximates the stochastic process whereas over time the error accumulates.

We can now define the system on the graph by the linear differential equation:

[TABLE]

where $x(t)=[x_{1}(t),...,x_{n}(t)]^{T}$ with $t\geq 0$, is the state of the system and the sparse state matrix $A$ is defined by the linearized spreading dynamics

[TABLE]

Here $E$ is the set of edges of the graph and therefore the off-diagonal entries $a_{ij}$ are determined by the spreading rate $\beta_{ij}$ and the adjacency matrix of the graph, whereas the diagonal entries are given by $-\delta_{i}$, the recovery rate. Because all off-diagonal entries $a_{ij}$ are assumed to be nonnegative, $A$ is Metzler and the system is positive, i.e. if $x_{i}(0)\geq 0$ for all $i$, then $x_{i}(t)\geq 0$ for all $t\geq 0$ [33].

The cost function $J$ associated with the system is defined as

[TABLE]

where $C=[c_{1},...,c_{n}]$ is a row vector defining the cost associated with each node $i$, with each $c_{i}\geq 0$. For example, in a wildfire the cost of a populated area burning is much higher than open grassland. The discount rate $r>0$ can be tuned to emphasize near-term cost over long-term cost. In surveillance problems discounting is necessary for the above integral to exist since the spreading process dynamics of interest are generally unstable.

The question now remains how to define the priority $p$ of the nodes in the graph. In this paper we use the property of positive systems that they have linear and hence separable Lyapunov and storage functions, which has previously been used for scalable methods of robust control and identification [35, 34, 36]. We note that this property also extends to certain classes of nonlinear systems [40, 41].

II-A Surveillance Problem

In the surveillance problem, the dynamics (1) are fixed and the aim is to identify the future cost of a fire spreading from each node. This then serves as a priority for surveillance. Using the properties of positive systems, we propose to find a priority vector $p$ such that

[TABLE]

where $c_{i}$ can be seen as the cost of each node burning, and $p_{i}x_{i}(0)$ as the discounted cost-to-go associated with each node $i$. The important point is the separable structure of this cost bound, i.e. it is the sum of terms each dependent on a single node $i$ without cross-terms, and hence these terms can serve as priority indicators for the nodes. In this paper, we focus on using $p_{i}$ as the priority for node $i$, but if state estimates are available then $p_{i}x_{i}(0)$ can also be used.

For any (unstable) matrix $A$, there exists a discount rate $r>0$ sufficiently large that $A-rI$ is Hurwitz-stable, i.e. all eigenvalues have negative real parts. Under this assumption, the node priorities can be calculated via the simple formula:

[TABLE]

The non-negativity of each $p_{i}$ follows from the fact that $rI-A$ is a positive-stable M-matrix, and hence has an elementwise non-negative inverse [33, p. 134], i.e. all elements of the matrix $(rI-A)^{-1}$ are non-negative, as are all elements of $C$ by construction. Then (5) implies (4) since integrating over the time interval $[0,\infty)$ both sides of

[TABLE]

gives (4), while expanding the left-hand-side gives

[TABLE]

which must hold for all $x(t)$ hence implies (5).

While (5) is suitable for calculating surveillance priorities, we also provide an equivalent formulation as a linear program (LP) which is suitable for extension to intervention problems. The equivalent LP is

[TABLE]

To show equivalence, let $\bar{p}$ be the priority vector calculated using (5), then $\bar{p}\geq 0$ and $p^{T}A-rp^{T}=-C$ hence $\bar{p}$ is feasible for the LP, but any other feasible $p$ has $p^{T}A-rp^{T}\leq-C$ and hence (via inverse-positivity) $p^{T}\geq=C(rI-A)^{-1}=p$, so $|p|_{1}=\sum p_{i}\geq\sum\bar{p}_{i}=|\bar{p}_{1}|$ and so $p=\bar{p}$ is optimal for the LP.

II-B Intervention Problem

In the intervention problem, we consider the scenario where we assume that the system dynamics can be altered to obtain a reduced cost-to-go vector $p$ and the new closed loop system is given by (6). The task is to reduce the discounted cost-to-go by finding a sparse set of locations for intervention within a budgeted constraint. Assuming a linear effect of intervention, the closed-loop dynamics become

[TABLE]

where $K$ is the control matrix and restricted to the total resource constraint $\sum k_{ij}\leq\Gamma$, where $\Gamma$ is the maximum budget available. Note that this imposes an $l^{1}$-type constraint on the elements of $K$ and hence encourages sparsity. A key difference here compared to [31] and [32] is that the control resources can be put on the nondiagonal entries of $K$, i.e. $k_{ij}\geq 0$ if $(i,j)\in E$ and $k_{ij}=0$ otherwise (instead of $k_{ij}\geq 0$ if $i=j$). Hence, the spreading rate is influenced instead of the recovery rate.

Extending the surveillance LP to this scenario we have

[TABLE]

However the multiplication of decision variables $p^{T}K$ means this is not a linear program. Therefore following [35] a new variable $q=p^{T}K$ is introduced, resulting in the updated constraint $p^{T}A-rp^{T}-q\leq-C$. However to take into account the restriction imposed upon the control matrix $K$ by the total resource constraint $\Gamma$ and the adjacency matrix, we take $P=\text{diag}(p)$ and $Q=PK$ to preserve the specified structure of K. This implies $\vec{1}P=p^{T}$ and $\vec{1}Q=q$, where $\vec{1}$ is the row vector with all ones of the appropriate dimension and the problem can be reformulated as

[TABLE]

where $\sum_{ij}|P_{0}^{-1}Q|\approx\sum_{ij}|K_{ij}|\leq\Gamma$ and $P_{0}$ is the initial guess of P taken as the solution of the surveillance problem. The optimization problem now consists of minimizing the trace of $P$, while updating $P_{0}$ to the obtained $P$ till the outcome converges, i.e. $P_{0}=P$.

Afterwards, the final control matrix $K=P^{-1}Q$. The optimization can be extended by including nodes as possible control resource options by changing the constraints on $k_{ij}$.

III SURVEILLANCE RESULTS

We look at two examples with different dynamics to evaluate the method presented in Section II-A and its response to changes in the different parameters.

III-A Example 1: 16-node grid

Let us consider a graph with $n=16$ nodes, connected as visualized in Fig. 2(a) with a spreading rate $\beta=0.5$ and recovery rate $\delta=0.2$, where the cost $c_{i}=1$ for $i=16$ and $c_{i}=0.1$ for every other node. This results in a system for which $a_{ij}=-0.2$ if $i=j$, $a_{ij}=0.5$ if $(i,j)\in E$ and $a_{ij}=0$ otherwise. This can be seen as a landscape with a city, which has a high cost of burning down, on node $16$, which is only connected to node $11$. The objective is stated in Section II-A and taking $r=2$ the normalized results are shown in Fig. 2(b).

As expected due to the connectivity, equal spreading and recovery rate and with one high cost node, Fig. 2(b) shows high priorities for the diagonal nodes $6$ and $11$ leading up to the high cost node $16$.

III-B Example 2: 1000 node fictional landscape

A more interesting problem is to look at the effect of different spreading rates $\beta$, which would be a more realistic representation of actual spreading processes. Therefore we consider the fictional landscape given in Fig. 3 consisting of three different vegetation types, a city and water. It can be represented as a network graph with $n=1000$ nodes. The adjacency matrix and therefore, the set of Edges, is based on a 8 node spreading direction grid, i.e., the landscape is seen as a grid where each node is connected to its direct horizontal, vertical and diagonal neighbors. The recovery rate $\delta=0.2$ for all nodes, whereas $\beta$ is generated using data from cellular automata bushfire models presented in [6] and [15], where the bushfire spreading probabilities are determined based on real bushfire observations. For the spreading dynamics the baseline $\beta=0.5$, where a correction factor for the vegetation is taken as $\beta_{veg}=0.4$, $\beta_{veg}=1$ and $\beta_{veg}=1.4$ for respectively desert, grassland and eucalyptus forest. For the city this is $\beta_{veg}=0.5$ and for the unburnable water areas $\beta_{veg}=0$ is taken. Finally, $\beta$ is corrected for spreading between diagonally connected nodes, following [6]. The cost of the city nodes is again $c_{i}=1$, wheras $c_{i}=0.01$ for all other nodes.

The spreading rate is furthermore adjusted for respectively a southwesterly, southeasterly, northeasterly and northwesterly wind with a speed of $V=4$ m/s. Taking a discount rate of $r=5$, this results in Fig. 4, where clearly the effects of wind direction on the node priority map can be seen. In general, the most important nodes are the nodes in the bottom left corner where the eucalyptic with grassland in between leads up to the city. Other possible hazards are the city itself due to its high cost and the smaller patches of eucalyptus closer to the city. However this area is smaller and restricted above by slow spreading desert and water. Therefore it poses less danger in case of fire igniting in the top part of the landscape. For better interpretation of the priority map, simulations were run starting fires in different priority nodes and the results can be seen in Fig. 5.

The cost that is taken for the landscape nodes also has an effect on the priority map as shown in Fig. 6. With a higher cost on the landscape nodes, they are prioritized over the city, because more landscape burning down will now result in a higher total cost.

III-C Discount Rate

In section II we discussed the discount rate $r$ in (3). For Example 2 with $c_{i}=0.1$ for landscape nodes, $c_{i}=1$ for city nodes and a southwesterly wind of $V=4$ m/s, this effect is visualized in Fig. 7, where the minimal discount rate equals $r=4.705$ before $rI-A$ becomes singular. A higher discount rate prioritizes the near future and hence, prioritizes what is happening in the city, i.e. the highest cost area of the landscape, over what possible outbreaks in other parts of the graph could lead to. In the previously and future examples the discount rate is taken close to the minimal discount rate, i.e. the positive real part of the largest eigenvalue of A, while still attributing some importance to the city, to sketch a more interesting and complete image of the situation.

IV INTERVENTION RESULTS

Let us consider the same examples as in Section III, however now the task is not to prioritize nodes, but the optimal control problem of lowering the total system cost by changing the systems dynamics in the most effective way as explained in Section II-B. Compared to Section III, the results do not directly represent the priority of the nodes, but give a priority on which links control resources should be applied to for the minimal cost of an outbreak in the system. The obtained map will therefore, henceforth be called a cost-to-go map to distinguish the difference with the previously obtained priority map.

IV-A Example 1

We consider the graph with $n=16$ nodes from Example 1, displayed in Fig. 2(a), with recovery rate $\delta=0.2$ and equal spreading rate of $\beta=0.5$ for connected nodes. In Fig. 8 the obtained control resource allocation is demonstrated by the color and width of the links. The node colors represent the updated priority or cost-to-go of the nodes, with the same color-scale as the case without intervention, i.e. $\Gamma=0$.

It is clear that increasing $\Gamma$ reduces the cost-to-go of the critical nodes. Furthermore, the allocation of resources is very sparse, which is well-suited to waterbombing or controlled-burning interventions. It can be seen that for $\Gamma=0.5$, applying all control resources to the link connecting diagonal node $11$ with the city node $16$ reduces the overall cost-to-go of the nodes significantly. This can be interpreted as the overall effect or outcome of one or more of these particular nodes being set on fire is reduced, i.e. the total cost will be less. For more resources available, more links are targeted but still in a sparse manner.

The effect of the discount rate $r$ on the distribution of the control resources is similar to the effect on the priority maps, where for long-term reduction in cost it is best to intervene at a central point, but for short-term focused optimisation it is best to intervene close to the city.

Now, taking the diagonal entries of K as nonzero and hence applying control resources to the nodes rather than edges, results in Fig. 9. The allocation of resources is similar to Fig. 8.

IV-B Example 2

To show the intervention problem is also scalable for larger examples, the intervention map for Example 2, with spreading dynamics as described in Section III-B, no wind present, $c_{i}=0.01$, $\Gamma=40$ and $r=11.5$, is displayed in Fig. 10. Again, it can be seen that the map of edges chosen for intervention is very sparse.

The intervention map behaves similar to Example 1 regarding the amount of resources, cost and discount rate. I.e. depending on the amount of resources available more resources are applied to the links connecting the city and to the eucalyptus forest to the southeast. With less resources available this forest is of less priority due to it being connected to the city by only a single link. Finally, taking a larger cost $c_{i}$ on the landscape nodes or a lower discount rate $r$ would result in a higher cost-to-go on the landscape and therefore more resources would be dedicated to the landscape at the expense of the links connecting the city.

The results of applying the interventions suggested in Fig. 10 are demonstrated in Fig. 11. In case the control resources are applied, the city does not burn down and the fire spreads significantly more slowly.

IV-C Computation and Convergence

The surveillance problem can be directly solved by using (5) and solves within 0.04 seconds for all given examples on a standard desktop computer.

For the intervention problem on the other hand, to solve the linear program and objective as stated in Section II-B, multiple iterations are required before $P_{0}=P$. For both Example 1 and Example 2 this happens quickly within 5 to 10 iterations. The problem is solved with MOSEK in Matlab and for Example 1, with $n=16$ nodes, takes on average 0.1 second per iteration, whereas for Example 2, $n=1000$ nodes, each iteration takes up to 35 seconds.

IV-D Application to UAV Path Planning

The main motivation for the developed methods is to obtain priority maps for UAV path planning purposes. Considering different constrained budgets (e.g. flight time, endurance, speed, sensor range) UAVs will have to make a trade-off where to fly if they are used for surveillance or intervention purposes. For the intervention problem, multiple targets have been identified that all have to be visited or hit. Therefore solving the traveling salesman problem is suitable and in Fig. 12 a possible UAV round trip is generated using binary integer programming for the intervention cost-to-go map for Example 2.

For surveillance problems with fuel-limited and turn-rate-limited UAVs it will also be interesting to consider path-planning algorithms such as [42, 10], which we leave as future work.

V CONCLUSIONS

In this paper we presented a method to obtain priority maps for spreading processes which can be used as an input for UAV path planning problems. Two different problems were considered: the surveillance problem and the intervention problem. Future work will include more testing with more realistic propagation models and sensing assumptions, and further integration with robot motion planning algorithms.